Manifold of m-by-n matrices with positive entries; scale invariant metric
function M = positivefactory(m)
function M = positivefactory(m, n)
A point X on the manifold M is represented as a matrix X of size mxn with
all individual entries real, strictly positive. By default, n = 1.
A tangent vector at X is represented as a matrix of the same size as X.
Entries of tangent vectors are free (in particular, not necessarily
positive.)
The Riemannian metric for each individual entry is the bi-invariant
metric for positive scalars, as a particular case of the bi-invariant
metric for positive definite matrices studied in Chapter 6 of the book
"Positive definite matrices" by Rajendra Bhatia,
Princeton University Press, 2007.
The Riemannian structure of M is obtained as the Cartesian product of the
geometry for mxn positive real numbers.
It should be stressed that matrices with one or more zero entries do not
belong to this manifold: they appear to be infinitely far away as a
result of the metric scaling like X.^(-1). Thus, if the solutions of an
optimization problem have entries equal to zero, these solutions are not
attainable on the manifold, which is likely to create serious numerical
issues. This geometry is best used when the solutions of the optimization
problem are indeed entry-wise positive, yet may have very different
scales (with some entries being very small, and some entries being very
large, relatively.)
See also: sympositivedefinitefactory0001 function M = positivefactory(m, n) 0002 % Manifold of m-by-n matrices with positive entries; scale invariant metric 0003 % 0004 % function M = positivefactory(m) 0005 % function M = positivefactory(m, n) 0006 % 0007 % A point X on the manifold M is represented as a matrix X of size mxn with 0008 % all individual entries real, strictly positive. By default, n = 1. 0009 % 0010 % A tangent vector at X is represented as a matrix of the same size as X. 0011 % Entries of tangent vectors are free (in particular, not necessarily 0012 % positive.) 0013 % 0014 % The Riemannian metric for each individual entry is the bi-invariant 0015 % metric for positive scalars, as a particular case of the bi-invariant 0016 % metric for positive definite matrices studied in Chapter 6 of the book 0017 % 0018 % "Positive definite matrices" by Rajendra Bhatia, 0019 % Princeton University Press, 2007. 0020 % 0021 % The Riemannian structure of M is obtained as the Cartesian product of the 0022 % geometry for mxn positive real numbers. 0023 % 0024 % It should be stressed that matrices with one or more zero entries do not 0025 % belong to this manifold: they appear to be infinitely far away as a 0026 % result of the metric scaling like X.^(-1). Thus, if the solutions of an 0027 % optimization problem have entries equal to zero, these solutions are not 0028 % attainable on the manifold, which is likely to create serious numerical 0029 % issues. This geometry is best used when the solutions of the optimization 0030 % problem are indeed entry-wise positive, yet may have very different 0031 % scales (with some entries being very small, and some entries being very 0032 % large, relatively.) 0033 % 0034 % See also: sympositivedefinitefactory 0035 0036 % This file is part of Manopt: www.manopt.org. 0037 % Original author: Bamdev Mishra, Dec 03, 2017. 0038 0039 if ~exist('n', 'var') || isempty(n) 0040 n = 1; 0041 end 0042 0043 M.name = @() sprintf('Element-wise positive %dx%d matrices', m, n); 0044 0045 M.dim = @() m*n; 0046 0047 % The metric is the scale invariant metric for scalars. 0048 M.inner = @myinner; 0049 function innerproduct = myinner(X, eta, zeta) 0050 innerproduct = (eta(:)./X(:))'*(zeta(:)./X(:)); 0051 end 0052 0053 M.norm = @(X, eta) sqrt(myinner(X, eta, eta)); 0054 0055 M.dist = @(X, Y) norm(log(Y./X), 'fro'); 0056 0057 M.typicaldist = @() sqrt(m*n); 0058 0059 M.egrad2rgrad = @egrad2rgrad; 0060 function eta = egrad2rgrad(X, eta) 0061 eta = X.*(eta).*X; 0062 end 0063 0064 M.ehess2rhess = @ehess2rhess; 0065 function Hess = ehess2rhess(X, egrad, ehess, eta) 0066 % Directional derivatives of the Riemannian gradient 0067 Hess = X.*(ehess).*X + 2*(eta.*(egrad).*X); 0068 0069 % Correction factor for the non-constant metric 0070 Hess = Hess - (eta.*(egrad).*X); 0071 end 0072 0073 % Since this manifold is an open subset of R^(nxm), the tangent space 0074 % at any X on M is all of R^(nxm). 0075 M.proj = @(X, eta) eta; 0076 0077 M.tangent = M.proj; 0078 M.tangent2ambient = @(X, eta) eta; 0079 0080 M.retr = @exponential; 0081 0082 M.exp = @exponential; 0083 function Y = exponential(X, eta, t) 0084 if nargin < 3 0085 t = 1.0; 0086 end 0087 % It is unclear whether this is the numerically most stable way to 0088 % implement this operation. If you run into trouble with this 0089 % factory, please get in touch on the forum. 0090 Y = (X.*(exp((t*eta)./X))); 0091 end 0092 0093 M.log = @logarithm; 0094 function H = logarithm(X, Y) 0095 % Same comment about numerical stability as for exp. 0096 H = (X.*(log(Y./X))); 0097 end 0098 0099 M.hash = @(X) ['z' hashmd5(X(:))]; 0100 0101 % Generate a random element-wise positive matrix following a 0102 % certain distribution. The particular choice of a distribution is of 0103 % course arbitrary, and specific applications might require different 0104 % ones. 0105 M.rand = @random; 0106 function X = random() 0107 X = exp(randn(m, n)); 0108 end 0109 0110 % Generate a uniformly random unit-norm tangent vector at X. 0111 M.randvec = @randomvec; 0112 function eta = randomvec(X) 0113 eta = randn(m, n).*X; 0114 nrm = M.norm(X, eta); 0115 eta = eta / nrm; 0116 end 0117 0118 M.lincomb = @matrixlincomb; 0119 0120 M.zerovec = @(X) zeros(m, n); 0121 0122 0123 M.transp = @(X1, X2, eta) eta; 0124 0125 % For reference, a proper vector transport is given here, following 0126 % work by Sra and Hosseini: "Conic geometric optimisation on the 0127 % manifold of positive definite matrices". 0128 % This is not used by default. To force the use of this transport, 0129 % execute "M.transp = M.paralleltransp;" on your M returned by the 0130 % present factory. 0131 M.paralleltransp = @parallel_transport; 0132 function zeta = parallel_transport(X, Y, eta) 0133 zeta = eta.*Y./X; 0134 end 0135 0136 % vec and mat are not isometries, because of the unusual inner metric. 0137 M.vec = @(X, U) U(:); 0138 M.mat = @(X, u) reshape(u, m, n); 0139 M.vecmatareisometries = @() true; 0140 0141 end