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

- grassmanngeneralizedfactory Returns a manifold struct of "scaled" vector subspaces.
- trustregions Riemannian trust-regions solver for optimization on manifolds.

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