geomstats.numerics package#

Submodules#

geomstats.numerics.bvp module#

Boundary value problem solvers implementation.

class geomstats.numerics.bvp.ScipySolveBVP(tol=0.001, max_nodes=1000, bc_tol=None, save_result=False)[source]#

Bases: object

Wrapper for scipy.integrate.solve_bvp.

integrate(fun, bc, x, y, fun_jac=None, bc_jac=None)[source]#

Solve a boundary value problem for a system of ODEs.

geomstats.numerics.finite_differences module#

Finite differences machinery.

geomstats.numerics.finite_differences.centered_difference(array, delta=None, axis=-1, endpoints=False)[source]#

Centered difference in a Euclidean space.

Points live in R^m, but are a k-dim embedding (e.g. a curve). Assumes points are in correspondence for cases higher than dim=1.

Parameters:

  • array (array-like) – Values of a function.

  • delta (float) – Spacing between points.

  • axis (int) – Axis in which perform the difference. Must be given backwards.

  • endpoints (bool) – If True, endpoints are computed by backward and forward differences, respectively.

Returns:

centered_diff (array-like) – Same shape as array.

geomstats.numerics.finite_differences.forward_difference(array, delta=None, axis=-1)[source]#

Forward difference in a Euclidean space.

Points live in R^m, but are a k-dim embedding (e.g. a curve). Assumes points are in correspondence for cases higher than dim=1.

Parameters:

  • array (array-like) – Values of a function.

  • delta (float) – Spacing between points.

  • axis (int) – Axis in which perform the difference. Must be given backwards.

Returns:

forward_diff (array-like) – Shape in the specified axis reduces by one.

geomstats.numerics.finite_differences.second_centered_difference(array, delta=None, axis=-1)[source]#

Second centered difference in a Euclidean space.

Points live in R^m, but are a k-dim embedding (e.g. a curve). Assumes points are in correspondence for cases higher than dim=1.

Parameters:

  • array (array-like) – Values of a function.

  • delta (float) – Spacing between points.

  • axis (int) – Axis in which perform the difference. Must be given backwards.

Returns:

second_centered_diff (array-like) – Shape in the specified axis reduces by two (endpoints).

geomstats.numerics.geodesic module#

Geodesic solvers implementation.

class geomstats.numerics.geodesic.ExpODESolver(space, integrator=None)[source]#

Bases: ExpSolver

Geodesic initial value problem solver.

Integrate geodesic equation.

Parameters:

  • space (Manifold) – Equipped manifold.

  • integrator (ODEIVPIntegrator) – Instance of ODEIVP integrator.

exp(tangent_vec, base_point)[source]#

Exponential map.

Parameters:

  • tangent_vec (array-like, shape=[…, *space.shape]) – Tangent vector at the base point.

  • base_point (array-like, shape=[…, *space.shape]) – Point on the manifold.

Returns:

end_point (array-like, shape=[…, *space.shape]) – Point on the manifold.

geodesic_ivp(tangent_vec, base_point)[source]#

Geodesic curve for initial value problem.

Parameters:

  • tangent_vec (array-like, shape=[…, *space.shape]) – Tangent vector at the base point.

  • base_point (array-like, shape=[…, *space.shape]) – Point on the manifold.

Returns:

path (callable) – Time parametrized geodesic curve. f(t).

property integrator#

An instance of ODEIVPIntegrator.

class geomstats.numerics.geodesic.ExpSolver(solves_ivp=False)[source]#

Bases: ABC

Abstract class for geodesic initial value problem solvers.

Parameters:

solves_ivp (bool) – Informs if solver is able to solve for geodesic at different t.

abstract exp(tangent_vec, base_point)[source]#

Exponential map.

Parameters:

  • tangent_vec (array-like, shape=[…, dim]) – Tangent vector at the base point.

  • base_point (array-like, shape=[…, dim]) – Point on the manifold.

Returns:

end_point (array-like, shape=[…, dim]) – Point on the manifold.

geodesic_ivp(tangent_vec, base_point)[source]#

Geodesic curve for initial value problem.

Parameters:

  • tangent_vec (array-like, shape=[…, dim]) – Tangent vector at the base point.

  • base_point (array-like, shape=[…, dim]) – Point on the manifold.

Returns:

path (callable) – Time parametrized geodesic curve. f(t).

class geomstats.numerics.geodesic.LogODESolver(space, n_nodes=10, integrator=None, initialization=None, use_jac=True)[source]#

Bases: _LogBatchMixins, LogSolver

Geodesic boundary value problem using an ODE solver.

Parameters:

  • space (Manifold) – Equipped manifold.

  • n_nodes (Number of mesh nodes.)

  • integrator (ScipySolveBVP) – Instance of ScipySolveBVP.

  • initialization (callable) – Function to provide initial solution. f( point, base_point). Defaults to linear initialization.

geodesic_bvp(point, base_point)[source]#

Geodesic curve for boundary value problem.

Parameters:

  • end_point (array-like, shape=[…, dim]) – Point on the manifold.

  • base_point (array-like, shape=[…, dim]) – Point on the manifold.

Returns:

path (callable) – Time parametrized geodesic curve. f(t). 0 <= t <= 1.

class geomstats.numerics.geodesic.LogShootingSolver(space, optimizer=None, initialization=None, flatten=True)[source]#

Bases: object

Geodesic boundary value problem solver using shooting.

Parameters:

  • space (Manifold) – Equipped manifold.

  • optimizer (ScipyMinimize) – Instance of ScipyMinimize.

  • initialization (callable) – Function to provide initial solution. f(point, base_point). Defaults to linear initialization.

  • flatten (bool) – If True, the optimization problem is solved together for all the batch points.

class geomstats.numerics.geodesic.LogSolver(solves_bvp=False)[source]#

Bases: ABC

Abstract class for geodesic boundary value problem solvers.

Parameters:

solves_bvp (bool) – Informs if solver is able to solve for geodesic at different t.

geodesic_bvp(point, base_point)[source]#

Geodesic curve for boundary value problem.

Parameters:

  • end_point (array-like, shape=[…, dim]) – Point on the manifold.

  • base_point (array-like, shape=[…, dim]) – Point on the manifold.

Returns:

path (callable) – Time parametrized geodesic curve. f(t).

abstract log(space, point, base_point)[source]#

Logarithm map.

Parameters:

  • end_point (array-like, shape=[…, dim]) – Point on the manifold.

  • base_point (array-like, shape=[…, dim]) – Point on the manifold.

Returns:

tangent_vec (array-like, shape=[…, dim]) – Tangent vector at the base point.

class geomstats.numerics.geodesic.MultiresPathStraightening(space, path_energy=None, n_nodes=None, n_levels=3, optimizer=None, initialization=None, symmetric=False, early_stop=False)[source]#

Bases: _DiscreteGeodesicBVPBatchMixins, PathBasedLogSolver

Geodesic boundary value problem with multiresolution path straightening.

Parameters:

  • space (Manifold) – Equipped manifold.

  • path_energy (callable) – Method to compute Riemannian path energy.

  • n_nodes (list) – Number of path discretization points at each resolution.

  • n_levels (int) – Number of resolutions to use. Sets number of nodes following a sequence. Ignored if n_nodes is not None.

  • optimizer (ScipyMinimize) – An optimizer to solve path energy minimization problem.

  • initialization (callable) – A method to get initial guess for optimization.

  • symmetric (bool) – If to use a symmetrized version of the energy.

  • early_stop (bool) – If to stop at a resolution if previous resolution was solved with one iteration.

class geomstats.numerics.geodesic.PathBasedLogSolver(space)[source]#

Bases: LogSolver, ABC

A geodesic BVP solver based on finding a discrete geodesic path.

Parameters:

space (Manifold) – Equipped manifold.

abstract discrete_geodesic_bvp(point, base_point)[source]#

Solve boundary value problem (BVP).

Given an initial point and an end point, solve the geodesic equation via minimizing the Riemannian path energy.

Parameters:

  • point (array-like, shape=[…, *point_shape]) – Point on the manifold.

  • base_point (array-like, shape=[…, *point_shape]) – Point on the manifold.

Returns:

discr_geod_path (array-like, shape=[…, n_nodes, *point_shape]) – Discrete geodesic.

geodesic_bvp(point, base_point)[source]#

Geodesic curve for boundary value problem.

Parameters:

  • end_point (array-like, shape=[…, *point_shape]) – Point on the manifold.

  • base_point (array-like, shape=[…, *point_shape]) – Point on the manifold.

Returns:

path (callable) – Time parametrized geodesic curve. f(t).

log(point, base_point)[source]#

Logarithm map.

Parameters:

  • end_point (array-like, shape=[…, *point_shape]) – Point on the manifold.

  • base_point (array-like, shape=[…, *point_shape]) – Point on the manifold.

Returns:

tangent_vec (array-like, shape=[…, *point_shape]) – Tangent vector at the base point.

path_n_nodes(path)[source]#

Get number of nodes of a path.

Parameters:

discr_path (array-like, shape=[…, n_nodes, *point_shape]) – Discrete path.

Returns:

n_nodes (int) – Number of path nodes.

class geomstats.numerics.geodesic.PathStraightening(space, path_energy=None, n_nodes=100, optimizer=None, initialization=None, symmetric=False)[source]#

Bases: _DiscreteGeodesicBVPBatchMixins, PathBasedLogSolver

Geodesic boundary value problem with path-straightening.

Parameters:

  • space (Manifold) – Equipped manifold.

  • path_energy (callable) – Method to compute Riemannian path energy.

  • n_nodes (int) – Number of path discretization points.

  • optimizer (ScipyMinimize) – An optimizer to solve path energy minimization problem.

  • initialization (callable) – A method to get initial guess for optimization.

  • symmetric (bool) – If to use a symmetrized version of the energy.

References

[HSKCB2022]

“Elastic shape analysis of surfaces with second-order Sobolev metrics: a comprehensive numerical framework”. arXiv:2204.04238 [cs.CV], 25 Sep 2022

geomstats.numerics.interpolation module#

Interpolation machinery.

class geomstats.numerics.interpolation.Interpolator[source]#

Bases: ABC

Abstract class for interpolator.

abstract interpolate(t)[source]#

Interpolate data.

Parameters:

t (array-like, shape=[n_time]) – Interpolation time.

Returns:

point (array-like, shape=[…, n_time, *point_shape])

class geomstats.numerics.interpolation.LinearInterpolator1D(times, data, point_ndim=1)[source]#

Bases: _LinearInterpolator1D

A 1D linear interpolator.

Assumes interpolation occurs in the unit interval.

Parameters:

  • times (array-like, [n_times]) – Times. Must be sorted.

  • data (array-like, […, *point_shape])

  • point_ndim (int) – Dimension of point.

class geomstats.numerics.interpolation.UniformUnitIntervalLinearInterpolator(data, point_ndim=1)[source]#

Bases: _LinearInterpolator1D

A 1D linear interpolator.

Assumes interpolation occurs in the unit interval and data is uniformly sampled.

Parameters:

  • data (array-like, […, *point_shape])

  • point_ndim (int) – Dimension of point.

geomstats.numerics.ivp module#

Initial value problem solvers implementation.

class geomstats.numerics.ivp.GSIVPIntegrator(n_steps=10, step_type='euler', save_result=False)[source]#

Bases: ODEIVPIntegrator

In-house ODE integrator.

Parameters:

  • n_steps (int) – Number of steps to perform.

  • step_type (str) – Type of integration step. Possible values are euler, rk2, rk4.

  • save_result (bool) – If True, result is stored after calling integrate or integrate_t.

integrate(force, initial_state, end_time=1.0)[source]#

Integrate force.

Parameters:

  • force (callable) – Function to integrate: f(state, t).

  • initial_state (array-like, shape=[…, n_vars, *point_shape]) – Initial state.

  • end_time (float or None) – Integration end time.

Returns:

result (OdeResult)

property step_type#

Integrator step type.

class geomstats.numerics.ivp.ODEIVPIntegrator(save_result=False, tchosen=False)[source]#

Bases: ABC

Abstract class for ode ivp solvers.

Parameters:

  • save_result (bool) – If True, result is stored after calling integrate or integrate_t.

  • tchosen (bool) – Informs about ability to solve at chosen times. If False, then does not implement integrate_t.

abstract integrate(force, initial_state, end_time)[source]#

Integrate force.

Parameters:

  • force (callable) – Function to integrate: f(state, t).

  • initial_state (array-like, shape=[…, n_vars, *point_shape]) – Initial state.

  • end_time (float or None) – Integration end time.

Returns:

result (OdeResult)

integrate_t(force, initial_state, t_eval)[source]#

Integrate force while choosing evaluating points.

Parameters:

  • force (callable) – Function to integrate: f(state, t).

  • initial_state (array-like, shape=[…, n_vars, *point_shape]) – Initial state.

  • t_eval (array-like) – Times at which to store the computed solution.

Returns:

result (OdeResult)

class geomstats.numerics.ivp.OdeResult[source]#

Bases: OptimizeResult

Bunch object (follows scipy).

Its purposes is to homogenize output of different integrators.

get_last_y()[source]#

Get value for last y.

Allows to have y represented as an gs.array or list[gs.array] (latter for cases where they have different shapes).

Assumes last t is the same.

class geomstats.numerics.ivp.ScipySolveIVP(method='RK45', save_result=False, point_ndim=1, **options)[source]#

Bases: ODEIVPIntegrator

Wrapper for scipy.integrate.solve_ivp.

Check https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html for additional options.

Parameters:

  • method (str) – Integration method.

  • save_result (bool) – If True, result is stored after calling integrate or integrate_t.

  • point_ndim = int – Dimension of array representing a point in the space.

integrate(force, initial_state, end_time=1.0)[source]#

Integrate force.

Parameters:

  • force (callable) – Function to integrate: f(state, t).

  • initial_state (array-like, shape=[…, n_vars, *point_shape]) – Initial state.

  • end_time (float or None) – Integration end time.

Returns:

result (OdeResult)

integrate_t(force, initial_state, t_eval)[source]#

Integrate force at t_eval points.

Parameters:

  • force (callable) – Function to integrate: f(state, t).

  • initial_state (array-like, shape=[…, n_vars, *point_shape]) – Initial state.

  • t_eval (array-like) – Times at which to store the computed solution.

Returns:

result (OdeResult)

geomstats.numerics.optimization module#

Optimizers implementations.

class geomstats.numerics.optimization.NewtonMethod(atol=1e-12, max_iter=100, damped=False)[source]#

Bases: RootFinder

Find a root of a vector-valued function with Newton’s method.

Check https://en.wikipedia.org/wiki/Newton%27s_method_in_optimization for details.

Parameters:

  • atol (float) – Tolerance to check algorithm convergence.

  • max_iter (int) – Maximum iterations.

  • damped (bool) – Whether to use a damped version. Check p358 of [N2018].

References

[N2018]

Yurii Nesterov. Lectures on Convex Optimization. Springer, 2018.

root(fun, x0, fun_jac)[source]#

Find a root of a vector-valued function.

Parameters:

  • fun (callable) – Vector-valued function.

  • x0 (array-like) – Initial guess.

  • fun_jac (callable) – Jacobian of fun.

class geomstats.numerics.optimization.RootFinder[source]#

Bases: ABC

Find a root of a vector-valued function.

abstract root(fun, x0, fun_jac=None)[source]#

Find a root of a vector-valued function.

Parameters:

  • fun (callable) – Vector-valued function.

  • x0 (array-like) – Initial guess.

  • fun_jac (callable) – Ignored if None.

Returns:

res (OptimizeResult)

class geomstats.numerics.optimization.ScipyMinimize(method='L-BFGS-B', autodiff_jac=False, autodiff_hess=False, bounds=None, constraints=(), tol=None, callback=None, options=None, save_result=False)[source]#

Bases: object

Wrapper for scipy.optimize.minimize.

Only autodiff_jac and autodiff_hess differ from scipy: if True, then automatic differentiation is used to compute jacobian and/or hessian.

minimize(fun, x0, fun_jac=None, fun_hess=None, hessp=None)[source]#

Minimize objective function.

Parameters:

  • fun (callable) – The objective function to be minimized.

  • x0 (array-like) – Initial guess.

  • fun_jac (callable) – Jacobian of fun.

  • fun_hess (callable) – Hessian of fun.

  • hessp (callable)

class geomstats.numerics.optimization.ScipyRoot(method='hybr', autodiff_jac=False, tol=None, callback=None, options=None, save_result=False)[source]#

Bases: RootFinder

Wrapper for scipy.optimize.root.

Only autodiff_jac differs from scipy: if True, then automatic differentiation is used to compute jacobian.

root(fun, x0, fun_jac=None)[source]#

Find a root of a vector-valued function.

Parameters:

  • fun (callable) – Vector-valued function.

  • x0 (array-like) – Initial guess.

  • fun_jac (callable) – Jacobian of fun. Ignored if None.

Returns:

res (OptimizeResult)

geomstats.numerics.path module#

Discrete-path related machinery.

class geomstats.numerics.path.UniformlySampledDiscretePath(path, interpolator=None, **interpolator_kwargs)[source]#

Bases: object

A uniformly-sampled discrete path.

Parameters:

  • path (array-like, […, *point_shape])

  • interpolator (Interpolator1D)

class geomstats.numerics.path.UniformlySampledPathEnergy(space)[source]#

Bases: object

Riemannian path energy of a uniformly-sampled path.

Parameters:

space (Manifold) – Equipped manifold.

energy(path)[source]#

Compute Riemannian path energy.

Parameters:

path (array-like, shape=[…, n_times, *point_shape]) – Piecewise linear path.

Returns:

energy (array-like, shape=[…,]) – Path energy.

energy_per_time(path)[source]#

Compute Riemannian path enery per time.

Parameters:

path (array-like, shape=[…, n_times, *point_shape]) – Piecewise linear path.

Returns:

energy (array-like, shape=[…, n_times - 1,]) – Stepwise path energy.

Module contents#