如何解决Sympy 中有没有办法简化矩阵代数?用伍德伯里身份说
我想指定一个 N x M 矩阵,其中 N >> M 并使用伍德伯里, 并让 Sympy 简化矩阵代数。
这是可能的还是类似的?
解决方法
Sympy 中有没有办法简化矩阵代数?用伍德伯里身份说
是的,可以应用 Woodbury's identity。在下面的代码中,标识应用于 (I + Phi^T * Phi)**(-1)
。
使用的特殊情况是 (I+UV)^-1 = I - U*(I+VU)^-1*V
from sympy.simplify.simplify import nc_simplify
from sympy import *
N,M,sigma = symbols("N M \sigma")
Phi = MatrixSymbol("Phi",N,M)
In,Im = Identity(N),Identity(M)
f = (Im + Phi.transpose()@Phi).inverse()
f = f @ Phi.transpose()
f = Phi @ f
f = In - sigma**(-2)*f
f = sigma**(-2)*Phi.transpose()@f
f = f@ Phi
def apply_woodburry_1(e,U,V):
# The Woodbury identity
#
# (I+UV)^-1 = I - U * (I+VU)^-1 * V
#
# Below,P will be LHS and Q will be RHS
#
UV_dim = (U * V).shape
I1 = Identity(UV_dim[0])
VU_dim = (V * U).shape
I2 = Identity(VU_dim[0])
P = (I1+U*V)**(-1)
Q = I1 - U * ( (I2+V*U)**(-1) ) * V
return e.replace(P,Q)
display(f)
f = apply_woodburry_1(f,Phi.transpose(),Phi)
display(f)
f = nc_simplify(f.expand())
display(f)
输出:
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。