Updating the eigensystem of modified symmetric matrices is an important task arising from certain fields of applications. The core of the problem is computing the eigenvalues and orthogonal eigenvectors of a diagonal matrix with symmetric low rank modifications, {\em i.e.} $D + UHU^T$. The eigenproblem of this type of matrix has long been studied since Golub et.al~\cite{pa} proposed theoretical considerations about it. Currently, there exist methods to compute eigenvalues accurately, but the difficulty remains in numerical stability and orthogonality of computed eigenvectors.
The main contribution of this thesis is a new method to compute all the eigenvalues and eigenvectors of a real diagonal matrix with a symmetric low rank perturbation. The algorithm computes an orthogonal matrix $Q = [q_1, q_2, \ldots, q_n]$ and a diagonal matrix $\Lambda = diag\{\lambda_1, \lambda_2, \ldots, \lambda_n\}$ such that $AQ = Q\Lambda$. Here the matrix $A = D + UHU^T$ has the special structure that $D \in\mathbb{R}^{n\times n}$ is a diagonal matrix, $U\in\mathbb{R}^{n\times r}$ is a column orthorgonal matrix and $H\in\mathbb{R}^{r\times r}$ is a symmetric matrix. $n$ is the dimension of $A$ and $r<
Aside from solving the eigensystem update problem mentioned above, our proposed method can also be used in the divide and conquer eigenvalue algorithm. Cuppen's divide and conquer algorithm~\cite{cuppen} solves a rank-one update of eigensystem in its merge step for a symmetric tri-diagonal matrix. A symmetric banded matrix will require solving the eigensystem of a low-rank perturbed diagonal matrix, {\em i.e.}, $D+UHU^T$. Efficient solution to this problem in the merge step can potentially enable application of divide and conquer algorithm directly on symmetric banded matrix.
In our proposed algorithm, eigenpairs are mostly computed by Rayleigh Quotient Iteration safe-guarded with bisection, with each eigenpair requiring $O(nr^2)$ flops to compute. Hence the overall computational complexity for our algorithm is $O(n^2r^2)$. This is an appealing quadratic algorithm since $r$ is usually considered a small constant relative to $n$. To ensure numerical stability, eigenvectors corresponding to eigenvalue clusters are computed through a special \textit{orthogonal deflation} method that completely avoids re-orthogonalization. In case of tight clusters, extended precision arithmetic is used for eigenvectors corresponding to close eigenvalues.
We present both theoretical analysis and numerical results to support our claim that the proposed algorithm is numerically stable.
