How to interprete an analytical eigendecomposition for a polynomial matrix with size nxn > 4x4 using the Symbolic Math Toolbox.

Hi,

I am using the Symbolic Math Toolbox to calculate the analytical eigenvalues of a polynomial matrix A(z), that is generated in my physical use-case. Depending on the input parameters, the matrix can scale in size from 1x1 to any nxn.

If I configure my matrix to have a size nxn > 4x4, the symbolic math toolbox still produces analytical eigenvalues (i.e. a polynomials without the "root" expression). How does Matlab calculate the analytical Eigenvalues in this particular case?

According to an answer to my previous question: https://www.mathworks.com/matlabcentral/answers/2132636-how-does-matlabs-symbolic-math-toolbox-always-finds-an-analytical-eigendecomposition-for-arbitrary-p#answer_1478061, Matlab uses numerical methods to find the roots of the characteristic polynomial.

However, from my experience the symbolic toolbox does not provide an analytic expression (i.e. a polynomial), if numerical methods where used and the solution is purely numerical. So my question is how Matlab arrives at this result?

1. Does Matlab indeed use a numerical approximation?
2. Does Matlab do some "tricks" in the background to reformulate my matrix into some analytically solvable matrix?

I also have attached a copy of my Matrix in a 6x6 configuration and the corresponding result.

Note: the variable "y" in my matrix represents 1/z, but this shouldn't make a difference.

My code does the following:

syms z y
A = A;
A_roots = eig(A);

[Torsten]

If I configure my matrix to have a size nxn > 4x4, the symbolic math toolbox still produces analytical eigenvalues (i.e. a polynomials without the "root" expression).

Can you show us an example where this happens ? It might produce the characteristic polynomial, but how should it be possible to get analytical expressions for the roots ? It would contradict Abel-Ruffini.

However, from my experience the symbolic toolbox does not provide an analytic expression (i.e. a polynomial), if numerical methods where used and the solution is purely numerical.

I don't understand what you mean here.

So my question is how Matlab arrives at this result?

Which result ? Matlab only provides the polynomial, but not its roots.

[Bastian Loß]

1. The attached file "Example_matrix.mat " contains a 6x6 polynomial matrix . Where y serves as a placeholder for . Calling eig(A) produces A_roots, which are the eigenvalue-polynomials of the matrix A in terms of z and y. Atleast from my understanding Matlab therefore did in fact obtained an analytical solution for this case.
2. As discussed in: https://www.mathworks.com/matlabcentral/answers/2132636-how-does-matlabs-symbolic-math-toolbox-always-finds-an-analytical-eigendecomposition-for-arbitrary-p#answer_1478061, it seems to me that Matlab does provide a "root(...)" expression in case it can not find analytical eigenvalues.
3. From my understanding it did provide the roots, as the characteristic polynomial would be and its roots are the eigenvalues, as these are only dependent on z and y.

[Torsten]

The attached file "Example_matrix.mat " contains a 6x6 polynomial matrix . Where y serves as a placeholder for . Calling eig(A) produces A_roots, which are the eigenvalue-polynomials of the matrix A in terms of z and y. Atleast from my understanding Matlab therefore did in fact obtained an analytical solution for this case.

eig(A) will produce a polynomial p of degree 6 by symbolically computing det(A-lambda*I) and return root(p). Is this a solution ? Formally yes, but not very helpful.

So if you consider the computation of this polynomial of degree 6 in lambda as "analytical solution", you are correct. Usually, the representation of the roots of this polynomial (that will depend on z and y) is referred to as "analytical solution" and these roots are the key of interest (because they are the eigenvalues).

[Bastian Loß]

Quoting my answer from below: If I use a random 6x6 matrix, the algorithm will produce a root expression - essentially telling: "I can't root this characteristic polynomial - here you have a placeholder".

But in my case: it does always return a polynomial in z and y and not a root expression. My matrix is constructed from the modelling of a physical process, and scales (depending on the inputs) to . No matter how I choose n (except n = 2), the algorithm will produce an actual polynomial. The attached example shows this - A_roots is the polynomial for the A input.

So my current understanding is that there are cases where an analytical solution for an matrix can be obtained , i.e. when exploiting the matrix structure. In this case Matlab returns an actual polynomial by computing det(A-lambda*I) and return root(p).

However, if there is no analytical solution, Matlab will just return a root-expression and essentially tell me: "I can't root this characteristic polynomial - here you have a placeholder".

Therefore, it looks like my matrix has some property, that Matlab can somehow exploit to produce the actual eigenvalues. Or why can Matlab produce a polynomial in my case, but not for a random matrix?

[Torsten]

I cannot run your code here because it takes too long.

As you can see below for a simpler example, eig(A) returns the roots of the characteristic polynomial for A, but again in a form that doesn't help. You need to specifiy numerical values for y and z in your code and use vpa to get the numerical roots of the resulting 5th order polynomial. @John D'Errico did this in his answer.

If you really get analytical expressions for the roots of your complicated 11x11 matrix A that depend on two other symbolic variables y and z, you are really in luck and this A- matrix must be very special.

syms x y
A = [x y 1 5 x*y;2 3 x y -2;1 x^2 y^4 -2 x*y^2;x^2+y^2 -2 1 9 10;3 -0.5 y^2 x^2 5];
eig(A)

ans = 

[Bastian Loß]

If you really get analytical expressions for the roots of your complicated 11x11 matrix A that depend on two other symbolic variables y and z, you are really in luck and this A- matrix must be very special.

Okay, that answers my question. Thanks for the effort @Torsten - I much appreciate your and @John D'Errico's help. If I find out which property that is, I will let you know.

[Torsten]

Maybe you could insert the 11 symbolic eigenvalues you get as a graphics - I'm quite interested.

[Bastian Loß]

While attempting to plot the results, I have checked my own Example_matrix.mat , and the solution there is in-fact not analytical. The solution I provided did also use the "root(...)" expression to represent the eigenvalues.

Therefore, everything makes sense now and fits into the in-depth explanation that @John D'Errico provided earlier.

I have re-run a few tests with my framework to check if I can find any case for > 4 matrix with a polynomial solution. However, all of them ended up with "root(...)" solutions.

Therefore: You were correct all along, I just misread the results for the 6x6 case. Again thanks for your patience in explaining the matter.

Answer 1

[John D'Errico]

Ok, do you understand that you CANNOT compute the roots of a polynomial, if the degree of that polynomial is greater than 4? Well, at least, unless the roots are trivial. For example, yes, this works.

syms x
solve(x^6 - 1 == 0)

ans = 

But in general, if the polynomial is non-trivial, then there is no general solution, and this is provably true? (Talk to some folks named Abel and Ruffini if that upsets you. They are both long dead of course, so you may need a Ouija board, or something as useful.)

Next, you need to understand that computing the eigenvalues of a matrix is mathematically equivalent to solving for the roots of a polynomial. So a 5x5 or 6x6 matrix has a chacteristic polynomial, and the roots of that polynomial are the same as the eigenvalues of the matrix. There is a complete correspondence between the two problems. If you can solve one, you can solve the other. But I just got done telling you that you CANNOT compute the roots of a general symbolic polynomial of degree 5 or higher, where the result has an algebraic solution. For degree 4 or below, the problem is solved, and the solutions have also existed for many years. That means, effectively, the eigenvalues of even a fully general 4x4 matrix can be found in algebraic form. They will be a little messy, but they are still trivially written in theory. For 5x5 or above, nothing can be done. NOTHING.

Do you see where this is going? Just wanting to solve a problem that is mathematically impossible to solve, where it has been proved to be the case. Well, the need or the desire is not sufficient. Sorry. This is not a question of MATLAB. It is a question of mathematical impossibility.

load Example_matrix
whos

  Name              Size            Bytes  Class             Attributes

  A                 6x6                 8  sym                         
  A_evs             6x1                 8  sym                         
  ans               6x1                 8  sym                         
  cmdout            1x33               66  char                        
  gdsCacheDir       1x14               28  char                        
  gdsCacheFlag      1x1                 8  double                      
  i                 0x0                 0  double                      
  managers          1x0                 0  cell                        
  managersMap       0x1                 8  containers.Map              
  x                 1x1                 8  sym                         

It appears your matrix is A. That is a good name.

A

A = 

And now you want to compute the eigenvalues of A, in an algebraic form, as a function of y AND of z. Sorry. Not gonna happen. EVER. Again, if you have a problem, all you can do is talk to either Abel of Ruffini. Maybe you can convince them that you really, very much need a solution, and that you need them to rewrite the proof. This may take a truly cosmic intervention though. Do you have the right connections for such a task? I certainly don't.

If you knew the values of both y and z. Well, then you can compute the numerical roots.

syms y z
vpa(eig(subs(A,[y,z],[3,17])))

ans = 

But you need to see at that point, the result is just a numerical matrix. Now a rootfinder has no problem handling what is essentially a polynomial rootfinding problem.

In terms of your original matrix A though, there is NO numerical solution that can be applied. There is no magic. You really cannot do more than that. I'm sorry. Mathematics can be a cruel mistress at times.

[Bastian Loß]

Hello John,

I am very well aware of the Abel-Ruffini theorem. But my question remains: How does Matlab then compute "eig(A)" in a case of a 6x6 matrix? Does it use check if the matrix can be somehow simplified or re-writte or does straight up resort to using numerics?

Update: Thanks for the update John, your Text was a joy to read :)

But I am still confused, Matlab does compute "something" i.e. running A_roots = eig(A); produces 6x polynomials. Assuming there is in fact no numerical solution and the problem is infeasible, what does Matlab compute then?

[John D'Errico]

How much does it try to simplify the matrix? Not even at all hard, if it does anything at all. For example,...

syms x
A = sym(rand(5));
eig(sym(x*A))

ans = 

So it is trivial to compute the eigenvalues of A (although computing them in an algebraic form is a mathematical impossibiity, as we should know.) Was it smart enough to factor out x, to know how the simple multiplicative factor x modifies the eigenvalues? Nope.

It is often the case that MATLAB almost looks like it is doing something, returning a rootof expression. For example:

syms a y
P = y^5 + a*y + 2;
ysol = solve(P == 0,y)

ysol = 

It tells us that IF it could solve the problem, which it apparently does not fully understand that mathematically, it cannot, the result would be one of the 5 roots of that polynomial. It is a non-answer though.

[Bastian Loß]

I have observed this behaviour as well, but this is essentially where my confusion comes from.

If I use a random 6x6 matrix, the algorithm will produce a root expression - essentially telling: "I can't root this characteristic polynomial - here you have a placeholder".

But in my case: it does always return a polynomial in z and y and not a root expression. My matrix is constructed from the modelling of a physical process, and scales (depending on the inputs) to . No matter how I choose n (except n = 2), the algorithm will produce an actual polynomial. The attached example shows this - A_roots is the polynomial for the A input.

Therefore, I thought Matlab would maybe use some "tricks" to transform my matrix beforehand.

Update: Or my Matrix maybe has some special properties that Matlab exploits to come up with an actual solution.