QR decomposition by LAPACK

Open discussion regarding features, bugs, issues, vendors, etc.

QR decomposition by LAPACK

Matlab use LAPACK to calculate the QR decomposition, the algorithm is based on Householder reflections.
In mathematics I know that the elementary Householder reflector used by the QR decomposition is: F = I - 2 (uu')/(u'u). " ' ": transpose operator and I: identity matrix.
Can anyone provide me with a mathematical form of the householder reflectors used by LAPACK to calculate the QR decomposition?
If it use a reflector Like the one mentioned, What is the mathematical form of u ?
BoureghdaMohammed_14

Posts: 2
Joined: Fri Feb 03, 2017 12:09 pm

Re: QR decomposition by LAPACK

The LAPACK subroutine to compute the Householder reflection vector (`u`) is DLARFG. To apply the reflector is DLARF.

So to answer your question, please look at DLARFG. This uses the standard numerical linear algebra formula. We do the +/- trick for stability reasons. (So this is either I - 2 (uu')/(u'u) or I + 2 (uu')/(u'u).) See any numerical linear algebra textbook for explanation. Since the scalar 2/(u'u) would need to be computed at each application of the reflection, we prefer to store it. This is `TAU`. So the formula is rewritten as I - TAU (uu'), where TAU = 2/(u'u). Also note that in the formula I - 2 (uu')/(u'u), the scaling of u is not specified. (Any scaling works.) The scaling used by LAPACK is such that the top entry in u is a `1`. (Why not.) So in practice we do not store the top entry of `u`. (Since we know it is `1`.)

For an example of use of DLARFG and DLARF, see DGEQR2.

You may also be interested in looking at DLARFGP. DLARFGP does not use the +/- trick for stability, it uses another trick (so called Parlett's trick).

Then there are subroutines to block the application of `NB` reflections. First to form the `T` matrix is DLARFT, (pre-processing,) and then to apply a sequence of `NB` reflections is DLARFB. This is interesting for performance reasons. For an example of use, please look at DGEQRF.

I am attaching a quick matlab code that reproduces dlarfg.m. Main difference is that: `sqrt ( alpha^2 + xnorm^2 )` is done in a unsafe way in the matlab code below.

Julien.

Code: Select all
function [tau,y] = dlarfg( y )

%     ..
%
%  Purpose
%  =======
%
%  DLARFG generates a real elementary reflector H of order n, such
%  that
%
%        H * y_in   = ( y_out(1) ),   H' * H = I.
%                     (   0      )
%
%  H is represented in the form
%
%        H = I - tau * (      1       ) * ( 1 y_out(2:n-1)' ) ,
%                      ( y_out(2:n-1) )
%
%  If the elements of y_in(2:n-1) are all zero, then tau = 0 and H is taken to be
%  the unit matrix.
%
%  Otherwise  1 <= tau <= 2.
%

n = size(y,1);
alpha = y(1);

if( n <= 1 ),
tau = 0;
return
end

xnorm = norm( y(2:n,1) ) ;

if( xnorm == 0 ),

tau = 0 ;
else

beta = - sign( alpha ) * sqrt ( alpha^2 + xnorm^2 ) ;
tau = ( beta - alpha ) / beta ;
y( 2:n ) = y( 2:n ) / ( alpha - beta );
y(1) = beta;

end

end
Julien Langou

Posts: 834
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Re: QR decomposition by LAPACK

Thank you so much, Julien for the reply, you helped me a lot.

Boureghda Mohammed.
BoureghdaMohammed_14

Posts: 2
Joined: Fri Feb 03, 2017 12:09 pm