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
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, sinxoverx(nrm_td)); 0120 0121 y = trnsp(y); 0122 end 0123 0124 M.log = @logarithm; 0125 function v = logarithm(x1, x2) 0126 x1 = trnsp(x1); 0127 x2 = trnsp(x2); 0128 0129 v = projection(x1, x2 - x1); 0130 dists = real(2*asin(.5*sqrt(sum((x1 - x2).^2, 1)))); 0131 norms = real(sqrt(sum(v.^2, 1))); 0132 factors = dists./norms; 0133 % For very close points, dists is almost equal to norms, but 0134 % because they are both almost zero, the division above can return 0135 % NaN's. To avoid that, we force those ratios to 1. 0136 factors(dists <= 1e-10) = 1; 0137 v = bsxfun(@times, v, factors); 0138 0139 v = trnsp(v); 0140 end 0141 0142 M.retr = @retraction; 0143 % Retraction on the oblique manifold 0144 function y = retraction(x, d, t) 0145 x = trnsp(x); 0146 d = trnsp(d); 0147 0148 if nargin < 3 0149 % t = 1; 0150 td = d; 0151 else 0152 td = t*d; 0153 end 0154 0155 y = normalize_columns(x + td); 0156 0157 y = trnsp(y); 0158 end 0159 0160 % Inverse retraction: see spherefactory.m for background 0161 M.invretr = @inverse_retraction; 0162 function d = inverse_retraction(x, y) 0163 0164 x = trnsp(x); 0165 y = trnsp(y); 0166 0167 d = bsxfun(@times, y, 1./sum(x.*y, 1)) - x; 0168 0169 d = trnsp(d); 0170 0171 end 0172 0173 0174 M.hash = @(x) ['z' hashmd5(x(:))]; 0175 0176 M.rand = @() trnsp(random(n, m)); 0177 0178 M.randvec = @(x) trnsp(randomvec(n, m, trnsp(x))); 0179 0180 M.lincomb = @matrixlincomb; 0181 0182 M.zerovec = @(x) trnsp(zeros(n, m)); 0183 0184 M.transp = @(x1, x2, d) M.proj(x2, d); 0185 0186 M.pairmean = @pairmean; 0187 function y = pairmean(x1, x2) 0188 y = trnsp(x1+x2); 0189 y = normalize_columns(y); 0190 y = trnsp(y); 0191 end 0192 0193 % vec returns a vector representation of an input tangent vector which 0194 % is represented as a matrix. mat returns the original matrix 0195 % representation of the input vector representation of a tangent 0196 % vector. vec and mat are thus inverse of each other. They are 0197 % furthermore isometries between a subspace of R^nm and the tangent 0198 % space at x. 0199 vect = @(X) X(:); 0200 M.vec = @(x, u_mat) vect(trnsp(u_mat)); 0201 M.mat = @(x, u_vec) trnsp(reshape(u_vec, [n, m])); 0202 M.vecmatareisometries = @() true; 0203 0204 end 0205 0206 % Given a matrix X, returns the same matrix but with each column scaled so 0207 % that they have unit 2-norm. 0208 function X = normalize_columns(X) 0209 % This is faster than norms(X, 2, 1) for small X, and as fast for large X. 0210 nrms = sqrt(sum(X.^2, 1)); 0211 X = bsxfun(@times, X, 1./nrms); 0212 end 0213 0214 % Orthogonal projection of the ambient vector H onto the tangent space at X 0215 function PXH = projection(X, H) 0216 0217 % Compute the inner product between each vector H(:, i) with its root 0218 % point X(:, i), that is, X(:, i)' * H(:, i). Returns a row vector. 0219 inners = sum(X.*H, 1); 0220 0221 % Subtract from H the components of the H(:, i)'s that are parallel to 0222 % the root points X(:, i). 0223 PXH = H - bsxfun(@times, X, inners); 0224 0225 % % Equivalent but slow code: 0226 % m = size(X, 2); 0227 % PXH = zeros(size(H)); 0228 % for i = 1 : m 0229 % PXH(:, i) = H(:, i) - X(:, i) * (X(:, i)'*H(:, i)); 0230 % end 0231 0232 end 0233 0234 % Uniform random sampling on the sphere. 0235 function x = random(n, m) 0236 0237 x = normalize_columns(randn(n, m)); 0238 0239 end 0240 0241 % Random normalized tangent vector at x. 0242 function d = randomvec(n, m, x) 0243 0244 d = randn(n, m); 0245 d = projection(x, d); 0246 d = d / norm(d(:)); 0247 0248 end