Repository: sympy/sympy

Issue:
Tensor product of sparse matrices failing
I'd like to form some tensor products of sparse matrices, but this does not seem to work:

    from sympy import diag
    from sympy.physics.quantum import TensorProduct
    from sympy.matrices import SparseMatrix
    spm = SparseMatrix(diag(1,0))
    tp = TensorProduct(spm,spm)

Here's the failure:

    Traceback (most recent call last):
      File "/Users/ddahl/anaconda3/envs/qiskit/lib/python3.8/site-packages/IPython/core/interactiveshell.py", line 3331, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
    File "<ipython-input-9-88417d1a4157>", line 1, in <module>
    tp = TensorProduct(spm,spm)
    File "/Users/ddahl/anaconda3/envs/qiskit/lib/python3.8/site-packages/sympy/physics/quantum/tensorproduct.py", line 123, in __new__
    c_part, new_args = cls.flatten(sympify(args))
    File "/Users/ddahl/anaconda3/envs/qiskit/lib/python3.8/site-packages/sympy/physics/quantum/tensorproduct.py", line 139, in flatten
    cp, ncp = arg.args_cnc()
    File "/Users/ddahl/anaconda3/envs/qiskit/lib/python3.8/site-packages/sympy/matrices/matrices.py", line 2339, in __getattr__
    raise AttributeError(
    AttributeError: ImmutableSparseMatrix has no attribute args_cnc.

When I look in sympy/physics/quantum/tensorproduct.py it looks to me
like there is an attempted test to determine whether the first argument
to the `TensorProduct(...)` function is a sparse SciPy matrix:

    def __new__(cls, *args):
        if isinstance(args[0], (Matrix, numpy_ndarray, scipy_sparse_matrix)):
            return matrix_tensor_product(*args)
        c_part, new_args = cls.flatten(sympify(args))
        c_part = Mul(*c_part)
        if len(new_args) == 0:
            return c_part
        elif len(new_args) == 1:
            return c_part * new_args[0]
        else:
            tp = Expr.__new__(cls, *new_args)
            return c_part * tp

But according to my debugger, the `isinstance(...)` call in the above code block returns false.
Looks to me like that sparse matrix should be of type scipy_sparse_matrix.

My environment is a freshly built Anaconda from two weeks ago and I'm running on macOS.

Discussion:

Patch:
diff --git a/sympy/matrices/dense.py b/sympy/matrices/dense.py
index 7840092da4..e16ccd9379 100644
--- a/sympy/matrices/dense.py
+++ b/sympy/matrices/dense.py
@@ -290,6 +290,8 @@ def _force_mutable(x):
 
 
 class MutableDenseMatrix(DenseMatrix, MatrixBase):
+    __hash__ = None
+
     def __new__(cls, *args, **kwargs):
         return cls._new(*args, **kwargs)
 
diff --git a/sympy/matrices/immutable.py b/sympy/matrices/immutable.py
index 4c4f390156..b559180b3a 100644
--- a/sympy/matrices/immutable.py
+++ b/sympy/matrices/immutable.py
@@ -1,4 +1,3 @@
-from typing import Callable
 from mpmath.matrices.matrices import _matrix
 
 from sympy.core import Basic, Dict, Integer, S, Tuple
@@ -21,7 +20,7 @@ def sympify_mpmath_matrix(arg):
 sympify_converter[_matrix] = sympify_mpmath_matrix
 
 
-class ImmutableDenseMatrix(DenseMatrix, MatrixExpr): # type: ignore
+class ImmutableDenseMatrix(DenseMatrix, MatrixExpr):
     """Create an immutable version of a matrix.
 
     Examples
@@ -49,7 +48,7 @@ class ImmutableDenseMatrix(DenseMatrix, MatrixExpr): # type: ignore
     def __new__(cls, *args, **kwargs):
         return cls._new(*args, **kwargs)
 
-    __hash__ = MatrixExpr.__hash__  # type: Callable[[MatrixExpr], int]
+    __hash__ = MatrixExpr.__hash__
 
     @classmethod
     def _new(cls, *args, **kwargs):
@@ -62,19 +61,15 @@ def _new(cls, *args, **kwargs):
         else:
             rows, cols, flat_list = cls._handle_creation_inputs(*args, **kwargs)
             flat_list = list(flat_list) # create a shallow copy
-        rows = Integer(rows)
-        cols = Integer(cols)
-        if not isinstance(flat_list, Tuple):
-            flat_list = Tuple(*flat_list)
 
-        return Basic.__new__(cls, rows, cols, flat_list)
-
-    @property
-    def _mat(self):
-        # self.args[2] is a Tuple.  Access to the elements
-        # of a tuple are significantly faster than Tuple,
-        # so return the internal tuple.
-        return self.args[2].args
+        obj = Basic.__new__(cls,
+            Integer(rows),
+            Integer(cols),
+            Tuple(*flat_list))
+        obj._rows = rows
+        obj._cols = cols
+        obj._mat = flat_list
+        return obj
 
     def _entry(self, i, j, **kwargs):
         return DenseMatrix.__getitem__(self, (i, j))
@@ -114,15 +109,15 @@ def _eval_extract(self, rowsList, colsList):
 
     @property
     def cols(self):
-        return int(self.args[1])
+        return self._cols
 
     @property
     def rows(self):
-        return int(self.args[0])
+        return self._rows
 
     @property
     def shape(self):
-        return tuple(int(i) for i in self.args[:2])
+        return self._rows, self._cols
 
     def as_immutable(self):
         return self
@@ -138,7 +133,7 @@ def is_diagonalizable(self, reals_only=False, **kwargs):
 ImmutableMatrix = ImmutableDenseMatrix
 
 
-class ImmutableSparseMatrix(SparseMatrix, Basic):
+class ImmutableSparseMatrix(SparseMatrix, MatrixExpr):
     """Create an immutable version of a sparse matrix.
 
     Examples
@@ -163,29 +158,41 @@ class ImmutableSparseMatrix(SparseMatrix, Basic):
     is_Matrix = True
     _class_priority = 9
 
+    def __new__(cls, *args, **kwargs):
+        return cls._new(*args, **kwargs)
+
+    __hash__ = MatrixExpr.__hash__
+
     @classmethod
     def _new(cls, *args, **kwargs):
         s = MutableSparseMatrix(*args)
-        rows = Integer(s.rows)
-        cols = Integer(s.cols)
-        mat = Dict(s._smat)
-        obj = Basic.__new__(cls, rows, cols, mat)
-        obj.rows = s.rows
-        obj.cols = s.cols
-        obj._smat = s._smat
+        rows, cols, smat = s.rows, s.cols, s._smat
+        obj = Basic.__new__(cls, Integer(rows), Integer(cols), Dict(smat))
+        obj._rows = rows
+        obj._cols = cols
+        obj._smat = smat
         return obj
 
-    def __new__(cls, *args, **kwargs):
-        return cls._new(*args, **kwargs)
-
     def __setitem__(self, *args):
         raise TypeError("Cannot set values of ImmutableSparseMatrix")
 
-    def __hash__(self):
-        return hash((type(self).__name__,) + (self.shape, tuple(self._smat)))
+    def _entry(self, i, j, **kwargs):
+        return SparseMatrix.__getitem__(self, (i, j))
 
     _eval_Eq = ImmutableDenseMatrix._eval_Eq
 
+    @property
+    def cols(self):
+        return self._cols
+
+    @property
+    def rows(self):
+        return self._rows
+
+    @property
+    def shape(self):
+        return self._rows, self._cols
+
     def as_immutable(self):
         return self
 
diff --git a/sympy/matrices/matrices.py b/sympy/matrices/matrices.py
index 076fa1f747..e9dfe8612f 100644
--- a/sympy/matrices/matrices.py
+++ b/sympy/matrices/matrices.py
@@ -768,9 +768,6 @@ class MatrixBase(MatrixDeprecated,
     zero = S.Zero
     one = S.One
 
-    # Mutable:
-    __hash__ = None  # type: ignore
-
     # Defined here the same as on Basic.
 
     # We don't define _repr_png_ here because it would add a large amount of
@@ -1268,10 +1265,12 @@ def cross(self, b):
         multiply
         multiply_elementwise
         """
-        if not is_sequence(b):
+        from sympy.matrices.expressions.matexpr import MatrixExpr
+
+        if not isinstance(b, MatrixBase) and not isinstance(b, MatrixExpr):
             raise TypeError(
-                "`b` must be an ordered iterable or Matrix, not %s." %
-                type(b))
+                "{} must be a Matrix, not {}.".format(b, type(b)))
+
         if not (self.rows * self.cols == b.rows * b.cols == 3):
             raise ShapeError("Dimensions incorrect for cross product: %s x %s" %
                              ((self.rows, self.cols), (b.rows, b.cols)))
diff --git a/sympy/printing/latex.py b/sympy/printing/latex.py
index 53e7b8f58c..fc9e61b3b1 100644
--- a/sympy/printing/latex.py
+++ b/sympy/printing/latex.py
@@ -1588,9 +1588,9 @@ def _print_MatrixBase(self, expr):
             out_str = r'\left' + left_delim + out_str + \
                       r'\right' + right_delim
         return out_str % r"\\".join(lines)
-    _print_ImmutableMatrix = _print_ImmutableDenseMatrix \
-                           = _print_Matrix \
-                           = _print_MatrixBase
+
+    _print_ImmutableDenseMatrix = _print_MatrixBase
+    _print_ImmutableSparseMatrix = _print_MatrixBase
 
     def _print_MatrixElement(self, expr):
         return self.parenthesize(expr.parent, PRECEDENCE["Atom"], strict=True)\


Test patch:
diff --git a/sympy/matrices/expressions/tests/test_matexpr.py b/sympy/matrices/expressions/tests/test_matexpr.py
index 1d9bfecc29..6f58ab32c3 100644
--- a/sympy/matrices/expressions/tests/test_matexpr.py
+++ b/sympy/matrices/expressions/tests/test_matexpr.py
@@ -112,11 +112,15 @@ def test_subs():
     C = MatrixSymbol('C', m, l)
 
     assert A.subs(n, m).shape == (m, m)
-
     assert (A*B).subs(B, C) == A*C
-
     assert (A*B).subs(l, n).is_square
 
+    A = SparseMatrix([[1, 2], [3, 4]])
+    B = Matrix([[1, 2], [3, 4]])
+    C, D = MatrixSymbol('C', 2, 2), MatrixSymbol('D', 2, 2)
+
+    assert (C*D).subs({C: A, D: B}) == MatMul(A, B)
+
 
 def test_ZeroMatrix():
     A = MatrixSymbol('A', n, m)
diff --git a/sympy/matrices/tests/test_immutable.py b/sympy/matrices/tests/test_immutable.py
index 7c53e50db6..5d8728f27f 100644
--- a/sympy/matrices/tests/test_immutable.py
+++ b/sympy/matrices/tests/test_immutable.py
@@ -1,29 +1,35 @@
 from itertools import product
 
 from sympy import (ImmutableMatrix, Matrix, eye, zeros, S, Equality,
-        Unequality, ImmutableSparseMatrix, SparseMatrix, sympify,
-        integrate)
+        Unequality, SparseMatrix, sympify, integrate)
+from sympy.matrices.immutable import \
+    ImmutableDenseMatrix, ImmutableSparseMatrix
 from sympy.abc import x, y
 from sympy.testing.pytest import raises
 
-IM = ImmutableMatrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
-ieye = ImmutableMatrix(eye(3))
+IM = ImmutableDenseMatrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
+ISM = ImmutableSparseMatrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
+ieye = ImmutableDenseMatrix(eye(3))
 
 
-def test_immutable_creation():
-    assert IM.shape == (3, 3)
-    assert IM[1, 2] == 6
-    assert IM[2, 2] == 9
+def test_creation():
+    assert IM.shape == ISM.shape == (3, 3)
+    assert IM[1, 2] == ISM[1, 2] == 6
+    assert IM[2, 2] == ISM[2, 2] == 9
 
 
 def test_immutability():
     with raises(TypeError):
         IM[2, 2] = 5
+    with raises(TypeError):
+        ISM[2, 2] = 5
 
 
 def test_slicing():
-    assert IM[1, :] == ImmutableMatrix([[4, 5, 6]])
-    assert IM[:2, :2] == ImmutableMatrix([[1, 2], [4, 5]])
+    assert IM[1, :] == ImmutableDenseMatrix([[4, 5, 6]])
+    assert IM[:2, :2] == ImmutableDenseMatrix([[1, 2], [4, 5]])
+    assert ISM[1, :] == ImmutableSparseMatrix([[4, 5, 6]])
+    assert ISM[:2, :2] == ImmutableSparseMatrix([[1, 2], [4, 5]])
 
 
 def test_subs():
@@ -41,11 +47,13 @@ def test_subs():
 
 
 def test_as_immutable():
-    X = Matrix([[1, 2], [3, 4]])
-    assert sympify(X) == X.as_immutable() == ImmutableMatrix([[1, 2], [3, 4]])
-    X = SparseMatrix(5, 5, {})
-    assert sympify(X) == X.as_immutable() == ImmutableSparseMatrix(
-            [[0 for i in range(5)] for i in range(5)])
+    data = [[1, 2], [3, 4]]
+    X = Matrix(data)
+    assert sympify(X) == X.as_immutable() == ImmutableMatrix(data)
+
+    data = {(0, 0): 1, (0, 1): 2, (1, 0): 3, (1, 1): 4}
+    X = SparseMatrix(2, 2, data)
+    assert sympify(X) == X.as_immutable() == ImmutableSparseMatrix(2, 2, data)
 
 
 def test_function_return_types():
diff --git a/sympy/physics/quantum/tests/test_tensorproduct.py b/sympy/physics/quantum/tests/test_tensorproduct.py
index c8955a08f1..0c47f0c2f5 100644
--- a/sympy/physics/quantum/tests/test_tensorproduct.py
+++ b/sympy/physics/quantum/tests/test_tensorproduct.py
@@ -1,4 +1,6 @@
-from sympy import I, symbols, Matrix
+from sympy import I, symbols
+from sympy.core.expr import unchanged
+from sympy.matrices import Matrix, SparseMatrix
 
 from sympy.physics.quantum.commutator import Commutator as Comm
 from sympy.physics.quantum.tensorproduct import TensorProduct
@@ -17,6 +19,11 @@
 mat2 = Matrix([[2*I, 3], [4*I, 2]])
 
 
+def test_sparse_matrices():
+    spm = SparseMatrix.diag(1, 0)
+    assert unchanged(TensorProduct, spm, spm)
+
+
 def test_tensor_product_dagger():
     assert Dagger(TensorProduct(I*A, B)) == \
         -I*TensorProduct(Dagger(A), Dagger(B))