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

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?

• Are the entries D all positive/negative, or are they mixed in sign? – rchilton1980 Sep 17 at 0:38

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.

• Isn't this basically the same thing as the Schur complement approach that OP wanted to get away from? $L_{22}U_{22}$ is the Schur complement. – Federico Poloni Sep 17 at 11:52
• Agreed with Federico, if the input $\mathbf H$ is ill conditioned, you might encounter cancellation problems when forming $\mathbf A - \mathbf B^T \mathbf D^{-1} \mathbf B$ – rchilton1980 Sep 17 at 13:01
• I actually did not know about the LU condition number trick! I have gotten the details from Golub and Van Loan. – Damien Sep 17 at 14:50
• I will have a look at some of the failure cases I have and see if there are any problems with the $L_{22} U_{22}$. – Damien Sep 17 at 14:53
• It seems OP is happy with this answer. But I'll mention a few avenues for improving things further. It may be possible to scale the original matrix to reduce the condition number before doing the LU factorization. See LAPACK's DGEEQU routine for an example. If the matrix is still close to singular, you could use a partial or full pivoting approach in the $L_{22} U_{22}$ step to help, and also reveal the rank of the matrix. – vibe Sep 17 at 15:23