Home > manopt > manifolds > symfixedrank > sympositivedefinitefactory.m

sympositivedefinitefactory

PURPOSE ^

Manifold of n-by-n symmetric positive definite matrices with

SYNOPSIS ^

function M = sympositivedefinitefactory(n)

DESCRIPTION ^

 Manifold of n-by-n symmetric positive definite matrices with
 the bi-invariant geometry.

 function M = sympositivedefinitefactory(n)

 A point X on the manifold is represented as a symmetric positive definite
 matrix X (nxn). Tangent vectors are symmetric matrices of the same size
 (but not necessarily definite).

 The Riemannian metric is the bi-invariant metric, described notably in
 Chapter 6 of the 2007 book "Positive definite matrices"
 by Rajendra Bhatia, Princeton University Press.


 The retraction / exponential map involves expm (the matrix exponential).
 If too large a vector is retracted / exponentiated (e.g., a solver tries
 to make too big a step), this may result in NaN's in the returned point,
 which most likely would lead to NaN's in the cost / gradient / ... and
 will result in failure of the optimization. For trustregions, this can be
 controlled by setting options.Delta0 and options.Delta_bar, to prevent
 too large steps.


 Note also that many of the functions involve solving linear systems in X
 (a point on the manifold), taking matrix exponentals and logarithms, etc.
 It could therefore be beneficial to do some precomputation on X (an
 eigenvalue decomposition for example) and store both X and the
 preprocessing in a structure. This would require modifying the present
 factory to work with such structures to represent both points and tangent
 vectors. We omit this in favor of simplicity, but it may be good to keep
 this in mind if efficiency becomes an issue in your application.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function M = sympositivedefinitefactory(n)
0002 % Manifold of n-by-n symmetric positive definite matrices with
0003 % the bi-invariant geometry.
0004 %
0005 % function M = sympositivedefinitefactory(n)
0006 %
0007 % A point X on the manifold is represented as a symmetric positive definite
0008 % matrix X (nxn). Tangent vectors are symmetric matrices of the same size
0009 % (but not necessarily definite).
0010 %
0011 % The Riemannian metric is the bi-invariant metric, described notably in
0012 % Chapter 6 of the 2007 book "Positive definite matrices"
0013 % by Rajendra Bhatia, Princeton University Press.
0014 %
0015 %
0016 % The retraction / exponential map involves expm (the matrix exponential).
0017 % If too large a vector is retracted / exponentiated (e.g., a solver tries
0018 % to make too big a step), this may result in NaN's in the returned point,
0019 % which most likely would lead to NaN's in the cost / gradient / ... and
0020 % will result in failure of the optimization. For trustregions, this can be
0021 % controlled by setting options.Delta0 and options.Delta_bar, to prevent
0022 % too large steps.
0023 %
0024 %
0025 % Note also that many of the functions involve solving linear systems in X
0026 % (a point on the manifold), taking matrix exponentals and logarithms, etc.
0027 % It could therefore be beneficial to do some precomputation on X (an
0028 % eigenvalue decomposition for example) and store both X and the
0029 % preprocessing in a structure. This would require modifying the present
0030 % factory to work with such structures to represent both points and tangent
0031 % vectors. We omit this in favor of simplicity, but it may be good to keep
0032 % this in mind if efficiency becomes an issue in your application.
0033 
0034 % This file is part of Manopt: www.manopt.org.
0035 % Original author: Bamdev Mishra, August 29, 2013.
0036 % Contributors: Nicolas Boumal
0037 % Change log:
0038 %
0039 %   March 5, 2014 (NB)
0040 %       There were a number of mistakes in the code owing to the tacit
0041 %       assumption that if X and eta are symmetric, then X\eta is
0042 %       symmetric too, which is not the case. See discussion on the Manopt
0043 %       forum started on Jan. 19, 2014. Functions norm, dist, exp and log
0044 %       were modified accordingly. Furthermore, they only require matrix
0045 %       inversion (as well as matrix log or matrix exp), not matrix square
0046 %       roots or their inverse.
0047 %
0048 %   July 28, 2014 (NB)
0049 %       The dim() function returned n*(n-1)/2 instead of n*(n+1)/2.
0050 %       Implemented proper parallel transport from Sra and Hosseini (not
0051 %       used by default).
0052 %       Also added symmetrization in exp and log (to be sure).
0053 %
0054 %   April 3, 2015 (NB):
0055 %       Replaced trace(A*B) by a faster equivalent that does not compute
0056 %       the whole product A*B, for inner product, norm and distance.
0057     
0058     symm = @(X) .5*(X+X');
0059     
0060     M.name = @() sprintf('Symmetric positive definite geometry of %dx%d matrices', n, n);
0061     
0062     M.dim = @() n*(n+1)/2;
0063     
0064     % Helpers to avoid computing full matrices simply to extract their trace
0065     vec     = @(A) A(:);
0066     trinner = @(A, B) vec(A')'*vec(B);  % = trace(A*B)
0067     trnorm  = @(A) sqrt(trinner(A, A)); % = sqrt(trace(A^2))
0068     
0069     % Choice of the metric on the orthonormal space is motivated by the
0070     % symmetry present in the space. The metric on the positive definite
0071     % cone is its natural bi-invariant metric.
0072     % The result is equal to: trace( (X\eta) * (X\zeta) )
0073     M.inner = @(X, eta, zeta) trinner(X\eta, X\zeta);
0074     
0075     % Notice that X\eta is *not* symmetric in general.
0076     % The result is equal to: sqrt(trace((X\eta)^2))
0077     % There should be no need to take the real part, but rounding errors
0078     % may cause a small imaginary part to appear, so we discard it.
0079     M.norm = @(X, eta) real(trnorm(X\eta));
0080     
0081     % Same here: X\Y is not symmetric in general.
0082     % Same remark about taking the real part.
0083     M.dist = @(X, Y) real(trnorm(real(logm(X\Y))));
0084     
0085     
0086     M.typicaldist = @() sqrt(n*(n+1)/2);
0087     
0088     
0089     M.egrad2rgrad = @egrad2rgrad;
0090     function eta = egrad2rgrad(X, eta)
0091         eta = X*symm(eta)*X;
0092     end
0093     
0094     
0095     M.ehess2rhess = @ehess2rhess;
0096     function Hess = ehess2rhess(X, egrad, ehess, eta)
0097         % Directional derivatives of the Riemannian gradient
0098         Hess = X*symm(ehess)*X + 2*symm(eta*symm(egrad)*X);
0099         
0100         % Correction factor for the non-constant metric
0101         Hess = Hess - symm(eta*symm(egrad)*X);
0102     end
0103     
0104     
0105     M.proj = @(X, eta) symm(eta);
0106     
0107     M.tangent = M.proj;
0108     M.tangent2ambient = @(X, eta) eta;
0109     
0110     M.retr = @exponential;
0111     
0112     M.exp = @exponential;
0113     function Y = exponential(X, eta, t)
0114         if nargin < 3
0115             t = 1.0;
0116         end
0117         % The symm() and real() calls are mathematically not necessary but
0118         % are numerically necessary.
0119         Y = symm(X*real(expm(X\(t*eta))));
0120     end
0121     
0122     M.log = @logarithm;
0123     function H = logarithm(X, Y)
0124         % Same remark regarding the calls to symm() and real().
0125         H = symm(X*real(logm(X\Y)));
0126     end
0127     
0128     M.hash = @(X) ['z' hashmd5(X(:))];
0129     
0130     % Generate a random symmetric positive definite matrix following a
0131     % certain distribution. The particular choice of a distribution is of
0132     % course arbitrary, and specific applications might require different
0133     % ones.
0134     M.rand = @random;
0135     function X = random()
0136         D = diag(1+rand(n, 1));
0137         [Q, R] = qr(randn(n)); %#ok<NASGU>
0138         X = Q*D*Q';
0139     end
0140     
0141     % Generate a uniformly random unit-norm tangent vector at X.
0142     M.randvec = @randomvec;
0143     function eta = randomvec(X)
0144         eta = symm(randn(n));
0145         nrm = M.norm(X, eta);
0146         eta = eta / nrm;
0147     end
0148     
0149     M.lincomb = @matrixlincomb;
0150     
0151     M.zerovec = @(X) zeros(n);
0152     
0153     % Poor man's vector transport: exploit the fact that all tangent spaces
0154     % are the set of symmetric matrices, so that the identity is a sort of
0155     % vector transport. It may perform poorly if the origin and target (X1
0156     % and X2) are far apart though. This should not be the case for typical
0157     % optimization algorithms, which perform small steps.
0158     M.transp = @(X1, X2, eta) eta;
0159     
0160     % For reference, a proper vector transport is given here, following
0161     % work by Sra and Hosseini: "Conic geometric optimisation on the
0162     % manifold of positive definite matrices", to appear in SIAM J. Optim.
0163     % in 2015; also available here: http://arxiv.org/abs/1312.1039
0164     % This will not be used by default. To force the use of this transport,
0165     % execute "M.transp = M.paralleltransp;" on your M returned by the
0166     % present factory.
0167     M.paralleltransp = @parallel_transport;
0168     function zeta = parallel_transport(X, Y, eta)
0169         E = sqrtm((Y/X));
0170         zeta = E*eta*E';
0171     end
0172     
0173     % vec and mat are not isometries, because of the unusual inner metric.
0174     M.vec = @(X, U) U(:);
0175     M.mat = @(X, u) reshape(u, n, n);
0176     M.vecmatareisometries = @() false;
0177     
0178 end

Generated on Sat 12-Nov-2016 14:11:22 by m2html © 2005