matrices: fixed MatPow.doit() so now == mat**x

[tom-pytel]

MatPow.doit() did not use _matrix_pow_by_jordan_blocks() where it could
and sometimes did not arrive at the same result as mat**x for matrices
where this was possible. Added the appropriate check an usage.

Example:

Matrix ([[1,0],[0,1]])**x
Matrix([
[1, 0],
[0, 1]])

MatPow (Matrix ([[1,0],[0,1]]), x).doit ()
Matrix([
[1, 0],
[0, 1]])**x

Now:

MatPow (Matrix ([[1,0],[0,1]]), x).doit ()
Matrix([
[1, 0],
[0, 1]])

References to other Issues or PRs

Fixes 17175

Brief description of what is fixed or changed

Other comments

Release Notes

• matrices
  • Added usage of _matrix_pow_by_jordan_blocks() in MatPow.doit() where aplicable.

[sympy-bot]

✅

Hi, I am the SymPy bot (v147). 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
  • Added usage of _matrix_pow_by_jordan_blocks() in MatPow.doit() where aplicable. (#17179 by @Pristine-Cat)

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

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.

MatPow.doit() did not use _matrix_pow_by_jordan_blocks() where it could
and sometimes did not arrive at the same result as mat**x for matrices
where this was possible. Added the appropriate check an usage.

Example:
```python
Matrix ([[1,0],[0,1]])**x
Matrix([
[1, 0],
[0, 1]])

MatPow (Matrix ([[1,0],[0,1]]), x).doit ()
Matrix([
[1, 0],
[0, 1]])**x
```
Now:
```python
MatPow (Matrix ([[1,0],[0,1]]), x).doit ()
Matrix([
[1, 0],
[0, 1]])
```
#### References to other Issues or PRs
Fixes 17175


#### Brief description of what is fixed or changed


#### Other comments


#### Release Notes

<!-- BEGIN RELEASE NOTES -->
* matrices
  * Added usage of `_matrix_pow_by_jordan_blocks()` in `MatPow.doit()` where aplicable.
<!-- END RELEASE NOTES -->

Update

The release notes on the wiki have been updated.

[codecov]

Codecov Report

Merging #17179 into master will decrease coverage by 0.011%.
The diff coverage is 100%.

@@              Coverage Diff              @@
##            master    #17179       +/-   ##
=============================================
- Coverage   74.529%   74.517%   -0.012%     
=============================================
  Files          623       623               
  Lines       161546    161545        -1     
  Branches     37910     37910               
=============================================
- Hits        120400    120380       -20     
- Misses       35811     35830       +19     
  Partials      5335      5335

[oscarbenjamin]

Is it valid to use the Jordan normal form like this for non-integer exponents?

[tom-pytel]

I can not say but the author of MatrixArithmetic.__pow__() could and that is where I copied the usage from under the same circumstances. All I did was make MatPow.doit() CONSISTENT with MatrixArithmetic.__pow__() and if the latter were broken I am guessing it should have been noticed by now? It does actually seem to work empirically testing matrices to symbol powers later substituted with actual values, even fractional, but I can not give you a mathematical proof of the correctness of this.

[asmeurer]

Can you add a test for this?

[tom-pytel]

Test for what specifically? The consistency of mat**x vs. MatPow (mat, x).doit()?

[tom-pytel]

Actually let me potentially reverse myself, it seems the condition in MatrixArithmetic.pow() may actually be incorrect as:

Matrix ([[0,1,0],[0,0,1],[0,0,1]])**Symbol('n', nonnegative=True)

raises a

ValueError: Matrix det == 0; not invertible

when it should return a MatPow object instead of trying to evaluate a zero determinant matrix?

[oscarbenjamin]

I can not say but the author of MatrixArithmetic.pow() could

Don't place too much faith in the existing code. You found an inconsistency which suggests a bug somewhere. Now it all needs closer inspection.

I'm not sure how to convince myself that something like this is correct:

In [2]: Matrix([[1, 1], [0, 1]]) ** I                                                                                             
Out[2]: 
⎡1  ⅈ⎤
⎢    ⎥
⎣0  1⎦

This one raises:

In [3]: Matrix([[0, 1], [0, 1]]) ** I                                                                                             
---------------------------------------------------------------------------
...
TypeError: Invalid comparison of complex I

which suggests that the author hadn't considered non-real exponents.

[tom-pytel]

|In [2]: Matrix([[1, 1], [0, 1]]) ** I |

Seems to be treating I as a generic Symbol. >> Matrix ([[1, 1], [0, 1]])**x Matrix([ [1, x], [0, 1]]) Which looks like it gives the correct answer for real x including negative exponents (inverses of positive powers). As well as dealing with more complicated cases: (Matrix([[1, 2], [3, 4]])**x).subs(x, 1.5) - Matrix([[1, 2], [3, 4]])**1.5 Matrix([ [0, 0], [0, 0]]) Also I could see an imaginary power of a matrix being valid because of: If a matrix is diagonalizable, then diagonalize it,A=PDP−1A=PDP−1and apply the power to the diagonal An=PDnP−1An=PDnP−1 The diagonal values are acted on individually. from: https://math.stackexchange.com/questions/732511/fractional-power-of-matrix

[oscarbenjamin]

If a matrix is diagonalizable

That's why I deliberately gave a non-diagonalisable example :)

The question is: what is the first principles definition that we are working to here?

[tom-pytel]

The question is: what is the first principles definition that we are working to here?

I see your point, and if you are asking if it even makes sense to raise a matrix to a complex power then yes its weird. Looking through the code for _matrix_pow_by_jordan_blocks() in matrices seems to imply that non-integer powers were never intended, there is a raise ValueError("Non-integer power cannot be evaluated") which is never reached. The fact that it works for some is probably coincidence.

[sylee957]

I think that imaginary numbers are subclasses of Expr, so the original check may not work.
The check in __pow__ may have to be changed to

        elif num.is_real:
            return jordan_pow(num)

[oscarbenjamin]

@sylee957 do you think that the results for imaginary exponents are incorrect? More generally what is supposed to happen for non-integer exponents?

Let's take M**(S(1)/2) which corresponds to the matrix square root. An NxN matrix has up to 2**N distinct square roots. If the matrix is positive definite then there is a unique positive definite square root. What should we do if the matrix is not positive definite? Do we use the Jordan form approach and take the principal root of each eigenvalue?

How do we generalise that to complex exponents? What about irrational exponents like M**sqrt(2) - should that be exp(sqrt(2)*log(M))?

[jksuom]

Would it be possible to try base**exp to see if it raises an exception? If it does, then the expression could be returned unmodified.

[oscarbenjamin]

Would it be possible to try base**exp to see if it raises an exception? If it does, then the expression could be returned unmodified.

It would but I feel like that using that now could be masking bugs.

I think what is needed is a definition of what should be the correct return value in SymPy and when the result should be able to evaluate. Here's an attempt:

1. For M**n where n is a non-negative integer the result is always well defined. It can be computed from the Jordan form including for symbolic n (that satisfies the assumptions). The trivial case n=0 always gives the identity.

2. If n is a negative integer then the matrix needs to be invertible but otherwise the result is as in case 1.

3. If n is an integer, not known to be positive or negative then it can evaluate if the matrix is invertible but should remain unevaluated otherwise. Substitution of a negative integer for n should raise. Substitution of a nonnegative integer should evaluate.

4. If n is rational p/q then this is the pth power of the qth root of M (or the qth root of the pth power?). As for integers negative n is only allowed for invertible M. The qth root is not unique and does not always exist e.g. [[0, 1], [0, 0]] has no square root so this can raise. Where the roots exist there can be up to q**N possibilities if M is NxN.

4a. If M is positive definite then I think it has a unique positive definite qth root so we use that.

4b. If M has negative or non-real eigenvalues then provided it is not defective we can just take the pth root of each eigenvalue (taking SymPy's definition of the pth root for scalars).

4c. If M is defective then the Jordan form method can still work using pth roots of the eigenvalues.

5. If n is irrational or non-real complex then M**n = exp(n*log(M)) - log(M) exists iff M is invertible and is unique if M is positive definite. If M is not positive definite then different choices to log(M) will not give equivalent results so raise? This can be taken as the definition of correctness although it might still be easiest to compute via the Jordan form method.

6. All of the above assumes that the entries of M are complex-valued. If M has matrices or operators or something as entries then only integer powers can be defined? The result could depend on whether the elements are invertible...

I guess a minimum requirement here is that all of the above needs to match SymPy scalars in the case of a 1x1 matrix.

Have I missed any important cases above?

[tom-pytel]

1. For |M**n| where n is a non-negative integer the result is always well defined. It can be computed from the Jordan form including for symbolic n (that satisfies the assumptions). The trivial case n=0 always gives the identity.

Consider disallowing symbolic n. It works because it winds up encoding the steps from _matrix_pow_by_jordan_blocks() within the matrix itself which later evaluates correctly only if you substitute in something which fits the assumptions (as you noted). But this can get ugly and slow to the point of never terminating, have a look at the output of Matrix([[a,b],[c,d]])**n (all Symbols). A 3x3 done like this never finishes (or maybe it does and I just didn't wait long enough).

[oscarbenjamin]

I don't think we should disallow symbolic n. Lots of things can lead to symbolic explosion, especially if you allow everything to be an arbitrary symbol (try solving an arbitrary quartic) but they can still produce intelligible results in more constrained cases e.g.:

In [5]: Matrix([[a, 0], [1, 1]])**n                                                                                               
Out[5]: 
⎡      n         ⎤
⎢     a         0⎥
⎢                ⎥
⎢   n            ⎥
⎢  a       1     ⎥
⎢───── - ─────  1⎥
⎣a - 1   a - 1   ⎦

That can be a reason for not evaluating by default though and requiring a doit from users.

Also it's always the case that you shouldn't substitute a value that violates the assumptions on a symbol. I don't know whether subs should check for that sort of thing.

[tom-pytel]

Warning, it is evident that the original MatrixArithmetic.pow() still has problems, all this change did was make MatPow.doit() consistent, but those problems will still pop up in some cases. It passes all the tests and is not common for the problems to pop up though.

[oscarbenjamin]

@smichr I didn't think this PR was finished yet. It has no tests for example.

We need a new issue to add tests for anything fixed here and to track the unresolved discussion around which cases should be defined...

[tom-pytel]

Ok, I am working to patch it up but I keep running into the problem that many of the tests are coded for the way this works now, expecting a specific exception or a symbolic representation in their own specific cases. Some tests seem to expect non-integers to work while others do not for example. Here is an example of two tests which seemingly expect a different result from the same evaluation (and get it in the current state of things):

A = Matrix([[0, 1, 0], [0, 0, 1], [0, 0, 1]])  # Nilpotent jordan block size 2
n = Symbol('n', integer=True)
assert isinstance(A**n, MatPow)

vs.

from sympy.abc import a, b, n
assert Matrix([[1, a], [0, 1]])**n == Matrix([[1, a*n], [0, 1]])

How to deal with such situation where cleaning up old code breaks tests?

Another example:
assert Matrix([[1, 0], [1, 1]])^0.5 == Matrix([[1.0, 0], [0.5, 1.0]])
vs.
raises(ValueError, lambda: A^2.1)

Should I test floats for rationality? And with what cutoff?

One more, this test does not even seem to make sense:

n = Symbol('n', integer=True, nonnegative=True)
raises(ValueError, lambda: A**n)

Any matrix to a nonnegative integer is always possible no? But the test expects it to fail (which it does currently).

[oscarbenjamin]

Here is an example of two tests

The first example is a non-invertible matrix so that's my case 3. above. The second example is invertible so it's different.

How to deal with such situation where cleaning up old code breaks tests?

First we need to define consistent rules for what should happen. Then if necessary the tests can be changed.

Let's take this discussion back to the original issue or a new PR since this PR is closed now.