How can cofactor matrix help find eigenvectors? Browsing around for math literature related to eigenvalues and eigenvectors I did not read yet, I came across the following algorithm. It seems to be functional for all 3x3 matrix I tried (symmetric, diagonal, triangular, etc. for a grand total of a few 10000 matrices), but I do not fully get the math logic behind...
Let $M$ be an arbitrary $3 \times 3$ matrix.
Let $\lambda_1$, $\lambda_2$, and $\lambda_3$ be the eigenvalues of $M$ (i.e.: the roots of its polynomial obtained by solving $det(M - \lambda I) = 0$).
Now, let
$
P = (M - \lambda_{i} I) = \begin{bmatrix}
a & b & c \\
d & e & f \\
g & h & i \\
\end{bmatrix}
$
The algorithm then creates three vectors $v_{ij}$ made of the determinant of the sub matrices $M[i,j;]$ as follow:
$
v_{01} = \begin{bmatrix}
   bf - ce  \\
 -(af - cd) \\
   ae - bd  \\
\end{bmatrix}
$
$
v_{02} = \begin{bmatrix}
   bi - ch  \\
 -(ai - cg) \\
   ah - bg  \\
\end{bmatrix}
$
$
v_{12} = \begin{bmatrix}
   ei - fh  \\
 -(di - fg) \\
   dh - eg  \\
\end{bmatrix}
$
From there, we find that the eigenvector associated to the current $\lambda_i$ is the $ \max \left ( |v_{01}|, |v_{021}|, |v_{12}| \right )$ dividing its corresponding vector to produce a normed eigenvector.
Q1: Could you explain and provide a bit more details about the math behind this algorithm?
I find some similarities with what I have done by substitution to express one unknown as a function of the others to inject it in a third equation and then solve the system by fixing a free unknown. The above algorithm looks, to me, like a great shortcut to obtain the eigenvectors. Furthermore, when I made my regular calculus there were some "corner cases" to check for in order to avoid, for instance, division by zero.. which seems not to be a problem here!
Q2: Although I doubt it, can this adapted to $2 \times 2$?
Q3: Can this be extended to higher matrix dimensions?
EDIT:
Staring at it long enough I realized that the vectors $v_{ij}$ form the rows of the cofactor matrix. I also found out that the odd lines of the cofactor matrix $P$ are timed by -1 (but this does not seem to change the final result).
Now I realize that this can be generalized to higher dimensions matrix (I went up to 5x5).
However, I still don't get the relation so Q1 is still a question..
Q4: why don't we always use this method to compute the eigenvectors? I mean, apart from the higher number of multiplications (which heavily impacts the performance though), this seems to be working for any types of matrix ..
NOTE: I renamed the title of the post accordingly.
 A: Q1. As we know, the eigenvector $\boldsymbol{v}$ must solve
$$
P\boldsymbol{v}=(M-\lambda I)\boldsymbol{v}=\begin{bmatrix}
a & b & c \\
d & e & f \\
g & h & i \\
\end{bmatrix}\boldsymbol{v}=\boldsymbol{0}\,.
$$
In other words, it must be orthogonal to all three vectors
$$
\boldsymbol{x}=\begin{bmatrix}
a  \\
b \\
c  \\
\end{bmatrix}\,,\quad
\boldsymbol{y}=\begin{bmatrix}
d  \\
e \\
f  \\
\end{bmatrix}\,,\quad
\boldsymbol{z}=\begin{bmatrix}
g  \\
h \\
i  \\
\end{bmatrix}\,.
$$
The existence of a non zero eigenvector is, as we know, equivalent to the matrix $P$ having rank at most two.
If we set $\boldsymbol{v}=\boldsymbol{x}\times\boldsymbol{y}$ (which is exactly what your algorithm is doing) then there are to possibilities:

*

*$\boldsymbol{x}$ and $\boldsymbol{y}$ are linearly independent. Then $\boldsymbol{v}$ is orthogonal to both, $\boldsymbol{x}$ and $\boldsymbol{y}$. Since $P$ has rank at most two this means that $\boldsymbol{z}$ must be a linear combination of $\boldsymbol{x}$ and $\boldsymbol{y}$. Then $\boldsymbol{v}$ is orthogonal to all three vectors and therefore an eigenvector. $\quad\Box$


*$\boldsymbol{x}$ and $\boldsymbol{y}$ are not linearly independent.
Then $\boldsymbol{v}=\boldsymbol{0}$ which is not an eigenvector that we are seeking.
We see that this algorithm does not always produce eigenvectors. The most simple case where it doesn't is $M=\operatorname{diag}(\lambda,\lambda,\lambda)$ in which case $P=0$ and $\boldsymbol{x}=\boldsymbol{y}=\boldsymbol{z}=\boldsymbol{0}\,.$
Q2. Q3. Don't think so. The algorithm relies heavily on the properties of the cross product in $\mathbb R^3$.
Q4. It should be obvious by now why we don't always use this method.
