C#
行列分解
行列分解とは、数値の正方行列を 2 つの正方行列に分割する手法を指し、連立方程式を効率的に解くための基本です。つまり、逆行列を計算する (行列を反転する) ための基本でもあります。この逆行列の計算は、多くの重要なアルゴリズムに含まれています。今回は、行列分解、逆行列の計算、連立方程式の解法、および関連する演算を実行する C# コードを具体的に示しながら説明します。
確かに行列分解は華々しい話題ではありませんが、一連の行列メソッドを開発者のコード ライブラリに備えておけば大きな効果を得られます。今回は、ニーズに合わせてソース コードを編集できるよう、これらのメソッドについて詳しく説明していきます。また、行列メソッドで使用している手法の一部は、他のコーディング シナリオにも再利用できます。
今回説明する内容の概要については、図 1 のスクリーンショットをご覧いただければ一目瞭然です。デモ プログラムでは、まず 4 x 4 正方行列を作成し、その値を表示します。次に、この行列を LUP 行列という行列に分解します。L は下三角行列、U は上三角行列を表します。P 部分 (P は置換を表します) は、値が {3,1,2,0} の配列で、分解処理中に行 0 と行 3 を交換したことを表します。分解では、トグル値 -1 も生成しています。この値は、交換した行数が奇数であることを表します。デモ プログラムは、分解を 2 つの方法で表示します。最初に結合した LU 行列として表示し、次に個別の L 行列と U 行列として表示します。続いて、元の行列の逆行列を計算し、表示します。この処理の内部では、LUP 行列を使用しています。デモ プログラムでは、再び分解を使用して、元の行列の行列式を計算します。その後、逆行列を使用して連立方程式を解き、最後に L 行列と U 行列を結合して元の行列に戻します。
.png)
図 1 行列分解のデモ
なぜ独自の行列分解メソッドと関連メソッドのライブラリを作成するのでしょう。スタンドアロンの行列ツールは多数存在しますが、他のアプリケーションやシステムとの統合が難しい場合があります。また、行列分解は根本的に重要な技法ですが、著作権で保護されていない無料の .NET コード実装はほとんど存在しません。存在していても、多くの場合は実際のコーディング シナリオに合わせて変更できるほど解説が充実していません。
ここでは、中級レベルの C# プログラミング スキルがあり、行列演算とその用語について少なくとも基礎知識があることを前提としています。主な C# コードはすべてここに掲載しています。また、MSDN のコード ダウンロード サイト (arichve.msdn.microsoft.com/mag201212matrix、英語) でもコードをダウンロードできます。
行列の定義
C# で行列を実装する方法はいくつかあります。従来の方法 (ここで使用する方法) では、ジャグ配列とも呼ばれる、配列の配列を使用します。たとえば、次のコードでは、3 行 2 列の行列を定義します。
double[][] m = new double[3][];
m[0] = new double[2];
m[1] = new double[2];
m[2] = new double[2];
m[2][1] = 5.0; // set row 2, col 1
ほとんどのプログラミング言語とは異なり、C# には組み込みの多次元配列型があるので、別の方法も利用できます。次に例を示します。
double[,] m = new double[3,2];
m[2,1] = 5.0;
C# で行列を実装する 3 つ目の方法は、次のように、配列インデックスと組み合わせた 1 つの配列を使用します。
int rows = 3;
int cols = 2;
double[] m = new double[rows * cols];
int i = 2;
int j = 1;
m[i * cols + j] = 5.0;
使用する格納方法に関係なく、行列の実装には OOP 手法も静的メソッド手法も使用できます。たとえば、OOP 手法では次のようになります。
public class MyMatrix
{
int m; // number rows
int n; // number columns
data[][]; // the values
...
}
行列の設計には、最適な設計が 1 つだけということはありません。最適な設計は、コードを運用する特定のコーディング シナリオと、コーディングに関する個人的な好みによって異なります。ここでは、最もわかりやすくリファクタリングが簡単なことから、静的メソッド手法を使用します。
行列に "配列の配列" 設計を使用する場合は各行を個別に割り当てる必要があるので、メモリ割り当てを実行するヘルパー メソッドを定義すると役立つことがよくあります。次に例を示します。
static double[][] MatrixCreate(int rows, int cols)
{
// creates a matrix initialized to all 0.0s
// do error checking here?
double[][] result = new double[rows][];
for (int i = 0; i < rows; ++i)
result[i] = new double[cols]; // auto init to 0.0
return result;
}
このメソッドは次のように呼び出せます。
double[][] m = MatrixCreate(3,2);
m[2][1] = 5.0;
このメソッドは、独自の行列メソッド ライブラリを作成する利点の 1 つを示しています。パフォーマンスを向上する場合は、入力パラメーターのエラー チェックを省略できますが、コードの呼び出しで例外が発生するリスクが高くなります (ここでは、説明を簡単にするためにほとんどのエラー チェックを取り除いています)。もう 1 つの利点は、ライブラリをカスタマイズして、実際のシナリオ向けに最適化できることです。独自のライブラリ作成の主な欠点は、既存のライブラリを使用するよりも時間がかかる可能性があることです。
カスタム行列ライブラリに追加する別の便利なメソッドは、行列を文字列として表示するメソッドです。これは、次のようになります。
static string MatrixAsString(double[][] matrix)
{
string s = "";
for (int i = 0; i < matrix.Length; ++i)
{
for (int j = 0; j < matrix[i].Length; ++j)
s += matrix[i][j].ToString("F3").PadLeft(8) + " ";
s += Environment.NewLine;
}
return s;
}
表示する桁数、列幅のパディング、またはその両方をパラメーターにすることをお勧めします。
行列乗算
Microsoft .NET Framework 4 以降には、行列乗算メソッドのパフォーマンスを大幅に向上できる便利な機能が用意されています。図 2 に、行列乗算を示します。
.png)
図 2 行列乗算
結果の各セル値の計算は結果に含まれる他のセル値に依存していないので、各計算は独立しており、複数のプロセッサを搭載したコンピューターなら並列に実行できます。行列乗算の標準手法は次のとおりです。
static double[][] MatrixProduct(double[][] matrixA,
double[][] matrixB)
{
int aRows = matrixA.Length; int aCols = matrixA[0].Length;
int bRows = matrixB.Length; int bCols = matrixB[0].Length;
if (aCols != bRows)
throw new Exception("Non-conformable matrices in MatrixProduct");
double[][] result = MatrixCreate(aRows, bCols);
for (int i = 0; i < aRows; ++i) // each row of A
for (int j = 0; j < bCols; ++j) // each col of B
for (int k = 0; k < aCols; ++k)
result[i][j] += matrixA[i][k] * matrixB[k][j];
return result;
}
行列乗算は O(n^3) 演算なので、パフォーマンスが問題になることがあります。たとえば、行列 A のサイズが 300 x 200 で、行列 B のサイズが 200 x 400 の場合、A と B の積を計算するには、300 * 200 * 400 = 24,000,000 回の乗算が必要です。.NET Framework 4 以降の System.Threading.Tasks 名前空間に含まれるタスク並列ライブラリ (TPL) を使用すると、単純な並列処理バージョンの行列乗算のコードを簡単に作成できます。たとえば、次のようになります。
static double[][] MatrixProduct(double[][] matrixA,
double[][] matrixB)
{
// error check, compute aRows, aCols, bCols
double[][] result = MatrixCreate(aRows, bCols);
Parallel.For(0, aRows, i =>
{
for (int j = 0; j < bCols; ++j)
for (int k = 0; k < aCols; ++k)
result[i][j] += matrixA[i][k] * matrixB[k][j];
}
);
return result;
}
このバージョンでは、計算を行単位に分割しています。内部の処理では、TPL によって、複数のプロセッサで計算を実行する複雑な同期プラミング コードがすべて生成されます。
整合性テスト
相互に関連するメソッドのライブラリに関する興味深い点の 1 つは、ライブラリ内のメソッドが整合性のある結果を生成するかどうか確認することで、多くの場合はメソッドをテストできる点です。たとえば、ランダムな行列を作成する次のようなメソッドがあるとします。
static double[][] MatrixRandom(int rows, int cols,
double minVal, double maxVal, int seed)
{
// return matrix with values between minVal and maxVal
Random ran = new Random(seed);
double[][] result = MatrixCreate(rows, cols);
for (int i = 0; i < rows; ++i)
for (int j = 0; j < cols; ++j)
result[i][j] = (maxVal - minVal) * ran.NextDouble() + minVal;
return result;
}
また、単位行列を作成する次のようなメソッドがあるとします。単位行列とは、主対角線上の値が 1.0 で残りの値が 0.0 の正方行列です。
static double[][] MatrixIdentity(int n)
{
double[][] result = MatrixCreate(n, n);
for (int i = 0; i < n; ++i)
result[i][i] = 1.0;
return result;
}
さらに、2 つの行列が等しいかどうか比較する次のようなメソッドがあるとします。
static bool MatrixAreEqual(double[][] matrixA,
double[][] matrixB, double epsilon)
{
// true if all values in A == corresponding values in B
int aRows = matrixA.Length;
int bCols = matrixB[0].Length;
for (int i = 0; i < aRows; ++i) // each row of A and B
for (int j = 0; j < aCols; ++j) // each col of A and B
if (Math.Abs(matrixA[i][j] - matrixB[i][j]) > epsilon)
return false;
return true;
}
値が double 型なので、MatrixAreEqual メソッドではセル値が正確に等しいかどうかを確認していないことに注意してください。代わりに、すべてのセル値どうしの差が非常に小さい (epsilon の値以下) かどうかを確認します。
正方行列 m と同じサイズの単位行列との積は、元の行列 m に等しいので、次のコードを使用して行列乗算メソッドをテストできます。
double[][] m = MatrixRandom(4, 4, -9.0, 9.0, 0);
double[][] i = MatrixIdentity(4);
double[][] mi = MatrixProduct(m, i);
if (MatrixAreEqual(m, mi, 0.00000001) == true)
Console.WriteLine("Consistent result");
else
Console.WriteLine("Something is wrong");
Consistency checking lends itself well to random input testing.
行列分解
行列分解では、正方行列 M を受け取り、掛け合わせると元の行列 M になる 2 つの新しい正方行列を計算します。この考え方は、通常の数の因数分解 (2 * 3 = 6 が成立するので、数 6 は因数 2 と 3 に分解できること) に似ています。行列分解は一見してほとんど意味がないように思えますが、実は行列分解を使用すると、非常に難しい逆行列の計算がかなり簡単になります。行列分解にはさまざまな種類があり、しかも各種類の分解計算に複数のアルゴリズムがあります。ここで紹介する手法は LUP 分解という手法で、部分ピボット選択を併用したドゥーリトル法を使用します。
LUP 分解を理解するには、比較的単純な LU 分解を理解しておくと役に立ちます。LU 分解は、著名な数学者アラン チューリングが発案した手法です。次のような、M という 4 x 4 行列があるとします。
9.000 5.000 3.000 4.000
4.000 8.000 2.000 5.000
3.000 5.000 7.000 1.000
2.000 6.000 0.000 8.000
M の LU 分解の一例として、L は次のようになります。
1.000 0.000 0.000 0.000
0.444 1.000 0.000 0.000
0.333 0.577 1.000 0.000
0.222 0.846 -0.219 1.000
この場合、U は次のとおりです。
9.000 5.000 3.000 4.000
0.000 5.778 0.667 3.222
0.000 0.000 5.615 -2.192
0.000 0.000 0.000 3.904
L * U = M となるので、この分解は成立します。下三角行列 L は、対角線上の値が 1.0 で、それより右上の値が 0.0 であることに注目してください。つまり、下三角行列で意味のあるセル値は、左下に存在します。同様に、上三角行列で意味のあるセル値は、主対角線上とその右上に存在します。
また、L と U の意味のあるセル値の場所がまったく重なっていないことにも注目してください。このため、通常の行列分解では、2 つの結果の行列 (L と U) を生成する代わりに、L と U の両方を保持する 1 つの行列に下三角行列と上三角行列の両方の結果を格納して、メモリ領域を節約します。
LUP 行列分解は、LU 分解とほぼ同じですが、重要な違いがあります。LUP 分解は、行列 M を受け取り、L 行列と U 行列だけでなく、P 配列も生成します。LUP における L と U の積は、元の行列 M そのものではなく、M の行を一部並べ替えた行列になります。行の並べ替え方は P 配列に格納されているため、この情報を使用して元の行列 M を再構築できます。
ここで紹介するドゥーリトル分解によく似た手法に、クラウトの分解という手法もあります。ドゥーリトル法とクラウト法の主な違いは、ドゥーリトル法では L 行列の主対角線上に 1.0 を配置するのに対し、クラウト法では U 行列の主対角線上に 1.0 を配置します。
LUP 分解が LU 分解よりも頻繁に使用されるのは些細な理由からです。後ほど説明するように、行列分解は、逆行列の計算に使用します。しかし、行列分解を逆行列計算のヘルパーとして使用する場合、LU 行列の主対角線上に値 0.0 が存在すると計算に失敗します。そのため、LUP 分解では、主対角線上に値 0.0 が発生する場合は 2 行を交換して値 0.0 を対角線上から除去し、交換した行を P 配列で追跡します。
図 3 に、行列分解メソッドを示します。
図 3 行列分解メソッド
static double[][] MatrixDecompose(double[][] matrix,
out int[] perm, out int toggle)
{
// Doolittle LUP decomposition.
// assumes matrix is square.
int n = matrix.Length; // convenience
double[][] result = MatrixDuplicate(matrix);
perm = new int[n];
for (int i = 0; i < n; ++i) { perm[i] = i; }
toggle = 1;
for (int j = 0; j < n - 1; ++j) // each column
{
double colMax = Math.Abs(result[j][j]); // largest val in col j
int pRow = j;
for (int i = j + 1; i < n; ++i)
{
if (result[i][j] > colMax)
{
colMax = result[i][j];
pRow = i;
}
}
if (pRow != j) // swap rows
{
double[] rowPtr = result[pRow];
result[pRow] = result[j];
result[j] = rowPtr;
int tmp = perm[pRow]; // and swap perm info
perm[pRow] = perm[j];
perm[j] = tmp;
toggle = -toggle; // row-swap toggle
}
if (Math.Abs(result[j][j]) < 1.0E-20)
return null; // consider a throw
for (int i = j + 1; i < n; ++i)
{
result[i][j] /= result[j][j];
for (int k = j + 1; k < n; ++k)
result[i][k] -= result[i][j] * result[j][k];
}
} // main j column loop
return result;
}
このメソッドは次のように呼び出すことができます。
double[][] m = MatrixRandom(4, 4, -9.0, 9.0, 0);
int[] perm;
int toggle;
double luMatrix = MatrixDecompose(m, out perm, out toggle);
MatrixDecompose メソッドは、入力として正方行列を受け取ります。メソッドの戻り値は 3 つあります。明示的な戻り値は、置換した LU 行列です。メソッドでは、出力パラメーターとして 2 つの値を返します。1 つ目の出力パラメーターは、行の置換方法に関する情報を保持する置換配列です。もう 1 つは、+1 または -1 のトグル値で、この値は交換した行数が偶数 (+1) か奇数 (-1) かに応じて変化します。トグル値は、逆行列計算には使用しませんが、行列分解を使用して行列式を計算する場合に必要です。
MatrixDecompose メソッドはかなり複雑ですが、コードを変更するために理解しておく必要がある細部はわずかです。次に示すバージョンでは、ヘルパー メソッド MatrixDuplicate を使用して、戻り値となる LU 行列の新しいメモリを割り当てます。
static double[][] MatrixDuplicate(double[][] matrix)
{
// assumes matrix is not null.
double[][] result = MatrixCreate(matrix.Length, matrix[0].Length);
for (int i = 0; i < matrix.Length; ++i) // copy the values
for (int j = 0; j < matrix[i].Length; ++j)
result[i][j] = matrix[i][j];
return result;
}
他には、メモリを節約するために、計算結果を入力行列に格納する方法もあります。C# のセマンティクスでは、この方法は matrix パラメーターを入力と出力の両方に使用するので、パラメーターの種類を ref パラメーターに変更することになります。この方法を使用すると、メソッド シグネチャは次のようになります。
static void MatrixDecompose(ref double[][] matrix, out int[] perm,
out int toggle)
また、明示的な戻り値がなくなったので、このメソッドを置換配列や交換のトグルにも使用できます。次に例を示します。
static int[] MatrixDecompose(ref double[][] matrix, out int toggle)
行列分解を行列式の計算に使用しない場合は、toggle パラメーターを削除してメソッド シグネチャを簡略化することをお勧めします。
変更をお勧めする MatrixDecompose の別の部分は、次のステートメントです。
if (Math.Abs(result[j][j]) < 1.0E-20)
return null;
文章にすると、このコードは、「主対角線上に値 0.0 が存在したために 2 つの行を入れ替えた後でも、主対角線上にまだ値 0.0 が存在する場合は、null を返す」ことになります。ゼロ チェックの等式に使用する任意の浮動小数点数を 1.0E-20 から変更して、システムの精度特性に基づいた値にすることをお勧めします。また、null を返す代わりに、例外をスローすることをお勧めします。これはこのメソッドを逆行列計算の一環として呼び出す場合、計算に失敗するためです。最後に、逆行列計算以外の目的で行列分解を使用する場合は、このステートメント全体を削除することをお勧めします。
逆行列計算
行列分解を使用して逆行列を計算する際に重要なのは、連立方程式を解くヘルパー メソッドを作成することです。この重要なヘルパー メソッドを図 4 に示します。
図 4 HelperSolve メソッド
static double[] HelperSolve(double[][] luMatrix,
double[] b)
{
// solve luMatrix * x = b
int n = luMatrix.Length;
double[] x = new double[n];
b.CopyTo(x, 0);
for (int i = 1; i < n; ++i)
{
double sum = x[i];
for (int j = 0; j < i; ++j)
sum -= luMatrix[i][j] * x[j];
x[i] = sum;
}
x[n - 1] /= luMatrix[n - 1][n - 1];
for (int i = n - 2; i >= 0; --i)
{
double sum = x[i];
for (int j = i + 1; j < n; ++j)
sum -= luMatrix[i][j] * x[j];
x[i] = sum / luMatrix[i][i];
}
return x;
}
HelperSolve メソッドは、LU 行列に乗算すると配列 b になるような、配列 x を見つけます。このメソッドは非常に高度なので、完全に理解するには、いくつかの例を追跡する必要があります。メソッドには 2 つのループがあります。最初のループでは、LU 行列内の下三角行列に対する前進代入を使用します。次のループでは、LU 行列内の上三角行列に対する後退代入を使用します。同様のメソッドは、一部の行列分解の実装で luBackSub などと呼ばれています。
このコードは短くても複雑ですが、このコードを変更する必要に迫られる状況はないはずです。HelperSolve メソッドは MatrixDecompose メソッドから LU 行列を受け取りますが、perm 出力パラメーターは使用しないことに注意してください。これは、HelperSolve メソッドが実際にヘルパー メソッドであり、連立方程式を解くには追加のラッパー コードが必要なことを表しています。ここで紹介するコードを OOP 設計にリファクタリングする場合は、HelperSolve メソッドをプライベート メソッドにすることをお勧めします。
HelperSolve メソッドの準備が整ったので、逆行列計算メソッドを実装できます (図 5 参照)。
図 5 MatrixInverse メソッド
static double[][] MatrixInverse(double[][] matrix)
{
int n = matrix.Length;
double[][] result = MatrixDuplicate(matrix);
int[] perm;
int toggle;
double[][] lum = MatrixDecompose(matrix, out perm, out toggle);
if (lum == null)
throw new Exception("Unable to compute inverse");
double[] b = new double[n];
for (int i = 0; i < n; ++i)
{
for (int j = 0; j < n; ++j)
{
if (i == perm[j])
b[j] = 1.0;
else
b[j] = 0.0;
}
double[] x = HelperSolve(lum, b);
for (int j = 0; j < n; ++j)
result[j][i] = x[j];
}
return result;
}
このコードも複雑です。逆行列を計算するアルゴリズムは、行列 M とその逆行列の積が単位行列になるという考え方を基にしています。MatrixInverse メソッドは、基本的に連立方程式 Ax = b を解きます。この連立方程式では、A は LU 行列分解です。また、定数の集合 b は、それぞれの値が 1.0 または 0.0 で、単位行列に対応しています。MatrixInverse では、MatrixDecompose の呼び出しで取得した perm 配列を使用していることに注意してください。
MatrixInverse メソッドは次のように呼び出すことができます。
double[][] m = MatrixRandom(4, 4, -9.0, 9.0, 0);
double[][] inverse = MatrixInverse(m);
Console.WriteLine(MatrixAsString(inverse));
まとめると、重要な行列演算の 1 つが逆行列計算ですが、これは非常に難しい演算です。逆行列計算手法の 1 つでは、まず元の行列を分解します。次に、前進代入と後退代入を実行するヘルパー解法メソッドを作成します。最後に、分解を LU 置換配列とヘルパー解法メソッドと併用して、逆行列を特定します。この手法は複雑に思えますが、通常は逆行列を直接計算するよりも効率的で簡単です。
行列式
行列分解メソッドがあると、行列式を計算する次のようなメソッドを簡単に作成できます。
static double MatrixDeterminant(double[][] matrix)
{
int[] perm;
int toggle;
double[][] lum = MatrixDecompose(matrix, out perm, out toggle);
if (lum == null)
throw new Exception("Unable to compute MatrixDeterminant");
double result = toggle;
for (int i = 0; i < lum.Length; ++i)
result *= lum[i][i];
return result;
}
実は、少し意外なことに、行列式とは行列の LUP 分解で主対角線上にある値の積を (トグル値に応じて) 正または負にしたものです。トグルの値を int から double に暗黙に変換していることに注意してください。MatrixDeterminant にエラー チェックを追加する以外に、入力行列のサイズが 1 (戻り値が単一の値) または 2 x 2 (戻り値が a*d - b*c) の場合を対象に、省略された戻り値を追加することをお勧めします。行列式メソッドは次のように呼び出すことができます。
double[][] m = MatrixRandom(4, 4, -9.0, 9.0, 0);
double det = MatrixDeterminant(m);
Console.WriteLine("Determinant = " + det.ToString("F2"));
連立方程式を解く
HelperSolve メソッドは、次のような連立一次方程式の解法に簡単に応用できます。
static double[] SystemSolve(double[][] A, double[] b)
{
// Solve Ax = b
int n = A.Length;
int[] perm;
int toggle;
double[][] luMatrix = MatrixDecompose(A,
out perm, out toggle);
if (luMatrix == null)
return null; // or throw
double[] bp = new double[b.Length];
for (int i = 0; i < n; ++i)
bp[i] = b[perm[i]];
double[] x = HelperSolve(luMatrix, bp);
return x;
}
図 1 のスクリーンショットを生成した、次の連立方程式を解くコードを示します。
3x0 + 7x1 + 2x2 + 5x3 = 49
x0 + 8x1 + 4x2 + 2x3 = 30
2x0 + x1 + 9x2 + 3x3 = 43
5x0 + 4x1 + 7x2 + x3 = 52
double[][] m = MatrixCreate(4, 4);
m[0][0] = 3.0; m[0][1] = 7.0; m[0][2] = 2.0; m[0][3] = 5.0;
m[1][0] = 1.0; m[1][1] = 8.0; m[1][2] = 4.0; m[1][3] = 2.0;
m[2][0] = 2.0; m[2][1] = 1.0; m[2][2] = 9.0; m[2][3] = 3.0;
m[3][0] = 5.0; m[3][1] = 4.0; m[3][2] = 7.0; m[3][3] = 1.0;
double[] b = new double[] { 49.0, 30.0, 43.0, 52.0 };
double[] x = SystemSolve(m, b);
Console.WriteLine("\nSolution is x = \n" + VectorAsString(x));
SystemSolve メソッドでは、MatrixDecompose から取得した perm 配列を使用して b 入力パラメーターを並べ替えてから、HelperSolve を呼び出していることに注意してください。
置換配列とは
図 1 のスクリーンショットの最後の数行では、元の行列を取得可能な方法で L 行列と U 行列を乗算できることを示しています。この計算方法について知っても実用的な行列問題を解けるようにはなりませんが、LUP 分解の P 部分を理解する手助けとなります。また、L 成分と U 成分から元の行列を再生成すると、行列ライブラリ メソッドの整合性をテストする場合に便利です。
LUP 分解の実行後に元の行列を再生成する方法の 1 つでは、次のように、L と U を掛け合わせてから、P 配列を使用して積の行を置換します。
// create matrix m
// call MatrixDecompose, returning lu and perm
// extract the lower part of lu as matrix 'lower'
// extract the upper part of lu as matrix 'upper'
double[][] lu = MatrixProduct(lower, upper);
double[][] orig = UnPermute(lu, perm);
if (MatrixAreEqual(orig, m, 0.000000001) == true)
Console.WriteLine("L and U unpermuted using perm array");
UnPermute メソッドは次のように記述できます。
static double[][] UnPermute(double[][] luProduct, int[] perm)
{
double[][] result = MatrixDuplicate(luProduct);
int[] unperm = new int[perm.Length];
for (int i = 0; i < perm.Length; ++i)
unperm[perm[i]] = i; // create un-perm array
for (int r = 0; r < luProduct.Length; ++r) // each row
result[r] = luProduct[unperm[r]];
return result;
}
LUP 分解から元の行列を再生成する 2 つ目の方法では、次のように、perm 配列を perm 行列に変換し、結合した LU 行列を perm 行列に乗算します。
// create matrix m
// call MatrixDecompose, returning lu and perm
// extract the lower part of lu as matrix 'lower'
// extract the upper part of lu as matrix 'upper'
double[][] permMatrix = PermArrayToMatrix(perm);
double[][] orig = MatrixProduct(permMatrix, lu);
if (MatrixAreEqual(orig, m, 0.000000001) == true)
Console.WriteLine("L and U unpermuted using perm matrix");
perm 行列は、各行と各列に値 1.0 が 1 回ずつ出現する正方行列です。perm 配列から perm 行列を作成するメソッドは次のように記述できます。
static double[][] PermArrayToMatrix(int[] perm)
{
// Doolittle perm array to corresponding matrix
int n = perm.Length;
double[][] result = MatrixCreate(n, n);
for (int i = 0; i < n; ++i)
result[i][perm[i]] = 1.0;
return result;
}
まとめ
連立一次方程式の解法、逆行列の特定、または行列式の特定を必要とするアルゴリズムは多数あります。行列分解の使用は、このような演算すべての実行に効果的な手法です。ここで紹介したコードは、コード ベースに外部依存関係が存在することを望まない場合や、パフォーマンス向上や機能変更を目的として演算をカスタマイズする必要がある場合に使用できます。行列の真実に目覚めましょう。
Dr. James McCaffrey は Volt Information Sciences Inc. に勤務し、ワシントン州レドモンドにあるマイクロソフト本社で働くソフトウェア エンジニアの技術トレーニングを管理しています。これまでに、Internet Explorer、MSN サーチなどの複数のマイクロソフト製品にも携わってきました。また、『.NET Test Automation Recipes: A Problem-Solution Approach』(Apress、2006 年) の著者でもあります。連絡先は jammc@microsoft.com (英語のみ) です。
この記事のレビューに協力してくれた技術スタッフの Paul Koch と Dan Liebling に心より感謝いたします。