8
$\begingroup$

I am interested in solving a sequence of shifted linear systems $(A+\sigma I)x = b$ for various values of $\sigma$. The matrix $A$ is sparse and not too large, so I have its LU factorization available. What is the best way to do this?

This is for calculating a frequency spectrum; $\sigma$ is roughly the frequency so it will vary smoothly in an interval. I expect that solutions to nearby $\sigma$ will be somewhat similar. My first thought was to use an iterative solver like GMRES with the LU factors as a preconditioner, but I'm wondering if there is any better method.

Improve this question
$\endgroup$
6
  • 1
    $\begingroup$ This question provides several approaches. I recommend that these questions be merged. $\endgroup$
    – Jed Brown
    May 31, 2012 at 0:19
  • $\begingroup$ I'm not sure that a simple merge would work well here. The top two answers for the question you linked to exploit positive semidefinite structure of $A$. Only your answer applies to the scenario in this question. $\endgroup$
    – Geoff Oxberry
    May 31, 2012 at 2:06
  • $\begingroup$ Fair, I withdraw my suggestion. $\endgroup$
    – Jed Brown
    May 31, 2012 at 3:24
  • $\begingroup$ @GeoffOxberry: Actually, nothing about my answer requires $A$ to be positive semi-definite. VictorLiu, Could you clarify what you mean by calculating the "frequency spectrum"? Are you seeking the Schur decomposition of $A$? $\endgroup$
    – Jack Poulson
    May 31, 2012 at 5:12
  • $\begingroup$ @JackPoulson: I apologize; I stopped reading your answer as soon as I got to the first mention of positive semidefiniteness. You're absolutely right; your answer does not require symmetry. $\endgroup$
    – Geoff Oxberry
    May 31, 2012 at 5:31

3 Answers 3

Reset to default
8
$\begingroup$

There is a method called Automated MultiLevel Substructuring (AMLS) which was originally designed for a similar problem in vibration analysis, where the solution of the linear system with a particular shift corresponds to the frequency response problem at a frequency which is the square-root of the shift. The basic idea is to use nested dissection in order to generate a tree of separators and substructures, and to use extensions of the low-frequency modes of substructures (which can be found cheaply) as a means of reducing the global model. It has a setup cost proportional to a single factorization, but the cost grows with the highest frequency that you would like to solve to.

The original paper is here, and an extension to Maxwell's equations is here.

Disclaimer: I used to work in the lab that invented AMLS

Improve this answer
$\endgroup$
1
  • $\begingroup$ The link to the paper is dead. $\endgroup$
    – asmeurer
    Mar 3, 2017 at 0:28
5
$\begingroup$

The solution $x$ (and hence any simple response computed from it) is an analytic function of $\sigma$, hence it can be well approximated locally by rational functions. The response has poles at the eigenvalues of $-A$, which is why rational interpolation is needed. The following interpolation technique should work quite well with a number of factorizations approximately equal to the number of peaks in the responses.

Compute factorizations for three values of $\sigma$, namely the endpoints and the midpoint of your interval, and then calculate the responses of interest and their first two (or more) derivatives with respect to $\sigma$. Note that given the factorization, you can easily get the derivatives by multiple backsolves.

Then use piecewise rational Hermite interpolation, first a single rational, then two rational pieces, one in each subinterval. Assuming that the responses of interest decay to zero for $\sigma\to\infty$ (if not one needs to change the degrees in the interpolant), $d$ derivatives at the endpoints of a subinterval allow one to interpolate by a quotient of a degree $d-1$ polynomial divided by a degree $d$ polynomial. For $d=2$, this captures a single Lorentz peak exactly, for $d=4$ it captures a superposition of two. (If this is not enough, one can also consider multipoint Hermite interpolation, making better use of the computed data, at the expense of more programming work.)

Test the accuracy of the current interpolant by an additional exact calculation at the point where the last two interpolants differed most, using a new factorization and a new rational fit in the two intervals created by the new point, to decide when you can quit.

Unless your spectrum is very complicated, you'll need few factorizations only.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Interesting approach, but I would tend to classify the spectrum of indefinite Maxwell (particularly at high frequency) as "very complicated" since several eigenvalues will likely cross the origin as you move between two values of $\sigma$. If you have many sample points lying between eigenvalues, I would expect this method to work well. Have you had different experience? $\endgroup$
    – Jed Brown
    May 31, 2012 at 15:23
  • $\begingroup$ @JedBrown: I updated the answer a bit. Only that part of the spectral density that gives rise to peaks in the responses is relevant for the complexity. I don't believe that there is any method that will work well with significantly fewer factorizations (or full iterative linear solves) than the number of peaks to be resolved. $\endgroup$
    – Arnold Neumaier
    May 31, 2012 at 15:48
3
$\begingroup$

You can use the direct solve with $A + \sigma_0 I$ as a preconditioner for $A + \sigma_i I$ as long as $\sigma_i$ is reasonably close to $\sigma_0$. This tends to work well if the spectrum is away from the origin. You can improve the approximation somewhat using the first order correction described in this answer. If the spectrum crosses the origin between the two shifts $\sigma_0$ and $\sigma_i$, as occurs for Helmholtz, I do not know a way to utilize the old factorization.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Unfortunately, the spectrum of A is already indefinite. You could say that A is the Laplacian in the Helmholtz equation, but possibly with negative weights, so it's not even positive semidefinite. $\endgroup$
    – Victor Liu
    May 31, 2012 at 9:44
  • $\begingroup$ As I suspected. Of course it is not surprising that a solve with one shift is not a good preconditioner for another since the solutions of the two problems may be very different (the difference between the two waves gives you beats because they have slightly different wavelength). $\endgroup$
    – Jed Brown
    May 31, 2012 at 12:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.