Python sympy做代数运算解Cholesky求逆的L和L逆矩阵示例

Posted 2022-12-08 飞凡可期

提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档

文章结构

• sympy包问题
• 一、代码
• 二、运行结果
• 总结

---

sympy包问题

在1.6.x的包中发现用“from sympy import *"遇上频繁的Warning问题？
解决方法：shell power >> pip install -U sympy 更新到了1.9系列，问题解决！

---

提示：sympy不仅仅有symbol系列的代数运算功能，其Matrix系列的矩阵数值运算也比较方便，比如Matrix（[[1,2],[3,4]]）这样对于习惯matlab运算的非常友好。

一、代码

"""
利用cholesky进行矩阵求逆，这里代数演算求解元素
"""
from sympy import *
import numpy as np

#1 alg ex1
x = Symbol('x')
y = Symbol('y')
a, b, c = symbols('a b c')
print(3*x+b*y**2)

#2 解方程组
expr = a*x -b
expr2 = b*y -c
s_expr = solve([expr, expr2], [x,y])
print(s_expr)

#3 matrix comb:利用numpy轻松做行、列list拼接
m = [1,2,3,4,5,6]
m2 = [7,8,9,10,11,12]
Mr = np.r_[m,m2]
Mc = np.c_[m,m2].T
print(Mr, " \\n mc = ", Mc)

#np.r_[]， np.c_[]注意是索引不是函数调用，拼接之后变为nparray矩阵 不再是list；
"""
求二元一次方程组：
2x + 9y = 81
3x + y = 34
"""
Coeff = Matrix([[2,9],[3,1]])
x,y = symbols('x y')
p_expr = Coeff * Matrix([x,y]) - Matrix([[81],[34]])
s_expr = solve(p_expr, [x,y])
print(s_expr)





####### 下面是cholesky分解求逆算法，一种针对共轭对称矩阵的快速、低复杂度求逆计算方法 ########
# A = L*D*L^H; L是下三角、D是对角阵；   A逆 = L^-H * D^-1 * L^-1
def printVar(exp_a, DX, x_exp):
    """
    :param exp_a: 公式矩阵
    :param DX: 未知变量
    :param x_exp: 求解结果[rst], dict
    :return:
    """
    print("Formula:\\n", exp_a, "\\n unKnow var:\\n", DX)
    print("Rst:\\n", x_exp)
    if type(x_exp) == list:
        for x, y in zip(DX, x_exp[0]):
            print(x, " = ", y)
    elif type(x_exp) == dict:
        for x, y in x_exp.items():
            print(x, " = ", y)

#P1 A = L * D * L.T, find L and D
#print(Matrix(a,b,c)) 2x2 矩阵求逆，代数解
a = symbols('a(0:2)(0:2)')
A = Matrix([[a[0], a[2]], [a[2], a[3]]]) # hermit sym matrix
d = [Symbol('d0'), Symbol('d1')]
D = Matrix([[d[0],0],[0, d[1]]])
l = [Symbol('l10')]
L = Matrix([[1,0],[l[0],1] ,])
exp_a = L*D*L.T - A
DX = d.copy()
DX.extend(l)
x_exp = solve(list(exp_a), DX)
print("L D result: ")
printVar(exp_a, DX, x_exp)

    #3-1 L*G = I,即求下三角矩阵L的逆
g = [symbols('g10')]
G = Matrix([[1,0], [g[0], 1] ])
lg_exp = L*G - eye(2)
g_exp = solve(list(lg_exp),g)
print("L^-1 result: ")
printVar(lg_exp, g, g_exp) #[list(g_exp.values())])


#P2 A = L * D * L.T, find L and D
#print(Matrix(a,b,c)) 2x2 矩阵求逆，代数解
a = symbols('a(0:4)(0:4)')
A = Matrix([[a[0], a[4], a[8], a[12]], [a[4], a[5], a[9], a[13],],
            [a[8], a[9], a[10], a[14], ], [a[12], a[13], a[14], a[15]]]) # hermit sym matrix
d = [Symbol('d0'), Symbol('d1'), Symbol('d2'), Symbol('d3')]
D = Matrix([[d[0],0, 0, 0],[0, d[1], 0, 0], [0, 0, d[2], 0], [0, 0, 0, d[3]]])
l = symbols('l10 l20 l21 l30 l31 l32')
L = Matrix([[1, 0, 0, 0],[l[0],1, 0,0] ,[l[1], l[2], 1, 0], [l[3], l[4], l[5], 1]])
exp_a = L*D*L.T - A
DX = d.copy()
DX.extend(l)
x_exp = solve(list(exp_a), DX)
print("L D result: ")
printVar(exp_a, DX, x_exp)

    #3-1 L*G = I,即求下三角矩阵L的逆
g =symbols('g10 g20 g21 g30 g31 g32')
G = Matrix([[1, 0, 0, 0], [g[0], 1, 0, 0], [g[1], g[2], 1, 0], [g[3], g[4], g[5], 1]])
lg_exp = L*G - eye(4)
g_exp = solve(list(lg_exp),g)
print("L^-1 result: ")
printVar(lg_exp, g, g_exp) #[list(g_exp.values())])

二、运行结果

by**2 + 3x
x: b/a, y: c/b
[ 1 2 3 4 5 6 7 8 9 10 11 12]
mc = [[ 1 2 3 4 5 6]
[ 7 8 9 10 11 12]]
x: 9, y: 7
L D result:
Formula:
Matrix([[-a00 + d0, -a10 + d0l10], [-a10 + d0l10, -a11 + d0l10**2 + d1]])
unKnow var:
[d0, d1, l10]
Rst:
[(a00, (a00a11 - a102)/a00, a10/a00)]
d0 = a00
d1 = (a00*a11 - a102)/a00
l10 = a10/a00
L^-1 result:
Formula:
Matrix([[0, 0], [g10 + l10, 0]])
unKnow var:
[g10]
Rst:
g10: -l10
g10 = -l10
L D result:
Formula:
Matrix([[-a00 + d0, -a10 + d0l10, -a20 + d0l20, -a30 + d0l30], [-a10 + d0l10, -a11 + d0l10**2 + d1, -a21 + d0l10l20 + d1l21, -a31 + d0l10l30 + d1l31], [-a20 + d0l20, -a21 + d0l10l20 + d1l21, -a22 + d0l202 + d1*l212 + d2, -a32 + d0l20l30 + d1l21l31 + d2l32], [-a30 + d0l30, -a31 + d0l10l30 + d1l31, -a32 + d0l20l30 + d1l21l31 + d2l32, -a33 + d0l30**2 + d1l312 + d2*l322 + d3]])
unKnow var:
[d0, d1, d2, d3, l10, l20, l21, l30, l31, l32]
Rst:
[(a00, (a00a11 - a10**2)/a00, (a00a11a22 - a00a212 - a102a22 + 2a10a20a21 - a11a20**2)/(a00a11 - a102), (a00a11a22a33 - a00a11*a322 - a00a21**2a33 + 2a00a21a31a32 - a00a22a312 - a102a22a33 + a102*a322 + 2a10a20a21a33 - 2a10a20a31a32 - 2a10a21a30a32 + 2a10a22a30a31 - a11a20**2a33 + 2a11a20a30a32 - a11a22a302 + a202a31**2 - 2a20a21a30a31 + a21**2a302)/(a00a11a22 - a00*a212 - a102a22 + 2a10a20a21 - a11*a202), a10/a00, a20/a00, (a00a21 - a10a20)/(a00a11 - a10**2), a30/a00, (a00a31 - a10a30)/(a00a11 - a102), (a00a11a32 - a00a21a31 - a102a32 + a10a20a31 + a10a21a30 - a11a20a30)/(a00a11a22 - a00a212 - a102a22 + 2a10a20a21 - a11a20**2))]
d0 = a00
d1 = (a00a11 - a102)/a00
d2 = (a00a11a22 - a00*a212 - a102a22 + 2a10a20a21 - a11*a202)/(a00a11 - a10**2)
d3 = (a00a11a22a33 - a00a11a322 - a00*a212a33 + 2a00a21a31a32 - a00a22a312 - a102a22a33 + a10**2a322 + 2a10a20a21a33 - 2a10a20a31a32 - 2a10a21a30a32 + 2a10a22a30a31 - a11*a202a33 + 2a11a20a30a32 - a11a22a302 + a202a312 - 2a20a21a30a31 + a212a30**2)/(a00a11a22 - a00a212 - a102a22 + 2a10a20a21 - a11a20**2)
l10 = a10/a00
l20 = a20/a00
l21 = (a00a21 - a10a20)/(a00a11 - a102)
l30 = a30/a00
l31 = (a00a31 - a10a30)/(a00*a11 - a102)
l32 = (a00a11a32 - a00a21a31 - a102a32 + a10a20a31 + a10a21a30 - a11a20a30)/(a00a11a22 - a00a212 - a102a22 + 2a10a20a21 - a11*a202)
L^-1 result:
Formula:
Matrix([[0, 0, 0, 0], [g10 + l10, 0, 0, 0], [g10l21 + g20 + l20, g21 + l21, 0, 0], [g10l31 + g20l32 + g30 + l30, g21l32 + g31 + l31, g32 + l32, 0]])
unKnow var:
(g10, g20, g21, g30, g31, g32)
Rst:
g10: -l10, g20: l10l21 - l20, g30: -l10l21l32 + l10l31 + l20l32 - l30, g21: -l21, g31: l21l32 - l31, g32: -l32
g10 = -l10
g20 = l10l21 - l20
g30 = -l10l21l32 + l10l31 + l20l32 - l30
g21 = -l21
g31 = l21l32 - l31
g32 = -l32
s_expr
Out[3]: x: 9, y: 7

总结

sympy非常强大，好用，对于小众的大部分解方程、代数演算、矩阵代数验算和多元方程联合求解都够用了！

• 熟练使用var = Symbol(‘v’)， arr = symbols(’(a b c)’) ， Mat = Matrix([var[0], var[1]])这几种代数定义和数据结构组织方式非常管用，类似list，array，dict在基本python的意义。
• 熟悉rst_exp = solve(program_exp, unknown varibale list）这个解代数的语句，熟悉其返回的dict格式，以及其它如可变格式（返回列表list，内嵌差异化格式）的差异。