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