Simplify on Non-Exprs

[Matthew Rocklin]

Should we allow

x = Matrix(....)

y = simplify(x)  ?

How about

x = ImmutableMatrix(...)

y = simplify(x)  ?

Historically the first has never worked but the second did work for the last few months because of an odd inheritance tree. I'm currently proposing a PR which removes the ImmutableMatrix inherits from Expr link so simplify(ImmutableMatrix) will cease to function. Also we have always had, (and will continue to have

ImmutableMatrix.simplify()

The core issue is that simplify was written with scalar expressions in mind. It (indirectly) calls methods like `as_numer_denom` that assumes that what you're working with is a scalar. This has worked on ImmutableMatrix in the past because ImmutableMatrix inherited from Expr in order to get a few nice features.  A lot of things sort of worked and covered up a few problems. In a new PR I cut this inheritance so this question opens up again. 

I see two options:

1. Do not allow simplify(x) syntax where x is not an Expr

2. Add a shortcut `if hasattr(x, '_simplify'): return x._simplify()` to the top of simplify

3. Someone else suggests something better than 1 or 2

Do we want the function simplify to be just for exprs or do we want it to be for all sorts of things? If we want it to be for all sorts of things what is the best way to handle the polymorphism?

[Julien Rioux]

Hi,

I certainly think I should be able to simplify(m) for m a matrix type. As I use sympy I expect that all objects interact together naturally. I should be able to combine objects (as long as they make mathematical sense) and be able to use the same functions on all objects. It used to be that simplify wasn't built to deal very well with noncommutative symbols, so e.g. simplify(A*B - B*A) would return 0. Those were considered bugs and were eventually fixed. I think the same approach should be taken about simplifying matrices, any issue should be considered a bug and fixed. I stopped using Matrix and have since used ImmutableMatrix because I could use most regular sympy functions on it. By the way, I wish I could do factor(m), too. Although it doesn't do anything useful at the moment, we should improve upon it.

I think the inheritance on Expr ensures that a lot of functions such as simplify will work (possibly for free) or at least they won't error out. It doesn't bother me too much how it's implemented as long as it's not overly complicated and it works.

Regards,
Julien

[Matthew Rocklin]

On Mon, Oct 15, 2012 at 2:16 PM, Julien Rioux <julien...@gmail.com> wrote:

Hi,

I certainly think I should be able to simplify(m) for m a matrix type. As I use sympy I expect that all objects interact together naturally. I should be able to combine objects (as long as they make mathematical sense) and be able to use the same functions on all objects.

This is a good ideal. Unfortunately most of the Expr code was written to only work with scalars. I.e. people often use 0 when they mean "additive identity" (I use this example a lot, please forgive me). I don't think that incremental bug fixes will solve this problem.

It used to be that simplify wasn't built to deal very well with noncommutative symbols, so e.g. simplify(A*B - B*A) would return 0. Those were considered bugs and were eventually fixed. I think the same approach should be taken about simplifying matrices, any issue should be considered a bug and fixed. 

I stopped using Matrix and have since used ImmutableMatrix because I could use most regular sympy functions on it. By the way, I wish I could do factor(m), too. Although it doesn't do anything useful at the moment, we should improve upon it.

I think that the easiest way to improve upon it is to have factor call ImmutableMatrix._factor. Trying to make one giant function to do factoring for all classes sounds difficult to me. It's hard to write a single function that factors all sorts of things. 

I think the inheritance on Expr ensures that a lot of functions such as simplify will work (possibly for free) or at least they won't error out. 

This is not my experience. It is very easy to stress the system and get into very weird states; or at least it was until I monkey patched some things. Now it is only moderately easy :) In these cases I would rather error out at a higher level than at a very low level. My experience is that the MatrixExpr(Expr) inheritance is very hard to reason about and has high potential for bugs. I've duct-taped things together so that they appear to usually work. It's very hard to develop new functionality with this though. I always end up running into Expr trouble. MatrixExpr can't grow effectively and will have many more bugs while it is still connected to Expr.

 

> It doesn't bother me too much how it's implemented as long as it's not overly complicated and it works.

I think the easiest way to make it work is to remove the Expr dependence and use classes for polymorphism (i.e. put `if hasattr(expr, '_simplify'): return expr._simplify()` in the publicly visible simplify funciton. )

[Aaron Meurer]

I think the real question here is, do we want most functions in sympy
to automatically vectorize over matrices (or similar objects)? And I
think the answer is no. Consider for example if we had an unevaluated
Matrix type that didn't multiply through by default (what
ImmutableMatrix used to be). We might want something like
simplify(A*B + A*C) to give A*(B + C), where A, B, and C are all
explicit unevaluated matrices.

It's similar to the reason that solve() does not automatically call
dsolve() on ODEs. It's useful to be able solve for the function
algebraically, i.e., both

>>> dsolve(f(x).diff(x, x) + f(x) - 1, f(x))
f(x) == C1*sin(x) + C2*cos(x) + 1

and

>>> solve(f(x).diff(x, x) + f(x) - 1, f(x))
[-Derivative(f(x), x, x) + 1]

can be useful.

Another good argument for matrices is exp(). We want exp(Matrix) to
give the exponential of a Matrix, which is not the matrix of
exponentials. And actually you can define f(Matrix) for any analytic
f. Won't it be a little confusing if math_function(Matrix) works
different from simplify_function(Matrix)? And again, what if you have
simplify_function(combination_of_math_functions(unevaluated_matrices))?

And aside from that, automatic vectorization would be a mess, because
only the most popular functions would support it, and the rest
wouldn't. So we would end up with a situation where function `a`
works as expected, but very similar function `b` doesn't, and it's not
clear why (and even if `b` doesn't support it, it doesn't mean that
Matrix defines `_b`). It's similar to automatic sympification of
strings that some but not all SymPy functions support.

Aaron Meurer

> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To post to this group, send email to sy...@googlegroups.com.
> To unsubscribe from this group, send email to
> sympy+un...@googlegroups.com.
> For more options, visit this group at
> http://groups.google.com/group/sympy?hl=en.

[Matthew Rocklin]

So what should we do with `simplify(ImmutableMatrix)`. What is the preferred behavior?

[Julien Rioux]

Hi,

On Tuesday, 16 October 2012 07:48:52 UTC+2, Aaron Meurer wrote:

I think the real question here is, do we want most functions in sympy
to automatically vectorize over matrices (or similar objects)?  And I
think the answer is no.

I'm not sure how you end up at this question, but fine. The answer is no. I'm fine with that.

 

> Consider for example if we had an unevaluated
Matrix type that didn't multiply through by default (what
ImmutableMatrix used to be).  We might want something like
simplify(A*B + A*C) to give A*(B + C), where A, B, and C are all
explicit unevaluated matrices.

Sure, we might want that. Simplifying a MatAdd or MatMul should be implemented, too. Of course what we mean by "simpler" is always ambiguous. There might be different "strategies", but simplify is the general one. The point is that it actually works right now to some extend:

In [9]: A = ImmutableMatrix(2, 2, [x*y + x*z, 0, 0, 0])

In [10]: B = ImmutableMatrix(2, 2, [x*y + x*z, 0, 0, 1])

In [11]: C = ImmutableMatrix(2, 2, [x*y + x*z, 1, 1, 0])

In [12]: A*B + A*C  # automatically evaluated
Out[12]:
[             2           ]
[2*(x*y + x*z)   x*y + x*z]
[                         ]
[      0             0    ]

In [13]: simplify(A*B + A*C)
Out[13]:
[   2        2           ]
[2*x *(y + z)   x*(y + z)]
[                        ]
[      0            0    ]

In [14]: MatAdd(MatMul(A, B), MatMul(A, C))  # unevaluated
Out[14]:
[x*y + x*z  0]*[x*y + x*z  0] + [x*y + x*z  0]*[x*y + x*z  1]
[            ] [            ]   [            ] [            ]
[    0      0] [    0      1]   [    0      0] [    1      0]

In [15]: simplify(MatAdd(MatMul(A, B), MatMul(A, C)))
Out[15]:
[x*(y + z)  0]*[x*(y + z)  0] + [x*(y + z)  0]*[x*(y + z)  1]
[            ] [            ]   [            ] [            ]
[    0      0] [    0      1]   [    0      0] [    1      0]

Should it work better? yes absolutely. We see that it doesn't factor out the identical part. But why should it be removed?
 

It's similar to the reason that solve() does not automatically call
dsolve() on ODEs.  It's useful to be able solve for the function
algebraically, i.e., both

>>> dsolve(f(x).diff(x, x) + f(x) - 1, f(x))
f(x) == C1*sin(x) + C2*cos(x) + 1

and

>>> solve(f(x).diff(x, x) + f(x) - 1, f(x))
[-Derivative(f(x), x, x) + 1]

can be useful.

So you mean that solving a differential equation for f(x) is ambiguous. At least here, we can identify two clearly distinct strategies, and define a function for each approach. It sure can be useful to have both. There are such different "strategies" in the simplify module. One could think of strategies specific to matrices, I guess. simplify itself has been the wrapper around those strategies.
 

Another good argument for matrices is exp().  We want exp(Matrix) to
give the exponential of a Matrix, which is not the matrix of
exponentials. And actually you can define f(Matrix) for any analytic
f.

To map mathematical functions to each element of the matrix does not make mathematical sense to me. If we were talking about an Array type, maybe.
 

 Won't it be a little confusing if math_function(Matrix) works
different from simplify_function(Matrix)?

At least not to me. These are vastly different things.
 

 And again, what if you have
simplify_function(combination_of_math_functions(unevaluated_matrices))?

I suppose what would happen is something akin to

In [16]: sin(x*y + x*z)
Out[16]: sin(x*y + x*z)

In [17]: simplify(sin(x*y + x*z))
Out[17]: sin(x*(y + z))

Well yes it is in fact what happens:

In [23]: sin(A)  
Out[23]:
sin/[x*y + x*z  0]\
   |[            ]|
   \[    0      0]/

In [24]: simplify(sin(A))
Out[24]:
sin/[x*(y + z)  0]\
   |[            ]|
   \[    0      0]/

Now I can't use sin(MatAdd), that does not currently work, but I would expect it to, and calling simplify on that would, you know, simplify the argument of the sin().

 

And aside from that, automatic vectorization would be a mess, because
only the most popular functions would support it, and the rest
wouldn't.  So we would end up with a situation where function `a`
works as expected, but very similar function `b` doesn't, and it's not
clear why (and even if `b` doesn't support it, it doesn't mean that
Matrix defines `_b`).  It's similar to automatic sympification of
strings that some but not all SymPy functions support.

Aaron Meurer

I don't want automatic map of a function to all elements of a matrix, either. You introduced that idea.

Regards,
Julien

[Matthew Rocklin]

I agree with Julien that vectorization is a different question. The question at hand is (I think) the following:

We have a function simplify that tries to reduce the complexity of a scalar expression (for example by turning  y* x / y into x). We have grown accustomed to asking SymPy to reduce the complexity of an object by calling `simplify(object)`. Should we extend the meaning of simplify to non-Exprs? Many classes know some way to reduce their complexity. For example matrices have a Matrix.simplify() method. Some sets like Union and Intersection have a Set.reduce() method. Should we expose these methods through simplify? I.e. do we want things like the following to work `simplify(Set)` ? 

--

You received this message because you are subscribed to the Google Groups "sympy" group.

To view this discussion on the web visit https://groups.google.com/d/msg/sympy/-/-eALL2LJk7kJ.

[Sergiu Ivanov]

Hello,

On Tue, Oct 16, 2012 at 08:15:12AM -0500, Matthew Rocklin wrote:
>
> I agree with Julien that vectorization is a different question. The
> question at hand is (I think) the following:
>
> We have a function simplify that tries to reduce the complexity of a scalar
> expression (for example by turning y* x / y into x). We have grown
> accustomed to asking SymPy to reduce the complexity of an object by calling
> `simplify(object)`. Should we extend the meaning of simplify to non-Exprs?
> Many classes know some way to reduce their complexity. For example matrices
> have a Matrix.simplify() method. Some sets like Union and Intersection have
> a Set.reduce() method. Should we expose these methods through simplify?
> I.e. do we want things like the following to work `simplify(Set)` ?

I think `simplify(Set)` should work.

I do believe however that the simplification functionality itself
should be provided by the object on which `simplify` is invoked,
rather than stuffed into `simplify`. (I think the majority agree that
this is way to go, so I'm just redundantly restating the fact.)

I think letting `simplify` work on all kinds of SymPy objects is the
behaviour a user would expect.

Sergiu

[Matthew Rocklin]

This is implemented in the last few commits of 

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

def simplify(expr, ...):

    try:

        return expr._eval_simplify(ratio=ratio, measure=measure)

    except AttributeError:

        pass

--
You received this message because you are subscribed to the Google Groups "sympy" group.

[Matthew Rocklin]

On this PR https://github.com/sympy/sympy/pull/1598 the question of sqrt came up. Do we want an _eval_sqrt method? Should sqrt work on Matrices? How far do we want to take this idea?