QR decomposition by LAPACK

[BoureghdaMohammed_14]

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 ?

[Julien Langou]

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


[BoureghdaMohammed_14]

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

Boureghda Mohammed.