Content

Least Square Solution

lapack::ls solves an overdetermined or underdetermined real linear systems involving an \(m \times n\) matrix \(A\), or its transpose, using a \(QR\) or \(LQ\) factorization of \(A\). It is assumed that \(A\) has full rank.

The following options are provided:

Case: Least Squares Solution of an Overdetermined System

In this case we compute the least squares solution of an overdetermined system, i.e., we solve the least squares problem: minimize \(\| B - AX \|\) where

On input the FLENS-LAPACK function lapack::ls receives \(A\) and \(B\). On exit the first \(n\) rows of \(B\) will be overwritten with the least square solution \(X\). Matrix (or vector) views can be helpful and convenient. Just let \(X\) be a matrix view (or vector) that references the first \(n\) rows (or elements) of \(X\).

Example Code

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

using namespace flens;
using namespace std;

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

    const Underscore<IndexType>  _;

    GeMatrix     A(43);
    DenseVector  b(4);

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

    b =  6,
        15,
        24,
        34;

    auto x = b(_(1,3));

    lapack::ls(NoTrans, A, b);

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

    return 0;
}

Comments on Example Code

Setup the matrix \(A\).

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

Setup vector \(b\)

    b =  6,
        15,
        24,
        34;

The least square solution \(x\) will be stored in the first three components of vector \(b\). So vectors view come in handy:

    auto x = b(_(1,3));

Solve \(\min\limits_{x} \|Ax - b \|\)

    lapack::ls(NoTrans, A, b);

Compile

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

Run

$shell> cd flens/examples                                                       
$shell> ./lapack-gels-case1                                                     
x = 
            1              1              1 

Case: Minimum Norm Solution of an Undetermined System

In this case we compute the minimum norm solution of the underdetermined system \(A^T X = B\) where

On input the first \(n\) rows of matrix \(X\) contain the matrix \(B\). Then the FLENS-LAPACK function lapack::ls receives \(A\) and \(X\). On exit \(X\) will be overwritten with the minimal norm solution. Again, matrix views can be convenient as \(X\) and \(B\) partially share the same data. More precise, matrix \(B\) is a matrix view that references the first \(n\) rows of \(X\).

Example Code

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

using namespace flens;
using namespace std;

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

    const Underscore<IndexType>  _;

    GeMatrix     A(43);
    DenseVector  x(4);

    A = 15,  9,
        2610,
        3711,
        4812;

    auto b = x(_(1,3));

    b = 30,
        70,
       110;

    lapack::ls(Trans, A, x);

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

    return 0;
}

Comments on Example Code

Setup the matrix \(A\).

    A = 15,  9,
        2610,
        3711,
        4812;

Initially the first the elements of \(x\) contain the right hand side vector \(b\). For convenience we define a vector view \(b\) that references these elements.

    auto b = x(_(1,3));

Setup vector \(b\).

    b = 30,
        70,
       110;

Find the minimal norm solution of \(A^T x = b\).

    lapack::ls(Trans, A, x);

Compile

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

Run

$shell> cd flens/examples                                                       
$shell> ./lapack-gels-case2                                                     
x = 
            1              2              3              4 

Case: Minimum Norm Solution of an Undetermined System

In this case we compute the minimum norm solution of the underdetermined system \(AX = B\) where

On input the first \(m\) rows of matrix \(X\) contain the matrix \(B\). The FLENS-LAPACK function lapack::ls receives \(A\) and \(X\). On exit \(X\) will be overwritten with the minimal norm solution. Again, matrix views can be convenient as \(X\) and \(B\) partially share the same data. More precise, matrix \(B\) is a matrix view that references the first \(m\) rows of \(X\).

Example Code

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

using namespace flens;
using namespace std;

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

    const Underscore<IndexType>  _;

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

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

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

    lapack::ls(NoTrans, A, x);

    cout << "x = " << x << 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;

Find the minimal norm solution of \(Ax = b\).

    lapack::ls(NoTrans, A, x);

Compile

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

Run

$shell> cd flens/examples                                                       
$shell> ./lapack-gels-case3                                                     
x = 
            1              2              3              4 

Case: Least Squares Solution of an Overdetermined System

In this case we compute the least squares solution of an overdetermined system, i.e., we solve the least squares problem: minimize \(\| B - A^T X \|\) where

On input the FLENS-LAPACK function lapack::ls receives \(A\) and \(B\). On exit the first \(m\) rows of \(B\) will be overwritten with the least square solution \(X\). Matrix (or vector) views can be helpful and convenient. Just let \(X\) be a matrix view (or vector) that references the first \(m\) rows (or elements) of \(X\).

Example Code

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

using namespace flens;
using namespace std;

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

    const Underscore<IndexType>  _;

    GeMatrix     A(34);
    DenseVector  b(4);

    A = 14710,
        25811,
        36913;

    b =  6,
        15,
        24,
        34;

    auto x = b(_(1,3));

    lapack::ls(Trans, A, b);

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

    return 0;
}

Comments on Example Code

Setup the matrix \(A\).

    A = 14710,
        25811,
        36913;

Setup vector \(b\)

    b =  6,
        15,
        24,
        34;

The least square solution \(x\) will be stored in the first three components of vector \(b\). So vectors view come in handy:

    auto x = b(_(1,3));

Solve \(\min\limits_{x} \|A^T x - b \|\)

    lapack::ls(Trans, A, b);

Compile

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

Run

$shell> cd flens/examples                                                       
$shell> ./lapack-gels-case4                                                     
x = 
            1              1              1 

Examples with Complex Numbers

Currently complex numbers for lapack::ls require an external LAPACK implementation as computational back-end.

Case: Least Squares Solution of an Overdetermined System

In this case we compute the least squares solution of an overdetermined system, i.e., we solve the least squares problem: minimize \(\| B - AX \|\) where

On input the FLENS-LAPACK function lapack::ls receives \(A\) and \(B\). On exit the first \(n\) rows of \(B\) will be overwritten with the least square solution \(X\). Matrix (or vector) views can be helpful and convenient. Just let \(X\) be a matrix view (or vector) that references the first \(n\) rows (or elements) of \(X\).

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 DenseVector<Array<Complex> >     DenseVector;
    typedef typename DenseVector::IndexType  IndexType;

    const Underscore<IndexType>  _;

    GeMatrix     A(43);
    DenseVector  b(4);

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

    b =  6,
        15,
        24,
        34;

    A *= I;
    b *= I;

    auto x = b(_(1,3));

    lapack::ls(NoTrans, A, b);

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

    return 0;
}

Comments on Example Code

Setup the matrix \(A\).

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

Setup vector \(b\)

    b =  6,
        15,
        24,
        34;

Make the matrix/vectors somehow complex

    A *= I;
    b *= I;

The least square solution \(x\) will be stored in the first three components of vector \(b\). So vectors view come in handy:

    auto x = b(_(1,3));

Solve \(\min\limits_{x} \|Ax - b \|\)

    lapack::ls(NoTrans, A, b);

Compile

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

Run

$shell> cd flens/examples                                                       
$shell> ./lapack-complex-gels-case1                                             
x = 
             (1,2.92971e-15)               (1,2.80786e-15)              (1,-4.66206e-15) 

Case: Minimum Norm Solution of an Undetermined System

In this case we compute the minimum norm solution of the underdetermined system \(A^T X = B\) where

On input the first \(n\) rows of matrix \(X\) contain the matrix \(B\). Then the FLENS-LAPACK function lapack::ls receives \(A\) and \(X\). On exit \(X\) will be overwritten with the minimal norm solution. Again, matrix views can be convenient as \(X\) and \(B\) partially share the same data. More precise, matrix \(B\) is a matrix view that references the first \(n\) rows of \(X\).

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 DenseVector<Array<Complex> >     DenseVector;
    typedef typename DenseVector::IndexType  IndexType;

    const Underscore<IndexType>  _;

    GeMatrix     A(43);
    DenseVector  x(4);

    A = 15,  9,
        2610,
        3711,
        4812;

    auto b = x(_(1,3));

    b = 30,
        70,
       110;

    A *= I;
    b *= 2.0*I;

    lapack::ls(ConjTrans, A, x);

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

    return 0;
}

Comments on Example Code

Setup the matrix \(A\).

    A = 15,  9,
        2610,
        3711,
        4812;

Initially the first the elements of \(x\) contain the right hand side vector \(b\). For convenience we define a vector view \(b\) that references these elements.

    auto b = x(_(1,3));

Setup vector \(b\).

    b = 30,
        70,
       110;

Make the matrix/vectors somehow complex

    A *= I;
    b *= 2.0*I;

Find the minimal norm solution of \(A^H x = b\).

    lapack::ls(ConjTrans, A, x);

Compile

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

Run

$shell> cd flens/examples                                                       
$shell> ./lapack-complex-gels-case2                                             
x = 
                      (-2,0)                        (-4,0)                        (-6,0)                        (-8,0) 

Case: Minimum Norm Solution of an Undetermined System

In this case we compute the minimum norm solution of the underdetermined system \(AX = B\) where

On input the first \(m\) rows of matrix \(X\) contain the matrix \(B\). The FLENS-LAPACK function lapack::ls receives \(A\) and \(X\). On exit \(X\) will be overwritten with the minimal norm solution. Again, matrix views can be convenient as \(X\) and \(B\) partially share the same data. More precise, matrix \(B\) is a matrix view that references the first \(m\) rows of \(X\).

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 DenseVector<Array<Complex> >     DenseVector;
    typedef typename DenseVector::IndexType  IndexType;

    const Underscore<IndexType>  _;

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

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

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

    A *= I;
    b *= 4.0*I;
    //b *= 3.0*I;

    lapack::ls(NoTrans, A, x);

    cout << "x = " << x << 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 the matrix/vectors somehow complex.

    A *= I;
    b *= 4.0*I;
    //b *= 3.0*I;

Find the minimal norm solution of \(Ax = b\).

    lapack::ls(NoTrans, A, x);

Compile

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

Run

$shell> cd flens/examples                                                       
$shell> ./lapack-complex-gels-case3                                             
x = 
                       (4,0)                         (8,0)                        (12,0)                        (16,0) 

Case: Least Squares Solution of an Overdetermined System

In this case we compute the least squares solution of an overdetermined system, i.e., we solve the least squares problem: minimize \(\| B - A^T X \|\) where

On input the FLENS-LAPACK function lapack::ls receives \(A\) and \(B\). On exit the first \(m\) rows of \(B\) will be overwritten with the least square solution \(X\). Matrix (or vector) views can be helpful and convenient. Just let \(X\) be a matrix view (or vector) that references the first \(m\) rows (or elements) of \(X\).

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 DenseVector<Array<Complex> >     DenseVector;
    typedef typename DenseVector::IndexType  IndexType;

    const Underscore<IndexType>  _;

    GeMatrix     A(34);
    DenseVector  b(4);

    A = 14710,
        25811,
        36913;

    b =  6,
        15,
        24,
        34;

    A *= I;
    b *= 4.0*I;

    auto x = b(_(1,3));

    lapack::ls(ConjTrans, A, b);

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

    return 0;
}

Comments on Example Code

Setup the matrix \(A\).

    A = 14710,
        25811,
        36913;

Setup vector \(b\)

    b =  6,
        15,
        24,
        34;

Make the matrix/vectors somehow complex

    A *= I;
    b *= 4.0*I;

The least square solution \(x\) will be stored in the first three components of vector \(b\). So vectors view come in handy:

    auto x = b(_(1,3));

Solve \(\min\limits_{x} \|A^H x - b \|\)

    lapack::ls(ConjTrans, A, b);

Compile

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

Run

$shell> cd flens/examples                                                       
$shell> ./lapack-complex-gels-case4                                             
x = 
            (-4,1.33126e-14)              (-4,8.04828e-15)             (-4,-1.70266e-14)