in [Q, R, E] = qr(A), how E is determined ?

Marco Pappalepore — Tue, 05 May 2009 17:51:10 -0400

Hi,

I'm trying to implement qr factorization, as it works in matlab, in another programming language. I've found an algorithm to do the qr factorization (but without ordering abs(diag(R)). That's the part I'm missing. On qr's help page, it is said:

"[Q,R,E] = qr(A) for full matrix A, produces a permutation matrix E, an upper triangular matrix R with decreasing diagonal elements, and a unitary matrix Q so that A*E = Q*R. The column permutation E is chosen so that abs(diag(R)) is decreasing."

I've looked also at the fortran implementation of qr in netlib at this address :

http://www.netlib.org/lapack/double/dgeqp3.f

there I can get the permutation (=pivot) matrix is received as an input parameter:

(...)
* JPVT (input/output) INTEGER array, dimension (N)
* On entry, if JPVT(J).ne.0, the J-th column of A is permuted
* to the front of A*P (a leading column); if JPVT(J)=0,
* the J-th column of A is a free column.
* On exit, if JPVT(J)=K, then the J-th column of A*P was the
* the K-th column of A.
(...)

it seems that the pivot strategy is determined outside of the fortran routine, somewhere in matlab internals. Now I'm able to find the correct permutation matrix by searching it with a brute-force approach, i.e. simply exploring all the permutations until I find one good, that makes abs(diag(R)) decreasing; but that could bring the elaboration time to exponential or even factorial dimensions. I'm guessing if, given a not ordered Q-R factorization, is there some transformation or formula (like the gauss algorithm for finding linear independent columns) that "understands" rapidly how the permutation matrix E must be. So, how is E determined in the qr routine ?

many thanks, Marco

Re: in [Q, R, E] = qr(A), how E is determined ?

Bruno Luong — Tue, 05 May 2009 18:12:01 -0400

QR is essentially a Gram-Smidth orthogonalization procedure. During the iteration #k, it is easy to show that if we pick among the set of remaining vectors (i.e., not yet included in the span) the vector that has the *largest* orthogonal component to the current subspace (generated by k first vectors), then the diagonal of R must decrease. Essentially, this can be easily prove using triangular inequality (norm of the projection <= norm of the vector). This procedure provides a permutation E and warranty the expected property without any extra cost in QR decomposition (certainly not exponential or factorial complexity).

Bruno

MATLAB Central Newsreader - in [Q, R, E] = qr(A), how E is determined ?

in [Q, R, E] = qr(A), how E is determined ?

Re: in [Q, R, E] = qr(A), how E is determined ?