# Sparse, underdetermined system of linear equations

I'm looking for an algorithm to solve the underdetermined system of linear equations $$\mathbf{A}\,\mathbf{x} = \mathbf{b}$$ with $\mathbf{A} \in \mathbb{R}^{n\times n}$, $\mathbf{b} \in \mathbb{R}^{n\times 1}$ and $rank(\mathbf{A}) < n$. Clearly, there exists an infinite number of solution vectors $\mathbf{x}$ (or no solution at all). However, if $\mathbf{A}$ is sparse, some components of $\mathbf{x} = [x_1,\ldots,x_n]$ might be determined uniquely. Example: $$\begin{bmatrix}2 & -1 & 0 & 0\\1&0&2&4\\1&0&2&4\\1 & 2 & 0 & 0\end{bmatrix} \,\begin{bmatrix}x_1 \\x_2 \\ x_3 \\ x_4\end{bmatrix} = \begin{bmatrix}0 \\ 1 \\ 1 \\ 5\end{bmatrix}$$ How do I find those components $x_i$ of $\mathbf{x}$ which have this property (in this example $x_1$ and $x_2$)? I'd be grateful for pointers to the relevant literature.

• How would you resolve cases like the one you have above, where $x_{1}$ and $x_{3}$ could be chosen as the components determined uniquely instead of $x_{1}$ and $x_{2}$? If such a situation occurs, do you want all possible sets of uniquely determined variables? – Geoff Oxberry Jan 9 '14 at 19:50
• @Geoff: If I'm not mistaken, the set of $x_i$, which can be determined uniquely, depend on the (given) values of $\mathbf{A}$ and $\mathbf{b}$; thus, there's just one set of uniquely determined $x_i$. – user71769 Jan 10 '14 at 11:05
• I think I misunderstood your question. Looking again at the example, you're right that $x_{1}$ and $x_{2}$ are uniquely determined. – Geoff Oxberry Jan 10 '14 at 19:35

The problem here is to determine where or not there are any components of $x$ that are uniquely determined by the linear system of equations.

Here's one approach: Compute a basis for the null space of $A$ by using a rank revealing QR factorization or by SVD or by any of a number of other methods. Let $U$ be the matrix of size $n$ by $p$ that contains as its columns the basis for $N(A)$. If row $i$ of $U$ is $0$, then $x_{i}$ is uniquely determined by the equations.

Unfortunately, this won't work very well in practice because of small round off errors in the determination of the basis for $N(A)$. For example, using MATLAB, I computed the null space for your example matrix:

A=[2 -1 0 0; 1 0 2 4; 1 0 2 4; 1 2 0 0]

A =

 2    -1     0     0
1     0     2     4
1     0     2     4
1     2     0     0


null(A)

ans =

        0


-5.5511e-17

8.9443e-01

-4.4721e-01

In this case, the null space is one dimensional. Looking at the basis vector for this null space, it's clear that the $x_{1}$ is uniquely determined by the equations. However, the component of the null space basis vector corresponding to $x_{2}$ is small but not quite exactly 0. Thus it isn't clear that $x_{2}$ is uniquely determined by the equations.

Another way of approaching this is to compute the RREF of the matrix $A$ in exact rational arithmetic and then obtain a basis for the null space of $A$. In this case, we get:

rref(A)

ans =

 1     0     0     0
0     1     0     0
0     0     1     2
0     0     0     0


Thus x=[0; 0; -2; 1] is a basis vector for the one dimensional $N(A)$. It's clear from this that $x_{1}$ and $x_{2}$ are uniquely determined by the equations. We were lucky in this case that the RREF computed by MATLAB produced an exact answer- for a much larger matrix we'd probably see some round-off error and then the use of higher (or infinite) precision arithmetic would be necessary.

Use a standard sparse reordering of the equations and variables. If all goes well, the matrix then decomposes into diagonal blocks (and possibly some remains) that can be easier examined for partial uniqueness of the solution.

The Scotch project provides "sequential and parallel sparse matrix block ordering"

The btf tool of SuiteSparse uses "maximal matching which maximizes the number of nonzeros on the diagonal" and "a method for finding the strongly connected components of a graph" to "give the permutation to block upper triangular form."

Google provides more links to papers discussing the theory and other stackexchange and overflow questions, like (https://mathoverflow.net/questions/68041/showing-block-diagonal-structure-of-matrix-by-reordering).

These algorithm interpret the sparsity structure of the matrix as a bipartite graph and construct edge weights from the entries. For instance log magnitude (bit length for integers) plus number of other row entries times number of other column entries. For this matrix it would give the weight matrix $$\begin{bmatrix} 4&1&-&-\\ 6&-&3&4\\ 6&-&3&4\\ 3&1&-&- \end{bmatrix}$$ where then the Hungarian algorithm or a variant is applied to find a max weight matching defining the reordered diagonal. Which can be seen here to be $(1,1),(2,3),(3,4),(4,2)$. Then semi-heuristic algorithms of clustering and partial graph ordering are applied to get as many entries as possible above the diagonal. Which would result here in moving the $(2:3, 3:4)$ block to the topleft diagonal position.

• Could you enhance your answer with a link to a description of these reordering methods? – Bill Barth Jan 9 '14 at 23:04
• Done. Also shortly explained what would happen with the example. – Lutz Lehmann Jan 10 '14 at 8:30

You might want to minimize $L_0$ norm $||x||_0$ s.t. $Ax=b$, as to find the the sparsest solution. Other $p$-norms also applicable depending on the case. Some resources to follow up are:

These methods should already cover the cases you like to solve.

• I don't think this helps. If you have the sparsest solution there may be some null-space vector that you can add to completely destroy the sparsity. Hence, the sparsest solution does not tell you which components are uniquely determined. You need to know the null space of the matrix (as @Brian Borchers explains). – Dirk Jan 9 '14 at 13:21
• At the begining of the question, he states: "I'm looking for an algorithm to solve the underdetermined system of linear equation" – Tolga Birdal Jan 9 '14 at 13:56
• As far as I see, the question is "How do I find those components $x_i$ of $x$ which have this property (in this example $x_1$ and $x_2$)?" Regarding just "find a solution": Why trying to find the sparsest solution (which is NP-hard in general) instead of calculating the minimum-norm solution (which is P (polynomial))? – Dirk Jan 9 '14 at 13:58