Solver freezes on matrix equation

[jbuckman]

I'm trying to use SymPy to symbolically solve some basic matrix equations. Here's an example, using it to compute the closed-form optimal coefficient vector for least-squares linear regression:

from sympy.solvers import solve

# n := number of data points
# d := dimension of input
n, d = symbols('n d', integer=True, positive=True)

# data matrices
X = MatrixSymbol('X', n, d)
y = MatrixSymbol('y', n, 1)
# weight matrix
β = MatrixSymbol('β', d, 1)

errors = y - X*β
loss = errors.T*errors
gradient = diff(loss, β)
ans = solve(gradient, β) # this line freezes

This is a pretty straightforward problem with a solution of (X.T*X)**(-1)*X.T*y. Everything looks fine up until the point that the solver is invoked, at which point it simply freezes up and never returns an answer.

Is this a bug? Is there something I need to do differently to make it go through? Or is this a known limitation of the solver?

[oscarbenjamin]

This shouldn't hang but to be clear solve does not have any capability to handle symbolic matrix equations.

Internally solve tries to convert this to a list and that hangs:

In [35]: list(gradient)
^C

Probably iterating over this should give an error because the shape is symbolic.

[jbuckman]

I see. Is there any plan to implement this functionality at some point in the future?

It also fails when I just define arbitrary non-symbolic dimensions n=5, d=2. It does not hang, it gives an error: TypeError: 'Zero' object is not subscriptable. I'll see whether I can debug this first.

[oscarbenjamin]

I think you might be hitting up against the same issue as #19162 #19217 #20386. There is a PR to fix this but it seems to have stalled on another problem: #20691

[EdwardRaff]

Following up in this thread that the issue still seems to occur, though all the issues you referenced seem to be closed now.

[oscarbenjamin]

The issue causing the freeze is basically this:

In [20]: n = Symbol('n')

In [21]: M = MatrixSymbol('M', n, 1)

In [22]: list(M)
^C---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-22-ccf281550fe7> in <module>
----> 1 list(M)

~/current/sympy/sympy.git/sympy/matrices/expressions/matexpr.py in __getitem__(self, key)
    303             i = key // cols
    304             j = key % cols
--> 305             if self.valid_index(i, j) != False:
    306                 return self._entry(i, j)
    307             else:

The diff here fixes it so that attempting to iterate over a matrix expression with symbolic shape gives TypeError:

diff --git a/sympy/matrices/expressions/matexpr.py b/sympy/matrices/expressions/matexpr.py
index 78db834..1c21b51 100644
--- a/sympy/matrices/expressions/matexpr.py
+++ b/sympy/matrices/expressions/matexpr.py
@@ -316,6 +316,11 @@ def _is_shape_symbolic(self) -> bool:
         return (not isinstance(self.rows, (SYMPY_INTS, Integer))
             or not isinstance(self.cols, (SYMPY_INTS, Integer)))
 
+    def __iter__(self):
+        if self._is_shape_symbolic():
+            raise TypeError("Matrix expression with symbolic shape is not iterable")
+        return (self[i] for i in range(self.rows*self.cols))
+
     def as_explicit(self):
         """
         Returns a dense Matrix with elements represented explicitly

Then the OP example gives:

   ...: ans = solve(gradient, β) # this line freezes
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-1-8b3ed4f39ee9> in <module>
     14 loss = errors.T*errors
     15 gradient = diff(loss, β)
---> 16 ans = solve(gradient, β) # this line freezes

~/current/sympy/sympy.git/sympy/solvers/solvers.py in solve(f, *symbols, **flags)
    920         if f[i].is_Matrix:
    921             bare_f = False
--> 922             f.extend(list(f[i]))
    923             f[i] = S.Zero
    924 

~/current/sympy/sympy.git/sympy/matrices/expressions/matexpr.py in __iter__(self)
    319     def __iter__(self):
    320         if self._is_shape_symbolic():
--> 321             raise TypeError("Matrix expression with symbolic shape is not iterable")
    322         return super().__iter__(self)
    323 

TypeError: Matrix expression with symbolic shape is not iterable