R(行列)

2023年12月15日 (金)

[R]行列に列を追加する

cbind関数を使う。最後の例のとおり、行数が異なる場合はエラーが発生する。

> mx1 <- matrix(1:9, nrow = 3)
> mx2 <- matrix(10:15, nrow = 3)
> print(mx1)
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
> print(mx2)
[,1] [,2]
[1,] 10 13
[2,] 11 14
[3,] 12 15
> mx3 <- cbind(mx1, mx2)
> print(mx3)
[,1] [,2] [,3] [,4] [,5]
[1,] 1 4 7 10 13
[2,] 2 5 8 11 14
[3,] 3 6 9 12 15
> mx4 <- matrix(16:19, nrow = 4)
> print(mx4)
[,1]
[1,] 16
[2,] 17
[3,] 18
[4,] 19
> cbind(mx1, mx4)
cbind(mx1, mx4) でエラー:
行列の行数は一致していなければなりません (2 番目の引数を参照)

2023年12月13日 (水)

[R]行列に行を追加する

rbind関数を使う。最後の例のとおり、列数が異なる場合はエラーが発生する。

> mx1 <- matrix(1:9, ncol = 3)
> mx2 <- matrix(10:15, ncol = 3)
> print(mx1)
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
> print(mx2)
[,1] [,2] [,3]
[1,] 10 12 14
[2,] 11 13 15
> mx3 <- rbind(mx1, mx2)
> print(mx3)
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
[4,] 10 12 14
[5,] 11 13 15
> mx4 <- matrix(16:19, ncol = 4)
> print(mx4)
[,1] [,2] [,3] [,4]
[1,] 16 17 18 19
> rbind(mx1, mx4)
rbind(mx1, mx4) でエラー:
行列の列数は一致していなければなりません (2 番目の引数を参照)

2023年10月13日 (金)

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

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

2023年2月13日 (月)

[R]連立一次方程式を解く

solve関数を使う。以下の連立一次方程式を解く(解はx=3,y=2)。

 x + 2y = 7
3x - 4y = 1

計算する。

> mxaa <- matrix(c(1, 2, 3, -4), 2, 2, byrow = TRUE)
> mxy <- matrix(c(7, 1), 2, 1)
> print(mxaa)
[,1] [,2]
[1,] 1 2
[2,] 3 -4
> print(mxy)
[,1]
[1,] 7
[2,] 1
> solve(mxaa, mxy)
[,1]
[1,] 3
[2,] 2

同じく以下の連立一次方程式を解く(解はx=2,y=-3,z=-14)。

 3x - 3y + z =  1
3x + 2y = 0
-1x - 5y + z = -1

計算する。

> mxaa <- matrix(c(3, 3, -1, -3, 2, -5, 1, 0, 1), 3, 3)
> mxy <- matrix(c(1, 0, -1), 3, 1)
> print(mxaa)
[,1] [,2] [,3]
[1,] 3 -3 1
[2,] 3 2 0
[3,] -1 -5 1
> print(mxy)
[,1]
[1,] 1
[2,] 0
[3,] -1
> solve(mxaa, mxy)
[,1]
[1,] 2
[2,] -3
[3,] -14

2023年2月 9日 (木)

[R]一般逆行列を求める

MASSパッケージのginv関数を使う。

> library(MASS)
> mx <- matrix(1:16, 4, 4, byrow = TRUE)
> print(mx)
[,1] [,2] [,3] [,4]
[1,] 1 2 3 4
[2,] 5 6 7 8
[3,] 9 10 11 12
[4,] 13 14 15 16
> det(mx)
[1] 4.733165e-30
> solve(mx)
solve.default(mx) でエラー:
システムは数値的に特異です: 条件数の逆数 = 2.95756e-18
> ginv(mx)
[,1] [,2] [,3] [,4]
[1,] -0.2850 -0.1450 -0.0050 0.1350
[2,] -0.1075 -0.0525 0.0025 0.0575
[3,] 0.0700 0.0400 0.0100 -0.0200
[4,] 0.2475 0.1325 0.0175 -0.0975

一般逆行列の定義を満たしているか否か確認。

> mxi <- ginv(mx)
> mx %*% mxi %*% mx
[,1] [,2] [,3] [,4]
[1,] 1 2 3 4
[2,] 5 6 7 8
[3,] 9 10 11 12
[4,] 13 14 15 16

当然、正方行列でなくとも求めることができる。

> mx <- matrix(1:8, 2, 4, byrow = TRUE)
> print(mx)
[,1] [,2] [,3] [,4]
[1,] 1 2 3 4
[2,] 5 6 7 8
> det(mx)
determinant.matrix(x, logarithm = TRUE, ...) でエラー:
'x' は正方行列でなければなりません
> solve(mx)
solve.default(mx) でエラー: 'a' (2 x 4) は正方行列である必要があります
> ginv(mx)
[,1] [,2]
[1,] -0.550 2.500000e-01
[2,] -0.225 1.250000e-01
[3,] 0.100 -2.081668e-17
[4,] 0.425 -1.250000e-01
> mxi <- ginv(mx)
> mx %*% mxi %*% mx
[,1] [,2] [,3] [,4]
[1,] 1 2 3 4
[2,] 5 6 7 8

ginv関数はムーア・ペンローズ一般逆行列を数値的に求める関数のため、求まった一般逆行列は近似値であることに注意。ムーア・ペンローズ一般逆行列の定義の4つの式について、比較をした例は以下のとおり。すべての要素がTRUEにならないことがわかる。

> mx %*% mxi %*% mx == mx
[,1] [,2] [,3] [,4]
[1,] FALSE FALSE TRUE TRUE
[2,] FALSE FALSE FALSE TRUE
> mxi %*% mx %*% mxi == mxi
[,1] [,2]
[1,] FALSE FALSE
[2,] FALSE FALSE
[3,] FALSE FALSE
[4,] FALSE FALSE
> t(mx %*% mxi) == mx %*% mxi
[,1] [,2]
[1,] TRUE FALSE
[2,] FALSE TRUE
> t(mxi %*% mx) == mxi %*% mx
[,1] [,2] [,3] [,4]
[1,] TRUE FALSE FALSE FALSE
[2,] FALSE TRUE FALSE FALSE
[3,] FALSE FALSE TRUE FALSE
[4,] FALSE FALSE FALSE TRUE

2023年2月 7日 (火)

[R]逆行列を求める

solve関数を使う。

> mx <- matrix(3, 1, 1)
> print(mx)
[,1]
[1,] 3
> solve(mx)
[,1]
[1,] 0.3333333
> mx <- matrix(1:4, 2, 2, byrow = TRUE)
> print(mx)
[,1] [,2]
[1,] 1 2
[2,] 3 4
> solve(mx)
[,1] [,2]
[1,] -2.0 1.0
[2,] 1.5 -0.5
> mx <- matrix(c(3, 3, -1, -3, 2, -5, 1, 0, 1), 3, 3)
> print(mx)
[,1] [,2] [,3]
[1,] 3 -3 1
[2,] 3 2 0
[3,] -1 -5 1
> solve(mx)
[,1] [,2] [,3]
[1,] 1.0 -1 -1.0
[2,] -1.5 2 1.5
[3,] -6.5 9 7.5

正則ではない行列(行列式の値が0)の場合は、solve関数はエラーを返す。

> mx <- matrix(1:9, 3, 3)
> print(mx)
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
> det(mx)
[1] 0
> solve(mx)
solve.default(mx) でエラー:
Lapack routine dgesv: システムは正確に特異です: U[3,3] = 0

2023年2月 6日 (月)

[R]転置行列を求める

t関数を使う。

> mx <- matrix(1:9, 3, 3, byrow = TRUE)
> print(mx)
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6
[3,] 7 8 9
> t(mx)
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
> mx <- matrix(1:6, 2, 3, byrow = TRUE)
> print(mx)
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6
> t(mx)
[,1] [,2]
[1,] 1 4
[2,] 2 5
[3,] 3 6

2023年2月 4日 (土)

[R]行列式を求める

det関数を使う。

> mx <- matrix(3, 1, 1)
> print(mx)
[,1]
[1,] 3
> det(mx)
[1] 3
> mx <- matrix(1:4, 2, 2, byrow = TRUE)
> print(mx)
[,1] [,2]
[1,] 1 2
[2,] 3 4
> det(mx)
[1] -2
> mx <- matrix(c(0, -4, 0, -1, 3, 6, -2, 0, 2, 1, 3, 0, -1, 5, -1, 2), 4, 4)
> print(mx)
[,1] [,2] [,3] [,4]
[1,] 0 3 2 -1
[2,] -4 6 1 5
[3,] 0 -2 3 -1
[4,] -1 0 0 2
> det(mx)
[1] 28

2023年2月 2日 (木)

[R]行列の掛け算

%*% 演算子を使う。最後の例のとおり、掛けられる行列の列数と掛ける行列の行数が同じでなければ計算することはできない。

> mx1 <- matrix(c(-1, 8, 0, 1, 3, -5), 2, 3)
> mx2 <- matrix(c(3, -2, 4, 1, 0, 2), 2, 3)
> mx3 <- matrix(c(4, 2, 1, -3, 0, 5, 7, -1, 0), 3, 3)
> print(mx1)
[,1] [,2] [,3]
[1,] -1 0 3
[2,] 8 1 -5
> print(mx2)
[,1] [,2] [,3]
[1,] 3 4 0
[2,] -2 1 2
> print(mx3)
[,1] [,2] [,3]
[1,] 4 -3 7
[2,] 2 0 -1
[3,] 1 5 0
> mx1 %*% mx3
[,1] [,2] [,3]
[1,] -1 18 -7
[2,] 29 -49 55
> mx2 %*% mx3
[,1] [,2] [,3]
[1,] 20 -9 17
[2,] -4 16 -15
> mx1 %*% mx2
mx1 %*% mx2 でエラー: 適切な引数ではありません

2022年12月24日 (土)

[R]余因子行列を求める

polyMatrixパッケージのadjoint関数を使う。以下、使用例。計算に用いた行列はここhttps://jp.mathworks.com/help/symbolic/adjoint.htmlを参考にしている。

> d <- c(8, 3, 4, 1, 5, 9, 6, 7, 2)
> mx <- matrix(d, 3)
> mx
[,1] [,2] [,3]
[1,] 8 1 6
[2,] 3 5 7
[3,] 4 9 2
> adjoint(mx)
[,1] [,2] [,3]
[1,] -53 52 -23
[2,] 22 -8 -38
[3,] 7 -68 37

余因子行列の性質 A A~ = A~ A ~ = |A|・En を確認する。

> print(mx %*% adjoint(mx))
[,1] [,2] [,3]
[1,] -3.600000e+02 2.842171e-13 -5.684342e-14
[2,] -7.815970e-14 -3.600000e+02 -5.684342e-14
[3,] -8.704149e-14 1.421085e-13 -3.600000e+02
> print(adjoint(mx) %*% mx)
[,1] [,2] [,3]
[1,] -3.600000e+02 2.842171e-14 -1.136868e-13
[2,] 5.684342e-14 -3.600000e+02 5.684342e-14
[3,] 8.526513e-14 1.705303e-13 -3.600000e+02
> print(det(mx) * diag(nrow(mx)))
[,1] [,2] [,3]
[1,] -360 0 0
[2,] 0 -360 0
[3,] 0 0 -360

上2つの計算結果の非対角成分は非常に小さいので0と見なしてよく、上記の性質を満たしていることが確認できた。