LoginSignup
4
4

More than 3 years have passed since last update.

Pythonによる線形代数入門:A=LU

Last updated at Posted at 2020-11-15

LU分解とは、正方行列 A を下三角行列 L と上三角行列 U の積、すなわち A = LU が成立するような L と Uに分解することです。連立方程式の解法、逆行列の計算、行列式の計算などに用いられます。


毎週木曜の7時から8時半まで 線形代数 オンライン勉強会
初心者大歓迎 xx回目 線形代数イントロダクション Gストラング
を行っています。奮ってご参加ください。


A=LU

3x3の行列$A$を消去により対角要素にピボットをもつ上三角行列$U$と下三角行列$L$に分解できることを確認します。

消去による対角要素にピボットをもつ上三角行列Uの取得

Uは対角要素にピボットをもつ上三角行列です。Lは下三角行列です。Aは消去の手順を通して上三角行列と下三角行列に分解することができます。

  A=\left(
    \begin{array}{cc}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      2 & 7 & 9 \\
    \end{array}
  \right)

とすると、$A$の3行1列目の要素($l_{31}$)2を消去するには3行目から1行目の2倍を引きます。

  \left(
    \begin{array}{cc}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      0 & 3 & 7 \\
    \end{array}
  \right)

つぎに3行目2列目の要素($l_{32}$)3は2行目を3倍して3行目から引くことで消去できます。

  \left(
    \begin{array}{cc}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      0 & 0 & 4 \\
    \end{array}
  \right)

この手順をまとめてかくと

  \left(
    \begin{array}{cc}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      2 & 7 & 9 \\
    \end{array}
  \right)
\ \rightarrow \ 
  \left(
    \begin{array}{cc}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      0 & 3 & 7 \\
    \end{array}
  \right)
\ \rightarrow \ 
  \left(
    \begin{array}{cc}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      0 & 0 & 4 \\
    \end{array}
  \right)

となります。

行列を使った消去

ある行列を$A$の左から掛けると消去の手順を行列を使って行うことができます。まず$A$の要素($l_{21}$)2を消去します。

  \left(
    \begin{array}{rr}
      1 & 0 & 0 \\
      0 & 1 & 0 \\
      -2 & 0 & 1 \\
    \end{array}
  \right)

  \left(
    \begin{array}{rr}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      2 & 7 & 9 \\
    \end{array}
  \right)
  =
  \left(
    \begin{array}{cc}
      1\times1+ 0\times0 +0\times2 & 1\times2+ 0\times1 +0\times7 & 1\times1+ 0\times1 +0\times9 \\
      0\times1+ 1\times0 +0\times2 & 0\times2+ 1\times1 +0\times7 & 0\times1+ 1\times1 +0\times9 \\
      -2\times1+0\times0 +1\times2 &-2\times2+ 0\times1 +1\times7 & -2\times1+ 0\times1 +1\times9 \\
    \end{array}
  \right)
  =
  \left(
    \begin{array}{rr}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      0 & 3 & 7 \\
    \end{array}
  \right)

となります。

この消去に用いた行列は消したい要素の正負を逆にした乗数$l_{31}=-2$を
単位行列$I$

  \left(
    \begin{array}{cc}
      1 & 0 & 0 \\
      0 & 1 & 0 \\
      0 & 0 & 1 \\
    \end{array}
  \right)

に加えたものです。$E_{31}$と書きます。これは

  E_{31}A=
  \left(
    \begin{array}{cc}
      1 & 0 & 0 \\
      0 & 1 & 0 \\
      -2 &  0 & 1 \\
    \end{array}
  \right)
  \left(
    \begin{array}{cc}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      2 & 7 & 9 \\
    \end{array}
  \right)
=
  \left(
    \begin{array}{cc}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      0 & 3 & 7 \\
    \end{array}
  \right)

となります。

つぎに右辺の要素($l_{32}$)3を消去します。したがって、$l_{32}=-3$の$E_{32}$を左から掛けると、

  \left(
    \begin{array}{cc}
      1 & 0 & 0 \\
      0 & 1 & 0 \\
      0 & -3 & 1 \\
    \end{array}
  \right)

  \left(
    \begin{array}{cc}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      0 & 3 & 7 \\
    \end{array}
  \right)
  =
  \left(
    \begin{array}{cc}
      1\times1+  0 \times0+0\times0 & 1\times2+  0\times1 +0\times3 & 1\times1+  0\times1 +0\times7 \\
      0\times1+  1 \times0+0\times0 & 0\times2+  1\times1 +0\times3 & 0\times1+  0\times1 +0\times7 \\
      0\times1+(-3)\times0+1\times0&0\times2+(-3)\times1+1\times3 & 0\times1+(-3)\times1+1\times7 \\
    \end{array}
  \right)
  =
  \left(
    \begin{array}{rr}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      0 & 0 & 4 \\
    \end{array}
  \right)

となります。

これは

E_{32}
  \left(
    \begin{array}{cc}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      0 & 3 & 7 \\
    \end{array}
  \right)
  =
  \left(
    \begin{array}{cc}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      0 & 0 & 4 \\
    \end{array}
  \right)
  =U

です。また

  E_{31}A=
  \left(
    \begin{array}{cc}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      0 & 3 & 7 \\
    \end{array}
  \right)

から

E_{32}E_{31}A=
E_{32}
  \left(
    \begin{array}{cc}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      0 & 3 & 7 \\
    \end{array}
  \right)
  =U

となります。これでピッボト付き上三角行列が得られました。

逆行列を使った表現

$E_{32}^{-1}$は$E_{32}$の逆行列を表し、$E_{32}^{-1}E_{32}=I$となります。逆行列の求め方は次節で示します。
$E_{32}E_{31}A=U$両辺に$E_{32}^{-1}$を掛けると

E_{32}^{-1}E_{32}E_{31}A=IE_{31}A=E_{31}A=E_{32}^{-1}U

となります。これは

E_{31}A=E_{32}^{-1}U

さらに両辺に$E_{31}^{-1}$を掛けると

E_{31}^{-1}E_{31}A=E_{31}^{-1}E_{32}^{-1}U

ガウスジョルダン法による逆行列の取得

つぎに$E_{31}^{-1}$と$E_{32}^{-1}$を求めます。これはガウスジョルダン法では、逆行列を求めたい行列を最初の3行3列に書き、その後の3行3列が単位行列からなる長方行列を作り、これを消去により

  [E_{31} e_1 e_2 e_3]=
  \left(
    \begin{array}{cc}
      1 & 0 & 0 &1&0&0\\
      0 & 1 & 0 &0&1&0\\
      -2 & 0 & 1&0&0&1 \\
    \end{array}
  \right)
  \rightarrow
  \left(
    \begin{array}{cc}
      1 & 0 & 0 &1&0&0\\
      0 & 1 & 0 &0&1&0\\
      0 & 0 & 1&2&0&1 \\
    \end{array}
  \right)

と計算していきます。右端の3x3の行列が逆行列になります。

  E_{31}^{-1}=
  \left(
    \begin{array}{cc}
      1&0&0\\
      0&1&0\\
      2&0&1 \\
    \end{array}
  \right)

が得らました。また、

  [E_{32} e_1 e_2 e_3]=
  \left(
    \begin{array}{cc}
      1 & 0 & 0 &1&0&0\\
      0 & 1 & 0 &0&1&0\\
      0 & -3 & 1&0&0&1 \\
    \end{array}
  \right)
  \rightarrow
  \left(
    \begin{array}{cc}
      1 & 0 & 0 &1&0&0\\
      0 & 1 & 0 &0&1&0\\
      0 & 0 & 1&0&3&1 \\
    \end{array}
  \right)

から

  E_{32}^{-1}=
  \left(
    \begin{array}{cc}
      1&0&0\\
      0&1&0\\
      0&3&1 \\
    \end{array}
  \right)

が得られます。

下三角行列の取得

つぎに

  E_{31}^{-1}E_{32}^{-1}=
  \left(
    \begin{array}{cc}
      1&0&0\\
      0&1&0\\
      2&0&1 \\
    \end{array}
  \right)

  \left(
    \begin{array}{cc}
      1&0&0\\
      0&1&0\\
      0&3&1 \\
    \end{array}
  \right)
  =
  \left(
    \begin{array}{cc}
      1&0&0\\
      0&1&0\\
      2&3&1 \\
    \end{array}
  \right)

これは下三角行列であるので

ゆえに

  A=LU
  =
  \left(
    \begin{array}{cc}
      1&0&0\\
      0&1&0\\
      2&3&1 \\
    \end{array}
  \right)
  \left(
    \begin{array}{cc}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      0 & 0 & 4 \\
    \end{array}
  \right)  

となります。$E_{31}, E_{32}$がわかれば下三角行列がわかります。

三角行列への分解

ピッボト付き上三角行列を対角行列を用いて、ピッボトなしの上三角行列に変換します。

  A=LU
  =
  \left(
    \begin{array}{cc}
      1&0&0\\
      0&1&0\\
      2&3&1 \\
    \end{array}
  \right)
  \left(
    \begin{array}{cc}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      0 & 0 & 4 \\
    \end{array}
  \right)  

は対角行列$D$を用いて

  A
  =
  \left(
    \begin{array}{cc}
      1&0&0\\
      0&1&0\\
      2&3&1 \\
    \end{array}
  \right)
  \left(
    \begin{array}{cc}
      1 & 0 & 0 \\
      0 & 1 & 0 \\
      0 & 0 & 4 \\
    \end{array}
  \right)  
  \left(
    \begin{array}{cc}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      0 & 0 & 1 \\
    \end{array}
  \right)  
  = LDU_u = LDU

と書き換えることができます。

対称行列の例

$A$を対称行列とすると$A=LDU=LDL^{T}$となることを確認します。

  A=\left(
    \begin{array}{rrr}
       1 & -1 & 2 \\
      -1 & 2 & -1 \\
       2 & -1 & 1 \\
    \end{array}
  \right)

これを行列により消去し、ピボット付き上三角行列を作ります。

  E_{32}E_{31}E_{21}A=\left(
    \begin{array}{rrr}
       1 & -1 & 2 \\
       0& 1 & 1 \\
       0 & 0 & -4 \\
    \end{array}
  \right)

これは

  U=DU_u=
  \left(
    \begin{array}{rrr}
       1 & 0 & 0 \\
       0& 1 & 0 \\
       0 & 0 & -4 \\
    \end{array}
  \right)
  \left(
    \begin{array}{rrr}
       1 & -1 & 2 \\
       0& 1 & 1 \\
       0 & 0 & 1 \\
    \end{array}
  \right)

に分解できます。

また、下三角行列は

  E_{32}^{-1}E_{31}^{-1}E_{21}^{-1}=\left(
    \begin{array}{rrr}
       1 & 0 & 0 \\
       -1& 1 & 0 \\
       2 & 1 & 1 \\
    \end{array}
  \right)
  = U_u^{T}

となり
$$A=LDL^T$$

が確認できます。

sciPy、numpyによる確認

numpy, scipyを用いて、上述の手順を確認します。

import numpy as np
from scipy.linalg import inv
A = np.array([[1, 2, 1], [0, 1, 1], [2, 7, 9]])
E31=np.array([[1, 0, 0], [0, 1, 0], [-2, 0, 1]])
E32=np.array([[1, 0, 0], [0, 1, 0], [0, -3, 1]])
U=E32@E31@A
L=inv(E32)@inv(E31)
print("U:\n{}\n".format(U))
print("L:\n{}\n".format(L))
print("A:\n{}\n".format(L@U))
UU=np.triu(U,1)+np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
D=np.diag(U)
D=np.diag(D)
print("UU:\n{}\n".format(UU))
print("D:\n{}\n".format(D))
AA=L@D@(UU)
print("A:\n{}\n".format(AA))

U:
[[1 2 1]
[0 1 1]
[0 0 4]]

L:
[[1. 0. 0.]
[0. 1. 0.]
[2. 3. 1.]]

A:
[[1. 2. 1.]
[0. 1. 1.]
[2. 7. 9.]]

UU:
[[1 2 1]
[0 1 1]
[0 0 1]]

D:
[[1 0 0]
[0 1 0]
[0 0 4]]

A:
[[1. 2. 1.]
[0. 1. 1.]
[2. 7. 9.]]

これで確認ができました。

つぎに、対称行列についても確認してみます。

from scipy.linalg import inv
A = np.array([[1, -1, 2], [-1, 2, -1], [2, -1, 1]])
print("A:\n{}\n".format(A))
E21=np.array([[1, 0, 0], [1, 1, 0], [0, 0, 1]])
E31=np.array([[1, 0, 0], [0, 1, 0], [-2, 0, 1]])
E32=np.array([[1, 0, 0], [0, 1, 0], [0, -1, 1]])
U=E32@E31@E21@A
L=inv(E21)@inv(E31)@inv(E32)
print("U:\n{}\n".format(U))
print("L:\n{}\n".format(L))
print("A:\n{}\n".format(L@U))
UU=np.triu(U,1)+np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
D=np.diag(U)
D=np.diag(D)
print("Uu:\n{}\n".format(UU))
print("D:\n{}\n".format(D))
AA=L@D@(UU)
print("LDUu:\n{}\n".format(AA))

A:
[[ 1 -1 2]
[-1 2 -1]
[ 2 -1 1]]

U:
[[ 1 -1 2]
[ 0 1 1]
[ 0 0 -4]]

L:
[[ 1. 0. 0.]
[-1. 1. 0.]
[ 2. 1. 1.]]

A:
[[ 1. -1. 2.]
[-1. 2. -1.]
[ 2. -1. 1.]]

Uu:
[[ 1 -1 2]
[ 0 1 1]
[ 0 0 1]]

D:
[[ 1 0 0]
[ 0 1 0]
[ 0 0 -4]]

LDUu:
[[ 1. -1. 2.]
[-1. 2. -1.]
[ 2. -1. 1.]]

これで確認ができました。

A=LUの別解法

Aの第一ピボット行を$u_1^$、第一列を$l_1$とします。Aを$l_1 u_1^$とそれ以外に分解します。

  A=
  \left(
    \begin{array}{ccc}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      2 & 7 & 9 \\
    \end{array}
  \right)
 =l_1 u_1^* +A_2=
  \left(
    \begin{array}{ccc}
      1 \times& 1行目& \\
      0 \times& 1行目& \\
      2 \times& 1行目& \\
    \end{array}
  \right)
+
  \left(
    \begin{array}{ccc}
      0 & 0 & 0 \\
      0 & 1 & 1 \\
      0 & 3 & 7 \\
    \end{array}
  \right)
 =
  \left(
    \begin{array}{ccc}
      1 & 2 & 1 \\
      0 & 0 & 0 \\
      2 & 4 & 2 \\
    \end{array}
  \right)
+
  \left(
    \begin{array}{ccc}
      0 & 0 & 0 \\
      0 & 1 & 1 \\
      0 & 3 & 7 \\
    \end{array}
  \right)

右辺の2番目の行列を$A_2$とします。$A$の2番目のピボット行を$u_2^*$としてそれと$A_2$の2列目を$l_2$としてそれらの積とそのあまりの行列に分解していきます。

 A=l_1 u_1^* +A_2=l_1 u_1^* +
  \left(
    \begin{array}{ccc}
      0 & 0 & 0 \\
      0 & 1 & 1 \\
      0 &  3 & 7 \\
    \end{array}
  \right)

= l_1 u_1^* + l_2 u_2^*
+
  \left(
    \begin{array}{ccc}
      0 & 0 & 0 \\
      0 & 0 & 0\\
      0 & 0 & 4 \\
    \end{array}
  \right)
 = l_1 u_1^* +
  \left(
    \begin{array}{ccc}
      0 \times& 2行目& \\
      1 \times& 2行目& \\
      3 \times& 2行目& \\
    \end{array}
  \right)
+
  \left(
    \begin{array}{ccc}
     0 & 0 & 0 \\
     0 & 0 & 0 \\
     0 & 0 & 4 \\
    \end{array}
  \right)
 = l_1 u_1^* +
  \left(
    \begin{array}{ccc}
      0 & 0 & 0 \\
      0 & 1 & 1 \\
      0 & 3 & 3 \\
    \end{array}
  \right)
+
  \left(
    \begin{array}{ccc}
     0 & 0 & 0 \\
     0 & 0 & 0 \\
     0 & 0 & 4 \\
    \end{array}
  \right)

右辺の3番目の行列を$A_3$としますが、これはピボット行しかないのでここで分解を終了とします。
したがって、

  A=
  \left(
    \begin{array}{ccc}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      2 & 7 & 9 \\
    \end{array}
  \right)
 =
  \left(
    \begin{array}{ccc}
      1 & 2 & 1 \\
      0 & 0 & 0 \\
      2 & 4 & 2 \\
    \end{array}
  \right)
+
  \left(
    \begin{array}{ccc}
      0 & 0 & 0 \\
      0 & 1 & 1 \\
      0 & 3 & 3 \\
    \end{array}
  \right)
+
  \left(
    \begin{array}{ccc}
      0 & 0 & 0 \\
      0 & 0 & 0 \\
      0 & 0 & 4 \\
    \end{array}
  \right)


 =l_1 u_1^*+l_2 u_2^* +l_3 u_3^* =
  \left(
    \begin{array}{ccc}
      1 & 0 & 0 \\
      0 & 1 & 0 \\
      2 & 3 & 1 \\
    \end{array}
  \right)
  \left(
    \begin{array}{cc}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      0 & 0 & 4 \\
    \end{array}
  \right)-------(a)

PA=LUについて

Aの第一ピボット行を1列目の最も大きな数値をもつ3行目とします。$u_1^* $、第一列の1/2を$l_1$とします。Aを$l_1 u_1^*$とそれ以外に分解します。

  A=
  \left(
    \begin{array}{ccc}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      2 & 7 & 9 \\
    \end{array}
  \right)
 =l_1 u_1^* +A_2=
  \left(
    \begin{array}{ccc}
      0.5 \times& 3行目& \\
      0 \times& 3行目& \\
      1 \times& 3行目& \\
    \end{array}
  \right)
+
  \left(
    \begin{array}{ccc}
      0 & -3/2 & -7/2 \\
      0 & 1 & 1 \\
      0 & 0 & 0 \\
    \end{array}
  \right)
 =
  \left(
    \begin{array}{ccc}
      1 & 3.5 & 4.5 \\
      0 & 0 & 0 \\
      2 & 7 & 9 \\
    \end{array}
  \right)
+
  \left(
    \begin{array}{ccc}
      0 & -1.5 & -3.5 \\
      0 & 1 & 1 \\  
      0 & 0 & 0 \\
    \end{array}
  \right) 

右辺の2番目の行列を$A_2$とします。$A_2$の2列目の一番大きな数値は2行目なので、2番目をピボット行$u_2^*$としてそれと$A_2$の2列目を$l_2$としてそれらの積とそのあまりの行列に分解していきます。

 A=l_1 u_1^* +A_2=l_1 u_1^* +
  \left(
    \begin{array}{ccc}
      0 & -3/2 & -7/2 \\
      0 & 1 & 1 \\
      0 &  0 & 0 \\
    \end{array}
  \right)

= l_1 u_1^* + l_2 u_2^*
+
  \left(
    \begin{array}{ccc}
      0 & 0 & -2 \\
      0 & 0 & 0\\
      0 & 0 & 0 \\
    \end{array}
  \right)
 = l_1 u_1^* +
  \left(
    \begin{array}{cc}
      -3/2 \times& 2行目& \\
      1 \times& 2行目& \\
      0 \times& 2行目& \\
    \end{array}
  \right)
+
  \left(
    \begin{array}{cc}
     0 & 0 & -2 \\
     0 & 0 & 0 \\
     0 & 0 & 0 \\
    \end{array}
  \right)
 = l_1 u_1^* +
  \left(
    \begin{array}{ccc}
      0 & -3/2 & -3/2 \\
      0 & 1 & 1 \\
      0 &  0 & 0 \\
    \end{array}
  \right)+
  \left(
    \begin{array}{cc}
     0 & 0 & -2 \\
     0 & 0 & 0 \\
     0 & 0 & 0 \\
    \end{array}
  \right)
=l_1 u_1^* + l_2 u_2^*+ l_3 u_3^*= 
  \left(
    \begin{array}{ccc}
      0.5 & -3/2 & 1 \\
      0 & 1 & 0 \\
      1 &  0 & 0 \\
    \end{array}
  \right)
  \left(
    \begin{array}{cc}
     2 & 7 & 9 \\
     0 & 1 & 1 \\
     0 & 0 & -2 \\
    \end{array}
  \right)

これでは右辺の右側の行列は上三角行列ですが左側は下三角行列になっていません。これは行を入れ替える必要があります。ここではAのピボット行を321としてしまいました。実際は123である必要があります。そこで

  A=
  \left(
    \begin{array}{ccc}
      2 & 7 & 9 \\
      1 & 2 & 1 \\
      0 & 1 & 1 \\
    \end{array}
  \right)
 =l_1 u_1^* +A_2=
  \left(
    \begin{array}{ccc}
      1 \times& 1行目& \\
      0.5 \times& 1行目& \\
      0 \times& 1行目& \\
    \end{array}
  \right)
+
  \left(
    \begin{array}{ccc}
      2 & 7 & 9 \\
      1 & 3.5 & 4.5 \\
      0 & 0 & 0 \\
    \end{array}
  \right)
 =
  \left(
    \begin{array}{ccc}
      2 & 7 & 9 \\
      1 & 3.5 & 4.5 \\
      0 & 0 & 0 \\
    \end{array}
  \right)
+
  \left(
    \begin{array}{ccc}
      0 & 0 & 0 \\
      0 & -1.5 & -3.5 \\
      0 & 1 & 1 \\
    \end{array}
  \right)
 =l_1 u_1^* +l_2 u_2^* +A_3=
  l_1 u_1^* +
  \left(
    \begin{array}{ccc}
      0 \times& 2行目& \\
      1 \times& 2行目& \\
      -1/1.5 \times& 2行目& \\
    \end{array}
  \right)
+
  \left(
    \begin{array}{ccc}
      0 & 0 & 0 \\
      0 & 0 & 0 \\
      0 & 0 & -1.33 \\
    \end{array}
  \right)
 = l_1 u_1^* +
  \left(
    \begin{array}{ccc}
      0 & 0 & 0 \\
      0 & -1.5 & -3.5 \\
      0 & 1 & -2.333 \\
    \end{array}
  \right)
+
  \left(
    \begin{array}{ccc}
      0 & 0 & 0 \\
      0 & 0 & 0 \\
      0 & 0 & -1.33 \\
    \end{array}
  \right)
 =
  \left(
    \begin{array}{ccc}
      1 & 0 & 0 \\
      0.5 & 1 & 0 \\
      0 & -0.66 & 1 \\
    \end{array}
  \right)
  \left(
    \begin{array}{ccc}
      2 & 7 & 9 \\
      0 & -1.5 & -3.5 \\
      0 & 0 & -1.33 \\
    \end{array}
  \right)

これを実際にScipy.linalgで実行してみます。

import numpy as np
from scipy.linalg import lu

A = np.array([[2, 7, 9],[1, 2, 1], [0, 1, 1]])

# PA=LU分解
P, L, U = lu(A,overwrite_a=True)

print("A:\n{}\n".format(A))
print("P:\n{}\n".format(P))
print("L:\n{}\n".format(L))
print("U:\n{}\n".format(U))
print("PLU:\n{}".format(P@L@U))
#
A:
[[2 7 9]
 [1 2 1]
 [0 1 1]]

P:
[[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]

L:
[[ 1.          0.          0.        ]
 [ 0.5         1.          0.        ]
 [ 0.         -0.66666667  1.        ]]

U:
[[ 2.          7.          9.        ]
 [ 0.         -1.5        -3.5       ]
 [ 0.          0.         -1.33333333]]

PLU:
[[2. 7. 9.]
 [1. 2. 1.]
 [0. 1. 1.]]

つぎに

  A=
  \left(
    \begin{array}{ccc}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      2 & 7 & 9 \\
    \end{array}
  \right)

について実行してみます。

A = np.array([[1, 2, 1], [0, 1, 1], [2, 7, 9]])

# PA=LU分解
P, L, U = lu(A,overwrite_a=True)

print("A:\n{}\n".format(A))
print("P:\n{}\n".format(P))
print("L:\n{}\n".format(L))
print("U:\n{}\n".format(U))
print("PLU:\n{}".format(P@L@U))
#
A:
[[1 2 1]
 [0 1 1]
 [2 7 9]]

P:
[[0. 1. 0.]
 [0. 0. 1.]
 [1. 0. 0.]]

L:
[[ 1.          0.          0.        ]
 [ 0.5         1.          0.        ]
 [ 0.         -0.66666667  1.        ]]

U:
[[ 2.          7.          9.        ]
 [ 0.         -1.5        -3.5       ]
 [ 0.          0.         -1.33333333]]

PLU:
[[1. 2. 1.]
 [0. 1. 1.]
 [2. 7. 9.]]

結果は(a)とは異なるものになりました。これはScipy.linalg.luがAの一列目の一番大きな数値の行を第一ピボット行として計算し始めるからだと思います。置換行列Pを見てみるとわかります。そのためにLは三角行列にはなりませんが、Lを下三角行列としているのでPLU=Aとなります。

4
4
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
4
4