Clean up of MatPow code, forbid non-square matrix powers

[jlherren]

Completely forbid matrix powers with non-square bases. Further this moves responsibility of evaluating matrix powers to the respective sub-classes via _eval_power().

Other comments

Also includes a fix for Inverse(Inverse(C)).doit()

Release Notes

• matrices
  • Powers of non-square matrices now always raise an error, even when constructed using MatPow(...).

[sympy-bot]

✅

Hi, I am the SymPy bot (v158). I'm here to help you write a release notes entry. Please read the guide on how to write release notes.

Your release notes are in good order.

Here is what the release notes will look like:

• matrices
  • Powers of non-square matrices now always raise an error, even when constructed using MatPow(...). (#19104 by @jlherren)

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

Note: This comment will be updated with the latest check if you edit the pull request. You need to reload the page to see it.

Click here to see the pull request description that was parsed.

Completely forbid matrix powers with non-square bases. Further this moves responsibility of evaluating matrix powers to the respective sub-classes via _eval_power().

#### Other comments
Also includes a fix for Inverse(Inverse(C)).doit()

#### Release Notes
<!-- BEGIN RELEASE NOTES -->
* matrices
  * Powers of non-square matrices now always raise an error, even when constructed using `MatPow(...)`.
<!-- END RELEASE NOTES -->

Update

The release notes on the wiki have been updated.

[codecov]

Codecov Report

Merging #19104 into master will decrease coverage by 0.016%.
The diff coverage is 92.592%.

@@              Coverage Diff              @@
##            master    #19104       +/-   ##
=============================================
- Coverage   75.786%   75.770%   -0.017%     
=============================================
  Files          650       650               
  Lines       168933    168931        -2     
  Branches     39831     39826        -5     
=============================================
- Hits        128029    128000       -29     
- Misses       35332     35352       +20     
- Partials      5572      5579        +7     

[asmeurer]

You can add a release notes entry for this.

[oscarbenjamin]

With this we have:

In [14]: MatPow(A, 1)                                                                                                             
Out[14]: 
 1
A 

In [15]: MatPow(A, 1).doit()                                                                                                      
Out[15]: 
 1
A 

In [16]: A**1                                                                                                                     
---------------------------------------------------------------------------
NonSquareMatrixError

I think that it's important that the outcome of MatPow(A, n).doit() is the same as A ** n (which would ideally be implemented simply as MatPow(...).doit()).

Is it possible to raise NonSquareMatrixError on MatPow(A, 0) and MatPow(A, 1) for non-square A or would that break things?

[oscarbenjamin]

Is it possible to raise NonSquareMatrixError on MatPow(A, 0) and MatPow(A, 1) for non-square A or would that break things?

Alternatively could we raise on doit?

[jlherren]

@oscarbenjamin Actually that was my original intent, but see https://groups.google.com/forum/#!topic/sympy/o622ArN_F1g

However raising an error in doit() is maybe an option.

[oscarbenjamin]

It would also be good to have tests showing that the equivalents of these with MatPow do not evaluate to anything e.g.:

from sympy.core.expr import unchanged

C2 = MatPow(C, 2)
assert unchanged(MatPow, C, 0)
assert unchanged(MatPow, C, 1)
assert unchanged(MatPow, C, -1)
assert unchanged(MatPow, C2, 2)

[oscarbenjamin]

It seems there are no tests at all like that in the MatPow tests. The idea of unchanged is that it asserts that the expression does not evaluate or simplify in any way which is the intention here e.g.:

In [7]: MatPow(Identity(2), 2)                                                                                                    
Out[7]: 
   2
(𝕀) 

In [8]: MatPow(ZeroMatrix(2, 2), 2)                                                                                               
Out[8]: 
   2
(𝟘) 

[oscarbenjamin]

I don't know whether this should raise:

In [9]: MatPow(ZeroMatrix(2, 2), -1)                                                                                              
Out[9]: 
   -1
(𝟘) 

For a non-square matrix it raises before calling doit:

In [4]: MatPow(ZeroMatrix(2, 3), -1)                                                                                              
---------------------------------------------------------------------------
NonSquareMatrixError

[oscarbenjamin]

Thanks @jlherren. This is good.

I still think that this #19104 (comment) needs addressing but I'll open a new issue.