Eigen()関数は、対称行列の固有値分解を実行します。固有値分解は、主成分分析と正準相関分析をはじめとするさまざまな統計的手法で用いられています。主成分分析と正準相関分析は、最大の固有値に対応した変換が分散を最大化する変換となっています。
Eigen()は行列のリストを戻します。戻されたリストの最初の行列は固有値の列ベクトルで、2番目の行列には列として固有ベクトルが含まれています。
A = [ 5 4 1 1, 4 5 1 1, 1 1 4 2, 1 1 2 4];
Eigen( A );
{evals, evecs} = Eigen( A );
固有値分解をすることで、n x nの行列Aに対する方程式Ax = λxがゼロベクトルでない解xを持つときの、すべてのλ(ラムダ)とベクトルxが求められます。このλを固有値、xベクトルを固有ベクトルと呼びます。これは(A - λI)x = 0を解くことと同じです。下のようなスクリプトで、固有値と固有ベクトルからAを再構成できます。
newA = evecs * Diag( evals ) * evecs`;
A * x = E` * Diag( M ) * E * x;
E = [ 5 4 1 1, 4 5 1 1, 1 1 4 2, 1 1 2 4];
L = Cholesky( E );
L*L`;
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行列です。
Cはm × mの対称行列、Vはn × m行列です。
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分解を実行する
exV = [1 1, 0 0, 0 1, 0 0, 0 0, 0 1];
// 計画行列に2つの列ベクトルを追加する
exC = [1 0, 0 -1];
/* 最初の列ベクトルは、計画行列の行の1つに足される。
2番目の列ベクトルは、計画行列の行の1つから引かれる */ */
exAnew = exS + exV * exC * exV`;
exAcholnew = Cholesky( exAnew );
// Cholesky分解を手動で更新する
手動で更新するのではなく、Chol Update()を使ってCholesky分解をより効率的に更新する手順は、次のとおりです。
exAcholnew_test = Chol Update( exAchol, exV, exC );
// Cholesky分解をより効率的に更新する
Show( exAcholnew_test );
// 結果は手動での手順の場合と同じ
Show( exAcholnew );
•
|
結果の行列Mでは、余計な0の対角要素は省かれています。
|
•
|
特異値分解によって、AはUSV’の形で表されます。ここで、
|
–
|
UとVは、互いに直交している列ベクトルを含んでいる行列です(「直交するベクトル」とは、「直角」もしくは「独立」なベクトルのことを指します)。
|
–
|
•
|
特異値分解は対応分析などで使われています。
|
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`;
Ortho()関数は、列の直交化を行い、さらにベクトルの大きさで割って正規化します。この関数はGram-Schmidt法を使用します。直交行列の列ベクトルは、大きさが単位長で互いに直交(内積が0)します。
B = Ortho( [1 -1, 1 0, 0 1] );
C = B` * B;
特に指定がない限り、ベクトルは正規化されます。つまり、ベクトルの大きさでベクトルそのものを割るスケーリングで、大きさ1のベクトルを作ります。Scaled(0)を指定して、スケーリングをしない場合は、次のようになります。
Ortho( [1 -1, 1 0, 0 1], Scaled( 0 ) );
全要素の和が0になるベクトルを作るには、オプションのCentered(1)を指定します。このオプションは対比の行列を作る場合に便利です。
result = Ortho( [1 -1, 1 0, 0 1], Centered( 1 ) );
J(1, 3) * result;
Ortho Poly()関数は、引数で指定された次元までのベクトルの直交多項式 を戻します。直交多項式は、多項式モデルを普通にあてはめると回帰係数の間に強い相関が生じてしまう場合に有用です。
Ortho Poly( [1 2 3], 2 );
QR()分解は、数値的に安定した行列の処理に便利です。QR()は2つの行列のリストを戻します。以下はその例です。
{Q, R} = QR( [11 22, 33 44] );
M`M行列の逆行列に行を追加または削除するには、Inv Update(S, X, w)関数を使用します。逆行列の更新は、1行除去による影響診断や候補計画の評価に役立ちます。
•
|
第1引数Sは更新する行列です。
|
•
|
第2引数Xは、追加または削除する行を含む行列です。
|
•
|
Inv Update(S, X, w)関数は、次の式を計算するのと同じです。
S - w * S * X` * Inv( I + w * X * S * X` ) * X * S