LU decomposition in Python

Introduction

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

Triangular matrix

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.

LU decomposition

Overview

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.)

How to find

Now, I will explain how to find it concretely, but before that, I will explain the submatrix.

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.

How to find the main subject

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.

Implementation

Memory-saving implementation

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.

Implementation subject

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())

Partial pivot selection

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()

Supplement

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.

reference

Numerical calculation for engineering

Finally

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

LU decomposition in Python
Introduction to Linear Algebra in Python: A = LU Decomposition
Quadtree in Python --2
CURL in python
Metaprogramming in Python
Python 3.3 in Anaconda
Geocoding in python
SendKeys in Python
Meta-analysis in Python
Unittest in python
Epoch in Python
Discord in Python
Sudoku in Python
DCI in Python
quicksort in python
nCr in python
N-Gram in Python
Programming in python
Plink in Python
Constant in python
Lifegame in Python.
FizzBuzz in Python
Sqlite in python
StepAIC in Python
N-gram in python
LINE-Bot [0] in Python
Csv in python
Disassemble in Python
Reflection in Python
Constant in python
nCr in Python.
format in python
Scons in Python3
Puyo Puyo in python
python in virtualenv
PPAP in Python
Quad-tree in Python
Reflection in Python
Chemistry in Python
Hashable in python
DirectLiNGAM in Python
LiNGAM in Python
Flatten in python
flatten in python
Sorted list in Python
Daily AtCoder # 36 in Python
Clustering text in Python
Daily AtCoder # 2 in Python
Implement Enigma in python
Daily AtCoder # 32 in Python
Daily AtCoder # 6 in Python
Edit fonts in Python
Singleton pattern in Python
File operations in Python
Read DXF in python
Daily AtCoder # 53 in Python
Key input in Python
Use config.ini in Python
Daily AtCoder # 33 in Python
Solve ABC168D in Python
Logistic distribution in Python