extfieldsofstudy:
Computer Science
Mathematics

text:
Dual optimization for convex constrained objectives without the gradient-Lipschitz assumption

The minimization of convex objectives coming from linear supervised learning problems, such as penalized generalized linear models, can be formulated as finite sums of convex functions. For such problems, a large set of stochastic first-order solvers based on the idea of variance reduction are available and combine both computational efficiency and sound theoretical guarantees (linear convergence rates). Such rates are obtained under both gradient-Lipschitz and strong convexity assumptions. Motivated by learning problems that do not meet the gradient-Lipschitz assumption, such as linear Poisson regression, we work under another smoothness assumption, and obtain a linear convergence rate for a shifted version of Stochastic Dual Coordinate Ascent (SDCA) that improves the current state-of-the-art. Our motivation for considering a solver working on the Fenchel-dual problem comes from the fact that such objectives include many linear constraints, that are easier to deal with in the dual. Our approach and theoretical findings are validated on several datasets, for Poisson regression and another objective coming from the negative log-likelihood of the Hawkes process, which is a family of models which proves extremely useful for the modeling of information propagation in social networks and causality inference.

Introduction
In the recent years, much effort has been made to minimize strongly convex finite sums with first order information. Recent developments, combining both numerical efficiency and sound theoretical guarantees, such as linear convergence rates, include SVRG [19], SAG [35], SDCA [36] or SAGA [13] to solve the following problem: where the functions ϕ i correspond to a loss computed at a sample i of the dataset, and g is a (eventually non-smooth) penalization. However, theoretical guarantees about these algorithms, such as linear rates guaranteeing a numerical complexity O(log(1/ε)) to obtain a solution ε-distant to the minimum, require both strong convexity of 1 n n i=1 ϕ i + λg and a gradient-Lipschitz property on each ϕ i , namely ϕ i (x) − ϕ i (y) ≤ L i x − y for any x, y ∈ R d , where · stands for the Euclidean norm on R d and L i > 0 is the Lipschitz constant. However, some problems, such as the linear Poisson regression, which is of practical importance in statistical image reconstruction among others (see [6] for more than a hundred references) do not meet such a smoothness assumption. Indeed, we have in this example ϕ i (w) = w x i − y i log(w x i ) for i = 1, . . . , n where x i ∈ R d are the features vectors and y i ∈ N are the labels, and where the model-weights must satisfy the linear constraints w x i > 0 for all i = 1, . . . , n.
Motivated by machine learning problems described in Section 4 below, that do not satisfy the gradient-Lipschitz assumption, we consider a more specific task relying on a new smoothness assumption. Given convex functions f i : D f → R with D f = (0, +∞) such that lim t→0 f i (t) = +∞, a vector ψ ∈ R d , features vectors x 1 , . . . , x n ∈ R d corresponding to the rows of a matrix X we consider the objective min w∈Π(X) where λ > 0, g : R d → R is a 1-strongly convex function and Π(X) is the open polytope that we assume to be non-empty. Note that the linear term ψ w can be included in the regularization g but the problem stands clearer if it is kept out.
Definition 1. We say that a function f : D f ⊂ R → R is L-log smooth, where L > 0, if it is a differentiable and strictly monotone convex function that satisfies for all x, y ∈ D f .
We detail this property and its specificities in Section 2. All along the chapter, we assume that the functions f i are L i -log smooth. Note also that the Poisson regression objective fits in this setting, where f i (x) = −y i log x is y i -log smooth and ψ = 1 n n i=1 x i . See Section 4.1 below for more details.
Related works. Standard first-order batch solvers (non stochastic) for composite convex objectives are ISTA and its accelerated version FISTA [5] and first-order stochastic solvers are mostly built on the idea of Stochastic Gradient Descent (SGD) [34]. Recently, stochastic solvers based on a combination of SGD and the Monte-Carlo technique of variance reduction [35], [36], [19], [13] turn out to be both very efficient numerically (each update has a complexity comparable to vanilla SGD) and very sound theoretically, because of strong linear convergence guarantees, that match or even improve the one of batch solvers. These algorithms involve gradient steps on the smooth part of the objective and theoretical guarantees justify such steps under the gradient-Lipschitz assumptions thanks to the descent lemma [7,Proposition A.24]. Without this assumption, such theoretical guarantees fall apart. Also, stochastic algorithms loose their numerical efficiency if their iterates are projected on the feasible set Π(X) at each iteration as Equation (2) requires. STORC [18] can deal with constrained objectives without a full projection but is restricted to compact sets of constraints which is not the case of Π(X). Then, a modified proximal gradient method from [40] provides convergence bounds relying on self-concordance [26] rather than the gradient-Lipschitz property. However, the convergence rate is guaranteed only once the iterates are close to the optimal solution and we observed in practice that this algorithm is simply not working (since it ends up using very small step-sizes) on the problems considered here. Recently, [21] has provided new descent lemmas based on relative-smoothness that hold on a wider set of functions including Poisson regression losses. This work is an extension of [4] that presented the same algorithm and detailed its application to Poisson regression losses. While this is more generic than our work, they only manage to reach sublinear convergence rates O(1/t) that applies only on positive solution (namely w * ∈ [0, ∞) d ) while we reach linear rates for any solution w * ∈ R d .
Our contribution. The first difficulty with the objective (2) is to remain in the open polytope Π(X). To deal with simpler constraints we rather perform optimization on the dual problem where for a function h, the Fenchel conjugate h * is given by h * (v) = sup u uv − h(u), and −D f * is the domain of the function x → n i=1 f * i (−x). This strategy is the one used by Stochastic Dual Coordinate Ascent (SDCA) [36]. The dual problem solutions are box-constrained to −D n f * which is much easier to maintain than the open polytope Π(X). Note that as all f i are strictly decreasing (because they are strictly monotone on (0, +∞) with lim t→0 f i (t) = +∞), their dual are defined on D f * ⊂ (−∞, 0). By design, this approach keeps the dual constraints maintained all along the iterations and the following proposition, proved in Section 7.1, ensures that the primal iterate converges to a point of Π(X). Proposition 1. Assume that the polytope Π(X) is non-empty, the functions f i are convex, differentiable, with lim t→0 f i (t) = +∞ for i = 1, . . . , n and that g is strongly convex. Then, strong duality holds, namely P (w * ) = D(α * ) and the Karush-Kuhn-Tucker conditions relate the two optima as ∀i ∈ {1, . . . , n}, α * i = −f * i (w * x i ) and w * = ∇g * 1 λn where w * ∈ Π(X) is the minimizer of P and α * ∈ −D n f * the maximizer of D.
In this chapter, we introduce the log smoothness property and its characteristics and then we derive linear convergence rates for SDCA without the gradient-Lipschitz assumption, by replacing it with log smoothness, see Definition 1. Our results provide a state-of-the-art optimization technique for the considered problem (2), with sound theoretical guarantees (see Section 3) and very strong empirical properties as illustrated on experiments conducted with several datasets for Poisson regression and Hawkes processes likelihood (see Section 5). We study also SDCA with importance sampling [42] under log smoothness and prove that it improves both theoretical guarantees and convergence rates observed in practical situations, see Sections 3.2 and 5. We provide also a heuristic initialization technique in Section 5.3 and a "mini-batch" [32] version of the algorithm in Section 5.4 that allows to end up with a particularly efficient solver for the considered problems. We motivate even further the problem considered in this chapter in Figure 1, where we consider a toy Poisson regression problem (with 2 features and 3 data points), for which L-BFGS-B typically fails while SDCA works. This illustrates the difficulty of the problem even on such an easy example.
Outline. We first introduce the log smoothness property in Section 2, relate it to self-concordance in Proposition 2 and translate it in the Fenchel conjugate space in Proposition 3. Then, we present the shifted SDCA algorithm in Section 3 and state its convergence guarantees in Theorem 4 under the log smoothness assumption. We also provide theoretical guarantees for variants of the algorithm, one using proximal operators [37] and the second using importance sampling [42] which leads to better convergence guarantees in Theorem 5. In Section 4 we focus on two specific problems, namely Poisson regression and Hawkes processes, and explain how they fit into the considered setting. Section 5 contains experiments that illustrate the benefits of our approach compared to baselines. This Section also proposes a very efficient heuristic initialization and numerical details allowing to optimize over several indices at each iteration, which is a trick to accelerate even further the algorithm.

A tighter smoothness assumption
To have a better overview of what log smoothness is, we formulate the following proposition giving an equivalent property for log smooth functions that are twice differentiable.
Proposition 2. Let f : D f ⊂ R → R be a convex strictly monotone and twice differentiable function. Then, This proposition is proved in Section 7.4 and we easily derive from it that x → −L log x on (0, +∞), x → Lx on R and x → L exp(x) on [0, +∞) are L-log smooth. This proposition is linked to the self-concordance property introduced by Nesterov [26] widely used to study losses involving logarithms. For the sake of clarity, the results will be presented for functions whose domain D f is a subset of R as this leads to lighter notations.
This property has been generalized [1,38] but always consists in controlling the third order derivative by the second order derivative, initially to bound the second order Taylor approximation used in the Newton descent algorithm [26]. While right hand sides of both properties (f (x) 2 and 2f (x) 3/2 ) might look arbitrarily chosen, in fact they reflect the motivating use case of the logarithmic barriers where f : t → − log(t) and for which the bound is reached. Hence, log smoothness is the counterpart of self-concordance but to control the second order derivative with the first order derivative. As it is similarly built, log smoothness shares the affine invariant property with self-concordance. It means that if f 1 is L log-smooth then f 2 : x → ax + b with a, b ∈ R is also L-log smooth with the same constant L. An extension to the multivariate case where D f ⊂ R d is likely feasible but is useless for our algorithm and hence beyond the scope of this paper. From the log smoothness property of a function f , we derive several characteristics for its Fenchel conjugate f * starting with the following proposition.
Proposition 3. Let f : D f ⊂ R → R be a strictly monotone convex function and f * be its twice differentiable Fenchel conjugate. Then, This proposition is proved in Section 7.5 and is the first step towards a series of convex inequalities for the Fenchel conjugate of log smooth functions. These inequalities, detailed in Section 7.6, bounds the Bregman divergence of such functions and are compared to what can be obtained with self-concordance or strong convexity (on a restricted set) in the canonical case where f : t → − log(t). It appears with log smoothness we obtain tighter bounds than what is achievable with other assumptions and that all these bounds are reached (and hence cannot be improved) in the canonical case (see Table 3).

The Shifted SDCA algorithm
The dual objective (4) cannot be written as a composite sum of convex functions as in the general objective (1), which is required for stochastic algorithms such as SRVG [19] or SAGA [13]. It is better to use a coordinate-wise approach to optimize this problem, which leads to SDCA [37], in which the starting point has been shifted by 1 λ ψ. This shift is induced by the relation linking primal and dual variables at the optimum: the second Karush-Kuhn-Tucker condition from Proposition 1, We first present the general algorithm (Algorithm 1), then its proximal alternative (Algorithm 2) and finally how importance sampling leads to better theoretical results. We assume that we know bounds (β i ) 1≤i≤n such that β i /α * i ≥ 1 for any i = 1, . . . , n, such bounds can be explicitly computed from the data in the particular cases considered in this chapter, see Section 4 for more details.
Sample i uniformly at random in {1, . . . , n} 4: The next theorem provides a linear convergence rate for Algorithm 1 where we assume that each f i is L i -log smooth (see Definition 1).
Theorem 4. Suppose that we known bounds β i ∈ −D f * such that R i = β i α * i ≥ 1 for i = 1, . . . , n and assume that all f i are L i -log smooth with differentiable Fenchel conjugates and g is 1-strongly convex. Then, Algorithm 1 satisfies where The proof of Theorem 4 is given in Section 7.7. It states that in the considered setting, SDCA achieves a linear convergence rate for the dual objective. The bounds β i are provided in Section 4 below for two particular cases: Poisson regression and likelihood Hawkes processes. Equipped with these bounds, we can compare the rate obtained in Theorem 4 with already known linear rates for SDCA under the gradient-Lipschitz assumption [36]. Indeed, we can restrict the domain of all f * i to (−β i , 0) on which Proposition 3 states that all f i are L i /(α * i 2 R 2 i )-strongly convex. Now, following carefully the proof in [36] leads to the convergence rate given in Equation (6) but with , Theorem 4 provides a faster convergence rate. The comparative gain depends on the values of ( x i 2 α * i 2 )/(λnL i ) and R i but it increases logarithmically with the value of R i . Table 1 below compares the explicit values of these linear rates on a dataset used in our experiments for Poisson regression.
Remark 1. Convergence rates for the primal objective are not provided since the primal iterate w (t) typically belongs to Π(X) only when it is close enough to the optimum. This would make most of the values of the primal objective P (w (t) ) undefined and therefore not comparable to P (w * ).

Proximal algorithm
Algorithm 1 maximizes the dual over one coordinate at Line 4 whose solution might not be explicit and requires inner steps to obtain α (t) i . But, whenever g can be written as where h is a convex, prox capable and possibly non-differentiable function, we use the same technique as Prox-SDCA [37] with a proximal lower bound that leads to see Section 7.2 for details. This leads to a proximal variant described in Algorithm 2 below, which is able to handle various regularization techniques and which has the same convergence guarantees as Algorithm 1 given in Theorem 4. Also, note that assuming that g can be written as (8) with a prox-capable function h is rather unrestrictive, since one can always add a ridge penalization term in the objective.

Importance sampling
Importance sampling consists in adapting the probabilities of choosing a sample i (which is by default done uniformly at random, see Line 3 from Algorithm 1) using the improvement which is expected by sampling it. Consider a distribution ρ on {1, . . . , n} with probabilities {ρ 1 , . . . , ρ n } such that ρ i ≥ 0 for any i and n i=1 ρ i = 1. The Shifted SDCA and Shifted Prox-SDCA with Sample i uniformly at random in {1, . . . , n} 4: w (t) ← prox h (v (t) ) 10: end for importance sampling algorithms are simply obtained by modifying the way i is sampled in Line 3 of Algorithms 1 and 2: instead of sampling uniformly at random, we sample using such a distribution ρ. The optimal sampling probability ρ is obtained in the same way as [42] and it also leads under our log smoothness assumption to a tighter convergence rate, as stated in Theorem 5 below.
. . , n and assume that all f i are L i -log smooth with differentiable Fenchel conjugates and g is 1-strongly convex. Consider σ defined by (7) and consider the distribution . . , n}. Then, Algorithm 1 and 2 where Line 3 is replaced by sampling i ∼ ρ satisfy The proof is given in Section 7.8. This convergence rate is stronger than the previous one from Theorem 4 since ( 1 Table 1 below compares the explicit values of these linear rates on a dataset used in our experiments for Poisson regression (facebook dataset). We observe that the log smooth rate with importance sampling is orders of magnitude better than the one obtained with the standard theory for SDCA which exploits only the

Applications to Poisson regression and Hawkes processes
In this Section we describe two important models that fit into the setting of this chapter. We precisely formulate them as in Equation (2) and give the explicit value of bounds β i such as α * i ≤ β i , where α * is the solution to the dual problem (4).

Linear Poisson regression
Poisson regression is widely used to model count data, namely when, in the dataset, each observation x i ∈ R d is associated an integer output y i ∈ N for i = 1, . . . , n. It aims to find a vector w ∈ R d strongly convex strongly convex with log smooth log smooth with importance sampling importance sampling Table 1: Theoretical convergence rates obtained on the facebook dataset (see Section 5.1) in four different settings: strongly convex (which is the rate obtained when all functions f i are considered with and without importance sampling [36,42] and the rate obtained in the setting considered in the chapter, with and without importance sampling. In this experiment, the maximum value for R i is 9062 and its average value is 308. As expected, the best rate is obtained by combining the log-smoothness property with importance sampling. such that for a given function φ : D φ ⊂ R → (0, +∞) + , y i is the realization of a Poisson random variable of intensity φ(w x i ). A convenient choice is to use exp for φ as it always guarantees that φ(w x i ) > 0. However, using the exponential function assumes that the covariates have a multiplicative effect that often cannot be justified. The tougher problem of linear Poisson regression, where φ(t) = t and D φ is the polytope Π(X), appears to model additive effects. For example, this applies in image reconstruction. The original image is retrieved from photons counts y i distributed as a Poisson distribution with intensity w x i , that are received while observing the image with different detectors represented by the vectors x i ∈ R d . This application has been extensively studied in the literature, see [16,4,40] and [6] for a review with a hundred references. Linear Poisson regression is also used in various fields such as survival analysis with additive effects [8] and web-marketing [9] where the intensity corresponds to an intensity of clicks on banners in web-marketing. To formalize, we consider a training dataset (x 1 , y 1 ), . . . , (x n 0 , y n 0 ) with x i ∈ R d and y i ∈ N and assume without loss of generality that y i > 0 for i ∈ {1, . . . , n} while y i = 0 for i ∈ {n + 1, . . . , n 0 } where n = #{i : y i > 0} ≤ n 0 (this simply means that we put first the samples corresponding to a label y i > 0). The negative log-likelihood of this model with a penalization function g can be written as where λ 0 > 0 corresponds to the level of penalization, with the constraint that w x i for i = 1, . . . , n. This corresponds to Equation (2) with f i (w) = −y i log(x i w) for i = 1, . . . , n, which are y i -log smooth functions, and with Note that the zero labeled observations can be safely removed from the sum and are fully encompassed in ψ. The algorithms and results proposed in Section 3 can therefore be applied for this model.

Hawkes processes
Hawkes processes are used to study cross causality that might occur in one or several events series. First, they were introduced to study earthquake propagation, the network across which the aftershocks propagate can be recovered given all tremors timestamps [30]. Then, they have been used in high frequency finance to describe market reactions to different types of orders [3]. In the recent years Hawkes processes have found many new applications including crime prediction [23] or social network information propagation [22]. A Hawkes process [17] is a multivariate point-process: it models timestamps {t i k } i≥1 of nodes i = 1, . . . , I using a multivariate counting process with a particular auto-regressive structure in its intensity. More precisely, we say that a multivariate counting process is a Hawkes process if the intensity of N i has the following structure: The µ i ≥ 0 are called baselines intensities, and correspond to the exogenous intensity of events from node i, and the functions φ ij for 1 ≤ i, j ≤ I are called kernels. They quantify the influence of past events from node j on the intensity of events from node i. The main parametric model for the kernels is the so-called exponential kernel, in which we consider 1≤i,j≤d is understood as an adjacency matrix, since entry A i,j quantifies the impact of the activity of node j on the activity of node i, while b u > 0 are memory parameters. We stack these parameters into a vector θ containing the baselines µ i and the self and cross-excitation parameters a ij u . Note that in this model the memory parameters b u are supposed to be given. The associated goodness-of-fit is the negative log-likelihood, which is given by the general theory of point processes (see [11]) as Let us define the following weights for i, j = 1, . . . I and u = 1, . . . , U , that can be computed efficiently for exponential kernels thanks to recurrence formulas (the complexity is linear with respect to the number of events of each node). Using the parametrization of the kernels from Equation (9) we can rewrite each term of the negative log-likelihood as To rewrite i in a vectorial form we define n i as the number of events of node i and the following vectors for i = 1, . . . , I: that are the model weights involved in i , and which correspond to the vector involved in the linear part of the primal objective (2) and finally for k = 1, . . . , n i which contains all the timestamps data computed in the weights computed in Equation (10). With these notations the negative log-likelihood for node i can be written as First, it shows that the negative log-likelihood can be separated into I independent sub-problems with goodness-of-fit − i (w i ) that corresponds to the intensity of node i with the weights x i,k carrying data from the events of the other nodes j. Each subproblem is a particular case of the primal objective (2), where all the labels y i are equal to 1. As a consequence, we can use the algorithms and results from Section 3 to train penalized multivariate Hawkes processes very efficiently.

Closed form solution and bounds on dual variables
In this Section with provide the explicit solution to Line 4 of Algorithm 2 when the objective corresponds to the linear Poisson regression or the Hawkes process goodness-of-fit. In Proposition 6 below we provide the closed-form solution of the local maximization step corresponding to Line 4 of Algorithm 2.
Proposition 6. For Poisson regression and Hawkes processes, Line 4 of Algorithm 2 has a closed form solution, namely This closed-form expression allows to derive a numerically very efficient training algorithm, as illustrated in Section 5 below. For these two use cases, the dual loss is given by for any α i > 0 (with y i = 1 for the Hawkes processes). For this specific dual loss, we can provide also upper bounds β i for all optimal dual variables α * i , as stated in the next Proposition.
Proposition 7. For Poisson regression and Hawkes processes, if g(w) = 1 2 w 2 and if x i x j ≥ 0 for all 1 ≤ i, j ≤ n, we have the following upper bounds on the dual variables at the optimum: The proofs of Propositions 6 and 7 are provided in Section 7.9. Note that the inner product assumption x i x j ≥ 0 from Proposition 7 is mild: it is always met for the Hawkes process with kernels given by (9) and it it met for Poisson regression whenever one applies for instance a min-max scaling on the features matrix.
Remark 2. The closed form solution from Proposition 6 is always lower than the generic bound β i , as explained in Section 7.11. Hence, we actually do not need to manually bound α

Experiments
To evaluate efficiently Shifted SDCA we have compared it with other optimization algorithms that can handle the primal problem (2) nicely, without the gradient-Lipschitz assumptions. We have discarded the modified proximal gradient method from [40] since most of the time it was diverging while computing the initial step with the Barzilai-Borwein method on the considered examples. We consider the following algorithms.
NoLips. This is a first order batch algorithm that relies on relative-smoothness [21] instead of the gradient Lipschitz assumption. Its application to linear Poisson regression has been detailed in [4] and its analysis provides convergence guarantees with a sublinear convergence rate in O(1/n). However, this method is by design limited to solutions with positive entries (namely w * ∈ [0, ∞) d ) and provides guarantees only in this case. Its theoretical step-size decreases linearly with 1/n and is too small in practice. Hence, we have tuned the step-size to get the best objective after 300 iterations.
SVRG. This is a stochastic gradient descent algorithm with variance reduction introduced in [19,41]. We used a variant introduced in [39], which uses Barzilai-Borwein in order to adapt the step-size, since gradient-Lipschitz constants are unavailable in the considered setting. We consider this version of variance reduction, since alternatives such as SAGA [13] and SAG [35] do not propose variants with Barzilai-Borwein type of step-size selection.
L-BFGS-B. L-BFGS-B is a limited-memory quasi-Newton algorithm [28,29]. It relies on an estimation of the inverse of the Hessian based on gradients differences. This technique allows L-BFGS-B to consider the curvature information leading to faster convergence than other batch first order algorithms such as ISTA and FISTA [5].
Newton algorithm. This is the standard second-order Newton algorithm which computes at each iteration the hessian of the objective to solve a linear system with it. In our experiments, the considered objectives are both log-smooth and self-concordant [26]. The self-concordant property bounds the third order derivative by the second order derivative, giving explicit control of the second order Taylor expansion [1]. This ensures supra-linear convergence guarantees and keeps all iterates in the open polytope (3) if the starting point is in it [27]. However, the computational cost of the hessian inversion makes this algorithm scale very poorly with the number of dimensions d (the size of the vectors x i ).
SDCA. This is the Shifted-SDCA algorithm, see Algorithm 2, without importance sampling. Indeed, the bounds given in Proposition 7 are not tight enough to improve convergence when used for importance sampling in the practical situations considered in this Section (despite the fact that the rates are theoretically better). A similar behavior was observed in [31].
SVRG and L-BFGS-B are almost always diverging in these experiments just like in the simple example considered in Figure 1. Hence, the problems are tuned to avoid any violation of the open polytope constraint (3), and to output comparable results between algorithms. Namely, to ensure that w x i > 0 for any iterate w, we scale the vectors x i so that they contain only non-negative entries, and the iterates of SVRG and L-BFGS-B are projected onto [0, +∞) d . This highlights two first drawbacks of these algorithms: they cannot deal with a generic feature matrix and their solutions contain only non-negative coefficients. For each run, the simply take λ = x/n where x = 1 n n i=1 x i 2 . This simple choice seemed relevant for all the considered problems.

Poisson regression
For Poisson regression we have processed our feature matrices to obtain coefficients between 0 and 1. Numerical features are transformed with a min-max scaler and categorical features are one hot encoded. We run our experiments on six datasets found on UCI dataset repository [20] and Kaggle 1 (see Table 2 for more details). These datasets are used to predict a number of interactions for a social post (news and facebook), the rating of a wine or a hotel (wine and vegas) or the number of hazards occurring in a property (property). The last one comes from simulated data which follows a Poisson regression. In Figure 2 we present the convergence speed of the five algorithms. As our algorithms follow quite different schemes, we measure this speed regarding to the computational time. In all runs, NoLips, SVRG and L-BFGS-B cannot reach the optimal solution as the problem minimizer contains negative values. This is illustrated in detail in Figure 3 for vegas dataset where it appears that all solvers obtain similar results for the positive values of w * but only Newton and SDCA algorithms are able to estimate the negatives values of w * . As expected, the Newton algorithm becomes very slow as the number of features d increases. SDCA is the only first order solver that reaches the optimal solution. It combines the best of both world, the scalability of a first order solver and the ability to reach solutions with negative entries.    Figure 4: Adjacency matrix A of a Hawkes process fitted on high-frequency financial data from the Bund market. This reproduces an experiment run in [3] where P u (resp. P d ) counts the number of upward (resp. downward) mid-price moves and T a (resp. T b ) counts the number of market orders at the best ask (resp. best bid) that do not move the price. SDCA detects inhibitive behaviors while L-BFGS-B cannot.

Hawkes processes
If the adjacency matrix A is forced to be entrywise positive, then no event type can have an inhibitive effect on another. This ability to exhibit inhibitive effect has direct implications on real life datasets especially in finance where these effects are common [2,3,33]. In Figure 4 we present the aggregated influence of the kernels obtained after training a Hawkes process on a finance dataset exploring market microstructure [3]. While L-BFGS-B (or SVRG, or NoLips) recovers only excitation in the adjacency matrix, SDCA also retrieves inhibition that one event type might have on another. It is expected that when stocks are sold (resp. bought) the price is unlikely to go up (resp. down) but this is retrieved by SDCA only. On simulated data this is even clearer and in Figure 5 we observe the same behavior when the ground truth contains inhibitive effects. Our experiment consists in two simulated Hawkes processes with 10 nodes and sum-exponential kernels with 3 decays. There are only excitation effects -all a ij u are positive -in the first case and we allow inhibitive effects in the second. Events are simulated according to these kernels that we try to recover. While it would be standard to compare the performances in terms of log-likelihood obtained on the a test sample, nothing ensures that the problem optimizer lies in the feasible set of the test set. Hence the results are compared by looking at the estimation error (RMSE) of the adjacency matrix A across iterations. Figure 5 shows that SDCA always converges faster towards its solution in both cases and that when the adjacency matrix contains inhibitive effects, SDCA obtains a better estimation error than L-BFGS-B.

Heuristic initialization
The default dual initialization in [37] (α (0) = 0 n ) is not a feasible dual point. Instead of setting arbitrarily α (0) to 1 n , we design, from three properties, a vector κ ∈ −D n f * that is linearly linked to α * and then rely on Proposition 9 to find a heuristic starting point α (0) from κ for Poisson regression and Hawkes processes.
This Proposition is proved in Section 7.12. It suggests to consider κ i ∝ 1/ x i for all i = 1, . . . n.  There is a linear link relating initial guess κ i to the dual optimum α * i on Poisson datasets but the amplitude is not adjusted yet.
for i = 1, . . . , n. Hence, α * and y are correlated (a change in y i only leads to a minor change in w * ), so we will consider κ i ∝ y i / x i .
Property 3: link with the features matrix. The inner product w * x i is positive and at the optimum, the Karush-Kuhn-Tucker Condition (5) (which links w * to x i through α * i ) tells that α * i is likely to be large if x i is poorly correlated to other features, i.e. if x i n j=1 x j is small. Finally, the choice for i = 1, . . . , n, takes these three properties into account. Figure 6 plots the optimal dual variables α * from the Poisson regression experiments of Section 5.1 against the the κ vector from Equation 11. We observe in these experiments a good correlation between the two, but κ is only a good guess for initialization α (0) up to a multiplicative factor that the following proposition aims to find. Proposition 9. For Poisson regression and Hawkes processes and g(w) = 1 2 w 2 , if we constraint the dual solution α * ∈ (0, +∞) n to be collinear with a given vector κ ∈ (0, +∞) n , i.e. α * =ᾱκ for someᾱ ∈ R, then the optimal value forᾱ is given bȳ Combined with the previous Properties, we suggest to consider as an initial point, where κ i is defined in Equation (11).
This Proposition is proved in Section 7.13. Figure 7 presents the values of α * i given its initial value α (0) i for i = 1, . . . , n and shows that the rescaling has worked properly. We validate this heuristic initialization by showing that it leads to a much faster convergence in Figure 8 below. Indeed, we observe that SDCA initialized with Equation (12) reaches optimal objective much faster than when initialization consists in setting all dual variables arbitrarily to 1.

Using mini batches
At each step t, SDCA [36] maximizes the dual objective by picking one index i t ∈ {1, . . . , n} and maximizing the dual objective over the coordinate i t of the dual vector α, and sets In some cases this maximization has a closed-form solution, such as for Poisson regression (see Proposition 6) or least-squares regression where f i (w) = (y i − w x i ) 2 leads to the explicit update In some other cases, such as logistic regression, this closed form solution cannot be exhibited and we must perform a Newton descent algorithm. Each Newton step consists in computing ∂D(α)/∂α i and ∂ 2 D(α)/∂α 2 i , which are one dimensional operations given x i 2 and w x i . Hence, in a large dimensional setting, when the observations x i have many non zero entries, the main cost of the steps resides mostly in computing x i 2 and w x i . Since x i 2 and w x i must also be computed when using a closed-form solution, using Newton steps instead of the closed-form is eventually not much more computationally expensive. So, in order to obtain a better trade-off between Newton steps and inner-products computations, we can consider more than a single index on which we maximize the dual objective. This is called the mini-batch approach, see Stochastic Dual Newton Ascent (SDNA) [32]. It consists in selecting a set I ⊂ {1, . . . , n} of p indices at each iteration t. The value of α t+1 i becomes in this case The two extreme cases are p = 1, which is the standard SDCA algorithm, and p = n for which we perform a full Newton algorithm. After computing the inner products w x i and x i x j for all (i, j) ∈ I 2 each iteration will simply performs up to 10 Newton steps in which the bottleneck is to solve a p × p linear system. This allows to better exploit curvature and obtain better convergence guarantees for gradient-Lipschitz losses [32]. We can apply this to Poisson regression and Hawkes processes where f i (x) = −y i log x. The maximization steps of Line 4 in Algorithm 2 is now performed on a set of coordinates I ⊂ {1, . . . n} and consists in finding where we denote by α I the sub-vector of α of size p containing the values of all indices in I. We initialize the vector α The gradient ∇D t I (α k I ) and the hessian ∇ 2 D t I (α k I ) have the following explicit formulas: and Note that D t I is a concave function hence −∇ 2 D t I (α k I ) will be positive semi-definite and the system in Equation (13) can be solved very efficiently with BLAS and LAPACK libraries. Let us explicit computations when p = 2. Suppose that I = {i, j} and put The gradient and the Hessian inverse are then given by This direct computation leads to even faster computations than using the dedicated libraries. We plot in Figure 9 the convergence speed for three sizes of batches 1, 2 and 10. Note that in all cases using a batch of size p = 2 is faster than standard SDCA. Also, in the last simulated experiment where d has been set on purpose to 1000, the solver using batches of size p = 10 is the fastest one. The bigger number of features d gets, the better are solvers using big batches.

About the pessimistic upper bounds
The generic upper bounds derived in Proposition 7 are general but pessimistic as they depend linearly on n. In fact this dependence is also observed in the NoLips algorithm [4] where the rate depends on a constant L = n i y i . Note that, for Nolips algorithms, L is involved in the step size definition and leads to step too small to be used in practice but that in our algorithm, these bounds  Figure 10 shows the evolution of the maximum optimal dual obtained (max i=1,...,n α * i ) for the six datasets considered in Section 5.1 for Poisson regression, on an increasing fraction of the dataset. These values are averaged over 20 samples and we provide the associated 95% confidence interval on this value. We observe that max i=1,...,n α * i has a much lower dependence in n than the bounds given by Proposition 7.

Conclusion
This work introduces the log-smoothness assumption in order to derive improved linear rates for SDCA, for objectives that do not meet the gradient-Lipschitz assumption. This provides, to the best of our knowledge, the first linear rates for a stochastic first order algorithm without the gradient-Lipschitz assumption. The experimental results also prove the efficiency of SDCA to solve such problems and its ability to deal with the open polytope constraints, improving the state-of-the-art. Finally, this work also presents several variants of SDCA and experimental heuristics to make the most of it on real world datasets. Future work could provide better linear rates under more specialized assumptions on the Gram matrix, as observed on numerical experiments. Also, to extend this work to more applications, we aim to find a generalization of the log smoothness assumption such as what [38] has done for self-concordance.

Acknowledgments
We would like to acknowledge support for this project from the Datascience Initiative of École polytechnique

Proofs
We start this Section by providing extra details on the derivation of the dual problem and the proximal version of SDCA. We then provides the proofs of all the results stated in the chapter, namely Proposition 2, Proposition 3, Theorem 4, Theorem 5, Proposition 6, Proposition 7, Remark 2, Proposition 8 and Proposition 9.

Duality and proof of Proposition 1
It is not straightforward to obtain strong duality for a convex problem with strict inequalities (such as the ones enforced by the polytope Π(X) from Equation 3). To bypass this difficulty we consider the same problem but constrained on a closed set and show how it relates to the original Problem (2). But first we formulate the two following lemmas.
Lemma 10. It exists ε > 0 such that the following problem min w∈Π |ε (X) has a solution w * that is also the solution of the original Problem (2).
Proof. First, notice that this problem has a unique solution which is unique as we minimize a convex function on a closed convex set Π(X). Also, as the function w → ψ w + λg(w) is strongly convex, since g is strongly convex, it is lower bounded. We denote by M ∈ R a lower bound of this function. Then, we consider w 0 ∈ Π(X), and choose ε sufficiently small such that for all i = 1, . . . , n, ∀t < ε, f i (t) > nP (w 0 ) − nM, this value of ε always exists since by assumption lim t→0 f i (t) = +∞ for all i = 1, . . . , n. For any w ε ∈ Π(X) \ Π ε (X) (so a w ε is such that ∃i ∈ {1, . . . , n}, w ε x i < ε), we thus have Hence, for such a value of ε, the solution to the original Problem (2) is necessarily in Π ε (X) and both the problems constrained on Π ε (X) and Π(X) share the same solution w * .
Lemma 11. For all i = 1, . . . , n, if we define by Likewise, the maximization step in the computation of f * i would share the same maximizer and f * We form the dual of the problem constrained on Π ε (X) where ε is such that Lemma 10 applies. We replace the inner products w x i by the scalars u i for i = 1, . . . , n and their equality is constrained to form the strictly equivalent problem: We maximize the Lagrangian to include the constraints. This introduces the vector of dual variables α ∈ R n as following: that leads to the corresponding dual problem: where −D f * |ε is the domain of all f * i|ε . The primal problem constrained on Π ε (X) verifies the Slater's conditions so strong duality holds and the maximizer of D |ε , α * |ε , is reached. As D |ε is concave, α * |ε is the only vector such that ∇D |ε (α * |ε ) = 0. Also, as strong duality holds, we can relate α * |ε to the primal optimum though the Karush-Kuhn-Tucker condition where w * is such that w * x i ≥ ε (see Lemma 10). Hence Lemma 11 applies and the dual formulation of the original Problem (2) that writes is such that ∇D(α * |ε ) = ∇D |ε (α * |ε ) = 0. Since D(α) is concave, this means that α * |ε = α * where α * is the solution of the dual formulation of the original Problem (2). Thus, the Karush-Kuhn-Tucker conditions that link the primal and dual optima of the problem constrained on Π |ε (X) also links the primal and dual optima of the original Problem (2). The first one is given in Equation (5) and the second one writes for any i ∈ {1, . . . n}. From the first we can define two functions linking vector w ∈ R d to α ∈ (−D f * ) n and such that w(α * ) = w * and

Proximal algorithm
Given that g * is smooth since its Fenchel conjugate is strongly convex, the gradient-Lipschitz property from Definition 4 below entails g * (v + ∆v) ≤ g * (v) + ∇g * (v) ∆v + 1 2 ∆v 2 . Hence, maximization step of Algorithm 1, namely, arg max x i − 1 λ ψ can be simplified by setting α t i such that it maximizes the lower bound (16) Setting w (t−1) = ∇g * (v (t−1) ) and discarding constants terms leads to the equivalent relation, While convergence speed is guaranteed for any 1-strongly convex g, to simplify the algorithm we will consider that g is not only 1-strongly convex but that that it can also be decomposed as where h is a prox capable function. With Proposition 12 below the relation between w t and v t becomes which is the proximal operator stated in Definition 6 below: w t = prox h (v(α t )).

Preliminaries for the proofs
Let us first recall some definitions and basic properties.
This is equivalent to This is equivalent to Definition 5. Fenchel conjugate. For a convex function f : D f → R we call Fenchel conjugate the function f * defined by Proposition 12. For a convex differentiable function f , the gradient of its differentiable Fenchel conjugate f * is the maximizing argument of (19): Proposition 13. For a convex differentiable function f and its differentiable Fenchel conjugate f * we have ∀u ∈ D f , f * (f (u)) = u and ∀v ∈ D f * , f (f * (v)) = v.

This leads to
Note that if f is γ-strongly convex (respectively L-smooth), then its Fenchel conjugate f * is 1/γ smooth (respectively 1/L strongly convex). We also recall results from [27] on self-concordant functions introduced in Definition 2. This concept is widely used to study losses involving logarithms. For the sake of clarity, the results will be presented for functions whose domain D f is a subset of R as this leads to lighter notations. Unlike smoothness and strong convexity, this property is affine invariant. From this definition, some inequalities are derived in [27]. Two of them provide lower bounds that are comparable to strong convexity inequalities: where ω(t) = t − log(1 + t), and Finally, we define the proximal operator used to apply the penalization g.
Definition 6. Proximal operator. For a convex function g : D g → R, the proximal operator associated to g is given by prox g (y) = arg min The proximal operator always exists and is uniquely defined as the minimizer of a strongly convex function. Before the proof of Theorem 4, we need to introduce new convex inequalities for L-log smooth functions. This class of function includes x → −L log x which is our function of interest in the Poisson and in the Hawkes cases.

Proof of Proposition 2
First order implies second order We start by showing that if f is a L log-smooth function then we can bound its second derivative by the square of its gradient. For any x ∈ D f , we set y = x + h in the Definition 1, which now writes Taking the limit of the previous inequality when h tends towards 0 leads to the desired inequality, Second order implies first order We now prove that if f is convex strictly monotone, twice differentiable and |f (x)| ≤ 1 L f (x) 2 then f is L-log smooth. If for all x ∈ D f , we denote by φ : x → 1 f (x) , (note that ∀x ∈ D f , f (x) = 0 as f is strictly monotone), then From this inequality, we limit the increasings of the function φ, that is the definition of a L-log smooth function for a convex strictly monotone function.

Proof of Proposition 3
We working by exhibiting several statements equivalent to log smoothness. First, we divide both sides of the log smoothness definition by f (x)f (y) > 0 since f is strictly monotone, Then we rewrite the equation in the dual space using Proposition 13, This can be rewritten into the following integrated form which becomes equivalent to the desired result with the fundamental theorem of calculus

Inequalities for log-smooth functions
The proof of SDCA [36] relies on the smoothness of the functions f i which implies strong convexity of their Fenchel conjugates f * i . Indeed, a γ strongly convex function f * satisfies the following inequality This inequality is not satisfied for L-log smooth functions. However, we can derive for such functions another inequality which can be compared to such inequalities based on self-concordance and strongly convex properties.
Let f : D f ⊂ R → R be a strictly monotone convex function and f * be its differentiable Fenchel conjugate. Then, This bound is an equality for f (x) = −L log x.
Proof. From log smoothness definition, we obtain by multiplying both sides by |f (x) − f (y)| and dividing by f (x)f (y) > 0 (since f is strictly monotone), Since f is a convex function, (f (x) − f (y))(x − y) ≥ 0 and using Proposition 13, we can rewrite the previous equivalence in the dual space: which concludes the proof.
Lemma 15. Let f : D f ⊂ R → R be a strictly monotone convex function and f * be its differentiable Fenchel conjugate. Then, This bound is an equality for f (x) = −L log x.
Proof. Let x, y ∈ D f * , we have the following on the one hand, On the other hand, applying Lemma 14 together with the fundamental theorem of calculus gives Finally, solving the integral leads to the desired result: Lemma 16. Assume that f is L-log smooth and f * is its differentiable Fenchel conjugate, then, for any y, x ∈ D f * and s ∈ [0, 1]. This bound is an equality for f (x) = −L log x.
Proof. Let x, y ∈ D f * and define for any s ∈ [0, 1], u(s) = sx + (1 − s)y. We apply Lemma 15 twice for x, u(s) and y, u(s): Combining s (23) and (1 − s)(24) leads to This Lemma which implies the barycenter u(s) = y + s(x − y) for s ∈ [0, 1] is the lower bound that we actually use in proof of Theorem 4. To compare this result with strong convexity and self-concordance assumptions, we will suppose that f * is twice differentiable and hence that Proposition 3 applies.
Comparison with self-concordant functions Instead of building our lower bounds on log smoothness, we rather exhibit what can be obtained with self-concordance combined with the lower bound f * (y) ≥ Ly −2 from Proposition 3. In this paragraph, we consider that 1 L f * is standard selfconcordant, such an hypothesis is verified for f : t → − 1 L log(t). Hence, using lower bound (21) on 1 L f and then Proposition 3, we obtain Since ∀x, y ∈ D f * , xy > 0, this lower bound is equivalent to Lemma 14 if |x| ≥ |y| but not as good otherwise. Lemma 15 can also be compared to what can be obtained applying Inequality (20) on 1 L f . Since ω : t → t − log(1 + t) is an increasing function, it leads to Again, this lower bound the same as Lemma 15 if |x| ≥ |y| but not as good otherwise. Finally, a bound equivalent to Lemma 16 for self-concordant functions is not easy to explicit in a clear form. However, it is numerically smaller than the lower bound stated in Lemma 16 for any s ∈ [0, 1] and any x, y ∈ D f * .
Comparison with strongly convex functions We cannot directly assume that f * is strongly convex as it would mean that f is gradient-Lipschitz. But, for fixed values of x and y ∈ D f * , we define on {u ∈ D f * , |u| < max(|x|, |y|)} the function f * {x,y} : u → f * (u) as the restriction of f * on this interval. The lower bound f * (y) ≥ Ly −2 from Proposition 3 implies that f * {x,y} is L/ max(x 2 , y 2 ) strongly-convex on its domain to which x and y belong. Equation (18) leads to the following inequality, valid for f * {x,y} and thus for f * , .
As soon as x = y, this lower bound is not as good as the one provided by Lemma 14. Following the same logic, we exhibit the two following lower bounds. The first one corresponds to Lemma 15 strongly convex self-concordant log smoothness Lemma 14 and is entailed by Equation (17), and the second to Lemma 16 and is entailed by Equation (22) ∀x, .
In both cases the reached lower bounds are not as tight as the ones stood in Lemmas 15 and 16. All these bounds are reported in Table 3 for an easy comparison.
Finally, two lemmas to lower bound Lemma 16 are needed as well. The analytical proof of these lemmas are very technical and not much informative. Thus, we rather illustrate them with the two following figures

Proof of Theorem 4
This proof is very similar to SDCA's proof [37] but it uses the new convex inequality on the Fenchel conjugate of log smooth functions from Lemma 16 to get a tighter inequality. We first prove the following lemma which is an equivalent of Lemma 6 from [37] but with convex functions f i that are L i -log smooth instead of being L i -gradient Lipschitz.
Lemma 19. Suppose that we known bounds β i ∈ −D f * such that R i = β i /α * i ≥ 1 for i = 1, . . . , n and assume that all f i are L i -log smooth with differentiable Fenchel conjugates and that g is 1strongly convex. Then, if α (t,i) is the value of α (t) when i is sampled at iteration t for Algorithms 1 and 2, we have for any s 1 , . . . , s n ∈ [0, 1], where Proof. At iteration t, if the dual vector is set to α (t) (and v (t) = v(α (t) ), see Equation (15)) the dual gain is where i is the index sampled at iteration t (see Line 3). For Algorithm 1, by the definition of α (t) i given on Lines 4 and 5 we have Using the smoothness inequality on g * which is 1-smooth as g is 1-strongly convex, where h(v (t−1) , ∆v) = g * (v (t−1) ) + ∇g * (v (t−1) ) ∆v + 1 2 ∆v 2 .
Hence setting ∆v to (λn) −1 (α i − α (t−1) i )x i , we can lower bound A i with For Algorithm 2, by definition of α (t) i stated at Lines 4 and 5 combined with the modified argmax relation (16), As both α (t−1) i and α * i belong to {α i ∈ −D f * , β i /α i ≥ 1}, for any s i ∈ [0, 1], the convex combination α i = (1 − s i )α (t−1) i + s i α * i belongs to it as well. Hence, for both algorithms, A i is higher than the previous quantity evaluated at this specific α i . Namely, We then use Lemma 16, in which −α * i ∈ D f * stands for x and −α (t−1) i ∈ D f * for y: This inequality is used instead of the strong convex inequality of the classic SDCA analysis [37]. If we plug this inequality into A i we obtain Hence, we retrieve B i and rewrite the previous inequality as We can sum over all possible sampled i and weight each entry with s −1 i to obtain Then since g * is convex, we obtain ≤ λn g * (v(α * )) − g * (v (t−1) ) , which can be injected in Equation (27) leading to Finally, since A i − B i = n(D(α (t,i) ) − D(α (t−1) )), we obtain D(α (t,i) ) − D(α (t−1) ) ≥ min j σ j n (D(α * ) − D(α (t−1) )), by multiplying Equation (26) with min j∈{1,...,n} σ j /n and removing the quantity G (t−1) ≥ 0. This leads to the following convergence speed after t iterations, E[D(α * ) − D(α (t) )] ≤ 1 − min j σ j n and since α t i is increasing with (α (t−1) i − λn x i 2 x i w (t−1) ), we obtain which concludes the proof.

Proof of Proposition 8
This proposition easily follows from the following computation ζ * = arg max where D is the original dual problem and c · ζ is the element wise product of the vectors c and ζ. Then, since arg max which remains valid in the multivariate case, we obtain ζ * i = α * i /c i for any i = 1, . . . , n.

Proof of Proposition 9
Using α * =ᾱκ the dual problem becomes one dimensional This problem is concave inᾱ and the optimalᾱ is obtained by setting the derivative to zero: This leads to the following second order equation which has a unique positive solution and concludes the proof.


text:
Distribution-aware Block-sparse Recovery via Convex Optimization

We study the problem of reconstructing a block-sparse signal from compressively sampled measurements. In certain applications, in addition to the inherent block-sparse structure of the signal, some prior information about the block support, i.e. blocks containing non-zero elements, might be available. Although many block-sparse recovery algorithms have been investigated in Bayesian framework, it is still unclear how to incorporate the information about the probability of occurrence into regularization-based block-sparse recovery in an optimal sense. In this work, we bridge between these fields by the aid of a new concept in conic integral geometry. Specifically, we solve a weighted optimization problem when the prior distribution about the block support is available. Moreover, we obtain the unique weights that minimize the expected required number of measurements. Our simulations on both synthetic and real data confirm that these weights considerably decrease the required sample complexity.

I. INTRODUCTION
C OMPRESSED Sensing (CS) has emerged in the past decade as a modern technique for recovering a sparse vector x P R n from compressed measurements (see [1], [2] for more explanations about this field). In this work, we consider signals x P R n that have a block-sparse structure, namely their non-zero entries appear in blocks. This property has been referred in the literature as block sparsity. It is common to use the following optimization problem for recovering the signal x from compressive measurements. P η 1,2 : min zPR n }z} 1,2 :" where A P R mˆn represents a fat measurement matrix with m ! n, tV b u q b"1 are the default disjoint blocks of size tk b u q b"1 that partition the set t1, ..., qu, y :" Ax`e P R m is the observation vector 1 , e is the noise term which is considered to be i.i.d. Gaussian with variance σ 2 , and η is an upper-bound for }e} 2 . Most of the earlier literature in block-sparse recovery is focused on the case of single constraint }Az´y} 2 ď η. However, in many applications such as DNA micro-arrays [3], [4], computational neuroscience [5] multi-band signal reconstruction, multiple measurement vector (MMV) problem S. Daei and F. Haddadi are with the School of Electrical Engineering, Iran University of Science & Technology. A. Amini is with EE department, Sharif University of Technology. 1 Our analysis holds for both real and complex-valued signals and measurements. [6], and the reconstruction of signals in union of subspaces [7], [8] [9], [10], there exist additional information (or alternatively additional constraints in P η 1,2 ) about the signal of interest. In this work, we explore the benefits of having access to extra information about the distribution of the block support (blocks containing non-zero elements) on the required number of measurements. To this end, we propose the optimization problem P η 1,2,w : min zPR n }z} 1,2,w :" where the quantities w b , b " 1, ..., q are some positive scalars, w :" rw 1 , ..., w q s T and M is some predefined model that restricts the feasible set of the solution. Specifically, we consider two new models for prior information that are of practical interest: ‚ Model 1 (Prior distribution) : We assume that the prior distribution of the block support is available. Under this setting, there are known probabilities associated with each block index b P t1, ..., qu. Namely, where bsuppp¨q returns the block support of a vector. ‚ Model 2 (Multiple Bayesian estimates): We consider L independent Bayesian estimates of bsupppxq i.e. tP i u L i"1 with accuracies We assume, without loss of generality, α i 's sum to unity. To each subset P i Ă t1, ..., qu, we assign a fixed weight ω i . In fact, it holds that w " where 1 C is the indicator function of the set C. In this work, we obtain the weights w˚P R q and ω˚P R L that minimize a threshold m 0 describing the expected number of required measurements for Models 1 and 2, respectively. Our approach is to find a suitable upper-bound for m 0 . The bound is not necessarily tight but leads to closed-form expressions for w P R q and ω P R L in Models 1 and 2, respectively. Availability of the parameters (p b 's and α i 's) associated with Models 1 and 2 is feasible in certain practical scenarios. For instance, they might be obtained from previous arXiv:1809.03005v1 [cs.IT] 9 Sep 2018 signal estimation in a dynamic scenario or from the statistic of previous training data. We further show that how Model 1 and 2 apply to brain functional magnetic resonance imaging ( fMRI) reconstruction and direction of arrival (DOA) estimation, respectively (see sections I-B,III). It often occurs in practice that the parameters p b 's and α i 's are not accurately accessible, but only an approximation is available. Our last contribution in this work is about this challenge. In fact, we show that our method of obtaining the optimal weights wå nd ω˚is robust to small changes of p b and α i as long as p b Ç 1 10 and α i Ç 1 10 . It is worth mentioning that in this work, two structures are used to construct the measurement matrix A mˆn :randomly subsampled discrete Fourier transform (DFT) 2 and i.i.d. standard normal distribution. The former lies in the structure of DOA estimation (see subsection I-B) while we use the latter for fMRI reconstruction. Although, at first look, might be different, these two distributions have a common property making them usable in our conic geometry approach [11,Theorem 1]; in fact, for both, null space of A mˆn is drawn uniformly at random from Grassmann manifold of subspaces of dimension n´m.

A. Related works and Key Differences
CS in presence of prior information has been studied in different signal models. While a large part of research (see for example [12]- [17]) deals with deterministic signal models, only a few works (see [18]- [20]) have investigated random signal models with Bayesian information. In the deterministic model, the ground truth signal has intersected with a few sets which called support estimates. The contributing level of each set to the support is available to the experimenter [12], [13]. This exact situation is investigated in [17]. They propose a non-uniform model for capturing deterministic prior information. The work [19] considers a probabilistic model where there is a continuous shape function describing the probability of contributing each element to the support. The authors obtain an upper-bound for failure probability of weighted 1 minimization. Their approach is based on calculating the internal and external angles of a weighted cross polytope. With a different approach, [20] has investigated a discrete measure for Bayesian information (a special case of Model 1 with k " 1). They relate the weights of weighted Fig. 1. Schematic diagram of DOA estimation of far-field sources. The angular half-space is divided into q " 30 angular clusters tV b u q b"1 with equal length. The associated parameters corresponding to Model 2 are L " 5,

B. Broadband DOA Estimation
Suppose that s far-field broadband signals tx i pt, f qu s i"1 in the frequency range f P rf L , f H s incident on an q-element uniform linear array (ULA). The received signal in m sensors at time t and l-th frequency bin can be expressed as: where a i pf l q :" r1, ..., e´j 2πf l pq´1qd sinpθiq c s T is the steering vector, c is the propagation velocity, d is the inter-sensor spacing and ept, f l q is a Gaussian noise term with variance σ 2 . In practice, one has to take several snapshots typt, f l qu k t"1 . This temporal redundancy is crucial in practice since the array size is limited due to physical constraints [21], [22]. Consequently, one may write Y pf l q :" ryp1, f l q, ..., ypk, f l qs " Apf l qXpf l q`Epf l q P C mˆk , where Xpf l q " rxp1, f l q, ..., xpq, f l qs P C qˆk , Apf l q " ra 1 pf l q, ..., a q pf l qs P C mˆq and Epf l q P C mˆk is defined similarly. If the sources is time-invariant over the period of snapshotting, then for all t " 1, ..., k the non-zero dominant peaks in xpt, f l q occur at the same locations corresponding to the ground truth DOAs. Hence, DOA estimation can be cast as recovering a joint sparse signal Xpf l q from Y pf l q. In addition, it is realistic for a radar engineer to know the probability of appearing the ground truth DOAs in some angular bands [23] (see Figure 1 for a schematic model of this scenario). Notation. Throughout, scalars are denoted by lowercase letters, vectors by lowercase boldface letters, and matrices by uppercase boldface letters. The ith element of a vector x is shown either by xpiq or x i . C˝denotes the polar of a cone C. We show sets (e.g. B) by calligraphic uppercase letters. B is used to represent the complement t1, ..., nuzB of a set B Ă t1, ..., nu. g P R n is a vector with standard normal distribution. The performance of an image x X P R n1ˆn2 is measured by peak signal to noise ratio (PSNR) given by: where X is the true image and }¨} 8 returns the maximum intensity.

II. MAIN RESULTS
Before stating our main results, we introduce two concepts from conic integral geometry that are used in our analysis.
Descent cone: Let x be a vector and f be its structure inducing convex function. Then, the set of decent directions forms a convex set defined by: Dpf, xq " ď tě0 tz P C n : f px`tzq ď f pxqu¨ (8) Statistical dimension: It is shown in [11] that statistical dimension of the above decent cone defined by δpDpf, xqq :" Edist 2 pg, Dpf, xq˝q, specifies the required number of measurements (i.e. m) that needs for perfect recovery. Define the expected number of measurements needed for P 1,2,w as s m :" E x δpDp}¨} 1,2,w , xqq.
Then, we call the weights that minimize s m optimal in sense of expected sample complexity. In the following propositions, we obtain closed-form solutions for optimal weights in case that x satisfies Model 1 and 2. The proofs are provided in Appendix A. Proposition 1. Let x P C n satisfies Model 1 with parameter p " rp 1 , ..., p q s T . Then, the optimal weights w˚" rw1 , ..., wq s T in P 1,2,w are obtained by solving the following equations simultaneously: Proposition 2. Let x P C n be decomposed into q blocks V b of equal size k. Assume that there exist L independent Bayesian estimates tP i u L i"1 of bsupppxq with parameter α " rα 1 , ..., α L s T . Then, the optimal weights ω˚" rω1 , ..., ωLs T in Model 2 are obtained by In applications, it is of important practical value to know how the inaccuracies of p b and α i in Model 1 and 2, respectively, affect the optimal weights. The following theorem is about this challenge. Theorem 1. Assume that p and p 1 be the true and approximate estimate of bsupppxq, respectively. Let w and w 1 be the Fig. 2. The middle image compares the success rate of P 0 1,2,w˚w ith that in P 0 1,2 when there exists a prior distribution p P R q (depicted in the left image) about the block indexes t1, ..., 50u. The optimal weight w˚is obtained by (12). The right image represents cpk b , p b q in Theorem 11 for different k b 's. In the right image, recovery with P η 1,2 is shown in the dashed line while recovery using P η 1,2,w˚i s shown in the solid line. The prior distribution p P R q that is used to calculate w˚via (12) is depicted in the left image. The used parameters are k b " 1 for all b and q " 128 2 " 16384. Selected true frames are also displayed in the right image, with green markers showing where these frames correspond to the PSNR of P η 1,2 .
optimal weights corresponding to p, and p 1 . Then, there exists a constant cpk b , p b q such that where γpa, zq :" ş 8 z u a´1 e´udu is incomplete gamma function and hp¨q : p b Ñ w b is the nonlinear function in (12).
From the above theorem and the right image of Figure 2, one can infer that the method of obtaining wb and ωi is robust to slight changes of p b and α i as long as p b Ç 1 10 and α i Ç 1 10 .

III. SIMULATIONS
In the first experiment, we construct a random blocksparse x P R 250 , whose building blocks V b have equal size k " 5. The probability of activating each block is taken from the vector p depicted in the left image of Figure 2 . Then, this signal is observed through a measurement matrix A P R mˆ250 , of which the elements are drawn from i.i.d. standard normal distribution. We obtain the optimal weights w˚corresponding to p by solving (12). In the middle image of Figure 2, we plot the success rate as a function of m. For each m, we average over 100 realizations of A and x. As turns out, P 1,2,w˚n eeds less measurements than P 1,2 . Fig. 4. Broadband DOA estimation for 10 sources using an m " 15 sensor ULA. The associated parameters are k " 10, q " 100, d " 5, c " 3ˆ10 8 , and σ " 1. The accuracies corresponding to three Bayesian estimates are α 1 " 4 5 , α 2 " 2 3 , and α 3 " 0. The top and bottom images correspond to the recovery using P η 1,2 and P η 1,2,Dω˚, respectively, where D :" r1 P 1 , 1 P 2 , 1 P 3 s and ω˚is calculated using (13).
In the second experiment, we test brain fMRI reconstruction of n 1ˆn2 images when there exists a two dimensional prior distribution (the left image in Figure 3). We utilize 27 frames of a MRI video scan, each of which contains 128ˆ128 pixels (k b " 1 @b and q " 128 2 " 16384), and consider only the gray scale components. The measurement matrix A P R 9830ˆ16384 for each frame is chosen from an i.i.d. standard normal distribution. Also, the additive noise e is chosen to be i.i.d. Gaussian with variance 1 100 . We implement two recovery programs P η 1,2 and P η 1,2,w˚( w˚is obtained by (12)) using TFOCS package [24]. From the right image of Figure 3, one can see that, for a fixed m " 9830 and with our optimal weighting strategy, P η 1,2,w˚c learly acts better than P η 1,2 (leading to PSNR values about 30dB higher than that in P η 1,2 ). In the last experiment, we test DOA estimation using broadband signals (see Subsection I-B and Figure 1). The angular half-space´r90˝, 90˝q is divided into q " 100 angular clusters (grids). For each frequency bin f l P r0, 5s GHz, we take k " 10 snapshots. We assume that there exist three (L " 3) Bayesian estimates P i with accuracies α 1 " 4 5 , α 2 " 2 3 , and α 3 " 0. Also, we use m " 15 sensors for recovering s " 10 ground truth sources and implement the optimization problems P η 1,2 and P η 1,2,w˚w here w˚" ř L i"1 ωi 1 Pi . ωi 's are obtained using (13). As turns out from Figure 4, while P η 1,2 detects many nonexisting sources with spurious DOAs, P η 1,2,w˚l ocates the ground truth sources correctly. This in turn suggests that our optimal weighting strategy considerably decreases the required number of sensors.
APPENDIX A. Proof of Propositions 1 and 2 Proof. Since δpDp}¨} 1,2,w , xqq only depends on bsupppxq and w, we show it by δpDpB, wqq. It holds that: We proceed by using an upper-bound for δpDpB, wqq a special case of which is obtained in [17,Lemma 1]. Lemma 1. The statistical dimension of descent cone of any vector x nˆ1 with bsupppxq " B satisfies: with φ B pz, kq :" Moreover, the minimum achieved at a unique t ě 0. The inequality is in fact equality in the asymptotic case (q Ñ 8) [17, Proposition 3].
Thus, by the aid of Lemma 1, it holds that By using p b :" Ppb P bsupppxqq " ř BQb q B and 1´p b " ř BSb q B , we have:

16)
Moreover, by the inequality I and (4), one may write where we used α i " PpP i Ă bsupppxqq " ř BĄPi q B and 1´α i " PpP i Ć bsupppxqq " p ř BČPi q B q. Since scaling is effectless on P 1,2,w˚, by minimizing the expressions in parentheses in (16) and (17) with respect to w b and ω i , one can get to the expressions (12) and (13).