Errors in determinant of 'large' symbolic matrices as opposed to working over better rings

[Linden Disney]

Attached is a jupyter notebook that runs Sage 9.1, a (slightly more) minimal example of a problem that I discovered. When calculating the determinant of a large (in the sense n>=9 I have currently found) symbolic matrix the answer is not correct. To see this, run the notebook with Qsimplify either True or False. When Qsimplify is false, the calculation is done when variables lie in the symbolic ring, when true a specially constructed ring is used instead. The output of the script shows a matrix and the resulting characteristic polynomial after some simplification has occured.  While the two matrices look the same regardless of Qsimplify, the characteristic polynomial changes. This error goes away for smaller matrices (it first turns up at rank=4, where the rank is of the Lie algebra involved in the calculation, but this just gives the basis the matrix is constructed from). General theory tells us that the answer when Qsimplify is true is the correct one.  

[Linden Disney]

As a follow up on this, it seems that sage implements the determinant by evaluating the characteristic polynomials at 0, and the characteristic polynomial is calculated by maxima. Is it possible to edit the maxima source code in sage? 

[Michael Orlitzky]

On 10/12/20 8:16 AM, Linden Disney wrote:
> Attached is a jupyter notebook that runs Sage 9.1, a (slightly more)
> minimal example of a problem that I discovered. When calculating the
> determinant of a large (in the sense n>=9 I have currently found)

> symbolic matrix the answer is not correct....

I've tried to translate this into a plain sage script (for the
command-line) to get a minimal example. Really, all we need is the
matrices, but I ran into some trouble with indices in the Qsimplify=True
branch. That could mean I don't know enough about Jupyter, but it also
might just mean that there's a bug in that part of the code.

Is this still a problem for you? If so, can you try to reproduce only
the two (equivalent) matrices whose determinants are different? The
sage_input() function might be of some help here. Ideally the example
would look like

<some variable definitions>
A = matrix(...)
B = matrix(...)
A.det() == B.det()

[Linden Disney]

Ok I've modified the code to plain sage to make it more useful and I've copied it below. Given that it's hard to compare the determinants of the raw matrices, as they are defined in terms of different variables, I have found the z^2 coefficient in each case and you can see they are different. I've tested this in Sage 9.1 command line and this has worked. 

rank = 4

N = 9

p = [var("p{}".format(i), latex_name="p_{}".format(i)) for i in (1..rank)]

w = var('w')

z = var('z')

Q = [var("Q{}".format(i), latex_name="Q_{}".format(i)) for i in (0..rank)]

ni = [1,1,2,2,2]

Qn = [Q[i]^ni[i] for i in (0..rank)]

L1 = Matrix([[-p1, Q1, 0, 0, 0, -4*Q0/z, 0, 0, 0],

            [Q1, p1-p2, Q2, 0, 4*Q0/z, 0, 0, 0, 0],

            [0, Q2, p2-p3, Q3, 0, 0, 0, 0, 0],

            [0, 0, Q3, p3-2*p4, 0, 0, 0, 0, 2*Q4],

            [0, Q0*z, 0, 0, p1, -Q1, 0, 0, 0],

            [-Q0*z, 0, 0, 0, -Q1, -p1+p2, -Q2, 0, 0],

            [0, 0, 0, 0, 0, -Q2, -p2+p3, -Q3, 0],

            [0, 0, 0, 0, 0, 0, -Q3, -p3+2*p4, -Q4],

            [0, 0, 0, Q4, 0, 0, 0, -2*Q4, 0]])

C1 = L1 - w*matrix.identity(N)

C1 = C1.det()

C1 = C1.subs({Q0:prod(ni)/prod(Qn[1:])}).simplify_full()

D1 = C1*z

print(D1.collect(z).coefficients(z)[2])

variables = ["p{}".format(i) for i in (1..rank)] + ["Q{}".format(i) for i in (0..rank)]

R = PolynomialRing(QQ,variables)

R.inject_variables()

Q = R.gens()[rank:2*rank+1]

Qn = [Q[i]^ni[i] for i in (0..rank)]

I = R.ideal(prod(Qn)-prod(ni))

Rq = QuotientRing(R,I,variables)

p = Rq.gens()[0:rank]

Q = Rq.gens()[rank:2*rank+1]

LR = LaurentPolynomialRing(Rq,['w','z'])

LR.inject_variables()

L2 = Matrix([[-p1, Q1, 0, 0, 0, -4*Q0/z, 0, 0, 0],

            [Q1, p1-p2, Q2, 0, 4*Q0/z, 0, 0, 0, 0],

            [0, Q2, p2-p3, Q3, 0, 0, 0, 0, 0],

            [0, 0, Q3, p3-2*p4, 0, 0, 0, 0, 2*Q4],

            [0, Q0*z, 0, 0, p1, -Q1, 0, 0, 0],

            [-Q0*z, 0, 0, 0, -Q1, -p1+p2, -Q2, 0, 0],

            [0, 0, 0, 0, 0, -Q2, -p2+p3, -Q3, 0],

            [0, 0, 0, 0, 0, 0, -Q3, -p3+2*p4, -Q4],

            [0, 0, 0, Q4, 0, 0, 0, -2*Q4, 0]])

C2 = L2 - w*matrix.identity(N)

C2 = C2.det()

D2 = C2*z

display(D2.exponents(), D2.coefficients()[4])

[Michael Orlitzky]

On 12/16/20 3:27 PM, Linden Disney wrote:
> Ok I've modified the code to plain sage to make it more useful and I've
> copied it below. Given that it's hard to compare the determinants of the
> raw matrices, as they are defined in terms of different variables, I
> have found the z^2 coefficient in each case and you can see they are
> different. I've tested this in Sage 9.1 command line and this has worked.
>

Great, thanks. I was able to simplify this even further... the
determinants still disagree at the L1/L2 stage, but look a lot nicer.
Maybe now someone can figure out what's going wrong. With these
definitions, you can even try (L1 == L2) and it returns True.

p1, p2, p3, p4, Q0, Q1, Q2, Q3, Q4, w, z = SR.var("p1, p2, p3, p4, Q0,
Q1, Q2, Q3, Q4, w, z")

L1 = Matrix([[-p1, Q1, 0, 0, 0, -4*Q0/z, 0, 0, 0],
[Q1, p1-p2, Q2, 0, 4*Q0/z, 0, 0, 0, 0],
[0, Q2, p2-p3, Q3, 0, 0, 0, 0, 0],
[0, 0, Q3, p3-2*p4, 0, 0, 0, 0, 2*Q4],
[0, Q0*z, 0, 0, p1, -Q1, 0, 0, 0],
[-Q0*z, 0, 0, 0, -Q1, -p1+p2, -Q2, 0, 0],
[0, 0, 0, 0, 0, -Q2, -p2+p3, -Q3, 0],
[0, 0, 0, 0, 0, 0, -Q3, -p3+2*p4, -Q4],
[0, 0, 0, Q4, 0, 0, 0, -2*Q4, 0]])

L1.det()


R = LaurentPolynomialRing(QQ, "p1, p2, p3, p4, Q0, Q1, Q2, Q3, Q4, w, z")
p1, p2, p3, p4, Q0, Q1, Q2, Q3, Q4, w, z = R.gens()

L2 = Matrix([[-p1, Q1, 0, 0, 0, -4*Q0/z, 0, 0, 0],
[Q1, p1-p2, Q2, 0, 4*Q0/z, 0, 0, 0, 0],
[0, Q2, p2-p3, Q3, 0, 0, 0, 0, 0],
[0, 0, Q3, p3-2*p4, 0, 0, 0, 0, 2*Q4],
[0, Q0*z, 0, 0, p1, -Q1, 0, 0, 0],
[-Q0*z, 0, 0, 0, -Q1, -p1+p2, -Q2, 0, 0],
[0, 0, 0, 0, 0, -Q2, -p2+p3, -Q3, 0],
[0, 0, 0, 0, 0, 0, -Q3, -p3+2*p4, -Q4],
[0, 0, 0, Q4, 0, 0, 0, -2*Q4, 0]])

L2.det()

[dmo...@deductivepress.ca]

That does indeed seem simple. Here is an even shorter version that only needs one matrix.

R = LaurentPolynomialRing(QQ, "p1, p2, p3, p4, Q0, Q1, Q2, Q3, Q4, w, z")
p1, p2, p3, p4, Q0, Q1, Q2, Q3, Q4, w, z = R.gens()

L = Matrix([[-p1, Q1, 0, 0, 0, -4*Q0/z, 0, 0, 0],

           [Q1, p1-p2, Q2, 0, 4*Q0/z, 0, 0, 0, 0],
           [0, Q2, p2-p3, Q3, 0, 0, 0, 0, 0],
           [0, 0, Q3, p3-2*p4, 0, 0, 0, 0, 2*Q4],
           [0, Q0*z, 0, 0, p1, -Q1, 0, 0, 0],
           [-Q0*z, 0, 0, 0, -Q1, -p1+p2, -Q2, 0, 0],
           [0, 0, 0, 0, 0, -Q2, -p2+p3, -Q3, 0],
           [0, 0, 0, 0, 0, 0, -Q3, -p3+2*p4, -Q4],
           [0, 0, 0, Q4, 0, 0, 0, -2*Q4, 0]])

print(f"Laurent: {L.det()}")
print(f"SR: {L.change_ring(SR).det().factor()}")

For me, the output is:

Laurent: 0

SR: 4*Q0*Q1*Q2^2*Q3^2*Q4^2*p1*(z + 1)*(z - 1)/z

I thought maybe the problem had something to do with L being singular, but adding the identity matrix to L doesn't solve the problem:

I = matrix.identity(9)
((L.change_ring(SR) + 42*I).det() - SR((L + 42*I).det())).factor()

Both individual determinants are a mess, but the difference I get is

4*Q0*Q1*Q2^2*Q3^2*Q4^2*(p1 - 42)*(z + 1)*(z - 1)/z

It is surprising to me that this is the same as the above determinant in SR, except for the 42.  And I get a similar answer if I replace 42 with a different constant, or even put in a variable by defining a polynomial ring over R.

I think a ticket should be opened.

[Linden Disney]

Currently making sage_9.3 on machine so I can have open a new ticket and have a deeper look, but I'm quite convinced it is somehow a problem due to the complexity of the calculation in the symbolic ring. 

Firstly note that the answer det(L)=0 is the correct one (you can either trust this from theory or read "A Lie Theoretic Galois Theory for Spectral Curves of an Integrable System I" by McDaniel, Smolinsky where this is considered, albeit in different notation). 

If we try 

M = Matrix([[-p1, Q1, 0, 0, -4*Q0/z, 0, 0],

           [Q1, p1-p2, Q2, 4*Q0/z, 0, 0, 0],

           [0, Q2, p2-p3, 0, 0, 0, 2*Q3],

           [0, Q0*z, 0, p1, -Q1, 0, 0],

           [-Q0*z, 0, 0, -Q1, -p1+p2, -Q2, 0],

           [0, 0, 0, 0, -Q2, -p2+2*p3, -Q3],

           [0, 0, Q3, 0, 0, -2*Q3, 0]])

bool(M.change_ring(SR).det()==M.det())

which is in some way the natural predecessor to this matrix L, we do not have a problem, both dets being equal. 

Alternatively, consider making substitutions for the parameters which shouldn't really the outcome e.g. 

print(f"Laurent: {L.subs({p2:0}).det()}")

print(f"SR: {L.subs({p2:0}).change_ring(SR).det().factor()}")

or 

print(f"Laurent: {L.subs({Q0:1,Q1:1}).det()}")

print(f"SR: {L.subs({Q0:1,Q1:1}).change_ring(SR).det().factor()}")

[dmo...@deductivepress.ca]

I created trac ticket #31077 to continue this discussion.  But first I simplified the example a bit more by making some substitutions to reduce the number of variables, and changing the names to a,b,c,d,e,z. Also, I eliminated the need to compare with the Laurent polynomial ring, by noticing that the determinant of z*L is 0, which is inconsistent with L having nonzero determinant.