Re: [eigen] Specialized QR

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

Thanks for your email, I haven't had time yet to look at your code but
I can already answer this:

2009/5/13, Andrea Arteaga <yo.eres@xxxxxxxxx>:
> The result is: the QR_Givens class computes the QR
> decomposition about 4 times faster than the QR class.

For what sizes?

> The relative errors
> are
> sometimes bigger than those of QR, but I never had errors bigger than 1e-13.
> (I tested if Q^t * Q is the identity and if Q*R = A.)

again, this makes more sense if you say which size. But yes, 1e-13 sounds nice.

>
> Furthermore, I noticed that Eigen QR class has (run-time) troubles with mxn
> matrices, where m<n (I become this error:
> "../eigen2/Eigen/src/Core/MapBase.h:166:
> Eigen::MapBase<Derived>::MapBase(const
> typename Eigen::ei_traits<T>::Scalar*, int, int) [with Derived =
> Eigen::Block<Eigen::Block<Eigen::Matrix<double, 50, 100, 0, 50, 100>, 50, 1,

Can you send us code reproducing the problem?
Also, 50 and 100 are HUGE as compile-time fixed sizes. For matrices of
this size, you should definitely use dynamic size.

> In this moment my class works only with fixed-sized matrices.

ah ok.

> As I have
> time, I
> will implement the dynamically-sized version and try with larger matrices.

Yes, that would be nice. This shouldn't be any more complicated and
aside from possible meta-unrollers for fixed-size, both cases should
be exactly the same code. You can have a look how Eigen's QR or LU
does, to  handle both cases at once. It's entirely implicit.

> I also implemented a rank() method, but this is not very stable yet.

See the other thread on LU precision tuning: QR is basically not rank
revealing unless you add pivoting to it and even then i don't know
exactly what it takes to be really rank revealing but Hauke found
references to read on this topic.

For a fixed-size QR, i'd say: nevermind rank-revealing QR, just focus
on performance and no rank().

> Another question: the QR class returns a quadratic matrix R and a
> rectangular
> matrix Q; QR_Givens instead returns a rectangular R and a "really
> orthogonal"
> quadratic matrix Q.

Very good catch. See the other thread (as above), we have also been
discussing that issue.

> Which behaviour is better?

I would like us to standardize everywhere in Eigen on saying that,
once and for all, "unitary" implies "square", so, unless otherwise
explicitly stated (e.g. in a ThinQR class), all "Q" matrices are
square.
So: R rectangular and Q "really unitary".
Notice that i say unitary instead of orthogonal because i want support
for complex matrices.

Cheers
Benoit