11:05 a.m.
Hi all,
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?)
3) In Lapack this is given here:
https://www.netlib.org/lapack/lug/node42.html, and references the
routine xGEQP3 (Level 3 BLAS).
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?
thanks
sven
4:46 p.m.
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
6:12 a.m.
Am 03.06.2022 um 22:46 schrieb Cottrell, Allin:
On Fri, Jun 3, 2022 at 11:05 AM Sven Schreiber
<svetosch(a)gmx.net> wrote:
> 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.
OK, but I'm thinking more about rank-deficient square matrices.
Note: we currently
assess rank using regular QR, by counting the R elements greater
than
some specified "tiny" value.
Aha? The doc for "rank" says: "numerically computed via the singular
value decomposition" - so that's not true anymore?
And with "the R elements" you mean the ones on the diagonal, or all of
them (plus some assessment of the columns)?
> 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.
OK, good to know.
> 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.
Well, let's keep that possibility in mind, then.
thanks
sven
8:40 a.m.
On Sat, 4 Jun 2022, Sven Schreiber wrote:
Am 03.06.2022 um 22:46 schrieb Cottrell, Allin:
> Note: we currently assess rank using regular QR, by counting the
> R elements greater than some specified "tiny" value.
Aha? The doc for "rank" says: "numerically computed via the singular
value decomposition" - so that's not true anymore?
The hansl rank() function does use SVD as advertised, but internally
we use QR to assess rank in some contexts, using the approach
illustrated in Jack's reply.
Allin
7:11 a.m.
On Fri, 3 Jun 2022, Cottrell, Allin wrote:
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.
in Hansl:
<hansl>
R = {}
Q = qrdecomp(mnormal(4)~zeros(4)~mnormal(4) , &R)
scalar r = sumc(abs(diag(R)).>1.0e-12)
print r
</hansl>
> 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.
Me too.
> 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.
I'm not really sure what the benefit would be from having column pivoting,
which is relatively easy to do in Hansl anyway (note that P is a
permutation matrix, so it can be more compactly expressed as a set of
indices), although of course numerically speaking lapack is unbeatable.
Could you provide us with a use case where pivoting would be useful?
-------------------------------------------------------
Riccardo (Jack) Lucchetti
Dipartimento di Scienze Economiche e Sociali (DiSES)
Università Politecnica delle Marche
(formerly known as Università di Ancona)
r.lucchetti(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------
9:51 a.m.
Am 04.06.2022 um 13:11 schrieb Riccardo (Jack) Lucchetti:
On Fri, 3 Jun 2022, Cottrell, Allin wrote:
> 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.
in Hansl:
<hansl>
R = {}
Q = qrdecomp(mnormal(4)~zeros(4)~mnormal(4) , &R)
scalar r = sumc(abs(diag(R)).>1.0e-12)
print r
</hansl>
OK, thanks, that's what I thought.
I'm noting that in the square-invertible input case the diagonal
elements of R are not necessarily positive with qrdecomp(). Maybe this
could be mentioned in the doc. (The decomposition would be unique if
non-negativity were imposed.)
Next, I'm noting that in the rectangular case m,n with m>n the output Q
is m,n and R is n,n. According to Wikipedia this is the "thin" or
"reduced" QR factorization, which might be mentioned as well. In the
Netlib terminology, this corresponds to A = Q_1 * R, if this source is
preferred. The Matlab "economy-size" convention could be the same thing.
(https://www.mathworks.com/help/matlab/ref/qr.html)
I'm not really sure what the benefit would be from having column
pivoting, which is relatively easy to do in Hansl anyway (note that P is
a permutation matrix, so it can be more compactly expressed as a set of
indices), although of course numerically speaking lapack is unbeatable.
Could you provide us with a use case where pivoting would be useful?
I'm not sure myself yet, let me think about it some more. What the
Matlab people note is that for sparse inputs the result is much sparser
with permutation (see page linked above).
thanks
sven
1:39 p.m.
Am 04.06.2022 um 15:51 schrieb Sven Schreiber:
Am 04.06.2022 um 13:11 schrieb Riccardo (Jack) Lucchetti:
> On Fri, 3 Jun 2022, Cottrell, Allin wrote:
>> 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.
>
> in Hansl:
>
> <hansl>
> R = {}
> Q = qrdecomp(mnormal(4)~zeros(4)~mnormal(4) , &R)
> scalar r = sumc(abs(diag(R)).>1.0e-12)
> print r
> </hansl>
OK, thanks, that's what I thought.
Hm, but plain QR is not rank-revealing in general, that's why the thing
with column-pivoting exists. Trivial example:
<hansl>
R = {}
A = zeros(4,2) ~ I(4,2)
print A
Q = qrdecomp(A , &R)
print rank(A) # correct: 2
scalar r = sumc(abs(diag(R)).>1.0e-12)
print r # wrong: 0
</hansl>
Now, I guess that in your internal usage this is quite unlikely to
happen, but still, if you don't want to use the SVD for that --and I
presume it's because of the computational expenses-- then maybe there's
a role for the pivoted QR thing?
> I'm not really sure what the benefit would be from having
column
> pivoting, which is relatively easy to do in Hansl anyway
AFAICS, it's not enough to wiggle around the columns of R ex-post-facto,
or am I missing something?
I'm not sure myself yet, let me think about it some more. What
the
Matlab people note is that for sparse inputs the result is much sparser
with permutation (see page linked above).
Addendum: Aiming for this sparsity is a deliberate choice, however, and
then the diagonal entries are not decreasing anymore. So it's a
different issue, I guess.
thanks
sven
10:15 a.m.
On Sat, 4 Jun 2022, Sven Schreiber wrote:
Am 04.06.2022 um 15:51 schrieb Sven Schreiber:
> Am 04.06.2022 um 13:11 schrieb Riccardo (Jack) Lucchetti:
>> On Fri, 3 Jun 2022, Cottrell, Allin wrote:
>>> 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.
>>
>> in Hansl:
>>
>> <hansl>
>> R = {}
>> Q = qrdecomp(mnormal(4)~zeros(4)~mnormal(4) , &R)
>> scalar r = sumc(abs(diag(R)).>1.0e-12)
>> print r
>> </hansl>
>
> OK, thanks, that's what I thought.
Hm, but plain QR is not rank-revealing in general, that's why the thing
with column-pivoting exists. Trivial example:
<hansl>
R = {}
A = zeros(4,2) ~ I(4,2)
print A
Q = qrdecomp(A , &R)
print rank(A) # correct: 2
scalar r = sumc(abs(diag(R)).>1.0e-12)
print r # wrong: 0
</hansl>
Now, I guess that in your internal usage this is quite unlikely to
happen, but still, if you don't want to use the SVD for that --and I
presume it's because of the computational expenses-- then maybe there's
a role for the pivoted QR thing?
Very nice counterxample, Sven, kudos for that. Let me think about it.
-------------------------------------------------------
Riccardo (Jack) Lucchetti
Dipartimento di Scienze Economiche e Sociali (DiSES)
Università Politecnica delle Marche
(formerly known as Università di Ancona)
r.lucchetti(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------
10:21 a.m.
On Sun, 5 Jun 2022, Riccardo (Jack) Lucchetti wrote:
Very nice counterxample, Sven, kudos for that. Let me think about it.
In fact this bug also affects the dropcoll() built-in function:
<hansl>
nulldata 4
series z = 0
series x = 1 | zeros($nobs - 1,1)
list X = z x
D = dropcoll(X)
list print D
</hansl>
So yeah, we need to give this proper consideration.
-------------------------------------------------------
Riccardo (Jack) Lucchetti
Dipartimento di Scienze Economiche e Sociali (DiSES)
Università Politecnica delle Marche
(formerly known as Università di Ancona)
r.lucchetti(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------
8:46 p.m.
On Sun, Jun 5, 2022 at 10:21 AM Riccardo (Jack) Lucchetti
<p002264(a)staff.univpm.it> wrote:
On Sun, 5 Jun 2022, Riccardo (Jack) Lucchetti wrote:
> Very nice counterxample, Sven, kudos for that. Let me think about it.
In fact this bug also affects the dropcoll() built-in function:
<hansl>
nulldata 4
series z = 0
series x = 1 | zeros($nobs - 1,1)
list X = z x
D = dropcoll(X)
list print D
</hansl>
So yeah, we need to give this proper consideration.
I tend to think these are different error conditions: (a) the input
contains exactly collinear data, and (b) one or more of the input
series/columns are all zero.
The second case is already caught for series input to most estimators:
we flag an error if any regressor is all zeros, before checking for
perfect collinearity.
I suppose dropcoll() could add a first check for any list element
being all-zero. Though in most cases this would be a waste of time.
Allin
10:48 p.m.
On Sun, Jun 5, 2022 at 8:46 PM Cottrell, Allin <cottrell(a)wfu.edu> wrote:
I suppose dropcoll() could add a first check for any list element
being all-zero. Though in most cases this would be a waste of time.
But also in most cases it won't cost much (since we should reject the
all-zeros condition quite quickly). That's now added in git.
Allin
6:47 a.m.
Am 06.06.2022 um 02:46 schrieb Cottrell, Allin:
On Sun, Jun 5, 2022 at 10:21 AM Riccardo (Jack) Lucchetti
<p002264(a)staff.univpm.it> wrote:
>
> On Sun, 5 Jun 2022, Riccardo (Jack) Lucchetti wrote:
>
>> Very nice counterxample, Sven, kudos for that. Let me think about it.
Full disclosure: I didn't invent it myself, but found it in the depths
of the net.
> In fact this bug also affects the dropcoll() built-in function:
...
> So yeah, we need to give this proper consideration.
I tend to think these are different error conditions: (a) the input
contains exactly collinear data, and (b) one or more of the input
series/columns are all zero.
But that was just an example, not the encompassing definition of the
problem.
Consider this variation without all-zero columns:
<hansl>
R = {}
A = I(4,1) ~ I(4,2)
print A
Q = qrdecomp(A , &R)
print rank(A) # correct: 2
print R
scalar r = sumc(abs(diag(R)).>1.0e-12)
print r # wrong: 1
</hansl>
cheers
sven
11:29 a.m.
On Mon, Jun 6, 2022 at 6:47 AM Sven Schreiber <svetosch(a)gmx.net> wrote:
Am 06.06.2022 um 02:46 schrieb Cottrell, Allin:
>
> I tend to think these are different error conditions: (a) the input
> contains exactly collinear data, and (b) one or more of the input
> series/columns are all zero.
But that was just an example, not the encompassing definition of the
problem.
Consider this variation without all-zero columns:
<hansl>
R = {}
A = I(4,1) ~ I(4,2)
print A
Q = qrdecomp(A , &R)
print rank(A) # correct: 2
print R
scalar r = sumc(abs(diag(R)).>1.0e-12)
print r # wrong: 1
</hansl>
Is there any general characterization of what it takes for plain QR
not to reveal rank?
It would appear that the presence of all-zero rows and/or columns is
necessary but not sufficient. With m > n you can construct mostly-zero
matrices where diag(R) does give the right answer. E.g.
<hansl>
function void rank_check (const matrix X)
R = {}
qrdecomp(X, &R)
printf "rank via SVD: %d\n", rank(X)
printf "rank via QR: %d\n\n", sumc(abs(diag(R)) .> 1.0e-12)
end function
matrix X = zeros(8,4)
X[6,1] = 1
X[7,4] = 1
X[8,4] = 1
print X
rank_check(X)
</hansl>
Allin
12:13 p.m.
Am 06.06.2022 um 17:29 schrieb Cottrell, Allin:
Is there any general characterization of what it takes for plain QR
not to reveal rank?
It would appear that the presence of all-zero rows and/or columns is
necessary but not sufficient. With m > n you can construct mostly-zero
matrices where diag(R) does give the right answer. E.g.
My general answer is: I don't know. There seems to be an elaborate
literature on this, and I haven't dealt with it in the past. However,
I'm kind of skeptical that there are easy rules, because otherwise they
would be expected to become part of the algorithms.
BTW, it is my (current) understanding that also the QR decomp with
pivoting is not as "safe" as SVD, so there seems to be a tradeoff as
well. But I guess this is always true with numerics.
What might be useful is to check out the database of singular matrices
at San Jose State U. (It's astounding what things exist in science!)
They also have a number of fully dense matrices, i.e. 100% non-zeros,
see here:
http://www.math.sjsu.edu/singular/matrices/html/list_by_nnz_percent_desce...
So I guess it would be interesting to run all those matrices through
your simple QR diag check and see what that yields.
thanks
sven
12:55 p.m.
On Mon, Jun 6, 2022 at 12:13 PM Sven Schreiber <svetosch(a)gmx.net> wrote:
Am 06.06.2022 um 17:29 schrieb Cottrell, Allin:
> Is there any general characterization of what it takes for plain QR
> not to reveal rank?
>
> It would appear that the presence of all-zero rows and/or columns is
> necessary but not sufficient. With m > n you can construct mostly-zero
> matrices where diag(R) does give the right answer. E.g.
My general answer is: I don't know. There seems to be an elaborate
literature on this, and I haven't dealt with it in the past. However,
I'm kind of skeptical that there are easy rules, because otherwise they
would be expected to become part of the algorithms.
OK, thanks.
What might be useful is to check out the database of singular
matrices
at San Jose State U. (It's astounding what things exist in science!)
I'll take a look.
Meanwhile (in git), you can now call QR with pivoting, the trigger
being a non-null third (matrix-pointer) argument. That gets you the
ordering of the columns of the input, which can be compared with
seq(1,n) to see if anything has moved. This is available for complex
matrices too.
R = {}
P = {}
Q = qrdecomp(A, &R, &P)
Allin
8:18 a.m.
Am 06.06.2022 um 18:55 schrieb Cottrell, Allin:
Meanwhile (in git), you can now call QR with pivoting, the trigger
being a non-null third (matrix-pointer) argument.
For m < n (more rows) this seems to crash gretl. (I don't need it, but I
kind of thought it might cause trouble.)
cheers
sven
11:55 a.m.
On Tue, 7 Jun 2022, Sven Schreiber wrote:
Am 06.06.2022 um 18:55 schrieb Cottrell, Allin:
>
> Meanwhile (in git), you can now call QR with pivoting, the trigger
> being a non-null third (matrix-pointer) argument.
For m < n (more rows) this seems to crash gretl. (I don't need it, but I
kind of thought it might cause trouble.)
There was a guard against that in place, but maybe it didn't come
soon enough. Now fixed in git and snapshots.
BTW, the qrdecomp() doc is now more explicit, and includes the
requirement that m >= n. (I gather that for short-and-fat input you
really want the LQ decomposition rather than QR, but I don't see a
need to implement that yet.)
Allin
4:33 p.m.
On Mon, Jun 6, 2022 at 12:13 PM Sven Schreiber <svetosch(a)gmx.net> wrote:
What might be useful is to check out the database of singular matrices
at San Jose State U. (It's astounding what things exist in science!)
Most of the matrices in that database are square. From the point of
view of libgretl I'm particularly interested in the case m > n, since
our internal usage of QR is mostly confined to such cases (matrices of
regressors or instruments, where we'd generally expect m >> n). The
"NYPA" group in the database contains such matrices. There's one case
of full "structural rank" and the other matrices are rank-deficient.
All of them have "numerical rank" (assessed via SVD) < "structural
rank".
Here's what I found: gretl's (SVD) rank() agrees with the numerical
rank values shown in the database. So do the results of counting
abs(R[i,i]) values from QR, provided the minimal value is in the range
1.0e-10 to 1.0e-8. Plain QR and QR with column pivoting produce the
same "revealed rank" in all cases. The illustrative min(R) of 1.0e-12
used upthread seems a bit too small: it gives a numerical rank
slightly greater than the "official" value for 3 of the 5 specimen
matrices. (Internally, we use R_DIAG_MIN = 1.0e-8.)
Allin
5:23 p.m.
On 06.06.22 22:33, Cottrell, Allin wrote:
Here's what I found: gretl's (SVD) rank() agrees with the
numerical
rank values shown in the database. So do the results of counting
abs(R[i,i]) values from QR, provided the minimal value is in the range
1.0e-10 to 1.0e-8. Plain QR and QR with column pivoting produce the
same "revealed rank" in all cases. The illustrative min(R) of 1.0e-12
used upthread seems a bit too small: it gives a numerical rank
slightly greater than the "official" value for 3 of the 5 specimen
matrices. (Internally, we use R_DIAG_MIN = 1.0e-8.)
Thanks, Allin, that's good to know. But the (repeated) unit vector case
does not seem totally irrelevant for econometrics and the dropcoll
function, thinking of impulse dummies for example. Maybe your hunch
provides a useful heuristic: only use column-pivoting if there's a row
of zeros. (Plus the precheck for all-zeros.)
OTOH, I ran some very crude and simple speed comparisons between
qrdecomp with and without pivoting (on Linux, current git), and
depending on the input the advantage could go either way it seemed,
which I found a little surprising. But then I also ran SVD on the same
input (and also grabbing the optional output), and that went even
faster, so now I'm puzzled... isn't SVD supposed to be the most
expensive thing of the candidates? What am I missing?
thanks
sven
5:56 p.m.
On Mon, Jun 6, 2022 at 5:24 PM Sven Schreiber <svetosch(a)gmx.net> wrote:
OTOH, I ran some very crude and simple speed comparisons between
qrdecomp with and without pivoting (on Linux, current git), and
depending on the input the advantage could go either way it seemed,
which I found a little surprising. But then I also ran SVD on the same
input (and also grabbing the optional output), and that went even
faster, so now I'm puzzled... isn't SVD supposed to be the most
expensive thing of the candidates? What am I missing?
Are you doing enough replications to get good resolution from the timer?
Here's what I see for a 100 x 10 random matrix with 10,000 replications:
SVD 0.324311s
QR 0.094924s
QRP 0.114874s
But as I mentioned earlier, big m favors SVD. Here's 800 x 20 with
1000 replications.
SVD 0.137257s
QR 0.231528s
QRP 0.248974s
Allin
3:23 a.m.
Am 06.06.2022 um 23:56 schrieb Cottrell, Allin:
On Mon, Jun 6, 2022 at 5:24 PM Sven Schreiber
<svetosch(a)gmx.net> wrote:
>
> OTOH, I ran some very crude and simple speed comparisons between
> qrdecomp with and without pivoting (on Linux, current git), and
> depending on the input the advantage could go either way it seemed,
> which I found a little surprising. But then I also ran SVD on the same
> input (and also grabbing the optional output), and that went even
> faster, so now I'm puzzled... isn't SVD supposed to be the most
> expensive thing of the candidates? What am I missing?
Are you doing enough replications to get good resolution from the timer?
Here's what I see for a 100 x 10 random matrix with 10,000 replications:
SVD 0.324311s
QR 0.094924s
QRP 0.114874s
I think I had 500x12 or so. But interesting that QRP is almost as fast
as QR.
But as I mentioned earlier, big m favors SVD. Here's 800 x 20
with
1000 replications.
SVD 0.137257s
QR 0.231528s
QRP 0.248974s
Apparently I missed that, this is something I was completely unaware of.
I thought that one tries to avoid SVD because of the costs. But here I
see no costs!? For mols() for example, SVD has to be forced explicitly.
OK, you could say that at least in macro T=100 is more common than
T=800. But maybe then a heuristic could help, so SVD is used for some
large enough T (or T-K, or T/K). (Only as a fallback from Cholesky, of
course.)
For the dropcoll() functions and related ones that currently use the R
diagonal internally, I guess a similar argument applies. The cost of QRP
seems negligible, and SVD takes over at some point. The advantage of
plain QR seems to be limited.
thanks
sven
894
days inactive
898
days old
20 comments
4 participants
participants (4)
-
Allin Cottrell -
Cottrell, Allin
-
Riccardo (Jack) Lucchetti -
Sven Schreiber