[R]行列関連の計算

R の初歩的な行列関連の計算。僕自身は行列操作は手計算できるようなものしかやったことがない。

行列の積
%*% を使う。なぜこの演算子が与えられたのかは、不明。

> A <- matrix(c(1,3,2,-4,-1,2), 2, 3) > A
     [,1] [,2] [,3]
[1,]    1    2   -1
[2,]    3   -4    2
> B <- matrix(c(2,1,4,1,-1,-3), 3, 2) > B
     [,1] [,2]
[1,]    2    1
[2,]    1   -1
[3,]    4   -3
> A %*% B
     [,1] [,2]
[1,]    0    2
[2,]   10    1
> B %*% A
     [,1] [,2] [,3]
[1,]    5    0    0
[2,]   -2    6   -3
[3,]   -5   20  -10

行列式(determinant)
det を使う

A <- matrix(c(1,0,0,1), 2, 2) > A
     [,1] [,2]
[1,]    1    0
[2,]    0    1
> det(A)
> det(A)
[1] 7

> B <- matrix(c(2,4,8,1,6,8,1,3,9), 3, 3) > B
     [,1] [,2] [,3]
[1,]    2    1    1
[2,]    4    6    3
[3,]    8    8    9
> det(B)
> det(B)
[1] 32

逆行列(invertible matrices)
solve を使う

> solve(A)
     [,1] [,2]
[1,]    1    0
[2,]    0    1

> B
     [,1] [,2] [,3]
[1,]    2    1    1
[2,]    4    6    3
[3,]    8    8    9

固有値・固有ベクトル(eigenvalues/eigenvectors)

> A <- matrix(c(3,2, 1,2), 2,2) > A
     [,1] [,2]
[1,]    3    1
[2,]    2    2
> B <- eigen(A)
> B
$values
[1] 4 1

$vectors
          [,1]       [,2]
[1,] 0.7071068 -0.4472136
[2,] 0.7071068  0.8944272

> (A - B$values[1] * diag(1, 2))
     [,1] [,2]
[1,]   -1    1
[2,]    2   -2
> B$vectors[,1]
[1] 0.7071068 0.7071068
> (A - B$values[1] * diag(1, 2)) %*% cbind(B$vectors[,1])
     [,1]
[1,]    0
[2,]    0

> (A - B$values[2] * diag(1, 2))
     [,1] [,2]
[1,]    2    1
[2,]    2    1
> B$vectors[,2]
[1] -0.4472136  0.8944272
> (A - B$values[2] * diag(1, 2)) %*% cbind(B$vectors[,2])
     [,1]
[1,]    0
[2,]    0

> A <- matrix(c(5,-1,4,1), 2, 2) > A
     [,1] [,2]
[1,]    5    4
[2,]   -1    1
> eigen(A)
$values
[1] 3 3

$vectors
           [,1]       [,2]
[1,]  0.8944272 -0.8944272
[2,] -0.4472136  0.4472136
> B <- eigen(A) > B
$values
[1] 3 3

$vectors
           [,1]       [,2]
[1,]  0.8944272 -0.8944272
[2,] -0.4472136  0.4472136

> x <- B$vectors[,1]
> (A - B$values[1] * diag(1, 2)) %*% cbind(x)
                 x
[1,] -2.220446e-16
[2,]  1.014271e-16

ランク(rank)

グーグると、直接求める方法は見当たらず、QR分解したときの rank 属性を利用するものが見つかった。(QR分解は数分前に名前を聞いただけなので、理論的なことは知らない)

> A
     [,1] [,2]
[1,]    3    1
[2,]    2    2
> qr(A)$rank
[1] 2
> C <- matrix(c(1,5,9,13, 2,6,10,14,3,7,11,15,4,8,12,16), 4, 4) > C
     [,1] [,2] [,3] [,4]
[1,]    1    2    3    4
[2,]    5    6    7    8
[3,]    9   10   11   12
[4,]   13   14   15   16
> qr(C)$rank
[1] 2

wikipedia にのっている Gram-Schmidt process の例を R で計算すると次のようになる

> A <- matrix(c(12, 6, -4, -51, 167, 24, 4, -68, -41), 3, 3) > A
     [,1] [,2] [,3]
[1,]   12  -51    4
[2,]    6  167  -68
[3,]   -4   24  -41

> foo <- qr(A) > foo
$qr
            [,1]         [,2] [,3]
[1,] -14.0000000  -21.0000000   14
[2,]   0.4285714 -175.0000000   70
[3,]  -0.2857143    0.1107692  -35

$rank [1] 3

$qraux
[1]  1.857143  1.993846 35.000000

$pivot
[1] 1 2 3

attr(,"class")
[1] "qr"
> Q <- qr.Q(foo) > R <- qr.R(foo) > qr.solve(foo)
             [,1]         [,2]         [,3]
[1,]  0.060816327 0.0232653061 -0.032653061
[2,] -0.006040816 0.0055510204 -0.009795918
[3,] -0.009469388 0.0009795918 -0.026938776
> qr.X(foo)
     [,1] [,2] [,3]
[1,]   12  -51    4
[2,]    6  167  -68
[3,]   -4   24  -41

▼Q は直行行列

> t(Q) %*% Q
              [,1]          [,2]          [,3]
[1,]  1.000000e+00 -2.541099e-19 -3.081845e-17
[2,] -2.541099e-19  1.000000e+00  2.467915e-17
[3,] -3.081845e-17  2.467915e-17  1.000000e+00

> Q %*% t(Q)
              [,1]          [,2]          [,3]
[1,]  1.000000e+00 -1.051439e-16 -2.217193e-17
[2,] -1.051439e-16  1.000000e+00  4.248378e-17
[3,] -2.217193e-17  4.248378e-17  1.000000e+00

▼次の3つは同じ結果が得られる

> t(Q) %*% A
              [,1]          [,2] [,3]
[1,] -1.400000e+01 -2.100000e+01   14
[2,] -2.220446e-16 -1.750000e+02   70
[3,]  0.000000e+00  1.776357e-15  -35
> t(Q) %*% Q %*% R
              [,1]          [,2] [,3]
[1,] -1.400000e+01 -2.100000e+01   14
[2,]  3.557538e-18 -1.750000e+02   70
[3,]  4.314583e-16 -3.671664e-15  -35
> R
     [,1] [,2] [,3]
[1,]  -14  -21   14
[2,]    0 -175   70
[3,]    0    0  -35
Tagged with: , , , , ,
Posted in algorithm, R

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Archives
  • RT @__apf__: How to write a research paper: a guide for software engineers & practitioners. docs.google.com/presentation/d… /cc @inwyrd 6 days ago
  • RT @HayatoChiba: 昔、自然と対話しながら数学に打ち込んだら何かを悟れるのではと思いたち、専門書1つだけ持ってパワースポットで名高い奈良の山奥に1週間籠ったことがある。しかし泊まった民宿にドカベンが全巻揃っていたため、水島新司と対話しただけで1週間過ぎた。 それ… 3 weeks ago
  • RT @googlecloud: Ever wonder what underwater fiber optic internet cables look like? Look no further than this deep dive w/ @NatAndLo: https… 3 weeks ago
  • @ijin UTC+01:00 な時間帯で生活しています、、、 6 months ago
  • RT @mattcutts: Google's world-class Site Reliability Engineering team wrote a new book: amazon.com/Site-Reliabili… It's about managing produc… 9 months ago
%d bloggers like this: