Home > examples > dominant_invariant_subspace.m

# dominant_invariant_subspace

## PURPOSE Returns an orthonormal basis of the dominant invariant p-subspace of A.

## SYNOPSIS function [X, info] = dominant_invariant_subspace(A, p)

## DESCRIPTION ``` 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.

## CROSS-REFERENCE INFORMATION This function calls:
• grassmannfactory Returns a manifold struct to optimize over the space of vector subspaces.
• trustregions Riemannian trust-regions solver for optimization on manifolds.
• hessianspectrum Returns the eigenvalues of the (preconditioned) Hessian at x.
This function is called by:

## SOURCE CODE ```0001 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 %
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);
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.
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```

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