Repository: sympy/sympy
Pull request: Fix Matrix Equality and MatrixExpr
<!-- Your title above should be a short description of what
was changed. Do not include the issue number in the title. -->

#### References to other Issues or PRs
Fixes #7842 
Fixes #16042
<!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
format, e.g. "Fixes #1234". See
https://github.com/blog/1506-closing-issues-via-pull-requests . Please also
write a comment on that issue linking back to this pull request once it is
open. -->


#### Brief description of what is fixed or changed
Test case for equality in matrices.

#### Other comments


#### Release Notes

<!-- Write the release notes for this release below. See
https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more information
on how to write release notes. The bot will check your release notes
automatically to see if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
*  matrices
   *  _eval_Eq added for MatExpr
<!-- END RELEASE NOTES -->


Problem statement:
`Equality` and `MatrixExpr`
There are a couple issues here, just not sure about the best way to fix them.

A couple things to setup the examples:

```
>>> A = MatrixSymbol('A', 3, 1)
>>> B = MatrixSymbol('B', 2, 1)
```

1.) A and B cannot be equal, but Symbol doesn't care:

```
>>> Eq(A, B)
A == B
```

I think this could be fixed by giving `MatrixSymbol` its own `_eval_Eq` method, as `ImmutableMatrix` does, but this method could also go in `MatrixExpr` itself. The second is probably preferrable, as it would work for anything containing `MatrixExpr`

2.) Individual elements should be able to be set equivalent (as you can with Symbol), but you can't:

```
>>> Eq(A[1,0], B[1, 0])
# Trimming long error
/home/jimmy/Code/sympy/sympy/core/expr.py in equals(self, other, failing_expression)
    554         # because if the expression ever goes to 0 then the subsequent
    555         # simplification steps that are done will be very fast.
--> 556         diff = factor_terms((self - other).simplify(), radical=True)
    557 
    558         if not diff:

/home/jimmy/Code/sympy/sympy/core/expr.py in simplify(self, ratio, measure)
   2882         from sympy.core.function import count_ops
   2883         measure = measure or count_ops
-> 2884         return simplify(self, ratio, measure)
   2885 
   2886     def nsimplify(self, constants=[], tolerance=None, full=False):

/home/jimmy/Code/sympy/sympy/simplify/simplify.py in simplify(expr, ratio, measure, fu)
   3692     _e = cancel(expr)
   3693     expr1 = shorter(_e, _mexpand(_e).cancel())  # issue 6829
-> 3694     expr2 = shorter(together(expr, deep=True), together(expr1, deep=True))
   3695 
   3696     if ratio is S.Infinity:

/home/jimmy/Code/sympy/sympy/simplify/simplify.py in shorter(*choices)
   3686         if not has_variety(choices):
   3687             return choices[0]
-> 3688         return min(choices, key=measure)
   3689 
   3690     expr = bottom_up(expr, lambda w: w.normal())

/home/jimmy/Code/sympy/sympy/core/function.py in count_ops(expr, visual)
   2263         while args:
   2264             a = args.pop()
-> 2265             if a.is_Rational:
   2266                 #-1/3 = NEG + DIV
   2267                 if a is not S.One:

AttributeError: 'str' object has no attribute 'is_Rational'
```

I have two thoughts on fixing this:
1. Make `MatrixElement` a subclass of Symbol. Perhaps this is a bad idea, but then it would work with the current `Eq` method.
2. Add `MatrixElement` to the `has` statement in `_eval_sides` method of `Relational` in sympy/core/Relational.py. This seems like a hackish thing though.

I think making MatrixElement a subclass of Symbol is the best bet, but doing that change results in tons of errors right now. Thoughts?

Hints from issue discussion:
Ping @mrocklin : Matrix Expressions seem to be your thing. What are your thoughts on this?

Regarding 1.  `Eq` subclasses `Expr` and so will be assumed by much of sympy to be a scalar, not a Matrix.  Historically we've had to reinvent all of the Expr terms, e.g. `Add -> MatAdd`.  One might have to do the same with `Eq -> MatEq`.  `_eval_Eq` seems reasonable as long as it doesn't return something that is an instance of `Expr`.

Regarding 2.  Exprs, and MatrixExprs don't mix well.  They assume different invariants.  In particular MatrixExprs contain non-Basic as args.  MatrixElement which is an Expr but contains MatrixExprs, often has issues as a result.

Here are my thoughts:

*\* warning, thought spew ahead **

## 1. `Eq` and `Matrices`:

I actually do think there should be a `MatEq`, but I also think that `Eq` should work with matrices (both `MatrixExpr` and `Matrix`, and all its variations). I think they should do different things though. Basically there are two kinds of statements that could be thought of to represent matrix equality:
1. Are the expressions equal (A == B, where they are both Matrix like). This is inherently a scalar value (True or False, either they are or they're not, or it remains unevaluated).
2. Are the terms equival? Similar to how other matrix libraries (numpy) work with relationals, will return a matrix of True, False, or unevaluated.

The main question is, what should do what? I have some inner conflicts on this, so I'll just lay them all out. Many disagree, as I'm not actually certain on what's best.

**Lots of thoughts:**
- `Eq` previously always returned a scalar, so `Eq(Matrix, Matrix)` and `Eq(Matrix, MatrixExpr)` should go for number 1.
- `MatrixExpr` operations in the past have always been enclosed in their own types. Anything with a `MatrixExpr` should return a `Mat*` thing.
- `Mat*` in the past have always referred to operations that were done _on_ matrices. Since `MatEq` should be a matrix operation, it should be attempting to return a matrix of `Eq`-like things, and remain unevaluated if it can't.
- Some `MatrixExpr` types actually do evaluate to scalars (or at least 1 x 1 matrices). Do we treat these as scalars? Should this work with `Eq` or `MatEq`?  
- Elementwise will never work between `MatrixExpr` types. They are either algebraically equivalent, or they aren't. This is a scalar value, and should be represented by a `Eq`.

I really am not sure how to handle this.

The main reason for having an elementwise instead of just an overall `Eq` is because if we support `Eq`, other relationals should be added as well. There are other interpretations for these as well, but most will assume it's elementwise (what other matrix libraries do). I do think it would be neat to support linear matrix inequalities though, but this is probably outside scope.

## 2. `Eq` and `MatrixElement`:

So subclassing `MatrixElement` from Symbol doesn't make technical sense. Still, I think more support for this type of operation needs to be looked into.  `MatrixElement`s are inherently scalar, symbolic references to a the value held in a single element of a matrix. I should be able to do anything with them that I can do with a Symbol. They aren't atomic though, as both the base matrix and the indices can be expressions themselves. I'm really not sure how to handle this.  I don't have immediate need for this feature, but it really should be supported.

Latest situation in master:

```pypy
>>> Eq(A, B)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "./sympy/core/relational.py", line 384, in __new__
    dif = lhs - rhs
  File "./sympy/matrices/expressions/matexpr.py", line 24, in __sympifyit_wrapper
    return func(a, b)
  File "./sympy/core/decorators.py", line 132, in binary_op_wrapper
    return func(self, other)
  File "./sympy/matrices/expressions/matexpr.py", line 97, in __sub__
    return MatAdd(self, -other, check=True).doit()
  File "./sympy/matrices/expressions/matadd.py", line 49, in __new__
    validate(*args)
  File "./sympy/matrices/expressions/matadd.py", line 89, in validate
    raise ShapeError("Matrices %s and %s are not aligned"%(A, B))
sympy.matrices.common.ShapeError: Matrices A and -B are not aligned
>>> Eq(A[1,0], B[1, 0])
Eq(A[1, 0], B[1, 0])
```

The first should be easy to fix by defining `_eval_Eq` on MatExpr. It would need to unconditionally return an Equality when it cannot return a bool, to avoid the code in the Eq constructor that does `lhs - rhs`. The `_eval_Eq` would be something like

```py
def _eval_Eq(self, other):
    if not isinstance(other, MatrixExpr:
        return False
    if self.shape != other.shape:
        return False
    if (self - other).is_ZeroMatrix:
        return True
    return Eq(self, other, evaluate=False)
```

`Eq(A[1,0], B[1, 0])` works so we should make sure it is tested, so it doesn't break again.

Commit message: fix-matexpr

Patch to review:
diff --git a/sympy/matrices/expressions/matexpr.py b/sympy/matrices/expressions/matexpr.py
index 124e9994e25e..18457d9a19c5 100644
--- a/sympy/matrices/expressions/matexpr.py
+++ b/sympy/matrices/expressions/matexpr.py
@@ -549,6 +549,15 @@ def applyfunc(self, func):
         from .applyfunc import ElementwiseApplyFunction
         return ElementwiseApplyFunction(func, self)
 
+    def _eval_Eq(self, other):
+        if not isinstance(other, MatrixExpr):
+            return False
+        if self.shape != other.shape:
+            return False
+        if (self - other).is_ZeroMatrix:
+            return True
+        return Eq(self, other, evaluate=False)
+
 
 def _matrix_derivative(expr, x):
     from sympy import Derivative
diff --git a/sympy/matrices/expressions/tests/test_matexpr.py b/sympy/matrices/expressions/tests/test_matexpr.py
index ace95c2802fb..6cfbf8765fca 100644
--- a/sympy/matrices/expressions/tests/test_matexpr.py
+++ b/sympy/matrices/expressions/tests/test_matexpr.py
@@ -360,3 +360,14 @@ def test_issue_2749():
 def test_issue_2750():
     x = MatrixSymbol('x', 1, 1)
     assert (x.T*x).as_explicit()**-1 == Matrix([[x[0, 0]**(-2)]])
+
+
+def test_matrix_equality():
+    from sympy import Eq
+    A = MatrixSymbol('A', 3, 1)
+    B = MatrixSymbol('B', 2, 1)
+    assert Eq(A, B) == False
+    assert Eq(A[1,0], B[1, 0]).func is Eq
+    A = ZeroMatrix(2, 3)
+    B = ZeroMatrix(2, 3)
+    assert Eq(A, B) == True


Code review comments:
diff_hunk: @@ -360,3 +360,14 @@ def test_issue_2749():
 def test_issue_2750():
     x = MatrixSymbol('x', 1, 1)
     assert (x.T*x).as_explicit()**-1 == Matrix([[x[0, 0]**(-2)]])
+
+
+def test_matrix_equality():
+    from sympy import Eq
path: sympy/matrices/expressions/tests/test_matexpr.py
text: @user1:
import Eq at top and please refer to issue number: as I noted on the closed PR, if this ever fails it gives an easy reference to the original issue whereas the subject of the test -- matrix equality -- is self evident from the tests added.

@author:
okay I'll do that.