Home > examples > generalized_eigenvalue_computation.m

# generalized_eigenvalue_computation

## PURPOSE

Returns orthonormal basis of the dominant invariant p-subspace of B^-1 A.

## SYNOPSIS

function [Xsol, Ssol] = generalized_eigenvalue_computation(A, B, p)

## DESCRIPTION

``` Returns orthonormal basis of the dominant invariant p-subspace of B^-1 A.

function [Xsol, Ssol] = generalized_eigenvalue_computation(A, B, p)

Input: A is a real, symmetric matrix of size nxn,
B is a symmetric positive definite matrix, same size as A
p is an integer such that p <= n.

Output: Xsol: a real, B-orthonormal matrix X of size nxp such that
trace(X'*A*X) is maximized, subject to X'*B*X = identity.
That is, the columns of X form a B-orthonormal basis of a
dominant subspace of dimension p of B^(-1)*A. These are thus
generalized eigenvectors associated with the largest generalized
eigenvalues of B^(-1)*A  (in no particular order). Sign is
important: 2 is deemed a larger eigenvalue than -5.
Ssol: the eigenvalues associated with the eigenvectors Xsol, in a
vector.

We intend to solve the homogeneous system A*X = B*X*S,
where S is a diagonal matrix of dominant eigenvalues of B^-1 A.

The optimization is performed on the generalized Grassmann manifold,
since only the space spanned by the columns of X matters in the
optimization problem.

The optimization problem that we are solving here is
maximize trace(X'*A*X) subject to X'*B*X = eye(p).
Consequently, the solutions remain invariant to transformation
X --> XQ, where Q is a p-by-p orthogonal matrix. The search space, in
essence, is set of equivalence classes
[X] = {XQ : X'*B*X = I and Q is orthogonal matrix}. This space is called
the generalized Grassmann manifold.
Before returning, Q is chosen such that Xsol = Xq matches the output one
would expect from eigs.

## CROSS-REFERENCE INFORMATION

This function calls:
• grassmanngeneralizedfactory Returns a manifold struct of "scaled" vector subspaces.
• norm NORM Norm of a TT/MPS tensor.
• norm NORM Norm of a TT/MPS block-mu tensor.
• trustregions Riemannian trust-regions solver for optimization on manifolds.
This function is called by:

## SOURCE CODE

```0001 function [Xsol, Ssol] = generalized_eigenvalue_computation(A, B, p)
0002 % Returns orthonormal basis of the dominant invariant p-subspace of B^-1 A.
0003 %
0004 % function [Xsol, Ssol] = generalized_eigenvalue_computation(A, B, p)
0005 %
0006 % Input: A is a real, symmetric matrix of size nxn,
0007 %        B is a symmetric positive definite matrix, same size as A
0008 %        p is an integer such that p <= n.
0009 %
0010 % Output: Xsol: a real, B-orthonormal matrix X of size nxp such that
0011 %         trace(X'*A*X) is maximized, subject to X'*B*X = identity.
0012 %         That is, the columns of X form a B-orthonormal basis of a
0013 %         dominant subspace of dimension p of B^(-1)*A. These are thus
0014 %         generalized eigenvectors associated with the largest generalized
0015 %         eigenvalues of B^(-1)*A  (in no particular order). Sign is
0016 %         important: 2 is deemed a larger eigenvalue than -5.
0017 %         Ssol: the eigenvalues associated with the eigenvectors Xsol, in a
0018 %         vector.
0019 %
0020 % We intend to solve the homogeneous system A*X = B*X*S,
0021 % where S is a diagonal matrix of dominant eigenvalues of B^-1 A.
0022 %
0023 %
0024 % The optimization is performed on the generalized Grassmann manifold,
0025 % since only the space spanned by the columns of X matters in the
0026 % optimization problem.
0027 %
0028 % The optimization problem that we are solving here is
0029 % maximize trace(X'*A*X) subject to X'*B*X = eye(p).
0030 % Consequently, the solutions remain invariant to transformation
0031 % X --> XQ, where Q is a p-by-p orthogonal matrix. The search space, in
0032 % essence, is set of equivalence classes
0033 % [X] = {XQ : X'*B*X = I and Q is orthogonal matrix}. This space is called
0034 % the generalized Grassmann manifold.
0035 % Before returning, Q is chosen such that Xsol = Xq matches the output one
0036 % would expect from eigs.
0037 %
0039
0040
0041 % This file is part of Manopt and is copyrighted. See the license file.
0042 %
0043 % Main author: Bamdev Mishra, June 30, 2015.
0044 % Contributors:
0045 % Change log:
0046 %
0047 %     Aug. 10, 2016 (NB): the eigenvectors Xsol are now rotated by Vsol
0048 %     before they are returned, to ensure the output matches what you would
0049 %     normally expect calling eigs.
0051
0052     % Generate some random data to test the function
0053     if ~exist('A', 'var') || isempty(A)
0054         n = 128;
0055         A = randn(n);
0056         A = (A+A')/2;
0057     end
0058     if ~exist('B', 'var') || isempty(B)
0059         n = size(A, 1);
0060         e = ones(n, 1);
0061         B = spdiags([-e 2*e -e], -1:1, n, n); % Symmetric positive definite
0062     end
0063
0064     if ~exist('p', 'var') || isempty(p)
0065         p = 3;
0066     end
0067
0068     % Make sure the input matrix is square and symmetric
0069     n = size(A, 1);
0070     assert(isreal(A), 'A must be real.')
0071     assert(size(A, 2) == n, 'A must be square.');
0072     assert(norm(A-A', 'fro') < n*eps, 'A must be symmetric.');
0073     assert(p <= n, 'p must be smaller than n.');
0074
0075     % Define the cost and its derivatives on the generalized
0076     % Grassmann manifold, i.e., the column space of all X such that
0077     % X'*B*X is identity.
0078     gGr = grassmanngeneralizedfactory(n, p, B);
0079
0080     problem.M = gGr;
0081     problem.cost  = @(X)    -trace(X'*A*X);
0083     problem.ehess = @(X, H) -2*(A*H); % Only Euclidean Hessian needed.
0084
0085     % An alternative way to compute the egrad and the ehess is to use
0086     % automatic differentiation provided in the deep learning toolbox (slower)
0087     % Notice that the function trace is not supported for AD so far.
0088     % Replace it with ctrace described in the file manoptADhelp.m
0089     % problem.cost = @(X)    -.5*ctrace(X'*A*X);
0090     % call manoptAD to automatically obtain the egrad and the ehess
0092
0093     % Execute some checks on the derivatives for early debugging.
0094     % These things can be commented out of course.
0096     % pause;
0097     % checkhessian(problem);
0098     % pause;
0099
0100     % Issue a call to a solver. A random initial guess will be chosen and
0101     % default options are selected except for the ones we specify here.
0102     options.Delta_bar = 8*sqrt(p);
0104     options.verbosity = 2; % set to 0 to silence the solver, 2 for normal output.
0105     [Xsol, costXsol, info] = trustregions(problem, [], options); %#ok<ASGLU>
0106
0107     % To extract the eigenvalues, solve the small p-by-p symmetric
0108     % eigenvalue problem.
0109     [Vsol, Dsol] = eig(Xsol'*(A*Xsol));
0110     Ssol = diag(Dsol);
0111
0112     % To extract the eigenvectors, rotate Xsol by the p-by-p orthogonal
0113     % matrix Vsol.
0114     Xsol = Xsol*Vsol;
0115
0116     % This quantity should be small.
0117     % norm(A*Xsol - B*Xsol*diag(Ssol));
0118
0119 end```

Generated on Fri 30-Sep-2022 13:18:25 by m2html © 2005