MFOPT - Optimization over Manifolds

Download the most current release: Version 1.0, 15/10/2013

This package contains the MFOPT MATLAB library that was used for the computations in the original ICCV 2013 paper.

Features: MFOPT can be used to approximately solve variational problems with total variation-base regularization and arbitrary data terms for functions that are constrained to values in a manifold. The library includes implementations for

  • Euclidean spaces R^n,
  • n-dimensional unit spheres S^n,
  • the Moebius band,
  • the three-dimensional rotation group SO(3).

Other manifolds can be easily added due to the modular structure.

Requirements: 64-bit version of MATLAB 2013a running under 64-bit Windows 7 or later.


  • Figures 3,4,6 in the original paper were generated using a different code base and are not included.
  • To clearly see the discretization effects as in Fig. 2, the problems have to be solved to a high accuracy. Low accuracy solutions that are not fully converged can actually look better (although they are in fact incorrect solutions). For the paper, a commercial solver was used to obtain high accuracy. The package contains a version that uses the simpler primal-dual method but gives much lower accuracy.

License: You may freely use and modify the source code for non-commercial research, as long as a reference is included in related publications. Please cite:

J. Lellmann, E. Strekalovskiy, S. Koetter, and D. Cremers: Total Variation Regularization for Functions with Values in a Manifold. In: International Conference of Computer Vision (ICCV), 2013. In press.

Quick Start

  • Download and unzip to a folder of your choice
  • run any of the paper_* scripts from within MATLAB to solve the optimization problems, or the vis_* scripts to visualize the results and generate the figures. A good starting point is paper_moebius.m.
  • Stopping criteria are set to high accuracy by default and will take a long time to compute. Computation time can be reduced by lowering the maximum iteration count term_maxiter or increasing the target accuracy term_relgap.


  • The paper_* scripts read the input data from the input folder, and solve the corresponding optimization problems.
  • All output is stored in separate folders within the results folder, named by the date and time the script was run.
  • The results_important folder contains copies of the exact results shown in the paper.
  • The vis_* scripts use the data in the results_important folder to generate the figures. The figures are stored in the results folder.
  • The actual minimization happens in solve_manifold(). If use_gpu = true (default), the actual iterations are carried out by the MEX function mf_iterate. Despite the field name, the open implementation is currently CPU only, so no GPU is required.

Adding your own manifold

All information about the manifold is passed to the solver in form of a struct that is created by one of the manifold_* functions. In order to add your own manifold, the best way is to modify one of the manifold_* functions.

The manifold structure has the following fields:

  • m.mdims - struct holding the manifold dimensions and various constants:
    • m.mdims.dmanifold - dimension of the manifold: n for R^n, m for the m-dimensional sphere S^m, 2 for the Moebius band, and so on.
    • m.mdims.gradients - number of “evaluation points” at which the gradient is evaluated. Denoted by m in the paper.
    • m.mdims.points - size of the neighbourhood N_j in (10). Each gradient computation in one of the evaluation points uses the value of u at m.mdims.points points. This is used together with m.P (see below).
    • m.mdims.vecs - the number of vectors (labels) used to discretize the manifold. Denoted by l in the paper.
  • m.v - the vectors/points on the manifold, embedded into some R^n, that correspond to the labels. Denoted by z^k in the paper.
  • m.b - the weights b used to integrate functions over the manifold as in (9).
  • m.P - defines the neighbourhood N_j for the evaluation points as in (10). For the j-th evaluation point y^j, the neighbourhood consists of all x^{i+1} with i in P(:,j). P must be zero-based and of int64 type. mf_iterate contains a hard-coded limit which requires size(P,1) ⇐ 25 if use_gpu = true is used.
  • m.A,m.B - These are the matrices A and B from equation (14) in the paper. They define the way gradients are computed on the manifold.
  • m.dist(v1,v2) - function handle, computes the geodesic distance between a pair of points x,y on the manifold. This uses the same vector representation/embedding into R^n as m.v.
  • m.mean(z,u) - function handle, computes the Fréchet mean (22) for the vectors z and weights u.
  • m.plot() - function handle, sets up a (2D or 3D) figure and axis and plots a view of the manifold into which individual points can be plotted later.
  • m.embed(v) - function handle, for given points v on the manifold returns a points in the same (2D or 3D) space as m.plot that are used to visualize the points v.

The easiest way to find out the correct dimensions is to check out manifold_moebius.m or manifold_Rn.m. The solver will display an error if any of the dimensions are inconsistent with the m.mdims struct.

If you have any questions or need help, feel free to send me an e-mail at Jan.Lellmann (at)


  • The computation of the area on the sphere uses parts of the SPHERE_DELAUNAY and SPHERE_CVT libraries written by John Burkardt.