Content |
QR Factorization with Column Pivoting
In this example we compute a \(QR\) factorization with column pivoting and use it for solving a system of linear equations. Another (more interesting) of this method are rank deficient least square problems.
For a general \(m \times n\) matrix \(A\) the FLENS-LAPACK function lapack::qp3 the factorization \(AP = QR\) where \(P\) is a permutation matrix, \(Q\) a orthogonal/unitary matrix and \(R\) an upper triangular or trapezoidal matrix.
Example Code
#include <flens/flens.cxx>
using namespace std;
using namespace flens;
int
main()
{
typedef GeMatrix<FullStorage<double> > GeMatrix;
typedef typename GeMatrix::IndexType IndexType;
typedef DenseVector<Array<IndexType> > IndexVector;
typedef DenseVector<Array<double> > DenseVector;
GeMatrix A(4,4);
IndexVector jPiv;
DenseVector tau;
DenseVector b(4);
//DenseVector work;
A = 2, 3, -1, 0,
-6, -5, 0, 2,
2, -5, 6, -6,
4, 6, 2, -3;
b = 20,
-33,
-43,
49;
cout << "A = " << A << endl;
cout << "b = " << b << endl;
lapack::qp3(A, jPiv, tau);
//lapack::qp3(A, jPiv, tau, work);
lapack::ormqr(Left, Trans, A, tau, b);
//lapack::ormqr(Left, Trans, A, tau, b, work);
blas::sv(NoTrans, A.upper(), b);
DenseVector x(4);
for (IndexType i=1; i<=x.length(); ++i) {
x(jPiv(i)) = b(i);
}
cout << "x = " << x << endl;
cout << "jPiv = " << jPiv << endl;
}
Comments on Example Code
Compute the factorization \(AP = QR\). As vector jPiv has initially length zero it gets internally resized to length \(n\) and initialized with Zero. This means that all columns are free columns. Also note that the workspace gets created implicitly and temporarily. So you might not want to do this inside a loop.
//lapack::qp3(A, jPiv, tau, work);
Compute \(\tilde{b} = Q^H b\). Vector \(b\) gets overwritten with \(\tilde{b}\).
//lapack::ormqr(Left, Trans, A, tau, b, work);
Solve \(R u = \tilde{b}\). Vector \(b\) gets overwritten with \(u =P^T x\).
Compute \(x = Pu\). Note that we need an extra vector here!
for (IndexType i=1; i<=x.length(); ++i) {
x(jPiv(i)) = b(i);
}
Output some additional information:
}
Compile
$shell> cd flens/examples $shell> g++ -std=c++11 -Wall -I../.. -o lapack-geqp3 lapack-geqp3.cc
Run
$shell> cd flens/examples $shell> ./lapack-geqp3 A = 2 3 -1 0 -6 -5 0 2 2 -5 6 -6 4 6 2 -3 b = 20 -33 -43 49 x = 1 9 9 9 jPiv = 2 4 1 3
Example with Complex Numbers
Example Code
#include <flens/flens.cxx>
using namespace std;
using namespace flens;
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;
GeMatrix A(4,4);
IndexVector jPiv;
DenseVector tau;
DenseVector b(4);
//DenseVector work;
A = 2, 3, -1, 0,
-6, -5, 0, 2,
2, -5, 6, -6,
4, 6, 2, -3;
b = 20,
-33,
-43,
49;
A *= I;
b *= 2.0*I;
cout << "A = " << A << endl;
cout << "b = " << b << endl;
lapack::qp3(A, jPiv, tau);
//lapack::qp3(A, jPiv, tau, work);
lapack::unmqr(Left, ConjTrans, A, tau, b);
//lapack::unmqr(Left, ConjTrans, A, tau, b, work);
blas::sv(NoTrans, A.upper(), b);
DenseVector x(4);
for (IndexType i=1; i<=x.length(); ++i) {
x(jPiv(i)) = b(i);
}
cout << "x = " << x << endl;
cout << "jPiv = " << jPiv << endl;
}
Comments on Example Code
Make that \(A\) and \(b\) somehow look complex ;-)
b *= 2.0*I;
Compute the factorization \(AP = QR\). As vector jPiv has initially length zero it gets internally resized to length \(n\) and initialized with Zero. This means that all columns are free columns. Also note that the workspace gets created implicitly and temporarily. So you might not want to do this inside a loop.
//lapack::qp3(A, jPiv, tau, work);
Compute \(\tilde{b} = Q^H b\). Vector \(b\) gets overwritten with \(\tilde{b}\).
//lapack::unmqr(Left, ConjTrans, A, tau, b, work);
Solve \(R u = \tilde{b}\). Vector \(b\) gets overwritten with \(u =P^T x\).
Compute \(x = Pu\). Note that we need an extra vector here!
for (IndexType i=1; i<=x.length(); ++i) {
x(jPiv(i)) = b(i);
}
Output some additional information:
}
Compile
$shell> cd flens/examples $shell> g++ -DUSE_CXXLAPACK -framework vecLib -std=c++11 -Wall -I../.. -o lapack-complex-geqp3 lapack-complex-geqp3.cc
Run
$shell> cd flens/examples $shell> ./lapack-complex-geqp3 A = (0,2) (0,3) (-0,-1) (0,0) (-0,-6) (-0,-5) (0,0) (0,2) (0,2) (-0,-5) (0,6) (-0,-6) (0,4) (0,6) (0,2) (-0,-3) b = (0,40) (-0,-66) (-0,-86) (0,98) x = (2,1.29776e-14) (18,8.39706e-16) (18,4.76582e-14) (18,5.05304e-14) jPiv = 2 4 1 3