function [X, info] = dominant_invariant_subspace_complex(A, p)
% Returns a unitary basis of the dominant invariant p-subspace of A.
%
% function X = dominant_invariant_subspace_complex(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

% This file is part of Manopt and is copyrighted. See the license file.
%
% Main author: Nicolas Boumal, June 30, 2015
% Contributors:
%
% Change log:
%    
%   Xiaowen Jiang Aug. 31, 2021
%       Added AD to compute the egrad and the ehess  


    % Generate some random data to test the function
    if ~exist('A', 'var') || isempty(A)
        A = randn(128) + 1i*randn(128);
        A = (A+A')/2;
    end
    if ~exist('p', 'var') || isempty(p)
        p = 3;
    end
    
    % Make sure the input matrix is Hermitian
    n = size(A, 1);
    assert(size(A, 2) == n, 'A must be square.');
    assert(norm(A-A', 'fro') < n*eps, 'A must be Hermitian.');
	assert(p<=n, 'p must be smaller than n.');
    
    % Define the cost and its derivatives on the complex Grassmann manifold
    Gr = grassmanncomplexfactory(n, p);
    problem.M = Gr;
    problem.cost  = @(X)    -real(trace(X'*A*X));
    problem.egrad = @(X)    -2*A*X;
    problem.ehess = @(X, H) -2*A*H;
    
    % An alternative way to compute the egrad and the ehess is to use 
    % automatic differentiation provided in the deep learning toolbox
    % (slower). AD does not support complex numbers if the Matlab version
    % is R2021a or earlier. The cost function should be defined differently
    % In this case. See complex_example_AD.m and manoptADhelp.m for more
    % information.
    % problem.cost = @cost_complex;
    %    function f = cost_complex(X)
    %        AX = cprod(A,X);
    %        Xtransp = ctransp(X);
    %        product = cprod(Xtransp,AX);
    %        f = -creal(ctrace(product));
    %    end
    % call manoptAD to automatically obtain the egrad and the ehess
    % problem = manoptAD(problem);
    
    % If the version of Matlab installed is R2021b or later, specify the 
    % cost function in the normal way and call manoptAD. Notice that
    % the function trace is not supported for AD so far. Replace it with 
    % ctrace described in the file manoptADhelp.m
    % problem.cost  = @(X)    -real(ctrace(X'*A*X));
    % problem = manoptAD(problem);

    % Execute some checks on the derivatives for early debugging.
    % These can be commented out.
    % checkgradient(problem);
    % pause;
    % checkhessian(problem);
    % pause;
    
    % Issue a call to a solver. A random initial guess will be chosen and
    % default options are selected except for the ones we specify here.
    options.Delta_bar = 8*sqrt(p);
    [X, costX, info, options] = trustregions(problem, [], options); %#ok<ASGLU>
    
    fprintf('Options used:\n');
    disp(options);
    
    % For our information, Manopt can also compute the spectrum of the
    % Riemannian Hessian on the tangent space at (any) X. Computing the
    % spectrum at the solution gives us some idea of the conditioning of
    % the problem. If we were to implement a preconditioner for the
    % Hessian, this would also inform us on its performance.
    %
    % Notice that (typically) all eigenvalues of the Hessian at the
    % solution are positive, i.e., we find an isolated minimizer. If we
    % replace the Grassmann manifold by the Stiefel manifold, hence still
    % optimizing over orthonormal matrices but ignoring the invariance
    % cost(XQ) = cost(X) for all Q orthogonal, then we see
    % dim O(p) = p(p-1)/2 zero eigenvalues in the Hessian spectrum, making
    % the optimizer not isolated anymore.
    if Gr.dim() < 512
        evs = hessianspectrum(problem, X);
        stairs(sort(evs));
        title(['Eigenvalues of the Hessian of the cost function ' ...
               'at the solution']);
        xlabel('Eigenvalue number (sorted)');
        ylabel('Value of the eigenvalue');
    end

end
