Eigenvalues over symbolic ring

asked 2019-02-28 01:50:15 +0200

angelolucia
23 ●1 ●6

updated 2019-02-28 20:35:34 +0200

I am trying to compute the eigenvalues of a 16x16 matrix whose eigenvalues are multivariate polynomials in 4 variables of degree at most 3 and integer coefficients. When I try to do so using the Symbolic Ring SR my machine quickly runs out of memory (~15GB available I believe).

TypeError: ECL says: Memory limit reached. Please jump to an outer pointer, quit program and enlarge the memory limits before executing the program again.

I am using sage 8.1 (and I could try to upgrade if that is the issue here). The example I am struggling with is included at the end of this question (I did not write this by hand, but I don't think that it is relevant how the matrix was generated. One should be able to copy-and-paste this into Sage)

I tried to change the ring to the fraction field of a Multivariate Polynomial Ring over Rational Field, and this allows to quickly compute the characteristic polynomial and its roots, without killing my machine. Unfortunately though some of the eigenvalues are not rational functions, so they do not show up. This is what I did

R = QQ['a,b,c,d']
N = matrix(R, 16, 16, M)
N.change_ring(R.fraction_field()).characteristic_polynomial().roots()

Mathematica on the same machine is able to compute the eigenvalues symbolically (quite quickly actually), and I can solve a similar problem in Sage if I only use 2 variables. Am I missing something, or is simply the algorithm used by Sage not efficient enough for 16x16 matrices in 4 variables?

a,b,c,d = var('a','b','c','d', domain='positive')
M = matrix(16,16,[
(3*c^2*d + 6*c + 1, 0, 0, (a - b)*c*d + a - b, 0, (a - b)*c*d + a - b, (a - b)*c*d + a - b, 0, 0, (a - b)*c*d + a - b, (a - b)*c*d + a - b, 0, (a - b)*c*d + a - b, 0, 0, 3*(a - b)^2*d),
(0, -3*c^2*d + 1, (a + b)*c*d + a + b, 0, (a + b)*c*d + a + b, 0, 0, -(a - b)*c*d + a - b, (a + b)*c*d + a + b, 0, 0, -(a - b)*c*d + a - b, 0, -(a - b)*c*d + a - b, 3*(a + b)*(a - b)*d, 0),
(0, (a + b)*c*d + a + b, -3*c^2*d + 1, 0, (a + b)*c*d + a + b, 0, 0, -(a - b)*c*d + a - b, (a + b)*c*d + a + b, 0, 0, -(a - b)*c*d + a - b, 0, 3*(a + b)*(a - b)*d, -(a - b)*c*d + a - b, 0),
((a - b)*c*d + a - b, 0, 0, 3*c^2*d - 2*c + 1, 0, -(a + b)*c*d + a + b, -(a + b)*c*d + a + b, 0, 0, -(a + b)*c*d + a + b, -(a + b)*c*d + a + b, 0, (2*(a + b)^2 + (a - b)^2)*d, 0, 0, (a - b)*c*d + a - b),
(0, (a + b)*c*d + a + b, (a + b)*c*d + a + b, 0, -3*c^2*d + 1, 0, 0, -(a - b)*c*d + a - b, (a + b)*c*d + a + b, 0, 0, 3*(a + b)*(a - b)*d, 0, -(a - b)*c*d + a - b, -(a - b)*c*d + a - b, 0),
((a - b)*c*d + a - b, 0, 0, -(a + b)*c*d + a + b, 0, 3*c^2*d - 2*c + 1, -(a + b)*c*d + a + b, 0, 0, -(a + b)*c*d + a + b, (2*(a + b)^2 + (a - b)^2)*d, 0, -(a + b)*c*d + a + b, 0, 0, (a - b)*c*d + a - b),
((a - b)*c*d + a - b, 0, 0, -(a + b)*c*d + a + b, 0, -(a + b)*c*d + a + b, 3*c^2*d - 2*c + 1, 0, 0, (2*(a + b)^2 + (a - b)^2)*d, -(a + b)*c*d + a + b, 0, -(a + b)*c*d + a + b, 0, 0, (a - b)*c*d + a - b),
(0, -(a - b)*c*d + a - b, -(a - b)*c*d + a - b, 0, -(a - b)*c*d + a - b, 0, 0, -3*c^2*d + 1, 3*(a + b)*(a - b)*d, 0, 0, (a + b)*c*d + a + b, 0, (a + b)*c*d + a + b, (a + b)*c*d + a + b, 0),
(0, (a + b)*c*d + a + b, (a + b)*c*d + a + b, 0, (a + b)*c*d + a + b, 0, 0, 3*(a + b)*(a - b)*d, -3*c^2*d + 1, 0, 0, -(a - b)*c*d + a - b, 0, -(a - b)*c*d + a - b, -(a - b)*c*d + a - b, 0),
((a - b)*c*d + a - b, 0, 0, -(a + b)*c*d + a + b, 0, -(a + b)*c*d + a + b, (2*(a + b)^2 + (a - b)^2)*d, 0, 0, 3*c^2*d - 2*c + 1, -(a + b)*c*d + a + b, 0, -(a + b)*c*d + a + b, 0, 0, (a - b)*c*d + a - b),
((a - b)*c*d + a - b, 0, 0, -(a + b)*c*d + a + b, 0, (2*(a + b)^2 + (a - b)^2)*d, -(a + b)*c*d + a + b, 0, 0, -(a + b)*c*d + a + b, 3*c^2*d - 2*c + 1, 0, -(a + b)*c*d + a + b, 0, 0, (a - b)*c*d + a - b),
(0, -(a - b)*c*d + a - b, -(a - b)*c*d + a - b, 0, 3*(a + b)*(a - b)*d, 0, 0, (a + b)*c*d + a + b, -(a - b)*c*d + a - b, 0, 0, -3*c^2*d + 1, 0, (a + b)*c*d + a + b, (a + b)*c*d + a + b, 0),
((a - b)*c*d + a - b, 0, 0, (2*(a + b)^2 + (a - b)^2)*d, 0, -(a + b)*c*d + a + b, -(a + b)*c*d + a + b, 0, 0, -(a + b)*c*d + a + b, -(a + b)*c*d + a + b, 0, 3*c^2*d - 2*c + 1, 0, 0, (a - b)*c*d + a - b),
(0, -(a - b)*c*d + a - b, 3*(a + b)*(a - b)*d, 0, -(a - b)*c*d + a - b, 0, 0, (a + b)*c*d + a + b, -(a - b)*c*d + a - b, 0, 0, (a + b)*c*d + a + b, 0, -3*c^2*d + 1, (a + b)*c*d + a + b, 0),
(0, 3*(a + b)*(a - b)*d, -(a - b)*c*d + a - b, 0, -(a - b)*c*d + a - b, 0, 0, (a + b)*c*d + a + b, -(a - b)*c*d + a - b, 0, 0, (a + b)*c*d + a + b, 0, (a + b)*c*d + a + b, -3*c^2*d + 1, 0),
(3*(a - b)^2*d, 0, 0, (a - b)*c*d + a - b, 0, (a - b)*c*d + a - b, (a - b)*c*d + a - b, 0, 0, (a - b)*c*d + a - b, (a - b)*c*d + a - b, 0, (a - b)*c*d + a - b, 0, 0, 3*c^2*d + 6*c + 1)]
)

edit retag flag offensive close merge delete

Comments

Welcome to Ask Sage! Thank you for your question! Please provide a specific matrix to illustrate the question and help others study it.

slelievre ( 2019-02-28 11:51:53 +0200 )edit

Thanks! I have tried to include the example, I hope it is at least copy-pastable (it's a 16x16 matrix after all)

angelolucia ( 2019-02-28 20:36:32 +0200 )edit

Thanks for providing an example. One way to provide the matrix would be to define

a, b, c, d = var('a b c d', domain='positive')
e, f, g, h, i = a + b, a - b, 3*c^2*d, f*c*d, e*c*d
A, B, C = g + 6*c + 1, -g + 1, g - 2*c + 1
D, E, F, G, H = 3*e*f*d, -i + e, i + e, -h + f, h + f
I, J = (2*e^2 + f^2)*d, 3*f^2*d

slelievre ( 2019-03-01 07:46:24 +0200 )edit

then

M = matrix(16, 16, [
(A, 0, 0, H, 0, H, H, 0, 0, H, H, 0, H, 0, 0, J),
(0, B, F, 0, F, 0, 0, G, F, 0, 0, G, 0, G, D, 0),
(0, F, B, 0, F, 0, 0, G, F, 0, 0, G, 0, D, G, 0),
(H, 0, 0, C, 0, E, E, 0, 0, E, E, 0, I, 0, 0, H),
(0, F, F, 0, B, 0, 0, G, F, 0, 0, D, 0, G, G, 0),
(H, 0, 0, E, 0, C, E, 0, 0, E, I, 0, E, 0, 0, H),
(H, 0, 0, E, 0, E, C, 0, 0, I, E, 0, E, 0, 0, H),
(0, G, G, 0, G, 0, 0, B, D, 0, 0, F, 0, F, F, 0),

slelievre ( 2019-03-01 07:51:13 +0200 )edit

(continued)

(0, F, F, 0, F, 0, 0, D, B, 0, 0, G, 0, G, G, 0),
(H, 0, 0, E, 0, E, I, 0, 0, C, E, 0, E, 0, 0, H),
(H, 0, 0, E, 0, I, E, 0, 0, E, C, 0, E, 0, 0, H),
(0, G, G, 0, D, 0, 0, F, G, 0, 0, B, 0, F, F, 0),
(H, 0, 0, I, 0, E, E, 0, 0, E, E, 0, C, 0, 0, H),
(0, G, D, 0, G, 0, 0, F, G, 0, 0, F, 0, B, F, 0),
(0, D, G, 0, G, 0, 0, F, G, 0, 0, F, 0, F, B, 0),
(J, 0, 0, H, 0, H, H, 0, 0, H, H, 0, H, 0, 0, A)])

(sorry, can't fit in one comment).

Or maybe you have a nicer way to produce the matrix.

slelievre ( 2019-03-01 07:52:01 +0200 )edit

see more comments

answered 2019-02-28 22:38:58 +0200

rburing

10644 ●4 ●77 ●211 https://www.rburing.nl/

updated 2019-02-28 22:40:02 +0200

You can change the base ring of the matrix to the ring of rational functions (as you did), factor the characteristic polynomial, and set the factors equal to zero symbolically again (and solve).

R = PolynomialRing(QQ, names='a,b,c,d')
M_rat = M.change_ring(R.fraction_field())
S.<x> = R[]
f = M_rat.characteristic_polynomial(x)
f_factors = f.factor()
linear_factors = [(g,e) for (g,e) in f_factors if g.degree() == 1]
rational_roots = [(solve(SR(g)==0,SR(x),solution_dict=True)[0][SR(x)],e) for (g,e) in linear_factors]
nonlinear_factors = [(g,e) for (g,e) in f_factors if g.degree() > 1]
nonrational_roots = sum([map(lambda eqn: (eqn.rhs(), e), solve(SR(g)==0,SR(x))) for (g,e) in nonlinear_factors], [])
eigenvalues = rational_roots + nonrational_roots

This runs in 33 seconds on my machine. The result is

[(3*(a^2 - b^2 + 2*b*c - c^2)*d + 6*a + 1, 1),
 (-3*(a^2 - b^2 - 2*a*c + c^2)*d + 6*b + 1, 1),
 (-3*(a^2 - 2*a*b + b^2 - c^2)*d + 6*c + 1, 1),
 ((3*a^2 + 2*a*b + 3*b^2 + 2*(a + b)*c + 3*c^2)*d - 2*a - 2*b - 2*c + 1, 2),
 ((3*a^2 - 3*b^2 - 2*b*c - 3*c^2)*d - 2*a + 1, 3),
 (-(3*a^2 - 3*b^2 + 2*a*c + 3*c^2)*d - 2*b + 1, 3),
 (-(3*a^2 + 2*a*b + 3*b^2 - 3*c^2)*d - 2*c + 1, 3),
 ((3*a^2 - 2*a*b + 3*b^2 - 2*(a + b)*c + 3*c^2)*d + 2*a + 2*b + 2*c - 4*sqrt((a^2*b^2 + (a^2 - a*b + b^2)*c^2 - (a^2*b + a*b^2)*c)*d^2 + a^2 - a*b + b^2 - (a + b)*c + c^2 + (a^2*b + a*b^2 + (a + b)*c^2 + (a^2 - 6*a*b + b^2)*c)*d) + 1,
  1),
 ((3*a^2 - 2*a*b + 3*b^2 - 2*(a + b)*c + 3*c^2)*d + 2*a + 2*b + 2*c + 4*sqrt((a^2*b^2 + (a^2 - a*b + b^2)*c^2 - (a^2*b + a*b^2)*c)*d^2 + a^2 - a*b + b^2 - (a + b)*c + c^2 + (a^2*b + a*b^2 + (a + b)*c^2 + (a^2 - 6*a*b + b^2)*c)*d) + 1,
  1)]

Maybe an algorithm like this could be included into Sage in the future.

(The hard part is to detect when this can be applied.)

edit flag offensive delete link

Comments

Thank you very much! This is indeed a pretty good idea. I have noticed that doing computations with the polynomial ring Sage uses Singular (while I guess it is using Maxima for the symbolic ring). Even computing the characteristic polynomial over the symbolic ring is too much for my machine, but the answer should be the same in both cases, right?

Regarding the applicability of this approach in general: I suppose a polynomial matrix could have a lot (all?) eigenvalues that are non-rational functions, so this approach would fail in that cases: but wouldn't it reduce to the original symbolic calculations in that case? In other words, can this actually be harmful?

angelolucia ( 2019-02-28 23:15:40 +0200 )edit

You're welcome! Yes, the characteristic polynomial can be computed in either ring, and when this can be applied it is never harmful.

The matrix Matrix(SR,[[0,0,0,0,1], [1,0,0,0,1], [0,1,0,0,0], [0,0,1,0,0], [0,0,0,1,0]]) is 'multiplication by $\alpha$' w.r.t. the basis $1,\alpha,\alpha^2,\alpha^3,\alpha^4$ of $\mathbb{Q}(\alpha)$ where $\alpha$ is a root of $x^5 - x - 1$. So its characteristic polynomial is $x^5 - x - 1$ with Galois group $S_5$ which is not solvable, so the roots cannot be expressed in terms of radicals.

This is a special case where there are no variables; then one should factor over $\mathbb{Q}$ instead of a polynomial ring. The polynomial won't factor in this case, and the subsequent attempt at a symbolic solution also won ...(more)

rburing ( 2019-03-01 00:17:39 +0200 )edit

add a comment

Eigenvalues over symbolic ring

Comments

1 Answer

Comments

Your Answer

Question Tools

Stats

Related questions

Eigenvalues over symbolic ring edit

Comments

1 Answer

Comments

Your Answer

Question Tools

Stats

Related questions

Eigenvalues over symbolic ring