PARTICLE STOCHASTIC DUAL COORDINATE ASCENT: EXPONENTIAL CONVERGENT ALGORITHM FOR MEAN FIELD NEURAL NETWORK OPTIMIZATION
Year: 2021

Abstract:
We introduce Particle-SDCA, a gradient-based optimization algorithm for two-layer neural networks in the mean field regime that achieves exponential convergence rate in regularized empirical risk minimization. The proposed algorithm can be regarded as an infinite dimensional extension of Stochastic Dual Coordinate Ascent (SDCA) in the probability space: we exploit the convexity of the dual problem, for which the coordinate-wise proximal gradient method can be applied. Our proposed method inherits advantages of the original SDCA, including (i) exponential convergence (with respect to the outer iteration steps), and (ii) better dependency on the sample size and condition number than the full-batch gradient method. One technical challenge in implementing the SDCA update is the intractable integral over the entire parameter space at every step. To overcome this limitation, we propose a tractable particle method that approximately solves the dual problem, and an importance re-weighting technique to reduce the computational cost. The convergence rate of our method is verified by numerical experiments.

heading: 1 INTRODUCTION
text: A major challenge in developing an optimization theory for neural network is to handle the nonconvexity of the loss landscape. Recent works showed that overparameterization is one of the key properties that makes global optimization possible. In particular, the effectiveness of overparameterization has been extensively investigated in the theory of Neural Tangent Kernel (NTK) (Jacot et al., 2018; Du et al., 2019; Allen-Zhu et al., 2019) and the mean field analysis (Nitanda & Suzuki, 2017; Mei et al., 2018; Chizat & Bach, 2018; Sirignano & Spiliopoulos, 2020).
The NTK analysis considers a linear approximation of neural networks and casts the gradient descent dynamics to that in a corresponding reproducing kernel Hilbert space (Jacot et al., 2018; Du et al., 2019; Allen-Zhu et al., 2019). Since the whole argument can be carried out essentially in a linear space, we can derive both global convergence and generalization guarantees. However, one of the biggest drawbacks of the NTK approach is that the fixed kernel fails to capture the feature training aspect of deep learning. Indeed, this feature learning ability of neural network is an important ingredient that separates deep learning from shallow models such as kernel methods (Suzuki, 2018; Yehudai & Shamir, 2019; Ghorbani et al., 2019).
In contrast, the mean field analysis considers a smaller scale of parameters than the NTK analysis, which allows the parameters to travel away from the initialization and learn features. However, the optimization problem is no longer convex and the analysis requires much more involved mathematical tools. A typical approach is to formulate the optimization of neural networks as a partial differential equation (PDE) of the distribution of parameters, for which global convergence can be established under certain conditions (Nitanda & Suzuki, 2017; Mei et al., 2018; Chizat & Bach, 2018; Rotskoff & Vanden-Eijnden, 2018; Sirignano & Spiliopoulos, 2020). A benefit of the PDE formulation is that the objective function becomes convex with respect to the probability measure of the parameters. However, the dynamics is essentially infinite dimensional, and thus most convergence results are shown in the continuous limit (infinite-width limit). Moreover, time discretization of the corresponding PDE is also challenging, so most existing works dealt with a continuous time dynamics. Indeed, Nitanda & Suzuki (2017); Chizat & Bach (2018); Rotskoff & Vanden-Eijnden (2018) proved convergence to a global optimal solution for essentially infinite width networks under
continuous time settings. On the other hand, global optimality may be shown under less restrictive settings if we consider the noisy gradient descent algorithm; the corresponding stochastic dynamics can be formulated as McKean–Vlasov dynamics, and its convergence to a global solution has been studied in Mei et al. (2018); Hu et al. (2019). However, most of these works establish convergence result for the discrete-time finite-width method by bounding the difference from its continuous-time and infinite-width limit counterpart; this usually incurs a large discretization error, which presents a limitation in its efficiency. Hence, it is important to develop a practical algorithm with optimization guarantees for finite width network and discrete time setting.
Recently, Nitanda et al. (2021) proposed a completely time-space discretized algorithm termed Particle Dual Averaging (PDA) that attains polynomial time convergence guarantee to the (regularized) global optimal solution. This method combines the dual averaging technique (Nesterov, 2005; 2009; Xiao, 2009) and the gradient Langevin dynamics (Vempala & Wibisono, 2019). More specifically, the algorithm realized time discretization by extending the dual averaging technique for finite dimensional convex optimization to the infinite dimensional setting, and realized the space discretization by sampling finite number of particles via the gradient Langevin dynamics. Eventually, the method achieves O(1/t) convergence with respect to the outer-loop iteration t. Although O(1/t) is optimal in an online optimization setting, there is still room for improvement for a finite sample size setting. Motivated by this observation, we propose a method that improves the convergence rate from the polynomial order (O(1/t)) to an exponential order (O(exp(−Ct))).
Our contribution. We propose Particle-SDCA, a novel optimization algorithm for two-layer neural networks in the mean field regime that achieves linear convergence (w.r.t. the outer loop steps) for regularized empirical risk minimization1. We adopt the stochastic dual coordinate ascent (SDCA) technique (Shalev-Shwartz & Zhang, 2013a; Takáč et al., 2013; Shalev-Shwartz & Zhang, 2013b) to achieve the linear convergence rate and much better dependency on the sample size n. Importantly, the integration with respect to the probability measure on the parameters is approximated by a finite number of particles. In contrast to existing methods, in which particles are sampled at each iteration, we employ a novel update rule of the weights on the particles so that the sampling procedure is only performed once in ñ-iterations, which provides significant computational speedup.
Furthermore, unlike existing analyses concerning the time-space discretized version of McKean–Vlasov dynamics (Mei et al., 2018; Bou-Rabee & Schuh, 2020; Bou-Rabee & Eberle, 2021), we do not couple the discrete time update with the continuous time counterpart by taking the small step size limit. Instead, we directly analyze the convergence of the discrete time update; hence, our method converges to the optimal solution with a fixed (non-vanishing) step size. Our contribution can be summarized as follows (for comparison with existing works see Table 1): • We propose a new optimization method for mean field representation of neural networks that is
efficient in the finite sample setting. By utilizing the SDCA technique, the algorithm achieves exponential convergence with respect to the outer iteration. Moreover, the global convergence is guaranteed without letting the step size converge to 0.
• The integral required in each update is approximated by an average over finite number of particles. The number of particles can be of polynomial order with respect to a required precision, and the computational cost for re-sampling can be of poly-log order.
• Thanks to a novel re-weighting scheme for particles, we only resample particles once in ñ updates, where ñ can be as large as the total sample size n; this significantly reduces the computation cost.
1We note that our goal is to propose a more efficient algorithm for optimizing mean field neural network with faster convergence guarantee, rather than to gain a better understanding on standard neural network training.

heading: 2 PRELIMINARIES
text: We first introduce the problem setting and some notations used in the following sections. We consider the empirical risk minimization problem in the supervised learning setting with a twolayer neural networks, where the input and output spaces are denoted by X ⊂ Rd and Y ⊂ R, respectively. Let hθ : X → Y be one neuron with the parameters θ ∈ Rd̃. For example, we may set hθ(x) = tanh(r)σ(w
>x) for θ = (r, w) ∈ R×Rd and an activation function σ (the tanh operation is merely to ensure boundedness of the output). In the mean field regime, we consider neural networks represented as an average of neurons: fΘ(x) = 1M ∑M m=1 hθm(x) where Θ = (θm) M m=1 ⊂ Rd̃ is a set of parameters. The continuous limit in the mean field regime is obtained by taking M →∞, and by an analogy to the law of large numbers, the function fΘ converges to the following integral representation of a two-layer neural network:
fp(x) = ∫ hθ(x)p(θ)dθ,
where p : Rd̃ → R is a probability density function on Rd̃ representing the weight of parameters. Here, we denote by P the set of probability density functions on Rd̃. As in the typical mean field analysis (Nitanda & Suzuki, 2017; Mei et al., 2018; Chizat & Bach, 2018), we aim to optimize the density function p ∈ P so that the neural network fp accurately predicts the output y ∈ Y from the input x ∈ X . To define the empirical risk and predictive risk, we let `(z, y) : Y × Y → R be a convex loss function, such as the squared loss `(z, y) = (z − y)2/2 for regression, and the logistic loss `(z, y) = log(1 + exp(−yz)) for classification. For each (xi, yi), we use the notation hi(θ) and `i(f(xi)) to indicate hθ(xi) and `(f(xi), yi) respectively. To estimate the density function p ∈ P , we minimize the regularized empirical risk defined by
min p∈P
1
n n∑ i=1 `i(fp(xi)) + λKL(p||N(0, σ2I)) (1)
where λ > 0 is a regularization parameter and KL(p,N(0, σ2I)) represents the KL divergence from the normal distributionN(0, σ2I) with mean 0 and covariance σ2I to the distribution with the density function p. This type of regularization naturally arises in training neural networks by the gradient Langevin dynamics (Hu et al., 2019; Chen et al., 2020b). In our framework, we consider directly optimizing the empirical risk with KL regularization, instead of running gradient Langevin dynamics (i.e., noisy gradient descent) on the unregularized objective. By decomposing the KL divergence in Eq. (1), we obtain the following equivalent representation of the objective:
min p∈P
P (p) := 1
n n∑ i=1 `i(fp(xi)) + λ1 ∫ p(θ)‖θ‖2dθ + λ2 ∫ p(θ) log(p(θ))dθ, (2)
where λ1, λ2 > 0 are the regularization parameters. Importantly, the objective is convex with respect to the density p ∈ P although it is non-convex with respect to each parameter θ. We make full use of this convexity in the following analysis.
The main difficulty to optimize (2) stems from the following two factors: (i) there is an integral with respect to p that should be approximated by a computationally tractable scheme, and (ii) there is a summation over n data points, which could lead to costly gradient evaluations. To overcome the former difficulty, we sample a finite number of particles from the density p. As for the later, we adopt the stochastic dual coordinate ascent (SDCA) technique, which we explain in the next section.

heading: 3 PROPOSED METHOD: PARTICLE-SDCA
text: We now explain our proposed Particle-Stochastic Dual Coordinate Ascent (P-SDCA) method.

heading: 3.1 DUAL PROBLEM AND ALGORITHM DESCRIPTION
text: The SDCA method relies on the dual problem. The dual of the optimization problem (2) can be derived from the following Fenchel’s duality theorem. Note that the convexity of P (p) with respect to the density function p is important in the derivation.
Proposition 1 (Fenchel’s duality theorem). Suppose that `i : R → R is a proper convex function, hi : Rd → R is bounded, and the primal problem satisfies infp∈P P (p) > −∞, then the following duality holds:
inf p∈P P (p) = sup g∈Rn D(g) := sup g∈Rn { − 1 n n∑ i=1 `∗i (gi)− λ2 log (∫ q[g](θ)dθ )} where q[g](θ) := exp { − 1λ2 ( 1 n ∑n i=1 hi(θ)gi + λ1‖θ‖2 )} and `∗i is the convex conjugate
2 of `i. Moreover, p∗ ∈ P and g∗ ∈ Rn are both optimal solutions of primal and dual problems if and only if fp∗(xi) ∈ ∂`∗i (g∗i ) and p[g∗] = p∗ , where p[g](θ) := q[g](θ)∫ q[g](θ′)dθ′
and ∂`∗i denotes the sub-differential of `∗i .
This proposition can be derived from Theorems 3 and 8 of Rockafellar (1967) (its finite dimensional version can be found in Corollary 31.2.1 of Rockafellar (1970)). See Appendix A for the complete proof. By the optimality condition, we can recover the optimal primal solution p∗ from the optimal dual solution g∗. Hence, we aim to optimize the dual problem as p∗ = p[g∗]. In the dual problem, updating one coordinate gi can be carried out by picking up one data points (xi, yi). The main strategy of SDCA is to randomly pick the sample index i ∈ [n] 3 and update the corresponding coordinate gi by maximizing the dual (i.e., dual coordinate ascent). Before we state the algorithm, we assume the following condition. Assumption 1.
(A1) `i : R→ R is (1/γ)-smooth4 for all i ∈ [n].
(A2) supθ∈Rd̃ |hi(θ)| ≤ 1 for all i ∈ [n].
The first assumption (A1) is satisfied, for example, by the squared loss `i(u) = (yi − u)2/2 with γ = 1 and the logistic loss `i(u) = log(1 + exp(−yu)) with γ = 1. As for the second assumption (A2), our analysis is valid for bounded neurons, and the upper-bound 1 is chosen just for simplicity and can be replaced by another value. In fact, if we assume supθ∈Rd̃ |hi(θ)| ≤ R instead, the analysis can be reduced to the setting of R = 1 by rescaling hi and the loss function, i.e., R can be absorbed into the smoothness parameter, which then becomes R/γ.
As stated above, the strategy to optimize the dual problem is to (i) randomly pick i ∈ [n] and (ii) update gi by approximately maximizing D(g) with respect to gi with other coordinates fixed. For the update of each coordinate, we employ the proximal gradient-type update where we optimize the lower bound of D(g) given by the following lemma (the proof can be found in Appendix C). Lemma 1. Under (A2) of Assumption 1, for g ∈ Rn and δg ∈ Rn, the second term of D(g) can be bounded as λ2 log (∫ q[g + δg](θ)dθ ) ≤ − 1n ∑n i=1 ∫ hi(θ)p[g](θ)dθδgi+ 1 2n2λ2 ‖δg‖21 +λ2 log (∫ q[g](θ)dθ ) .
Using this inequality, we can derive the update rule as a maximizer of an lower bound of the objective as follows. When we update only the i-th component (i.e., δgj = 0 (∀j 6= i)) for randomly chosen index i, the dual objective can be bounded as follows:
D(ǧ + δg) ≥ D(ǧ) + 1n (` ∗ i (ǧi)− `∗i (ǧi + δgi)) + 1n ∫ hi(θ)p[ǧ](θ)dθδgi − 12n2λ2 δg 2 i .
Hence, by substituting ǧ ← ǧ(t) and maximizing the right hand side with respect to δgi = gi − ǧ(t)i , we obtain the following ideal version of the update rule of ǧ(t):
ǧ (t+1) i ← argmax
gi
{ −`∗i (gi) + ∫ hi(θ)p[ǧ (t)](θ)dθ(gi − ǧ(t)i )− 1
2nλ2 (gi − ǧ(t)i ) 2
} , (3)
while the other components are not updated: ǧ(t+1)j = ǧ (t) j (j 6= i). Note that this update requires only one data point while the full gradient descent requires the whole n data points. Along with the updated dual variable ǧ(t), we update the primal variable as p̌(t) = p[ǧ(t)]. The whole optimization procedure based on this update rule is summarized in Algorithm 1.
2The convex conjugate f∗ of a convex function f : R→ R is defined by f∗(u) := supx∈R{ux− f(x)}. 3We write [N ] := {1, . . . , N} for an integer N . 4A function f : RD → R is L-smooth if and only if f is differentiable and ‖∇f(x)−∇f(y)‖ ≤ L‖x− y‖
for all x, y ∈ RD .
Algorithm 1 Dual Coordinate Descent in the continuous limit Require: training data {(xi, yi)}ni=1 and number of iterations tend
1: Choose ǧ(0)i s.t. |`∗i ′(ǧ (0) i )| ≤ 1 (i = 1, . . . , n) and `∗i (ǧ (0) i ) ≤ `∗i (0) 2: for t = 0, 1, . . . , tend do 3: Randomly choose it from {1, 2, . . . , n} 4: ǧ(t+1)it ← argmaxgit∈R { −`∗it(git) + ∫ hit(θ)p̌ (t)(θ)dθ(git − ǧ (t) it )− 12nλ2 |git − ǧ (t) it |2 } . 5: Update the primal solution as p̌(t+1) = p[ǧ(t+1)]. 6: end for 7: return g(tend)

heading: 3.2 OUR PROPOSAL: PARTICLE UPDATE METHOD
text: The update (3) requires an integration with respect to p[ǧ(t)]. We approximate the integral by a weighted sum over M particles:∫
hi(θ)p[ǧ (t)](θ)dθ ≈ ∑M m=1 rmhi(θm)∑M
m=1 rm , (4)
where (θm)Mm=1 are particles that approximately cover the support of p[ǧ (t)], and (rm)Mm=1 are the weights of each particle. Therefore, the key point of designing an efficient algorithm is to construct a “good” set of particles (θm)Mm=1 and their weights (rm) M m=1.
The construction of particle approximation consists of two stages: (1) the sampling stage, and (2) the re-weighting stage. We perform the sampling stage once every ñ ∈ [n] iterations, and perform the re-weighting stage in the other iterations. The detailed description of each stage is given below.
(1) Sampling stage. In the sampling stage, we sample M i.i.d. particles (θm)Mm=1 from the density function p[g(t)] by using a Monte Carlo sampler, and place an even weight r(t)m = 1/M (m = 1, . . . ,M) on the sampled particles (θm)Mm=1. As for the sampler, we may employ the unadjusted gradient Langevin algorithm (ULA), or the Metropolis-adjusted Langevin algorithm (MALA), both described in the following subsection. Note that the distribution of particles obtained from the Monte Carlo samplers as listed above may be biased away from the target distribution. We tolerate this inaccuracy by only requiring the particles to obey an approximated density function p(t) that is in (t)C distance from the target density p[g(t)] measured by the total variation distance5. Our convergence rate analysis will take this error into account.
(2) Re-weighting stage. Since the sampling procedure requires relatively heavy computation, it is desirable to perform this operation once in a while. We consider performing the sampling step once in ñ iterations. In between the sampling stages, we instead perform the re-weighting stage, in which we iteratively update the weight of each particle, rather than resampling them. Recall that p[g(t)](θ) ∝ q[g(t)](θ) = exp { − 1λ2 ( 1 n ∑n i=1 hi(θ)g (t) i + λ1‖θ‖2 )} which yields q[g(t+1)] = q[g(t)] exp[−hit(θ)(g (t+1) it
− g(t)it )/(nλ2)] for the index it chosen at the t-th iteration. We adapt this relation to update the weight r(t)m of particles as
r (t+1) m = r (t) m exp[−hit(θm)(g (t+1) it − g(t)it )/(nλ2)],
during the re-weighting stage because r(t)m represents “importance” of θm which is proportional to p[g(t)](θm)/p[g(ñT )](θm) so that the particle approximation (4) is (nearly) unbiased. We allow the number of particles to depend on the outer-loop iteration T , and write M (ñT ) as the number of particles in the T -th outer-loop. The entire procedure is summarized in Algorithm 2.
Finally, we output the solution of the last iterate. We call this output as Option (A) and write its solution as g(A)out . Alternatively, there is another choice Option (B), which randomly selects the stopping time t′end from the final n-iterations, i.e., that is randomly chosen from {ñTend − n + 1, . . . , ñTend}, and returns the solution g(B)out = g
(t′end). This technique improves the accuracy up to O(n) factor. Therefore, we can reduce the number of particles by O(n) times.
5The total variation (TV) distance between two probability measures with densities p, q is defined as TV(p||q) := 1
2
∫ |p− q|dθ.
Algorithm 2 Dual Coordinate Descent with the particle method Require: training data {(xi, yi)}ni=1 and numbers of inner-loop iterations ñ and outer-loop iterations
Tend, 1: Choose g(0)i s.t. |`∗i ′(g (0) i )| ≤ 1 (i = 1, . . . , n) and `∗i (g (0) i ) ≤ `∗i (0) 2: for T = 0, 1, . . . , Tend − 1 do 3: Randomly (approximately) draw i.i.d. parameters θm (m = 1, . . . ,M (ñT )) from p(ñT )(θ)dθ that satisfies TV(p(ñT )||p[g(ñT )]) ≤ (ñT )C . 4: r(ñT )m ← 1M(ñT ) (m = 1, . . . ,M
(ñT )) 5: for t = ñT, ñT + 1, . . . , ñT + ñ− 1 do 6: Randomly choose it from {1, 2, . . . , n}
7: g(t+1)it ← argmax git∈R
{ −`∗it(git) + ∑M(ñT ) m=1 r (t) m hit(θm)∑M(ñT )
m=1 r (t) m
(git − g (t) it )− 1 2nλ2 (git − g (t) it )2
} .
8: r(t+1)m ← r(t)m exp ( − 1nλ2hit(θm)(g (t+1) it − g(t)it ) )
(m = 1, . . . ,M (ñT )). 9: end for
10: end for 11: return Option (A): g(A)out = g(ñTend); Option (B): g (B) out = g
(t′end) for t′end that is randomly chosen from ñ(Tend − dnñe) + 1, . . . , ñTend}.
Sampling algorithms. Here, we present two examples that can be applied to sampling the particles (θm) M m=1 from p[g](θ) ∝ exp(−U(θ)) where U(θ) is the potential function defined as U(θ) :=
1 λ2 ( 1 n ∑n i=1 hi(θ)gi + λ1‖θ‖2 ) .
(1) ULA (unadjusted gradient Langevin algorithm): The algorithm generates samples from the following noisy gradient descent update:
θk+1 = θk − η∇U(θk) + √
2ηξk,
where η > 0 is the step-size and ξk ∼ N(0, I) is a random noise generated from the d̃-dimensional standard normal distribution. The sampling efficiency of ULA has been extensively studied in the literature (Dalalyan, 2014; Durmus & Moulines, 2017; Cheng & Bartlett, 2017; Vempala & Wibisono, 2019; Erdogdu & Hosseinzadeh, 2020).
(2) MALA (Metropolis-adjusted Langevin algorithm): The algorithm combines ULA with an acceptance-rejection step according to the Metropolis–Hastings algorithm. That is, we generate a proposal θ̃k+1 by ULA, i.e., θ̃k+1 = θk − η∇U(θk) + √ 2ηξk, and this proposal is accepted with probability α given by
α = min { 1, U(θ̃k+1)q(θ k|θ̃k+1)
U(θk)q(θ̃k+1|θk)
} ,
where q(θ′|θ) ∝ exp(−‖θ′ − θ − η∇U(θ)‖2/4η). If θ̃k+1 is not accepted, we set θ̃k+1 = θk. The convergence rate of MALA has been studied by Bou-Rabee & Hairer (2013); Ma et al. (2019); Chen et al. (2020a), to name a few. Due to the Metropolis-Hastings step, MALA allows for larger step size than ULA, which leads to an improved iteration complexity with respect to the desired accuracy C .

heading: 4 CONVERGENCE ANALYSIS
text: In this section, we give the convergence analysis of our proposed method. We first introduce the following assumptions. Assumption 2.
(A3) infx `i(x) > −∞ holds for all i ∈ [n]. There exists a constant γ′ such that for any i ∈ [n], `∗i is γ′-smooth in the interval {y ∈ dom(`∗i ) | ∂`∗i (y) ∩ [−2, 2] 6= ∅}. B1 := sup {|`i(x)− infx′ `i(x′)| | |x| ≤ 1, i ∈ [n]} and B2 := sup {|`′i(x)| | |x| ≤ 1, i ∈ [n]} are both finite.
(A4) hi is L-smooth for all i = 1, . . . , n.
It is straightforward to verify the finiteness of B1, B2 and infx `i(x) in (A3) for several practical loss functions such as the squared loss and logistic loss. The smoothness assumption (γ′) on `∗i is
required to fill in the gap between g(t) and ǧ(t), and is also satisfied by the loss functions listed above: the squared loss satisfies it with γ′ = 1 and the logistic loss satisfies it with γ′ = (e2 + 1)2/e2. The second assumption (A4) is required to prove convergence of the sampling procedure, and can be omitted if we employ a particle sampler that does not require smoothness of the density function. This assumption is satisfied, for example, by a (truncated) two-layer neuron hθ(x) = tanh(r)σ(w>x) for θ = (r, w) ∈ R× Rd (which also satisfies the boundedness assumption (A2)). We are now ready to present the convergence rate of our proposed method Algorithm 2. Here, we define s̃ def= ñγλ22(1+nγλ2) .
Theorem 1. Suppose that Assumptions 1 and 2 hold. Let tend = ñTend. There exist constants Ĉ1, Ĉ2 and Ĉ3 which can depend on s̃ and λ2 such that, if
(ñT ) C ≤ Ĉ1 exp(−s̃T/2), M (ñT ) ≥ Ĉ 2 2
( (ñT ) C )
2 log(4nTend/δ) (∀T ∈ [Tend]), (5)
for 0 < δ < 1 and λ2 ≥ 2B/n, then there exists an event E6 satisfying P (E) ≥ 1− δ in which the conditional expectation of the duality gap can be bounded by E [ P (p[g
(A) out ])−D(g (A) out ) ∣∣ E] ≤ P as long as the number of iterations tend satisfies the following conditions:
(Option A) tend ≥ 2 ( n+ 1
γλ2
) log (( n+ 1
γλ2 ) Ĉ3 P ) , (6)
(Option B) tend ≥ 2 ( n+ 1
γλ2
) log (( 1 + 1
nγλ2 ) Ĉ3 P ) + ⌈n ñ ⌉ . (7)
Moreover, the constants have the following dependency on n, ñ, λ2, s̃: Ĉ2 = exp( 4B1ñλ2n ), Ĉ −1 1 = O ( λ −1/2 2 Ĉ2 exp ( ñ(2Ĉ2+1) nγλ2+1 )) and Ĉ3 = O(max{s̃−1, λ−3/22 (1 + s̃−1/2)}).
The proof and detailed descriptions of each constant are given in Appendix C. The theorem indicates that our proposed method achieves exponential convergence with respect to the number of outer iterations. This coincides with the convergence rate of the finite dimensional version of SDCA (Shalev-Shwartz & Zhang, 2013a). In contrast, a (non-stochastic) full-batch gradient descent requires O((n/λ2γ) log(1/ P )) total gradient evaluations to optimize a convex objective with a condition number 1/(λ2γ). Importantly, SDCA (and our method) “decouples” the sample size n and the condition number 1/(λ2γ) yielding only O((n+ 1/λ2γ) log(1/ P )) iterations, which is beneficial in a setting of big data and a small regularization parameter.
The difference between options (A) and (B) is only an O(log(n)) factor in terms of the number of iterations. However, this improvement affects the number of particles given by Eq. (5). Indeed, the log(n) factor improvement of tend allows us to take M (ñTend) ≥ Ω(Ĉ22 Ĉ3/(Ĉ21 P )) in the last inner loop, whereas option (A) requires M (ñTend) ≥ Ω(nĈ22 Ĉ3/(Ĉ21 P )). Either way, the number of particles can be of polynomial order in terms of the target precision P .
Remark 1. Depending on the regularization strength λ2, the re-sampling interval ñ should be chosen differently. For constant order (non-vanishing) λ2, we may take the re-sampling interval as ñ = n so that all constants are bounded. On the other hand, in the small λ2 regime, say λ2 = O(1/n), we have Ĉ2 = exp(O(ñ)), and there also appears an exponential dependency on Ĉ1 w.r.t. ñ. In this case, ñ = O(1) is a better choice to avoid the large particle complexity. But for such choice of ñ, we should re-sample at almost every iteration, which is computationally demanding; this presents a trade-off. In practice, re-sampling can be executed in parallel and thus efficient training (using GPU) is possible. Importantly, to obtain the usual O(1/ √ n)-generalization error, λ2 = 1/ √ n is sufficient
(Nitanda et al., 2021) and thus we may take ñ = Ω( √ n) in that situation. We however conjecture that this trade-off is unavoidable without additional assumptions.
Sampling complexity. In Theorem 1, it is assumed that the density p(ñt) from which the particles are generated is within (ñT )C from the target density p[g (ñT )] measured by the total variation distance.
6Intuitively, the event E represents an event in which “good” particles are sampled from the sampler to approximate the integral.
The next proposition provides an upper bound of the number of iterations required to obtain the (ñT )C accurate distribution by the Monte-Carlo samplers such as ULA and MALA. The key property of the target density that determines the convergence rate of these methods is the log-Sobolev inequality with a constant cLS (see Definition 2 for its precise definition). Proposition 2. Under Assumptions 1 and 2, the target density p[g(ñT )] at each outer iteration satisfies the log-Sobolev inequality with constant cLS = 2λ1λ2 exp(−4B2/λ2). Accordingly, the number of iterations K(T )∗ to obtain the (ñT ) C -accurate distribution can be evaluated as follows:
(ULA) K(T )∗ ≤ 4L 2
c2LS max
{ 1, 4d̃
( (ñT ) C ) 2
} log [ 4B2/(λ2 (ñT ) C ) ] ,
(MALA) K(T )∗ ≤ O [ 1
c 5/2 LS
( B2(1+L)+λ1
λ2 + log(1/ (ñT ) C )
)3/2 d̃ ] .
The proof and a more detailed statement can be found in Appendix B. Roughly speaking, this proposition states that the number of iterations for ULA can be simplified as K(T )∗ ≥ Ω ( (cLS (ñT ) C ) −2 log(n) ) . On the other hand, MALA requires only poly-log order complexity with respect to (ñT )C as K (T ) ∗ = Ω ( c −5/2 LS log(1/ (ñT ) C ) 3/2 ) , which is therefore more preferable to obtain higher accuracy. Combining the required number of particles M (ñT ) and the number of iterations K (T ) ∗ for the sampler, we see that the total computational complexity is still of polynomial order.
Remark 2. Note that the log-Sobolev constant cLS also depends exponentially on 1/λ2. While the exponential factor can be avoided in the outer loop (Theorem 1) by appropriately setting ñ, such dependence is known to be unavoidable for cLS in the most general setting (Menz & Schlichting, 2014). An interesting direction is to explore additional assumptions that remove this dependence.
Total complexity. By summarizing the results we obtained above, the total computational complexity of our algorithm (Algorithm 2) with the MALA sampler can be roughly bounded as∑Tend
T=1M (ñT )(1 + nK (T ) ∗ /ñ)ñ = O ( n log(n) ñ P ( n+ 1λ2γ )[(( 1 + 1λ2 ) log( −1P ) cLS )5/2]) ,
where we assumed Ĉ−11 , Ĉ2 = O(1) (see Remark 1 on how this is achieved). Detailed evaluation of the computational complexity with the ULA sampler can be found in Appendix B.2.

heading: 5 ADDITIONAL RELATED WORKS
text: Hu et al. (2019); Sirignano & Spiliopoulos (2020) showed a linear convergence rate of the McKeanVlasov dynamics. However, the theory only handles the infinite width limit and continuous time dynamics, and it is from trivial to covert their results to a discrete time and finite width setting. Mei et al. (2018) also analyzed space-time discretization of a similar dynamics. However, a quantitative convergence rate is not provided, and the time discretization error grows exponentially with the time horizon. Compared to these results, we deal with general (1/γ)-smooth loss functions and made the dependency on the parameters λ2 and γ explicit in our convergence rate. Bou-Rabee & Schuh (2020); Bou-Rabee & Eberle (2021) showed an exponential convergence of gradient Langevin dynamics with an interaction potential, which includes certain mean field models. When specialized to mean field neural networks, our formulation covers a much wider range of objectives and loss functions. Furthermore, they considered only the Hamiltonian Monte Carlo sampler, and the analysis relied on the coupling between the discrete time dynamics and the continuous time counterpart. Consequently, their outer-loop iteration complexity is larger than O(log(1/ P )/ P ), whereas our bound requires only poly-log order cost due to the flexible choice of sampling algorithm. Jabir et al. (2021) also considered spatial-time discretization of mean field dynamics. However, it contains the time discretization lag h in their convergence rate, while our method is purely discrete time and is not affected by such a time discretization error.

heading: 6 NUMERICAL EXPERIMENTS
text: We verify our convergence rate analysis on regression problems with the squared loss. We consider teacher-student setup in which teacher and student are both two-layer neural networks of different
width. The training dataset is generated by the teacher model:
yi = σ(w ∗>xi + b ∗) + i,
where xi and i are independently drawn from Gaussian distributions. We optimize a tanh neural network (hθ(x) = tanh(w>x+ b) with θ = (w, b)) of width Ms as a student using Algorithm 2.
Convergence rate of P-SDCA. We verify the linear convergence rate of P-SDCA as suggested by our theoretical analysis in a simple setting. To verify the convergence rate, we run the method with 1000-steps of ULA. As for the teacher model, we use the ReLU activation for σ and (note that even though this is a simple setting, the teacher model is “not” included in the student model). The number of examples and dimensionality of the input space are n = 500 and d = 50, respectively. The other hyperparameters are set to λ1 = 10−4, λ2 = 10−4, and Ms = 200. The left side of Figure 1 depicts the convergence of the primal objective to the optimal value. The x-axis is the number of outer-iteration and the y-axis is the excess primal objective P (·)− infp P (p) including the regularization terms on a training set and a test set in log-scale. The integral with respect to the primal solution was computed by the particle average and the entropy is estimated by using the k-NN entropy estimator (Kozachenko & Leonenko, 1987; Brodersen, 2020). Since the y-axis is log-scaled, we can clearly observe the linear convergence of P-SDCA.
Comparison of predictive accuracy with SGD and PDA. We compare the empirical performance of P-SDCA against SGD and PDA in terms of predictive accuracy. In this experiment, we employ a sign activation for the teacher model (σ(·) = sign(·)), and the student is an overparameterized tanh network of width Ms = 200. The number of examples and input dimensionality are n = 1000 and d = 50. The regularization parameters are set to λ1 = 10−3, λ2 = 10−4 or 10−5. As for SGD, we cannot apply the entropic regularization and thus the term is omitted in the optimization. We also used the neural tangent kernel (NTK) scaling (fΘ(x) = 1√M ∑M m=1 hθm(x)) for SGD optimization in addition to the mean field scaling. The test error with respect to the number of gradient evaluations is displayed in the right side of Figure 1. We make the following observations: (i) P-SDCA achieves a better test error than other methods; (ii) SGD achieves minimal test error at the early stage but does not perform well afterwards, which is due to the absence of the entropic regularization; (iii) the mean field parameterization outperforms the NTK parameterization in terms of test loss.

heading: 7 CONCLUSION
text: In this paper, we proposed a new stochastic optimization technique for two-layer neural networks with the mean field representation. Our algorithm can be seen as an infinite dimensional extension of SDCA that maintains its advantage, namely, that it converges linearly with respect to the outeriteration, and is much more efficient in terms of the sample size than non-stochastic methods. In the algorithm, we approximate the integral required in each update by an average over finite number of particles. Remarkably, the number of particles can be of polynomial order with respect to the required precision, and its sampling cost can be of poly-log order. This is significantly different from existing numerical approximation approaches of McKean-Vlasov equations, which typically require either vanishing step size or exponentially large particle size. We validated the linear convergence of P-SDCA in numerical experiments, and demonstrated that it outperforms other existing methods.

heading: 1 Introduction 1

heading: 2 Preliminaries 3

heading: 3 Proposed method: Particle-SDCA 3
text: 3.1 Dual problem and algorithm description . . . . . . . . . . . . . . . . . . . . . . . 3
3.2 Our proposal: Particle update method . . . . . . . . . . . . . . . . . . . . . . . . 5

heading: 4 Convergence analysis 6

heading: 5 Additional related works 8

heading: 6 Numerical experiments 8

heading: 7 Conclusion 9

heading: A Fenchel’s duality theorem (Proof of Proposition 1) 14

heading: B Computational complexity of sampling algorithms (Proof of Proposition 2) 15
text: B.1 Log-Sobolev inequality of the target function . . . . . . . . . . . . . . . . . . . . 15
B.2 Sampling complexity of ULA . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
B.3 Sampling complexity of MALA . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

heading: C Proof of Theorem 1 19
text: C.1 Convergence of the ideal update ǧ(t) (Algorithm 1) . . . . . . . . . . . . . . . . . 21
C.2 Convergence of the particle-version update g(t) (Algorithm 2) . . . . . . . . . . . 25
C.2.1 Convergence of auxiliary variables g̃(t) and ḡ(t) . . . . . . . . . . . . . . . 29
C.2.2 Particle approximation, and error analysis of g(t) . . . . . . . . . . . . . . 30
C.2.3 Integration of two sources of discretization error . . . . . . . . . . . . . . 34
C.2.4 Discrete Gronwall’s lemma . . . . . . . . . . . . . . . . . . . . . . . . . . 40
D Additional experiments and details of experimental settings 40
——Appendix——

heading: A FENCHEL’S DUALITY THEOREM (PROOF OF PROPOSITION 1)
text: To prove Proposition 1, we utilize the infinite dimensional Fenchel’s duality theorem by Rockafellar (1967). Two real vector spaces E and E∗ with locally convex Hausdorff topology are called topologically paired vector spaces if all the elements of each space can be identified with all the continuous linear functionals on the other. Therefore, for a continuous linear functional y : X → R, there exists y ∈ E∗ such that we can write y(x) = 〈y, x〉E∗,E as a bilinear form 〈·, ·〉E∗,E . Definition 1 (Fenchel-Legendre convex conjugate). Let X,X∗ be topologically paired vector spaces. For f : X → R ∪ {∞}, its Fenchel-Legendre convex conjugate f∗ : X∗ → R ∪ {∞} is defined by
f∗(y) = sup{〈y, x〉X∗,X − f(x) | x ∈ X}.
Remark 3. When X = X∗ = Rd equipped with the usual inner product, the definition of the Fenchel-Legendre convex conjugate coincides with the definition introduced in the main text. Lemma 2 (Theorems 3 and 8 of Rockafellar (1967)). Let (X,X∗) and (Y, Y ∗) be two pairs of topologically paired vector spaces. Let f : X → R ∪ {∞} and g : Y → R ∪ {∞} are proper lower semi-continuous convex functions, and let A : X → Y be a bounded linear operator. Assume that ∃x ∈ X such that f(x) is finite and g is finite and continuous at Ax. Then,
inf{g(Ax) + f(x) | x ∈ X} = sup{−g∗(y)− f∗(−A∗y) | y ∈ Y ∗},
where A∗ is the adjoint operator of A such that 〈y,Ax〉Y ∗,Y = 〈A∗y, x〉X∗,X . Moreover, x̂ and ŷ are the optimal solutions of the primal and dual problems respectively if and only if Ax̂ ∈ ∂g∗(ŷ) and A∗ŷ ∈ ∂f(x̂).
Now, we are ready to prove Proposition 1.
Proof of Proposition 1. We utilize Lemma 2 to show the assertion. We let Y = Rn. Then, we know that Y is self-dual, that is, its topological dual Y ∗ is Y itself where the linear functional is given by the standard inner product 〈u, v〉Y ∗,Y = ∑n i=1 uivi for u = (ui) n i=1 ∈ Y and v = (vi) n i=1 ∈ Y ∗. We let X be the set of L1-integral functions with respect to the Lebesgue measure, i.e., X = L1(Rd̃) = L1(Rd̃;µ), where µ is the Lebesgue measure on (Rd̃,B(Rd̃)). It is known that its topological dual is the set of essentially bounded functions equipped with L∞-norm, i.e., X∗ = L∞(Rd̃)(= L∞(Rd̃;µ)) = {f : Rd̃ → R | f is B(Rd̃)-measurable, ‖f‖∞ < ∞}, where ‖f‖∞ = inf{M > 0 | µ({x | |f(x)| ≥ M}) = 0} for the Lebesgue measure µ. The bilinear functional 〈·, ·〉X∗,X is given by 〈g, f〉X∗,X = ∫ fgdµ. Then, X,X∗ become topologically paired vector spaces by equipping X∗ with its weak∗ topology with respect to the inner product 〈·, ·〉X∗,X .
We let A : L1(Rd̃) → Rn be Af = ( ∫ hi(θ)f(θ)dθ) n i=1 ∈ Rn for f ∈ L1(Rd̃). Then, we can see that A is a bounded linear functional by the boundedness of hi(·). It is also known that the KL divergence is jointly lower semi-continuous with respect to the weak convergence topology (Posner, 1975, Theorem 1). Since the L1-norm between probability density functions is equivalent to the total variation norm between the corresponding probability distributions, and the total variation norm gives a stronger topology than the weak topology, KL divergence is lower semi-continuous with respect to the L1-norm of the corresponding density functions. It is also well known that KL divergence is jointly convex and it is easy to check that p ∈ P 7→ KL(p,N(0, σ2I)) is proper, i.e., ∃p ∈ P such that KL(p,N(0, σ2I)) <∞.
Under these settings, if we set g(u) = 1n ∑n i=1 `i(ui) for u = (ui) n i=1 ∈ Y = Rd̃ and
f(p) = { λ2KL(p||N(0, λ1/λ2I)) (p ∈ P), ∞ (otherwise),
for p ∈ L1(Rd̃), then we can see that P (p) = g(Ap) + f(p) for p ∈ P and infp∈P P (p) = infh∈X{g(Ah) + f(h)}. We can also verify the condition in Lemma 2. Indeed, if we let p0 ∈ P be the density function of N(0, λ1/λ2I), then f(p0) = 0 < ∞ and g is continuous at Ap0 because a
convex function is continuous when its effective domain is the whole space Rn (Rockafellar, 1976, Theorem 10.1), and the domain of g is the whole space Rn (remember that the range of `i is R).
Therefore, by applying Lemma 2 and going through a simple calculation to obtain f∗, we obtain the desired assertion.

heading: B COMPUTATIONAL COMPLEXITY OF SAMPLING ALGORITHMS (PROOF OF
text: PROPOSITION 2)
B.1 LOG-SOBOLEV INEQUALITY OF THE TARGET FUNCTION
In each outer iteration, we sample M particles from the density p[g(ñT )] so that the particles are sampled from an approximated density p(ñT ) with TV(p(ñT )||p[g(ñT )]) ≤ (T )S . To estimate the convergence of Monte-Carlo methods such as ULA and MALA, the log-Sobolev inequality associated with the target density plays the central role.
Definition 2 (Log-Sobolev inequality). Let p(θ) be a smooth probability density function on Rd̃. Then, we say p(θ) (or its corresponding probability measure on (R,B(Rd̃))) satisfies the log-Sobolev inequality (LSI) with a constant cLS > 0 if and only if, for any smooth function φ : Rp → R with Ep[φ2] <∞, it holds that
Ep[φ2 log(φ2)]− Ep[φ2] log(Ep[φ2]) ≤ 2
cLS Ep[‖∇φ‖22].
The following lemma shows that the target density p[g(ñT )] satisfies the LSI.
Lemma 3. Suppose maxi |g(ñT )i | ≤ B (this is true by Lemma 6 below with B = B2) and assume (A2) of Assumption 1, then p[g(ñT )] satisfies the LSI with the constant
cLS = 2λ1 λ2 exp(−4B/λ2).
Proof. Remember that p[g(ñT )] is given as
p[g(ñT )](θ) ∝ exp ( − 1 λ2 1 n n∑ i=1 hi(θ)g (ñT ) i − 1 2 2λ1 λ2 ‖θ‖2 ) .
Therefore, p[g(ñT )](θ) can be written as
p[g(ñT )](θ) = exp(f(θ))q(θ)
Eq[exp(f(·))] , (8)
where f(θ) = − 1λ2 1 n ∑n i=1 hi(θ)g (ñT ) i and q(θ) is the probability density function of a Gaussian measure N(0, λ22λ1 I). Note that ‖f‖∞ ≤ B/λ2 by the assumption (A2) and maxi |g (ñT ) i | ≤ B. Then, it is well known that q satisfies the LSI with a LSI constant cLS(q) = 2λ1λ2 and Holley & Stroock (1987) showed that the bounded perturbation as in Eq. (8) preserves the LSI condition where the LSI constant of p[g(ñT )] is bounded by cLS = cLS(q)/ exp(4‖f‖∞) ≥ 2λ1λ2 exp(−4B/λ2). This concludes the proof.
By Lemma 3, we show that existing convergence rate analyses of ULA and MALA such as Vempala & Wibisono (2019); Ma et al. (2019) can be applied to our target density p[g(ñT )]. From now on, we denote our target density p[g(ñT )] by p∗, and the negative log density of p[g(ñT )] by U , that is, p∗(θ) = exp(−U(θ)). In addition, let pk be the marginal distribution of θk (sample in the k-th iteration of ULA or MALA).
B.2 SAMPLING COMPLEXITY OF ULA
First, we apply the convergence rate analysis of ULA developed by Vempala & Wibisono (2019) to our setting. We utilize the following proposition given by Vempala & Wibisono (2019). Proposition 3 (Theorem 1 of Vempala & Wibisono (2019)). Assume that p∗(θ)dθ satisfies the log-Sobolev inequality with a constant cLS and U is L-smooth. Under these assumptions, for the step size η satisfying η ≤ cLS4L2 , we have that
KL (pk ‖ p∗) ≤ exp(−cLSηk)KL (p1 ‖ p∗) + 8ηd̃L2
cLS .
Proposition 3 implies that if the step size satisfies η ≤ cLS4L2 min { 1, 4d̃ } and the number of iterations satisfies
k ≥ 1 cLSη
log 2KL (p0 ‖ p∗) ,
then we obtain an -accurate distribution pk in terms of the KL-divergence.
If we employ p0(θ) ∝ exp(−λ1λ2 ‖θ‖ 2) (a Gaussian measure), then 2KL (p0 ‖ p∗) is bounded by 4B2λ2 , from Lemma 6:
KL (p0 ‖ p∗) = ∫ p0(θ) log p0(θ)
p∗(θ) dθ =− ∫ p0(θ) [log p ∗(θ)− log p0(θ)] dθ =− ∫ p0(θ) log exp ( − ∑n i=1 hi(θ)gi−λ1‖θ‖ 2 2 λ2 ) ∫
exp ( − ∑n i=1 hi(θ
′)gi−λ1‖θ′‖22 λ2
) dθ′ − log p0(θ) dθ ≤− ∫ p0(θ) log exp ( − supi |gi|λ2 ) exp ( −λ1‖θ‖ 2 2 λ2 ) exp ( supi |gi| λ2 ) ∫ exp ( −λ1‖θ‖ 2 2 λ2 ) dθ − log p0(θ)
 dθ ≤− ∫ p0(θ) [ log { exp ( −2 ·B2
λ2
) · p0(θ) } − log p0(θ) ] dθ
≤2B2 λ2
Moreover, using Pinsker’s inequality bounds total variation norm by KL-divergence: ‖pk − p∗‖TV ≤ √ 2KL (pk ‖ p∗).
Therefore, in order to have an C -accurate distribution in total variation, it suffices to let the step size satisfy η ≤ cLS4L2 min { 1, 2C 8d̃ } and the number of iterations satisfy
k ≥ 2 cLSη log 4B2 λ2 C .
Hence, the total computational complexity can be bounded by O ( d̃L2
c2LS 2 C
log B2λ2 C
) .
Total computational complexity of sampling by ULA According to Theorem 1, sampling accuracy in total variation (ñT )C should be smaller than O(exp ( − s̃T2 ) ) in each sampling step T . Therefore, the total computational complexity for sampling can be calculated as follows:
Algorithm 2 with Option A: Suppose we run outer-loop up to Tend ≥ 1 s̃ log (( n+ 1γλ2 ) (B1+B2+ C0 1−exp(−s̃) ) P ) times. Indeed it is sufficient to have an P -
approximated solution, according to Theorem 3, which gives a tighter and more detailed
bound than in Theorem 1. Then, the total required computational complexity for sampling is
Tend∑ T=1 O
( d̃L2
c2LS( (ñT ) C )
2 log
B2
λ2 (ñT ) C
)
≤O
( d̃L2
Ĉ21c 2 LS
log B2
λ2 (ñTend) C ) Tend∑ T=1 exp (s̃T )
=O
( d̃L2
Ĉ21c 2 LS
log B2
λ2 (ñTend) C
) exp (s̃Tend)
1− e−s̃
=Õ  ( n+ 1γλ2 )2 ñ P d̃L2 c2LS log B2(B1 +B2 + C0) λ2 + 1 nγ2  , where Õ hides log log terms.
Algorithm 2 with Option B: According to Theorem 1, it is sufficient to run outer-loop up to Tend ≥ 1 s̃ log (( 1 + 1nγλ2 ) Ĉ3 P ) times to acquire an P -approximated solution. You can see the
detailed definition of Ĉ3 = O(max{s̃−1, λ−3/22 s̃−1/2}) in Eqs. (52a)-(52c) of Theorem 3 and in Proof of Theorem 1. Therefore, the total required computational complexity for sampling is
Tend∑ T=1 O  d̃L2 c2LS ( (ñT ) C )2 log B2 λ2 (ñT ) C  exp (s̃T ) ≤O ( d̃L2
Ĉ21c 2 LS
log B2
λ2 (ñTend) C ) Tend∑ T=1 exp (s̃T )
=O
( d̃L2
Ĉ21c 2 LS
log B2
λ2 (ñTend) C
) exp (s̃Tend)
1− e−s̃
=Õ max{λ−3/22 ,(nγλ2+1ñγλ2 )1/2} n3/2 ( 1 + 1nγλ2 )5/2 ñ3/2 P d̃L2 c2LS log B2(B1 +B2 + C0) λ2 + 1 nγ2  , where Õ hides log log terms as well.
B.3 SAMPLING COMPLEXITY OF MALA
Here, we adopt the convergence rate analysis of MALA obtained by Ma et al. (2019). They showed the convergence under the setting where U is µ-strongly convex outside of a compact domain and is L-smooth. Although our situation does not necessarily satisfy this condition, the proof of Ma et al. (2019) essentially utilizes it to show the LSI of the target density. Fortunately, we have already seen that the target density p∗ satisfies the LSI in Lemma 3. Thus, the analysis of Ma et al. (2019) can be applied to our setting.
From now on, we denote p[g(ñT )] by p∗. Let R(s) be a positive real such that P ∗(‖θ‖ ≥ R(s)) ≤ s for 0 < s < 1 where P ∗ is the probability measure corresponding to p∗ and L̃ > 0 be the smoothness parameter of log(p∗): ‖∇ log(p∗(θ))−∇ log(p∗(θ′))‖ ≤ L̃‖θ− θ′‖. Then, Ma et al. (2019) showed the following proposition.
Proposition 4 (Lemma 9 of Ma et al. (2019) with modification). Let the initial density p0 satisfy
sup θ∈Rp
p0(θ) p∗(θ) ≤ β
for β > 0. Then, for the step size η > 0 satisfying
η = O ( min { 1
L̃3/2R( /2β) ,
1
L̃d̃ ,
1
L̃4/3R( /2β)2/3d̃1/3
}) ,
and the number k of iterations satisfying
k = Ω
( 1
c2LSη log
( 2β )) ≥ Ω ( 1
c2LS log
( 2β ) max{L̃3/2R( /2β), L̃d̃, L̃4/3R( /2β)2/3d̃1/3} ) , (9)
it holds that ‖p∗ − pk‖TV ≤ .
We apply this proposition to our setting. From now on, we assume Assumptions 1 and 2 hold. Then, it can be easily seen that
L̃ = LB2 + 2λ1
λ2 .
If we employ p0(θ) ∝ exp(−λ1λ2 ‖θ‖ 2) (a Gaussian measure), then we may set β ≤ exp(2B2/λ2) because
p0(θ) p∗(θ) =
( 2λ1/λ2
2π )p/2 ∫ exp(− 1λ2n∑ni=1 hi(θ′)g(ñT )i − λ1λ2 ‖θ′‖2) dθ′ exp ( − 1λ2n ∑n i=1 hi(θ)g (ñT ) i
) ≤ (
2λ1/λ2 2π
)p/2 exp (2B2/λ2) ∫ exp ( −λ1λ2 ‖θ ′‖2 ) dθ′
= exp (2B2/λ2)
( 2λ1/λ2
2π
)p/2 ∫ exp ( −λ1λ2 ‖θ ′‖2 ) dθ′
= exp (2B2/λ2) . (10)
Moreover, the proof of Lemma 16 of Ma et al. (2019) showed that we may set
R(s) = √√√√2 log( d̃ s ) d̃ cLS + √ Ep∗ [‖θ‖2].
The second term in the right hand side can be evaluated as Ep∗ [‖θ‖2] = ∫ ‖θ‖2p∗(θ)dθ = ∫ ‖θ‖2p0(θ) p∗(θ)
p0(θ) dθ ≤ exp(2B2/λ2) ∫ ‖θ‖2p0(θ)dθ
= exp(2B2/λ2) λ2 2λ1 d̃,
where the inequality in the second line can be shown in the same manner as Eq. (10). Then, we obtain
R(s) = √√√√2 log( d̃ s ) d̃ cLS + √ exp ( 2B2 λ2 ) λ2 2λ1 d̃.
Therefore, substituting these evaluations into Eq. (9) yields that
k ≥Ω
( 1
c2LS ( 2B2 λ2 + log ( 1 )) ×
max { L̃3/2 √√√√(log(2d̃ ) +
2B2 λ2
) d̃
cLS , L̃d̃, L̃4/3
( log ( 2d̃ ) +
2B2 λ2
)1/3 d̃2/3
c 1/3 LS
})
Algorithm 3 Definition of g̃(t)
1: p̃(ñT ) ← p[g(ñT )]
2: g̃(ñT+1)iñT ← argmax g̃iñT
−l∗iñT (g̃iñT )+ ∫ hiñT (θ)p̃ (ñT )(θ)dθ(g̃iñT −g (ñT ) iñT
)− 1 2nλ2 |g̃iñT −g (ñT ) iñT |2
3: for t = ñT + 1, . . . , ñ(T + 1)− 1 do 4: p̃(t) ← p[g̃(t)]
5: g̃(t+1)it ← argmax g̃it
−l∗it(g̃it) + ∫ hit(θ)p̃ (t)(θ)dθ(g̃it − g̃ (t) it
)− 1 2nλ2 |g̃it − g̃ (t) it |2
6: end for
is sufficient to ensure ‖p∗ − pk‖TV ≤ , where cLS = 2λ1λ2 exp(−4B2/λ2) and L̃ = LB2+2λ1 λ2 . We can see that the right hand side can be bounded by
O
[ 1
c 5/2 LS
( B2(1 + L) + λ1
λ2 + log(1/ )
)3/2 d̃ ] ,
which gives the assertion of Proposition 2.

heading: C PROOF OF THEOREM 1
text: In this section, we give the convergence rate analysis of our algorithm. We first show the convergence rate Algorithm 1 (the continuous limit version), then we prove Theorem 1 by bounding the discretization error. Before we begin the proofs, we need some auxiliary variables that interpolate between the continuous limit solution and the discretized one.
The key ingredient of our proof is to show the discrepancy between the solution with particle discretization and the one with the continuous limit is sufficiently small throughout each inner loop from t = ñT to t = ñ(T + 1). To rigorously prove this intuition, we introduce g̃(t) which follows the continuous limit update but is initialized from the discretized one g(ñT ) at the every starting point of the inner loop (see its definition in Algorithm 3). To fill in the gap between the continuous limit g̃(t) and the discretized version g(t), we also define ḡ(t) (t = ñT, . . . , T (ñ+ 1)) which lies in between (see Algorithm 4). Basically, ḡ(t) and g̃(t) are updated by the same rule, but are initialized in different ways. Indeed, the primal variables corresponding to these two sequences are updated as p[g̃(t+1)](θ) ∝ p[g̃(t)](θ)×exp ( − 1nλ2hit(θ)δg̃ (t) ) and p̄(t+1) ∝ p̄(t)(θ)×exp ( − 1nλ2hit(θ)δḡ (t) ) , respectively. However, different initial primal variables p̃(nT ) = p[g(nT )] and p̄(nT ) = p(nT ) are employed for the two schemes (remember that p(nT ) is an approximation of p[g(nT )] by the sampling scheme). We denote difference of the updated coordinates before and after updates as
δg(t) def = g (t+1) it − g(t)it , δg̃ (t+1) def= g̃(t+1) − g̃(t), δḡ(t) = ḡ(t+1) − ḡ(t) (t = ñT, . . . , T (ñ+ 1)− 1),
respectively.
At the same time, we define some variants of r(t)j , in order to fill the gap more finely between ḡ (t) and g(t). We define r̄(t)j as a direct discretization of p̄ (t) using the sampling point θj (see Algorithm
5). Let r̂(t)j be a numerator of r̄ (t) j (see Algorithm 6). We expect that r̂ (t) j∑M
j=1 r̂ (t) j′ plays a role of
interpolation between p̄(t) and r (t) j∑M
j=1 r (t) j′ .
In the convergence rate analysis, the smoothness of the dual objective plays an important role. The following lemma shows the smoothness of the dual of the regularization term.
Algorithm 4 Definition of ḡ(t)
1: p̄(ñT ) ← p(ñT )
2: ḡ(ñT+1)iñT ← argmax ḡiñT
−l∗iñT (ḡiñT )+ ∫ hiñT (θ)p̄ (ñT )(θ)dθ(ḡiñT −g (ñT ) iñT
)− 1 2nλ2 |ḡiñT −g (ñT ) iñT |2
3: for t = ñT + 1, . . . , ñ(T + 1)− 1 do
4: p̄(t)(θ)← p(ñT )(θ) exp
( − 1nλ2 ∑t−1 τ=ñT hiτ (θj)δḡ (τ) iτ ) ∫ p(ñT )(θ′) exp ( − 1nλ2 ∑t−1 τ=ñT hiτ (θ ′)δḡ (τ) iτ ) dθ′
5: ḡ(t+1)it ← argmax ḡit
−l∗it(ḡit) + ∫ hit(θ)p̄ (t)(θ)dθ(ḡit − ḡ (t) it
)− 1 2nλ2 |ḡit − ḡ (t) it |2
6: end for
Algorithm 5 Definition of r̄(t)j
1: r̄(ñT )j ← 1M (j = 1, . . . ,M) 2: if ñ > 1 then 3: for t = ñT, ñT + 1, . . . , ñ(T + 1)− 2 do
4: r̄(t+1)j ← 1 M exp
( − 1nλ2 ∑t τ=ñT hiτ (θj)δḡ (τ) iτ ) ∫ p(nT )(θ′) exp ( − 1nλ2 ∑t τ=ñT hiτ (θ ′)δḡ (τ) iτ ) dθ′
5: end for 6: end if
Algorithm 6 Definition of r̂(t)j
1: r̂(ñT )j ← 1M (j = 1, . . . ,M) 2: if ñ > 1 then 3: for t = ñT, ñT + 1, . . . , ñ(T + 1)− 2 do
4: r̂(t+1)j ← 1
M exp ( − 1 nλ2 t∑ τ=ñT hiτ (θj)δḡ (τ) iτ ) 5: end for 6: end if
Lemma 4. The first and second order derivatives of the second term of the dual objective D(g) can be evaluated as[
∇gλ2 log [∫ q[g](θ)dθ ]] i = − 1 n ∫ hi(θ)p[g](θ)dθ, (11)[
∇2gλ2 log [∫ q[g](θ)dθ ]] i,j =
1
n2λ2
{∫ hi(θ)hj(θ)p[g](θ)dθ − ∫ hi(θ)p[g](θ)dθ ∫ hj(θ)p[g](θ)dθ } . (12)
In particular, the dual of the regularization term can be bounded by a quadratic function as follows: (a) L1-norm bound,
λ2 log
[∫ q[g + δg](θ)dθ ] ≤ λ2 log [∫ q[g](θ)dθ ] − 1 n n∑ i=1 ∫ hi(θ)p[g](θ)dθδgi + 1 2n2λ2 ‖δg‖21,
(13)
and (b) L2-norm bound,
λ2 log
[∫ q[g + δg](θ)dθ ] ≤ λ2 log [∫ q[g](θ)dθ ] − 1 n n∑ i=1 ∫ hi(θ)p[g](θ)dθδgi + 1 2nλ2 ‖δg‖22.
(14)
The L1-norm bound is useful to derive the dual coordinate ascent method (Algorithms 1 and 2). On the other hand, the L2-norm bound is also useful in the following convergence rate analysis.
Proof. First of all, it is easy to see that λ2 log [∫ q[g](θ)dθ ] is two times differentiable by gi. Since a direct calculation yields
∂
∂gi
∫ q[g](θ)dθ = − 1
nλ2
∫ hi(θ) exp (− 1n∑ni=1 hi(θ)(gi)− λ1‖θ‖22 λ2 ) dθ, (15)
we immediately obtain Eq. (13):
∂
∂gi λ2 log
[∫ q[g](θ)dθ ] = − 1n ∫ hi(θ)q[g](θ)dθ∫ q[g](θ)dθ′ = − 1 n ∫ hi(θ)p[g](θ)dθ.
By differentiating the right hand side by gj again and using Eq. (15), we get the desired equality (12) for the Hessian.
Note that the right hand side of the equation (12) for the Hessian evaluation is given by the covariance between hi and hj with respect to the probability distribution p[g](θ). Hence, each element is bounded by 1n2λ2 due to the boundedness assumption of hi (Assumption (A2)), which yields Eq. (13). As for the L2-norm bound (14), it is easily obtained by the relation ‖δg‖1 ≤ √ n‖δg‖2 .
C.1 CONVERGENCE OF THE IDEAL UPDATE ǧ(t) (ALGORITHM 1)
First, we prove the convergence of the ideal update ǧ(t). Convergence analysis in the continuous limit is the basis for that of discretized updates; more concretely, we can use results and proofs for continuous limit in this section as a reference to derive Lemma 8, which is one of the key lemmas for Theorem 1 (convergent proof for discretized updates). Lemma 5. Suppose that `i is 1/γ-smooth for all i ∈ {1 . . . , n} (Assumption (A1)). Consider the update of ǧ(t) given in Algorithm 1 (continuous limit update), then, for any iteration t and any s ∈ [0, 1], it holds that
Eit [ D(ǧ(t+1))−D(ǧ(t)) ] ≥ s n (P (p̌(t))−D(ǧ(t)))− ( s n )2 G(t) 2λ2 , (16)
where p̌(t) = p[ǧ(t)],
G(t) def = ( 1− γ(1− s)λ2n
s
) Eit [ (u (t) it − ǧ(t)it ) 2 ]
(17)
and u(t)i def = `′i (∫ p̌(t)(θ)hi(θ)dθ ) . Here, the expectation was taken with respect to the choice of the coordinate it. Accordingly, we also have the same bound for the ideal update in Algorithm 2 (g(ñT ) and g̃(t) (t = ñT + 1, . . . , ñ(T + 1))), considering the definition of g̃(t).
Proof. Since we will focus on the update of one step, it is omitted as i in the following. At each update, the improvement in the dual objective can be written as
n[D(ǧ(t+1))−D(ǧ(t))]
=− `∗i (ǧ (t+1) i ) + ` ∗ i (ǧ (t) i )− λ2 log
[∫ q[ǧ(t+1)](θ)dθ ] + λ2 log [∫ q[ǧ(t)](θ)dθ ] ≥− `∗i (ǧ (t+1) i ) + ∫ hi(θ)p̌ (t)(θ)dθ (ǧ (t+1) i − g (t) i )− 1
2nλ2 |gi − ǧ(t)i | 2 + `∗i (ǧ (t) i )
≥− `∗i (ǧ (t) i + s(u (t) i − ǧ (t) i )) + s
∫ h(θ)p̌(t)(θ)dθ (u
(t) i − ǧ (t) i )−
s2
2nλ2 |u(t)i − ǧ (t) i | 2 + `∗i (ǧ (t) i ),
(18)
where we used Eq. (13) in Lemma 4 for the first inequality, and the second inequality follows from the definition of ǧ(t+1).
Since `∗i is γ-strongly convex (γ = 0 is also allowed), it holds that
`∗i (ǧ (t) i + s(u (t) i − ǧ (t) i )) ≤ s` ∗ i (u (t) i ) + (1− s)` ∗ i (ǧ (t) i )−
γ 2 s(1− s)|u(t)i − ǧ (t) i | 2. (19)
From the property of the Legendre transform, u(t)i = ` ′ i (∫ p̌(t)(θ)hi(θ)dθ ) implies that
`i (∫ p̌(t)(θ)hi(θ)dθ ) + `∗i (u (t) i ) = ∫ hi(θ)p̌ (t)(θ)dθu (t) i . (20)
Similarly, we also have the following relation from the duality of the regularization term:
λ1 ∫ p̌(t)(θ)||θ||22dθ + λ2 ∫ p̌(t)(θ) log p(t)(θ)dθ + λ2 log [∫ q[ǧ(t)](θ)dθ ] = − 1
n n∑ i=1 ∫ hi(θ)p̌ (t)(θ)dθǧ (t) i . (21)
Combining Eqs. (18), (19) and (20), we have that
n[D(ǧ(t+1))−D(ǧ(t))] ≥− s ( `∗i (u (t) i )− ∫ hi(θ)p̌ (t)(θ)dθ u (t) i ) + s
2
( γ(1− s)− s
nλ2
) |u(t)i − ǧ (t) i | 2
+ s ( `∗i (ǧ (t) i )− ǧ (t) i ∫ hi(θ)p̌ (t)(θ)dθ ) =s [ `i (∫ p̌(t)(θ)hi(θ)dθ ) + `∗i (ǧ (t) i )− ∫ hi(θ)p̌ (t)(θ)dθǧ (t) i
] + s
2
( γ(1− s)− s
nλ2
) |u(t)i − ǧ (t) i | 2. (22)
By taking expectation of both sides of (22) with respect to i, it follows that E [ n[D(ǧ(t+1))−D(ǧ(t))] ] ≥s [ 1
n n∑ i=1 `i (∫ p̌(t)(θ)hi(θ)dθ ) + 1 n n∑ i=1 `∗i (ǧ (t) i )− 1 n n∑ i=1 ∫ hi(θ)p̌ (t)(θ)dθǧ (t) i ]
+ s
2n
( γ(1− s)− s
nλ2 ) n∑ i=1 |u(t)i − ǧ (t) i | 2
=s
[ 1
n n∑ i=1 `i (∫ p̌(t)(θ)hi(θ)dθ ) + 1 n n∑ i=1 `∗i (ǧ (t) i )
+ λ1 ∫ p̌(t)(θ)||θ||22dθ + λ2 ∫ p̌(t)(θ) log p̌(t)(θ)dθ +λ2 log (∫ q[ǧ(t)](θ)dθ )] + s
2n
( γ(1− s)− s
nλ2 ) n∑ i=1 |u(t)i − ǧ (t) i | 2
=s(P (p̌(t))−D(ǧ(t)))− s 2
n
G(t) 2λ2 , (23)
where we used Eq. (21) for the first equality. Dividing the both sides of Eq. (23) by n, we obtain the desired inequality (16).
Lemma 6. ǧ(t)i in Algorithm 1 is uniformly bounded by B2, that is,
|ǧ(t)i | ≤ B2 (∀t = 0, 1, . . . , nTend, ∀i ∈ {1, . . . , n}).
In the same way, g(t)i , g̃ (t) i and ḡ (t) i in Algorithm 2 satisfy
|g(t)i |, |g̃ (t) i |, |ḡ (t) i | ≤ B2 (∀t = 0, 1, . . . , nTend, ∀i ∈ {1, . . . , n}).
Proof. Let p be an arbitrary probability distribution. Let us consider the update
z = argmax z′∈R
{ −`∗i (z′) + (z′ − y) ∫ hi(θ)p(θ)dθ } , (24)
for an arbitrary i ∈ {1, . . . , n}. Since hi is bounded by 1 (Assumption (A2)), the sub-differential of `∗i at z satisfies that ∂`∗i (z) ∩ [−1, 1] 6= ∅. By the assumption (A3), i.e.,
B2 = sup i,x
{ |`i′(x)| ∣∣ |x| ≤ 1} = sup i,y {|y| | ∂`∗i (y) ∩ [−1, 1] 6= ∅} <∞,
we have that z in Eq. (24) is bounded by B2.
Next, we consider its proximal-version update:
z = argmax z′∈R
{ −`∗i (z′) + ∫ hi(θ)p(θ)dθ(z
′ − y)− 1 2nλ2
|z′ − y|2 } . (25)
In this proximal-version update, we added the term − 12nλ2 (z − y) 2 whose maximizer with respect to z is given by z = y. Since the maximum point of the sum of two concave functions is located between the maximizers of them, z in Eq. (25) satisfies that
|z| ≤ max {B2, |y|} .
From this and the initialization |ǧ(0)i | ≤ B2, we have that |ǧ (t) i | in Algorithm 1 is uniformly bounded by B2 by induction.
In the same way, g(t)i , g̃ (t) i and ḡ (t) i in Algorithm 2 minimize the sum of the average of hi(θj) and quadratic term as well. Note that the absolute value of the the average of hi(θj) is bounded with 1, regardless of methods of averaging, that is, weighted particle for g(t)i and probability distribution for g̃
(t) i and ḡ (t) i . Thus, g (t) i , g̃ (t) i and ḡ (t) i are uniformly bounded with B2, as well as ǧ (t)i in Algorithm 1.
Remark 4. In the Section C.2, we sometimes face the situation where we should upper-bound increment of the coordinates δg(t), δg̃(t) and δḡ(t).
Lemma 6 directly implies that δg(t), δg̃(t) and δḡ(t) are all bounded by 2B2 for all t = 0, 1, . . . , ñTend. Thus, for future convenience, we define
C1 def = exp ( 2κB2 λ2 ) with which exp ( 2κ|δg(t)|
λ2
) , exp ( 2κ|δg̃(t)|
λ2
) , exp ( 2κ|δḡ(t)|
λ2
) ≤ C1 hold, where κ is defined as
κ = ñn .
Note that, however, we could have much tighter evaluations of δg(t), δg̃(t) and δḡ(t), which are dependent of residual D(g∗)−D(g(t)) or duality gap P [p(g(t))]−D(g(t)), since the dual problem is strongly convex due to the smoothness of primal problem.
Lemma 7. Under Assumptions (A2) and (A3), it holds that
P (p[ǧ(0)])−D(ǧ(0)) ≤ B1 +B2,
for Algorithm 1 and
P (p[g(0)])−D(g(0)) ≤ B1 +B2,
in Algorithm 2. Proof. By Assumptions (A2) and (A3), we can see that ∣∣`i (∫ p(0)(θ)hi(θ)dθ)− infx `i(x)∣∣ is bounded by B1. Since infx `i(x) > − inf (Assumption (A3)), we also have `∗i (0) = supx{0 · x− `i(x)} = − infx `i(x) <∞. In Algorithm 1, every initial coordinate ǧ(0)i is bounded by B2. Also, g
(0) i in Algorithm 2 is similarly bounded by B2. Therefore, in the following, g will be used without distinguishing between g and ǧ. Moreover, in a similar manner to Eq. (21) in Lemma 5, we have that
λ1 ∫ p[g(t)](θ)||θ||22dθ + λ2 ∫ p[g(t)](θ) log p[g(t)](θ)dθ + λ2 log [∫ q[g(t)](θ)dθ ] = − 1
n n∑ i=1 ∫ hi(θ)p[g (t)](θ)dθg (t) i ,
from the duality.
Combining these observations, we have that
P (p[g(0)])−D(g(0))
= 1
n n∑ i=1 `i (∫ p(0)(θ)hi(θ)dθ ) + 1 n n∑ i=1 `∗i (g (0) i ) + λ1 ∫ p(0)(θ)||θ||22dθ
+ λ2 ∫ p(0)(θ) log p(0)(θ)dθ + λ2 log [∫ q[g(0)](θ)dθ ] = 1
n n∑ i=1 `i (∫ p(0)(θ)hi(θ)dθ ) + 1 n n∑ i=1 `∗i (g (0) i )− 1 n n∑ i=1 ∫ hi(θ)p (0)(θ)dθ · g(0)i
≤ 1 n n∑ i=1 [ `i (∫ p(0)(θ)hi(θ)dθ ) − inf x `i(x) ] + 1 n n∑ i=1 |g(0)i |
≤ 1 n n∑ i=1 B1 +B2 = B1 +B2,
where we used Lemma 6 in the last inequality. Theorem 2. Let P (> 0) be a desired accuracy. To obtain an expected duality gap E [ P (p̌(tend))−D(ǧ(tend)) ] ≤ P in Algorithm 1, it suffices to have the total number Tend of iterations as
tend ≥ ( n+ 1
γλ2
) log (( n+ 1
γλ2 ) B1 +B2 P ) .
Furthermore, if we randomly choose (p̌(t), ǧ(t)) from t = tend − n+ 1, . . . , tend in Algorithm 1, then tend ≥ ( n+ 1
γλ2
) log (( 1 + 1
nγλ2 ) B1 +B2 P ) + n
is sufficient to achieve P duality gap in expectation.
Proof. In Lemma 4, choosing s = nγλ21+nγλ2 ∈ [0, 1], then we have that
1− (1− s)nγλ2 s ≤ 1− 1 = 0.
Hence, we have G(t) ≤ 0 for all t. Thus, Lemma 5 implies that E [ D(ǧ(t+1))−D(ǧ(t)) ] ≥ s n E [ P (p̌(t))−D(ǧ(t)) ] . (26)
Therefore, we have that s n E [ D(g∗)−D(ǧ(t)) ] ≤ s n E [ P (p̌(t))−D(ǧ(t)) ] ≤ s n E [ D(ǧ(t+1))−D(ǧ(t)) ] = E [ D(g∗)−D(ǧ(t)) ] − E [ D(g∗)−D(ǧ(t+1)) ] .
This implies that E [ D(g∗)−D(ǧ(t+1)) ] ≤ (
1− s n
) E [ D(g∗)−D(ǧ(t)) ] . (27)
By applying this inequality repeatedly, we obtain that E [ D(g∗)−D(ǧ(t)) ] ≤ (
1− s n
)t E [ D(g∗)−D(ǧ(0)) ] ≤ (B1 +B2) exp ( −st n ) , (28)
where the last inequality is obtained by E [ D(g∗)−D(ǧ(0)) ] ≤ E [ P (p̌(0))−D(ǧ(0)) ] and Lemma 7.
Moreover, making use of Eq. (26), the duality gap can be bounded as follows: E [ P (p̌(t))−D(ǧ(t)) ] ≤ n s E [ D(ǧ(t+1))−D(ǧ(t)) ] ≤ n s E [ D(g∗)−D(ǧ(t)) ] . (29)
Hence, combining Eqs. (28) and (29) yields E [ P (p̌(t))−D(ǧ(t)) ] ≤ n (B1 +B2)
s exp
( − γλ2t
1 + nγλ2
) .
Therefore, in order to have an P -approximated solution in expectation, it suffices to let tend ≥ ( n+ 1
γλ2
) log (( n+ 1
γλ2 ) B1 +B2 P ) .
Finally, we show the second part. By summing up Eq. (26) over t = tend − n+ 1, . . . , tend, we have that
E
[ 1
n tend∑ t=tend−n+1 ( P (p̌(t))−D(ǧ(t)) )] ≤1 s E [ D(ǧ(tend+1))−D(ǧ(tend−n)) ] ≤1 s E [ D(g∗)−D(ǧ(tend−n))
] ≤B1 +B2
1 + 1nγλ2 exp
( −γλ2(tend − n)
1 + nγλ2
) .
Note that although we do not actually calculate ǧ(tend+1), we can virtually define ǧ(tend+1) as if the algorithm continued.
Therefore, in order to have an P -approximated solution, it suffices to let tend ≥ ( n+ 1
γλ2
) log (( 1 + 1
nγλ2 ) B1 +B2 P ) + n,
when we randomly choose (p̌(t), ǧ(t)) from t = tend − n+ 1, . . . , tend.
C.2 CONVERGENCE OF THE PARTICLE-VERSION UPDATE g(t) (ALGORITHM 2)
Next, we consider the convergence in Algorithm 2.
In order to prove Theorem 1 and give a convergence proof for Algorithm 2, we adopt the following strategy. First we prove Lemma 8, by rewriting Theorem 2 to allow the discretization error in each
series of inner loops. Then, the assumptions of Theorem 2 is verified, by Theorem 3. There are two major streams leading up to Theorem 3. The first is about the error originated from the sampling accuracy (Langevin algorithms), and the second is to evaluate the error accumulated during the re-weighting stage. The former is considered in Section C.2.1 and the latter is in Section C.2.2, both of which are combined in Section C.2.3 and finally Theorem 3 is proved.
Apart from the sources of the error, we bound the discrepancy of the duality gap between continuous limit and discretized updates (P (g(t)) − D(g(t))) − (P (g̃(t)) − D(g̃(t))) with the following two parts,
(t) A def = E
[ D(g̃(t))−D(g(t)) ] ,
(t) B def = E
[ P (p[g(t)])− P (p[g̃(t)]) ] ,
each of which corresponds to the discrepancy in the primal problem and dual problem, respectively. As mentioned above, these two terms are attributed to the two sources, that is, sampling and reweighting. These variables (t)A , (t) B represent how much the two sequences g
(t) and g̃(t) diverge during the updates of the inner-loop. Intuitively, in each inner-loop, the absolute values of (t)A and (t) B are increasing sequences between t = ñT + 1 and ñ(T + 1). Then, they are “reset” at t = ñ(T + 1) according to the initialization rule of g̃(ñ(T+1)+1) and again will increase from t = ñ(T + 1) + 1 to ñ(T + 2).
The following lemma states that under the assumption about (t)A , (t) B , Algorithm 2 achieves linear convergence with regard to the number of coordinate update. From now on, we fix s = nγλ21+nγλ2 . Also, we define s̃ = κ · s2 = ñγλ2 1+nγλ2 , where κ = ñn .
Lemma 8. Suppose that Assumptions 1 and 2 hold. If (ñT )A is smaller than C0e −s̃T for ev-
ery T = 1, . . . , Tend with a constant C0 you can arbitrarily take, an expected duality gap E [ P (p[g(nTend)])−D(g(nTend)) ] ≤ P can be achieved by a total number Tend of outer iterations given by
Tend ≥ 2 n
ñ
( 1 + 1
nγλ2
) log (( n+ 1
γλ2
) (B1 +B2 + C0 1−exp(−s̃) )
P
) . (30)
Furthermore, if (t)A and (t) B are uniformly bounded by P 4 for all t = ñ(Tend − d 1 κe) + 1, . . . , nTend in addition to the assumptions above, then it also suffices to set
Tend ≥ 2 n
ñ
( 1 + 1
nγλ2
) log (( 2 + 2
nγλ2
) (B1 +B2 + C0 1−exp(−s̃) )
P
) + ⌈ 1 κ ⌉ , (31)
to achieve P -duality gap when we randomly choose output (p(t), g(t)) from t = ñ(Tend − d 1κe) + 1, . . . , ñTend. Remark 5. Comparing the speed of convergence with regard to coordinate update t, it is exp(−st/ñ) for Theorem 2 and exp(−s̃T ) ∝ exp(− s2 t/ñ) for Lemma 8, which can be interpreted as follows. In the descretized version, Lemma 8 allows discretization error to appear, at the cost of a slight decrease in convergence speed. Indeed, you can choose arbitrary s̃ from the open interval (0, κs). Our choice of s̃ = κs2 is just for simplicity. If you choose another s̃ much closer to κs, Tend can be smaller. Remark 6. If you change the length of inner loop in the middle, the convergence rate after the change is obtained by simply replacing B1 +B2 with the duality gap at the time of the change.
Proof. Since the update of g̃(t) is the same as that of ǧ(t) in Algorithm 1 during the inner loop, it holds that
E [ D(g∗)−D(g̃(ñT+1)) ] ≤ (
1− s n
) E [ D(g∗)−D(g(ñT )) ] ,
for T = 0, 1, . . . , Tend − 1. We also have E [ D(g∗)−D(g̃(t+1)) ] ≤ (
1− s n
) E [ D(g∗)−D(g̃(t)) ] ,
for t = ñT + 1, . . . , ñ(T + 1)− 1, in the same way as the derivation of Eq. (27). By applying these inequalities above iteratively, we have that
E [ D(g∗)−D(g̃(ñ(T+1))) ] ≤ (
1− s n
)ñ E [ D(g∗)−D(g(ñT )) ] ≤ e−2s̃E [ D(g∗)−D(g(ñT )) ] . (32)
Bounding D(g̃(ñ(T+1))) by D(g(ñ(T+1)))− (ñ(T+1))A , we have that E [ D(g∗)−D(g(ñ(T+1))) ] ≤ e−2s̃E [ D(g∗)−D(g(ñT )) ] +
(ñ(T+1)) A .
Therefore, by the Gronwall’s lemma (“vanilla version” of Lemma 18), we obtain that
E [ D(g∗)−D(g(ñT )) ] ≤ e−2s̃TE [ D(g∗)−D(g(0)) ] + e−2s̃T T∑ T ′=1 e2s̃T ′ (nT ′) A .
Applying the assumption (nT )A ≤ C0e−s̃T , this implies that
E [ D(g∗)−D(g(ñT )) ] ≤ e−sTE [ D(g∗)−D(g(0)) ] + C0 e−s̃T − e−2s̃T
1− e−s̃ ≤e−s̃T ( E [ P (g(0))−D(g(0)) ] +
C0 1− e−s̃ ) ≤e−s̃T ( B1 +B2 +
C0 1− e−s̃
) , (33)
where we used Lemma 7 for the last inequality.
Moreover, the same reasoning as Eq. (29) yields that E [ P (p[g(ñT )])−D(g(ñT )) ] ≤ n s E [ D(g∗)−D(g(ñT )) ] . (34)
Therefore, combining Eqs. (33) and (34), we can bound the duality gap as E [ P (p[g(ñT )])−D(g(ñT )) ] ≤ e−s̃T n
s
( B1 +B2 +
C0 1− e−s̃
) .
This inequality provides the sufficient number of iterations to attain the desired duality gap. In order to have an P -approximated solution, it is sufficient to let
Tend ≥ 2 n
ñ
( 1 + 1
nγλ2
) log (( n+ 1
γλ2
) B1 +B2 + C0 1−exp(−s̃)
P
) .
Next, we show the convergence of the averaging version. Similarly to Theorem 2, we have that for k = 1, . . . , d 1κe,
E  1 ñ ñ(Tend−k+1)∑ t=ñ(Tend−k)+1 ( P (p[g̃(t)])−D(g̃(t)) ) ≤1 s E [ D(g∗)−D(g̃(ñ(Tend−k)+1)) ] ≤1 s e−s̃(Tend−d1/κe) ( B1 +B2 + C0 1− e−s̃ ) .
Note that, if Eq. (31) holds, we have that 1se −s̃(Tend−d1/κe) ( B1 +B2 + C0 1−e−s̃ ) ≤ P2 . Hence, if
(t) A and (t) B are uniformly bounded by P 4 for all t = ñ (Tend − d1/κe) + 1, . . . , ñTend, it holds that
E  1 d1/κeñ nTend∑ t=n(Tend−d1/κe)+1 ( P (p[g(t)])−D(g(t)) )
≤ d1/κe∑ k=1 1 d1/κe E  1 ñ ñ(Tend−k+1)∑ t=ñ(Tend−k)+1 ( P (p[g̃(t)])−D(g̃(t)) ) + 1
d1/κeñ ñTend∑ t=ñ(Tend−d1/κe)+1 (t) A +
1
d1/κeñ ñTend∑ t=ñ(Tend−d1/κe)+1 (t) B
≤ P 2 + P 4 + P 4 = P .
This concludes the proof.
In the following, we derive the sufficient conditions for the assumptions in Lemma 8 by considering the error induced by the sampling stage (Section C.2.1) and the error induced at the re-weighting stage (Section C.2.2), respectively. For that purpose, we first prepare lemmas that will be used in the both proofs.
The following Lemma 9 evaluates how much the difference of two solutions in the previous step incurs difference of them in the next update in Algorithm 2. This will be used to bound the accumulated error through each inner loop. Lemma 9. For arbitrary i ∈ [n], consider the following two updates:
z1 = argmax z′1∈R
{ −`∗i (z′1) +H1(z′1 − y1)− 1
2nλ2 |z′1 − y1|2
} ,
z2 = argmax z′2∈R
{ −`∗i (z′2) +H2(z′2 − y2)− 1
2nλ2 |z′2 − y2|2
} .
Then, we have that
|z2 − z1| ≤ 1
γ + 1nλ2
( |H2 −H1|+ 1
nλ2 |y2 − y1|
) .
Proof. From the definition of z1, it holds that
−`∗i ′(z1) +H1 −
1
nλ2 (z1 − y1) = 0, (35)
and similarly,
−`∗i ′(z2) +H2 −
1
nλ2 (z2 − y2) = 0,
where `∗i ′(z) is an element in the sub-differential ∂`∗i (z). From (35), we have that
−`∗i ′(z1) +H2 −
1
nλ2 (z1 − y2) = H2 −H1 +
1
nλ2 (y2 − y1). Since `∗i is γ-strongly convex,−`∗i (z)+H2(z−y2)− 12nλ2 |z−y2| 2 is ( γ + 1nλ2 ) -strongly concave. Thus, we obtain that |z2 − z1|
≤ 1 γ + 1nλ2
∣∣∣∣∣∇ ( −`∗i (z) +H2(z − y2)− 1 2nλ2 |z − y2|2 )∣∣∣∣ z=z2
−∇ ( −`∗i (z) +H2(z − y2)− 1
2nλ2 |z − y2|2 )∣∣∣∣ z=z1 ∣∣∣∣∣ ≤ 1 γ + 1nλ2 ∣∣∣∣H2 −H1 + 1nλ2 (y2 − y1) ∣∣∣∣
≤ 1 γ + 1nλ2
( |H2 −H1|+ 1
nλ2 |y2 − y1|
) ,
where ∇f(z)|z=z′ is an element of the sub-differential of f at z = z′. This concludes the proof.
Now we are ready to get into the main proofs which evaluate the discretization error.
C.2.1 CONVERGENCE OF AUXILIARY VARIABLES g̃(t) AND ḡ(t)
We now fix some T , and consider the sequances {g(t)}ñ(T+1)t=ñT , {ḡ (t)}ñ(T+1)t=ñT+1, {g̃ (t)}ñ(T+1)t=ñT+1. Lemma 10. Suppose that ∫ ∣∣p(ñT )(θ)− p[g(ñT )](θ)∣∣dθ ≤ (ñT )C holds in Algorithm 2. Under these assumptions, it holds that
|ḡ(t)i − g̃ (t) i | ≤ C2 (ñT ) C , (36)
for t = ñT + 1, . . . , ñ(T + 1), where C2 is defined as follows:
C2 def = 2C21 γ + 1nλ2 exp
( ñ(2C21 + 1)
nγλ2 + 1
) .
Proof. By Lemma 9, we have that
|ḡ(t+1)it − g̃ (t+1) it |
≤ 1 γ + 1nλ2 (∣∣∣∣∫ hit(θ)p̄(t)(θ)dθ − ∫ hit(θ)p[g̃(t)](θ)dθ∣∣∣∣+ 1nλ2 |ḡ(t)it − g̃(t)it | )
≤ 1 γ + 1nλ2 (∫ ∣∣∣p̄(t)(θ)dθ − p[g̃(t)](θ)∣∣∣dθ + 1 nλ2 |ḡ(t)it − g̃ (t) it | ) , (37)
so we evaluate each terms of Eq. (37).
First, we bound the first term of the right hand side of Eq. (37). Using Lemma 6, we have the following bounds:
1
C1 ≤exp ( − 1 nλ2 t−1∑ τ=ñT hiτ (θ)δḡ (τ) iτ ) ≤C1, 1 C1 ≤exp ( − 1 nλ2 t−1∑ τ=ñT hiτ (θ)δg̃ (τ) iτ ) ≤C1. (38)
Using (38), we have that∫ ∣∣∣∣∣p(ñT )(θ) exp ( − 1 nλ2 t−1∑ τ=ñT hiτ (θ)δḡ (τ) iτ ) −p[g(ñT )](θ) exp ( − 1 nλ2 t−1∑ τ=ñT hiτ (θ)δg̃ (τ) iτ )∣∣∣∣∣dθ ≤C1
∫ ∣∣∣p(ñT )(θ)− p[g(ñT )](θ)∣∣∣dθ + ∫ p[g(ñT )](θ) ∣∣∣∣∣exp ( − 1 nλ2 t−1∑ τ=ñT hiτ (θ)δḡ (τ) iτ ) − exp ( − 1 nλ2 t−1∑ τ=ñT hiτ (θ)δg̃ (τ) iτ
)∣∣∣∣∣dθ ≤C1 ( (ñT ) C + 1
nλ2 t−1∑ τ=ñT |ḡ(τ+1)iτ − g̃ (τ+1) iτ |
) . (39)
Thus, combining Eqs. (38) and (39) yields we have the following bound for the first term of Eq. (37):∫ ∣∣∣p̄(t)(θ)− p[g̃(t)](θ)∣∣∣dθ = ∫ ∣∣∣∣ p(ñT )(θ) exp ( − 1nλ2 ∑t−1 τ=ñT hiτ (θ)δḡ (τ) iτ ) ∫ p(ñT )(θ′) exp ( − 1nλ2 ∑t−1 τ=ñT hiτ (θ ′)δḡ (τ) iτ ) dθ′ − p[g(ñT )](θ) exp ( − 1nλ2 ∑t−1 τ=ñT hiτ (θ)δg̃ (τ) iτ ) ∫ p[g(ñT )](θ′) exp ( − 1nλ2 ∑t−1 τ=ñT hiτ (θ ′)δg̃ (τ) iτ ) dθ′
∣∣∣∣ dθ ≤ ∫ ∣∣∣∣ p(ñT )(θ) exp ( − 1nλ2 ∑t−1 τ=ñT hiτ (θ)δḡ (τ) iτ ) ∫ p(ñT )(θ′) exp ( − 1nλ2 ∑t−1 τ=ñT hiτ (θ ′)δḡ (τ) iτ ) dθ′ − p(ñT )(θ) exp ( − 1nλ2 ∑t−1 τ=ñT hiτ (θ)δḡ (τ) iτ ) ∫ p[g(ñT )](θ′) exp ( − 1nλ2 ∑t−1 τ=ñT hiτ (θ ′)δg̃ (τ) iτ ) dθ′
∣∣∣∣ dθ + ∫ ∣∣∣∣p(ñT )(θ) exp ( − 1nλ2 ∑t−1 τ=ñT hiτ (θ)δḡ (τ) iτ ) −p[g(ñT )](θ) exp ( − 1nλ2 ∑t−1 τ=ñT hiτ (θ)δg̃ (τ) iτ ) ∫ p[g(ñT )](θ′) exp ( − 1nλ2 ∑t−1 τ=ñT hiτ (θ ′)δg̃ (τ) iτ ) dθ′
∣∣∣∣ dθ ≤ ∣∣∣∫ p[g(ñT )](θ′) exp(− 1nλ2 ∑t−1τ=ñT hiτ (θ′)δg̃(τ)iτ )dθ′−∫ p(ñT )(θ′) exp(− 1nλ2 ∑t−1τ=ñT hiτ (θ′)δḡ(τ)iτ )dθ′∣∣∣∫
p[g(ñT )](θ′) exp ( − 1nλ2 ∑t−1 τ=ñT hiτ (θ ′)δg̃ (τ) iτ ) dθ′
+ ∫ ∣∣∣p(ñT )(θ) exp(− 1nλ2 ∑t−1τ=ñT hiτ (θ)δḡ(τ)iτ )−p[g(ñT )](θ) exp(− 1nλ2 ∑t−1τ=ñT hiτ (θ)δg̃(τ)iτ )∣∣∣dθ∫ p[g(ñT )](θ′) exp ( − 1nλ2 ∑t−1 τ=ñT hiτ (θ ′)δg̃ (τ) iτ ) dθ′
≤2C1 ∣∣∣∣∣p(ñT )(θ) exp ( − 1 nλ2 t−1∑ τ=ñT hiτ (θ)δḡ (τ) iτ ) − p[g(ñT )](θ) exp ( − 1 nλ2 t−1∑ τ=ñT hiτ (θ)δg̃ (τ) iτ )∣∣∣∣∣ ≤2C21 ( (ñT ) C + 1
nλ2 t−1∑ τ=ñT |ḡ(τ+1)iτ − g̃ (τ+1) iτ |
) . (40)
where we used the boundedness of each term following Eq. (38) for the second from last inequality and (39) for the last.
Finally, we bound |ḡ(t+1)it − g̃ (t+1) it | by applying (40) to (37):
|ḡ(t+1)it − g̃ (t+1) it |
≤ 1 γ + 1nλ2
( 2C21 ( (ñT ) C + 1
nλ2 t−1∑ τ=ñT |ḡ(τ+1)iτ − g̃ (τ+1) iτ |
) + 1
nλ2 |ḡ(τ+1)iτ − g̃ (τ+1) iτ |
)
≤ 1 γ + 1nλ2
( 2C21 (ñT ) C + 2C21 + 1
nλ2 t−1∑ τ=ñT |ḡ(τ+1)iτ − g̃ (τ+1) iτ |
) . (41)
Moreover, we have that
|ḡ(ñT+1)iñT − g̃ (ñT+1) iñT | ≤ 1 γ + 1nλ2 ∣∣∣∣∫ hiñT (θ)p(ñT )(θ)dθ − ∫ hiñT (θ)p[g(ñT )](θ)dθ∣∣∣∣ ≤ 1 γ + 1nλ2 ∫ |p(ñT )(θ)− p[g(ñT )](θ)|dθ ≤ (ñT ) C γ + 1nλ2 . (42)
Combining Eqs. (41) and (42), we can apply Gronwall’s lemma (“summation version” of Lemma 18): we have that, for t = ñT, ñT + 1, . . . , ñ(T + 1)− 1,
|ḡ(t+1)it − g̃ (t+1) it |
≤  2C21 γ + 1nλ2 + 2C21 + 1(
γ + 1nλ2
)2 nλ2 (1 + 2C21 + 1 nγλ2 + 1 )ñ−2 (ñT ) C
≤ 2C 2 1
γ + 1nλ2
( 1 + 2C21 + 1
nγλ2 + 1
)ñ−1 (ñT ) C
≤ 2C 2 1
γ + 1nλ2 exp
( ñ(2C21 + 1)
nγλ2 + 1
) (ñT ) C ,
where we used C1 ≥ 1 for the second inequality. Remember that we defined C2 as C2 = 2C21 γ+ 1nλ2 exp ( ñ(2C21+1) nγλ2+1 ) . By the definition of g̃(t) and ḡ(t), ḡ(t)i = ḡ (t) i holds for t = ñT +
1, . . . , ñ(T + 1) unless the coordinate i is chosen during the inner loop. Otherwise, we have∣∣∣g̃(t)i − ḡ(t)i ∣∣∣ ≤ C2 (nT )C for all t ≥ τ + 1 where τ is the first time at which the coordinate i is chosen as we have seen above. Thus, Eq. (36) holds for all t = ñT + 1, . . . , ñ(T + 1) and i = 1, . . . , n.
C.2.2 PARTICLE APPROXIMATION, AND ERROR ANALYSIS OF g(t)
From now on, we bound the error induced by particle sampling in the re-weighting stage. The following three Lemmas relate, in three steps, the update of the continuous limit to the approximation by particle weighting.
Lemma 11. With probability 1− δ, it holds that∣∣∣∣∣∣ M∑ j=1 hit(θj)r̄ (t) j − ∫ hit(θ)p̄ (t)(θ)dθ ∣∣∣∣∣∣ ≤ 7C21 √ 1 M log ( 2n δ ) ,
uniformly over all t = ñT, . . . , ñ(T + 1)− 1 and the choice of the coordinates (iτ )ñ(T+1)−1τ=ñT . .
Proof. For t = ñT + 1, . . . , ñ(T + 1)− 1, let
H := { fI,a(·) = hit(·) exp ( − 1 nλ2 t−1∑ τ=ñT hiτ (·)aτ ) C−1I,a
| it ∈ [n], I = (iñT , . . . , it−1) ∈ [n]t−ñT , |aτ | ≤ 2B2 (τ ∈ {ñT, . . . , t− 1}) } ,
where CI,a := ∫ p(ñT )(θ′) exp ( − 1nλ2 ∑t−1 τ=ñT hiτ (θ ′)aτ ) dθ′ for I = (iñT , . . . , it−1) and a = (añT , . . . , at−1). We also define
Hi := { fI,a(·) = hi(·) exp ( − 1 nλ2 t−1∑ τ=ñT hiτ (·)aτ ) C−1I,a
| I = (iñT , . . . , it−1) ∈ [n]t−ñT , |aτ | ≤ 2B2 (τ ∈ {ñT, . . . , t− 1}) } .
Then, we can see thatH = ∪i∈[n]Hi. Then, ∑M j=1 hit(θj)r̄ (t) j can be rewritten as
M∑ j=1 hit(θj)r̄ (t) j = 1 M M∑ j=1 fI,a(θj).
In a similar way, we also have that∫ hit(θ)p̄ (t)(θ)dθ = Eθ[fI,a(θ)].
To check this, first note that (Mhit(θj)r̄ (t) j ) M j=1 are i.i.d. random variables for each t = ñT + 1, . . . , ñ(T + 1)− 1. Indeed, r̄(t)j is determined by ḡ(t) which is updated by using a continuous limit and thus is independent of the choice of particles (θj)Mj=1. Their expectations can be evaluated as
Eθj [ Mhit(θj)r̄ (t) j ] = ∫ hit(θ)
p(ñT )(θ) exp ( − 1nλ2 ∑t−1 τ=ñT hiτ (θj)δḡ (τ) iτ ) ∫ p(ñT )(θ′) exp ( − 1nλ2 ∑t−1 τ=ñT hiτ (θ ′)δḡ (τ) iτ ) dθ′ dθ
= ∫ hit(θ)p̄ (t)(θ)dθ.
Moreover, they are uniformly bounded as
|Mhit(θj)r̄ (t) j | ≤ |hit(θj)|
exp ( − 1nλ2 ∑t−1 τ=ñT hiτ (θj)δḡ (τ) iτ ) ∫ p(ñT )(θ′) exp ( − 1nλ2 ∑t−1 τ=ñT hiτ (θ ′)δḡ (τ) iτ ) dθ′ ≤ C21 . (43)
The same results are also true for t = ñT .
Note that it suffices to bound Z = supfI,a∈H | 1 M ∑M j=1 fI,a(θj)− Eθ[fI,a(θ)]| to obtain the assertion. For that purpose, we define the Rademacher complexity of a function classH′ as
R(H′) := E  sup f∈H′ 1 M M∑ j=1 εjf(θj)  , where (εj)Mj=1 is an i.i.d. sequence of the Rademacher random variables (P (εj = 1) = P (εj = −1) = 1/2) and the expectation is taken with respect to both (εj)Mj=1 and (θj)Mj=1. Then, the
Rademacher concentration inequality (Mohri et al., 2012, Theorem 3.1) tells that the uniform bound on a function classH′ can be obtained as
P  sup f ′∈H′ ∣∣∣∣∣∣ 1M M∑ j=1 (f(θj)− Eθ[f(θ)]) ∣∣∣∣∣∣ ≥ 2R(H′) + supf∈H′ ‖f‖∞ √ 1 2M log ( 2 δ ) ≤ δ, (44) for 0 < δ < 1. Here, let Zi = supf∈Hi | 1 M ∑M j=1 f(θj) − Eθ[f(θ)]| and then we also see that Z = maxi∈[n] Zi. Hence, we have that
P (Z ≥ s) ≤ n∑
it=1
P (Zit ≥ s),
for s > 0. If we substitute s ← maxi∈[n] 2R(Hi) + supf∈Hi ‖f‖∞ √ 1 2M log ( 2n δ ) , then by Eq. (44) we can see that the right hand side can be bounded by δ. Note that we have already seen that ‖f‖∞ ≤ C21 for any f ∈ H (see Eq. (43)). Therefore, we just need to bound R(Hit). By the contraction property of the Rademacher complexity ((Boucheron et al., 2013, Theorem 11.6) or (Ledoux & Talagrand, 1991, Theorem 4.12) and its proof) and remembering the fact that the Rademacher complexity of a function class is same as that of its convex hull (Mohri et al., 2012, Theorem 6.2), we obtain
R(Hit)
≤ E sup I,a 1 M M∑ j=1 εj exp ( − 1 nλ2 t−1∑ τ=ñT hiτ (θj)aτ − logCI,a ) (∵ |hit(·)| ≤ 1) ≤ C1 λ2 E sup I,a 1 M M∑ j=1 εj ( − 1 n t−1∑ τ=ñT hiτ (θj)aτ − λ2 logCI,a
) ≤ 2C1B2
λ2 E  max i∈[n],s∈{±1},ξ:|ξ|≤log(C1) 1 M M∑ j=1 εj ( shi(θj) + λ2 2B2 ξ ) ≤ 4C1B2
λ2
√ 2
M log(2n),
where we used the contraction properties in the first and second inequalities, we used the convex hull argument in the third inequality, and the last inequality is by Massart’s Lemma (Mohri et al., 2012, Theorem 3.3). Summarizing these evaluations, we obtain that
P ( Z ≥ 8C1B2
λ2
√ 2
M log(2n) + C21
√ 1
2M log
( 2n
δ
)) ≤ δ.
The same bound also holds for t = ñT .
Then, by taking the uniform bound over t = ñT, . . . , ñ(T + 1)− 1 and simpifying the bound using ñ ≤ n, we yield that∣∣∣∣∣∣ M∑ j=1 hit(θj)r̄ (t) j − ∫ hit(θ)p̄ (t)(θ)dθ ∣∣∣∣∣∣ ≤ 7C21 √ 1 M log ( 2n δ
) uniformly over t = ñT, . . . , ñ(T + 1) − 1 and the choice of the sequence (iτ )ñ(T+1)−1τ=ñT with probability 1− δ, which yields the assertion.
Lemma 12. With probability 1− δ, it holds that∣∣∣∣∣∣ M∑ j=1 hi(θj) r̂ (t) j∑M j′=1 r̂ (t) j′ − M∑ j=1 hi(θj)r̄ (t) j ∣∣∣∣∣∣ ≤ 7C21 √ 1 M log ( 2n δ ) ,
uniformly over all t = ñT, . . . , ñ(T + 1) − 1, i = 1, . . . , n and the choice of coordinates (iτ ) ñ(T+1)−1 τ=ñT .
Proof. Let Cp = ∫ p(ñT )(θ) exp ( − 1nλ2 ∑t−1 τ=ñT hiτ (θ)δḡ (τ) ) dθ. Then, noticing that r̂(t)j = r̂ (t) j /Cp and ‖hi‖∞ ≤ 1, we have that∣∣∣∣∣∣ M∑ j=1 hi(θj) r̂ (t) j∑M j′=1 r̂ (t) j′ − M∑ j=1 hi(θj)r̄ (t) j
∣∣∣∣∣∣ =
∣∣∣∣∣∣ M∑ j=1 hi(θj) r̂ (t) j∑M j′=1 r̂ (t) j′ − M∑ j=1 hi(θj) r̂ (t) j Cp ∣∣∣∣∣∣ =
∣∣∣∣∣∣ M∑ j=1 hi(θj) r̂ (t) j∑M j′=1 r̂ (t) j′ ∣∣∣∣∣∣ ∣∣∣∣∣1− ∑M j=1 r̂ (t) j Cp ∣∣∣∣∣ ≤
∣∣∣∣∣∣1− M∑ j=1 r̄ (t) j ∣∣∣∣∣∣ . We apply the same argument as Lemma 11 to evaluate the far right hand side
∣∣∣1−∑Mj=1 r̄(t)j ∣∣∣ by replacing hit(θ) with 1 in the definition ofH. Then, we obtain the assertion.
From now on, we define D as
D def = 7
√ 1
M log
( 4n
δ
) .
Then, we have the following bound. Lemma 13. We have that, with probability 1− δ with respect to realization of particles (θj)Mj=1,
|ḡ(t)i − g (t) i | ≤ C2C 2 1 D (45)
for all t = ñT +1, . . . , ñ(T +1), i = 1, . . . , n and the choice of coordinates (iτ ) ñ(T+1)−1 τ=ñT uniformly.
Proof. From Lemma 11 and 12, with probability 1− δ, we have that∣∣∣∣∣∣ M∑ j=1 hit(θj)r̄ (t) j − ∫ hit(θ)p̄ (t)(θ)dθ ∣∣∣∣∣∣ ≤ C21 D, and that ∣∣∣∣∣∣ M∑ j=1 hit(θj) r̂ (t) j∑M j′=1 r̂ (t) j′ − M∑ j=1 hit(θj)r̄ (t) j
∣∣∣∣∣∣ ≤ C21 D, for all t = ñT, . . . , ñ(T + 1)− 1. Therefore, by Lemma 8, we have that
|g(t+1)it − ḡ (t+1) it |
≤ 1 γ + 1nλ2 ∣∣∣∣∣∣ M∑ j=1 hi(θj) r (t) j∑M j′=1 r (t) j′ − ∫ hi(θ)p̄ (t)(θ)dθ ∣∣∣∣∣∣+ 1nλ2 |g(t)i − ḡ(t)i | 
≤ 1 γ + 1nλ2 2C21 D + ∣∣∣∣∣∣ M∑ j=1 hi(θj) ( r (t) j∑M j′=1 r (t) j′ − r̂ (t) j∑M j′=1 r̂ (t) j′ )∣∣∣∣∣∣+ 1nλ2 |g(t)i − ḡ(t)i |  . (46)
The second term of the right hand side of Eq. (46) can be further bounded as∣∣∣∣∣∣ M∑ j=1 hi(θj) ( r (t) j∑M j′=1 r (t) j′ − r̂ (t) j∑M j′=1 r̂ (t) j′ )∣∣∣∣∣∣
≤ M∑ j=1 ∣∣∣r(t)j − r̂(t)j ∣∣∣∑M j′=1 r (t) j′ + M∑ j=1 ∣∣∣∣∣ r̂ (t) j∑M j′=1 r (t) j′ − r̂ (t) j∑M j′=1 r̂ (t) j′ ∣∣∣∣∣ ≤C1
M∑ j=1 ∣∣∣r(t)j − r̂(t)j ∣∣∣+ M∑ j=1 r̂ (t) j∑M j′=1 r̂ (t) j′ ∣∣∣∣∣ ∑M j′=1 r (t) j′ − ∑M j′=1 r̂ (t) j′∑M j′=1 r (t) j′ ∣∣∣∣∣ ≤C1
M∑ j=1 ∣∣∣r(t)j − r̂(t)j ∣∣∣+ C1 ∣∣∣∣∣∣ M∑ j′=1 r (t) j′ − M∑ j′=1 r̂ (t) j′ ∣∣∣∣∣∣ ≤C1
M∑ j=1 ∣∣∣r(t)j − r̂(t)j ∣∣∣+ C1 M∑ j′=1 ∣∣∣r(t)j′ − r̂(t)j′ ∣∣∣ ≤2C1
M∑ j=1 ∣∣∣r(t)j − r̂(t)j ∣∣∣ ≤2C1
M∑ j=1 ∣∣∣∣∣ 1M exp ( − 1 nλ2 t∑ τ=ñT hiτ (θj)δg (τ) iτ ) − 1 M exp ( − 1 nλ2 t∑ τ=ñT hiτ (θj)δḡ (τ) iτ )∣∣∣∣∣ ≤2C 2 1
nλ2 t−1∑ τ=ñT |g(τ+1)iτ − ḡ (τ+1) iτ | (47)
for t = ñT, . . . , ñ(T + 1)− 1. This and Eq. (46) yields
|g(t+1)it − ḡ (t+1) it |
≤ 1 γ + 1nλ2
( 2C21 D + 2C21 + 1
nλ2 t−1∑ τ=ñT |g(τ+1)iτ − ḡ (τ+1) iτ | ) for t = ñT + 1, . . . , ñ(T + 1)− 1. Also, it holds that
|g(ñT+1)i − ḡ (ñT+1) i | ≤ C21 γ + 1nλ2 D,
from Lemma11.
Therefore, applying Gronwall’s lemma (“summation version” of Lemma 18) yields that for t = ñT, . . . , ñ(T + 1)− 1,
|g(t+1)it − ḡ (t+1) it | ≤  2C21 γ + 1nλ2 + 2C21 + 1(
γ + 1nλ2
)2 nλ2 (1 + 2C21 + 1 nγλ2 + 1 )ñ−2 C21 D
≤ 2C 2 1
γ + 1nλ2
( 1 + 2C21 + 1
nγλ2 + 1
)ñ−1 C21 D
≤ 2C 2 1
γ + 1nλ2 exp
( ñ(2C21 + 1)
nγλ2 + 1
) C21 D
≤C2C21 D,
which concludes the proof of Eq. (45).
C.2.3 INTEGRATION OF TWO SOURCES OF DISCRETIZATION ERROR
From now on, we integrate the convergence analysis from Section C.2.1 and Section C.2.2. Then, we extend our analysis from that of each inner loop to that of the outer loop over T = 0, . . . , Tend. So far, we have been arguing only within one inner loop, while we may choose a different number
M of particles in each outer-loop T . Accordingly, we denote the number of particles at the T -th outer-iteration by M (ñT ) and add a subscript to D as
(ñT ) D def = 7
√ 1
M (ñT ) log
( 4n(Tend)
δ
) .
Moreover, we assume that the bound of Lemma 13 holds for every T = 0, . . . , Tend−1. We denote by E this even. We know that P (E) ≥ 1−δ. In the following, we analyze the convergence under the event E where the bounds in Lemma 13 hold, so that |ḡ(t)i −g (t) i | ≤ (C2 (ñT ) C +C2C 2 1 (ñT ) D ) (T = d t ñe−1) holds owing to Lemma 10 and Lemma 13. Therefore, the expectation with respect to the choice of coordinates is taken under the condition of this event. We note that the bound in Lemma 13 holds uniformly over the choice of the coordinates and thus the distribution of the coordinates is not affected by being conditioned by the event.
In order to evaluate the effect of the discrepancy of coordinate updates on the difference of A, we present the following Lemma 14 on the local smoothness of dual space.
Lemma 14. Suppose that n ≥ 2B2λ2 and Assumption 1 holds. Then, D(g) is 1 nλ2 (λ2γ ′ + 1)-smooth with respect to g ∈ Rn on the set A defined by
A = { g ∈ Rn ∣∣∣ min{min t g (t) i ,mint g̃ (t) i } ≤ gi ≤ max { max t g (t) i ,maxt g̃ (t) i } (i = 1, . . . , n) } ,
where t takes t = 0, 1, . . . , ñTend.
Proof. Suppose that the coordinate i is chosen at the t-th iteration. By the definition of g(t), we have that
|`∗i ′(g (t) i )| ≤ M∑ j=1 |hi(θj)| r (t) j∑M j′=1 r (t) j′ + 1 nλ2 |g(t)i − g (t−1) i |, (48)
where `∗i ′ is an element of the sub-differential. The first term of the right hand side can be bounded by 1 as follows:
M∑ j=1 |hi(θj)| r (t) j∑M j′=1 r (t) j′ ≤ max j |hi(θj)| M∑ j=1 r (t) j∑M j′=1 r (t) j′ ≤ 1.
Also, Lemma 6 implies that the second term can be bounded by 2B2nλ2 . Therefore, under the assumption n ≥ 2B2λ2 , we have that
|`∗i ′(g (t) i )| ≤ 2.
Applying the same argument, we also have that
|`∗i ′(g̃ (t) i )| ≤ 2.
Therefore, for i ∈ [n] that has already chosen in the algorithm until the t-th step, we see that `∗i is γ′-smooth on the interval between g(t)i and g̃ (t) i . For the coordinates i that have not yet chosen, we have that g(t)i = g̃ (t) i = g (0) i and |`∗i ′(g (0) i )| ≤ 1. Therefore `∗i is γ′-smooth on the box A. By summarizing this argument and applying Lemma 4, we see that D(g) is 1nλ2 (λ2γ ′ + 1)-smooth.
Combining the above results, we obtain the evalution of A and sufficient conditions for the first statement of Lemma 8. Lemma 15. Suppose that n ≥ 2B2λ2 and fix an arbitrary T . if E[D(g̃
(ñT ′)) − D(g(ñT ′))] ≤ C0 exp (−s̃T ′) is satisfied for all T ′ ≤ T , then we have that
E[D(g̃(t))−D(g(t))] ≤ C0 2 exp (−s̃(T + 1)) + 1 (2− µ)µ
( γ′ + 1
λ2
)( C2 (nT ) C + C2C 2 1 (nT ) D )2 ,
(49)
for all t = ñT + 1, . . . , ñ(T + 1), where µ = 2 1+ 2e s̃
C0
( B1+B2+
C0 1−e−s̃ ) . In particular, if
(ñT ) C ≤ C3 exp
( − s̃T
2
) ,
(ñT ) D ≤ C3 C21 exp
( − s̃T
2 ) are satisfied for all T = 0, . . . , Tend, where
C3 ≤ √ C0(2− µ)µ
C2 √ 2es̃(γ′ + 1/λ2) , (50)
then it holds that for all t = 0, 1, . . . , ñTend,
E[D(g̃(t))−D(g(t))] ≤ C0 exp (−s̃T ) , where T = ⌈ t ñ ⌉ .
Proof. Under the assumptions that n ≥ 2B2λ2 , Lemma 14 holds. Fix some T and we will show the first assertion for t = ñT + 1, . . . , ñ(T + 1).
Let At = { sg̃(t) + (1− s)g(t) ∣∣ s ∈ [0, 1]}, then Lemma 15 implies that D(g) is 1nλ2 (λ2γ′ + 1)- smooth on the set. Hence, for µ = 2
1+ 2e s̃
C0
( B1+B2+
C0 1−e−s̃ ) , the Cauchy-Schwarz inequality yields D(g̃(t))−D(g(t)) ≤ sup g∈At ‖∇D(g)‖2 ‖g̃ (t) − g(t)‖2
≤1 2
[ µ
nλ2 λ2γ′ + 1 ( sup g∈At ‖∇D(g)‖2 )2 + 1 µ λ2γ ′ + 1 nλ2 ‖g̃(t) − g(t)‖22 ]
≤1 2 [ µ sup g∈At [D(g∗)−D(g)] + 1 µ ( γ′ + 1 λ2 )( C2 (ñT ) C + C2C 2 1 (ñT ) D )2] .
If D(g̃(t))−D(g(t)) ≥ 0, then we have that
D(g̃(t))−D(g(t))
≤1 2
[ µ ( D(g∗)−D(g(t)) ) + 1
µ
( γ′ + 1
λ2
)( C2 (ñT ) C + C2C 2 1 (ñT ) D )2] .
Rearranging the terms we have that
D(g̃(t))−D(g(t))
≤ µ 2− µ
( D(g∗)−D(g̃(t)) ) +
1
(2− µ)µ
( γ′ + 1
λ2
)( C2 (ñT ) C + C2C 2 1 (ñT ) D )2 ,
This holds as well when D(g̃(t))−D(g(t)) ≤ 0.
Here, Lemma 5 implies that E[D(g̃(t))] is monotonically non-decreasing within the same inner loop, and thus we have that E [ D(g∗)−D(g̃(t)) ] ≤ E [ D(g∗)−D(g(ñT )) ] . Moreover, the right
hand side E [ D(g∗)−D(g(ñT )) ] can be further bounded by ( B1 +B2 + C0 1−e−s̃ ) exp (−s̃T ) = C0 2 2−µ µ exp(−s̃(T + 1)) by using the assumption and applying Eq. (33) of Lemma 8. Therefore,
E [ D(g̃(t))−D(g(t)) ] ≤C0
2 exp (−s̃(T + 1)) + 1 (2− µ)µ
( γ′ + 1
λ2
)( C2 (ñT ) C + C2C 2 1 (ñT ) D )2 .
This concludes the first assertion.
We will show the second assertion by induction. It is obvious that this is true at initialization. Suppose that Eq. (50) holds for t = 0, 1, . . . ñT . Then the assumption of the first statement of Lemma 8 holds for T ′ = 0, . . . , T, which leads to Eq. (49) for t = ñT + 1, . . . , ñ(T + 1). Then we have that Eq. (50) also holds for t = ñT + 1, . . . , ñ(T + 1), since we chose C3 as
C3 =
√ C0(2− µ)µ
C2 √ 2es̃(γ′ + 1/λ2) ,
so that the right hand side of (49) is bounded by C0 exp (−s̃(T + 1)). Thus, the proof is completed by induction.
Although The proof for the first statement of Lemma 8 is almost done, we have some proof left for the second claim. The next two Lemmas measure the effect of discrepancy of coordinates in the dual space on the value of the main problem. Lemma 16. P (p[g]) is 1
nλ22 (γ + 3λ2 + 9B2)-smooth with respect to g on the set A.
Proof. We show the assertion by a direct calculation. Let Ep[g] [hi] def = ∫ hi(θ)p[g](θ)dθ. First, the second derivative of the loss function can be calculated as ∂2
∂gjgk `i ( Ep[g] [hi] ) = 1
nλ2
∂ ∂gk `i ′ (Ep[g] [hi]) (Ep[g] [hi]Ep[g] [hj ]− Ep[g] [hihj ])
= 1
n2λ22
[ `i ′′ (Ep[g] [hi]) (Ep[g] [hi]Ep[g] [hj ]− Ep[g] [hihj ]) (Ep[g] [hi]Ep[g] [hk]− Ep[g] [hihk])
+ `i ′ (Ep[g] [hi]) · ((Ep[g] [hi]Ep[g] [hk]− Ep[g] [hihk])Ep[g] [hj ]
+Ep[g] [hi] ( Ep[g] [hj ]Ep[g] [hk]− Ep[g] [hjhk] ) − ( Ep[g] [hihj ]Ep[g] [hk]− Ep[g] [hihjhk] ))] ≤ 1 n2λ22 (γ + 3B2) ,
where we used smoothness of `i for the last inequality.
The regularization term can be split into the following two terms:
λ1 ∫ p[g](θ)‖θ‖22dθ + λ2 ∫ p[g](θ) log p[g](θ)dθ
=− 1 n n∑ i=1 ∫ hi(θ)p[g](θ)dθgi + λ2 log [∫ q[g](θ)dθ ] We already confirmed 1nλ2−smoothness of the second term in Lemma 4. As for the first term, we have that
− ∂ 2
∂gjgk
( 1
n n∑ i=1 Ep[g] [hi] gi
)
=− 1 nλ2
∂2
∂gk
( λ2Ep[g] [hj ] + 1
n n∑ i=1 ( Ep[g] [hi]Ep[g] [hj ]− Ep[g] [hihj ] ) gi
)
=− 1 n2λ22
[ 2λ2 ( Ep[g] [hj ]Ep[g] [hk]− Ep[g] [hjhl] ) + 1
n n∑ i=1 (( Ep[g] [hi]Ep[g] [hk]− Ep[g] [hihk] ) Ep[g] [hj ]
+Ep[g] [hi] ( Ep[g] [hj ]Ep[g] [hk]− Ep[g] [hjhk] ) − ( Ep[g] [hihj ]Ep[g] [hk]− Ep[g] [hihjhk] )) gi ]
≤ 1 n2λ22 (2λ2 + 6B2) , where we used Lemma 6 for the last inequality to bound gi.
Combining the two inequalities and the result from Lemma 4, we obtain the assertion.
Lemma 17. Suppose Eq. (50) of Lemma 15 and its assumption hold for t = 0, 1, . . . , ñTend. Then, we have that
E [ P (p[g(t)])− P (p[g̃(t)]) ] ≤2γ + 6λ2 + 18B2
λ 3 2 2 √ γ
√ B1 +B2 + C0 +
C0 1− e−s̃ √ C0es̃(2− µ)µ√ 2(γ′λ2 + 1) exp(−s̃T ),
where T = ⌈ t n ⌉ .
Proof. By using the smoothness of P (p[g]) by Lemma 16, we have that
P (p[g(t)])− P (p[g̃(t)]) ≤ sup g∈At ‖∇P (p[g])‖2 ‖g̃ (t) − g(t)‖2
≤
√ nλ22
γ + 3λ2 + 9B2 { sup g ‖∇P (p[g])‖2 }2√ γ + 3λ2 + 9B2 nλ22 ‖g̃(t) − g(t)‖2
≤ √
sup g P (p[g])− P (p[g∗])
√ γ + 3λ2 + 9B2
λ22
( C2 (ñT ) C + C2C 2 1 (ñT ) C ) .
where T = ⌈ t n ⌉ − 1. if P (p[g(t)])− P (p[g̃(t)]) ≥ 0, then we have that
P (p[g(t)])− P (p[g̃(t)])
≤ √ P (p[g(t)])− P (p[g∗]) √ γ + 3λ2 + 9B2
λ22
( C2 (ñT ) C + C2C 2 1 (ñT ) C ) ≤γ + 3λ2 + 9B2
λ22 √ n
‖g(t) − g∗‖2 ( C2 (ñT ) C + C2C 2 1 (ñT ) C ) .
Since D is γn -strongly convex, we have that‖g (t) − g∗‖2 ≤ √ n γ ( D(g∗)−D(g(t)) ) . Therefore, it holds that
P (p[g(t)])− P (p[g̃(t)])
≤γ + 3λ2 + 9B2 λ22 √ γ
√( D(g∗)−D(g(t)) ) ( C2 (ñT ) C + C2C 2 1 (ñT ) C ) .
which holds as well when P (p[g(t)])− P (p[g̃(t)]) ≤ 0. By taking expectations, it yields that
E [ P (p[g(t)])− P (p[g̃(t)]) ] ≤γ + 3λ2 + 9B2
λ22 √ γ
E [√( D(g∗)−D(g(t)) )] (
C2 (ñT ) C + C2C 2 1 (ñT ) C ) ≤γ + 3λ2 + 9B2
λ22 √ γ
√ E [ D(g∗)−D(g(t)) ] ( C2 (ñT ) C + C2C 2 1 (ñT ) C ) ≤γ + 3λ2 + 9B2
λ22 √ γ
√ E [ D(g∗)−D(g̃(t)) +D(g̃(t))−D(g(t)) ] ( C2 (ñT ) C + C2C 2 1 (ñT ) C ) ≤2γ + 6λ2 + 18B2
λ22 √ γ
√ B1 +B2 +
C0 1− e−s̃ + C0C2C3 exp (−s̃T ) ,
where we used Lemma 8 and Lemma 15 for the last inequality. Furthermore, substituting C2C3 = √ C0(2−µ)µ√
2es̃(γ′+1/λ2) , we have the assertion.
Finally, we will prove the Theorem 3, which gives sufficient conditions about sampling accuracy and number of particles for Theorem 1.
Theorem 3. Suppose that Assumptions 1 and 2 hold. When we choose an Option (A), suppose that we set Tend so that
Tend ≥ 2 n
ñ
( 1 + 1
nγλ2
) log (( n+ 1
γλ2
) (B1 +B2 + C0 1−exp(−s̃) )
P
) , (51)
as in Eq. (30) in Lemma 8. Furthermore, when we choose an Option (B), where we choose a solution randomly from t = ñ(Tend − d 1κe), . . . , ñTend, suppose that we set Tend so that
Tend ≥ 2 n
ñ
( 1 + 1
nγλ2
) log ( Cq P ) + ⌈ 1 κ ⌉ ,
where
Cq = max
{( 2 + 2
nγλ2
)( B1 +B2 +
C0 1− exp (−s̃)
) , (52a)
4C0, (52b)( (8γ+24λ2+72B2) √ B1+B2+C0+
C0 1−e−s̃
√ C0es̃(2−µ)µ
λ 3 2 2 √ 2γ(γ′λ2+1)
)} . (52c)
Then, the condition imposed on (t)A and (t) B in Lemma 8 can be satisfied for both options under the event E if the following inequalities are satisfied: n ≥ 2B2 λ2 , (nT ) C ≤ C3 exp ( − s̃T 2 ) ,
(nT ) D ≤ C3 C21 exp
( − s̃T
2
) ⇔ M (nT ) ≥ 49C 4 1
C23 log ( 4nTend δ ) exp (s̃T ) ,
(53)
for T = 0, 1, . . . , Tend − 1. That is, under the conditions above, E [P (p[ĝ])−D(ĝ))|E ] ≤ P is satisfied for the solution ĝ of either Option A or B.
Proof. As for the Option (A), according to Lemma 15, Eqs. (51) and (53) are the sufficient conditions for the first statement of Lemma 8 in the event E . Next, we will show the convergence of Option (B). In the same way, Eqs. (52a) and (53) assure that (ñT )A is smaller than C0e
−s̃T for every T = 1, . . . , Tend, according to Lemma 15. Furthermore, Eqs. (52b) and (52c) states that both (t)A and (t) B are uniformly bounded by P 4 for all t = ñ(Tend − d 1κe) + 1, . . . , nTend. Note that the evaluation of (t) A follows from Lemma 15, and that of (t) B is derived in Lemma 17.
Therefore, all assertions are proved.
Theorem 3 essentially shows Theorem 1. More detailed explanations can be found in the following.
Proof of Theorem 1. Theorem 3 gives the sufficient conditions for Lemma 8. Thus, we have that, for Option (A), with Eqs. (51) and (53), an expected duality gap E [ P (p[g(ñTend)])−D(g(ñTend))|E ] ≤
P can be achieved, conditioned by the event E whose probability is no less than 1 − δ. In the same way, for Option (B), with Eqs. (52a)-(52c) and (53), an conditional expectation of duality gap E [ P (p[g(t ′ end)])−D(g(t′end))|E ] ≤ P can be achieved. The sufficient condition of tend is then
obtained by noticing tend = ñTend and setting Ĉ1 = C3/7, Ĉ2 = C21 and Ĉ3 = 2Cq
1+ 1nγλ2 .
Next, we evaluate the constants as follows. First, we note that Ĉ2 = C21 by their definitions. Moreover, C−13 = √ 2es̃(γ′+1/λ2)√ C0(2−µ)µ C2 = O ( λ −1/2 2 C 2 1 exp ( ñ(2C21+1) nγλ2+1 )) = O ( λ −1/2 2 Ĉ2 exp ( ñ(2Ĉ2+1) nγλ2+1 )) where we used C2 =
2C21 γ+ 1nλ2 exp ( ñ(2C21+1) nγλ2+1 ) and C21 = Ĉ2. Finally, the evaluation of Ĉ3 can be
obtained by bounding Cq using 1/(1 − exp(−s̃)) = O(s̃−1) and 1 ≤ 1 + 1nγλ2 ≤ 1 + 1 2B2γ = O(1).
C.2.4 DISCRETE GRONWALL’S LEMMA
In the proofs of Lemmas 8, 10 and 13, we used a discrete version of the Gronwall’s lemma. We give its formal statement in the following lemma.
Lemma 18 (Discrete Gronwall’s lemma).
(1) (Vanilla version) Suppose that a sequence of real numbers (un) satisfies
un+1 ≤ aun + bn+1 (∀n ≥ 0),
for a positive real a > 0 and a sequence of real numbers (bn). Then, it holds that
un ≤ anu0 + an n∑ k=1 a−kbk.
(2) (Summation version) Suppose that a sequence of real numbers (un) satisfies
un ≤ n−1∑ k=0 auk +B (∀n ≥ 1),
for a positive real a > 0 and a real number B. Then, it holds that
un ≤ (1 + a)n−1(B + au0).
Proof. The assertion is just an application of the standard discrete version of the Gronwall’s lemma. However, we give a proof for completeness.
The first assertion is a direct consequence of Lemma 5.1 of Mischler (2019). The second assertion is also proven by Lemma 5.2 of Mischler (2019). Indeed, Lemma 5.2 of Mischler (2019) yields that
un ≤ B + n−1∑ k=1 a(1 + a)n−1−kB + a(1 + a)n−1u0.
The right hand side is equivalent to
B [ 1 + a
n−2∑ k=0 (1 + a)k
] + a(1 + a)n−1u0
=B [ 1 + a
(1 + a)n−1 − 1 1 + a− 1
] + a(1 + a)n−1u0
=(1 + a)n−1(B + au0),
which yields the second assertion.

heading: D ADDITIONAL EXPERIMENTS AND DETAILS OF EXPERIMENTAL SETTINGS
text: Here, we give additional experiments and more detailed explanations of our experimental settings. All experiments were conducted under Google Colaboratory. All codes of the experiments in the main text and the following additional experiments are provided in the supplementary material.
Additional experiments for a single neuron teacher network We consider a setting where the teacher network consists of a single neuron:
yi = σ(w ∗> 1 xi) + i,
where the input dimensionality is d = 10, w∗1 ∈ Rd is generated from the uniform distribution on the unit sphere Sd−1 in the d-dimensional Euclidean space, xi ∈ Rd is distributed from a standard normal distribution on Rd, i ∈ R obeys the Gaussian distribution with mean 0 and standard deviation 0.5 and the activation function σ is ReLU. We generated n = 100 trainig data points from this model and trained a student network by using P-SDCA, PDA, and a naive SGD. Here, as the number of particles, we employed M = 200 and the student network is the tanh neural network (hθ(x) = tanh(w>x + b) with θ = (w, b)). We employed the squared loss `(f, y) = (y − u)2/2 for the loss function. We conducted experiments with three patterns of regularization parameters: λ1 = λ2 = λ ∈ {0.01, 0.001, 0.0001}. The naive SGD was trained to minimize the primal objective for the student network with M = 200 without the entropy regularization, i.e.,
1
n n∑ i=1 `
( 1
M M∑ m=1 hθm(xi), yi ) + λ1 M M∑ m=1 ‖θm‖2.
We employed the ULA for sampling the particles for PDA and P-SDCA where we ran the algorithm with 200 steps in each re-sampling stage. The step-sizes η used for sampling (ULA) and SGD were set as 10−3.
We plotted the excess primal objective value against the number of gradient evaluations for each method in Figure 2. The excess primal objective value is computed by subtracting the minimum objective value throughout the optimization and all optimization methods from the current primal objective. Since we are approximating the primal variable by a finite particles, we cannot determine the entropy term exactly. To handle this issue, we employed the k-NN entropy estimator (Kozachenko & Leonenko, 1987; Brodersen, 2020) to estimate the entropy term from the finite particles. This is also applied to SGD to compute the primal objective with the regularization term for the solution obtained by SGD. From Figure 2, we can see that P-SDCA and PDA behave almost identically for large regularization parameter λ. In that regime, SGD once decreases the primal objective rapidly but increases it afterward. This is because SGD is minimizing an objective without the entropy regularization. For λ = 10−2, P-SDCA and PDA does not decrease the excess primal objective below 10−2. This would be due to the instability of the entropy estimator. For small λ, we see that P-SDCA outperforms other methods after some iterations because it yields linear convergence. This observation supports our theoretical analyses.
Additional experiments for a teacher network with multiple neurons We consider a setting where the teacher network consists of multiple neurons:
yi = 1
Mt M∑ m=1 σ(w∗>m xi) + i,
where the width of the teacher network is Mt = 10, the input dimensionality is d = 5, the activation function of the teacher network is σ(·) = tanh(·), the noise is drawn from Normal distribution and
its Signal to Noise Ratio (SNR) is 7.0. For each method, we selected the optimal step size: 10−4 for ULA in P-SDCA and PDA, and 10−5 for SGD. In the following, “step size of P-SDCA or PDA” refers to the step size for ULA in the sampling scheme. In addition to that, we plotted SGD with its step size 10−4 to see how the step size affects the convergence of SGD. The purpose to include SGD as well as P-SDCA and PDA is to see how the entropy term will behave when we do not explicitly include the entropic regularization term into the objective function.
We compared our method with SGD and PDA with the regularization parameters λ1 = 0.01 (fixed) and λ2 = 0.01, 0.001, 0.0001. We repeated the experiments 10 times and took the average of primal objective values, which includes regularization terms. The primal objective value is plotted against the number of gradient evaluations in Figure 3 where the error-bar represents the standard error over the 10 repetitions.
For λ2 = 0.01, we can see that particle methods (P-SDCA and PDA) showed faster convergence than that of SGD. Next, as for λ2 = 0.001, we can see that P-SDCA converged faster than SGD and PDA. Although PDA and SGD with step size η = 10−5 also converged eventually, the behavior of PDA in the intermediate stage was quite unstable and SGD with η = 10−5 showed quite slow convergence. These observations support the superiority of our method. Finally, as for λ2 = 0.0001, SGD with step size η = 10−5 was the fastest due to the small regularization parameter λ2 for the particle methods which makes the convergence of them slower. However, we still see that P-SDCA outperformed PDA in terms of both convergence speed and stability. The convergence of SGD is strongly affected by its step size. In fact, SGD with larger step size (η = 10−4) could decrease the objective rapidly at the early stage, but it eventually diverged. This would be partly for the same reason as in the single neuron setting, that is, due to the lack of the entropic regularization in the objective of SGD. In summary, we can see that P-SDCA properly minimizes the corresponding objective function in wide range of λ2 and shows faster and more stable convergence than PDA as indicated by our theory.
Experiments on MNIST dataset We conduct numerical experiments to illustrate the “feature learning” aspect of mean field neural network. For this purpose, we run the algorithm in a binary classification task that separates “2” and “4” in the MNIST dataset. We subsampled n = 1000 training examples and trained the two layer tanh network with width M = 2500. The step size for ULA in P-SDCA was set as η = 10−5 and the regularization parameters were set as λ1 = 10−2, λ2 = 10−4.
To visualize the adaptivity of the feature map, we first plot the evolution of singular values of the first layer’s weight matrix, W = [w1, w2, . . . , wM ] ∈ Rd×M , until convergence, when the solution achieves 99.8% training accuracy. We can see that the singular values of the network trained by P-SDCA drastically changed during optimization. This indicates that the mean-field neural network does not ”freeze” and can adaptively learn features during optimization by P-SDCA. We also plotted the same quantity for SGD updates in the NTK regime. In contrast, we observe only small change of the singular values. This is because the first layer’s parameters do not change much during optimization in the NTK regime; in other words, SGD in the NTK regime does not perform adaptive feature extraction. In that sense, we expect that training in the mean field regime is effective in a setting where adaptive feature extraction is required (e.g., see Daniely & Malach (2020)).
In addition to the singular value evolution, we evaluate how well the extracted features are “aligned” to the label. For that purpose, we define a kernel function kW (x, x′) = ∑M m=1 σ(w > mx)σ(w > mx ′) for a weight matrix W = [w1, w2, . . . , wM ] ∈ Rd×M . Then, we define the kernel alignment (Cristianini
et al., 2002) on the training data (xi, yi)ni=1 as
A(kW ) := 〈KW , yy>〉F√
〈KW ,KW 〉F〈yy>, yy>〉F where 〈A,B〉F := ∑n i=1 ∑n j=1AijBij , KW = (kW (xi, xj) n,n i=1,j=1) and y = (y1, . . . , yn)
>. We also define the kernel alignment on the test data in the same manner. We can see that the kernel alignment represents how strongly the kernel function defined by the features after the first layer is aligned to the target signal.
The kernel alignment of the solutions obtained by each method (P-SDCA, SGD in NTK regime, SGD in mean field regime) is depicted in Figure 5. We can see that P-SDCA (and SGD in the mean field regime) properly improves the kernel alignment, which means that P-SDCA can extract informative features adaptively depending on the data. On the other hand, SGD in the NTK regime does not increase the kernel alignment due to the parameters almost “frozen” at initialization. This experiment highlights how well the mean field neural networks can execute feature learning.
Evaluation of effect of the negative entropy term Here, we evaluate the effect of the negative entropy in a similar setting to Appendix F.3 of Nitanda et al. (2021). We consider a synthetic 1D dataset with n = 15. Figure 6 plots the estimated curve by a two-layer tanh network with M = 2000 neurons trained by P-SDCA and SGD where the squared loss is employed. We used the same L2-regularization parameter λ1 = 10−3 for both methods, and employed λ2 = 10−3 for the negative entropy regularization in P-SDCA. We observe that the solution by P-SDCA is less complicated than that by SGD. This indicates that the negative entropy regularization can serve as a smoothness penalty, which may be useful to avoid overfitting in a noisy dataset. Similar findings have also been reported in (Nitanda et al., 2021).

Review summaries:
The paper considers optimization of MF two-layer neural nets (the infinite-width version with entropic regularization). Specifically the paper proposes to do dual coordinate ascent, which allows for exponential convergence rate, together with a particle approximation scheme. The paper shows that in the discrete-time setting, an improved convergence rate is obtained (w.r.t. desired error $\epsilon_P$, the number of data points $n$, the regularization $\lambda_2$).

The problem of training an infinitely wide neural network with one hidden layer can be written as an optimization problem in the space of probability distributions, where each atom in the distribution corresponds to a hidden neuron. The paper proposes to solve instead the dual problem (of an entropy regularized variant), which in the case of finite data is a finite-dimensional optimization problem.  The main contribution is the conception and analysis of an efficient dual coordinate ascent scheme to solve the dual, which in theory and practice requires fewer iterations to find an epsilon-accurate solution than existing methods.  

This work solves the entropy regularized mean-field model for a two-layer neural network. As the Fenchel–Rockafellar dual of that problem is a finite dimension problem, they consider an SDCA-type scheme to solve the dual problem, which still contains an integral term and they use the Langevin iterate procedure to solve the subproblem approximately. In contrast to the existing analysis for SGD in the mean-field regime, the new algorithm has a fixed length dual step size and has a fast linear convergence rate.

This paper presents an optimization alogrithm for mean-field shallow neural networks. The algorithm optimizes the parameter measure using a particle-based implementation of the stochastic dual coordinate ascent algorithm. The authors establish convergence results, namely an exponential convergence rate owing to the convexity of the dual problem.