Home > examples > dominant_invariant_subspace_complex.m

dominant_invariant_subspace_complex

PURPOSE ^

Returns a unitary basis of the dominant invariant p-subspace of A.

SYNOPSIS ^

function [X, info] = dominant_invariant_subspace_complex(A, p)

DESCRIPTION ^

 Returns a unitary basis of the dominant invariant p-subspace of A.

 function X = dominant_invariant_subspace(A, p)

 Input: A complex, Hermitian matrix A of size nxn and an integer p < n.
 Output: A complex, unitary matrix X of size nxp such that trace(X'*A*X)
         is maximized. That is, the columns of X form a unitary basis
         of a dominant subspace of dimension p of A.

 The optimization is performed on the complex Grassmann manifold, since
 only the space spanned by the columns of X matters.

 See dominant_invariant_subspace for more details in the real case.

 See also: dominant_invariant_subspace grassmanncomplexfactory

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [X, info] = dominant_invariant_subspace_complex(A, p)
0002 % Returns a unitary basis of the dominant invariant p-subspace of A.
0003 %
0004 % function X = dominant_invariant_subspace(A, p)
0005 %
0006 % Input: A complex, Hermitian matrix A of size nxn and an integer p < n.
0007 % Output: A complex, unitary matrix X of size nxp such that trace(X'*A*X)
0008 %         is maximized. That is, the columns of X form a unitary basis
0009 %         of a dominant subspace of dimension p of A.
0010 %
0011 % The optimization is performed on the complex Grassmann manifold, since
0012 % only the space spanned by the columns of X matters.
0013 %
0014 % See dominant_invariant_subspace for more details in the real case.
0015 %
0016 % See also: dominant_invariant_subspace grassmanncomplexfactory
0017 
0018 % This file is part of Manopt and is copyrighted. See the license file.
0019 %
0020 % Main author: Nicolas Boumal, June 30, 2015
0021 % Contributors:
0022 %
0023 % Change log:
0024     
0025     % Generate some random data to test the function
0026     if ~exist('A', 'var') || isempty(A)
0027         A = randn(128) + 1i*randn(128);
0028         A = (A+A')/2;
0029     end
0030     if ~exist('p', 'var') || isempty(p)
0031         p = 3;
0032     end
0033     
0034     % Make sure the input matrix is Hermitian
0035     n = size(A, 1);
0036     assert(size(A, 2) == n, 'A must be square.');
0037     assert(norm(A-A', 'fro') < n*eps, 'A must be Hermitian.');
0038     assert(p<=n, 'p must be smaller than n.');
0039     
0040     % Define the cost and its derivatives on the complex Grassmann manifold
0041     Gr = grassmanncomplexfactory(n, p);
0042     problem.M = Gr;
0043     problem.cost  = @(X)    -real(trace(X'*A*X));
0044     problem.egrad = @(X)    -2*A*X;
0045     problem.ehess = @(X, H) -2*A*H;
0046     
0047     % Execute some checks on the derivatives for early debugging.
0048     % These can be commented out.
0049     % checkgradient(problem);
0050     % pause;
0051     % checkhessian(problem);
0052     % pause;
0053     
0054     % Issue a call to a solver. A random initial guess will be chosen and
0055     % default options are selected except for the ones we specify here.
0056     options.Delta_bar = 8*sqrt(p);
0057     [X, costX, info, options] = trustregions(problem, [], options); %#ok<ASGLU>
0058     
0059     fprintf('Options used:\n');
0060     disp(options);
0061     
0062     % For our information, Manopt can also compute the spectrum of the
0063     % Riemannian Hessian on the tangent space at (any) X. Computing the
0064     % spectrum at the solution gives us some idea of the conditioning of
0065     % the problem. If we were to implement a preconditioner for the
0066     % Hessian, this would also inform us on its performance.
0067     %
0068     % Notice that (typically) all eigenvalues of the Hessian at the
0069     % solution are positive, i.e., we find an isolated minimizer. If we
0070     % replace the Grassmann manifold by the Stiefel manifold, hence still
0071     % optimizing over orthonormal matrices but ignoring the invariance
0072     % cost(XQ) = cost(X) for all Q orthogonal, then we see
0073     % dim O(p) = p(p-1)/2 zero eigenvalues in the Hessian spectrum, making
0074     % the optimizer not isolated anymore.
0075     if Gr.dim() < 512
0076         evs = hessianspectrum(problem, X);
0077         stairs(sort(evs));
0078         title(['Eigenvalues of the Hessian of the cost function ' ...
0079                'at the solution']);
0080         xlabel('Eigenvalue number (sorted)');
0081         ylabel('Value of the eigenvalue');
0082     end
0083 
0084 end

Generated on Sat 12-Nov-2016 14:11:22 by m2html © 2005