I have a matrix and I would like to know if it is diagonalizable. How do I do this in the R programming language?
+2
A:
You can implement the full algorithm to check if the matrix reduces to a Jordan form or a diagonal one (see e.g., this document). Or you can take the quick and dirty way: for an n-dimensional square matrix, use eigen(M)$values and check that they are n distinct values. For random matrices, this always suffices: degeneracy has prob.0.
P.S.: based on a simple observation by JD Long below, I recalled that a necessary and sufficient condition for diagonalizability is that the eigenvectors span the original space. To check this, just see that eigenvector matrix has full rank (no zero eigenvalue). So here is the code:
diagflag = function(m,tol=1e-10){
x = eigen(m)$vectors
y = min(abs(eigen(x)$values))
return(y>tol)
}
# nondiagonalizable matrix
m1 = matrix(c(1,1,0,1),nrow=2)
# diagonalizable matrix
m2 = matrix(c(-1,1,0,1),nrow=2)
> m1
[,1] [,2]
[1,] 1 0
[2,] 1 1
> diagflag(m1)
[1] FALSE
> m2
[,1] [,2]
[1,] -1 0
[2,] 1 1
> diagflag(m2)
[1] TRUE
gappy
2009-07-25 02:47:38
Don't you need to worry about exactly what distinct means when calculated with floating point math?
hadley
2009-07-26 14:09:00
Yes, and that is where knowledge of the problem structure helps. If the matrix has certain random structures, e.g. gaussian iid, the one can claim certain distributional properties of eigenvalues (typically, uniform on the unit circle), and then implement simple tests for the null hypothesis that eigenvalues are identical. This paper is basically all I know about random matrices: http://tinyurl.com/eigenvalues. If nothing is known about the matrix, I'd try a few more heuristics..P.S.: I love ggplot2!
gappy
2009-07-27 03:36:58
Thanks for taking what little I knew and stretching it into a much better answer! Good job. This is exactly why I like StackOverflow; the good stuff floats to the top much more obviously than with a pure threaded discussion. (ggplot does, indeed, rock!)
JD Long
2009-07-27 13:33:26
+1
A:
If you have a given matrix, m, then one way is the take the eigen vectors times the diagonal of the eigen values times the inverse of the original matrix. That should give us back the original matrix. In R that looks like:
m <- matrix( c(1:16), nrow = 4)
p <- eigen(m)$vectors
d <- diag(eigen(m)$values)
p %*% d %*% solve(p)
m
so in that example [p %% d %% solve(p)] should be the same as m
and yes, I totally stole that from the R Wiki.
JD Long
2009-07-26 01:14:56
That's a good point. To answer the original question (to check for diagonalizability of matrix m): 1) compute the eigenvector matrix p = eigen(m)$vectors2) check that the eigenvector matrix p is invertible (i.e., solve(p) does not result in error). If yes, then it's diagonalizable, and the unitary action is P. If not, then it is not diagonalizable. Why? Because invertibility means that the dimension of the direct product of the eigenspaces is the same as the dimension of the original matrix, and this condition is necessary and sufficient for diagonalizability.
gappy
2009-07-27 04:02:58