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.

 See also dominant_invariant_subspace nonlinear_eigenspace

CROSS-REFERENCE INFORMATION ^

This function calls: 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 %
0038 % See also dominant_invariant_subspace nonlinear_eigenspace
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.
0050     
0051     % Generate some random data to test the function
0052     if ~exist('A', 'var') || isempty(A)
0053         n = 128;
0054         A = randn(n);
0055         A = (A+A')/2;
0056     end
0057     if ~exist('B', 'var') || isempty(B)
0058         n = size(A, 1);
0059         e = ones(n, 1);
0060         B = spdiags([-e 2*e -e], -1:1, n, n); % Symmetric positive definite
0061     end
0062     
0063     if ~exist('p', 'var') || isempty(p)
0064         p = 3;
0065     end
0066     
0067     % Make sure the input matrix is square and symmetric
0068     n = size(A, 1);
0069     assert(isreal(A), 'A must be real.')
0070     assert(size(A, 2) == n, 'A must be square.');
0071     assert(norm(A-A', 'fro') < n*eps, 'A must be symmetric.');
0072     assert(p <= n, 'p must be smaller than n.');
0073     
0074     % Define the cost and its derivatives on the generalized
0075     % Grassmann manifold, i.e., the column space of all X such that
0076     % X'*B*X is identity.
0077     gGr = grassmanngeneralizedfactory(n, p, B);
0078     
0079     problem.M = gGr;
0080     problem.cost  = @(X)    -trace(X'*A*X);
0081     problem.egrad = @(X)    -2*(A*X); % Only Euclidean gradient needed.
0082     problem.ehess = @(X, H) -2*(A*H); % Only Euclidean Hessian needed.
0083     
0084     % Execute some checks on the derivatives for early debugging.
0085     % These things can be commented out of course.
0086     % checkgradient(problem);
0087     % pause;
0088     % checkhessian(problem);
0089     % pause;
0090     
0091     % Issue a call to a solver. A random initial guess will be chosen and
0092     % default options are selected except for the ones we specify here.
0093     options.Delta_bar = 8*sqrt(p);
0094     options.tolgradnorm = 1e-7;
0095     options.verbosity = 2; % set to 0 to silence the solver, 2 for normal output.
0096     [Xsol, costXsol, info] = trustregions(problem, [], options); %#ok<ASGLU>
0097     
0098     % To extract the eigenvalues, solve the small p-by-p symmetric
0099     % eigenvalue problem.
0100     [Vsol, Dsol] = eig(Xsol'*(A*Xsol));
0101     Ssol = diag(Dsol);
0102     
0103     % To extract the eigenvectors, rotate Xsol by the p-by-p orthogonal
0104     % matrix Vsol.
0105     Xsol = Xsol*Vsol;
0106     
0107     % This quantity should be small.
0108     % norm(A*Xsol - B*Xsol*diag(Ssol));
0109   
0110 end

Generated on Mon 10-Sep-2018 11:48:06 by m2html © 2005