首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >自定义求解器之Cholesky分解法

自定义求解器之Cholesky分解法

作者头像
fem178
发布2021-04-16 15:38:54
1.1K0
发布2021-04-16 15:38:54
举报

对称正定矩阵

A

可以分解为

A=LL^T

,这种分解被称为Cholesky分解,是LU分解的一个重要特例,可以显著降低计算量。在计算机程序中常常用到这种方法解线性代数方程组。它的优点是节省存储量,得到的L矩阵可以覆盖原来的A矩阵。

对于方程组

Ax=b

,可以写成

LL^Tx=b

,令

L^Tx=y

,则

Ly=b

考察一个3X3矩阵:

\begin{bmatrix} A_{11} & A_{12} & A_{13}\\ A_{21} & A_{22} & A_{23}\\ A_{31} & A_{32} & A_{33}\\ \end{bmatrix} = \begin{bmatrix} L_{11} & 0 & 0\\ L_{21} & L_{22} & 0\\ L_{31} & L_{32} & L_{33}\\ \end{bmatrix} \begin{bmatrix} L_{11} & L_{21} & L_{31}\\ 0 & L_{22} & L_{32}\\ 0 & 0 & L_{33}\\ \end{bmatrix}
= \begin{bmatrix} L^2_{11} & L_{11}L_{21} & L_{11}L_{31}\\ L_{11}L_{21} & L^2_{21}+L^2_{22} & L_{21}L_{31}+L_{22}L_{32}\\ L_{11}L_{31} & L_{21}L_{31}+L_{22}L_{32} &L^2_{31}+L^2_{32}+L^2_{33} \end{bmatrix}

A_{11}=L^2_{11},L_{11} = \sqrt {A_{11}}
A_{21}=L_{11}L_{21}, L_{21}=A_{21}/L_{11}
A_{31}=L_{11}L_{31}, L_{31}=A_{31}/L_{11}
A_{22}=L^2_{21}+L^2_{22}, L_{22}=\sqrt{A_{22}-L^2_{21} }
A_{32}=L_{21}L_{31}+L_{22}L_{32}, L_{32}=(A_{32}-L_{21}L_{31})/L_{22}
A_{33}=L^2_{31}+L^2_{32}+L^2_{33}, L_{33}=\sqrt{A_{33}-L^2_{31}-L^2_{32}}
LL^T

分解的算法:

L_{11} = \sqrt {A_{11}},L_{i1} = A_{i1}/L_{11} ,i=2,3,...,n
L_{jj} = \sqrt {A_{jj}-\sum_{k=1}^{j-1}L^2_{jk}} ,j = 2,3,...,n
L_{ij}=(A_{ij}-\sum_{k=1}^{j-1})/L_{ik} L_{jk} ,j=2,3,...,n-1,i=j+1,j+2,...,n
import numpy as np
import math

class LinerSolver:
    def __init__(self, A, b):
        self.A = A
        self.b = b

    def CholeskiSolver(self):
        n = len(self.A)
        # 分解 [A] = [L] [L^T]
        for k in range(n):
            self.A[k,k] = math.sqrt(self.A[k,k] - np.dot(self.A[k,0:k], self.A[k,0:k]))
            for i in range(k+1,n):
                self.A[i,k] = (self.A[i,k] - np.dot(self.A[i,0:k], self.A[k,0:k])) / self.A[k,k]
        for k in range(1,n): 
            self.A[0:k,k] = 0.0
        
        # 求解 [L]{y} = {b} 
        for k in range(n):
            self.b[k] = (self.b[k] - np.dot(self.A[k,0:k], self.b[0:k])) / self.A[k,k]
        # 求解 [L^T]{x} = {y} 
        for k in range(n-1,-1,-1):
            self.b[k] = (self.b[k] - np.dot(self.A[k+1:n,k], self.b[k+1:n])) / self.A[k,k]
        
        return self.b



A = np.array([  [ 4,   -1,     1],
                [-1,  4.25,  2.75],
                [1,   2.75,   3.5] ])



b = np.array([4, 6, 7.25])

cls =  LinerSolver(A, b) #创建一个求解器的实例cls

x = cls.CholeskiSolver() #调用Choleski法求解
print(x)


与高斯消去法相比,LL分解的优点在于,一旦A被分解,我们就可以对任意多个常量向量b求解Ax=b。每个附加解决方案的成本相对较小,因为前向和后向替换操作比分解过程花费的时间少得多.

本文参与 腾讯云自媒体分享计划,分享自微信公众号。
原始发表:2021-04-06,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 数值分析与有限元编程 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体分享计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档