python - 在 sympy 中求解正规方程

Question

我正在尝试学习 SymPy，我想弄清楚如何完成一项很酷的任务，派生 normal equations象征性地用于最小二乘问题。

from sympy import *
init_session()

x, y, b = Matrix(), Matrix(), Matrix()
sqNorm = (y - x*b).dot(y- x*b)
solve(diff(sqNorm, b), b)

当我运行它时，我得到

Traceback (most recent call last):
  File "<console>", line 1, in <module>
  File "/usr/local/lib/python2.7/dist-packages/sympy/core/function.py", line 1638, in diff
    return Derivative(f, *symbols, **kwargs)
  File "/usr/local/lib/python2.7/dist-packages/sympy/core/function.py", line 1005, in __new__
    if v._diff_wrt:
  File "/usr/local/lib/python2.7/dist-packages/sympy/matrices/matrices.py", line 3084, in __getattr__
    "%s has no attribute %s." % (self.__class__.__name__, attr))
AttributeError: ImmutableMatrix has no attribute _diff_wrt.

我希望得到类似 (x'x)^{-1}x'y 之类的结果。这在 SymPy 中可能吗？

Answer 1

不，SymPy 没有内置这种级别的抽象矩阵微积分。为了能够区分矩阵，它们必须具有特定的大小并填充可以区分的东西:我在下面举了一个例子。也就是说，您可能对 this project 感兴趣，它通过公理化某些规则在 SymPy 中实现抽象矩阵微积分的元素。

这是您可以使用 SymPy 在符号最小二乘法方面实际执行的示例。可以使用 symarray(使用 a_0_0 表示法)用符号变量填充矩阵。然后计算残差、微分和求解。

from sympy import *
m = 3
n = 2
x = Matrix(symarray('x', (m, n)))
y = Matrix(symarray('y', (m, 1)))
b = Matrix(symarray('b', (n, 1)))
z = (y-x*b).dot(y-x*b)
vars = [b[i,0] for i in range(n)]
eqs = [z.diff(b) for b in vars]
print solve(eqs, vars)

输出虽然正确，但并不具有启发性:

{b_0_0: ((x_0_1**2 + x_1_1**2 + x_2_1**2)*(x_0_0*y_0_0 + x_1_0*y_1_0 + x_2_0*y_2_0) - (x_0_0*x_0_1 + x_1_0*x_1_1 + x_2_0*x_2_1)*(x_0_1*y_0_0 + x_1_1*y_1_0 + x_2_1*y_2_0))/((x_0_0**2 + x_1_0**2 + x_2_0**2)*(x_0_1**2 + x_1_1**2 + x_2_1**2) - (x_0_0*x_0_1 + x_1_0*x_1_1 + x_2_0*x_2_1)**2), b_1_0: ((x_0_0**2 + x_1_0**2 + x_2_0**2)*(x_0_1*y_0_0 + x_1_1*y_1_0 + x_2_1*y_2_0) - (x_0_0*x_0_1 + x_1_0*x_1_1 + x_2_0*x_2_1)*(x_0_0*y_0_0 + x_1_0*y_1_0 + x_2_0*y_2_0))/((x_0_0**2 + x_1_0**2 + x_2_0**2)*(x_0_1**2 + x_1_1**2 + x_2_1**2) - (x_0_0*x_0_1 + x_1_0*x_1_1 + x_2_0*x_2_1)**2)}