Sparse, underdetermined system of linear equations

Question

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.

Answer 1

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.

Answer 2

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.

Answer 3

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:

http://www.math.vanderbilt.edu/~foucart/FL08.pdf http://ftp.math.uga.edu/~mjlai/papers/surveySS.pdf

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