Returns an orthonormal basis of the dominant invariant p-subspace of A.
function X = dominant_invariant_subspace(A, p)
Input: A real, symmetric matrix A of size nxn and an integer p < n.
Output: A real, orthonormal matrix X of size nxp such that trace(X'*A*X)
is maximized. That is, the columns of X form an orthonormal basis
of a dominant subspace of dimension p of A. These are thus
eigenvectors associated with the largest eigenvalues of A (in no
particular order). Sign is important: 2 is deemed a larger
eigenvalue than -5.
The optimization is performed on the Grassmann manifold, since only the
space spanned by the columns of X matters. The implementation is short to
show how Manopt can be used to quickly obtain a prototype. To make the
implementation more efficient, one might first try to use the caching
system, that is, use the optional 'store' arguments in the cost, grad and
hess functions. Furthermore, using egrad2rgrad and ehess2rhess is quick
and easy, but not always efficient. Having a look at the formulas
implemented in these functions can help rewrite the code without them,
possibly more efficiently.
See also: dominant_invariant_subspace_complex0001 function [X, info] = dominant_invariant_subspace(A, p) 0002 % Returns an orthonormal basis of the dominant invariant p-subspace of A. 0003 % 0004 % function X = dominant_invariant_subspace(A, p) 0005 % 0006 % Input: A real, symmetric matrix A of size nxn and an integer p < n. 0007 % Output: A real, orthonormal matrix X of size nxp such that trace(X'*A*X) 0008 % is maximized. That is, the columns of X form an orthonormal basis 0009 % of a dominant subspace of dimension p of A. These are thus 0010 % eigenvectors associated with the largest eigenvalues of A (in no 0011 % particular order). Sign is important: 2 is deemed a larger 0012 % eigenvalue than -5. 0013 % 0014 % The optimization is performed on the Grassmann manifold, since only the 0015 % space spanned by the columns of X matters. The implementation is short to 0016 % show how Manopt can be used to quickly obtain a prototype. To make the 0017 % implementation more efficient, one might first try to use the caching 0018 % system, that is, use the optional 'store' arguments in the cost, grad and 0019 % hess functions. Furthermore, using egrad2rgrad and ehess2rhess is quick 0020 % and easy, but not always efficient. Having a look at the formulas 0021 % implemented in these functions can help rewrite the code without them, 0022 % possibly more efficiently. 0023 % 0024 % See also: dominant_invariant_subspace_complex 0025 0026 % This file is part of Manopt and is copyrighted. See the license file. 0027 % 0028 % Main author: Nicolas Boumal, July 5, 2013 0029 % Contributors: 0030 % 0031 % Change log: 0032 % 0033 % NB Dec. 6, 2013: 0034 % We specify a max and initial trust region radius in the options. 0035 % NB Jan. 20, 2018: 0036 % Added a few comments regarding implementation of the cost. 0037 0038 % Generate some random data to test the function 0039 if ~exist('A', 'var') || isempty(A) 0040 A = randn(128); 0041 A = (A+A')/2; 0042 end 0043 if ~exist('p', 'var') || isempty(p) 0044 p = 3; 0045 end 0046 0047 % Make sure the input matrix is square and symmetric 0048 n = size(A, 1); 0049 assert(isreal(A), 'A must be real.') 0050 assert(size(A, 2) == n, 'A must be square.'); 0051 assert(norm(A-A', 'fro') < n*eps, 'A must be symmetric.'); 0052 assert(p<=n, 'p must be smaller than n.'); 0053 0054 % Define the cost and its derivatives on the Grassmann manifold 0055 Gr = grassmannfactory(n, p); 0056 problem.M = Gr; 0057 problem.cost = @(X) -.5*trace(X'*A*X); 0058 problem.grad = @(X) -Gr.egrad2rgrad(X, A*X); 0059 problem.hess = @(X, H) -Gr.ehess2rhess(X, A*X, A*H, H); 0060 0061 % Notice that it would be more efficient to compute trace(X'*A*X) via 0062 % the formula sum(sum(X .* (A*X))) -- the code above is written so as 0063 % to be as close as possible to the familiar mathematical formulas, for 0064 % ease of interpretation. Also, the product A*X is needed for both the 0065 % cost and the gradient, as well as for the Hessian: one can use the 0066 % caching capabilities of Manopt (the store structures) to save on 0067 % redundant computations. 0068 0069 % Execute some checks on the derivatives for early debugging. 0070 % These can be commented out. 0071 % checkgradient(problem); 0072 % pause; 0073 % checkhessian(problem); 0074 % pause; 0075 0076 % Issue a call to a solver. A random initial guess will be chosen and 0077 % default options are selected except for the ones we specify here. 0078 options.Delta_bar = 8*sqrt(p); 0079 [X, costX, info, options] = trustregions(problem, [], options); %#ok<ASGLU> 0080 0081 fprintf('Options used:\n'); 0082 disp(options); 0083 0084 % For our information, Manopt can also compute the spectrum of the 0085 % Riemannian Hessian on the tangent space at (any) X. Computing the 0086 % spectrum at the solution gives us some idea of the conditioning of 0087 % the problem. If we were to implement a preconditioner for the 0088 % Hessian, this would also inform us on its performance. 0089 % 0090 % Notice that (typically) all eigenvalues of the Hessian at the 0091 % solution are positive, i.e., we find an isolated minimizer. If we 0092 % replace the Grassmann manifold by the Stiefel manifold, hence still 0093 % optimizing over orthonormal matrices but ignoring the invariance 0094 % cost(XQ) = cost(X) for all Q orthogonal, then we see 0095 % dim O(p) = p(p-1)/2 zero eigenvalues in the Hessian spectrum, making 0096 % the optimizer not isolated anymore. 0097 if Gr.dim() < 512 0098 evs = hessianspectrum(problem, X); 0099 stairs(sort(evs)); 0100 title(['Eigenvalues of the Hessian of the cost function ' ... 0101 'at the solution']); 0102 xlabel('Eigenvalue number (sorted)'); 0103 ylabel('Value of the eigenvalue'); 0104 end 0105 0106 end