Linear Algebra 1: Computing canonical forms in exact linear algebra - PowerPoint PPT Presentation

Algorithms: main types Three ways to group operations: 1. simple iterative ◮ Apply the standard Gaussian elimination in dimension n ◮ Main loop for i=1 to n 2. block algorithms 2.1 block iterative (Tile) ◮ Apply Gaussian elimination in dimension n / k over blocks of size k ◮ Main loop: for i=1 to n/k 2.2 block recursive ◮ Apply Gaussian elimination in dimension 2 recursively on blocks of size n / 2 i ◮ Main loop: for i=1 to 2

Type of algorithms Data locality: prefer block algorithms ◮ cache aware: block iterative ◮ cache oblivious: block recursive Base case efficieny: simple iterative Asymptotic time complexity: block recursive Parallelization: block iterative

Block recursive gaussian elimination Author Year Computation Requirement Strassen 69 Inverse gen. rank prof. Bunch, Hopcroft 74 LUP gen. row rank prof Ibarra, Moran, Hui 82 LSP , LQUP none Schönage, Keller-Gerig 85 StepForm none Storjohann 00 Echelon, RedEch none here 11 CUP ,PLE,PLUQ none

Block recursive gaussian elimination Author Year Computation Requirement Strassen 69 Inverse gen. rank prof. Bunch, Hopcroft 74 LUP gen. row rank prof Ibarra, Moran, Hui 82 LSP , LQUP none Schönage, Keller-Gerig 85 StepForm none Storjohann 00 Echelon, RedEch none here 11 CUP ,PLE,PLUQ none Comparison according to: ◮ No requirement on the input matrix

Block recursive gaussian elimination Author Year Computation Requirement Strassen 69 Inverse gen. rank prof. Bunch, Hopcroft 74 LUP gen. row rank prof Ibarra, Moran, Hui 82 LSP , LQUP none Schönage, Keller-Gerig 85 StepForm none Storjohann 00 Echelon, RedEch none here 11 CUP ,PLE,PLUQ none Comparison according to: ◮ No requirement on the input matrix ◮ Rank sensitive complexity

Block recursive gaussian elimination Author Year Computation Requirement Strassen 69 Inverse gen. rank prof. Bunch, Hopcroft 74 LUP gen. row rank prof Ibarra, Moran, Hui 82 LSP , LQUP none Schönage, Keller-Gerig 85 StepForm none Storjohann 00 Echelon, RedEch none here 11 CUP ,PLE,PLUQ none Comparison according to: ◮ No requirement on the input matrix ◮ Rank sensitive complexity ◮ Memory allocations ◮ Constant factor in the time complexity

Memory requirements: Definition In place = output overrides the input and computation does not need extra memory (considering Matrix multiplication C ← C + AB as a black box) Remark: a unit lower triangular and an upper triangular matrix can be stored on the same m × n storage!

Preliminaries TRSM: TRiangular Solve with Matrix � − 1 � D � A � � A − 1 � � I � � I � � D � − B B = C − 1 C E I I E

Preliminaries TRSM: TRiangular Solve with Matrix � − 1 � D � A � � A − 1 � � I � � I � � D � − B B = C − 1 C E I I E Compute F = C − 1 E (Recursive call) Compute G = D − BF (MM) Compute H = A − 1 G (Recursive call) � H � Return F A B D C E

Preliminaries TRSM: TRiangular Solve with Matrix � − 1 � D � A − 1 � A � � � I � � I � � D � B − B = C − 1 C E I I E Compute F = C − 1 E (Recursive call) Compute G = D − BF (MM) Compute H = A − 1 G (Recursive call) � � H Return F A B D C F

Preliminaries TRSM: TRiangular Solve with Matrix � − 1 � D � A − 1 � A � � � I � � I � � D � B − B = C − 1 C E I I E Compute F = C − 1 E (Recursive call) Compute G = D − BF (MM) Compute H = A − 1 G (Recursive call) � � H Return F A B G C F

Preliminaries TRSM: TRiangular Solve with Matrix � − 1 � D � A − 1 � A � � � I � � I � � D � B − B = C − 1 C E I I E Compute F = C − 1 E (Recursive call) Compute G = D − BF (MM) Compute H = A − 1 G (Recursive call) � � H Return F A B H C F

Preliminaries TRSM: TRiangular Solve with Matrix � − 1 � D � A − 1 � A � � � I � � I � � D � B − B = C − 1 C E I I E Compute F = C − 1 E (Recursive call) Compute G = D − BF (MM) Compute H = A − 1 G (Recursive call) � � H Return F A B H ◮ O ( n ω ) C ◮ In place F

The CUP decomposition 1. Split A Row-wise A1 A2

The CUP decomposition 1. Split A Row-wise 2. Recursive call on A 1 U V 1 C 1 A A 21 22

The CUP decomposition 1. Split A Row-wise 2. Recursive call on A 1 U V 3. G ← A 21 U − 1 1 ( trsm ) C 1 1 A G 22

The CUP decomposition 1. Split A Row-wise 2. Recursive call on A 1 U V 3. G ← A 21 U − 1 1 ( trsm ) C 1 1 4. H ← A 22 − G × V ( MM ) H G

The CUP decomposition 1. Split A Row-wise 2. Recursive call on A 1 U 1 3. G ← A 21 U − 1 ( trsm ) C 1 1 4. H ← A 22 − G × V ( MM ) U G C 2 2 5. Recursive call on H

The CUP decomposition 1. Split A Row-wise 2. Recursive call on A 1 U 1 3. G ← A 21 U − 1 ( trsm ) C 1 U 1 2 4. H ← A 22 − G × V ( MM ) G C 2 5. Recursive call on H 6. Row permutations

Memory: LSP vs LQUP vs PLUQ vs CUP Decomposition In place LSP N LQUP N PLUQ Y CUP Y

Echelon forms From CUP to ColumnEchelon form U = C P A � − 1 � U 1 U 2 P T = Y U I n − r = C P A Id U − 1 − U − 1 � � 1 U 2 P T = 1 I n − r A Y = C Additional operations: − U − 1 U 2 trsm (triangular system solve) in-place

Echelon forms From CUP to ColumnEchelon form U = C P A � − 1 � U 1 U 2 P T = Y U I n − r = C P A Id U − 1 − U − 1 � � 1 U 2 P T = 1 I n − r A Y = C Additional operations: − U − 1 U 2 trsm (triangular system solve) in-place U − 1 1 : trtri (triangular inverse)

Echelon forms From CUP to ColumnEchelon form U = C P A � − 1 � U 1 U 2 P T = Y U I n − r = C P A Id U − 1 − U − 1 � � 1 U 2 P T = 1 I n − r A Y = C Additional operations: − U − 1 U 2 trsm (triangular system solve) in-place U − 1 1 : trtri (triangular inverse) in-place

From CUP to Column Echelon form TRTRI : triangular inverse � − 1 � U − 1 − U − 1 1 U 2 U − 1 � U 1 � U 2 = 1 3 U − 1 U 3 3 1: if n = 1 then U ← U − 1 2: 3: else U 2 ← U − 1 3 U 2 4: TRSM U 2 ← − U 2 U − 1 5: TRSM 3 U 1 ← U − 1 6: TRTRI 1 U 3 ← U − 1 7: TRTRI 3 8: end if

Reduced Echelon forms From Col. Echelon form to Reduced Col. Echelon form A Y = C M A Y = Q � − 1 � M Z = Y M Id I n − r A Y = Q Id C A Z = Similarly, from PLE to RowEchelon form Again reduces to: U − 1 X : TRSM , in-place U − 1 : TRTRI , in-place

Reduced Echelon forms From Col. Echelon form to Reduced Col. Echelon form A Y = C M A Y = Q � − 1 � M Z = Y M Id I n − r A Y = Q Id C A Z = Similarly, from PLE to RowEchelon form Again reduces to: U − 1 X : TRSM , in-place U − 1 : TRTRI , in-place UL : TRTRM ,

Reduced Echelon forms From Col. Echelon form to Reduced Col. Echelon form A Y = C M A Y = Q � − 1 � M Z = Y M Id I n − r A Y = Q Id C A Z = Similarly, from PLE to RowEchelon form Again reduces to: U − 1 X : TRSM , in-place U − 1 : TRTRI , in-place UL : TRTRM ,

Reduced Echelon forms From Col. Echelon form to Reduced Col. Echelon form A Y = C M A Y = Q � − 1 � M Z = Y M Id I n − r A Y = Q Id C A Z = Similarly, from PLE to RowEchelon form Again reduces to: U − 1 X : TRSM , in-place U − 1 : TRTRI , in-place UL : TRTRM , in-place

From Echelon to Reduced Echelon TRTRM : triangular triangular multiplication � U 1 � � L 1 � � U 1 L 1 + U 2 L 2 � U 2 U 2 L 3 = U 3 L 2 L 3 U 3 L 2 U 3 L 3 1: X 1 ← U 1 L 1 TRTRM 2: X 1 ← X 1 + U 2 L 2 MM 3: X 2 ← U 2 L 3 TRMM 4: X 3 ← U 3 L 2 TRMM 5: X 4 ← U 3 L 3 TRTRM

From Echelon to Reduced Echelon TRTRM : triangular triangular multiplication � U 1 � � L 1 � � U 1 L 1 + U 2 L 2 � U 2 U 2 L 3 = U 3 L 2 L 3 U 3 L 2 U 3 L 3 1: X 1 ← U 1 L 1 TRTRM 2: X 1 ← X 1 + U 2 L 2 MM 3: X 2 ← U 2 L 3 TRMM 4: X 3 ← U 3 L 2 TRMM 5: X 4 ← U 3 L 3 TRTRM U U L 1 1 2 U L L 3 2 3

From Echelon to Reduced Echelon TRTRM : triangular triangular multiplication � U 1 � � L 1 � � U 1 L 1 + U 2 L 2 � U 2 U 2 L 3 = U 3 L 2 L 3 U 3 L 2 U 3 L 3 1: X 1 ← U 1 L 1 TRTRM 2: X 1 ← X 1 + U 2 L 2 MM 3: X 2 ← U 2 L 3 TRMM 4: X 3 ← U 3 L 2 TRMM 5: X 4 ← U 3 L 3 TRTRM X 1 U 2 U L L 3 2 3

From Echelon to Reduced Echelon TRTRM : triangular triangular multiplication � U 1 � � L 1 � � U 1 L 1 + U 2 L 2 � U 2 U 2 L 3 = U 3 L 2 L 3 U 3 L 2 U 3 L 3 1: X 1 ← U 1 L 1 TRTRM 2: X 1 ← X 1 + U 2 L 2 MM 3: X 2 ← U 2 L 3 TRMM 4: X 3 ← U 3 L 2 TRMM 5: X 4 ← U 3 L 3 TRTRM X 1 U 2 U L L 3 2 3

From Echelon to Reduced Echelon TRTRM : triangular triangular multiplication � U 1 � � L 1 � � U 1 L 1 + U 2 L 2 � U 2 U 2 L 3 = U 3 L 2 L 3 U 3 L 2 U 3 L 3 1: X 1 ← U 1 L 1 TRTRM 2: X 1 ← X 1 + U 2 L 2 MM 3: X 2 ← U 2 L 3 TRMM 4: X 3 ← U 3 L 2 TRMM 5: X 4 ← U 3 L 3 TRTRM X 1 X 2 U L L 3 2 3

From Echelon to Reduced Echelon TRTRM : triangular triangular multiplication � U 1 � � L 1 � � U 1 L 1 + U 2 L 2 � U 2 U 2 L 3 = U 3 L 2 L 3 U 3 L 2 U 3 L 3 1: X 1 ← U 1 L 1 TRTRM 2: X 1 ← X 1 + U 2 L 2 MM 3: X 2 ← U 2 L 3 TRMM 4: X 3 ← U 3 L 2 TRMM 5: X 4 ← U 3 L 3 TRTRM X 1 X 2 U X L 3 3 3

From Echelon to Reduced Echelon TRTRM : triangular triangular multiplication � U 1 � � L 1 � � U 1 L 1 + U 2 L 2 � U 2 U 2 L 3 = U 3 L 2 L 3 U 3 L 2 U 3 L 3 1: X 1 ← U 1 L 1 TRTRM 2: X 1 ← X 1 + U 2 L 2 MM 3: X 2 ← U 2 L 3 TRMM 4: X 3 ← U 3 L 2 TRMM 5: X 4 ← U 3 L 3 TRTRM X 1 X ◮ O ( n ω ) 2 ◮ In place X 4 X 3

Example: in place matrix inversion A

Example: in place matrix inversion U A L CUP decomposition A = LU

Example: in place matrix inversion −1 U U A L L CUP decomposition Echelon AU − 1 = L

Example: in place matrix inversion −1 U U −1 A A L L CUP decomposition Echelon Reduced Echelon A ( U − 1 L − 1 ) = I

Experiments ./test-invert 65521 A1000 1 518,996,125,000 bytes x ms bytes 16M x804AA7E:_Z10read_fieldI 14M x8049A98:main 12M 10M x804EDC0:void%FFLAS::Win 8M 6M x804EDEE:void%FFLAS::Win 4M x804EDD4:void%FFLAS::Win 2M 0M 0.0 5000.0 10000.0 15000.0 20000.0 25000.0 30000.0 ms ./test-redechelon 65521 A1000 1 280,663,687,500 bytes x ms bytes 8M x804A47E:_Z10read_fieldI 6M x804B0EC:void%FFLAS::Win 4M x804B11A:void%FFLAS::Win 2M x804B100:void%FFLAS::Win 0M 1.0 5001.0 10001.0 15001.0 20001.0 25001.0 30001.0 ms

Direct computation of the Reduced Echelon form ◮ Strassen 69: inverse of generic matrices ◮ Storjohann 00: Gauss-Jordan generalization for any rank profile Matrix Inversion [Strassen 69] � − 1 � A � A − 1 � � I � � I � � � B − B I = ( D − CA − 1 B ) − 1 CA − 1 C D I I I

Direct computation of the Reduced Echelon form ◮ Strassen 69: inverse of generic matrices ◮ Storjohann 00: Gauss-Jordan generalization for any rank profile Matrix Inversion [Strassen 69] � − 1 � A � A − 1 � � I � � I � � � B − B I = ( D − CA − 1 B ) − 1 CA − 1 C D I I I 1: Compute E = A − 1 (Recursive call) 2: Compute F = D − CEB (MM) 3: Compute G = F − 1 (Recursive call) 4: Compute H = − EB (MM) 5: Compute J = HG (MM) 6: Compute K = CE (MM) 7: Compute L = E + JK (MM) 8: Compute M = GK (MM) � E � J 9: Return M G

Strassen-Storjohann’s Gauss-Jordan elimination Problem Needs to perform operations of the form A ← AB ⇒ not doable in place by a usual matrix multiplication algorithm

Strassen-Storjohann’s Gauss-Jordan elimination Problem Needs to perform operations of the form A ← AB ⇒ not doable in place by a usual matrix multiplication algorithm Workaround [Storjohann]: 1. Decompose B = LU LU 2. A ← AL trmm 3. A ← AU trmm

Rank sensitive time complexity Fact Algorithms LSP , CUP , LQUP , PLUQ, ... have a rank sensitive mnr ω − 2 � � computation time: O n=3000, PIII−1.6Ghz, 512Mb RAM 140 Block algorithm Iterative algorithm 120 100 Time (s) 80 60 40 20 0 0 500 1000 1500 2000 2500 3000 Rank

Time complexity: comparing constants O ( n ω ) = C ω n 3 Algorithm Constant C ω in-place C 3 C log 2 7 C ω 2 6 × MM C ω 1 4 ∨ TRSM 2 ω − 1 − 2 C ω 1 8 3 ≈ 0 . 33 5 = 1 . 6 ∨ TRTRI ( 2 ω − 1 − 2 )( 2 ω − 1 − 1 ) TRTRM , CUP C ω C ω 2 14 2 ω − 1 − 2 − 3 ≈ 0 . 66 5 = 2 . 8 ∨ PLUQ LQUP , 2 ω − 2 C ω 3 C ω 22 2 ω − 2 − 1 − 1 5 ≈ 4 . 4 ∨ Echelon 2 ω − 2 C ω ( 2 ω − 1 + 2 ) 44 2 5 = 8 . 8 ∨ RedEchelon ( 2 ω − 1 − 2 )( 2 ω − 1 − 1 ) 5 C ω C ω 76 2 ω − 1 − 1 + 4 5 = 15 . 2 × StepForm ( 2 ω − 1 − 1 )( 2 ω − 2 − 1 ) C ω GJ ∗ 2 8 × 2 ω − 2 − 1 ∗ : GJ : GaussJordan alg of [Storjohann00] computing the reduced echelon form

Applications to standard linalg problems Problem Using C ω C 3 C log 2 7 In place Rank RankProfile C ω 2 8 × GJ 2 ω − 2 − 1 IsSingular C ω C ω 0 . 66 2 . 8 ∨ 2 ω − 1 − 2 − CUP 2 ω − 2 Det Solve C ω 2 × GJ 8 2 ω − 2 − 1 Inverse C ω ( 2 ω − 1 + 2 ) 2 8 . 8 ∨ CUP ( 2 ω − 1 − 2 )( 2 ω − 1 − 1 )

Summary - Det - Echelon - Reduced - Rank - Nullspace basis Echelon form - RankProfile - Inverse - Solve -1 U -1 LUDecomp U TRTRI U TRTRI TRTRM -1 A A -1 L L L 2/3 1/3 1/3 2/3 LU decomp Echelon form decomp Reduced Echelon A=LU Gauss (TA=E) form decomp TA=R 1 GaussJordan 2

Outline Reduced Echelon forms and Gaussian elimination Gaussian elimination based matrix decompositions Relations between decompositions Algorithms Hermite normal form Micciancio & Warinschi algorithm Double Determinant AddCol Frobenius normal form Krylov method Algorithm Reduction to matrix multiplication

Computing Hermite Normal form Equivalence over a ring: H = UA , where det ( U ) = ± 1  p 1 ∗ x 1 , 2 ∗ ∗ x 1 , 3 ∗  ∗ ∗ x 2 , 3 ∗ p 2 Hermite normal form: H =    , with 0 ≤ x ∗ , j < p j   ∗ p 3  Reduced Echelon form over a Ring Improving Micciancio-Warinshi algorihm n 5 log � A � n 3 log � A � ◮ O ˜ � � � � (heuristically: O ˜ ) n 2 log � A � ◮ space: O � � ⇒ Good on random matrices, common in num. theory, crypto.

Computing Hermite Normal form Equivalence over a ring: H = UA , where det ( U ) = ± 1  p 1 ∗ x 1 , 2 ∗ ∗ x 1 , 3 ∗  ∗ ∗ x 2 , 3 ∗ p 2 Hermite normal form: H =    , with 0 ≤ x ∗ , j < p j   ∗ p 3  Reduced Echelon form over a Ring Improving Micciancio-Warinshi algorihm n 5 log � A � n 3 log � A � ◮ O ˜ � � � � (heuristically: O ˜ ) n 2 log � A � ◮ space: O � � ⇒ Good on random matrices, common in num. theory, crypto. Implementation, reduction to building blocks: ◮ LinSys over Z , ◮ CUP and MatMul over Z p

Naive algorithm begin 1 foreach i do 2 ( g , t i , . . . , t n ) = xgcd ( A i , i , A i + 1 , i , . . . , A n , i ) ; 3 L i ← � n j = i + 1 t j L j ; 4 5 for j = i + 1 . . . n do A j , i L j ← L j − g L i ; /* eliminate */ 6 for j = 1 . . . i − 1 do 7 A j , i L j ← L j − ⌊ g ⌋ L i ; /* reduce */ 8 end 9   p 1 ∗ x 1 , 2 ∗ ∗ x 1 , 3 ∗ p 2 ∗ ∗ x 2 , 3 ∗     p 3 ∗  

Computing modulo the determinant [Domich & Al. 87] Property For A non-singular: max i � j H ij ≤ det H Example     − 5 8 − 3 − 9 5 5 1 0 3 237 − 299 90 − 2 8 − 2 − 2 8 5 0 1 1 103 − 130 40     A =  , H =     7 − 5 − 8 4 3 − 4 0 0 4 352 − 450 135    1 − 1 6 0 8 − 3 0 0 0 486 − 627 188 det A = 1944

Computing modulo the determinant [Domich & Al. 87] Property For A non-singular: max i � j H ij ≤ det H Example     − 5 8 − 3 − 9 5 5 1 0 3 237 − 299 90 − 2 8 − 2 − 2 8 5 0 1 1 103 − 130 40     A =  , H =     7 − 5 − 8 4 3 − 4 0 0 4 352 − 450 135    1 − 1 6 0 8 − 3 0 0 0 486 − 627 188 det A = 1944 Moreover, every computation can be done modulo d = det A : � A � � H � U ′ = dI n I n I n n 4 log � A � � n 3 � � � ⇒O × M ( n ( log n + log � A � )) = O ˜

Computing modulo the determinant ◮ Pessimistic estimate on the arithmetic size ◮ d large but most coefficients of H are small ◮ On the average : only the last few columns are large ⇒ Compute H ′ close to H but with small determinant

Computing modulo the determinant ◮ Pessimistic estimate on the arithmetic size ◮ d large but most coefficients of H are small ◮ On the average : only the last few columns are large ⇒ Compute H ′ close to H but with small determinant [Micciancio & Warinschi 01]   B b c T A = a n − 1 , n   d T a n , n �� B �� B �� �� d 1 = det , d 2 = det c T d T g = gcd ( d 1 , d 2 ) = sd 1 + td 2 Then �� �� B det = g sc T + td T

Micciancio & Warinschi algorithm begin 1 �� B �� B �� �� Compute d 1 = det , d 2 = det ; /* Double Det */ c T d T 2 ( g , s , t ) = xgcd ( d 1 , d 2 ) ; 3 � � B Compute H 1 the HNF of mod g ; /* Modular HNF */ sc T + td T 4 � � B b Recover H 2 the HNF of ; /* AddCol */ sc T + td T sa n − 1 , n + ta n , n 5   B b  ; c T Recover H 3 the HNF of /* AddRow */ a n − 1 , n  d T a n , n 6 end 7

Micciancio & Warinschi algorithm begin 1 �� B �� B �� �� Compute d 1 = det , d 2 = det ; /* Double Det */ c T d T 2 ( g , s , t ) = xgcd ( d 1 , d 2 ) ; 3 � � B Compute H 1 the HNF of mod g ; /* Modular HNF */ sc T + td T 4 � � B b Recover H 2 the HNF of ; /* AddCol */ sc T + td T sa n − 1 , n + ta n , n 5   B b  ; c T Recover H 3 the HNF of /* AddRow */ a n − 1 , n  d T a n , n 6 end 7

Double Determinant First approach: LU mod p 1 , . . . , p k + CRT ◮ Only one elimination for the n − 2 first rows ◮ 2 updates for the last rows (triangular back substitution) i = 1 p i > n n log � A � n / 2 ◮ k large such that � k

Double Determinant First approach: LU mod p 1 , . . . , p k + CRT ◮ Only one elimination for the n − 2 first rows ◮ 2 updates for the last rows (triangular back substitution) i = 1 p i > n n log � A � n / 2 ◮ k large such that � k Second approach: [Abbott Bronstein Mulders 99] ◮ Solve Ax = b . ◮ δ = lcm ( q 1 , . . . , q n ) s.t. x i = p i / q i Then δ is a large divisor of D = det A . ◮ Compute D /δ by LU mod p 1 , . . . , p k + CRT i = 1 p i > n n log � A � n / 2 /δ ◮ k small , such that � k

Double Determinant : improved Property � � Let x = [ x 1 , . . . , x n ] be the solution of x = d . Then A c x n , . . . , − x n − 1 y = [ − x 1 x n , 1 � � x n ] is the solution of y = c . A d ◮ 1 system solve ◮ 1 LU for each p i ⇒ d 1 , d 2 computed at about the cost of 1 déterminant

AddCol Problem Find a vector e such that � � B b � � = U H 1 e sc T + td T sa n − 1 , n + ta n , n � � b = e U sa n − 1 , n + ta n , n � − 1 � � � B b = H 1 sc T + td T sa n − 1 , n + ta n , n ⇒ Solve a system. ◮ n − 1 first rows are small ◮ last row is large

AddCol Idea: replace the last row by a random small one w T . � B � � � b y = w T a n − 1 , n − 1 Let k be a basis of the kernel of B . Then x = y + α k . where α = a n − 1 , n − 1 − ( sc T + td T ) · y ( sc T + td T ) · k ⇒ limits the expensive arithmetic to a few dot products

Outline Reduced Echelon forms and Gaussian elimination Gaussian elimination based matrix decompositions Relations between decompositions Algorithms Hermite normal form Micciancio & Warinschi algorithm Double Determinant AddCol Frobenius normal form Krylov method Algorithm Reduction to matrix multiplication

Krylov Method Definition (degree d Krylov matrix of one vector v ) � A d − 1 v � K = . . . v Av Property   ∗ 0 1 ∗   A × K = K ×   ...   ∗   ∗ 1

Krylov Method Definition (degree d Krylov matrix of one vector v ) � A d − 1 v � K = . . . v Av Property   ∗ 0 1 ∗   A × K = K ×   ...   ∗   ∗ 1 ⇒ if d = n , K − 1 AK = C P A car

Krylov Method Definition (degree d Krylov matrix of one vector v ) � A d − 1 v � K = . . . v Av Property   ∗ 0 1 ∗   A × K = K ×   ...   ∗   ∗ 1 ⇒ if d = n , K − 1 AK = C P A car ⇒ [Keller-Gehrig, alg. 2] computes K in O ( n ω log n )

Definition (degree k Krylov matrix of several vectors v i ) � A k − 1 v 1 A k − 1 v 2 A k − 1 v l � K = v 1 . . . v 2 . . . . . . v l . . . Property ≤ k k k 0 1 1 0 1 1 A × K = K × 0 1 1

Hessenberg poly-cyclic form Fact If ( d 1 , . . . d l ) is lexicographically maximal such that A d 1 − 1 v 1 A d l − 1 v l � � K = v 1 . . . . . . v l . . . is non-singular, then d 1 d 2 dl 0 1 1 0 1 A × K = K × 1 0 1 1

Principle k -shifted form: k k ≤ k 0 1 1 0 1 1 0 1 1