>  Docs Center  >  Libraries  >  Markwardt  >  QRFAC






  Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770


  Perform QR decomposition of a rectangular matrix

Major Topics

  Linear Systems

Calling Sequence

  QRFAC, A, R, [ IPVT, /PIVOT, QMATRIX=qmatrix ]


  Given an MxN matrix A (M>N), the procedure QRFAC computes the QR
  decomposition (factorization) of A. This factorization is useful
  in least squares applications solving the equation, A # x = B.
  Together with the procedure QRSOLV, this equation can be solved in
  a least squares sense.
  The QR factorization produces two matrices, Q and R, such that
    A = Q ## R
  where Q is orthogonal such that TRANSPOSE(Q)##Q equals the identity
  matrix, and R is upper triangular. This procedure does not compute
  Q directly, but returns the more-compact Householder reflectors,
  which QRSOLV applies in constructing the solution.
  Pivoting can be performed by setting the PIVOT keyword. Rows with
  the largest L2-norm are pivoted into the top positions of the
  matrix. The permutation matrix is returned in the IPVT parameter.


  A - upon input, an MxN matrix ( =XARRAY(M,N) ) to be factored,
      where M is greater than N.
      Upon output, the upper triangular MxN matrix of Householder
      reflectors used in reconstructing Q. Obviously the original
      matrix A is destroyed upon output.
      Note that the dimensions of A in this routine are the
      *TRANSPOSE* of the conventional appearance in the least
      squares matrix equation.
  R - upon ouptut, an upper triangular NxN matrix
  IPVT - upon output, the permutation indices used in partial
          pivoting. If pivoting is used, this array should be passed
          to the PIVOTS keyword of QRSOLV. If the PIVOT keyword is
          not set, then IPVT returns an unpermuted array of indices.

Keyword Parameters

  PIVOT - if set, then partial pivoting is performed, to bring the
          rows with the largest norm to the top of the matrix.
  QMATRIX - upon return, the fully explicit "Q" matrix is returned.
            This matrix is optional since the Householder vectors
            needed to solve QR problems, and to compute QMAT, are
            also stored in A. This square matrix can be used to
            perform explicit matrix multiplication (although not
            super efficiently).
  Upon return, A is in standard parameter order; A(*,IPVT) is in
  permuted order. RDIAG and QMATRIX are in permuted order upon
  return. QRSOLV accounts for these facts at the solution stage.


  Decompose the 3x2 matrix [[9.,2.,6.],[4.,8.,7.]]
    aa = [[9.,2.,6.],[4.,8.,7.]]
    qrfac, aa, r, ipvt
    IDL> print, aa
          1.81818 0.181818 0.545455
        XXXXXXXXX 1.90160 0.432573
    (position marked with Xs is undefined)
  Construct the matrix Q by expanding the Householder reflectors
  returned in AA. ( M = 3, N = 2 ) This same procedure is
  accomplished by using the QMATRIX keyword.
    ident = fltarr(m,m) ;; Construct an identity matrix
    ident(lindgen(m),lindgen(m)) = 1
    q = ident
    for i = 0, n-1 do begin
    v = aa(*,i) & if i GT 0 then v(0:i-1) = 0 ;; extract reflector
    q = q ## (ident - 2*(v # v)/total(v * v)) ;; generate matrix
  Verify that Q ## R returns to the original AA
    print, q(0:1,*) ## r
        9.00000 4.00000
        2.00000 8.00000
        6.00000 7.00000
  See example in QRSOLV to solve a least squares problem.


  More', Jorge J., "The Levenberg-Marquardt Algorithm:
    Implementation and Theory," in *Numerical Analysis*, ed. Watson,
    G. A., Lecture Notes in Mathematics 630, Springer-Verlag, 1977.

Modification History

  Written (taken from MPFIT), CM, Feb 2002
  Added usage message, error checking, CM 15 Mar 2002
  Corrected error in EXAMPLE, CM, 10 May 2002
  Now returns Q matrix explicitly if requested, CM, 14 Jul 2002
  Documented QMATRIX keyword, CM, 22 Jul 2002
  Corrected errors in computations of R and Q matrices when
    pivoting, CM, 21 May 2004
  Small correction to documentation, CM, 05 Oct 2007
  Documentation, CM, 17 Dec 2007

© 2020 Harris Geospatial Solutions, Inc. |  Legal
My Account    |    Store    |    Contact Us