Numerical continuation of all eigenvalues of small, dense matrices

Question

Consider a one-parameter family of matrices $A(q)\in \mathbb{R}^{n\times n}$, $q\in\mathbb{R}$. For my applications, $n$ is typically between $5$ and $50$ and $A(q)$ is generally dense, so direct methods for obtaining the entire spectrum of $A(q)$ are applicable and not too expensive. I would like to analyze the eigenvalues of $A(q)$ as I vary $q$, but I do not know which eigenvalues are interesting (changes sign, multiplicity, real to complex conjugate pairs, etc.) a priori, so I need to see how all the eigenvalues of $A(q)$ vary.

The problem I am running into is that in general, the eigenvalues returned by some eigensolver will not be ordered the same, so it is difficult to make a correspondence between $\sigma(A(q_1))$ and $\sigma(A(q_2))$, even when $q_1$ and $q_2$ are close. I know that there will be issues when multiplicity changes, but even for simpler cases this is still an issue.

I have come up with two methods that seem to work on some test problems, but they seem either inefficient or hacky.

1. Compute a single eigenvalue at one point $\lambda_1(q_0)$ and numerically continue it in $q$; repeat for all eigenvalues of initial matrix $\lambda_i(q_0)$. This will certainly get be the eigenvalue loci I want, but it requires potentially many sequential eigenvalue decompositions. Additionally, I would prefer not to use sophisticated continuation techniques and treat the eigensolver purely as a black box, as I eventually would like to extend this to a parameterized nonlinear eigenvalue problem.

2. Compute all eigenvalues for a grid of $q$: $\sigma(A(q_0))$, $\sigma(A(q_1))$, $\sigma(A(q_3))$, etc., where $q_{i+1} = q_i+\Delta q$. Assign each of these to the appropriate locus via nearest neighbors between successive sets of eigenvalues. That is, assign $$\lambda_k(q_{i+1}) = \mathrm{argmin}_{\mu\in\sigma(A(q_{i+1}))}\{|\mu - \lambda_k(q_i)|\}.$$ This feels hacky and will almost certainly get me into trouble around changing multiplicities when the loci cross, but it allows for all of the eigendecompositions to be computed fully in parallel.

Are there any established methods for performing this kind of computation efficiently?

Answer 1

I think a good direction would be to look at continuation of eigenvectors/invariant subspace, not just the eigenvalues. You could proceed as follows:

• calculate schur form: $A(q_i) = Q(q_i)T(q_i)Q(q_i)^*$
• apply Q(q_i) to new matrix: $Q(q_i)^*A(q_{i+1})Q(q_i) $, this should give you a matrix that is almost (quasi-)upper-triangular, which you can then refine to give you the true Schur form. To do so, you could use iterative refinement as described by Bujanovic et al. [1].

[1] Bujanović, Zvonimir; Kressner, Daniel; Schröder, Christian, Iterative refinement of Schur decompositions, Numer. Algorithms 92, No. 1, 247-267 (2023). ZBL1507.65076.