Inverse of high condition number matrices

Question

I am working on estimating the parameters of a linear system which is as follows:

$Ax = b$, $A \in \mathbb{R}^{n \times k} $, $x \in \mathbb{R}^{k \times 1} $ and $ b \in \mathbb{R}^{n \times 1} $

where $A$ and $b$ are given and $x$ has to be computed. The problem is: $A$ has a very high condition number ($~10e6 - 10e7$) which is affecting my results of $x$ even with a small noise($~10e-3 - 10e-2$) in $b$.

I have come across literature where such problem was tackled. They have broken down the data matrix($A$) column-wise into multiple sub-matrices such that the sub-matrices are well-conditioned. Please refer to pages $3-4$ and equations $9-15$.

$ A = [A_{1} \quad A_{2} \quad ... \quad A_{l}$]

The parameter vector x is also broken down accordingly which gives

$\Sigma_{i =1}^{l} A_{i}x_{i} = b$

Now, the author solves for the following set of equations iteratively for $i = 1: l$ to compute the the value of x

$x_{i} = A_{i}^{+}b_{i}$

where ${\{.\}}^+$ denotes the pseudo-inverse

My doubt is: How to calculate $b_{i}$. We only know $b$, but how do we compute the component of $b$ contributed by $x_{i}$

Answer 1

If $$ A=[A_1,\ldots,A_\ell] \quad \text{and} \quad x=\begin{bmatrix}x_1\\\vdots\\ x_\ell\end{bmatrix} $$ are conforming partitions of $A$ and $x$ in $Ax=b$, then from $$ Ax=b \quad \Leftrightarrow \quad \sum_{i=1}^\ell A_ix_i=b $$ we have $$\tag{1} x_i=A_i^+b_i, \quad b_i=b-\sum_{j=1\\j\neq i}^\ell A_jx_j, $$ which is a recipe for a fixed point iteration.

Let $x^{(k)}$ be an approximation of $x$ at step $k$. A simple way to implement an iteration on (1) is for $k=1,2,\ldots$, do $$\tag{2} x_i^{(k+1)}=A_i^+\left(b-\sum_{j=1\\j\neq i}^\ell A_jx_j^{(k)}\right), \quad i=1,\ldots,\ell. $$

Since at step $i$ of the inner loop you already have values of $x_1^{(k+1)},\ldots,x_{i-1}^{(k+1)}$, you can also improve (2) by using the already computed components on the right-hand side as $$\tag{3} x_i^{(k+1)}=A_i^+\left(b-\sum_{j=1\\j<i}^\ell A_jx_j^{(k+1)}-\sum_{j=1\\j>i}^\ell A_jx_j^{(k)}\right), \quad i=1,\ldots,\ell. $$

Note that the step from (2) to (3) is similar to how Gauss-Seidel is derived from Jacobi.

I don't know how this method is called or what are its convergence properties. It reminds me a bit of the Cimmino's method which works by splitting the rows instead of columns.