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.

 See also: spherefactory obliquecomplexfactory

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

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 %
0018 % See also: spherefactory obliquecomplexfactory
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.
0087     M.egrad2rgrad = M.proj;
0088     
0089     M.ehess2rhess = @ehess2rhess;
0090     function rhess = ehess2rhess(X, egrad, ehess, U)
0091         X = trnsp(X);
0092         egrad = trnsp(egrad);
0093         ehess = trnsp(ehess);
0094         U = trnsp(U);
0095         
0096         PXehess = projection(X, ehess);
0097         inners = sum(X.*egrad, 1);
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