C++で任意の行列のLU分解を行うプログラムの実装方法

行列のLU分解とは、ある行列を下三角行列（L）と上三角行列（U）の積として表現する手法です。名前の「LU」は「Lower（下三角）」「Upper（上三角）」の頭文字に由来しています。

LU分解は、連立一次方程式の求解や逆行列の計算など、数値計算の分野で広く活用されている重要なアルゴリズムです。

LU分解の例

以下に、3×3の行列をLU分解した具体例を示します。

元の行列 A：
1 1 0
2 1 3
3 1 1

下三角行列 L：
1 0 0
2 -1 0
3 -2 -5

上三角行列 U：
1 1 0
0 1 -3
0 0 1

このとき、L × U = A が成立します。この分解を実際に行うC++プログラムは以下の通りです。

C++プログラムの全体コード

#include<iostream>
using namespace std;

void LUdecomposition(float a[10][10], float l[10][10], float u[10][10], int n) {
    int i = 0, j = 0, k = 0;
    for (i = 0; i < n; i++) {
        // 下三角行列 L の計算
        for (j = 0; j < n; j++) {
            if (j < i)
                l[j][i] = 0;
            else {
                l[j][i] = a[j][i];
                for (k = 0; k < i; k++) {
                    l[j][i] = l[j][i] - l[j][k] * u[k][i];
                }
            }
        }
        // 上三角行列 U の計算
        for (j = 0; j < n; j++) {
            if (j < i)
                u[i][j] = 0;
            else if (j == i)
                u[i][j] = 1;
            else {
                u[i][j] = a[i][j] / l[i][i];
                for (k = 0; k < i; k++) {
                    u[i][j] = u[i][j] - ((l[i][k] * u[k][j]) / l[i][i]);
                }
            }
        }
    }
}

int main() {
    float a[10][10], l[10][10], u[10][10];
    int n = 0, i = 0, j = 0;

    cout << "Enter size of square matrix : "<<endl;
    cin >> n;

    cout<<"Enter matrix values: "<endl;
    for (i = 0; i < n; i++)
        for (j = 0; j < n; j++)
            cin >> a[i][j];

    LUdecomposition(a, l, u, n);

    cout << "L Decomposition is as follows..."<<endl;
    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            cout<<l[i][j]<<" ";
        }
        cout << endl;
    }

    cout << "U Decomposition is as follows..."<<endl;
    for (i = 0; i < n; i++) {
        for (j = 0; j < n; j++) {
            cout<<u[i][j]<<" ";
        }
        cout << endl;
    }
    return 0;
}

実行結果

上記プログラムを実行すると、以下のような出力が得られます。

Enter size of square matrix : 3
Enter matrix values:
1 1 0
2 1 3
3 1 1
L Decomposition is as follows...
1 0 0
2 -1 0
3 -2 -5
U Decomposition is as follows...
1 1 0
0 1 -3
0 0 1

プログラムの解説

LUdecomposition関数の処理内容

上記プログラムでは、LUdecomposition関数が与えられた行列に対してL行列とU行列を求めます。具体的には、ネストされたforループを用いて各要素を順番に計算し、その結果を配列 l[][] と u[][] に格納していきます。

下三角行列Lを計算する部分のコードは以下の通りです。

for (i = 0; i < n; i++) {
    for (j = 0; j < n; j++) {
        if (j < i)
            l[j][i] = 0;
        else {
            l[j][i] = a[j][i];
            for (k = 0; k < i; k++) {
                l[j][i] = l[j][i] - l[j][k] * u[k][i];
            }
        }
    }

対角成分より下側（j < i）は0を設定し、それ以外の要素については、元の行列の値からすでに計算済みのL・Uの積を差し引くことで求めています。

続いて、上三角行列Uを計算する部分です。

    for (j = 0; j < n; j++) {
        if (j < i)
            u[i][j] = 0;
        else if (j == i)
            u[i][j] = 1;
        else {
            u[i][j] = a[i][j] / l[i][i];
            for (k = 0; k < i; k++) {
                u[i][j] = u[i][j] - ((l[i][k] * u[k][j]) / l[i][i]);
            }
        }
    }
}

こちらでは、対角成分より上側（j > i）の要素を計算し、対角成分には1を設定します。これにより、対角成分がすべて1となる「単位上三角行列」としてUが構成されます（ドーリトル法と呼ばれる手法です）。

main関数の処理内容

main()関数では、まずユーザーから行列のサイズと各要素の値を入力として受け取ります。

cout << "Enter size of square matrix : "<<endl;
cin >> n;
cout<<"Enter matrix values: "<endl;
for (i = 0; i < n; i++)
    for (j = 0; j < n; j++)
        cin >> a[i][j];

その後、LU分解関数を呼び出し、計算結果であるL行列とU行列を画面に表示します。

LUdecomposition(a, l, u, n);
cout << "L Decomposition is as follows..."<<endl;
for (i = 0; i < n; i++) {
    for (j = 0; j < n; j++) {
        cout<<l[i][j]<<" ";
    }
    cout << endl;
}
cout << "U Decomposition is as follows..."<<endl;
for (i = 0; i < n; i++) {
    for (j = 0; j < n; j++) {
        cout<<u[i][j]<<" ";
    }
    cout << endl;
}

注意点

このアルゴリズムでは、対角成分 l[i][i] で除算を行うため、ピボット（対角要素）が0に近い場合には数値的な不安定さが生じる可能性があります。実際の数値計算ライブラリでは、行の入れ替え（ピボット選択）を組み合わせた「ピボット付きLU分解」が一般的に使用されます。