Issue regarding the result from symbolic computation

Hello,

Currently, I am trying to use symbolic computation in the process of solving the eigenvalue problem.

I need to do this because I need the solution of the eigendecomposition problem to be expressed using ingredients of the original matrix.

However, the result of algebraic expression from symbolic computation is not reliable.

For instance,

syms  a11 a12 a22 
A = [a11 a12 0; a12 a22 0; 0 0 0];   

[R, EV] = eig(A);  % Calculate eigenvalues and eigenvectors
[EV,ind] = sort(diag(EV),'descend');
R2= R(:,ind);

gives

R2 = [(a11/2 + a22/2 + (a11^2 - 2*a11*a22 + 4*a12^2 + a22^2)^(1/2)/2)/a12 - a22/a12 (a11/2 + a22/2 - (a11^2 - 2*a11*a22 + 4*a12^2 + a22^2)^(1/2)/2)/a12 - a22/a12 0;
1 1 0
;
0 0 1;];

However, this is not correct as if I give it a11 = 1, a12 = 0, a22 = 0, the matrix is

R2 =

   NaN   NaN     0
     1     1     0
     0     0     1

Although, the correct result is

 R2 =

     1     0     0
     0     1     0
     0     0     1

It is obvious that I get NaN values as the algebraic expression has the term with zero division.

Moreover, in the case a11 = 1/2, a12 = 1/2, a22 = 1/2, the result from symbolic computation is

 R2 =

     1    -1     0
     1     1     0
     0     0     1

Although, the correct result is

   R2 =

    0.7071   -0.7071         0
    0.7071    0.7071         0
         0         0    1.0000

Does anyone know the solution to circumvent this issue?

Answer 1

[Paul]

Hi ArxIv,

Expectations might be too high for symbolic solutions of eigenvectors.

Consider a simpler, but very closely related problem:

a11 = 5;a22 = 10;
syms a12
A = [a11,a12;a12,a22];
[V,D] = eig(A);
diag(D)
disp(char(ans))

[15/2 - (4*a12^2 + 25)^(1/2)/2; (4*a12^2 + 25)^(1/2)/2 + 15/2]

V = simplify(V)
disp(char(V))

[-((4*a12^2 + 25)^(1/2) + 5)/(2*a12), ((4*a12^2 + 25)^(1/2) - 5)/(2*a12); 1, 1]

I believe that the eigenvalues will yield the correct result for any value of a12, as will the eigenvectors except for the problematic case if we attempt to evaluate them at a12 = 0. The question is if there exists a general expression for V that works for all possible values of a12. If not, then the ideal solution would be for eig to return a piecewise expression of some sort for V taking into account all possible special cases. Other functions, like int also typically do not identify and return piecewise solutions for special cases. I'd imagine it may be difficult to find all special cases for more complicated problems. So I guess eig does the best it can.

Sometimes, we can get the solution by taking the limit as a12 ->0, which works in this case for V(:,2)

limit(V(:,2),a12,0),disp(char(ans))

[0; 1]

But that "trick" doesn't work as well for V(:,1)

limit(V(:,1),a12,0),disp(char(ans))

[NaN; 1]

though it kind of does if we realize that V(1,1) is blowing up to infinity, and so the normalized eigenvector is, in fact, approaching [1;0].

Also, using sort on symbolic expressions is problematic.

syms a11 a22 real
sort([a11,a22],'descend'),disp(char(ans))

[a22, a11]

Is there any mathematical reason to expect a22 >= a11? What happens if we sub values for a11 and a22 into that result when the latter is smaller than the former?

In addition, sort apparently doesn't respect assumptions

assumeAlso(a11 > a22);
sort([a11,a22],'descend'),disp(char(ans))

[a22, a11]

[John D'Errico]

@Paul has made some very valid points in terms of your expectations for the eigenvectors.

[ArxIv]

@John D'Errico @Paul

Thank you very much for your answers.

I learned a lot.

Answer 2

[John D'Errico]

If you have this:

a11 = 1/2; a12 = 1/2; a22 = 1/2;
A = [a11 a12 0; a12 a22 0; 0 0 0]

A = 3×3

0.5000 0.5000 0 0.5000 0.5000 0 0 0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

The eigenvectors were found as:

[V,D] = eig(A)

V = 3×3

-0.7071 0 0.7071 0.7071 0 0.7071 0 1.0000 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

D = 3×3

0 0 0 0 0 0 0 0 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Vs,Ds] = eig(sym(A))
Vs =
[-1, 0, 1]
[ 1, 0, 1]
[ 0, 1, 0]
Ds =
[0, 0, 0]
[0, 0, 0]
[0, 0, 1]

Is the set of eigenvectors still a valid set? YES. Eigenvectors can be arbitrarily scaled, so that is v is an eigenvector, then k*v is still just as validly an eigenvector. Note that the symbolic eig likes to scale its eigenvectors differently from the numerical version. Complain to the author, not to us.

Regardless, you should see the result from the symbolic compiutation is as valid as the result from the numerical one.

In the fully symbolic case, doing exactly as you said, we see:

syms  a11 a12 a22 
A = [a11 a12 0; a12 a22 0; 0 0 0];   
[Vs,Ds] = eig(A)

Now substitute

subs(Vs,[a11,a12,a22],[1,0,0])

Error using sym/subs (line 156)
Division by zero.

and we get a 0/0 result, which is properly a NaN. However, if I plug in those values for a11,a12, and a22, first, then we get a reasonable result.

[Vs,Ds] = eig(subs(A,[a11,a12,a22],[1 0 0]))
Vs =
[0, 0, 1]
[1, 0, 0]
[0, 1, 0]
Ds =
[0, 0, 0]
[0, 0, 0]
[0, 0, 1]

[John D'Errico]

Note that Answers is bugged as it often is lately, so symbolic results are not properly displayed.