Returns a manifold struct defining complex matrices w/ unit-norm columns. function M = obliquecomplexfactory(n, m) function M = obliquecomplexfactory(n, m, transposed) Oblique manifold: deals with complex matrices of size n x m such that each column has unit 2-norm, i.e., is a point on the unit sphere in C^n. The geometry is a product geometry of m unit spheres in C^n. For the metric, C^n is treated as R^(2n), so that the real part and imaginary parts are treated separately as 2n real coordinates. As such, the complex oblique manifold is a Riemannian submanifold of (R^2)^(n x m), with the usual metric <u, v> = real(u'*v). 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. In transposed form, a point Y is such that Y*Y' is a Hermitian, positive semidefinite matrix of size m and of rank at most n, such that all the diagonal entries are equal to 1. Note: obliquecomplexfactory(1, n, true) is equivalent to (but potentially slower than) complexcirclefactory(n). See also: spherecomplexfactory complexcirclefactory obliquefactory
0001 function M = obliquecomplexfactory(n, m, transposed) 0002 % Returns a manifold struct defining complex matrices w/ unit-norm columns. 0003 % 0004 % function M = obliquecomplexfactory(n, m) 0005 % function M = obliquecomplexfactory(n, m, transposed) 0006 % 0007 % Oblique manifold: deals with complex matrices of size n x m such that 0008 % each column has unit 2-norm, i.e., is a point on the unit sphere in C^n. 0009 % The geometry is a product geometry of m unit spheres in C^n. For the 0010 % metric, C^n is treated as R^(2n), so that the real part and imaginary 0011 % parts are treated separately as 2n real coordinates. As such, the complex 0012 % oblique manifold is a Riemannian submanifold of (R^2)^(n x m), with the 0013 % usual metric <u, v> = real(u'*v). 0014 % 0015 % If transposed is set to true (it is false by default), then the matrices 0016 % are transposed: a point Y on the manifold is a matrix of size m x n and 0017 % each row has unit 2-norm. It is the same geometry, just a different 0018 % representation. 0019 % 0020 % In transposed form, a point Y is such that Y*Y' is a Hermitian, positive 0021 % semidefinite matrix of size m and of rank at most n, such that all the 0022 % diagonal entries are equal to 1. 0023 % 0024 % Note: obliquecomplexfactory(1, n, true) is equivalent to (but potentially 0025 % slower than) complexcirclefactory(n). 0026 % 0027 % See also: spherecomplexfactory complexcirclefactory obliquefactory 0028 0029 % This file is part of Manopt: www.manopt.org. 0030 % Original author: Nicolas Boumal, Sep. 3, 2014. 0031 % Contributors: 0032 % Change log: 0033 % 0034 % Oct. 21, 2016 (NB) 0035 % Formatted for inclusion in Manopt release. 0036 % 0037 % July 20, 2017 (NB) 0038 % Distance function is now accurate for close-by points. See notes 0039 % inside the spherefactory file for details. Also improves distances 0040 % computation as part of the log function. 0041 0042 0043 if ~exist('transposed', 'var') || isempty(transposed) 0044 transposed = false; 0045 end 0046 0047 if transposed 0048 trnsp = @(X) X.'; 0049 else 0050 trnsp = @(X) X; 0051 end 0052 0053 M.name = @() sprintf('Complex oblique manifold COB(%d, %d)', n, m); 0054 0055 M.dim = @() (2*n-1)*m; 0056 0057 M.inner = @(x, d1, d2) real(d1(:)'*d2(:)); 0058 0059 M.norm = @(x, d) norm(d(:)); 0060 0061 M.dist = @(x, y) norm(real(2*asin(.5*sqrt(sum(trnsp(abs(x - y).^2), 1))))); 0062 0063 M.typicaldist = @() pi*sqrt(m); 0064 0065 M.proj = @(X, U) trnsp(projection(trnsp(X), trnsp(U))); 0066 0067 M.tangent = M.proj; 0068 0069 % For Riemannian submanifolds, converting a Euclidean gradient into a 0070 % Riemannian gradient amounts to an orthogonal projection. 0071 M.egrad2rgrad = M.proj; 0072 0073 M.ehess2rhess = @ehess2rhess; 0074 function rhess = ehess2rhess(X, egrad, ehess, U) 0075 X = trnsp(X); 0076 egrad = trnsp(egrad); 0077 ehess = trnsp(ehess); 0078 U = trnsp(U); 0079 0080 PXehess = projection(X, ehess); 0081 inners = sum(real(conj(X).*egrad), 1); 0082 rhess = PXehess - bsxfun(@times, U, inners); 0083 0084 rhess = trnsp(rhess); 0085 end 0086 0087 M.exp = @exponential; 0088 % Exponential on the complex oblique manifold 0089 function y = exponential(x, d, t) 0090 x = trnsp(x); 0091 d = trnsp(d); 0092 0093 if nargin == 2 0094 % t = 1; 0095 td = d; 0096 else 0097 td = t*d; 0098 end 0099 0100 nrm_td = sqrt(sum(real(td).^2 + imag(td).^2, 1)); 0101 0102 y = bsxfun(@times, x, cos(nrm_td)) + ... 0103 bsxfun(@times, td, sinxoverx(nrm_td)); 0104 0105 y = trnsp(y); 0106 end 0107 0108 M.log = @logarithm; 0109 function v = logarithm(x1, x2) 0110 x1 = trnsp(x1); 0111 x2 = trnsp(x2); 0112 0113 v = projection(x1, x2 - x1); 0114 dists = real(2*asin(.5*sqrt(sum(trnsp(abs(x - y).^2), 1)))); 0115 norms = sqrt(sum(real(v).^2 + imag(v).^2, 1)); 0116 factors = dists./norms; 0117 % For very close points, dists is almost equal to norms, but 0118 % because they are both almost zero, the division above can return 0119 % NaN's. To avoid that, we force those ratios to 1. 0120 factors(dists <= 1e-10) = 1; 0121 v = bsxfun(@times, v, factors); 0122 0123 v = trnsp(v); 0124 end 0125 0126 M.retr = @retraction; 0127 % Retraction on the oblique manifold 0128 function y = retraction(x, d, t) 0129 x = trnsp(x); 0130 d = trnsp(d); 0131 0132 if nargin < 3 0133 td = d; 0134 else 0135 td = t*d; 0136 end 0137 0138 y = normalize_columns(x + td); 0139 0140 y = trnsp(y); 0141 end 0142 0143 M.hash = @(x) ['z' hashmd5([real(x(:)) ; imag(x(:))])]; 0144 0145 M.rand = @() trnsp(random(n, m)); 0146 0147 M.randvec = @(x) trnsp(randomvec(n, m, trnsp(x))); 0148 0149 M.lincomb = @matrixlincomb; 0150 0151 M.zerovec = @(x) trnsp(zeros(n, m)); 0152 0153 M.transp = @(x1, x2, d) M.proj(x2, d); 0154 0155 M.pairmean = @pairmean; 0156 function y = pairmean(x1, x2) 0157 y = trnsp(x1+x2); 0158 y = normalize_columns(y); 0159 y = trnsp(y); 0160 end 0161 0162 % vec returns a vector representation of an input tangent vector which 0163 % is represented as a matrix. mat returns the original matrix 0164 % representation of the input vector representation of a tangent 0165 % vector. vec and mat are thus inverse of each other. They are 0166 % furthermore isometries between a subspace of R^2nm and the tangent 0167 % space at x. 0168 vect = @(X) X(:); 0169 M.vec = @(x, u_mat) [vect(real(trnsp(u_mat))) ; ... 0170 vect(imag(trnsp(u_mat)))]; 0171 M.mat = @(x, u_vec) trnsp(reshape(u_vec(1:(n*m)), [n, m])) + ... 0172 1i*trnsp(reshape(u_vec((n*m+1):end), [n, m])); 0173 M.vecmatareisometries = @() true; 0174 0175 end 0176 0177 % Given a matrix X, returns the same matrix but with each column scaled so 0178 % that they have unit 2-norm. 0179 function X = normalize_columns(X) 0180 norms = sqrt(sum(real(X).^2 + imag(X).^2, 1)); 0181 X = bsxfun(@times, X, 1./norms); 0182 end 0183 0184 % Orthogonal projection of the ambient vector H onto the tangent space at X 0185 function PXH = projection(X, H) 0186 0187 % Compute the inner product between each vector H(:, i) with its root 0188 % point X(:, i), that is, real(X(:, i)' * H(:, i)). 0189 % Returns a row vector. 0190 inners = real(sum(conj(X).*H, 1)); 0191 0192 % Subtract from H the components of the H(:, i)'s that are parallel to 0193 % the root points X(:, i). 0194 PXH = H - bsxfun(@times, X, inners); 0195 0196 end 0197 0198 % Uniform random sampling on the sphere. 0199 function x = random(n, m) 0200 0201 x = normalize_columns(randn(n, m) + 1i*randn(n, m)); 0202 0203 end 0204 0205 % Random normalized tangent vector at x. 0206 function d = randomvec(n, m, x) 0207 0208 d = randn(n, m) + 1i*randn(n, m); 0209 d = projection(x, d); 0210 d = d / norm(d(:)); 0211 0212 end