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_eigenspace0001 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 % Aug. 20, 2021 (XJ): Added AD to compute the egrad and the ehess 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); 0082 problem.egrad = @(X) -2*(A*X); % Only Euclidean gradient needed. 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 0091 % problem = manoptAD(problem); 0092 0093 % Execute some checks on the derivatives for early debugging. 0094 % These things can be commented out of course. 0095 % checkgradient(problem); 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); 0103 options.tolgradnorm = 1e-7; 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