There are many articles about LU decomposition, but there was no article that mentioned the feature of LU decomposition that calculations can be performed with memory saving, so I tried to explain LU decomposition with a little more attention to memory saving. think.
After studying React, I created a site that can easily perform LU decomposition. However, for the convenience of rendering, the memory-saving implementation described below is not implemented. LU decomposition experience
Before discussing LU decomposition, let's talk about the triangular matrix. The goal of LU decomposition is to ** represent the coefficients of simultaneous equations with only a triangular matrix **. I'll explain why. First, consider a simple system of equations.
\left\{
\begin{array}{rrrrrrr}
2x_1 & & & & &=& 2 \\
x_1 &+& 4x_2 & & &=& 9 \\
4x_1 &-& 3x_2 &+& 3x_3 &=& -5
\end{array}
\right.
The solution of this simultaneous equation can be easily obtained by solving it sequentially.
\begin{align}
x_1 &= \frac{2}{2} = 1 \\
x_2 &= \frac{9 - x_1}{4} = \frac{9-1}{4} = 2 \\
x_3 &= \frac{-5 - 4x_1 + 3x_2}{3} = \frac{-5 - 4 + 6}{3} = -1
\end{align}
Here, if this simultaneous equation is expressed by a matrix,
\left(
\begin{matrix}
2 & 0 & 0 \\
1 & 4 & 0 \\
4 & -3 & 3
\end{matrix}
\right)
\left(
\begin{matrix}
x_1 \\
x_2 \\
x_3
\end{matrix}
\right)
=
\left(
\begin{matrix}
2 \\
9 \\
-5
\end{matrix}
\right) \tag{1}
It looks like. You can see that the coefficients of the simultaneous equations are represented by a ** lower triangular matrix ** whose elements above the diagonal are $ 0 $. In this case, when you find $ x_i (i \ in 2, 3) $, you already know $ x_ {i-1} (i \ in 1,2) $, so you can solve them in order. This method is called ** forward substitution **. Next, consider the following simultaneous equations.
\left\{
\begin{array}{rrrrrrr}
3x_1 &-& 3x_2 &+& 4x_3 &=& -5 \\
& & 4x_2 &+& x_3 &=& 9 \\
& & & & 2x_3 &=& 2
\end{array}
\right. \\
\left(
\begin{matrix}
3 & -3 & 4\\
0 & 4 & 1 \\
0 & 0 & 2
\end{matrix}
\right)
\left(
\begin{matrix}
x_1 \\
x_2 \\
x_3
\end{matrix}
\right)
=
\left(
\begin{matrix}
-5 \\
9 \\
2
\end{matrix}
\right)
In this case, the coefficients of the simultaneous equations are represented by the ** upper triangular matrix **, where the components below the diagonal are $ 0 $. In this case, when finding $ x_i (i \ in 2, 1) $, we know $ x_ {i + 1} (i \ in 3, 2) $, so the reverse order in the case of the lower triangular matrix. Can be solved. This method is called ** backward substitution **. In this way, when the coefficients of simultaneous equations are represented by a triangular matrix, they can be solved by a simple algorithm. Here, I will define the words. The expression $ (1) $ is expressed as follows.
A \boldsymbol{x} = \boldsymbol{b} \\
A =
\left(
\begin{matrix}
3 & -3 & 4\\
0 & 4 & 1 \\
0 & 0 & 2
\end{matrix}
\right), \
\boldsymbol{x} =
\left(
\begin{matrix}
x_1 \\
x_2 \\
x_3
\end{matrix}
\right), \
\boldsymbol{b} =
\left(
\begin{matrix}
2 \\
9 \\
-5
\end{matrix}
\right)
Here, $ A $ is called the coefficient matrix, $ \ boldsymbol {x} $ is called the solution vector, and $ \ boldsymbol {b} $ is called the right-hand side vector. In the following, this notation is generalized, where $ A $ is the $ n \ times n $ matrix and $ \ boldsymbol {x}, \ \ boldsymbol {b} $ is the $ n $ dimensional vector.
As mentioned above, if the coefficients of simultaneous equations can be represented by a triangular matrix, it is almost like solving the equations. Now consider expressing the coefficient matrix $ A $ as the product of the lower triangular matrix $ L $ and the upper triangular matrix $ U $.
LU \boldsymbol{x} = \boldsymbol{b} \\
If $ U \ boldsymbol {x} = \ boldsymbol {y} $, then $ \ boldsymbol {y} $ can be obtained by forward substitution, and by using the result, $ \ boldsymbol {x} $ can be obtained by backward substitution. You can see that it is required.
\begin{align}
L \boldsymbol{y} &= \boldsymbol{b} \\
U \boldsymbol{x} &= \boldsymbol{y}
\end{align}
So how do you find $ L $ and $ U $? Before getting into the main subject, let's focus on the degree of freedom. There are $ n ^ 2 $ elements for $ A $ and $ n (n + 1) / 2 $ elements for $ L and \ U $, respectively, for a total of $ n ^ 2 + n $ elements. .. Therefore, after decomposition, the degrees of freedom are $ n $ larger, which means that all the diagonal components of $ L $ are $ 1 $, so the degrees of freedom are $ n ^ 2 $, which is $ A $. On the other hand, $ L and \ U $ can be uniquely determined. (By the way, it doesn't matter whether the diagonal component of $ L $ is $ 1 $ or the diagonal component of $ U $ is $ 1 $. It depends on the person.)
Now, I will explain how to find it concretely, but before that, I will explain the submatrix.
Consider the following matrix $ M $.
M = \left(
\begin{matrix}
m_{11} & m_{12} &|& m_{13} & m_{14} & m_{15} & m_{16} \\
m_{21} & m_{22} &|& m_{23} & m_{24} & m_{25} & m_{26} \\
- & - & - & - & - & - & - \\
m_{31} & m_{32} &|& m_{33} & m_{34} & m_{35} & m_{36} \\
m_{41} & m_{42} &|& m_{43} & m_{44} & m_{45} & m_{46} \\
m_{51} & m_{52} &|& m_{53} & m_{54} & m_{55} & m_{56} \\
m_{61} & m_{62} &|& m_{63} & m_{64} & m_{65} & m_{66}
\end{matrix}
\right)
Separate this matrix with a dotted line and give each block (** submatrix **) a name.
M_{11} =
\left(
\begin{matrix}
m_{11} & m_{12} \\
m_{21} & m_{22}
\end{matrix}
\right)
, \
M_{12} =
\left(
\begin{matrix}
m_{13} & m_{14} & m_{15} & m_{16} \\
m_{23} & m_{24} & m_{25} & m_{26}
\end{matrix}
\right) \\
M_{21} =
\left(
\begin{matrix}
m_{31} & m_{32} \\
m_{41} & m_{42} \\
m_{51} & m_{52} \\
m_{61} & m_{62}
\end{matrix}
\right)
, \
M_{22} =
\left(
\begin{matrix}
m_{33} & m_{34} & m_{35} & m_{36} \\
m_{43} & m_{44} & m_{45} & m_{46} \\
m_{53} & m_{54} & m_{55} & m_{56} \\
m_{63} & m_{64} & m_{65} & m_{66}
\end{matrix}
\right)
Then the matrix $ M $ is expressed as:
M =
\left(
\begin{matrix}
M_{11} & M_{12} \\
M_{21} & M_{22}
\end{matrix}
\right)
Now, consider the $ 6 \ times 6 $ matrix $ N $, and divide it in the same way as $ M $.
N =
\left(
\begin{matrix}
N_{11} & N_{12} \\
N_{21} & N_{22}
\end{matrix}
\right)
And the product of $ M $ and $ N $ can actually be expressed as:
\begin{align}
MN &=
\left(
\begin{matrix}
M_{11} & M_{12} \\
M_{21} & M_{22}
\end{matrix}
\right)
\left(
\begin{matrix}
N_{11} & N_{12} \\
N_{21} & N_{22}
\end{matrix}
\right) \\
&=
\left(
\begin{matrix}
N_{11}N_{11} + N_{12}N_{21} & N_{11}N_{12} + N_{12}N_{22} \\
N_{21}N_{11} + N_{22}N_{21} & N_{21}N_{12} + N_{22}N_{22}
\end{matrix}
\right) \\
\end{align}
In other words, you can treat the submatrix like a matrix element. (You can see it by actually calculating it.) However, each submatrix must be divided into $ M $ and $ N $ so that it can be calculated consistently.
The introduction has been lengthened, but here are the steps to finally find the matrices $ L $ and $ U $. First, divide $ L $ and $ U $ into submatrixes as shown below. Where $ O $ is a matrix with all elements $ 0 $.
L =
\left(
\begin{matrix}
1 & O \\
L_{21} & L_{22}
\end{matrix}
\right), \
U =
\left(
\begin{matrix}
u_{11} & U_{12} \\
O & U_{22}
\end{matrix}
\right)
$ L_ {22} and U_ {22} $ are the lower and upper triangular matrices of $ (n-1) \ times (n-1) $, respectively. Also, divide the coefficient matrix $ A $ in the same way.
A =
\left(
\begin{matrix}
a_{11} & A_{12} \\
A_{21} & A_{22}
\end{matrix}
\right)
Since $ A = LU $,
\begin{align}
\left(
\begin{matrix}
a_{11} & A_{12} \\
A_{21} & A_{22}
\end{matrix}
\right)
&=
\left(
\begin{matrix}
1 & O \\
L_{21} & L_{22}
\end{matrix}
\right)
\left(
\begin{matrix}
u_{11} & U_{12} \\
O & U_{22}
\end{matrix}
\right) \\
&=
\left(
\begin{matrix}
u_{11} & U_{12} \\
L_{21}u_{11} & L_{21}U_{12} + L_{22}U_{22}
\end{matrix}
\right)
\end{align}
Is expressed as. From this relationship
u_{11} = a_{11}, \ U_{12} = A_{12}
To get Here, assuming $ a_ {11} \ neq 0 $,
L_{21} = \frac{A_{21}}{a_{11}}
And $ L_ {21} $ are also determined. This will give you the first column of $ L $ and the first row of $ U $. Well, I think it's probably easy to understand so far. The question is how to handle $ A_ {22} = L_ {21} U_ {12} + L_ {22} U_ {22} $. Here, as mentioned earlier, $ L_ {21} $ and $ U_ {12} $ have already been obtained, so the first term on the right side has already been obtained. Next, regarding the second term on the right side, how to find it As mentioned at the beginning of the main subject, $ L_ {22} $ is the lower triangular matrix and $ U_ {22} $ is the upper triangular matrix. Therefore, by transposing the first term,
\hat{A}_{22} = A_{22} - L_{21}U_{12}
given that,
A^\prime_{22} = L_{22}U_{22}
Can be expressed as. This means that the $ (n-1) \ times (n-1) $ matrix $ A ^ \ prime_ {22} $ is decomposed into a lower triangular matrix and an upper triangular matrix. Therefore, by applying the same procedure as above for this $ A ^ \ prime_ {22} $, the problem of LU decomposition of the matrix of $ (n-2) \ times (n-2) $ is reduced. In addition, it can be solved sequentially, such as the problem of LU decomposition of the matrix of the $ (n-3) \ times (n-3) $ matrix. Once you understand this, all you have to do is put it in your code.
Now, I will write the code to actually find the matrices $ L $ and $ U $, but if I don't think about anything, a two-dimensional array for the coefficient matrix $ A $ and a two-dimensional array for the lower triangular matrix $ L $ , You might want to have three variables in the two-dimensional array for the upper triangular matrix $ U $. In fact, in an article that introduces LU decomposition (for example, I wrote the LU decomposition method programmatically), I defined variables for each matrix. doing. However, although it is a destructive method, if you take the method of overwriting $ A $, only one variable is required. Specifically, consider using the fact that the diagonal component of $ L $ is $ 1 $ and managing $ L $ and $ U $ in one matrix $ V $ as follows.
L =
\left(
\begin{matrix}
1 & 0 & 0 \\
l_{11} & 1 & 0 \\
l_{21} & l_{22} & 1 \\
\end{matrix}
\right), \
U =
\left(
\begin{matrix}
u_{11} & u_{12} & u_{13} \\
0 & u_{22} & u_{22} \\
0 & 0 & u_{33} \\
\end{matrix}
\right) \\
V =
\left(
\begin{matrix}
u_{11} & u_{12} & u_{13} \\
l_{11} & u_{22} & u_{22} \\
l_{21} & l_{22} & u_{33} \\
\end{matrix}
\right)
By considering such $ V $ and overwriting $ V $ over $ A $, the memory required to save the array is only $ A $, and memory can be saved.
The introduction has become quite a volume, but I will introduce an implementation example in Python. As I've explained at length, the implementation itself is surprisingly easy.
import numpy as np
def main():
n = int(input())
A = np.array([np.array(list(map(float, input().split()))) for i in range(n)])
for i in range(n):
for j in range(i + 1, n):
A[j][i] /= A[i][i]
for k in range(i + 1, n):
A[j][k] -= A[j][i] * A[i][k]
return A
if __name__ == "__main__":
print(main())
Commentary will be added. I will put only the code.
import numpy as np
def main():
n = int(input())
A = np.array([np.array(list(map(float, input().split()))) for i in range(n)])
P = np.identity(n)
for i in range(n):
maxidx = np.argmax(abs(A[i:n, i])) + i
A[i, :], A[maxidx, :] = A[maxidx, :], A[i, :].copy()
P[i, :], P[maxidx, :] = P[maxidx, :], P[i, :].copy()
for j in range(i + 1, n):
A[j][i] /= A[i][i]
for k in range(i + 1, n):
A[j][k] -= A[j][i] * A[i][k]
return A, P
if __name__ == "__main__":
A, P = main()
In fact, the calculation may not be stable as it is. That's because we're assuming $ a_ {ii} \ neq 0 $. To work around this issue, you need to do a ** partial pivot selection ** that swaps the lines. I will add this explanation when I have time.
Numerical calculation for engineering
If you have any questions, please ask! I'm not very familiar with mathematics, so I may not be able to deal with it.
Recommended Posts