R-Matrix

行列の作成

R では以下の手順で行列を作成する．

行列の要素をベクトル（配列やリスト，行列でも可）で用意する
関数 matrix(ベクトル, 行数, 列数) でベクトルから行列に変換する

例として，ベクトル（1，2，3，4，5，6）を変換して，行列 $\begin{pmatrix}1$ を作成する．

 
 matrix(1:6, nrow=2, ncol=3)           # nrow で行数，ncol で列数を指定する 
                                       # matrix(1:6, 2, 3) と略記しても良い 
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

行列を作成する場合，指定された要素は左の列から順に (上から下へ) 埋められる．要素を上の行から順に (左から右へ) 埋める場合は，引数 byrow=T を指定する．

 
 matrix(1:6, nrow=2, ncol=3, byrow=T)  # matrix(1:6, 2, 3, byrow=T)でも可 
 
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6

行列要素へのアクセス

行列の中の数を「要素」と呼び，行は左から順に 1, 2, … と，列は上から順に 1, 2, … と番号が振られている．以下では行列 x について要素にアクセスする方法を一覧表で示している．

基本は，[ と ] の間にカンマで区切って行の位置（番号）と列の位置（番号）を指定ればよい．

 
 ( x <- matrix(1:6, nrow=2, ncol=3) )
 
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6 
 
 x[1, 2]
 
[1] 3 
 
 x[-1, c(T,F,T)]                       # 1 行目を非表示，2 行目 1,3 成分を指定 
 
[1] 2 6

以下に指定例の一覧を示す．

コマンド	機能
x[2, ]	2 行目を取り出す
x[ ,2]	2 列目を取り出す
x[1, 2]	1 行 2 列目の成分を取り出す
x[c(1,2), 2]	1,2 行 2 列目の成分を取り出す
x[c(1,2), c(1,3)]	1, 2 行目と 1, 3 列を取り出す
x[ ,-c(1, 3)]	1, 3列を除外した行列を取り出す
x[ , c(T,F,T)]	1, 3列を取り出す
rowSums()	行の総和
colSums()	列の総和
rowMeans()	行の平均
colMeans()	列の平均

行列の計算

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

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

行列計算用関数の一覧

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 乗操作 A² は成分毎に 2 乗した行列を与えるだけで，行列のべき乗 A %*% A は与えてくれない．そこで，A のべき乗 Aⁿ を計算する場合は，

A の固有値分解 A = V %*% D %*% t(V)

を用いて

A^n = V %*% D^n %*% t(V)

を計算すればよい．D は対角行列なので Dⁿと 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

Back to R