行列のLU分解は、Pythonで簡単に実施することができます。
本記事では、Pythonを使用した行列のLU分解の方法について、詳しくご説明します。
- Python初心者の人
- Pythonを使用した行列のLU分解の方法について学びたい人
LU分解
LU分解（LU decomposition）とは、正方行列\(A\)を下三角行列\(L\)と上三角行列\(U\)の積に分解することを言います。
三角行列については、下記記事を参考にしてください。
任意の正方行列\(A\)に対し、三角行列を生成するために、置換行列\(P\)を使用した分解方法を、\(PA=LU\)分解と言います。
PA形式の場合、必ずLU分解できますので、使い勝手が良いです。
本記事では、上式両辺に左側から\(P^{-1}\)を掛け、\(P^{-1}\)を\(P\)と表現した、\(A=PLU\)で分解する方法をご紹介します。
scipy.linalg.lu()
SciPyのscipy.linalg.lu()を使用すると、引数に指定した正方行列をA=PLU分解することができます。
#input
import numpy as np
np.set_printoptions(precision=3)
from scipy.linalg import lu
np.random.seed(7)
A = np.random.randint(0, 10, (5, 5))
P, L, U = lu(A)
print("A:\n{}\n".format(A))
print("P:\n{}\n".format(P))
print("L:\n{}\n".format(L))
print("U:\n{}".format(U))#output
A:
[[4 9 6 3 3]
[7 7 9 7 8]
[9 8 7 6 4]
[0 7 0 7 6]
[3 5 8 8 7]]
P:
[[0. 0. 0. 1. 0.]
[0. 0. 0. 0. 1.]
[1. 0. 0. 0. 0.]
[0. 1. 0. 0. 0.]
[0. 0. 1. 0. 0.]]
L:
[[1. 0. 0. 0. 0. ]
[0. 1. 0. 0. 0. ]
[0.333 0.333 1. 0. 0. ]
[0.444 0.778 0.51 1. 0. ]
[0.778 0.111 0.627 0.107 1. ]]
U:
[[ 9. 8. 7. 6. 4. ]
[ 0. 7. 0. 7. 6. ]
[ 0. 0. 5.667 3.667 3.667]
[ 0. 0. 0. -6.98 -5.314]
[ 0. 0. 0. 0. 2.489]]式が成り立っているかどうか、確認してみます。
#input
PA = np.dot(P,A)
LU = np.dot(L,U)
print("PA:\n{}\n".format(PA))
print("LU:\n{}\n".format(LU))#output
A:
[[4 9 6 3 3]
[7 7 9 7 8]
[9 8 7 6 4]
[0 7 0 7 6]
[3 5 8 8 7]]
PLU:
[[4. 9. 6. 3. 3.]
[7. 7. 9. 7. 8.]
[9. 8. 7. 6. 4.]
[0. 7. 0. 7. 6.]
[3. 5. 8. 8. 7.]]scipy.linalg.lu_factor()
SciPyのscipy.linalg.lu_factor()を使用すると、引数に指定した正方行列をA=PLU分解し、\(L+U\)から\(L\)の主対角成分を除いた行列\(LU\)と置換行列\(P\)を表すピボットインデックスを出力することができます。
#input
import numpy as np
np.set_printoptions(precision=3)
from scipy.linalg import lu_factor
np.random.seed(8)
A = np.random.randint(0, 10, (5, 5))
LU, piv = lu_factor(A)
print("A:\n{}\n".format(A))
print("LU:\n{}".format(LU))
print("piv:\n{}".format(piv))#output
A:
[[3 4 1 9 5]
[8 3 8 0 5]
[1 3 9 2 2]
[6 8 9 3 4]
[5 5 7 9 2]]
LU:
[[ 8. 3. 8. 0. 5. ]
[ 0.75 5.75 3. 3. 0.25 ]
[ 0.125 0.457 6.63 0.63 1.261]
[ 0.375 0.5 -0.528 7.833 3.666]
[ 0.625 0.543 0.056 0.936 -4.763]]
piv:
[1 3 2 3 4]scipy.linalg.lu_solve()
SciPyのscipy.linalg.lu_solve()を使用すると、scipy.linalg.lu_factor()の戻り値A=PLUを使用して、連立方程式\(Ax=b\)を解くことができます。
#input
import numpy as np
np.set_printoptions(precision=3)
from scipy.linalg import lu_factor, lu_solve
np.random.seed(9)
A = np.random.randint(0, 10, (5, 5))
b = np.random.randint(0, 10, 5)
x = lu_solve(lu_factor(A),b)
print("A:\n{}\n".format(A))
print("b:\n{}".format(b))
print("x:\n{}".format(x))#output
A:
[[5 6 8 6 1]
[6 4 8 1 8]
[5 1 0 8 8]
[8 2 6 8 1]
[8 3 5 3 6]]
b:
[7 9 0 8 1]
x:
[-1.326 -2.268 2.797 0.749 0.362]出力が正しいか、確かめてみます。
#input
print("Ax:\n{}\n".format(np.dot(A,x)))
print("b:\n{}".format(b))#output
Ax:
[7.000e+00 9.000e+00 4.441e-16 8.000e+00 1.000e+00]
b:
[7 9 0 8 1]\(Ax=b\)が成り立っていることが確認できました。
まとめ
この記事では、Pythonを使用した、行列のLU分解の方法について、ご説明しました。
本記事を参考に、ぜひ試してみて下さい。
参考
Python学習用おすすめ教材
Pythonの基本を学びたい方向け
統計学基礎を学びたい方向け
Pythonの統計解析を学びたい方向け
おすすめプログラミングスクール
Pythonをはじめ、プログラミングを学ぶなら、TechAcademy（テックアカデミー）がおすすめです。
私も入っていますが、好きな時間に気軽にオンラインで学べますので、何より楽しいです。
現役エンジニアからマンツーマンで学べるので、一人では中々続かない人にも、向いていると思います。
無料体験ができますので、まずは試してみてください！
