# What is a good way to solve the following linear system? (repeatedly)

Let $n,m\in \mathbb{N}$ be such that $m\ge n$.

Let $M_1\in \mathbb{R}^{n\times n}$, $\{M_{2},M_{3} \} \subset \mathbb{R}^{m\times m}$ be symmetric positive definite and computationally cheap to invert. You can think of them as a diagonal matrix. Let $S\in \mathbb{R}^{n\times m}$ be full rank. Let $B\in \mathbb{R}^m$. I am interested in solving the following linear system over a range of $\omega \in \mathbb{R}$. \begin{align*} \left\lbrack \begin{array}{cc} -i \omega M_1& S\\ -S^\mathrm{T} & -i\omega M_2+ M_3 \end{array} \right\rbrack \left\lbrack \begin{array}{c} U\\Q \end{array} \right\rbrack = \left\lbrack \begin{array}{c} 0\\B \end{array} \right\rbrack. \end{align*}

If $M_1$ and $M_2$ were both $0$ this would be a saddle point problem and if $M_3$ was zero this would be a Helmholtz problem. I am interested in solving this equation over a range of values of $\omega$, so reusing decompositions might be useful. This system comes from a Telegrapher's equation.

I am writing to inquire if there is a body of research on these types of problem. I know the research for the saddle point problems and Helmholtz problem is pretty extensive. Any ideas or suggestions are appreciated.

Edit.

This is what I am planning on doing, without getting a better idea from here.

Writing $U= U_R + i U_I$ and $Q= Q_R + i Q_I$ and rewriting the system by equating real and imaginary parts we get

$$\left\lbrack \begin{array}{cccc} -\omega M_1 & 0 & 0 & S\\ 0 & \omega M_1& S&0\\ 0& -S^\mathrm{T}& -\omega M_2 & M_3\\ -S^{\mathrm{T}}&0&M_3& \omega M_2 \end{array} \right\rbrack \left\lbrack \begin{array}{c} U_R\\ U_I\\ Q_R\\ Q_I \end{array} \right\rbrack = \left\lbrack \begin{array}{c} 0\\ 0\\ 0\\ B \end{array} \right\rbrack.$$

Eliminating $U_I$ and $U_R$ from these we get

$$\left\lbrack \begin{array}{cc} \omega M_3 & S^\mathrm{T} M_1^{-1} S - \omega^2 M_2\\ S^\mathrm{T} M_1^{-1} S - \omega^2 M_2& \omega M_3 \end{array} \right\rbrack \left\lbrack \begin{array}{c} Q_R\\ Q_I \end{array} \right\rbrack = \left\lbrack \begin{array}{c} \omega B\\0 \end{array} \right\rbrack.$$

Now we have a standard saddle point problem. I plan on solving it using a Schur complement (inner/outer CG solves) method. Note that I will use the $\omega M_3$ matrices as "pivots" because there exists $\omega$ for which $S^\mathrm{T} M_1^{-1} S - \omega^2 M_2$ is singular. Having obtained $Q_I$ and $Q_R$ we can use the eliminated equations to recover $U_I$ and $U_R$.

What I don't like about this approach is that there is no reuse. I mean the method doesn't achieve any economies of scale even though it will solve many closely related problems.

• How big are $n$ and $m$ likely to be? Have you tried a block factorization? That should let you reuse factorizations of all of the submatrices, though you will have to refactor $i\omega M_2+M_3$ each time you change $\omega$ unless there's some very special structure in $M_2$ and $M_3$. – Bill Barth Oct 1 '15 at 21:14
• The matrices will be as can big as can fit on the computer's RAM. (normal workstation not a super computer). The $M$ matrices are actually very sparse and almost diagonal so inverting them for a given RHS is a trivial operation. $m$ will be $3$ times $n$ for $3$- dimensional problems and and $2$ times $n$ for $2$ dimensional problems. As for as factorizations, I really only think it would be worthwhile to to factor $S$, maybe into some sort of sparse SVD. – fred Oct 1 '15 at 21:49
• – Bill Barth Oct 2 '15 at 2:52
• "Do block-wise Gaussian elimination" -- which is also known as forming the "Schur complement" in this context. – Wolfgang Bangerth Oct 2 '15 at 2:53
• @BillBarth So nothing glaringly obviously wrong or inappropriate as far as you see it? Thanks. – fred Oct 2 '15 at 14:33

Sorry this is late, but I've been thinking about an iterative solution of this which might be an alternative. So you're solving $(\omega A+B)x = b$ for $x$ with $A$ easily invertible. Can you make an iterative version as $$\omega A x^{(i+1)} + Bx^{(i)} = b,$$ or equivalently $$x^{(i+1)} = \omega^{-1}A^{-1}b - \omega^{-1}A^{-1}Bx^{(i)}.$$ You can precalculate $A^{-1}b$ and $A^{-1}B$, and as @BillBarth said, you could try the $x$ for the previous $\omega$ as your starting guess. Maybe I'm missing something.
• Hi, I think this is valid. I will check and try to get an idea of the spectrum of the $A^{-1}B$ operator to see how fast this will converge. It looks like for large $\omega$ it will converge fast. – fred Oct 5 '15 at 14:03
• Also, since you're only doing the inversion once, could check out $$x^{(i+1)} = B^{-1}b - \omega B^{-1}Ax^{(i)}$$ – Steve Oct 5 '15 at 14:36
• This will only work for $\omega$ large enough. For your first example, we need that $\left\|\omega^{-1} A^{-1}B \right\| <1$. This is a pretty drastic restriction on $\omega$ in my application since $\left\| A^{-1}B \right\| \to \infty$ as the mesh gets smaller. There might be a way to fix this though. Thanks for the idea. – fred Oct 5 '15 at 15:04