Given matrix $A$ (symmetric and positive definite) and its Cholesky factors $A = LL^\top$, I want to obtain a new Cholesky factor $L_{new}$ such that
$$L_{new}L_{new}^\top = A + uv^\top + vu^\top$$
where $u,v$ are vectors.
Question: How can I do this efficiently?
Sub-optimal solution
I know rank-1 up/downdates can be done efficiently, so I could do
\begin{align*}
L_{new}L_{new}^\top = A + uv^\top + vu^\top = A + (u+v)(u+v)^\top – uu^\top – vv^\top
\end{align*}
so the required rank-2 update can be achieved via 3 separate rank-1 updates. Do I really need to do 3 rank-1 updates?
Best Answer
The update can be completed using two rank-1 updates subject to certain conditions. Let $B$ denote the perturbation given by $$B = uv^T + vu^T.$$ This is a symmetric matrix of rank at most $2$. Hence it can be written as $$B = \lambda_1 x_1x_1^T + \lambda_2 x_2x_2^T.$$ However, without further restrictions there is no guarantee that $B$ or even $A+B$ is symmetric positive semidefinite. In order to generate a counterexample, consider $v=-tu$ for a large value of $t > 0$. Then $$A + B = A - 2tuu^T$$ and $$u^T(A+B)u = u^T A u - 2 t \|u\|^2 \leq (\lambda_{\max}(A) - 2 t)\|u\|^2.$$ It remains to compute the eigenvalues and eigenvectors of $B$. Given the definition of $B$ we search for eigenvectors of the form $\alpha u + \beta v$ for real scalars $\alpha$ and $\beta$. Specifically, we wish to find $\alpha, \beta$ and $\lambda$ such that $$B (\alpha u + \beta v) = \lambda (\alpha u + \beta v).$$ If we insert the definition of $B$, we discover that this is equivalent to the 2-by-2 eigenvalue problem given by $$ C \begin{bmatrix} \alpha \\ \beta \end{bmatrix} := \begin{bmatrix} v^T u & v^T v \\ u^T u & u^T v \end{bmatrix} \begin{bmatrix} \alpha \\\beta \end{bmatrix} = \lambda \begin{bmatrix} \alpha \\ \beta \end{bmatrix}.$$ We observe that the matrix $C$ is not necessarily symmetric, so in order to eliminate the possibility of complex eigenpairs, we need to study $C$ further. The determinant is $$\text{det} (C) = (u^Tv)^2 - \|u\|_2^2 \|v\|_2^2.$$ By the Cauchy-Schwartz inequality, we conclude that $\text{det}(C) \leq 0$. We can now conclude that the eigenvalues of $C$ cannot be complex! Why is that? The matrix $C$ is real and has dimension $2$, so if the eigenvalues are complex, then they are necessarily a pair of complex conjugate numbers, i.e., $$\lambda_1 = \lambda \quad \text{and} \quad \lambda_2 = \bar{\lambda}$$ for some $\lambda \in \mathbb{C}$. However, the determinant is also the product of the eigenvalues, i.e., $$\text{det}(C) = \lambda_1 \lambda_2 = |\lambda|^2,$$ so $\text{C} \leq 0$ forces $\lambda = 0$ and the eigenvalues of $C$ must necessarily be real. Given the possible eigenvalues we can now compute the eigenvectors. This demonstrates that the problem can potentially be reduced to two rank-1 updates.
There remains the very practical problem of solving the 2-by-2 eigenvalue problem. You may experience overflow in extreme cases and if this is an issue, then this would be the basis for a good new question.
Be mindful that you do not equate speed with computational efficiency. The rank-1 update procedure for a Cholesky factorization is not efficient in the sense that all level-1 BLAS operations have low arithmetic intensity, i.e., few arithmetic operations per memory operation. Therefore, their performance is bounded by memory bandwidth/latency and they only run at a small fraction of the peak flop rate. While you obviously will save time and energy by doing 2 rather than 3 rank-1 updates, and the wall-time may be sufficiently short, the resource utilization will be rather poor. This has been an issue since the 1990s, because the peak flop rate has grown faster than the memory bandwidth.
To the best of my knowledge there is no algorithm for your problem that can be considered efficient by today's standards.
EDIT: The remark by @LutzLehmann simplifies the analysis considerably. We literally have $$ uv^T + vu^T = \frac{1}{2} \left[ (u+v)(u+v)^T - (u-v)(u-v)^T \right],$$ so there is no need to detour to the 2-by-2 eigenvalue problem.