Analytic formula for leading eigenvector of $uu^T + vv^T$?

Question

Let $u$ and $v$ be nonzero column vectors of size $n$ and consider the $n \times n$ positive-definite matrix $A:=uu^T + vv^T$. In this post https://math.stackexchange.com/a/112201/168758, the eigenvalues of $A$ were computed analytically.

Question

I wonder whether there is an analytic formula for the eigenvectors of $A$.

Observation

This answer https://math.stackexchange.com/a/112197/87355 shows how to compute the eigenvalues of $A$ via 2 iterations of Gram-schmidt. The method also contains a blueprint for computing the eigenvectors, in principle. The issue is that this method gives an eigendecomposition in a basis $\mathcal B=\{u_1,u_2,\ldots\}$ which is different from the standard basis. With all of this in hand, how to get eigenvectors of $A$ in the standard basis efficiently (for example, without computing the other $n-2$ vectors $u_3,\ldots,u_n$ of the basis $\mathcal B$, and then doing a change of basis formula).

Update

Fibonatic's answer along with the refs above, here is some python code which solves the problem (hope it helps someone else):

def special_eig(u, v, tol=1e-9):
    """
    Computes leading eigenvalue and eigenvector of uu^T + vv^T

    Notes
    =====
    Many of the computations can become numerically unstable. Most of the code
    is to overcome these potential issues.
    """
    from math import sqrt
    u2 = u.dot(u)
    v2 = v.dot(v)
    uv = u.dot(v)

    # check for linear dependence
    if u2 <= tol:
        return v2, v
    if v2 <= tol:
        return u2, u.copy()
    tmp = uv ** 2 / (u2 * v2)
    if abs(tmp - 1) <= tol:
        return u2 + v2, u.copy()

    # at this point, u and v are linear independent
    disc = sqrt((u2 - v2) ** 2 + 4 * (uv) ** 2)
    eigval = .5 * (u2 + v2 + disc)
    if abs(uv) <= tol:
        x = uv / (eigval - v ** 2)
    else:
        x = (eigval - u2) / uv
    eigvec = u + x * v
    return eigval, eigvec

Answer 1

One can start by writing the eigenvalue problem

$$ \left(A - \lambda_i\,I\right) w_i = 0, \tag{1} $$

with $\lambda_i$ one of the two eigenvalues and $w_i$ its eigenvector. By using the definition of $A$ $(1)$ can also be written as

$$ \langle u, w_i\rangle u + \langle v, w_i\rangle v = \lambda_i\,w_i. \tag{2} $$

This implies that $w_i$ has to be a linear combination of $u$ and $v$, so $w_i = \alpha_i\,u + \beta_i\,v$. Plugging this into $(2)$ yields

$$ \langle u, \alpha_i\,u + \beta_i\,v\rangle u + \langle v, \alpha_i\,u + \beta_i\,v\rangle v = \lambda_i (\alpha_i\,u + \beta_i\,v). \tag{3} $$

In the case that $u$ and $v$ are linearly dependent (so $v = x\,u$ with $x\in\mathbb{R}$) the matrix can also be written as $A = (1+x^2) u\,u^\top$. This means that it has only one non-zero eigenvalue and the corresponding eigenvector would just be $u$.

In the case that $u$ and $v$ are linearly independent $(3)$ can be factored into the following two scalar equation which are multiplied by $u$ and $v$ respectively

$$ \langle u, \alpha_i\,u + \beta_i\,v\rangle = \lambda_i\,\alpha_i, \tag{4a} $$ $$ \langle v, \alpha_i\,u + \beta_i\,v\rangle = \lambda_i\,\beta_i, \tag{4b} $$

which should be equivalent equations, because $|w_i|$ should remain undetermined. For example when defining $x_i = \beta_i/\alpha_i$ it can be solved for using $(4a)$ and $(4b)$ respectively

$$ x_i = \frac{\lambda_i - |u|^2}{\langle u,v\rangle}, \tag{5a} $$ $$ x_i = \frac{\langle u,v\rangle}{\lambda_i - |v|^2}. \tag{5b} $$

Choosing $\alpha_i=1$ gives $w_i = u + x_i\,v$ (one could normalize $w_i$ if desired). In the case that $\langle u,v\rangle = 0$ $(5a)$ is not well defined and similarly for $(5b)$ when $\lambda_i = |v|^2$. However, in both cases $u$ and $v$ are eigenvectors of $A$, since if $\langle u,v\rangle = 0$ then $A\,v = |v|^2 v$ and $A\,u = |u|^2 u$.

It can be shown that $\lambda_i=u^2$, $\lambda_j=v^2$ ($i\neq j$) and $\langle u,v\rangle=0$ are all equivalent. For this I assume that $|u|^2 \geq |v|^2$ (if this is not the case one could always swap the definitions of $u$ and $v$) in which case there should always exist a $x\geq0$ such that

$$ \sqrt{(|u|^2 - |v|^2)^2 + 4 \langle u,v\rangle^2} = \pm (|u|^2 - |v|^2 + 2\,x), $$

with $x=0$ only when $\langle u,v\rangle=0$. Substituting this into the equation for the eigenvalues yields

$$ \lambda_i = \frac{|u|^2 - |v|^2 \pm (|u|^2 - |v|^2 + 2\,x)}{2}, $$

which can be simplified to $\lambda_1 = |u|^2 + x$ and $\lambda_2 = |v|^2 - x$. However, since $|u|^2 \geq |v|^2$ and $x\geq0$ implies that $|u|^2$ can only be an eigenvalue if $x=0$ and thus $|v|^2$ has to be an eigenvalue as well and $\langle u,v\rangle=0$ and vice versa.

One side note, namely this does demonstrate an analytical solution for the eigenvectors, but it might not always yield accurate results when implementing this directly numerically. For example when $\lambda_i$ is very close to but not exactly equal to $|u|^2$ or $|v|^2$ in $(5a)$ and $(5b)$ respectively, then one could lose quite some numerical accuracy.