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 improvies distances
0056 %       computation 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     M.hash = @(x) ['z' hashmd5(x(:))];
0165     
0166     M.rand = @() trnsp(random(n, m));
0167     
0168     M.randvec = @(x) trnsp(randomvec(n, m, trnsp(x)));
0169     
0170     M.lincomb = @matrixlincomb;
0171     
0172     M.zerovec = @(x) trnsp(zeros(n, m));
0173     
0174     M.transp = @(x1, x2, d) M.proj(x2, d);
0175     
0176     M.pairmean = @pairmean;
0177     function y = pairmean(x1, x2)
0178         y = trnsp(x1+x2);
0179         y = normalize_columns(y);
0180         y = trnsp(y);
0181     end
0182 
0183     % vec returns a vector representation of an input tangent vector which
0184     % is represented as a matrix. mat returns the original matrix
0185     % representation of the input vector representation of a tangent
0186     % vector. vec and mat are thus inverse of each other. They are
0187     % furthermore isometries between a subspace of R^nm and the tangent
0188     % space at x.
0189     vect = @(X) X(:);
0190     M.vec = @(x, u_mat) vect(trnsp(u_mat));
0191     M.mat = @(x, u_vec) trnsp(reshape(u_vec, [n, m]));
0192     M.vecmatareisometries = @() true;
0193 
0194 end
0195 
0196 % Given a matrix X, returns the same matrix but with each column scaled so
0197 % that they have unit 2-norm.
0198 function X = normalize_columns(X)
0199     % This is faster than norms(X, 2, 1) for small X, and as fast for large X.
0200     nrms = sqrt(sum(X.^2, 1));
0201     X = bsxfun(@times, X, 1./nrms);
0202 end
0203 
0204 % Orthogonal projection of the ambient vector H onto the tangent space at X
0205 function PXH = projection(X, H)
0206 
0207     % Compute the inner product between each vector H(:, i) with its root
0208     % point X(:, i), that is, X(:, i).' * H(:, i). Returns a row vector.
0209     inners = sum(X.*H, 1);
0210     
0211     % Subtract from H the components of the H(:, i)'s that are parallel to
0212     % the root points X(:, i).
0213     PXH = H - bsxfun(@times, X, inners);
0214 
0215     % % Equivalent but slow code:
0216     % m = size(X, 2);
0217     % PXH = zeros(size(H));
0218     % for i = 1 : m
0219     %     PXH(:, i) = H(:, i) - X(:, i) * (X(:, i)'*H(:, i));
0220     % end
0221 
0222 end
0223 
0224 % Uniform random sampling on the sphere.
0225 function x = random(n, m)
0226 
0227     x = normalize_columns(randn(n, m));
0228 
0229 end
0230 
0231 % Random normalized tangent vector at x.
0232 function d = randomvec(n, m, x)
0233 
0234     d = randn(n, m);
0235     d = projection(x, d);
0236     d = d / norm(d(:));
0237 
0238 end

Generated on Fri 08-Sep-2017 12:43:19 by m2html © 2005