# Solution to underdetermined linear equations

I have a set of numbers $x_i$ and I know sums of certain subsets $y_i=\sum x_{\sigma_k}$. All $x_i>0$ and I'm looking for a simple solution.

With some internet research I found that this might be related to problems in signal processing. So basically I have given a vector $\mathbf{y}$ and a matrix $\mathbf{A}$ with $y_i>0$ and $A_{ij}\in\{0,1\}$. I'm looking for a solution to the vector $\mathbf{x}$ ($x_i\geq 0$) with

$\mathbf{A}\mathbf{x}=\mathbf{y}$

where this linear equation is underdetermined.

Apparently to complete this problem several norms to minimize on $\mathbf{x}$ are possible. For my particular task it's not clear whether I need L0, L1 or L2 norm, so any solution will do - as long as it's simple. Approximate solution like iterative approaches are also fine.

Can you suggest a way to solve this problem? I'm looking for a reference to an algorithm which I can understand as a non-mathematician. Even better would be an open source implementation that I can download. And it would be perfect if it were a Python solution.

-

## migrated from dsp.stackexchange.comJan 6 '13 at 22:29

This question came from our site for practitioners of the art and science of signal, image and video processing.

Hey! From the way you stated your question, it is not obvious how it's related to signal processing, and why it would be interesting to the community. If you can provide context that puts it in to the signal processing field, we'll be glad to help. –  penelope Dec 16 '12 at 18:22
So, you want to minimize x, subject to the constraints Ax=y (equality constraint) and x >= 0 (inequality constraint). Is that correct? I think that's called a "quadratic programming" problem. I don't know a Python implementation of QP, though. Maybe you can solve it using the scipy.optimize. –  nikie Dec 17 '12 at 9:22

You can solve the system in a least-squares sense:

$$\mathbf{Ax}=\mathbf{y}$$ $$\mathbf{A^{T}Ax}=\mathbf{A^{T}y}$$ $$\mathbf{Jx}=\mathbf{r}$$ $$\mathbf{x}=\mathbf{J^{-1}r}$$

where $\mathbf{J=A^{T}A}$ and $\mathbf{r=A^{T}y}$.

Note that $\mathbf{J^{-1}r=(\mathbf{A^{T}A)^{-1}A^{T}y}}$ which is application of left pseudoinverse of $\mathbf{A}$ - this obtain the least-squares solution in $\mathbf{x}$ if $\mathbf{A}$ is overdetermined or have full rank - but this may no be our case.

The $\mathbf{J}$ is $n\times n$ and is possibly rank-deficient (underdetermined solution).

The SVD of $\mathbf{J}$ is then

$$\mathbf{J}=USV^{T}$$

where

$U$ is $n\times n$ orthogonal.

$V$ is $n\times n$ orthogonal.

$S$ is $n\times n$ diagonal, with diagonal elements $\sigma_{1} \geq \sigma_{2} \geq \cdots \geq \sigma_{n} > 0$.

The solution of your linear system is given by

$$\mathbf{x}=\mathbf{J^{-1}r}=\left(USV^{T}\right)^{-1}=VS^{-1}U^{T}\mathbf{y}$$

or more specifically:

$$\mathbf{x}=\sum_{i=1}^{n}\frac{u_{i}^{T}\mathbf{r}}{\sigma_{i}}v_{i}$$

where $u_{i}\in \mathbb{R}^{m}$ and $v_{i}\in \mathbb{R}^{n}$ are i-th columns of $U$ and $V$, respectively.

We can extend the above sum for rank-deficient cases:

$$\mathbf{x}=\sum_{\sigma_{i}\neq 0}\frac{u_{i}^{T}\mathbf{r}}{\sigma_{i}}v_{i}+\sum_{\sigma_{i}=0}\tau_{i}v_{i}$$

where $\tau_{i}$ are arbitrary coefficients (any choice of $\tau_{i}$ satisfies your linear system).

It should be noted that by choosing $\tau_{i}=0$ yield minimum-norm solution, which is usually the most desirable one in undetermined and ill-conditioned systems (where singular values are almost zero).

Source: Nocedal, Wright: "Numerical Optimization, Second Edition", chapter 10.2 Linear Least-Squares Problems, p. 250

-
Is it correct to go to the $\mathbf{x}=\mathbf{A}^{-1}\mathbf{y}$ step? I'd think you can not write that since there doesn't exist a matrix that would cancel $\mathbf{A}$ at the $\mathbf{x}$. –  Gerenuk Dec 16 '12 at 12:30
Yes, the $\mathbf{A^{-1}}$ may not exist, but SVD of $\mathbf{A}$ always exists. When you take inverse of the SVD of $\mathbf{A}$, The $V$ and $U$ matrices are orthogonal and their inverses are guaranteed to exists, but $S^{-1}$ does not exists when $\mathbf{A}$ is not invertible (because of zero signular values - division by zero). The last formula in the answers deals with the zero signular values without explicitly computing the inverse. This allows you to get over the problem of non-invertible matrices. Another approach to undetermined systems is regularization. –  Libor Dec 16 '12 at 18:19
I have updated the answer by adding a link to Wikipedia article showing the formula I have used in relation to regularization of ill-posed problems and the SVD. Hope it will help. –  Libor Dec 16 '12 at 19:36
Yes, sure the SVD exists. But how do would proceed to get your equation $\mathbf{x}=\mathbf{A}^{-1}\mathbf{y}$ from my previous one? You would premuliply the inverse of that complete SVD to both sides. But that inverse doesnt exists. So you cannot say $\mathbf{A}^{-1}\mathbf{A}\mathbf{x}=\mathbf{x}$ for the left side, or can you? –  Gerenuk Dec 17 '12 at 8:56
$\mathbf{x}=\mathbf{A^{-1}y}$ is equivalent to your first equation and all the conditions for its solution are the same, so we can "afford" doing the step. –  Libor Dec 17 '12 at 12:55

If you have MATLAB simply use the "\" or "/" operator will do the trick. It solves the equation Ax = y in a least squared sense if it's under-determined. It determines the solution by solving |Ax-y|^2 = min. If you don't have Matlab you can do it manually by defining your error E = (Ax-y) * (Ax-y)' and then calculating the partial differentials dE/dyi setting them to zero. This give and system of linear equations that matches the number of elements and can therefore be solved using a regular linear equation solver.

-
It comes close, but I'm afraid it doesn't guarantee a positive solution. –  Gerenuk Dec 16 '12 at 12:17