QRFAC Name
QRFAC
Author
Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
craigm@lheamail.gsfc.nasa.gov Purpose
Perform QR decomposition of a rectangular matrix
Major Topics
Linear Systems
Calling Sequence
QRFAC, A, R, [ IPVT, /PIVOT, QMATRIX=qmatrix ]
Description
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.
Parameters
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).
IMPLEMENTATION NOTE:
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.
Example
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
endfor
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
(transposed)
See example in QRSOLV to solve a least squares problem.
References
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