# What strategies / decompositions would be useful to solve the following linear system repeatedly if I only care about time to solution?

Let $A\in \mathbb{R}^{n\times n}$ be symmetric positive semidefinite, and $B\in \mathbb{R}^{n\times n}$ be symmetric positive definite. Suppose $B$ is block diagonal so it is easy to invert. (We have access to the entries of $B^{-1}$ in the computer). Let $\omega \in \mathbb{R}$. Suppose I want to solve the following symmetric and positive definite system $$((A-\omega^2 B)B^{-1}(A-\omega^2 B)+\omega^2 B)x = b,$$

for a range of values of $\omega$. Note that the RHS, $b$, may depend on $\omega$.

I am interested in strategies to compute $x(\omega)$ repeatedly. As an idealization, suppose I only care about correctness and the time it takes to produce a solution from a given $\omega$. I mean, any time spent performing decompositions and preprocessing will be considered free.

What decompositions / strategies would be useful?

# Edit

I have been playing around with this today. This is my proposed solution.

We use a combination of a Cholesky decomposition and an eigen decomposition. Let $C\in \mathbb{R}^{n\times n}$ be a lower triangular matrix such that $CC^\mathrm{T} = B$. Then

\begin{align*} ((A&-\omega^2 B)B^{-1}(A-\omega^2 B)+\omega^2 B)\\ &= \cdots \text{just algebra}\\ &= C\left((C^{-1} A C^{-\mathrm{T}} - \omega^2 I)(C^{-1} A C^{-\mathrm{T}} - \omega^2 I)+ \omega^2 I) \right)C^{\mathrm{T}} \end{align*}

Now we take an eigen decomposition of $C^{-1} A C^{-\mathrm{T}}$, i.e.

$C^{-1}A C^{-\mathrm{T}} = Q\Lambda Q^\mathrm{T}$ for $Q\in \mathrm{R}^{n\times n}$, orthonormal, $\Lambda\in \mathbb{R}^{n\times n}$, diagonal.

Then \begin{align*} C\big((C^{-1} &A C^{-\mathrm{T}} - \omega^2 I)(C^{-1} A C^{-\mathrm{T}} - \omega^2 I)+ \omega^2 I) \big)C^{\mathrm{T}} \\ &= \cdots \text{(just algebra)}\\ &= CQ \big((\Lambda - \omega^2 I)(\Lambda- \omega^2I)+ \omega^2 I\big) Q^\mathrm{T} C^\mathrm{T} \end{align*}

Now, this system can be trivially inverted, and none of the decompositions depends on $\omega$.

The real question is if any of this can be used to design a preconditioner for a system like this that is too large to be factored in this way. Any ideas on this would be much appreciated.

• How large is $n$? – Kirill Oct 14 '15 at 22:04
• For the purposes of the this question, we can say that $n$ is on the order of 100. Eigen decompositions, LU, cholesky are all feasible, while brute force Cramer's rule is not. I am mainly interested in solving this idealized problem as a way to build intuition to build a preconditioner. – fred Oct 15 '15 at 2:41
• I think you already got the best solution in your edit. The key decomposition here is that two symmetric matrices can be simultaneously diagonalized by congruence, as you have already figured out by yourself. – Federico Poloni Oct 15 '15 at 9:09
• @FredericoPoloni To reveal my ignorance, I don't know what it means to be diagonalized by congruence. Could you give me a link? – fred Oct 15 '15 at 13:13
• Given $A$ spd and $B$ symmetric, there is an invertible $M$ such that $M^*AM$ and $M^*BM$ are both diagonal (and actually a stronger statement holds, $M^*AM=I$). The operation $M^*AM$ is called congruence. See for instance Horn and Johnson, Matrix analysis, Corollary 4.6.12 (but you have basically already proved it yourself above, it is enough to take $M=CQ$). – Federico Poloni Oct 15 '15 at 13:39