20. 行列計算


行列の計算

行列の和,差,積 などの基本的な行列演算は普通の数値と同じコマンドで行うことが出来る.

 
 a <- matrix(1:4, 2, 2)                # 2 * 2 行列
 b <- matrix(0:3, 2, 2)                # 2 * 2 行列
 a + b - a %*% b                       # 和・差・積
 a * b                                 # 要素ごとの掛け算
 a %o% b                               # 外積(関数 outer(x,y) と同じ)
 a %x% b                               # クロネッカー積

クロネッカー積は関数 kronecker() でも計算できる.

行列の積は %*% 演算子によって求める点に注意.行列に対して * 演算子を用いると要素毎の積を計算してしまう.

関数 outer(x, y, FUN) の引数 FUN に 2 変数関数を指定することで,外積を求める際の関数を指定することが出来る.

 
 x <- seq(-3,3,0.1)
 y <- seq(-3,3,0.1)
 f <- function(x,y) 1/(2*pi)*exp(-(x^2+y^2)/2)
 z <- outer(x, y, f)                   # この後,persp(x,y,z,expand=0.5) などで3次元のグラフが描ける

演算子 * については,定数部分にベクトルを持ってくることも出来る.また,似たような働きで / を用いることで,要素ごとの逆数を求めることも出来る.

 
 a <- matrix(c(1,2,3,4), 2, 2)         # 2 * 2 行列
 a * (1:2)                             # (1:2) * a も同じ結果

     [,1] [,2]
[1,]    1    3
[2,]    4    8

 1/a

     [,1]      [,2]
[1,]  1.0 0.3333333
[2,]  0.5 0.2500000

行列計算用関数の一覧

R には行列操作を行う関数が多数用意されている.以下に一覧表を示す.

コマンド

機能

matrix(0, nrow=2, ncol=3)

2 行 3 列のゼロ行列を作る.matrix(1,nrow=2,ncol=3) と応用することも可.

diag(3)

3 × 3 の単位行列を作る.diag(1, ncol=2,  nrow=2),diag(1, 3),diag(rep(1, 3)),diag(c(1,1)) などでも可.

diag(1:3)

対角成分が (1 2 3) の 3 × 3 対角行列を作る.

diag(X) <- 1

行列 X の対角成分を全て 1 にする.

diag(X) <- c(1, 2)

2 × 2 行列 X の対角成分を (1 2) にする.

t(X)

行列 X を転置する.

X[upper.tri(X)] <- 0

上三角成分(対角成分含まず)を全て 0 にする.対角成分も 0 にする場合は upper.tri(X,diag=TRUE) とする.

X[lower.tri(X)] <- 0

下三角成分(対角成分含まず)を全て 0 にする.

sum(X^2)

行列成分の二乗.sum( diag( t(X) %*% X ) ) でも可.

solve(X)

行列 X の逆行列を求める.

ginv(X)

行列 X のムーア・ペンローズ型一般化逆行列を求める( library(MASS) を実行する必要あり).

crossprod(X)

行列 X 自身のクロス積を求める(t(X) %*% X).

crossprod(X, Y)

行列 X と Y のクロス積を求める(t(X) %*% Y).

eigen(X)

行列 X の固有値と固有ベクトルを求める.

qr(X)

行列 X の QR 分解を行う.

chol(X)

正値対称行列(エルミート行列)X のコレスキー分解を行う.すなわち,X = t(A) %*% A を満たす上三角行列 A を求める.

svd(X)

行列 X の特異値分解を行なう.すなわち,X = U %*% D %*% t を満たす「直交行列 U 」,「X の特異値を対角成分とする対角行列 D 」,「直交行列 V 」を求める.

det(X)

行列 X の行列式を求める.prod( svd(X)$d ) でも可.

一覧表では伝え切れなかった事項を以下で述べる.

対称行列

三角行列から対称行列を生成することが出来る.

 
 ( x <- matrix(c(1,2,0,3), 2, 2) )

     [,1] [,2]
[1,]    1    0
[2,]    2    3

 y <- x+t(x)             # 三角行列+三角行列の転置行列
 diag(y) <- diag(y)/2    # 対角成分を元に戻す
 y                       # y に対称行列が出来ている

     [,1] [,2]
[1,]    1    2
[2,]    2    3

連立方程式の解

A.x = b なる形の x についての連立方程式の解は関数 solve() を使って解く事が出来る.

 
 a <- matrix(c(0,1,2,3,4,5,6,7,9), 3,3)       #      3y + 6z =  1
 b <- matrix(c(1,0,-2))                       #  x + 4y + 7z =  0
 solve(a,b)                                   # 2x + 5y + 9z = -2

          [,1]
[1,] -2.333333
[2,]  2.333333
[3,] -1.000000

係数行列が上三角行列の場合は関数 backsolve() を使う.backsolve() は係数行列の下三角部分を無視するので下三角部分が 0 でない行列を係数行列に指定することも出来るが,この場合に連立一次方程式の係数として使われるのは上三角の部分のみとなる.

 
 r <- rbind(c(1,2,3),                         # x + 2y + 3z = 8
            c(0,1,1),                         #      y +  z = 4
            c(0,0,2))                         #          2z = 2
 ( y <- backsolve(r, x <- c(8,4,2)) )         #
 r %*% y                                      # 結果は (8   4  2)
 ( y2 <- backsolve(r, x, transpose = TRUE))   # 結果は (8 -12 -5)
 all(t(r) %*% y2 == x)
 all(y  == backsolve(t(r), x, upper = FALSE, transpose = TRUE) )
 all(y2 == backsolve(t(r), x, upper = FALSE, transpose = FALSE))

下三角行列の場合は関数 forwardsolve() を用いる.

固有値と固有ベクトル

行列の固有値  ( 実正方行列の固有値) と固有ベクトルは関数 eigen() で求められる.このとき,行列のスペクトル分解がリストの成分として返される.

 
 X <- matrix(c(1,2,3,4), 2, 2)
 z <- eigen(X)                                # X の固有値を求める
 z$values

[1]  5.3722813 -0.3722813

 z$vectors[,1]                                # 固有値  5.3722813 に対する固有ベクトル

[1] -0.5657675 -0.8245648

 z$vectors[,2]                                # 固有値 -0.3722813 に対する固有ベクトル

[1] -0.9093767  0.4159736

行列の平方根

関数 svd() を用いると,行列の平方根を求めることが出来る.

 
 A <- array(c(2,1,1,1,2,1,1,1,2), dim=c(3,3))
 U <- svd(A)$u
 V <- svd(A)$v
 D <- diag(sqrt(svd(A)$d))
 B <- U %*% D %*% t(V)              # 行列 A の平方根
 A                                  # 行列 A を表示

     [,1] [,2] [,3]
[1,]    2    1    1
[2,]    1    2    1
[3,]    1    1    2

 B                                  # 行列 B を表示

          [,1]      [,2]      [,3]
[1,] 1.3333333 0.3333333 0.3333333
[2,] 0.3333333 1.3333333 0.3333333
[3,] 0.3333333 0.3333333 1.3333333

 B%*%B                              # 検算(行列 A と一致している)

     [,1] [,2] [,3]
[1,]    2    1    1
[2,]    1    2    1
[3,]    1    1    2

正方行列のべき乗

例えば正方行列 A に対する 2 乗操作 A2 は成分毎に 2 乗した行列を与えるだけで,行列のべき乗 A %*% A は与えてくれない.そこで,A のべき乗 An を計算する場合は,

を用いて

を計算すればよい.D は対角行列なので Dn と D の行列としての n 乗は一致する.以下では行列の指数乗を計算する関数 matpow() を定義している(関数 matpow() は 『RjpWiki』 の記事より引用した.).

 
 matpow <- function(x, pow=2) {
   y <- eigen(x)
   y$vectors %*% diag( (y$values)^pow ) %*% t(y$vectors)
 }
 matpow(diag(1:2))         # 行列 ((1, 0) , (0, 2)) の 2 乗

     [,1] [,2]
[1,]    1    0
[2,]    0    4

この概念を応用すると行列の指数乗が計算出来る.以下では行列の指数乗を計算する関数 matexp() を定義している.

 
 matexp <- function(x) {
   y <- eigen(x)
   y$vectors %*% diag( exp(y$values) ) %*% t(y$vectors)
 }
 matexp(matrix(0, 2, 2))   # ゼロ行列の指数乗は単位行列

     [,1] [,2]
[1,]    1    0
[2,]    0    1

 matexp(diag(2))           # exp(E) = ((e, 0) , (0, e))

         [,1]     [,2]
[1,] 2.718282 0.000000
[2,] 0.000000 2.718282