Home > manopt > manifolds > sphere > spherefactory.m

# spherefactory

## PURPOSE

Returns a manifold struct to optimize over unit-norm vectors or matrices.

## SYNOPSIS

function M = spherefactory(n, m, gpuflag)

## DESCRIPTION

``` Returns a manifold struct to optimize over unit-norm vectors or matrices.

function M = spherefactory(n)
function M = spherefactory(n, m)
function M = spherefactory(n, m, gpuflag)

Manifold of n-by-m real matrices of unit Frobenius norm.
By default, m = 1, which corresponds to the unit sphere in R^n. The
metric is such that the sphere is a Riemannian submanifold of the space
of nxm matrices with the usual trace inner product, i.e., the usual
metric.

Set gpuflag = true to have points, tangent vectors and ambient vectors
stored on the GPU. If so, computations can be done on the GPU directly.

## CROSS-REFERENCE INFORMATION

This function calls:
• norm NORM Norm of a TT/MPS tensor.
• norm NORM Norm of a TT/MPS block-mu tensor.
• factorygpuhelper Returns a manifold struct to optimize over unit-norm vectors or matrices.
• hashmd5 Computes the MD5 hash of input data.
• matrixlincomb Linear combination function for tangent vectors represented as matrices.
• sinxoverx Computes sin(x) ./ x entrywise, such that sin(0) / 0 yields 1.
This function is called by:
• using_counters Manopt example on how to use counters during optimization. Typical uses,
• basic_example_AD A basic example that shows how to apply automatic differentiation to
• realtest_AD1 Test AD for a real optimization problem on a product manifold (struct)
• realtest_AD2 Test AD for a real optimization problem on a power manifold (cell)
• realtest_AD3 Test AD for a real optimization problem on a manifold which is stored in

## SOURCE CODE

```0001 function M = spherefactory(n, m, gpuflag)
0002 % Returns a manifold struct to optimize over unit-norm vectors or matrices.
0003 %
0004 % function M = spherefactory(n)
0005 % function M = spherefactory(n, m)
0006 % function M = spherefactory(n, m, gpuflag)
0007 %
0008 % Manifold of n-by-m real matrices of unit Frobenius norm.
0009 % By default, m = 1, which corresponds to the unit sphere in R^n. The
0010 % metric is such that the sphere is a Riemannian submanifold of the space
0011 % of nxm matrices with the usual trace inner product, i.e., the usual
0012 % metric.
0013 %
0014 % Set gpuflag = true to have points, tangent vectors and ambient vectors
0015 % stored on the GPU. If so, computations can be done on the GPU directly.
0016 %
0018
0019 % This file is part of Manopt: www.manopt.org.
0020 % Original author: Nicolas Boumal, Dec. 30, 2012.
0021 % Contributors:
0022 % Change log:
0023 %
0024 %   Oct. 8, 2016 (NB)
0025 %       Code for exponential was simplified to only treat the zero vector
0026 %       as a particular case.
0027 %
0028 %   Oct. 22, 2016 (NB)
0029 %       Distance function dist now significantly more accurate for points
0030 %       within 1e-7 and less from each other.
0031 %
0032 %   July 20, 2017 (NB)
0033 %       Following conversations with Bruno Iannazzo and P.-A. Absil,
0034 %       the distance function is now even more accurate.
0035 %
0036 %   Sep. 7, 2017 (NB)
0037 %       New isometric vector transport available in M.isotransp,
0038 %       contributed by Changshuo Liu.
0039 %
0040 %   April 17, 2018 (NB)
0041 %       ehess2rhess: Used to compute projection of ehess, then subtract a
0042 %       multiple of u (which is assumed tangent.) Now, similarly to what
0043 %       happens in stiefelfactory, we first subtract the multiple of u from
0044 %       ehess, then we project. Mathematically, these operations are the
0045 %       same. Numerically, the former version used to be better because tCG
0046 %       in trustregions had some drift near fine convergence. Now that the
0047 %       drift in tCG has been fixed, it is reasonable to apply the
0048 %       projection last, to ensure best tangency of the output.
0049 %
0050 %   July 18, 2018 (NB)
0051 %       Added the inverse retraction (M.invretr) for the sphere.
0052 %
0053 %   Aug. 3, 2018 (NB)
0054 %       Added GPU support: just set gpuflag = true.
0055 %
0056 %   Jan. 8, 2021 (NB)
0058
0059
0060     if ~exist('m', 'var') || isempty(m)
0061         m = 1;
0062     end
0063     if ~exist('gpuflag', 'var') || isempty(gpuflag)
0064         gpuflag = false;
0065     end
0066
0067     % If gpuflag is active, new arrays (e.g., via rand, randn, zeros, ones)
0068     % are created directly on the GPU; otherwise, they are created in the
0069     % usual way (in double precision).
0070     if gpuflag
0071         array_type = 'gpuArray';
0072     else
0073         array_type = 'double';
0074     end
0075
0076
0077     if m == 1
0078         M.name = @() sprintf('Sphere S^%d', n-1);
0079     else
0080         M.name = @() sprintf('Unit F-norm %dx%d matrices', n, m);
0081     end
0082
0083     M.dim = @() n*m-1;
0084
0085     M.inner = @(x, d1, d2) d1(:)'*d2(:);
0086
0087     M.norm = @(x, d) norm(d, 'fro');
0088
0089     M.dist = @dist;
0090     function d = dist(x, y)
0091
0092         % The following code is mathematically equivalent to the
0093         % computation d = acos(x(:)'*y(:)) but is much more accurate when
0094         % x and y are close.
0095
0096         chordal_distance = norm(x - y, 'fro');
0097         d = real(2*asin(.5*chordal_distance));
0098
0099         % Note: for x and y almost antipodal, the accuracy is good but not
0100         % as good as possible. One way to improve it is by using the
0101         % following branching:
0102         % % if chordal_distance > 1.9
0103         % %     d = pi - dist(x, -y);
0104         % % end
0105         % It is rarely necessary to compute the distance between
0106         % almost-antipodal points with full accuracy in Manopt, hence we
0107         % favor a simpler code.
0108
0109     end
0110
0111     M.typicaldist = @() pi;
0112
0113     M.proj = @(x, d) d - x*(x(:)'*d(:));
0114
0115     M.tangent = M.proj;
0116
0117     M.tangent2ambient_is_identity = true;
0118     M.tangent2ambient = @(X, U) U;
0119
0120     % For Riemannian submanifolds, converting a Euclidean gradient into a
0121     % Riemannian gradient amounts to an orthogonal projection.
0123
0124     M.ehess2rhess = @ehess2rhess;
0125     function rhess = ehess2rhess(x, egrad, ehess, u)
0126         rhess = M.proj(x, ehess - (x(:)'*egrad(:))*u);
0127     end
0128
0129     M.exp = @exponential;
0130
0131     M.retr = @retraction;
0132     M.invretr = @inverse_retraction;
0133
0134     M.log = @logarithm;
0135     function v = logarithm(x1, x2)
0136         v = M.proj(x1, x2 - x1);
0137         di = M.dist(x1, x2);
0138         % If the two points are "far apart", correct the norm.
0139         if di > 1e-6
0140             nv = norm(v, 'fro');
0141             v = v * (di / nv);
0142         end
0143     end
0144
0145     M.hash = @(x) ['z' hashmd5(x(:))];
0146
0147     M.rand = @() random(n, m, array_type);
0148
0149     M.randvec = @(x) randomvec(n, m, x, array_type);
0150
0151     M.zerovec = @(x) zeros(n, m, array_type);
0152
0153     M.lincomb = @matrixlincomb;
0154
0155     M.transp = @(x1, x2, d) M.proj(x2, d);
0156
0157     % Isometric vector transport of d from the tangent space at x1 to x2.
0158     % This is actually a parallel vector transport, see Ch. 5 in
0159     % http://epubs.siam.org/doi/pdf/10.1137/16M1069298
0160     % "A Riemannian Gradient Sampling Algorithm for Nonsmooth Optimization
0161     %  on Manifolds", by Hosseini and Uschmajew, SIOPT 2017
0162     M.isotransp = @(x1, x2, d) isometricTransp(x1, x2, d);
0163     function Td = isometricTransp(x1, x2, d)
0164         v = logarithm(x1, x2);
0165         dist_x1x2 = norm(v, 'fro');
0166         if dist_x1x2 > 0
0167             u = v / dist_x1x2;
0168             utd = u(:)'*d(:);
0169             Td = d + (cos(dist_x1x2)-1)*utd*u ...
0170                     -  sin(dist_x1x2)  *utd*x1;
0171         else
0172             % x1 == x2, so the transport is identity
0173             Td = d;
0174         end
0175     end
0176
0177     M.pairmean = @pairmean;
0178     function y = pairmean(x1, x2)
0179         y = x1+x2;
0180         y = y / norm(y, 'fro');
0181     end
0182
0183     M.vec = @(x, u_mat) u_mat(:);
0184     M.mat = @(x, u_vec) reshape(u_vec, [n, m]);
0185     M.vecmatareisometries = @() true;
0186
0187
0188     % Automatically convert a number of tools to support GPU.
0189     if gpuflag
0190         M = factorygpuhelper(M);
0191     end
0192
0193
0194 end
0195
0196 % Exponential on the sphere
0197 function y = exponential(x, d, t)
0198
0199     if nargin == 2
0200         % t = 1
0201         td = d;
0202     else
0203         td = t*d;
0204     end
0205
0206     nrm_td = norm(td, 'fro');
0207     y = x*cos(nrm_td) + td*sinxoverx(nrm_td);
0208
0209 end
0210
0211 % Retraction on the sphere
0212 function y = retraction(x, d, t)
0213
0214     if nargin == 2
0215         % t = 1;
0216         td = d;
0217     else
0218         td = t*d;
0219     end
0220
0221     y = x + td;
0222     y = y / norm(y, 'fro');
0223
0224 end
0225
0226 % Given x and y two points on the manifold, if there exists a tangent
0227 % vector d at x such that Retr_x(d) = y, this function returns d.
0228 function d = inverse_retraction(x, y)
0229
0230     % Since
0231     %   x + d = y*||x + d||
0232     % and x'd = 0, multiply the above by x' on the left:
0233     %   1 + 0 = x'y * ||x + d||
0234     % Then solve for d:
0235
0236     d = y/(x(:)'*y(:)) - x;
0237
0238 end
0239
0240 % Uniform random sampling on the sphere.
0241 function x = random(n, m, array_type)
0242
0243     x = randn(n, m, array_type);
0244     x = x / norm(x, 'fro');
0245
0246 end
0247
0248 % Random normalized tangent vector at x.
0249 function d = randomvec(n, m, x, array_type)
0250
0251     d = randn(n, m, array_type);
0252     d = d - x*(x(:)'*d(:));
0253     d = d / norm(d, 'fro');
0254
0255 end```

Generated on Fri 30-Sep-2022 13:18:25 by m2html © 2005