LU分解とは、ひとつの正方行列を、下三角行列L と上三角行列U の２つの積に分解するアルゴリズムです。
２つの行列に分解することにより、連立一次方程式の解や行列式の計算が容易になるメリットがあります。
例えば、3×3 行列であれば、次のような分解を行います。
$$
\begin{pmatrix}
a_{11} & a_{12} & a_{13} \\
a_{21} & a_{22} & a_{23} \\
a_{31} & a_{32} & a_{33}
\end{pmatrix}
=
\begin{pmatrix}
1 & 0 & 0 \\
l_{21} & 1 & 0 \\
l_{31} & l_{32} & 1
\end{pmatrix}
\begin{pmatrix}
u_{11} & u_{12} & u_{13} \\
0 & u_{22} & u_{23} \\
0 & 0 & u_{33}
\end{pmatrix}
$$
本記事では、LU分解の実装例を紹介します。
目次
Matrix クラス
まずはアルゴリズムで利用するための、Matrix クラスを掲載します。
以前のガウス消去法の記事で掲載したものをブラッシュアップしたものになります。
一部、本記事では利用しないメソッドも定義があります。
// matrix.h
#ifndef MATRIX_H_
#define MATRIX_H_
#include <cstddef>
#include <iostream>
#include <iomanip>
#include <vector>
class Matrix
{
public:
using Index = std::pair<std::size_t, std::size_t>;
Matrix(const Matrix &) = delete;
Matrix &operator=(const Matrix &) = delete;
Matrix(Matrix &&) = default;
Matrix &operator=(Matrix &&) = default;
Matrix clone() const { return Matrix(n_rows_, n_cols_, std::vector(v_)); }
std::size_t n_rows() const { return n_rows_; }
std::size_t n_cols() const { return n_cols_; }
double &operator[](Index i) { return v_[i.first * n_cols_ + i.second]; }
double operator[](Index i) const { return v_[i.first * n_cols_ + i.second]; }
void swap_elems(Index i1, Index i2);
void swap_rows(std::size_t i1, std::size_t i2);
void swap_columns(std::size_t i1, std::size_t i2);
void print() const;
static Matrix new_matrix(std::size_t n_rows, std::size_t n_cols, std::vector<double> &&v);
static Matrix new_square_matrix(std::size_t n, std::vector<double> &&v);
static Matrix new_zero_matrix(std::size_t n_rows, std::size_t n_cols);
static Matrix new_identity_matrix(std::size_t n);
private:
Matrix(std::size_t n_rows, std::size_t n_cols, std::vector<double> &&v)
: n_rows_(n_rows), n_cols_(n_cols), v_(std::move(v)) {}
std::vector<double> v_;
std::size_t n_rows_;
std::size_t n_cols_;
};
Matrix operator-(const Matrix &lhs, const Matrix &rhs);
Matrix operator*(const Matrix &lhs, const Matrix &rhs);
#endif // MATRIX_H_
// matrix.cc
#include "matrix.h"
// ２つの要素を交換する
void Matrix::swap_elems(Index i1, Index i2)
{
double tmp = (*this)[i1];
(*this)[i1] = (*this)[i2];
(*this)[i2] = tmp;
}
// ２つの行を交換する
void Matrix::swap_rows(std::size_t i1, std::size_t i2)
{
for (size_t k = 0; k < n_cols_; ++k)
{
swap_elems({i1, k}, {i2, k});
}
}
// ２つの列を交換する
void Matrix::swap_columns(std::size_t i1, std::size_t i2)
{
for (size_t k = 0; k < n_rows_; ++k)
{
swap_elems({k, i1}, {k, i2});
}
}
// 行列を標準出力する
void Matrix::print() const
{
std::cout << std::fixed;
std::cout << std::setprecision(2);
for (size_t i = 0; i < n_rows_; ++i)
{
if (i == 0)
{
std::cout << "[ ";
}
else
{
std::cout << " ";
}
for (size_t j = 0; j < n_cols_; ++j)
{
std::cout << std::setw(5) << (*this)[{i, j}] << " ";
}
if (i == n_rows_ - 1)
{
std::cout << "]\n";
}
else
{
std::cout << "\n";
}
}
}
// 任意の行列を生成
Matrix Matrix::new_matrix(std::size_t n_rows, std::size_t n_cols, std::vector<double> &&v)
{
v.resize(n_rows * n_cols);
return Matrix(n_rows, n_cols, std::move(v));
}
// 任意の正方行列を生成
Matrix Matrix::new_square_matrix(std::size_t n, std::vector<double> &&v)
{
v.resize(n * n);
return Matrix(n, n, std::move(v));
}
// 零行列を生成
Matrix Matrix::new_zero_matrix(std::size_t n_rows, std::size_t n_cols)
{
return Matrix(n_rows, n_cols, std::vector<double>(n_rows * n_cols));
}
// 単位行列を生成
Matrix Matrix::new_identity_matrix(std::size_t n)
{
auto a = new_square_matrix(n, {});
for (std::size_t i = 0; i < n; ++i)
{
a[{i, i}] = 1.0;
}
return a;
}
Matrix operator-(const Matrix &lhs, const Matrix &rhs)
{
std::size_t n_rows = lhs.n_rows();
std::size_t n_cols = rhs.n_cols();
auto out = Matrix::new_zero_matrix(n_rows, n_cols);
for (std::size_t i = 0; i < n_rows; ++i)
{
for (std::size_t j = 0; j < n_cols; ++j)
{
out[{i, j}] = lhs[{i, j}] - rhs[{i, j}];
}
}
return out;
}
Matrix operator*(const Matrix &lhs, const Matrix &rhs)
{
std::size_t n_rows = lhs.n_rows();
std::size_t n_cols = rhs.n_cols();
std::size_t n_inner = lhs.n_cols();
auto out = Matrix::new_zero_matrix(n_rows, n_cols);
for (std::size_t i = 0; i < n_rows; ++i)
{
for (std::size_t j = 0; j < n_cols; ++j)
{
for (std::size_t k = 0; k < n_inner; ++k)
{
out[{i, j}] += lhs[{i, k}] * rhs[{k, j}];
}
}
}
return out;
}
LU分解の実装例
LU分解アルゴリズムの実装例を掲載します。
Matrix クラスの定義さえできていれば、分解アルゴリズムそのものは比較的簡易なものとなります。
#include "matrix.h"
std::pair<Matrix, Matrix> lu_decomposition(const Matrix &a)
{
std::size_t n = a.n_rows();
auto lower_matrix = Matrix::new_identity_matrix(n);
auto upper_matrix = Matrix::new_zero_matrix(n, n);
for (std::size_t i = 0; i < n; ++i)
{
for (std::size_t j = 0; j < n; ++j)
{
auto x = a[{i, j}];
for (std::size_t k = 0; k < j; ++k)
{
x -= lower_matrix[{i, k}] * upper_matrix[{k, j}];
}
if (j < i)
{
x /= upper_matrix[{j, j}];
lower_matrix[{i, j}] = x;
}
else
{
upper_matrix[{i, j}] = x;
}
}
}
return {std::move(lower_matrix), std::move(upper_matrix)};
}LU分解の実行
3×3 の簡潔な行列
行列 \(
A=\begin{pmatrix}
1 & 1 & 1 \\
1 & 2 & 2 \\
2 & 3 & 2
\end{pmatrix}
\) に対して、LU分解を実行します。
２つの行列 \(
\begin{pmatrix}
1 & 0 & 0 \\
1 & 1 & 0 \\
2 & 1 & 1
\end{pmatrix}
\begin{pmatrix}
1 & 1 & 1 \\
0 & 1 & 1 \\
0 & 0 & -1
\end{pmatrix}
\) に分解されます。
また、\( A – LU \) の計算結果が零行列であることを出力することで、正しく分解が行われていることを確認します。
int main()
{
auto a = Matrix::new_square_matrix(3, {1, 1, 1, 1, 2, 2, 2, 3, 2});
std::cout << "入力行列：" << std::endl;
a.print();
auto [lower, upper] = lu_decomposition(a);
std::cout << "下三角行列L：" << std::endl;
lower.print();
std::cout << "上三角行列U：" << std::endl;
upper.print();
std::cout << "結果確認：" << std::endl;
(a - lower * upper).print();
return 0;
}入力行列：
[ 1.00 1.00 1.00
1.00 2.00 2.00
2.00 3.00 2.00 ]
下三角行列L：
[ 1.00 0.00 0.00
1.00 1.00 0.00
2.00 1.00 1.00 ]
上三角行列U：
[ 1.00 1.00 1.00
0.00 1.00 1.00
0.00 0.00 -1.00 ]
結果確認：
[ 0.00 0.00 0.00
0.00 0.00 0.00
0.00 0.00 0.00 ]5×5 の乱数行列
5×5 の乱数行列についても同様の確認を実施します。
\( A – LU \) の計算結果が零行列になることが確認できます。
#include <algorithm>
#include <random>
int main()
{
size_t n = 5;
std::vector<double> v(n * n);
std::random_device seed_gen;
std::mt19937 engine(seed_gen());
std::uniform_real_distribution<double> dist{-1.0, 1.0};
auto rng = [&engine, &dist]()
{
return dist(engine);
};
std::generate(v.begin(), v.end(), rng);
auto a = Matrix::new_square_matrix(n, std::move(v));
std::cout << "入力行列：" << std::endl;
a.print();
auto [lower, upper] = lu_decomposition(a);
std::cout << "下三角行列L：" << std::endl;
lower.print();
std::cout << "上三角行列U：" << std::endl;
upper.print();
std::cout << "結果確認：" << std::endl;
(a - lower * upper).print();
return 0;
}入力行列：
[ -0.57 0.68 0.14 -0.65 -0.85
0.89 -0.28 -0.79 -0.60 0.24
-0.38 0.06 0.64 0.91 0.25
-0.25 0.37 0.65 -0.02 0.76
0.74 -0.82 -0.14 0.22 -0.28 ]
下三角行列L：
[ 1.00 0.00 0.00 0.00 0.00
-1.57 1.00 0.00 0.00 0.00
0.67 -0.51 1.00 0.00 0.00
0.43 0.10 2.50 1.00 0.00
-1.31 0.09 0.38 0.78 1.00 ]
上三角行列U：
[ -0.57 0.68 0.14 -0.65 -0.85
0.00 0.78 -0.56 -1.63 -1.09
0.00 0.00 0.26 0.52 0.27
0.00 0.00 0.00 -0.87 0.56
0.00 0.00 0.00 0.00 -1.84 ]
結果確認：
[ 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 ]
コメント