Eigenvalue Problems and Singular Value Decomposition Sanzheng Qiao - PowerPoint PPT Presentation

EIG Singular Value Decomposition Software Summary Eigenvalue Problems and Singular Value Decomposition Sanzheng Qiao Department of Computing and Software McMaster University August, 2012

EIG Singular Value Decomposition Software Summary Outline Eigenvalue Problems 1 Singular Value Decomposition 2 Software Packages 3

EIG Singular Value Decomposition Software Summary Outline Eigenvalue Problems 1 Singular Value Decomposition 2 Software Packages 3

EIG Singular Value Decomposition Software Summary Problem setting Eigenvalue problem: A x = λ x , λ : eigenvalue x : right eigenvector. y H A = λ y H , y left eigenvector.

EIG Singular Value Decomposition Software Summary Canonical forms Decomposition: A = SBS − 1 where B is in a canonical (simple) form, whose eigenvalues and eigenvectors can be easily obtained. A and B have the same eigenvalues. (They are similar.) If x is an eigenvector of B , then S x is the eigenvector of A corresponding to the same eigenvalue.

EIG Singular Value Decomposition Software Summary Jordan canonical form A = SJS − 1 , J = diag ( J n 1 ( λ 1 ) , ..., J n k ( λ k ))   λ i 1 0   ... ...   J n i ( λ i ) =   .   ...   1 0 λ i

EIG Singular Value Decomposition Software Summary Jordan canonical form The algebraic multiplicity of λ i is n i . A Jordan block has one right eigenvector [ 1 , 0 , ..., 0 ] T and one left eigenvector [ 0 , ..., 0 , 1 ] T . If all n i = 1, then J is diagonal, A is called diagonalizable; otherwise, A is called defective. An n -by- n defective matrix has fewer than n eigenvectors.

EIG Singular Value Decomposition Software Summary Example In practice, confronting defective matrices is a fundamental fact. Mass-spring problem k m b

EIG Singular Value Decomposition Software Summary Mass-spring problem Newton’s law F = ma implies m ¨ x ( t ) = − kx ( t ) − b ˙ x ( t ) . Let � ˙ � x ( t ) y ( t ) = , x ( t ) we transform the second order ODE into a system of the first order ODEs � � − b − k y ( t ) = m m y ( t ) =: A y ( t ) . ˙ 1 0

EIG Singular Value Decomposition Software Summary Mass-spring problem The characteristic polynomial of A is λ 2 + b m λ + k m and the eigenvalues are �� � 2 − b b − 4 k � � � m ± b 1 − 4 km m m λ ± = = − 1 ± . 2 m b 2 2 When 4 km / b 2 = 1, critically damped, two equal eigenvalues, A is not diagonalizable.

EIG Singular Value Decomposition Software Summary Jordan canonical form A = SJS − 1 � � − b − k A = m m 4 km = b 2 , , 1 0 � − b � � � − b b 1 1 − J = S = 2 m 2 m 2 m , − b 0 1 1 2 m

EIG Singular Value Decomposition Software Summary Jordan canonical form It is undesirable to compute Jordan form, because Jordan block is discontinuous � 0 � ǫ � � 1 1 J ( 0 ) = J ( ǫ ) = , 2 ǫ 0 0 0 while J ( 0 ) has an eigenvalue of multiplicity two, J ( ǫ ) has two simple eigenvalues. In general, computing Jordan form is unstable, that is, S − 1 = A + E for a small E . S � J � there is no guarantee that �

EIG Singular Value Decomposition Software Summary Schur canonical form A = QTQ H Q : unitary T : upper triangular The eigenvalues of A are the diagonal elements of T .

EIG Singular Value Decomposition Software Summary Schur canonical form A = QTQ H Q : unitary T : upper triangular The eigenvalues of A are the diagonal elements of T . Real case A = QTQ T Q : orthogonal T : quasi-upper triangular, 1-by-1 or 2-by-2 blocks on the diagonal.

EIG Singular Value Decomposition Software Summary Conditioning Let λ be a simple eigenvalue of A with unit right eigenvector x and left eigenvector y . λ + ǫ be the corresponding eigenvalue of A + E , then ǫ = y H E x y H x + O ( � E � 2 ) or 1 | y H x |� E � + O ( � E � 2 ) . | ǫ | ≤ Condition number for finding a simple eigenvalue 1 / | y H x |

EIG Singular Value Decomposition Software Summary Computing the Schur decomposition A = QTQ T , T : quasi-upper triangular Step 1: Reduce A to upper Hessenberg A = Q 1 HQ T h ij = 0 , i > j + 1 1 , Step 2: Compute the Schur decomposition of H H = Q 2 TQ T 2

EIG Singular Value Decomposition Software Summary Introducing zeros into a vector Householder transformation H = I − 2 uu T with u T u = 1 H is symmetric and orthogonal ( H 2 = I ). Goal: H a = α e 1 .

EIG Singular Value Decomposition Software Summary Introducing zeros into a vector Householder transformation H = I − 2 uu T with u T u = 1 H is symmetric and orthogonal ( H 2 = I ). Goal: H a = α e 1 . Choose u = a ± � a � 2 e 1

EIG Singular Value Decomposition Software Summary A geometric interpretation u a a u e 1 b (a) (b) Figure (a) shows the image b = ( I − 2 uu T ) a for an arbitrary u , in figure (b), u = a − � a � 2 e 1 .

EIG Singular Value Decomposition Software Summary Computing Householder transformations Given a vector x , it computes scalars σ , α , and vector u such that ( I − σ † uu T ) x = − α e 1 where σ † = 0 if σ = 0 and σ − 1 otherwise α = sign x 1 � x � 2 u = x + α e 1 2 = 2 ( α 2 + α x 1 ) = 2 α u 1 � u � 2

EIG Singular Value Decomposition Software Summary house.m function [u,sigma,alpha] = house(x) u = x; alpha = sign(u(1))*norm(u); u(1) = u(1) + alpha; sigma = alpha*u(1);

EIG Singular Value Decomposition Software Summary Reducing A to upper Hessenberg n = length(A(1,:)); Q = eye(n); for j=1:(n-2) [u, sigma, alpha] = myhouse(A(j+1:n, j)); if sigma ˜= 0.0 for k = j:n A(j+1:n, k) = A(j+1:n, k) - ((u’*A(j+1:n, k))/sigma)*u; end %for k for i=1:n A(i, j+1:n) = A(i, j+1:n) - ((A(i, j+1:n)*u)/sigma)*u’; Q(i, j+1:n) = Q(i, j+1:n) - ((Q(i, j+1:n)*u)/sigma)*u’; end %for i end %if end %for j

EIG Singular Value Decomposition Software Summary Computing eigenvalues and eigenvectors Suppose A has distinct eigenvalues λ i , i = 1 , ..., n , where | λ 1 | > | λ 2 | ≥ ... ≥ | λ n | , and x i are the eigenvectors (linear independent). An arbitrary vector u can be expressed as u = µ 1 x 1 + µ 2 x 2 + · · · + µ n x n If µ 1 � = 0, A k u has almost the same direction as x 1 when k and large and thus ( λ i /λ 1 ) k ( i > 1) is small. Thus the Rayleigh quotient ( A k u ) T A ( A k u ) ( A k u ) T ( A k u ) ≈ λ 1 .

EIG Singular Value Decomposition Software Summary Power method Initial u 0 ; i = 0; repeat v i + 1 = A u i ; u i + 1 = v i + 1 / � v i + 1 � 2 ; λ i + 1 = u T i + 1 A u i + 1 ; ˜ i = i + 1 until convergence

EIG Singular Value Decomposition Software Summary Power method Initial u 0 ; i = 0; repeat v i + 1 = A u i ; u i + 1 = v i + 1 / � v i + 1 � 2 ; λ i + 1 = u T i + 1 A u i + 1 ; ˜ i = i + 1 until convergence Problems Computes only ( λ 1 , x 1 ) Converges slowly when | λ 1 | ≈ | λ 2 | Does not work when | λ 1 | = | λ 2 |

EIG Singular Value Decomposition Software Summary Inverse power method Suppose that µ is an estimate of λ k , then ( λ k − µ ) − 1 is the dominant eigenvalue of ( A − µ I ) − 1 . Applying the power method to ( A − µ I ) − 1 , we can compute x k and λ k . Example. Eigenvalues of A : − 1, − 0 . 2, 0 . 5, 1 . 5 Shift µ : − 0 . 8 Eigenvalues of ( A − µ I ) − 1 : − 5 . 0, 1 . 7, 0 . 78, 0 . 43 Very effective when we have a good estimate for an eigenvalue.

EIG Singular Value Decomposition Software Summary QR method Goal: Generate a sequence A 0 = A , A 1 , ..., A k + 1 � B � u A i + 1 = Q T i AQ i = s T µ where s is small and Q i is orthogonal, i.e., Q T i = Q − 1 (so A k + 1 i and A have the same eigenvalues). Since s is small, µ is an approximation of an eigenvalue of A k + 1 ( A ); Deflate A k + 1 and repeat the procedure on B when s is sufficiently small. The problem size is reduced by one.

EIG Singular Value Decomposition Software Summary QR method What does Q k look like? If the last column of Q k is a left eigenvector y of A , then � P T � Q T k AQ k k A [ P k y ] = y T � P T � k [ AP k A y ] = y T � B � u = 0 T λ

EIG Singular Value Decomposition Software Summary QR method How do we get an approximation of a left eigenvector y of A ( y T A = λ y T )?

EIG Singular Value Decomposition Software Summary QR method How do we get an approximation of a left eigenvector y of A ( y T A = λ y T )? One step of the inverse power method: Solve for q in ( A − µ I ) T q = e n , where µ is an estimate for an eigenvalue of A . (How? Later.)

EIG Singular Value Decomposition Software Summary QR method How do we get an approximation of a left eigenvector y of A ( y T A = λ y T )? One step of the inverse power method: Solve for q in ( A − µ I ) T q = e n , where µ is an estimate for an eigenvalue of A . (How? Later.) How do we construct an orthogonal Q whose last column is q ?