LU разложение

  • главные миноры $\neq 0$

Тогда $A = LU$

  • $L$ - нижняя треугольная
  • $U$ - верхняя треугольная

Вычисление: https://courses.physics.illinois.edu/cs357/sp2020/notes/ref-9-linsys.html

$$ \begin{bmatrix}
a_{11} & a_{12}\
a_{21} & A_{22} \end{bmatrix}

\begin{bmatrix}
1 & 0\
l_{21} & L_{22} \end{bmatrix} \times \begin{bmatrix}
u_{11} & u_{12}\
0 & U_{22} \end{bmatrix}

\begin{bmatrix}
u_{11} & u_{12}\
u_{11} l_{21} & (l_{21} u_{12} + L_{22} U_{22}) \end{bmatrix} $$

Отсюда

$$ \begin{align}
u_{11} &= a_{11}\ u_{12} &= a_{12}\ l_{21} &= \frac{1}{u_{11}}a_{21}\ L_{22}U_{22} &= A_{22} - a_{21} \frac{1}{a_{11}} a_{12} = A_{22} - u_{12} l_{21} \end{align} $$

import numpy as np
def lu_decomp(A):
    """(L, U) = lu_decomp(A) is the LU decomposition A = L U
       A is any matrix
       L will be a lower-triangular matrix with 1 on the diagonal, the same shape as A
       U will be an upper-triangular matrix, the same shape as A
    """
    n = A.shape[0]
    if n == 1:
        L = np.array([[1]])
        U = A.copy()
        return (L, U)

    A11 = A[0,0]
    A12 = A[0,1:]
    A21 = A[1:,0]
    A22 = A[1:,1:]

    L11 = 1
    U11 = A11

    L12 = np.zeros(n-1)
    U12 = A12.copy()

    L21 = A21.copy() / U11
    U21 = np.zeros(n-1)

    S22 = A22 - np.outer(L21, U12)
    (L22, U22) = lu_decomp(S22)

    L = np.block([[L11, L12], [L21, L22]])
    U = np.block([[U11, U12], [U21, U22]])
    return (L, U)