[Math] Solving a complex, sparse, linear system using the Schur complement

linear algebramatricesnumerical linear algebraschur-complementsparse matrices

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.

Best Answer

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.

Related Question