# Solving a complex, sparse, linear system using the Schur complement

## Solution method

I am repetitively solving sparse linear systems (for the need of ARNOLDI iterations) of the type: $$\underbrace{\begin{bmatrix} J_1 & J_2 \\ J_3 & J_4 \end{bmatrix}}_J \underbrace{\begin{bmatrix} x \\ y \end{bmatrix}}_z= \underbrace{\begin{bmatrix} b_x \\ b_y \end{bmatrix}}_b$$

I build the Schur-complement ($J_4$ is non-singular): $$J_D=J_1-J_2J_4^{-1}J_3$$ Then, I solve for $x$: $$J_Dx=b_x-J_2J_4^{-1}b_y$$ and then for $y$: $$J_4y=b_y-J_3x$$ At this point I have a full solution of $z$.

## Intermediate calculations

To calculate $J_2J_4^{-1}J_3$, I first compute $\tilde{J_2}=J_2J_4^{-1}$ as follows: $$\tilde{J_2}^T=(J_2J_4^{-1})^T=(J_4^{-1})^TJ_2^T$$ Thus, I get $\tilde{J_2}$ by solving the system: $$J_4^T\tilde{J_2}^T=J_2^T$$ and tacking the transpose.

Edit: The above is implemented in Matlab and I use the conjugate transpose ' (mathworks.com/help/matlab/ref/ctranspose.html) for both the real and the complex. Also I use dot (mathworks.com/help/matlab/ref/dot.html) for the dot products.

## Problem

This way of solving works perfectly when $J$ is real (even if $b$ is complex). That is, solving $Jz=b$ directly or by building the Schur-complement gives the same $z$ as expected. However, when $J_4$ becomes complex, then I get different results... I suspect that I did something really wrong when taking the transpose for $\tilde{J_2}$ or doing the dot product for $J_2J_4^{-1}b_y$. Yet, I fail to detect my error.

Note: The reason I solve using the Schur-complement is that matrices $J_2,J_3,J_4$ have some particular structures that allow for good parallelization of the operations involved. E.g. $J_4$ is block diagonal and non-singular.

• Interesting! Do you have any reference of how this method works? – mathreadler Sep 18 '15 at 15:59
• Hi, there is nothing really complicated, just linear algebra manipulations. At the wiki article (en.wikipedia.org/wiki/Schur_complement) there is a section "Application to solving linear equations". If you need more help, let me know! – electrique Sep 18 '15 at 17:32

Your approach looks good, so you may have a bug somewhere in your implementation. As a validation, I ran the following simple test case involving a randomly generated dense complex matrix $J$ and a randomly generated real right-hand side $b$.

N = 100;
M = 30;
J = exp(2i*pi*rand(N));
b = rand(N,1);
Jd = (J(M+1:end,M+1:end)' \ J(1:M,M+1:end)')';
x = b(1:M) - Jd*b(M+1:end);
Jd = J(1:M,1:M) - Jd*J(M+1:end,1:M);
x = Jd \ x;
y = b(M+1:end) - J(M+1:end,1:M)*x;
y = J(M+1:end,M+1:end) \ y;
z = [x; y];
ztrue = J \ b;
err = norm(ztrue - z) / norm(ztrue);
res = norm(J*z - b) / norm(b);
restrue = norm(J*ztrue - b) / norm(b);


Running this code multiple times the approximate averaged values for $\texttt{err}$, $\texttt{res}$, and $\texttt{restrue}$ were $1\times 10^{-13}$, $1\times 10^{-13}$, and $4\times 10^{-14}$, respectively.

• Thank you! Indeed it was a bug and your program showed me the problem. In several places of the algorithm, I was using the dot product instead of *. Changing the latter fixed my problem! – electrique Sep 18 '15 at 17:15