Advice regarding solving ill-conditioned systems

condition numberlinear algebranumerical linear algebrasystems of equations

I start from a well-conditioned linear system $AX=B$ where $A$, $X$ and $B$ are complex square matrices of the same shape ($M$,$M$) and cond(A)=1

Then I need to solve linear systems that derive from the one above, where I replace some columns of A and B as follows: given a bit string $b$ of length $M$, if b[i] is 1 I replace the $i$-th column of $A$ with minus the $i$-th column of the identity matrix and if it's 0 I replace the i-th column of B with the same.

Example for $M=3$ (here the entries of A and B are integers for illustration):

$$
A = \begin{pmatrix}
1&2&3\\
4&5&6\\
7&8&9\\
\end{pmatrix}
\qquad
B = \begin{pmatrix}
1&2&3\\
4&5&6\\
7&8&9\\
\end{pmatrix}
\qquad
b = [0,0,1]
$$

Using this bit string new system is $AX=B$ with
$$
A = \begin{pmatrix}
1&2&0\\
4&5&0\\
7&8&-1\\
\end{pmatrix}
\qquad
B = \begin{pmatrix}
-1&0&3\\
0&-1&6\\
0&0&9\\
\end{pmatrix}
$$

The issue is that almost always the derived system is badly conditioned (cond(A) is anywhere between 1 and 100 in my tests) and eventually the algorithm I'm developing which relies on solving for X for all bit strings gives incorrect results.

I've tried solving with QR decomposition, but it didn't help much. Given the particular nature of the derived system, what else could I try?

Best Answer

I'll consider the case where b is a sequence of 1's followed by 0's so that $A$ has its first columns removed, and $B$ has its last columns removed; the approach can be generalized from there. The resulting matrices can be written as $$ A_1 = \pmatrix{-I & A_{12}\\0 & A_{22}}, \quad B_1 = \pmatrix{B_{11} & 0\\B_{21} & -I}. $$ We can multiply $A_1$ by a matrix from the right to get a new system: $$ \left[\pmatrix{-I & A_{12}\\0 & A_{22}} \pmatrix{I & A_{12}\\0 & I}\right] \left[\pmatrix{I & A_{12}\\0 & I}^{-1}X\right] = \pmatrix{B_{11} & 0\\B_{21} & -I} \implies\\ \underbrace{\pmatrix{-I & 0\\0 & A_{22}}}_{A_2} \underbrace{\left[\pmatrix{I & A_{12}\\0 & I}^{-1}X\right]}_{X_2} = \pmatrix{B_{11} & 0\\B_{21} & -I} $$ From there, the matrix $A_2$ should be relatively well conditioned (cf. this post for instance), and the solution process therefore more straightforward. However, we can actually use the structure of $A_2$ to make the solution process easier. Breaking $X_2$ up vertically, we have $$ \pmatrix{-I & 0\\0 & A_{22}} \pmatrix{X_{21}\\X_{22}} = \pmatrix{B_{11} & 0\\ B_{21} & -I} \implies\\ \pmatrix{-X_{21}\\ A_{22}X_{22}} = \pmatrix{B_{11} & 0\\B_{21} & -I}\implies\\ \begin{cases} X_{21} = -\pmatrix{B_{11} & 0}\\ A_{22}X_{22} = \pmatrix{B_{21} & -I}, \end{cases} $$ which leaves you with only one, smaller system to solve in order to build $X_2$. From there, we have $$ X_2 = \pmatrix{I & A_{12}\\0 & I}^{-1}X \implies\\ X = \pmatrix{I & A_{12}\\0 & I} X_2 = \pmatrix{X_{21} + A_{12}X_{22}\\ X_{22}}. $$


Here's how we can generalize this approach. Let $i_1<i_2<\cdots<i_k$ denote the indices $i$ for which $b[i] = 1$, and let $j_1<j_2<\cdots<j_{n-k}$ denote the rest. Let $P$ denote the matrix whose columns are the columns of the identity matrix, but in the order $i_1,i_2,\dots,i_k,j_1,j_2,\dots,j_{n-k}$.

I claim that the matrix $P$ is a permutation matrix (so that $P^{-1} = P^T$) designed so that $AP$ and $BP$ are of the correct form for the above approach. With that, we can solve the system $$ (AP)\underbrace{(P^{-1}XP)}_{X_3} = BP. $$ for $X_3$, then recover $X = PX_3P^T$.

Related Question