Solving a singular system of linear equations with smallest solution?

Question

I have a problem where I am solving a system of linear equations but sometimes the system results in a singular matrix which cannot be easily solved. In this case I would like that those rows for which system is ill-defined would be computed as 0. How can I assure that?

An example matrix:

$$ \left[ \begin{array}{ccccc} 1 & 0 & -1 & 0 \\ 0 & 1 & -0.3 & -0.7 \\ -1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1\\ \end{array} \right] \boldsymbol{x} = \left[ \begin{array}{c} 0 \\ 0 \\ 0 \\ 1 \\ \end{array} \right] $$

I would like that for $\boldsymbol{x}$ the result is:

$$ \left[ \begin{array}{c} 0 \\ 0.7 \\ 0 \\ 1 \\ \end{array} \right] $$

I tried to use least-squares approach but it is not really necessary that I get a solution which has zeros for rows which are ill-defined because it is trying to minimize $|A\boldsymbol{x} - \boldsymbol{b}|$ and not $|\boldsymbol{x}|$ itself.

Some properties of the matrix above: there is always a diagonal with elements $1$, other elements are from $[-1,0]$, sum of each row is always or 1 or 0. The vector on the right can contain only 0 or 1 elements.

Answer 1

If you are willing to minimize the 2-norm of the solution rather than minimizing the number of non-zeros, you can use the pseudo-inverse. To solve A*x = b in MATLAB, that would be pinv(A)*b. In your example, this produces 2-norm of solution of 1.2120, vs. 2-norm of your preferred solution of 1.2207.

You can come closer to your objective by minimizing the one-norm of x, subject to A*x=b. That is a Linear Programming problem. In this case, it produces your preferred result, but is not guaranteed to minimize the number of non-zeros in general (unless something about the special properties of your A and b wind up doing so?).

What you really want is to minimize the 0-norm of x, subject to A*x=b. You can solve this as a Mixed Integer Linear Programming problem by using a big M approach to minimize the number of elements of x which are more than some numerical tolerance, say 1e-6, from 0. This is easy to do in YALMIP under MATLAB using YALMIP"s "implies", but is more computationally intensive than the previously mentioned approaches.

Answer 2

I think I found a solution. I think the following algorithm can do the trick.

1. Define a set of restricted indexes $M$.
2. Populate $M$ with indexes of all rows which have $1$ only on a diagonal (equal to: all rows for which the sum of its elements is $1$; or, for which all non-diagonal elements are $0$).
3. Remove all rows and columns corresponding to those indexes.
4. Add (original) index of every row for which the sum (of its elements in this reduced matrix) is $\ne 0$ to $M$.
5. If you added such indexes, go to step 3.
6. Otherwise you are left with rows which are making the matrix ill-defined.

In the original matrix we can now change all rows with indexes which are missing in $M$ such that for those rows only diagonal elements are 1 and the rest of elements 0. Then we solve the original problem.

For the above example we start with:

$$ \left[ \begin{array}{cccc} 1 & 0 & -1 & 0 \\ 0 & 1 & -0.3 & -0.7 \\ -1 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1\\ \end{array} \right] $$

Removing the last row and last column:

$$ \left[ \begin{array}{ccc} 1 & 0 & -1 \\ 0 & 1 & -0.3 \\ -1 & 0 & 1 \\ \end{array} \right] $$

Removing the second row and the second column:

$$ \left[ \begin{array}{cc} 1 & -1 \\ -1 & 1 \\ \end{array} \right] $$

The rest are resulting rows. We modify them in the original matrix:

$$ \left[ \begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & -0.3 & -0.7 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1\\ \end{array} \right] $$

With this we can go and normally solve the matrix equation.