exception while simplifying product over matrices

[albertz]

a, b = symbols("a b")
k = Symbol("k", integer=True)
n = Symbol("n", integer=True)
expr = Product(Matrix([[a, b], [b, a]]), (k, 1, n))

Now expr.simplify() results in the exception TypeError: cannot add <class 'sympy.matrices.immutable.ImmutableDenseMatrix'> and <class 'int'>.

[oscarbenjamin]

Looks like a bug in Product that should be fixed. This line assumes that the expression in the Product is a number rather than a matrix:

sympy/sympy/concrete/products.py

Line 299 in aaf43fc

if (term - 1).is_zero:

For your example SymPy can compute it as a power:

In [48]: M = Matrix([[a, b], [b, a]])                                                                                                          

In [49]: M                                                                                                                                     
Out[49]: 
⎡a  b⎤
⎢    ⎥
⎣b  a⎦

In [50]: M**n                                                                                                                                  
Out[50]: 
⎡        n          n            n          n⎤
⎢ (a - b)    (a + b)      (a - b)    (a + b) ⎥
⎢ ──────── + ────────   - ──────── + ────────⎥
⎢    2          2            2          2    ⎥
⎢                                            ⎥
⎢         n          n          n          n ⎥
⎢  (a - b)    (a + b)    (a - b)    (a + b)  ⎥
⎢- ──────── + ────────   ──────── + ──────── ⎥
⎣     2          2          2          2     ⎦

A side note here might be that perhaps Product should recognise when a product is just a power. I'm not sure what other matrix cases could be handled in general... (for symbolic limits)

[sylee957]

I anticipate that there are lots of problems regarding sympy functional used with matrix, because Basic Expr do not have the necessary attributes like shape.
I'd expect that there should be lots of matrix expression classes should be introduced, to fix these sort of problems, like MatSum, MatProd, MatIntegral, MatSin, MatCos, ...

[oscarbenjamin]

I'm not sure we need special MatSum etc. At some point we will want something else like VecSum and so on. If we have a Cartesian product of things: operators (Add, Mul, ...) and objects (Number, Matrix, Vector, Tensor, Quaternion, ...) then making a class for every pair of the Cartesian product is not a great idea.

I can understand introducing MatAdd and MatMul etc but it is not ultimately a scalable solution to implement everything in that way. A MatSum would just be a Sum that is matrix-valued and has code to handle sums of matrices using applyfunc. We can add the code to handle sums of matrices without a new class so the only thing that MatSum solves is making the objects as a whole appear as a Matrix so we don't get:

In [1]: M = Matrix([[x, y], [y, x]])

In [2]: S = Sum(M, (n, 1, m))

In [3]: S
Out[3]: 
  m         
 ____       
 ╲          
  ╲         
   ╲  ⎡x  y⎤
   ╱  ⎢    ⎥
  ╱   ⎣y  x⎦
 ╱          
 ‾‾‾‾       
n = 1       

In [4]: A = eye(3)

In [5]: A
Out[5]: 
⎡1  0  0⎤
⎢       ⎥
⎢0  1  0⎥
⎢       ⎥
⎣0  0  1⎦

In [6]: S*A
Out[6]: 
⎡  m                                     ⎤
⎢ ____                                   ⎥
⎢ ╲                                      ⎥
⎢  ╲                                     ⎥
⎢   ╲  ⎡x  y⎤                            ⎥
⎢   ╱  ⎢    ⎥       0             0      ⎥
⎢  ╱   ⎣y  x⎦                            ⎥
⎢ ╱                                      ⎥
⎢ ‾‾‾‾                                   ⎥
⎢n = 1                                   ⎥
⎢                                        ⎥
⎢                m                       ⎥
⎢               ____                     ⎥
⎢               ╲                        ⎥
⎢                ╲                       ⎥
⎢                 ╲  ⎡x  y⎤              ⎥
⎢     0           ╱  ⎢    ⎥       0      ⎥
⎢                ╱   ⎣y  x⎦              ⎥
⎢               ╱                        ⎥
⎢               ‾‾‾‾                     ⎥
⎢              n = 1                     ⎥
⎢                                        ⎥
⎢                              m         ⎥
⎢                             ____       ⎥
⎢                             ╲          ⎥
⎢                              ╲         ⎥
⎢                               ╲  ⎡x  y⎤⎥
⎢     0             0           ╱  ⎢    ⎥⎥
⎢                              ╱   ⎣y  x⎦⎥
⎢                             ╱          ⎥
⎢                             ‾‾‾‾       ⎥
⎣                            n = 1       ⎦

What I think we need is a way to distinguish different classes of object (matrix vs number vs set etc) that is not inherently tied to the the class hierarchy.