R入門

行列計算

最終更新：2026年06月23日 23:56

r-intro

- view

管理者のみ編集可

連立方程式

前進代入法で連立方程式を解く

forwardsolve関数を使う。以下は、連立方程式を行列で表記し、上三角行列の係数行列をA、定数項行列をy、求める解をxとおき、x = Ayと考えてxを求めた例。関数だけではなく、solve関数を使った例、逆行列を求めて正規方程式を解いた例も示してあるが、得られた解は同じであることがわかる。

> set.seed(7)
> mxaa <- matrix(trunc(runif(16, 1, 10)), 4, 4)
> mxaa[upper.tri(mxaa)] <- 0
> print(mxaa)
     [,1] [,2] [,3] [,4]
[1,]    9    0    0    0
[2,]    4    8    0    0
[3,]    2    4    2    0
[4,]    1    9    3    1
> mxy <- matrix(trunc(runif(4, 1, 10)), 4, 1)
> print(mxy)
     [,1]
[1,]    6
[2,]    1
[3,]    9
[4,]    3
> forwardsolve(mxaa, mxy)
           [,1]
[1,]  0.6666667
[2,] -0.2083333
[3,]  4.2500000
[4,] -8.5416667
> solve(mxaa, mxy)
           [,1]
[1,]  0.6666667
[2,] -0.2083333
[3,]  4.2500000
[4,] -8.5416667
> print(solve(t(mxaa) %*% mxaa) %*% (t(mxaa) %*% mxy))
           [,1]
[1,]  0.6666667
[2,] -0.2083333
[3,]  4.2500000
[4,] -8.5416667

forwardsolve関数は、与えられた行列の下三角の成分だけを使って解く。以下のとおり、すべての成分を使って解くsolve関数とは得られた解が違うことがわかる。

> set.seed(7)
> mxaa <- matrix(trunc(runif(16, 1, 10)), 4, 4)
> mxy <- matrix(trunc(runif(4, 1, 10)), 4, 1)
> print(mxaa)
     [,1] [,2] [,3] [,4]
[1,]    9    3    2    7
[2,]    4    8    5    1
[3,]    2    4    2    5
[4,]    1    9    3    1
> print(mxy)
     [,1]
[1,]    6
[2,]    1
[3,]    9
[4,]    3
> forwardsolve(mxaa, mxy)
           [,1]
[1,]  0.6666667
[2,] -0.2083333
[3,]  4.2500000
[4,] -8.5416667
> solve(mxaa, mxy)
            [,1]
[1,] -0.97516556
[2,]  0.06125828
[3,]  0.49337748
[4,]  1.94370861

後退代入法で連立方程式を解く

backsolve関数を使う。以下は連立方程式を行列で表記し、上三角行列の係数行列をA、定数項行列をy、求める解をxとおき、x = Ayと考えてxを求めた例。backsolve関数だけではなく、solve関数を使った例、逆行列を求めて正規方程式を解いた結果も示したが、得られた解は同じであることがわかる。

> set.seed(6)
> mxaa <- matrix(trunc(runif(16, 1, 10)), 4, 4)
> mxaa[lower.tri(mxaa)] <- 0
> print(mxaa)
     [,1] [,2] [,3] [,4]
[1,]    6    8    5    1
[2,]    0    9    1    3
[3,]    0    0    6    7
[4,]    0    0    0    3
> mxy <- matrix(trunc(runif(4, 1, 10)), 4, 1)
> print(mxy)
     [,1]
[1,]    5
[2,]    7
[3,]    2
[4,]    7
> backsolve(mxaa, mxy)
           [,1]
[1,]  2.0812757
[2,]  0.2654321
[3,] -2.3888889
[4,]  2.3333333
> solve(mxaa, mxy)
           [,1]
[1,]  2.0812757
[2,]  0.2654321
[3,] -2.3888889
[4,]  2.3333333
> print(solve(t(mxaa) %*% mxaa) %*% (t(mxaa) %*% mxy))
           [,1]
[1,]  2.0812757
[2,]  0.2654321
[3,] -2.3888889
[4,]  2.3333333

backsolve関数は、与えられた行列の上三角の成分だけを使って解く。以下のとおり、すべての成分を使って解くsolve関数とは得られた解が違うことがわかる。

> set.seed(6)
> mxaa <- matrix(trunc(runif(16, 1, 10)), 4, 4)
> mxy <- matrix(trunc(runif(4, 1, 10)), 4, 1)
> print(mxaa)
     [,1] [,2] [,3] [,4]
[1,]    6    8    5    1
[2,]    9    9    1    3
[3,]    3    9    6    7
[4,]    4    7    9    3
> print(mxy)
     [,1]
[1,]    5
[2,]    7
[3,]    2
[4,]    7
> backsolve(mxaa, mxy)
           [,1]
[1,]  2.0812757
[2,]  0.2654321
[3,] -2.3888889
[4,]  2.3333333
> solve(mxaa, mxy)
           [,1]
[1,]  1.8801598
[2,] -1.3715047
[3,]  0.8322237
[4,]  0.5299601

分解

コレスキー分解する

chol関数を使う。

> mx0 <- matrix(c(1, 0, 1, 0, 2, 0, 1, 0, 3), 3, byrow = TRUE)
> mx0
     [,1] [,2] [,3]
[1,]    1    0    1
[2,]    0    2    0
[3,]    1    0    3
> mx <- chol(mx0)
> mx
     [,1]     [,2]     [,3]
[1,]    1 0.000000 1.000000
[2,]    0 1.414214 0.000000
[3,]    0 0.000000 1.414214

得られた結果が正しいか、確認する。

> t(mx) %*% mx
     [,1] [,2] [,3]
[1,]    1    0    1
[2,]    0    2    0
[3,]    1    0    3

上記と同じ計算（A^T A）はcrossprod関数を使うと同じ結果を得ることができる。

> crossprod(mx)
     [,1] [,2] [,3]
[1,]    1    0    1
[2,]    0    2    0
[3,]    1    0    3

コレスキー分解を行う

chol関数を使う。コレスキー分解とは、既知の正定値対称行列Aに対してA = U^T Uを満たす上三角行列Uを求めること。

> aa <- matrix(c(1, 1, 1, 1, 2, 3, 1, 3, 6), 3, 3)
> print(aa)
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    1    2    3
[3,]    1    3    6
> det(aa)
[1] 1
> uu <- chol(aa)
> print(uu)
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    0    1    2
[3,]    0    0    1
> print(t(uu) %*% uu)
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    1    2    3
[3,]    1    3    6

特異値分解を行う

svd関数を使う。特異値分解の等式は以下のとおり。

A = U D V^T

UとVは直交行列。DはAの特異値を対角成分にもつ対角行列。

> a <- matrix(1:9, 3, 3)
> a
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9
> r <- svd(a)
> r
$d
[1] 1.684810e+01 1.068370e+00 5.543107e-16
$u
           [,1]        [,2]       [,3]
[1,] -0.4796712  0.77669099  0.4082483
[2,] -0.5723678  0.07568647 -0.8164966
[3,] -0.6650644 -0.62531805  0.4082483
$v
           [,1]       [,2]       [,3]
[1,] -0.2148372 -0.8872307  0.4082483
[2,] -0.5205874 -0.2496440 -0.8164966
[3,] -0.8263375  0.3879428  0.4082483

svd関数の戻り値はリスト。要素dに特異値が代入されている。

> r$d
[1] 1.684810e+01 1.068370e+00 5.543107e-16

上記の恒等式における、行列U（戻り値では要素u）と行列V（戻り値では要素v）が直交行列であることを確認する。

> t(r$u) %*% r$u
             [,1]          [,2]          [,3]
[1,] 1.000000e+00  1.110223e-16  5.551115e-17
[2,] 1.110223e-16  1.000000e+00 -1.110223e-16
[3,] 5.551115e-17 -1.110223e-16  1.000000e+00
> t(r$v) %*% r$v
              [,1]         [,2]          [,3]
[1,]  1.000000e+00 0.000000e+00 -1.110223e-16
[2,]  0.000000e+00 1.000000e+00  1.110223e-16
[3,] -1.110223e-16 1.110223e-16  1.000000e+00

上記の恒等式が成り立っているか確認する。

> r$u %*% diag(r$d) %*% t(r$v)
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9

行列Aの特異値は、A^T Aの正の固有値の正の平方根でも求められる。

> sqrt(eigen(t(a) %*% a)$values)
[1] 1.684810e+01 1.068370e+00 7.942757e-08

3つ目は実質0ということだが、高精度に計算されていないことがわかる。

LU分解を行う

Matrixパッケージのlu関数を使う。LU分解したそれぞれの行列を取り出すには、同パッケージに含まれているexpand関数を使う。

n次の正方行列Aが与えられたとき、下三角行列Lと上三角行列Uを用いてA = L Uと分解することを、行列AのLU分解という。lu関数は置換行列Pを用いてP A = L Uと分解する。戻り値はLが下三角行列、Uが上三角行列、Pが置換行列。以下は4次の正方行列のパスカル行列をLU分解した例。最後に、元の行列に戻るかどうか計算している。

> n <- 4
> aa <- matrix(1, n, n)
> for (i in 2:n) for (j in 2:n) aa[i, j] <- aa[i, j - 1] + aa[i - 1, j]
> print(aa)
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    1    2    3    4
[3,]    1    3    6   10
[4,]    1    4   10   20
> library(Matrix)
> r <- expand(lu(aa))
> print(r$L)
4 x 4 Matrix of class "dtrMatrix" (unitriangular)
     [,1]      [,2]      [,3]      [,4]     
[1,] 1.0000000         .         .         .
[2,] 1.0000000 1.0000000         .         .
[3,] 1.0000000 0.6666667 1.0000000         .
[4,] 1.0000000 0.3333333 1.0000000 1.0000000
> print(r$U)
4 x 4 Matrix of class "dtrMatrix"
     [,1]       [,2]       [,3]       [,4]      
[1,]  1.0000000  1.0000000  1.0000000  1.0000000
[2,]          .  3.0000000  9.0000000 19.0000000
[3,]          .          . -1.0000000 -3.6666667
[4,]          .          .          .  0.3333333
> print(r$P)
4 x 4 sparse Matrix of class "pMatrix"
[1,] | . . .
[2,] . . . |
[3,] . . | .
[4,] . | . .
> print(solve(r$P) %*% r$L %*% r$U)
4 x 4 Matrix of class "dgeMatrix"
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    1    2    3    4
[3,]    1    3    6   10
[4,]    1    4   10   20

LU分解を行う

matrixcalcパッケージのlu.decomposition関数を使う。n次の正方行列Aが与えられたとき、下三角行列Lと上三角行列Uを用いてA = L Uと分解することを、行列AのLU分解という。戻り値はLが下三角行列、Uが上三角行列。lu.decomposition関数はピボット選択は行わないため、戻り値に置換行列はない。以下は4次の正方行列のパスカル行列をLU分解した例。最後に、元の行列に戻るかどうか計算している。

> n <- 4
> aa <- matrix(1, n, n)
> for (i in 2:n) for (j in 2:n) aa[i, j] <- aa[i, j - 1] + aa[i - 1, j]
> print(aa)
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    1    2    3    4
[3,]    1    3    6   10
[4,]    1    4   10   20
> library(matrixcalc)
> r <- lu.decomposition(aa)
> print(r$L)
     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    1    1    0    0
[3,]    1    2    1    0
[4,]    1    3    3    1
> print(r$U)
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    0    1    2    3
[3,]    0    0    1    3
[4,]    0    0    0    1
> print(r$L %*% r$U)
     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    1    2    3    4
[3,]    1    3    6   10
[4,]    1    4   10   20

QR分解をする

行列XをQR分解で行列Qと行列Rに分解することを考える。

X = Q R

行列QはQ^T Q=I（Iは単位行列）となる直交行列。行列Rは上三角行列（対角成分はすべて0以外で、対角成分より下の要素はすべて0）。

qr.Q関数、qr.R関数を使うと、それぞれQR分解した際の行列Qと行列Rが得られる。引数にはqrオブジェクトを与える必要があり、qr関数を使ってあらかじめ分解したい行列Xからqrオブジェクトを作成しておく必要がある。

> d <- c(1, 2, 3, 4, 5, 6, 9, 9, 8)
> x <- matrix(d, 3, 3, byrow = TRUE)
> x
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
[3,]    9    9    8
> objqr <- qr(x)
> q <- qr.Q(objqr)
> q
           [,1]       [,2]       [,3]
[1,] -0.1010153 -0.7184092 -0.6882472
[2,] -0.4040610 -0.6025367  0.6882472
[3,] -0.9091373  0.3476173 -0.2294157
> r <- qr.R(objqr)
> r
          [,1]       [,2]        [,3]
[1,] -9.899495 -10.404571 -10.0005102
[2,]  0.000000  -1.320946  -2.9895090
[3,]  0.000000   0.000000   0.2294157

Q^T Q = I となるか試す。

> t(q) %*% q
              [,1]          [,2]          [,3]
[1,]  1.000000e+00 -1.110223e-16 -1.665335e-16
[2,] -1.110223e-16  1.000000e+00 -9.714451e-17
[3,] -1.665335e-16 -9.714451e-17  1.000000e+00

qr.X関数を使うことで、qrオブジェクトの元の行列を取り出すことができる。

> qr.X(objqr)
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
[3,]    9    9    8

QR分解を行う

qr関数を使う。分解した行列はqr関数の戻り値からそれぞれqr.Q関数、qr.R関数を使うことで得ることができる。

m行n列の行列A(m >= n)が与えられたとき、Q^T Q = Iを満たす行列Qと上三角行列（対角成分より下の成分はすべて0の行列）Rを用いてA = Q Rと分解することを、行列AのQR分解という。QR分解はcomplete型とreduced型の2種類があり、それぞれ分解した行列の行数と列数が異なる。

complete型 A(m×n) = Q(m×m) R(m×n)
reduced型 A(m×n) = Q(m×n) R(n×n)

以下は適当な行列を作成してQR分解した例。Q^T Q = Iを満たしていることと、Q R = Aと元に戻ることも計算している。qr.Q関数とqr.R関数にそれぞれcompleteオプションにTRUEを指定すると、complete型の計算を行う。何も指定しなければreduced型の計算結果が返される。

> m <- 4
> n <- 3
> aa <- matrix(0.0, m, n)
> for (i in 1:m) for (j in 1:n) aa[i, j] <- (-1) ^ i * choose(2 * i + j, i)
> print(aa)
     [,1] [,2] [,3]
[1,]   -3   -4   -5
[2,]   10   15   21
[3,]  -35  -56  -84
[4,]  126  210  330
> 
> # complete型 A(m×n) = Q(m×m) R(m×n)
> r <- qr(aa)
> qq <- qr.Q(r, complete = TRUE)  # m×m
> rr <- qr.R(r, complete = TRUE)  # m×n
> print(qq)
            [,1]       [,2]       [,3]      [,4]
[1,] -0.02286814  0.3344769  0.7843075 0.5219808
[2,]  0.07622713 -0.5474374 -0.2854080 0.7829713
[3,] -0.26679495  0.7243109 -0.5399946 0.3355591
[4,]  0.96046183  0.2526086 -0.1086731 0.0434984
> print(rr)
         [,1]       [,2]        [,3]
[1,] 131.1869 217.872381 341.0782902
[2,]   0.0000   2.936931   9.3501598
[3,]   0.0000   0.000000  -0.4176769
[4,]   0.0000   0.000000   0.0000000
> print(t(qq) %*% qq)
              [,1]          [,2]          [,3]          [,4]
[1,]  1.000000e+00 -8.326673e-17 -2.220446e-16 -7.632783e-17
[2,] -8.326673e-17  1.000000e+00  9.714451e-17  2.775558e-17
[3,] -2.220446e-16  9.714451e-17  1.000000e+00  1.257675e-16
[4,] -7.632783e-17  2.775558e-17  1.257675e-16  1.000000e+00
> print(qq %*% rr)
     [,1] [,2] [,3]
[1,]   -3   -4   -5
[2,]   10   15   21
[3,]  -35  -56  -84
[4,]  126  210  330
> 
> # reduced型 A(m×n) = Q(m×n) R(n×n)
> r <- qr(aa)
> qq <- qr.Q(r)  # m×n
> rr <- qr.R(r)  # n×n
> print(qq)
            [,1]       [,2]       [,3]
[1,] -0.02286814  0.3344769  0.7843075
[2,]  0.07622713 -0.5474374 -0.2854080
[3,] -0.26679495  0.7243109 -0.5399946
[4,]  0.96046183  0.2526086 -0.1086731
> print(rr)
         [,1]       [,2]        [,3]
[1,] 131.1869 217.872381 341.0782902
[2,]   0.0000   2.936931   9.3501598
[3,]   0.0000   0.000000  -0.4176769
> print(t(qq) %*% qq)
              [,1]          [,2]          [,3]
[1,]  1.000000e+00 -8.326673e-17 -2.220446e-16
[2,] -8.326673e-17  1.000000e+00  9.714451e-17
[3,] -2.220446e-16  9.714451e-17  1.000000e+00
> print(qq %*% rr)
     [,1] [,2] [,3]
[1,]   -3   -4   -5
[2,]   10   15   21
[3,]  -35  -56  -84
[4,]  126  210  330

qr関数にLAPACKオプションでTRUEを指定すると（デフォルトはFALSE)、計算にLAPACKを使いピボット選択が行われる（Using LAPACK (including in the complex case) uses column pivoting and does not attempt to detect rank-deficient matrices.）。このピボットの情報は、qr関数の戻り値のpivotというベクトルに含まれている。

> print(aa)
     [,1] [,2] [,3]
[1,]   -3   -4   -5
[2,]   10   15   21
[3,]  -35  -56  -84
[4,]  126  210  330
> r <- qr(aa, LAPACK = TRUE)
> print(r$pivot)
[1] 3 1 2
> qq <- qr.Q(r)
> rr <- qr.R(r)
> print(qq)
            [,1]       [,2]       [,3]
[1,] -0.01465387  0.2996579  0.7984525
[2,]  0.06154627 -0.5360454 -0.3095536
[3,] -0.24618510  0.7547242 -0.5071335
[4,]  0.96715574  0.2307639 -0.0972919
> print(rr)
         [,1]       [,2]        [,3]
[1,] 341.2067 131.137526 217.8708796
[2,]   0.0000  -3.598527  -3.0434558
[3,]   0.0000   0.000000   0.1310637
> print(t(qq) %*% qq)
             [,1]          [,2]          [,3]
[1,] 1.000000e+00  5.551115e-17  1.387779e-16
[2,] 5.551115e-17  1.000000e+00 -9.020562e-17
[3,] 1.387779e-16 -9.020562e-17  1.000000e+00
> print(aa[, r$pivot])
     [,1] [,2] [,3]
[1,]   -5   -3   -4
[2,]   21   10   15
[3,]  -84  -35  -56
[4,]  330  126  210
> print(qq %*% rr)
     [,1] [,2] [,3]
[1,]   -5   -3   -4
[2,]   21   10   15
[3,]  -84  -35  -56
[4,]  330  126  210

固有値と固有ベクトルを求める

eigen関数を使う。

> # 例1
> d <- c(1, 3, -2, -4)
> mx <- matrix(d, 2, 2)
> print(mx)
     [,1] [,2]
[1,]    1   -2
[2,]    3   -4
> eig <- eigen(mx)
> print(eig)
eigen() decomposition
$values
[1] -2 -1
$vectors
          [,1]      [,2]
[1,] 0.5547002 0.7071068
[2,] 0.8320503 0.7071068
> # 固有値（合計2つ）
> print(eig$value)
[1] -2 -1
> # 固有値 -2 のときの固有ベクトル
> print(eig$vector[, 1])
[1] 0.5547002 0.8320503
> # 固有値 -1 のときの固有ベクトル
> print(eig$vector[, 2])
[1] 0.7071068 0.7071068
> 
> # 例2
> d <- c(1, 1, 1, 0, 1, 0, 2, 1, 0)
> mx <- matrix(d, 3, 3)
> print(mx)
     [,1] [,2] [,3]
[1,]    1    0    2
[2,]    1    1    1
[3,]    1    0    0
> eig <- eigen(mx)
> print(eig)
eigen() decomposition
$values
[1]  2  1 -1
$vectors
          [,1] [,2]       [,3]
[1,] 0.5345225    0 -0.7071068
[2,] 0.8017837    1  0.0000000
[3,] 0.2672612    0  0.7071068
> # 固有値（合計3つ）
> print(eig$value)
[1]  2  1 -1
> # 固有値 2 のときの固有ベクトル
> print(eig$vector[, 1])
[1] 0.5345225 0.8017837 0.2672612
> # 固有値 1 のときの固有ベクトル
> print(eig$vector[, 2])
[1] 0 1 0
> # 固有値 -1 のときの固有ベクトル
> print(eig$vector[, 3])
[1] -0.7071068  0.0000000  0.7071068

固有値は数値的に求めているため、基本的に戻り値は近似値である。そのため、解析的には整数で得られる場合でも戻り値が複素数になる場合がある。この場合、固有ベクトルも複素数になる。以下の例では、解析的には固有値は3と-1（二重根）と求められるが、戻り値は複素数になっている。虚部が無視できるくらい小さいため、虚部を切り捨てて実数ととして扱ってかまわないが、ケースバイケースであることに注意。

> d <- c(1, 2, -1, 2, 1, -1, -4, 4, -1)
> mx <- matrix(d, 3, 3)
> eig <- eigen(mx)
> print(eig)
eigen() decomposition
$values
[1]  3+0i -1+0i -1-0i
$vectors
              [,1]                        [,2]                        [,3]
[1,] -0.9045340+0i -7.071068e-01+0.000000e+00i -7.071068e-01+0.000000e+00i
[2,] -0.3015113+0i  7.071068e-01-0.000000e+00i  7.071068e-01+0.000000e+00i
[3,]  0.3015113+0i  0.000000e+00+2.637113e-09i  0.000000e+00-2.637113e-09i

固有値と固有ベクトルを求める

eigen関数を使う。戻り値はリストで、順番にvaluesに固有値が、vectorsにその固有値のときの固有ベクトルが含まれる。

> mxaa <- matrix(c(3, 1, 1, 1, 2, 0, 1, 0, 2), nrow = 3)
> print(mxaa)
     [,1] [,2] [,3]
[1,]    3    1    1
[2,]    1    2    0
[3,]    1    0    2
> lam <- eigen(mxaa)
> print(lam)
eigen() decomposition
$values
[1] 4 2 1
$vectors
          [,1]       [,2]       [,3]
[1,] 0.8164966  0.0000000  0.5773503
[2,] 0.4082483 -0.7071068 -0.5773503
[3,] 0.4082483  0.7071068 -0.5773503

この例では固有値と固有ベクトルは3つあり、例えば3番目の固有値と固有ベクトルを取り出すには、以下のようにする。

> lam$values[3]
[1] 1
> lam$vectors[, 3]
[1]  0.5773503 -0.5773503 -0.5773503

「行列計算」をウィキ内検索

最近更新されたスレッド

+新規スレッドを作成する

人気記事ランキング

最近更新されたページ

急上昇Wikiランキング

急上昇中のWikiランキングです。今注目を集めている話題をチェックしてみよう！

新規Wikiランキング

最近作成されたWikiのアクセスランキングです。見るだけでなく加筆してみよう！

人気Wikiランキング

atwikiでよく見られているWikiのランキングです。新しい情報を発見してみよう！

全体ページランキング

最近アクセスの多かったページランキングです。話題のページを見に行こう！

行列計算

目次

連立方程式

前進代入法で連立方程式を解く

後退代入法で連立方程式を解く

分解

コレスキー分解する

コレスキー分解を行う

特異値分解を行う

LU分解を行う

LU分解を行う

QR分解をする

QR分解を行う

固有値と固有ベクトルを求める

固有値と固有ベクトルを求める

メニュー

更新履歴