# Rank Deficient Least Square Problems

Computes the minimum-norm solution to a real linear least squares problem: minimize $$\| A X - B \|$$ using a complete orthogonal factorization of $$A$$. $$A$$ is an $$m \times n$$ matrix which may be rank-deficient. The rank of $$A$$ gets determined using a incremental condition estimation.

The routine first computes a $$QR$$ factorization with column pivoting:

$A P = Q \begin{pmatrix} R_{11} & R_{12} \\ \mathbb{0} & R_{22} \end{pmatrix}$

with $$R_{11}$$ defined as the largest leading submatrix whose estimated condition number is less than $$1/\text{rCond}$$ (where $$\text{rCond}$$ is a input parameter). The order of $$R_{11}$$ is the effective rank of $$A$$.

Then, $$R_{22}$$ is considered to be negligible, and $$R_{12}$$ is annihilated by unitary transformations from the right, arriving at the complete orthogonal factorization:

$A P = Q \begin{pmatrix} T_{11} & \mathbb{0} \\ \mathbb{0} & \mathbb{0} \end{pmatrix} Z = \begin{pmatrix} Q_1 & Q_2 \end{pmatrix} \begin{pmatrix} T_{11} & \mathbb{0} \\ \mathbb{0} & \mathbb{0} \end{pmatrix} Z$

The minimum-norm solution is then

$X = P Z^T \begin{pmatrix} T_{11}^{-1} Q_1^T B \\ \mathbb{0} \end{pmatrix}$

where $$Q_1$$ consists of the first $$\text{rank}$$ columns of $$Q$$.

## Example Code

#include <flens/flens.cxx>
#include <iostream>

using namespace flens;
using namespace std;

int
main()
{
typedef GeMatrix<FullStorage<double> >   GeMatrix;
typedef typename GeMatrix::IndexType     IndexType;
typedef DenseVector<Array<IndexType> >   IndexVector;
typedef DenseVector<Array<double> >      DenseVector;

const Underscore<IndexType>  _;

GeMatrix     A(34);
DenseVector  x(4);
IndexVector  jPiv(4);

A =  1,  2,  3,  4,
5,  6,  7,  8,
9101112;

auto b = x(_(1,3));
b =  30,
70,
110;

jPiv = 0;

double    rCond = 1E-8;

IndexType rank = lapack::lsy(A, x, jPiv, rCond);
cout << "x = " << x << endl;

cout << endl << endl << endl
<< "Some additional information:"
<< endl;

cout << "jPiv = " << jPiv << endl;
cout << "rank = " << rank << endl;

return 0;
}

## Comments on Example Code

Setup the $$3 \times 4$$ matrix $$A$$.

A =  1,  2,  3,  4,
5,  6,  7,  8,
9101112;

The right hand side vector $$b$$ will be stored in vector $$x$$.

auto b = x(_(1,3));
b =  30,
70,
110;

All columns are free columns.

jPiv = 0;

We want that the condition of the leading submatrix $$R_{11}$$ is less than 1/rCond.

double    rCond = 1E-8;

Find the minimal norm solution of $$Ax = b$$.

IndexType rank = lapack::lsy(A, x, jPiv, rCond);
cout << "x = " << x << endl;

Output of permutation $$P$$ and the effective rank.

cout << endl << endl << endl
<< "Some additional information:"
<< endl;

## Compile

$shell> cd flens/examples$shell> g++ -std=c++11 -Wall -I../.. -o lapack-gelsy lapack-gelsy.cc


## Run

$shell> cd flens/examples$shell> ./lapack-gelsy
x =
1              2              3              4
Some additional information:
jPiv =
4              1              3              2
rank = 2


## Example with Complex Numbers

### Example Code

#include <flens/flens.cxx>
#include <iostream>

using namespace flens;
using namespace std;

int
main()
{
typedef complex<double>                  Complex;
const Complex                            I(0,1);

typedef GeMatrix<FullStorage<Complex> >  GeMatrix;
typedef typename GeMatrix::IndexType     IndexType;
typedef DenseVector<Array<IndexType> >   IndexVector;
typedef DenseVector<Array<Complex> >     DenseVector;

const Underscore<IndexType>  _;

GeMatrix     A(34);
DenseVector  x(4);
IndexVector  jPiv(4);

A =  1,  2,  3,  4,
5,  6,  7,  8,
9101112;

auto b = x(_(1,3));
b =  30,
70,
110;

A *= I;
b *= 3.0*I;

jPiv = 0;

double    rCond = 1E-8;

IndexType rank = lapack::lsy(A, x, jPiv, rCond);

cout << "x = " << x << endl;

cout << endl << endl << endl
<< "Some additional information:"
<< endl;

cout << "jPiv = " << jPiv << endl;
cout << "rank = " << rank << endl;

return 0;
}

### Comments on Example Code

Setup the $$3 \times 4$$ matrix $$A$$.

A =  1,  2,  3,  4,
5,  6,  7,  8,
9101112;

The right hand side vector $$b$$ will be stored in vector $$x$$.

auto b = x(_(1,3));
b =  30,
70,
110;

Make $$A$$ and $$b$$ somehow complex ;-)

A *= I;
b *= 3.0*I;

All columns are free columns.

jPiv = 0;

We want that the condition of the leading submatrix $$R_{11}$$ is less than 1/rCond.

double    rCond = 1E-8;

Find the minimal norm solution of $$Ax = b$$.

IndexType rank = lapack::lsy(A, x, jPiv, rCond);

Output of permuation $$P$$ and the effective rank.

cout << endl << endl << endl
<< "Some additional information:"
<< endl;

### Compile

$shell> cd flens/examples$shell> g++ -DUSE_CXXLAPACK -framework vecLib -std=c++11 -Wall -I../.. -o lapack-complex-gelsy lapack-complex-gelsy.cc


### Run

$shell> cd flens/examples$shell> ./lapack-complex-gelsy
x =
(3,6.62737e-15)               (6,1.23591e-15)              (9,-4.75902e-15)             (12,-1.18941e-15)
Some additional information:
jPiv =
4              1              3              2
rank = 2