SymPy中矩阵的逆？

Question

我在想，怎么用Python里的SymPy库来创建一个矩阵并计算它的逆矩阵呢？

比如说，对于这个符号矩阵：

Answer 1

对于符号矩阵，sympy提供的方法速度非常慢。逐个元素的操作要快得多。我在一篇名为“部分逆”的论文中发现了这个问题：doi.org/10.1088/1402-4896/ad298a。

这里有一些代码和关于(2,2)和(3,3)维度符号矩阵的时间性能图表。（.inv()这个方法在我尝试更高维度时无法完成编译，而“部分逆”这个方法的速度相对较快。）

from sympy import Matrix, Symbol, simplify

def sp_partial_inversion(m, *cells):
    ''' Partial inversion algorithm.
    ARGUMENTS
        m <sympy.Matrix> : symbolic matrix
        *cells <tuple>   : 2-tuples with matrix indices to perform partial inversion on.
    RETURNS
        <sympy.Matrix>   : matrix m partially-inverted at indices *cells
    '''
    # Partial inversion algorithm
    M = m.copy()
    for cell in cells:
        i,k = cell
        z = M[i,k]
        newmat = []
        for r in range(m.shape[0]):
            newrow = []
            for c in range(m.shape[1]):
                if i == r:
                    if k == c:  # Pivot cell
                        newrow.append( 1/z )
                    else:       # Pivot row
                        newrow.append( -M[r,c]/z )
                else:
                    if k == c:  # Pivot col
                        newrow.append( M[r,c]/z )
                    else:       # Elsewhere
                        newrow.append( M[r,c] - M[i,c]*M[r,k]/z )
            newmat.append(newrow)
        M =  Matrix(newmat)
    #
    return M

def SymbolicMatrix(n):
    "Generates symbolic matrix for tests"
    return Matrix( [ Symbol( f'm_{i}' ) for i in range(n**2) ] ).reshape(n,n)

def FullInversion(m):
    "Full matrix inversion is partial inversion at all i==j."
    cells = [(i,i) for i in range(m.shape[0])]
    return sp_partial_inversion(m, *cells)

# Test 1 --- Z should be matrix of zeros
M = SymbolicMatrix(3)
#Z = simplify( FullInversion(M) - M.inv() )

# Test 2 ---
M = SymbolicMatrix(4)
M_ = simplify( FullInversion(M) )
M_

看看这些图表：

这是关于1000个符号矩阵从(2,2)到(3,3)的时间性能图。

这是关于100个符号矩阵从(2,2)到(4,4)的时间性能图。

这个方法也可以用于数值矩阵，但在我的测试中，它的速度并没有比numpy默认的矩阵求逆快。

Answer 2

这里有一个例子，展示了我们如何计算一个符号矩阵的逆矩阵（用问题中的那个矩阵）：

import sympy as sym


# Not necessary but gives nice-looking latex output
# More info at: http://docs.sympy.org/latest/tutorial/printing.html
sym.init_printing()

sx, sy, rho = sym.symbols('sigma_x sigma_y rho')
matrix = sym.Matrix([[sx ** 2, rho * sx * sy], 
                     [rho * sx * sy, sy ** 2]])

现在打印逆矩阵 matrix.inv() 会得到：

这个结果可以进一步简化，使用 sym.simplify(matrix.inv())：

Answer 3

如果你的问题是：如何在sympy中计算矩阵M的逆，那么：

M_inverse = M.inv()

关于如何创建一个矩阵：

M = Matrix(2,3, [1,2,3,4,5,6])

这将给你一个如下的2行3列的矩阵：

1 2 3

4 5 6

查看详细信息：http://docs.sympy.org/0.7.2/modules/matrices/matrices.html

SymPy中矩阵的逆？

3 个回答

撰写回答