Home > manopt > manifolds > hyperbolic > poincareballfactory.m

# poincareballfactory

## PURPOSE

Factory for matrices whose columns live in the Poincare ball manifold

## SYNOPSIS

function M = poincareballfactory(k, n, gpuflag)

## DESCRIPTION

``` Factory for matrices whose columns live in the Poincare ball manifold

function M = poincareballfactory(k)
function M = poincareballfactory(k, n)

Manifold of k-by-n real matrices whose columns live in the Poincare ball.
By default, n = 1, which corresponds to a single Poincare ball.
The metric is such that each ball has constant sectional curvature -1.

This manifold is an open submanifold of R^{kxn}, so that tangent vectors
and vectors in the embedding space are represented as real matrices of
size kxn, without any restrictions. Points are likewise represented as
real matrices of size kxn such that each column has (Euclidean 2-norm)
strictly less than 1. The embedding space is endowed with its usual
and ehess2rhess thus expect to be given Euclidean gradients and Hessians.

Set gpuflag = true to have points, tangent vectors and ambient vectors
stored on the GPU. If so, computations can be done on the GPU directly.

## CROSS-REFERENCE INFORMATION

This function calls:
• factorygpuhelper Returns a manifold struct to optimize over unit-norm vectors or matrices.
• hashmd5 Computes the MD5 hash of input data.
• matrixlincomb Linear combination function for tangent vectors represented as matrices.
This function is called by:

## SOURCE CODE

```0001 function M = poincareballfactory(k, n, gpuflag)
0002 % Factory for matrices whose columns live in the Poincare ball manifold
0003 %
0004 % function M = poincareballfactory(k)
0005 % function M = poincareballfactory(k, n)
0006 %
0007 % Manifold of k-by-n real matrices whose columns live in the Poincare ball.
0008 % By default, n = 1, which corresponds to a single Poincare ball.
0009 % The metric is such that each ball has constant sectional curvature -1.
0010 %
0011 % This manifold is an open submanifold of R^{kxn}, so that tangent vectors
0012 % and vectors in the embedding space are represented as real matrices of
0013 % size kxn, without any restrictions. Points are likewise represented as
0014 % real matrices of size kxn such that each column has (Euclidean 2-norm)
0015 % strictly less than 1. The embedding space is endowed with its usual
0016 % Euclidean structure (with the trace inner product): the tools egrad2rgrad
0017 % and ehess2rhess thus expect to be given Euclidean gradients and Hessians.
0018 %
0019 % Set gpuflag = true to have points, tangent vectors and ambient vectors
0020 % stored on the GPU. If so, computations can be done on the GPU directly.
0021 %
0023
0024 % This file is part of Manopt: www.manopt.org.
0025 % Original author: Quentin Rebjock, Sep. 28, 2020
0026 % Contributors: NB
0027 % Change log:
0028
0029     if ~exist('n', 'var') || isempty(n)
0030         n = 1;
0031     end
0032
0033     if ~exist('gpuflag', 'var') || isempty(gpuflag)
0034         gpuflag = false;
0035     end
0036
0037     % If gpuflag is active, new arrays (e.g., via rand, randn, zeros, ones)
0038     % are created directly on the GPU; otherwise, they are created in the
0039     % usual way (in double precision).
0040     if gpuflag
0041         array_type = 'gpuArray';
0042     else
0043         array_type = 'double';
0044     end
0045
0046
0047     if n == 1
0048         M.name = @() sprintf('Poincare ball B_%d', k);
0049     else
0050         M.name = @() sprintf('Poincare ball B_%d^%d', k, n);
0051     end
0052
0053     M.dim = @() k * n;
0054
0055     M.conformal_factor = @(x) 2 ./ (1 - sum(x .* x, 1));
0056
0057     M.inner = @(x, u, v) sum(sum(u .* v, 1) .* (M.conformal_factor(x).^2));
0058
0059     M.norm = @(x, d) sqrt(M.inner(x, d, d));
0060
0061     M.dist = @dist;
0062     function d = dist(x, y)
0063         norms2x = sum(x .* x, 1);
0064         norms2y = sum(y .* y, 1);
0065         norms2diff = sum((x - y) .* (x - y), 1);
0066         d = sqrt(sum(acosh(1 + 2 * norms2diff ./ (1 - norms2x) ./ (1 - norms2y)) .^ 2));
0067     end
0068
0069     M.typicaldist = @() M.dim() / 8;
0070
0071     % Identity map since the embedding space is the tangent space.
0072     M.proj = @(x, d) d;
0073
0074     M.tangent = M.proj;
0075
0076     % The Poincare ball is not a Riemannian submanifold hence the Euclidean
0077     % gradient is not just a projection of the Euclidean gradient.
0080         factor = M.conformal_factor(x);
0082     end
0083
0084     M.ehess2rhess = @ehess2rhess;
0085     function rhess = ehess2rhess(x, egrad, ehess, u)
0086         factor = M.conformal_factor(x);
0087         rhess = ( u .* sum(egrad .* x, 1) - ...
0088                   egrad .* sum(u .* x, 1) - ...
0089                   x .* sum(u .* egrad, 1) + ...
0090                   ehess ./ factor ...
0091                 ) ./factor;
0092     end
0093
0095     function res = mobius_addition(x, y)
0096         sp = sum(x .* y, 1);
0097         norm2x = sum(x .* x, 1);
0098         norm2y = sum(y .* y, 1);
0099         res = ( x .* (1 + 2 .* sp + norm2y) + y .* (1 - norm2x) ) ...
0100                                        ./ (1 + 2 .* sp + norm2x .* norm2y);
0101     end
0102
0103     M.exp = @exponential;
0104     M.log = @logarithm;
0105
0106     M.retr = M.exp;
0107     M.invretr = M.log;
0108
0109     % This is not a parallel transport.
0110     M.transp = @(x1, x2, v) v;
0111
0112     M.hash = @(x) ['z' hashmd5(x(:))];
0113
0114     % Columns are sampled uniformly at random in the unit ball.
0115     M.rand = @() sample_ball_uniformly(k, n, array_type);
0116
0117     M.randvec = @randvec;
0118     function v = randvec(x)
0119         v = randn(k, n, array_type);
0120         v = v / M.norm(x, v);
0121     end
0122
0123     M.zerovec = @(x) zeros(k, n, array_type);
0124
0125     M.lincomb = @matrixlincomb;
0126
0127     M.pairmean = @pairmean;
0128     function y = pairmean(x1, x2)
0129         y = M.exp(x1, M.log(x1, x2) / 2);
0130     end
0131
0132     M.vec = @vec;
0133     function u_vec = vec(x, u_mat)
0134         u_vec = bsxfun(@times, u_mat, M.conformal_factor(x));
0135         u_vec = u_vec(:);
0136     end
0137     M.mat = @mat;
0138     function u_mat = mat(x, u_vec)
0139         u_mat = reshape(u_vec, [k, n]);
0140         u_mat = bsxfun(@times, u_mat, 1./M.conformal_factor(x));
0141     end
0142     M.vecmatareisometries = @() true;
0143
0144
0145     % Automatically convert a number of tools to support GPU.
0146     if gpuflag
0147         M = factorygpuhelper(M);
0148     end
0149
0150 end
0151
0152 function z = mobius_addition(x, y)
0153     inner = sum(x .* y, 1);
0154     norms2x = sum(x .* x, 1);
0155     norms2y = sum(y .* y, 1);
0156     z = ((1 + 2 * inner + norms2y) .* x + (1 - norms2x) .* y ) ./ (1 + 2 * inner + norms2x .* norms2y);
0157 end
0158
0159 % Exponential on the PoincarĂ© ball.
0160 function y = exponential(x, d, t)
0161     if nargin == 2
0162         % t = 1
0163         td = d;
0164     else
0165         td = t*d;
0166     end
0167
0168     normstd = vecnorm(td);
0169     factor = (1 - sum(x .* x, 1));
0170     % Avoid dividing by zero.
0171     w = td .* (tanh(normstd ./ factor) ./ (normstd + (normstd == 0)));
0173 end
0174
0175 function v = logarithm(x, y)