2012/08/10

文系のための「擬逆行列」

さて、前回の話の最後に、「逆行列」というをやった。
逆行列」というのは、スカラーの「逆数」に対応し、
元の行列に、逆行列を掛けると、「単位行列」が出てくるのであった。

覚えていない人は、もう一度、逆行列の投稿を参照すること。
この投稿の最後に、「逆行列」が「正方行列」でないといけない、と言った。

ところが、この「正方行列」という制限は、色々と不便なのである。
しかし、方法が無いかと言われると、そうでもない。
近似的に「逆行列」を求める方法がある。
っと、いうことを今日は書いてみる。

では、早速、「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」で計算するには、以下のようにする。

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」になっている。逆行列になっている。

このように、擬似的に求めた逆行列のことを「擬逆行列」と呼ぶ。
これは、「正方行列」でなくても使うことができる。これは便利である。

ちなみに、今回の例では、n > p の行列であったが、
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()関数では、「特異値分解」を使って計算しているのである。

まだ、「特異値分解」の説明を見ていない人は、
とにかく、先にそっちを確認すること。

さて、特異値分解を行うと、ある行列 は次のように分解されるのであった。


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 件のコメント:

コメントを投稿