Finding null vectors of a parameter-dependent matrix

Question

I have dense complex matrices $M(z)$ in which each element $M_{ij} = M_{ij}(z)$ depends on a complex parameter $z$. I need to find $z$ such that the matrix $M$ gets singular, i.e. I am looking for null vectors $\vec v$ which satisfy $M(z) \vec k = \vec 0$.

So far, I use a Newton method. As the wanted $M(z)$ is singular, I need $M$'s determinant to vanish. Starting with an initial guess $z_0$, I iterate $$z_{i+1} = z_i - \frac{g(z_i)}{g'(z_i)}$$ with $g(z) := \det M(z)$. To avoid computing the determinant from definition, I use LU factorization. The iteration is stopped when $|z_{i+1} - z_{i}|$ is sufficiently small.

Using the matrix identity $\frac{\mathrm{d}}{\mathrm{d}z}\ln \det M(z) = \mathrm{trace}\,(M^{-1}(z) M'(z))$, the reciprocal logarithmic derivative in the iteration formula can replaced to yield $$z_{i+1} = z_i - \frac{1}{\mathrm{trace}\,(M^{-1}(z) M'(z))}\,,$$ which is what I use in the end.

A few details on steps involved: so far, I compute all derivatives $g'(z)$ from a forward finite difference scheme $g'(z) = \frac{g(z + h) - g(z)}{h} - \mathrm{i}\frac{g(z + \mathrm{i}h) - g(z)}{h}$, ($h$ real), likewise for derivatives of $M$. Also note I'm assuming that singular vectors have multiplicity $m = 1$ (although adaptation of the formulas above is possible using $g=(\det M(z))^m$.

The matrix $M$ has no special structure in global (i.e. not Hermitian). Depending on the problem, it may be rectangular or square. In any case, it is dense and well-conditioned. Typical sizes range from ~ $100\times100$ to $10k \times 10k$. $M$'s origin is in the boundary element method. For open systems (that may have resonances), the null vectors I am looking are the resonance wavefunctions on the boundaries of the domains under consideration.

I'd like to learn about alternative methods of adjusting the parameter $z$ such that $M$ gets singular. Do you have any other ideas? Or any comments on the method I described? Although the results are fine, I don't like the approach too much as both iterations described here involve a step that is essentially $O(N^3)$ (LU or the inverse). Furthermore, for the second scheme, I need to compute the inverse of near singular matrices.

Thank you for your input!

Answer 1

I have used the method you described for a similar problem. I found that a combination of your Newton iteration, along with Brent's method on minimizing the smallest singular value gave pretty good results. I believe I used Newton's method first for a few iterations, then switched over to Brent's. This paper of mine contains the description for the analogous iteration on the SVD.

Your problem is generally called a "nonlinear eigenvalue problem". Applying Newton's method as you have done is one standard way of solving it. The literature desribes a lot of methods for it (for example). Besides these kinds of iterations to find singular eigenvalues, you can also perform contour integration around the region in the complex plane where you want to compute eigenvalues. There is also a fairly extensive set of literature on this (see this or this). These links are more or less the first things I found on Google, so please don't think that they are the best examples out there.

Edit: Note that your formula for the finite difference approximation of $g'$ does not need the imaginary correction (just the first term is enough) if your function is analytic. If you can compute the derivatives analytically or automatically (with automatic differentiation), then that would be much better. However, I would guess this is not possible if your matrix comes from a BEM discretization.

Answer 2

The problem of finding singular points of a parameter-dependent matrix is called the nonlinear eigenvalue problem in the field of numerical linear algebra. I think you could benefit from using developed methods. It is a quite active area of research. You can consider solving it with the Julia software package NEP-PACK (I am a lead developer of NEP-PACK). NEP-PACK contains an implementation of the method of Beyn cited in the previous post. There are some methods also in SLEPc.