Give more precise answer on symbolic power of a matrix #36845

RuchitJagodara · 2023-12-09T21:44:43Z

Changed _matrix_power_symbolic function to consider mk=0 case
Changed the _matrix_power_symbolic function to handle all the cases
Created tests covering the changes

This PR improves the answer given by _matrix_power_symbolic() function, which involves finding the symbolic power of a matrix, and also gives more precise answer on evalution. eg, before, when we were evaluating the zeroth power of a matrix sometimes it was giving wrong answer like below,

        sage: n = var("n")
        sage: A = matrix(ZZ, [[1,-1], [-1,1]])
        sage: An = A^(n) 
        [ 1/2*2^n -1/2*2^n]
        [-1/2*2^n  1/2*2^n]
        sage: An({n:0})
        [ 1/2 -1/2]
        [-1/2  1/2]

Here, for n=0, it is not giving an identity matrix which is wrong, but after the change it will give answer as below,

        sage: n = var("n")
        sage: A = matrix(ZZ, [[1,-1], [-1,1]])
        sage: An = A^(n) 
        [ 1/2*2^(2*n + 1) + 1/2*kronecker_delta(0, 2*n + 1) -1/2*2^(2*n + 1) + 1/2*kronecker_delta(0, 2*n + 1)]
        [-1/2*2^(2*n + 1) + 1/2*kronecker_delta(0, 2*n + 1)  1/2*2^(2*n + 1) + 1/2*kronecker_delta(0, 2*n + 1)]
        sage: An({n:0})
        [1 0]
        [0 1]

Also, This patch fixes #36838 . The function, _matrix_power_symbolic() was giving symbolic power
of a nilpotent matrix as a null matrix but the correct answer should be represented
in the form of kronecker_delta() function. In the case of mk=0, simplifying only binomial
term should work because afterall the form 0^(n-i) should be simplifed to kronecker_delta() function and no further simplification is needed.
And in every other case it is simplifying whole term so it will handle all other cases.

📝 Checklist

The title is concise, informative, and self-explanatory.
The description explains in detail what this PR is about.
I have linked a relevant issue or discussion.
I have created tests covering the changes.

⌛ Dependencies

None

phul-ste · 2023-12-09T23:57:00Z

src/sage/matrix/matrix2.pyx

+            vk = [(binomial(n, i) * mk._pow_(n-i)).simplify_full()
+                  for i in range(nk)]
+        else:
+            vk = [(binomial(n, i).simplify_full() * mk._pow_(n-i))


In fact, a more suitable result would be given by
vk = [(binomial(n, i).simplify_full() * kronecker_delta(n,i)) for i in range(nk)]
This even allows to evaluate with small values of n

phul-ste · 2023-12-10T00:04:04Z

src/sage/matrix/matrix2.pyx

+        n
+        sage: An = A^n; An
+        [        0^n 0^(n - 1)*n]
+        [          0         0^n]


According to the comment on line 18678,
We would then have something like
[ kronecker_delta(0, n) n*kronecker_delta(1, n)]
[ 0 kronecker_delta(0, n)]
that can then be evaluated
sage: An({n:0})
[1 0]
[0 1]

phul-ste

I checked and it works perfectly to me.

sage: An = _matrix_power_symbolic(A,n); An

[ kronecker_delta(0, n) n*kronecker_delta(1, n)]

[ 0 kronecker_delta(0, n)]

sage: print(all(A^k == An.subs({n: k}) for k in range(8)))
#True

phul-ste · 2023-12-11T21:31:27Z

src/sage/matrix/matrix2.pyx

+        if mk:
+            vk = [(binomial(n, i) * mk._pow_(n-i)).simplify_full()
+                  for i in range(nk)]
+        else:


In fact, I realize that after the computation, we can get a value even for negative n (or even also for non-integer n).

I think that, somehow, the system should raise an error if mk=0 but the evaluation is made at bad values of n. I have no idea on how to do that properly.

This can be done by declaring variable with domain="integer", here when you do not give any domain to a variable it will consider it in complex plane by default.

But we can check whether the variable have an integer value of not at the time of evaluation.

I think that, somehow, the system should raise an error if mk=0 but the evaluation is made at bad values of n.

I didn't get that, why system should give an error if mk=0 ?

I'll implement this in a new PR as this topic is not concentrated on nilpotent matrices.

This can be done by declaring variable with domain="integer", here when you do not give any domain to a variable it will consider it in complex plane by default.

great!

But we can check whether the variable have an integer value of not at the time of evaluation.

I think that, somehow, the system should raise an error if mk=0 but the evaluation is made at bad values of n.

I didn't get that, why system should give an error if mk=0 ?

It should raise an error if mk=0 and n negative since the matrix is then non-invertible. For n positive non-integer, I don't know if it makes sense either.

I'll implement this in a new PR as this topic is not concentrated on nilpotent matrices.

By the way, I already opened an other issue #36863 for general simplification of 0^n.

But the problem for matrices is somehow different since:

the exponent is generally assumed to be an integer

the general formula for the power of a Jordan block with eigenvalue 0 involves some expression of the form P(n)*0^{n-i} (P polynomial) where the meaning of 0^{n-i}for 0<i <=nshould be understood as kronecker_delta(i,n) in this special case (assuming that n is a *non-negative integer)

In fact, I realize that after the computation, we can get a value even for negative n (or even also for non-integer n).
I think that, somehow, the system should raise an error if mk=0 but the evaluation is made at bad values of n. I have no idea on how to do that properly.

I checked what I can do to determine whether we are evaluating an integer power or not. I found that when evaluating the matrix, we only replace the variables with their values, and we cannot impose a condition that the variable should be an integer for obvious reasons because we are not sure that the variable is the symbolic power only. The only place where we can make a change is when calculating the power of the matrix (i.e., A^n). However, if we do so, it will not be ideal because we can define a real power of a matrix, and imposing a condition would not be appropriate. So, I think implementing a system for this is not practical. So it is better to use domain=integer while defining a variable.

RuchitJagodara · 2023-12-13T19:04:51Z

    Testing exponentiation in the integer ring::

        sage: A = matrix(ZZ, [[1,-1], [-1,1]])
        sage: A^(2*n+1)                                                                 # needs sage.symbolic
        [ 1/2*2^(2*n + 1) -1/2*2^(2*n + 1)]
        [-1/2*2^(2*n + 1)  1/2*2^(2*n + 1)]

@dimpase, Can I change the above doctest which was there before my changes, because here, according to my change, I think the answer should be

    [ 1/2*2^(2*n + 1) + 1/2*kronecker_delta(0, 2*n + 1) -1/2*2^(2*n + 1) + 1/2*kronecker_delta(0, 2*n + 1)]
    [-1/2*2^(2*n + 1) + 1/2*kronecker_delta(0, 2*n + 1)  1/2*2^(2*n + 1) + 1/2*kronecker_delta(0, 2*n + 1)]

But the answer given in doctest is different. And to support my answer, I'll say that this is working for 2n+1=0 where we should get an indentity matrix but the answer given in doctest is false as it will not give identity matrix for 2n+1=0. And becuase of this my PR is failing that doctest check.

phul-ste · 2023-12-13T19:43:52Z

@dimpase, Can I change the above doctest which was there before my changes, because here, according to my change, I think the answer should be
    [ 1/2*2^(2*n + 1) + 1/2*kronecker_delta(0, 2*n + 1) -1/2*2^(2*n + 1) + 1/2*kronecker_delta(0, 2*n + 1)]
    [-1/2*2^(2*n + 1) + 1/2*kronecker_delta(0, 2*n + 1)  1/2*2^(2*n + 1) + 1/2*kronecker_delta(0, 2*n + 1)]
But the answer given in doctest is different. And to support my answer, I'll say that this is working for n=0 where we should get an indentity matrix but the answer given in doctest is false as it will not give identity matrix for n=0. And becuase of this my PR is failing that doctest check.

A precision (in case your argument is misunderstood by others), you mean 0 as the power of the matrix (and not for n=0 when we consider the power 2n+1)

RuchitJagodara · 2023-12-13T19:46:58Z

Ohh, ya that's right, thank you for pointing out.

dimpase · 2023-12-14T00:19:11Z

It makes sense to talk about changing the doctest result to match the output of the changed function.

By the way, I don't quite understand the comment

# Return mk^(n-i) instead of mk**(n-i) if mk=0

Do you mean ^ as in Python? That's confusing. Better use a maths way to explain what you mean.

RuchitJagodara · 2023-12-14T08:24:45Z

@dimpase I have updated the results of doctest and also updated the comment. Can you please review this? And I have modified code style for this file according to sage's developer guide which can be found in #36877. After merging that PR, I think Lint check for this PR will be passed.

RuchitJagodara · 2023-12-15T02:57:59Z

In Build & Test / build (pull_request) it is showing that doctests are failing on some files. But when I checked it in my machine was working fine and passing all doctests for all files. Like for below it is showing that 130 doctests are failed but when I ran the same command in my machine it is showing that all doctests are passes.

sage -t --long --random-seed=286735480429121101562228604801325644303 sage/libs/coxeter3/coxeter_group.py  # 130 doctests failed

This fixes sagemath#36838. Here, I have added a condition to check whether mk=0 or not. Because whenever mk=0 we should give mk^(n-i) (i.e. 0^(n-i)) instead of only 0 considering (n-i) can be equal to zero and in this case 0^(n-i) will be more accurate than only 0.

This change handles all the errors which were occuring in previous commit and this commit handles the case of mx=0 very effectively.

Created tests covering the changes and checking whether the trac:`36838` is fixed or not.

Instead of returning 0^(n-i), it would be more precise if we reutrn value of kroncker_delta function, it will be helpful when we try to evaluate the final matrix using the value of n.

Changed the final answer given by test of this PR.

Changed the answer of doctest according to new changes.

Corrected the doctest and improved some code styles, which were not correct according to guidelines for python coding for sage.

Updated the old doctests which were failing before and updated the comment to make it more readable.

dimpase · 2023-12-17T15:06:03Z

Please fix the linting errors

/sage/matrix/matrix2.pyx:3185:30: W291 trailing whitespace
./sage/matrix/matrix2.pyx:3191:29: W291 trailing whitespace
./sage/matrix/matrix2.pyx:3199:44: W291 trailing whitespace
./sage/matrix/matrix2.pyx:18665:106: W291 trailing whitespace
./sage/matrix/matrix2.pyx:18666:106: W291 trailing whitespace

github-actions · 2023-12-17T18:19:34Z

Documentation preview for this PR (built with commit ed4e712; changes) is ready! 🎉

dimpase

looks good, thanks

- Changed _matrix_power_symbolic function to consider mk=0 case - Changed the _matrix_power_symbolic function to handle all the cases - Created tests covering the changes This PR improves the answer given by _matrix_power_symbolic() function, which involves finding the symbolic power of a matrix, and also gives more precise answer on evalution. eg, before, when we were evaluating the zeroth power of a matrix sometimes it was giving wrong answer like below, ```python sage: n = var("n") sage: A = matrix(ZZ, [[1,-1], [-1,1]]) sage: An = A^(n) [ 1/2*2^n -1/2*2^n] [-1/2*2^n 1/2*2^n] sage: An({n:0}) [ 1/2 -1/2] [-1/2 1/2] ``` Here, for n=0, it is not giving an identity matrix which is wrong, but after the change it will give answer as below, ```python sage: n = var("n") sage: A = matrix(ZZ, [[1,-1], [-1,1]]) sage: An = A^(n) [ 1/2*2^(2*n + 1) + 1/2*kronecker_delta(0, 2*n + 1) -1/2*2^(2*n + 1) + 1/2*kronecker_delta(0, 2*n + 1)] [-1/2*2^(2*n + 1) + 1/2*kronecker_delta(0, 2*n + 1) 1/2*2^(2*n + 1) + 1/2*kronecker_delta(0, 2*n + 1)] sage: An({n:0}) [1 0] [0 1] ``` Also, This patch fixes sagemath#36838 . The function, _matrix_power_symbolic() was giving symbolic power of a nilpotent matrix as a null matrix but the correct answer should be represented in the form of kronecker_delta() function. In the case of mk=0, simplifying only binomial term should work because afterall the form 0^(n-i) should be simplifed to kronecker_delta() function and no further simplification is needed. And in every other case it is simplifying whole term so it will handle all other cases.      ### 📝 Checklist     - [x] The title is concise, informative, and self-explanatory. - [x] The description explains in detail what this PR is about. - [x] I have linked a relevant issue or discussion. - [x] I have created tests covering the changes. ### ⌛ Dependencies   None URL: sagemath#36845 Reported by: Ruchit Jagodara Reviewer(s): Dima Pasechnik, phul-ste, Ruchit Jagodara

RuchitJagodara changed the title ~~Fix wrong answer on symbolic power of nilpotent matrix~~ Fix answer on symbolic power of nilpotent matrix Dec 9, 2023

RuchitJagodara changed the title ~~Fix answer on symbolic power of nilpotent matrix~~ Give more precise answer on symbolic power of nilpotent matrix Dec 9, 2023

phul-ste reviewed Dec 9, 2023

View reviewed changes

phul-ste reviewed Dec 10, 2023

View reviewed changes

RuchitJagodara force-pushed the power_of_nilpotent_matrix branch from 4757770 to 3042d2a Compare December 10, 2023 04:24

RuchitJagodara requested a review from phul-ste December 10, 2023 04:27

phul-ste approved these changes Dec 10, 2023

View reviewed changes

phul-ste reviewed Dec 11, 2023

View reviewed changes

dimpase force-pushed the power_of_nilpotent_matrix branch from ab1ec91 to fc7526d Compare December 13, 2023 15:42

RuchitJagodara changed the title ~~Give more precise answer on symbolic power of nilpotent matrix~~ Give more precise answer on symbolic power of a matrix Dec 13, 2023

RuchitJagodara force-pushed the power_of_nilpotent_matrix branch from fc7526d to 700f425 Compare December 14, 2023 08:22

RuchitJagodara added 8 commits December 17, 2023 14:25

Changed the _matrix_power_symbolic function to handle all the cases

dc79089

This change handles all the errors which were occuring in previous commit and this commit handles the case of mx=0 very effectively.

Created tests covering the changes

7139bb8

Created tests covering the changes and checking whether the trac:`36838` is fixed or not.

Give more precise answer by using kroncker_delta function

42454c0

Instead of returning 0^(n-i), it would be more precise if we reutrn value of kroncker_delta function, it will be helpful when we try to evaluate the final matrix using the value of n.

Rewriting tests for the PR

96442fb

Changed the final answer given by test of this PR.

Correct answers of the doctest for _matrix_power_symbolic function

d40d348

Changed the answer of doctest according to new changes.

Corrected the doctest

935a420

Corrected the doctest and improved some code styles, which were not correct according to guidelines for python coding for sage.

Modified the doctest and changed the comment

7955a51

Updated the old doctests which were failing before and updated the comment to make it more readable.

dimpase force-pushed the power_of_nilpotent_matrix branch from 700f425 to 7955a51 Compare December 17, 2023 14:25

Corrected lint errors

ed4e712

dimpase approved these changes Dec 18, 2023

View reviewed changes

dimpase added c: symbolics s: positive review labels Dec 18, 2023

Merge branch 'sagemath:develop' into power_of_nilpotent_matrix

f95793b

vbraun merged commit 58e899c into sagemath:develop Dec 26, 2023
1 check passed

github-actions bot removed the s: positive review label Dec 26, 2023

mkoeppe added this to the sage-10.3 milestone Dec 26, 2023

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Give more precise answer on symbolic power of a matrix #36845

Give more precise answer on symbolic power of a matrix #36845

RuchitJagodara commented Dec 9, 2023 •

edited

phul-ste Dec 9, 2023

phul-ste Dec 10, 2023 •

edited

phul-ste left a comment

phul-ste Dec 11, 2023

RuchitJagodara Dec 12, 2023

RuchitJagodara Dec 12, 2023

phul-ste Dec 12, 2023

phul-ste Dec 12, 2023 •

edited

RuchitJagodara Dec 13, 2023

RuchitJagodara commented Dec 13, 2023 •

edited

phul-ste commented Dec 13, 2023

RuchitJagodara commented Dec 13, 2023

dimpase commented Dec 14, 2023

RuchitJagodara commented Dec 14, 2023 •

edited

RuchitJagodara commented Dec 15, 2023

dimpase commented Dec 17, 2023

github-actions bot commented Dec 17, 2023

dimpase left a comment

Give more precise answer on symbolic power of a matrix #36845

Give more precise answer on symbolic power of a matrix #36845

Conversation

RuchitJagodara commented Dec 9, 2023 • edited

📝 Checklist

⌛ Dependencies

phul-ste Dec 9, 2023

Choose a reason for hiding this comment

phul-ste Dec 10, 2023 • edited

Choose a reason for hiding this comment

phul-ste left a comment

Choose a reason for hiding this comment

[ kronecker_delta(0, n) n*kronecker_delta(1, n)]

[ 0 kronecker_delta(0, n)]

phul-ste Dec 11, 2023

Choose a reason for hiding this comment

RuchitJagodara Dec 12, 2023

Choose a reason for hiding this comment

RuchitJagodara Dec 12, 2023

Choose a reason for hiding this comment

phul-ste Dec 12, 2023

Choose a reason for hiding this comment

phul-ste Dec 12, 2023 • edited

Choose a reason for hiding this comment

RuchitJagodara Dec 13, 2023

Choose a reason for hiding this comment

RuchitJagodara commented Dec 13, 2023 • edited

phul-ste commented Dec 13, 2023

RuchitJagodara commented Dec 13, 2023

dimpase commented Dec 14, 2023

RuchitJagodara commented Dec 14, 2023 • edited

RuchitJagodara commented Dec 15, 2023

dimpase commented Dec 17, 2023

github-actions bot commented Dec 17, 2023

dimpase left a comment

Choose a reason for hiding this comment

RuchitJagodara commented Dec 9, 2023 •

edited

phul-ste Dec 10, 2023 •

edited

phul-ste Dec 12, 2023 •

edited

RuchitJagodara commented Dec 13, 2023 •

edited

RuchitJagodara commented Dec 14, 2023 •

edited