LU分解大法好！

LU分解

之前提过的

Randomized SVD

的实现中用到了LU factorization, 所谓LU factorization就是把矩阵分解为一个下三角矩阵 LL 和一个上三角矩阵 UU 的乘积。

高斯消元法

先手动算一发高斯消元：

A=⎛⎝⎜⎜⎜⎜13−1−3−2−92−6−20426−3−972⎞⎠⎟⎟⎟⎟A=(1−2−2−33−90−9−1247−3−6262)

对 AA 做高斯消元法得到上三角矩阵 UU 和下三角矩阵 LL :

LU=⎡⎣⎢⎢⎢⎢13−1−30104001−20001⎤⎦⎥⎥⎥⎥⋅⎡⎣⎢⎢⎢⎢1000−2−300−2620−3041⎤⎦⎥⎥⎥⎥LU=[10003100−1010−34−21]⋅[1−2−2−30−36000240001]

Gaussian Elimination transform a linear system into an upper triangular one by applying linear transformations on the left. It is triangular triangularization.

换个角度看：将 AA 做高斯消元得到 UU 可以看成是对 AA 连续做基础行变换。

Lm−1…L2L1A=ULm−1…L2L1A=U

LL is unit lower-triangular: LL 的对角线上的元素全为1！

def LU(A):
U = np.copy(A)
m, n = A.shape
# 注意L一开始是Identity
L = np.eye(n)

# 最后一列(n)不用处理
for k in range(n-1):
# 第1行不用处理
for j in range(k+1,n):
L[j,k] = U[j,k]/U[k,k]
U[j,k:n] -= L[j,k] * U[k,k:n]
# print(U)会发现U就是在做高斯消元
print(U)
return L, U

高斯消元法复杂度

Work for Gaussian Elimination: 2⋅13n32⋅13n3.

Memory

Above, we created two new matrices, LL and UU. However, we can store the values of LL and UU in our matrix A (overwriting the original matrix). Since the diagonal of LL is all 11s, it doesn’t need to be stored. Doing factorizations or computations in-place is a common technique in numerical linear algebra to save memory.

Note: you wouldn’t want to do this if you needed to use your original matrix AA again in the future.

LU分解的用处

分解为 LL 和 UU 之后，先解 Ly=bLy=b 再解 Ux=yUx=y 会变得简单许多。

Solving Ax=bAx=b becomes LUx=bLUx=b :

1. find A=LUA=LU

2. solve Ly=bLy=b

3. solve Ux=yUx=y

Speed up LU Factorization

Parallelized LU Decomposition LU decomposition can be fully parallelized

Randomized LU Decomposition (2016 article): The randomized LU is fully implemented to run on a standard GPU card without any GPU–CPU data transfer.