A Riemannian low-rank method for optimization over semidefinite matrices with block-diagonal constraints

Introduction

and the set of $Y$’s obtained by stacking as

This paper is concerned with solving optimization problems of the form

with the additional constraint $\operatorname{rank}(X)\leq p$. As often, the rank constraint is the culprit. Indeed, continuing with the Max-Cut example (linear $f$), optimization over $\mathcal{C}$ without the rank constraint is a semidefinite program, which can be solved to arbitrary precision in polynomial time .

This is motivation to study the relaxation obtained by ignoring the rank constraint (for linear $f$, it is also the dual of the dual of ($\textrm{RP}_{p}$)):

For linear $f$, solutions of (P) can of course be computed using standard SDP solvers, such as interior point methods (IPM). Unfortunately, as demonstrated in Section 5, IPM’s do not scale well. The main reason for it is that, as the name suggests, IPM’s iterate inside the interior of the search space $\mathcal{C}$. The latter is formed by full-rank, dense matrices of size $n$: this quickly becomes unmanageable.

The full-rank operations seem even more wasteful considering that, still for linear $f$, problem (P) always admits a solution of rank at most

Indeed, this follows a general result of Shapiro , Barvinok and Pataki regarding extreme points of spectrahedra,The name spectrahedron for the search space of a semidefinite program echoes the name polyhedron for the search space of a linear program. that is, intersections of the positive semidefinite cone with an affine subspace—the geometry of $\mathcal{C}$ is discussed in Section 2.2. This prompted Burer and Monteiro to propose SDPLR, a generic SDP solver which exploits the low-rank phenomenon. Applying SDPLR to our problem amounts to computing a local minimizer $Y$ of ($\textrm{RP}_{p}$) for some small $p$, using classical nonlinear optimization algorithms and penalizing for the constraints in a Lagrangian way. Then, $p$ is increased as needed until $YY^{\top}\!$ can be certified as a solution to the SDP.

SDPLR is powerful and generic, and the theory accompanying the algorithm brings great insight into the problem. But it also has some downsides we want to improve on in the context of (P). First, it is not an easy matter to guarantee convergence to (even local) optimizers in the nonlinear subproblems. Furthermore, since constraints are enforced by penalization, they are not accurately satisfied by the returned solution. Finally, we would like to allow for nonlinear $f$. Nevertheless, Section 5 shows SDPLR improves significantly upon IPM’s.

In both and , one of the keys to practical efficiency is (well-justified) optimism: ($\textrm{RP}_{p}$) is first solved for small values of $p$, and $p$ is increased only as needed. In both papers, it is observed that, in practice, it often suffices to reach $p$ just above the rank of the target solution of (P), which may be quite small; but there is no theory to confirm this. We do not prove such a strong result either, but we give some nontrivial, deterministic bounds on “how high one must lift”, refining certain results of .

We then turn our attention to computing Karush-Kuhn-Tucker (KKT) points for (P). These are points that satisfy first-order necessary optimality conditions. If $f$ is convex, the conditions are also sufficient. Our goal is to compute KKT points via the computation of second-order critical points of ($\textrm{RP}_{p}$), which is lower-dimensional. A key property that makes this possible is the availability of an explicit dual matrix $S(X)$ (21) which intervenes in both sets of conditions.

Using this dual matrix, we show that rank-deficient second-order critical points $Y$ reveal KKT points $X=YY^{\top}\!$. Furthermore, when a computed second-order critical point is full rank, it is shown how to use it as a warm-start for the computation of a second-order critical point of ($\textrm{RP}_{p}$) with a larger value of $p$. It is guaranteed that if $p$ is allowed to grow up to $n$, then all second-order critical points reveal KKT points, so that the procedure terminates. This is formalized in Algorithm 1, which we call the Riemannian Staircase, as it lifts ($\textrm{RP}_{p}$) to (P) step by step, instead of all at once.

The above points rest extensively on work discussed earlier in this introduction , and improve upon those along the lines announced in the same. In particular, we do not require the computation of local optimizers of ($\textrm{RP}_{p}$), and we avoid the geometry breakdown tied to the quotient approach in . We also stress that the latter reference only covers $d=1$, and SDPLR only covers linear $f$.

We further take particular interest in understanding how large $p$ may grow in the staircase algorithm. We view this part as our principal theoretical contribution. This investigation calls for inspection of the convex geometry of $\mathcal{C}$, with particular attention to its faces and their dimension. To this effect, we use results by Pataki to describe the face of $\mathcal{C}$ which contains a given $X$ in its relative interior, and we quote the lower-bound on the dimension of that face as a function of $\operatorname{rank}(X)$. We further argue that this bound is almost always tight, and we give an essentially tight upper bound on the dimension of a face, generalizing a result of Laurent and Poljak to $d>1$.

Using this facial description of $\mathcal{C}$, we establish that for strongly concave $f$, for $p>p^{*}$ (3), all second-order critical points of ($\textrm{RP}_{p}$) reveal KKT points. Also, for concave $f$, we show the same for $p>\frac{d+1}{d+3}n$ (Corollary 3.16), and argue that $p>p^{*}$ is sufficient under an additional condition we believe to be mild. Hence,

For linear $f$, above a certain threshold for $p$, all second-order critical points of ($\textrm{RP}_{p}$) are global optimizers.

The condition is stronger than the one proposed in , and the statement is about second-order critical points, rather than about local optimizers of ($\textrm{RP}_{p}$). There are no similar results for convex $f$, as then solutions can have any rank.

We close the paper with numerical experiments showing the efficiency of the staircase algorithm to solve (P) on certain synchronization problems involving rotations and permutations, as compared to IPM’s and SDPLR.

Note that, up to a linear change of variable, problem (P) also encompasses constraints of the form $X_{ii}=B_{i}$ where each $B_{i}$ is positive definite. We assume all diagonal blocks have identical size $d$ as this simplifies exposition, but the proposed method can easily accommodate inhomogeneous sizes, and many of the developments go through for complex matrices as well.

2 Applications

Problem ($\textrm{RP}_{p}$) and its relaxation (P) appear in numerous applications. Many of those belong to the class of synchronization problems, which consist in estimating group elements from measurements of pairwise ratios. Further applications are also described, e.g., in .

can be modeled in ($\textrm{RP}_{p}$) with $d=p=1$. A seminal example is Max-Cut: the problem of clustering a graph in two classes, so as to maximize the sum of weights of edges joining the two classes. The cost $f$ is linear, determined by the graph’s adjacency matrix. Its relaxation to (P) is the subject of an influential analysis by Goemans and Williamson , which helped popularize the type of lifts considered here. See [31, eq. (3)] for a recent application of Max-Cut to genomics. The same setup, but with different linear costs, appears in the stochastic block model , in community detection , in maximum a posteriori (MAP) inference in Markov random fields with binary variables and pairwise interactions and in robust PCA [46, Alg. 1]. All of these study the effects of the relaxation on the final outcome, mostly under random data models. Their linear cost matrices are often structured (sparse or low-rank), which is easily exploited here.

is an important biomedical imaging instance of ($\textrm{RP}_{p}$), where orthonormal matrices are to be estimated with $d=2,p=3$ .

3 Related work

Problem ($\textrm{RP}_{p}$) is an instance of optimization on manifolds . Optimization over orthonormal matrices is also studied in, e.g., . Being equivalent to (P) with a rank constraint, ($\textrm{RP}_{p}$) also falls within the scope of optimization over matrices with bounded rank , where the latter is also an extension of . The particular case of optimization over bounded-rank positive semidefinite matrices with linear constraints was already addressed in . The same without positive semidefiniteness constraint is studied recently in , also with a discussion of global optimality of second-order critical points. With a linear cost $f$, problem ($\textrm{RP}_{p}$) (which then has a quadratic cost $g$) is a subclass of quadratically constrained quadratic programming (QCQP). QCQP’s and their SDP relaxations have been extensively studied, notably in , with particular attention to approximation ratios. For (P), these approximation ratios can be found in .

In part owing to the success of (P) with linear $f$ in adequately solving a myriad of hard problems, there has been strong interest in developing fast, large-scale SDP solvers. The present paper is one example of such a solver, restricted to the class of problems (P). SDPLR is a more generic such solver . See also for a review, and for a recent low-complexity example with precise convergence results, but which does not handle constraints.

Much of this paper is concerned with characterizing the rank of solutions of (P), especially with respect to how large $p$ must be allowed to grow in ($\textrm{RP}_{p}$) to solve (P). There is also considerable value in determining under what conditions (P) admits solutions of the desired rank for a specific application, that is: when is the relaxation tight? This question is partially answered in for the closely related phase synchronization problem, under a stochastic model for the data. See for a proof in the stochastic block model, and for a study of phase transitions in random convex programs. There also exist deterministic tightness results, typically relying on special structure in a graph underlying the problem data. See for example . See also Appendix A for a deterministic proof of tightness in the case of single-cycle synchronization of rotations. The proof rests on the availability of a closed-form expression for the dual matrix $S$ (21), and for the solution to be certified. With the same ingredients, it is easy to show, for example, that (P) is tight for Max-Cut when the graph is bipartite.

Semidefinite relaxations in the form of (P) with additional constraints have also appeared in the literature. In particular, this occurs in estimation of rotations, with explicit care for the determinant constraints: Saunderson et al. explicitly constrain off-diagonal blocks to belong to the convex hull of the rotation group; this is not necessary for the orthogonal group—see Proposition 2.1. Similarly, for synchronization of permutations in joint shape matching, off-diagonal blocks are restricted to be doubly stochastic . Finally, in recent work, Bandeira et al. study a more powerful class of synchronization problems with additional linear constraints of various forms. An example with an additional nonlinear constraint appears in , which imposes an upperbound on the spectral norm of $X$. All of these are motivation to generalize the framework studied here, in future work.

We mention in passing that the MaxBet and MaxDiff problems do not fall within the scope of this paper. Indeed, although they also involve estimating orthonormal matrices as in ($\textrm{RP}_{p}$), their cost function has a different type of invariance, which would also lead to a different type of relaxation.

4 Notation

Geometry

Both search spaces of ($\textrm{RP}_{p}$) and (P) enjoy rich geometry, which leads to efficient analytical and numerical tools for the study of these optimization problems. The former is a smooth Riemannian manifold, while the latter is a compact convex set. Figure 1 depicts the two.

This allows for a simple definition of the manifold via an equality constraint as

The total computational cost of a projection is thus $\mathcal{O}(m\cdot d^{2}p)=\mathcal{O}(ndp)$ flops.

Optimization algorithms on Riemannian manifolds typically are iterative. As such, they require a means of moving away from a point $Y$ along a prescribed tangent direction $\dot{Y}$, to reach a new point on the manifold: the next iterate. Since $Y+\dot{Y}$ does not, in general, belong to the manifold, extra operations are required. Retractions achieve exactly this [4, § 4.1]. One possible retraction for (1) is as follows. For each $d\times p$ “slice” $i$ in $\{1,\ldots,m\}$,

where $U_{i}\Sigma_{i}V_{i}^{\top}\!$ is a thin singular value decomposition of the $i$th slice $Y_{i}+\dot{Y}_{i}$. This retraction projects each slice of $Y+\dot{Y}$ to the closest orthonormal matrix. Consequently, this is even a second-order retraction . The total cost of computing a retraction is $\mathcal{O}(m\cdot(p^{2}d+d^{3}))=\mathcal{O}(np^{2})$ flops.

2 Convex geometry of the full search space

The optimization problem (P) is defined over the compact convex set

Since $X$ is positive semidefinite, for all $i\neq j$, the submatrix formed by the blocks $X_{ii},X_{ij},X_{ji}$ and $X_{jj}$ is positive semidefinite. By Schur, this holds if and only if $X_{ij}^{\top}\!X_{ij}\preceq I_{d}$, which in turn happens if and only if all singular values of $X_{ij}$ are at most 1. The set of such matrices is the convex hull of all orthogonal matrices of size $d$ .

The set $\mathcal{C}$ may be decomposed into faces of various dimensions.

A face of $\mathcal{C}$ is a convex subset $\mathcal{F}$ of $\mathcal{C}$ such that every (closed) line segment in $\mathcal{C}$ with a relative interior point in $\mathcal{F}$ has both endpoints in $\mathcal{F}$. The empty set and $\mathcal{C}$ itself are faces of $\mathcal{C}$.

By [56, Thm. 18.2], the collection of relative interiorsThe relative interior of a singleton is the singleton. of the non-empty faces forms a partition of $\mathcal{C}$. That is, each $X\in\mathcal{C}$ is in the relative interior of exactly one face of $\mathcal{C}$, called $\mathcal{F}_{X}$ . Furthermore, all faces of $\mathcal{C}$ are exposed [55, Cor. 1], that is, for every face $\mathcal{F}$, there exists a linear function $f$ such that $\mathcal{F}$ is the set of solutions of (P). Of particular interest are the zero-dimensional faces of $\mathcal{C}$ (singletons), called its extreme points.

$X\in\mathcal{C}$ is an extreme point of $\mathcal{C}$ if there does not exist $X^{\prime},X^{\prime\prime}\in\mathcal{C}\backslash\{X\}$ and $0<\lambda<1$ such that $X=\lambda X^{\prime}+(1-\lambda)X^{\prime\prime}$. $X$ is an exposed point of $\mathcal{C}$ if there exists $C$ such that $X$ is the unique maximizer of $\left\langle{C},{X}\right\rangle$ in $\mathcal{C}$.

In other words, $X$ is extreme if it does not lie on an open line segment included in $\mathcal{C}$. Since $\mathcal{C}$ is compact, it is the convex hull of its extreme points [56, Cor. 18.5.1]. Extreme points are of interest notably because they often arise as the solution of optimization problems. Specifically, if $f$ is a concave function (in particular, if $f$ is linear), then $f$ attains its minimum on $\mathcal{C}$ at one of its extreme points [56, Cor. 32.3.2].

The dimension of $\mathcal{F}_{X}$ is the dimension of the kernel of $\mathcal{L}_{X}$. The rank-nullity theorem gives a lowerbound (see also Theorem 3.15 for an upperbound):

It follows that extreme points $X$ (i.e., points such that $\dim\mathcal{F}_{X}=0$) have small rank:

Note that $\Delta\geq 0$ when $p\geq p^{*}$.

For linear $f$, (P) admits an extreme point as global optimizer, so that ($\textrm{RP}_{p}$) and (P) have the same optimal value as soon as $p\geq p^{*}$ (16). In other words: for linear $f$, ($\textrm{RP}_{p}$) is not NP-hard if $p\geq p^{*}$.

Not all feasible $X$’s with rank as in (16) are extreme. For example, setting $d=1$ and $m\geq 3$ as in Figure 1, select two distinct, admissible matrices of rank 1, $X_{0}$ and $X_{1}$. For all $0<\lambda<1$, the matrix $X_{\lambda}=\lambda X_{1}+(1-\lambda)X_{0}$, lying on the open line segment between $X_{0}$ and $X_{1}$, is admissible and has rank 2. Thus, $X_{\lambda}$ satisfies (16), but it is not an extreme point, by construction. Notwithstanding, the expectation that $\mathcal{L}_{X}$ is generically of full rank suggests that almost all feasible $X$’s satisfying (16) should be extreme; an intuition that is supported by Figure 1. More generally, in Theorem B.1, we prove for $d=1$ that $\dim\mathcal{F}_{X}=\Delta$ for almost all $X$ of rank $p$.

Many applications look for solutions of rank $d$. All $X$’s of rank $d$ are exposed (hence extreme), meaning they can all be recovered as unique solutions of (P).

Let $\sigma_{1}(X_{ij})\geq\cdots\geq\sigma_{d}(X_{ij})$ denote the singular values of $X_{ij}$. By Proposition 2.1, $\sigma_{k}(X_{ij})\leq 1$ for all $i,j,k$. Hence,

The upperbound is attained if and only if $\sigma_{k}(X_{ij})=1$ for all $i,j,k$, thus, if and only all $X_{ij}$’s are orthogonal. By Proposition 2.1, this is the case if and only if $\operatorname{rank}(X)=d$. Now consider $X$ has rank $d$. We show it is exposed (and hence extreme):

The second inequality follows by Cauchy-Schwarz, and equality is attained if and only if $\hat{X}=X$, which is in $\mathcal{C}$. Thus, both max problems admit $X$ as unique solution, confirming that $X$ is an exposed extreme point. ∎

Proposition 2.3 is an extension of [45, Thm. 1] to the case $d>1$. We note that the proof in that reference does not generalize to $d>1$.

From second-order critical points Y𝑌Y to KKT points X𝑋X

In this section, we show that rank-deficient second-order critical points $Y$ of ($\textrm{RP}_{p}$) yield KKT points $X=YY^{\top}\!$ of (P). Furthermore, when $Y$ is second-order critical but $X$ is not KKT, it is shown how to escape the saddle point by increasing $p$. If $p$ increases all the way to $n$, then all second-order critical points reveal KKT points. The proofs parallel those in . The main novelty is explicit bounds on $p$ such that all second-order critical points of ($\textrm{RP}_{p}$) reveal KKT points. The proofs bring us to consider the facial structure of $\mathcal{C}$.

A key ingredient for all proofs in this section is the availability of an explicit matrix $S(X)$ (21) which is positive semidefinite if and only if $X$ is KKT (Theorem 3.3). The formula for $S$ is simply read off from the first-order optimality conditions of ($\textrm{RP}_{p}$), owing to smoothness of the latter.

If $X$ is a local optimizer for (P), then $X$ is a KKT point. If $f$ is convex, all KKT points are global optimizers.

Apply Theorems 3.25 and 3.34, and Example 3.36 in . KKT conditions are necessary since Slater’s condition holds: $I_{n}$ is feasible for (P) and it is strictly positive definite. ∎

If $Y$ is a local optimizer for ($\textrm{RP}_{p}$), then it is a second-order critical point.

Lemmas 3.1 and 3.2 suggest the definition of an (as yet merely tentative) formula for the dual certificate $\hat{S}$, based on (19):

$X\in\mathcal{C}$ is a KKT point for (P) if and only if $S$ (21) is positive semidefinite. If so, $\hat{S}=S$ is the unique dual certificate for Lemma 3.1.

We show the if and only if parts of the first statement separately.

Since $X$ is a KKT point for (P), there exist $\hat{S}\succeq 0$ and $\hat{\Lambda}$ symmetric, block-diagonal satisfying the conditions in Lemma 3.1. In particular, $\hat{S}X=0$ and $\nabla f(X)=\hat{S}-\hat{\Lambda}$. Thus, $\nabla f(X)X=-\hat{\Lambda}X$ and $\operatorname{symblockdiag}\!\left({\nabla f(X)X}\right)=-\operatorname{symblockdiag}\!\left({\hat{\Lambda}X}\right)=-\hat{\Lambda}$. Here, we used both the fact that $\hat{\Lambda}$ is symmetric, block-diagonal and the fact that $X_{ii}=I_{d}$. Consequently, $S=\nabla f(X)-\operatorname{symblockdiag}\!\left({\nabla f(X)X}\right)=\hat{S}-\hat{\Lambda}+\hat{\Lambda}=\hat{S}\succeq 0$.

The last point also shows there exists only one pair $(\hat{S},\hat{\Lambda})$ certifying $X$ is a KKT point. ∎

Notice how the Riemannian structure underlying problem (P) made it possible to simply read off an analytical expression for a dual certificate from the necessary optimality conditions of ($\textrm{RP}_{p}$). This smooth geometry also leads to uniqueness of the dual certificate (this is connected to nondegeneracy [5, Thm. 7]). Theorem 3.3 makes for an unusually comfortable situation and will be helpful throughout the paper.

For convex $f$, we can make the following statement regarding uniqueness of the solution.

Assume $f$ is convex. If $X\in\mathcal{C}$ is an extreme point for (P) (which is true in particular if $\operatorname{rank}(X)=d$), and $S\succeq 0$, and $\operatorname{rank}(X)+\operatorname{rank}(S)=n$ (strict complementarity), then $X$ is the unique global optimizer of (P).

From Theorem 3.3, it is clear that $X$ is a global optimizer. We prove by contradiction that it is unique. Let $X^{\prime}\neq X$ be another global optimizer. Since (P) is a convex problem in this setting, $f$ is constant over the whole (optimal) segment $t\mapsto X+t(X^{\prime}-X)$ for $t\in$. Hence, the directional derivative of $f$ at $X$ along $\dot{X}=X^{\prime}-X$ is zero:

Since both $S$ and $X^{\prime}$ are positive semidefinite, it ensues that $SX^{\prime}=0$. (Note that for linear $f$, this shows $S$ is the dual certificate for all global optimizers of (P), not only for $X$.) Hence, $S\dot{X}=0$.

In general, the condition that $X$ be an extreme point in the previous theorem cannot be removed. Indeed, if $f(X)\equiv 0$, then all admissible $X$’s are globally optimal and $S(X)\equiv 0\succeq 0$. In particular, $X=I_{n}$ satisfies strict complementarity, but if $m>1$, it is not extreme, and it is not a unique global optimizer. Likewise, any rank-$d$ admissible $X$ is extreme and globally optimal, but does not satisfy strict complementarity. Similar examples can be built with nonzero linear costs $f(X)=\left\langle{C},{X}\right\rangle$, where the sparsity pattern of $C$ corresponds to a disconnected graph.

Continuing with convex cost functions, KKT points of (P) coincide with global optimizers. This and the fact that ($\textrm{RP}_{p}$) is a relaxation of (P) lead to the following summary regarding global optimality conditions. For ($\textrm{RP}_{p}$), these sufficient conditions are conclusive whenever the relaxation is tight.

Since ($\textrm{RP}_{p}$) is essentially equivalent to (P) with the additional constraint $\operatorname{rank}(X)\leq p$, assuming $Y$ is optimal for ($\textrm{RP}_{p}$), we expect $X$ to be a KKT point at least if either of the following holds: (1) if $Y$ is rank deficient, since then the extra constraint is not active, meaning it is “as if” we were solving (P); or (2) if $p=n$, since then the extra constraint is vacuous. The two following theorems show this still holds for second-order critical points $Y$.

If $Y$ is a rank-deficient, second-order critical point for ($\textrm{RP}_{p}$), then $X=YY^{\top}\!$ is a KKT point for (P).

If $Y$ is square ($p=n$) and it is a second-order critical point for $(\textrm{\emph{RP}}_{n})$, then $X=YY^{\top}\!$ is a KKT point for (P). If $Y$ is full-rank, it needs only be first-order critical for $X$ to be a KKT point.

If $Y$ is rank deficient, then the result follows from Theorem 3.7. If $Y$ is full rank, then it is invertible. Since $Y$ is also a critical point, first-order optimality conditions (19) imply $SY=0$, hence $S=0$. This completes the proof, as per Theorem 3.3. ∎

In particular, if $f$ is convex and (P) has a unique solution of rank $r$, then all second-order critical points of ($\textrm{RP}_{p}$) have rank either $r$ or $p$. Thus, if $p$ is larger than the rank of a solution of (P), we may hope that minimizing ($\textrm{RP}_{p}$) until we reach a second-order critical point will result in a rank-deficient $Y$, revealing a KKT point. Unfortunately, in general, we cannot guarantee rank deficiency beforehand. For those cases, the following theorem and corollary provide a means of escaping unsatisfactory critical points.

By Theorem 3.3, $u$ exists because $X$ is not a KKT point. Since $Y\dot{Y}^{\top}\!=0$, $\dot{Y}$ is indeed tangent at $Y$ (6). From (22), it follows that

Since $u^{\top}\!Su<0$, there exists $t_{0}>0$ such that $\phi(t)<\phi(0)$ for all $t$ in $]0,t_{0}[$.∎

Since $Y_{+}Y_{+}^{\top}\!=YY^{\top}\!$, $Y_{+}$ is indeed feasible for $(\text{RP}_{p_{+}})$ and $g(Y_{+})=g(Y)$. Since $Y$ is critical, $SY=0$, so $SY_{+}=0$: $Y_{+}$ is a critical point. The rest follows from Theorem 3.9. ∎

Later in this section, we show how to escape full-rank points without increasing the rank, under additional assumptions—see Proposition 3.17.

An important question remains: for moderate $p$, should we expect to encounter second-order critical points $Y$ that do not correspond to KKT points of (P)? We provide partial answers for concave $f$ below. In particular, this covers the important case of linear costs. The stronger results do not include strictly convex functions, as for these (P) can have solutions of arbitrary rank.

The result below shows that $Y$ can only be a critical point of ($\textrm{RP}_{p}$) if $YY^{\top}\!$ is a critical point for (P) restricted to the face $\mathcal{F}_{X}$; and similarly for second-order critical points. This brings a useful corollary.

Thus, $X$ is a second-order critical point for the face-restricted optimization problem. ∎

The following corollary can be put in perspective with [20, Thm. 3.4]. The latter states a similar result for (P) with general linear equality constraints, for linear $f$. Their result characterizes local optimizers of ($\textrm{RP}_{p}$), whereas the following result characterizes first- and second-order critical points (computationally more manageable objects).

If $f$ is linear and $Y$ is critical, then $f$ is constant over $\mathcal{F}_{X}$. If furthermore $Y$ is second-order critical and $p>\lfloor p^{*}\rfloor$ (16), then either $X$ is globally optimal for (P), or ($Y$ has full rank and) the face $\mathcal{F}_{X}$ has positive dimension (15) and is suboptimal.

If $f$ is convex (resp., strictly convex) and $Y$ is critical, then $X$ is optimal (resp., the unique optimizer) for (P) restricted to $\mathcal{F}_{X}$.

For linear $f$, Corollary 3.12 is not quite sufficient to determine how large $p$ must be to exclude “bad” second-order critical points. Paraphrasing the comment following [20, Thm. 3.4], the latter showed that, for linear $f$ and $p>p^{*}$ (essentially), local optima of ($\textrm{RP}_{p}$) are global optima, with the caveat that positive-dimensional faces (over which $f$ must be constant) may harbor non-global local optima. In the literature, this has sometimes been quoted as saying that local optima are global optima if $f$ is not constant over any proper face of $\mathcal{C}$ [see, e.g., 46, footnote 3], but there is no indication that this is a mild condition.In fact, we found in numerical experiments (not reported) that, for $d=1$, $\Delta=1$ (thus, almost all faces have dimension 1 and $p>p^{*}$) and a random linear cost $\langle{C},{X}\rangle$, we could easily find a face $\mathcal{F}_{X}$ of dimension 1 over which the cost is constant but not optimal.

We thus set out to further refine the implications of second-order criticality of $Y$. We do so by leveraging the tight relationship (22) between the Hessian of the cost on ($\textrm{RP}_{p}$) and the dual certificate $S$.

negative eigenvalues. Hence, if $X$ is not a KKT point for (P), then it has rank $p$ and $\dim\mathcal{F}_{X}\geq\Delta+p$.

In particular, since $M\succeq 0$, we have $0\leq\mu_{0}\leq\lambda_{\lfloor(np-k^{\prime})/p\rfloor}$. It remains to determine $k^{\prime}$.

Combine with $\lambda_{\lfloor(np-k^{\prime})/p\rfloor}\geq 0$ and the definition of $\Delta$ (15) to conclude. ∎

Theorem 3.13 is particularly meaningful for linear $f$, considering the intuition that for $p\geq p^{*}$, generically, $\dim\mathcal{F}_{X}=\Delta\geq 0$ (15) (Theorem B.1 gives a proof for $d=1$). Thus, for such $p$, a second-order critical point is either globally optimal, or it maps to a face of abnormally high dimension, over which $f$ must be constant and suboptimal. We could not produce an example of the latter. We summarize this in a corollary, followed by a question.

Use Corollary 3.12, Theorem 3.13 and Theorem 3.7. ∎

For linear $f(X)=\langle{C},{X}\rangle$ and $p>p^{*}$, the question is the following: if $C$ is sampled uniformly at random from the unit-norm symmetric matrices, what is the probability that $f$ is constant over a face $\mathcal{F}_{X}$ of dimension $\Delta+p$ or larger, with $\operatorname{rank}(X)=p$? If it is zero, then almost surely all second-order critical points of ($\textrm{RP}_{p}$) are global optimizers. We do not answer this question here, but refer to Theorem B.1 to argue that there are few such faces.

Theorem 3.13 is motivation to investigate upper-bounds on the dimensions of faces of $\mathcal{C}$. The following result extends [42, Thm. 3.1(i)] to $d\geq 1$.

If $X\in\mathcal{C}$ has rank $p$, then the face $\mathcal{F}_{X}$ (13) has dimension bounded as:

If $p$ is an integer multiple of $d$, the upperbound is attained for some $X$.

The constraint $\mathcal{L}_{X}(A)=0$ means $\left\langle{A},{E_{ij}}\right\rangle=0$ for each $k$ and for each $i,j\in s_{k}$. To establish the theorem, we need to extract a subset $T$ of at least $p(d+1)/2$ of these $md(d+1)/2$ constraint matrices, and guarantee their linear independence. To this end, let

That is, for each slice $k$, $T$ includes all constraints of that slice which involve at least one of the selected rows. For each slice $k$, there are $|c_{k}|d-\frac{|c_{k}|(|c_{k}|-1)}{2}$ such constraints—note the correction for double-counting the $E_{ij}$’s where both $i$ and $j$ are in $c_{k}$. Thus, using $|c_{1}|+\cdots+|c_{m}|=p$, the cardinality of $T$ is:

We first show matrices in $T$ are linearly independent. Then, we show $|T|$ is large enough.

It remains to lowerbound (28). To this effect, use $|c_{k}|\leq d$ to obtain:

Indeed, the maximum is attained by making as many of the entries of $x$ as large as possible—this can be verified using KKT conditions. In combination with (28), this confirms at least $p(d+1/2)-pd/2=p(d+1)/2$ linearly independent constraints act on $A$, thus upperbounding $\dim\mathcal{F}_{X}$.

To conclude, we argue that the proposed upperbound is essentially tight. Indeed, build $Y$ by repeating $m$ times the $d$ first rows of $I_{p}$, then by replacing its $p$ first rows with $I_{p}$ (to ensure $Y$ is full-rank). If $p/d$ is an integer, then exactly the $p/d$ first slices each contribute $d(d+1)/2$ independent constraints, i.e., $\dim\mathcal{F}_{YY^{\top}\!}=p(p+1)/2-p(d+1)/2$. ∎

Theorems 3.3, 3.13 and 3.15 combined give a sufficient condition on $p$ to ensure all second-order critical points of ($\textrm{RP}_{p}$) correspond to KKT points of (P).

Since $\operatorname{rank}(X)\leq p$, we have $\dim\mathcal{F}_{X}-\Delta\leq(n-p)\frac{d+1}{2}$. Theorem 3.13 then gives $S\succeq 0$ (21) if $(n-p)(d+1)<2p$, which is the case. Apply Theorem 3.3 to conclude. ∎

In particular, for the Max-Cut SDP ($d=1$, $f$ linear), this shows that computing a second-order critical point of ($\textrm{RP}_{p}$) with $p=\left\lfloor n/2\right\rfloor+1$ certainly solves (P). This is an interesting and new result, but of course, in practice, it is desirable (and empirically sufficient) to take $p=\left\lfloor p^{*}\right\rfloor+1$ (much smaller). In the unlikely event we would encounter a “bad” second-order critical point with such $p$, the following theorem provides an escape route (for concave $f$) which does not require increasing the rank. It proceeds by moving inside a face.

Recall the definitions of $\mathcal{F}_{X}$ (13) and $\mathcal{L}_{X}$ (14). All follows from $\mathcal{H}=\mathcal{L}_{X}^{*}\mathcal{L}_{X}$, where $\mathcal{L}_{X}^{*}$ is the adjoint of $\mathcal{L}_{X}$. ∎

The latter proposition suggests an explicit numerical method to compute $A$, by computing a minimal eigenvector of $\mathcal{H}$. Applying $\mathcal{H}$ costs $\mathcal{O}(m(d^{2}p+p^{2}d))$ flops. Assuming $p=\lfloor p^{*}\rfloor+1=\Theta(d\sqrt{m})$ and that up to $p(p+1)/2$ applications are necessary, this brings the cost of computing $A$ to $\mathcal{O}(d^{2}n^{3})$ flops.

The Riemannian staircase algorithm

The above results suggest a simple algorithm to compute KKT points of (P): for some small value of $p\geq d+1$, find a second-order critical point $Y$ of ($\textrm{RP}_{p}$). If $Y$ is rank deficient, then Theorem 3.7 guarantees $X=YY^{\top}\!$ is KKT for (P). Otherwise, increase $p$ and find a second-order critical point of $(\text{RP}_{p_{+}})$, possibly warm-starting as suggested by Corollary 3.10. Iterating this procedure, the worst-case scenario is when $p$ increases all the way to $n$, in which case any second-order critical point of $(\text{RP}_{n})$ yields a KKT point of (P), as per Theorem 3.8. Specific results pertaining to classes of functions $f$ limit how large $p$ could grow. We call this the Riemannian Staircase, listed as Algorithm 1. Of course, the hope is that the algorithm returns for some small $p$, and in practice we find that it is often sufficient to take $p$ just above the rank of a solution.

Inside the else-block, the augmented $Y_{i}$ (with additional columns of zeros) is (usually) a saddle point. Although the second-order procedure RiemannianOptimization should be able to escape it, we make this step explicit via the procedure EscapeDirection. The latter can be implemented using Corollary 3.10, which indicates how computing an eigenvector of $S$ (21) associated to its smallest eigenvalue, combined with a line-search, allows to escape the saddle with strict cost decrease (unless that eigenvalue is nonnegative, in which case $Z=0$ and $Y_{i}$ is returned with $Y_{i}Y_{i}^{\top}\!$ being KKT).

For all sufficiently smooth $f$, taking $p_{k}=n$ guarantees Algorithm 1 returns $Y$ such that $YY^{\top}\!$ is a KKT point. For convex $f$, KKT points may have arbitrary rank, so that allowing large $p_{k}$ seems necessary in general. For strongly concave $f$, it is sufficient to take $p_{k}=\lfloor p^{*}\rfloor+1$ (Corollary 3.12) ; for concave (and linear) $f$, it is sufficient to take $p_{k}=\lfloor\frac{d+1}{d+3}n\rfloor+1$ (Corollary 3.16), and it is expected that $p_{k}=\lfloor p^{*}\rfloor+1$ should be sufficient (Corollary 3.14 and discussion).

In the latter case, in the unlikely event that Algorithm 1 terminates with $Y$ of size $n\times p_{k}$, $p_{k}\geq\lfloor p^{*}\rfloor+1$, full-rank and second-order critical such that $X=YY^{\top}\!$ is not a KKT point of (P), it is possible to further optimize without increasing the rank. Indeed, since $\dim\mathcal{F}_{X}>0$, Proposition 3.17 shows how to compute $Y^{\prime}$ such that $X^{\prime}=Y^{\prime}(Y^{\prime})^{\top}\!$ is on the boundary of $\mathcal{F}_{X}$. Since $f$ is concave, $f(X^{\prime})\leq f(X)$ (Lemma 3.11). $Y^{\prime}$ is critical and rank-deficient. If $Y^{\prime}$ is second-order critical, $X^{\prime}$ is KKT. Otherwise, Theorem 3.9 shows how to escape with a strict cost decrease. Iterating this procedure as needed, the cost decreases strictly (no cycling), and the rank $p$ never exceeds $p_{k}$. We expect this procedure to terminate since (P) admits KKT points of rank at most $\lfloor p^{*}\rfloor$, but we do not prove this.

In practice, for the RiemannianOptimization procedure, we use the Riemannian trust-region method (RTR) , through the Manopt toolbox . RTR is a descent method. It converges toward critical points regardless of the initial iterate (global convergence).If the local optimizers of $g$ were isolated, we could also guarantee local convergence at a quadratic rate, but $g(Y)=g(YQ)$ for all orthogonal $Q$, so this is never the case. In practice though, we do observe a characteristically superlinear convergence. Furthermore, the stable fixed points of RTR are local optimizers, thus making convergence to points which are not second-order critical unlikely (but not impossible). Should this happen, Theorem 3.9 shows how to escape. Admittedly, it is unclear how many times this might have to be repeated in the worst case.

Ideally, one would modify the RTR algorithm itself to ensure global convergence to second-order critical points. To the best of our knowledge, algorithms with such properties have not yet been described in the Riemannian setting. Nevertheless, we are hopeful that this should be possible, in the light of recent work by Cartis et al. . These authors indeed describe a modification of the classical trust-region method and guarantee polynomial-time convergence to approximate second-order critical points. Encouragingly, Sun et al. achieved a strong result in this vein for dictionary learning with RTR on a sphere, hinting to a possible generalization on manifolds.

Special case: linear cost function

Let $f(X)=\langle{C},{X}\rangle$ be linear and let $f^{*}$ denote the optimal value of (P). Then, for all $X\in\mathcal{C}$,

Algorithm 1 solves the SDP by optimizing $g$ in ($\textrm{RP}_{p}$), whose differentials are:

As an illustrative example, we here apply Algorithm 1 and competing SDP solvers to random instances of the orthogonal synchronization problem . In this setting, one wishes to estimate $m$ orthogonal matrices $Q_{1},\ldots,Q_{m}$, based on noisy measurements of the relative transformations $Q_{i}Q_{j}^{\top}\!$. See the introduction for applications.

In this benchmark, for increasing values of $m$, target matrices of size $d=3$ are generated uniformly at random. The measurements of relative rotations are $H_{ij}=Q_{i}Q_{j}^{\top}\!+\sigma N_{ij}$ ($i<j$), where $\sigma=0.3$ is the noise level and the $N_{ij}$’s are independent random noise matrices with i.i.d. normal entries. We also set $H_{ji}=H_{ij}^{\top}\!$ and $H_{ii}=I_{d}$. To estimate the $Q_{i}$’s from the $H_{ij}$’s, we set $C=-H/(nm)$ and solve (P). If the solution has rank $d$, this is equivalent to solving the maximum likelihood problem:

Remarkably, for all instances generated, (P) admits a rank $d$ solution, thus revealing the true maximum likelihood estimator: a hard quantity to compute, in general. This serendipitous phenomenon is partly explained in .

Figure 2 shows how much time it takes various solvers to find this solution of rank $d$ (they all do). Algorithm 1 runs RTR once on ($\textrm{RP}_{p}$) with $p=d+1$, with a random initial guess, and returns with an optimal rank $d$ solution. We compare against interior point methods SeDuMi , SDPT3 and Mosek (the latter two via CVX ) as well as against SDPLR with and without forcing the search rank to $d+1$ (the forced version is labeled SDPLR*). We also depict how much time it takes to simply compute the top $d$ eigenvectors of $H$, which, after projection, reveal an (empirically) equally good estimator for this problem, but with weaker guarantees .

Special case: Pseudo-Huber loss cost function

In the previous section, orthogonal synchronization is considered with Gaussian noise on the relative measurements. Maximum likelihood estimation then naturally leads to the minimization of a quadratic cost in $Y$, which simplifies to a linear cost in $X$.

Tools in this paper do not directly apply to the LUD cost, because it is nonsmooth.Recent work on nonsmooth optimization on manifolds may prove useful in this regard. Unfortunately, in our experiments we also found that smoothing the LUD cost typically leads to higher rank solutions, at a significant computational premium. We formalize this observation in the following theorem. The assumptions on $H$ are not restrictive: they require just the slightest inconsistency in the measurements.

In view of these results, we take interest in minimizing the related smoothed cost:

We may still compute a KKT point for (P), but there is no guarantee that such a point will be even a local minimizer anymore. On the bright side, Corollary 3.12 states that, by strong concavity of $f$, all KKT points of (P) are extreme points (thus they have rank at most $\lfloor p^{*}\rfloor$ (16)) and, for $p>p^{*}$, all second-order critical points of ($\textrm{RP}_{p}$) reveal KKT points of (P). The numerical experiment below shows that, empirically, even for $\varepsilon>0$, the proposed algorithm typically converges to a rank-$d$ KKT point of excellent quality. Furthermore, as $\varepsilon$ is decreased, the quality of the found KKT point increases (with warm-starting).

We now use the proposed robust formulation of orthogonal synchronization to situations where the sought matrices are in fact permutationsPermutation matrices are binary matrices with exactly one 1 on each row and column. They are orthogonal. (without modifying the algorithms). Synchronization of permutations notably arises in image association problems in computer vision .

Let the $Q_{i}$’s be permutations to estimate and let the $H_{ij}$’s be measurements of the relative permutations $Q_{i}Q_{j}^{\top}\!$. A subset of the measurements of a given size is selected uniformly at random and replaced by uniformly random permutations (outliers). The other measurements are correct. If perfect recovery of the $Q_{i}$’s is achieved, then the permutations are recovered.

Figure 3 exhibits the perfect recovery phenomenon hinted by Wang and Singer . We say “hinted” as the chosen scenario does not exactly fit the assumptions of these authors. Even in the face of many outliers, the true permutations are recovered, showing the applicability of the proposed methods to permutation estimation.

Conclusions and perspectives

We proposed a novel algorithm to compute KKT points for optimization problems over a class of spectrahedra that come up in relaxations of various problems involving orthonormal matrices. Our approach consists in exploiting the smooth geometry of bounded-rank subsets of those spectrahedra, to reduce the problem to Riemannian optimization. This effectively allows one to control how much lifting (dimension increase) is involved in the relaxation. An investigation of both the convex and the Riemannian geometries of the total and the bounded-rank problem showed that, under certain conditions, it is only necessary to compute second-order critical points on a low-dimensional portion of the boundary of the spectrahedron. Numerical experiments confirm the usefulness of this observation.

The present work triggers a number of questions for future investigation.

Which spectrahedra are such that their elements of bounded rank form a smooth manifold? (Journée et al. cover a number of such sets.) When the search space is of such form with additional constraints, can those be accommodated efficiently? This would be useful to address the SDP’s in, e.g., .

What is the computational complexity of obtaining a second-order critical point of a sufficiently smooth function on a Riemannian manifold, up to a given accuracy? This might be answered by following work in . When $Y$ is only approximately second-order critical and rank deficient, is $YY^{\top}\!$ approximately KKT, in a certain sense? For linear $f$, Proposition 5.1 offers a positive answer.

For nonconvex $f$, the set of KKT points includes the local optimizers of (P), as well as a number of uninteresting points. All KKT points give rise to critical points, but not necessarily second-order critical points. Can this be used to improve guarantees? Could we compute second-order KKT points instead, thus possibly excluding even more spurious points? A starting point might be [57, Thm. 3.45] and .

Assuming linear $f$, if (P) admits a unique solution of rank $r$ (see, e.g., ), is it sufficient to explore ($\textrm{RP}_{p}$) with $p=r+1$? Under the noise model of Section 5, this is observed empirically. Perhaps, this could be investigated via the expected size of the attraction basin of the global optimizers, similarly to in the context of dictionary learning.

Finally, regarding Corollary 3.14 and its attached question: for a random linear cost function and $p>p^{*}$, what is the probability that ($\textrm{RP}_{p}$) admits second-order critical points which are not global optimizers?

The author thanks P.-A. Absil, A. d’Aspremont, A. Bandeira, X. Cheng, B. Gerencsér, Y. Khoo, B. Mishra, A. Singer and B. Vandereycken for fruitful discussions. Parts of this research were conducted while N.B. was a research fellow with the FNRS in Belgium, and while generously supported by a Research in Paris grant, by the “Fonds Spéciaux de Recherche” (FSR) at UCLouvain and by the Chaire Havas “Chaire Economie et gestion des nouvelles données”, the ERC Starting Grant SIPA and a Research in Paris grant in France.

References

Appendix A Tightness for synchronization on a cycle

with the indexing convention that $Y_{m+1}\equiv Y_{1}$ and $H_{m,m+1}\equiv H_{m,1}$. Under a permissive condition on the measurements, Sharp et al. , Peters et al. exhibit an explicit formula for the solution (they restrict their attention to rotation matrices, that is, orthogonal matrices of determinant +1). We show that under that same condition, the corresponding SDP relaxation (P) with

is tight: there exists a unique solution of rank $d$ which reveals the global optimum.

The proof rests on two key ingredients: (a) we have an explicit formula for the solution $X$ to certify, and (b) we have an explicit formula for the dual certificate $S(X)$ to check. It seems reasonable to expect that the result should carry over to connected graphs whose cycles have disjoint edges.

Recently, Zhang and Singer gave a powerful tightness result for such SDP’s: as stated, their result is restricted to cycles of length 3, but it nicely accommodates non-orthogonal blocks in the data matrix $C$.

Part 1: guessing $X$. When the measurements are perfectly consistent, $P=I_{d}$, and it is easy to construct $X$: set $Y_{m}=I_{d}$ and $Y_{i}=H_{i,i+1}Y_{i+1}$ for $i=1\ldots m-1$; then $X_{ij}=Y_{i}Y_{j}^{\top}\!$. This construction does not use $H_{m,1}$ but still achieves $g(Y)=0$ owing to $P=I_{d}$, hence $X$ is optimal. When the cycle is inconsistent, it is reasonable to guess that the least-squares criterion will attempt to spread the inconsistency evenly over each edge . One $m$th of the error is represented by an offset $P^{1/m}$—taking this principal matrix root requires $P$ not to have negative eigenvalues. We build $Y$ by incorporating part of the error at each step, appropriately aligned. First define this recurrence: $Q_{m}=H_{m,1}$ and $Q_{i}=H_{i,i+1}Q_{i+1}$ for $i=(m-1)\ldots 1$ (note that $Q_{1}=P$). Then, $Y_{m}=I_{d}$ and $Y_{i}=Q_{i}P^{-1/m}Q_{i}^{\top}\!H_{i,i+1}Y_{i+1}$ for $i=(m-1)\ldots 1$. As previously, $X_{ij}=Y_{i}Y_{j}^{\top}\!$. It is not hard to check that $X_{ij}=Q_{i}P^{(i-j)/m}Q_{j}^{\top}\!$. Of course, $X$ is admissible for (P) with rank $d$.

All blocks of (32) are diagonal, so that its rows and columns may be permuted (without affecting its spectrum) to make it block diagonal, with $k$th block of size $m$ given by

In particular, $\lambda_{2}(A),\ldots,\lambda_{m}(A)>0$. Since the vector $[e^{i(1/m)\theta_{k}},e^{i(2/m)\theta_{k}},\ldots,e^{i(m/m)\theta_{k}}]^{*}$ is in the kernel of $A$, it must be that $\lambda_{1}(A)=0$. This concludes the proof. ∎

In general, the condition on the eigenvalues of $P$ is necessary. Indeed, for $d=1$ and $m=3$, choose the measurements such that $P=-1$ (for example, $+1$, $+1$ and $-1$) and verify that none of the 4 admissible rank-1 matrices are optimal. For the frequent case where the measurements are rotation matrices (that is, orthogonal with determinant +1), the condition is not too restrictive: even if they were distributed uniformly at random, $P$ would satisfy the condition almost surely.

Appendix B Generic face dimension

For $d=1$, the following theorem shows almost all faces of $\mathcal{C}$ have minimal dimension, as per the bound (15) . Key parts of the proof are due to Xiuyuan Cheng and Balázs Gerencsér.

We expect that this result remains valid for $d\geq 1$, but we are missing a proof.