Define your own manifold

Define your own manifold#

Although CDOpt provides various pre-defined manifolds in cdopt.manifold_np and cdopt.manifold_torch, in many scenarios the users may only know the expression of the constraints, and it is difficult to determine all the essential materials for Riemannian optimization. Existing Riemannian optimization packages cannot deal with that issue, since they require geometrical materials (e.g., retractions and their inverse, vector transports, exponential mappings, logarithmic mappings, …) to describe the Riemannian manifold.

CDOpt provides an easy and direct way to define customized manifolds. For generalized constraints \(c(x) = 0\), the corresponding manifold can be directly constructed by the cdopt.manifold.basic_manifold structure, or its dependent class such as cdopt.manifold.basic_manifold_np and cdopt.manifold.basic_manifold_torch. one only needs to provide the expression of \(c(x)\), then all the other essential materials can be automatically generated by the selected AD packages.

A basic example#

Consider the Riemannian manifold

\[ \{ X \in \mathbb{R}^{m\times s}: X^\top XT = I_s \}, \]

where \(T \in \mathbb{R}^{s\times s}\) is a positive definite matrix (i.e., all its eigenvalues are positive). Such manifold have not been defined by existing Riemannian optimization packages, as far as we known in June, 2022. In CDOpt, we can manually define it through cdopt.manifold.basic_manifold and autograd package. In this way, we only need to provide several basic functions:

  1. basic_manifold.C: The expression of \(c\).

  2. basic_manifold.v2m: Converting flattened vector to the variables in autograd.

  3. basic_manifold.m2v: Flattening the variables in autograd.

  4. cdopt.manifold.basic_manifold.Init_point: The function that generated initial points for solvers.

import numpy as np
import cdopt
from cdopt.manifold import basic_manifold
from cdopt.core.problem import problem
from cdopt.core.backbone_autograd import backbone_autograd
class my_manifold(basic_manifold):
    def __init__(self, var_shape, T):
        
        m = var_shape[0]
        s = var_shape[1]
        self.T = T
        self.Is = np.eye(s)
        
        super().__init__('custom_manifold',(m,s), (s,s),  regularize_value = 0.01, backbone = backbone_autograd)


    def C(self, X):
        return (X.T @ X)@ self.T  - self.Is
    
    def v2m(self,x):
        return np.reshape(x, self.var_shape)

    def m2v(self,X):
        return X.flatten()

    def Init_point(self, Xinit = None):
        if Xinit == None:
        	Xinit = np.random.randn(*self.var_shape)
        return Xinit

Then we can set the parameters \(n\), \(p\) and randomly generate \(T\), and initialize the Riemannian manifold in the following code.

m = 50
s = 8
T = np.random.randn(s,s)
T = T @ T.T 
T = T/np.linalg.norm(T,2) + np.eye(s)

M = my_manifold((m, s), T)

Next, we define the objective function. Notice that previously we set backbone = 'autograd' in my_manifold. As a result, we need to construct the objective function that adapts to the autograd package, and set backbone = 'autograd' in initializing the Problem structure.

import autograd.numpy as anp
def obj_fun(X):
    return -0.5 * anp.sum(X * X) 

problem_test = problem(M, obj_fun, beta = 1000)

Finally, we retrieve the gradients and hessians of CDF function from problem.test and apply scipy.optimize.lbfgs solver to minimize CDF.

cdf_fun_np = problem_test.cdf_fun_vec_np   
cdf_grad_np = problem_test.cdf_grad_vec_np
cdf_hvp_np = problem_test.cdf_hvp_vec_np


import scipy as sp
from scipy.optimize import fmin_bfgs, fmin_cg, fmin_l_bfgs_b, fmin_ncg
# Implement L-BFGS solver from scipy.optimize
Xinit = problem_test.Xinit_vec_np
out_msg = sp.optimize.minimize(cdf_fun_np, Xinit,method='L-BFGS-B',jac = cdf_grad_np)

Furthermore, one can provide more information in define the my_manifold. For example, its constraint dissolving mapping can be chosen as $\( \mathcal{A}(X) = X - \frac{1}{2}X( X^TX T - I_p ). \)$ Then we can add such information to the definition of my_manifold,

import numpy as np
import cdopt
from cdopt.manifold import basic_manifold
from cdopt.core.problem import Problem
class my_manifold_with_A(basic_manifold):
    def __init__(self, var_shape, T):
        
        m = var_shape[0]
        s = var_shape[1]
        self.T = T
        self.Is = np.eye(s)
        super().__init__('custom_manifold',(m,s), (s,s),  regularize_value = 0.01, backbone = 'autograd')


    def C(self, X):
        return (X.T @ X)@ self.T  - self.Is
    
    def A(self, X):
        return X - 0.5 * X@((X.T @ X)@ self.T  - self.Is)
    
    def v2m(self,x):
        return np.reshape(x, self.var_shape)

    def m2v(self,X):
        return X.flatten()

    def Init_point(self, Xinit = None):
        if Xinit == None:
            Xinit = np.random.randn(*self.var_shape) 
          
        return Xinit

That would significantly accelerate the optimization. Generally speaking, the more one manually provides in defining the manifold, the faster those optimization algorithms could become. If one pursues the conveniences in defining the manifolds, he/she can only provide the expression of constraints \(c(x)\). On the other hand, if one pursues the efficiency in optimization, he/she may need to provide all the essential materials in defining the manifolds.

Higher-level APIs#

CDOpt provides the following higher-level APIs for defining customized manifolds,

  • cdopt.manifold_np.basic_manifold_np: The manifold class that defines the Riemannian manifold based on numpy and autograd.

  • cdopt.manifold_torch.basic_manifold_torch: The manifold class that defines the Riemannian manifold based on pytorch.

In these higher-level APIs, defining customized manifolds only requires the expression of constraints \(c(x) = 0\). Here is an example:

import numpy as np
import cdopt
from cdopt.manifold_np import basic_manifold_np
from cdopt.core.problem import Problem
from cdopt.core.backbone_autograd import backbone_autograd
class my_manifold_np(basic_manifold_np):
    def __init__(self, var_shape, T):
        
        m = var_shape[0]
        s = var_shape[1]
        self.T = T
        self.Is = np.eye(s)
        super().__init__('custom_manifold',(m,s), (s,s),  regularize_value = 0.01)


    def C(self, X):
        return (X.T @ X)@ self.T  - self.Is

Define manifolds through PyTorch#

Defining the manifolds based on PyTorch is slightly different with those on Numpy and autograd. We strongly suggest that all the parameters should be stored in basic_manifold._parameters, which is an ordered dictionary (OrderedDict). Moreover, as PyTorch supports the numerical computations in different data types and devices, we could specify the device and dtype in the initialization of my_manifold_torch.

from cdopt.manifold_torch import basic_manifold_torch

class my_manifold_torch(basic_manifold_torch):
    def __init__(self, var_shape, T, device = torch.device('cpu'), dtype = torch.float64):
        
        m = var_shape[0]
        s = var_shape[1]

        super().__init__('custom_manifold',(m,s), (s,s),  regularize_value = 0.01,  device= device ,dtype= dtype)
        
        self._parameters["T"] = T.to(device = device, dtype = dtype)
        self._parameters["Is"] = torch.eye(s).to(device = device, dtype = dtype)
        
        self.T = self._parameters["T"]
        self.Is = self._parameters["Is"]
        
        


    def C(self, X):
        return (X.T @ X)@ self.T  - self.Is

Then we can test how it works.

m = 50
s = 8
T = np.random.randn(s,s)
T = T @ T.T 
T = T/np.linalg.norm(T,2) + np.eye(s)
T = torch.as_tensor(T)

M = my_manifold_torch((m, s), T, device = torch.device('cpu'), dtype = torch.float64) 

def obj_fun(X):
    return -0.5 * torch.sum(X * X) 

problem_test = problem(M, obj_fun, beta = 1000)


cdf_fun_np = problem_test.cdf_fun_vec_np   
cdf_grad_np = problem_test.cdf_grad_vec_np
cdf_hvp_np = problem_test.cdf_hvp_vec_np


import scipy as sp
from scipy.optimize import fmin_bfgs, fmin_cg, fmin_l_bfgs_b, fmin_ncg
# Implement L-BFGS solver from scipy.optimize
Xinit = problem_test.Xinit_vec_np
out_msg = sp.optimize.minimize(cdf_fun_np, Xinit, method='L-BFGS-B', jac = cdf_grad_np)
Contents