Julia equivalent to Matlab's [Q, R, P] = qr(A)?

[Stuart Brorson]

Hello Julia users,

Matlab has a variant of the QR decomposition invoked like this:

[Q,R,P] = qr(A)

This variant of qr() returns matrix R where the diagonal elements are
sorted from largest to smallest magnitude as you go from upper left to
lower right. The matrix P is the permutation matrix which permutes
the rows/cols of A to give this ordering of the diagonal elements of
R. That is, Q*R = A*P.

I tried doing the naive, analogous thing in Julia, but get an error:

julia> A = rand(3,3)
3x3 Array{Float64,2}:
0.243071 0.454947 0.89657
0.112843 0.802457 0.375417
0.154241 0.0182734 0.992542

julia> Q,R,P = qr(A)
ERROR: BoundsError: attempt to access (
3x3 Array{Float64,2}:
-0.786117 0.0985642 -0.610168
-0.364946 -0.870763 0.329523
-0.498833 0.481723 0.720492,

3x3 Array{Float64,2}:
-0.309204 -0.659611 -1.33693
0.0 -0.645106 0.2396
0.0 0.0 0.29177)
at index [3]
in indexed_next at tuple.jl:21

My question: What's the best way to get the equivalent of Matlab's
[Q,R,P] = qr(A) in Julia? Should I write my own qr() (not too
difficult)? Or just do some row/col permutation of the output of
regular qr()?

Thanks for any advice,

Stuart

p.s. I am using Version 0.4.3-pre+6 (2015-12-11 00:38 UTC)

[Michele Zaffalon]

I think you should use qrfact(A,value{true}): http://docs.julialang.org/en/release-0.4/stdlib/linalg/?highlight=qrfact#Base.qrfact

[Chris Rackauckas]

Use qrfact() (or the in-place version). It doesn't return [Q,R,P] because Julia is smarter than that. It returns a specially designed type for the QR return. It holds all of the those parts so you can get Q, R, and P out, but most of the time you won't need to. For example, if A is the return of qrfact, then A\b will automatically dispatch to do the fast algorithm for a QR-factored matrix (which is the most common use for accessing Q and R). 

[Stuart Brorson]

Thank you, Michele & Chris for your suggestion to use qrfact. Good
idea. Unfortunately, I neglected to mention that my matrix is a
BigFloat array, and qrfact with pivoting fails like this:

julia> Z = qrfact(A, true)
ERROR: MethodError: `qrfact` has no method matching
qrfact(::Array{BigFloat,2}, ::Bool)

Hmmmm.... I guess I'll need to figure out how to implement the
QR-with-pivot algorithm in Golub and Van Loan (ugh!). Or is there a
clear reference to the algorithm Julia uses?

Stuart

[Andreas Noack]

We use LAPACK's QR with column pivoting. See http://www.netlib.org/lapack/lug/node42.html. LAPACK uses blocking for BLAS3 but that is not necessary for a generic implementation. So it's the task is just to sort the columns by norm at each step. If you want to try an implementation you can look for the reflector! and reflectorApply! functions in LinAlg.

[Andreas Noack]

My memory is too short. I just realized that I implemented a generic pivoted QR last fall so if you try out the prerelease of Julia 0.5 then you'll be able to compute the pivoted QR of a BigFloat matrix.

[Steven G. Johnson]

On Wednesday, September 7, 2016 at 7:53:41 AM UTC-4, Stuart Brorson wrote:

Thank you, Michele & Chris for your suggestion to use qrfact.  Good
idea.  Unfortunately, I neglected to mention that my matrix is a
BigFloat array, and qrfact with pivoting fails like this:

julia> Z = qrfact(A, true)

Do qrfact(A, Val{true}) in Julia 0.5