An interesting linear algebra problem from signal processing. 

This problem was presented on https://stackoverflow.com/questions/52373413/solve-underdetermined-system-of-equations-for-sparse-solution/52404876#52404876 
It appears that a discussion devoted to it might be interesting, in order to analyse possible approaches to solution. 

We approach the problem by transforming the equation into a form more "recognizable" from the classic linear algebra approach to the solution of systems of linear equations. 

The last expression, essentially, represents a set of m systems of linear equations with one and the same matrix B^T and different right-hand sides, which are the columns of C^T. In order to proceed further, we need to perform Gaussian elimination for the matrix B^T , simultaneously with all the right-hand sides, using some pivoting algorithm. As the matrix is singular, the resulting system of equations, for each particular column j from C^T and, correspondingly, for each particular column j from A^T, will have the following form: 

Since the matrix of the system is singular, k<n , where k=rank(B^T) . The vector of unknowns consists of two smaller vectors: 

 of the size k , containing the components from 1st to k-th of the column j of the A^T

                matrix;

  of the size n-k, containing the components from (k+1)th to n-th of the same                           column of  A^T. 

Correspondingly:

  of the size k, containing (permuted, due to the Gaussian elimination with column                    pivoting performed) components of the j-th column of C^T ;

  of the size n-k, containing the remaining components  of the j-th column of C^T . 

Therefore, upon completion of the Gaussian elimination, the last few equations in the system are trivial:

Therefore, two possibilities only exist: 


1) All of the components of the lower part of the right-hand side vector are zeros:

In that case, the system has infinite number of solutions, since any alpha-vector would satisfy the subsystem. Since our goal is to find a solution of the original system as "sparse" as possible, we can safely set all components of the alpha-vector to zeros, and solve the remaining equations by backward substitution, thus obtaining a unique solution. 

2) At least one component of the right-hand side (ie. of the chi-vector) of the above subsystem is different from zero. In that case, the subsystem has no solutions at all, and, therefore, the original system has no solutions as well. In such case, the only way to proceed is to try to find an approximation to the solution in terms of the least squares. 

We will perform the analysis from the very beginning, using a well-established method of Singular Value Decomposition, or SVD for short. Given a system of linear equations, say, Fx=g , matrix F can be represented as the product of multiplication of three matrices (we limit ourselves by case of square matrix F , i.e. in the context of the current problem): 

In the context of the problem which is being discussed here, the diagonal matrix does contain at least one zero element on the diagonal, since the original matrix is supposed to be singular. Here, we will consider a simple SLAE of the size 5x5 where 

- two rows of the matrix are linearly dependent, and

- the system is incompatible, so that no exact solution exists, and therefore we will try to find a least squares approximation to the solution.

The 5-th row of the matrix equals the product of multiplication of the 3rd row by the factor of -2. Now, if we select the right-hand side of the system in such a way that the 5th element of the vector g does not equal the 3rd element multiplied by -2 , the whole SLAE is incompatible, and we can only try to find a least squares approximation to the solution. 

Lets put the right-hand side vector g as follows: g = [ 1 ,  2,  3,  4,  5 ]^T 

The SLAE is, obviously, incompatible. 

Now, what the SVD algorithm delivers:

The elements of the diagonal matrix, named here as singular values, contain one zero element, as expected. The SVD-algorithm prescribes the following procedure to obtain a least-squares solution of the system:

1) Introduce vector d=U^Tg .

2) Calculate vector z as follows:

   

                                z_j = d_j / s_j  , if s_j is a non-zero value;

                                z_j can be chosen arbitrarily, if s_j = 0. 

3) The solution is obtained as the following product: x=Vz . 

And here comes the beauty of the SVD-method: regardless of the chosen value of z_j for which  s_j = 0, the program will deliver a solution vector from the R⁵ space such that the norm of the residuals vector will be the same, and it will be the least possible value for this system. (Obviously, if the system was compatible, the residuals would all equal zero). 

We provide just three cases to demonstrate that for our example. Namely, we set z₅ to 0., to 1., and to a randomly selected value, and see what happens. 

As we can see, each particular choice of z₅ leads to a different solution vector. Nevertheless, the residuals vector is one and the same in all cases, which means that all three solutions are equally valid in terms of least squares. 

Since there is an infinite number of least-squares solutions, each particular one corresponding to one particular value of chosen z₅, we can probably consider an opposite task: set one component of the solution vector to zero, and find another components such that the vector of the residuals will stay the same (obviously, some particular value of z will be calculated along the way). The more singular values s_i  equal zero, the more components of the solution vector can theoretically be made equal zero. We should recall that our original goal was to deliver as sparse matrix A as possible. We consider that as possible subject of further algorithm development.