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}$


1 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.

  • $\begingroup$ Used it on a regressor matrix with approx. 550 condition number and with slight noise in the data (~ 1e-4). I get the same result as $A^{+}b$ by using this method. $\endgroup$
    – bluebird
    Sep 16, 2017 at 10:27
  • $\begingroup$ @bluebird That is good or bad? :) I would not expect it to converge to anything else. In addition, condition number of 550 is not really high. $\endgroup$ Sep 18, 2017 at 5:15
  • $\begingroup$ Since the paper used an iterative method by computing the inverse of column submatrices rather than $A^{+}b$, I thought the earlier would yield a better result. Your method seems to be very close to the author's description. Saing that, I would like to ask if you have anything else in mind regarding the implementation in the paper? Also, thanks for informing that such a condition number is not really high. Even after trajectory optimization, I could land only there. Probably, I will have to look into preconditioning methods. $\endgroup$
    – bluebird
    Sep 19, 2017 at 7:17

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .