Home > manopt > manifolds > oblique > obliquecomplexfactory.m

obliquecomplexfactory

PURPOSE ^

Returns a manifold struct defining complex matrices w/ unit-norm columns.

SYNOPSIS ^

function M = obliquecomplexfactory(n, m, transposed)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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 improvies 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, sin(nrm_td) ./ nrm_td);
0104         
0105         % For those columns where the step is 0, replace y by x
0106         exclude = (nrm_td == 0);
0107         y(:, exclude) = x(:, exclude);
0108 
0109         y = trnsp(y);
0110     end
0111 
0112     M.log = @logarithm;
0113     function v = logarithm(x1, x2)
0114         x1 = trnsp(x1);
0115         x2 = trnsp(x2);
0116         
0117         v = projection(x1, x2 - x1);
0118         dists = real(2*asin(.5*sqrt(sum(trnsp(abs(x - y).^2), 1))));
0119         norms = sqrt(sum(real(v).^2 + imag(v).^2, 1));
0120         factors = dists./norms;
0121         % For very close points, dists is almost equal to norms, but
0122         % because they are both almost zero, the division above can return
0123         % NaN's. To avoid that, we force those ratios to 1.
0124         factors(dists <= 1e-10) = 1;
0125         v = bsxfun(@times, v, factors);
0126         
0127         v = trnsp(v);
0128     end
0129 
0130     M.retr = @retraction;
0131     % Retraction on the oblique manifold
0132     function y = retraction(x, d, t)
0133         x = trnsp(x);
0134         d = trnsp(d);
0135         
0136         if nargin < 3
0137             td = d;
0138         else
0139             td = t*d;
0140         end
0141 
0142         y = normalize_columns(x + td);
0143         
0144         y = trnsp(y);
0145     end
0146 
0147     M.hash = @(x) ['z' hashmd5([real(x(:)) ; imag(x(:))])];
0148     
0149     M.rand = @() trnsp(random(n, m));
0150     
0151     M.randvec = @(x) trnsp(randomvec(n, m, trnsp(x)));
0152     
0153     M.lincomb = @matrixlincomb;
0154     
0155     M.zerovec = @(x) trnsp(zeros(n, m));
0156     
0157     M.transp = @(x1, x2, d) M.proj(x2, d);
0158     
0159     M.pairmean = @pairmean;
0160     function y = pairmean(x1, x2)
0161         y = trnsp(x1+x2);
0162         y = normalize_columns(y);
0163         y = trnsp(y);
0164     end
0165 
0166     % vec returns a vector representation of an input tangent vector which
0167     % is represented as a matrix. mat returns the original matrix
0168     % representation of the input vector representation of a tangent
0169     % vector. vec and mat are thus inverse of each other. They are
0170     % furthermore isometries between a subspace of R^2nm and the tangent
0171     % space at x.
0172     vect = @(X) X(:);
0173     M.vec = @(x, u_mat) [vect(real(trnsp(u_mat))) ; ...
0174                          vect(imag(trnsp(u_mat)))];
0175     M.mat = @(x, u_vec)    trnsp(reshape(u_vec(1:(n*m)),     [n, m])) + ...
0176                         1i*trnsp(reshape(u_vec((n*m+1):end), [n, m]));
0177     M.vecmatareisometries = @() true;
0178 
0179 end
0180 
0181 % Given a matrix X, returns the same matrix but with each column scaled so
0182 % that they have unit 2-norm.
0183 function X = normalize_columns(X)
0184     norms = sqrt(sum(real(X).^2 + imag(X).^2, 1));
0185     X = bsxfun(@times, X, 1./norms);
0186 end
0187 
0188 % Orthogonal projection of the ambient vector H onto the tangent space at X
0189 function PXH = projection(X, H)
0190 
0191     % Compute the inner product between each vector H(:, i) with its root
0192     % point X(:, i), that is, real(X(:, i)' * H(:, i)).
0193     % Returns a row vector.
0194     inners = real(sum(conj(X).*H, 1));
0195     
0196     % Subtract from H the components of the H(:, i)'s that are parallel to
0197     % the root points X(:, i).
0198     PXH = H - bsxfun(@times, X, inners);
0199 
0200 end
0201 
0202 % Uniform random sampling on the sphere.
0203 function x = random(n, m)
0204 
0205     x = normalize_columns(randn(n, m) + 1i*randn(n, m));
0206 
0207 end
0208 
0209 % Random normalized tangent vector at x.
0210 function d = randomvec(n, m, x)
0211 
0212     d = randn(n, m) + 1i*randn(n, m);
0213     d = projection(x, d);
0214     d = d / norm(d(:));
0215 
0216 end

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