Home > manopt > manifolds > euclidean > euclideansubspacefactory.m

euclideansubspacefactory

PURPOSE ^

Returns a manifold struct to optimize over a subspace of a linear manifold

SYNOPSIS ^

function M = euclideansubspacefactory(E, proj, dim)

DESCRIPTION ^

 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:

 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 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function M = euclideansubspacefactory(E, proj, dim)
0002 % Returns a manifold struct to optimize over a subspace of a linear manifold
0003 %
0004 % The factory produces a manifold description M for a linear subspace S of
0005 % a linear space E.
0006 %
0007 % Points and tangent vectors on M are represented in memory in exactly the
0008 % same way as they are for E. Furthermore, E is considered the embedding
0009 % space of S. In particular, the Euclidean gradient of a function is
0010 % understood as the gradient of that function in E, whereas the gradient of
0011 % that function in S is here called the Riemannian gradient. Similar
0012 % considerations hold for the Hessian.
0013 %
0014 % The linear space E is described by a structure obtained from a factory.
0015 % For example, E can be the space of vectors, matrices or nD-arrays of real
0016 % or complex numbers, or a subspace of those. E is equipped with a (real)
0017 % inner product, making it a (real) Euclidean space.
0018 %
0019 % The subspace S is described by an orthogonal projector proj from E to E,
0020 % whose image (i.e., span, range, ...) is S. Orthogonality is judged with
0021 % respect to the inner product on E, and S inherits that inner product,
0022 % making it a Euclidean space itself, and a Riemannian submanifold of E.
0023 %
0024 % The dimension of the subspace S must ideally also be specified.
0025 %
0026 % Inputs:
0027 %   E is the factory for a linear space; for example:
0028 %       E = euclideanfactory(n)
0029 %       E = euclideanfactory(n, m)
0030 %       E = euclideanfactory([n1, n2, n3])
0031 %       E = euclideancomplexfactory(n) % and other dimensions too
0032 %       ...
0033 %
0034 %   proj is a function handle that takes as input a point u of E and
0035 %       returns another point of E: the orthogonal projection of u to S.
0036 %       Orthogonality is understood with respect to the inner product on E,
0037 %       that is, E.inner(x, u, v) = <u, v>. Like any orthogonal projector,
0038 %       it must be linear, and it must satisfy:
0039 %           <u, Pv> = <Pu, v>    and    PPw = Pw
0040 %       for all points u, v, w in E.
0041 %       To check these properties on random vectors, call M.checkproj();
0042 %
0043 %   dim is the dimension of the subspace S (an integer). As always in
0044 %       Manopt, we consider linear spaces over the real numbers, hence,
0045 %       this is the dimension of S as a real linear space. For example, the
0046 %       dimension of R^n is n, and the dimension of C^n is 2n.
0047 %
0048 % Output:
0049 %   M is a factory for optimization over the subspace S.
0050 %
0051 %
0052 % Example:
0053 %
0054 % Among complex vector of length n, consider the subspace S of such vectors
0055 % that can be the discrete Fourier transform of a real vector of length n
0056 % (that is, the subspace generated by feeding all real vectors of length n
0057 % to fft). The factory M below describes that subspace.
0058 %
0059 %   n = 17;
0060 %   E = euclideancomplexfactory(n);
0061 %   proj = @(u) (u + conj(u([1 ; (n:-1:2)'])))/2;
0062 %   M = euclideansubspacefactory(E, proj, n);
0063 %   M.checkproj(); % for debugging
0064 %
0065 %
0066 % Note 1: this factory is designed to work with linear subspaces only: it
0067 % does not work for affine subspaces. Explicitly: S must contain the origin
0068 % of E. If you need support for affine subspaces, let us know on the Manopt
0069 % forum and we can help (or share your improved code :)).
0070 %
0071 % Note 2: the linear space E can itself be a linear subspace of another.
0072 % For example, we can have:
0073 %       E = symmetricfactory(n, k)
0074 %       E = skewfactory(n, k)
0075 %       E = euclideansubspacefactory(...) % recursive nesting
0076 % For these use-cases, bear in mind that the embedding space of M is E,
0077 % hence, the Euclidean gradient and Hessian are expected to be given in E,
0078 % not in the 'bigger' linear space E possibly lives in.
0079 %
0080 % See also: euclideanfactory euclideancomplexfactory
0081 
0082 % This file is part of Manopt: www.manopt.org.
0083 % Original author: Nicolas Boumal, April 25, 2019.
0084 % Contributors:
0085 % Change log:
0086 
0087     M = E;
0088     
0089     M.name = @() ['Subspace of ', E.name()];
0090     
0091     M.proj = @(x, u) proj(u);
0092     
0093     M.egrad2rgrad = @(x, eg) proj(eg);
0094     
0095     M.ehess2rhess = @(x, eg, eh, d) proj(eh);
0096     
0097     M.tangent = @(x, u) proj(u);
0098     
0099     M.rand = @() proj(E.rand());
0100     
0101     M.randvec = @randvec;
0102     function v = randvec(x)
0103         v = proj(E.randvec(x));
0104         v = v / E.norm(x, v);
0105     end
0106     
0107     if exist('dim', 'var')
0108         M.dim = @() dim;
0109         M.typicaldist = @() sqrt(dim);
0110     else
0111         M = rmfield(M, 'dim');
0112         warning('manopt:subspacedim', ...
0113                ['Since the dimension of the subspace was not specified' ...
0114                 ', M.dim() is not available in the returned factory.']);
0115     end
0116     
0117     M.checkproj = @() checkproj(E, proj);
0118 
0119 end
0120 
0121 % Tool to check that proj is indeed an orthogonal projector to a subspace
0122 % of E, with respect to the inner product on the linear space E.
0123 function checkproj(E, proj)
0124     x = E.rand();
0125     u = E.randvec(x);
0126     v = E.randvec(x);
0127     a = randn();
0128     b = randn();
0129     % Check that P is linear.
0130     Paubv = proj(E.lincomb(x, a, u, b, v));
0131     aPvbPv = E.lincomb(x, a, proj(u), b, proj(v));
0132     fprintf(['Is proj linear? This number should be zero up to ' ...
0133              'machine precision:\n\t%.16e\n'], ...
0134                 E.norm(x, E.lincomb(x, 1, Paubv, -1, aPvbPv)));
0135     % Check that PPw = Pw.
0136     fprintf(['Is it a projector? This number should be zero up to ' ...
0137              'machine precision:\n\t%.16e\n'], ...
0138                 E.norm(x, E.lincomb(x, 1, proj(u), -1, proj(proj(u)))));
0139     % Check that <u, Pv> = <Pu, v>.
0140     fprintf(['Is it self-adjoint? These two numbers should be equal ' ...
0141              'up to machine precision:\n\t%.16e\n\t%.16e\n'], ...
0142                 E.inner(x, u, proj(v)), ...
0143                 E.inner(x, proj(u), v));
0144 end

Generated on Sun 05-Sep-2021 17:57:00 by m2html © 2005