[eigen] decompositions: how to deal with unitaries

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

Hi,

this mail is to:
- tell the list about some new stuff in the devel branch
- get some feedback if i'm getting something wrong

*** New stuff ***

By no means a comprehensive overview of recent developments --- just
some picks relevant to what I want to discuss.

- new Householder module:

What's done: provide general, reusable code for Householder transformations
To do: block Householder transformations

We've also been discussing with Gael that we'd like to be able to
treat a sequence of Householder transformations (typically stored in a
triangular part of a matrix) as a "special matrix" that can be used as
an ordinary matrix except that it takes advantage of the packed
Householder storage, e.g. operator* would apply the sequence of
householder transformations, etc.

- new Jacobi module:

What's done: provide general, reusable code for "planar rotations"
---> Givens and Jacobi.
To do: find a better name for the module, perhaps PlanarRotations... or not

No need to worry about "block versions" here, block Jacobi is done at
the level of the algorithm using it, the Jacobi rotations stay the
same.

- block version

- QR module:

What's done:
- HouseholderQR "modernized" after recent mailing list discussions
- new ColPivotingHouseholderQR and FullPivotingHouseholderQR allow
rank determination and are more numerically stable.

To do:

- block versions (requires block householder transformations)
- currently matrixQ() returns a square matrix. More on that below.

- SVD module:

What's done:

New JacobiSVD class, computes the SVD by Jacobi rotations.
Theoretically slower than "normal" SVD but:
- has guaranteed precision (inherent in the algorithm)
- is faster for very small sizes
Moreover, compared to our existing SVD class, it has the advantage of
supporting all recangular sizes (also m<n) and supporting complex
numbers, and allowing the user to skip the computation of U or V or
both. And it is very fast for rectangular sizes (n<<m or m<<n) as it
does a full-pivoting QR to reduce to the smaller dimension.

To do:

- consider block algorithm
- consider research paper, www.netlib.org/lapack/lawnspdf/lawn170.pdf
- currently U and V, if required, are square. More on that below.

*** Stuff open for discussion ***

What I want to discuss is how all these unitaries should be handled.

Here there's a fundamental difference between Householder algorithms
and Jacobi/Givens algorithms.

Householder algorithms just store the sequence of Householder vectors,
from which the unitary can be retrieved later.

For Jacobi it's not a good idea, since the number of Jacobi rotations
is big, and isn't even known beforehand.

For Givens .... i still don't know. It not a good idea for large
sizes, for the reasons discussed in the GivensQR thread, but GivensQR
is mostly useful for small sizes... will write in the other thread.

So anyway here's the plan:
- for HouseholderQR: as said above, with Gael we have a plan for
HouseholderSequence special matrix. So the HouseholderQR classes would
have a matrixQ() method that returns a HouseholderSequence, and then
doing square_matrix = matrixQ() evaluates the whole Q, and it just
occured to me that we could have logic to make
small_rectangular_matrix = matrixQ() evaluate only the "economy Q" (as
HouseholderQR used to do until the recent changes).

- for JacobiSVD: here the unitaries must be computed simultaneously
with the decomposition. So the options concerning them need to be
options of JacobiSVD itself. Currently we only have an option to skip
them (SkipU, SkipV) but if they are computed then it's as a full
unitary matrix. We need an option too for computing only the "compact
version" that's all what's needed in most cases --- U and V as small
rectangular matrices. Since the user will never need both versions at
once (anyway he could get the compact version from the full version)
that should be an option of the SVD decomposition class. So: CompactU
and CompactV. Concretely:

what you can already do:

JacobiSVD<MatrixXd,SkipU|SkipV> svd(matrix);

what you'll be able to do:

JacobiSVD<MatrixXd,CompactU|CompactV> svd(matrix);

Cheers,
Benoit