euclideansubspacefactory.m

function M = euclideansubspacefactory(E, proj, dim)
% Returns a manifold struct to optimize over a subspace of a linear manifold
%
% The factory produces a manifold description M for a linear subspace S of
% a linear space E.
%
% Points and tangent vectors on M are represented in memory in exactly the
% same way as they are for E. Furthermore, E is considered the embedding
% space of S. In particular, the Euclidean gradient of a function is
% understood as the gradient of that function in E, whereas the gradient of
% that function in S is here called the Riemannian gradient. Similar
% considerations hold for the Hessian.
%
% The linear space E is described by a structure obtained from a factory.
% For example, E can be the space of vectors, matrices or nD-arrays of real
% or complex numbers, or a subspace of those. E is equipped with a (real)
% inner product, making it a (real) Euclidean space.
%
% The subspace S is described by an orthogonal projector proj from E to E,
% whose image (i.e., span, range, ...) is S. Orthogonality is judged with
% respect to the inner product on E, and S inherits that inner product,
% making it a Euclidean space itself, and a Riemannian submanifold of E.
%
% The dimension of the subspace S must ideally also be specified.
%
% Inputs:
%   E is the factory for a linear space; for example:
%       E = euclideanfactory(n)
%       E = euclideanfactory(n, m)
%       E = euclideanfactory([n1, n2, n3])
%       E = euclideancomplexfactory(n) % and other dimensions too
%       ...
%
%   proj is a function handle that takes as input a point u of E and
%       returns another point of E: the orthogonal projection of u to S.
%       Orthogonality is understood with respect to the inner product on E,
%       that is, E.inner(x, u, v) = <u, v>. Like any orthogonal projector,
%       it must be linear, and it must satisfy:
%           <u, Pv> = <Pu, v>    and    PPw = Pw
%       for all points u, v, w in E.
%       To check these properties on random vectors, call M.checkproj();
%
%   dim is the dimension of the subspace S (an integer). As always in
%       Manopt, we consider linear spaces over the real numbers, hence,
%       this is the dimension of S as a real linear space. For example, the
%       dimension of R^n is n, and the dimension of C^n is 2n.
%
% Output:
%   M is a factory for optimization over the subspace S.
%
%
% Example 1:
%
% The set of upper-triangular matrices of size n x n can be handled as so:
%
%   n = 5;
%   E = euclideanfactory(n, n);
%   M = euclideansubspacefactory(E, @triu, n*(n+1)/2);
%
% Indeed, Matlab's built-in triu function is the orthogonal projector from
% ths space of n x n matrices to the subspace of upper-triangular matrices
% of that same size, and the dimension of that subspace is n*(n+1)/2.
%
% For strictly upper-triangular matrices, create M as follows instead:
%
%   M = euclideansubspacefactory(E, @(X) triu(X, 1), n*(n-1)/2);
%
%
% Example 2:
%
% Among complex vector of length n, consider the subspace S of such vectors
% that can be the discrete Fourier transform of a real vector of length n
% (that is, the subspace generated by feeding all real vectors of length n
% to fft). The factory M below describes that subspace.
%
%   n = 17;
%   E = euclideancomplexfactory(n);
%   proj = @(u) (u + conj(u([1 ; (n:-1:2)'])))/2;
%   M = euclideansubspacefactory(E, proj, n);
%   M.checkproj(); % for debugging
%
%
% Note 1: this factory is designed to work with linear subspaces only: it
% does not work for affine subspaces. Explicitly: S must contain the origin
% of E. If you need support for affine subspaces, let us know on the Manopt
% forum and we can help (or you can share your improved code :)).
%
% Note 2: the linear space E can itself be a linear subspace of another.
% For example, we can have:
%       E = symmetricfactory(n, k)
%       E = skewfactory(n, k)
%       E = euclideansubspacefactory(...) % recursive nesting
% For these use-cases, bear in mind that the embedding space of M is E,
% hence, the Euclidean gradient and Hessian are expected to be given in E,
% not in the 'bigger' linear space E possibly lives in.
%
% See also: euclideanfactory euclideancomplexfactory symmetricfactory

% This file is part of Manopt: www.manopt.org.
% Original author: Nicolas Boumal, April 25, 2019.
% Contributors: 
% Change log: 

    M = E;
    
    M.name = @() ['Subspace of ', E.name()];
    
    M.proj = @(x, u) proj(u);
    
    M.egrad2rgrad = @(x, eg) proj(eg);
    
    M.ehess2rhess = @(x, eg, eh, d) proj(eh);
    
    M.tangent = @(x, u) proj(u);
    
    M.rand = @() proj(E.rand());
    
    M.randvec = @randvec;
    function v = randvec(x)
        v = proj(E.randvec(x));
        v = v / E.norm(x, v);
    end
    
    if exist('dim', 'var')
        M.dim = @() dim;
        M.typicaldist = @() sqrt(dim);
    else
        M = rmfield(M, 'dim');
        warning('manopt:subspacedim', ...
               ['Since the dimension of the subspace was not specified' ...
                ', M.dim() is not available in the returned factory.']);
    end
    
    M.checkproj = @() checkproj(E, proj);

end

% Tool to check that proj is indeed an orthogonal projector to a subspace
% of E, with respect to the inner product on the linear space E.
function checkproj(E, proj)
    x = E.rand();
    u = E.randvec(x);
    v = E.randvec(x);
    a = randn();
    b = randn();
    % Check that P is linear.
    Paubv = proj(E.lincomb(x, a, u, b, v));
    aPvbPv = E.lincomb(x, a, proj(u), b, proj(v));
    fprintf(['Is proj linear? This number should be zero up to ' ...
             'machine precision:\n\t%.16e\n'], ...
                E.norm(x, E.lincomb(x, 1, Paubv, -1, aPvbPv)));
    % Check that PPw = Pw.
    fprintf(['Is it a projector? This number should be zero up to ' ...
             'machine precision:\n\t%.16e\n'], ...
                E.norm(x, E.lincomb(x, 1, proj(u), -1, proj(proj(u)))));
    % Check that <u, Pv> = <Pu, v>.
    fprintf(['Is it self-adjoint? These two numbers should be equal ' ...
             'up to machine precision:\n\t%.16e\n\t%.16e\n'], ...
                E.inner(x, u, proj(v)), ...
                E.inner(x, proj(u), v));
end