行列の分解と正規化
この節では、固有値および固有ベクトルを計算する関数と、行列を分解する関数について説明します。
固有値
Eigen()関数は、対称行列の固有値分解を実行します。固有値分解は、主成分分析と正準相関分析をはじめとするさまざまな統計的手法で用いられています。主成分分析と正準相関分析は、最大の固有値に対応した変換が分散を最大化する変換となっています。
Eigen()は行列のリストを戻します。戻されたリストの最初の行列は固有値の列ベクトルで、2番目の行列には列として固有ベクトルが含まれています。
A = [ 5 4 1 1, 4 5 1 1, 1 1 4 2, 1 1 2 4];
Eigen( A );
{[10, 5, 2, 1][0.632455532033676 - 0.316227766016838 - 2.77555756156289e-16 -0.707106781186547, 0.632455532033676 - 0.316227766016837 - 1.66533453693773e-16 0.707106781186547, 0.316227766016838 0.632455532033676 0.707106781186548 0, 0.316227766016837 0.632455532033676 - 0.707106781186547 0]}
この関数は行列のリストを結果として戻します。結果を受け取るには、一方には固有値のベクトル、もう一方には固有ベクトルの行列が代入されるように、2つのグローバル変数のリストに代入してください。
{evals, evecs} = Eigen( A );
固有値分解をすることで、n x nの行列Aに対する方程式Ax = λxがゼロベクトルでない解xを持つときの、すべてのλ（ラムダ）とベクトルxが求められます。このλを固有値、xベクトルを固有ベクトルと呼びます。これは(A - λI)x = 0を解くことと同じです。下のようなスクリプトで、固有値と固有ベクトルからAを再構成できます。
newA = evecs * Diag( evals ) * evecs`;
[5.00000000000001 4 1 1,
4 5 1 1,
1 1 4 2,
1 1 2 4]
固有値および固有ベクトルに関しては、次のことに注意してください。
• 固有ベクトルの行列は直交行列で、転置行列が逆行列になります（
）。
• 固有値が一意であるときのみ、固有ベクトルも一意に求められます。
• 固有値がゼロの行列は、特異行列です。
• 固有値の逆行列と固有ベクトルにより、逆行列も求められます。Moore-Penroseの一般化逆行列は、ゼロ以外の固有値の逆行列から求められます（GInverseを参照）。
メモ: 非常に小さい固有値を0とみなすかどうかを判断する必要があります。
• 固有値分解をすることで、正方行列の掛け算を次のように考えることができます。
– 回転（直交行列を掛ける）
– スケーリング（対角行列を掛ける）
– 逆回転（直交行列の逆行列、つまり転置行列を掛ける）
A * x = E` * Diag( M ) * E * x;
Eで回転、Diag( M )でスケーリング、E`で逆回転します。
Cholesky分解
Cholesky()関数はCholesky分解を行います。半正定値行列Aを、非特異下三角行列Lとその転置行列の積の形、L*L' = Aに書き直します。
E = [ 5 4 1 1, 4 5 1 1, 1 1 4 2, 1 1 2 4];
L = Cholesky( E );
[2.23606797749979 0 0 0,
1.788854381999832 1.341640786499874 0 0,
0.447213595499958 0.149071198499986 1.9436506316151 0,
0.447213595499958 0.149071198499986 0.914659120760047 1.71498585142509]
結果を確認するには、次のように入力します。
L*L`;
[5 4 1 1,
4 5 1 1,
1 1 4 2,
1 1 2 4]
Cholesky分解について
Cholesky()は、行列の式を扱いやすい形式に再構成するのに便利です。たとえば、JMPでは固有値を求めるには、行列が対称行列である必要があります。したがって、対称行列の積ABの固有値は直接求められない場合があります（AとBが対称行列であっても、その積が対称行列であるとは限らないからです）。しかし、この問題は、L′BLの固有値の問題に置き換えることができます。ここで、LはAのCholesky根です。L`BLはABと同じ固有値を持ちます。
Cholesky()の他の使用法としては、Trace()式で行列の対角和を求めるときに、行列の順序を変えるということが考えられます。Trace(A*B*A`)などの式で、Aの行数が多い場合、大きな計算量を要することになります。しかし、BがLL′とCholesky分解できるのであれば、元の式はTrace(A*L*L`*A`)と書き直せます。これは、Trace(L`A`*AL)と等しく、したがってALの要素の平方和を求めるだけでよくなるので、演算量はかなり少なくなります。
Chol Update()関数を使用すると、Cholesky分解を効率的に更新できます。n × n行列AのCholesky根をLとした場合、cholUpdate(L, C, V)を呼び出すと、A+V*C*V'のCholesky根が算出されます。Cはm × mの対称行列、Vはn × m行列です。
例
Cholesky分解を手動で更新する手順は次のとおりです。
exS = [16 1 0 11 -1 12, 1 11 -1 1 -1 1, 0 -1 12 -1 1 0, 11 1 -1 11 -1 9, -1 -1 1 -1 9 -1, 12 1 0 9 -1 12];
exAchol = Cholesky( exS ); // Cholesky分解を実行する
// 計画行列に2つの列ベクトルを追加する
exV = [1 1, 0 0, 0 1, 0 0, 0 0, 0 1];
/* 最初の列ベクトルは、計画行列の行の1つに足される。
2番目の列ベクトルは、計画行列の行の1つから引かれる */
exC = [1 0, 0 -1];
exAnew = exS + exV * exC * exV`;
// Cholesky分解を手動で更新する
exAcholnew = Cholesky( exAnew );
手動で更新するのではなく、Chol Update()を使ってCholesky分解をより効率的に更新する手順は、次のとおりです。
// Cholesky分解をより効率的に更新する
exAcholnew_test = Chol Update( exAchol, exV, exC );
Show( exAcholnew_test );
exAcholnew_test =
[ 4 0 0 0 0 0,
0.25 3.30718913883074 0 0 0 0, ...]
// 結果は手動での手順の場合と同じ
Show( exAcholnew );
exAcholnew=
[ 4 0 0 0 0 0,
0.25 3.30718913883074 0 0 0 0, ...]
特異値分解
SVD()関数を使うと、行列の特異値分解()が求められます。つまり、行列Aについて、SVD()はU*diag(M)*V`=Aを満たす3つの行列、U、M、Vのリストを戻します。
メモ:
• 結果の行列Mでは、余計な0の対角要素は省かれています。
• 特異値分解によって、AはUSV′の形で表されます。ここで、
– UとVは、互いに直交している列ベクトルを含んでいる行列です（「直交するベクトル」とは、「直角」もしくは「独立」なベクトルのことを指します）。
– Sはn × n対角行列で、Aの特異値、つまりA′Aの固有値の負でない平方根を要素として持っています。
• 特異値分解は対応分析などで使われています。
例
A = [1 2 1 0, 2 3 0 1, 1 0 1 5, 0 1 5 1];
{U, M, V} = SVD( A );
newA = U * Diag( M ) * V`;
[1 2 0.999999999999997 -2.99456640040496e-15,
2 3 -1.17505831453979e-15 1,
0.999999999999997 -2.16851276518826e-15 0.999999999999999 5,
2.22586706011274e-15 1 5 0.999999999999997]
正規直交化
Ortho()関数は、列の直交化を行い、さらにベクトルの大きさで割って正規化します。この関数はGram-Schmidt法を使用します。直交行列の列ベクトルは、大きさが単位長で互いに直交（内積が0）します。
B = Ortho( [1 -1, 1 0, 0 1] );
[0.408248290463863 -0.707106781186548, 0.408248290463863 0.707106781186548, -0.816496580927726 3.14018491736755e-16]
各列ベクトルが直交していることは、BとBの転置行列を掛け合わせると単位行列になることから確かめられます。
C = B` * B;
[1 -3.119061760824e-16, -3.119061760824e-16 1]
特に指定がない限り、ベクトルは正規化されます。つまり、ベクトルの大きさでベクトルそのものを割るスケーリングで、大きさ1のベクトルを作ります。Scaled(0)を指定して、スケーリングをしない場合は、次のようになります。
Ortho( [1 -1, 1 0, 0 1], Scaled( 0 ) );
[0.408248290463863 -0.353553390593274, 0.408248290463863 0.353553390593274, -0.816496580927726 1.57009245868377e-16]
全要素の和が0になるベクトルを作るには、オプションのCentered(1)を指定します。このオプションは対比の行列を作る場合に便利です。
result = Ortho( [1 -1, 1 0, 0 1], Centered( 1 ) );
[0.408248290463863 -0.707106781186548, 0.408248290463863 0.707106781186548, -0.816496580927726 3.14018491736755e-16]
各列の要素の和が0になることは、要素がすべて1のベクトルを前から掛けて列の要素の和を求めることで確認できます。
J(1, 3) * result;
[1.11022302462516e-16 2.02996189274239e-16]
直交多項式
Ortho Poly()関数は、引数で指定された次元までのベクトルの直交多項式を戻します。を戻します。直交多項式は、多項式モデルを普通にあてはめると回帰係数の間に強い相関が生じてしまう場合に有用です。
Ortho Poly( [1 2 3], 2 );
[-0.707106781186548 0.408248290463862, 0 -0.816496580927726, 0.707106781186548 0.408248290463864]
次元はベクトルの次元より小さいものである必要があります。正規直交化の項にあるように、Scaled(1)を指定することで単位長のベクトルを作れます。
QR分解
QR()分解は、数値的に安定した行列の処理に便利です。QR()は2つの行列のリストを戻します。以下はその例です。
{Q, R} = QR( [11 22, 33 44] );
QとRに結果が入ります。m × n行列Xに対して、QR()はX = Q*Rと分解します。ここで、Qは、m × mの直交行列、Rはm × nの上三角行列です。
逆行列の更新
M`M行列の逆行列に行を追加または削除するには、Inv Update(S, X, w)関数を使用します。逆行列の更新は、1行除去による影響診断や候補計画の評価に役立ちます。
メモ:
• 第1引数Sは更新する行列です。
• 第2引数Xは、追加または削除する行を含む行列です。
• 第3引数wには、行を追加する場合は1、削除する場合は-1を指定します。
• 複数の行を追加または削除できます。
Inv Update(S, X, w)関数は、次の式を計算するのと同じです。
S - w * S * X` * Inv( I + w * X * S * X` ) * X * S
ここで、Iは単位行列、Inv(A)はAの逆行列です。