MADMM: a generic algorithm for non-smooth optimization on manifolds

Introduction

A wide range of problems in machine learning, pattern recognition, computer vision, and signal processing is formulated as optimization problems where the variables are constrained to lie on some Riemannian manifold. For example, optimization on the Grassman manifold comes up in multi-view clustering and matrix completion . Optimization on the Stiefel manifold arises in a plethora of applications ranging from classical ones such as eigenvalue problems, assignment problems, and Procrustes problems , to more recent ones such as 1-bit compressed sensing . Problems involving products of Stiefel manifolds include shape correspondence , manifold learning , sensor localization , structural biology , and structure from motion recovery . Optimization on the sphere is used in principle geodesic analysis , a generalization of the classical PCA to non-Euclidean domains. Optimization over the manifold of fixed-rank matrices arises in maxcut problems , sparse principal component analysis , regression , matrix completion , and image classification . Oblique manifolds are encountered in problems such as independent component analysis and joint diagonalization , blind source separation , and prediction of stock returns .

Though some instances of manifold optimization such as eigenvalues problems have been treated extensively in the distant past, the first general purpose algorithms appeared only in the 1990s . With the emergence of numerous applications during the last decade, especially in the machine learning community, there has been an increased interest in general-purpose optimization on different manifolds , leading to several manifold optimization algorithms such as conjugate gradients , trust regions , and Newton . Boumal et al. released the MATLAB package Manopt, as of today the most complete generic toolbox for smooth optimization on manifolds, including a variety of manifolds and solvers.

In this paper, we are interested in manifold-constrained minimization of non-smooth functions. Recent applications of such problems include several formulations of robust PCA , the computation of compressed modes in Euclidean and non-Euclidean spaces, robust Euclidean embedding , synchronization of rotation matrices , and functional correspondence .

Broadly speaking, optimization methods for non-smooth functions break into three classes of approaches. First, smoothing methods replace the non-differentiable objective function with its smooth approximation . Such methods typically suffer from a tradeoff between accuracy (how far is the smooth approximation from the original objective) and convergence speed (less smooth functions are usually harder to optimize). A second class of methods use subgradients as a generalization of derivatives of non-differentiable functions. In the context of manifold optimization, several subgradient approaches have been proposed . The third class of methods are splitting approaches. Such methods have been studied mostly for problems involving the minimization of matrix functions with orthogonality constraints. Lai and Osher proposed the method of splitting orthogonal constraints (SOC) based on the Bregman iteration . Neumann et al. used a different splitting scheme for the same class of problems.

Contributions

In this paper, we propose Manifold alternating direction method of multipliers (MADMM), an extension of the classical ADMM scheme for manifold-constrained non-smooth optimization problems. The core idea is a splitting into a smooth problem with manifold constraints and a non-smooth unconstrained optimization problem. Our method has a number of advantages common to ADMM approaches. First, it is very simple to grasp and implement. Second, it is generic and not limited to a specific manifold, as opposed e.g. to developed for the Stiefel manifold, or developed for the oblique manifold. Third, it makes very few assumptions about the properties of the objective function. Fourth, in some settings, our method lends itself to parallelization on distributed computational architectures . Finally, our method demonstrates faster convergence than previous methods in a broad range of applications.

Manifold optimization

The term manifold- or manifold-constrained optimization refers to a class of problems of the form

A conceptual gradient descent-like manifold optimization is presented in Algorithm 1. For a comprehensive introduction to manifold optimization, the reader is referred to .

Manifold ADMM

Let us now consider general problems of the form

where $f$ and $g$ are smooth and non-smooth real-valued functions, respectively, $A$ is a $k\times m$ matrix, and the rest of the notation is as in problem (1). Examples of $g$ often used in machine learning applications are nuclear-, $L_{1}$-, or $L_{2,1}$-norms. Because of non-smoothness of the objective function, Algorithm 1 cannot be used directly to minimize (1).

In this paper, we propose treating this class of problems using the alternating directions method of multipliers (ADMM). The key idea is that problem (2) can be equivalently formulated as

by introducing an artificial variable $Z$ and a linear constraint. The method of multipliers , applied to only the linear constraints in (3), leads to the minimization problem

Note that MADMM is extremely simple and easy to implement. The $X$-step is the setting of Algorithm 1 and can be carried out using any standard smooth manifold optimization method. Similarly to common implementation of ADMM algorithms, there is no need to solve the $X$-step problem exactly; instead, only a few iterations of manifold optimization are done. Furthermore, for some manifolds and some functions $f$, the $X$-step has a closed-form solution. The implementation of the $Z$-step depends on the non-smooth function $g$, and in many cases has a closed-form expression: for example, when $g$ is the $L_{1}$-norm, the $Z$-step boils down to simple shrinkage, and when $g$ is nuclear norm, the $Z$-step is performed by singular value shrinkage. $\rho 0$ is the only parameter of the algorithm and its choice is not critical for convergence. In our experiments, we used a rather arbitrary fixed value of $\rho$, though in the ADMM literature it is common to adapt $\rho$ at each iteration, e.g. using the strategy described in .

Results and Applications

In this section, we show experimental results providing a numerical evaluation of our approach on several challenging applications from the domains of dimensionality reduction, data analysis, pattern recognition, and manifold learning. All our experiments were implemented in MATLAB; we used the conjugate gradients and trust regions solvers from the Manopt toolbox for the $X$-step. Time measurements were carried out on a PC with Intel Xeon 2.4 GHz CPU.

Our first application is the computation of compressed modes, an approach for constructing localized Fourier-like orthonormal bases recently introduced in . Let us be given a manifold $\mathcal{S}$ with a Laplacian $\Delta$, where in this context, ‘manifold’ can refer to both continuous or discretized manifolds of any dimension, represented as graphs, triangular meshes, etc., and should not be confused with the matrix manifolds we have discussed so far referring to manifold-constrained optimization problems. Here, we assume that the manifold is sampled at $n$ points and the Laplacian is represented as an $n\times n$ sparse symmetric matrix. In many machine learning applications such as spectral clustering , non-linear dimensionality reduction, and manifold learning , one is interested in finding the first $k$ eigenvectors of the Laplacian $\Delta\Phi=\Phi\Lambda$, where $\Phi$ is the $n\times k$ matrix of the first eigenvectors arranged as columns, and $\Lambda$ is the diagonal $k\times k$ matrix of the corresponding eigenvalues.

The first $k$ eigenvectors of the Laplacian can be computed by minimizing the Dirichlet energy

with orthonormality constraints. Laplacian eigenfunctions form an orthonormal basis on the Hilbert space $L^{2}(\mathcal{S})$ with the standard inner product, and are a generalization of the Fourier basis to non-Euclidean domains. The main disadvantage of such bases is that its elements are globally supported. Ozoliņš et al. proposed a construction of localized quasi-eigenbases by solving

where $\mu>0$ is a parameter. The $L_{1}$-norm (inducing sparsity of the resulting basis) together with the Dirichlet energy (imposing smoothness of the basis functions) lead to orthogonal basis functions, referred to as compressed modes that are localized and approximately diagonalize $\Delta$. Lai and Osher and Neumann et al. proposed two different splitting methods for solving problem (6). Due to lack of space, the reader is referred to for a detailed description of both methods.

Solution

Results

To study the behavior of ADMM, we used a simple 1D problem with a Euclidean Laplacian constructed on a line graph with $n$ vertices. Figure 2 (top left) shows the convergence of MADMM with different random initializations. We observe that the method converges globally independently of the initialization. Figure 2 (top right) shows the convergence of MADMM using different solvers and number of iterations in the $X$-step. We did not observe any significant change in the behavior. Figure 2 (bottom left) studies the scalability of different algorithms, speaking clearly in favor of MADMM compared to the methods of . Figure 1 shows compressed modes computed on a triangular mesh of a human sampled at 8K vertices, using the cotangent weights formula to discretize the Laplacian. Figure 2 (bottom right) shows the convergence of different methods on this dataset. MADMM shows the best performance among the compared methods.

2 Functional correspondence

Our second problem is coupled diagonalization, which is used for finding functional correspondence between manifolds and multi-view clustering . Let us consider a collection of $L$ manifolds $\{\mathcal{S}_{i}\}_{i=1}^{L}$, each discretized at $n_{i}$ points and equipped with a Laplacian $\Delta_{i}$ represented as an $n_{i}\times n_{i}$ matrix. The functional correspondence between manifolds $\mathcal{S}_{i}$ and $\mathcal{S}_{j}$ is an $n_{j}\times n_{i}$ matrix $T_{ij}$ mapping functions from $L^{2}(\mathcal{S}_{i})$ to $L^{2}(\mathcal{S}_{j})$. It can be efficiently approximated using the first $k$ Laplacian eigenvectors as $T_{ij}\approx\Phi_{j}X_{ij}\Phi_{i}^{\top}$, where $X_{ij}$ is the $k\times k$ matrix translating Fourier coefficients from basis $\Phi_{i}$ to basis $\Phi_{j}$, represented as $n_{i}\times k$ and $n_{j}\times k$ matrices, respectively. Imposing a further assumption that $T_{ij}$ is volume-preserving, $X_{ij}$ must be an orthonormal matrix , and thus can be represented as a product of two orthonormal matrices $X_{ij}=X_{i}X_{j}^{\top}$. For each pair of manifolds $\mathcal{S}_{i},\mathcal{S}_{j}$, we assume to be given a set of $q_{ij}$ functions in $L^{2}(\mathcal{S}_{i})$ arranged as columns of an $n_{i}\times q_{ij}$ matrix $F_{ij}$ and the corresponding functions in $L^{2}(\mathcal{S}_{j})$ represented by the $n_{j}\times q_{ij}$ matrix $G_{ij}$. The correspondence between all the manifolds can be established by solving the problem

Solution

Results

We computed functional correspondences between $L=6$ human 3D shapes from the TOSCA dataset using $k=25$ basis functions and $q=25$ seeds as correspondence data, contaminated by 16% outliers. Figure 3 (left) analyzes the resulting correspondence quality using the Princeton protocol , plotting the percentage of correspondences falling within a geodesic ball of increasing radius w.r.t. the groundtruth correspondence. For comparison, we show the results of a least-squares solution used in (see Figure 4). Figure 3 (right) shows the convergence of MADMM in a correspondence problem with $L=2$ shapes. For comparison, we show the convergence of a smoothed version of the $L_{2,1}$-norm $\|A\|_{2,1}\approx\sum_{j}\left(\sum_{i}a_{ij}^{2}+\epsilon\right)^{1/2}$ in (7) for various values of the smoothing parameter $\epsilon$.

3 Robust Euclidean embedding

known as classical MDS or classical scaling, which has a closed form solution by means of eigendecomposition of $HDH$.

The main disadvantage of classical MDS is the fact that noise in a single entry of the distance matrix $D$ is spread over entire column/row by the double centering transformation. To cope with this problem, Cayton and Dasgupta proposed an $L_{1}$ version of the problem,

where the use of the $L_{1}$-norm efficiently rejects outliers. The authors proposed two solutions for problem (9): a semi-definite programming (SDP) formulation and a subgradient descent algorithm (the reader is referred to for a detailed description of both methods).

Solution

Here, we consider (9) as a non-smooth optimization of the form (2) on the manifold of fixed-rank positive semi-definite matrices and solve it using MADMM. Note that in this case, we have only the non-smooth function $g$ and $f\equiv 0$. The $X$-step of the MADMM algorithm is manifold optimization of a quadratic function, carried out using two iterations of manifold conjugate gradients solver. The $Z$-step is performed by shrinkage. In our experiments, all the compared methods were initialized with the classical MDS solution and the value $\rho=10$ was used for MADMM. SDP approach was implemented using MATLAB CVX toolbox .

Results

Figure 5 shows an example of 2D Euclidean embedding of the distances between 500 US cities, contaminated by sparse noise (the distance between two cities was doubled, as in a similar experiment in ). The robust embedding is insensitive to such outliers, while the classical MDS result is completely ruined. Figure 6 (right) shows an example of convergence of the proposed MADMM method and the subgradient descent of on the same dataset. We observed that our algorithm outperforms the subgradient method in terms of convergence speed. Furthermore, the subgradient method appears to be very sensitive to the initial step size $c$; choosing too small a step leads to slower convergence, and if the step is too large the algorithm may fail to converge. Figure 6 (left) studies the scalability of the subgradient-, SDP-, and MADMM-based solutions for the REE problem, plotting the complexity of a single iteration as function of the problem size on random data. The typical number of iterations was of the order of 20 for SDP, 50 for MADMM, and 500 for the subgradient method. We see that MADMM scaled better than other approaches, and that SDP is not applicable to large problems.

Discussion and Conclusions

We presented MADMM, a generic algorithm for optimization of non-smooth functions with manifold constraints, and showed that it can be efficiently used in many important problems from the domains of machine learning, computer vision and pattern recognition, and data analysis. Among the key advantages of our method is its remarkable simplicity and lack of parameters to tune - in all our experiments, it worked entirely out-of-the-box. We believe that MADMM will be very useful in many other applications in this community. In our experiments, we observed that MADMM converged independently of the initialization; a theoretical study of convergence properties is an important future direction.

Acknowledgement

The work is supported by the ERC Starting grant No. 307047.

References