Modified Gram-Schmidt QR-decomposition

Describe the feature you would like to be implemented.

Modified Gram-Schmidt (MGS) QR decomposition for dense matrices should be added. So far you only have Householder (HH) QR decompositions.

In particular, if I provide matrix A of size m x n, where m>=n, and matrix R of size n x n, then there should be an inplace MGS procedure

qr(A,R)

that yields

A_before = A_after * R

such that cond(A_after)=1 and R upper triangular. (A is to be overwritten.)

Would such a feature be useful for other users? Why?

MGS has several advantages.

  1. Storage:

Though HH can be computed in-place, there is one additional vector needed of size n for storing the reflectors' first values v(1). Thus, if I want to use one factorization routine for several matrices A of different dimensions then each time eigen will interiorly a new vector of different size n. This can ruin memory management.

  1. Performance:

In principle, factoring A into QR requires 2mnn flops, regardless whether we use MGS or HH. HH will only drop that cost to 4/3mnn when Q is not factored. However, I have applications where I will compute BQ=AQQ and BR=AR/R after computation of the factorization.

In that case, applying the factorization of HH is way more tedious because I either have to use some phrase like "matrixQR.apply_Q_from_right(AQ)", whereas in MGS I can literally write BQ=AQ*Q, which is much much faster.

Lastly, if I want to compute many products like BQ=AQ*Q, or when QB has many rows, then probably I want to compute Q even in the case of HH. Now in that case, this is impossible to do in-place on the matrix A that holds the QR factorization. Notice that the comparision with two matrices of MGS is not unfair because the matrix R typically requires much less storage than A when m>>n.

  1. Accuracy:

The matrix Q, computed by MGS, is more accurate than that of HH. The accuracy of MGS can be even enhanced further by using iterative refinement. HH cannot do that. I agree HH has been analyzed thoroughly, but it remains limited.

  1. Comfort:

It is more user-friendly in terms of storage and computations when they can supply two matrices of appropriate sizes and obtain the factors Q and R. They can then use Q and R in subsequent computations without the need to involve specially crafted Eigen-syntax-sensitive functions for HH; when Q is a thing that should not need to have anything to do with HH.

Any hints on how to implement the requested feature?

Yes.

[m,n]=size(A);

for i=1:n

for j=1:(n-1)

R(j,i) = A(:,j).dot( A(:,i) );

A(:,i)-= A:,j) * R(j,i);

end

R(i,i) = A(:,j).norm();

A(:,i)/=R(i,i)

end

Additional resources

wikipedia :P

Edited by Martin Neuenhofen