【Python入門】行列のLU分解|scipy.linalg.lu()を学ぶ

行列のLU分解|scipy.linalg.lu()を学ぶ_アイキャッチ プログラミング

※ 当サイトはアフィリエイト広告を利用しています。

levtech-ad
スポンサーリンク

行列のLU分解は、Pythonで簡単に実施することができます。

本記事では、Pythonを使用した行列のLU分解の方法について、詳しくご説明します。

こんな人に読んでほしい
  • Python初心者の人
  • Pythonを使用した行列のLU分解の方法について学びたい人
levtech-ad

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(テックアカデミー)がおすすめです。

私も入っていますが、好きな時間に気軽にオンラインで学べますので、何より楽しいです。

現役エンジニアからマンツーマンで学べるので、一人では中々続かない人にも、向いていると思います。

無料体験ができますので、まずは試してみてください!

\まずは無料体験!/
スポンサーリンク
タイトルとURLをコピーしました