Faster and Simpler Width-Independent Parallel Algorithms for Positive Semidefinite Programming

Introduction

Semidefinite programming (SDP), alongside linear programming (LP), is an important tool in approximation algorithms, optimization, and discrete mathematics. In the context of approximation algorithms alone, it has emerged as a key technique which underlies a number of impressive results that substantially improve the approximation ratios. To solve a semidefinite program, algorithms from the linear programming literature such as Ellipsoid or interior-point algorithms can be applied to derive near exact solutions. But they are often costly. As a result, finding efficient approximations to such problems is a critical step in making them more practical.

From a parallel algorithms standpoint, both LPs and SDPs are P-complete even to approximate to any constant accuracy, suggesting that it is unlikely that they have a polylogarithmic depth algorithm. For linear programs, however, the special case of positive linear programs, first studied by Luby and Nisan , has an algorithm that finds a $(1+\varepsilon)$-approximate solution in $O(\textrm{poly}(\tfrac{1}{\varepsilon}\log{n}))$ iterations. This weaker approximation guarantee is still sufficient for approximation algorithms (e.g., solutions to vertex cover and set cover via randomized rounding), spurring interest in studying these problems in both sequential and parallel contexts (see, e.g., ).

The importance of problems such as MaxCut and Sparsest Cut has led to the identification and study of positive SDPs. Our work is motivated by a result by Jain and Yao that gave the first positive SDP algorithm whose work and depth are independent of the width parameter (commonly known as width-independent algorithms). As there are substantial differences in the analysis in this version compared to the conference version , we address the relation between our paper and other works in Section 1.1.

We present a simple algorithm that offers guarantees similar to but has less work-depth complexity. Each iteration of our algorithm involves only simple matrix operations and computing the trace of the product of a matrix exponential and a positive semidefinite matrix. The input consists of an accuracy parameter $\varepsilon>0$ and a positive semidefinite program (PSDP) in the following standard primal form:

where the matrices $\mathbf{C},\mathbf{A}_{1},\dots,\mathbf{A}_{n}$ are $m$-by-$m$ symmetric positive semidefinite matrices, $\bullet$ denotes the pointwise dot product between matrices (see Section 2), and the scalars $b_{1},\dots,b_{n}$ are nonnegative reals. This is a subclass of SDPs where the matrices and scalars are “positive” in their respective settings. We also assume, as is standard, that the SDP has strong duality. Our main result is as follows:

Given a primal positive SDP involving $m\times m$ matrices with $n$ constraints and an accuracy parameter $\varepsilon>0$, there is an algorithm approxPSDP that produces a $(1+\varepsilon)$-approximation in $O(\frac{1}{\varepsilon^{3}}\log^{3}n)$ iterations, where each iteration involves computing matrix sums and a special primitive that computes $\exp(\mathbf{\Phi})\bullet\mathbf{A}$ in the case when $\mathbf{\Phi}$ and $\mathbf{A}$ are both positive semidefinite.

The theorem quantifies the cost of our algorithm in terms of the number of iterations. The work and depth bounds implied by this theorem vary with the format of the input and how the matrix exponential is computed in each iteration. As we will discuss in Section 4, with input given in a suitable form, our algorithm runs in nearly-linear work and polylogarithmic depth.

The first definition of positive SDPs was due to Klein and Lu , who used it to characterize the MaxCut SDP. The MaxCut SDP can be viewed as a direct generalization of positive (packing) LPs. More recent work defines the notion of positive packing SDPs , which captures problems such as MaxCut, sparse PCA, and coloring; and the notion of covering SDPs , which captures the ARV relaxation of Sparsest Cut among others. Works this area tend to focus on developing fast sequential algorithms for finding a $(1+\varepsilon)$-approximation, leading to a series of sequential algorithms (e.g., ). The iteration count for these algorithms, however, depends on the so-called “width” parameter of the input program or some parameter of the spectrum of the input program. In some instances, the width parameter can be as large as $\Omega(n)$, making it a bottleneck in the depth of direct parallelization. This is even more so in the case of the Sparsest Cut SDP even though the problem has been labeled as a covering SDP .

Over the past few years, width-independent algorithms for positive SDPs have received much attention. Jain and Yao gave the first result in this direction: a polylog depth algorithm based on the first width independent linear programming algoirthm by Luby and Nisan . Their algorithm is based on updating the primal matrix. This leads to an intricate analysis based on carefully analyzing the eigenspaces of a particular matrix before and after each update. It takes $O(\frac{1}{\varepsilon^{13}}\log^{13}m\log n)$ iterations, each of which involves computing spectral decompositions using least $\Omega(m^{\omega})$ work.

Our algorithm follows a different approach. It updates the dual program in a way motivated by the width-independent positive LP algorithm by Young . Concurrently, Jain and Yao gave a similar algorithm for positive SDPs. Their algorithm solves a class of SDPs which contains both packing and diagonal covering constraints. Since matrix packing conditions between diagonal matrices are equivalent to point-wise conditions of the diagonal entries, these constraints are closer to a generalization of positive covering LP constraints. We believe that removing this restriction on diagonal packing matrices would greatly widen the class of problems included in this class of SDPs and discuss possibilities in this direction in Section 5.

The convergence analyses of both of these routines go through the matrix multiplicative weights update framework for solving semidefinite programs . However, a crucial algebraic piece is missing for adapting the iteration count bound from Young’s algorithm. Unlike scalar exponentials, matrix exponentials is not a monotonic function under standard notions of matrix ordering such as the Loewner partial order.

Most recently, Allen-Zhu et al. gave the first rigorous adaptation of Young’s algorithm to positive SDPs . Motivated by this result, as well as its precursor in positive LPs , we complete the analysis of our original algorithm. Specifically, we show that the dual generation scheme in these routines provide good bounds on the iteration count. This leads to a different analysis of Young’s algorithm that omits phases. Our modified analysis is for a simplified pseudocode of the algorithm from that removes these phases. However, the phase-based version can be analyzed similarly.

Compared to the algorithm by Allen-Zhu et al. , our analysis uses more elementary linear algebraic techniques, and is closer to multiplicative weights update schemes. We believe the dynamic bucketing method for obtaining better dependencies on $\varepsilon$ is also applicable to our analysis. Finally, as many recent results on faster positive LP / SDP algorithms take optimization based views, we believe our approach is also of independent interest.

2 Overview

We derive our algorithm by generalizing Young’s algoirthm . In place of the “soft max” function for bounding the maximum of a set of linear constraints, we use matrix exponential and the matrix multiplicative weights update (MMWU) mechanism. Moreover, besides standard operations on (sparse) matrix, the only other primitive needed is the matrix dot product $\exp(\mathbf{\Phi})\bullet\mathbf{A}$, where $\mathbf{\Phi}$ and $\mathbf{A}$ are positive semidefinite.

Intuitions. For intuition about packing SDPs and the matrix multiplicative weights update method for finding approximate solutions, a useful analogy of the decision problem is that of packing a (fractional) amount of ellipses into the unit ball. Figure 1 provides an example involving $3$ matrices (ellipses) in $2$ dimensions. Note that $\mathbf{A}_{1}$ and $\mathbf{A}_{2}$ are axis-aligned; their sum is also axis-aligned in this case. In fact, positive linear programs in the broader context corresponds exactly to the restriction of all ellipsoids being axis-aligned. In this setting, the algorithm of can be viewed as creating a penalty function by weighting the length of the axises using an exponential function. Then, ellipsoids with sufficiently small penalty subject to this function have their weights increased.

However, once we allow general ellipsoids such as $\mathbf{A}_{3}$, the resulting sum will no longer be axis-aligned. In this setting, a natural extension is to take the exponential of the semimajor axises of the resulting ellipsoid instead, leading to the matrix multiplicative weights scheme. Our analysis then focuses on showing that Young’s algorithm still has a width-independent iteration count in this setting.

Work and Depth. We now discuss the work and depth bounds of our algorithm. The main cost of each iteration of our algorithm comes from computing the dot product between a matrix exponential and a PSD matrix. Like in the sequential setting , we need to compute for each iteration the product $\mathbf{A}_{i}\bullet\exp(\mathbf{\Phi})$, where $\mathbf{\Phi}$ is some PSD matrix. The cost of our algorithm therefore depends on how the input is specified. When the input is given prefactored—that is, the $m$-by-$m$ matrices $\mathbf{A}_{i}$’s are given as $\mathbf{A}_{i}=\mathbf{Q}_{i}\mathbf{Q}_{i}^{\top}$ and the matrix $\mathbf{C}^{-1/2}$ is given, then Theorem 4.1 can be used to compute matrix exponential in $O(\frac{1}{\varepsilon^{3}}(m+q)\log n\log q\log(\nicefrac{{1}}{{\varepsilon}}))$ work and $O(\frac{1}{\varepsilon}\log n\log q\log(\nicefrac{{1}}{{\varepsilon}}))$ depth, where $q$ is the number of nonzero entries across $\mathbf{Q}_{i}$’s and $\mathbf{C}^{-1/2}$. This is because the matrix $\mathbf{\Phi}$ that we exponentiate has $\|{\mathbf{\Phi}}\|_{2}\leq O(\frac{1}{\varepsilon}\log n)$, as shown in Lemma 3.5. Therefore, as a corollary to the main theorem, we have the following cost bounds:

The algorithm approxPSDP has in $\widetilde{O}(\frac{1}{\epsilon^{6}}(n+m+q))$ work and $O(\frac{1}{\epsilon^{4}}\log^{O(1)}(n+m+q))$ depth.

If, however, the input program is not given in this form, we can add a preprocessing step that factors each $\mathbf{A}_{i}$ into $\mathbf{Q}_{i}\mathbf{Q}_{i}^{\top}$ since $\mathbf{A}_{i}$ is positive semidefinite. In general, this preprocessing requires at most $O(m^{4})$ work and $O(\log^{3}m)$ depth using standard parallel QR factorization . Furthermore, these matrices often have certain structure that makes them easier to factor. Similarly, we can factor and invert $\mathbf{C}$ with the same cost bound, and can do better if it also has specialized structure.

Background and Notation

where $\mathbf{v}_{1},\mathbf{v}_{2},\dots,\mathbf{v}_{m}$ are the eigenvectors of $\mathbf{A}$ with eigenvalues $\lambda_{1}\geq\dots\geq\lambda_{m}$ respectively. We will use $\lambda_{1}(\mathbf{A})$, $\lambda_{2}(\mathbf{A}),\dots,\lambda_{m}(\mathbf{A})$ to represent the eigenvalues of $\mathbf{A}$ in decreasing order and also use $\lambda_{\max}(\mathbf{A})$ to denote $\lambda_{1}(\mathbf{A})$. Notice that positive semidefiniteness induces a partial ordering on matrices. We write $\mathbf{A}\preccurlyeq\mathbf{B}$ if $\mathbf{B}-\mathbf{A}\succcurlyeq\mathbf{0}$.

The trace of a matrix $\mathbf{A}$, denoted $\mathop{\mathsf{Tr}}\left[\mathbf{A}\right]$, is the sum of the matrix’s diagonal entries: $\mathop{\mathsf{Tr}}\left[\mathbf{A}\right]=\sum_{i}A_{i,i}$. Alternatively, the trace of a matrix can be expressed as the sum of its eigenvalues, so $\mathop{\mathsf{Tr}}\left[\mathbf{A}\right]=\sum_{i}\lambda_{i}(\mathbf{A})$. Furthermore, we define

It follows that $\mathbf{A}$ is positive semidefinite if and only if $\mathbf{A}\bullet\mathbf{B}\geq 0$ for all PSD $\mathbf{B}$.

where, again, $\mathbf{v}_{i}$ is the eigenvector corresponding to the eigenvalue $\lambda_{i}$. It is not difficult to check that for $\exp(\mathbf{A})$, this definition coincides with $\exp(\mathbf{A})=\sum_{i\geq 0}\frac{1}{i!}\mathbf{A}^{i}$.

Our algorithm relies on a matrix multiplicative weights (MMW) algorithm, which can be summarized as follows. For a fixed $\varepsilon_{0}\leq\frac{1}{2}$ and $\mathbf{W}^{(1)}=\mathbf{I}$, we play a “game” a number of times, where in iteration $t=1,2,\dots$, the following steps are performed:

Produce a “probability” matrix $\mathbf{P}^{(t)}=\mathbf{W}^{(t)}/\mathop{\mathsf{Tr}}\left[\mathbf{W}^{(t)}\right]$;

Incur a gain matrix $\mathbf{M}^{(t)}$; and

Like in the standard setting of multiplicative weights algorithms, the gain matrix is chosen by an external party, possibly adversarially. In our algorithm, the gain matrix is chosen to reflect the step we make in the iteration. Arora and Kale shows that the MMW algorithm has the following guarantees (restated for our setting):

For $\varepsilon_{0}\leq\frac{1}{2}$, if $\mathbf{M}^{(t)}$’s are all PSD and $\mathbf{M}^{(t)}\preccurlyeq\mathbf{I}$, then after $T$ iterations,

2 Reduction to Bounded Decision Version

Our algorithm works with normalized primal/dual programs shown in Figure 2.

By using binary search and appropriately scaling the input program, such an SDP can be approximated using the following decision problem:

or a PSD matrix $\mathbf{Y}$ (a primal solution) such that

This reduction can also ensure bounded trace on all $\mathbf{A}_{i}$’s. The following lemma summarizes key properties of such a reduction:

For $0<\varepsilon<1$, a positive packing semidefinite program can be approximated to a relative error to $\varepsilon$ using $O(\log{n})$ calls to the $\varepsilon$-decision problem. Furthermore, each $\mathbf{A}_{i}$ supplied to the decision problem has $\mathop{\mathsf{Tr}}\left[\mathbf{A}_{i}\right]\leq O(n^{3})$.

The reduction is common to most positive linear and semidefinite program solvers ; we briefly sketch the idea for completeness.

(Sektch) First, transform to the normalized form by “dividing through” by $\mathbf{C}$ (see Appendix A). Since $\mathop{\mathsf{Tr}}\left[\sum_{i}x_{i}\mathbf{A}_{i}\right]$ is within a factor of $n$ of the maximum eigenvalue of $\sum_{i}x_{i}\mathbf{A}_{i}$, $\frac{1}{\min_{i}\mathop{\mathsf{Tr}}\left[\mathbf{A}_{i}\right]}$ gives a value that is within a factor $n$ of the optimum. Therefore, the optimization version can be solved by binary searching on the objective a total of at most $O(\log(\frac{n}{\varepsilon}))$ iterations.

In each of these decision instances, we can rescale $\mathbf{A}_{i}$’s so that the threshold in question is $1$. In this setting, the total sum of $x_{i}$ with $\mathop{\mathsf{Tr}}\left[\mathbf{A}_{i}\right]\geq n^{3}$ is at most $1/n$. By not using these $x_{i}$’s, the optimum changes by an additive value of less than $\varepsilon$.

Solving Positive SDPs

In this section, we describe a parallel algorithm for solving the decision version of positive packing SDPs, inspired by Young’s algorithm for positive LPs. The following theorem presents the guarantees of our algorithm.

Let $0<\varepsilon<1$. There is an algorithm decisionPSDP that given a positive SDP, solves the $\varepsilon$-decision problem in $O(\varepsilon^{-3}\log^{2}n)$ iterations.

Presented in Algorithm 3.1 is an algorithm that we will show to satisfy the theorem.

The algorithm is a multiplicative weights update routine, which proceeds in several rounds. The starting solution is $\mathbf{x}^{(0)}_{i}=\frac{1}{n\mathop{\mathsf{Tr}}\left[\mathbf{A}_{i}\right]}$. This solution is chosen to be small so that $\sum_{i}x_{i}^{(0)}\mathbf{A}_{i}\preccurlyeq\mathbf{I}$, hence respecting the dual constraint. Each subsequent update is a multiple of the current solution, so this $\mathbf{x}^{(0)}$ is also chosen to ensure that these updates make rapid progress.

In each iteration following that, the algorithm computes

Our presentation follows the multiplicative weights update framework from Arora-Kale and Kale . Several intermediate variables are helpful for further discussion. Define the cumulative sum corresponding to $\mathbf{x}^{(t)}$ as

This allows us to write the exponential as $\mathbf{W}^{(t)}=\exp(\mathbf{\Psi}^{(t)})$. We will also specifically define the probability matrix by which we use to pick the update coordinates:

This matrix is easier to work with because it has trace $1$:

Using this “probability” matrix $\mathbf{P}^{(t)}$, the algorithm identifies which $\mathbf{x}$ coordinates to update. These are the coordinates $x_{i}$’s for which their contributions with respect to the $\mathbf{A}_{i}$’s are still small—$\mathbf{P}^{(t)}\bullet\mathbf{A}_{i}\leq 1+\varepsilon$. Each of these $\mathbf{x}$ coordinates will be incremented by $\delta^{(t)}\stackrel{{\scriptstyle{\smash{\textsf{def}}}}}{{=}}\alpha\cdot\mathbf{x}_{i}$, where $\alpha=\frac{\varepsilon/K}{1+10\varepsilon}$. Therefore, in terms of $\delta^{(t)}$, we have

The main loop in Algorithm 3.1 terminates when $\|{\mathbf{x}^{(t)}}\|_{1}>K$ or the number of iterations exceeds a preset threshold $R=O(\varepsilon^{-3}\log^{2}n)$. For a desired accuracy parameter $\varepsilon>0$, we set $K$ to $O(\frac{1}{\varepsilon}(1+\ln n))$ so that the $\ln n$ additive term from Theorem 2.1 can be absorbed into the relative error. This additive term comes from the starting point $\mathbf{x}^{(0)}$.

The choice of $\alpha$ may seem mysterious at this point; it is chosen to prevent us from taking a step that is too big from the current solution while still making substantial progress. It ensures that

The update has small width, aka. $\sum_{i}\delta_{i}^{(t)}\mathbf{A}_{i}\preccurlyeq\varepsilon\mathbf{I}$, and

We cannot overshoot by much when exiting from the while loop, $\mathbf{1}^{\top}\delta^{(t)}\leq\varepsilon$.

Our analysis of decisionPSDP and in turn our proof of Theorem 3.1 hinge on two complementary processess: show that after $R$ steps, $\overline{\mathbf{Y}}$ is indeed a feasible primal solution, or that if the the algorithm terminates because $\|{\mathbf{x}^{(t)}}\|_{1}>K$, then $\widehat{\mathbf{x}}$ is a feasible dual solution. In the latter case, we have

For this $\widehat{\mathbf{x}}$ to be a dual solution, we still need to show that it satisfies $\sum_{i}\widehat{\mathbf{x}}_{i}\mathbf{A}_{i}\preccurlyeq\mathbf{I}$.

Let $T$ be the final iteration count (i.e., the final $t$). To meet the requirement above, we only need to show that $\frac{1}{(1+10\varepsilon)K}\mathbf{\Psi}^{(T)}\preccurlyeq\mathbf{I}$. We prove the following spectrum bound:

We prove this lemma by induction on the iteration number, resorting to properties of the MMW algorithm (Theorem 2.1), which relates the final spectral values to the “gain” derived at each intermediate step.

$\mathbf{\Psi}^{(0)}\preccurlyeq\mathbf{I}$. In other words, $\lambda_{\max}\left(\mathbf{\Psi}^{(0)}\right)=\lambda_{\max}\left(\sum_{i=1}^{n}x_{i}^{(0)}\mathbf{A}_{i}\right)\leq 1$.

Now for any iteration $t\leq T$, we can rewrite $\mathbf{\Psi}^{(t)}$ as

since both sums yield positive semidefinite matrices and the $\lambda_{\max}$ of the first sum is at most $1$ by Claim 3.3.

Proof of Lemma 3.2: We will prove (3.5) by (strong) induction on $t$. The base case of $t=0$ is true by Claim 3.3. For a given $t$, if we inductively assume that $\mathbf{\Psi}^{(\tau)}\preccurlyeq(1+10\varepsilon)K\mathbf{I}$ for all $\tau<t$, then for each $1\leq\tau<t$,

This makes Theorem 2.1 applicable, which gives

which allows us to conclude that $\mathbf{\Psi}^{(t)}\preccurlyeq(1+10\varepsilon)K\mathbf{I}$, as desired. $\blacksquare$

It remains to show the claims about the algorithm utilized in the above proof.

Our choice of $\mathbf{x}^{(0)}$ guarantees that for all $i=1,\dots,n$,

Summing across $i=1,\dots,n$ gives the desired bound.

Now every $i\in B^{(t)}$, though $B^{(t)}$ can be empty, has the property that $\mathbf{A}_{i}\bullet\mathbf{P}^{(t)}\leq(1+\varepsilon)$, so

It remains to examine the case where the algorithm returns a primal solution: Equation (3.3) gives

Furthermore, this $\overline{\mathbf{Y}}$ satisfies the primal constraints:

If $\|{\mathbf{x}^{(T)}}\|_{1}\leq K$—i.e, the algorithm exits the while-loop because it reaches $R$ iterations—then for all $i=1,\dots,n$, $\mathbf{A}_{i}\bullet\overline{\mathbf{Y}}\geq 1$.

Assume for a contradiction that there is an $i\in[n]$ such that $\mathbf{A}_{i}\bullet\overline{\mathbf{Y}}<1$. This means

Let $U=\{\tau:\mathbf{P}^{(\tau)}\bullet\mathbf{A}_{i}<1+\varepsilon\}$ be the iterations in which the $i$-th coordinate of $\mathbf{x}$ is incremented. By Markov’s inequality, the number of such iterations is bounded by $|U|<\frac{\varepsilon}{1+\varepsilon}T$. But then, every time the $i$-th coordinate changes, it increases by a factor of $1+\alpha$, so

Because $\|{\mathbf{x}^{(T)}}\|_{1}\leq K$, the algorithm exits the while-loop with $T=R$. Therefore, for $0<\varepsilon<1$,

as $\mathop{\mathsf{Tr}}\left[\mathbf{A}_{i}\right]\leq O(n^{3})$ by Lemma 2.2. This is a contradiction to $\|{\mathbf{x}^{(T)}}\|_{1}\leq K$, which proves the lemma.

Proof of Theorem 3.1: The algorithm terminates after at most $R$ iterations. Notice that $B^{(t)}$ may be empty in some iterations but this does not harm the algorithm nor the proof. It is standard to check that $R=O(\varepsilon^{-3}\log^{2}n)$. If we do run this many iterations, then Lemma 3.6 gives that we terminate with a primal solution.

Otherwise, Lemma 3.2 gives that at any point in the algorithm, the solution vector $\mathbf{x}^{(t)}$ satisfies $\sum_{i}x^{(t)}_{i}\mathbf{A}_{i}\preccurlyeq(1+10\varepsilon)K\mathbf{I}$. Together with Equation (3.4), we know that any $\widehat{\mathbf{x}}$ returned $\|{\mathbf{x}^{*}}\|_{1}\geq 1-10\varepsilon$ and

Thus, $\widehat{\mathbf{x}}$ is indeed a dual solution with value at least $1-10\varepsilon$. Replacing $\varepsilon$ with $\varepsilon/10$ then meets the requirements of the decision problem. $\blacksquare$

Matrix Exponential Evaluation

We describe a fast algorithm for computing the matrix dot product of a positive semidefinite matrix and the matrix exponential of another positive semidefinite matrix.

There is an algorithm bigDotExp that when given a $m$-by-$m$ matrix $\mathbf{\Phi}$ with $p$ non-zero entries, $\kappa\geq\max\{1,\|{\mathbf{\Phi}}\|_{2}\}$, and $m$-by-$m$ matrices $\mathbf{A}_{i}$ in factorized form $\mathbf{A}_{i}=\mathbf{Q}_{i}\mathbf{Q}_{i}^{\top}$ where the total number of nonzeros across all $\mathbf{Q}_{i}$ is $q$; $\mbox{{bigDotExp}}(\mathbf{\Phi},\{\mathbf{A}_{i}=\mathbf{Q}_{i}\mathbf{Q}_{i}^{\top}\}_{i=1}^{n})$ computes $(1\pm\varepsilon)$ approximations to all $\exp{(\mathbf{\Phi})}\bullet\mathbf{A}_{i}$ in $O(\kappa\log{m}\log(\nicefrac{{1}}{{\varepsilon}}))$ depth and $O(\frac{1}{\varepsilon^{2}}(\kappa\log(\nicefrac{{1}}{{\epsilon}})p+q)\log{m})$ work.

The idea behind Theorem 4.1 is to approximate the matrix exponential using a low-degree polynomial because evaluating matrix exponentials exactly is costly. For this, we will apply the following lemma, reproduced from Lemma 6 in :

If $\mathbf{B}$ is a PSD matrix such that $\|{\mathbf{B}}\|_{2}\leq\kappa$, then the operator

The given factorization of each $\mathbf{A}_{i}$ allows us to write $\exp{(\mathbf{\Phi})}\bullet\mathbf{A}_{i}$ as the $2$-norm of a vector:

By Lemma 4.2, it suffices to evaluate $\widehat{\mathbf{B}}\bullet\mathbf{A}_{i}$ where $\widehat{\mathbf{B}}$ is an approximation to $\mathbf{B}=\exp(\tfrac{1}{2}\mathbf{\Phi})$. To further reduce the work, we can apply the Johnson-Lindenstrauss transformation to reduce the length of the vectors to $O(\log{m})$; specifically, we find a $O(\frac{1}{\varepsilon^{2}}\log{m})\times m$ Gaussian matrix $\mathbf{\Pi}$ and evaluate

Since $\mathbf{\Pi}$ only has $O(\frac{1}{\varepsilon^{2}}\log{m})$ rows, we can compute $\mathbf{\Pi}\widehat{\mathbf{B}}$ using $O(\log{m})$ evaluations of $\widehat{\mathbf{B}}$. The work/depth bounds follow from doing each of the evaluations of $\widehat{\mathbf{B}}\Pi_{i}$, where $\Pi_{i}$ denotes the $i$-th column of $\mathbf{\Pi}$, and matrix-vector multiplies involving $\mathbf{\Phi}$ in parallel.

Conclusion

We presented a simple NC parallel algorithm for packing SDPs that requires $O(\frac{1}{\varepsilon^{4}}\log^{4}n\log(\frac{1}{\varepsilon}))$ iterations, where each iteration involves only simple matrix operations and computing the trace of the product of a matrix exponential and a positive semidefinite matrix. When a positive SDP is given in a factorized form, we showed how the dot product with matrix exponential can be implemented in nearly-linear work, leading to an algorithm with $\widetilde{O}(m+n+q)$ work, where $n$ is the number of constraint matrices, $m$ is the dimension of these matrices, and $q$ is the total number of nonzero entries in the factorization.

Compared to the situation with positive LPs, the classification of positive SDPs is much richer because packing/covering constraints can take many forms, either as matrices (e.g. $\sum_{i=1}^{n}x_{i}\mathbf{A}_{i}\preccurlyeq\mathbf{I}$ for packing, $\sum_{i=1}^{n}x_{i}\mathbf{A}_{i}\succcurlyeq\mathbf{I}$ for covering) or as dot products between matrices (e.g. $\mathbf{A}_{i}\bullet\mathbf{Y}\leq 1$ for packing, $\mathbf{A}_{i}\bullet\mathbf{Y}\geq 1$ for covering). The positive SDPs studied in and our paper should be compared with the closely related notion of covering SDPs studied by Iyengar et al ; however, among the applications they examine, only the beamforming SDP relaxation discussed in Section 2.2 of falls completely within the framework of packing SDPs as defined in 2.10. Problems such as MaxCut and SparsestCut require additional matrix-based packing constraints. We believe extending these algorithms to solve mixed packing/covering SDPs is an interesting direction for future work.

Acknowledgments

This work is partially supported by the National Science Foundation under grant numbers CCF-1018463, CCF-1018188, and CCF-1016799 and by generous gifts from IBM, Intel, and Microsoft. Richard Peng was partly supported by a Microsoft Research PhD. Fellowship.

We thank the SPAA reviewers, as well as Guy Blelloch and Gary Miller for suggestions that helped improve this paper. While making this revision, we benefitted greatly from discussions with Jon Kelner and Di Wang.

References

Appendix A Normalized Positive SDPs

This is the same transformation that Jain and Yao presented ; we only present it here for easy reference.

Consider the primal program in (1.4). It suffices to show that it can be transformed into the following program without changing the optimal value:

We can make the following assumptions without loss of generality: First, $b_{i}>0$ for all $i=1,\dots,m$ because if $b_{i}$ were , we could have thrown it away. Second, all $\mathbf{A}_{i}$’s are the support of $\mathbf{C}$, or otherwise we know that the corresponding dual variable must be set to and therefore can be removed right away. Therefore, we will treat $C$ as having a full-rank, allowing us to define

It is not hard to verify that the normalized program (A.4) has the same optimal value as the original SDP (1.4).

Note that if we’re given factorization of $\mathbf{A}_{i}$ into $\mathbf{Q}_{i}\mathbf{Q}_{i}^{\top}$, then $\mathbf{B}_{i}$ can also be factorized as:

Furthermore, it can be checked that the dual of the normalized program is the same as the dual in Equation 2.10.