さて、前回の話の最後に、「逆行列」というをやった。
「逆行列」というのは、スカラーの「逆数」に対応し、
元の行列に、逆行列を掛けると、「単位行列」が出てくるのであった。
覚えていない人は、もう一度、逆行列の投稿を参照すること。
この投稿の最後に、「逆行列」が「正方行列」でないといけない、と言った。
ところが、この「正方行列」という制限は、色々と不便なのである。
しかし、方法が無いかと言われると、そうでもない。
近似的に「逆行列」を求める方法がある。
っと、いうことを今日は書いてみる。
では、早速、「R」を使って、この問題に挑戦してみる。
まずは、4行×3列の行列を準備する。つまり、n=4、p=3 の行列である。
「R」を立ち上げたら、次のように入力する。
A <- matrix(c(50,0.5,50,5,7,0.2,0,6,75,3,20,50), nrow=4, ncol=3)
今回は、語呂でもなく、特に意味の無いデータであるが、
とにかく、この行列をつかって、演習を行う。
何のエラーも出なければ成功。これで、Aという変数に行列が代入された。
どのようなデータであるのか、変数名を入力してエンターキーを押すと、
変数に格納されているデータが表示される。この例では、変数は「A」である。
> A
[,1] [,2] [,3]
[1,] 50.0 7.0 75
[2,] 0.5 0.2 3
[3,] 50.0 0.0 20
[4,] 5.0 6.0 50
ふむ。一応、逆行列が計算できるか、確認しみる。
逆行列のコマンドは何であったか？
> solve(A)
Error in solve.default(A) : only square matrices can be inverted
この式を「R」で計算するには、以下のようにする。「逆行列」というのは、スカラーの「逆数」に対応し、
元の行列に、逆行列を掛けると、「単位行列」が出てくるのであった。
覚えていない人は、もう一度、逆行列の投稿を参照すること。
この投稿の最後に、「逆行列」が「正方行列」でないといけない、と言った。
ところが、この「正方行列」という制限は、色々と不便なのである。
しかし、方法が無いかと言われると、そうでもない。
近似的に「逆行列」を求める方法がある。
っと、いうことを今日は書いてみる。
では、早速、「R」を使って、この問題に挑戦してみる。
まずは、4行×3列の行列を準備する。つまり、n=4、p=3 の行列である。
「R」を立ち上げたら、次のように入力する。
A <- matrix(c(50,0.5,50,5,7,0.2,0,6,75,3,20,50), nrow=4, ncol=3)
今回は、語呂でもなく、特に意味の無いデータであるが、
とにかく、この行列をつかって、演習を行う。
何のエラーも出なければ成功。これで、Aという変数に行列が代入された。
どのようなデータであるのか、変数名を入力してエンターキーを押すと、
変数に格納されているデータが表示される。この例では、変数は「A」である。
> A
[,1] [,2] [,3]
[1,] 50.0 7.0 75
[2,] 0.5 0.2 3
[3,] 50.0 0.0 20
[4,] 5.0 6.0 50
ふむ。一応、逆行列が計算できるか、確認しみる。
逆行列のコマンドは何であったか？
> solve(A)
Error in solve.default(A) : only square matrices can be inverted
やはり、エラーが出る。「正方行列」でないと計算できない。
困ったことになった。どうしようも無いのか？
仕方が無いので、少し、物の見方というのを変えてみる。
何とか、元の行列を正方行列に変換できないのか？
つまり、という見方である。そのようにすれば、計算できるかもしれない。
ふむ。このような方法は使えるだろうか？
元の行列を転置し、その行列に元の行列を掛けてみる。
説明が解りにくい。これは、文系であっても、数式で示した方が理解しやすい。
そういった事はよくある。完全に数式を避けるべきではない。まぁ、それは良い。
要するに、
ということである。
要するに、
ということである。t(A) %*% A
以下は計算結果である。
> t(A) %*% A
[,1] [,2] [,3]
[1,] 5025.25 380.10 5001.5
[2,] 380.10 85.04 825.6
[3,] 5001.50 825.60 8534.0
なるほど、「正方行列」の形になっている。
これを使って「逆行列」を「擬似的」に計算するには、
以下のように計算する。
したがって、「R」では以下のようになる。
> solve(t(A)%*%A)%*%t(A)
[,1] [,2] [,3] [,4]
[1,] 0.07946341 -0.2671138 -0.04841192 -0.08380351
[2,] 1.54547248 -5.3787729 -1.34597554 -1.45709213
[3,] -0.18729532 0.6772539 0.16092918 0.19593608
ふむ、なるほど。よく見れば、上記の式と対応関係にあることが分かる。
これが「擬似的」に求めた「逆行列」である...本当か？疑いたくなる。
どうすれば、確認ができるのだろうか？確か、何かを掛けると、何かになった。
逆行列が逆数に対応するのだったことを思い出す。
> solve(t(A)%*%A)%*%t(A) %*% A
[,1] [,2] [,3]
[1,] 1.000000e+00 -6.095124e-14 -7.318590e-13
[2,] -7.414513e-12 1.000000e+00 -8.988366e-12
[3,] 2.988720e-13 -3.286260e-14 1.000000e+00
つまり、元の行列を掛けると、「単位行列」が出てくるのであった。
なるほど、確かに、対角成分が全て「1」となっていて、
それ以外の成分が全て「0」になっている。逆行列になっている。
このように、擬似的に求めた逆行列のことを「擬逆行列」と呼ぶ。
これは、「正方行列」でなくても使うことができる。これは便利である。
p > n の場合には、式が次のようになるので注意が必要。
行列の掛け算は、掛ける方と掛けられる方の順番が重要なのであった。
n > p の場合：
n ＜ p の場合：
実を言うと、「R」では、「擬逆行列」を計算する方法がいくつかある。
これを利用するには、「MASS」というパッケージを読み込む方法がある。
とにかく、次のようにする。
library(MASS)
これで、MASSというライブラリを読み込むことができた。
「R」というソフトウェアは、このようにライブラリと呼ばれるものを読み込み、
様々な機能を追加することができるのである。
このMASSというライブラリを読み込むと擬逆行列を計算する関数が使える。
その関数とは、ginv()関数である。とにかく、使ってみる。
ginv(A)
実行すると、次のようになる。
> ginv(A)
[,1] [,2] [,3] [,4]
[1,] 0.07946341 -0.2671138 -0.04841192 -0.08380351
[2,] 1.54547248 -5.3787729 -1.34597554 -1.45709213
[3,] -0.18729532 0.6772539 0.16092918 0.19593608
先に、手計算した結果と見比べてみると、なるほど、同じである。
回りくどいことせずとも、逆行列が計算できる。便利である。これを使おう。
この計算が、内部的にされているのだろう...
と思ってはいけない。実は、少々、手の混んだことをやっている。
実は、「擬逆行列」というのは、「特異値分解」からも導くことができ、
ginv()関数では、「特異値分解」を使って計算しているのである。
まだ、「特異値分解」の説明を見ていない人は、
とにかく、先にそっちを確認すること。
さて、特異値分解を行うと、ある行列 X は次のように分解されるのであった。
さて、特異値分解を行うと、ある行列 X は次のように分解されるのであった。
Uは、「左特異ベクトル」と呼ばれる行列であり、
Dは、「特異値行列」と呼ばれる対角行列。そして、
V'は、「右特異ベクトル」と呼ばれる行列、である。
あまり、深く考えてはいけない。スカラーの「因数分解」のようなものである。
そういう分解をしているのである。つまり、無秩序な任意の行列を、
ある「規則的な行列同士」の「掛け算」に分解しているのである。
複雑な事象を、本質を変えずに、見通しの良い状態にしている、と言っても良い。
ということがニュアンスとして理解できていたら、ここでは良いだろう。
とにかく、この分解された要素の組み合わせ次第で、様々な、難問が解ける。
まさに、数学の魔法が詰まった分解なのである。
擬逆行列の計算というものは、特異値分解の魔法を応用した一例にすぎないのだが、
分解された要素を次のように組み替えれば、擬逆行列となる。
では、「R」で行うとどのようになるのか、やってみる。
まず、変則的な命名方法であるが、A.svd という変数をつくり、
この変数に特異値分解の結果を格納する。
A.svd <- svd(A)
特異値分解は、svd()関数を使えば簡単に計算できる。
また、計算結果がどのように格納されているかは、次の関数で確認できる。
str(A.svd)
つまり、構造（structure）を確認する関数であり、str() というのは、
英語の最初の三文字を取っている。プログラミングでは、このような省略を用いる。
長い変数名は、バグの遠因であるし、また、コーディングも大変である。
実は、変数の命名法にもいくつか流儀というのがあるのだが、この話はいずれ。
すぐに、話が逸れてしまうので困る。私の性格の問題ではあるが。
できれば、愛嬌とか、ユーモアで片付けてもらいたい。まぁ、良いのだが。
つまり、str()関数 の実行結果は次のようになる。
> str(A.svd)
List of 3
$ d: num [1:3] 110.209 38.706 0.167
$ u: num [1:4, 1:3] -0.82 -0.0249 -0.409 -0.3995 -0.0647 ...
$ v: num [1:3, 1:3] -0.5758 -0.0739 -0.8142 0.8161 -0.1118 ...
そして、A.svd という変数に格納されている d, u, v の3つの要素は以下のように取り出す。
> # 特異値行列の抽出
> A.svd$d
[1] 110.2092350 38.7064165 0.1669015
> # 左特異ベクトルの抽出
> A.svd$u
[,1] [,2] [,3]
[1,] -0.82003552 -0.06468901 0.2601641
[2,] -0.02491042 -0.03398276 -0.9059111
[3,] -0.40900495 0.76122378 -0.2263649
[4,] -0.39954495 -0.64435926 -0.2457616
> # 右特異ベクトルの抽出
> A.svd$v
[,1] [,2] [,3]
[1,] -0.5758338 0.8160908 0.04910397
[2,] -0.0738822 -0.1117586 0.99098508
[3,] -0.8142216 -0.5670148 -0.12464897
ふむ。では、まずは、以下のようにしてみよう。
D <- diag(A.svd$d)
V <- A.svd$v
U <- A.svd$u
こうして、3つの変数、D, V, U ができた。
なお、Dに関しては、ベクトルとして出力されているため、
対角行列を作るための diag()関数 を使って、行列に変換している。
さて、ここまで遠回りが続いているが、
公式に基づいて、「擬逆行列」を計算すると以下のようになる。
> V%*%solve(D)%*%t(U)
[,1] [,2] [,3] [,4]
[1,] 0.07946341 -0.2671138 -0.04841192 -0.08380351
[2,] 1.54547248 -5.3787729 -1.34597554 -1.45709213
[3,] -0.18729532 0.6772539 0.16092918 0.19593608
見事に、擬逆行列が計算されているのが分かる。
特異値分解は、後に主成分分析や対応分析にも登場するので、
今の段階で、どういった特徴を持つ行列であるか、理解しておきたい。
0 件のコメント:
コメントを投稿