Re: [eigen] RFC: BasicSVD class |

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

*To*: eigen@xxxxxxxxxxxxxxxxxxx*Subject*: Re: [eigen] RFC: BasicSVD class*From*: Benoit Jacob <jacob.benoit.1@xxxxxxxxx>*Date*: Mon, 10 Dec 2012 11:18:23 -0500*Dkim-signature*: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=20120113; h=mime-version:in-reply-to:references:date:message-id:subject:from:to :content-type; bh=JXQpuMkejO555v7RfUn3KFvLKhCGzifC9RwI4j5DelU=; b=y0PjffxsSAKVgkykYnyyjRgTdg6TpIdnqTsixnguHLFeztwad4JI1LTI5gNdudvjew Xu5QZAOsJ2FZ6u/DUdOHRHW2Sq3NZmiPVIKiPswTRy2bNzltw3RUdKcsMd8KCVRCs5TP kI+3B64EUzv7gVS6o96uHeQU2Q4vKKHwuWcPVMEM/xckf+oCIwQmsMcFgtA1K6Y5fY/I kyI0fAoMpK/Xdy/zy6a9zwLW3lxb1lSCdX3BgGxFq6UWmzSIfx9bCck/rPE6JIzIyZQ2 iP5ju1Oy8Kaxwu3RTG8uRD5qHa95paOCApZCojVfFxFTmgfEhcDknzL0tiXieGdBJ2UE siaA==

I've drifted far away from the project now, but I just want to say that I'm very happy that people are taking care of our SVD situation. Clearly JacobiSVD leaves much to be desired in terms of speed. There is room for a fast-but-less-reliable SVD as I understand that BasicSVD is. I would just suggest a more specific name than "BasicSVD": use a name that is somehow specific to the particular method that your SVD employs.

Do you know how BasicSVD compares to Householder-bidiagonalization SVD's such as the LAPACK methods you mention? Those are what I thought we wanted (http://eigen.tuxfamily.org/bz/show_bug.cgi?id=67 ); but it's been so long now, I'd rather have another thing than nothing at all.

Benoit

2012/12/10 Mark Borgerding <mark@xxxxxxxxxxxxxx>

Thanks for your assessment of "not useless" ;). To specify your statement below: the algorithm struggleswith getting the smallest singular valuesfrom ill-conditioned matrices. This is the harder problem, to be sure, than finding the dominant dimensions. It is impressive that JacobiSVD can do such a great job with pulling out dimensions that are 10 orders of magnitude down.

However ...

Getting theweakestdimensions is not whatIneed an SVD for. I want the first few dominant dimensions. I suspect this is a common requirement (hence matlab eigs and svds functions). I have been frustated in the past by attempts to use JacobiSVD because of thehundred-fold speed differencecompared to the LAPACK functions ( sgesdd, cgesdd, dgesdd, zgesdd )

I'm not suggesting BasicSVD replace JacobiSVD. I'm offering it as an alternative.

On the subject of power iteration. Note that what I implemented is not textbook power iteration. I use feedback guided by the maximization of the Rayleigh Quotient. If you want to disable the feedback, comment out the lines with "damp" in them.

On 12/10/2012 09:12 AM, Jitse Niesen wrote:On Fri, 7 Dec 2012, Gael Guennebaud wrote:

Indeed, this simple algorithm seems pretty fast and it could even be even

faster once we have fully optimized SelfAdjointEigenSolver. However

I'm skeptical about its accuracy since it is squaring the entries that is

typically something you want to avoid when using a SVD solver.

The algorithm indeed seems to struggle with ill-conditioned matrices, which have both small and large singular values. The test program below uses a 6-by-6 matrix with doubles and gives the following output:

-------------

Singular values (exact):

1 0.01 0.0001 1e-06 1e-08 1e-10

Singular values (with JacobiSVD):

1 0.01 0.0001 1e-06 1e-08 1e-10

Singular values (with BasicSVD):

1 0.01 0.0001 9.99981e-07 1.05497e-08 1.43198e-09

-------------

Of course this does not mean that it's useless; there are probably circumstances in which BasicSVD is good enough as Mark indicates, and it would make a nice addition. The power method by itself would already be very nice.

Jitse

-----------------------------------

#include "BasicSVD.h"

#include <iostream>

using namespace Eigen;

int main()

{

// The singular values

const int size = 6;

VectorXd singvals(size);

singvals << 1, 1e-2, 1e-4, 1e-6, 1e-8, 1e-10;

// Generate two random orthogonal matrices

MatrixXd Uexact = MatrixXd::Random(size, size).householderQr().householderQ();

MatrixXd Vexact = MatrixXd::Random(size, size).householderQr().householderQ();

// Generate matrix whose SVD is to be computed

MatrixXd A = Uexact * singvals.asDiagonal() * Vexact;

std::cout << "Singular values (exact):\n";

std::cout << singvals.transpose() << "\n\n";

// Compute SVD using standard algorithm

VectorXd singvalsJac = A.jacobiSvd().singularValues();

std::cout << "Singular values (with JacobiSVD):\n";

std::cout << singvalsJac.transpose() << "\n\n";

// Compute SVD using new algorithm

BasicSVD<MatrixXd> svdBasic(A, size);

VectorXd singvalsBas = svdBasic.singularValues();

std::cout << "Singular values (with BasicSVD):\n";

std::cout << singvalsBas.transpose() << "\n";

}

**References**:**[eigen] RFC: BasicSVD class***From:*Mark Borgerding

**Re: [eigen] RFC: BasicSVD class***From:*Mark Borgerding

**Re: [eigen] RFC: BasicSVD class***From:*Gael Guennebaud

**Re: [eigen] RFC: BasicSVD class***From:*Jitse Niesen

**Re: [eigen] RFC: BasicSVD class***From:*Mark Borgerding

**Messages sorted by:**[ date | thread ]- Prev by Date:
**Re: [eigen] RFC: BasicSVD class** - Next by Date:
**[eigen] Transform class performance and inconsistencies** - Previous by thread:
**Re: [eigen] RFC: BasicSVD class** - Next by thread:
**[eigen] Displaying an Eigen object during debugging?**

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