Home > manopt > manifolds > oblique > obliquefactory.m

obliquefactory

PURPOSE

Returns a manifold struct to optimize over matrices w/ unit-norm columns.

SYNOPSIS

function M = obliquefactory(n, m, transposed)

DESCRIPTION

``` Returns a manifold struct to optimize over matrices w/ unit-norm columns.

function M = obliquefactory(n, m)
function M = obliquefactory(n, m, transposed)

Oblique manifold: deals with matrices of size n x m such that each column
has unit 2-norm, i.e., is a point on the unit sphere in R^n. The metric
is such that the oblique manifold is a Riemannian submanifold of the
space of nxm matrices with the usual trace inner product, i.e., the usual
metric.

If transposed is set to true (it is false by default), then the matrices
are transposed: a point Y on the manifold is a matrix of size m x n and
each row has unit 2-norm. It is the same geometry, just a different
representation.

CROSS-REFERENCE INFORMATION

This function calls:
• hashmd5 Computes the MD5 hash of input data.
• matrixlincomb Linear combination function for tangent vectors represented as matrices.
This function is called by:
• elliptope_SDP Solver for semidefinite programs (SDP's) with unit diagonal constraints.
• packing_on_the_sphere Return a set of points spread out on the sphere.
• thomson_problem Simple attempt at computing n well distributed points on a sphere in R^d.

SOURCE CODE

```0001 function M = obliquefactory(n, m, transposed)
0002 % Returns a manifold struct to optimize over matrices w/ unit-norm columns.
0003 %
0004 % function M = obliquefactory(n, m)
0005 % function M = obliquefactory(n, m, transposed)
0006 %
0007 % Oblique manifold: deals with matrices of size n x m such that each column
0008 % has unit 2-norm, i.e., is a point on the unit sphere in R^n. The metric
0009 % is such that the oblique manifold is a Riemannian submanifold of the
0010 % space of nxm matrices with the usual trace inner product, i.e., the usual
0011 % metric.
0012 %
0013 % If transposed is set to true (it is false by default), then the matrices
0014 % are transposed: a point Y on the manifold is a matrix of size m x n and
0015 % each row has unit 2-norm. It is the same geometry, just a different
0016 % representation.
0017 %
0019
0020 % This file is part of Manopt: www.manopt.org.
0021 % Original author: Nicolas Boumal, Dec. 30, 2012.
0022 % Contributors:
0023 % Change log:
0024 %
0025 %    July 16, 2013 (NB)
0026 %       Added 'transposed' option, mainly for ease of comparison with the
0027 %       elliptope geometry.
0028 %
0029 %    Nov. 29, 2013 (NB)
0030 %       Added normalize_columns function to make it easier to exploit the
0031 %       bsxfun formulation of column normalization, which avoids using for
0032 %       loops and provides performance gains. The exponential still uses a
0033 %       for loop.
0034 %
0035 %    April 4, 2015 (NB)
0036 %       Log function modified to avoid NaN's appearing for close by points.
0037 %
0038 %    April 13, 2015 (NB)
0039 %       Exponential now without for-loops.
0040 %
0041 %   Oct. 8, 2016 (NB)
0042 %       Code for exponential was simplified to only treat the zero vector
0043 %       as a particular case.
0044 %
0045 %  Oct. 21, 2016 (NB)
0046 %       Bug caught in M.log: the function called v = M.proj(x1, x2 - x1),
0047 %       which internally applies transp to inputs and outputs. But since
0048 %       M.log had already taken care of transposing things, this introduced
0049 %       a bug (which only triggered if using M.log in transposed mode.)
0050 %       The code now calls "v = projection(x1, x2 - x1);" since projection
0051 %       assumes the inputs and outputs do not need to be transposed.
0052 %
0053 %   July 20, 2017 (NB)
0054 %       Distance function is now accurate for close-by points. See notes
0055 %       inside the spherefactory file for details. Also improves distance
0056 %       computations as part of the log function.
0057
0058
0059     if ~exist('transposed', 'var') || isempty(transposed)
0060         transposed = false;
0061     end
0062
0063     if transposed
0064         trnsp = @(X) X';
0065     else
0066         trnsp = @(X) X;
0067     end
0068
0069     M.name = @() sprintf('Oblique manifold OB(%d, %d)', n, m);
0070
0071     M.dim = @() (n-1)*m;
0072
0073     M.inner = @(x, d1, d2) d1(:)'*d2(:);
0074
0075     M.norm = @(x, d) norm(d(:));
0076
0077     M.dist = @(x, y) norm(real(2*asin(.5*sqrt(sum(trnsp(x - y).^2, 1)))));
0078
0079     M.typicaldist = @() pi*sqrt(m);
0080
0081     M.proj = @(X, U) trnsp(projection(trnsp(X), trnsp(U)));
0082
0083     M.tangent = M.proj;
0084
0085     % For Riemannian submanifolds, converting a Euclidean gradient into a
0086     % Riemannian gradient amounts to an orthogonal projection.
0088
0089     M.ehess2rhess = @ehess2rhess;
0090     function rhess = ehess2rhess(X, egrad, ehess, U)
0091         X = trnsp(X);
0093         ehess = trnsp(ehess);
0094         U = trnsp(U);
0095
0096         PXehess = projection(X, ehess);
0098         rhess = PXehess - bsxfun(@times, U, inners);
0099
0100         rhess = trnsp(rhess);
0101     end
0102
0103     M.exp = @exponential;
0104     % Exponential on the oblique manifold
0105     function y = exponential(x, d, t)
0106         x = trnsp(x);
0107         d = trnsp(d);
0108
0109         if nargin < 3
0110             % t = 1;
0111             td = d;
0112         else
0113             td = t*d;
0114         end
0115
0116         nrm_td = sqrt(sum(td.^2, 1));
0117
0118         y = bsxfun(@times, x, cos(nrm_td)) + ...
0119             bsxfun(@times, td, sin(nrm_td) ./ nrm_td);
0120
0121         % For those columns where the step is 0, replace y by x
0122         exclude = (nrm_td == 0);
0123         y(:, exclude) = x(:, exclude);
0124
0125         y = trnsp(y);
0126     end
0127
0128     M.log = @logarithm;
0129     function v = logarithm(x1, x2)
0130         x1 = trnsp(x1);
0131         x2 = trnsp(x2);
0132
0133         v = projection(x1, x2 - x1);
0134         dists = real(2*asin(.5*sqrt(sum((x1 - x2).^2, 1))));
0135         norms = real(sqrt(sum(v.^2, 1)));
0136         factors = dists./norms;
0137         % For very close points, dists is almost equal to norms, but
0138         % because they are both almost zero, the division above can return
0139         % NaN's. To avoid that, we force those ratios to 1.
0140         factors(dists <= 1e-10) = 1;
0141         v = bsxfun(@times, v, factors);
0142
0143         v = trnsp(v);
0144     end
0145
0146     M.retr = @retraction;
0147     % Retraction on the oblique manifold
0148     function y = retraction(x, d, t)
0149         x = trnsp(x);
0150         d = trnsp(d);
0151
0152         if nargin < 3
0153             % t = 1;
0154             td = d;
0155         else
0156             td = t*d;
0157         end
0158
0159         y = normalize_columns(x + td);
0160
0161         y = trnsp(y);
0162     end
0163
0164     % Inverse retraction: see spherefactory.m for background
0165     M.invretr = @inverse_retraction;
0166     function d = inverse_retraction(x, y)
0167
0168         x = trnsp(x);
0169         y = trnsp(y);
0170
0171         d = bsxfun(@times, y, 1./sum(x.*y, 1)) - x;
0172
0173         d = trnsp(d);
0174
0175     end
0176
0177
0178     M.hash = @(x) ['z' hashmd5(x(:))];
0179
0180     M.rand = @() trnsp(random(n, m));
0181
0182     M.randvec = @(x) trnsp(randomvec(n, m, trnsp(x)));
0183
0184     M.lincomb = @matrixlincomb;
0185
0186     M.zerovec = @(x) trnsp(zeros(n, m));
0187
0188     M.transp = @(x1, x2, d) M.proj(x2, d);
0189
0190     M.pairmean = @pairmean;
0191     function y = pairmean(x1, x2)
0192         y = trnsp(x1+x2);
0193         y = normalize_columns(y);
0194         y = trnsp(y);
0195     end
0196
0197     % vec returns a vector representation of an input tangent vector which
0198     % is represented as a matrix. mat returns the original matrix
0199     % representation of the input vector representation of a tangent
0200     % vector. vec and mat are thus inverse of each other. They are
0201     % furthermore isometries between a subspace of R^nm and the tangent
0202     % space at x.
0203     vect = @(X) X(:);
0204     M.vec = @(x, u_mat) vect(trnsp(u_mat));
0205     M.mat = @(x, u_vec) trnsp(reshape(u_vec, [n, m]));
0206     M.vecmatareisometries = @() true;
0207
0208 end
0209
0210 % Given a matrix X, returns the same matrix but with each column scaled so
0211 % that they have unit 2-norm.
0212 function X = normalize_columns(X)
0213     % This is faster than norms(X, 2, 1) for small X, and as fast for large X.
0214     nrms = sqrt(sum(X.^2, 1));
0215     X = bsxfun(@times, X, 1./nrms);
0216 end
0217
0218 % Orthogonal projection of the ambient vector H onto the tangent space at X
0219 function PXH = projection(X, H)
0220
0221     % Compute the inner product between each vector H(:, i) with its root
0222     % point X(:, i), that is, X(:, i)' * H(:, i). Returns a row vector.
0223     inners = sum(X.*H, 1);
0224
0225     % Subtract from H the components of the H(:, i)'s that are parallel to
0226     % the root points X(:, i).
0227     PXH = H - bsxfun(@times, X, inners);
0228
0229     % % Equivalent but slow code:
0230     % m = size(X, 2);
0231     % PXH = zeros(size(H));
0232     % for i = 1 : m
0233     %     PXH(:, i) = H(:, i) - X(:, i) * (X(:, i)'*H(:, i));
0234     % end
0235
0236 end
0237
0238 % Uniform random sampling on the sphere.
0239 function x = random(n, m)
0240
0241     x = normalize_columns(randn(n, m));
0242
0243 end
0244
0245 % Random normalized tangent vector at x.
0246 function d = randomvec(n, m, x)
0247
0248     d = randn(n, m);
0249     d = projection(x, d);
0250     d = d / norm(d(:));
0251
0252 end```

Generated on Mon 10-Sep-2018 11:48:06 by m2html © 2005