Rで線形代数 Rで数学の色々をいじってみる

• ベクトル、行列、行列の多次元化→array

v<-floor(runif(10)*10)
# 行列を作るには、行数と列数を指定する
n.row<-3
n.col<-4
m<-matrix(floor(runif(n.row*n.col)*10),n.row,n.col)
# arrayを作るには、各軸の長さを指定する
arr.dim<-c(2,3,4)
a<-array(floor(runif(prod(arr.dim))*10),arr.dim)
v
m
a

> v
 [1] 8 3 5 6 4 7 3 6 2 8
> m
     [,1] [,2] [,3] [,4]
[1,]    5    6    1    7
[2,]    9    2    2    0
[3,]    1    3    6    3
> a
, , 1

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

, , 2

     [,1] [,2] [,3]
[1,]    8    5    2
[2,]    8    1    5

, , 3

     [,1] [,2] [,3]
[1,]    9    7    7
[2,]    2    0    8

, , 4

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

• ベクトルの内積と外積

# 内積
v%*%v
# 外積
v%o%v

> v%*%v
     [,1]
[1,]  312
> v%o%v
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]   64   24   40   48   32   56   24   48   16    64
 [2,]   24    9   15   18   12   21    9   18    6    24
 [3,]   40   15   25   30   20   35   15   30   10    40
 [4,]   48   18   30   36   24   42   18   36   12    48
 [5,]   32   12   20   24   16   28   12   24    8    32
 [6,]   56   21   35   42   28   49   21   42   14    56
 [7,]   24    9   15   18   12   21    9   18    6    24
 [8,]   48   18   30   36   24   42   18   36   12    48
 [9,]   16    6   10   12    8   14    6   12    4    16
[10,]   64   24   40   48   32   56   24   48   16    64

• 行列の処理

# 転置
> m
     [,1] [,2] [,3] [,4]
[1,]    5    6    1    7
[2,]    9    2    2    0
[3,]    1    3    6    3
> t(m)
     [,1] [,2] [,3]
[1,]    5    9    1
[2,]    6    2    3
[3,]    1    2    6
[4,]    7    0    3

• 正方行列、上三角行列、下三角行列、対角成分、対角行列

n.row<-3
n.col<-3
m2<-matrix(floor(runif(n.row*n.col)*10),n.row,n.col)
m3<-m2
m3[!upper.tri(m3)]<-0
m3
m3<-m2
m3[!lower.tri(m3)]<-0
m3
m3<-m2
diag(m3)
diag(diag(m3))

• 逆行列、固有値分解

m4<-solve(m2)
m4%*%m2
m2%*%m4
e.out<-eigen(m2)
e.out
v<-e.out[[2]]
s<-diag(e.out[[1]])
v%*%s%*%solve(v)
m2

> n.row<-3
> n.col<-3
> m2<-matrix(floor(runif(n.row*n.col)*10),n.row,n.col)
> m3<-m2
> m3[!upper.tri(m3)]<-0
> m3
     [,1] [,2] [,3]
[1,]    0    6    2
[2,]    0    0    2
[3,]    0    0    0
> m3<-m2
> m3[!lower.tri(m3)]<-0
> m3
     [,1] [,2] [,3]
[1,]    0    0    0
[2,]    6    0    0
[3,]    1    5    0
> m3<-m2
> diag(m3)
[1] 5 9 6
> diag(diag(m3))
     [,1] [,2] [,3]
[1,]    5    0    0
[2,]    0    9    0
[3,]    0    0    6
> 
> m4<-solve(m2)
> 
> m4%*%m2
             [,1]          [,2]          [,3]
[1,] 1.000000e+00 -1.540434e-15 -3.608225e-16
[2,] 4.024558e-16  1.000000e+00  1.942890e-16
[3,] 0.000000e+00 -4.440892e-16  1.000000e+00
> m2%*%m4
              [,1]          [,2]         [,3]
[1,]  1.000000e+00 -6.661338e-16 0.000000e+00
[2,] -2.220446e-16  1.000000e+00 5.551115e-17
[3,]  1.110223e-16  0.000000e+00 1.000000e+00
> 
> e.out<-eigen(m2)
> e.out
$values
[1] 14.7877109  4.3001970  0.9120921

$vectors
           [,1]        [,2]       [,3]
[1,] -0.5299298 -0.07540878 -0.6546323
[2,] -0.7097543 -0.30740064  0.5995079
[3,] -0.4641370  0.94858756 -0.4604854

> 
> v<-e.out[[2]]
> s<-diag(e.out[[1]])
> 
> v%*%s%*%t(v)
         [,1]     [,2]     [,3]
[1,] 4.568092 5.303678 3.604534
[2,] 5.303678 8.183490 3.365697
[3,] 3.604534 3.365697 7.248417
> m2
     [,1] [,2] [,3]
[1,]    5    6    2
[2,]    6    9    2
[3,]    1    5    6
> v%*%s%*%solve(v)
     [,1] [,2] [,3]
[1,]    5    6    2
[2,]    6    9    2
[3,]    1    5    6
> m2
     [,1] [,2] [,3]
[1,]    5    6    2
[2,]    6    9    2
[3,]    1    5    6