Powers of non-square matrices

[jlh]

Hello everyone!

I noticed that powers of non-square matrices (with MatPow) do not raise an error, e.g. if A = MatrixSymbol('A', 2, 3) then MatPow(A, 10) does not raise. I would assume that this is just an oversight, but then I see explicit unit tests for MatPow(A, 1), which make me wonder if perhaps there is a reason this is supported. However the behavior is not consistent, since A ** 1 does raise an error anyway. But then again ZeroMatrix(2, 3) ** 1 does not raise.

As far as I know most textbooks do not define A ** n at all for any value of n when A is not square. This might also be of interest: https://math.stackexchange.com/questions/101590/what-is-zero-power-of-a-non-square-matrix

So question: Would it be an acceptable breaking change to completely forbid powers of non-square matrices, including the barely useful MatPow(A, 0) and MatPow(A, 1)? It doesn't seem to break any tests (other than test_matpow.py).

Background: I plan to do some clean-up in the MatPow code to make it follow the same convention used in other places, e.g. move the argument validation from MatrixExpr.__pow__() to MatPow.__new__() and move some specialized code from MatPow.doit() to the _eval_power() methods of respective sub-classes.

Thanks!

[Oscar Benjamin]

Yes, that sounds like a reasonable change.

> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/8b817360-def8-43c9-84a0-bfd59b3c6651%40googlegroups.com.

[Aaron Meurer]

IMO A**1 does make sense, so I don't see any reason to disallow it. I
agree the others are just an oversight.

> Background: I plan to do some clean-up in the MatPow code to make it follow the same convention used in other places, e.g. move the argument validation from MatrixExpr.__pow__() to MatPow.__new__() and move some specialized code from MatPow.doit() to the _eval_power() methods of respective sub-classes.

That sounds good, but be aware that the design of things in the matrix
module is a bit different, in particular, the matrix expressions
remain completely unevaluated by default and only evaluate on doit().
Your redesign might not be affected by that, but it's worth noting
that just because something happens some way in the core, it shouldn't
necessarily happen the same way in the matrix expressions.

Aaron Meurer

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAHVvXxS3try%3DYC14%3DnAto1%2B3j7P4qPuvAiS4CpSC2R8%2BeqSWWg%40mail.gmail.com.

[jlh]

On Friday, April 10, 2020 at 6:55:32 PM UTC+2, Aaron Meurer wrote:

IMO A**1 does make sense, so I don't see any reason to disallow it. I 
agree the others are just an oversight. 

Technically even A ** 0 makes sense, which would just return Identity(A.rows) (and funnily return a matrix of a different shape than the power's base).

My reasoning for removing were:

- It's not consistent: A ** 1 vs MatPow(A, 1)

- Probably nobody relies on that behavior

- Makes the code more complex for hardly any gain

I can keep that behavior in if desired of course, but then it should at least be made consistent, which means it should work for A ** 1 as well as for Matrix([[1, 2]]) ** 1. And A ** 0 should work too.

That sounds good, but be aware that the design of things in the matrix
module is a bit different, in particular, the matrix expressions
remain completely unevaluated by default and only evaluate on doit().
Your redesign might not be affected by that, but it's worth noting
that just because something happens some way in the core, it shouldn't
necessarily happen the same way in the matrix expressions.  

Yeah I've seen they don't evaluate when constructed. They evaluate when calling doit(), which also happens when operators are used with those objects.

[Aaron Meurer]

On Fri, Apr 10, 2020 at 12:39 PM jlh <j...@gmx.ch> wrote:
>
> On Friday, April 10, 2020 at 6:55:32 PM UTC+2, Aaron Meurer wrote:
>>
>> IMO A**1 does make sense, so I don't see any reason to disallow it. I
>> agree the others are just an oversight.
>
>
> Technically even A ** 0 makes sense, which would just return Identity(A.rows) (and funnily return a matrix of a different shape than the power's base).

Why rows and not columns though? I think A**0 should give an error.

>
> My reasoning for removing were:
> - It's not consistent: A ** 1 vs MatPow(A, 1)
> - Probably nobody relies on that behavior
> - Makes the code more complex for hardly any gain

Well I wouldn't be surprised if it breaks something somewhere if A**1
stops working, since something might assume that A**1 is always the
same as A. I guess you can try removing it and seeing what tests fail.
But to me, mathematically, A**1 does make sense, even if none of the
other powers do.

Aaron Meurer

>
> I can keep that behavior in if desired of course, but then it should at least be made consistent, which means it should work for A ** 1 as well as for Matrix([[1, 2]]) ** 1. And A ** 0 should work too.
>
>> That sounds good, but be aware that the design of things in the matrix
>> module is a bit different, in particular, the matrix expressions
>> remain completely unevaluated by default and only evaluate on doit().
>> Your redesign might not be affected by that, but it's worth noting
>> that just because something happens some way in the core, it shouldn't
>> necessarily happen the same way in the matrix expressions.
>
>
> Yeah I've seen they don't evaluate when constructed. They evaluate when calling doit(), which also happens when operators are used with those objects.
>

> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/c697efff-ded3-4412-b727-a0fc69150e3b%40googlegroups.com.

[jlh]

On Friday, April 10, 2020 at 8:48:26 PM UTC+2, Aaron Meurer wrote:

On Fri, Apr 10, 2020 at 12:39 PM jlh <j...@gmx.ch> wrote:
> Technically even A ** 0 makes sense, which would just return Identity(A.rows) (and funnily return a matrix of a different shape than the power's base).

Why rows and not columns though? I think A**0 should give an error. 

Actually it should be GenericIdentity(), which is also what is returned from MatMul() (with no arguments), which is pretty much equivalent to A ** 0. I don't think it should be an error, it's really the same as n ** 0 being 1.

 

Well I wouldn't be surprised if it breaks something somewhere if A**1
stops working, since something might assume that A**1 is always the
same as A. I guess you can try removing it and seeing what tests fail.
But to me, mathematically, A**1 does make sense, even if none of the
other powers do.

I tried it, it didn't break any tests other than the ones in test_matpow.py explicitly testing for that.

[Samuel Lelièvre]

Fri 2020-04-10 20:40:07 UTC, jlh:

I would not say A**0 or A**1 make sense. I would say

- the context for powers is a multiplicative semigroup

- nonsquare matrices therefore don't qualify

 

[S.Y. Lee]

But I think that the conservative change is making MatPow(A, 1) and MatPow(A, 0) unevaluated objects that doesn't canonicalize to anything.

And I also find the idea of generic identity confusing.