pymanopt

In Pymanopt, cost functions, gradients, and Hessian-vector products (hvps) must be defined as Python callables. To enable automatic differentiation, you must annotate these callables with a backend decorator.

Decorating a callable wraps it in a pymanopt.autodiff.Function instance. This provides a backend-agnostic API to the pymanopt.core.problem.Problem class, allowing Pymanopt to compute derivatives automatically.

Important: Signature Matching The signature of your decorated callable must match the point layout of the manifold it is defined on. For example, if using the pymanopt.manifolds.fixed_rank.FixedRankEmbedded manifold, points are represented via singular value decomposition rather than a single matrix. Therefore, your cost function must accept three arguments (u, s, and vt) instead of one matrix. Always check the specific manifold's documentation to determine the correct argument signature for your cost function.

# Example conceptual usage:
# @backend_decorator
# def cost_function(u, s, vt):
#     ...
#     return value

Pymanopt follows a modular four-step workflow to solve optimization problems on manifolds:

1. Instantiate a manifold: Select a manifold from the pymanopt.manifolds package (e.g., Sphere).
2. Define a cost function: Create a function $f: \mathcal{M} \to \mathbb{R}$ to minimize. Use a backend decorator from pymanopt.function (like @pymanopt.function.autograd) to enable automatic differentiation.
3. Create a Problem: Instantiate a pymanopt.Problem by tying the manifold and the cost function together.
4. Run an optimizer: Instantiate an optimizer from pymanopt.optimizers and call .run(problem).

import autograd.numpy as anp
import pymanopt

# 1. Instantiate manifold
manifold = pymanopt.manifolds.Sphere(3)

# 2. Define cost function with backend decorator
@pymanopt.function.autograd(manifold)
def cost(point):
    return -point @ matrix @ point

# 3. Create problem
problem = pymanopt.Problem(manifold, cost)

# 4. Run optimizer
optimizer = pymanopt.optimizers.SteepestDescent()
result = optimizer.run(problem)

If your manifold is a submanifold of a Euclidean space, subclass RiemannianSubmanifold. This class provides automatic implementations for converting Euclidean gradients and Hessians to their Riemannian counterparts, provided you implement the necessary geometric components.

Requirements

To use the automatic conversions, you must implement:

1. projection(point, vector): To project vectors onto the tangent space.
2. weingarten(point, tangent_vector, normal_vector): The Weingarten map (shape operator), which takes a tangent vector and a normal vector to produce a tangent vector.

Automatic Methods Provided

• euclidean_to_riemannian_gradient(point, euclidean_gradient): Uses projection to convert the gradient.
• euclidean_to_riemannian_hessian(point, euclidean_gradient, euclidean_hessian, tangent_vector): Uses both projection and weingarten to convert the Hessian.

from pymanopt.manifolds import RiemannianSubmanifold

class MySubmanifold(RiemannianSubmanifold):
    def __init__(self):
        super().__init__(name="Submanifold", dimension=2)

    # Implement abstract Manifold methods (inner_product, projection, etc.)
    
    def weingarten(self, point, tangent_vector, normal_vector):
        # Implementation of the Weingarten map
        return tangent_vector

To implement a new Riemannian manifold in Pymanopt, you must create a subclass of the Manifold base class.

Initialization

When initializing your subclass, provide:

• name: A string representation.
• dimension: The dimension of the tangent spaces.
• point_layout: An integer or sequence of integers describing how points are represented. For standard numpy arrays, use 1. For manifolds representing points as a tuple/list of $n$ arrays, use n.

Required Implementations

Every custom manifold must implement these abstract methods:

• inner_product(point, tangent_vector_a, tangent_vector_b): The Riemannian inner product.
• projection(point, vector): Projects an ambient space vector onto the tangent space.
• norm(point, tangent_vector): Computes the norm of a tangent vector.
• random_point(): Returns a random point on the manifold.
• random_tangent_vector(point): Returns a random tangent vector at a point.
• zero_vector(point): Returns the origin of the tangent space.

Optimizer-Specific Requirements

Depending on the optimizer you intend to use, you may need to implement additional methods:

• First-order optimizers (e.g., SteepestDescent, ConjugateGradient): Require euclidean_to_riemannian_gradient.
• Second-order optimizers (e.g., TrustRegions): Require euclidean_to_riemannian_hessian.

from pymanopt.manifolds import Manifold
import numpy as np

class MyCustomManifold(Manifold):
    def __init__(self):
        super().__init__(name="MyManifold", dimension=3, point_layout=1)

    def inner_product(self, point, tangent_vector_a, tangent_vector_b):
        return np.dot(tangent_vector_a, tangent_vector_b)

    def projection(self, point, vector):
        # Implementation of projection
        return vector

    def norm(self, point, tangent_vector):
        return np.linalg.norm(tangent_vector)

    def random_point(self):
        return np.array([0.0, 0.0, 0.0])

    def random_tangent_vector(self, point):
        return np.array([1.0, 0.0, 0.0])

    def zero_vector(self, point):
        return np.array([0.0, 0.0, 0.0])

In MoG models, a common issue is a Gaussian component collapsing onto a single data point, causing singular covariance matrices. You can mitigate this by implementing a custom line search rule that detects np.linalg.LinAlgError during the retraction step and resets the problematic component to a random point on the manifold.

To implement this, create a class that implements a search method (similar to BackTrackingLineSearcher) and a helper method (e.g., _newxnewf) that wraps the manifold.retraction call in a try-except block.

class LineSearchMoG:
    def __init__(self, contraction_factor=0.5, optimism=2, sufficient_decrease=1e-4, max_iterations=25, initial_step_size=1):
        # ... initialization ...
        pass

    def search(self, objective, manifold, x, d, f0, df0):
        # ... backtracking logic ...
        pass

    def _newxnewf(self, x, d, objective, manifold):
        newx = manifold.retraction(x, d)
        try:
            newf = objective(newx)
        except np.linalg.LinAlgError:
            # Reset component to a random point if matrix is singular
            replace = np.asarray([
                np.linalg.matrix_rank(newx[0][k, :, :]) != newx[0][0, :, :].shape[0]
                for k in range(newx[0].shape[0])
            ])
            x[0][replace, :, :] = manifold.random_point()[0][replace, :, :]
            return x, objective(x), True
        return newx, newf, False