Multireference Alignment using Semidefinite Programming

Introduction

This problem has a vast list of applications. Alignment is directly used in structural biology [Dia92] [TS12]; radar [ZvdHGG03] [PZAF05]; crystalline simulations [SSK13]; and image registration in a number of important contexts, such as in geology, medicine, and paleontology [DM98] [FZB02]. Various methods to solve this problem are used in these communities (see Appendix A).

A naïve approach to estimate the shifts in (1) would be to fix one of the observations, say $y_{i}$, as a reference template and align every other $y_{j}$ with it by the shift $\rho_{ij}$ minimizing their distance

This solution works well at a high signal-to-noise ratio (SNR), but performs poorly at low SNR. A more democratic approach would be to calculate all of the pairwise relative shift estimates $\rho_{ij}$ before attempting to recover the shifts $\{l_{i}\}$. This can be done by solving the minimization problem

This problem is known as angular synchronization [Sin11, BSS12] and the solution can be approximated via a SDP-based relaxation.

The issue with attempting to solve the alignment problem using either (2) or (3) is that one is evaluating the performance of a given choice of $l_{i},l_{j}$ for a pair $(i,j)$ by how far $l_{i}-l_{j}$ is from $\rho_{ij}$, but not taking into account the cost associated with other possible relative shifts. Relating $R_{-l_{i}}y_{i}$ and $R_{-l_{i}}y_{j}$, for example, would take into account information about all possible shifts instead of just the best one. The quasi maximum likelihood estimator (Section 2) attempts to do exactly that by solving the minimization problem:

Finding the MLE (4) is a non-trivial computational task because the parameter space is of exponential size, and the likelihood function is non-convex. While one can apply optimization methods such as gradient descent, simulated annealing, and expectation-maximization (EM), these are only guaranteed to find local minima of (4), but not the global minimum.

In this paper we take a different approach and propose a semidefinite relaxation for the quasi maximum likelihood problem (4). This particular SDP is inspired by an approximation algorithm designed to solve the Unique Games problem [CMM06] (Section 3).

Convex relaxations of hard combinatorial problems have seen many successes in applied mathematics. They became particularly popular in the last decade with the introduction of Compressed Sensing, in the seminal work of Donoho, Candes, Tao, and others [CRT06, Don06]. This idea has since been applied to a vast list of problems. Semidefinite programming (SDP) has served as a convex surrogate for problems arising in applications such as low-rank matrix completion [CR09], phase retrieval [CSV11], Robust PCA [LMTZ12], multiple-input multiple-output (MIMO) channel detection [MCS10], and many others. In many of these applications the same phenomenon is present: for typical instances, solving the convex problem is often equivalent to solving the original combinatorial problem [ALMT13].

Convex relaxations (and, in particular SDP based relaxations) also play a central role in the design of approximation algorithms in theoretical computer science. Almost two decades ago, Goemans and Williamson [GW95] proposed a SDP based approximation algorithm for the MAX-CUT problem with approximation ratio $\alpha^{\textsf{GW}}\approx 0.878$. That is, for any instance of the problem, the computed solution is guaranteed to provide performance (in this case, a cut) at least $\alpha^{\textsf{GW}}$ of the optimum. Many semidefinite relaxations have since been proposed as approximation algorithms for a long list of NP-hard problems [WS11].

In order to better understand the theoretical limitations of approximation algorithms, substantial work has been done to establish limits on the approximation ratios achievable by poly-time algorithms for certain NP-hard problems (hardness of approximation). The Unique Games conjecture by Khot [Kho02] is central to many recent developments: For $\delta,\varepsilon>0$, it is impossible for a polynomial-time algorithm to distinguish between $\delta$-satisfiable and $(1-\varepsilon)$-satisfiable Unique-Games instances. A Unique-Games instance consists of a graph along with a permutations for each edge. The problem is to choose the best assignment of labels to each vertex such that as many of the edge permutations are satisfied. The validity of the UGC would imply the optimality of certain poly-time approximation ratios, in particular the Goemans-Williamson constant $\alpha_{\textsf{GW}}$ for the MAX-CUT problem [GW95, KKMO07].

The best known polynomial time approximation to the unique games problem [CMM06] is based on an SDP relaxation of a formulation that uses indicator variables. It is quite different from SDPs normally used in applications (such as those described above). In particular, the variable matrix has size $NL\times NL$ and $\Omega\left(N^{2}L^{2}\right)$ constraints. We adapt this SDP to approximate the quasi maximum likelihood problem (4).

As we show that it is Unique-Games hard to approximate (4) within any constant (see Section 2) it is hopeless to aim for good guarantees for general instances. However worst case analysis is often too pessimistic and not indicative of performance observed in practice. In fact, under the random noise model we have for the observations, numerical simulations suggest that the SDP relaxation performs remarkably well, seeming to outperform existing methods. In an attempt to explain this phenomenon we show that our SDP is stable at high SNR levels and that it is tight at extremely high SNR levels. By stability, we mean that with high probability the solution to the SDP does not deviate much from the true solution (see Theorem 4.2). The stability for this SDP is particularly interesting as it is more challenging to analyze than random instances of Unique-Games [AKK+08, KT07], since our noise model is on the vertices, which translates into dependent noise on the edges. Still, these results fall short of properly explaining the remarkable performance that we see in simulations and more research is needed towards understanding the typical behavior of this SDP.

In order to simplify the SDP we also study a version with fewer constraints. Interestingly, this weaker SDP can be solved explicitly and is equivalent to the pairwise alignment method called phase correlation [HG84]. This method does not take into account information between all pairs of measurements, which suggests that the full complexity of the Unique-Games SDP [CMM06] is needed to obtain a good approximation to (4).

The fact that a global shift does not affect the solution to (4) creates symmetries in our SDP relaxation. In fact, we leverage such structure by using symmetry reduction techniques from matrix representation theory to simplify the analysis and computation of the SDP, greatly decreasing its computational cost. This is quite useful given the high computational cost of semidefinite programming.

Contributions: Our main contribution is applying techniques from theoretical Computer Science to a problem in applied Math. We introduce an Unique Games style SDP relaxation for the alignment problem that is novel for the applied Math community. From the theoretical Computer Science point of view, we introduce a new problem that has a similar flavor to the Unique Games problem – in fact we show that the worst case version is at least as hard as Unique Games. We introduce a natural average case version of this alignment problem - aligning several shifted copies of a signal corrupted by independent Gaussian noise. Existing analyses of semi-random models of Unique Games do not seem to apply to this problem. We show that for sufficiently high SNR, the SDP solution is close to an integer solution – this is a first step to establishing a signal recovery result which we leave as an open problem. We believe that future investigations into this problem will yield interesting insights into the Unique Games SDP and on dealing with correlated noise in average case analysis.

Quasi Maximum Likelihood Estimator

The log likelihood function for the model (1) is

Maximizing $\mathcal{L}$ is equivalent to minimizing the sum of squared residuals $\sum\|R_{-l_{i}}y_{i}-x\|^{2}.$ Fixing the $l_{i}$’s, the minimal value of $\mathcal{L}$ occurs at the average $x=\frac{1}{N}\sum_{i=1}^{N}R_{-l_{i}}y_{i}$. Making the tame assumption that $\|x\|^{2}$ is estimable (the norm is shift-invariant), maximizing (5) is thus equivalent to maximizing the sum of the inner products $\left\langle R_{-l_{i}}y_{i},\,R_{-l_{j}}y_{j}\right\rangle$ across all pairs $(i,j)$. Thus we consider the estimator

Unfortunately, the search space for this optimization problem has exponential size. Indeed, assuming no model for the vectors $\{y_{i}\}$, it is NP-hard to find the shifts which maximize (6), or even estimate it within a close constant factor.

It is NP-hard (under randomized reductions) to find a set of labels approximating (6) within $16/17+\varepsilon$ of its optimum value. It is UG-hard (under randomized reductions) to approximate (6) within any constant factor.

Proof. (outline) We give a randomized reduction from the class of $\Gamma$-MAX-2LIN($q$) instances consisting of a set of $2$ variable linear equations of the form $x_{i}-x_{j}\equiv c_{ij}\pmod{q}$, with the goal of choosing an assignment for the variables which maximizes the number of satisfied equations. We construct a vector $y_{k}$ for every variable $x_{k}$ such that shifts of $y_{k}$ correspond to an assignment to $x_{k}$. We pick a random vector $z_{ij}$ corresponding to a constraint on variables $x_{i},x_{j}$ and place a copy of $z_{ij}$ at specific locations in $y_{i}$ and $y_{j}$. Shifts of $y_{i},y_{j}$ corresponding to satisfying assignments of the constraint results in a superposition of the copies of $z_{ij}$. We choose parameters so that the only non-trivial contributions to the objective function (6) come from such superposition. The value of the objective is (within small error) a scaled version of the number of constraints of the $\Gamma$-MAX-2LIN($q$) instance satisfied by the assignment corresponding to the shifts. Thus hardness results for $\Gamma$-MAX-2LIN($q$) directly translate to hardness results for the alignment problem. The details are given in Appendix C.

The discrete optimization problem (6) may be formulated using indicator variables as an integer programming problem

The data Gram matrix $C$ with entries $C_{ik;jl}=\left\langle R_{-k}y_{i},R_{-l}y_{j}\right\rangle$ satisfies:

$C\succeq 0$ and has rank $L$, with non-zero eigenvalues $\lambda_{k}=L\sum_{i=1}^{N}\left|\mathcal{F}(Y_{i},k)\right|^{2}.$

Unfortunately, this spectral gap will not be apparent for a large class of signals. As long as a single power spectra $|\mathcal{F}(x,k)|^{2}$ of $x$ is near zero, the corresponding eigenvalue $\lambda_{k}$ will separate less from the small eigenvalues of $C$, and hence the space of the top $L$ eigenvectors of $C$ will contain less relevant information. For this reason, it would only be meaningful for us to characterize the SNR by the spectral gap $\min_{k}|\mathcal{F}(x,k)|^{2}$. Furthermore, our simulations suggest that recovery from a spectral relaxation performs worse than the semidefinite relaxation we are about to propose.

Semidefinite relaxation

This formulation attempts to count the number of satisfied edge constraints for an Unique-Games instance. One can treat the SDP (8) as an instance of $\Gamma$-MAX-2LIN($L$) on a weighted complete graph, with each cyclic permutation weighted by $C_{ik;jl}=\langle R_{-l_{i}}y_{i},R_{-l_{j}}y_{j}\rangle$. In this context, the matrix $C$ is dubbed the label-extended adjacency matrix [Kol10]. Thus the significant body of literature conducted on Unique-Games may be useful in understanding (8). Another common feature of the aligment SDP with $\Gamma$-MAX-2LIN($L$) instances is that the assigned labels may be chosen upto cyclic symmetry. This induces a block circulant symmetry in the semidefinite program (see Lemma D.5). For example, this symmetry will allow us to presume that $V_{ik;ik}=1/L$ from the constraint $\sum_{k}V_{ik;ik}=1$.

A major feature of (8) which does not feature in the study of Unique-Games is the structure of the data coefficient matrix $C$. While Unique-Games specifies constraints on edges of a graph (there are $N^{2}$ pieces of information), the alignment problem only specifies information on its vertices ($N$ pieces of information). This does assist our understanding of the semidefinite program, since it enables us to apply more symmetry conditions, but it also significantly complicates our analysis.

which is the concatenation of phase correlation vectors (see Appendix A) between $y_{1}$ and $y_{i}$.

are respectively equivalent to the Fourier side constraints

Since $\mathcal{V}_{k}\succeq 0$, the absolute value of each entry is bounded by the maximum of its diagonal entries, so each entry has magnitude at most $1/L$. Hence

with equality occuring when $\mathcal{V}_{k}$ has the entries of $\mathcal{C}_{k}$, normalized to magnitude $1/L$. Then $(\mathcal{V}_{k})_{ij}=\frac{\mathcal{F}(y_{i},l)\mathcal{F}^{*}(y_{j},l)}{|\mathcal{F}(y_{i},l)\mathcal{F}^{*}(y_{j},l)|}$, and a basis of the non-trivial eigenspace of $V$ is given by

From the SDP solution, one must round back to a solution in the original search space. There is a significant body of literature on the topic of rounding the solutions to various SDPs for Unique Games. The analysis and guarantees for these rounding schemes are in terms of the number of constraints satisfied by the solution produced and do not immediately give a result about signal recovery in our setting. In their study of semi-random instances of Unique-Games, the authors of [KMM11] give a rounding technique that uses both SDP and LP solutions when the SDP is known to be somewhat sparse – this is similar to a condition we obtain in the following stability section. Exploiting these ideas to establish an exact signal recovery guarantee is an interesting open problem.

Stability

$D_{ij}$ is a measure of how much the SDP would weight shift preferences other than the ground truth. When all $D_{ij}=0$ we obtain the integral instance.

For convenience, we make the definitions $\xi_{0}=2x$ and $\eta_{ij}=2\max_{l}|\langle R_{l}\xi_{i},\xi_{j}\rangle|\geq 0$ for $i,j=0,1,\ldots,N$. The following lemma demonstrates that the $\eta_{ij}$’s can be treated as worst-case bounds on the noise terms of the entries of the data Gram matrix $C$.

If $\operatorname{tr}(CV)\geq\operatorname{tr}(CV^{int})$, then $\sum_{i\neq j}(\Delta-\eta_{ij}-\eta_{j0})D_{ij}\leq 0$.

Proof. We can find a block circulant matrix which attains the same SDP objective value as $V$, so without loss of generality presume $V$ is block circulant. Hence $\sum_{k}V_{ik;j0}=1/L$ and $D_{ij}=1-LV_{i0;j0}$. Expanding,

The second inequality requires the SDP constraint $V\geq 0$. The pessimistic bound

and rearrangement gives the desired inequality.

With probability $1-e^{-N+o(N)}$, the solution to the SDP satisfies

Proof. For sufficiently high SNR, we would expect that the inequality in Lemma 4.1 would fail to hold. Indeed, by Lemma E.2, the inequality

holds with probability at least $1-e^{-N+\mathcal{O}(\log N)}$. It arises from tail bounds on the sum of slightly dependent random variables, and is independent of the structure of the SDP (as opposed to Lemma 4.1). Combined with Lemma 4.1, we can obtain a guarantee on the deviation between the SDP solution and the integral instance. The full proof may be found following Theorem E.3.

This result may be strengthened in other manners. For constants $0<\delta,\varepsilon\ll 1$, we can attain a tighter concentration condition of the form

instead of the current condition $\sum D_{ij}\geq N^{2}/SNR$. The argument is based on the analysis of the SDP for adversarial semi-random Unique-Games instances by [KMM11]. However, the proof must be more nuanced in our case, due to correlations in the noise model of $C$, caused by the smaller source of randomness available to us. Hence we omit the full details, but provide a sketch below.

Instead of finding a global tail bound in Lemma E.2, we can derive a local tail bound for each $\alpha_{JL}$-ball of $D_{ij}^{\varphi}$’s, each tail bound holding with probability at least $1-2Ne^{-t^{2}N}$. For $0<\delta,\,\varepsilon=\varepsilon^{\prime}\ll 1$ constants, the number of $N$-tuples of vectors in $\mathfrak{N}$ is of size $\exp(\mathcal{O}(N))$. With a sufficiently large constant $t$, the local tail bound may be union bounded across all balls of vectors in $\mathfrak{N}^{N}$, and thus will hold for all SDP-feasible $V$. With some care, this argument would yield a concentration condition of the form (14).

Numerical Results

We implemented several baseline methods discussed thus far, and plotted their average error performance across $500$ iterations in Figure 1. For each iteration, we chose a signal $x$ randomly from the distribution $\mathcal{N}(0,I_{L})$, as well as $N$ i.i.d noise vectors $\xi_{i}\sim\mathcal{N}(0,\sigma^{2}I_{L})$, and apply each of our methods. These simulations confirm our intuition that the UG-based SDP performs better than other benchmark methods. In particular, they suggest that the UG-based SDP is highly stable around integral instances.

The implementation using bispectrum-like invariants is discussed in Appendix A. For each of the other procedures, we construct a $NL\times L$ matrix $W$ recording alignment preferences. For cross- or phase correlation, let $W_{ik;l}$ be the $(k-l)$th entry of the cross- or phase correlation vector between $y_{1}$ and $y_{i}$. For the spectral rounding off the Gram matrix $C$, and the solution of the semidefinite program $V$, let $W$ be the top $L$ eigenvectors of the respective matrix. The shifts are read off this matrix, and the un-shifted $y_{i}$ are averaged to produce an estimate for $x$. The first plot (a) shows the difference between this estimate and $x$.

Generalizations and Future Work

It would be interesting to see how the performance of the SDP changes when we study the alignment problem across more difficult groups, especially non-abelian ones, which arise in several applications.

The numerical simulations in Section 5 suggest that the UG-based SDP achieves exact recovery with high probability for sufficiently high SNR. That is, the resulting SDP matrix is integral, so by solving the SDP we are indeed obtaining the solution to the quasi-ML estimator. Indeed, as the SNR decreases, there appears to be a phase transition during which the SDP almost always recovers an integral solution. Our analysis of the stability of the semidefinite program does not fully explain this phenomenon. The authors believe this to be an interesting direction for future work, especially since guarantees of exact recovery are attainable in high SNR settings for a few semidefinite relaxations, for example for the MIMO problem [MCS10].

Another important question is to understand the sample complexity of our approach to the alignment problem. Since the objective is to recover the underlying signal $x$, having a larger number of observations $N$ should give us better recovery. The question can be formalized as: for a given value of $L$ and $\sigma$, how large does $N$ need to be in order to have reasonably accurate recovery? The sample complexity of methods like the bispectral invariants would be expected to require $N=\Omega(\sigma^{2}L^{2}\log L)$ observations. We would hypothesize on the strength of our numerical results that the UG-based SDP requires fewer observations for meaningful recovery, and establishing this is an interesting open problem.

Along with expanding the domain of the alignment problem, it would be interesting to attempt the style of analysis discussed in this paper for other maximum likelihood problems. Maximum likelihood estimators play an important role in many estimation problems, but often (as in our problem) computing or approximating the MLE is a challenging problem and semidefinite programming could provide a tractable alternative in an average case setting.

Acknowledgements

The authors thank Yutong Chen for valuable assistance with the implementation of our algorithm. A. S. Bandeira was supported by AFOSR Grant No. FA9550-12-1-0317. M. Charikar was supported by NSF grants CCF 0832797, AF 0916218 and AF 1218687. A. Singer was partially supported by Award Number FA9550-12-1-0317 and FA9550-13-1-0076 from AFOSR, by Award Number R01GM090200 from the NIGMS, and by Award Number LTR DTD 06-05-2012 from the Simons Foundation. Parts of this work have appeared in A. Zhu’s undergraduate thesis at Princeton University.

References

Appendix A Phase correlation and the Bispectrum

Many scientific fields have isolated discussions and solutions for the alignment problem, occasionally with slight context-specific adjustments. Such methods include iterative template aligning [KDM95], zero phase representations [ZvdHGG03], angular synchronization [SSK13], and machine learning [PZAF05]. Unfortunately, most of these are either do not fairly weight each observation, do not make full use of the available information, or do not have rigorous performance discussion. We outline just a couple of very well-studied alignment methods in this section.

In the case of a pair of noisily shifted vectors ($N=2$), several long-studied methods assign scores to each possible shift between the vectors, and estimate the shift as the one with the highest score. The most natural is to take the inner product between each possible shift of the vectors $y_{i},y_{j}$, which gives the cross-correlation vector $v_{l}^{\textsf{cross}}=\langle y_{i},R_{-l}y_{j}\rangle$, with maximal entry $\widehat{l}=\operatorname*{argmax}v_{l}^{\textsf{cross}}$. Another frequently used in practice is the phase correlation vector

For human-friendly images, phase correlation tends to be more appealing, since it tends to have a single sharp peak. The relation between cross- and phase correlation can be seen from the convolution theorem:

Convolution theorem: $\left\{\langle y_{i},R_{-l}y_{j})\rangle\right\}_{l=0}^{L-1}=\sqrt{L}\cdot\mathcal{F}\left\{\mathcal{F}(y_{i},k)\mathcal{F}^{*}(y_{j},k)\right\}_{k=0}^{L-1}.$

Proof. Let $M_{l}$ denote the modulation operator $(M_{l}x)_{k}=x_{k}e(kl/L)$, defined such that $\mathcal{F}R_{l}x=M_{l}\mathcal{F}x$. By Parseval’s Theorem,

Cross- and phase correlation can be used directly for the problem of multireference alignment, say by aligning all of the observations against the first one. However, this is certainly not a robust method.

its shift invariance can be seen since it is the Fourier transform of the $d$-th order autocorrelation function

by the Wiener-Khinchin theorem. How much information do these invariants capture about the signal $x$? Sadler and Giannakis give iterative and least squares algorithms in [SG92] for reconstructing the Fourier phases from full knowledge of the bispectrum, and Kakarala in [KI93] show the uniqueness of real function given all of its higher order moment spectra. These arguments generally tend to be of a more number-theoretic flavor and may be tedious in practice.

A simple special case occurs when we take $k_{1}=k_{2}=\ldots=k_{d}=1$ for $d=1,\ldots,L$. This gives us the set of bispectral invariants $\mathcal{F}(x,d)\mathcal{F}^{*}(x,1)^{d}$. This quantity, as a shift invariant, can be consistently estimated from our observations $y_{i}$. With an estimate for the phase of the first Fourier coefficient $\mathcal{F}(x,1)$ (for example we can sample over a discretization of the unit circle), we can recover the remaining Fourier phases $\mathcal{F}(x,d)$ and thus the signal $x$.

Appendix B Maximum Likelihood

The data Gram matrix $C$ with entries $C_{ik;jl}=\left\langle R_{-k}y_{i},R_{-l}y_{j}\right\rangle$ satisfies the following properties:

$C\succeq 0$ and has rank $L$, with non-zero eigenvalues $\lambda_{k}=L\sum_{i=1}^{N}\left|\mathcal{F}(y_{i},k)\right|^{2}.$

so $C\succeq 0$ has rank $L$. Since $\langle R_{-k}y_{i},R_{-l}y_{j}\rangle=\langle R_{-k+m}y_{i},R_{-l+m}y_{j}\rangle$, $C$ is composed of $N\times N$ circulant blocks of size $L\times L$. After permuting the rows and columns of $C$, one may write $C$ as a block circulant matrix

is block diagonal. These block diagonals are componentwise Discrete Fourier Transforms of the “vector” $\left(C_{0},C_{1},\ldots,C_{L-1}\right)$:

by Lemma A.1. Note also that each $\mathcal{C}_{k}$ is Hermitian and rank $1$ ($C$ is of rank $L$, and each of the blocks $\mathcal{C}_{k}$ has positive rank). The unique nonzero eigenvalue $\lambda_{k}$ of $\mathcal{C}_{k}$ is given by

Appendix C NP-hardness of Quasi MLE

We demonstrate this by a poly-time approximation preserving reduction from a special class of MAX-2LIN$(q)$ instances called $\Gamma$-MAX-2LIN$(q)$. Consider a (connected) MAX-2LIN$(q)$ instance where each constraint has the form $x_{i}-x_{j}\equiv c_{ij}\pmod{q}$, of which at most $\rho^{*}$ are satisfiable. Representing each variable as a vertex of a graph and each constraint as an edge, the instance is associated with a connected graph $G=(V(G),E(G)),$ where $V(G)=[N]$ and $|E(G)|=M$.

Let $\textsf{poly}(M)$ be the space of integer functions which are bounded by polynomial order, i.e. $f\in\textsf{poly}(M)$ iff there are constants $C,k$ such that $f(M)\leq CM^{k}$ for all $M\geq 1$. We say that an event occurs w.h.p if it occurs with probability $1-\epsilon(q,G)$, where $\frac{1}{\epsilon(q,G)}\notin\textsf{poly}(q\cdot|E(G)|)$. Notice under this definition that if $\textsf{poly}(qM)$ events occur w.h.p, then by an union bound the event that all occur will also be w.h.p.

For each edge constraint $x_{i}-x_{j}\equiv c_{ij}$, let $z_{ij}$ be a vector uniformly at random chosen from $\{\pm 1\}^{s}$. Assign the entries $(0,\{i,j\},\cdot)$ of $y_{i}$ to $z_{ij}$, and the entries $(c_{ij},\{i,j\},\cdot)$ of $y_{j}$ to $z_{ij}$. Assign all of the remaining entries of the $y_{i}$’s to i.i.d Rademacher random variables ($\{\pm 1\}$ with probability $1/2$). Intuitively, a relative shift of $c_{ij}\cdot qM$ between $y_{i}$ and $y_{j}$ will produce a large inner product due to the overlapping of the $z_{ij}$’s, while any other shift between them would produce low inner products.

Proof. By independence, each inner product is the sum of $\gamma$ independent Bernoulli random variables which take values $\pm 1$ with probability $1/2$. Hoeffding’s inequality indicates

Assuming the Unique Games conjecture, it is NP-hard to approximate $\Gamma$-MAX-2LIN$(q)$ to any constant factor [KKMO07]. Hence it is UG-hard (under randomized reductions) to approximate ALIGN within any constant factor. MAX-CUT is a special case of $\Gamma$-MAX-2LIN$(q)$. Since it is NP-hard to approximate MAX-CUT within $\frac{16}{17}+\varepsilon$, it is NP-hard (under randomized reductions) to approximate ALIGN within $\frac{16}{17}+\varepsilon$.

Appendix D ∗*-matrix algebras

In this appendix section, we briefly sketch the ideas behind symmetry reduction in semidefinite programs. Most of the details are deferred to other sources. For example, [DdK11] is a good reference on this topic. While the application of symmetry to the SDPs we are interested in is relatively straightforward and can be explained without reference to $*$-matrix algebras, the notation and ideas behind the theory of $*$-matrix algebras make it especially clear to express.

A coherent configuration is a set of matrices $\{A_{1},\ldots,A_{d}\}\subset\{0,1\}^{n\times n}$ satisfying:

Some subset of the $A_{i}$’s sum to $I_{n}$.

$A_{i}A_{j}$ is a linear combination of $A_{1},\ldots,A_{d}$.

The span of the $A_{i}$’s is a $*$-matrix algebra $\mathcal{A}$. If the $A_{i}$’s commute and contain the identity matrix $I_{n}$, then this configuration is called an association scheme, and $\mathcal{A}$ is also known as a Bose-Mesner algebra. Treating $\mathcal{A}$ as a matrix vector space, a coherent configuration forms a basis for $\mathcal{A}$.

A $*$-homomorphism $\phi$ is a linear mapping of $*$-matrix algebras $\phi:\mathcal{A}\to\mathcal{B}$ which additionally satisfies $\phi(A^{*})=\phi(A)^{*}$, $\phi(AB)=\phi(A)\phi(B)$, and $\phi(I_{\mathcal{A}})=I_{\mathcal{B}}$ if the identity matrix $I_{\mathcal{A}}\in\mathcal{A}$.

If $\phi$ is an injective $*$-homomorphism $\phi:\mathcal{A}\to\mathcal{B}$, then $A\in\mathcal{A}$ and $\phi(A)$ share the same set of eigenvalues (ignoring multiplicity), and $A\succeq 0\Longleftrightarrow\phi(A)\succeq 0$.

To leverage symmetry in semidefinite programs, it will be useful to take the images of the objective semidefinite matrix under dimension-reducing $*$-homomorphisms. Two such $*$-homomorphisms which may be systematically constructed are given by Artin-Wedderburn’s Theorem and regular $*$-representations [Mur07]. Consider a primal semidefinite problem of the form

there is an optimal solution $V$ to the primal SDP problem contained in $\mathcal{A}_{SDP}$

$\text{Re}(V)\in\mathcal{A}_{SDP}$ is also an optimal solution to the primal SDP problem.

Proof. Refer to ([DdK11], Corollary 2.5.2).

A special case of this is the result applies for the $*$-algebra of $G$-invariant matrices.

In the semidefinite program (22), suppose the constraint matrices $C_{0},C_{1},\ldots,C_{m}$ are $G$-invariant. Then there is a solution $V$ to the SDP which is also $G$-invariant.

is $G$-invariant, feasible, and has the same objective value as $V^{\prime}$. The averaging map $\mathcal{R}_{G}$ is known as the Reynolds operator.

Appendix E Stability

The lemma follows by taking $T=\log L$, since \eta_{ij}=2\sigma\|\xi_{i}\|\cdot\eta\big{(}\xi_{i}/\|\xi_{i}\|\big{)} conditional on $\xi_{i}$.

Let $t>0$. With probability at least $1-2Ne^{-t^{2}N}$,

Proof. Define the random variables $Z_{ij}=(\eta_{ij}+\eta_{j0})D_{ij}\leq\eta_{ij}+\eta_{j0}$. For fixed $i$, the $\eta_{ij}$’s are independent positive random variables. The one-sided tail bound for independent positive random variables by Maurer [Mau03] gives

Choose $t_{i}=t\sqrt{N}$. By union bound, with probability at least $1-2Ne^{-t^{2}N}$,

With probability $1-e^{-N+o(N)}$, the solution to the SDP satisfies

Proof. Let $V$ be a solution matrix to the SDP, so that $\operatorname{tr}(CV)\geq\operatorname{tr}(CV^{int})$. It follows from Lemma 4.1 that

The $\chi^{2}$-tail bound of Laurent-Massart [LM00] states

with probability at least $1-e^{-NL}$. Union bounding with the tail bound of Lemma E.2, and applying Lemma E.1,

with probability at least $1-e^{-N+o(N)}$.

Let $0<\delta,\varepsilon,\varepsilon^{\prime}\ll 1$ be small constants, and define $\alpha_{JL}(x)=(1+\delta)x+\varepsilon$. There is a set of unit vectors $\mathfrak{N}$ of size at most

such that for any set of unit vectors $\{v_{i}\}$, there is a map $\varphi:\{v_{i}\}\to\mathfrak{N}$ satisfying the inequality

for at least $1-\varepsilon^{\prime}$ fraction of the pairs $(i,j)\in[N]\times[N]$.

Proof. This lemma appears frequently in SDP literature, and in particular is used to analyze adversarial semi-random Unique-Games instances in [KMM11]. Notice that the size of the set $\mathfrak{N}$ is independent of the number of vectors in the set $\{v_{i}\}$.

Construct a $\varepsilon/32$-net $\mathfrak{N}$ of the unit hypersphere in a $\mathcal{O}(\delta^{2}\log(1/\varepsilon^{\prime}))$-dimensional space $\mathfrak{L}$. By the (strong version of the) Johnson-Lindenstrauss lemma, there is a randomized mapping $\varphi^{\prime}:\{v_{i}\}\to\mathfrak{L}$ satisfying

with probability at least $1-\varepsilon^{\prime}$. Define $\varphi(v_{i})$ to be the closest point to $\varphi^{\prime}(v_{i})$ in $\mathfrak{N}$, and observe that

for all $x>0$, so $\varphi$ satisfies the conditions of the lemma.