cholesky 分解物理意义

转载

clghxq 2025-01-02 13:05:36

文章标签 cholesky 分解物理意义 python 线性代数 for循环嵌套 文章分类 架构后端开发

Cholesky分解定理：就是把一个正定的矩阵分解成两个对称的三角阵

$cholesky 分解物理意义_python$

若要 $cholesky 分解物理意义_线性代数_02$ 为单位下三角阵，则有

$cholesky 分解物理意义_线性代数_03$ （改进的平方根法）

cholesky 分解物理意义_线性代数_04

cholesky 分解物理意义_for循环_05

下面是 $cholesky 分解物理意义_线性代数_03$ 的代码实现

import numpy as np
def factory(A):
   
    n = len(A)
    L = np.eye((n))
    D = np.ones(n)
    mid_l = 0
    
    D[0] = A[0][0]
    for i in range(1,n):
        L[i][0] = A[i][0]/D[0]
    
    
    for j in range(1,n-1):
        D[j]=A[j][j]
        for i in range(0,j):
            D[j] = D[j] - (L[j][i]**2)*D[i]
        for i in range(j+1,n):
            for k in range(j):
                mid_l = mid_l + L[i][k]*D[k]*L[j][k] 
                L[i][j] = (A[i][j]-mid_l)/D[j]
                mid_l = 0
    
    D[n-1] = A[n-1][n-1]
    for i in range(n-1):
        D[n-1] = D[n-1] - (L[n-1][i]**2)*D[i]
        
    print(L)
    print(D)

但是很快就发现问题了！

cholesky 分解物理意义_python_07

对比三种算法的运行速度，Cholesky分解居然是最慢的！这显然不符合我们的初衷，于是参考CSDN愚者吃鱼博主写的算法，发现他并没有用很多个for循环嵌套来实现该算法，反而用了numpy库中的dot（点积），这使得该算法变得十分整洁，不像我写的那么累赘

#为方便示意不考虑时间空间代价
def Cholesky_plus(matrix):
    w = matrix.shape[0]
    L = np.zeros((w,w))
    for i in range(w):
        L[i,i] = 1
    D = np.zeros((w,w))
    for i in range(w):
        D[i,i] = matrix[i,i] - np.dot(np.dot(L[i,:i],D[:i,:i]),L[i,:i].T)
        for j in range(i+1,w):
            L[j,i] = (matrix[j,i] - np.dot(np.dot(L[j,:i],D[:i,:i]),L[i,:i].T))/D[i,i]
    return L,D

同时我对比了一下运行速度，发现我的比他的要慢很多。

cholesky 分解物理意义_嵌套_08

cholesky 分解物理意义_线性代数_09

所以得出了一个结论：numpy里面的dot比python原本的for循环更快，并且能使代码更加整洁！

之后修改了自己的代码：

import numpy as np
def Cholesky(A):

    n = len(A)
    L = np.eye((n))
    D = np.ones(n)
    D[0] = L[0,0]
    
    for j in range(n-1):
        for i in range(j+1,n):
            D[j] = A[j,j] - np.dot(L[j,:j]*L[j,:j],D[:j])
            L[i,j] = (A[i,j] - np.dot(L[i,:j]*L[j,:j],D[:j]))/D[j]
            
    D[n-1] = A[n-1,n-1] - np.dot(L[n-1,:n-1]*L[n-1,:n-1],D[:n-1])
    
    return L,D

这样写不仅使我不再因为4个for循环的嵌套而头疼，而且还提升了运行速度。按照这样的思想，我对之前写的Gauss分解和PGuass分解算法进行了改进：

“”“Gauss消元法”
import numpy as np
def Gauss(A):
    n = len(A)
    for j in range(n-1):
            A[j+1:,j] = A[j+1:,j]/A[j,j]
            A[j+1:,j+1:]=A[j+1:,j+1:] - np.dot(A[j,j+1:],A[j+1:,j])

    return A

"""PGauss消元法"""
import numpy as np
def PGauss(A):
    n = len(A)
    for j in range(n-1):
            
        B = np.argmax(A[j:,j:]**2, axis=0)
        k = B[0]+j
        tmp = np.copy(A[k])
        A[k] = A[j]
        A[j] = tmp
            
        if A[j,j] != 0:
            A[j+1:,j] = A[j+1:,j]/A[j,j]
            A[j+1:,j+1:]=A[j+1:,j+1:] - np.dot(A[j,j+1:],A[j+1:,j])
        else:
            print("该矩阵奇异，无法进行LU分解")

    return A

用向量的思想去编程，减少了for的使用，同时使这两个算法的速度也大大提升！

但是，对比三个算法的运行速度，我又迷惑了…

cholesky 分解物理意义_cholesky 分解物理意义_10