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.

- hashmd5 Computes the MD5 hash of input data.
- matrixlincomb Linear combination function for tangent vectors represented as matrices.

- positive_definite_karcher_mean Computes a Karcher mean of a collection of positive definite matrices.

- function eta = egrad2rgrad(X, eta)
- function Hess = ehess2rhess(X, egrad, ehess, eta)
- function Y = exponential(X, eta, t)
- function H = logarithm(X, Y)
- function X = random()
- function eta = randomvec(X)
- function zeta = parallel_transport(X, Y, eta)

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