eigenvectors when A-lx is singular with no solutio

问题:

How is R able to find eigenvectors for the following matrix? Eigenvalues are 2,2 so eigenvectors require solving solve(matrix(c(0,1,0,0),2,2)) which is singular matrix with no solution.

> eigen(matrix(c(2,1,0,2),2,2))
$values
[1] 2 2
$vectors
[,1]          [,2] 
[1,]    0  4.440892e-16
[2,]    1 -1.000000e+00

> solve(matrix(c(0,1,0,0),2,2))
Error in solve.default(matrix(c(0, 1, 0, 0), 2, 2)) : 
Lapack routine dgesv: system is exactly singular

Both the routines essentially do the same thing. They find x such that (A-lambdaI)x = 0 without finding the inverse of A-lambdaI. Clearly (0 1) is a solution but how I can't understand why solve did not come up with it and how do I manually solve it.

回答1:

Maybe it's using one of the algorithms listed here:

http://en.wikipedia.org/wiki/List_of_numerical_analysis_topics#Eigenvalue_algorithms

?

According to http://stat.ethz.ch/R-manual/R-devel/library/base/html/eigen.html, eigen seems to use the LAPACK routine at http://netlib.org/lapack/double/dgeev.f (if you have a square matrix which is not symmetric).

Note: you're right that A - lambda * I is singular if lambda is an eigenvalue but in order to find eigenvectors, one does need invert A - lambda * I or solve an equation y = (A - lambda * I) * x (with y not being the null vector). It is sufficient to find non-zero vectors x which satisfy

(A - lambda * I) * x = 0

One strategy is to find a non-singular transformation matrix T such that (A - lambda * I) * T is an upper triangular matrix (i.e. all elements below the diagonal are zero). Because A-lambda*I is singular, T can be constructed such that the last element on the diagonal (or even more diagonal elements if the multiplicity of the eigenvalue is larger than one) is zero.

A vector z which only has it's last element equal to a non-zero value (i.e. z = (0,....,0,1) ) will then give the zero vector when multiplied with (A-lambda *I) * T. So one has:

0 = ((A - lambda * I) * T) * z

or in other words, T*z is an eigenvector of A.

回答2:

You asked for an eigen decomposition, you got an eigen decomposition.

Had you asked for rcond(), the condition number of the matrix, or for kappa(), you would have gotten the appropriate response.

For your second example:

> mat <- matrix(c(0,1,0,0), 2, 2)
> kappa(mat)
[1] Inf
> 
> rcond(mat)
[1] 0
>

For your first example, there is actually no problem:

> mat <- matrix(c(2,1,0,2), 2, 2)
> kappa(mat)
[1] 1.772727
> 
> rcond(mat)
[1] 0.5714286
> 
> 

See e.g. this previous question on SO for more.

回答3:

They are not doing the same. Just take a look at the wikipedia link. And look in code for dgeev.f and you will see that a simple solve of (A-lambda*I)x=0 is not performed.

回答4:

Clearly (0 1) is a solution but how I can't understand why solve did not come up with it and how do I manually solve it.

You try to get eigenvectors by solving (A - λI)x = 0. Function solve will fail because it does LU decomposition and requires the system to be full-rank.

See section "How to find eigenvectors using textbook method" of my this answer and do the following.

## your matrix
A <- matrix(c(2,1,0,2), 2, 2)

## read my linked answer for function "NullSpace"
NullSpace(A - diag(2, nrow(A)))
#     [,1]
#[1,]    0
#[2,]    1