Re: [eigen] Specialized QR

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

Hi,

2009/4/28 Andrea Arteaga <yo.eres@xxxxxxxxx>:
> I did some work on QR module to do some specialized version for small
> matrices.
> My idea is to introduce a class qr_compute, which does all the work (QR
> decomposition, matrixQ, matrixR, rank, solve,...) whit respect to the
> choosen algorithm (Householder, Gram-Schmidt or Givens). The QR class were
> only a front-end to this class. In this way, we don't need to specialize the
> QR, but only the backend qr_compute class.
>
> [......]
>
> The QR class should have a private member, e.g.:
>
> qr_compute<MatrixType, qr_decide<MatrixType::RowsAtCompileTime,
> MatrixType::ColsAtCompileTime>::Algorithm > compute_;

Don't worry, we can handle that meta-programming part. What's really
useful is your work below:

> The logic to decide which algorithm to use were done by a qr_decide class in
> this way:
>
> template<int rows, int cols> struct qr_decide {
> enum { Algorithm = rows < 3 ? 1 : ((rows < 15 && cols < 15) ? 2 : 3) };
> };
>
> Number 1 is Givens (for matrices with only 2 rows), 2 Gram-Schmidt (3x2
> upper to 15x15), 3 Householder.

We don't need specializations for sizes bigger than 8x8. So I'd say
that you can skip case 3.

So what's really interesting to me here is that you have identified
that Givens was best for 2 rows and that Gram-Schmidt was best in all
other cases in which we want to specialize, and that actually
fixed-size specializations /never/ should use Householder.

I think the big issue to discuss now is how the QR decomposition
should be stored.

The generic QR we have does Householder. This allows for compact
storage of the QR decomposition in a single matrix, the householder
vectors being stored in a triangular part, LAPACK-compatible. But this
is specific to Householder: when we store the QR in this way, we
hardcode the fact that it was computed using Householder reflections.

This means that your fixed-size specializations will have to use
another kind of storage, incompatible with this one. This means a
different API, as in Householder QR the QR matrix should be exposed in
the API (for interoperability with other libs like LAPACK, and because
btw for Eigen 2.1 i'm suggesting we remove all these methods returning
Part's and let the user do that himself -- so I'm saying here we
should remove the matrixR method in QR and have a matrixQR instead in
Eigen 2.1).

This means that actually, the fixed-size specializations are so
different from Householder QR that they probably should be a separate
class!

I really didn't see this coming, it all comes from your finding that
fixed-size QR should not use Householder. (And I never looked closely
at QR before, i didn't contribute to it in Eigen 2.0).

Now this makes it all look nicer. This TinyQR class can assume right
away that the Q and the R are stored in two distinct matrices.
Sub-optimal memory usage and locality is not an issue as this is only
for small matrices. So the only distinction to make is whether to use
Givens or Gram-Schmidt, but that only affects the algorithm, not the
API, and not the other methods.

You can then do a meta-selector (like you suggested) but only for the
compute() method: that's far nicer than having to do a big
meta-selector class _compute specializing half of the QR class!

Does this make sense? as I said i never look closely at QR.

Cheers,
Benoit