Provide support for QR decomposition

[Klaus Iglberger]

Description

The QR decomposition is often used to solve the linear least squares problem, and is the basis for a particular eigenvalue algorithm, the QR algorithm. Unfortunately, this feature is still missing in Blaze. As it is an often required algorithm, Blaze should provide support for it.

Tasks

• design the interface for the QR decomposition
• provide an implementation of the QR decomposition that ...
  • ... is robust
  • ... is as efficient as possible
  • ... can deal with errors in a graceful way
• provide a shared-memory parallel implementation (if possible)
• provide a full documentation of the feature, including ...
  • ... an explanation of the algorithm
  • ... limitations of the algorithms
  • ... examples of use
  • ... example computations
• ensure compatibility to all existing matrix types
• add appropriate test cases for the feature

[Klaus Iglberger]

The Blaze library already provides LAPACK wrapper functions for the QR decomposition of a dense matrix:

namespace blaze {

void geqrf( int m, int n, float* A, int lda, float* tau, float* work, int lwork, int* info );
void geqrf( int m, int n, double* A, int lda, double* tau, double* work, int lwork, int* info );
void geqrf( int m, int n, complex<float>* A, int lda, complex<float>* tau, complex<float>* work, int lwork, int* info );
void geqrf( int m, int n, complex<double>* A, int lda, complex<double>* tau, complex<double>* work, int lwork, int* info );

template< typename MT, bool SO >
void geqrf( DenseMatrix<MT,SO>& A, typename MT::ElementType* tau );

} // namespace blaze

These functions immediately available via cloning the Blaze repository and will be officially released in Blaze 2.6.

[Klaus Iglberger]

Summary

The feature has been implemented, tested, optimized, and documented as required. It is immediately available via cloning the Blaze repository and will be officially released in Blaze 2.6.

QR Decomposition

The QR decomposition of a dense matrix can be computed via the qr() function:

blaze::DynamicMatrix<double,blaze::rowMajor> A;
// ... Resizing and initialization

blaze::DynamicMatrix<double,blaze::columnMajor> Q;
blaze::DynamicMatrix<double,blaze::rowMajor> R;

qr( A, Q, R );  // QR decomposition of a row-major matrix

assert( A == Q * R );

The function works for both rowMajor and columnMajor matrices and the three matrices A, Q and R can have any storage order.

Furthermore, qr() can be used with adaptors. For instance, the following example demonstrates the QR decomposition of a symmetric matrix into a general matrix and an upper triangular matrix:

blaze::SymmetricMatrix< blaze::DynamicMatrix<double,blaze::columnMajor> > A;
// ... Resizing and initialization

blaze::DynamicMatrix<double,blaze::rowMajor> Q;
blaze::UpperMatrix< blaze::DynamicMatrix<double,blaze::columnMajor> > R;

qr( A, Q, R );  // QR decomposition of A

Note that the QR decomposition can only be used for dense matrices with float, double, complex<float> or complex<double> element type. The attempt to call the function with matrices of any other element type or with a sparse matrix results in a compile time error!

Also note that the qr() function decomposes a dense matrix by means of LAPACK kernels. Thus the function can only be used if the fitting LAPACK library is available and linked to the executable. Otherwise a linker error will be created.

[Mikhail Katliar]

According to Blaze documentation, the qr() function returns matrix Q which is “a general m-by-min(m,n) matrix”. This does not match the definition of QR decomposition, which implies that the Q matrix must be orthogonal (and hence square): “QR decomposition, also known as a QR factorization, is a decomposition of a matrix A into a product A = QR of an orthogonal matrix Q and an upper triangular matrix R.” The same is stated in the documentation of SGEQRF LAPACK function, which is used to implement blaze::qr():

A is REAL array, dimension (LDA,N)
On entry, the M-by-N matrix A.
On exit, the elements on and above the diagonal of the array
contain the min(M,N)-by-N upper trapezoidal matrix R (R is
upper triangular if m >= n); the elements below the diagonal,
with the array TAU, represent the orthogonal matrix Q as a
product of min(m,n) elementary reflectors (see Further
Details).

Is there an easy way to recover the orthogonal Q and an m-by-n R from the matrices returned by blaze::qr()? If no, I suggest changing blaze::qr() function such that it returns an orthogonal m-by-m Q and an m-by-n R.

[Mikhail Katliar]

qr() calls LAPACK orgqr function which recovers the first columns of the Q matrix. In order to recover the entire Q matrix, there is another LAPACK function orgr2.