## 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;   endend`
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