cholesky分解计算多项式

百度：https://baike.baidu.com/item/cholesky%E5%88%86%E8%A7%A3/10640395?fr=aladdin

知乎：https://zhuanlan.zhihu.com/p/112091443

csdn应用及理论：https://blog.csdn.net/qq_30142403/article/details/80482581

性能测试：https://www.5axxw.com/questions/content/cegbf2

代码参考：https://vimsky.com/examples/usage/python-scipy.linalg.cho_factor.html

官方说明：https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.diag.html

Cholesky 分解是把一个对称正定的矩阵表示成一个下三角矩阵L和其转置的乘积的分解。它要求矩阵的所有特征值必须大于零，故分解的下三角的对角元也是大于零的。Cholesky分解法又称平方根法，是当A为实对称正定矩阵时，LU三角分解法的变形。

Cholesky分解主要被用于求解Ax=b。。如果A是对称正定的，我们可以先求出A=LL'，随后通过后向替代法 Ly=b求得y，再通过前向替代法L'x=y求得x。实际应用中，Cholesky分解及其LDL变形在求解线性方程组中的效率约两倍于LU分解。

基本思想还是高斯消元：

要求$Ax=b$

一般来说是使用高斯消元法将A转化为行阶梯矩阵，然后就可以快速解出整个多元函数方程。

但是如果A这个矩阵比较特殊，是一个实对称矩阵，那么我们可以先使用cholesky分解将矩阵A分解为一个下三角矩阵及其转置，即:

$$ A = L*L^T $$

那么此时，我们就不需要使用复杂的高斯消元来计算A的行阶梯矩阵了，而是直接得到了行阶梯矩阵：

$$ A = L*L^T $$

$$ Ax=b $$

$$ L*L^Tx=b $$

$$ L*(L^Tx)=b $$

令

$$ y=(L^Tx) $$

$$ Ly=b $$

此时只需要接两个简单的多元函数即可

python的代码实现如下：

但是有的时候还不如直接求逆矩阵

from scipy.linalg import cho_factor, cho_solve
import numpy as np
import time

A = np.array([[9, 3, 1, 5], [3, 7, 5, 1], [1, 5, 9, 2], [5, 1, 2, 6]])

b = np.array([1,1,1,1])

def cholesky_calc(A,b):
    c, low = cho_factor(A)
    x = cho_solve((c, low), b)
    return x

def normal_inv_calc(A,b):
    x = np.dot(np.linalg.inv(A), b)
    return x

if __name__ == '__main__':
    x = cholesky_calc(A,b)
    print(x)
    x = normal_inv_calc(A,b)
    print(x)

    n = 1000000
    start_time = time.time()
    for i in range(n):
        x = cholesky_calc(A, b)
    end_time = time.time()

    cost_time = end_time-start_time
    print("cost_time: {:.2f}s".format(cost_time))

    start_time = time.time()
    for i in range(n):
        x = normal_inv_calc(A, b)
    end_time = time.time()

    cost_time = end_time - start_time
    print("cost_time: {:.2f}s".format(cost_time))