[eigen] decompositions: how to deal with unitaries |

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]

*To*: eigen <eigen@xxxxxxxxxxxxxxxxxxx>*Subject*: [eigen] decompositions: how to deal with unitaries*From*: Benoit Jacob <jacob.benoit.1@xxxxxxxxx>*Date*: Thu, 3 Sep 2009 10:35:03 -0400*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=gamma; h=domainkey-signature:mime-version:received:date:message-id:subject :from:to:content-type; bh=l9YYpJtGnKCrIjMVSaZIeu+g4Ms0RswEIkcUTMp7DEQ=; b=srXUShNRTc2a5dIoAXG/kY2+sRAtJ/TAO0KtDgcPe5QEtgeTspFMYXnmF3aY65/km9 F7nW/KD/ewnDse7CiRnfyQX2dPSpMC895Dhr1nKCSTuq9sXe+QFw2BikrDe7iEZAVLWg GF5XAW3gW3gwXJmWPE7wdpkxGBy18mFEEmufc=*Domainkey-signature*: a=rsa-sha1; c=nofws; d=gmail.com; s=gamma; h=mime-version:date:message-id:subject:from:to:content-type; b=pyniNm3jLYhaxaZAZDY2bW2tCnVmrC4yfl+CMOrAMomn/oKN0tyrtlczt7pB0FyyXD YGgHXfHeklZ0B415h+/u8UHs1mGAtKyu4dcSgTEnBSRVDrgeQWCGPKJdL81Hlcb25WqN OM8euL2EiXyezUnZZXAW/Wp6KLZXrMMtbY6T8=

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

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] Cost of a map operation** - Next by Date:
**Re: [eigen] GivensQR** - Previous by thread:
**Re: [eigen] Cost of a map operation** - Next by thread:
**Re: [eigen] GivensQR**

Mail converted by MHonArc 2.6.19+ | http://listengine.tuxfamily.org/ |