Manifold of n-by-n symmetric positive definite matrices with the
Bures-Wassterstein geometry.
function M = sympositivedefiniteBWfactory(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 Euclidean embedding space is the set of symmetric matrices of size n
with their usual trace inner product (Frobenius norm). In particular,
this means egrad and ehess (for Euclidean gradients and Hessians), if
implemented, must return symmetric matrices.
Please cite the Manopt paper as well as the research paper:
@article{malago2018wasserstein,
title={Wasserstein Riemannian geometry of Gaussian densities},
author={Malag{\`o}, Luigi and Montrucchio, Luigi and Pistone, Giovanni},
journal={Information Geometry},
volume={1},
number={2},
pages={137--179},
year={2018},
publisher={Springer}
}0001 function M = sympositivedefiniteBWfactory(n) 0002 % Manifold of n-by-n symmetric positive definite matrices with the 0003 % Bures-Wassterstein geometry. 0004 % 0005 % function M = sympositivedefiniteBWfactory(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 Euclidean embedding space is the set of symmetric matrices of size n 0012 % with their usual trace inner product (Frobenius norm). In particular, 0013 % this means egrad and ehess (for Euclidean gradients and Hessians), if 0014 % implemented, must return symmetric matrices. 0015 % 0016 % 0017 % Please cite the Manopt paper as well as the research paper: 0018 % @article{malago2018wasserstein, 0019 % title={Wasserstein Riemannian geometry of Gaussian densities}, 0020 % author={Malag{\`o}, Luigi and Montrucchio, Luigi and Pistone, Giovanni}, 0021 % journal={Information Geometry}, 0022 % volume={1}, 0023 % number={2}, 0024 % pages={137--179}, 0025 % year={2018}, 0026 % publisher={Springer} 0027 % } 0028 0029 % This file is part of Manopt: www.manopt.org. 0030 % Original author: Bamdev Mishra, January 23, 2020. 0031 % Contributors: 0032 % Change log: 0033 0034 symm = @(X) .5*(X+X'); 0035 0036 M.name = @() sprintf('Symmetric positive definite geometry of %dx%d matrices with the Bures-Wasserstein metric', n, n); 0037 0038 M.dim = @() n*(n+1)/2; 0039 0040 % Helpers to avoid computing full matrices simply to extract their trace 0041 vec = @(A) A(:); 0042 trAB = @(A, B) vec(A')'*vec(B); % = trace(A*B) 0043 0044 % Choice of the metric on the orthonormal space is motivated by the 0045 % symmetry present in the space. The metric on the positive definite 0046 % cone is the Bures-Wasserstein metric. 0047 M.inner = @myinner; 0048 function ip = myinner(X, eta, zeta) 0049 ip = 0.5*trAB(symm(lyapunov_symmetric(X, eta)), zeta); % BM: okay 0050 end 0051 0052 M.norm = @(X, eta) real(sqrt(myinner(X, eta, eta))); 0053 0054 M.dist = @mydist; 0055 function d = mydist(X, Y) 0056 Xhalf = sqrtm(X); 0057 d = real(sqrt(trace(X) + trace(Y) - 2*trace(symm(sqrtm(Xhalf*Y*Xhalf))))); 0058 end 0059 0060 M.typicaldist = @() sqrt(n*(n+1)/2); % BM: okay 0061 0062 M.egrad2rgrad = @egrad2rgrad; 0063 function eta = egrad2rgrad(X, eta) 0064 eta = 4*symm(eta*X); 0065 end 0066 0067 M.ehess2rhess = @ehess2rhess; 0068 function Hess = ehess2rhess(X, egrad, ehess, eta) 0069 % Directional derivatives of the Riemannian gradient 0070 Hess = 4*symm(ehess*X) + 4*symm(egrad*eta); 0071 0072 % Correction factor for the non-constant BW metric 0073 rgrad = egrad2rgrad(X, egrad); 0074 rgrad1 = lyapunov_symmetric(X, rgrad); 0075 eta1 = lyapunov_symmetric(X, eta); 0076 Hess = Hess ... 0077 - symm(rgrad1 * eta) ... 0078 - symm(rgrad * eta1) ... 0079 + 2*symm(X*symm(rgrad1 * eta1)); 0080 end 0081 0082 M.proj = @(X, eta) symm(eta); 0083 0084 M.tangent = M.proj; 0085 M.tangent2ambient = @(X, eta) eta; 0086 0087 M.exp = @exponential; 0088 function Y = exponential(X, eta, t) 0089 if nargin < 3 0090 t = 1.0; 0091 end 0092 teta = t*eta; 0093 teta1 = symm(lyapunov_symmetric(X, teta)); 0094 Y = X + teta + teta1*X*teta1; 0095 end 0096 0097 M.retr = @exponential; 0098 0099 function ABhalf = myhalf(A, B) 0100 Ahalf = sqrtm(A); 0101 ABhalf = (Ahalf*symm(sqrtm(Ahalf*B*Ahalf)))/Ahalf; 0102 end 0103 0104 M.log = @logarithm; 0105 function H = logarithm(X, Y) 0106 H = symm(myhalf(X, Y) + myhalf(Y, X) - 2*X); 0107 end 0108 0109 M.hash = @(X) ['z' hashmd5(X(:))]; 0110 0111 % Generate a random symmetric positive definite matrix following a 0112 % certain distribution. The particular choice of a distribution is of 0113 % course arbitrary, and specific applications might require different 0114 % ones. 0115 M.rand = @random; 0116 function X = random() 0117 D = diag(1+rand(n, 1)); 0118 [Q, R] = qr(randn(n)); %#ok 0119 X = Q*D*Q'; 0120 end 0121 0122 % Generate a uniformly random unit-norm tangent vector at X. 0123 M.randvec = @randomvec; 0124 function eta = randomvec(X) 0125 eta = symm(randn(n)); 0126 nrm = M.norm(X, eta); 0127 eta = eta / nrm; 0128 end 0129 0130 M.lincomb = @matrixlincomb; 0131 0132 M.zerovec = @(X) zeros(n); 0133 0134 % Poor man's vector transport: exploit the fact that all tangent spaces 0135 % are the set of symmetric matrices, so that the identity is a sort of 0136 % vector transport. It may perform poorly if the origin and target (X1 0137 % and X2) are far apart though. This should not be the case for typical 0138 % optimization algorithms, which perform small steps. 0139 M.transp = @(X1, X2, eta) eta; 0140 0141 % vec and mat are not isometries, because of the unusual inner metric. 0142 M.vec = @(X, U) U(:); 0143 M.mat = @(X, u) reshape(u, n, n); 0144 M.vecmatareisometries = @() false; 0145 0146 end