[PPT] - Faster algorithms for the characteristic polynomial Clment P ERNET PowerPoint Presentation

SLIDE 1

State of the art A new algorithm The new algorithm into practice

Faster algorithms for the characteristic polynomial

Clément PERNET and Arne STORJOHANN

Symbolic Computation Group University of Waterloo, Canada.

ISSAC 2007, Waterloo, July 30

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 2

State of the art A new algorithm The new algorithm into practice

Problem Compute the characteristic polynomial of a dense matrix

ver a field
C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 3

State of the art A new algorithm The new algorithm into practice

Problem Compute the characteristic polynomial of a dense matrix

ver a field

Result Randomized Las-Vegas algorithm in O (nω) field operations for large fields (#F > 2n2).

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 4

State of the art A new algorithm The new algorithm into practice

Problem Compute the characteristic polynomial of a dense matrix

ver a field

Result Randomized Las-Vegas algorithm in O (nω) field operations for large fields (#F > 2n2). Improves previous complexity by a log n factor, Optimal reduction to Matrix multiplication.

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 5

State of the art A new algorithm The new algorithm into practice

Problem Compute the characteristic polynomial of a dense matrix

ver a field

Result Randomized Las-Vegas algorithm in O (nω) field operations for large fields (#F > 2n2). Improves previous complexity by a log n factor, Optimal reduction to Matrix multiplication. Practical efficiency. E.g. over Z547 909: n 500 5000 15 000 LinBox 0.91s 4m44s 2h20m magma-2.13 1.27s 15m32s 7h28m

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 6

State of the art A new algorithm The new algorithm into practice

Outline

1

State of the art

2

A new algorithm Shifted forms Principle of the new algorithm Complexity

3

The new algorithm into practice

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 7

State of the art A new algorithm The new algorithm into practice

Outline

1

State of the art

2

A new algorithm Shifted forms Principle of the new algorithm Complexity

3

The new algorithm into practice

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 8

State of the art A new algorithm The new algorithm into practice

Pre-Strassen age

Leverrier 1840: trace of powers of A, and Newton’s formula improved/rediscovered by Souriau, Faddeev, Frame and Csanky O

n4

, based on Matrix multiplication Suited for parallel computation model

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 9

State of the art A new algorithm The new algorithm into practice

Pre-Strassen age

Leverrier 1840: trace of powers of A, and Newton’s formula improved/rediscovered by Souriau, Faddeev, Frame and Csanky O

n4

, based on Matrix multiplication Suited for parallel computation model Danilevskii 1937: elementary row/column operations ⇒O

n3
C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 10

State of the art A new algorithm The new algorithm into practice

Pre-Strassen age

Leverrier 1840: trace of powers of A, and Newton’s formula improved/rediscovered by Souriau, Faddeev, Frame and Csanky O

n4

, based on Matrix multiplication Suited for parallel computation model Danilevskii 1937: elementary row/column operations ⇒O

n3

Hessenberg 1942: transformation to quasi-upper triangular and determinant expansion formula. ⇒O

n3
C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 11

State of the art A new algorithm The new algorithm into practice

Post-Strassen age

Preparata & Sarwate 1978: Update Csanky with fast matrix multiplication ⇒O

nω+1
C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 12

State of the art A new algorithm The new algorithm into practice

Post-Strassen age

Preparata & Sarwate 1978: Update Csanky with fast matrix multiplication ⇒O

nω+1

Keller-Gehrig 1985, alg.1: computes (A2i)i=1... log2 n to form a Krylov basis. O (nω log n) the best complexity up to now

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 13

State of the art A new algorithm The new algorithm into practice

Post-Strassen age

Preparata & Sarwate 1978: Update Csanky with fast matrix multiplication ⇒O

nω+1

Keller-Gehrig 1985, alg.1: computes (A2i)i=1... log2 n to form a Krylov basis. O (nω log n) the best complexity up to now Keller-Gehrig 1985, alg.2: inspired by Danilevskii, block

perations

O (nω) but only valid with generic matrices

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 14

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Outline

1

State of the art

2

A new algorithm Shifted forms Principle of the new algorithm Complexity

3

The new algorithm into practice

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 15

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Definition (degree d Krylov matrix of one vector v) K =

v

Av . . . Ad−1v

Property

A × K = K ×

2 6 6 6 4 ∗ 1 ∗ ... ∗ 1 ∗ 3 7 7 7 5 | {z } CPA,v

min

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 16

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Definition (degree d Krylov matrix of one vector v) K =

v

Av . . . Ad−1v

Property

A × K = K ×

2 6 6 6 4 ∗ 1 ∗ ... ∗ 1 ∗ 3 7 7 7 5 | {z } CPA,v

min

⇒if d = n, K −1AK = CPA

car

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 17

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Definition (degree d Krylov matrix of one vector v) K =

v

Av . . . Ad−1v

Property

A × K = K ×

2 6 6 6 4 ∗ 1 ∗ ... ∗ 1 ∗ 3 7 7 7 5 | {z } CPA,v

min

⇒if d = n, K −1AK = CPA

car

[Keller-Gehrig, alg. 2] : K −1AK in O (nω) for A generic

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 18

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Definition (degree k Krylov matrix of several vectors vi)

K =

v1

. . . Ak−1v1 v2 . . . Ak−1v2 . . . vl . . . Ak−1vl

Property

1 1 1 1 1 1 k ≤ k

A × K = K×

k

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 19

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Fact (Shift Hessenberg form) If (d1, . . . dl) is lexicographically maximal such that K =

v1

. . . Ad1−1v1 . . . vl . . . Adl−1vl

is non-singular, then

1 1 1 1 1 1

A × K = K×

d2 d1 dl

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 20

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Principle

k-shifted form:

1 1 1 1 1 1 k k ≤ k

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 21

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Principle

k + 1-shifted form:

1 1 1 1 1 1 k + 1 k + 1 dl

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 22

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Principle

Compute iteratively from 1-shifted form to d1-shifted form

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 23

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Principle

Compute iteratively from 1-shifted form to d1-shifted form each completed block appears in the increasing degree

rder
C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 24

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Principle

Compute iteratively from 1-shifted form to d1-shifted form each completed block appears in the increasing degree

rder

until the shifted Hessenberg form is obtained:

1 1 1 1 1 1

d2 d1 dl

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 25

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Example

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 26

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Example

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 27

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Example

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 28

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Example

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 29

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Example

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 30

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Example

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 31

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Example

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 32

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Example

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 33

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Example

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 34

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Example

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 35

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Example

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 36

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Example

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 37

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Example

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 38

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Lemma If #F > 2n2, the transformation will succeed with high

probability. Failure is detected.
C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 39

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Lemma If #F > 2n2, the transformation will succeed with high

probability. Failure is detected.

How to use fast matrix arithmetic ?

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 40

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Permutations: compressing the dense columns

1 1 1 1 1 1 1 1

×P Ak = = Q×

c2 c3 c3 c1 c1 c2

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 41

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Permutations: compressing the dense columns

1 1 1 1 1 1 1 1

×P Ak = = Q×

c2 c3 c3 c1 c1 c2 1 1 1 1 1 1 1 1

×P′ K = = Q′×

c2 c1 c2 c1

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 42

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Reduction to Matrix multiplication

Similarity transformation:

K −1AK = Q′T I ∗ ∗

P′TQ

I ∗ ∗

PQ′

I ∗ ∗

P′
C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 43

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Reduction to Matrix multiplication

Similarity transformation:

K −1AK = Q′T I ∗ ∗ P′TQ I ∗ ∗ PQ′ I ∗ ∗

P′
C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 44

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Reduction to Matrix multiplication

Similarity transformation:

K −1AK = Q′T I ∗ ∗ P′TQ I ∗ ∗ PQ′ I ∗ ∗

P′

⇒O

k

n

k

ω

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 45

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Reduction to Matrix multiplication

Similarity transformation:

K −1AK = Q′T I ∗ ∗ P′TQ I ∗ ∗ PQ′ I ∗ ∗

P′

⇒O

k

n

k

ω Rank profile: derived from LQUP ⇒O

k

n

k

ω

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 46

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

Reduction to Matrix multiplication

Similarity transformation: ⇒O

k

n

k

ω Rank profile: derived from LQUP ⇒O

k

n

k

ω

n

k=1

k n k ω = nω

n

k=1

1 kω−1 = O (nω)

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 47

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

A new type of reduction

xIn − A det(xIn − A)

dimension = 1 degree = n degree = 1 dimension = n

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 48

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

A new type of reduction

xIn − A det(xIn − A)

dimension = 1 degree = n degree = 1 dimension = n

Keller-Gehrig 2

1 1 1 1

dimension =

n 2i

degree = 2i

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 49

State of the art A new algorithm The new algorithm into practice Shifted forms Principle of the new algorithm Complexity

A new type of reduction

xIn − A det(xIn − A)

dimension = 1 degree = n degree = 1 dimension = n

Keller-Gehrig 2

1 1 1 1

dimension =

n 2i

degree = 2i

New algorithm

1 1 1 1 1 1 1 1 1

dimension = n

k

degree = k

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 50

State of the art A new algorithm The new algorithm into practice

Outline

1

State of the art

2

A new algorithm Shifted forms Principle of the new algorithm Complexity

3

The new algorithm into practice

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 51

State of the art A new algorithm The new algorithm into practice

Improving the preconditioning

The preconditioning phase: A ← U−1AU for a random matrix U. (reminds [Kaltofen, Krish- namoorthy, Saunders 87])

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 52

State of the art A new algorithm The new algorithm into practice

Improving the preconditioning

The preconditioning phase: A ← U−1AU for a random matrix U. (reminds [Kaltofen, Krish- namoorthy, Saunders 87]) Instead, use a block Krylov precon- ditioning: A ← V −1AV, V =

W

AW . . . Ac−1W

for a random n × n/c matrix W.
C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 53

State of the art A new algorithm The new algorithm into practice

Improving the preconditioning

The preconditioning phase: A ← U−1AU for a random matrix U. (reminds [Kaltofen, Krish- namoorthy, Saunders 87]) Instead, use a block Krylov precon- ditioning: A ← V −1AV, V =

W

AW . . . Ac−1W

for a random n × n/c matrix W.

Property V −1AV is in c shifted form.

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 54

State of the art A new algorithm The new algorithm into practice

Efficiency balancing parameter

c small: full square matrix multiplications, but more ops c large: tends to matrix-vector products, but less ops

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 55

State of the art A new algorithm The new algorithm into practice

Efficiency balancing parameter

c small: full square matrix multiplications, but more ops c large: tends to matrix-vector products, but less ops ⇒parameter c balances efficiency

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 56

State of the art A new algorithm The new algorithm into practice

Efficiency balancing parameter

c small: full square matrix multiplications, but more ops c large: tends to matrix-vector products, but less ops ⇒parameter c balances efficiency

120 130 140 150 160 170 180 190 200 210 220 230 50 100 150 200 Time (s) Preconditionning parameter c Finding the optimal preconditionning paramater, n=5000 1 block of order 5000 5 blocks of order 1000 10 blocks of order 500

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 57

State of the art A new algorithm The new algorithm into practice

Experiments

n LU-Krylov New algorithm 200 0.024 0.032 300 0.06s 0.088s 500 0.248s 0.316s 750 1.084s 1.288s 1000 2.42s 2.296s 5000 267.6s 153.9s 10 000 1827s 991s 20 000 14 652s 7097s 30 000 48 887s 24 928s Computation time for 1 Frobenius block matrices, Itanium2-64 1.3Ghz, 192Gb

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 58

State of the art A new algorithm The new algorithm into practice

Experiments

0.01 0.1 1 10 100 1000 10000 100000 100 1000 10000 100000 Time (s) Matrix order Comparison for 1 frobenius block matrices over Z/(547909) New algorithm LU−Krylov

Timing comparison between the new algorithm and LU-Krylov, logarithmic scales, Itanium2-64 1.3Ghz, 192Gb

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 59

State of the art A new algorithm The new algorithm into practice

Comparison to Magma and previous LinBox

5000 10000 15000 20000 25000 30000 2000 4000 6000 8000 10000 12000 14000 16000 Time (s) Matrix order Comparison for 1 frobenius block matrices Magma 2.13 LU−Krylov New algorithm

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 60

State of the art A new algorithm The new algorithm into practice

Conclusion and perspectives

Results: Las Vegas reduction to matrix multiplication, The Frobenius normal form is easily derivable in O (nω)... ...but no transformation matrix Adaptive combination with block Krylov in practice.

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial

SLIDE 61

State of the art A new algorithm The new algorithm into practice

Conclusion and perspectives

Results: Las Vegas reduction to matrix multiplication, The Frobenius normal form is easily derivable in O (nω)... ...but no transformation matrix Adaptive combination with block Krylov in practice. Still to be done: Condition on the size of the field is a limitation. Eberly’s algorithm ? Ideally: derandomization ? (deterministic) Unification with matrix polynomial algorithms

C. PERNET and A. STORJOHANN

Faster algorithms for the characteristic polynomial