4

I'm trying to learn SymPy, and I'd like to figure out how to do a cool task, deriving the normal equations for a least-squares problem symbolically.

from sympy import *
init_session()

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

When I run that, I get

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.

I am hoping for a result like (x'x)^{-1}x'y or something. Is that possible in SymPy?

Hatshepsut
  • 5,962
  • 8
  • 44
  • 80
  • `Matrix()` is the [empty matrix](https://en.wikipedia.org/wiki/Matrix_(mathematics)#Empty_matrices), which is not what you want. You should be able to use MatrixSymbol, but it looks like it doesn't support differentiation yet. – asmeurer Jan 13 '16 at 15:30
  • See also https://github.com/sympy/sympy/issues/5858 – asmeurer Jan 13 '16 at 15:31
  • Might be related: https://stackoverflow.com/questions/21166771/is-there-a-vectorized-way-to-calculate-the-gradient-in-sympy?rq=1 – Rufus Sep 26 '19 at 01:07

1 Answers1

1

No, SymPy doesn't come with this level of abstract matrix calculus built-in. To be able to differentiate matrices, they have to have specific size and be filled with things you can differentiate in: I give an example below. That said, you may be interested in this project, which implements elements of abstract matrix calculus in SymPy by axiomatizing certain rules.

Here is a sample of you can actually do in terms of symbolic least squares with SymPy. Filling matrices with symbolic variables can be done with symarray (which uses a_0_0 notation). Then compute the residuals, differentiate, and solve. 

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)

The output, though correct, is not enlightening: 

{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)}