Matrices, Jordan Normal Forms, and Spectral ... - Semantic Scholar

Report 3 Downloads 55 Views
Matrices, Jordan Normal Forms, and Spectral Radius Theory∗ Ren´e Thiemann and Akihisa Yamada August 23, 2015

Abstract Matrix interpretations are useful as measure functions in termination proving. In order to use these interpretations also for complexity analysis, the growth rate of matrix powers has to examined. Here, we formalized an important result of spectral radius theory, namely that the growth rate is polynomially bounded if and only if the spectral radius of a matrix is at most one. To formally prove this result we first studied the growth rates of matrices in Jordan normal form, and partially prove the result that every complex matrix has a Jordan normal form: we are restricted to upper-triangular matrices since we have not yet formalized the Schur decomposition. The whole development is based on a new abstract type for matrices, which is also executable by a suitable setup of the code generator. It completely subsumes our former AFP-entry on executable matrices [6], and its main advantage is its close connection to the HMArepresentation which allowed us to easily adapt existing proofs on determinants. All the results have been applied to improve CeTA [7, 1], our certifier to validate termination and complexity proof certificates.

Contents 1 Introduction

3

2 Missing Unsorted

5

3 Missing Ring

6

4 Missing Polynomial

8

5 Missing Permutations ∗

10

Supported by FWF (Austrian Science Fund) project Y757.

1

6 Missing Abstract-Rewriting

13

7 Haskell Serialization for IArrays

13

8 Vectors and Matrices 14 8.1 Vectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 8.2 Matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 8.3 Block Vectors and Matrices . . . . . . . . . . . . . . . . . . . 34 9 Code Generation for Basic Matrix Operations

39

10 Gauss-Jordan Algorithm 42 10.1 Row Operations . . . . . . . . . . . . . . . . . . . . . . . . . . 42 10.2 Gauss-Jordan Elimination . . . . . . . . . . . . . . . . . . . . 45 11 Determinants

50

12 Code Equations for Determinants 12.1 Properties of triangular matrices . . . . . . 12.2 Algorithms for Triangulization . . . . . . . 12.3 Finding Non-Zero Elements . . . . . . . . . 12.4 Determinant Preserving Growth of Triangle 12.5 Recursive Triangulization of Columns . . . 12.6 Triangulization . . . . . . . . . . . . . . . . 12.7 Divisor will not be 0 . . . . . . . . . . . . . 12.8 Determinant Preservation Results . . . . . . 12.9 Determinant Computation . . . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

54 54 55 57 58 59 59 60 61 62

13 Converting Matrices to Strings

62

14 Characteristic Polynomial

62

15 Jordan Normal Form

64

16 Elementary Column Operations

68

17 Computing Jordan Normal Forms from Upper-Triangular Matrices 17.1 Pseudo Code Algorithm . . . . . . . . . . . . . . . . . . . . . 17.2 Real Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . 17.3 Algorithm to extract JNF vector from matrix in JNF . . . . . 17.4 Preservation of Dimensions . . . . . . . . . . . . . . . . . . . 17.5 Properties of Auxiliary Algorithms . . . . . . . . . . . . . . . 17.6 Proving Similarity . . . . . . . . . . . . . . . . . . . . . . . . 17.7 Invariants for Proving that Result is in JNF . . . . . . . . . .

71 71 72 75 75 77 78 79

2

17.8 Alternative Characterization of identify-blocks in Presence of ev-blocks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17.9 Proving the Invariants . . . . . . . . . . . . . . . . . . . . . . 17.10A Combination with JNF Represented as Vectors . . . . . . . 17.11Application for Complexity . . . . . . . . . . . . . . . . . . .

83 83 86 87

18 Code Equations for All Algorithms

87

19 Comparison of Matrices

88

20 Derivation Bounds

99

21 Matrix Conversions

100

22 Complexity Carrier

102

23 Converting Arctic Numbers to Strings

103

24 Application: Complexity of Matrix Orderings 104 24.1 Locales for Carriers of Matrix Interpretations and Polynomial Orders . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 24.2 The Integers as Carrier . . . . . . . . . . . . . . . . . . . . . . 105 24.3 The Rational and Real Numbers as Carrier . . . . . . . . . . 105 24.4 The Arctic Numbers as Carrier . . . . . . . . . . . . . . . . . 106 24.5 Estimations of Matrix Powers . . . . . . . . . . . . . . . . . . 106

1

Introduction

The spectral radius of a square, complex valued matrix A is defined as the largest norm of some eigenvalue c with eigenvector v. It is a central notion to estimate how the values in An for increasing n. If the spectral radius is larger than 1, clearly the values grow exponentially, since then An · v = cn · v becomes exponentially large. The other results, namely that the values in An are bounded by a constant, if the spectral radius is smaller than 1, and that there is a polynomial bound if the spectral radius is exactly 1 are only immediate for matrices which have an eigenbasis, a precondition which is not satisfied by every matrix. However, these results are derivable via Jordan normal forms (JNFs): If J is a JNF of A, then the growth rates of An and J n are related by a constant as A and J are similar matrices. And for the values in J n there is a closed formula which gives the desired complexity bounds. To be more precise, the values in J n are bounded by O(|c|n · nk−1 ) where k is the size of the largest block of an eigenvalue c which has maximal norm w.r.t. the set of all eigenvalues. And since every complex matrix has a JNF, we can 3

derive the polynomial (resp. constant bounds), if the spectral radius is 1 (resp. smaller than 1). These results are already applied in current complexity tools, and the motivation of this development was to extend our certifier CeTA to be able to validate corresponding complexity proofs. To this end, we formalized the following main results: • an algorithm to compute the characteristic polynomial, since the eigenvalues are exactly the roots of this polynomial; • the complexity bounds for JNFs; and • an algorithm which computes JNFs for every upper-triangular matrix by means of elementary row- and column-transformation: it can be used as such, and also serves as part of the statement that every complex matrix has a JNF As future work it remains to formalize Schur-decomposition which shows that every complex matrix has a similar upper-triangular matrix. Since CeTA is generated from Isabelle/HOL via code-generation, all the algorithms and results need to be available at code-generation time. Especially there is no possibility to create types on the fly which are chosen to fit the matrix dimensions of the input. To this end, we cannot use the matrix-representation of HOL multivariate analysis (HMA). Instead, we provide a new matrix library which is based on HOL-algebra with its explicit carriers. In contrast to our earlier development [6], we do not immediately formalize everything as lists of lists, but use a more mathematical notion as triples of the form (dimension, dimension, characteristicfunction). This makes reasoning very similar to HMA, and a suitable implementation type can be chosen afterwards: we provide one via immutable arrays (we use IArray’s from the HOL library), but one can also think of an implementation for sparse matrices, etc. Even the infinite carrier itself is executable where we rely upon Lochbihler’s container framework [4] to have different set representations at the same time. As a consequence of not using HMA, we could not directly reuse existing algorithms which have been formalized for this representation. For instance, we formalized our own version of Gauss-Jordan elimination which is not very different to the one of Divas´ on and Aransay in [2]: both define row-echelon form and apply elementary row transformations. Whereas Gauss-Jordan elimination has been developed from scratch as a case-study to see how suitable our matrix representation is, in other cases we often just copied and adjusted existing proofs from HMA. For instance, most of the library for determinants has been copied from the Isabelle distribution and adapted to our matrix representation.

4

As a result of our formalization, CeTA is now able to check polynomial bounds for matrix interpretations [3], provided the interpretations are restricted to upper-triangular matrices: it applies the algorithm to compute JNFs, checks that the spectral radius is at most 1, and determines the size of largest Jordan blocks of an eigenvalue with norm 1.

2

Missing Unsorted

This theory contains several lemmas which might be of interest to the Isabelle distribution. For instance, we prove that bn · nk is bounded by a constant whenever 0 < b < 1. theory Missing-Unsorted imports Archimedean-Field Real Rat begin lemma bernoulli-inequality: assumes x : −1 ≤ (x :: 0a :: linordered-field ) shows 1 + of-nat n ∗ x ≤ (1 + x ) ˆ n hproof i context fixes b :: 0a :: archimedean-field assumes b: 0 < b b < 1 begin private lemma pow-one: b ˆ x ≤ 1 hproof i lemma pow-zero: 0 < b ˆ x hproof i lemma exp-tends-to-zero: assumes c: c > 0 shows ∃ x . b ˆ x ≤ c hproof i lemma linear-exp-bound : ∃ p. ∀ x . b ˆ x ∗ of-nat x ≤ p hproof i lemma poly-exp-bound : ∃ p. ∀ x . b ˆ x ∗ of-nat x ˆ deg ≤ p hproof i end lemma listprod-replicate[simp]: listprod (replicate n a) = a ˆ n hproof i lemma set-upt-Suc: {0 ..< Suc i } = insert i {0 ..< i } hproof i Q lemma setprod-pow [simp]: ( i = 0 .. 1 int-mono hproof i

24.3

The Rational and Real Numbers as Carrier

definition delta-bound :: 0a :: floor-ceiling ⇒ 0a ⇒ nat where delta-bound d x = nat (ceiling (x ∗ of-int (ceiling (1 / d )))) lemma delta-complexity: assumes d0 : d > 0 and d1 : d ≤ def shows mono-matrix-carrier (delta-gt d ) def (delta-bound d ) delta-mono hproof i

105

lemma delta-weak-complexity-carrier : assumes d0 : def > 0 shows weak-complexity-linear-poly-order-carrier op > def delta-mono hproof i

24.4

The Arctic Numbers as Carrier

lemma arctic-delta-weak-carrier : weak-SN-both-mono-ordered-semiring-1 weak-gt-arctic-delta 1 pos-arctic-delta hproof i lemma arctic-weak-carrier : weak-SN-both-mono-ordered-semiring-1 op > 1 pos-arctic hproof i

24.5

Estimations of Matrix Powers

We connect the strict ordering mat-gt with norm-bound where the latter can be approximated via compute-jnf-complexity and counting-ones-complexity. Since compute-jnf-complexity provides sharper bounds, we only use this one. definition mat-estimate-complexity :: nat ⇒ 0a :: large-real-ordered-semiring-1 mat ⇒ nat option where mat-estimate-complexity n A = (let B = mat IR A in (if A ∈ carrier m n n ∧ upper-triangular A then let jnf = triangular-to-jnf-vector B in if (∀ na ∈ set jnf . norm (snd na) ≤ 1 ) then Some (max-list (map fst (filter (λ na. norm (snd na) = 1 ) jnf )) − 1 ) else None else None)) lemma mat-estimate-complexity-norm-bound : assumes ∗: mat-estimate-complexity n A = Some d shows ∃ c1 c2 . ∀ k . norm-bound (mat IR A ˆ m k ) (c1 + c2 ∗ of-nat k ˆ d ) hproof i context fixes A B C :: 0a :: large-real-ordered-semiring-1 mat and n d :: nat assumes ∗: mat-estimate-complexity n A = Some d and A: A ∈ carrier m n n and B : B ∈ carrier m n n and C : C ∈ carrier m n n begin lemma mat-estimate-complexity-norm-bound-prod : ∃ c. ∀ k . k > 0 −→ norm-bound (mat IR B ⊗m (mat IR A ˆ m k ⊗m mat IR C )) (c ∗ of-nat k ˆ d ) hproof i

This is the main result for real valued matrices. lemma mat-estimate-complexity-norm-mat-sum-prod : ∃ c. ∀ k . k > 0 −→ norm (mat-sum (mat IR B ⊗m (mat IR A ˆ m k ⊗m mat IR C ))) ≤ (c ∗ of-nat k ˆ d ) hproof i

106

And via conversion, it also holds for matrices over the intended semiring. lemma mat-estimate-complexity-mat-sum-prod : ∃ c. ∀ k . k > 0 −→ mat-sum (B ⊗m (A ˆ m k ⊗m C )) ≤ (c ∗ of-nat k ˆ d ) hproof i end

end

References [1] M. Avanzini, C. Sternagel, and R. Thiemann. Certification of complexity proofs using CeTA. In Proc. RTA 2015, LIPIcs 36, pages 23–39, 2015. [2] J. Divas´ on and J. Aransay. Gauss-jordan algorithm and its applications. Archive of Formal Proofs, Sept. 2014. http://afp.sf.net/entries/Gauss Jordan.shtml, Formal proof development. [3] J. Endrullis, J. Waldmann, and H. Zantema. Matrix Interpretations for Proving Termination of Term Rewriting. Journal of Automated Reasoning, 40(2-3):195–220, 2008. [4] A. Lochbihler. Light-weight containers. Archive of Formal Proofs, Apr. 2013. http://afp.sf.net/entries/Containers.shtml, Formal proof development. [5] R. Piziak and P. L. Odell. Matrix theory: from generalized inverses to Jordan form. CRC Press, 2007. [6] C. Sternagel and R. Thiemann. Executable matrix operations on matrices of arbitrary dimensions. Archive of Formal Proofs, June 2010. http://afp.sf.net/entries/Matrix.shtml, Formal proof development. [7] R. Thiemann and C. Sternagel. Certification of termination proofs using CeTA. In Proc. TPHOLs’09, LNCS 5674, pages 452–468, 2009.

107