![[RjpWiki]](image/../Rlogo.png "[RjpWiki]")
行列に関する Tips 大全
行列に関する操作は、R をマスターする基本です。関連するTips を脈絡なくできるだけ集めたいと思います。お気づきの正統派・裏技テクニックをお寄せください。一部重複はむしろ好ましいと思います。
> m <- matrix(x, nrow=i, ncol=j, byrow=TRUE, dirnames=z) # 一般形
nrow, ncol の一方だけを与えると(可能な限り) x のサイズから、もう一方が暗黙のうちに決定される
x の長さが行列のサイズに足りないと、最初からりサイクル使用される
> mm <- matrix(1:12, nrow=3, ncol=4)
> mm <- matrix(1:12, nrow=3) # 自動的に ncol=4 とされる
> mm <- matrix(1:12, ncol=4) # 自動的に ncrow=3 とされる
> mm <- matrix(1:12, 3, 4)
> mm <- matrix(1:12, nc=4, nr=3, b=F) # オプション引数名は(一意的に判別可能な範囲で)省略可能
> mm <- array(1:12, c(3,4)) # 行列は length(dim)=2 の行列
> mm # 上は同じ行列 mm を作る
[,1] [,2] [,3] [,4]
[1,] 1 4 7 10 # 既定の byrow=FALSE では要素は列主導で並べられる
[2,] 2 5 8 11
[3,] 3 6 9 12
> mm <- matrix(1:12, 3, 4, byrow=TRUE) # 要素を行主導で行列化
> mm
[,1] [,2] [,3] [,4]
[1,] 1 2 3 4
[2,] 5 6 7 8
[3,] 9 10 11 12
注意：array は必ず列主導で要素が並べられる
> b <- matrix(1:12,nc=4, b=T) # 行主導
> b
[,1] [,2] [,3] [,4]
[1,] 1 2 3 4
[2,] 5 6 7 8
[3,] 9 10 11 12
> d <- t(array(1:12, c(4,3)) # array を使うならこうする c(3,4) でないことを注意
> d
[,1] [,2] [,3] [,4]
[1,] 1 2 3 4
[2,] 5 6 7 8
[3,] 9 10 11 12
> M <- matrix(1:6, , 3, byrow=TRUE) # ncol=3 (自動的に nrow=2 になる)
> M
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6
> M <- matrix(1:6,,3) # matrix(1:6,,3,byrow=FALSE) と同じ
> M
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
> M <- matrix(1:3, 3, 2)
[,1] [,2]
[1,] 1 1 # ベクトル 1:3 がリサイクル使用されている
[2,] 2 2
[3,] 3 3
次の方法が間違いがなくお勧め(rbind 関数で行ベクトルを縦に並べる)
> x <- rbind(c(0, 0.57, 0.57, 0.57, 0.57),
c(0.46, 0, 0, 0, 0 ),
c(0, 0.77, 0, 0, 0 ),
c(0, 0, 0.82, 0, 0 ),
c(0, 0, 0, 0.91, 0.65) )
> x
[,1] [,2] [,3] [,4] [,5]
[1,] 0.00 0.57 0.57 0.57 0.57
[2,] 0.46 0.00 0.00 0.00 0.00
[3,] 0.00 0.77 0.00 0.00 0.00
[4,] 0.00 0.00 0.82 0.00 0.00
[5,] 0.00 0.00 0.00 0.91 0.65
代わりに cbind 関数を使うと、列ベクトルとして横に並べられる
> x<-cbind(c(0,0.57,0.57,0.57,0.57), c(0.46,0,0,0,0), c(0,0.77,0,0,0),
c(0,0,0.82,0,0), c(0,0,0,0.91,0.65) )
> x
[,1] [,2] [,3] [,4] [,5]
[1,] 0.00 0.46 0.00 0.00 0.00
[2,] 0.57 0.00 0.77 0.00 0.00
[3,] 0.57 0.00 0.00 0.82 0.00
[4,] 0.57 0.00 0.00 0.00 0.91
[5,] 0.57 0.00 0.00 0.00 0.65
行列はベクトルとしても操作可能
> mm <- matrix(1:16, ncol=4) > mm[3] # ベクトルとしてアクセス [1] 3 > mm[3:5] [1] 3 4 5
逆にベクトルに次元属性を与えれば行列になる
> x <- 1:16 # ベクトル
> x
[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
> dim(x) <- c(4,4) # x に次元属性 dim(x)=c(4,4) を与える
> x # x は行列になった
[,1] [,2] [,3] [,4]
[1,] 1 5 9 13
[2,] 2 6 10 14
[3,] 3 7 11 15
[4,] 4 8 12 16
> dim(x) <- NULL # x の次元属性を取り去る
> x
[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
> a <- matrix(1:6, byrow=TRUE, nc=3)
> rownames(a) <- c("male", "female") # 行の名前
> colnames(a) <- c("left", "middle", "right") # 列の名前
> a
left middle right
male 1 2 3
female 4 5 6
> a["male","right"] # 名前をつけるとこうした分かり易い操作が出来る利点
[1] 3
rownames(a) <- NULL # 行列 x の行に付けられた名前をすべて落とす colnames(a) <- NULL # 列に付けられた名前をすべて落とす a <- unname(a) # とすれば，行名・列名を一遍に落とす dimnames(a) <- NULL # でも同じ
だが，元の行列を残してたまま，名前を取り除いた別の行列を作りたいときには unname の方が一度でできる。
x <- unname(a) # 1 ステップでできる y <- a # 2 ステップ掛かってしまう dimnames(y) <- NULL
> x <- matrix(1:16, 4, 4)
> x
[,1] [,2] [,3] [,4]
[1,] 1 5 9 13
[2,] 2 6 10 14
[3,] 3 7 11 15
[4,] 4 8 12 16
> dimnames(x) # 次元名無し
NULL
> names(dimnames(x)) # 次元名の名前も無し
NULL
# dimnames は次元名からなるリスト
> dimnames(x) <- list(letters[1:4], c("Y", "E", "A", "R"))
> dimnames(x)
[[1]]
[1] "a" "b" "c" "d"
[[2]]
[1] "Y" "E" "A" "R"
> colnames(x) # colnames(x) は dimnames(x)[[2]]
[1] "Y" "E" "A" "R"
> rownames(x) # rownames(x) は dimnames(x)[[1]]
[1] "a" "b" "c" "d"
> x
Y E A R
a 1 5 9 13
b 2 6 10 14
c 3 7 11 15
d 4 8 12 16
# 次元名に更に名前（成分名）をつける
> names(dimnames(x)) <- c("First", "Second")
> dimnames(x)
$First
[1] "a" "b" "c" "d"
$Second
[1] "Y" "E" "A" "R"
> x
Second
First Y E A R
a 1 5 9 13
b 2 6 10 14
c 3 7 11 15
d 4 8 12 16
> dimnames(x)$First # rownames(x) と同じ
[1] "a" "b" "c" "d"
> x1 <- 1:5 # 二つのベクトル
> x2 <- 2:6
> xname <- c("x1", "x2") # そのオブジェクト名の文字列ベクトル
> sapply(xname, get) # x1,x2 を行列にバインド(オブジェクト名が列ラベルになる)
> x1 x2
> [1,] 1 2
> [2,] 2 3
> [3,] 3 4
> [4,] 4 5
> [5,] 5 6
> A <- matrix(1:12, 3, 4) > dim(A) # 行列の大きさ（次元属性） [1] 3 4 > ncol(A) # 列数。 dim(A)[2] と同じ [1] 4 > nrow(A) # 行数。 dim(A)[1] と同じ [1] 3
> DF <- data.frame(XX= runif(20), YY=runif(20)) > colSums(DF)
colSums関数の中を倫理型ベクトルにする事で、Excelのcountif関数のようにも使える
> colSums(DF > 0.5) XX YY 11 9 > colSums(DF > -Inf) XX YY 20 20
> DF <- data.frame(XX= runif(20), YY=runif(20)) > colMeans(DF) > rowMeans(DF)
> d <- matrix(1:6, nr=2)
> d
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
転置行列を作る関数は t
> e <- t(d)
> e
[,1] [,2]
[1,] 1 2
[2,] 3 4
[3,] 5 6
> diag(1, ncol=3, nrow=3) # 値 1 がリサイクル使用される
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 1
> diag(1, 3) # 単位行列ならこのように省略もできる > diag(3) # 単位行列ならこれが最も簡単。上の例は正方対角行列の最も簡単な方法。 > diag(5, 3) # 対角成分が5の3×3行列（2004/08/14）
> diag(rep(1, 3)) # diag(c(1,1,1)) と同じ
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 1
> x <- matrix(0, ncol=3, nrow=3) # 全成分が 0 の行列
> diag(x) <- 1 # 対角成分を 1 にする
> x
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 1
> diag(x) <- c(1, 2, 3) # 対角成分を 1, 2, 3 にする
> x # diag(1:3) で直接作れる
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 2 0
[3,] 0 0 3
> x
[,1] [,2] [,3] [,4]
[1,] 1 5 9 13
[2,] 2 6 10 14
[3,] 3 7 11 15
[4,] 4 8 12 16
> diag(x)
[1] 1 6 11 16
> x = matrix(1:12, 3,4)
> x
[,1] [,2] [,3] [,4]
[1,] 1 4 7 10
[2,] 2 5 8 11
[3,] 3 6 9 12
> diag(x)
[1] 1 5 9
> diag(diag(x), 3,4)
[,1] [,2] [,3] [,4]
[1,] 1 0 0 0
[2,] 0 5 0 0
[3,] 0 0 9 0
上の応用でもいいが，もっと簡単になる
> i <- diag(4) # 引数が一つだけなら，それは単位行列の大きさを表すと解釈される
> i
[,1] [,2] [,3] [,4]
[1,] 1 0 0 0
[2,] 0 1 0 0
[3,] 0 0 1 0
[4,] 0 0 0 1
> x <- matrix(0, nr=3, nc=3)
> x
[,1] [,2] [,3]
[1,] 0 0 0
[2,] 0 0 0
[3,] 0 0 0
> x <- diag(0,3) # 正方行列ならこれが簡単（2004/08/14）
> x
[,1] [,2] [,3]
[1,] 0 0 0
[2,] 0 0 0
[3,] 0 0 0
> y <- diag(0,3,4) # 正方行列以外ならもう一個引数を追加（2004/08/14）
> y
[,1] [,2] [,3] [,4]
[1,] 0 0 0 0
[2,] 0 0 0 0
[3,] 0 0 0 0
例とする正方行列を作る
> b <- matrix(1:9, 3)
> b
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
下三角行列（対角を含む）を作る（上三角行列を0にする）
> b[upper.tri(b)] <- 0
> b
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 2 5 0
[3,] 3 6 9
対角を含めない場合は以下のように
> b[upper.tri(b,diag=TRUE)] <- 0
> b
[,1] [,2] [,3]
[1,] 0 0 0
[2,] 2 0 0
[3,] 3 6 0
上三角行列を作るときには lower.tri を使う
ファイルから読み込んで作成する場合，たとえば次の tab 区切りで0値のないテキストでは，"testdata.txt"
A0 A B C A1 A2 16.81 A3 78.80 14.4
以下のようにread.tableのオプションfillを使用する．
> x <- read.table("testdata.txt", header=TRUE, fill=TRUE)
A0 A B C
1 A1 NA NA NA
2 A2 16.81 NA NA
3 A3 78.80 14.4 NA
あとはis.na(x)，as.matrix, row.namesなどで整形する．
> x <- matrix(c(1,2,3,4, 0,5,6,7, 0,0,8,9, 0,0,0,10),4)
> x
[,1] [,2] [,3] [,4]
[1,] 1 0 0 0
[2,] 2 5 0 0
[3,] 3 6 8 0
[4,] 4 7 9 10
> y <- x+t(x)
> diag(y) <- diag(y)/2
> y
[,1] [,2] [,3] [,4]
[1,] 1 2 3 4
[2,] 2 5 6 7
[3,] 3 6 8 9
[4,] 4 7 9 10
x ＞＞ 三角行列の要素を含むベクトル（対角成分を含む）
n ＞＞ 正方行列のサイズ
byrow ＞＞ x が，行優先で使用されるとき TRUE（例を参照）
lower ＞＞ x が下三角行列のとき TRUE（例を参照）
関数の本体3行目は，if 文で使用する関数名を選択している
> tri.mat <- function(x, n, byrow=TRUE, lower=TRUE)
+ {
+ stopifnot(length(x) == n*(n+1)/2)
+ mat <- diag(n)
+ mat[(if (xor(byrow, lower)) lower.tri else upper.tri)(mat, diag=TRUE)] <- x
+ mat <- t(mat)+mat
+ diag(mat) <- diag(mat)/2
+ mat
+ }
> tri.mat(1:15, 5, byrow=TRUE, lower=TRUE)
[,1] [,2] [,3] [,4] [,5]
[1,] 1 2 4 7 11
[2,] 2 3 5 8 12
[3,] 4 5 6 9 13
[4,] 7 8 9 10 14
[5,] 11 12 13 14 15
> tri.mat(1:15, 5, byrow=TRUE, lower=FALSE)
[,1] [,2] [,3] [,4] [,5]
[1,] 1 2 3 4 5
[2,] 2 6 7 8 9
[3,] 3 7 10 11 12
[4,] 4 8 11 13 14
[5,] 5 9 12 14 15
> tri.mat(1:15, 5, byrow=FALSE, lower=TRUE)
[,1] [,2] [,3] [,4] [,5]
[1,] 1 2 3 4 5
[2,] 2 6 7 8 9
[3,] 3 7 10 11 12
[4,] 4 8 11 13 14
[5,] 5 9 12 14 15
> tri.mat(1:15, 5, byrow= FALSE, lower=FALSE)
[,1] [,2] [,3] [,4] [,5]
[1,] 1 2 4 7 11
[2,] 2 3 5 8 12
[3,] 4 5 6 9 13
[4,] 7 8 9 10 14
[5,] 11 12 13 14 15
> x <- matrix(1:6, nrow=2)
> row(x)
[,1] [,2] [,3]
[1,] 1 1 1
[2,] 2 2 2
> x <- matrix(1:6, nrow=2)
> col(x)
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 1 2 3
> x <- matrix(runif(9), ncol=3, nrow=3)
> x
[,1] [,2] [,3]
[1,] 0.8296189 0.9144157 0.60020004
[2,] 0.9440765 0.5413652 0.03560994
[3,] 0.0623679 0.6051179 0.57717385
> x[1, 3] # (1,3)成分
[1] 0.6002
> x[, 3] # 第3列 (ベクトルになる)
[1] 0.60020004 0.03560994 0.57717385
> x[1,] # 第1列 (ベクトルになる)
[1] 0.8296189 0.9144157 0.6002000
> x[, 3, drop=FALSE] #第3列 (3x1 行列になる)
[,1]
[1,] 0.60020004
[2,] 0.03560994
[3,] 0.57717385
> x[1, , drop=FALSE] # 第1行 (1x3行列になる)
[,1] [,2] [,3]
[1,] 0.8296189 0.9144157 0.6002
> x[c(1, 2), c(3, 2)] # 2x2副行列の取り出し、c(3,2)の順序に注意
[,1] [,2]
[1,] 0.60020004 0.9144157
[2,] 0.03560994 0.5413652
> x[-1,] # 第2,3行からなる2x3行列
[,1] [,2] [,3]
[1,] 0.9440765 0.5413652 0.03560994
[2,] 0.0623679 0.6051179 0.57717385
> x[x >= 0.5] # 0.5より大きな成分は？
[1] 0.8296189 0.9440765 0.9144157 0.5413652 0.6051179 0.6002000 0.5771739
> x[x >= 0.5] <- 0 # 0.5より大きな成分を0に置き換え
> x
[,1] [,2] [,3]
[1,] 0.0000000 0 0.00000000
[2,] 0.0000000 0 0.03560994
[3,] 0.0623679 0 0.00000000
> i <- which(x != 0) # 0でない成分の添字
> i # xをベクトルと考えた時の添字
[1] 3 8
> x[i] # x[x != 0] と同じ
[1] 0.06236790 0.03560994
> which(x != 0, arr.ind=TRUE) # 0でない添字を行列として取り出す
row col
[1,] 3 1 # つまり x[3,1]=0.06236790
[2,] 2 3
注意！ ベクトルや行列にNAが入っている場合は、which を使わない添字取り出し(NA が入る！)と、which を使った添字取り出しで結果が異なるので注意が必要です。
入力
x <- c(1, 2, NA) x[x == 1] x[which(x == 1)]
結果
> x<-c(1, 2, NA) > x[x == 1] [1] 1 NA > x[which(x == 1)] [1] 1
8バイト実数から4バイト整数に変換(必要メモリが半分になる
storage.mode(Array) <- "integer" Array[] <- as.integer(Array)
両者とも dim, dimnames 等の属性を保 (注：Array 要素は実際は整数型に強制変換されている)
> x <- matrix(runif(5*5), ncol=5, nrow=5) > y <- x[cbind(c(1, 2, 3), c(5, 4, 3))] # y はベクトル c(x[1, 5], x[2, 4], x[3, 3]) になる
注：cbind(1:3, 5:3) は実際は次の行列をあらわす
[,1] [,2] [1,] 1 5 [2,] 2 4 [3,] 3 3
特に x[cbind(1:5), cbind(1:5)] は x の対角成分からなるベクトル diag(x) と一致する
x[cbind(1:5, 1:5)] も同じ
> x <- matrix(1:16, 4, 4)
> x[(x == 1) | (x == 2)] <- 100
> x
[,1] [,2] [,3] [,4]
[1,] 100 5 9 13
[2,] 100 6 10 14
[3,] 3 7 11 15
[4,] 4 8 12 16
> x <- matrix(1:16, 4, 4)
> x[x %in% c(1, 3, 5, 7)] <- 100
> x
[,1] [,2] [,3] [,4]
[1,] 100 100 9 13
[2,] 2 6 10 14
[3,] 100 100 11 15
[4,] 4 8 12 16
> a <- matrix(1:25, nrow=5, byrow=TRUE)
> a
[,1] [,2] [,3] [,4] [,5]
[1,] 1 2 3 4 5
[2,] 6 7 8 9 10
[3,] 11 12 13 14 15
[4,] 16 17 18 19 20
[5,] 21 22 23 24 25
上のような行列から，1，2，4 行と２，４列を抽出した行列を作る
> a[c(1, 2, 4), c(2, 4)]
[,1] [,2]
[1,] 2 4
[2,] 7 9
[3,] 17 19
> a[-c(3, 5), -c(1, 3, 5)] # 上と同じ結果を得る（ "-" は除外を意味する）
> x <- matrix(1:100, 10, 10)
> x
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1 11 21 31 41 51 61 71 81 91
[2,] 2 12 22 32 42 52 62 72 82 92
[3,] 3 13 23 33 43 53 63 73 83 93
[4,] 4 14 24 34 44 54 64 74 84 94
[5,] 5 15 25 35 45 55 65 75 85 95
[6,] 6 16 26 36 46 56 66 76 86 96
[7,] 7 17 27 37 47 57 67 77 87 97
[8,] 8 18 28 38 48 58 68 78 88 98
[9,] 9 19 29 39 49 59 69 79 89 99
[10,] 10 20 30 40 50 60 70 80 90 100
> head(x) # 既定では先頭６行分を取り出す
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1 11 21 31 41 51 61 71 81 91
[2,] 2 12 22 32 42 52 62 72 82 92
[3,] 3 13 23 33 43 53 63 73 83 93
[4,] 4 14 24 34 44 54 64 74 84 94
[5,] 5 15 25 35 45 55 65 75 85 95
[6,] 6 16 26 36 46 56 66 76 86 96
> head(x, n=4) # 先頭4行分を取り出す
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1 11 21 31 41 51 61 71 81 91
[2,] 2 12 22 32 42 52 62 72 82 92
[3,] 3 13 23 33 43 53 63 73 83 93
[4,] 4 14 24 34 44 54 64 74 84 94
> tail(x) # 既定では末尾６行分を取り出す
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 5 15 25 35 45 55 65 75 85 95
[2,] 6 16 26 36 46 56 66 76 86 96
[3,] 7 17 27 37 47 57 67 77 87 97
[4,] 8 18 28 38 48 58 68 78 88 98
[5,] 9 19 29 39 49 59 69 79 89 99
[6,] 10 20 30 40 50 60 70 80 90 100
> tail(x, n=3) # 末尾３行分を取り出す
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 8 18 28 38 48 58 68 78 88 98
[2,] 9 19 29 39 49 59 69 79 89 99
[3,] 10 20 30 40 50 60 70 80 90 100
NA を含まない，完全な行列を作るということ
> x <- c(1, 3, 2, 1, NA, 4)
> y <- c(2, NA, 3, 2, 5, 4)
> z <- c(5, 3, 2, 3, 2, 1)
> mat <- cbind(x, y, z)
> mat
x y z
[1,] 1 2 5
[2,] 3 NA 3
[3,] 2 3 2
[4,] 1 2 3
[5,] NA 5 2
[6,] 4 4 1
> matc <- subset(mat, complete.cases(mat))
> matc
x y z
[1,] 1 2 5
[2,] 2 3 2
[3,] 1 2 3
[4,] 4 4 1
> matc2 <- mat[complete.cases(mat),]
> matc2
x y z
[1,] 1 2 5
[2,] 2 3 2
[3,] 1 2 3
[4,] 4 4 1
> na.omit(mat) # NA を除くなら na.omit も使える
x y z
[1,] 1 2 5
[2,] 2 3 2
[3,] 1 2 3
[4,] 4 4 1
attr(,"na.action")
[1] 5 2
attr(,"class")
[1] "omit"
> x <- matrix(1:64, ncol=8)
> x
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,] 1 9 17 25 33 41 49 57
[2,] 2 10 18 26 34 42 50 58
[3,] 3 11 19 27 35 43 51 59
[4,] 4 12 20 28 36 44 52 60
[5,] 5 13 21 29 37 45 53 61
[6,] 6 14 22 30 38 46 54 62
[7,] 7 15 23 31 39 47 55 63
[8,] 8 16 24 32 40 48 56 64
> x2 <- x
> dim(x2) <- c(2, 4, 2, 4) # 2x2 副行列を系統的に取り出すために一旦 4 次元配列化
> x2[, 1, , 1]
[,1] [,2]
[1,] 1 9
[2,] 2 10
> x2[, 2, , 1]
[,1] [,2]
[1,] 3 11
[2,] 4 12
> x2[, 1, , 2]
[,1] [,2]
[1,] 17 25
[2,] 18 26
> x4 <- x
> dim(x4) <- c(4, 2, 4, 2) # 4x4 副行列を系統的に取り出すために一旦 4 次元配列化
> x4[, 1, , 1]
[,1] [,2] [,3] [,4]
[1,] 1 9 17 25
[2,] 2 10 18 26
[3,] 3 11 19 27
[4,] 4 12 20 28
> invisible(apply(x4, c(2, 4), print)) # 一度に出力する
# apply 関数の既定の(この場合余計な)出力をしないため invisible 関数を使ていることに注意
[,1] [,2] [,3] [,4]
[1,] 1 9 17 25
[2,] 2 10 18 26
[3,] 3 11 19 27
[4,] 4 12 20 28
[,1] [,2] [,3] [,4]
[1,] 5 13 21 29
[2,] 6 14 22 30
[3,] 7 15 23 31
[4,] 8 16 24 32
[,1] [,2] [,3] [,4]
[1,] 33 41 49 57
[2,] 34 42 50 58
[3,] 35 43 51 59
[4,] 36 44 52 60
[,1] [,2] [,3] [,4]
[1,] 37 45 53 61
[2,] 38 46 54 62
[3,] 39 47 55 63
[4,] 40 48 56 64
> mat <- rbind(c(1, 2, 3), c(4, 5, 6)) # 行ベクトルをあたえて行列を作る
> mat
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6
> cbind(c(1, 2, 3), c(4, 5, 6)) # 列ベクトルをあたえて行列を作る
[,1] [,2]
[1,] 1 4
[2,] 2 5
[3,] 3 6
> mat <- rbind(c(1, 2, 3), c(4, 5, 6)) # 2x3 行列
> mat
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6
縦ベクトル t(c(1, 2, 3, 4, 5, 6)) にするには
> vec <- as.vector(t(mat)) > vec [1] 1 2 3 4 5 6
とする。行列は列主導で配置されるから転置演算 t(mat) が必要。転置しないと
> vec <- as.vector(mat) > vec [1] 1 4 2 5 3 6
もしベクトルでなく 1x6 行列が必要なら次のようにする： (行列は dim 属性 c(nrow,ncol) を持つベクトルに他ならないことを注意)：
> matrix(t(mat), ncol=1)
[,1]
[1,] 1
[2,] 2
[3,] 3
[4,] 4
[5,] 5
[6,] 6
行列 X に対して関数
> Idx <- function(i, j, X) dim(X)[1]*(j-1)+i
を定義すると、行列成分 X[i, j] はベクトル成分 X[Idx(i, j, X)] としてアクセスできる。
> X <- matrix(runif(200), 20, 10) > X[5,7]==X[Idx(5, 7, X)] [1] TRUE > X[15,8]==X[Idx(15, 8, X)] [1] TRUE
> x <- matrix(1:12, ncol=3)
> x
[,1] [,2] [,3]
[1,] 1 5 9
[2,] 2 6 10
[3,] 3 7 11
[4,] 4 8 12
> rowSums(x)
[1] 15 18 21 24
> colSums(x)
[1] 10 26 42
この他に
> drop(x %*% rep(1, ncol(x))) # 行和 (drop は 1x4 行列をベクトル化する) [1] 15 18 21 24 > drop(t(x) %*% rep(1,nrow(x))) # 列和 [1] 10 26 42
もう一つのやり方
> apply(x, 2, sum) # 列和 [1] 10 26 42 > apply(x, 1, sum) # 行和 [1] 15 18 21 24
があります。
> x <- 1:12 > dim(x) <- c(4, 3) > y <- x[c(2, 4),] > dim(x[c(2, 4),]) [1] 2 3 > dim(x) <- c(12, 1) # x は 12x1 行列 > dim(x[c(2, 4),]) NULL # 結果はもはや行列でない！ > dim(x[c(2, 4),,drop=FALSE]) # オプション drop=FALSE [1] 2 1 # 結果は 2x1 の行列！
行列を画像と考え image 関数で表示する。option の col で色調を変えることが出来る (col の例: rainbow,heat.colors,topo.colors,terrain.colors, gray)
> x <- y <- 1:10 > z <- matrix(rnorm(100), ncol=10) > image(x, y, z, col=topo.colors(12)) # 地形図色調を12段階で使う
> x <- matrix(c(24,22,23,21, 13,12,11,14, 33,32,31,34), nr=4)
> x
[,1] [,2] [,3]
[1,] 24 13 33
[2,] 22 12 32
[3,] 23 11 31
[4,] 21 14 34
> order(x[, 2]) # 2列目の数値の順序
[1] 3 2 1 4
> x[order(x[, 2]),] # 2列目の数値により並べ替え
[,1] [,2] [,3]
[1,] 23 11 31
[2,] 22 12 32
[3,] 24 13 33
[4,] 21 14 34
> order(x[3,]) # 3行目の数値の順序
[1] 2 1 3
> x[, order(x[3,])] # 3 行目の数値により並べ替え
[,1] [,2] [,3]
[1,] 13 24 33
[2,] 12 22 32
[3,] 11 23 31
[4,] 14 21 34
max.at <- function (x) # 関数にするほどでもないが
{
o <- order(x)
cbind(row(x)[o], col(x)[o], x[o])
}
> x
[,1] [,2] [,3] [,4]
[1,] 2 5 11 8
[2,] 1 4 10 7
[3,] 3 6 12 9
> max.at(x)
[,1] [,2] [,3]
[1,] 2 1 1
[2,] 1 1 2
[3,] 3 1 3
[4,] 2 2 4
[5,] 1 2 5
[6,] 3 2 6
[7,] 2 4 7
[8,] 1 4 8
[9,] 3 4 9
[10,] 2 3 10
[11,] 1 3 11
[12,] 3 3 12
> x <- matrix(1:12, 3, 4)
> x
[,1] [,2] [,3] [,4]
[1,] 1 4 7 10
[2,] 2 5 8 11
[3,] 3 6 9 12
> apply(x,2,rev) #上下反転
[,1] [,2] [,3] [,4]
[1,] 3 6 9 12
[2,] 2 5 8 11
[3,] 1 4 7 10
> apply(x,1,rev) #反時計回り90度回転
[,1] [,2] [,3]
[1,] 10 11 12
[2,] 7 8 9
[3,] 4 5 6
[4,] 1 2 3
> t(apply(x,1,rev)) #左右反転
[,1] [,2] [,3] [,4]
[1,] 10 7 4 1
[2,] 11 8 5 2
[3,] 12 9 6 3
> t(apply(x,2,rev)) #時計回り90度回転
[,1] [,2] [,3]
[1,] 3 2 1
[2,] 6 5 4
[3,] 9 8 7
[4,] 12 11 10
別解
> x[nrow(x):1,] #上下反転
[,1] [,2] [,3] [,4]
[1,] 3 6 9 12
[2,] 2 5 8 11
[3,] 1 4 7 10
> t(x[,ncol(x):1]) #反時計回り90度回転
[,1] [,2] [,3]
[1,] 10 11 12
[2,] 7 8 9
[3,] 4 5 6
[4,] 1 2 3
> x[,ncol(x):1] #左右反転
[,1] [,2] [,3] [,4]
[1,] 10 7 4 1
[2,] 11 8 5 2
[3,] 12 9 6 3
> t(x[nrow(x):1,]) #時計回り90度回転
[,1] [,2] [,3]
[1,] 3 2 1
[2,] 6 5 4
[3,] 9 8 7
[4,] 12 11 10
行列のサイズがほどほどの場合には後者が速いが，速いと行ってもオーダーレベルの話で，実際の実行時間の差はない。行列サイズが大きいとほとんど差がなくなる。ほとんど，好みの問題と言っていいレベルか...
> x <- matrix(rnorm(1000000), 1000)
> system.time(a <- apply(x,2,rev)) # 上下反転
ユーザ システム 経過
0.059 0.006 0.083
> system.time(b <- x[nrow(x):1,])
ユーザ システム 経過
0.015 0.000 0.014
> all(a==b)
[1] TRUE
> system.time(a <- apply(x,1,rev)) # 反時計回り90度回転
ユーザ システム 経過
0.044 0.004 0.062
> system.time(b <- t(x[,ncol(x):1]))
ユーザ システム 経過
0.023 0.001 0.023
> all(a==b)
[1] TRUE
> system.time(a <- t(apply(x,1,rev))) # 左右反転
ユーザ システム 経過
0.052 0.003 0.072
> system.time(b <- x[,ncol(x):1]) # 左右反転
ユーザ システム 経過
0.016 0.001 0.016
> all(a==b)
[1] TRUE
> system.time(a <- t(apply(x,2,rev))) # 時計回り90度回転
ユーザ システム 経過
0.046 0.003 0.057
> system.time(b <- t(x[nrow(x):1,]))
ユーザ システム 経過
0.022 0.000 0.022
> all(a==b)
[1] TRUE
ある条件を満たす行列の特定要素の行列添字の全部を取り出すには which関数をオプション arr.ind=TRUE で使用する。既定の arr.ind=FALSEでは行列をベクトル化した上で該当要素の添字ベクトルを返す
> x <- matrix(1:12, 3, 4)
> x
[,1] [,2] [,3] [,4]
[1,] 1 4 7 10
[2,] 2 5 8 11
[3,] 3 6 9 12
> y <- which(x%%3==0, arr.ind=TRUE) # 3 の倍数である要素の行列添字を取り出す
> y
row col
[1,] 3 1
[2,] 3 2
[3,] 3 3
[4,] 3 4
> x[y[1, 1], y[1, 2]]
[1] 3
> x[y[4, 1], y[4, 2]]
[1] 12
単に該当要素を全部取り出すならベクトル化した方がわかりやすい。R の行列・配列は次元属性を持つベクトルであり、ベクトルとしてアクセスできる。
> x[x%%3 == 0] # 3の倍数である x の全要素 [1] 3 6 9 12 # 行列 x の偶数要素を全て 0 に置き換え(行列をベクトルとして添字操作) > x[x%%2 == 0] <- 0 > x [,1] [,2] [,3] [,4] [1,] 1 0 7 0 [2,] 0 5 0 11 [3,] 3 0 9 0
> m <- matrix(c(1,2,3, 5,6,2, 7,5,9), nr=3)
> m
[,1] [,2] [,3]
[1,] 1 5 7
[2,] 2 6 5
[3,] 3 2 9
> det(m)
[1] -69
> kappa()
正方行列の逆行列を solve() 関数で求めることができる。しかしながら逆行列を求めることは数値精度の点で問題を生じることが多く、また時間を食うという意味で"本当に必要なときだけ"使うことが数値計算の常識(実際逆行列そのものが必要になることは殆んど無いはず)
> A <- array(runif(9), c(3, 3))
> B <- solve(A)
> A%*%B # 検査(数値精度の範囲内で単位行列
[,1] [,2] [,3]
[1,] 1.000000e+00 -1.103176e-17 -2.745742e-17
[2,] 1.083863e-16 1.000000e+00 -1.039106e-16
[3,] 3.718813e-17 1.237075e-16 1.000000e-00
> B%*%A # 検査(数値精度の範囲内で単位行列
[,1] [,2] [,3]
[1,] 1.000000e+00 5.573382e-17 -2.075354e-16
[2,] -2.778268e-17 1.000000e+00 2.326698e-16
[3,] -7.285839e-17 -5.995638e-17 1.000000e+00
> x <- array(runif(8), c(2, 4)) > sum(x^2) # x^2 は行列 x*x [1] 3.185785 > sum(diag(t(x)%*%x)) # 同じことを複雑に(線形代数に詳しい人なら納得) [1] 3.185785 > sum(diag(x%*%t(x))) [1] 3.185785
2×3行列を作る
> a <- matrix(1:6, byrow=TRUE, nc=3)
> a
[,1] [,2] [,3]
[1,] 1 2 3
[2,] 4 5 6
行列の各要素のスカラー倍
> a*2
[,1] [,2] [,3]
[1,] 2 4 6
[2,] 8 10 12
行列の各要素に (1,2,3) という 3 つの要素を持つベクトルを「列優先で」掛ける。不足する場合には繰り返し適用する（リサイクルされるという）。
> a*1:3 # : という演算子は一番強い。a*(1:3) と書くと誤解がないが。
[,1] [,2] [,3]
[1,] 1 6 6
[2,] 8 5 18
行列の各要素に (1, 2, 3) という 3 つの要素を持つベクトルを「行優先で」掛けるには，以下のようにする。つまり，転置して掛けてから転置して元に戻す。
> t(t(a)*1:3)
[,1] [,2] [,3]
[1,] 1 4 9
[2,] 4 10 18
> b <- matrix(11:16, byrow=TRUE, nc=3)
> a+b # 要素どうしの足し算
[,1] [,2] [,3]
[1,] 12 14 16
[2,] 18 20 22
> 1/a # このような記述は，要するに要素ごとの演算（これは，要素の逆数）
[,1] [,2] [,3]
[1,] 1.00 0.5 0.3333333
[2,] 0.25 0.2 0.1666667
> sum(1/a) # 要素の逆数の総和
[1] 2.45
qr(X)$rank
singular valueのどの値を0とみなすかは toleranceで決定する。 svd() を利用してsingular valueの数を数えても求まる。
ベクトル化されている関数はまた行列化されている
> A <- array(1:9, c(3, 3))
> A
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
> log(A)
[,1] [,2] [,3]
[1,] 0.000000 1.38629 1.94591
[2,] 0.693147 1.60944 2.07944
[3,] 1.098612 1.79176 2.19722
> exp(A)
[,1] [,2] [,3]
[1,] 2.71828 54.5982 1096.63
[2,] 7.38906 148.4132 2980.96
[3,] 20.08554 403.4288 8103.08
> library(Matrix)
関数 qr は行列の QR 分解を計算する。これは LAPACK のルーチンDQRDC, DGEQP3 そして ZGEQP3 (複素行列の場合)へのインタフェイスである。
一般書式
qr(x, tol = 1e-07 , LAPACK = FALSE)
qr.coef(qr, y)
qr.qy(qr, y)
qr.qty(qr, y)
qr.resid(qr, y)
qr.fitted(qr, y, k = qr$rank)
qr.solve(a, b, tol = 1e-7)
solve(a, b, ...)
is.qr(x)
as.qr(x)
引数
x: QR 分解を計算すべき行列
tol: x の列間の線形従属性を検出するための数値的許容度。
LINPACK が真の時だけ使われる。
qr: qr() により計算された形式の QR 分解
y, b: 線形方程式の右辺のベクトル、もしくは行列
a: QR 分解、もしくは (`qr.solve' だけ)三角行列
k: 有効ランク
LAPACK: 論理値。実の x に対しては LAPACK を使うか LINPACK を使うか指示
...: その他の引数で、他のメソッドへ、又はから、引き渡される引数
返り値:
LINPACK もしくは LAPACK が計算した行列の QR 分解
返り値中の成分は DQRDC/DGEQP3/ZGEQP3 が返す値に直接対応
qr: x と同じ次元の行列。
上三角は分解の R 部分を含む
下三角部分は分解の Q 部分に関する情報を圧縮して含む
圧縮法は DQRDC と DGEQP3 では異なることを注意
qraux: 長さ ncol(x) のベクトルで、Q に関する補助情報を含む。
rank: 分解により計算された x のランク。LAPACK ケースでは常にフルランク
pivot: 分解の過程で使われたピボット選択に関する情報
LAPACK が計算した非複素 QR オブジェクトは値が TRUE の属性 useLAPACK を持つ
詳細：QR 分解は多くの統計理論で重要な役目を果たす。特に、与えられた行列 A、ベクトル b に対し方程式 Ax=b を解く。回帰係数の計算や、ニュートン・ラフソンアルゴリズムで有用である。
関数 qr.coef, qr.resid そして qr.fitted は、QR 分解された行列に y を当てはめた際の、係数、残差、そして当てはめ値を返す。qr.qy とqr.qty は QR 行列 Q に対する Q %*% y と t(Q) %*% y を計算する。これらのすべての関数は、もしそれらが存在すれば、 x と y の dimanmes (そして names)を保存する。
solve.qr は qr オブジェクトに対する solve のメッソドである。qr.solve は QR 分解を用いて線形方程式系を解く。もし a が QR 分解ならば、これは solve.qr と一致するが、もし a が長方形行列なら最初にその QR 分解が計算される。どちらも過剰、過小方程式系を扱え、最小長解、又は適当ならば最小自乗解を与える。
is.qr は、もし x が成分名 qr, rank, qraux を持つリストなら TRUE を返し、さもなければ FALSE を返す。
オブジェクトをモード qr に強制変換することはできない。オブジェクトは QR 分解であるか、無いかのいずれかである。
注意；行列の行列式を計算する(本当に必要かどうか要検討)には、QR 分解は固有値を計算する(eigen)より相当効率的である。関数 det を参照。複素ケースを含む LAPACK の使用は列ピボットを用い、行列のランク落ちを検出しようとしない。
行列の再構成には qr.Q, qr.R, qr.X を使う。
例:
hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, "+") }
h9 <- hilbert(9); h9
qr(h9)$rank # ランクはたった 7
qrh9 <- qr(h9, tol = 1e-10) # 数値的許容度を厳しくする
qrh9$rank # ランクは 9
## 線形方程式系 H %*% x = y を解く:
y <- 1:9/10
x <- qr.solve(h9, y, tol = 1e-10) # または同値であるが :
x <- qr.coef(qrh9, y) # 同値であるが次よりははるかに良い
#-- solve(h9) %*% y
h9 %*% x # y に等しい
svd(x, nu = min(n, p), nv = min(n, p), LINPACK = FALSE) # 書式
La.svd(x, nu = min(n, p), nv = min(n, p), method = c("dgesdd", "dgesvd"))
引数
x: 特異値分解を求めたい行列
nu: 計算したい左特異値の数。`method = "dgesdd"' でない限り
`0', `nrow(x)' そして `ncol(x)' のいずれ
nv: 計算したい右特異値の数。`0' か `ncol(x)' のいずれ
LINPACK: 論理値。LINPACK を使うか ( R < 1.7.0 との両立性のため)
method: 実数値の場合に使うべき LAPACK ルーティン
返り値。次の成分を持つリス
d: x の特異値を含むベクトル
u: X の左特異ベクトルからなる列を含む行列。nu > 0 の時与えられる
v: X の右特異ベクトルからなる列を含む行列。nv > 0 の時与えられる
La.svd は v の代わりにその転置である vt (複素行列の場合は共役) を返す
svd と La.svd は少し異なるインタフェイスを用いる。主な関数は LAPACK のルーチンDGESDD と ZGESVD である。`svd(LINPACK=TRUE)' は単に過去との一貫性のために存在しLINPACK ルーチン DSVDC へのインタフェイスである。
La.svd は LAPACK ルーチン DGESVD と DGESDD へのインタフェイス。後者は、特異値が必要な場合、普通かなり高速。効果は最適化 BLAS を用いる事で目に見えるようになる`method="dgesdd"' は IEEE 754 算術を必要とする。使用のプラットフォームがこれをサポートしていないと、警告とともに `method="dgesvd"' が使われる
LINPACK における行列 X の特異値分解とは
X = U D V' (V' は V の転置行列) U, V 直交行列 D 対角行列。 D[i, i] は i 番目の特異値。D = U' X V となる
例 (example(svd))
> hilbert <- function(n) { # ヒルベルト行列を作る関数
i <- 1:n
1/outer(i - 1, i, "+")
}
> str(X <- hilbert(9)[, 1:6])
num [1:9, 1:6] 1.000 0.500 0.333 0.250 0.200 ...
> str(s <- svd(X)) # 特異値I分解を実行し、内容を見る
List of 3
$ d: num [1:6] 1.67e+00 2.77e-01 2.22e-02 1.08e-03 3.24e-05 ...
$ u: num [1:9, 1:6] -0.724 -0.428 -0.312 -0.248 -0.206 ...
$ v: num [1:6, 1:6] -0.736 -0.443 -0.327 -0.263 -0.220 ...
> D <- diag(s$d) # 特異値ベクトルを対角成分とする対角行列
> s$u %*% D %*% t(s$v) # 数値誤差範囲内で X と一致するはず
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1.0000000 0.5000000 0.33333333 0.25000000 0.20000000 0.16666667
[2,] 0.5000000 0.3333333 0.25000000 0.20000000 0.16666667 0.14285714
[3,] 0.3333333 0.2500000 0.20000000 0.16666667 0.14285714 0.12500000
[4,] 0.2500000 0.2000000 0.16666667 0.14285714 0.12500000 0.11111111
[5,] 0.2000000 0.1666667 0.14285714 0.12500000 0.11111111 0.10000000
[6,] 0.1666667 0.1428571 0.12500000 0.11111111 0.10000000 0.09090909
[7,] 0.1428571 0.1250000 0.11111111 0.10000000 0.09090909 0.08333333
[8,] 0.1250000 0.1111111 0.10000000 0.09090909 0.08333333 0.07692308
[9,] 0.1111111 0.1000000 0.09090909 0.08333333 0.07692308 0.07142857
t(s$u) %*% X %*% s$v # 数値誤差の範囲内で対角行列 D になっている
[,1] [,2] [,3] [,4] [,5]
[1,] 1.668433e+00 4.229676e-16 8.466264e-17 -3.968180e-17 6.133874e-17
[2,] -5.922285e-17 2.773727e-01 -5.047300e-17 3.649696e-17 -3.833671e-17
[3,] -2.718785e-17 -4.455393e-17 2.223722e-02 5.586606e-18 -2.382407e-17
[4,] 4.463291e-17 4.649194e-17 -3.204432e-19 1.084693e-03 1.337568e-17
[5,] -6.266287e-18 1.087336e-18 -3.472332e-18 2.902897e-18 3.243788e-05
[6,] 2.591748e-17 -3.282344e-18 1.366355e-17 -3.892934e-18 2.680517e-18
[,6]
[1,] 7.870630e-17
[2,] 1.825525e-17
[3,] 4.753761e-18
[4,] 1.035196e-17
[5,] 6.299358e-18
[6,] 5.234864e-07
要約
> eigen(Sm) #固有値分解 > eigen(Sm)$val #固有値 > eigen(Sm)$vec # 固有ベクトル
固有値分解は正方行列 m を m = V D V' と分解 V' は V の転置行列 D は対角行列で、各成分は m の固有値。つまり m %*% V[, i] = D[i, i] %*% V[, i]
m は次のような対称行列とする
0.4015427 0.08903581 -0.2304132 0.08903581 1.60844812 0.9061157 -0.23041322 0.9061157 2.9692562
> sm <- eigen(m, sym=TRUE) # 関数 eigen を適用 > sm$values [1] 3.4311626 1.1970027 0.3510817 # 三つの固有値
>sm$vectors # それぞれに対応する固有ベクトルを並べた行列
[,1] [,2] [,3]
[1,] -0.05508142 -0.2204659 0.9738382
[2,] 0.44231784 -0.8797867 -0.1741557
[3,] 0.89516533 0.4211533 0.1459759
> V <- sm$vectors
> t(V) %*% V # (計算誤差内で) V は直交行列
[,1] [,2] [,3]
[1,] 1.000000e+00 -1.665335e-16 -5.551115e-17
[2,] -1.665335e-16 1.000000e+00 2.428613e-16
[3,] -5.551115e-17 2.428613e-16 1.000000e+00
> V %*% diag(sm$values) %*% t(V) # (計算誤差内で) m が復元されている
[,1] [,2] [,3]
[1,] 0.40154270 0.08903581 -0.2304132
[2,] 0.08903581 1.60844812 0.9061157
[3,] -0.23041320 0.90611570 2.9692562
R'R=x となる行列 R を計算する
一般書式
chol(x, pivot = FALSE, LINPACK = pivot) La.chol(x)
引数
x: (実正方対称)正定値符号行列 pivot: ピボット操作(x の行・列の置換)を行う LINPACK: 論理値。LINPACK を使うかどうか(R < 1.7.0 との互換性のため)
返り値
コレスキ分解の上三角因子、つまり R'R = x となる行列 R もしピボット操作が用いられると、二つの追加属性 pivot と rank が返される
説明:
chol(pivot = TRUE) は LINPACK ルーティン DCHDC へのインタフェイス La.chol は LAPACK ルーティン DPOTRF へのインタフェイス x の上三角行列部分だけが使われる(x が対称と仮定) pivot = FALSE で x が正定値符号ならばエラーとなる x が非負値定符号(ゼロの固有値を持つ)ならば、数値的許容度のせいで、やはりエラーが起きる もし pivot = TRUE ならば、非負値定符号行列 x のコレスキ分解を計算できる x のランクは attr(Q, "rank") で計算できるが、数値的エラーの可能性 ピボットは attr(Q, "pivot") で返される。必ずしも t(Q) %*% Q = x とはならない pivot <- attr(Q, "pivot") かつ oo <- order(pivot) とすれば、t(Q[, oo]) %*% Q[, oo] = x もしくは t(Q) %*% Q = x[pivot, pivot] となる
注意: x の対称性の検査は行わない。もし pivot = TRUE で x が非負値定符号でなければ、警告はでず、結果は無意味になる。もし x が非負値定符号ならば pivot = TRUE を使うこと
関連関数
chol2inv: (ピボット操作を行わない)逆行列を計算 backsolve: 上三角部分を持つ線形システムを解く
例 (example(chol))
m <- matrix(c(5, 1, 1, 3), 2, 2) cm <- chol(m) t(cm) %*% cm # m に等しい crossprod(cm) # m に等しい # 非負値定符号行列のケース x <- matrix(c(1:5, (1:5)^2), 5, 2) x <- cbind(x, x[, 1] + 3*x[, 2]) m <- crossprod(x) qr(m)$rank # 期待されるように 2 になる # 数値丸め誤差により chol() は失敗することもある # qr() と異なり chol() は数値的許容度を用いない try(chol(m)) Q <- chol(m, pivot = TRUE)) # 誤ったランクに注意(数値誤差による) ## これは次のように使う pivot <- attr(Q, "pivot") oo <- order(pivot) t(Q[, oo]) %*% Q[, oo] # m が復元された
> x %*% A %*% x
但し x はベクトル。もし x が length(x) 行、1 列の行列なら次のようにしないとエラー
> t(x) %*% A %*% x
機能：Matrix product，別表現：なし
A <- matrix(1:12, 4, 3) B <- matrix(runif(15), 3, 5) A %*% B # (4x3)%*%(3x5) = 4x5 行列になる
機能：Outer product，別表現：outer()
z <- outer(x, y, FUN="afun") は次元 c(dim(x), dim(y)) で、要素 z[i, j] が afun(x[i], y[j]) である行列になる。関数 afun は任意で、ベクトルにタグ名があれば保存される。x%o%y は outer(x, y, FUN="*") の省略形
> x <- runif(3)
> y <- 1:3
> outer(x, y, FUN="^")
[,1] [,2] [,3]
[1,] 0.56709640 0.321598331 0.1823772568 # [1,2] 成分は 0.56709640^2
[2,] 0.06480775 0.004200045 0.0002721955
[3,] 0.06183690 0.003823802 0.0002364520
機能：Kronecker product，別表現：kronecker()
> A <- matrix(10^(1:4), 2, 2)
> A
[,1] [,2]
[1,] 10 1000
[2,] 100 10000
> B <- matrix(runif(9), 3, 3)
> A%x%B # 小行列 A[i,j]*B を並べたサイズ (2*3)x(2*3) の大きな行列になる
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1.223105 9.172152 5.125222 122.3105 917.2152 512.5222
[2,] 6.986403 6.852122 1.428127 698.6403 685.2122 142.8127
[3,] 8.840513 3.958830 1.759327 884.0513 395.8830 175.9327
[4,] 12.231050 91.721521 51.252222 1223.1050 9172.1521 5125.2222
[5,] 69.864034 68.521222 14.281271 6986.4034 6852.1222 1428.1271
[6,] 88.405128 39.588303 17.593270 8840.5128 3958.8303 1759.3270
> a + b # いずれも行列 a,b のサイズが一致する必要 > a - c > a * b > a / b
Rの行列・配列は添字が1から始まるが、場合により0等の添字を使いたいこともある。アドオンパッケージ Oarrayは任意のオフセットの行列・配列を定義することを可能にする。主な使い方は添字に0を含める、添字が例えばA[i,j], 3<=i,j<=10, しか必要ないとき不要なメモリーを消費しない、等が考えられる。Oarray を使うにはあらかじめパッケージをインストールし、library(Oarray) で読み込んでおく必要がある
注意: R では負の添字は「除外」を意味するので、問題が起きる。真の負の添字を使いたければ drop.negative=FALSE にする。base パッケージ中の関数 as.array は as.array.Oarray() 関数で置き換えられる
一般形
Oarray(data=NA, dim=length(data), dimnames=NULL, offset=NA,
drop.negative=TRUE)
as.Oarray(x, offset=NA, drop.negative=TRUE)
as.array(x)
print(x)
x[i]
x[i, j, ...]
x[i, j, ..., drop=TRUE]
引数
data, dim, dimnames: array 関数と同じ意味 offset: 各添字組に対する最初の添字値のベクトル (既定値はすべて1) drop.negative: 論理値。負の添字が「除外」を意味するかどうか (既定値はTRUEで除外) x: 配列で、クラス Oarray でも良い i, j, ...: x への添字引数
返り値 典型的にはクラス Oarray を持つ、もしくは持たない配列。副配列を抽出すると単なる配列になる。Oarray オブジェクトを指定すると Oarray に再びなる。print 関数は Oarray オブジェクトを適切に表示する
> A <- Oarray(1:16, c(4, 4), offset=c(-2, -1))
> A
[,-1] [, 0] [, 1] [, 2]
[-2,] 1 5 9 13
[-1,] 2 6 10 14
[ 0,] 3 7 11 15
[ 1,] 4 8 12 16
> A[-2, 0] # 負の添字が除外を意味しないことを注意
[1] 5
> A[-2, -1]
[1] 1
> A[(-2):0, -1]
[1] 1 2 3
> A[(-2):0, (-1):1] # 副配列を取り出すと通常の配列になる (drop=FALSE をつけても同様 ?)
[,1] [,2] [,3]
[1,] 1 5 9
[2,] 2 6 10
[3,] 3 7 11
MASSパッケージ所収。以下の４つの条件（Penrose conditions）を満たす一般逆行列「ムーア・ペンローズ逆行列」を求める。（注：一般逆行列 generalized inverse of matrix)
Penrose conditions
(i) AGA = A (一般逆行列の性質) (ii) GAG = G (反射型の性質) (iii) (GA)' = GA (ノルム最小型の性質) (iv) (AG)' = AG (最小二乗型の性質)
一般逆行列にはこの他に、反射型一般逆行列（(i)(ii)を満たす）、ノルム最小型一般逆行列、最小二乗型一般逆行列がある。ginv()では特異値分解に基づき一般逆行列を計算している。
> library(MASS) > A <- rbind(c(1, -1, 0), c(0, 1, 1)) # 行列Aは、singular matrix > solve(A) # 普通の逆行列は求まらない > G <- ginv(A) # 行列Gは、ムーア・ペンローズ逆行列 > A %*% G %*% A # 条件(i)の確認 > G %*% A %*% G # 条件(ii)の確認 > t(G %*% A) # 条件(iii)の確認 > t(A %*% G) # 条件(iv)の確認
> x1 <- 1:5
> x2 <- runif(5)
> x3 <- rnorm(5)
> (x <- rbind(x1, x2, x3))
[,1] [,2] [,3] [,4] [,5]
x1 1.000000 2.0000000 3.0000000 4.00000000 5.00000000
x2 0.494798 0.7612015 0.5567693 0.01437583 0.76684719
x3 1.430874 -0.7378893 0.8044608 0.96198322 0.09665524
> max.col(x)
[1] 5 5 1
> c(which.max(x1), which.max(x2), which.max(x3)) # ベクトルの最大値位置は which.max
[1] 5 5 1
> apply(x, 1, max) # 最大値そのものを見つけるには例えばこうする
x1 x2 x3
5.0000000 0.7668472 1.4308745
ベクトルを次元属性を与えて行列化するように、リストも行列属性を与えれば行列としてアクセスできる。
> tm <- list(1, 2, "a", "b")
> ttm <- as.matrix(tm) # 行列化する
> dim(ttm) <- c(2, 2) # 2x2 行列にする
> ttm
[,1] [,2]
[1,] 1 "a"
[2,] 2 "b"
> is.list(ttm) # リストである
[1] TRUE
> is.matrix(ttm) # 行列でもある
[1] TRUE
> ttm[[1]] # リストとしてアクセス
[1] 1
> ttm[1, 1] # 行列としてアクセス
[[1]]
[1] 1
> ttm[1, 1][[1]]
[1] 1
> ttm[2, 2]
[[1]]
[1] "b"
> ttm[2, 2][[1]]
[1] "b"
行列 x を x = (eigen(x)$vectors) %*% (eigen(x)$values) %*% t(eigen(x)$vectors) と対角化した上で、(数値的許容度範囲で)負の対角成分を(数値的許容度範囲で)零に置き換えた xx を返す関数 (使い道は？)
> make.positive.definite <-
function(x, tol=1e-6) { # tol は数値的許容度
eig <- eigen(x, symmetric=TRUE)
rtol <- tol * eig$values[1]
if(min(eig$values) < rtol) {
vals <- eig$values
vals[vals < rtol] <- rtol
srev <- eig$vectors %*% (vals * t(eig$vectors))
dimnames(srev) <- dimnames(x)
return(srev)
} else {return(x)}
}
> x <- matrix(runif(4),ncol=2)
> x <- (x+t(x))/2 # 対称行列を作る
> y <- eigen(x)
> y$values
[1] 0.6038513 -0.3303502 # 固有値の一つが負だから正定値行列でない
> xx <- make.positive.definite(x)
> xx
[,1] [,2]
[1,] 0.2421127 0.2959413
[2,] 0.2959413 0.3617392
> yy <- eigen(xx)
# xx の正固有値は x のそれと一致、x の負固有値は(数値許容度範囲内で)零で置き換え
> yy$values
[1] 6.038513e-01 6.038513e-07
> eigen(x)$vector-eigen(xx)$vector # x,xx の固有ベクトル行列は数値許容度範囲内で一致
[,1] [,2]
[1,] 1.110223e-16 0.000000e+00
[2,] 0.000000e+00 -1.110223e-16
> z <- x %*% t(x) # 正定値符合行列を作る
> eigen(z)$values # 固有値は当然皆正
[1] 0.3646363 0.1091312
> z - make.positive.definite(z) # 正定値符合行列は不変
[,1] [,2]
[1,] 0 0
[2,] 0 0
正方行列 A に対する巾乗操作 A^2 は成分毎に巾乗した行列を与え、行列そのものの巾乗 A%*%A を与えない。A の任意巾乗 A^n を計算するには A の固有値分解 A = V %*% D %*% t(V) を用いて、A^n = V %*% D^n %*% t(V) と計算する。D は対角行列だから D^n と D の行列としての n 乗巾は一致する。例えば次の関数を定義して使う。
matpow <- function(A, pow=2) {
y <- eigen(A)
y$vectors %*% diag( (y$values)^pow ) %*% t(y$vectors)
}
発展問題として行列の指数関数を定義するには次の関数を定義すれば良い。
> matexp <- function(A) {
y <- eigen(x)
y$vectors %*% diag( exp(y$values) ) %*% t(y$vectors)
}
> matexp(x) %*% matexp(-x) # チェック１ (数値誤差の範囲で単位行列)
[,1] [,2] [,3] [,4]
[1,] 1.000000e+00 -6.788190e-16 -1.338827e-15 1.399597e-15
[2,] 4.406265e-16 1.000000e-00 -1.397808e-15 9.696020e-16
[3,] 4.826733e-17 -7.704883e-16 1.000000e+00 8.359199e-17
[4,] 1.034085e-15 -7.592668e-16 -1.682140e-15 1.000000e+00
> matexp(x) %*% matexp(x) - matexp(2*x) # チェック２
[,1] [,2] [,3] [,4]
[1,] -7.105427e-15 -5.329071e-15 -3.552714e-15 -3.552714e-15
[2,] -1.776357e-15 -1.776357e-15 -3.552714e-15 0.000000e+00
[3,] -3.552714e-15 -3.552714e-15 -3.552714e-15 -1.776357e-15
[4,] -3.552714e-15 0.000000e+00 -3.552714e-15 -7.105427e-15
> matexp(diag(rep(0,4))) # 零行列の指数関数は単位行列
[,1] [,2] [,3] [,4]
[1,] 1 0 0 0
[2,] 0 1 0 0
[3,] 0 0 1 0
[4,] 0 0 0 1
この定義は A が対称でなくても意味を持つが、結果は(おそらく)無意味。
例えば、X の４乗を (X %*% X) %*% (X %*% X) として計算し高速化している。
matPower <- function(X,n)
## Function to calculate the n-th power of a matrix X
{
if(n != round(n)) {
n <- round(n)
warning("rounding exponent `n' to", n)
}
phi <- diag(nrow = nrow(X))
pot <- X # the first power of the matrix.
while (n > 0)
{
if (n %% 2)
phi <- phi %*% pot
n <- n %/% 2
pot <- pot %*% pot
}
return(phi)
}
二つの行列の全体としての同等性を検査するには all と any を使う。返り値は単一の論理値なので if 文中でも使える。
> x <- y <- matrix(1:16, 4, 4)
> all(x == y) # 要素がすべて等しいか
[1] TRUE
> y[1,2] <- 4 # 一箇所変更
> all(x == y) # 要素がすべて等しいか
[1] FALSE
> any(x == y) # 要素が一箇所でも等しいか
[1] TRUE
> x == y # この構文は要素毎の同等性を与える行列を返す。
[,1] [,2] [,3] [,4]
[1,] TRUE FALSE TRUE TRUE
[2,] TRUE TRUE TRUE TRUE
[3,] TRUE TRUE TRUE TRUE
[4,] TRUE TRUE TRUE TRUE
> if(x == y) cat ("x == y?n") # 間違い易い構文(x[1,1]==y[1,1] だけをチェック)
x == y
Warning message:
the condition has length > 1 and only the first element will be used in:
if (x == y) cat("x == y?n")
ifelse(x==y, T, F) はどのような結果を返すか（考えてから実行してみよ）
Q：以下のような行列を作成する。(R-helpメーリングリスト，Subject: [R] question about matrix，Cynthia He氏からの質問)
0^1 0^2 ......... 0^n 1^1 1^2 ......... 1^n 2^1 2^2 ......... 2^n ................................ ................................ q^1 q^2 ......... q^n
> outer(0:5, 1:4, "^")
[,1] [,2] [,3] [,4] [1,] 0 0 0 0 [2,] 1 1 1 1 [3,] 2 4 8 16 [4,] 3 9 27 81 [5,] 4 16 64 256 [6,] 5 25 125 625
> x <- diag(0, 6, 4) # 必要な大きさの行列を作る
> (row(x)-1)^col(x) # 上の方にある row と col を使う
[,1] [,2] [,3] [,4]
[1,] 0 0 0 0
[2,] 1 1 1 1
[3,] 2 4 8 16
[4,] 3 9 27 81
[5,] 4 16 64 256
[6,] 5 25 125 625
Q：以下のような行列を作成する。(R-helpメーリングリスト，Subject: [R] Questions about Matrix，duoduo chen氏からの質問)
(0-1)^n (0-2)^n ,,,, (0-m)^n (1-1)^n (1-2)^n ,,,, (1-m)^n (2-1)^n (2-2)^n ,,,, (2-m)^n .. ... ,,,, ....... (t-1)^n (t-2)^n ,,,, (t-m)^n
> whatisit <- function(t, m, n) outer(0:t, 1:m, function(x, y) (x-y)^n)
> whatisit(3, 2, 4)
[,1] [,2]
[1,] 1 16
[2,] 0 1
[3,] 1 0
[4,] 16 1
> x <- diag(0, 4, 2)
> (row(x)-1-col(x))^4
[,1] [,2]
[1,] 1 16
[2,] 0 1
[3,] 1 0
[4,] 16 1
以下のような行列を作ろう（答えはコメントアウトしてある）。
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1 1 1 1 1 1 1 1 1 1
[2,] 1 2 2 2 2 2 2 2 2 2
[3,] 1 2 3 3 3 3 3 3 3 3
[4,] 1 2 3 4 4 4 4 4 4 4
[5,] 1 2 3 4 5 5 5 5 5 5
[6,] 1 2 3 4 5 6 6 6 6 6
[7,] 1 2 3 4 5 6 7 7 7 7
[8,] 1 2 3 4 5 6 7 8 8 8
[9,] 1 2 3 4 5 6 7 8 9 9
[10,] 1 2 3 4 5 6 7 8 9 10
----------
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1 2 3 4 5 6 7 8 9 10
[2,] 2 2 3 4 5 6 7 8 9 10
[3,] 3 3 3 4 5 6 7 8 9 10
[4,] 4 4 4 4 5 6 7 8 9 10
[5,] 5 5 5 5 5 6 7 8 9 10
[6,] 6 6 6 6 6 6 7 8 9 10
[7,] 7 7 7 7 7 7 7 8 9 10
[8,] 8 8 8 8 8 8 8 8 9 10
[9,] 9 9 9 9 9 9 9 9 9 10
[10,] 10 10 10 10 10 10 10 10 10 10
アドオンパッケージ Matrix には Shur 分解や、疎行列(要素の大部分が 0 である巨大行列)を扱うクラスが定義されている。以下にそのインデックスの一覧を紹介する：
Cholesky-class Cholesky decompositions Hilbert Generate a Hilbert matrix LU-class LU decompositions Matrix Construct a Classed Matrix Matrix-class Class "Matrix" Virtual class of Matrices Schur Schur Decomposition of a Matrix ScotsSec Scottish secondary school scores coef<- Assign parameters corFactor Extract the correlation factor corMatrix Generate the correlation matrix corrmatrix-class Class "corrmatrix" cscBlocked-class Class "cscBlocked" compressed, sparse, blocked Matrix cscMatrix-class Class "cscMatrix" compressed, sparse, column Matrix expand Expand a Decomposition into Factors facmul Multiplication by Decomposition Factors geMatrix-class General S4 Matrix class lmeRep-class lme models representation lu Triangular Decomposition of a Square Matrix matrix<- Assign the matrix mm A sample sparse model matrix norm Norm of a Matrix pdFactor Square-root factor pdMatrix Return the positive-definite matrix pdfactor-class Factors of positive-definite matrices pdmatrix-class Positive-definite matrices poMatrix-class Positive semi-definite, dense, real matrices rcond Estimate the Reciprocal Condition Number sscChol-class Cholesky decompositions of sscMatrix objects sscCrosstab Create pairwise crosstabulation sscCrosstab-class Pairwise crosstabulations sscMatrix-class Symmetric, compressed, sparse column matrices ssclme-class lme models as sparse, symmetric matrices syMatrix-class Class "syMatrix" of symmetric matrices trMatrix-class Triangular, dense, non-packed, real matrices tripletMatrix-class Class "tripletMatrix" sparse matrices in triplet form tscMatrix-class Triangular, compressed, sparse column matrices unpack Full Storage Representation of Packed Matrices y A sample response vector
このような大行列の逆行列(もしくは線形方程式)を解くことなど、普通の人には無縁でしょうが、R でやればどうなるかは一寸興味がありますね。Matrix パッケージを使うと「何度も計算する場合は」格段に早いようです。
On a modest desktop system (Athlon XP 2500+, approx 1.8 GHz) it took about 20 seconds for the inversion and about 8 seconds for the LU decomposition. This is using the Matrix package. > mm <- Matrix(rnorm(2900*2900), ncol = 2900) # 2900x2900 行列 > object.size(mm) [1] 67280636 > system.time(rcond(mm)) # rcond 適用の際密かに LU 分解され、情報として付け加わる # (これ以降は mm そのものではなく、そのLU分解が使われる) [1] 7.73 0.10 7.96 0.00 0.00 > object.size(mm) [1] 134573052 > system.time(minv <- solve(mm)) [1] 19.93 0.10 20.04 0.00 0.00 # Matrix パッケージを使うと 20 秒かかる > rcond(mm) [1] 2.466351e-07 The initial call to rcond calculates and stores the LU decomposition. That is why the size of the object approximately doubles. I get similar results without the matrix package except that the decomposition must be recomputed for each such operation on the matrix. > m1 <- matrix(rnorm(2900*2900), ncol = 2900) > system.time(minv1 <- solve(m1)) [1] 30.08 0.86 32.32 0.00 0.00 # Matrix パッケージを使わないと 30 秒
