Spectral Decomposition
## function (A)
## {
## InputOK = TRUE
## if (!is.matrix(A))
## InputOK = FALSE
## else {
## n = dim(A)[1]
## if (n != dim(A)[2])
## InputOK = FALSE
## for (i in 2:n) {
## for (j in 1:(i - 1)) if (A[i, j] - A[j, i] > 2 *
## .Machine$double.eps)
## InputOK = FALSE
## }
## }
## if (InputOK == FALSE) {
## warning("This is for real symmetric matrix only!")
## return(NULL)
## }
## V = diag(rep(1, n))
## MAXITER = 30
## for (Iteration in 1:MAXITER) {
## for (i in 1:(n - 1)) {
## for (j in (i + 1):n) {
## rotn = TRUE
## p = 0.5 * (A[i, j] + A[j, i])
## q = A[i, i] - A[j, j]
## t = sqrt(4 * p * p + q * q)
## if (t == 0)
## rotn = FALSE
## else {
## if (q >= 0) {
## oki = okj = FALSE
## if (abs(A[i, i]) == abs(A[i, i]) + 100 *
## abs(p))
## oki = TRUE
## if (abs(A[j, j]) == abs(A[j, j]) + 100 *
## abs(p))
## okj = TRUE
## if (oki == TRUE & okj == TRUE)
## rotn = FALSE
## else rotn = TRUE
## if (rotn == TRUE) {
## c = sqrt((t + q)/(2 * t))
## s = p/(t * c)
## }
## }
## else {
## rotn = TRUE
## s = sqrt((t - q)/(2 * t))
## if (p < 0)
## s = -s
## c = p/(t * s)
## }
## if (1 + abs(s) == 1)
## rotn = FALSE
## }
## if (rotn == TRUE) {
## for (k in 1:n) {
## q = A[i, k]
## A[i, k] = c * q + s * A[j, k]
## A[j, k] = -s * q + c * A[j, k]
## }
## for (k in 1:n) {
## q = A[k, i]
## A[k, i] = c * q + s * A[k, j]
## A[k, j] = -s * q + c * A[k, j]
## q = V[k, i]
## V[k, i] = c * q + s * V[k, j]
## V[k, j] = -s * q + c * V[k, j]
## }
## }
## }
## }
## }
## Res = list(diag(A), V)
## names(Res) = c("values", "vectors")
## return(Res)
## }
## <bytecode: 0x000001e2b4bc0960>
## <environment: namespace:math>
M = matrix(c(3, 2, 1, 2, 3, 2, 1, 2, 3), nrow=3) ; M
## [,1] [,2] [,3]
## [1,] 3 2 1
## [2,] 2 3 2
## [3,] 1 2 3
## eigen() decomposition
## $values
## [1] 6.3723 2.0000 0.6277
##
## $vectors
## [,1] [,2] [,3]
## [1,] 0.5418 -0.7071 0.4544
## [2,] 0.6426 0.0000 -0.7662
## [3,] 0.5418 0.7071 0.4544
## $values
## [1] 6.3723 2.0000 0.6277
##
## $vectors
## [,1] [,2] [,3]
## [1,] 0.5418 -7.071e-01 0.4544
## [2,] 0.6426 -5.969e-17 -0.7662
## [3,] 0.5418 7.071e-01 0.4544
## Loading required package: microbenchmark
microbenchmark(
eigen(M),
evJacobi(M)
)
## Unit: microseconds
## expr min lq mean median uq max neval
## eigen(M) 173.4 178.8 188.2 186.0 189.4 427.4 100
## evJacobi(M) 129.1 134.7 142.0 141.5 145.6 232.5 100
M2 = matrix(c(0.2, 0.1, 0.1, 0.1, 0.2, 0.1, 0.1, 0.1, 0.2), nrow=3) ; M2
## [,1] [,2] [,3]
## [1,] 0.2 0.1 0.1
## [2,] 0.1 0.2 0.1
## [3,] 0.1 0.1 0.2
## eigen() decomposition
## $values
## [1] 0.4 0.1 0.1
##
## $vectors
## [,1] [,2] [,3]
## [1,] 0.5774 0.0000 0.8165
## [2,] 0.5774 -0.7071 -0.4082
## [3,] 0.5774 0.7071 -0.4082
## $values
## [1] 0.4 0.1 0.1
##
## $vectors
## [,1] [,2] [,3]
## [1,] 0.5774 -0.2113 -0.7887
## [2,] 0.5774 0.7887 0.2113
## [3,] 0.5774 -0.5774 0.5774
M3 = matrix(c(0, 1, 1, 1, 0, 1, 1, 1, 0), nrow=3) ; M3
## [,1] [,2] [,3]
## [1,] 0 1 1
## [2,] 1 0 1
## [3,] 1 1 0
## eigen() decomposition
## $values
## [1] 2 -1 -1
##
## $vectors
## [,1] [,2] [,3]
## [1,] 0.5774 0.8165 0.0000
## [2,] 0.5774 -0.4082 -0.7071
## [3,] 0.5774 -0.4082 0.7071
## $values
## [1] 2 -1 -1
##
## $vectors
## [,1] [,2] [,3]
## [1,] 0.5774 -0.2113 -0.7887
## [2,] 0.5774 0.7887 0.2113
## [3,] 0.5774 -0.5774 0.5774
TestEV = function(n)
{
x = runif(n*(n+1)/2)
M = matrix(nrow=n, ncol=n)
for (i in 1:n) {
for (j in 1:i) {
M[i,j] = x[i] * (x[i] - 1) / 2 + j
M[j,i] = M[i,j]
}
}
print(M)
print(eigen(M))
print(evJacobi(M))
}
TestEV(2)
## [,1] [,2]
## [1,] 0.9475 0.9377
## [2,] 0.9377 1.9377
## eigen() decomposition
## $values
## [1] 2.5030 0.3823
##
## $vectors
## [,1] [,2]
## [1,] 0.5163 -0.8564
## [2,] 0.8564 0.5163
##
## $values
## [1] 2.5030 0.3823
##
## $vectors
## [,1] [,2]
## [1,] 0.5163 -0.8564
## [2,] 0.8564 0.5163
## [,1] [,2] [,3]
## [1,] 0.9585 0.8951 0.9887
## [2,] 0.8951 1.8951 1.9887
## [3,] 0.9887 1.9887 2.9887
## eigen() decomposition
## $values
## [1] 4.9475 0.5742 0.3206
##
## $vectors
## [,1] [,2] [,3]
## [1,] -0.3163 0.8309 0.4578
## [2,] -0.5812 0.2117 -0.7857
## [3,] -0.7498 -0.5146 0.4160
##
## $values
## [1] 4.9475 0.5742 0.3206
##
## $vectors
## [,1] [,2] [,3]
## [1,] 0.3163 -0.8309 0.4578
## [2,] 0.5812 -0.2117 -0.7857
## [3,] 0.7498 0.5146 0.4160
## [,1] [,2] [,3] [,4]
## [1,] 0.9066 0.9634 0.9044 0.8762
## [2,] 0.9634 1.9634 1.9044 1.8762
## [3,] 0.9044 1.9044 2.9044 2.8762
## [4,] 0.8762 1.8762 2.8762 3.8762
## eigen() decomposition
## $values
## [1] 7.9150 1.0569 0.4084 0.2703
##
## $vectors
## [,1] [,2] [,3] [,4]
## [1,] 0.2160 0.5473300 0.56388 -0.5795
## [2,] 0.4281 0.6083434 -0.07168 0.6644
## [3,] 0.5798 -0.0001828 -0.68124 -0.4469
## [4,] 0.6587 -0.5747592 0.46132 0.1516
##
## $values
## [1] 7.9150 1.0569 0.4084 0.2703
##
## $vectors
## [,1] [,2] [,3] [,4]
## [1,] 0.2160 -0.5473300 0.56388 -0.5795
## [2,] 0.4281 -0.6083434 -0.07168 0.6644
## [3,] 0.5798 0.0001828 -0.68124 -0.4469
## [4,] 0.6587 0.5747592 0.46132 0.1516