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

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

---

毎週木曜の７時から８時半まで　線形代数　オンライン勉強会
[初心者大歓迎 ｘｘ回目　線形代数イントロダクション Gストラング](https://kinyuzaimu.connpass.com/)
を行っています。奮ってご参加ください。

---

# A=LU

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

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

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

```math
  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倍を引きます。

```math
  \left(
    \begin{array}{cc}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      0 & 3 & 7 \\
    \end{array}
  \right)
```
つぎに3行目2列目の要素($l_{32}$)3は2行目を3倍して3行目から引くことで消去できます。

```math
  \left(
    \begin{array}{cc}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      0 & 0 & 4 \\
    \end{array}
  \right)
```
この手順をまとめてかくと

```math
  \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を消去します。

```math
  \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)
```

```math
  =
  \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)
```

```math
  =
  \left(
    \begin{array}{rr}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      0 & 3 & 7 \\
    \end{array}
  \right)
```
となります。

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

```math
  \left(
    \begin{array}{cc}
      1 & 0 & 0 \\
      0 & 1 & 0 \\
      0 & 0 & 1 \\
    \end{array}
  \right)
```
に加えたものです。$E_{31}$と書きます。これは

```math
  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}$を左から掛けると、

```math
  \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)
```

```math
  =
  \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)
```

```math
  =
  \left(
    \begin{array}{rr}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      0 & 0 & 4 \\
    \end{array}
  \right)
```
となります。

これは

```math
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
```
です。また

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

```math
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}$を掛けると

```math
E_{32}^{-1}E_{32}E_{31}A=IE_{31}A=E_{31}A=E_{32}^{-1}U
```
となります。これは

```math
E_{31}A=E_{32}^{-1}U
```
さらに両辺に$E_{31}^{-1}$を掛けると

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


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

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

```math
  [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の行列が逆行列になります。

```math 
  E_{31}^{-1}=
  \left(
    \begin{array}{cc}
      1&0&0\\
      0&1&0\\
      2&0&1 \\
    \end{array}
  \right)
```
が得らました。また、  

```math 
  [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)
```
から

```math 
  E_{32}^{-1}=
  \left(
    \begin{array}{cc}
      1&0&0\\
      0&1&0\\
      0&3&1 \\
    \end{array}
  \right)
```
が得られます。

## 下三角行列の取得

つぎに

```math
  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)
```
これは下三角行列であるので

ゆえに

```math
  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}$がわかれば下三角行列がわかります。

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

```math
  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$を用いて

```math
  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}$となることを確認します。

```math
  A=\left(
    \begin{array}{rrr}
       1 & -1 & 2 \\
      -1 & 2 & -1 \\
       2 & -1 & 1 \\
    \end{array}
  \right)
```
これを行列により消去し、ピボット付き上三角行列を作ります。

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

```math
  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)
```
に分解できます。

また、下三角行列は

```math
  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を用いて、上述の手順を確認します。

```Python
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.]]


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

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

```Python
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^*$とそれ以外に分解します。

```math
  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)
```
右辺の２番目の行列を$A_2$とします。$A$の2番目のピボット行を$u_2^*$としてそれと$A_2$の２列目を$l_2$としてそれらの積とそのあまりの行列に分解していきます。

```math
　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)
```
```math
　= 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)
```

```math
　= 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)
```
右辺の３番目の行列を$A_3$としますが、これはピボット行しかないのでここで分解を終了とします。
したがって、

```math
  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)

```

```math

　=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の第一ピボット行を１列目の最も大きな数値をもつ３行目とします。$u_1^\* $、第一列の1/2を$l_1$とします。Aを$l_1 u_1^\*$とそれ以外に分解します。

```math
  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)
```

```math
　=
  \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) 
```

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

```math
　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)
```

```math
　= 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)
```

```math
　= 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)
```

```math
=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である必要があります。そこで

```math
  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)
```
```math
　=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)
```
```math
　=
  \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で実行してみます。

```Python
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.]]
```

つぎに

```math
  A=
  \left(
    \begin{array}{ccc}
      1 & 2 & 1 \\
      0 & 1 & 1 \\
      2 & 7 & 9 \\
    \end{array}
  \right)
```
について実行してみます。

```Python
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となります。