[Gretl-devel] Re: qrdecomp with column pivoting?
On Fri, Jun 3, 2022 at 11:05 AM Sven Schreiber <svetosch(a)gmx.net> wrote:
I'm wondering whether it would be useful (or more precisely, whether the
cost-benefit calculation would be net positive...) to generalize the
qrdecomp() function to allow column pivoting. Here are a couple of
thoughts and remarks on that:
1) First of all, qrdecomp currently does not return an ordered result.
For example, when doing this:
Q = qrdecomp(mnormal(4)~zeros(4)~mnormal(4) , &R)
then the middle column of R is zero, just following the input. Of
course, this is fully according to spec, nothing else is promised currently.
2) A QR decomposition with pivoting would provide a rank-revealing
operation. The natural workaround and alternative way to do this is SVD.
It is my understanding that SVD would be noticeably slower, however.
(Correct?)
Maybe, though I gather that for a matrix with a lot more rows than
columns the execution time is not that different. Note: we currently
assess rank using regular QR, by counting the R elements greater than
some specified "tiny" value.
3) In Lapack this is given here:
https://www.netlib.org/lapack/lug/node42.html, and references the
routine xGEQP3 (Level 3 BLAS).
There's an interface to that in libgretl's gretl_matrix.c, namely
gretl_matrix_QR_pivot_decomp(), though it's not hooked up to any hansl
function at present. I find it kinda confusing trying to deal with
results after pivoting. You need to pay careful attention to the
reordering to avoid producing nonsense.
4) In terms of a possible interface, the immediate idea would be to
have
a third and optional boolean switch which defaults to the old behavior,
and if set to TRUE then it would use column pivoting. However, perhaps
one also wants to retrieve the resulting permutation matrix P in the
relationship AP=QR (where A is the input). Then perhaps the third
optional argument would be another pointer-to-matrix, just like the
second one. Then column pivoting would only be done if a matrix pointer
is actually provided in the third slot, capturing the P result.
So, what do you think, how difficult would it be to implement this, and
does it make sense?
It "makes sense" OK, and it's already part-way implemented. Whether
it's worthwhile -- I'm agnostic.
Allin