eigenvectors obtained with the example function not expected

[hellow...@yeah.net]

Thanks for the great package. I am testing it with a generalised eigenvalue problem using the provided example function:
[Xsol, Ssol] = generalized_eigenvalue_computation(A, B, p)
I got totally different eigenvectors from those I got with eigs from matlab. Anyone can explain this? what are the relations between those different results? Any comments will be appreciated.

[hellow...@yeah.net]

When I set the number of eigenvalues as 1 which means the largest/1st-dominant will be calculated, the results from two methods are exactly the same, but when the number is set greater than 1, the results become different and I cannot figure out any relations.

[Nicolas Boumal]

Hello,

You are right, thank you for pointing this out!

I tried this code:

% Generate test data

A = randn(5); A = A*A';

B = randn(5); B = B*B';

% How many eigenvalue/eigenvector pairs must be computed

p = 2;

% Check Matlab's built-in result

[V1, D1] = eigs(A, B, p)

(A*V1) ./ (B*V1)            % we get constant columns corresponding to the generalized eigenvalues indeed

% Check Manopt's example

[V2, D2] = generalized_eigenvalue_computation(A, B, p) % we get the same eigenvalues

(A*V2) ./ (B*V2)  % NOT constant columns

To make up for that, in the example file, please replace everything that occurs after the call to trustregions by the following code:

    % To extract the eigenvalues, solve the small p-by-p symmetric 

    % eigenvalue problem.

    [Vsol, Dsol] = eig(Xsol'*(A*Xsol));

    Ssol = diag(Dsol);

    

    % To extract the eigenvectors, rotate Xsol by the p-by-p orthogonal

    % matrix Vsol.

    Xsol = Xsol*Vsol;

 

This will be modified in the next release (which we hope will be this summer.)

Best, and thanks again for letting us know!

Nicolas