使用 SymPy 对符号矩阵表达式进行指数化

Question

我正在尝试使用SymPy以符号形式处理多元正态分布。但是，它似乎无法意识到我想要幂的矩阵表达式的计算结果为标量。代码如下：

from sympy import *
# Sample size, number of covariates
n, k = symbols('n k')
# Data
X = MatrixSymbol('X', n, k)
y = MatrixSymbol('y', n, 1)
# Parameters
beta = MatrixSymbol('beta', k, 1)
Sigma = MatrixSymbol('Sigma', n, n)
# Exponent expression, OK:
(-Rational(1, 2)*(y - X*beta).T*Sigma*(y - X*beta))
# Trying to use as exponent, not OK:
exp(-Rational(1, 2)*(y - X*beta).T*Sigma*(y - X*beta))

错误消息如下：

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-76-62c854b75962> in <module>()
----> 1 exp((-Rational(1, 2)*(y - X*beta).T*Sigma*(y - X*beta)))
D:ProgrammingAnaconda3libsite-packagessympycorefunction.py in __new__(cls, *args, **options)
425 
426         evaluate = options.get('evaluate', global_evaluate[0])
--> 427         result = super(Function, cls).__new__(cls, *args, **options)
428         if not evaluate or not isinstance(result, cls):
429             return result
D:ProgrammingAnaconda3libsite-packagessympycorefunction.py in __new__(cls, *args, **options)
248 
249         if evaluate:
--> 250             evaluated = cls.eval(*args)
251             if evaluated is not None:
252                 return evaluated
D:ProgrammingAnaconda3libsite-packagessympyfunctionselementaryexponential.py in eval(cls, arg)
300 
301         elif arg.is_Matrix:
--> 302             return arg.exp()
303 
304     @property
AttributeError: 'MatMul' object has no attribute 'exp'

Answer 1

指数函数不支持矩阵表达式;它需要一个显式矩阵，这意味着具有明确条目的确定大小的矩阵(可以是符号(。若要从矩阵表达式获取显式矩阵，请使用as_explicit方法。

exp((-Rational(1, 2)*(y - X*beta).T*Sigma*(y - X*beta)).as_explicit())

返回一个 1 x 1 矩阵，即

Matrix([[exp(Sum(-y[_k, 0]*Sum(Sigma[_k, _k]*y[_k, 0], (_k, 0, n - 1)) + y[_k, 0]*Sum(Sigma[_k, _k]*Sum(X[_k, _k]*beta[_k, 0], (_k, 0, k - 1)), (_k, 0, n - 1)) + Sum(Sigma[_k, _k]*y[_k, 0], (_k, 0, n - 1))*Sum(X[_k, _k]*beta[_k, 0], (_k, 0, k - 1)) - Sum(Sigma[_k, _k]*Sum(X[_k, _k]*beta[_k, 0], (_k, 0, k - 1)), (_k, 0, n - 1))*Sum(X[_k, _k]*beta[_k, 0], (_k, 0, k - 1)), (_k, 0, n - 1))/2)]])

或者你可以做

exp((-Rational(1, 2)*(y - X*beta).T*Sigma*(y - X*beta)).as_explicit()[0, 0])

从 1 x 1 矩阵下降到标量。