Parallel Coordinate Descent Methods for Big Data Optimization

Keywords:

Parallel coordinate descent, big data optimization, partial separability, huge-scale optimization, iteration complexity, expected separable over-approximation, composite objective, convex optimization, LASSO.

Introduction

Recently there has been a surge in interest in the design of algorithms suitable for solving convex optimization problems with a huge number of variables . Indeed, the size of problems arising in fields such as machine learning , network analysis , PDEs , truss topology design and compressed sensing usually grows with our capacity to solve them, and is projected to grow dramatically in the next decade. In fact, much of computational science is currently facing the “big data” challenge, and this work is aimed at developing optimization algorithms suitable for the task.

Coordinate descent methods.

Parallelization.

We wish to point out that for truly huge-scale problems it is absolutely necessary to parallelize. This is in line with the rise and ever increasing availability of high performance computing systems built around multi-core processors, GPU-accelerators and computer clusters, the success of which is rooted in massive parallelization. This simple observation, combined with the remarkable scalability of serial CDMs, leads to our belief that the study of parallel coordinate descent methods (PCDMs) is a very timely topic.

Research Idea.

The work presented in this paper was motivated by the desire to answer the following question:

Under what natural and easily verifiable structural assumptions on the objective function does parallelization of a coordinate descent method lead to acceleration?

Our starting point was the following simple observation. Assume that we wish to minimize a separable function $F$ of $n$ variables (i.e., a function that can be written as a sum of $n$ functions each of which depends on a single variable only). For simplicity, in this thought experiment, assume that there are no constraints. Clearly, the problem of minimizing $F$ can be trivially decomposed into $n$ independent univariate problems. Now, if we have $n$ processors/threads/cores, each assigned with the task of solving one of these problems, the number of parallel iterations should not depend on the dimension of the problemFor simplicity, assume the distance from the starting point to the set of optimal solutions does not depend on the dimension.. In other words, we get an $n$-times speedup compared to the situation with a single processor only. Note that any parallel algorithm of this type can be viewed as a parallel coordinate descent method. Hence, a PCDM with $n$ processors should be $n$-times faster than a serial one. If $\tau$ processors are used instead, where $1\leq\tau\leq n$, one would expect a $\tau$-times speedup.

By extension, one would perhaps expect that optimization problems with objective functions which are “close to being separable” would also be amenable to acceleration by parallelization, where the acceleration factor $\tau$ would be reduced with the reduction of the “degree of separability”. One of the main messages of this paper is an affirmative answer to this. Moreover, we give explicit and simple formulae for the speedup factors.

As it turns out, and as we discuss later in this section, many real-world big data optimization problems are, quite naturally, “close to being separable”. We believe that this means that PCDMs is a very promising class of algorithms when it comes to solving structured big data optimization problems.

Minimizing a partially separable composite objective.

where $f$ is a (block) partially separable smooth convex function and $\Omega$ is a simple (block) separable convex function. We allow $\Omega$ to have values in $\mathbf{R}\cup\{\infty\}$, and for regularization purposes we assume $\Omega$ is proper and closed. While (1) is seemingly an unconstrained problem, $\Omega$ can be chosen to model simple convex constraints on individual blocks of variables. Alternatively, this function can be used to enforce a certain structure (e.g., sparsity) in the solution. For a more detailed account we refer the reader to . Further, we assume that this problem has a minimum ($F^{*}>-\infty$). What we mean by “smoothness” and “simplicity” will be made precise in the next section.

Let us now describe the key concept of partial separability. Let $x\in\mathbf{R}^{N}$ be decomposed into $n$ non-overlapping blocks of variables $x^{(1)},\dots,x^{(n)}$ (this will be made precise in Section 2). We assume throughout the paper that $f:\mathbf{R}^{N}\to\mathbf{R}$ is partially separable of degree $\omega$, i.e., that it can be written in the form

where $\mathcal{J}$ is a finite collection of nonempty subsets of ${[n]}\stackrel{{\scriptstyle\text{def}}}{{=}}\{1,2,\dots,n\}$ (possibly containing identical sets multiple times), $f_{J}$ are differentiable convex functions such that $f_{J}$ depends on blocks $x^{(i)}$ for $i\in J$ only, and

Clearly, $1\leq\omega\leq n$. The PCDM algorithms we develop and analyze in this paper only need to know $\omega$, they do not need to know the decomposition of $f$ giving rise to this $\omega$.

Examples of partially separable functions.

Many objective functions naturally encountered in the big data setting are partially separable. Here we give examples of three loss/objective functions frequently used in the machine learning literature and also elsewhere. For simplicity, we assume all blocks are of size 1 (i.e., $N=n$). Let

where $m$ is the number of examples, $x\in\mathbf{R}^{n}$ is the vector of features, $(A_{j},y_{j})\in\mathbf{R}^{n}\times\mathbf{R}$ are labeled examples and ${\cal L}$ is one of the three loss functions listed in Table 1. Let $A\in\mathbf{R}^{m\times n}$ with row $j$ equal to $A_{j}^{T}$.

Often, each example depends on a few features only; the maximum over all features is the degree of partial separability $\omega$. More formally, note that the $j$-th function in the sum (4) in all cases depends on $\|A_{j}\|_{0}$ coordinates of $x$ (the number of nonzeros in the $j$-th row of $A$) and hence $f$ is partially separable of degree

All three functions of Table 1 are smooth (based on the definition of smoothness in the next section). We refer the reader to for more examples of interesting (but nonsmooth) partially separable functions arising in graph cuts and matrix completion.

Brief literature review.

Several papers were written recently studying the iteration complexity of serial CDMs of various flavours and in various settings. We will only provide a brief summary here, for a more detailed account we refer the reader to .

Classical CDMs update the coordinates in a cyclic order; the first attempt at analyzing the complexity of such a method is due to . Stochastic/randomized CDMs, that is, methods where the coordinate to be updated is chosen randomly, were first analyzed for quadratic objectives , later independently generalized to $L_{1}$-regularized problems and smooth block-structured problems , and finally unified and refined in . The problems considered in the above papers are either unconstrained or have (block) separable constraints. Recently, randomized CDMs were developed for problems with linearly coupled constraints .

A greedy CDM for $L_{1}$-regularized problems was first analyzed in ; more work on this topic include . A CDM with inexact updates was first proposed and analyzed in . Partially separable problems were independently studied in , where an asynchronous parallel stochastic gradient algorithm was developed to solve them.

When writing this paper, the authors were aware only of the parallel CDM proposed and analyzed in . Several papers on the topic appeared around the time this paper was finalized or after . Further papers on various aspects of the topic of parallel CDMs, building on the work in this paper, include .

Contents.

We start in Section 2 by describing the block structure of the problem, establishing notation and detailing assumptions. Subsequently we propose and comment in detail on two parallel coordinate descent methods. In Section 3 we summarize the main contributions of this paper. In Section 4 we deal with issues related to the selection of the blocks to be updated in each iteration. It will involve the development of some elementary random set theory. Sections 5-6 deal with issues related to the computation of the update to the selected blocks and develop a theory of Expected Separable Overapproximation (ESO), which is a novel tool we propose for the analysis of our algorithms. In Section 7 we analyze the iteration complexity of our methods and finally, Section 8 reports on promising computational results. For instance, we conduct an experiment with a big data (cca 350GB) LASSO problem with a billion variables. We are able to solve the problem using one of our methods on a large memory machine with 24 cores in 2 hours, pushing the difference between the objective value at the starting iterate and the optimal point from $10^{22}$ down to $10^{-14}$. We also conduct experiments on real data problems coming from machine learning.

Parallel Block Coordinate Descent Methods

In Section 2.1 we formalize the block structure of the problem, establish notationTable 8 in the appendix summarizes some of the key notation used frequently in the paper. that will be used in the rest of the paper and list assumptions. In Section 2.2 we propose two parallel block coordinate descent methods and comment in some detail on the steps.

Some elements of the setup described in this section was initially used in the analysis of block coordinate descent methods by Nesterov (e.g., block structure, weighted norms and block Lipschitz constants).

The block structure of (1) is given by a decomposition of $\mathbf{R}^{N}$ into $n$ subspaces as follows. Let $U\in\mathbf{R}^{N\times N}$ be a column permutationThe reason why we work with a permutation of the identity matrix, rather than with the identity itself, as in , is to enable the blocks being formed by nonconsecutive coordinates of $x$. This way we establish notation which makes it possible to work with (i.e., analyze the properties of) multiple block decompositions, for the sake of picking the best one, subject to some criteria. Moreover, in some applications the coordinates of $x$ have a natural ordering to which the natural or efficient block structure does not correspond. of the $N\times N$ identity matrix and further let $U=[U_{1},U_{2},\dots,U_{n}]$ be a decomposition of $U$ into $n$ submatrices, with $U_{i}$ being of size $N\times N_{i}$, where $\sum_{i}N_{i}=N$.

Any vector $x\in\mathbf{R}^{N}$ can be written uniquely as

where $x^{(i)}\in\mathbf{R}^{N_{i}}$. Moreover, $x^{(i)}=U_{i}^{T}x$.

Noting that $UU^{T}=\sum_{i}U_{i}U_{i}^{T}$ is the $N\times N$ identity matrix, we have $x=\sum_{i}U_{i}U_{i}^{T}x$. Let us now show uniqueness. Assume that $x=\sum_{i}U_{i}x_{1}^{(i)}=\sum_{i}U_{i}x_{2}^{(i)}$, where $x_{1}^{(i)},x_{2}^{(i)}\in\mathbf{R}^{N_{i}}$. Since

for every $j$ we get $0=U_{j}^{T}(x-x)=U_{j}^{T}\sum_{i}U_{i}(x_{1}^{(i)}-x_{2}^{(i)})=x_{1}^{(j)}-x_{2}^{(j)}$. ∎

In view of the above proposition, from now on we write $x^{(i)}\stackrel{{\scriptstyle\text{def}}}{{=}}U_{i}^{T}x\in\mathbf{R}^{N_{i}}$, and refer to $x^{(i)}$ as the $i$-th block of $x$. The definition of partial separability in the introduction is with respect to these blocks. For simplicity, we will sometimes write $x=(x^{(1)},\dots,x^{(n)})$.

For $S\subset{[n]}$ and $x\in\mathbf{R}^{N}$ we write

That is, given $x\in\mathbf{R}^{N}$, $x_{[S]}$ is the vector in $\mathbf{R}^{N}$ whose blocks $i\in S$ are identical to those of $x$, but whose other blocks are zeroed out. In view of Proposition 1, we can equivalently define $x_{[S]}$ block-by-block as follows

Inner products.

The standard Euclidean inner product in spaces $\mathbf{R}^{N}$ and $\mathbf{R}^{N_{i}}$, $i\in{[n]}$, will be denoted by $\langle\cdot,\cdot\rangle$. Letting $x,y\in\mathbf{R}^{N}$, the relationship between these inner products is given by

For any $w\in\mathbf{R}^{n}$ and $x,y\in\mathbf{R}^{N}$ we further define

For vectors $z=(z_{1},\dots,z_{n})^{T}\in\mathbf{R}^{n}$ and $w=(w_{1},\dots,w_{n})^{T}\in\mathbf{R}^{n}$ we write $w\odot z\stackrel{{\scriptstyle\text{def}}}{{=}}(w_{1}z_{1},\dots,w_{n}z_{n})^{T}$.

Norms.

Spaces $\mathbf{R}^{N_{i}}$, $i\in{[n]}$, are equipped with a pair of conjugate norms: $\|t\|_{(i)}\stackrel{{\scriptstyle\text{def}}}{{=}}\langle B_{i}t,t\rangle^{1/2}$, where $B_{i}$ is an $N_{i}\times N_{i}$ positive definite matrix and $\|t\|_{(i)}^{*}\stackrel{{\scriptstyle\text{def}}}{{=}}\max_{\|s\|_{(i)}\leq 1}\langle s,t\rangle=\langle B_{i}^{-1}t,t\rangle^{1/2}$, $t\in\mathbf{R}^{N_{i}}$. For $w\in\mathbf{R}^{n}_{++}$, define a pair of conjugate norms in $\mathbf{R}^{N}$ by

Note that these norms are induced by the inner product (9) and the matrices $B_{1},\dots,B_{n}$. Often we will use $w=L\stackrel{{\scriptstyle\text{def}}}{{=}}(L_{1},L_{2},\dots,L_{n})^{T}\in\mathbf{R}^{n}$, where the constants $L_{i}$ are defined below.

Smoothness of f𝑓f.

We assume throughout the paper that the gradient of $f$ is block Lipschitz, uniformly in $x$, with positive constants $L_{1},\dots,L_{n}$, i.e., that for all $x\in\mathbf{R}^{N}$, $i\in{[n]}$ and $t\in\mathbf{R}^{N_{i}}$,

where $\nabla_{i}f(x)\stackrel{{\scriptstyle\text{def}}}{{=}}(\nabla f(x))^{(i)}=U^{T}_{i}\nabla f(x)\in\mathbf{R}^{N_{i}}$. An important consequence of (11) is the following standard inequality :

Separability of ΩΩ\Omega.

We assume thatFor examples of separable and block separable functions we refer the reader to . For instance, $\Omega(x)=\|x\|_{1}$ is separable and block separable (used in sparse optimization); and $\Omega(x)=\sum_{i}\|x^{(i)}\|$, where the norms are standard Euclidean norms, is block separable (used in group lasso). One can model block constraints by setting $\Omega_{i}(x^{(i)})=0$ for $x\in X_{i}$, where $X_{i}$ is some closed convex set, and $\Omega_{i}(x^{(i)})=+\infty$ for $x\notin X_{i}$. $\Omega:\mathbf{R}^{N}\to\mathbf{R}\cup\{+\infty\}$ is (block) separable, i.e., that it can be decomposed as follows:

where the functions $\Omega_{i}:\mathbf{R}^{N_{i}}\to\mathbf{R}\cup\{+\infty\}$ are convex and closed.

Strong convexity.

In one of our two complexity results (Theorem 20) we will assume that either $f$ or $\Omega$ (or both) is strongly convex. A function $\phi:\mathbf{R}^{N}\to\mathbf{R}\cup\{+\infty\}$ is strongly convex with respect to the norm $\|\cdot\|_{w}$ with convexity parameter $\mu_{\phi}(w)\geq 0$ if for all $x,y\in\operatorname{dom}\phi$,

where $\phi^{\prime}(x)$ is any subgradient of $\phi$ at $x$. The case with $\mu_{\phi}(w)=0$ reduces to convexity. Strong convexity of $F$ may come from $f$ or $\Omega$ (or both); we write $\mu_{f}(w)$ (resp. $\mu_{\Omega}(w)$) for the (strong) convexity parameter of $f$ (resp. $\Omega$). It follows from (14) that

The following characterization of strong convexity will be useful. For all $x,y\in\operatorname{dom}\phi$ and $\lambda\in$,

It can be shown using (12) and (14) that $\mu_{f}(w)\leq\tfrac{L_{i}}{w_{i}}$.

2 Algorithms

In this paper we develop and study two generic parallel coordinate descent methods. The main method is PCDM1; PCDM2 is its “regularized” version which explicitly enforces monotonicity. As we will see, both of these methods come in many variations, depending on how Step 3 is performed.

Let us comment on the individual steps of the two methods.

Step 3. At the beginning of iteration $k$ we pick a random set ($S_{k}$) of blocks to be updated (in parallel) during that iteration. The set $S_{k}$ is a realization of a random set-valued mapping $\hat{S}$ with values in $2^{[n]}$ or, more precisely, it the sets $S_{k}$ are iid random sets with the distribution of $\hat{S}$. For brevity, in this paper we refer to such a mapping by the name sampling. We limit our attention to uniform samplings, i.e., random sets having the following property: $\mathbf{P}(i\in\hat{S})$ is independent of $i$. That is, the probability that a block gets selected is the same for all blocks. Although we give an iteration complexity result covering all such samplings (provided that each block has a chance to be updated, i.e., $\mathbf{P}(i\in\hat{S})>0$), there are interesting subclasses of uniform samplings (such as doubly uniform and nonoverlapping uniform samplings; see Section 4) for which we give better results.

Step 4. For $x\in\mathbf{R}^{N}$ we defineA similar map was used in (with $\Omega\equiv 0$ and $\beta=1$) and (with $\beta=1$) in the analysis of serial coordinate descent methods in the smooth and composite case, respectively. In loose terms, the novelty here is the introduction of the parameter $\beta$ and in developing theory which describes what value $\beta$ should have. Maps of this type are known as composite gradient mapping in the literature, and were introduced in .

and $\beta>0$, $w=(w_{1},\dots,w_{n})^{T}\in\mathbf{R}^{n}_{++}$ are parameters of the method that we will comment on later. Note that in view of (5), (10) and (13), $H_{\beta,w}(x,\cdot)$ is block separable;

Consequently, we have $h(x)=(h^{(1)}(x),\cdots,h^{(n)}(x))\in\mathbf{R}^{N}$, where

We mentioned in the introduction that besides (block) separability, we require $\Omega$ to be “simple”. By this we mean that the above optimization problem leading to $h^{(i)}(x)$ is “simple” (e.g., it has a closed-form solution). Recall from (8) that $(h(x_{k}))_{[S_{k}]}$ is the vector in $\mathbf{R}^{N}$ identical to $h(x_{k})$ except for blocks $i\notin S_{k}$, which are zeroed out. Hence, Step 4 of both methods can be written as follows:

Parameters $\beta$ and $w$ depend on $f$ and $\hat{S}$ and stay constant throughout the algorithm. We are not ready yet to explain why the update is computed via (17) and (18) because we need technical tools, which will be developed in Section 4, to do so. Here it suffices to say that the parameters $\beta$ and $w$ come from a separable quadratic overapproximation of $\mathbf{E}[f(x+h_{[\hat{S}]})]$, viewed as a function of $h\in\mathbf{R}^{N}$. Since expectation is involved, we refer to this by the name Expected Separable Overapproximation (ESO). This novel concept, developed in this paper, is one of the main tools of our complexity analysis. Section 5 motivates and formalizes the concept, answers the why question, and develops some basic ESO theory.

Section 6 is devoted to the computation of $\beta$ and $w$ for partially separable $f$ and various special classes of uniform samplings $\hat{S}$. Typically we will have $w_{i}=L_{i}$, while $\beta$ will depend on easily computable properties of $f$ and $\hat{S}$. For example, if $\hat{S}$ is chosen as a subset of ${[n]}$ of cardinality $\tau$, with each subset chosen with the same probability (we say that $\hat{S}$ is $\tau$-nice) then, assuming $n>1$, we may choose $w=L$ and $\beta=1+\tfrac{(\omega-1)(\tau-1)}{n-1}$, where $\omega$ is the degree of partial separability of $f$. More generally, if $\hat{S}$ is any uniform sampling with the property $|\hat{S}|=\tau$ with probability 1, then we may choose $w=L$ and $\beta=\min\{\omega,\tau\}$. Note that in both cases $w=L$ and that the latter $\beta$ is always larger than (or equal to) the former one. This means, as we will see in Section 7, that we can give better complexity results for the former, more specialized, sampling. We analyze several more options for $\hat{S}$ than the two just described, and compute parameters $\beta$ and $w$ that should be used with them (for a summary, see Table 4).

Step 5. The reason why, besides PCDM1, we also consider PCDM2, is the following: in some situations we are not able to analyze the iteration complexity of PCDM1 (non-strongly-convex $F$ where monotonicity of the method is not guaranteed by other means than by directly enforcing it by inclusion of Step 5). Let us remark that this issue arises for general $\Omega$ only. It does not exist for $\Omega=0$, $\Omega(\cdot)=\lambda\|\cdot\|_{1}$ and for $\Omega$ encoding simple constraints on individual blocks; in these cases one does not need to consider PCDM2. Even in the case of general $\Omega$ we sometimes get monotonicity for free, in which case there is no need to enforce it. Let us stress, however, that we do not recommend implementing PCDM2 as this would introduce too much overhead; in our experience PCDM1 works well even in cases when we can only analyze PCDM2.

Smmary of Contributions

In this section we summarize the main contributions of this paper (not in order of significance).

Problem generality. We give the first complexity analysis for parallel coordinate descent methods for problem (1) in its full generality.

Complexity. We show theoretically (Section 7) and numerically (Section 8) that PCDM accelerates on its serial counterpart for partially separable problems. In particular, we establish two complexity theorems giving lower bounds on the number of iterations $k$ sufficient for one or both of the PCDM variants (for details, see the precise statements in Section 7) to produce a random iterate $x_{k}$ for which the problem is approximately solved with high probability, i.e., $\mathbf{P}(F(x_{k})-F^{*}\leq\epsilon)\geq 1-\rho$. The results, summarized in Table 2, hold under the standard assumptions listed in Section 2.1 and the additional assumption that $f,\hat{S},\beta$ and $w$ satisfy the following inequality for all $x,h\in\mathbf{R}^{N}$:

This inequality, which we call Expected Separable Overapproximation (ESO), is the main new theoretical tool that we develop in this paper for the analysis of our methods (Sections 4-6 are devoted to the development of this theory).

The main observation here is that as the average number of block updates per iteration increases (say, $\hat{\tau}=\operatorname{\mathbf{E}}[|\hat{S}|]$), enabled by the utilization of more processors, the leading term in the complexity estimate, $n/\hat{\tau}$, decreases in proportion. However, $\beta$ will generally grow with $\hat{\tau}$, which has an adverse effect on the speedup. Much of the theory in this paper goes towards producing formulas for $\beta$ (and $w$), for partially separable $f$ and various classes of uniform samplings $\hat{S}$. Naturally, the ideal situation is when $\beta$ does not grow with $\hat{\tau}$ at all, or if it only grows very slowly. We show that this is the case for partially separable functions $f$ with small $\omega$. For instance, in the extreme case when $f$ is separable ($\omega=1$), we have $\beta=1$ and we obtain linear speedup in $\hat{\tau}$. As $\omega$ increases, so does $\beta$, depending on the law governing $\hat{S}$. Formulas for $\beta$ and $\omega$ for various samplings $\hat{S}$ are summarized in Table 4.

Algorithm unification. Depending on the choice of the block structure (as implied by the choice of $n$ and the matrices $U_{1},\dots,U_{n}$) and the way blocks are selected at every iteration (as given by the choice of $\hat{S}$), our framework encodes a family of known and new algorithmsAll the methods are in their proximal variants due to the inclusion of the term $\Omega$ in the objective. (see Table 3).

In particular, PCDM is the first method which “continuously” interpolates between serial coordinate descent and gradient (by manipulating $n$ and/or $\mathbf{E}[|\hat{S}|]$).

Partial separability. We give the first analysis of a coordinate descent type method dealing with a partially separable loss / objective. In order to run the method, we need to know the Lipschitz constants $L_{i}$ and the degree of partial separability $\omega$. It is crucial that these quantities are often easily computable/predictable in the huge-scale setting. For example, if $f(x)=\tfrac{1}{2}\|Ax-b\|^{2}$ and we choose all blocks to be of size $1$, then $L_{i}$ is equal to the squared Euclidean norm of the $i$-th column of $A$ and $\omega$ is equal to the maximum number of nonzeros in a row of $A$. Many problems in the big data setting have small $\omega$ compared to $n$.

Choice of blocks. To the best of our knowledge, existing randomized strategies for paralleling gradient-type methods (e.g., ) assume that $\hat{S}$ (or an equivalent thereof, based on the method) is chosen as a subset of $[n]$ of a fixed cardinality, uniformly at random. We refer to such $\hat{S}$ by the name nice sampling in this paper. We relax this assumption and our treatment is hence much more general. In fact, we allow for $\hat{S}$ to be any uniform sampling. It is possible to further consider nonuniform samplingsRevision note: See ., but this is beyond the scope of this paper.

In particular, as a special case, our method allows for a variable number of blocks to be updated throughout the iterations (this is achieved by the introduction of doubly uniform samplings). This may be useful in some settings such as when the problem is being solved in parallel by $\tau$ unreliable processors each of which computes its update $h^{(i)}(x_{k})$ with probability $p_{b}$ and is busy/down with probability $1-p_{b}$ (binomial sampling).

Uniform, doubly uniform, nice, binomial and other samplings are defined, and their properties studied, in Section 4.

ESO and formulas for $\beta$ and $w$. In Table 4 we list parameters $\beta$ and $w$ for which ESO inequality (19) holds. Each row corresponds to a specific sampling $\hat{S}$ (see Section 4 for the definitions). The last 5 samplings are special cases of one or more of the first three samplings. Details such as what is $\nu,\gamma$ and “monotonic” ESO are explained in appropriate sections later in the text. When a specific sampling $\hat{S}$ is used in the algorithm to select blocks in each iteration, the corresponding parameters $\beta$ and $w$ are to be used in the method for the computation of the update (see (17) and (18)).

En route to proving the iteration complexity results for our algorithms, we develop a theory of deterministic and expected separable overapproximation (Sections 5 and 6) which we believe is of independent interest, too. For instance, methods based on ESO can be compared favorably to the Diagonal Quadratic Approximation (DQA) approach used in the decomposition of stochastic optimization programs .

Parallelization speedup. Our complexity results can be used to derive theoretical parallelization speedup factors. For several variants of our method, in case of a non-strongly convex objective, these are given in Section 7.1 (Table 5). For instance, in the case when all block are updated at each iteration (we later refer to $\hat{S}$ having this property by the name fully parallel sampling), the speedup factor is equal to $\tfrac{n}{\omega}$. If the problem is separable ($\omega=1$), the speedup is equal to $n$; if the problem is not separable ($\omega=n$), there may be no speedup. For strongly convex $F$ the situation is even better; the details are given in Section 7.2.

Relationship to existing results. To the best of our knowledge, there are just two papers analyzing a parallel coordinate descent algorithm for convex optimization problems. In the first paper all blocks are of size $1$, $\hat{S}$ corresponds to what we call in this paper a $\tau$-nice sampling (i.e., all sets of $\tau$ coordinates are updated at each iteration with equal probability) and hence their algorithm is somewhat comparable to one of the many variants of our general method. While the analysis in works for a restricted range of values of $\tau$, our results hold for all $\tau\in{[n]}$. Moreover, the authors consider a more restricted class of functions $f$ and the special case $\Omega=\lambda\|x\|_{1}$, which is simpler to analyze. Lastly, the theoretical speedups obtained in , when compared to the serial CDM method, depend on a quantity $\sigma$ that is hard to compute in big data settings (it involves the computation of an eigenvalue of a huge-scale matrix). Our speedups are expressed in terms of natural and easily computable quantity: the degree $\omega$ of partial separability of $f$. In the setting considered by , in which more structure is available, it turns out that $\omega$ is an upper boundRevision note requested by a reviewer: In the time since this paper was posted to arXiv, a number of follow-up papers were written analyzing parallel coordinate descent methods and establishing connections between a discrete quantity analogous to $\omega$ (degree of partial/Nesterov separability) and a spectral quantity analogous to $\sigma$ (largest eigenvalue of a certain matrix), most notably . See also , which uses a spectral quantity, which can be directly compared to $\omega$. on $\sigma$. Hence, we show that one can develop the theory in a more general setting, and that it is not necessary to compute $\sigma$ (which may be complicated in the big data setting). The parallel CDM method of the second paper only allows all blocks to be updated at each iteration. Unfortunately, the analysis (and the method) is too coarse as it does not offer any theoretical speedup when compared to its serial counterpart. In the special case when only a single block is updated in each iteration, uniformly at random, our theoretical results specialize to those established in .

Computations. We demonstrate that our method is able to solve a LASSO problem involving a matrix with a billion columns and 2 billion rows on a large memory node with 24 cores in 2 hours (Section 8), achieving a $20\times$ speedup compared to the serial variant and pushing the residual by more than 30 degrees of magnitude. While this is done on an artificial problem under ideal conditions (controlling for small $\omega$), large speedups are possible in real data with $\omega$ small relative to $n$. We also perform additional experiments on real data sets from machine learning (e.g., training linear SVMs) to illustrate that the predictions of our theory match reality.

Code. The open source code with an efficient implementation of the algorithm(s) developed in this paper is published here: http://code.google.com/p/ac-dc/.

Block Samplings

In Step 3 of both PCDM1 and PCDM2 we choose a random set of blocks $S_{k}$ to be updated at the current iteration. Formally, $S_{k}$ is a realization of a random set-valued mapping $\hat{S}$ with values in $2^{{[n]}}$, the collection of subsets of $[n]$. For brevity, in this paper we refer to $\hat{S}$ by the name sampling. A sampling $\hat{S}$ is uniquely characterized by the probability mass function

that is, by assigning probabilities to all subsets of ${[n]}$. Further, we let $p=(p_{1},\dots,p_{n})^{T}$, where

In Section 4.1 we describe those samplings for which we analyze our methods and in Section 4.2 we prove several technical results, which will be useful in the rest of the paper.

A sampling is proper if $p_{i}>0$ for all blocks $i$. That is, from the perspective of PCDM, under a proper sampling each block gets updated with a positive probability at each iteration. Clearly, PCDM can not converge for a sampling that is not proper.

A sampling $\hat{S}$ is uniform if all blocks get updated with the same probability, i.e., if $p_{i}=p_{j}$ for all $i,j$. We show in (33) that, necessarily, $p_{i}=\tfrac{\mathbf{E}[|\hat{S}|]}{n}$. Further, we say $\hat{S}$ is nil if $\mathbf{P}(\emptyset)=1$. Note that a uniform sampling is proper if and only if it is not nil.

All our iteration complexity results in this paper are for PCDM used with a proper uniform sampling (see Theorems 19 and 20) for which we can compute $\beta$ and $w$ giving rise to an inequality (we we call “expected separable overapproximation”) of the form (43). We derive such inequalities for all proper uniform samplings (Theorem 12) as well as refined results for two special subclasses thereof: doubly uniform samplings (Theorem 15) and nonoverlapping uniform samplings (Theorem 13). We will now give the definitions:

Doubly Uniform (DU) samplings. A DU sampling is one which generates all sets of equal cardinality with equal probability. That is, $\mathbf{P}(S^{\prime})=\mathbf{P}(S^{\prime\prime})$ whenever $|S^{\prime}|=|S^{\prime\prime}|$. The name comes from the fact that this definition postulates a different uniformity property, “standard” uniformity is a consequence. Indeed, let us show that a DU sampling is necessarily uniform. Let $q_{j}=\mathbf{P}(|\hat{S}|=j)$ for $j=0,1,\dots,n$ and note that from the definition we know that whenever $S$ is of cardinality $j$, we have $\mathbf{P}(S)=q_{j}/{n\choose j}$. Finally, using this we obtain

It is clear that each DU sampling is uniquely characterized by the vector of probabilities $q$; its density function is given by

Nonoverlapping Uniform (NU) samplings. A NU sampling is one which is uniform and which assigns positive probabilities only to sets forming a partition of ${[n]}$. Let $S^{1},S^{2},\dots,S^{l}$ be a partition of ${[n]}$, with $|S^{j}|>0$ for all $j$. The density function of a NU sampling corresponding to this partition is given by

Note that $\mathbf{E}[|\hat{S}|]=\tfrac{n}{l}$.

Let us now describe several interesting special cases of DU and NU samplings:

Nice sampling. Fix $1\leq\tau\leq n$. A $\tau$-nice sampling is a DU sampling with $q_{\tau}=1$.

Interpretation: There are $\tau$ processors/threads/cores available. At the beginning of each iteration we choose a set of blocks using a $\tau$-nice sampling (i.e., each subset of $\tau$ blocks is chosen with the same probability), and assign each block to a dedicated processor/thread/core. Processor assigned with block $i$ would compute and apply the update $h^{(i)}(x_{k})$. This is the sampling we use in our computational experiments.

Independent sampling. Fix $1\leq\tau\leq n$. A $\tau$-independent sampling is a DU sampling with

where $c_{1}=\left(\tfrac{1}{n}\right)^{\tau}$ and $c_{k}=\left(\tfrac{k}{n}\right)^{\tau}-\sum_{i=1}^{k-1}{k\choose i}c_{i}$ for $k\geq 2$.

Interpretation: There are $\tau$ processors/threads/cores available. Each processor chooses one of the $n$ blocks, uniformly at random and independently of the other processors. It turns out that the set $\hat{S}$ of blocks selected this way is DU with $q$ as given above. Since in one parallel iteration of our methods each block in $\hat{S}$ is updated exactly once, this means that if two or more processors pick the same block, all but one will be idle. On the other hand, this sampling can be generated extremely easily and in parallel! For $\tau\ll n$ this sampling is a good (and fast) approximation of the $\tau$-nice sampling. For instance, for $n=10^{3}$ and $\tau=8$ we have $q_{8}=0.9723$, $q_{7}=0.0274$, $q_{6}=0.0003$ and $q_{k}\approx 0$ for $k=1,\dots,5$.

Binomial sampling. Fix $1\leq\tau\leq n$ and $0<p_{b}\leq 1$. A $(\tau,p_{b})$-binomial sampling is defined as a DU sampling with

Notice that $\mathbf{E}[|\hat{S}|]=\tau p_{b}$ and $\mathbf{E}[|\hat{S}|^{2}]=\tau p_{b}(1+\tau p_{b}-p_{b})$.

Interpretation: Consider the following situation with independent equally unreliable processors. We have $\tau$ processors, each of which is at any given moment available with probability $p_{b}$ and busy with probability $1-p_{b}$, independently of the availability of the other processors. Hence, the number of available processors (and hence blocks that can be updated in parallel) at each iteration is a binomial random variable with parameters $\tau$ and $p_{b}$. That is, the number of available processors is equal to $k$ with probability $q_{k}$.

Case 1 (explicit selection of blocks): We learn that $k$ processors are available at the beginning of each iteration. Subsequently, we choose $k$ blocks using a $k$-nice sampling and “assign one block” to each of the $k$ available processors.

Case 2 (implicit selection of blocks): We choose $\tau$ blocks using a $\tau$-nice sampling and assign one to each of the $\tau$ processors (we do not know which will be available at the beginning of the iteration). With probability $q_{k}$, $k$ of these will send their updates. It is easy to check that the resulting effective sampling of blocks is $(\tau,p_{b})$-binomial.

Serial sampling. This is a DU sampling with $q_{1}=1$. Also, this is a NU sampling with $l=n$ and $S^{j}=\{j\}$ for $j=1,2,\dots,l$. That is, at each iteration we update a single block, uniformly at random. This was studied in .

Fully parallel sampling. This is a DU sampling with $q_{n}=1$. Also, this is a NU sampling with $l=1$ and $S^{1}={[n]}$. That is, at each iteration we update all blocks.

The following simple result says that the intersection between the class of DU and NU samplings is very thin. A sampling is called vacuous if $\mathbf{P}(\emptyset)>0$.

There are precisely two nonvacuous samplings which are both DU and NU: i) the serial sampling and ii) the fully parallel sampling.

Assume $\hat{S}$ is nonvacuous, NU and DU. Since $\hat{S}$ is nonvacuous, $\mathbf{P}(\hat{S}=\emptyset)=0$. Let $S\subset{[n]}$ be any set for which $\mathbf{P}(\hat{S}=S)>0$. If $1<|S|<n$, then there exists $S^{\prime}\neq S$ of the same cardinality as $S$ having a nonempty intersection with $S$. Since $\hat{S}$ is doubly uniform, we must have $\mathbf{P}(\hat{S}=S^{\prime})=\mathbf{P}(\hat{S}=S^{\prime})>0$. However, this contradicts the fact that $\hat{S}$ is non-overlapping. Hence, $\hat{S}$ can only generate sets of cardinalities $1$ or $n$ with positive probability, but not both. One option leads to the fully parallel sampling, the other one leads to the serial sampling. ∎

2 Technical results

For a given sampling $\hat{S}$ and $i,j\in{[n]}$ we let

The following simple result has several consequences which will be used throughout the paper.

Let $\emptyset\neq J\subset{[n]}$ and $\hat{S}$ be any sampling. If $\theta_{i}$, $i\in{[n]}$, and $\theta_{ij}$, for $(i,j)\in{[n]}\times{[n]}$ are real constants, thenSum over an empty index set will, for convenience, be defined to be zero.

Proof. We prove the first statement, proof of the remaining statements is essentially identical:

The consequences are summarized in the next theorem and the discussion that follows.

Let $\emptyset\neq J\subset{[n]}$ and $\hat{S}$ be an arbitrary sampling. Further, let $a,h\in\mathbf{R}^{N}$, $w\in\mathbf{R}^{n}_{+}$ and let $g$ be a block separable function, i.e., $g(x)=\sum_{i}g_{i}(x^{(i)})$. Then

Moreover, the matrix $P\stackrel{{\scriptstyle\text{def}}}{{=}}(p_{ij})$ is positive semidefinite.

Proof. Noting that $|J\cap\hat{S}|=\sum_{i\in J\cap\hat{S}}1$, $|J\cap\hat{S}|^{2}=(\sum_{i\in J\cap\hat{S}}1)^{2}=\sum_{i\in J\cap\hat{S}}\sum_{j\in J\cap\hat{S}}1$, $\langle a,h_{[\hat{S}]}\rangle_{w}=\sum_{i\in\hat{S}}w_{i}\langle a^{(i)},h^{(i)}\rangle$, $\|h_{[\hat{S}]}\|_{w}^{2}=\sum_{i\in\hat{S}}w_{i}\|h^{(i)}\|_{(i)}^{2}$ and

all five identities follow directly by applying Lemma 3. Finally, for any $\theta=(\theta_{1},\dots,\theta_{n})^{T}\in\mathbf{R}^{n}$,

The above results hold for arbitrary samplings. Let us specialize them, in order of decreasing generality, to uniform, doubly uniform and nice samplings.

Uniform samplings. If $\hat{S}$ is uniform, then from (28) using $J={[n]}$ we get

Plugging (33) into (28), (30), (31) and (32) yields

Doubly uniform samplings. Consider the case $n>1$; the case $n=1$ is trivial. For doubly uniform $\hat{S}$, $p_{ij}$ is constant for $i\neq j$:

Substituting (38) and (33) into (29) then gives

Nice samplings. Finally, if $\hat{S}$ is $\tau$-nice (and $\tau\neq 0$), then $\mathbf{E}[|\hat{S}|]=\tau$ and $\mathbf{E}[|\hat{S}|^{2}]=\tau^{2}$, which used in (39) gives

Moreover, assume that $\mathbf{P}(|J\cap\hat{S}|=k)\neq 0$ (this happens precisely when $0\leq k\leq|J|$ and $k\leq\tau\leq n-|J|+k$). Then for all $i\in J$,

Expected Separable Overapproximation

Recall that given $x_{k}$, in PCDM1 the next iterate is the random vector $x_{k+1}=x_{k}+h_{[\hat{S}]}$ for a particular choice of $h\in\mathbf{R}^{N}$. Further recall that in PCDM2,

again for a particular choice of $h$. While in Section 2 we mentioned how $h$ is computed, i.e., that $h$ is the minimizer of $H_{\beta,w}(x,\cdot)$ (see (17) and (18)), we did not explain why is $h$ computed this way. The reason for this is that the tools needed for this were not yet developed at that point (as we will see, some results from Section 4 are needed). In this section we give an answer to this why question.

Given $x_{k}\in\mathbf{R}^{N}$, after one step of PCDM1 performed with update $h$ we get $\mathbf{E}[F(x_{k+1})\;|\;x_{k}]=\mathbf{E}[F(x_{k}+h_{[\hat{S}]})\;|\;x_{k}]$. On the the other hand, after one step of PCDM2 we have

So, for both PCDM1 and PCDM2 the following estimate holds,

A good choice for $h$ to be used in the algorithms would be one minimizing the right hand side of inequality (42). At the same time, we would like the minimization process to be decomposable so that the updates $h^{(i)}$, $i\in\hat{S}$, could be computed in parallel. However, the problem of finding such $h$ is intractable in general even if we do not require parallelizability. Instead, we propose to construct/compute a “simple” separable overapproximation of the right-hand side of (42). Since the overapproximation will be separable, parallelizability is guaranteed; “simplicity” means that the updates $h^{(i)}$ can be computed easily (e.g., in closed form).

From now on we replace, for simplicity and w.l.o.g., the random vector $x_{k}$ by a fixed deterministic vector $x\in\mathbf{R}^{N}$. We can thus remove conditioning in (42) and instead study the quantity $\mathbf{E}[F(x+h_{[\hat{S}]})]$. Further, fix $h\in\mathbf{R}^{N}$. Note that if we can find $\beta>0$ and $w\in\mathbf{R}^{n}_{++}$ such that

we indeed find a simple separable overapproximation of $\mathbf{E}[F(x+h_{[\hat{S}]})]$:

where we recall from (18) that $H_{\beta,w}(x,h)=f(x)+\langle\nabla f(x),h\rangle+\tfrac{\beta}{2}\|h\|_{w}^{2}+\Omega(x+h)$.

That is, (44) says that the expected objective value after one parallel step of our methods, if block $i\in\hat{S}$ is updated by $h^{(i)}$, is bounded above by a convex combination of $F(x)$ and $H_{\beta,w}(x,h)$. The natural choice of $h$ is to set

Note that this is precisely the choice we make in our methods. Since $H_{\beta,w}(x,0)=F(x)$, both PCDM1 and PCDM2 are monotonic in expectation.

The above discussion leads to the following definition.

Let $\beta>0$, $w\in\mathbf{R}^{n}_{++}$ and let $\hat{S}$ be a proper uniform sampling. We say that $f:\mathbf{R}^{N}\to\mathbf{R}$ admits a $(\beta,w)$-ESO with respect to $\hat{S}$ if inequality (43) holds for all $x,h\in\mathbf{R}^{N}$. For simplicity, we write $(f,\hat{S})\sim ESO(\beta,w)$.

Inflation. If $(f,\hat{S})\sim ESO(\beta,w)$, then for $\beta^{\prime}\geq\beta$ and $w^{\prime}\geq w$, $(f,\hat{S})\sim ESO(\beta^{\prime},w^{\prime})$.

Reshuffling. Since for any $c>0$ we have $\|h\|_{cw}^{2}=c\|h\|_{w}^{2}$, one can “shuffle” constants between $\beta$ and $w$ as follows:

Strong convexity. If $(f,\hat{S})\sim ESO(\beta,w)$, then

Indeed, it suffices to take expectation in (14) with $y$ replaced by $x+h_{[\hat{S}]}$ and compare the resulting inequality with (43) (this gives $\beta\|h\|_{w}^{2}\geq\mu_{f}(w)\|h\|_{w}^{2}$, which must hold for all $h$).

Recall that Step 5 of PCDM2 was introduced so as to explicitly enforce monotonicity into the method as in some situations, as we will see in Section 7, we can only analyze a monotonic algorithm. However, sometimes even PCDM1 behaves monotonically (without enforcing this behavior externally as in PCDM2). The following definition captures this.

Assume $(f,\hat{S})\sim ESO(\beta,w)$ and let $h(x)$ be as in (45). We say that the ESO is monotonic if $F(x+(h(x))_{[\hat{S}]})\leq F(x)$, with probability 1, for all $x\in\operatorname{dom}F$.

The following theorem will be useful in deriving ESO for uniform samplings (Section 6.1) and nonoverlapping uniform samplings (Section 6.2). It will also be useful in establishing monotonicity of some ESOs (Theorems 12 and 13).

Assume $f$ is partially separable (i.e., it can be written in the form (2)). Letting $\operatorname{Supp}(h)\stackrel{{\scriptstyle\text{def}}}{{=}}\{i\in{[n]}\;:\;h^{(i)}\neq 0\}$, for all $x,h\in\mathbf{R}^{N}$ we have

Proof. Let us fix $x$ and define $\phi(h)\stackrel{{\scriptstyle\text{def}}}{{=}}f(x+h)-f(x)-\langle\nabla f(x),h\rangle$. Fixing $h$, we need to show that $\phi(h)\leq\frac{\theta}{2}\|h\|_{L}^{2}$ for $\theta=\max_{J\in\mathcal{J}}\theta^{J}$, where $\theta^{J}\stackrel{{\scriptstyle\text{def}}}{{=}}|J\cap\operatorname{Supp}(h)|$. One can define functions $\phi^{J}$ in an analogous fashion from the constituent functions $f_{J}$, which satisfy

Now, since $\phi^{J}$ depends on the intersection of $J$ and the support of its argument only, we have

The argument in the last expression can be written as a convex combination of $1+\theta^{J}$ vectors: the zero vector (with weight $\tfrac{\theta-\theta^{J}}{\theta}$) and the $\theta^{J}$ vectors $\{\theta U_{i}h^{(i)}:i\in J\cap\operatorname{Supp}(h)\}$ (with weights $\tfrac{1}{\theta}$):

Finally, we plug (53) into (52) and use convexity and some simple algebra:

Besides the usefulness of the above result in deriving ESO inequalities, it is interesting on its own for the following reasons.

Block Lipschitz continuity of $\nabla f$. The DSO inequality (48) is a generalization of (12) since (12) can be recovered from (48) by choosing $h$ with $\operatorname{Supp}(h)=\{i\}$ for $i\in{[n]}$.

Global Lipschitz continuity of $\nabla f$. The DSO inequality also says that the gradient of $f$ is Lipschitz with Lipschitz constant $\omega$ with respect to the norm $\|\cdot\|_{L}$:

Indeed, this follows from (48) via $\max_{J\in\mathcal{J}}|J\cap\operatorname{Supp}(h)|\leq\max_{J\in\mathcal{J}}|J|=\omega$. For $\omega=n$ this has been shown in ; our result for partially separable functions appears to be new.

Tightness of the global Lipschitz constant. The Lipschitz constant $\omega$ is “tight” in the following sense: there are functions for which $\omega$ cannot be replaced in (54) by any smaller number. We will show this on a simple example. Let $f(x)=\tfrac{1}{2}\|Ax\|^{2}$ with $A\in\mathbf{R}^{m\times n}$ (blocks are of size 1). Note that we can write $f(x+h)=f(x)+\langle\nabla f(x),h\rangle+\tfrac{1}{2}h^{T}A^{T}Ah$, and that $L=(L_{1},\dots,L_{n})=\operatorname{diag}(A^{T}A)$. Let $D=\operatorname{Diag}(L)$. We need to argue that there exists $A$ for which $\sigma\stackrel{{\scriptstyle\text{def}}}{{=}}\max_{h\neq 0}\tfrac{h^{T}A^{T}Ah}{\|h\|_{L}^{2}}=\omega$. Since we know that $\sigma\leq\omega$ (otherwise (54) would not hold), all we need to show is that there is $A$ and $h$ for which

Since $f(x)=\sum_{i=1}^{m}(A_{j}^{T}x)^{2}$, where $A_{j}$ is the $j$-th row of $A$, we assume that each row of $A$ has at most $\omega$ nonzeros (i.e., $f$ is partially separable of degree $\omega$). Let us pick $A$ with the following further properties: a) $A$ is a 0-1 matrix, b) all rows of $A$ have exactly $\omega$ ones, c) all columns of $A$ have exactly the same number ($k$) of ones. Immediate consequences: $L_{i}=k$ for all $i$, $D=kI_{n}$ and $\omega m=kn$. If we let $e_{m}$ be the $m\times 1$ vector of all ones and $e_{n}$ be the $n\times 1$ vector of all ones, and set $h=k^{-1/2}e_{n}$, then

establishing (55). Using similar techniques one can easily prove the following more general result: Tightness also occurs for matrices $A$ which in each row contain $\omega$ identical nonzero elements (but which can vary from row to row).

2 ESO for a convex combination of samplings

Let $\hat{S}_{1},\hat{S}_{2},\dots,\hat{S}_{m}$ be a collection of samplings and let $q\in\mathbf{R}^{m}$ be a probability vector. By $\sum_{j}q_{j}\hat{S}_{j}$ we denote the sampling $\hat{S}$ given by

This procedure allows us to build new samplings from existing ones. A natural interpretation of $\hat{S}$ is that it arises from a two stage process as follows. Generating a set via $\hat{S}$ is equivalent to first choosing $j$ with probability $q_{j}$, and then generating a set via $\hat{S}_{j}$.

Let $\hat{S}_{1},\hat{S}_{2},\dots,\hat{S}_{m}$ be arbitrary samplings, $q\in\mathbf{R}^{m}$ a probability vector and $\kappa:2^{[n]}\to\mathbf{R}$ any function mapping subsets of ${[n]}$ to reals. If we let $\hat{S}=\sum_{j}q_{j}\hat{S}_{j}$, then

$\mathbf{E}[\kappa(\hat{S})]=\sum_{j=1}^{m}q_{j}\mathbf{E}[\kappa(\hat{S}_{j})]$,

$\mathbf{E}[|\hat{S}|]=\sum_{j=1}^{m}q_{j}\mathbf{E}[|\hat{S}_{j}|]$,

$\mathbf{P}(i\in\hat{S})=\sum_{j=1}^{m}q_{j}\mathbf{P}(i\in\hat{S}_{j})$, for any $i=1,2,\dots,n$,

If $\hat{S}_{1},\dots,\hat{S}_{m}$ are uniform (resp. doubly uniform), so is $\hat{S}$.

Statement (i) follows by writing $\mathbf{E}[\kappa(\hat{S})]$ as

Statement (ii) follows from (i) by choosing $\kappa(S)=|S|$, and (iii) follows from (i) by choosing $\kappa$ as follows: $\kappa(S)=1$ if $i\in S$ and $\kappa(S)=0$ otherwise. Finally, if the samplings $\hat{S}_{j}$ are uniform, from (33) we know that $\mathbf{P}(i\in\hat{S}_{j})=\mathbf{E}[|\hat{S}_{j}|]/n$ for all $i$ and $j$. Plugging this into identity (iii) shows that $\mathbf{P}(i\in\hat{S})$ is independent of $i$, which shows that $\hat{S}$ is uniform. Now assume that $\hat{S}_{j}$ are doubly uniform. Fixing arbitrary $\tau\in\{0\}\cup{[n]}$, for every $S\subset{[n]}$ such that $|S|=\tau$ we have

As the last expression depends on $S$ via $|S|$ only, $\hat{S}$ is doubly uniform. ∎

If we fix $S\subset{[n]}$ and define $k(S^{\prime})=1$ if $S^{\prime}=S$ and $k(S^{\prime})=0$ otherwise, then statement (i) of Lemma8 reduces to (56).

All samplings arise as a combination of elementary samplings, i.e., samplings whose all weight is on one set only. Indeed, let $\hat{S}$ be an arbitrary sampling. For all subsets $S_{j}$ of ${[n]}$ define $\hat{S}_{j}$ by $\mathbf{P}(\hat{S}_{j}=S_{j})=1$ and let $q_{j}=\mathbf{P}(\hat{S}=S_{j})$. Then clearly, $\hat{S}=\sum_{j}q_{j}\hat{S}_{j}$.

All doubly uniform samplings arise as convex combinations of nice samplings.

Often it is easier to establish ESO for a simple class of samplings (e.g., nice samplings) and then use it to obtain an ESO for a more complicated class (e.g., doubly uniform samplings as they arise as convex combinations of nice samplings). The following result is helpful in this regard.

Let $\hat{S}_{1},\dots,\hat{S}_{m}$ be uniform samplings satisfying $(f,\hat{S}_{j})\sim ESO(\beta_{j},w_{j})$ and let $q\in\mathbf{R}^{m}$ be a probability vector. If $\sum_{j}q_{j}\hat{S}_{j}$ is not nil, then

First note that from part (iv) of Lemma 8 we know that $\hat{S}\stackrel{{\scriptstyle\text{def}}}{{=}}\sum_{j}q_{j}\hat{S}_{j}$ is uniform and hence it makes sense to speak about ESO involving this sampling. Next, we can write

It now remains to use (43) and part (ii) of Lemma 8:

where $w=\sum_{j}q_{j}\mathbf{E}[|\hat{S}_{j}|]\beta_{j}w_{j}$. In the third step we have also used the fact that $\mathbf{E}[|\hat{S}|]>0$ which follows from the assumption that $\hat{S}$ is not nil. ∎

3 ESO for a conic combination of functions

We now establish an ESO for a conic combination of functions each of which is already equipped with an ESO. It offers a complementary result to Theorem 9.

If $(f_{j},\hat{S})\sim ESO(\beta_{j},w_{j})$ for $j=1,\dots,m$, then for any $c_{1},\dots,c_{m}\geq 0$ we have

Proof. Letting $f=\sum_{j}c_{j}f_{j}$, we get

Expected Separable Overapproximation (ESO) of Partially Separable Functions

Here we derive ESO inequalities for partially separable smooth functions $f$ and (proper) uniform (Section 6.1), nonoverlapping uniform (Section 6.2), nice (Section 6.3) and doubly uniform (Section 6.4) samplings.

Consider an arbitrary proper sampling $\hat{S}$ and let $\nu=(\nu_{1},\dots,\nu_{n})^{T}$ be defined by

Proof. Let us use Theorem 7 with $h$ replaced by $h_{[\hat{S}]}$. Note that $\max_{J\in\mathcal{J}}|J\cap\operatorname{Supp}(h_{[\hat{S}]})|\leq\max_{J\in\mathcal{J}}|J\cap\hat{S}|\leq\min\{\omega,|\hat{S}|\}$. Taking expectations of both sides of (48) we therefore get