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.


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.

  • $\begingroup$ Interesting! Do you have any reference of how this method works? $\endgroup$ – mathreadler Sep 18 '15 at 15:59
  • $\begingroup$ 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! $\endgroup$ – 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.

  • 1
    $\begingroup$ 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! $\endgroup$ – electrique Sep 18 '15 at 17:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.