Repository: sympy/sympy

Issue:
Matrix expressions should be less strict about shape checks
I think the fundamental problem in the matrix expressions is that it validates dimensions using symbolic, structural equality. In other words, behavior like this:

```
>>> var('n m k l', integer=True)
(n, m, k, l)
>>> MatrixSymbol("A", n, m)*MatrixSymbol("B", k, l)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/aaronmeurer/Documents/Python/sympy/sympy/./sympy/matrices/expressions/matexpr.py", line 27, in __sympifyit_wrapper
    return func(a, b)
  File "/Users/aaronmeurer/Documents/Python/sympy/sympy/./sympy/core/decorators.py", line 106, in binary_op_wrapper
    return func(self, other)
  File "/Users/aaronmeurer/Documents/Python/sympy/sympy/./sympy/matrices/expressions/matexpr.py", line 127, in __mul__
    return MatMul(self, other).doit()
  File "/Users/aaronmeurer/Documents/Python/sympy/sympy/./sympy/matrices/expressions/matmul.py", line 60, in __new__
    validate(*matrices)
  File "/Users/aaronmeurer/Documents/Python/sympy/sympy/./sympy/matrices/expressions/matmul.py", line 249, in validate
    raise ShapeError("Matrices %s and %s are not aligned"%(A, B))
sympy.matrices.common.ShapeError: Matrices A and B are not aligned
```

The problem here is that this expression is completely valid if `m` and `k` are equal. There's no reason mathematically to assume they aren't. 

We could instead allow any multiplication of any shapes unless the inner dimensions are explicitly not equal, like `(n x k) * (k + 1 x m)`. If we did this, there would be no need for GenericIdentity or GenericZeroMatrix or GenericOneMatrix. We could just use `n = Dummy('n', integer=True); Identity(n)`.

Of course, we would still want a strict validation function where inner dimensions are required to be exactly equal, but this should be something like an optional flag to `doit`, rather than something that is automatically called in the constructor. 

One problem with this approach is the assumptions system is still not very good at handling generic constraints like `Ne(k, m)`, although in principle the new assumptions machinery is capable of representing them. This isn't just for user constraints: you would need this to validate dimensions across multiple operations, e.g., `(n x n) * (m x m) * (n + 1 x n + 1)` should never validate, because it would imply that `n = m = n + 1`.

_Originally posted by @asmeurer in https://github.com/sympy/sympy/pull/24141#discussion_r998771206_
      

Discussion:
sympy-bot: :white_check_mark:

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot) (v167). I'm here to help you write a release notes entry. Please read the [guide on how to write release notes](https://github.com/sympy/sympy/wiki/Writing-Release-Notes).



Your release notes are in good order.

Here is what the release notes will look like:
* matrices
  - Runtime type, shape checks in `MatAdd`, `MatMul`, `MatPow` and `HadamardProduct` are disabled. ([#24168](https://github.com/sympy/sympy/pull/24168) by [@sylee957](https://github.com/sylee957))

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.12.

<details><summary>Click here to see the pull request description that was parsed.</summary>

    <!-- 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
    <!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
    format, e.g. "Fixes #1234" (see
    https://tinyurl.com/auto-closing for more information). Also, please
    write a comment on that issue linking back to this pull request once it is
    open. -->

    Fixes #24149

    #### Brief description of what is fixed or changed

    I attempt to demonstrate a way to fix #24149 and introduce potentially a new direction to solve these kinds of problems reasonably.
    We are often bothered with asserting mathematically correctness expressions beforehand defining expressions or not.
    However, I think that this has serious tradeoffs for runtime performance.

    However, nowadays I don't think that we may bother with runtime type checking for asserting some objects are mathematically true. Because I believe that when we define logic correct only for mathematically correct objects, it won't have problem for deriving mathematically wrong facts.
    And anything that has invalid definitions before can give arbitrary result like [vacuous truth](https://en.wikipedia.org/wiki/Vacuous_truth).

    Also, I believe that it also aligns with EAFP principle. (https://docs.python.org/3/glossary.html)
    And LBYL can rather be asserted in static time type checking or theorem proving by developers before making code,
    or become the responsibility for users to implement logic and then they will have option to choose reasonable decision procedures for it, so it would not be a responsibility here.

    However, I think that this can be a bit of controversial change because, and also it can make some code for end users more harder to debug.
    I also want to same question for how many libraries had done with it, for example, if any library features partial functions like division, should it bother with checking the dividend is zero and reporting it to users, or should it just propagate the `ZeroDivisionError` for it, or if we can try to catch deeper errors and wrap errors more elegantly for users such that we can guide them with more correct semantic information about how to solve them than giving errors like `object has no attribute 'as_coeff_mul'`, or whether this can be implemented in more higher level frameworks for exceptions.
    There are other articles I introduce to read https://stackoverflow.com/questions/29282819/coq-qarith-division-by-zero-is-zero-why, to have some context about how it is no problem.

    #### Other comments


    #### Release Notes

    <!-- Write the release notes for this release below between the BEGIN and END
    statements. The basic format is a bulleted list with the name of the subpackage
    and the release note for this PR. For example:

    * solvers
      * Added a new solver for logarithmic equations.

    * functions
      * Fixed a bug with log of integers.

    or if no release note(s) should be included use:

    NO ENTRY

    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
      - Runtime type, shape checks in `MatAdd`, `MatMul`, `MatPow` and `HadamardProduct` are disabled.
    <!-- END RELEASE NOTES -->


</details><p>

sympy-bot: ### 🟠

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot) (v167). I've noticed that some of your commits add or delete files. Since this is sometimes done unintentionally, I wanted to alert you about it.

This is an experimental feature of SymPy Bot. If you have any feedback on it, please comment at https://github.com/sympy/sympy-bot/issues/75.

The following commits **add new files**:
* 4020a24a191568e9089eab63234d5a47336becb1:
  - `Untitled.ipynb`

If these files were added/deleted on purpose, you can ignore this message.

sylee957: I have different demo to use `Eq` and `And` evaluation for possibly detecting `m = m + 1` False in `doit`.
I have not tried using the equivalent logic old assumption check yet.
sylee957: I have changed the shape check to `Optional[bool]` type, however, I think that there are some technical difficulties about it because of the compatibility with the existing `GenericZeroMatrix, GenericOneMatrix`.
There was also minor difficulty with the type because `MatrixBase.shape` returns `int` than `Integer`.

I've also implemented a second level of shape checks that uses more stronger `Eq` than the old assumptions check. It is able to infer $M_{n, n+1}$ is not square or such, but there can be tradeoff for speed. I've implemented it only inside `doit` logic, but I'd question if it will be necessary in the future there is another argument that if makes things slower.

Also, I think that we need to undefine `GenericOneMatrix ** -1` or should not guarantee any meaningful result for it, because it complicates the implementation.
github-actions[bot]: Benchmark results from GitHub Actions

Lower numbers are good, higher numbers are bad. A ratio less than 1
means a speed up and greater than 1 means a slowdown. Green lines
beginning with `+` are slowdowns (the PR is slower then master or
master is slower than the previous release). Red lines beginning
with `-` are speedups.

Significantly changed benchmark results (PR vs master)
```diff
       before           after         ratio
     [5e50b55c]       [a66d0874]
-      68.4±0.2μs      24.0±0.05μs     0.35  matrices.TimeMatrixExpression.time_MatAdd
+     9.87±0.01ms      16.2±0.04ms     1.64  matrices.TimeMatrixExpression.time_MatAdd_doit
-       105±0.1μs      27.3±0.05μs     0.26  matrices.TimeMatrixExpression.time_MatMul
+     7.51±0.01ms      16.4±0.03ms     2.18  matrices.TimeMatrixExpression.time_MatMul_doit

```
Significantly changed benchmark results (master vs previous release)
```diff
       before           after         ratio
     [41d90958]       [5e50b55c]
     <sympy-1.11.1^0>                 
-         960±3μs          621±3μs     0.65  solve.TimeSparseSystem.time_linear_eq_to_matrix(10)
-     2.80±0.01ms         1.16±0ms     0.42  solve.TimeSparseSystem.time_linear_eq_to_matrix(20)
-     5.62±0.02ms      1.72±0.01ms     0.31  solve.TimeSparseSystem.time_linear_eq_to_matrix(30)

```
Full benchmark results can be found as artifacts in GitHub Actions
(click on checks at the top of the PR).
sylee957: The benchmark results indicates that there is slowdown for `doit` when we use more stronger logic to test out the shape.
oscarbenjamin: I was expecting at least one of these to raise:
```python
In [7]: n, m = symbols('n, m')

In [8]: M = MatrixSymbol('M', m, m)

In [9]: N = MatrixSymbol('N', n, n)

In [10]: M1 = MatrixSymbol('M1', m+1, m)

In [11]: MatAdd(M, N, M1)
Out[11]: M + N + M₁

In [12]: MatAdd(M, M1, N)
Out[12]: M + M₁ + N
```
I'm not sure what would be the best way of handling these comparisons without all to all checking.
oscarbenjamin: > I was expecting at least one of these to raise:

To be clear I don't think it is completely necessary that it should raise. I'm just surprised that it does not because it least this equality evaluates:
```python
In [17]: Eq(m, m+1)
Out[17]: False
```
In general though I am not sure what is an efficient way of handling a situation where many symbolic expressions are expected to be equal without doing the `N^2` comparisons.
sylee957: > I am not sure what is an efficient way of handling a situation where many symbolic expressions are expected to be equal without doing the N^2 comparisons.

The problem is such worst time complexity is achieved every time when the input supplied by user is correct

Patch:
diff --git a/sympy/matrices/expressions/_shape.py b/sympy/matrices/expressions/_shape.py
new file mode 100644
index 0000000000..0e43410648
--- /dev/null
+++ b/sympy/matrices/expressions/_shape.py
@@ -0,0 +1,102 @@
+from sympy.core.relational import Eq
+from sympy.core.expr import Expr
+from sympy.core.numbers import Integer
+from sympy.logic.boolalg import Boolean, And
+from sympy.matrices.expressions.matexpr import MatrixExpr
+from sympy.matrices.common import ShapeError
+from typing import Union
+
+
+def is_matadd_valid(*args: MatrixExpr) -> Boolean:
+    """Return the symbolic condition how ``MatAdd``, ``HadamardProduct``
+    makes sense.
+
+    Parameters
+    ==========
+
+    args
+        The list of arguments of matrices to be tested for.
+
+    Examples
+    ========
+
+    >>> from sympy import MatrixSymbol, symbols
+    >>> from sympy.matrices.expressions._shape import is_matadd_valid
+
+    >>> m, n, p, q = symbols('m n p q')
+    >>> A = MatrixSymbol('A', m, n)
+    >>> B = MatrixSymbol('B', p, q)
+    >>> is_matadd_valid(A, B)
+    Eq(m, p) & Eq(n, q)
+    """
+    rows, cols = zip(*(arg.shape for arg in args))
+    return And(
+        *(Eq(i, j) for i, j in zip(rows[:-1], rows[1:])),
+        *(Eq(i, j) for i, j in zip(cols[:-1], cols[1:])),
+    )
+
+
+def is_matmul_valid(*args: Union[MatrixExpr, Expr]) -> Boolean:
+    """Return the symbolic condition how ``MatMul`` makes sense
+
+    Parameters
+    ==========
+
+    args
+        The list of arguments of matrices and scalar expressions to be tested
+        for.
+
+    Examples
+    ========
+
+    >>> from sympy import MatrixSymbol, symbols
+    >>> from sympy.matrices.expressions._shape import is_matmul_valid
+
+    >>> m, n, p, q = symbols('m n p q')
+    >>> A = MatrixSymbol('A', m, n)
+    >>> B = MatrixSymbol('B', p, q)
+    >>> is_matmul_valid(A, B)
+    Eq(n, p)
+    """
+    rows, cols = zip(*(arg.shape for arg in args if isinstance(arg, MatrixExpr)))
+    return And(*(Eq(i, j) for i, j in zip(cols[:-1], rows[1:])))
+
+
+def is_square(arg: MatrixExpr, /) -> Boolean:
+    """Return the symbolic condition how the matrix is assumed to be square
+
+    Parameters
+    ==========
+
+    arg
+        The matrix to be tested for.
+
+    Examples
+    ========
+
+    >>> from sympy import MatrixSymbol, symbols
+    >>> from sympy.matrices.expressions._shape import is_square
+
+    >>> m, n = symbols('m n')
+    >>> A = MatrixSymbol('A', m, n)
+    >>> is_square(A)
+    Eq(m, n)
+    """
+    return Eq(arg.rows, arg.cols)
+
+
+def validate_matadd_integer(*args: MatrixExpr) -> None:
+    """Validate matrix shape for addition only for integer values"""
+    rows, cols = zip(*map(lambda x: x.shape, args))
+    if len(set(filter(lambda x: isinstance(x, (int, Integer)), rows))) > 1:
+        raise ShapeError(f"Matrices have mismatching shape: {rows}")
+    if len(set(filter(lambda x: isinstance(x, (int, Integer)), cols))) > 1:
+        raise ShapeError(f"Matrices have mismatching shape: {cols}")
+
+
+def validate_matmul_integer(*args: MatrixExpr) -> None:
+    """Validate matrix shape for multiplication only for integer values"""
+    for A, B in zip(args[:-1], args[1:]):
+        i, j = A.cols, B.rows
+        if isinstance(i, (int, Integer)) and isinstance(j, (int, Integer)) and i != j:
+            raise ShapeError("Matrices are not aligned", i, j)
diff --git a/sympy/matrices/expressions/determinant.py b/sympy/matrices/expressions/determinant.py
index af9cd2feb4..db40424b17 100644
--- a/sympy/matrices/expressions/determinant.py
+++ b/sympy/matrices/expressions/determinant.py
@@ -27,7 +27,7 @@ def __new__(cls, mat):
         if not mat.is_Matrix:
             raise TypeError("Input to Determinant, %s, not a matrix" % str(mat))
 
-        if not mat.is_square:
+        if mat.is_square is False:
             raise NonSquareMatrixError("Det of a non-square matrix")
 
         return Basic.__new__(cls, mat)
diff --git a/sympy/matrices/expressions/hadamard.py b/sympy/matrices/expressions/hadamard.py
index 4f9721d9c1..08e610e1db 100644
--- a/sympy/matrices/expressions/hadamard.py
+++ b/sympy/matrices/expressions/hadamard.py
@@ -5,8 +5,8 @@
 from sympy.core.expr import ExprBuilder
 from sympy.core.sorting import default_sort_key
 from sympy.functions.elementary.exponential import log
-from sympy.matrices.common import ShapeError
 from sympy.matrices.expressions.matexpr import MatrixExpr
+from sympy.matrices.expressions._shape import validate_matadd_integer as validate
 from sympy.matrices.expressions.special import ZeroMatrix, OneMatrix
 from sympy.strategies import (
     unpack, flatten, condition, exhaust, rm_id, sort
@@ -33,11 +33,9 @@ def hadamard_product(*matrices):
     """
     if not matrices:
         raise TypeError("Empty Hadamard product is undefined")
-    validate(*matrices)
     if len(matrices) == 1:
         return matrices[0]
-    else:
-        return HadamardProduct(*matrices).doit()
+    return HadamardProduct(*matrices).doit()
 
 
 class HadamardProduct(MatrixExpr):
@@ -69,19 +67,18 @@ def __new__(cls, *args, evaluate=False, check=None):
         if len(args) == 0:
             # We currently don't have a way to support one-matrices of generic dimensions:
             raise ValueError("HadamardProduct needs at least one argument")
+
+        if not all(isinstance(arg, MatrixExpr) for arg in args):
+            raise TypeError("Mix of Matrix and Scalar symbols")
+
         if check is not None:
             sympy_deprecation_warning(
                 "Passing check to HadamardProduct is deprecated and the check argument will be removed in a future version.",
                 deprecated_since_version="1.11",
                 active_deprecations_target='remove-check-argument-from-matrix-operations')
 
-        if check in (True, None):
+        if check is not False:
             validate(*args)
-        else:
-            sympy_deprecation_warning(
-                "Passing check=False to HadamardProduct is deprecated and the check argument will be removed in a future version.",
-                deprecated_since_version="1.11",
-                active_deprecations_target='remove-check-argument-from-matrix-operations')
 
         obj = super().__new__(cls, *args)
         if evaluate:
@@ -100,10 +97,11 @@ def _eval_transpose(self):
         return HadamardProduct(*list(map(transpose, self.args)))
 
     def doit(self, **hints):
-        expr = self.func(*[i.doit(**hints) for i in self.args])
+        expr = self.func(*(i.doit(**hints) for i in self.args))
         # Check for explicit matrices:
         from sympy.matrices.matrices import MatrixBase
         from sympy.matrices.immutable import ImmutableMatrix
+
         explicit = [i for i in expr.args if isinstance(i, MatrixBase)]
         if explicit:
             remainder = [i for i in expr.args if i not in explicit]
@@ -164,15 +162,6 @@ def _eval_derivative_matrix_lines(self, x):
         return lines
 
 
-def validate(*args):
-    if not all(arg.is_Matrix for arg in args):
-        raise TypeError("Mix of Matrix and Scalar symbols")
-    A = args[0]
-    for B in args[1:]:
-        if A.shape != B.shape:
-            raise ShapeError("Matrices %s and %s are not aligned" % (A, B))
-
-
 # TODO Implement algorithm for rewriting Hadamard product as diagonal matrix
 # if matmul identy matrix is multiplied.
 def canonicalize(x):
@@ -381,15 +370,12 @@ class HadamardPower(MatrixExpr):
     def __new__(cls, base, exp):
         base = sympify(base)
         exp = sympify(exp)
+
         if base.is_scalar and exp.is_scalar:
             return base ** exp
 
-        if base.is_Matrix and exp.is_Matrix and base.shape != exp.shape:
-            raise ValueError(
-                'The shape of the base {} and '
-                'the shape of the exponent {} do not match.'
-                .format(base.shape, exp.shape)
-                )
+        if isinstance(base, MatrixExpr) and isinstance(exp, MatrixExpr):
+            validate(base, exp)
 
         obj = super().__new__(cls, base, exp)
         return obj
diff --git a/sympy/matrices/expressions/inverse.py b/sympy/matrices/expressions/inverse.py
index a0cbc1b7da..7cb8436915 100644
--- a/sympy/matrices/expressions/inverse.py
+++ b/sympy/matrices/expressions/inverse.py
@@ -39,7 +39,7 @@ def __new__(cls, mat, exp=S.NegativeOne):
         exp = _sympify(exp)
         if not mat.is_Matrix:
             raise TypeError("mat should be a matrix")
-        if not mat.is_square:
+        if mat.is_square is False:
             raise NonSquareMatrixError("Inverse of non-square matrix %s" % mat)
         return Basic.__new__(cls, mat, exp)
 
diff --git a/sympy/matrices/expressions/kronecker.py b/sympy/matrices/expressions/kronecker.py
index 9c8b836115..e12909422c 100644
--- a/sympy/matrices/expressions/kronecker.py
+++ b/sympy/matrices/expressions/kronecker.py
@@ -77,7 +77,6 @@ def kronecker_product(*matrices):
     """
     if not matrices:
         raise TypeError("Empty Kronecker product is undefined")
-    validate(*matrices)
     if len(matrices) == 1:
         return matrices[0]
     else:
diff --git a/sympy/matrices/expressions/matadd.py b/sympy/matrices/expressions/matadd.py
index 1504c1ed74..6402743a08 100644
--- a/sympy/matrices/expressions/matadd.py
+++ b/sympy/matrices/expressions/matadd.py
@@ -5,13 +5,13 @@
 from sympy.core.add import add, Add, _could_extract_minus_sign
 from sympy.core.sorting import default_sort_key
 from sympy.functions import adjoint
-from sympy.matrices.common import ShapeError
 from sympy.matrices.matrices import MatrixBase
 from sympy.matrices.expressions.transpose import transpose
 from sympy.strategies import (rm_id, unpack, flatten, sort, condition,
     exhaust, do_one, glom)
 from sympy.matrices.expressions.matexpr import MatrixExpr
 from sympy.matrices.expressions.special import ZeroMatrix, GenericZeroMatrix
+from sympy.matrices.expressions._shape import validate_matadd_integer as validate
 from sympy.utilities.iterables import sift
 from sympy.utilities.exceptions import sympy_deprecation_warning
 
@@ -45,6 +45,9 @@ def __new__(cls, *args, evaluate=False, check=None, _sympify=True):
         if _sympify:
             args = list(map(sympify, args))
 
+        if not all(isinstance(arg, MatrixExpr) for arg in args):
+            raise TypeError("Mix of Matrix and Scalar symbols")
+
         obj = Basic.__new__(cls, *args)
 
         if check is not None:
@@ -52,15 +55,9 @@ def __new__(cls, *args, evaluate=False, check=None, _sympify=True):
                 "Passing check to MatAdd is deprecated and the check argument will be removed in a future version.",
                 deprecated_since_version="1.11",
                 active_deprecations_target='remove-check-argument-from-matrix-operations')
-        if check in (True, None):
-            if not any(isinstance(i, MatrixExpr) for i in args):
-                return Add.fromiter(args)
+
+        if check is not False:
             validate(*args)
-        else:
-            sympy_deprecation_warning(
-                "Passing check=False to MatAdd is deprecated and the check argument will be removed in a future version.",
-                deprecated_since_version="1.11",
-                active_deprecations_target='remove-check-argument-from-matrix-operations')
 
         if evaluate:
             obj = cls._evaluate(obj)
@@ -109,14 +106,6 @@ def _eval_derivative_matrix_lines(self, x):
 
 add.register_handlerclass((Add, MatAdd), MatAdd)
 
-def validate(*args):
-    if not all(arg.is_Matrix for arg in args):
-        raise TypeError("Mix of Matrix and Scalar symbols")
-
-    A = args[0]
-    for B in args[1:]:
-        if A.shape != B.shape:
-            raise ShapeError("Matrices %s and %s are not aligned"%(A, B))
 
 factor_of = lambda arg: arg.as_coeff_mmul()[0]
 matrix_of = lambda arg: unpack(arg.as_coeff_mmul()[1])
diff --git a/sympy/matrices/expressions/matexpr.py b/sympy/matrices/expressions/matexpr.py
index fb45210bc7..5358f2e97a 100644
--- a/sympy/matrices/expressions/matexpr.py
+++ b/sympy/matrices/expressions/matexpr.py
@@ -171,8 +171,13 @@ def cols(self):
         return self.shape[1]
 
     @property
-    def is_square(self):
-        return self.rows == self.cols
+    def is_square(self) -> bool | None:
+        rows, cols = self.shape
+        if isinstance(rows, Integer) and isinstance(cols, Integer):
+            return rows == cols
+        if rows == cols:
+            return True
+        return None
 
     def _eval_conjugate(self):
         from sympy.matrices.expressions.adjoint import Adjoint
@@ -258,7 +263,7 @@ def T(self):
         return self.transpose()
 
     def inverse(self):
-        if not self.is_square:
+        if self.is_square is False:
             raise NonSquareMatrixError('Inverse of non-square matrix')
         return self._eval_inverse()
 
diff --git a/sympy/matrices/expressions/matmul.py b/sympy/matrices/expressions/matmul.py
index 314474a6da..baaf8a21ce 100644
--- a/sympy/matrices/expressions/matmul.py
+++ b/sympy/matrices/expressions/matmul.py
@@ -7,9 +7,10 @@
 from sympy.functions import adjoint
 from sympy.strategies import (rm_id, unpack, typed, flatten, exhaust,
         do_one, new)
-from sympy.matrices.common import ShapeError, NonInvertibleMatrixError
+from sympy.matrices.common import NonInvertibleMatrixError
 from sympy.matrices.matrices import MatrixBase
 from sympy.utilities.exceptions import sympy_deprecation_warning
+from sympy.matrices.expressions._shape import validate_matmul_integer as validate
 
 from .inverse import Inverse
 from .matexpr import MatrixExpr
@@ -56,13 +57,8 @@ def __new__(cls, *args, evaluate=False, check=None, _sympify=True):
                 deprecated_since_version="1.11",
                 active_deprecations_target='remove-check-argument-from-matrix-operations')
 
-        if check in (True, None):
+        if check is not False:
             validate(*matrices)
-        else:
-            sympy_deprecation_warning(
-                "Passing check=False to MatMul is deprecated and the check argument will be removed in a future version.",
-                deprecated_since_version="1.11",
-                active_deprecations_target='remove-check-argument-from-matrix-operations')
 
         if not matrices:
             # Should it be
@@ -184,17 +180,18 @@ def _eval_determinant(self):
         return factor**self.rows * Mul(*list(map(Determinant, square_matrices)))
 
     def _eval_inverse(self):
-        try:
-            return MatMul(*[
+        if all(arg.is_square for arg in self.args if isinstance(arg, MatrixExpr)):
+            return MatMul(*(
                 arg.inverse() if isinstance(arg, MatrixExpr) else arg**-1
-                    for arg in self.args[::-1]]).doit()
-        except ShapeError:
-            return Inverse(self)
+                    for arg in self.args[::-1]
+                )
+            ).doit()
+        return Inverse(self)
 
     def doit(self, **hints):
         deep = hints.get('deep', True)
         if deep:
-            args = [arg.doit(**hints) for arg in self.args]
+            args = tuple(arg.doit(**hints) for arg in self.args)
         else:
             args = self.args
 
@@ -241,16 +238,8 @@ def _eval_derivative_matrix_lines(self, x):
 
 mul.register_handlerclass((Mul, MatMul), MatMul)
 
-def validate(*matrices):
-    """ Checks for valid shapes for args of MatMul """
-    for i in range(len(matrices)-1):
-        A, B = matrices[i:i+2]
-        if A.cols != B.rows:
-            raise ShapeError("Matrices %s and %s are not aligned"%(A, B))
 
 # Rules
-
-
 def newmul(*args):
     if args[0] == 1:
         args = args[1:]
diff --git a/sympy/matrices/expressions/matpow.py b/sympy/matrices/expressions/matpow.py
index 3f075764a6..ce380b0f30 100644
--- a/sympy/matrices/expressions/matpow.py
+++ b/sympy/matrices/expressions/matpow.py
@@ -1,4 +1,3 @@
-from sympy.matrices.common import NonSquareMatrixError
 from .matexpr import MatrixExpr
 from .special import Identity
 from sympy.core import S
@@ -7,15 +6,16 @@
 from sympy.core.power import Pow
 from sympy.core.sympify import _sympify
 from sympy.matrices import MatrixBase
+from sympy.matrices.common import NonSquareMatrixError
 
 
 class MatPow(MatrixExpr):
-
     def __new__(cls, base, exp, evaluate=False, **options):
         base = _sympify(base)
         if not base.is_Matrix:
             raise TypeError("MatPow base should be a matrix")
-        if not base.is_square:
+
+        if base.is_square is False:
             raise NonSquareMatrixError("Power of non-square matrix %s" % base)
 
         exp = _sympify(exp)
@@ -59,7 +59,7 @@ def _entry(self, i, j, **kwargs):
 
     def doit(self, **hints):
         if hints.get('deep', True):
-            base, exp = [arg.doit(**hints) for arg in self.args]
+            base, exp = (arg.doit(**hints) for arg in self.args)
         else:
             base, exp = self.args
 
diff --git a/sympy/matrices/expressions/special.py b/sympy/matrices/expressions/special.py
index 776d019031..fb7eeb7904 100644
--- a/sympy/matrices/expressions/special.py
+++ b/sympy/matrices/expressions/special.py
@@ -195,6 +195,10 @@ def cols(self):
     def shape(self):
         raise TypeError("GenericIdentity does not have a specified shape")
 
+    @property
+    def is_square(self):
+        return True
+
     # Avoid Matrix.__eq__ which might call .shape
     def __eq__(self, other):
         return isinstance(other, GenericIdentity)
diff --git a/sympy/matrices/expressions/trace.py b/sympy/matrices/expressions/trace.py
index 048c2ec26a..02ec6650f6 100644
--- a/sympy/matrices/expressions/trace.py
+++ b/sympy/matrices/expressions/trace.py
@@ -37,7 +37,7 @@ def __new__(cls, mat):
         if not mat.is_Matrix:
             raise TypeError("input to Trace, %s, is not a matrix" % str(mat))
 
-        if not mat.is_square:
+        if mat.is_square is False:
             raise NonSquareMatrixError("Trace of a non-square matrix")
 
         return Basic.__new__(cls, mat)
diff --git a/sympy/physics/units/__init__.py b/sympy/physics/units/__init__.py
index 361dbfdbcb..18fee8c601 100644
--- a/sympy/physics/units/__init__.py
+++ b/sympy/physics/units/__init__.py
@@ -164,7 +164,6 @@
     josephson_constant,
     von_klitzing_constant,
     Da, dalton, amu, amus, atomic_mass_unit, atomic_mass_constant,
-    me, electron_rest_mass,
     gee, gees, acceleration_due_to_gravity,
     u0, magnetic_constant, vacuum_permeability,
     e0, electric_constant, vacuum_permittivity,
@@ -399,7 +398,6 @@ def find_unit(quantity, unit_system="SI"):
     'josephson_constant',
     'von_klitzing_constant',
     'Da', 'dalton', 'amu', 'amus', 'atomic_mass_unit', 'atomic_mass_constant',
-    'me', 'electron_rest_mass',
     'gee', 'gees', 'acceleration_due_to_gravity',
     'u0', 'magnetic_constant', 'vacuum_permeability',
     'e0', 'electric_constant', 'vacuum_permittivity',
diff --git a/sympy/physics/units/definitions/__init__.py b/sympy/physics/units/definitions/__init__.py
index ab5a0e1c2f..f60eb838be 100644
--- a/sympy/physics/units/definitions/__init__.py
+++ b/sympy/physics/units/definitions/__init__.py
@@ -82,7 +82,6 @@
     josephson_constant,
     von_klitzing_constant,
     Da, dalton, amu, amus, atomic_mass_unit, atomic_mass_constant,
-    me, electron_rest_mass,
     gee, gees, acceleration_due_to_gravity,
     u0, magnetic_constant, vacuum_permeability,
     e0, electric_constant, vacuum_permittivity,
@@ -214,7 +213,6 @@
     'josephson_constant',
     'von_klitzing_constant',
     'Da', 'dalton', 'amu', 'amus', 'atomic_mass_unit', 'atomic_mass_constant',
-    'me', 'electron_rest_mass',
     'gee', 'gees', 'acceleration_due_to_gravity',
     'u0', 'magnetic_constant', 'vacuum_permeability',
     'e0', 'electric_constant', 'vacuum_permittivity',
diff --git a/sympy/physics/units/definitions/unit_definitions.py b/sympy/physics/units/definitions/unit_definitions.py
index cf51b51f45..4939fa65ac 100644
--- a/sympy/physics/units/definitions/unit_definitions.py
+++ b/sympy/physics/units/definitions/unit_definitions.py
@@ -122,9 +122,6 @@
 t = metric_ton = tonne = Quantity("tonne", abbrev="t")
 tonne.set_global_relative_scale_factor(mega, gram)
 
-# Electron rest mass
-me = electron_rest_mass = Quantity("electron_rest_mass", abbrev="me")
-
 
 # Common length units
 
diff --git a/sympy/physics/units/systems/si.py b/sympy/physics/units/systems/si.py
index 2bfa780587..1517273ce1 100644
--- a/sympy/physics/units/systems/si.py
+++ b/sympy/physics/units/systems/si.py
@@ -25,7 +25,7 @@
     coulomb, volt, ohm, siemens, farad, henry, tesla, weber, dioptre, lux,
     katal, gray, becquerel, inch, liter, julian_year, gravitational_constant,
     speed_of_light, elementary_charge, planck, hbar, electronvolt,
-    avogadro_number, avogadro_constant, boltzmann_constant, electron_rest_mass,
+    avogadro_number, avogadro_constant, boltzmann_constant,
     stefan_boltzmann_constant, Da, atomic_mass_constant, molar_gas_constant,
     faraday_constant, josephson_constant, von_klitzing_constant,
     acceleration_due_to_gravity, magnetic_constant, vacuum_permittivity,
@@ -229,10 +229,6 @@
 SI.set_quantity_dimension(vacuum_impedance, impedance)
 SI.set_quantity_scale_factor(vacuum_impedance, u0 * c)
 
-# Electron rest mass
-SI.set_quantity_dimension(electron_rest_mass, mass)
-SI.set_quantity_scale_factor(electron_rest_mass, 9.1093837015e-31*kilogram)
-
 # Coulomb's constant:
 SI.set_quantity_dimension(coulomb_constant, force * length ** 2 / charge ** 2)
 SI.set_quantity_scale_factor(coulomb_constant, 1/(4*pi*vacuum_permittivity))
@@ -368,7 +364,7 @@
     'planck_angular_frequency', 'ohm', 'pound', 'planck_pressure', 'G', 'psi',
     'dHg0', 'von_klitzing_constant', 'planck_length', 'avogadro_number',
     'mole', 'acceleration', 'information', 'planck_energy_density',
-    'mebibyte', 's', 'acceleration_due_to_gravity', 'electron_rest_mass',
+    'mebibyte', 's', 'acceleration_due_to_gravity',
     'planck_temperature', 'units', 'mass', 'dimsys_MKSA', 'kelvin', 'kPa',
     'boltzmann', 'milli_mass_unit', 'planck_impedance', 'electric_constant',
     'derived_dims', 'kg', 'coulomb', 'siemens', 'byte', 'magnetic_flux',


Test patch:
diff --git a/sympy/matrices/expressions/tests/test_applyfunc.py b/sympy/matrices/expressions/tests/test_applyfunc.py
index 8500f7d191..d98732e275 100644
--- a/sympy/matrices/expressions/tests/test_applyfunc.py
+++ b/sympy/matrices/expressions/tests/test_applyfunc.py
@@ -7,8 +7,6 @@
 from sympy.matrices.expressions.matexpr import MatrixSymbol
 from sympy.matrices.expressions.matmul import MatMul
 from sympy.simplify.simplify import simplify
-from sympy.testing.pytest import raises
-from sympy.matrices.common import ShapeError
 
 
 X = MatrixSymbol("X", 3, 3)
@@ -56,8 +54,6 @@ def test_applyfunc_matrix():
     expr = ElementwiseApplyFunction(exp, Z)*ElementwiseApplyFunction(exp, Z.T)
     assert expr.shape == (2, 2)
 
-    raises(ShapeError, lambda: ElementwiseApplyFunction(exp, Z)*ElementwiseApplyFunction(exp, Z))
-
     M = Matrix([[x, y], [z, t]])
     expr = ElementwiseApplyFunction(sin, M)
     assert isinstance(expr, ElementwiseApplyFunction)
@@ -79,8 +75,6 @@ def test_applyfunc_matrix():
     expr3 = expr*M
     assert expr3.shape == (k, t)
 
-    raises(ShapeError, lambda: M*expr)
-
     expr1 = ElementwiseApplyFunction(lambda x: x+1, Xk)
     expr2 = ElementwiseApplyFunction(lambda x: x, Xk)
     assert expr1 != expr2
diff --git a/sympy/matrices/expressions/tests/test_hadamard.py b/sympy/matrices/expressions/tests/test_hadamard.py
index 1c664d0a75..c61d2ac899 100644
--- a/sympy/matrices/expressions/tests/test_hadamard.py
+++ b/sympy/matrices/expressions/tests/test_hadamard.py
@@ -1,10 +1,11 @@
 from sympy.matrices.dense import Matrix, eye
+from sympy.matrices.common import ShapeError
 from sympy.matrices.expressions.matadd import MatAdd
-from sympy.matrices.expressions.special import (Identity, OneMatrix, ZeroMatrix)
+from sympy.matrices.expressions.special import Identity, OneMatrix, ZeroMatrix
 from sympy.core import symbols
 from sympy.testing.pytest import raises, warns_deprecated_sympy
 
-from sympy.matrices import ShapeError, MatrixSymbol
+from sympy.matrices import MatrixSymbol
 from sympy.matrices.expressions import (HadamardProduct, hadamard_product, HadamardPower, hadamard_power)
 
 n, m, k = symbols('n,m,k')
@@ -17,7 +18,6 @@
 def test_HadamardProduct():
     assert HadamardProduct(A, B, A).shape == A.shape
 
-    raises(ShapeError, lambda: HadamardProduct(A, B.T))
     raises(TypeError,  lambda: HadamardProduct(A, n))
     raises(TypeError,  lambda: HadamardProduct(A, 1))
 
@@ -63,17 +63,13 @@ def test_hadamard():
     m, n, p = symbols('m, n, p', integer=True)
     A = MatrixSymbol('A', m, n)
     B = MatrixSymbol('B', m, n)
-    C = MatrixSymbol('C', m, p)
     X = MatrixSymbol('X', m, m)
     I = Identity(m)
-    with raises(TypeError):
-        hadamard_product()
+
+    raises(TypeError, lambda: hadamard_product())
     assert hadamard_product(A) == A
     assert isinstance(hadamard_product(A, B), HadamardProduct)
     assert hadamard_product(A, B).doit() == hadamard_product(A, B)
-    with raises(ShapeError):
-        hadamard_product(A, C)
-        hadamard_product(A, I)
     assert hadamard_product(X, I) == HadamardProduct(I, X)
     assert isinstance(hadamard_product(X, I), HadamardProduct)
 
@@ -133,3 +129,13 @@ def test_hadamard_power_explicit():
         Matrix([
             [A[0, 0]**B[0, 0], A[0, 1]**B[0, 1]],
             [A[1, 0]**B[1, 0], A[1, 1]**B[1, 1]]])
+
+
+def test_shape_error():
+    A = MatrixSymbol('A', 2, 3)
+    B = MatrixSymbol('B', 3, 3)
+    raises(ShapeError, lambda: HadamardProduct(A, B))
+    raises(ShapeError, lambda: HadamardPower(A, B))
+    A = MatrixSymbol('A', 3, 2)
+    raises(ShapeError, lambda: HadamardProduct(A, B))
+    raises(ShapeError, lambda: HadamardPower(A, B))
diff --git a/sympy/matrices/expressions/tests/test_inverse.py b/sympy/matrices/expressions/tests/test_inverse.py
index 342a47e910..baf1289544 100644
--- a/sympy/matrices/expressions/tests/test_inverse.py
+++ b/sympy/matrices/expressions/tests/test_inverse.py
@@ -1,6 +1,6 @@
 from sympy.core import symbols, S
 from sympy.matrices.expressions import MatrixSymbol, Inverse, MatPow, ZeroMatrix, OneMatrix
-from sympy.matrices.common import NonSquareMatrixError, NonInvertibleMatrixError
+from sympy.matrices.common import NonInvertibleMatrixError, NonSquareMatrixError
 from sympy.matrices import eye, Identity
 from sympy.testing.pytest import raises
 from sympy.assumptions.ask import Q
@@ -45,11 +45,7 @@ def test_inverse():
     assert isinstance(OneMatrix(n, n).I, Inverse)
 
 def test_inverse_non_invertible():
-    raises(NonSquareMatrixError, lambda: Inverse(A))
-    raises(NonSquareMatrixError, lambda: Inverse(A*B))
-    raises(NonSquareMatrixError, lambda: ZeroMatrix(n, m).I)
     raises(NonInvertibleMatrixError, lambda: ZeroMatrix(n, n).I)
-    raises(NonSquareMatrixError, lambda: OneMatrix(n, m).I)
     raises(NonInvertibleMatrixError, lambda: OneMatrix(2, 2).I)
 
 def test_refine():
@@ -59,3 +55,8 @@ def test_refine():
 def test_inverse_matpow_canonicalization():
     A = MatrixSymbol('A', 3, 3)
     assert Inverse(MatPow(A, 3)).doit() == MatPow(Inverse(A), 3).doit()
+
+
+def test_nonsquare_error():
+    A = MatrixSymbol('A', 3, 4)
+    raises(NonSquareMatrixError, lambda: Inverse(A))
diff --git a/sympy/matrices/expressions/tests/test_matadd.py b/sympy/matrices/expressions/tests/test_matadd.py
index d13890557c..2c544a8e83 100644
--- a/sympy/matrices/expressions/tests/test_matadd.py
+++ b/sympy/matrices/expressions/tests/test_matadd.py
@@ -1,5 +1,6 @@
 from sympy.matrices.expressions import MatrixSymbol, MatAdd, MatPow, MatMul
 from sympy.matrices.expressions.special import GenericZeroMatrix, ZeroMatrix
+from sympy.matrices.common import ShapeError
 from sympy.matrices import eye, ImmutableMatrix
 from sympy.core import Add, Basic, S
 from sympy.core.add import add
@@ -46,3 +47,12 @@ def test_zero_matrix_add():
 @XFAIL
 def test_matrix_Add_with_scalar():
     raises(TypeError, lambda: Add(0, ZeroMatrix(2, 2)))
+
+
+def test_shape_error():
+    A = MatrixSymbol('A', 2, 3)
+    B = MatrixSymbol('B', 3, 3)
+    raises(ShapeError, lambda: MatAdd(A, B))
+
+    A = MatrixSymbol('A', 3, 2)
+    raises(ShapeError, lambda: MatAdd(A, B))
diff --git a/sympy/matrices/expressions/tests/test_matexpr.py b/sympy/matrices/expressions/tests/test_matexpr.py
index 87a76441d5..36b3846c2e 100644
--- a/sympy/matrices/expressions/tests/test_matexpr.py
+++ b/sympy/matrices/expressions/tests/test_matexpr.py
@@ -12,8 +12,9 @@
 from sympy.functions import sin, cos, tan, sqrt, cbrt, exp
 from sympy.simplify import simplify
 from sympy.matrices import (ImmutableMatrix, Inverse, MatAdd, MatMul,
-        MatPow, Matrix, MatrixExpr, MatrixSymbol, ShapeError,
-        SparseMatrix, Transpose, Adjoint, NonSquareMatrixError, MatrixSet)
+        MatPow, Matrix, MatrixExpr, MatrixSymbol,
+        SparseMatrix, Transpose, Adjoint, MatrixSet)
+from sympy.matrices.common import NonSquareMatrixError
 from sympy.matrices.expressions.determinant import Determinant, det
 from sympy.matrices.expressions.matexpr import MatrixElement
 from sympy.matrices.expressions.special import ZeroMatrix, Identity
@@ -50,8 +51,7 @@ def test_matrix_symbol_creation():
 
 def test_matexpr_properties():
     assert A.shape == (n, m)
-    assert (A*B).shape == (n, l)
-    raises(ShapeError, lambda: B*A)
+    assert (A * B).shape == (n, l)
     assert A[0, 1].indices == (0, 1)
     assert A[0, 0].symbol == A
     assert A[0, 0].symbol.name == 'A'
@@ -114,14 +114,12 @@ def test_addition():
     assert (A + B).shape == A.shape
     assert isinstance(A - A + 2*B, MatMul)
 
-    raises(ShapeError, lambda: A + B.T)
     raises(TypeError, lambda: A + 1)
     raises(TypeError, lambda: 5 + A)
     raises(TypeError, lambda: 5 - A)
 
     assert A + ZeroMatrix(n, m) - A == ZeroMatrix(n, m)
-    with raises(TypeError):
-        ZeroMatrix(n,m) + S.Zero
+    raises(TypeError, lambda: ZeroMatrix(n, m) + S.Zero)
 
 
 def test_multiplication():
@@ -130,10 +128,7 @@ def test_multiplication():
     C = MatrixSymbol('C', n, n)
 
     assert (2*A*B).shape == (n, l)
-
     assert (A*0*B) == ZeroMatrix(n, l)
-
-    raises(ShapeError, lambda: B*A)
     assert (2*A).shape == A.shape
 
     assert A * ZeroMatrix(m, m) * B == ZeroMatrix(n, l)
diff --git a/sympy/matrices/expressions/tests/test_matmul.py b/sympy/matrices/expressions/tests/test_matmul.py
index e3126c8ae4..cca168bb0d 100644
--- a/sympy/matrices/expressions/tests/test_matmul.py
+++ b/sympy/matrices/expressions/tests/test_matmul.py
@@ -1,6 +1,7 @@
 from sympy.core import I, symbols, Basic, Mul, S
 from sympy.core.mul import mul
 from sympy.functions import adjoint, transpose
+from sympy.matrices.common import ShapeError
 from sympy.matrices import (Identity, Inverse, Matrix, MatrixSymbol, ZeroMatrix,
         eye, ImmutableMatrix)
 from sympy.matrices.expressions import Adjoint, Transpose, det, MatPow
@@ -12,7 +13,7 @@
 from sympy.assumptions.refine import refine
 from sympy.core.symbol import Symbol
 
-from sympy.testing.pytest import XFAIL
+from sympy.testing.pytest import XFAIL, raises
 
 n, m, l, k = symbols('n m l k', integer=True)
 x = symbols('x')
@@ -177,3 +178,9 @@ def test_issue_23519():
     I = Identity(N)
     z = (M2 + 2 * (M2 + I) * M1 + I)
     assert z.coeff(M1) == 2*I + 2*M2
+
+
+def test_shape_error():
+    A = MatrixSymbol('A', 2, 2)
+    B = MatrixSymbol('B', 3, 3)
+    raises(ShapeError, lambda: MatMul(A, B))
diff --git a/sympy/matrices/expressions/tests/test_special.py b/sympy/matrices/expressions/tests/test_special.py
index 794fb785e6..beeaf1d76a 100644
--- a/sympy/matrices/expressions/tests/test_special.py
+++ b/sympy/matrices/expressions/tests/test_special.py
@@ -6,7 +6,6 @@
 from sympy.concrete.summations import Sum
 from sympy.functions.elementary.complexes import im, re
 from sympy.functions.elementary.piecewise import Piecewise
-from sympy.matrices.common import NonSquareMatrixError, ShapeError
 from sympy.matrices.immutable import ImmutableDenseMatrix
 from sympy.matrices.expressions.matexpr import MatrixSymbol
 from sympy.matrices.expressions.matadd import MatAdd
@@ -130,13 +129,6 @@ def test_ZeroMatrix():
     assert im(Z) == Z
 
     assert ZeroMatrix(n, n)**0 == Identity(n)
-    with raises(NonSquareMatrixError):
-        Z**0
-    with raises(NonSquareMatrixError):
-        Z**1
-    with raises(NonSquareMatrixError):
-        Z**2
-
     assert ZeroMatrix(3, 3).as_explicit() == ImmutableDenseMatrix.zeros(3, 3)
 
 
@@ -151,7 +143,6 @@ def test_ZeroMatrix_doit():
 def test_OneMatrix():
     n, m = symbols('n m', integer=True)
     A = MatrixSymbol('A', n, m)
-    a = MatrixSymbol('a', n, 1)
     U = OneMatrix(n, m)
 
     assert U.shape == (n, m)
@@ -163,14 +154,6 @@ def test_OneMatrix():
     assert im(U) == ZeroMatrix(n, m)
 
     assert OneMatrix(n, n) ** 0 == Identity(n)
-    with raises(NonSquareMatrixError):
-        U ** 0
-    with raises(NonSquareMatrixError):
-        U ** 1
-    with raises(NonSquareMatrixError):
-        U ** 2
-    with raises(ShapeError):
-        a + U
 
     U = OneMatrix(n, n)
     assert U[1, 2] == 1
diff --git a/sympy/physics/units/tests/test_quantities.py b/sympy/physics/units/tests/test_quantities.py
index 962677d4d0..0d6784decf 100644
--- a/sympy/physics/units/tests/test_quantities.py
+++ b/sympy/physics/units/tests/test_quantities.py
@@ -18,7 +18,7 @@
     day, foot, grams, hour, inch, kg, km, m, meter, millimeter,
     minute, quart, s, second, speed_of_light, bit,
     byte, kibibyte, mebibyte, gibibyte, tebibyte, pebibyte, exbibyte,
-    kilogram, gravitational_constant, electron_rest_mass)
+    kilogram, gravitational_constant)
 
 from sympy.physics.units.definitions.dimension_definitions import (
     Dimension, charge, length, time, temperature, pressure,
@@ -278,9 +278,6 @@ def test_issue_quart():
     assert convert_to(4 * quart / inch ** 3, meter) == 231
     assert convert_to(4 * quart / inch ** 3, millimeter) == 231
 
-def test_electron_rest_mass():
-    assert convert_to(electron_rest_mass, kilogram) == 9.1093837015e-31*kilogram
-    assert convert_to(electron_rest_mass, grams) == 9.1093837015e-28*grams
 
 def test_issue_5565():
     assert (m < s).is_Relational
@@ -307,11 +304,11 @@ def test_find_unit():
         'deciliter', 'centiliter', 'deciliters', 'milliliter',
         'centiliters', 'milliliters', 'planck_volume']
     assert find_unit('voltage') == ['V', 'v', 'volt', 'volts', 'planck_voltage']
-    assert find_unit(grams) == ['g', 't', 'Da', 'kg', 'me', 'mg', 'ug', 'amu', 'mmu', 'amus',
-                                'gram', 'mmus', 'grams', 'pound', 'tonne', 'dalton', 'pounds',
-                                'kilogram', 'kilograms', 'microgram', 'milligram', 'metric_ton',
-                                'micrograms', 'milligrams', 'planck_mass', 'milli_mass_unit', 'atomic_mass_unit',
-                                'electron_rest_mass', 'atomic_mass_constant']
+    assert find_unit(grams) == ['g', 't', 'Da', 'kg', 'mg', 'ug', 'amu', 'mmu', 'amus',
+                                'gram', 'mmus', 'grams', 'pound', 'tonne', 'dalton',
+                                'pounds', 'kilogram', 'kilograms', 'microgram', 'milligram',
+                                'metric_ton', 'micrograms', 'milligrams', 'planck_mass',
+                                'milli_mass_unit', 'atomic_mass_unit', 'atomic_mass_constant']
 
 
 def test_Quantity_derivative():