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

flens/examples/lapack-gelsy.cc

#include <cxxstd/iostream.h>

#include <flens/flens.cxx>

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(3, 4);

    DenseVector  x(4);

    IndexVector  jPiv(4);

    A =  1,  2,  3,  4,

         5,  6,  7,  8,

         9, 10, 11, 12;

    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,

         9, 10, 11, 12;

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