[sympy] Remove auto-evaluation in MatrixExpr (#1670)

Ronan Lamy

unread,

Nov 23, 2012, 3:01:55 PM11/23/12

to sympy/sympy

After the lengthy discussions surrounding sympy.rules and sympy.unify, I've come to the conclusion that the only clean way for them to work is to avoid any form of auto-evaluation in their constructors.

For now, I've only dealt with Adjoint and Transpose. The rest will follow.

NB: this depends on #1667.

You can merge this Pull Request by running:

  git pull https://github.com/rlamy/sympy MatExpr

Or view, comment on, or merge it at:

https://github.com/sympy/sympy/pull/1670

Commit Summary

Make Matrix sympifiable.
Make Transpose and Adjoint not evaluate their argument

File Changes

M sympy/core/decorators.py (6)
M sympy/matrices/dense.py (7)
M sympy/matrices/expressions/adjoint.py (38)
M sympy/matrices/expressions/blockmatrix.py (22)
M sympy/matrices/expressions/matadd.py (11)
M sympy/matrices/expressions/matexpr.py (50)
M sympy/matrices/expressions/matmul.py (12)
M sympy/matrices/expressions/tests/test_adjoint.py (17)
M sympy/matrices/expressions/tests/test_matrix_exprs.py (29)
M sympy/matrices/expressions/transpose.py (36)
M sympy/matrices/immutable.py (17)
M sympy/matrices/matrices.py (11)
M sympy/matrices/tests/test_interactions.py (16)
M sympy/printing/tests/test_latex.py (4)

Patch Links

—
Reply to this email directly or view it on GitHub.

Matthew Rocklin

unread,

Nov 23, 2012, 4:00:50 PM11/23/12

to sympy/sympy

Great. I look forward to seeing the result.

Aaron Meurer

unread,

Nov 24, 2012, 1:04:09 AM11/24/12

to sympy/sympy

SymPy Bot Summary: Passed after merging rlamy/MatExpr (fa1d10d) into master (4ce940a).
Python 2.5.0-final-0: pass
Python 2.6.6-final-0: pass
Python 2.7.2-final-0: pass
Python 2.6.8-final-0: pass
Python 2.7.3-final-0: pass
PyPy 1.9.1-dev-0; 2.7.3-final-42: pass
Python 3.2.2-final-0: pass
Python 3.3.0-final-0: pass
Python 3.2.3-final-0: pass
Python 3.3.0-final-0: pass
Python 3.3.0-final-0: pass
**Sphinx 1.1.3:** pass

Julien Rioux

unread,

Nov 24, 2012, 6:40:30 AM11/24/12

to sympy/sympy

Looks good to me. I suppose we might also want trace to evaluate and Trace would remain unevaluated.

Ronan Lamy

unread,

Nov 24, 2012, 12:16:57 PM11/24/12

to sympy/sympy

Yes, that's part of the idea, but trace is a bit of a mess. Quantum code uses sympy.core.trace and AFAIK there's actually no trace() function anywhere. I started some work on it, but I put it on hold. I'll come back to it later.

Matthew Rocklin

unread,

Nov 24, 2012, 1:19:33 PM11/24/12

to sympy/sympy

Yeah, the Tr/Trace issue isn't great. I'm of the opinion that we shouldn't try to merge them though. A general trace function that branches between the two based on the input might be useful.

Julien Rioux

unread,

Nov 24, 2012, 2:38:43 PM11/24/12

to sympy/sympy

SymPy Bot Summary: Passed after merging rlamy/MatExpr (c4908e0) into master (4ce940a).
PyPy 2.0.0-beta-1; 2.7.3-final-42: pass

Python 2.7.2-final-0: pass

Python 3.2.1-final-0: pass
Sphinx 1.1.3: pass

Julien Rioux

unread,

Nov 25, 2012, 7:57:13 AM11/25/12

to sympy/sympy

SymPy Bot Summary: Passed after merging rlamy/MatExpr (cdfd8d6) into master (6084fa8).

PyPy 2.0.0-beta-1; 2.7.3-final-42: pass
Python 2.7.2-final-0: pass
Python 3.2.1-final-0: pass
Sphinx 1.1.3: pass

—

Reply to this email directly or view it on GitHub.

Julien Rioux

unread,

Nov 27, 2012, 9:07:16 AM11/27/12

to sympy/sympy

This establishes a very nice distinction between lower-case dothing and upper-case Thing. I started adding more tests, that I will submit as a PR to you. How far did you get on reworking trace/Trace/Tr?

Ronan Lamy

unread,

Nov 27, 2012, 11:31:27 AM11/27/12

to sympy/sympy

Thanks for the tests. I didn't actually do much about the trace - I tried to make the obvious changes in Trace.__new__ but didn't fix the resulting issues.

Ronan Lamy

unread,

Dec 5, 2012, 3:11:42 PM12/5/12

to sympy/sympy

This should be ready for merging now. In addition to the auto-evaluation thing, I have also tried to ensure that obj.__class__(*obj.args) == obj holds for all MatrixExprs, but I didn't deal with the trace.

Matthew Rocklin

unread,

Dec 5, 2012, 3:19:24 PM12/5/12

to sympy/sympy

In sympy/matrices/expressions/inverse.py:

>      B^-1*A^-1
> +    >>> (A*B).inverse() == Inverse(A*B)
> +    False

While this is important to understand I think that this line will be confusing to non-experts. I think this sort of thing belongs in documentation but not in docstrings. If you really think it should be here then I think that this issue should be explained.

Matthew Rocklin

unread,

Dec 5, 2012, 3:25:53 PM12/5/12

to sympy/sympy

In sympy/matrices/expressions/adjoint.py:

> @@ -1,33 +1,36 @@
> -from matexpr import MatrixExpr
> -from sympy import Basic
> -from sympy.functions.elementary.complexes import conjugate
> +from sympy.core import Basic
> +from sympy.functions import adjoint, conjugate, transpose

These inherit from Expr - does this make sense?

Matthew Rocklin

unread,

Dec 5, 2012, 3:31:23 PM12/5/12

to sympy/sympy

Are you sure that canonicalization runs deep?

In [21]: X = MatrixSymbol('X', 3, 3)
In [22]: transpose(X+X)
Out[22]: X'⋅2

Matthew Rocklin

unread,

Dec 5, 2012, 3:33:20 PM12/5/12

to sympy/sympy

In sympy/matrices/expressions/blockmatrix.py:

>  
> -    def inv(self, expand=False):

This was inv to match the Matrix API. I'm fine with the new version but some users might find it disconcerting.

Matthew Rocklin

unread,

Dec 5, 2012, 3:39:06 PM12/5/12

to sympy/sympy

In sympy/matrices/expressions/inverse.py:

>  
>      >>> from sympy import MatrixSymbol, Inverse
>      >>> A = MatrixSymbol('A', 3, 3)
>      >>> B = MatrixSymbol('B', 3, 3)
>      >>> Inverse(A)
>      A^-1
> -    >>> A.I

I prefer the A.I syntax. It keeps to the numpy standard

Matthew Rocklin

unread,

Dec 5, 2012, 3:44:33 PM12/5/12

to sympy/sympy

In sympy/matrices/expressions/tests/test_matrix_exprs.py:

> @@ -341,3 +177,14 @@ def test_matmul_simplify():
>      A = MatrixSymbol('A', 1, 1)
>      assert simplify(MatMul(A, ImmutableMatrix([[sin(x)**2 + cos(x)**2]]))) == \
>          MatMul(A, ImmutableMatrix([[1]]))
> +
> +def test_invariants():
> +    A = MatrixSymbol('A', n, m)
> +    B = MatrixSymbol('B', m, l)
> +    X = MatrixSymbol('X', n, n)
> +    objs = [Identity(n), ZeroMatrix(m, n), A, MatMul(A, B), MatAdd(A, A),
> +            Transpose(A), Adjoint(A), Inverse(X), MatPow(X, 2), MatPow(X, -1),
> +            MatPow(X, 0)]
> +    for obj in objs:
> +        print obj

print

Matthew Rocklin

unread,

Dec 5, 2012, 3:46:45 PM12/5/12

to sympy/sympy

In sympy/matrices/expressions/tests/test_hadamard.py:

>      assert isinstance(expr2, MatrixSymbol)
> +
> +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)
> +    with raises(TypeError):
> +        hadamard_product()
> +    assert hadamard_product(A) == A
> +    assert isinstance(hadamard_product(A, B), HadamardProduct)
> +    assert hadamard_product(A, B).doit() == hadamard_product(A, B)
> +    assert hadamard_product(A, B) == hadamard_product(B, A)
> +    with raises(ShapeError):
> +        hadamard_product(A, C)

Neat.

Matthew Rocklin

unread,

Dec 5, 2012, 3:48:58 PM12/5/12

to sympy/sympy

In sympy/matrices/expressions/matadd.py:

> @@ -18,21 +19,12 @@ class MatAdd(MatrixExpr):
>      is_MatAdd = True
>  
>      def __new__(cls, *args, **kwargs):
> -        evaluate = kwargs.get('evaluate', True)
> -        check    = kwargs.get('check'   , True)
> -
> -        # TODO: This is a kludge
> -        # We still use Matrix + 0 in a few places. This removes it
> -        # In particular see matrix_multiply
> -        args = [x for x in args if not x == 0]

Did you run across any issues when removing this? It was actually necessary at some point because (if I recall correctly) BlockMatrix used Matrix routines which assumed that they could sum elements onto scalar 0.

Matthew Rocklin

unread,

Dec 5, 2012, 4:07:20 PM12/5/12

to sympy/sympy

doit doesn't seem quite appropriate here. Often doit doesn't do anything. I wonder if there isn't a better choice of method name.

Ronan Lamy

unread,

Dec 5, 2012, 4:59:15 PM12/5/12

to sympy/sympy

fixed

Ronan Lamy

unread,

Dec 5, 2012, 5:53:51 PM12/5/12

to sympy/sympy

In sympy/matrices/expressions/inverse.py:

>      B^-1*A^-1
> +    >>> (A*B).inverse() == Inverse(A*B)
> +    False

Well, non-experts should not look at Inverse in the first place, but just use .I/.inverse(). If they do look at it, this is really the most important thing they should know. What do you want me to add to the paragraph above?

Ronan Lamy

unread,

Dec 5, 2012, 6:00:38 PM12/5/12

to sympy/sympy

In sympy/matrices/expressions/inverse.py:

>  
>      >>> from sympy import MatrixSymbol, Inverse
>      >>> A = MatrixSymbol('A', 3, 3)
>      >>> B = MatrixSymbol('B', 3, 3)
>      >>> Inverse(A)
>      A^-1
> -    >>> A.I

The problem with A.I is that it uses attribute syntax but can trigger expensive computations, so it doesn't follow the recommendation that properties should perform simple and fast operations. I didn't remove it but I think that we should not encourage its use.

Ronan Lamy

unread,

Dec 5, 2012, 6:06:35 PM12/5/12

to sympy/sympy

In sympy/matrices/expressions/matadd.py:

> @@ -18,21 +19,12 @@ class MatAdd(MatrixExpr):
>      is_MatAdd = True
>  
>      def __new__(cls, *args, **kwargs):
> -        evaluate = kwargs.get('evaluate', True)
> -        check    = kwargs.get('check'   , True)
> -
> -        # TODO: This is a kludge
> -        # We still use Matrix + 0 in a few places. This removes it
> -        # In particular see matrix_multiply
> -        args = [x for x in args if not x == 0]

Yes, that was the problem, but fixing it was just a matter of handling empty matrices separately. Cf. df17dac.

Ronan Lamy

unread,

Dec 5, 2012, 6:14:42 PM12/5/12

to sympy/sympy

In sympy/matrices/expressions/adjoint.py:

> @@ -1,33 +1,36 @@
> -from matexpr import MatrixExpr
> -from sympy import Basic
> -from sympy.functions.elementary.complexes import conjugate
> +from sympy.core import Basic
> +from sympy.functions import adjoint, conjugate, transpose

Well, they are only used as functions (in the Python sense), and I think they should actually be turned into functions, some day.

Ronan Lamy

unread,

Dec 5, 2012, 6:23:13 PM12/5/12

to sympy/sympy

In sympy/matrices/expressions/blockmatrix.py:

>  
> -    def inv(self, expand=False):

Generic MatrixExprs don't have .inv() so having it here wasn't really consistent. Also the distinction between .inv() and .inverse() was confusing.

Julien Rioux

unread,

Dec 7, 2012, 1:17:11 AM12/7/12

to sympy/sympy

SymPy Bot Summary: Passed after merging rlamy/MatExpr (5044db6) into master (02bc3d5).

PyPy 2.0.0-beta-1; 2.7.3-final-42: pass
Python 2.7.2-final-0: pass
Python 3.2.1-final-0: pass
Sphinx 1.1.3: pass

—

Reply to this email directly or view it on GitHub.

Julien Rioux

unread,

Dec 11, 2012, 4:08:01 PM12/11/12

to sympy/sympy

In sympy/matrices/expressions/adjoint.py:

>  
> -    Represents the Adjoint of a matrix expression.
> +class Adjoint(MatrixExpr):
> +    """
> +    The Hermitian adjoint of a matrix expression.

Could you please add the same "warning" as in the Inverse and Transpose docstrings i.e. that this is unevaluated.

Julien Rioux

unread,

Dec 11, 2012, 4:09:10 PM12/11/12

to sympy/sympy

In sympy/matrices/expressions/hadamard.py:

>  class HadamardProduct(MatrixExpr):
> -    """Elementwise Product of Matrix Expressions
> +    """
> +    Elementwise product of matrix expressions

Another unevaluated "warning".

Julien Rioux

unread,

Dec 11, 2012, 4:10:54 PM12/11/12

to sympy/sympy

In sympy/matrices/expressions/inverse.py:

>      B^-1*A^-1
> +    >>> (A*B).inverse() == Inverse(A*B)
> +    False

More useful might be to actually show how they differ.

>>> Inverse(A*B)
(A*B)^-1
>>> (A*B).inverse()
B^-1*A^-1

Julien Rioux

unread,

Dec 11, 2012, 4:15:46 PM12/11/12

to sympy/sympy

This looks about ready to merge.

@mrocklin What do you mean with your canonicalization comment?

Matthew Rocklin

unread,

Dec 11, 2012, 4:29:46 PM12/11/12

to sympy/sympy

What do you mean with your canonicalization comment?

Current behavior is

In [21]: X = MatrixSymbol('X', 3, 3)
In [22]: transpose(X+X)
Out[22]: X'⋅2

Ideal behavior is

In [21]: X = MatrixSymbol('X', 3, 3)
In [22]: transpose(X+X)
Out[22]: 2⋅X'

Christopher Smith

unread,

Dec 11, 2012, 4:36:56 PM12/11/12

to sympy/sympy

Current behavior in master is

Christopher Smith

unread,

Dec 11, 2012, 4:37:11 PM12/11/12

to sympy/sympy


>>> X = MatrixSymbol('X', 3, 3)
>>> transpose(X+X)
2*X'

Julien Rioux

unread,

Dec 11, 2012, 4:50:17 PM12/11/12

to sympy/sympy

In sympy/matrices/expressions/matmul.py:

>  
>      def _eval_transpose(self):
> -        from transpose import Transpose
> -        return MatMul(*[Transpose(arg) for arg in self.args[::-1]])
> +        return MatMul(*[transpose(arg) for arg in self.args[::-1]])

This is the line that causes the behavior

>>> transpose(2*X)
X'*2

something better could be done, maybe using as_coeff_matrices.

Christopher Smith

unread,

Dec 12, 2012, 12:15:19 AM12/12/12

to sympy/sympy

In sympy/matrices/expressions/matmul.py:

>  
>      def _eval_transpose(self):
> -        from transpose import Transpose
> -        return MatMul(*[Transpose(arg) for arg in self.args[::-1]])
> +        return MatMul(*[transpose(arg) for arg in self.args[::-1]])

On Wed, Dec 12, 2012 at 3:35 AM, Julien Rioux <notifi...@github.com>wrote: In sympy/matrices/expressions/matmul.py: > > def _eval_transpose(self): > - from transpose import Transpose > - return MatMul(*[Transpose(arg) for arg in self.args[::-1]]) > + return MatMul(*[transpose(arg) for arg in self.args[::-1]]) This is the line that causes the behavior >>> transpose(2*X) X'*2 something better could be done, maybe using as_coeff_matrices.

hmmm...it's traversing in reversed order. So if you separate off non-matrix args first, that should do it, something like ``` h, t = self.as_independent(Matrix, as_Add=False) return h*MatMul(*[transpose(a) for a in Mul.make_args(t)]) ```

Julien Rioux

unread,

Dec 16, 2012, 7:38:56 PM12/16/12

to sympy/sympy

SymPy Bot Summary: Passed after merging rlamy/MatExpr (748db07) into master (501534b).

PyPy 2.0.0-beta-1; 2.7.3-final-42: pass
Python 2.7.2-final-0: pass
Python 3.2.1-final-0: pass
Sphinx 1.1.3: pass

Docs build command: make clean && make html-errors && make latex && cd _build/latex && xelatex sympy-*.tex

Julien Rioux

unread,

Dec 19, 2012, 1:38:43 PM12/19/12

to sympy/sympy

SymPy Bot Summary: Passed after merging rlamy/MatExpr (e339511) into master (d503614).

PyPy 2.0.0-beta-1; 2.7.3-final-42: pass
Python 2.7.2-final-0: pass
Python 3.2.1-final-0: pass
Sphinx 1.1.3: pass
Docs build command: make clean && make html-errors && make latex && cd _build/latex && xelatex sympy-*.tex

Ronan Lamy

unread,

Dec 19, 2012, 11:56:33 PM12/19/12

to sympy/sympy

In sympy/matrices/expressions/inverse.py:

>      B^-1*A^-1
> +    >>> (A*B).inverse() == Inverse(A*B)
> +    False

done

Ronan Lamy

unread,

Dec 19, 2012, 11:56:49 PM12/19/12

to sympy/sympy

In sympy/matrices/expressions/adjoint.py:

>  
> -    Represents the Adjoint of a matrix expression.
> +class Adjoint(MatrixExpr):
> +    """
> +    The Hermitian adjoint of a matrix expression.

Ronan Lamy

unread,

Dec 19, 2012, 11:57:03 PM12/19/12

to sympy/sympy

In sympy/matrices/expressions/hadamard.py:

>  class HadamardProduct(MatrixExpr):
> -    """Elementwise Product of Matrix Expressions
> +    """
> +    Elementwise product of matrix expressions

Ronan Lamy

unread,

Dec 20, 2012, 12:02:17 AM12/20/12

to sympy/sympy

Are there any other comments? I think this is ready now.

Matthew Rocklin

unread,

Dec 20, 2012, 8:57:35 AM12/20/12

to sympy/sympy

Seems good to me. If you wait on merging until the weekend I can give it another look. I'd be fine merging now though.

Ronan Lamy

unread,

Dec 20, 2012, 8:40:11 PM12/20/12

to sympy/sympy

Thanks! I'd rather merge now, so if you see anything to change, please send another PR.

Ronan Lamy

unread,

Dec 20, 2012, 8:40:27 PM12/20/12

to sympy/sympy

Merged #1670.

Reply all

Reply to author

Forward