固有値分解と行列の階乗

一般線形最小二乗法を使うために行列を固有値分解して、-0.5乗などの行列を求める必要があり、このあたりをどうするかを調べました。 Rには、eigenというFunctionがあります。

> A = matrix(c(1,2,3,2,3,4,3,4,5), ncol = 3)
> A
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    2    3    4
[3,]    3    4    5
> eigen(A)
eigen() decomposition
$values
[1]  9.623475e+00 -1.887379e-15 -6.234754e-01

$vectors
           [,1]       [,2]       [,3]
[1,] -0.3850898  0.4082483  0.8276709
[2,] -0.5595102 -0.8164966  0.1424137
[3,] -0.7339306  0.4082483 -0.5428436

この出力を使って、

この形に分解してくれます。

> P = eigen(A)$vector
> D = diag(eigen(A)$values)
> invP = solve(P)
> P %*% D %*% invP
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    2    3    4
[3,]    3    4    5

なお、Aが分散共分散行列などのような対称行列であれば、

なので、solve(A) を t(A)にすることで計算量は減ります。

とにかく、固有値分解をする意味は、行列の階乗が簡単にできることです。例えばAの3乗は

> A %*% A %*% A
     [,1] [,2] [,3]
[1,]  132  192  252
[2,]  192  279  366
[3,]  252  366  480

> P %*% D^3 %*% invP
     [,1] [,2] [,3]
[1,]  132  192  252
[2,]  192  279  366
[3,]  252  366  480

で計算できます。下記のようなFunctionをつくっておくのも良いかもしれません。

> "%^%" <- function(x, n){
+   with(eigen(x), vectors %*% (values^n * solve(vectors)))
+ }
> A %^% 3
     [,1] [,2] [,3]
[1,]  132  192  252
[2,]  192  279  366
[3,]  252  366  480

この関数で、Aの半乗（0.5乗）を求めました。

> A = matrix(c(2,2,3,4,5,3,2,3,5), ncol = 3)
> A
     [,1] [,2] [,3]
[1,]    2    4    2
[2,]    2    5    3
[3,]    3    3    5
> B = A %^% 0.5
> B
          [,1]      [,2]      [,3]
[1,] 1.0577617 1.2406827 0.3642050
[2,] 0.4495838 2.0309631 0.6886054
[3,] 0.8878226 0.4609286 2.0878823
> B %*% B
     [,1] [,2] [,3]
[1,]    2    4    2
[2,]    2    5    3
[3,]    3    3    5