ycc ycc - 2 months ago 10
R Question

why solve(eigenvalues, eigenvectors) of cov are different for 2d and 3d matrix?

I'm learning eigenvalues, eigenvectors of cov by using R commands, I tried with 2d matrix first and the results seems right:

x <- c(60, 35, 74, 30, 80)
y <- c(58, 40, 68, 40, 70)
m <- cbind(x, y)
cv <- cov(m)
eig <- eigen(cv)
det(cv %*% eig$vectors) # 895.4
det(eig$values * diag(2) %*% eig$vectors) # 895.4
cv %*% eig$vectors
# [,1] [,2]
# x -604.7476 0.6751113
# y -391.0320 -1.0440884
solve(eig$vectors, eig$values*diag(2))
# [,1] [,2]
# [1,] -604.7476 -0.6751113
# [2,] 391.0320 -1.0440884


Then, I tried with 3d matrix:

x <- c(60, 35, 74, 30, 80)
y <- c(58, 40, 68, 40, 70)
z <- c(25, 75, 50, 60, 50)
m <- cbind(x, y, z)
cv <- cov(m)
eig <- eigen(cv)
det(cv %*% eig$vectors) # -178424.4
det(eig$values * diag(3) %*% eig$vectors) # -178424.4
cv %*% eig$vectors
# [,1] [,2] [,3]
# x -639.0223 78.10751 0.56623319
# y -416.7254 43.54886 -0.89605389
# z 389.0068 174.95930 -0.02974941
solve(eig$vectors, eig$values*diag(3))
# [,1] [,2] [,3]
# [1,] -639.0223 -95.61750 0.48169216
# [2,] 340.4124 43.54886 0.94419557
# [3,] 457.2808 -166.03865 -0.02974941


My question is: why "cv %% eig$vectors" and "solve(eig$vectors, eig$valuesdiag(n))" result the same numbers for 2d matrix, but result just the same diagonal numbers for 3d matrix?

Answer Source

cv %% eig$vectors and solve(eig$vectors, eig$valuesdiag(n)) are not equal even in the 2d case as the output in the question shows.

The first equals eig$vectors %*% diag(eig$values) and the second equals t(eig$vectors) %*% diag(eig$values) -- note the transpose.

To see this expand cv into eig$vectors %*% diag(eig$values) %*% t(eig$vectors) and note that eig$vectors is orthogonal (because cv is symmetric) so its transpose equals its inverse.

Note that det(eig$vectors) is 1 (it must be 1 or -1 because it is orthognoal) and if A and B are any two conformable square matrices and v is a vector then:

  • det(A %*% B) = det(A) * det(B)
  • det(t(A)) = det(A)
  • det(diag(v)) = prod(v)

so we have:

det(cv %*% eig$vectors)
= det(eig$vectors %*% diag(eig$values) %*% t(eig$vectors) %*% eig$vectors)
= det(eig$vectors) * det(diag(eig$values)) * det(t(eig$vectors)) * det(eig$vectors)
= 1 * prod(eig$values) * 1 * 1
= prod(eig$values)

det(eig$values * diag(2) %*% eig$vectors)
= det(diag(eig$values) %*% eig$vectors)
= det(diag(eig$values)) * det(eig$vectors)
= prod(eig$values) * 1  
= prod(eig$values)