Efficient solution to a structured symmetric linear system with condition number estimation

Question

I have a real-valued linear system $Hx = b$ where $H$ is symmetric matrix** (not necessarily positive/negative definite) with a very particular structure:

$$ H = \begin{bmatrix} D && B \\ B^T && A\end{bmatrix} $$

Where:

• $A$ is dense with a small fixed dimension (e.g. 3x3)
• $D$ is diagonal with a (relatively) large dimension that changes per-problem (e.g. 200 x 200)
• $B$ is dense with dimension (200 x 3)

(Note**: $A$, $B$ and $D$ are stored separately. These blocks can be re-arranged if it helps.)

The solution is currently via an explicit inverse using a Schur complement, which is problematic if H has a large conditions number (which arises from time to time).

I would like to compute the solution $x$ in an efficient way that exploits this peculiar problem structure, and provides an estimate of the condition number to evaluate whether we should trust the solution.

My thought was to tridiagonalize away $B$ using Householder transformations, which will make it easy to compute the eigenvalues and hence the condition number. The idea is that since $D$ is diagonal, there should be much less work involved.

Would there be any better approaches? For example, some other structure that exploits the large block diagonal component?

Answer 1

If we assume that $D$ is nonsingular, then there is a relatively straightforward (and efficient) solution based on an $LU$ decomposition. If we write $$ \pmatrix{D & B \\ B^T & A} = \pmatrix{ L_{11} & \\ L_{21} & L_{22}} \pmatrix{U_{11} & U_{12} \\ & U_{22}} = \pmatrix{L_{11} U_{11} & L_{11} U_{12} \\ L_{21} U_{11} & L_{21} U_{12} + L_{22} U_{22}} $$ where we have chosen the same partitioning of the $L$ and $U$ matrices, then we have the following four subproblems:

(1) $D = L_{11} U_{11} \rightarrow L_{11} = I, U_{11} = D$

(2) $B = L_{11} U_{12} \rightarrow U_{12} = B$

(3) $B^T = L_{21} U_{11} \rightarrow L_{21} = B^T D^{-1}$

(4) $A = L_{21} U_{12} + L_{22} U_{22} \rightarrow A - B^T D^{-1} B = L_{22} U_{22}$

So the only real effort here is to solve a 3x3 LU decomposition problem, $$ A - B^T D^{-1} B = L_{22} U_{22} $$ which can be done with any standard library. Once you have all the $L$ and $U$ factors, you can easily solve the linear system with backward/forward substitution. There also exist standard library routines to compute the condition number of a matrix in $LU$ form - see for example the LAPACK DGECON routine.

EDIT: the backward/forward substitution step can (and should) also be optimized for this problem. Once we have $L$ and $U$, we need to solve two problems, \begin{align} Lz &= b \\ Ux &= z \end{align} I will examine the first equation and leave the second for you to work out. We have $$ \pmatrix{I & \\ B^T D^{-1} & L_{22}} \pmatrix{z_1 \\ z_2} = \pmatrix{b_1 \\ b_2} $$ So we immediately see $z_1 = b_1$ and $$ L_{22} z_2 = b_2 - B^T D^{-1} b_1 $$ This equation can be solved with the TRSV BLAS call.