Eigenvalue/Eigenvector Sensitivity 

Hakan Tiftikci
Turkey
hakan.tiftikci@yahoo.com.tr 

Introduction 

In this worksheet, a method to compute sensitivities of eigenvalues/eigenvectors of a matrix with respect to parameters of the matrix is presented. Method is based on the Cayley-Hamilton theorem, which states that, the matrix itself satisfies its characteristic polynomial, i.e. 

 

 

 

where 

 

 

 

is the characteristic polynomial of matrix . Factored form of characteristic polynomial is given by 

 

 

 

which have a connection to algebraic form by symmetric functions of eigenvalues by 

 

 

 

Roots are the eigenvalues of matrix . 

 

Cayley-Hamilton theorem can be expressed using the factored form of characteristic polynomials, viz. 

 

 

 

Taking partial derivative of charactristic polynomial with respect to a parameter that the matrix depends on it, 

 

 

 

or if written more concisely 

 

 

 

Denoting 

 

 

 

the last equation becomes 

 

 

 

expanding and isolating eigenvalue derivatives on one side, following final matrix equation is obtained 

 

 

 

which is equations in unknown eigenvalue derivatives Selecting appropriate (well-formed/nonsingular) equations from a total of , eigenvalue derivatives can be obtained. For the derivatives of eigenvectors , the defining equation 

 

 

 

 

 

is differentiated with respect to to yield 

 

 

 

 

 

Initialization 

Implementation and Application 

In the first subsection, method is first illustrated step-by-step for a simple system and then in second part is put to procedural form.  

 

Application to a simple system (Spring-Mass-Damper) 

When spring-mass-damper system is converted to state-space form by matrix in linear representation is given by following 

>	A := <<0, -omega^2>|<1, -2*omega*zeta>>;

 

(3.1.1)

 

Characteristic polynomial is given by 

>	CharacteristicPolynomial(A, 's');

 

(3.1.2)

 

Eigenvalues are 

>	Lambda := Eigenvalues(A);

 

(3.1.3)

 

Matrix satisfies Cayley-Hamilton theorem,  

 

 

 

as shown by following computation 

>	J2 := Matrix(2,2,shape=identity);

 

(3.1.4)

 

>	testCayleyHamilton := A.A+MatrixScalarMultiply(A, 2*omega*zeta)+omega^2*J2;

 

(3.1.5)

 

Sensitivities of eigenvalues (by direct computation) are 

>	Mu1 := map(diff, Lambda, omega); Mu2 := map(diff, Lambda, zeta);

 

 

	(3.1.6)

 

Derivatives of the matrix with respect to its parameters 

>	A1 := map(diff,A,omega); A2 := map(diff,A,zeta);

 

 

	(3.1.7)

 

Derivative with respect to first parameter  

Denote derivatives of by , then left and right hand side of  E.1 and the matrix equation are 

>	LHS := mu1*(A-Lambda[2]*J2)+mu2*(A-Lambda[1]*J2): RHS := A1.(A-Lambda[2]*J2)+(A-Lambda[1]*J2).A1: EQU := LHS-RHS;

 

(3.1.1.1)

 

solving by elements of matrix equation 

>	sol_11_22 := solve({EQU[1, 1], EQU[2, 2]}, {mu1, mu2});

 

(3.1.1.1.1)

 

compare last result with directly computed results 

>	eval(<mu1,mu2>,sol_11_22)-Mu1;

 

(3.1.1.1.2)

 

solving by  

>	sol_11_22 := solve({EQU[1, 1], EQU[2, 2]}, {mu1, mu2});

 

(3.1.1.2.1)

 

compare last result with directly computed results 

>	eval(<mu1,mu2>,sol_11_22)-Mu1;

 

(3.1.1.2.2)

 

Performing computations for the and elements,  equivalence of solution obtained from different element selections of matrix equation (E.1) are demostrated (just) for the case. 

 

Derivative with respect to second parameter  

Denote derivatives of by , then left and right hand side of  E.1 and the matrix equation are 

>	LHS := mu1*(A-Lambda[2]*J2)+mu2*(A-Lambda[1]*J2): RHS := A2.(A-Lambda[2]*J2)+(A-Lambda[1]*J2).A2: EQU := LHS-RHS;

 

(3.1.2.1)

 

solving by and elements of matrix equation 

>	sol_11_12 := solve({EQU[1, 1], EQU[1, 2]}, {mu1, mu2});

 

(3.1.2.2)

 

compare last result with directly computed results 

>	eval(<mu1,mu2>,sol_11_12)-Mu2; simplify(%);

 

 

	(3.1.2.3)

 

Sensitivity Procedure  

In this section, the method whose step are explained in above example is defined as a Maple procedure. 

>

EigenSensitivity := proc(A,q,L) local LHS,RHS,N,i,Aq,AP,MATEQ,sol; description "computes the partial derivatives of eigenvalues 'L' of matrix 'A' with respect to parameter 'q'"; # if not given, compute the eigenvalues #if evalb(L=NULL) then # L := Eigenvalues(A); #end if; N := RowDimension(A);  # system size # derivative of matrix with respect to parameter "q" Aq := map(diff,A,q); # define auxiliary product function AP := proc(a,b) local i,result,J; description "auxiliary product function for defining matrix equation"; J := Matrix(N,N,shape=identity); result := Matrix(N,N,shape=identity); for i from a to b do result := result . (A-L[i]*J); end do; result; end proc; # construct the left and right side of matrix equation (E.1) LHS := Matrix(N,N,0); RHS := Matrix(N,N,0); for i from 1 to N do LHS := LHS + AP(1,i-1) . Aq . AP(i+1,N); RHS := RHS + MatrixScalarMultiply( AP(1,i-1) . AP(i+1,N), (mu||i)); end do; # form matrix equation MATEQ := LHS-RHS; # solve matrix equation for eigenvalue derivatives sol := solve( {seq(MATEQ[i,i],i=1..N)}, {mu||(1..N)}); # as default, linear system is selected from diagonals of matrix equation # return (simplified) result as vector #map(simplify,eval(<mu||(1..N)>,sol)); eval(<mu||(1..N)>,sol); end proc:

 

>	#untrace(EigenSensitivity);

 

 

Sensitivity Procedure Test 

Test the procedure on Spring-Mass-Damper 

>	EigenSensitivity(A,omega,Lambda); EigenSensitivity(A,zeta,Lambda); simplify(%%-Mu1); simplify(%%-Mu2);

 

 

 

 

	(3.3.1.1)

 

Test the procedure on general 2x2 matrix 

>	A2 := <<a11,a21>|<a12,a22>>;

 

(3.3.2.1)

 

>	Lambda2 := Eigenvalues(A2);

 

(3.3.2.2)

 

derivatives with respect to  

compute sensitivities directly 

>	Mu_11_direct := map(diff, Lambda2, a11);

 

(3.3.2.1.1)

 

compute sensitivities by the method 

>	Mu_11_method := EigenSensitivity(A2,a11,Lambda2);

 

(3.3.2.1.2)

 

>	simplify(Mu_11_direct-Mu_11_method);

 

(3.3.2.1.3)

 

derivatives with respect to  

compute sensitivities directly 

>	Mu_12_direct := map(diff, Lambda2, a12);

 

(3.3.2.2.1)

 

compute sensitivities by the method 

>	Mu_12_method := EigenSensitivity(A2,a12,Lambda2);

 

(3.3.2.2.2)

 

>	simplify(Mu_12_direct-Mu_12_method);

 

(3.3.2.2.3)

 

Test the procedure on simplified 3x3 matrix 

For this case eigenvalues used in obtaining their derivatives are kept as symbols until last step in which results are compared to those of direct computation. Otherwise computational resources are quickly consumed. For the same reason, in this case, matrix is simplified by fixing some elements and reducing number of parameters. 

>	A3 := Matrix(3,3, (i,j) -> a||i||j, shape=symmetric): A3[2,1]:=1: A3[3,1]:=2: A3[3,2]:=2*a22: A3;

 

(3.3.3.1)

 

>	Lambda3 := Eigenvalues(A3):

 

derivatives with respect to  

compute sensitivities directly 

>	Mu_22_direct := map(diff, Lambda3, a22):

 

compute sensitivities by the method 

>	Mu_22_method := EigenSensitivity(A3,a22,<L3_1,L3_2,L3_3>):

 

>	mustBeZero := subs([L3_1=Lambda3[1],L3_2=Lambda3[2],L3_3=Lambda3[3]],Mu_22_method)-Mu_22_direct:

 

>	simplify(mustBeZero);

 

(3.3.3.1.1)

 

Eigenvalue Jacobian and a Related Application 

Jacobian of eigenvalue vector is obtained by repeated application of the method ("EigenSensitivity"). A more efficient version may be obtained by identifying the repeated multiplications and/or common factors.

Eigenvalue Jacobian can be used to determine required direction/rate in parameter space for attaining a direction/rate in eigenvalue space, as formulated by following expressions 

 

 

 

or 

 

 

 

 

if it is desired to direct/move current eigenvalue vector towards a target eigenvalue vector then the choice  

 

 

 

serves to this purpose. 

>

EigenJacobian := proc(A,Q,L) local result,i,Lq; description "computes the Jacobian of matrix 'A' with respect to parameters 'Q'"; result := Matrix( RowDimension(A), Dimension(Q)); for i from 1 to Dimension(Q) do Lq := EigenSensitivity(A,Q[i],L); result[1..RowDimension(A),i] := Lq; end do; result; end proc:

 

 

Application to Spring-Mass-Damper 

Above formulation is applied to move the parameters from an arbitrary starting point to a desired target point. First Jacobian is computed as 

>	J := EigenJacobian(A,<omega,zeta>,Lambda);

 

(3.4.1.1)

 

>	Jinv := MatrixInverse(J);

 

(3.4.1)

 

define gain (positive definite) and target eigenvalues (yet as symbols) 

>	K := <<k11,k21>|<k12,k22>>; Ltarget := <L1,L2>;

 

 

	(3.4.2)

 

express differential vector in parameter space that corresponds to direction in eigenvalue space. 

>	dp := Jinv . K . (Ltarget-Lambda);

 

(3.4.3)

 

study-1 

Eigenvalues are moved to  

>	SL := [L1=-3,L2=-4,k11=1,k12=0,k21=0,k22=1,zeta=zeta(t),omega=omega(t)];

 

(3.4.2.1)

 

>	DEQ := diff(omega(t),t)=subs(SL,dp[1]), diff(zeta(t),t)=subs(SL,dp[2]);

 

(3.4.2.2)

 

solve differential system numerically (note complex option) 

>	dsol := dsolve({DEQ,zeta(0)=1/10,omega(0)=1},{zeta(t),omega(t)},'numeric',complex=true);

 

(3.4.2.3)

 

Plot locus in parameter space 

>	zetaplot := plots[odeplot](dsol,[Re(zeta(t)),Im(zeta(t))],t=0..5,style=POINT,color=red): omegaplot := plots[odeplot](dsol,[Re(omega(t)),Im(omega(t))],t=0..5,style=POINT,color=blue):

 

>	plots[display]({zetaplot,omegaplot},axes=boxed,labels=["Re","Im"],title="parameters in complex plane");

 

 

plot history of eigenvalues (compare to target values) 

>	solutionEigen1 := subs(omega=omega(t),zeta=zeta(t),Lambda[1]); solutionEigen2 := subs(omega=omega(t),zeta=zeta(t),Lambda[2]);

 

 

	(3.4.2.4)

 

>

plots[odeplot](dsol,[t,Re(solutionEigen1)],t=0..4,numpoints=64,style=POINT,color=blue,axes=frame,labels=["t","Re(L1)"]); plots[odeplot](dsol,[t,Im(solutionEigen1)],t=0..4,numpoints=64,style=POINT,color=blue,axes=frame,labels=["t","Im(L1)"]); plots[odeplot](dsol,[Re(solutionEigen1),Im(solutionEigen1)],t=0..4,numpoints=64,style=POINT,color=blue,axes=frame,labels=["Re(L1)","Im(L1)"]);

 

 

 

 

>

plots[odeplot](dsol,[t,Re(solutionEigen2)],t=0..4,numpoints=64,style=POINT,color=blue,axes=frame,labels=["t","Re(L2)"]); plots[odeplot](dsol,[t,Im(solutionEigen2)],t=0..4,numpoints=64,style=POINT,color=blue,axes=frame,labels=["t","Im(L2)"]); plots[odeplot](dsol,[Re(solutionEigen2),Im(solutionEigen2)],t=0..4,numpoints=64,style=POINT,color=blue,axes=frame,labels=["Re(L2)","Im(L2)"]);

 

 

 

 

study-2 

eigenvalues are moved to  

>	SL := [L1=-3+2*I,L2=-4,k11=1,k12=0,k21=0,k22=1,zeta=zeta(t),omega=omega(t)];

 

(3.4.3.1)

 

>	DEQ := diff(omega(t),t)=subs(SL,dp[1]), diff(zeta(t),t)=subs(SL,dp[2]);

 

(3.4.3.2)

 

>	dsol := dsolve({DEQ,zeta(0)=1/10,omega(0)=1},{zeta(t),omega(t)},'numeric',complex=true);

 

(3.4.3.3)

 

>	zetaplot := plots[odeplot](dsol,[Re(zeta(t)),Im(zeta(t))],t=0..1.2,numpoints=64,style=POINT,color=red): omegaplot := plots[odeplot](dsol,[Re(omega(t)),Im(omega(t))],t=0..1.2,numpoints=64,style=POINT,color=blue):

 

>	plots[display]({zetaplot,omegaplot},axes=boxed,labels=["Re","Im"],title="parameters in complex plane");

 

 

>	solutionEigen1 := subs(omega=omega(t),zeta=zeta(t),Lambda[1]); solutionEigen2 := subs(omega=omega(t),zeta=zeta(t),Lambda[2]);

 

 

	(3.4.3.4)

 

>

plots[odeplot](dsol,[t,Re(solutionEigen1)],t=0..8,numpoints=64,style=POINT,color=blue,axes=frame,labels=["t","Re(L1)"]); plots[odeplot](dsol,[t,Im(solutionEigen1)],t=0..8,numpoints=64,style=POINT,color=blue,axes=frame,labels=["t","Im(L1)"]); plots[odeplot](dsol,[Re(solutionEigen1),Im(solutionEigen1)],t=0..8,numpoints=64,style=POINT,color=blue,axes=frame,labels=["Re(L1)","Im(L1)"]);

 

 

 

 

>

plots[odeplot](dsol,[t,Re(solutionEigen2)],t=0..8,numpoints=64,style=POINT,color=blue,axes=frame,labels=["t","Re(L2)"]); plots[odeplot](dsol,[t,Im(solutionEigen2)],t=0..8,numpoints=64,style=POINT,color=blue,axes=frame,labels=["t","Im(L2)"]); plots[odeplot](dsol,[Re(solutionEigen2),Im(solutionEigen2)],t=0..8,numpoints=64,style=POINT,color=blue,axes=frame,labels=["Re(L2)","Im(L2)"]);

 

 

 

 

Note that eigenvalue locus is straight line as imposed while developing differential relations. 

References 

[1] Cayley-Hamilton@Wikipedia 

[2] Cayley-Hamilton@PlanetMath 

Conclusions 

A general method to compute sensitivities of eigenvalues/eigenvectors is described. Method is applied to simple system and simplified system. Even though, the method described does not necessarily rely on symbolic computation, CAS environment provides robust way of implementing and verifying results. 

 

The described method, since involves many matrix multiplications,  is computationally intensive, and for the symbolic case, eigenvalues (if available) are generally so complex that, resulting matrix equation bears all the complexity of eigenvalues involved, causing deficiency in computer resources even though equations are linear in eigenvalue derivatives to be determined. For such cases eigenvalues may be left as symbols until their values are needed (as illustrated in test case) 

 

Additionaly,  following properties of method are observed: 

• All the eigenvalues must be available. Hence, pure symbolic solutions are limited by degree 5 (Galois) for the general case.

 

• Formulation has the structure of explicit ODE . Thus, for continuous change of parameters, corresponding eigenvalues can be traced by integration of .

 

• Derivatives of the matrices with respect to parameters must be available. In pure numeric implementation, this requirement poses another problem. Either numerical (finite difference,..), Automatic-Differentiation (AD), or CAS generated code for differentiation may be used.

 

 

Legal Notice: The copyright for this application is owned by the author(s). Neither Maplesoft nor the author are responsible for any errors contained within and are not liable for any damages resulting from the use of this material. This application is intended for non-commercial, non-profit use only. Contact the author for permission if you wish to use this application in for-profit activities.