# 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$? Oct 14, 2015 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, 2015 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. Oct 15, 2015 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, 2015 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$). Oct 15, 2015 at 13:39