LU decomposition means decomposing a square matrix A into L and U such that the product of the lower triangular matrix L and the upper triangular matrix U, that is, A = LU holds. It is used for solving simultaneous equations, calculating inverse matrices, and calculating determinants. There are many solutions, such as analytical solutions, Gaussian elimination, and recursive methods. The pair of L and U seems to be uniquely determined, but it depends on the conditions.
November 19 Linear algebra online study session Beginners welcome 19th Linear Algebra Introduction G Strang
to hold. Please join us.
Confirm that the 3x3 matrix $ A $ can be decomposed into the upper triangular matrix $ U $ and the lower triangular matrix $ L $ with pivots on the diagonal elements.
U is an upper triangular matrix with pivots on diagonal elements. L is the lower triangular matrix. A can be decomposed into an upper triangular matrix and a lower triangular matrix through the elimination procedure.
A=\left(
\begin{array}{cc}
1 & 2 & 1 \\
0 & 1 & 1 \\
2 & 7 & 9 \\
\end{array}
\right)
Then, to erase element 2 in the 3rd row and 1st column of $ A $, subtract twice from the 3rd row to the 1st row.
\left(
\begin{array}{cc}
1 & 2 & 1 \\
0 & 1 & 1 \\
0 & 3 & 7 \\
\end{array}
\right)
Next, element 3 in the 3rd row and 2nd column can be deleted by multiplying the 2nd row by 3 and subtracting it from the 3rd row.
\left(
\begin{array}{cc}
1 & 2 & 1 \\
0 & 1 & 1 \\
0 & 0 & 4 \\
\end{array}
\right)
To summarize this procedure
\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)
Will be.
By multiplying a matrix from the left of $ A $, the elimination procedure can be performed using the matrix. First, delete element 2 of $ A $.
\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)
The matrix used for this erasure is the multiplier $ l_ {31} =-2 $, which is the reverse of the positive and negative elements to be erased, added to the identity matrix $ I $, and is written as $ E_ {31} $. this is
E_{31}A=
\left(
\begin{array}{cc}
1 & 2 & 1 \\
0 & 1 & 1 \\
0 & 3 & 7 \\
\end{array}
\right)
Can be written.
Next, I want to erase the element $ 3 $ on the right side. Therefore, multiplying $ E_ {32} $ of $ l_ {32} =-3 $ from the left,
\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)
this is
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
Will be. Also
E_{31}A=
\left(
\begin{array}{cc}
1 & 2 & 1 \\
0 & 1 & 1 \\
0 & 3 & 7 \\
\end{array}
\right)
From
E_{32}
\left(
\begin{array}{cc}
1 & 2 & 1 \\
0 & 1 & 1 \\
0 & 3 & 7 \\
\end{array}
\right)
=E_{32}E_{31}A=U
Call. This gives an upper triangular matrix with pivots.
Multiply both sides by $ E_ {32} ^ {-1} $
E_{32}^{-1}E_{32}E_{31}A=E_{32}^{-1}U
E_{32}^{-1}E_{32}E_{31}A=IE_{31}A=E_{31}A=E_{32}^{-1}U
this is
E_{31}A=E_{32}^{-1}U
If you multiply both sides by $ E_ {31} ^ {-1} $
E_{31}^{-1}E_{31}A=E_{31}^{-1}E_{32}^{-1}U
this is
A=E_{31}^{-1}E_{32}^{-1}U
Also call.
Next, find $ E_ {31} ^ {-1} $ and $ E_ {32} ^ {-1} $. In the Gauss-Jordan method, the matrix for which the inverse matrix is to be obtained is written in the first 3 rows and 3 columns, and then a rectangular matrix consisting of 3 rows and 3 columns is created by elimination.
[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)
I will calculate.
E_{31}^{-1}=
\left(
\begin{array}{cc}
1&0&0\\
0&1&0\\
2&0&1 \\
\end{array}
\right)
was gotten. Also,
[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)
From
E_{32}^{-1}=
\left(
\begin{array}{cc}
1&0&0\\
0&1&0\\
0&3&1 \\
\end{array}
\right)
Is obtained.
Next
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)
Because this is a lower triangular matrix
therefore
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)
Disassembly is complete.
Convert an upper triangular matrix with pivots to an upper triangular matrix without pivots using a diagonal matrix.
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)
Uses the diagonal matrix $ 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
Can be rewritten as.
Confirm that $ A = LDU = LDL ^ {T} $ when $ A $ is a symmetric matrix.
A=\left(
\begin{array}{rrr}
1 & -1 & 2 \\
-1 & 2 & -1 \\
2 & -1 & 1 \\
\end{array}
\right)
Eliminate this with a matrix to create a pivoted upper triangular matrix.
E_{32}E_{31}E_{21}A=\left(
\begin{array}{rrr}
1 & -1 & 2 \\
0& 1 & 1 \\
0 & 0 & -4 \\
\end{array}
\right)
this is
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)
Can be disassembled into.
Also, the lower triangular matrix
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}
Can be confirmed.
Check the above procedure using numpy and 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.]]
I was able to confirm this.
Next, let's check the symmetric matrix.
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.]]
I was able to confirm this.
Next, let's decompose using scipy's lu function and lul function. Note that the pivot may be attached to the lower triangular matrix. Moreover, it can be seen that the LU decomposition is not uniquely determined.
scipy.linalg.lu(a, permute_l=False, overwrite_a=False, check_finite=True)
Compute the pivot LU decomposition of the matrix.
Disassembly
A = P L U
P is the permutation matrix, The diagonal element of L is 1.
import numpy as np
from scipy.linalg import lu
#np.random.seed(10)
#Generate a matrix with random elements
#A = np.random.randint(1, 5, (5, 5))
A = np.array([[1, 2, 1], [0, 1, 1], [2, 7, 9]])
# PA=LU decomposition
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{}\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.]]
scipy.linalg.ldl(A, lower=True, hermitian=True, overwrite_a=False, check_finite=True) Compute the LDLt or Bunch-Kaufman decomposition of a symmetric / Hermitian matrix.
This function returns a block diagonal matrix D consisting of blocks of up to 2x2 and a lower triangular matrix L that holds the decomposition A = L D L ^ H or A = L D L ^ T and may require permutation. If lower is False, the upper triangular matrix (which may require permutation) returns as the return value.
You can use the permutation matrix to triangulate the return value by simply shuffling the rows. This corresponds to multiplication using the permutation matrix P.dot (lu), where P is the permutation matrix I [:, perm].
from scipy.linalg import ldl
A = np.array([[1, -1, 2], [-1, 2, -1], [2, -1, 1]])
lu, d, perm = ldl(A)
print("lu:\n{}\n".format(lu[perm,:]))
print("d:\n{}\n".format(d))
print("lu.t:\n{}\n".format(lu[perm,:].T))
print("A:\n{}\n".format(lu.dot(d).dot(lu.T)))
lu: [[ 1. 0. 0. ] [ 0. 1. 0. ] [-0.33333333 -0.33333333 1. ]]
d: [[1. 2. 0. ] [2. 1. 0. ] [0. 0. 1.33333333]]
lu.t: [[ 1. 0. -0.33333333] [ 0. 1. -0.33333333] [ 0. 0. 1. ]]
A: [[ 1. -1. 2.] [-1. 2. -1.] [ 2. -1. 1.]]