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?


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.

  • $\begingroup$ How large is $n$? $\endgroup$
    – Kirill
    Commented Oct 14, 2015 at 22:04
  • $\begingroup$ 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. $\endgroup$
    – fred
    Commented Oct 15, 2015 at 2:41
  • 2
    $\begingroup$ 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. $\endgroup$ Commented Oct 15, 2015 at 9:09
  • $\begingroup$ @FredericoPoloni To reveal my ignorance, I don't know what it means to be diagonalized by congruence. Could you give me a link? $\endgroup$
    – fred
    Commented Oct 15, 2015 at 13:13
  • 2
    $\begingroup$ 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$). $\endgroup$ Commented Oct 15, 2015 at 13:39


Your Answer

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