---
abstract: 'Non-convex optimization is ubiquitous in machine learning. The Majorization-Minimization (MM) procedure systematically optimizes non-convex functions through an iterative construction and optimization of upper bounds on the objective function. The bound at each iteration is required to *touch* the objective function at the optimizer of the previous bound. We show that this touching constraint is unnecessary and overly restrictive. We generalize MM by relaxing this constraint, and propose a new optimization framework, named Generalized Majorization-Minimization ([G-MM]{}) that is more flexible compared to MM. For instance, it can incorporate application-specific biases into the optimization procedure without changing the objective function. We derive [G-MM]{}algorithms for several latent variable models and show empirically that they consistently outperform their MM counterparts in optimizing non-convex objectives. In particular, [G-MM]{}algorithms appear to be less sensitive to initialization.'
author:
- |
    Sobhan Naderi Parizi sobhan@alumni.brown.edu\
    School of Engineering, Brown University\
    Providence, RI 02912, USA Kun He hekun@bu.edu\
    Dept. of Computer Science, Boston University\
    Boston, MA 02215, USA Stan Sclaroff sclaroff@bu.edu\
    Dept. of Computer Science, Boston University\
    Boston, MA 02215, USA Pedro Felzenszwalb pff@brown.edu\
    School of Engineering, Brown University\
    Providence, RI 02912, USA
bibliography:
- 'gmm\_bibfile.bib'
title: 'Generalized Majorization-Minimization'
---

majorization-minimization, non-convex optimization, latent variable models, expectation maximization

Introduction
============

Non-convex optimization is ubiquitous in machine learning. For example, data clustering [@kmeans; @kmeans:plusplus], training classifiers with latent variables [@yu09; @dpm; @Pirsiavash14; @glvm],and training visual object detectors from weakly labeled data [@Hyun14; @Rastegari15; @Ries15] all lead to non-convex optimization problems.

Majorization-Minimization (MM) [@mm] is an optimization framework for designing well-behaved optimization algorithms for non-convex functions. MM algorithms work by iteratively optimizing a sequence of easy-to-optimize surrogate functions that bound the objective. Two of the most successful instances of MM algorithms are Expectation-Maximization (EM) [@em] and the Concave-Convex Procedure (CCCP) [@cccp]. However, both have a number of drawbacks in practice, such as sensitivity to initialization and lack of uncertainty modeling for latent variables. This has been noted in works such as [@gen_em; @dpm; @naderi12; @kumar12; @ping14].

We propose a new procedure, *Generalized Majorization-Minimization ([G-MM]{})*, for optimizing non-convex objective functions. Our approach is inspired by MM, but we generalize the bound construction process. Specifically, we relax the strict bound construction method in MM to allow for a *set* of valid bounds to be used, while still maintaining algorithmic convergence. This relaxation gives us more freedom in bound selection and can be used to design better optimization algorithms.

In training latent variable models and in clustering problems, MM algorithms such as CCCP and [$k$-means]{}are known to be sensitive to the initial values of the latent variables or cluster memberships. We refer to this problem as *stickiness* of the algorithm to the initial latent values. Our experimental results show that [G-MM]{}leads to methods that tend to be less sticky to initialization. We demonstrate the benefit of using [G-MM]{}on multiple problems, including [$k$-means]{}clustering and applications of [Latent Structural SVM]{}[s]{} to image classification with latent variables.

Related Work
------------

Perhaps the most famous iterative algorithm for non-convex optimization in machine learning and statistics is the EM algorithm [@em]. EM is best understood in the context of maximum likelihood estimation in the presence of missing data, or *latent variables*. EM is a bound optimization algorithm: in each E-step, a lower bound on the likelihood is constructed, and the M-step maximizes this bound.

Countless efforts have been made to extend the EM algorithm since its introduction. In [@gen_em] it is shown that, while both steps in EM involve optimizing some functions, it is not necessary to fully optimize the functions in each step; in fact, each step only needs to [“make progress”]{}. This relaxation can potentially avoid sharp local minima and even speed up convergence.

In another attempt, @mm proposed the Majorization-Minimization (MM) framework. MM generalizes methods like EM by [“transferring”]{} the optimization to a sequence of surrogate functions (bounds) on the original objective function. The Concave-Convex Procedure (CCCP) [@cccp] is another widely-used instance of MM, where the surrogate function is obtained by *linearizing* the concave part of the objective function. Many successful machine learning algorithms employ CCCP, the Latent SVM [@dpm]. Other instances of MM algorithms include iterative scaling [@pietra97], and non-negative matrix factorization [@nmf]. Another related line of research concerns the Difference-of-Convex (DC) programming [@dc_prog], which can be shown to reduce to CCCP under certain conditions. Convergence properties of such general “bound optimization” algorithms have been discussed by @salakhutdinov02.

Despite widespread success, MM (and CCCP in particular) has a number of drawbacks, some of which have motivated our work. In practice, CCCP often exhibits *stickiness* to initialization, which necessitates expensive initialization or multiple trials [@naderi12; @Hyun14; @multifold]. In optimizing latent variable models, CCCP lacks the ability to incorporate application-specific information without making modifications to the objective function, such as prior knowledge or side information [@xing02; @yu12], latent variable uncertainty [@kumar12; @ping14], and posterior regularization [@ganchev10]. Our framework can deal with these drawbacks. Our key observation is that we can relax the constraint enforced by MM that requires the bounds to touch the objective function, and this relaxation gives us the ability to better avoid sensitivity to initialization, and to incorporate side information.

A closely related work to ours is the pseudo-bound optimization framework by @tang14, which generalizes CCCP by using [“pseudo-bounds”]{} that may intersect the objective function. In contrast, our framework still uses valid bounds, but only relaxes the touching requirement. Also, the pseudo-bound optimization framework is not as general as ours in that it is designed specifically for optimizing binary energies in MRFs, and it restricts the form of surrogate functions to parametric max-flow.

The generalized variants of EM proposed and analyzed by @gen_em and @gam are related to our work when we restrict our attention to probabilistic models and the EM algorithm. EM can be viewed as a bound optimization procedure where the likelihood function involves both the model parameters $\theta$ and a distribution $q$ over the latent variables, denoted by $F(\theta,q)$. Choosing $q$ to be the posterior leads to a lower bound on $F$ that is tight at the current estimate of $\theta$. Generalized versions of EM, such as those given by @gen_em, use distributions other than the posterior in an alternating optimization of $F$. This fits into our framework, as we use the exact same objective function, ans only change the bound construction step (which amounts picking the distribution $q$ in EM). We propose both *stochastic* and *deterministic* strategies for bound construction, and demonstrate that they lead to higher quality solutions and less sensitivity to initialization than other EM-like methods.

Proposed Optimization Framework
===============================

.5

------------------------------------------------------------------------

height 1pt

.5

.5

------------------------------------------------------------------------

height 0.5pt

.5

\[alg:our\_framework\]

$v_0 := F(w_0)$ select ${b}_t \in \mathcal{B}_t = \mathcal{B}(w_{t-1}, v_{t-1})$ as in (\[eq:valid\_bounds\_eta\]) \[alg:our\_framework:construct\_bound\] $w_t := \operatorname*{argmin}_w {b}_t(w)$ \[alg:our\_framework:optimize\_bound\] $d_t := {b}_t(w_t) - F(w_t)$ $v_t := {b}_t(w_t) - \eta d_t$ \[alg:our\_framework:vt\] **if** $d_t < \epsilon$ **break**

.5

------------------------------------------------------------------------

height 0.5pt

.5

![image](images/CCCP_vs_ours){width="55.00000%"}

We consider minimization of functions that are bounded from below. The extension to maximization is trivial. Let $F(w) : \mathbb{R} ^d\rightarrow \mathbb{R}$ be a lower-bounded function that we wish to minimize. We propose an iterative procedure that generates a sequence of solutions $w_0,w_1,\ldots$ until it converges. The solution at iteration $t\geq 1$ is obtained by minimizing an upper bound ${b}_t(w)$ to the objective function $w_t = \operatorname*{argmin}_w {b}_t(w)$. The bound at iteration $t$ is chosen from a set of [“valid”]{} bounds $\mathcal{B}_t$ (see Figure \[fig:cccp\_vs\_ours\]). In practice, we assume that the members of $\mathcal{B}_t$ come from a family $\mathcal{F}$ of functions that upper bound $F$ and can be optimized efficiently, such as quadratic functions, or quadratic functions with linear constraints. Algorithm \[alg:our\_framework\] gives the outline of the approach.

This general scheme is used in both MM and [G-MM]{}. However, as we shall see in the rest of this section, MM and [G-MM]{}have key differences in the way they measure progress and the way they construct new bounds.

Progress Measure
----------------

MM measures progress with respect to the objective values. To guarantee progress over time MM requires that the bound at iteration $t$ must touch the objective function at the previous solution, leading to the following constraint: $$\begin{aligned}
    \textbf{MM constraint:}
    \hspace{1cm}
    b_t(w_{t-1}) = F(w_{t-1}). 
\label{eq:MMconstraint}\end{aligned}$$ This touching constraint, together with the fact that $w_t$ minimizes $b_t$ leads to $F(w_t) \le F(w_{t-1})$. That is, the value of the objective function is non-increasing over time.

The touching requirement is restrictive and, in practice, can make it hard to avoid local minima. In particular, this can make MM algorithms such as CCCP sensitive to initialization [@naderi12; @multifold; @Hyun14]. The touching constraint also eliminates the possibility of using bounds that do not touch the objective function but may have other desirable properties.

In [G-MM]{}, we measure progress with respect to the bound values. It allows us to relax the touching constraint of MM, stated in (\[eq:MMconstraint\]), and require instead that, $$\begin{aligned}
\textbf{{G-MM\xspace}constraints:}
\hspace{0.3cm}
\begin{cases}
    b_1(w_0) = F(w_0) \\
    b_t(w_{t-1}) \le b_{t-1}(w_{t-1}).
\label{eq:GMMconstraint}
\end{cases}\end{aligned}$$ Note that the [G-MM]{}constraints are weaker than MM. In fact, (\[eq:MMconstraint\]) implies (\[eq:GMMconstraint\]) because $b_{t-1}$ is an upper bound on $F$.

The weaker constraints used in [G-MM]{}do not imply $F(w_t) \le
F(w_{t-1})$, but, are sufficient to ensure that $\forall t, F(w_t) \le F(w_0)$. This follows from the fact that $b_t$ is an upper bound of $F$, $w_t$ is a minimizer of $b_t(w)$, and (\[eq:GMMconstraint\]): $$\begin{aligned}
    F(w_t) & \le b_t(w_t) \le b_t(w_{t-1}) \leq \dots \notag \\
    & \le b_1(w_1) \le b_1(w_0) = F(w_0).
    \label{eq:bound_sequence_ineq}\end{aligned}$$ $$\begin{aligned}
    F(w_t) \le b_t(w_t) \le b_t(w_{t-1}) \leq \dots \le b_1(w_1) \le b_1(w_0) = F(w_0).
    \label{eq:bound_sequence_ineq}\end{aligned}$$

Bound Construction {#sec:bound_construction}
------------------

This section describes step \[alg:our\_framework:construct\_bound\] of Algorithm \[alg:our\_framework\]. To construct new bounds, [G-MM]{} considers a [“valid”]{} subset of the upper bounds in $\mathcal{F}$ that satisfy (\[eq:GMMconstraint\]). We denote the set of valid bounds at iteration $t$ by $\mathcal{B}_t$. We consider two scenarios for constructing bounds in [G-MM]{}. The first scenario is when we have a bias function $g:\mathcal{B}_t \times \mathbb{R}^d \rightarrow \mathbb{R}$ over the set of valid bounds. The function $g$ takes in a bound $b \in \mathcal{B}_t$ and a current solution $w \in \mathbb{R}^d$ and returns a scalar indicating the goodness of the bound. In this scenario we select the bound with the largest bias value ${b}_t = \operatorname*{argmax}_{{b}\in \mathcal{B}_t} g({b}, w_{t-1})$.

In the second scenario we propose to choose one of the valid bounds from $\mathcal{B}_t$ at random. Thus, we have both a deterministic (the $1^{st}$ scenario) and a stochastic (the $2^{nd}$ scenario) bound construction mechanism. MM algorithms, such as EM and CCCP, fall into the first scenario.

In general MM algorithms are instances of [G-MM]{}that use a specific bias function $g(b,w) = -b(w)$. The bound that touches $F$ at $w$ maximizes this bias function. It also maximizes the amount of progress with respect to the previous bound at $w$. By choosing bounds that maximize progress, MM algorithms tend to rapidly converge to a nearby local minimum. For instance, at iteration $t$, the CCCP bound for latent SVMs is obtained by fixing latent variables in the concave part of the objective function to the best values according to the previous solution $w_{t-1}$, making the new solution $w_t = \operatorname*{argmin}_w {b}_t(w)$ attracted to $w_{t-1}$. Similarly, in the E-step, EM sets the posterior distribution of the latent variables conditioned on the data according to the model from the previous iteration. Thus, in the next maximization step, the model is updated to [“match”]{} the fixed posterior distributions. This explains one reason why MM algorithms are observed to be sticky to initialization.

[G-MM]{}offers a more flexible bound construction scheme than MM. In Section \[sec:experiments\] we show empirically that picking bounds uniformly at random from the set of valid bounds is less sensitive to initialization and leads to better results compared to CCCP and EM. We also show that using good bias functions can further improve performance of the learned models.

To guarantee convergence, we restrict the set of valid bounds $\mathcal{B}_t$ to those that are below a threshold $v_{t-1}$ at the previous solution $w_{t-1}$: $$\begin{aligned}
    \mathcal{B}_t &= \mathcal{B}(w_{t-1}, v_{t-1}) \notag \\
    \mathcal{B}(w, v) &= \{ b \in \mathcal{F} \mid b(w) \leq v \}.
    \label{eq:valid_bounds_eta}\end{aligned}$$

Initially, we set $v_0 = F(w_0)$ to ensure that the first bound *touches* $F$. In the next iterations, we set $v_t$ using a *progress coefficient* ${  \mbox{$\medmuskip=0mu\displaystyle\eta \in (0,1]$}}$ that ensures making, at least, $\eta d_t$ progress where ${  \mbox{$\medmuskip=0mu\displaystyled_t = b_t(w_t) - F(w_t)$}}$ is the gap between the bound and the true objective value at $w_t$.

Note that when ${  \mbox{$\medmuskip=0mu\displaystyle\eta=1$}}$ all valid bounds touch $F$ at $w_{t-1}$ (see (\[eq:valid\_bounds\_eta\])), corresponding to the MM requirement. Small $\eta$ values allow for gradual exploratory progress in each step while large $\eta$ values greedily select bounds that guarantee immediate progress.

Convergence of [G-MM]{}
=======================

We showed in (\[eq:bound\_sequence\_ineq\]) that, in any iteration $t$, the [G-MM]{}solution is no worse than the initial solution in terms of the objective value, $F(w_t) \leq F(w_0)$. In this section we also show that the gap between the [G-MM]{}bounds and the true objective converges to zero at the minimizers of the bounds (Theorem \[thrm:gmm\_bounds\_touch\_F\]). Moreover, we show, under mild conditions, [G-MM]{}converges to a solution (Theorem \[thrm:gmm\_sol\_converges\]) that is a local extremum or a saddle point of $F$ (Theorem \[thrm:GMM-converges-to-local-min\]). Other convergence properties of [G-MM]{}depend on the structure of the objective function $F$, the family of the bounds $\mathcal{F}$, and the details of the bound selection strategy.

In the limit as $t \rightarrow \infty$, [G-MM]{}bounds and $F$ **touch** at the minimizer of the bounds, $\lim_{t \rightarrow \infty} d_t = 0$, where $d_t=b_t(w_t)-F(w_t)$. \[thrm:gmm\_bounds\_touch\_F\]

We have $b_t(w_t) \leq b_t(w_{t-1}) \leq v_{t-1}$, where $v_t$ is defined as in line \[alg:our\_framework:vt\] of Algorithm \[alg:our\_framework\]. The first inequality holds because $w_t$ minimizes $b_t$ and the second inequality follows from (\[eq:valid\_bounds\_eta\]). Summing over $t$ gives $$\begin{aligned}
        \sum_{t=1}^T b_t(w_t) \leq \sum_{t=1}^T v_{t-1} = v_0 + \sum_{t=1}^{T-1} b_t(w_t) - \eta d_t\\
        \Rightarrow \; \eta \sum_{t=1}^T d_t \leq v_0 - b_T(w_T) = F(w_0) - b_T(w_T).
        \label{eq:touching_inequality}
    \end{aligned}$$ Let $w^* \in \operatorname*{argmin}_w F(w)$ be an unknown *global* minimum of $F$. Using $F(w^*) \leq F(w_T) \leq b_T(w_T)$ and (\[eq:touching\_inequality\]), we have $$\begin{aligned}
        \eta \sum_{t=1}^T d_t \leq F(w_0) - F(w^*)  
        \; \Rightarrow \; \lim_{t \rightarrow \infty} d_t = 0.
        \label{eq:dt_converges_to_zero}
    \end{aligned}$$ If a global minimizer of $F$ does not exist, we can use the assumption that $F$ is bounded from below and replace $F(w^*)$ in (\[eq:dt\_converges\_to\_zero\]) with the value of the lower bound.

If $\mathcal{F}$ is a family of $m$-strongly convex functions where $m>0$, then $\exists w^\dagger$ [*s.t. *]{}$\lim_{t \rightarrow \infty} w_t = w^\dagger$; [G-MM]{}converges to a solution. \[thrm:gmm\_sol\_converges\]

Let $f$ be an $m$-strongly convex function, and let $x, y$ be two points in the domain of $f$, also let $g \in \partial f$ be a subgradient of $f$. We have $$\begin{aligned}
        f(y) \geq f(x) + g^T (y-x) + \frac{m}{2} || x - y ||^2 .
        \label{eq:functions_with_bounded_curvature}
    \end{aligned}$$ After substituting $f = b_t$, $x = w_t$, and $y = w_{t-1}$ in (\[eq:functions\_with\_bounded\_curvature\]), and noting that the zero vector is a subgradient of $b_t$ at $w_t$, $0 \in \partial b_t(w_t)$, we get: $$\begin{aligned}
        \frac{m}{2} \left( b_t(w_{t-1}) - b_t(w_t) \right) \geq || w_t - w_{t-1} ||^2 .
        \label{eq:bounds_with_bounded_curvature}
    \end{aligned}$$ We can replace $b_t(w_{t-1})$ by $v_{t-1}$ in (\[eq:bounds\_with\_bounded\_curvature\]) due to (\[eq:valid\_bounds\_eta\]). By adding the inequalities for $T$ iterations, similar to what we did in the proof of Theorem \[thrm:gmm\_bounds\_touch\_F\], we get: $$\begin{aligned}
        \frac{m}{2} \left( F(w_0) - F(w^*) \right) \geq \sum_{t=1}^T || w_t - w_{t-1} ||^2
        \label{eq:gmm_converges}\\
        \Rightarrow \exists w^\dagger \text{ {\emph{s.t.\ }}} \lim_{t \rightarrow \infty} w_t = w^\dagger.
    \end{aligned}$$ Similar to the proof of Theorem \[thrm:gmm\_bounds\_touch\_F\], if a global minimum does not exist, we can replace $F(w^*)$ with the value of the lower bound of $F$ in (\[eq:gmm\_converges\]).

Let $F$ be continuously differentiable, and $\mathcal{F}$ be a family of strongly convex and smooth functions. Further, if $\exists \infty>M>m>0, \forall b \in \mathcal{F}, M I \succeq \nabla^2 b(w) \succeq m I$, where $I$ is the identity matrix, then $\nabla F(w^\dagger) = 0$; that is, [G-MM]{}converges to a local extremum or saddle point of $F$. \[thrm:GMM-converges-to-local-min\]

We postpone the proof of Theorem \[thrm:GMM-converges-to-local-min\] to the appendix.

Derived Optimization Algorithms {#sec:derived_algorithms}
===============================

For simplicity and ease of exposition, we primarily focus on demonstrating [G-MM]{}on *latent variable models* where bound construction naturally corresponds to [imputing]{} latent variables in the model. In this section we derive [G-MM]{}algorithms for two models widely used in machine learning, namely, [$k$-means]{}clustering and [Latent Structural SVM]{}, both belonging to the category of latent variable models. It is worth mentioning that [G-MM]{}is fully capable of handling more general non-convex problems.

[**$\boldsymbol{k}$-means**]{}Clustering
----------------------------------------

Let $\{x_1, \dots, x_n\}$ denote a set of sample points and $w=(\mu_1, \dots, 
\mu_k)$ denote a set of cluster centers. We use $z_i \in \{1, \dots, k\}$ to denote the index of the cluster assigned to the $i$-th sample point. The objective function in [$k$-means]{}clustering is defined as follows, $$\begin{aligned}
    F(w) = \sum_{i=1}^n \min_{z_i = 1}^k || x_i - \mu_{z_i} ||^2, 
    \;\; w=(\mu_1, \dots, \mu_k). 
\label{eq:kmeans_objective}\end{aligned}$$ **Bound construction:** We obtain a convex upper bound on $F$ by fixing the latent variables $(z_1, \dots, z_n)$ to certain values instead of minimizing over these variables. Bounds constructed this way are quadratic convex functions of $w=(\mu_1,\ldots,\mu_k)$, $$\begin{aligned}
    \mathcal{F} = \left\{\sum_{i=1}^n || x_i - \mu_{z_i} ||^2 \;\middle|\; \forall i,\; z_i \in \{1,\ldots, k\}^n \right\}.
    \label{eq:boundFamily_kmeans}\end{aligned}$$ The [$k$-means]{}algorithm is an instance of MM methods. The algorithm repeatedly assigns each example to its nearest center to construct a bound, and then updates the centers by optimizing the bound. We can set $g({b}, w) = -{b}(w)$ in [G-MM]{}to obtain the [$k$-means]{}algorithm. We can also define $g$ differently to obtain a [G-MM]{}algorithm that exhibits other desired properties. For instance, a common issue in clustering is cluster starvation. One can design a bias function that encourages balanced clusters by selecting $g$ appropriately. We can select a random bound from $\mathcal{B}_t$ by sampling a latent configuration ${  \mbox{$\medmuskip=0mu\displaystylez=(z_1, \dots, z_n)$}}$ uniformly from the set of configurations leading to valid bounds. Specifically, we start from a valid initial configuration ([$k$-means]{}solution) and do a random walk on a graph whose nodes are latent configurations defining valid bounds. The neighbors of a latent configuration $z$ are the latent configurations that can be obtained by changing the value of one of the $n$ latent variables in $z$.

**Bound optimization:** Optimization of a bound ${b}\in 
\mathcal{F}$ can be done in closed form by setting $\mu_j$ to be the mean of all examples assigned to cluster $j$: $$\begin{aligned}
    \mu_j = \frac{\sum_{i \in I_j} x_i}{|I_j|}, \quad I_j =\{1 \leq i \leq n \mid z_i = j\}.
    \label{eq:kmeans_update}\end{aligned}$$

[Latent Structural SVM]{}
-------------------------

A [Latent Structural SVM]{}([LS-SVM]{}) [@yu09] defines a structured output classifier with latent variables. It extends the Structural SVM [@joachims09] by introducing latent variables.

Let $\{ (x_1, y_1), \dots, (x_n, y_n) \}$ denote a set of labeled examples with $x_i \in \mathcal{X}$ and $y_i \in \mathcal{Y}$. We assume that each example $x_i$ has an associated latent value $z_i \in \mathcal{Z}$. Let $\phi(x, y, z):\mathcal{X} \times 
\mathcal{Y} \times \mathcal{Z} \rightarrow \mathbb{R}^d$ denote a *feature map*. A vector $w \in \mathbb{R}^d$ defines a classifier $\hat{y}:\mathcal{X} \rightarrow \mathcal{Y}$, $$\begin{aligned}
  \hat{y}(x) = \operatorname*{argmax}_y ( \max_z w \cdot \phi(x,y,z) ).\end{aligned}$$ The LS-SVM training objective is defined as follows, $$\begin{aligned}
    F(w) = \frac{\lambda}{2} ||w||^2 + \frac{1}{n}\sum_{i=1}^n \hspace{-0.05cm} \Big( \hspace{-0.05cm} & \max_{y, z} \left( w \cdot \phi(x_i, y, z) + \Delta(y, y_i) \right) - \max_z w \cdot \phi(x_i, y_i, z) \Big),
    \label{eq:LSSVM_objective}\end{aligned}$$ where $\lambda$ is a hyper-parameter that controls regularization and $\Delta(y, y_i)$ is a non-negative loss function that penalizes the prediction $y$ when the ground truth label is $y_i$.

**Bound construction:** As in the case of [$k$-means]{}, a convex upper bound on the LS-SVM objective can be obtained by imputing latent variables. Specifically, for each example $x_i$, we fix $z_i \in \mathcal{Z}$, and replace the maximization in the last term of the objective with a linear function $w \cdot \phi(x_i, y_i, z_i)$. This forms a family of convex piecewise quadratic bounds, $$\begin{aligned}
    \hspace{-0.23cm}\mathcal{F} \hspace{-0.1cm} = \hspace{-0.1cm} \Bigg\{ \hspace{-0.1cm} \frac{\lambda}{2} ||w||^2 + \frac{1}{n} \hspace{-0.05cm} \sum_{i=1}^n & \max_{y, z}\left( w \cdot \phi(x_i, y, z) + \Delta(y, y_i) \right) 
    - w \cdot \phi(x_i, y_i, z_i) \Bigg| \forall i, z_i \in \mathcal{Z} \hspace{-0.05cm} \Bigg\}.
    \label{eq:boundFamily_lssvm}\end{aligned}$$ The CCCP algorithm for [LS-SVM]{}selects the bound $b_t$ defined by $z^t_i=\operatorname*{argmax}_{z_i} w_{t-1} \cdot \phi(x_i, y_i, z_i)$. This particular choice is a special case of [G-MM]{}when $g({b}, w) = -{b}(w)$.

To generate random bounds from $\mathcal{B}_t$ we use the same approach as in the case of [$k$-means]{}clustering. We perform a random walk in a graph where the nodes are latent configurations leading to valid bounds, and the edges connect latent configurations that differ in a single latent variable.

**Bound optimization:** Optimization of a bound ${b}\in 
\mathcal{F}$ corresponds to a convex quadratic program and can be solved using different techniques, including gradient based methods (SGD) and the cutting-plane method [@joachims09]. We use the cutting-plane method in our experiments.

### Bias Function for Multi-fold MIL {#sec:bias_for_mfmil}

The multi-fold MIL algorithm [@multifold] was introduced for training latent SVMs for weakly supervised object localization, to deal with stickiness issues in training with CCCP. It modifies how latent variables are updated during training. @multifold divide the training set into $K$ folds, and updates the latent variables in each fold using a model trained on the other $K-1$ folds. This algorithm does not have a formal convergence guarantee. By defining a suitable bias function, we can derive a [G-MM]{} algorithm that mimics the behavior of multi-fold MIL, and yet, is convergent.

Consider training an [LS-SVM]{}. Let $S=(1,\ldots,n)$ be an ordered sequence of training sample indices. Also, let ${  \mbox{$\medmuskip=0mu\displaystylez_i \in \mathcal{Z}$}}$ denote the latent variable associated to training example $(x_i, y_i)$. Let $I$ denote a subsequence of $S$. Also, let $z_{I}^t$ denote the fixed latent variable values of training examples in $I$ in iteration $t$, and let $w(I, z_{I}^t)$ be the model trained on $\{(x_i,y_i) | i\in I\}$ with latent variables fixed to $z_{I}^t$ in the last maximization of (\[eq:LSSVM\_objective\]).

We assume access to a loss function $\ell(w,x,y,z)$. For example, for the binary latent SVM where $y\in\{-1,1\}$, $\ell$ is the hinge loss: $\ell(w,x,y,z) = \max\{0,1-y\, w\cdot\phi(x,z)\}$.

We start by considering the Leave-One-Out (LOO) setting, $K=n$, and call the algorithm of @multifold LOO-MIL in this case. The update rule of LOO-MIL in iteration $t$ is to set $$\begin{aligned}
z_i^t=\operatorname*{argmin}_{z\in\mathcal{Z}}~\ell\left(w(S\backslash i,z_{S\backslash i}^{t-1}),x_i,y_i,z\right), ~~\forall i\in S.
\label{eq:update_rule}\end{aligned}$$ After updating the latent values for all training examples, the model $w$ is retrained by optimizing the resulting bound.

Now let us derive the bias function for a [G-MM]{}algorithm that mimics the behavior of LOO-MIL. Recall from Equation \[eq:boundFamily\_lssvm\] that each bound $b \in \mathcal{B}_t$ is associated with a joint latent configuration $z(b)=(z_1,\ldots,z_n)$. We use the following bias function: $$\begin{aligned}
g(b,w) = -\sum_{i\in S} \ell\left(w(S\backslash i,z_{S\backslash i}^{t-1}),x_i,y_i,z_i\right).
\label{eq:loomil}\end{aligned}$$ Note that picking a bound according to (\[eq:loomil\]) is equivalent to the LOO-MIL update rule of (\[eq:update\_rule\]) except that in (\[eq:loomil\]) only *valid* bounds are considered; that is bounds that make at least $\eta$-progress.

For the general multi-fold case ($K<n$), the bias function can be derived similarly.

**Remarks:** Let us denote the data generating distribution as $D$ and regard the training set as i.i.d. samples from $D$, $\{(x_i,y_i)\}_{i=1}^n\overset{iid}{\sim} D$. Then, $-g(b,w)$ is an unbiased estimate of an expected generalization risk when the training set size is $n-1$, defined with respect to loss $\ell$: $$\begin{aligned}
&\mathbb{E}_{S\sim D^{n-1}}\left[R_D^{\ell}\left(w(S,z_S^{t-1})\right)\right]\label{eq:gen_risk}\\
\text{where} ~~~ & R_D^{\ell}(w)=\mathbb{E}_{(x,y)\sim D}\left[\min_{z}\ell(w,x,y,z)\right]\end{aligned}$$ This gives theoretical support for using the bias function defined in (\[eq:loomil\]). However, leave-one-out estimates usually have high variance and are time-consuming to compute, and in practice we found that using fewer folds ($K=10$) offers a better trade-off.

Experiments {#sec:experiments}
===========

We evaluate [G-MM]{}and MM algorithms on [$k$-means]{}clustering and [LS-SVM]{}training on various datasets. Recall from (\[eq:valid\_bounds\_eta\]) that the progress coefficient $\eta$ defines the set of valid bounds $\mathcal{B}_t$ in each step. CCCP and standard [$k$-means]{}bounds correspond to setting ${  \mbox{$\medmuskip=0mu\displaystyle\eta = 1$}}$, thus taking maximally large steps towards a local minimum of the true objective.

[**$\boldsymbol{k}$-means**]{}Clustering {#sec:experiments:clustering}
----------------------------------------

We conduct experiments on four different clustering datasets: Norm-25 [@kmeans:plusplus], D31 [@kmeans:D31], Cloud [@kmeans:plusplus], and GMM-200. Norm-25, D31, and GMM-200 are synthetic and Cloud is from real data. See the references for details about the datasets. GMM-200 was created by us. It is a Gaussian mixture model on $2$-D data with 200 mixture components. Each component is a Gaussian distribution with $\sigma=1.0$ and $\mu$ on a square of size $70 \times 70$. The means are placed at least $2.5 \sigma$ apart from each other. The dataset contains 50 samples per mixture component.

  -- -- ------------ --------------------------------------------------------- ------------------------------------------------- --------------------------------------------------------- ------------------------------------------------- --------------------------------------------------------- -----------------------------------------------
          **opt.**                                                                                                                                                                                                                                                                                     
         **method**                      **avg $\pm$ std**                                         **best**                                          **avg $\pm$ std**                                         **best**                                          **avg $\pm$ std**                                        **best**
          hard-EM     ${  \mbox{$\medmuskip=0mu\displaystyle1.9e5 \pm 2e5$}}$   ${  \mbox{$\medmuskip=0mu\displaystyle7.0e4$}}$   ${  \mbox{$\medmuskip=0mu\displaystyle5.8e5 \pm 3e5$}}$   ${  \mbox{$\medmuskip=0mu\displaystyle2.2e5$}}$   ${  \mbox{$\medmuskip=0mu\displaystyle5.3e3 \pm 9e3$}}$   ${  \mbox{$\medmuskip=0mu\displaystyle1.5$}}$
          [G-MM]{}    ${  \mbox{$\medmuskip=0mu\displaystyle9.7e3 \pm 1e4$}}$    ${  \mbox{$\medmuskip=0mu\displaystyle1.5$}}$     ${  \mbox{$\medmuskip=0mu\displaystyle2.0e4 \pm 0$}}$    ${  \mbox{$\medmuskip=0mu\displaystyle2.0e4$}}$   ${  \mbox{$\medmuskip=0mu\displaystyle4.5e3 \pm 8e3$}}$   ${  \mbox{$\medmuskip=0mu\displaystyle1.5$}}$
          hard-EM                         $1.69 \pm 0.03$                                           $1.21$                                           $52.61 \pm 47.06$                                          $4.00$                                            $1.55 \pm 0.17$                                          $1.10$  
          [G-MM]{}                        $1.43 \pm 0.15$                                           $1.10$                                            $1.21 \pm 0.05$                                           $1.10$                                            $1.45 \pm 0.14$                                          $1.10$
          hard-EM                         $1929 \pm 429$                                            $1293$                                           $44453 \pm 88341$                                          $3026$                                             $1237 \pm 92$                                           $1117$
          [G-MM]{}                         $1465 \pm 43$                                            $1246$                                             $1470 \pm 8$                                             $1444$                                             $1162 \pm 95$                                           $1067$
          hard-EM                         $2.25 \pm 0.10$                                           $2.07$                                           $11.20 \pm 0.63$                                           $9.77$                                            $2.12 \pm 0.07$                                          $1.99$
          [G-MM]{}                        $2.04 \pm 0.09$                                           $1.90$                                            $1.85 \pm 0.02$                                           $1.80$                                            $1.98 \pm 0.06$                                          $1.89$
  -- -- ------------ --------------------------------------------------------- ------------------------------------------------- --------------------------------------------------------- ------------------------------------------------- --------------------------------------------------------- -----------------------------------------------

We compare results from three different initializations: **forgy** selects $k$ training examples uniformly at random without replacement to define initial cluster centers, **random partition** assigns training samples to cluster centers randomly, and [**$\boldsymbol{k}$-means++**]{}uses the algorithm in [@kmeans:plusplus]. In each experiment we run the algorithm 50 times and report the mean, standard deviation, and the best objective value (Equation \[eq:kmeans\_objective\]). Table \[table:kmeans\] shows the results using [$k$-means]{} (hard-EM) and [G-MM]{}. We note that the variance of the solutions found by [G-MM]{}is typically smaller than [$k$-means]{}. Moreover, the best and the average solutions found by [G-MM]{}are always better than (or the same as) those found by [$k$-means]{}. This trend generalizes over different initialization schemes as well as different datasets.

Although *random partition* seems to be a very bad initialization for [$k$-means]{}on all datasets, [G-MM]{} recovers from it. In fact, on D31 and GMM-200 datasets, [G-MM]{}initialized by *random partition* performs better than when it is initialized by other methods (including [$k$-means++]{}). Also, the variance of the best solutions (across different initialization methods) in [G-MM]{}is smaller than that of [$k$-means]{}. These suggest that the [G-MM]{}optimization is less sticky to initialization than [$k$-means]{}.

Figure \[fig:kmeans\_visualization\] visualizes the result of [$k$-means]{}and [G-MM]{} (with random bounds) on the D-31 dataset [@kmeans:D31], from the same initialization. [G-MM]{}finds a near perfect solution, while in [$k$-means]{}many clusters get merged incorrectly or die off. Dead clusters are those which do not get any points assigned to them. The update rule of Equation  collapses the dead clusters on to the origin.

Figure \[fig:different\_progress\_coeffs\] shows the effect of the progress coefficient on the quality of the solution found by [G-MM]{}. Different colors correspond to different initialization schemes. The solid line indicates the average objective over 50 iterations, the shaded area covers one standard deviation from the average, and the dashed line indicates the best solution over the 50 trials. Smaller progress coefficients allow for more extensive exploration, and hence, smaller variance in the quality of the solutions. On the other hand, when the progress coefficient is large [G-MM]{}is more sensitive to initialization (is more sticky) and, thus, the quality of the solutions over multiple runs is more diverse. However, despite the greater diversity, the best solution is worse when the progress coefficient is large. [G-MM]{}reduces to [$k$-means]{}if we set the progress coefficient to 1 (the largest possible value).

[0.32]{} ![Visualization of clustering solutions on the D31 dataset [@kmeans:D31] from identical initializations. Random partition initialization scheme is used. (a) color-coded ground-truth clusters. (b) solution of [$k$-means]{}. (c) solution of [G-MM]{}. The white crosses indicate location of the cluster centers. Color codes match up to a permutation.[]{data-label="fig:kmeans_visualization"}](images/D31_ground-truth.pdf "fig:"){width="\textwidth"}

[0.32]{} ![Visualization of clustering solutions on the D31 dataset [@kmeans:D31] from identical initializations. Random partition initialization scheme is used. (a) color-coded ground-truth clusters. (b) solution of [$k$-means]{}. (c) solution of [G-MM]{}. The white crosses indicate location of the cluster centers. Color codes match up to a permutation.[]{data-label="fig:kmeans_visualization"}](images/D31_hard-EM.pdf "fig:"){width="\textwidth"}

[0.32]{} ![Visualization of clustering solutions on the D31 dataset [@kmeans:D31] from identical initializations. Random partition initialization scheme is used. (a) color-coded ground-truth clusters. (b) solution of [$k$-means]{}. (c) solution of [G-MM]{}. The white crosses indicate location of the cluster centers. Color codes match up to a permutation.[]{data-label="fig:kmeans_visualization"}](images/D31_generalized-MM.pdf "fig:"){width="\textwidth"}

[0.49]{} ![image](images/D31_different-progress-coeffs.png){height="2in"}

[0.49]{} ![image](images/cloud_different-progress-coeffs.png){height="2in"}

[Latent Structural SVM]{}for Image Classification and Object Detection
----------------------------------------------------------------------

We consider the problem of training an [LS-SVM]{}classifier on the mammals dataset  [@mammals]. The dataset contains images of six mammal categories with image-level annotation. Locations of the objects in these images are not provided, and therefore, treated as latent variables in the model. Specifically, let $x$ be an image and $y$ be a class label ($y \in \{1,\ldots,6\}$ in this case), and let $z$ be the latent location of the object in the image. We define $\phi(x, y, z)$ to be a feature function with $6$ blocks; one block for each category. It extracts features from location $z$ of image $x$ and places them in the $y$-th block of the output and fills the rest with zero. We use the following multi-class classification rule: $$\begin{aligned}
    y(x) = \operatorname*{argmax}_{y, z} w \cdot \phi(x, y, z), \; w = (w_1, \ldots, w_6).\end{aligned}$$ In this experiment we use a setup similar to that in [@kumar12]: we use Histogram of Oriented Gradients (HOG) for the image feature $\phi$, and the 0-1 classification loss for $\Delta$. We set $\lambda=0.4$ in Equation \[eq:LSSVM\_objective\]. We report 5-fold cross-validation performance. Three initialization strategies are considered for the latent object locations: *image center*, *top-left corner*, and *random locations*. The first is a reasonable initialization since most objects are at the center in this dataset; the second initialization strategy is somewhat adversarial.

We try a stochastic as well as a deterministic bound construction method. For the stochastic method, in each iteration $t$ we uniformly sample a subset of examples $S_t$ from the training set, and update their latent variables using $z^t_i = \operatorname*{argmax}_{z_i} 
w_{t-1} \cdot \phi(x_i, y_i, z_i)$. Other latent variables are kept the same as the previous iteration. We increase the size of $S_t$ across iterations.

For the deterministic method, we use the bias function that we described in Section \[sec:bias\_for\_mfmil\]. This is inspired by the multi-fold MIL idea [@multifold] and is shown to reduce stickiness to initialization, especially in high dimensions. We set the number of folds to ${  \mbox{$\medmuskip=0mu\displaystyleK=10$}}$ in our experiments.

  -------------------- ----------------- ---------------- ----------------- ----------------- ----------------- ---------------- -- -- -- --
                                                                                                                                          
                         **objective**    **test error**    **objective**    **test error**     **objective**    **test error**           
  **CCCP**              1.21 $\pm$ 0.03   22.9 $\pm$ 9.7   1.35 $\pm$ 0.03   42.5 $\pm$ 4.6    1.47 $\pm$ 0.03   31.8 $\pm$ 2.6           
  **[G-MM]{}random**    0.79 $\pm$ 0.03   17.5 $\pm$ 3.9   0.91 $\pm$ 0.02   31.4 $\pm$ 10.1   0.85 $\pm$ 0.03   19.6 $\pm$ 9.2           
  **[G-MM]{}biased**    0.64 $\pm$ 0.02   16.8 $\pm$ 3.2   0.70 $\pm$ 0.02   18.9 $\pm$ 5.0    0.65 $\pm$ 0.02   14.6 $\pm$ 5.4           
  -------------------- ----------------- ---------------- ----------------- ----------------- ----------------- ---------------- -- -- -- --

![image](images/topleft_cccp_5folds){width="0.32\linewidth"} ![image](images/topleft_srl_5folds){width="0.32\linewidth"} ![image](images/topleft_mfold10_5folds){width="0.32\linewidth"}

Table \[table:mammal\] shows results on the mammals dataset. Both variants of [G-MM]{}consistently outperform CCCP in terms of training objective and test error. We observed that CCCP rarely updates the latent locations, under all initializations. On the other hand, both variants of [G-MM]{}significantly alter the latent locations, thereby avoiding the local minima close to the initialization. Figure \[fig:latent\_locations\] visualizes this for *top-left* initialization. Since objects rarely occur at the top-left corner in the mammals dataset, a good model is expected to significantly update the latent locations. Averaged over five cross-validation folds, about $90\%$ of the latent variables were updated in [G-MM]{}after training whereas this measure was $2.4\%$ for CCCP. This is consistent with the better training objectives and test errors of [G-MM]{}.

We provide additional experimental results on the mammals dataset. Figure \[fig:mammals\] shows example training images and the final imputed latent object locations by three algorithms: CCCP (red), G-MM random (blue), and G-MM biased (green). The initialization is *top-left*.

In most cases CCCP fails to update the latent locations given by initialization. The two G-MM variants, however, are able to update them significantly and often localize the objects in training images correctly. This is achieved *only* with image-level object category annotations, and with a very bad (even adversarial) initialization.

![image](images/mammals/train_bison1){height="6.8em"} ![image](images/mammals/train_bison2){height="6.8em"} ![image](images/mammals/train_bison3){height="6.8em"} ![image](images/mammals/train_bison4){height="6.8em"}

![image](images/mammals/train_elephant1){height="6.8em"} ![image](images/mammals/train_elephant4){height="6.8em"} ![image](images/mammals/train_elephant6){height="6.8em"} ![image](images/mammals/train_elephant5){height="6.8em"}

![image](images/mammals/train_giraffe1){height="6.8em"} ![image](images/mammals/train_giraffe2){height="6.8em"} ![image](images/mammals/train_giraffe3){height="6.8em"} ![image](images/mammals/train_giraffe5){height="6.8em"}

![image](images/mammals/train_deer2){height="6.8em"} ![image](images/mammals/train_deer3){height="6.8em"} ![image](images/mammals/train_deer6){height="6.8em"} ![image](images/mammals/train_deer7){height="6.8em"}

![image](images/mammals/train_llama1){height="6.8em"} ![image](images/mammals/train_llama2){height="6.8em"} ![image](images/mammals/train_llama3){height="6.8em"} ![image](images/mammals/train_llama4){height="6.8em"}

![image](images/mammals/train_rhino1){height="6.8em"} ![image](images/mammals/train_rhino2){height="6.8em"} ![image](images/mammals/train_rhino3){height="6.8em"}

[Latent Structural SVM]{}for Scene Recognition {#sec:experiments:scene_classification}
----------------------------------------------

We implement the reconfigurable model of @naderi12 to do scene classification on MIT-Indoor dataset [@mit_indoor], which has images from 67 indoor scene categories. We segment each image into a ${  \mbox{$\medmuskip=0mu\displaystyle10 \times 10$}}$ regular grid and treat the grid cells as image regions. We train a model with 200 shared parts. All parts can be used to describe the data in an image region. We use the pre-trained hybrid ConvNet from @places to extract features from image regions. We record from the 4096 neurons at the penultimate layer of the network and use PCA to reduce the dimensionality of these features to 240.

The reconfigurable model is an instance of [LS-SVM]{}models. The latent variables are the assignments of parts to image regions and the output structure is the multi-valued category label predictions. [LS-SVM]{}[s]{} are known to be sensitive to initialization ([*a.k.a.*]{}the stickiness issue). @naderi12 cope with this issue by training a generative version of the model first and using it to initialize the discriminative model. Generative models are typically less sticky but perform worse in practice. To validate the hypothesis regarding stickiness of [LS-SVM]{}[s]{} we run training using several initialization strategies.

Initializing training entails the assignment of parts to image regions setting $z_i$’s in Equation \[eq:boundFamily\_lssvm\] to define the first bound. To this end we first discover 200 parts that capture discriminative features in the training data. We then run graph cut on each training image to obtain part assignments to image regions. Each cell in the ${  \mbox{$\medmuskip=0mu\displaystyle10 \times 10$}}$ image grid is a node in the graph. Two nodes in the graph are connected if their corresponding cells in the image grid are next to each other. Unary terms in the graph cut are the dot product scores between the feature vector extracted from an image region and a part filter plus the corresponding region-to-part assignment score. Pairwise terms in the graph cut implement a *Potts* model that encourages coherent labelings. Specifically, the penalty of labeling two neighboring nodes differently is $\lambda$ and it is zero otherwise. $\lambda$ controls the coherency of the initial assignments. We experiment using $\lambda \in \{ 0, 0.25, 0.5, 1 \}$. We also experiment with random initialization, which corresponds to assigning $z_i$’s randomly. This is the simplest form of initialization and does not require discovering initial part filters.

[| l “ c|c ” c|c “ c|c ” c|c " c|c |]{} & & & & &\
& **Acc.% $\pm$ std** & **O.F.** & **Acc. %** & **O.F.** & **Acc. %** & **O.F.** & **Acc. %** & **O.F.** & **Acc. %** & **O.F.**\
**CCCP** & $41.94 \pm 1.1$ & 15.20 & $40.88$ & 14.81 & $43.99$ & 14.77 & $45.60$ & 14.72 & $46.62$ & 14.70\
**[G-MM]{}random** & $47.51 \pm 0.7$ & 14.89 & $43.38$ & 14.71 & $44.41$ &14.70 & $47.12$ & 14.66 & $49.88$ & 14.58\
**[G-MM]{}biased** & $49.34 \pm 0.9$ & 14.55 & $44.83$ & 14.63 & $48.07$ & 14.51 & $53.68$ & 14.33 & $56.03$ & 14.32\

We do [G-MM]{}optimization using both random and biased bounds. For the latter we use a bias function $g({b}, w)$ that measures coherence of the labeling from which the bound was constructed. Recall from Equation \[eq:boundFamily\_lssvm\] that each bound in ${b}\in B_t$ corresponds to a labeling of the image regions. We denote the labeling corresponding to the bound ${b}$ by $z({b})=(z_1, \dots, z_n)$ where $z_i=(z_{i, 1}, \dots, z_{i, 100})$ specifies part assignments for all the 100 regions in the $i$-th image. Also, let $E(z_i)$ denote a function that measures coherence of the labeling $z_i$. In fact, $E(z_i)$ is the Potts energy function on a graph whose nodes are $z_{i, 1}, \dots, z_{i, 100}$. The graph respects a 4-connected neighborhood system (recall that $z_{i, r}$ corresponds to the $r$-th cell in the ${  \mbox{$\medmuskip=0mu\displaystyle10 \times 10$}}$ grid defined on the $i$-th image). If two neighboring nodes $z_{i, r}$ and $z_{i, s}$ get different labels the energy $E(z_i)$ increases by 1. For biased bounds we use the following bias function which favors bounds that correspond to more coherent labelings: $$g({b}, w) = - \sum_{i=1}^n E(z_i), \hspace{0.5cm} z({b}) = (z_1, \dots, z_n).
    \label{eq:bias_function:LSSVM}$$

Table \[table:lssvm\] compares performance of models trained using CCCP and [G-MM]{}with random and biased bounds. For [G-MM]{}with random bounds we repeat the experiment five times and report the average over these five trials. Also, for random initialization, we do five trials using different random seeds and report the mean and standard deviation of the results. [G-MM]{}does better than CCCP under *all* initializations. It also converges to a solution with lower training objective value than CCCP. Our results show that picking bounds uniformly at random from the set of valid bounds is slightly (but consistently) better than committing to the CCCP bound. We get a remarkable boost in performance when we use a reasonable prior over bounds (the bias function of Equation \[eq:bias\_function:LSSVM\]). With $\lambda = 1$, CCCP attains accuracy of 46.6%, whereas G-MM attains 49.9%, and 56.0% accuracy with random and biased initialization respectively. Moreover, [G-MM]{}is less sensitive to initialization.

Running Time
------------

[G-MM]{}bounds make a fraction of the progress made in each bound construction step compared to the MM bound. Therefore, we would expect [G-MM]{}to take more steps before it converges. We report the number of iterations that MM and [G-MM]{}take to converge in Table \[table:CCCPvsGMM\_cntIters\]. The results for [G-MM]{}depend on the value of the progress coefficient $\eta$ which is set to match the experiments in the paper; ${  \mbox{$\medmuskip=0mu\displaystyle\eta=0.02$}}$ for the clustering experiment (Section \[sec:experiments:clustering\]) and ${  \mbox{$\medmuskip=0mu\displaystyle\eta=0.10$}}$ for the scene recognition experiment (Section \[sec:experiments:scene\_classification\]).

The overhead of the bound construction step depends on the application. For example, in the scene recognition experiment, optimizing the bounds takes orders of magnitude more than sampling them (a couple of hours vs. a few seconds). In the clustering experiment, however, the optimization step is solved in closed form whereas sampling a bound involves performing a random walk on a large graph which can take a couple of minutes to run.

  ------------- ----------------- ------------------- ------------ ------------
                                                                   
                                                       **random**   **biased**
      scene       $\lambda=0.0$           145             107           87
   recognition    $\lambda=1.0$           65               69          138
      data            forgy         $35.76 \pm 7.8$                
   clustering      rand. part.     $114.98 \pm 12.9$               
    (GMM-200)    [$k$-means++]{}    $32.92 \pm 5.8$                
      data            forgy        $37.18 \pm 12.1$                
   clustering      rand. part.     $65.14 \pm 18.7$                
     (Cloud)     [$k$-means++]{}    $21.3 \pm 4.1$                 
  ------------- ----------------- ------------------- ------------ ------------

  : Comparison of the number of iterations it takes for MM and [G-MM]{}to converge in the scene recognition and the data clustering experiment with different initializations and/or datasets. The numbers reported for the clustering experiment are the average and standard deviation over 50 trials.[]{data-label="table:CCCPvsGMM_cntIters"}

Conclusion
==========

We have introduced Generalized Majorization-Minimization ([G-MM]{}), a generic iterative bound optimization framework that generalizes upon Majorization-Minimization (MM). Our key observation is that MM enforces an overly-restrictive touching constraint in its bound construction mechanism, which is inflexible and can lead to sensitivity to initialization. By adopting a different measure of progress, in [G-MM]{}, this constraint is relaxed to allow for more freedom in bound construction. Specifically, we propose deterministic and stochastic ways of selecting bounds from a set of valid ones. This generalized bound construction process tends to be less sensitive to initialization, and enjoys the ability to directly incorporate rich application-specific priors and constraints, without modifications to the objective function. In experiments with several latent variable models, [G-MM]{} algorithms are shown to significantly outperform their MM counterparts.

Future work includes applying [G-MM]{}to a wider range of problems and theoretical analysis, such as convergence rate. We also note that, although [G-MM]{}is more conservative than MM in moving towards nearby local minima, it still requires *making progress* in every step. Another interesting research direction is to enable [G-MM]{}to occasionally pick bounds that do not make progress with respect to the solution of the previous bound, thereby making it possible to get out of local minima, while still maintaining the convergence guarantees of the method.

Proof of Theorem \[thrm:GMM-converges-to-local-min\]
====================================================

We have proved two theorems in the main paper. In Theorem 1, we showed that $\lim_{t \rightarrow \infty} d_t = 0$, where $d_t = b_t(w_t) - F(w_t)$ is the gap between the [G-MM]{}bound and the objective function at iteration $t$ at the minimizer of the bound. In Theorem 2, we showed that [G-MM]{}converges to a solution, or $\exists w^\dagger$ [*s.t. *]{}$\lim_{t \rightarrow \infty} w_t = w^\dagger$.

In this appendix we show that under under mild conditions, our [G-MM]{} framework converges to a local extremum or a saddle point of the objective function. Theorem \[thrm:GMM-converges-to-local-min\] roughly states that, if $F$ is continuously differentiable, and the family of bounds are smooth and have bounded curvature, then [G-MM]{}converges to a stationary point; that is, $\nabla F(w^\dagger) = 0$.

We briefly overview the sketch of the proof first before explaining the steps in detail. We prove the theorem by contradiction by showing that if $\nabla 
    F(w^\dagger) \neq 0$ then we can construct a lower bound on $F$ (stated in (\[eq:lower-bound\_gz\])) that, at some point $w$, goes above an upper bound of $b_t$ (stated in (\[eq:upper-bound\_bt\_v4\])). This implies that $\exists w \text{ {\emph{s.t.\ }}} F(w) > b_t(w)$ which contradicts with the assumption that $b_t$ is an upper bound of $F$. The proof has three steps: 1) constructing the lower bound on $F$, 2) constructing the upper bound on $b_t$, and 3) finding a point $w$ [*s.t. *]{}$F(w) > b_t(w)$. Details are presented below.

**Step 1:** Assume ${  \mbox{$\medmuskip=0mu\displaystyle||\nabla F(w^\dagger)||_2 \neq 0$}}$. This can be stated as $$\begin{aligned}
        \exists {A}> 0 \text{ {\emph{s.t.\ }}} ||\nabla F(w^\dagger)||_2 = 2{A}.
        \label{eq:wrong_assumption}
    \end{aligned}$$ Consider value of $F$ along the gradient direction at $w^\dagger$. We denote this by $g(z)$ where $z \in \mathbb{R}$ and define it as $$\begin{aligned}
        g(z) &= F(w^\dagger + zu), \quad 
        u = \frac{\nabla F(w^\dagger)}{|| \nabla F(w^\dagger) ||_2}. \label{eq:gz}
    \end{aligned}$$ In what follows, we only consider the case where $z\geq 0$. We can write the following equalities: $$\begin{aligned}
        g'(z) &= \frac{\partial g}{\partial z}(z) 
        = \nabla F(w^\dagger + zu) \cdot u \label{eq:partial_gz} \\
        g'(0) &= ||\nabla F(w^\dagger)||_2 \label{eq:partial_g0} = 2 {A}\\
        g(0) &= F(w^\dagger). \label{eq:g0}
    \end{aligned}$$

In Figure \[fig:F\_along\_gradient:F\] we show a non-convex function $F$ in $\mathbb{R}^2$. The gradient direction at a particular point $p$ is also shown in the figure. Figure \[fig:F\_along\_gradient:F-top\] shows the same function from the top view, together with the gradient direction at $p$. Figure \[fig:F\_along\_gradient:F-along-u\] shows how the value of $F$ changes along the tangent line as a function of the step size $z$ along the gradient direction $u$ (see $g(z)$ in Equation \[eq:gz\]).

[.35]{} ![image](images/proof/F){width="1\linewidth"}

[.3]{} ![image](images/proof/F_top){width="1\linewidth"}

[.3]{} ![image](images/proof/F_along-u_grid_new){width="1\linewidth"}

Since $F$ is differentiable, we can use the mean value theorem to rewrite $g$ for $z \geq 0$ as $$\begin{aligned}
        \forall z \geq 0, \exists k \in [0, z] \text{ {\emph{s.t.\ }}} g(z) = g(0)+g'(k) z.
        \label{eq:Fwdag-zu_from_mean_value_theorem}
    \end{aligned}$$

Further, since we assume $F$ to be continuously differentiable, so is $g$. Continuity of $g'$ implies that in a small neighborhood around $0$, $g'$ is larger than ${A}$: $$\begin{aligned}
        \exists \delta > 0 \text{ {\emph{s.t.\ }}} \forall z \in [0, \delta], \; 
        g'(z) > {A}.
        \label{eq:gradient_value_at_ball_around_wdagger}
    \end{aligned}$$ We can use this inequality in (\[eq:Fwdag-zu\_from\_mean\_value\_theorem\]) and get a lower bound on $g$ in the neighborhood of $0$, or equivalently, a lower bound on $F$ in the vicinity of $w^\dagger$ and in the direction of $u$: $$\begin{aligned}
        & \exists \delta>0, \forall z \in [0, \delta], \;  g(z) > g(0) + {A}z  \notag\\
        \Rightarrow~ & \exists \delta>0, \forall z \in [0, \delta], \;  F(w^\dagger + z u) > F(w^\dagger) + {A}z.
        \label{eq:lower-bound_gz}
    \end{aligned}$$ This lower bound is shown in Figure \[fig:F\_along\_gradient:F-along-u\] in green.

**Step 2:** To get an upper bound on $b_t$, we note that $\nabla b_t(w_t) = 0$, since $w_t = \operatorname*{argmin}_w b_t(w)$ is the minimizer of $b_t$ and $b_t$ is differentiable. From this and the assumption that the curvature of all bounds in all directions is finite ($\forall b \in \mathcal{F}, \nabla^2 b(w) \preceq M I$) we can upper bound $b_t$ with a quadratic function with curvature $M$ whose minimum occurs at $w = w_t$ and has value $b_t(w_t)$, as follows: $$\begin{aligned}
         b_t(w)\leq (w - w_t)^T M I (w - w_t) + b_t(w_t), \forall w.
        \label{eq:upper-bound_bt}
    \end{aligned}$$ Similar to step 1, we only focus on the points $w$ that lie on the line $w^\dagger + zu$ and get the following bound on $b_t$: $$\begin{aligned}
        b_t(w^\dagger + zu) & \leq M || w^\dagger - w_t +zu ||^2 + b_t(w_t) \\
        & \leq 2M || w^\dagger - w_t ||^2 + 2M z^2 + b_t(w_t), \forall z\in\mathbb{R}.
        \label{eq:upper-bound_bt_v2}
    \end{aligned}$$ From Young’s inequality in multiple dimensions we have that $\forall a, b \in \mathbb{R}^d, ||a+b||^2 \leq 2||a||^2 + 2||b||^2$. We can set $a = w^\dagger - w_t$ and $b = zu$ to get $$\begin{aligned}
        2M || w^\dagger - w_t ||^2 + 2M z^2 + b_t(w_t) \geq b_t(w^\dagger + zu).
        \label{eq:upper-bound_bt_v3}
    \end{aligned}$$ Finally, from Theorem 2 we have that $\forall \epsilon_1 > 0$, $\exists T_1$ [*s.t. *]{}$\forall t > T_1$, $||w^\dagger - w_t||_2 < \epsilon_1$. We can use this in (\[eq:upper-bound\_bt\_v2\]) to get the following strict upper bound: $$\begin{aligned}
         b_t(w^\dagger + zu) < 2M \epsilon_1^2 + 2M z^2 + b_t(w_t) ,\forall t>T_1, \forall z\in\mathbb{R}.
        \label{eq:upper-bound_bt_v4}
    \end{aligned}$$

**Step 3:** We show that $\exists z \in [0, \delta]$ for which the right hand side of (\[eq:lower-bound\_gz\]) is larger than the right hand side of (\[eq:upper-bound\_bt\_v4\]), which leads to the contradiction that for $w = w^\dagger + zu$, $F(w) > b_t(w)$. In other words, we need to find $z \in [0, \delta]$ such that $F(w^\dagger) + {A}z > 2M \epsilon_1^2 + 2M z^2 + b_t(w_t)$, or $$\begin{aligned}
        2M z^2 - {A}z + 2M \epsilon_1^2 + b_t(w_t) - F(w^\dagger) < 0.
        \label{eq:find_z}
    \end{aligned}$$

${  \mbox{$\medmuskip=0mu\displaystyle\forall \epsilon_2 > 0, \exists T_2 \text{ {\emph{s.t.\ }}} \forall t > T_2, 
        b_t(w_t) - F(w^\dagger) < \epsilon_2$}}$. \[lem:bt-wt\_gets-close-to\_F-wdagger\]

To prove the lemma we first show that $$\begin{aligned}
            \lim_{t \rightarrow \infty} |b_t(w_t) - F(w^\dagger)| = 0,
            \label{eq:bt_converges-to_Fwdagger}
        \end{aligned}$$ and since $b_1(w_1), b_2(w_2), \dots$ is a non-increasing sequence, $\forall w$, $b_t(w_t) \geq F(w^\dagger)$, thus we can drop the absolute value.

To show (\[eq:bt\_converges-to\_Fwdagger\]), we can write $$\begin{aligned}
            {  \mbox{$\medmuskip=0mu\displaystyle|b_t(w_t) - F(w^\dagger)|$}} 
                        &\leq |b_t(w_t) - F(w_t)| + |F(w_t) - F(w^\dagger)|.
            \notag         \end{aligned}$$ The first term on the right hand side can get arbitrarily close to zero due to Theorem 1. The second term, also, can get arbitrarily close to zero due to Theorem 2 and that $F$ is Lipschitz continuous (because it is continuously differentiable).

Applying lemma \[lem:bt-wt\_gets-close-to\_F-wdagger\] to inequality (\[eq:find\_z\]), now we only need to show that $$\begin{aligned}
        \exists z\in [0,\delta],~ 2Mz^2 - {A}z + 2M \epsilon_1^2 + \epsilon_2 < 0.
        \label{eq:find_z_v2}
    \end{aligned}$$ Note that $\epsilon_1$ and $\epsilon_2$ can be chosen arbitrarily. In particular, we can set them so that the discriminant of the quadratic function in (\[eq:find\_z\_v2\]) is positive: $$\begin{aligned}
        \Delta = {A}^2 - 4 \left( 2M \right) \left( 2M\epsilon_1^2 + \epsilon_2 \right) > 0 
        \Rightarrow 2M \epsilon_1^2 + \epsilon_2 < \frac{{A}^2}{8M}. 
        \label{eq:discriminant}
    \end{aligned}$$ This guarantees that the quadratic function in (\[eq:find\_z\_v2\]) has two [distinct]{} roots. If $\epsilon_1=\epsilon_2=0,2M\epsilon_1^2 + \epsilon_2 = 0$, the two roots are $z_1 = 0 $ and $z_2 = \frac{{A}}{2M} > 0$. For $\epsilon_1,\epsilon_2>0$, only the constant $2M\epsilon_1^2 + \epsilon_2$ is affected which shifts the quadratic function up. We can now control $\epsilon_1$ and $\epsilon_2$ to drive the constant arbitrarily close to $0$, causing $z_1 > 0$ to get arbitrarily close to zero, and thus fall within the interval $[0,\delta]$. Therefore (\[eq:find\_z\]) can be satisfied.

To summarize, if we assume $\|\nabla F(w^\dagger)\|_2=2{A}>0$, then we have a lower bound on $F$ along the gradient direction $u$ at $w^\dagger$ (step 1): $\exists \delta>0,\forall z\in [0,\delta],F(w^\dagger+zu)>F(w^\dagger)+{A}z$. Along the same direction, $b_t$ can be upper bounded for sufficiently large $t$ (step 2): $\forall \epsilon_1>0,\exists T_1 ~s.t.~ \forall t>T_1, b_t(w^\dagger+zu)<2M\epsilon_1^2+2Mz^2+b_t(w_t)$. Finally, we can find $z_1\in [0,\delta]$ such that $F(w^\dagger)+{A}z_1>2M\epsilon_1^2+2Mz_1^2+b_t(w_t)$ (step 3), leading to $F(w^\dagger + z_1u) > b_t(w^\dagger + z_1u)$. However, this violates the assumption that $b_t$ is an upper bound of $F$, therefore $|| \nabla F(w^\dagger) ||_2 = 0$.

![The quadratic function in the left hand side of (\[eq:find\_z\_v2\]). When $\epsilon_1 = \epsilon_2 = 0$, the quadratic function has two roots, $z_1=0$ and $z_2=\frac{A}{2M}>0$ (shown in blue). When $\epsilon_1,\epsilon_2>0$ but properly chosen, the quadratic function is shifted up but still has two distinct roots (shown in red). []{data-label="fig:quadratic_z"}](images/proof/quadratic){width="0.5\linewidth"}
