Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

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

Merged
merged 10 commits into from
Dec 26, 2023
Merged

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

Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
afa6a7c
Changed _matrix_power_symbolic function to condier mk=0 case
RuchitJagodara Dec 9, 2023
dc79089
Changed the _matrix_power_symbolic function to handle all the cases
RuchitJagodara Dec 9, 2023
7139bb8
Created tests covering the changes
RuchitJagodara Dec 9, 2023
42454c0
Give more precise answer by using kroncker_delta function
RuchitJagodara Dec 10, 2023
96442fb
Rewriting tests for the PR
RuchitJagodara Dec 10, 2023
d40d348
Correct answers of the doctest for _matrix_power_symbolic function
RuchitJagodara Dec 13, 2023
935a420
Corrected the doctest
RuchitJagodara Dec 13, 2023
7955a51
Modified the doctest and changed the comment
RuchitJagodara Dec 14, 2023
ed4e712
Corrected lint errors
RuchitJagodara Dec 17, 2023
f95793b
Merge branch 'sagemath:develop' into power_of_nilpotent_matrix
RuchitJagodara Dec 19, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
47 changes: 30 additions & 17 deletions src/sage/matrix/matrix2.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -216,10 +216,10 @@ cdef class Matrix(Matrix1):
from sage.matrix.constructor import matrix
if self.is_sparse():
return matrix({ij: self[ij].subs(*args, **kwds) for ij in self.nonzero_positions()},
nrows=self._nrows, ncols=self._ncols, sparse=True)
nrows=self._nrows, ncols=self._ncols, sparse=True)
else:
return matrix([a.subs(*args, **kwds) for a in self.list()],
nrows=self._nrows, ncols=self._ncols, sparse=False)
nrows=self._nrows, ncols=self._ncols, sparse=False)

def solve_left(self, B, check=True):
"""
Expand Down Expand Up @@ -3183,14 +3183,12 @@ cdef class Matrix(Matrix1):
"""

# Validate assertions
#
if not self.is_square():
raise ValueError("self must be a square matrix")

from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing

# Extract parameters
#
cdef Matrix M = <Matrix> self
n = M._ncols
R = M._base_ring
Expand All @@ -3199,7 +3197,6 @@ cdef class Matrix(Matrix1):
# Corner cases
# N.B. We already tested for M to be square, hence we do not need to
# test for 0 x n or m x 0 matrices.
#
if n == 0:
return S.one()

Expand All @@ -3215,7 +3212,6 @@ cdef class Matrix(Matrix1):
#
# N.B. The documentation is still 1-based, although the code, after
# having been ported from Magma to Sage, is 0-based.
#
from sage.matrix.constructor import matrix

F = [R.zero()] * n
Expand All @@ -3226,30 +3222,25 @@ cdef class Matrix(Matrix1):
for t in range(1,n):

# Set a(1, t) to be M(<=t, t)
#
for i in range(t+1):
a.set_unsafe(0, i, M.get_unsafe(i, t))

# Set A[1, t] to be the (t)th entry in a[1, t]
#
A[0] = M.get_unsafe(t, t)

for p in range(1, t):

# Set a(p, t) to the product of M[<=t, <=t] * a(p-1, t)
#
for i in range(t+1):
s = R.zero()
for j in range(t+1):
s = s + M.get_unsafe(i, j) * a.get_unsafe(p-1, j)
a.set_unsafe(p, i, s)

# Set A[p, t] to be the (t)th entry in a[p, t]
#
A[p] = a.get_unsafe(p, t)

# Set A[t, t] to be M[t, <=t] * a(p-1, t)
#
s = R.zero()
for j in range(t+1):
s = s + M.get_unsafe(t, j) * a.get_unsafe(t-1, j)
Expand All @@ -3262,7 +3253,7 @@ cdef class Matrix(Matrix1):
F[p] = s - A[p]

X = S.gen(0)
f = X ** n + S(list(reversed(F)))
f = X**n + S(list(reversed(F)))

return f

Expand Down Expand Up @@ -3634,7 +3625,7 @@ cdef class Matrix(Matrix1):
# using the rows of a matrix.) Also see Cohen's first GTM,
# Algorithm 2.2.9.

cdef Py_ssize_t i, m, n,
cdef Py_ssize_t i, m, n
n = self._nrows

cdef Matrix c
Expand Down Expand Up @@ -18607,8 +18598,8 @@ def _matrix_power_symbolic(A, n):

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)]
[ 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)]

Check if :trac:`23215` is fixed::

Expand All @@ -18618,12 +18609,25 @@ def _matrix_power_symbolic(A, n):
-1/2*I*(a + I*b)^k + 1/2*I*(a - I*b)^k,
1/2*I*(a + I*b)^k - 1/2*I*(a - I*b)^k,
1/2*(a + I*b)^k + 1/2*(a - I*b)^k]

Check if :trac:`36838` is fixed:Checking symbolic power of
nilpotent matrix::

sage: A = matrix([[0,1],[0,0]]); A
[0 1]
[0 0]
sage: n = var('n'); n
n
sage: B = A^n; B
[ kronecker_delta(0, n) n*kronecker_delta(1, n)]
[ 0 kronecker_delta(0, n)]
"""
from sage.rings.qqbar import AlgebraicNumber
from sage.matrix.constructor import matrix
from sage.functions.other import binomial
from sage.symbolic.ring import SR
from sage.rings.qqbar import QQbar
from sage.functions.generalized import kronecker_delta

got_SR = A.base_ring() == SR

Expand Down Expand Up @@ -18656,8 +18660,17 @@ def _matrix_power_symbolic(A, n):
# D^i(f) / i! with f = x^n and D = differentiation wrt x
if hasattr(mk, 'radical_expression'):
mk = mk.radical_expression()
vk = [(binomial(n, i) * mk**(n-i)).simplify_full()
for i in range(nk)]


# When the variable "mk" is equal to zero, it is advisable to employ the Kronecker delta function
# instead of utilizing the numerical value zero. This choice is made to encompass scenarios where
# the power of zero is also equal to zero.
if mk:
vk = [(binomial(n, i) * mk._pow_(n-i)).simplify_full()
for i in range(nk)]
else:
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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 ?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link

@phul-ste phul-ste Dec 12, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

vk = [(binomial(n, i).simplify_full() * kronecker_delta(n,i))
for i in range(nk)]

# Form block Mk and insert it in M
Mk = matrix(SR, [[SR.zero()]*i + vk[:nk-i] for i in range(nk)])
Expand Down