# Step by Step for Generalized Eigenvectors?

When I do

 WolframAlpha["Eigensystem[{{1,-3},{3,-5}}]", PodStates -> {"Step-by-step solution"}]


It only shows steps for non-generalized eigenvectors.

Is there any way to get the step by step from WA when you have generalized eigenvectors?

Note: I know how to do this by hand, so no solution needed for this specific problem, I was just curious if there was any way to get it out of MMA through WA or other. I would like to see it for an $$nxn$$ matrix (small number of eigenvalues, say up to $$5$$).

Maybe there is a more general approach that does not require a call to WA, so that is welcome if it is possible.

Update:

As mentioned in comments, Eigensystem does not return generalized eigenvectors, so is there a general way to get those with steps? Certainly I can use things like LinearSolve to find eigenvectors, but that leads to a solution with no steps. If I were going to do that, I would just use JordanDecomposition[A], which is much easier but suffers from the same issue of no steps.

• This is something I was never clear on: Do you happen to know if when using WolframAlpha command within a copy of Mathematica, if it will use the PRO version of Wolfram alpha or the free version? If the free version, then there is a pro version of wolfram alpha also to try. But my understanding of step-by-step solutions, is that it is done for basic problems at the level of high school or first year college level, and not the more advanced ones. Try ChatGPT, may be AI will do it :) Commented May 9, 2023 at 20:44
• @Nasser: Thanks for the feedback. I believe we get WA Pro because that is the only thing that reruns results at all. I can see solutions for all sorts of things - including DEQs. I always found it strange that it can return regular eigenvalues but not generalized ones, but maybe the internal CAS algorithms do not find them that way we do by hand and can thusly not give them. I was hoping that maybe the manual method can be mechanized and displayed since it is almost algorithmic for generalized eigenvectors.
– Moo
Commented May 9, 2023 at 20:52
• Are you sure Eigensystem is supposed to return generalized eigenvectors for defective matrix? I think it will not. Googling around, one needs to add extra code to solve for these themselves. Here is a link showing how to do this where it says To find the generalized eigenvector corresponding to λ=1, we use the following Mathematica command so the steps the above command returned (and it does return step-by-step on V 13.2.1) is all what you will get, because it is not meant to return generalized V anyway. right? Commented May 9, 2023 at 21:16
• Eigensystem does not return the generalized eigenvectors. Commented May 9, 2023 at 21:42
• @Nasser: That is the point of the post, that Eigensystem does not return generalized eigenvectors, so is there a a general way to get those with steps. Certainly I can do things like LinearSolve, but that leads to a solution with no steps. If I were going to do that, I would just use JordanDecomposition[A], which is much easier but suffers from the same issue of no steps.
– Moo
Commented May 10, 2023 at 1:47

This is an update to make it handle geometric multiplicity>1.

When an eigenvalue has geometric multiplicity=1, then this will use the basic method of solving $$(A-\lambda I) v_2 = v_1$$ to find the generalized eigenvector $$v_2$$ given $$v_1$$, and repeats this as many times as needed per each defective $$\lambda$$.

To find the first generalized $$v_2$$, in the above $$v_1$$ will be the normal eigenvector, which is found by Eigensystem. If more generalized eigenvectors are needed (because the eigenvalue could have large defect), then the next step will use the now found $$v_2$$ in order to find $$v_3$$ and so on. This is done for each defective $$\lambda$$.

When an eigenvalue has geometric multiplicity>1 however, the generalized eigenvector is found from the null space of the transpose of the matrix made up of the current set of known eigenvectors. We just pick one random row from the nullspace and use it (after transposing) as new generalized eigenvector and add it to the set of known eigenvectors for that $$\lambda$$. If we need more generalized eigenvectors, we repeat this process. The function NullSpace is used in Mathematica for this.

## Example 1

mat = {{4, 1, 1}, {-2, 1, -2}, {1, 1, 4}};
findV[mat, V];


Gives

## Example 2

mat = {{3, 1}, {0, 3}}
findV[mat, V];


## Example 3

mat = {{2, 1}, {-1, 4}};
findV[mat, V];


## Example 4

mat = {{1, 0, 0, 0, 0}, {3, 1, 0, 0, 0}, {6, 3, 2, 0, 0}, {10, 6, 3,
2, 0}, {15, 10, 6, 3, 2}};
findV[mat, V];


## Example 5

mat = {{5, 1, -4}, {4, 3, -5}, {3, 1, -2}};
findV[mat, V];


## Example 6

mat = {{5, 1, 3, 2}, {0, 5, 0, -3}, {0, 0, 5, 1}, {0, 0, 0, 5}};
findV[mat, V];


## Example 7

mat = {{39, 8, -16}, {-36, -5, 16}, {72, 16, -29}};
findV[mat, V];


## Example 8

mat = {{1, 0, 0}, {1, 3, 1}, {-2, -4, -1}};
findV[mat, V];


## Example 9

mat = {{3, 1, 0, 0}, {0, 3, 1, 0}, {0, 0, 3, 0}, {0, 0, 0, 3}};
findV[mat, V];


Jordan form is unique up to permutation of its Jordan blocks.

When there is one block, like in example 9 above, since there is only one $$\lambda$$, the arrangement of columns of $$P$$ matrix (the matrix whose columns are the eigenvectors of this $$\lambda$$) is arbitrary. Only condition is that they are all linearly independent of each others. So you might not get same Jordan form as Mathematica gives if the order of eigenvectors in $$P$$ happened to be different than what Mathematica used. Here is for the above.

(A = {{3, 1, 0, 0}, {0, 3, 1, 0}, {0, 0, 3, 0}, {0, 0, 0, 3}}) // MatrixForm


JordanDecomposition[A][[2]] // MatrixForm


Here is the $$P$$ matrix of the eigenvectors found for $$\lambda$$

(p = {{0, 1, 0, 0}, {0, 0, 0, 1}, {0, 0, 1, 0}, {1, 0, 0, 0}}) // MatrixForm


Now

 (Inverse[p] . A . p) // MatrixForm


This does not have same arrangement what Mathematica gave.

But if we simply swap the 3rd column with the 4th, now it will give same Joprdan matrix as Mathematica

(p = {{0, 1, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}, {1, 0, 0, 0}}) // MatrixForm


(Inverse[p] . A . p) // MatrixForm


## Code

Bug reports are welcome. In case code gets messed up when pasting here as I use code cell, here is link to the notebook also. Used V 13.2.1

(*Change history *)
(*may 10, 2023. Initial Version *)
(*may 10, 2023. Added check for geometric multiplicity >1 for now until supported*)
(*may 11, 2023. Added support for geometric multiplicity >1, improved formatting*)

findV[mat_,V_]:=Module[{data,keys},
{data,keys}=generateInitialTable[mat];
processTable[mat,data,keys,V]
]

(*****************************)

Print@Text[Style[
"We first use Mathematica to find the eigenvalues and the norma\n"<>
"eigenvectors and examine the result to determine which eigenvalues\n"<>
"are defective and which are not. This results in this table ",TextAlignment->Left]];

MapAt[MatrixForm[Transpose[#]]&,data,{All,4}]
],
Frame->All]]
]
(******************************)
{lambda,v}=Eigensystem[mat];
Print["The Matrix is ",MatrixForm@mat];
Print["its eigenvalues are ",MatrixForm@lambda];
Print["its normal Eigenvectors are ",MatrixForm@Transpose[v]];

p=PositionIndex[lambda];
keys=Keys[p];

(*this table will hold summary of Eigensystem finding. One row per each lambda*)
(*1=lambda,2=algebraic multiplicity,3=geometric multiplicity,4=normal V,5=defective?*)

data=Table[{0,0,0,0,0},Length@keys];

(*Now build the table *)
Do[
k=keys[[n]];
data[[n,1]]=k;
data[[n,2]]=Length@p[k];(*algebraic multiplicity*)
idx=p[k];
theV=v[[idx[[1]];;idx[[-1]]]];
theV=Cases[theV,x_/;Not@AllTrue[x,#==0&]];
data[[n,3]]=Length@theV;(*geometric multiplicity*)
data[[n,4]]=theV;
data[[n,5]]=If[data[[n,2]]>data[[n,3]],True,False]
,
{n,Length@keys}
];

printTable[data];
{data,keys}
]
(******************************)
(*geometric multiplicity =1, easy case *)
handleCaseOne[data_,mat_,V_,n_]:=Module[{howMuchDefective,currentV,g,eq,m,size=Length@mat,sol,setOfV},
howMuchDefective=data[[n,2]]-data[[n,3]];
Print["Since eigenvalue ",data[[n,1]]," is defective, and the degree of defect is ",
howMuchDefective," then we need to find ",howMuchDefective,
" additional generalized eigenvectors for this eigenvalue."];
currentV=data[[n,4]]; (*normal V found*)
setOfV=currentV;

Do[
g=Table[V[i],{i,1,size}];
eq=mat.g==currentV;
Print["Now we solve for generalized eigevector ",MatrixForm[g],
" Using (A-\[Lambda]I).g=v1. This gives"];
Print[MatrixForm[mat],"-(",data[[n,1]],")",
MatrixForm[IdentityMatrix[size]],".",MatrixForm[g],"=",MatrixForm[Transpose@currentV]];

Print["The above equation becomes\n",MatrixForm[mat-data[[n,1]]*IdentityMatrix[size]],
".",MatrixForm[g],"=",MatrixForm[Transpose@currentV]];

sol=LinearSolve[(mat-data[[n,1]]*IdentityMatrix[size]),currentV[[1]]];
Print["Using LinearSolve gives ",MatrixForm[g],"=",Style[MatrixForm@sol,Bold]];
currentV={sol};
setOfV=Join[setOfV, currentV]
,
{m,1,howMuchDefective}
];

setOfV
]

(******************************)
(* when geometric multiplicity >1, hard case *)
handleCaseTwo[data_,mat_,V_,n_]:=Module[{howMuchDefective,currentV,g,eq,m,size=Length@mat,sol,newV,i},
Print["This case of geometric multiplicity >1."];

howMuchDefective=data[[n,2]]-data[[n,3]];
Print["Since eigenvalue ",data[[n,1]]," is defective, and the\n",
"degree of defect is ",howMuchDefective," then we need to find ",howMuchDefective,"\n",
"additional generalized eigenvectors for this eigenvalue."];
currentV=data[[n,4]]; (*normal V's found. More than one*)
g=Table[V[i],{i,1,size}];

Do[
Print["Now we solve need to find a generlized eigenvector ",MatrixForm@g,"\n",
"which is linearly independent to the set of eigenvectors already found ",
MatrixForm[Transpose[currentV]],"\n",
"Which we do by looking at the NullSpace of the transpose of the above as a matrix\n",
"Which gives ",MatrixForm[NullSpace[currentV]],"\n",
"And we need to just pick one of these rows as our new eigenvector.\n"
];
newV=NullSpace[currentV];
newV=newV[[1,All]];
currentV=Join[currentV,{newV}];
Print["Taking the first row, and transposing it gives the generalized eigenvector as ",MatrixForm@g," = ",
Style[MatrixForm[Transpose[{newV}]],Bold],"\nwhich we now add to the set of known eigenvectors found giving\n",
MatrixForm[Transpose@currentV]
];

,
{m,1,howMuchDefective}
];

currentV
]
(******************************)
processTable[mat_,data_,keys_,V_]:=Module[{n,setOfV},
Do[
If[data[[n,5]]==True,(*i.e. defective eigenvalue*)
(*check if geometric multiplicity=1 or not. These cases are handled differently*)

If[data[[n,3]]==1,(*standard case, where one normal vector is found*)

(*geometric multiplicity =1 case*)
setOfV=handleCaseOne[data,mat,V,n];
Print["Hence the set of eigenvectors associated with eigenvalue ",data[[n,1]]," is ",
Style[MatrixForm[Transpose@setOfV],Bold]
]
,
(*geometric multiplicity >1 case*)
setOfV=handleCaseTwo[data,mat,V,n];
Print["Hence the set of eigenvectors associated with eigenvalue ",data[[n,1]]," is ",
Style[MatrixForm[Transpose@setOfV],Bold]
]
]
,
Print["Since eigenvalue ",data[[n,1]]," is not defective, then we already have its normal "];
Print["eigenvectors from the table above,  which was found by Eigensystem. Nothing to do."]
]
,
{n,Length@keys}
]
]

• That is awesome! Here is a failed case, see example 2 of math.purdue.edu/~neptamin/303Au21/Handouts/High_defect.pdf . It looks like we get two normal eigenvectors and that seems to be an issue (maybe just needs to use one?). Here is the input: mat = {{3, 1, 0, 0}, {0, 3, 1, 0}, {0, 0, 3, 0}, {0, 0, 0, 3}}; {data, keys} = doPhaseOne[mat]; data = doPhaseTwo[mat, data, keys, V]; Example 3 also fails for the same reason.
– Moo
Commented May 10, 2023 at 10:09
• Additional failure examples (same as above) for test are 4 and 5, mathcs.holycross.edu/~spl/old_courses/304_fall_2008/handouts/…
– Moo
Commented May 10, 2023 at 10:30
• Might need to solve using a linear combination of the eigenvectors with undetermined coefficients. Commented May 10, 2023 at 13:58
• @Moo, Yes, I will fix it today. I forgot to check for the case where one defective $\lambda$ can have more than one normal eigevector. Yes, I need to do linear combination in this case before solving for the missing generalized ones. Will do this later and will also try to improve formatting more. Commented May 10, 2023 at 14:18
• @Moo I think it is correct. Jordan form is unique up to permutation of its blocks. But here we have one block since one $\lambda$. So It depends on how the eigenvectors of the one λ are ordered in P . This order is arbitrary. By changing the order, you can get same form as Mathematica. Will add note to my answer showing this. If you think there is still a problem, please let me know. Commented May 11, 2023 at 13:02