atwiki-logo
  • 新規作成
    • 新規ページ作成
    • 新規ページ作成(その他)
      • このページをコピーして新規ページ作成
      • このウィキ内の別ページをコピーして新規ページ作成
      • このページの子ページを作成
    • 新規ウィキ作成
  • 編集
    • ページ編集
    • ページ編集(簡易版)
    • ページ名変更
    • メニュー非表示でページ編集
    • ページの閲覧/編集権限変更
    • ページの編集モード変更
    • このページにファイルをアップロード
    • メニューを編集
    • 右メニューを編集
  • バージョン管理
    • 最新版変更点(差分)
    • 編集履歴(バックアップ)
    • アップロードファイル履歴
    • ページ操作履歴
  • ページ一覧
    • ページ一覧
    • このウィキのタグ一覧
    • このウィキのタグ(更新順)
    • このページの全コメント一覧
    • このウィキの全コメント一覧
    • おまかせページ移動
  • RSS
    • このウィキの更新情報RSS
    • このウィキ新着ページRSS
  • ヘルプ
    • ご利用ガイド
    • Wiki初心者向けガイド(基本操作)
    • このウィキの管理者に連絡
    • 運営会社に連絡(不具合、障害など)
ページ検索 メニュー
R入門
  • ウィキ募集バナー
  • 目安箱バナー
  • 操作ガイド
  • 新規作成
  • 編集する
  • 全ページ一覧
  • 登録/ログイン
ページ一覧
R入門
  • ウィキ募集バナー
  • 目安箱バナー
  • 操作ガイド
  • 新規作成
  • 編集する
  • 全ページ一覧
  • 登録/ログイン
ページ一覧
R入門
ページ検索 メニュー
  • 新規作成
  • 編集する
  • 登録/ログイン
  • 管理メニュー
管理メニュー
  • 新規作成
    • 新規ページ作成
    • 新規ページ作成(その他)
      • このページをコピーして新規ページ作成
      • このウィキ内の別ページをコピーして新規ページ作成
      • このページの子ページを作成
    • 新規ウィキ作成
  • 編集
    • ページ編集
    • ページ編集(簡易版)
    • ページ名変更
    • メニュー非表示でページ編集
    • ページの閲覧/編集権限変更
    • ページの編集モード変更
    • このページにファイルをアップロード
    • メニューを編集
    • 右メニューを編集
  • バージョン管理
    • 最新版変更点(差分)
    • 編集履歴(バックアップ)
    • アップロードファイル履歴
    • ページ操作履歴
  • ページ一覧
    • このウィキの全ページ一覧
    • このウィキのタグ一覧
    • このウィキのタグ一覧(更新順)
    • このページの全コメント一覧
    • このウィキの全コメント一覧
    • おまかせページ移動
  • RSS
    • このwikiの更新情報RSS
    • このwikiの新着ページRSS
  • ヘルプ
    • ご利用ガイド
    • Wiki初心者向けガイド(基本操作)
    • このウィキの管理者に連絡
    • 運営会社に連絡する(不具合、障害など)
  • atwiki
  • R入門
  • 行列

R入門

行列

最終更新:2025年05月26日 20:33

r-intro

- view
管理者のみ編集可

目次

  • 目次
  • 基本操作
    • 行列を作る
    • 行列を作る
    • 行列から成分を取り出す
    • 行列から特定の行か特定の列のみ取り出す
    • 行列をベクトルに変換する
    • 行列から非対角成分だけを取り出す
    • 行列に行を追加する
    • 行列に列を追加する
  • 線形代数
    • 転置行列を求める
    • 単位行列を作る
    • 対角行列を作成する
    • 行列のトレース(固有和)を求める
    • 行列式を求める
    • 行列式を求める
    • 行列の掛け算
    • 逆行列を求める
    • 一般逆行列を求める
    • 連立一次方程式を解く
  • 計算
    • QR分解をする
    • コレスキー分解する
    • 固有値と固有ベクトルを求める
    • 固有値と固有ベクトルを求める
    • 特異値分解を行う
    • 余因子を求める
    • 余因子行列を求める
    • ヘッセ行列を求める
    • 勾配ベクトルを求める

基本操作

行列を作る

matrix関数を使う。行列の要素を格納したベクトルを第一引数に与えて、nrowオプションで行数を、ncolオプションに列数を与える。デフォルトでは列方向に順番に要素が格納される。これを行方向に順に格納するにはbyrowオプションにTRUEを指定する。

> d <- 1:12
> mx <- matrix(d, nrow = 3, ncol = 4)
> print(mx)
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12
> mx <- matrix(d, ncol = 4)
> print(mx)
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12
> mx <- matrix(d, ncol = 5)
 警告メッセージ: 
 matrix(d, ncol = 5) で: 
   データ長 [12] が列数 [5] を整数で割った、もしくは掛けた値ではありません 
> mx <- matrix(d, ncol = 4, byrow = TRUE)
> print(mx)
     [,1] [,2] [,3] [,4]
[1,]    1    2    3    4
[2,]    5    6    7    8
[3,]    9   10   11   12
> mx <- matrix(1:4, nrow = 2, ncol = 4, byrow = TRUE)
> print(mx)
     [,1] [,2] [,3] [,4]
[1,]    1    2    3    4
[2,]    1    2    3    4
> mx <- matrix(0., 2, 5)
> print(mx)
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    0    0    0    0
[2,]    0    0    0    0    0

上記の例のとおりに、特に何も指定しなければ第二引数が行列、第三引数が列数になる。ベクトルの要素数で行列どちらかが決まればもう片方も決まるため、行または列の省略は可能。ただし、ベクトルの要素数に応じた正しい値でなければならない。

最後の例のとおりに、指定した行数×列数の要素数に対してベクトルの要素数が足らないときは、そのベクトルの要素数が行数×列数の約数であれば、自動的に繰り返し使われる。

行列を作る

行列はmatrix関数を使用して、ベクトルから作ることができる。以下は、2行3列の行列を新規に作成している。

> d <- 1:6
> mx <- matrix(d, 2, 3)
> mx
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

matrix関数の第一引数には、成分(並べられている各々の数)に代入する値が含まれているベクトルを指定し、特に指定をしなければ第二引数には行数を、第三引数には列数を指定する。行数はnrowオプション、列数はncolオプションで明示的に指定することができ(順不同)、例えば、以下のようにしても同じ。

> mx <- matrix(d, ncol = 3, nrow = 2)
> mx
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

特にオプションを指定しない場合、第一引数に与えたベクトルの要素について、列方向(縦方向)を優先に値を埋めるように代入する。これは見た目の直感となかなか合わない。列方向ではなく、行方向(横方向)を優先に値を埋めるように代入させるには、byrowオプションにTRUEを与える。

> mx <- matrix(d, 2, 3, byrow = TRUE)
> mx
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6

行列から成分を取り出す

行ごと列ごとに角括弧(ブラケット)をつかってインデックス(1~)を指定する。インデックスは1から始まることに注意。

> mx <- matrix(1:6, 3, 2)
> print(mx)
     [,1] [,2]
[1,]    1    4
[2,]    2    5
[3,]    3    6
> mx[1, 2]
[1] 4
> mx[3, 2]
[1] 6
> mx[4, 2]
 mx[4, 2] でエラー:  添え字が許される範囲外です

最後の例のように、行列の各長さより大きい値を指定すると、エラーが発生する。

行列から特定の行か特定の列のみ取り出す

インデックスを指定する際に、数値を与えずに空欄のままにする。以下の例ではそれぞれ2行目のみ、3列目のみを取り出している。

> mx <- matrix(1:9, 3, 3)
> print(mx)
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9
> mx[2, ]
[1] 2 5 8
> mx[, 3]
[1] 7 8 9
> is.vector(mx[, 3])
[1] TRUE

最後のとおりに、1行(列)だけ取り出したものはベクトルになる。これを行列のままにする場合はdropオプションにFALSEを与える。

> mx[2, , drop = FALSE]
     [,1] [,2] [,3]
[1,]    2    5    8

複数行(列)取り出したい場合は、インデックスをベクトルで与える。以下は2行目と3行目のみを取り出した例。

> mx[c(2, 3), ]
     [,1] [,2] [,3]
[1,]    2    5    8
[2,]    3    6    9

行列をベクトルに変換する

as.vector関数を使用する。

> mx1 <- matrix(6, nrow = 1, ncol = 1)
> mx2 <- matrix(1:9, nrow = 3, ncol = 3)
> mx1
     [,1]
[1,]    6
> as.vector(mx1)
[1] 6
> mx2
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9
> as.vector(mx2)
[1] 1 2 3 4 5 6 7 8 9
> as.vector(t(mx2))
[1] 1 4 7 2 5 8 3 6 9

1列目の要素すべて、2列目の要素すべて、・・・、という順にベクトルの要素になる。これを1行目の要素すべて、2行目の要素すべて、・・・、としたい場合は、最後の例のとおりにt関数を使って元の行列の転置行列を指定する。

行列から非対角成分だけを取り出す

row関数とcol関数を使うと、ベクトル形式で簡単に取り出すことができる。

> mx <- matrix(1:16, 4, 4)
> mx
     [,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    2    6   10   14
[3,]    3    7   11   15
[4,]    4    8   12   16
> mx[row(mx) != col(mx)]
 [1]  2  3  4  5  7  8  9 10 12 13 14 15

応用として、上の例で2、7、12(対角成分のそれぞれ1行下の成分)を取り除くには、以下のようにする。

> mx[row(mx) != col(mx) + 1]
 [1]  1  3  4  5  6  8  9 10 11 13 14 15 16

行列に行を追加する

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 番目の引数を参照)

行列に列を追加する

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 番目の引数を参照)

線形代数

転置行列を求める

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

単位行列を作る

diag関数を使う。引数を一つだけ指定した場合は、n次の単位行列を作成する。

> diag(3)
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1
> diag(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

対角行列を作成する

diag関数を使う。第1引数に対角要素の値、第2引数に求める正方行列の次数を指定する。

diag(4, 2)

    [,1] [,2]

[1,] 4 0 [2,] 0 4

diag(5, 3)

    [,1] [,2] [,3]

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

diag(c(2, 4), 4)

    [,1] [,2] [,3] [,4]

[1,] 2 0 0 0 [2,] 0 4 0 0 [3,] 0 0 2 0 [4,] 0 0 0 4 }}

行列のトレース(固有和)を求める

Rには行列Aのトレース(固有和)Tr(A)を求めるための関数が標準では組み込まれていない。diag関数とsum関数を組み合わせることで、計算を実現することができる。以下の例では、Tr(A) = 1+5+9=15である。

> mx <- matrix(1:9, nrow = 3, byrow = TRUE)
> mx
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
[3,]    7    8    9
> tr <- sum(diag(mx))
> tr
[1] 15

行列式を求める

det関数を使う。

> d <- c(1, 2, 3, 4)
> mx <- matrix(d, 2, 2, byrow = TRUE)
> mx
     [,1] [,2]
[1,]    1    2
[2,]    3    4
> det(mx)
[1] -2

行列式を求める

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

行列の掛け算

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

> 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 でエラー:  適切な引数ではありません

逆行列を求める

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

一般逆行列を求める

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

連立一次方程式を解く

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

計算

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

コレスキー分解する

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

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

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

特異値分解を行う

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ということだが、高精度に計算されていないことがわかる。

余因子を求める

polyMatrixパッケージのcofactor関数を使う。以下、使用例。計算に用いた行列はここを参考にしており、先に余因子行列をadjoint関数を使用して求めている。

> d <- c(8, 3, 4, 1, 5, 9, 6, 7, 2)
> mx <- matrix(d, 3)
> print(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
> cofactor(mx, 1, 2)
[1] 22
> cofactor(mx, 2, 3)
[1] -68

余因子行列を求める

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と見なしてよく、上記の性質を満たしていることが確認できた。

ヘッセ行列を求める

numDerivパッケージのhessian関数を使う。第一引数に目的関数を、第二引数に変数の値をベクトルで与える。

> library(numDeriv)
> f1 <- function(x) x[1] ^ 2 + 4 * x[2] ^ 2
> hessian(f1, c(0, 0))
              [,1]          [,2]
[1,]  2.000000e+00 -3.308722e-16
[2,] -3.308722e-16  8.000000e+00
> f2 <- function(x) x[1] ^ 2 + 2 * x[2] ^ 2
> hessian(f2, c(0, 0))
             [,1]         [,2]
[1,] 2.000000e+00 3.308722e-16
[2,] 3.308722e-16 4.000000e+00
> f3 <- function(x) 2 * x[1] - 4 * x[2] + x[1] ^ 2 + 2 * x[2] ^ 2 + 2 * x[1] * x[2]
> hessian(f3, c(0, 0))
     [,1] [,2]
[1,]    2    2
[2,]    2    4

勾配ベクトルを求める

numDerivパッケージのjacobian関数を使う。第一引数に目的関数を、第二引数に変数の値をベクトルで与える。

> library(numDeriv)
> f1 <- function(x) {6 -  0.1 * (x[1] - 5) ^ 2} * {6 - 0.04 * (x[2] - 5) ^ 2}
> jacobian(f1, c(8.25, 5.25))
          [,1]      [,2]
[1,] -3.898375 -0.098875
> f2 <- function(x) 2 * x[1] - 4 * x[2] + x[1] ^ 2 + 2 * x[2] ^ 2 + 2 * x[1] * x[2]
> jacobian(f2, c(1, 1))
     [,1] [,2]
[1,]    6    2
「行列」をウィキ内検索
LINE
シェア
Tweet
R入門
記事メニュー

メニュー

  • トップページ
  • Rとは
  • Rを使ってみる
  • 画面出力と入力
  • 変数とオブジェクト
  • ベクトル
  • 演算子と制御構文
  • 数と式
  • 文字と文字列
  • 正規表現
  • 日付と時刻
  • データフレーム
  • tidyverse
  • リスト
  • ファイルの入出力
  • ファイル・ディレクトリ
  • 実行とデバッグ
  • 数学
  • 行列
  • 関数
  • 統計学
  • 統計解析
  • 数値計算
  • 応用
  • 時系列解析
  • 地理空間情報
  • Rの操作
  • データベース
  • パッケージ

  • 環境
  • インターネット
  • 作図
  • 画像

  • 本の計算を再現
  • セイバーメトリクス
  • その他


ここを編集
記事メニュー2

更新履歴

取得中です。

ここを編集
人気記事ランキング
  1. インターネット
  2. その他
  3. tidyverse
  4. 応用
  5. 情報量統計学(共立出版)
  6. パッケージ
  7. 正規表現
もっと見る
最近更新されたページ
  • 14時間前

    その他
  • 3日前

    tidyverse
  • 19日前

    データフレーム
  • 22日前

    リスト
  • 46日前

    行列
  • 87日前

    インターネット
  • 98日前

    日付と時刻
  • 98日前

    応用
  • 98日前

    文字と文字列
  • 115日前

    変数とオブジェクト
もっと見る
人気記事ランキング
  1. インターネット
  2. その他
  3. tidyverse
  4. 応用
  5. 情報量統計学(共立出版)
  6. パッケージ
  7. 正規表現
もっと見る
最近更新されたページ
  • 14時間前

    その他
  • 3日前

    tidyverse
  • 19日前

    データフレーム
  • 22日前

    リスト
  • 46日前

    行列
  • 87日前

    インターネット
  • 98日前

    日付と時刻
  • 98日前

    応用
  • 98日前

    文字と文字列
  • 115日前

    変数とオブジェクト
もっと見る
ウィキ募集バナー
新規Wikiランキング

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

  1. MadTown GTA (Beta) まとめウィキ
  2. R.E.P.O. 日本語解説Wiki
  3. AviUtl2のWiki
  4. 機動戦士ガンダム EXTREME VS.2 INFINITEBOOST wiki
  5. シュガードール情報まとめウィキ
  6. ソードランページ @ 非公式wiki
  7. SYNDUALITY Echo of Ada 攻略 ウィキ
  8. シミュグラ2Wiki(Simulation Of Grand2)GTARP
  9. 星飼いの詩@ ウィキ
  10. ドラゴンボール Sparking! ZERO 攻略Wiki
もっと見る
人気Wikiランキング

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

  1. アニヲタWiki(仮)
  2. ストグラ まとめ @ウィキ
  3. ゲームカタログ@Wiki ~名作からクソゲーまで~
  4. 初音ミク Wiki
  5. 機動戦士ガンダム バトルオペレーション2攻略Wiki 3rd Season
  6. MadTown GTA (Beta) まとめウィキ
  7. 検索してはいけない言葉 @ ウィキ
  8. 発車メロディーwiki
  9. オレカバトル アプリ版 @ ウィキ
  10. Grand Theft Auto V(グランドセフトオート5)GTA5 & GTAオンライン 情報・攻略wiki
もっと見る
全体ページランキング

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

  1. 参加者一覧 - ストグラ まとめ @ウィキ
  2. 魔獣トゲイラ - バトルロイヤルR+α ファンフィクション(二次創作など)総合wiki
  3. ロスサントス警察 - ストグラ まとめ @ウィキ
  4. ダギ・イルス - 機動戦士ガンダム バトルオペレーション2攻略Wiki 3rd Season
  5. 光の黄金櫃(遊戯王OCG) - アニヲタWiki(仮)
  6. 召喚 - PATAPON(パタポン) wiki
  7. ステージ - PATAPON(パタポン) wiki
  8. 美食神アカシア - アニヲタWiki(仮)
  9. 可愛い逃亡者(トムとジェリー) - アニヲタWiki(仮)
  10. 箱入り娘(パズル) - アニヲタWiki(仮)
もっと見る

  • このWikiのTOPへ
  • 全ページ一覧
  • アットウィキTOP
  • 利用規約
  • プライバシーポリシー

2019 AtWiki, Inc.