Home > manopt > manifolds > symfixedrank > sympositivedefiniteBWfactory.m

# sympositivedefiniteBWfactory

## PURPOSE

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

## SYNOPSIS

function M = sympositivedefiniteBWfactory(n)

## DESCRIPTION

``` 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,
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}
}```

## CROSS-REFERENCE INFORMATION

This function calls:
• hashmd5 Computes the MD5 hash of input data.
• lyapunov_symmetric Solves AX + XA = C when A = A', as a pseudo-inverse.
• matrixlincomb Linear combination function for tangent vectors represented as matrices.
This function is called by:

## SOURCE CODE

```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
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
0075         eta1 = lyapunov_symmetric(X, eta);
0076         Hess = Hess ...
0077             - symm(rgrad1 * eta) ...
0078             - symm(rgrad * 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```

Generated on Fri 30-Sep-2022 13:18:25 by m2html © 2005