[eigen] QR precision tuning

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

Hi,

This is mostly so that these results get stored somewhere one can find
them back ;) don't feel forced to read.

i'm doing rank-revealing QR these days. There were some questions to answer:
- what does full pivoting bring as compared to column pivoting?
theoretically it is required to make householder transformations
stable.
There is also this paper: http://citeseer.ist.psu.edu/old/518014.html
- since unitary transformations preserve the L2 norm, can we use the
L2 norm of columns to determine the pivots, avoiding to recompute that
every time, making pivoting much less costly? Or is that too
optimistic (not accurate enough)?
- what cutoff should be used? The same formula epsilon*size as in LU and LLT ?

To answer this, i measured the accuracy of full-pivoting QR, on the
operations of solving and of rank determination (on random matrices
with 50% singular values 0 and 50% singular values 1), using
- 2 pivot strategies: "traditional" (at each step look for the biggest
coeff) and "optimistic" (as described above)
- 2 cutoff formulas: "epsilon*size" and "epsilon".

Results:

First, with no pivoting at all (and hence no cutoff):

*** Matrix3f ***
 repeat = 2000000, size = 3
time:         1.37
solve error: 1.82e-07
*** Matrix4d ***
 repeat = 2000000, size = 4
time:         2.54
solve error: 3.67e-15
*** MatrixXf ***
 repeat = 1000, size = 100
time:         2.4
solve error: 1.81e-06
*** MatrixXd ***
 repeat = 1000, size = 100
time:         3.92
solve error: 2.98e-14
*** MatrixXf ***
 repeat = 1, size = 1000
time:         2.5
solve error: 0.000114
*** MatrixXd ***
 repeat = 1, size = 1000
time:         4.4
solve error: 3.76e-14


Now with full pivoting, with cutoff = epsilon*size:

traditional pivoting...
*** Matrix3f ***
 repeat = 2000000, size = 3
time:         1.44
rank error:   0
solve error: 1.05e-07
*** Matrix4d ***
 repeat = 2000000, size = 4
time:         2.07
rank error:   0
solve error: 2.86e-16
*** MatrixXf ***
 repeat = 1000, size = 100
time:         3.05
rank error:   0
solve error: 4.4e-07
*** MatrixXd ***
 repeat = 1000, size = 100
time:         4.42
rank error:   0
solve error: 8.01e-16
*** MatrixXf ***
 repeat = 1, size = 1000
time:         3.08
rank error:   0
solve error: 2.44e-06
*** MatrixXd ***
 repeat = 1, size = 1000
time:         5.01
rank error:   0
solve error: 5.02e-15

optimistic pivoting...
*** Matrix3f ***
 repeat = 2000000, size = 3
time:         1.5
rank error:   0
solve error: 5.12e-07
*** Matrix4d ***
 repeat = 2000000, size = 4
time:         2.04
rank error:   0
solve error: 2.1e-16
*** MatrixXf ***
 repeat = 1000, size = 100
time:         2.16
rank error:   0
solve error: 6.45e-07
*** MatrixXd ***
 repeat = 1000, size = 100
time:         3.51
rank error:   0
solve error: 3.43e-13
*** MatrixXf ***
 repeat = 1, size = 1000
time:         2.16
rank error:   0
solve error: 0.000813
*** MatrixXd ***
 repeat = 1, size = 1000
time:         3.79
rank error:   0
solve error: 3.35e-14


finally with full pivoting and cutoff=epsilon. Notice how this cutoff
is inappropriate for rank determination; but i still tried it to check
what solve() accuracy it gives.


traditional pivoting...
*** Matrix3f ***
 repeat = 2000000, size = 3
time:         1.46
rank error:   0
solve error: 1.05e-07
*** Matrix4d ***
 repeat = 2000000, size = 4
time:         2.07
rank error:   0
solve error: 2.86e-16
*** MatrixXf ***
 repeat = 1000, size = 100
time:         3.54
rank error:   44
solve error: 7.02e-07
*** MatrixXd ***
 repeat = 1000, size = 100
time:         5.09
rank error:   40
solve error: 1.14e-15
*** MatrixXf ***
 repeat = 1, size = 1000
time:         3.52
rank error:   497
solve error: 8.28e-06
*** MatrixXd ***
 repeat = 1, size = 1000
time:         5.63
rank error:   499
solve error: 2.12e-14

optimistic pivoting...
*** Matrix3f ***
 repeat = 2000000, size = 3
time:         1.5
rank error:   0
solve error: 5.12e-07
*** Matrix4d ***
 repeat = 2000000, size = 4
time:         2
rank error:   0
solve error: 2.1e-16
*** MatrixXf ***
 repeat = 1000, size = 100
time:         2.53
rank error:   45
solve error: 9.76e-07
*** MatrixXd ***
 repeat = 1000, size = 100
time:         4.06
rank error:   47
solve error: 1.99e-15
*** MatrixXf ***
 repeat = 1, size = 1000
time:         2.44
rank error:   499
solve error: 1.91e-05
*** MatrixXd ***
 repeat = 1, size = 1000
time:         4.29
rank error:   492
solve error: 1.54e-14



Conclusions:
 - when optimal accuracy is wanted, nothing equals traditional
pivoting. Then the cutoff epsilon*size is good and solve precision
isn't significantly enhanced by using a smaller cutoff.
- when speed is wanted, optimistic pivoting is better; using a smaller
cutoff increases precision but that's not worth the hassle (losing
rank determination) since if the goal is precision then traditional
pivoting is better anyway.

So I'll implement full-pivoting QR using traditional pivoting
(actually that's already in Hg), and  column-pivoting QR using
optimistic pivoting. In all cases, using the cutoff epsilon*size.

Cheers,
Benoit