Convergence and complexity of block coordinate descent with diminishing radius for nonconvex optimization

Block coordinate descent (BCD), also known as nonlinear Gauss-Seidel, is a simple iterative algorithm for nonconvex optimization that sequentially minimizes the objective function in each block coordinate while the other coordinates are held fixed. We propose a version of BCD that, for block multi-convex and smooth objective functions under constraints, is guaranteed to converge to the stationary points with worst-case rate of convergence of $O((\log n)^{2}/n)$ for $n$ iterations, and a bound of $O(\epsilon^{-1}(\log \epsilon^{-1})^{2})$ for the number of iterations to achieve an $\epsilon$-approximate stationary point. Furthermore, we show that these results continue to hold even when the convex sub-problems are inexactly solved if the optimality gaps are uniformly summable against initialization. A key idea is to restrict the parameter search within a diminishing radius to promote stability of iterates. As an application, we provide an alternating least squares algorithm with diminishing radius for nonnegative CP tensor decomposition that converges to the stationary points of the reconstruction error with the same robust worst-case convergence rate and complexity bounds. We also experimentally validate our results with both synthetic and real-world data and demonstrate that using auxiliary search radius restriction can in fact improve the rate of convergence.

INTRODUCTION
Throughout this paper, we are interested in the minimization of a continuous function f : R I 1 ×· · ·× R I m → [0, ∞) on a cartesian product of convex sets Θ = Θ (1) × · · · × Θ (m) : θ * ∈ arg min θ=[θ 1 ,...,θ m ]∈Θ f (θ 1 , . . . , θ m ). (1) We are interested in the case when f is block multi-convex, that is, when f is convex in each block coordinate θ i . When the objective function f is nonconvex, the convergence of any algorithm for solving (1) to a globally optimal solution can hardly be expected. Instead, global convergence to stationary points of the objective function is desired.
Block coordinate descent (BCD), also known as nonlinear Gauss-Seidel [Ber99,Wri15], is a simple iterative algorithm that sequentially minimizes the objective function in each block coordinate while the other coordinates are held fixed. Due to its simplicity, BCD has been widely used in various optimization problems such as nonnegative matrix or tensor factorization [LS99,LS01,KB09]. More precisely, fix an initial estimate θ 0 ∈ Θ, and obtain a new estimate θ n from the previous one θ n−1 by successively updating the i th block coordinate using the following rule for i = 1, . . . , m: Then (2) generates a sequence of esimates θ n = [θ (1) n , . . . , θ (m) n ] for n ≥ 1. Convergence properties of BCD for convex f has been extensively studied [H + 57, SS73,LT92,Tse91]. It is known that BCD does not always converge to the stationary points of the non-convex objective function that is convex in each block coordinate [Pow73], but such global convergenece is guaranteed under additional assumptions: Two-block (m = 2) or strict quasiconvexity for m − 2 blocks [GS99,GS00] and uniqueness of minimizer per block [Ber97,Sec. 2.7].
One way to guarantee convergence of BCD without requiring additional reguality assumptions is to use proximal regularization, which in the simplest setting takes the following form (see [GS00,XY13, The codes for the main algorithm and simulations are provided in https://github.com/HanbaekLyu/BCD-DR.
While there exists an extensive amount of work on the rate of convergence of BCD when the objective function f is convex on the whole domain, the rate of convergence in the general block multiconvex case is not very well understood. Namely, various sublinear bounds on the rate of convergence is known for convex problems assuming local error bounds [LT93], projected gradient descent [BT13], and randomized and accelerated coordinate descent [Nes12].
In this work, we introduce a new modification of the standard BCD (2) and achieve the following results: (1) Global convergence to stationary points for block multi-convex and differentiable objectives under constraints; (2) Worst-case rate of convergence to stationary points of order (log n) 2 /n; (3) Worst-case bound of O(ε −1 (log ε −1 ) 2 ) on the number of iterations to achieve ε-approximate stationary points; (4) Robustness of the aformentioned results under inexact execution of the algorithm.
To our best knowledge, we believe our work provides the first result on the global rate of convergence and iteration complexity of BCD for the general block multi-convex case, especially with the additional robustness result. We do not impose any additional assumptions such as the Kurdyka-Łojasiewicz property. For gradient descent methods with unconstrained nonconvex objective, it is known that such rate of convergence cannot be faster than O(1/n) [CGT10], so our rate bound matches the optimal result up to a (log n) 2 factor. We experimentally validate our results with both synthetic and real-world data and demonstrate that using auxiliary search radius restriction can in fact improve the rate of convergence over the standard BCD (2). However, we have not proved that our algorithm converges faster than the standard BCD in general.

3:
Update estimate θ n = [θ (1) n , · · · , θ (m) n ] by 4: For i = 1, · · · , m do: Our key idea is to use an additional (but rather natural) constraint of "diminishing radius", by which we mean restricting the range of parameter search within a diminishing radius from the previous estimation. Namely, instead of adding the proximal regularizer θ − θ (i ) n−1 2 F to the objective function, we impose an additional constraint θ − θ (i ) n−1 F ≤ c w n , where c > 0 is a fixed constant and (w n ) n≥1 is a non-increasing sequence of weights in (0, 1]. See Algorithm 1. As pointed out in [Ber99,Prop. 2.7.1], a necessary property that implies convergence of BCD to stationary points is the "stability of estimates", that the change made in each block update (2) vanishes. Our diminishing radius condition directly bakes such property in the algorithm in the form of an auxiliary constraint. Moreover, we remark that using proximal regularization in (3) only ensures square-summability of iterates, ∞ n=1 θ n − θ n−1 2 F < ∞ (see, e.g., [XY13, Lem 2.2]), while our diminishing radius constraint gives a finer control of individual iterates, θ n − θ n−1 F ≤ c w n . This advantage is crucial in obtaining the rate of convergence in our analysis. As a trade-off, we devote a substantial amount of work to show the auxiliary diminishing radius constraint vanishes along every convergent subsequence.

STATEMENT OF MAIN RESULTS
Throughout this paper, we assume the following two mild conditions: In (A1), the assumption that each Θ (i ) ⊆ R I i contains an open ball in R I i does not lose any generality. Indeed, otherwise, Θ (i ) is contained in some hyperplane and after a change of coordinates, we can reduce the full dimension of the ambient spcae R I i . The smoothness assumption is standard on the literature of BCD (see, e.g., [XY13,LX17]). It is important to note that under (A1)-(A2), the subproblem of block minimization (4) is a convex minimization problem over a convex set n F ≤ c w n+1 } (being the intersection of two convex sets), so Algorithm 1 can be readily executed using standard convex optimization procedures (see, e.g., [Ber99]). (We do cover the case when these convex sub-problems are only approximately solved in Theorem 2.2.) Recall that we say θ * ∈ Θ is a stationary point of f over Θ if where 〈·, ·〉 denotes the dot project on R I 1 +···+I m ⊇ Θ. This is equivalent to saying that −∇ f (θ * ) is in the normal cone of Θ at θ * . If θ * is in the interior of Θ, then it implies ∇ f (θ * ) F = 0. For iterative algorithms, such a first-order optimality condition may hardly be satisfied exactly in a finite number of iterations, so it is more important to know how the worst-case number of iterations required to achieve an ε-approximate solution scales with the desired precision ε. More precisely, we say θ * ∈ Θ is an ε-approxiate stationary point of f over Θ if This notion of ε-approximate solution is consistant with the corresponding notion for unconstrained problems. Indeed, if θ * is an interior point of Θ, then (5) reduces to ∇ f (θ * ) 2 F ≤ ε. Next, for each ε > 0 we define the worst-case iteration complexity N ε of an algorithm for solving (1) as where (θ n ) n≥0 is a sequence of estimates produced by the algorithm with initial estimate θ 0 . Note that N ε gives the worst-case bound on the number of iterations for an algorithm to achieve an εapproximate solution due to the supremum over the initialization θ 0 in (6). Now we state the first main result, Theorem 2.1. To our best knowledge, this gives the first worstcase rate of convergence and iteration complexity of BCD-type algorithms in the literature.
It is worth noting that, imposing an aditional assumption that all stationary points of f of Θ are isolated, then global convergence to a single stationary point can be easily deduced from Theorem 2.1 (i). However, requiring global convergence to a single stationary point in the general case (especially when the stationary point may accumulate) seems to be too strong. Indeed, the standard BCD does not warrant such a strong global convergence (e.g., when the objective function is globally minimized at every point in a small ball of factors). For our algorithm, θ n will converge to a single point trivially if the radii are summable (i.e., ∞ n=1 w n < ∞). However, in this case convergence to stationary points is not guaranteed in general, as all stationary points could be farther than the sum of radii from initialization.
Theorem 2.1 (ii) gives a worst-case rate of convergence of Algorithm 1. Namely, for each n ≥ 1, uniformly among all initialization θ 0 , there exists k ∈ {1, . . . , n} such that θ k is an ε-approximate solution with ε = c/ n i =1 w i . If we consider the weight w n = n −γ for some γ > 0, then the hypothesis of Theorem 2.1 requires that γ ∈ (1/2, 1], and the bound in (7) is optimized when γ 1/2. More precisely, consider the weights w n = 1/( n log n). Then Theorem 2.1 holds and we have for some constant c > 0 for all n ≥ 1. In particular, if all stationary points of f in Θ are in the interior of Θ, then (7) implies Hence we obtain the worst-case rate of convergence in gradient norm squared of order O((log n) 2 /n). We remark that in the literature of stochastic gradient descent for unconstrained problems, the rate of convergence is often measured by the worst-case gradient norm squared [SSY18, XYY + 19, WWB19,KW18]. Hence we can interpret (8) as giving a worst-case rate of convergence of order O((log n) 2 /n) for constrained nonconvex problems.
Our second set of results in this work considers robustness of our former result in Theorem 2.1 when the convex sub-problems in Algorithm 1 are solved only approximately. We show that Theorem 2.1 continues to hold under inexact compution of (4) if the errors we make in solving each convex sub-problem are "uniformly summable". To make a precise statement, we introduce the following notion of "optimality gap". Let (θ n ) n≥0 be an inexct output of Algorithm 1. Writeθ n = [θ (1) n , . . . ,θ (m) n ] for n ≥ 1 and denote the i th block-marginal cost function for i = 1, . . . , m and n ≥ 1. Eachθ (i ) n would be a minimizer of f (i ) n over the convex set Θ (i ) if (θ n ) n≥0 were an exact output of Algorithm 1. Hence we define the optimality gap at iteration n as Theorem 2.2. Assume (A1)-(A2). Suppose weights (w n ) n≥1 satisfy ∞ n=1 w 2 n < ∞. Let (θ n ) n≥0 be an inexct output of Algorithm 1. Then the following hold: F , thenθ n converges to the set of stationary points of f over Θ.
Note that since the sub-problems (4) are convex so that they can be solved approximately up to arbitrary prescribed precision. Hence summability of the optimality gap can easily be satisfied. For instance, the optimality gap shrinks linearly fast if the sub-problems (4) are strictly convex by using coordinate descent algorithms or using gradient descent algorithms, where also superlinear convergence is available using accelerated methods under additional assumptions [Wri15,BCN18].
Lastly, we apply our method and main results to Nonnegative CANDECMP/PARAFAC Decomposition (NCPD) [SH05,Zaf09]. The usual Alternating Least Squares (ALS) algorithm for solving NCPD [KB09] is the vanila BCD (2) for NCPD, which does not have a general convergence guarantee when the number of modes exceeds two. With additional diminishing raidus condition, the ALS is modified as Algorithm 2, which now enjoys all theoretical convergence guarentees and bounds on the rate of convergence and complexity stated in Theorems 2.1. Moreover, by Theorem 2.2 these results are robust under inexact computation of each nonnegative least squares computation. These results are stated in Corollary 5.1 and are new are new to the literature on NCPD.
2.1. Orgnization. This paper is organized as follows. In Section 3, we prove Theorem 2.1. Modifying the arguments in the exact case, we prove Theorem 2.2 in Section 4. Then we apply our theory to NCPD and drive Corollary 5.1 in Section 5. Finally, we provide some experimental results for NCPD problems in Section 6.
We start by recalling a classical lemma on the first-order approximation of functions with Lipshitz gradients.
Lemma 3.1 (First-order approximation of functions with Lipschitz gradient). Let f : R p → R be differentiable and ∇ f be L-Lipschitz continuous. Then for each θ, θ ∈ R p , Proof. This is a classical lemma. See [Nes98, Lem 1.2.3].

Proposition 3.2 (Forward monotonicity). Suppose (A1) and (A2). Then f
Recall that θ (i ) n is a minimizer of the convex function f (i ) n over the convex set Θ (i ) n that contains θ (i ) n−1 (see (11)). Hence each summand in the last expression above is nonnegative. This shows f (θ n−1 ) − f (θ n ) ≥ 0, as desired.

Remark 3.4. Proposition 3.3 still holds if we replace the assumption (A1) with uniform boundedness of f over Θ.
The following is a simple observation based on (A1).
Proof. Let R min denote the minimum radius of all open balls in Θ (i ) for i = 1, . . . , m (recall (A1)). We will show the assertion holds for Next, we show that the block coordinate descent we use to obtain θ n+1 should always give the optimal first-order descent up to a small additive error. Proposition 3.6 (Asymptotic first-order optimality). Assume (A1)-(A2) and w n = o(1). Fix a sequence (b n ) n≥1 such that 0 < b n ≤ w n for all n ≥ 1. Then there exists constant c 1 > 0 that do not depend on θ 0 ∈ Θ such that for all n ≥ 1, Proof. Let R 0 > 0 denote the constant as in Proposition 3.5. Since w n = o(1), there exists N 0 ≥ 1 such that c w n ≤ R 0 for all n ≥ N 0 . Such N 0 does not depend on θ 0 . We will prove the assertion for all n ≥ N 0 . Since all quantities in the assertion up to iteration N 0 are uniformly bounded (using compactness in (A1)), the assertion for all n ≥ 1 then follows. Recall which holds for all θ ∈ Θ with θ − θ n F = c b n+1 . Adding up these inequalities for i = 1, . . . , m, we get Noting that Θ is convex by (A1), we see that (13) also holds for all θ ∈ Θ with θ − θ n F ≥ c b n+1 . Indeed, for such θ, let θ be the point in the secant line between θ and θ n such that θ −θ n F ≥ c b n+1 . Then θ ∈ Θ and the above inequality holds for θ replaced with θ . However, the right hand side is unchanged when replacing θ with any point on the line passing through θ and θ n .
To finish the argument, it remains to argue that (13) holds for all θ ∈ Θ for all n ≥ N 0 with possibly different constant c 1 > 0. Let ∂Θ denote the boundary of Θ. Since Θ is convex and (13) is invariant under radial scaling of the vector θ − θ n , it suffices to verify (13) for boundary points of Θ that are within c b n+1 from θ n . To this end, fix θ ∈ ∂Θ with θ −θ n F < c b n+1 . Using Proposition 3.5, choose θ n ∈ Θ such that θ n − θ n F = c b n+1 . Now by a change of coordinates, we can view the graph of the boundary ∂Θ near θ n as a convex function. Then by a Taylor expansion, there exists a constant c 2 > 0 independent of θ 0 ∈ Θ such that for all n ≥ 1, Note that the constant c 2 above can be chosen to be independent of θ . Since (13) holds for θ = θ and ∇ f is L-Lipschitz, combining (13) and (14) gives Since this holds for all θ ∈ ∂Θ with θ − θ n F < c b n+1 , this completes the proof.
Recall that during the update θ n−1 → θ n each block coordinate of θ n−1 changes by at most c w n in Frobenius norm. For each n ≥ 1, we say θ n is a long point if none of the block coordinates of θ n−1 change by c w n in Frobenius norm and a short point otherwise. Observe that if θ n is a long point, then imposing the search radius restriction in (4) has no effect and θ n is obtained from θ n−1 by a single cycle of block coordinate descent on f over Θ.
Proof. By Proposition 3.3, we have Hence the former condition implies the latter condition in (21). Thus it suffices to show that this latter condition implies the assertion. Assume this condition, and let (t k ) k≥1 be a subsequence of (n k ) k≥1 for which the liminf in (21) is achieved. By taking a further subsequence, we may assume that θ ∞ = lim k→∞ θ t k exists. Now suppose for contradiction that θ ∞ is not a stationary point of f over Θ. Then there exists θ ∈ Θ and δ > 0 such that Noting that θ n −θ n−1 F = O(w n ) = o(1), we see that the right hand side goes to zero as k → ∞. Hence for all sufficiently large k ≥ 1, we have Noting that θ n+1 − θ n F ≤ mc w n+1 , applying Proposition 3.6 with b n = θ n+1 − θ n F /mc gives which is a contradiction. This shows the assertion.
The following proposition is crucial in establishing global convergence to stationary poitns (Theorem 2.1 (i)) and it is the only place where we use non-summability of the weights w n in the proof of Theorem 2.1. Proof. Let A denote the event that there exists a non-stationary limit point θ ∞ of Λ. Let B denote the event that there exists an ε-neighborhood B ε (θ ∞ ) of θ ∞ that does not contain any long points of Λ. We claim that A implies B . Suppose for contrary that A \ B is true. Then for all ε > 0, B ε (θ ∞ ) contains some long point of Λ. Then there exists a sequence of long points of A converging to θ ∞ , but then θ ∞ has to be stationary by Proposition 3.7. This contradicts to A. Hence A implies B .
Lastly, from the earlier results, on the event A, we can choose ε > 0 such that B ε (θ ∞ ) has no long point of Λ and also satisfies (a). We will show that B ε/2 (θ ∞ ) satisfies (b). Then B ε/2 (θ ∞ ) satisfies (a)-(b), as desired. Suppose for contradiction there are only finitely many θ n 's outside of B ε/2 (θ ∞ ). Then there exists an integer M ≥ 1 such that θ n ∈ B ε/2 (θ ∞ ) for all t ≥ M . Then each θ n for t ≥ M is a short point of Λ. By definition, it follows that θ n −θ n−1 F ≥ c w n for all t ≥ M . Then since ∞ n=1 w n = ∞, by Proposition 3.8 there exists a subsequence (s k ) k≥1 such that θ ∞ := lim k→∞ θ n k exists and is stationary. But since θ ∞ ∈ B ε (θ), this contradicts (a) for B ε (θ). This shows the assertion.
Lemma 3.10. Let (a n ) n≥0 and (b n ) n≥0 be sequences of nonnegative real numbers such that ∞ n=0 a n b n < ∞. Then Proof. The assertion follows from noting that n k=1 a k min We are now ready to give a proof of Theorem 2.1.
Proof of Theorem 2.1. Assume (A1), (A2) and ∞ n=1 w 2 n < ∞. We first show (ii). First note that by Proposition 3.6 with b k = w k , we have for some constants c 1 , c 2 > 0 for all k ≥ 1. From Proposition 3.3 and since ∞ n=1 w 2 n < ∞, for some constant M > 0 that does not depend on the initialization θ 0 . Thus by Lemma 3.10, we have Hence by (17), we conclude Since M does not depend on θ 0 , we can take the supremum over θ 0 ∈ Θ on both sides above and finish the proof of (ii). Next, (iii) is a direct consequence of (ii). Indeed, if w n = 1/( n log n), then the uniform upper bound on the rate of convergence in (ii) is of order O((log n) 2 /n). Then one can conclude by using the fact that n ≥ 2ε −1 (log ε −1 ) 2 implies (log n) 2 /n ≤ ε for all sufficiently small ε > 0.
Finally, we assume ∞ n=1 w n = ∞ and show (i). Suppose for contradiction that there exists a nonstationary limit point θ ∞ ∈ Θ of Λ. By Proposition 3.9, we may choose ε > 0 such that B ε (θ ∞ ) satisfies the conditions (a)-(b) of Proposition 3.9. Choose M ≥ 1 large enough so that w n < ε/4 whenever t ≥ M . We call an integer interval I := [ , ] a crossing if θ ∈ B ε/3 (θ ∞ ), θ ∈ B 2ε/3 (θ ∞ ), and no proper subset of I satisfies both of these conditions. By definition, two distinct crossings have empty intersection. Fix a crossing I = [ , ], it follows that by triangle inequality, Note that since θ ∞ is a limit point of Λ, θ n visits B ε/3 (θ ∞ ) infinitely often. Moreover, by condition (a) of Proposition 3.9, θ n also exits B ε (θ ∞ ) infinitely often. It follows that there are infinitely many crossings. Let t k denote the k th smallest integer that appears in some crossing. Then t k → ∞ as k → ∞, and by (18), Then by Proposition 3.8, there exists a subsequence (s k ) k≥1 of (t k ) k≥1 such that θ ∞ := lim k→∞ θ s k exists and is a stationary of f over Θ. However, since θ t k ∈ B 2ε/3 (θ ∞ ) for all k ≥ 1, we have θ ∞ ∈ B ε (θ ∞ ). This contradicts the condition (b) of Proposition 3.9 for B ε (θ ∞ ) that it cannot contain any stationary point of f over Θ. This shows that there is no non-stationary limit point θ ∞ ∈ Θ of Λ.

PROOF OF THEOREM 2.2: THE INEXACT CASE
In this section, we extend our analysis in the main text to the case where the convex sub-problems (4) in the main algorithm (Algorithm 1) is solved inexactly, thereby establishing Theorem 2.2. Throughout this section, (θ n ) n≥1 is assumed to be an inexact solution of Algorithm 1 and (∆) n≥1 denote the associated optimality gap defined in (10). We writeθ n = [θ (1) n , . . . ,θ (m) n ] for each n ≥ 1 and let f (i ) n denote the i th block marginal cost function at iteration n as defined in (12). Also, we define Λ and Θ (i ) n similarly as in the beginning of Section 3. Namely, Λ := {θ n | t ≥ 1} ⊆ Θ and for each n ≥ 1 and i ∈ {1, . . . , m}, we denote In order to prove Theorem 2.2, we will proceed by making modifications of the statements in the main text we proved for the exact trajectory of Algorithm 1. We first need the following approximate monotonicity result.
Let θ (i ) n be a minimizer of the convex function θ → f (i ) n (θ) over the convex set Θ (i ) n . By the definition of the optimality gap ∆ n in (10), we have Hence we have Then the assertion follows by summing over the above inequality for i = 1, . . . , m.

Remark 4.2. We remark that f
n is an estimate of the minimizer of f (i ) n over Θ (i ) n , which is a ball in Θ (i ) centered atθ (i ) n−1 , this is an entirely trivial condition to satisfy, e.g., by using a line search in each iteration. With this additional assumption, the numerator M + 2m ∞ n=1 ∆ n in the upper bound in Theorem 2.2 can be reduced to M + m ∞ n=1 ∆ n , see the proof of Theorem 2.2.
Proof. According to (A1) and (A2), it follows that ∇ f overθ is Lipschitz with some uniform Lipshitz constant L > 0. Hence by Lemma 3.1, for all t ≥ 1, Using Proposition 4.1, it follows that for n ≥ 1. Note that θ n −θ n+1 Then we can conclude since the first two terms in the middle of the equations above can be uniformly bounded from above with respect to θ 0 .
Next, we show an analogue of Proposition 3.6 in the inexact case.
Then it follows that where the first inequality above follows by the definition of the optimality gap ∆ n+1 in (10). Write this as Recall is Lipschitz with uniform Lipschitz constant L > 0 by (A2). Hence by Lemma 3.1, we get where ∇ 1 denotes the partial gradient with respect to the first block-coordinate. Applying a similar argument for all other blocks, it follows that ∇ f (θ ∞ ), (θ −θ ∞ ) ≥ 0 for all θ ∈ Θ. This shows the assertion.
Hence the former condition implies the latter condition in (21). Thus it suffices to show that this latter condition implies the assertion. Assume this condition, and let (t k ) k≥1 be a subsequence of (n k ) k≥1 for which the liminf in (21) is achieved. By taking a further subsequence, we may assume that θ ∞ = lim k→∞θt k exists. Now suppose for contradiction thatθ ∞ is not a stationary point of f over Θ. Then there exists θ ∈ Θ and δ > 0 such that ∇ f (θ ∞ ), (θ −θ ∞ ) < −δ < 0. By triangle inequality, write Noting that θ n −θ n−1 F = O(w n ) = o(1), we see that the right hand side goes to zero as k → ∞. Hence for all sufficiently large k ≥ 1, we have ∇ f (θ t k ), (θ −θ t k ) < −δ /2.
Noting that θ n+1 −θ n F ≤ mc w n+1 , applying Proposition 4.4 with b n = θ n+1 −θ n F /mc , for some constants c 1 , c 2 > 0 for all n ≥ 1. Note that here we have used crucially the assumption The rest of the proof is similar to the proof of Proposition 3.8. Proof. Proof is similar to that of Proposition 3.9 using the inexact counterparts of the earlier results shown in this section.
We now prove Theorem 2.2.
Lastly, we assume ∞ n=1 w n = ∞ and show (i). For this, suppose ∆ n = O( θ n −θ n−1 2 F ). (Since θ n −θ n−1 F = O(w n ) and ∞ n=1 w 2 n < ∞, this is a stronger condition than ∞ n=1 ∆ n < ∞.) Under this assumption, the exact counterparts of Proposition 3.9, Proposition 4.7, holds. Then an identical argument as in the proof of Theorem 2.1 shows thatθ n converges to the set of all stationary points of f over Θ.

APPLICATIONS TO CONSTRAINED MATRIX AND TENSOR FACTORIZATION
In this section, we discuss applications of our main result in the setting of nonnegative matrix and tensor factorization problems [LS99,LS01,KB09]. Nonnegative matrix factorization (NMF) has found applications in text analysis, image reconstruction, medical imaging, bioinformatics, and many other scientific fields more generally [SGH02, BB05, BBL + 07, CWS + 11, TN12, BMB + 15, RPZ + 18]. As matrix factorization is for unimodal data, nonnegative tensor factorization (NTF) provides a powerful and versatile tool that can extract useful latent information out of multi-model data tensors. As a result, tensor factorization methods have witnessed increasing popularity and adoption in modern data science [SH05, Zaf09, ZVB + 16, SLLC17,RLH20]. The use of nonnegativity constraint in both NMF and NTF is crucial in obtaining a "parts-based" representation of the input signals [LS99]. Given a d × T data matrix X ∈ R d ×T and an integer parameter r ≥ 1, consider the following constrained matrix factorization problem arg min where the two factors W and H are called dictionary and code matrices of X , respectively, and λ ≥ 0 is a 1 -regularizer for the code matrix H . Depending on the application contexts, we may impose some constraints Θ (1) and Θ (2) on the dictionary and code matrices, respectively, such as nonnegativity or some other convex constraints. An interpretation of this approximate factorization is that the r columns of W give an approximate basis for spanning the n columns of X , where the columns of H give suitable linear coefficients for each approximation. In dictionary learning problems [OF97, EAH99, LS00, EA06, LHK05] , one seeks for a sparse representation of the columns of X with respect to and an over-complete dictionary W , for which one can take r > d and λ > 0.
More generally, we can formulate a similar constrained factorization problem in a tensor setting. Suppose a data tensor X ∈ R I 1 ×···×I m is given and fix an integer R ≥ 1. In the CANDECOMP/PARAFAC (CP) decomposition of X [KB09], we would like to find R loading matrices U (i ) ∈ R I i ×R for i = 1, . . . , m such that the sum of the outer products of their respective columns approximate X: where U (i ) [:, k] denotes the k th column of the I i × R loading matrix matrix U (i ) and denotes the outer product. As an optimization problem, the above CP decomposition model can be formulated as the following the constrained tensor factorization problem: where Θ (i ) ⊆ R I i ×R denotes a compact and convex constraint set and λ i ≥ 0 is a 1 -regularizer for the i th loading matrix U (i ) for i = 1, . . . , m. In particular, by taking λ i = 0 and Θ (i ) to be the set of nonnegative I i × R matrices with bounded norm for i = 1, . . . , m, (24) reduces to the nonnegative CP decomposition (NCPD) [SH05,Zaf09]. Also, It is easy to see that (24) is equivalent to arg min U (1) ∈Θ (1) ,...,U (m) ∈Θ (i ) X − Out(U (1) , . . . ,U (m−1) ) × m (U (m) ) T 2 which is a CP-dictionary-learning problem introduced in [SLN20]. Here × m denotes the mode-m product (see [KB09]) the outer product of loading matrices U (1) , . . . ,U (m) is defined as Namely, we can think of the m-mode tensor X as I m observations of (m − 1)-mode tensors, and the R rank-1 tensors in Out(U (1) , . . . ,U (m) ) serve as dictionary atoms, whereas the tranpose of the last loading matrix U (m) can be regarded as the code matrix (see Figure 1). In particular, assuming m = 2, (25) becomes the constrained matrix factorization problem in (23).
The constrained tensor factorization problem (24) falls under the framework of BCD in (1), where the objective function f CTF is convex in each loading matrix U (i ) for i = 1, . . . , m. Indeed, BCD is a popular approach for both NMF and NTF problems [KHP14]. Namely, when we apply BCD (2) for (24), each block update amounts to solve a quadratic problem under convex constraint. BCD for (24) is known as the form of alternating least squares (ALS). For NMF, ALS (or vanila BCD in (2)) is known to converge to stationary points [GS00]. However, for NTF with m ≥ 3 modes, global convergence to stationary points of ALS is not guaranteed in general [KB09] and requires some additional regularity conditions [Ber97,GS00]. Applying our BCD with diminishing radius, we derive the following ALS with diminishing radius for (24) in Algorithm 2.

3:
Update loading matrices θ n = [U (1) n , · · · ,U (m) n ] by 4: For i = 1, · · · , m do: (unfold(·, i ) denotes the mode-i tensor unfolding (see [KB09])) 5: end for 6: end for 7: output: θ N An immediate corollary of our Theorem 2.1 (i) is that Algorithm 2 always converges to the stationary points of the objective function for the constrained tensor factorization problem (24). Moreover, Theorem 2.1 (ii)-(iii) imply that the worst-case rate of convergence is of O((log n) 2 /n) with worst-case iteration complexity of O(ε −1 (log ε −1 ) 2 ). Furthermore, by Theorem 2.2, these results hold even when the sub-problems of nonnegative least squares are solved approximately as long as the optimality gaps are uniformly bounded. These results, which we summarize in Corollary 5.1, are new especially in the context of NCPD.
Corollary 5.1. Let θ n := [U (1) n , . . . ,U (m) n ], n ≥ 0 denote the output of Algorithm 2 with ∞ n=1 w n = ∞ and n=1 w 2 n < ∞. Then the following hold: (i) θ n converges to the stationary points of the objective function f CTF in (24) for every data tensor X and initial estimate θ 0 .
(ii) If w n = 1/( n log n), then the worst-case convergence rate up to iteration n is of order (log n) 2 /n in the sense of Theorem 2.1 (ii). Furthermore, the worst-case iteration complexity N ε for Algorithm 2 satisfies N ε = O(ε −1 (log ε −1 ) 2 ).
(iii) Suppose θ n for n ≥ 1 are computed approximately from Algorithm 2 such that the optimality gap is uniformly summable. Then (i)-(ii) still hold.
Remark 5.2. In a recent joint work with Strohmeier and Needell [SLN20], we have developed an algorithm for an online version of constrained tensor factorization and convergence analysis under mild assumptions. The analysis we provide in this work plays an important role in the reference.

EXPERIMENTAL VALIDATION
In this section, we compare the performance of our proposed BCD-DR algorithm (Algorithm 1, implemented as Algorithm 2 for NCPD) for the task of NCPD (24) (with no L 1 -regularization, i.e., λ 1 = λ 2 = λ 3 = 0) against two most popular approaches in practice: 1) The vanila BCD (2) (known as Alternating Least Squares for NCPD); and 2) Multiplicative Update (MU) (see [SH05]). Recall that we recover the standard ALS from our Algorithm 2 of ALS with diminishing radius by setting c w n = ∞ for n ≥ 1. We consider one synthetic and two real-world tensor data as follow: (1) X synth ∈ R 300×300×300 ≥0 equals Out(V 1 ,V 2 ,V 3 ), where the true loading matrices V 1 ∈ R 300×50
(2) X Twitter ∈ R 90×5000×1000 ≥0 (days×words×tweets) is the anonymized Twitter text data related to the COVID-19 pandemic from Feb. 1 to May 1 of 2020, which is recently obtained and analyzed in [KKL + 20]. The three modes correspond to days, words, and tweets, in order. In each day, the top 1000 most retweeted English tweets are collected and encoded as a 5000-dimensional tf-idf vector [RU11].
(3) X Headlines ∈ R 203×7000×700 ≥0 (months × words × headlines) is processed in [KKL + 20] from "A Million News Headlines" data, which consists of news headlines published over 203 months years sourced from the Australian news source ABC [Kul18]. The three modes correspond to months, words, and headlines. FIGURE 2. Comparison of the performance of BCD-DR (Algorithm 1) for the nonnegative CPdecomposition (NCPD) problem against BCD and Multiplicative Update (MU). BCD-DR is implmented as Algorithm 2 with c = X F /10 5 (X denoting the data tensor) and weights w n = n −β / log n for n ≥ 1. BCD is implmented as Algorithm 2 with c w n = ∞ for n ≥ 1 (equivalent to Alternating Least Sqaures). For each data tensor, we apply each algorithm ten times from the same random initialization to find nonnegative loading matrices U 1 ,U 2 ,U 3 with R = 5 columns. The average reconstruction error with 1 standard deviation are shown by the solid lines and shaded regions of respective colors.
For each data tensor X of shape I 1 × I 2 × I 3 , we used all three algorithms to learn three loading matrices U i ∈ R I i ×R for i = 1, 2, 3 with R = 5 columns. We plot the objective function in (24) against elapsed time (CPU time) in both cases in Figure 2. Each algorithm is used 10 times for the same data, and the plot shows the average relative reconstruction error (the square-root of f CTF in (24) divided by X F ) together with their standard deviation in shades.
In all experiments, BCD and BCD-DR with various hyperparameters show similar performance, all significantly outperforming MU. Interestingly, for the synthetic data (Figure 2 left), BCD-DR converges faster with stronger radius restriction and β = 5 (in blue) gives the best performance. BCD (in black) converges slower than all BCD-DR variants. Note that for β = 5, the weights w n are summable so the worst-case rate of convergence in Theorem 2.1 (i) is not meaningful. Hence this experimental result indicates that, in some cases, the average speed of convergence of BCD-DR could be significantly faster than in the worst-case. Clearly, if w n decays so fast that there is no stationary point within the 'total search range' ∞ n=1 w n , then no subsequence can converge to a stationary point (i.e., the left hand side of Theorem 2.1 (ii) cannot vanish). From this, we may conjecture that an optimal average rate of convergence of BCD-DR may be achieved when w n 's are summable but decays slow enough to be able to reach stationary points within the total search range. Further investigation on this matter is subject for a future study.


extfieldsofstudy:
Computer Science
Mathematics

text:
Neural Ensemble Search for Uncertainty Estimation and Dataset Shift

Ensembles of neural networks achieve superior performance compared to stand-alone networks in terms of accuracy, uncertainty calibration and robustness to dataset shift. \emph{Deep ensembles}, a state-of-the-art method for uncertainty estimation, only ensemble random initializations of a \emph{fixed} architecture. Instead, we propose two methods for automatically constructing ensembles with \emph{varying} architectures, which implicitly trade-off individual architectures' strengths against the ensemble's diversity and exploit architectural variation as a source of diversity. On a variety of classification tasks and modern architecture search spaces, we show that the resulting ensembles outperform deep ensembles not only in terms of accuracy but also uncertainty calibration and robustness to dataset shift. Our further analysis and ablation studies provide evidence of higher ensemble diversity due to architectural variation, resulting in ensembles that can outperform deep ensembles, even when having weaker average base learners. To foster reproducibility, our code is available: \url{https://github.com/automl/nes}

discussed in Section 2, while ensembles with varying architectures has already been explored in the literature, variation in architectures is typically limited to just varying depth and/or width, in contrast to more complex variations, such as changes in the topology of the connections and operations used, as considered in our work. More generally, automatic ensemble construction is well-explored in AutoML [17,37,46,44]. Our work builds on this by demonstrating that, in the context of uncertainty estimation, automatically constructed ensembles with varying architectures outperform deep ensembles that use state-of-the-art, or even optimal, architectures ( Figure 1). Studied under controlled settings, we assess the ensembles by various measures, including predictive performance, uncertainty estimation and calibration, base learner performance and two ensemble diversity metrics, showing that architectural variation is beneficial in ensembles. Ensemble diagnostics for all 6 metrics smaller (i.e. closer to center) is better Figure 1: A comparison of a deep ensemble with the best architecture (out of 15,625 possible architectures) and ensembles constructed by our method NES on ImageNet-16-120 over NAS-Bench-201. Performance metrics (smaller is better): error, negative loglikelihood (NLL) and expected calibration error (ECE). We also measure average base learner NLL and two metrics for ensemble diversity (see Section 3): oracle NLL and [1−predictive disagreement]; small means more diversity for both metrics. NES ensembles outperform the deep ensemble, despite the latter having a significantly stronger average base learner.
Note that, a priori, it is not obvious how to find a set of diverse architectures that work well as an ensemble. On the one hand, optimizing the base learners' architectures in isolation may yield multiple base learners with similar architectures (like a deep ensemble). On the other hand, selecting the architectures randomly may yield numerous base learners with poor architectures harming the ensemble. Moreover, as in neural architecture search (NAS), we face the challenge of needing to traverse vast architectural search spaces. We address these challenges in the problem of Neural Ensemble Search (NES), an extension of NAS that aims to find a set of complementary architectures that together form a strong ensemble. In summary, our contributions are as follows: 1. We present two NES algorithms for automatically constructing ensembles with varying base learner architectures. As a first step, we present NES with random search (NES-RS), which is simple and easily parallelizable. We further propose NES-RE inspired by regularized evolution [49], which evolves a population of architectures yielding performant and robust ensembles.   We refer to Appendix A for a proof. Proposition 3.1 suggests that it can be beneficial for ensembles to not only have strong average base learners (smaller upper bound), but also more diversity in their predictions (smaller lower bound). There is extensive theoretical work relating strong base learner performance and diversity with the generalization properties of ensembles [24,62,6,31,3,21]. In Section 5, the two metrics we use for measuring diversity are oracle ensemble loss and (normalized) predictive disagreement, defined as the average pairwise predictive disagreement amongst the base learners, normalized by their average error [18].

Visualizing Similarity in Base Learner Predictions
The fixed architecture used to build deep ensembles is typically chosen to be a strong stand-alone architecture, either hand-crafted or found by NAS. However, optimizing the base learner's architecture and then constructing a deep ensemble can neglect diversity in favor of strong base learner performance. Having base learner architectures vary allows more diversity in their predictions. We provide empirical evidence for this intuition by visualizing the base learners' predictions. Fort et al. [18] found that base learners in a deep ensemble explore different parts of the function space by means of applying dimensionality reduction to their predictions. Building on this, we uniformly sample five architectures from the DARTS search space [40], train 20 initializations of each architecture on CIFAR-10 and visualize the similarity among the networks' predictions on the test dataset using t-SNE [51]. Experiment details are available in Section 5 and Appendix B.
As shown in Figure 2a, we observe clustering of predictions made by different initializations of a fixed architecture, suggesting that base learners with varying architectures explore different parts of the function space. Moreover, we also visualize the predictions of base learners of two ensembles, each of size M = 30, where one is a deep ensemble and the other has varying architectures (found by NES-RS as presented in Section 4). Figure 2b shows more diversity in the ensemble with varying architectures than in the deep ensemble. These qualitative findings can be quantified by measuring diversity: for the two ensembles shown in Figure 2b, we find the predictive disagreement to be 94.6% for the ensemble constructed by NES and 76.7% for the deep ensemble (this is consistent across independent runs). This indicates higher predictive diversity in the ensemble with varying architectures, in line with the t-SNE results.

Neural Ensemble Search
In this section, we define neural ensemble search (NES). In summary, a NES algorithm optimizes the architectures of base learners in an ensemble to minimize ensemble loss.
Given a network f : R D → R C , let L(f, D) = (x,y)∈D (f (x), y) be the loss of f over dataset D.   Figure 3: Illustration of one iteration of NES-RE. Network architectures are represented as colored bars of different lengths illustrating different layers and widths. Starting with the current population, ensemble selection is applied to select parent candidates, among which one is sampled as the parent. A mutated copy of the parent is added to the population, and the oldest member is removed.
to explore fully. Second, a larger search space also increases the risk of overfitting the ensemble loss to D val . A possible approach here is to consider the ensemble as a single large network to which we apply NAS, but joint training of an ensemble through a single loss has been empirically observed to underperform training base learners independently, specially for large networks [52]. Instead, our general approach to solve Eq. 1 consists of two steps: 1. Pool building: build a pool P = {f θ1,α1 , . . . , f θ K ,α K } of size K consisting of potential base learners, where each f θi,αi is a network trained independently on D train . Step 1 reduces the options for the base learner architectures, with the intention to make the search more feasible and focus on strong architectures; step 2 then selects a performant ensemble. This procedure also ensures that the ensemble's base learners are trained independently. We use ensemble selection without replacement [7] for step 2. More specifically, this is forward step-wise selection; that is, given the set of networks P, we start with an empty ensemble and add the network from P which minimizes ensemble loss on D val . We repeat this without replacement until the ensemble is of size M . ForwardSelect(P, D val , M ) denotes the resulting set of M base learners selected from P.
Note that selecting the ensemble from P is a combinatorial optimization problem; a greedy approach such as ForwardSelect is nevertheless effective as shown by Caruana et al. [7], while keeping computational overhead low, given the predictions of the networks on D val . We also experimented with various other ensemble selection algorithms, including weighted averaging, as discussed in Section 5.2 and Appendix C.8, finding ForwardSelect to perform best.
We have not yet discussed the algorithm for building the pool in step 1; we present two approaches, NES-RS (Section 4.1) and NES-RE (Section 4.2). NES-RS is a simple random search based algorithm, while NES-RE is based on regularized evolution [49], a state-of-the-art NAS algorithm. Note that while gradient-based NAS methods have recently become popular, they are not naively applicable in our setting as the base learner selection component ForwardSelect is typically non-differentiable.

NES with Random Search
In NAS, random search (RS) is a competitive baseline on carefully designed architecture search spaces [38,58,59]. Motivated by its success and simplicity, NES with random search (NES-RS) builds the pool P by independently sampling architectures uniformly with replacement from the search space A (and training them). Since the architectures of networks in P vary, applying ensemble selection is a simple way to exploit diversity, yielding a performant ensemble. Importantly, NES-RS is easy to parallelize, exactly like deep ensembles. See Algorithm 2 in Appendix B.3 for pseudocode.

NES with Regularized Evolution
A more guided approach for building the pool P is using regularized evolution (RE) [49]. While RS has the benefit of simplicity by sampling architectures uniformly, the resulting pool might contain many weak architectures, leaving few strong architectures for ForwardSelect to choose between. RE is an evolutionary algorithm used on NAS spaces. It explores the search space by evolving (via mutations) a population of architectures. We first briefly describe RE as background before NES-RE. RE starts with a randomly initialized fixed-size population of architectures. At each iteration, a subset of size m of the population is sampled, from which the best network by validation loss is selected as the parent. A mutated copy (e.g. changing an operation in the network, see Appendix
// α stays in p. 6 Apply mutation to α, yielding child architecture β. 7 Train β using D train and add the trained network f θ,β to p and P. 8 Remove the oldest member in p. // as done in RE [49]. 9 Select base learners {f θ * for examples of mutations) of the parent architecture, called the child, is trained and added to the population, and the oldest member of the population is removed, preserving the population size. This is iterated until the computational budget is reached, returning the history, i.e. all the networks evaluated during the search, from which the best model is chosen by validation loss.
Building on RE for NAS, we propose NES-RE to build the pool of potential base learners. NES-RE starts by randomly initializing a population p of size P . At each iteration, we first apply ForwardSelect to the population to select an ensemble of size m, then we uniformly sample one base learner from this ensemble to be the parent. A mutated copy of the parent is added to p and the oldest network is removed, as in usual regularized evolution. This process is repeated until the computational budget is reached, and the history is returned as the pool P. See Algorithm 1 for pseudocode and Figure 3 for an illustration.
Also, note the distinction between the population and the pool in NES-RE: the population is evolved, whereas the pool is the set of all networks evaluated during evolution (i.e., the history) and is used post-hoc for selecting the ensemble. Moreover, ForwardSelect is used both for selecting m parent candidates (line 4 in NES-RE) and choosing the final ensemble of size M (line 9 in NES-RE). In general, m = M .

Ensemble Adaptation to Dataset Shift
Deep ensembles offer (some) robustness to dataset shift relative to training data. In general, one may not know the type of shift that occurs at test time. By using an ensemble, diversity in base learner predictions prevents the model from relying on one base learner's predictions which may not only be incorrect but also overconfident.
We assume that one does not have access to data points with test-time shift at training time, but one does have access to some validation data D shift val with a validation shift, which encapsulates one's belief about test-time shift. Crucially, test and validation shifts are disjoint. To adapt NES-RS and NES-RE to return ensembles robust to shift, we propose using D shift val instead of D val whenever applying ForwardSelect to select the final ensemble. In Algorithms 1 and 2, this is in lines 9 and 3, respectively. Our experiments show that this is highly effective against shift. Baselines. One of our main baselines is deep ensembles with fixed, optimized architectures. We consider various optimized architectures, indicated as "DeepEns (arch)": • On the DARTS search space, we consider the architectures found by: 1. the DARTS algorithm (DeepEns (DARTS)), 2. regularized evolution (DeepEns (AmoebaNet)), 3. random search, with the same number of networks trained as NES algorithms (DeepEns (RS)). • On the NAS-Bench-201 search space, in addition to DeepEns (RS), we consider: 1. the architecture found by GDAS 2 [12] (DeepEns (GDAS)), 2. the best architecture in the search space by validation loss (DeepEns (best arch.)). 3 We also compare to anchored ensembles [48], a recent technique for approximately Bayesian ensembles, which explicitly regularizes each base learner towards a fresh initialization sample and aims to induce more ensemble diversity. We use the DARTS architecture, and our implementation is described in Appendix B.5.
NES shows improved predictive uncertainty (NLL) and robustness to dataset shift ( Figure 4). Figure 4a shows the NLL achieved by NES-RS, NES-RE and the baselines as functions of the ensemble size M without dataset shift. We find that NES algorithms consistently outperform deep ensembles, with NES-RE usually outperforming NES-RS. Next, we evaluate the robustness of the ensembles to dataset shift in Figure 4b. In our setup, all base learners are trained on D train without data augmentation of shifted examples. However, as explained in Section 4.3, we use a shifted validation dataset, D shift val , and evaluate on a shifted test dataset, D shift test . The types of shifts appearing in D shift val are disjoint from those in D shift test . The severity of the shift varies between 1-5. We refer to Appendix B and Hendrycks & Dietterich [27] for details. The fixed architecture used in the baseline DeepEns (RS) is selected based on its loss over D shift val , but the DARTS and AmoebaNet are architectures from the literature. As shown in Figure 4b, ensembles picked by NES-RS and NES-RE are significantly more robust to dataset shift than the baselines, highlighting the effectiveness of applying ForwardSelect with D shift val . Unsurprisingly, AnchoredEns (DARTS) and DeepEns (DARTS/AmoebaNet) perform poorly compared to the other methods, as they are not optimized to deal with dataset shift here. The results of our experiments on Fashion-MNIST are in Appendix C.1.1. In line with the results in this section, both NES algorithms outperform deep ensembles with NES-RE performing best.
Better uncertainty calibration versus dataset shift and classification error ( Figure 5, Table 1). We also assess the ensembles by error and ECE. In short, ECE measures the mismatch between the model's confidence and accuracy. Figure 5 shows the ECE achieved by the ensembles at  Oracle ensemble NLL ImageNet-16-120 Figure 6: Ensemble, average base learner and oracle ensemble NLL versus budget K on ImageNet-16-120 over the NAS-Bench-201 search space. Ensemble size is M = 3.
varying dataset shift severities, noting that uncertainty calibration is especially important when models are used during dataset shift. Ensembles found with NES tend to exhibit better uncertainty calibration and are either competitive with or outperform anchored and deep ensembles for most shift severities. Notably, on CIFAR-10, ECE is reduced by up to 40% relative to the baselines. In terms of classification error, we find that ensembles constructed by NES outperform deep ensembles, with reductions of up to 7 percentage points in error, shown in Table 1. As with NLL, NES-RE tends to outperforms NES-RS. Ensembles found by NES tend to be more diverse ( Table 2). We measure ensemble diversity using the two metrics predictive disagreement (larger means more diverse) and oracle ensemble NLL (smaller means more diverse) as defined in Section 3 and average base learner NLL. Table  2 contrasts NES with the baselines in terms of these metrics. In terms of both diversity metrics, ensembles constructed by NES tend to be more diverse than anchored and deep ensembles. Ranking of the methods is largely consistent for different ensemble sizes M (see Appendix C). Unsurprisingly, note that the average base learner in NES is not always the best: by optimizing the fixed architecture first and then ensembling it, deep ensembles end up with a strong average base learner at the expense of less ensemble diversity. Despite higher average base learner NLL, NES ensembles perform better (recall Figure 4), highlighting once again the importance of diversity.
NES outperforms the deep ensemble of the best architecture in the NAS-Bench-201 search space ( Figure 6). Next, we compare NES to deep ensembles over the NAS-Bench-201 search space, which has two benefits: we demonstrate that our findings are not specific to the DARTS search space, and NAS-Bench-201 is an exhaustively evaluated search space for which all architectures' trained weights are available (three initializations per architecture), allowing us to compare NES to the deep ensemble of the best architecture by validation loss. Results shown in Figure 6 compare the losses of the ensemble, average base learner and oracle ensemble versus the number of networks evaluated K. Interestingly, although DeepEns (best arch.) has a significantly stronger average base learner than the other methods, its lack of diversity, as indicated by higher oracle ensemble loss ( Figure 6) and lower predictive disagreement (Figure 1), yields a weaker ensemble than both NES algorithms. Also, NES-RE outperforms NES-RS with a 6.6x speedup as shown in Figure 6-left.

Analysis and Ablations
Why does NES work? What if deep ensembles use ensemble selection over initializations? (Figure 7). NES algorithms differ from deep ensembles in two important ways: the ensembles use  varying architectures and NES utilizes ensemble selection (ForwardSelect) to pick the base learners. On Tiny ImageNet over the DARTS search space, we conduct an experiment to explore whether the improvement offered by NES over deep ensembles is only due to ensemble selection. The baselines "DeepEns + ES" operate as follows: we optimize a fixed architecture for the base learners, train K random initializations of it to form a pool and apply ForwardSelect to select an ensemble of size M . Figure 7-left shows that both NES algorithms (each shown for two computational budgets K = 200, 400) outperform all DeepEns + ES baselines. Figure 7-right contains the cost of each method in terms of the number of networks trained. Note that DeepEns + ES (RS) is the most competitive of the deep ensemble baselines, and, at an equal budget of 400, it is outperformed by both NES algorithms. However, even at half the cost (200), NES-RE outperforms DeepEns + ES (RS) while NES-RS performs competitively. As expected, deep ensembles with ensemble selection consistently perform better than without ensemble selection at the expense of higher computational cost, but do not close the gap with NES algorithms. We also re-emphasize that comparisons between NES and deep ensembles in our experiments always fix the base learner training routine and method for composing the ensemble, so variations in ensemble performance are only due to the architecture choices. Therefore, to summarize, combined with the finding that NES outperforms deep ensembles even when NES' average base learner is weaker (as in e.g. Figure 6 and Table 2), we find architectural variation in ensembles to be important for NES' performance gains. Computational cost and parallelizability. The primary computational cost involved in NES is training K networks to build the pool P. In our experiments on the DARTS search space, we set K to be 400 for CIFAR-10/100 and 200 for Tiny ImageNet (except Figure 7 which additionally considers K = 400). Figure 7-right gives an example of costs for each method. The cost of a deep ensemble with an optimized architecture stems from the initial architecture search and the subsequent training of M random initializations to build ensemble. We refer to Appendix B for details and discussion of computation cost, including training times. We note that apart from DeepEns (DARTS) and AnchoredEns (DARTS), NES has a lower computational cost than the baselines in our experiments. Similar to deep ensembles, training the pool for NES-RS is embarrassingly parallel. Our implementation of NES-RE is also parallelized as described in Appendix B.
Comparison of different ensemble selection algorithms. Both NES algorithms utilize ForwardSelect as the ensemble selection algorithm (ESA). We experimented with various other choices of ESAs, including weighted averaging and explicit diversity regularization, as shown in Figure 9. A detailed description of each ESA is in Appendix C.8. In summary, our choice of ForwardSelect performs better than or at par with all ESAs considered. Moreover, weighted averaging using stacking and Bayesian model averaging has   a very minor impact since the weights end up being close to uniform. Explicit diversity regularization, as shown in Figure 9b, appears to slightly improve performance in some cases, provided that the diversity regularization strength hyperparameter is appropriately tuned.
NES is insensitive to the size of validation dataset D val . We study the sensitivity of NES to the size of D val . Specifically, we measure test loss of the ensembles selected using D val of different sizes (with as few as 10 validation samples). The results in Figure 8 indicate that NES is insensitive to validation size for different levels of dataset shift severity, achieving almost the same loss when using 10× fewer validation data. We provide more details in Appendix C.5 and we also discuss why overfitting to D val is averted during ensemble selection.
Mutations of an architecture are similar in function space and by performance. To explore how NES-RE traverses the search space, we analyzed how mutations affect an architecture as follows. We sampled five random architectures from the DARTS space, and for each one, we applied a single random mutation twenty times, yielding a "family" of parent-children architectures, which are then trained on CIFAR-10. Figure 10-left shows the result of t-SNE applied to the test predictions of these architectures, where each color corresponds to one "family", of which the "star" is the parent architecture and the circles are the children architectures. The clustering demonstrates that architectures which differ by only a single mutation are similar in function space after training. Figure 10-right shows the NLL and error achieved by these architectures. Again, similar clustering shows that architectures differing by a single mutation also perform similarly w.r.t. NLL and error. This confirms that mutations allow for locally exploring the search space.
Further experiments. We also provide additional experiments in the Appendix. This includes a comparison of ensembles built by averaging logits vs. probabilities (Appendix C.7), a comparison of NES to ensembles with other hyperparameters being varied (either width/depth or training hyperparameters similar to Wenzel et al. [54]) (Appendix C.4) and using a weight-sharing model [2] as a proxy to accelerate the search phase in NES (Appendix C.9).

Conclusion, Limitations & Broader Impact
We presented Neural Ensemble Search for automatically constructing ensembles with varying architectures and demonstrated that the resulting ensembles outperform state-of-the-art deep ensembles in terms of uncertainty estimation and robustness to dataset shift. Our work highlights the benefit of ensembling varying architectures. In future work, we aim to address the limitation posed by the computational cost of NES due to building the pool P. An interesting approach in this direction could be the use of differentiable NAS methods to simultaneously optimize base learner architectures within a one-shot model and reduce cost [40,57,8]. More generally, we also hope to explore what other hyperparameters can be varied to improve ensemble performance and how best to select them.
Our work can readily be applied to many existing ensemble-based deep learning systems to improve their predictive performance. This work also focuses on improving uncertainty estimation in neural networks, which is a key problem of growing importance with implications for safe deployment of such systems. We are not aware of any direct negative societal impacts of our work, since NES is task-agnostic and its impact depends on its applications. Taking the loss function to be NLL, we have (f (x), y)) = − log [f (x)] y , where [f (x)] y is the probability assigned by the network f of x belonging to the true class y, i.e. indexing the predicted probabilities f (x) with the true target y. Note that t → − log t is a convex and decreasing function.

Acknowledgments and Disclosure of Funding
That is, f θ k assigns the highest probability to the correct class y for input x. Since − log is a decreasing function, we have We apply Jensen's inequality in its finite form for the second inequality. Jensen's inequality states that for a real-valued, convex function ϕ with its domain being a subset of R and numbers t 1 , . . . , t n in its domain, ϕ( 1

B Experimental and Implementation Details
We describe details of the experiments shown in Section 5 and Appendix C and include Algorithms 2 and 3 describing NES-RS and ForwardSelect, respectively. Note that unless stated otherwise, all sampling over a discrete set is done uniformly in the discussion below.
// without replacement. 5 return E // set of selected base learners.

B.1 Architecture Search Spaces
DARTS search space.
The first architecture search space we consider in our experiments is the one from DARTS [40]. We search for two types of cells: normal cells, which preserve the spatial dimensions, and reduction cells, which reduce the spatial dimensions. See Figure 11b for an illustration. The cells are stacked in a macro architecture, where they are repeated and connected using skip connections (shown in Figure 11a). Each cell is a directed acyclic graph, where nodes represent feature maps in the computational graph and edges between them correspond to operation choices (e.g. a convolution operation). The cell parses inputs from the previous and previous-previous cells in its 2    information coming from 2 previous nodes in the cell and finally an output node that concatenates the output of all intermediate nodes across the channel dimension. AmoebaNet contains one more intermediate node, making that a deeper architecture. The set of possible operations (eight in total in DARTS) that we use for each edge in the cells is the same as DARTS, but we leave out the "zero" operation since that is not necessary for non-differentiable approaches such as random search and evolution. Specifically, this leaves us with the following set of seven operations: 3 × 3 and 5 × 5 separable convolutions, 3 × 3 and 5 × 5 dilated separable convolutions, 3 × 3 average pooling, 3 × 3 max pooling and identity. Randomly of architectures is done by sampling the structure of the cell and the operations at each edge. The total number of architectures contained in this space is ≈ 10 18 . We refer the reader to Liu et al. [40] for more details.
NAS-Bench-201 search space. NAS-Bench-201 [13] is a tabular NAS benchmark, i.e. all architectures in the cell search space are trained and evaluated beforehand so one can query their performance (and weights) from a table quickly. Since this space is exhaustively evaluated, its size is also limited to only normal cells containing 4 nodes in total (1 input, 2 intermediate and 1 output node) and 5 operation choices on every edge connecting two nodes. This means that there are only 15,625 possible architecture configurations in this space. The networks are constructed by stacking 5 cells with in-between fixed residual blocks for reducing the spacial resolution. Each of them is trained for 200 epochs 3 times with 3 different seeds on 3 image classification datasets. For more details, please refer to Dong & Yang [13].

B.2 Datasets
Fashion-MNIST [56]. Fashion-MNIST consists of a training set of 60k 28×28 grayscale images and a test set of 10k images. The number of total labels is 10 classes. We split the 60k training set images to 50k used to train the networks and 10k used only for validation.
CIFAR-10/100 [32]. CIFAR-10 and CIFAR-100 both consist of 60k 32×32 colour images with 10 and 100 classes, respectively. We use 10k of the 60k training images as the validation set. We use the 10k original test set for final evaluation.
Tiny ImageNet [35]. Tiny Imagenet has 200 classes and each class has 500 training, 50 validation and 50 test colour images with 64×64 resolution. Since the original test labels are not available, we split the 10k validation examples into 5k for testing and 5k for validation.
Note that the test data points are only used for final evaluation. The data points for validation are used by the NES algorithms and DeepEns + ES baselines during ensemble selection and by DeepEns (RS) for picking the best architecture from the pool to use in the deep ensemble. Note that when considering dataset shift for CIFAR-10, CIFAR-100 and Tiny ImageNet, we also apply two disjoint sets of "corruptions" (following the terminology used by Hendrycks & Dietterich [27]) to the validation and test sets. We never apply any corruption to the training data. More specifically, out of the 19 different corruptions provided by Hendrycks & Dietterich [27], we randomly apply one from {Speckle Noise, Gaussian Blur, Spatter, Saturate} to each data point in the validation set and one from {Gaussian Noise, Shot Noise, Impulse Noise, Defocus Blur, Glass Blur, Motion Blur, Zoom Blur, Snow, Frost, Fog, Brightness, Contrast, Elastic Transform, Pixelate, JPEG compression} to each data point in the test set. This choice of validation and test corruptions follows the recommendation of Hendrycks & Dietterich [27]. Also, as mentioned in Section 5, each of these corruptions has 5 severity levels, which yields 5 corresponding severity levels for D shift val and D shift test .

B.3 Training Routine & Time
The macro-architecture we use has 16 initial channels and 8 cells (6 normal and 2 reduction) and was trained using a batch size of 100 for 100 epochs for CIFAR-10/100 and 15 epochs for Fashion-MNIST. For Tiny ImageNet, we used 36 initial channels and a batch size of 128 for 100 epochs. Training a single network took roughly 40 minutes and 3 GPU hours 4 for CIFAR-10/100 and Tiny ImageNet, respectively. The networks are optimized using SGD with momentum set to 0.9. We used a learning rate of 0.025 for CIFAR-10/100 and Fashion-MNIST, and 0.1 for Tiny ImageNet. Unlike DARTS, we do not use any data augmentation procedure during training, nor any additional regularization such as ScheduledDropPath [65] or auxiliary heads, except for the case of Tiny ImageNet, for which we used ScheduledDropPath, gradient clipping and standard data augmentation as default in DARTS. All other hyperparameter settings are exactly as in DARTS [40].
All results containing error bars are averaged over multiple runs (at least five) with error bars indicating a 95% confidence interval. We used a budget K = 400 for CIFAR-10/100 (corresponding to 267 GPU hours). For Tiny ImageNet on the DARTS search space, unless otherwise stated, we used K = 200 (corresponding to 600 GPU hours). Note that only in Figure 7 and Table 30 we also tried NES algorithms using a higher budget of K = 400 on Tiny ImageNet for an equal-cost comparison with DeepEns + ES (RS). For ImageNet-16-120 on the NAS-Bench-201 search space, we used K = 1000 (compute costs are negligible since it is a tabular benchmark). Note that NES algorithms can be parallelized, especially NES-RS which is embarrassingly parallel like deep ensembles, therefore training time can be reduced easily when using multiple GPUs.

B.4 Implementation Details of NES-RE
Parallization. Running NES-RE on a single GPU requires evaluating hundreds of networks sequentially, which is tedious. To circumvent this, we distribute the "while |P| < K" loop in Algorithm 1 over multiple GPUs, called worker nodes. We use the parallelism scheme provided by the hpbandster [16] codebase. 5 In brief, the master node keeps track of the population and history (lines 1, 4-6, 8 in Algorithm 1), and it distributes the training of the networks to the individual worker nodes (lines 2, 7 in Algorithm 1). In our experiments, we always use 20 worker nodes and evolve a population p of size P = 50 when working over the DARTS search space. Over NAS-Bench-201, we used one worker since it is a tabular NAS benchmark and hence is quick to evaluate on. During iterations of evolution, we use an ensemble size of m = 10 to select parent candidates.
Mutations. We adapt the mutations used in RE to the DARTS search space. As in RE, we first pick a normal or reduction cell at random to mutate and then sample one of the following mutations: • identity: no mutation is applied to the cell.
• op mutation: sample one edge in the cell and replace its operation with another operation sampled from the list of operations described in Appendix B.1.
• hidden state mutation: sample one intermediate node in the cell, then sample one of its two incoming edges. Replace the input node of that edge with another sampled node, without altering the edge's operation.
For example, one possible mutation would be changing the operation on an edge in the cell (as shown in e.g. Figure 16) from say 5 × 5 separable convolution to 3 × 3 max pooling. See Real et al. [49] for further details and illustrations of these mutations. Note that for NAS-Bench-201, following Dong & Yang [13] we only use op mutation.
Adaptation of NES-RE to dataset shifts. As described in Section 4.3, at each iteration of evolution, the validation set used in line 4 of Algorithm 1 is sampled uniformly between D val and D shift val when dealing with dataset shift. In this case, we use shift severity level 5 for D shift val . Once the evolution is complete and the pool P has been formed, then for each severity level s ∈ {0, 1, . . . , 5}, we apply ForwardSelect with D shift val of severity s to select an ensemble from P (line 9 in Algorithm 1), which is then evaluated on D shift test of severity s. (Here s = 0 corresponds to no shift.) This only applies to CIFAR-10, CIFAR-100 and Tiny ImageNet, as we do not consider dataset shift for Fashion-MNIST and ImageNet-16-120 .

B.5 Implementation details of anchored ensembles.
Anchored ensembles [48] are constructed by independently training M base learners, each regularized towards a new initialization sample (called the anchor point). Specifically, the k-th base learner f θ k is trained to minimize the following loss (using the notation in Section 3.1): is the cross-entropy loss for the i-th datapoint. Moreover, the anchor point • θ k is an independently sampled initialization, and λ is the regularization strength. Defining p to be the number of parameters, i.e. θ ∈ R p , and letting σ 2 i be the variance of the initialization for the i-th parameter (θ k ) i , Γ ∈ R p×p is the diagonal matrix with Γ ii = 1/(2σ 2 i ). Although Pearce et al. [48] do not include a tune-able regularization parameter λ (they set λ = 1), we found anchored ensembles performed poorly without tuning λ. For CIFAR-10/100, we used λ = 0.4 and for Tiny ImageNet, we used λ = 0.1 (tuned by grid search). We followed the recommendation of Pearce et al. [48] in using a different initialization for optimization than the anchor point • θ k . Moreover, we turned off weight decay when training anchored ensembles, and we did not apply the regularization to the trainable parameters in batch normalization layers (which are initialized deterministically).

C Additional Experiments
In this section we provide additional results for the experiments conducted in Section 5. Note that, as with all results shown in Section 5, all evaluations are made on test data unless stated otherwise. Notice that, apart from the case when ensemble size M = 30, NES-RS and NES-RE find ensembles with both stronger and more diverse base learners (smaller losses in Figures 13 and 14, respectively). While it is expected that the oracle ensemble loss is smaller for NES-RS and NES-RE compared to DeepEns (RS), it initially appears surprising that DeepEns (RS) has a larger average base learner loss considering that the architecture for the deep ensemble is chosen to minimize the base learner loss. We found that this is due to the loss having a sensitive dependence not only on the architecture but also the initialization of the base learner networks. Therefore, re-training the best architecture by validation loss to build the deep ensemble yields base learners with higher losses due to the use of different random initializations. Fortunately, NES algorithms are not affected by this, since they simply select the ensemble's base learners from the pool without having to re-train anything which allows them to exploit good architectures as well as initializations. Note that, for CIFAR-10-C experiments, this was not the case; base learner losses did not have as sensitive a dependence on the initialization as they did on the architecture.
In Table 3, we compare the classification error and expected calibration error (ECE) of NES algorithms with the deep ensembles baseline for various ensemble sizes on Fashion-MNIST. Similar to the loss, NES algorithms also achieve smaller errors, while ECE remains approximately the same for all methods.

C.1.2 Entropy on out-of-distribution inputs
To assess how well models respond to completely out-of-distribution (OOD) inputs (inputs which do not belong to one of the classes the model can predict), we investigate the entropy of the predicted probability distribution over the classes when the input is OOD. Higher entropy of the predicted probabilities indicates more uncertainty in the model's output. For CIFAR-10 on the DARTS search space, we compare the entropy of the predictions made by NES ensembles with deep ensembles on two types of OOD inputs: images from the SVHN dataset and Gaussian noise. In Figure 15, we notice that NES ensembles indicate higher uncertainty when given inputs of Gaussian noise than deep ensembles but behave similarly to deep ensembles for inputs from SVHN.

C.1.3 Additional Results on the DARTS search space
In this section, we provide additional experimental results on CIFAR-10, CIFAR-100 and Tiny ImageNet on the DARTS search space, complementing the results in Section 5 as shown in Figures 17- Figure 16 for Tiny ImageNet, showcasing variation in architectures.  Results on CIFAR for larger models. In additional to the results on CIFAR-10 and CIFAR-100 on the DARTS search space using the settings described in Appendix B.3, we also train larger models (around 3M parameters) by scaling up the number of stacks cells and initial channels in the network. We run NES and other baselines similarly as done before and plot results in Figure 27 and 28 for NLL and classification test error with budget K = 90. As shown, NES algorithms tend to outperform or be competitive with the baselines. Note, more runs are needed including error bars for conclusive results in this case.  Ensemble size (M)    Ensemble size (M)

C.3 What if deep ensembles use ensemble selection over initializations?
We provide additional experimental results in this section for comparing NES to deep ensembles with ensemble selection building on those shown in Figure 7 in Section 5. The table in Figure  30-right compares the methods with respect to classification error as well as loss at different dataset shift severities. In the absence of any dataset shift, we find that NES-RE outperforms all baselines with respect to both metrics, specially at the higher budget of K = 400 and NES-RS performs competitively. At shift severity 5, NES algorithms also outperform the baselines, with NES-RS performing slightly better than NES-RE. The plot in Figure 30-left shows similar results on CIFAR-10 over the DARTS search space for the same experiment as the one conducted on Tiny ImageNet and shown in Figure 7.

C.4 Comparing NES to ensembles with other varying hyperparameters
Since varying the architecture in an ensemble improves predictive performance and uncertainty estimation as demonstrated in Section 5, it is natural to ask what other hyperparameters should be 5 10 15  Figure 31: Plots show NLL vs. ensemble sizes comparing NES to the baselines introduced in Appendix C.4 on CIFAR-10 and CIFAR-100 with and without respective dataset shifts [27].
varied in an ensemble. It is also unclear which hyperparameters might be more important than others. Note that concurrent work by Wenzel et al. [54] has shown that varying hyperparameters such as L 2 regularization strength, dropout rate and label smoothing parameter also improves upon deep ensembles. While these questions lie outside the scope of our work and are left for future work, we conduct preliminary experiments to address them.   In this section, we consider two additional baselines working over the DARTS search space on CIFAR-10/100: The results shown in Figures 31 and Table 32b compare the two baselines above to DeepEns (DARTS), NES-RS and NES-RE. 6 As shown in Figure 31, NES-RE tends to outperform the baselines, though is at par with HyperEns on CIFAR-100 without dataset shift (Figure 31a). Under the presence of dataset shift ( Figures  31b and 31c), both NES algorithms substantially outperform all baselines. Note that both HyperEns and NES-RS (depth, width) follow the same protocol as NES-RS and NES-RE: ensemble selection uses a shifted validation dataset when evaluating on a shifted test dataset. In terms of classification error, the observations are similar as shown in Table  32b. Lastly, we view the diversity of the ensembles from the perspective of oracle ensemble loss in Figure 32a. As in Section 5, results here also suggest that NES agorithms tend to find more diverse ensembles despite having higher average base learner loss.  In this section, we consider how sensitive the performance of NES is to changes in the size of the validation dataset D val , and we discuss why overfitting to D val is not a concern in our experiments. Recall that D val is used by NES algorithms during ensemble selection from the pool of base learners as outlined in Algorithm 3, and we use D val or D shift val depending on whether dataset shift is expected at test time as described in Section 4.3. Over the DARTS search space for CIFAR-100 and Tiny ImageNet, we measure the test loss of NES as a function of different validation dataset sizes. Specifically, we measure test loss of the ensembles selected using validation datasets of different sizes (with as few as 10 validation samples). Figure 33b shows this for different levels of dataset shift severity. The results indicate that NES is insensitive to the validation set size, achieving almost the same loss when using 10× fewer validation data.
NES also did not overfit to the validation data in our experiments. This can be seen in Figure 33a for CIFAR-100 over the DARTS search space, which shows that both the test and validation losses decrease over time (i.e. as the pool size K increases). This finding is consistent across datasets and search spaces in our experiments. We also remark that overfitting is unexpected, because ensemble selection "fits a small number of parameters" since it only selects which base learners are added to the ensemble.  In practice, one might choose to re-trained the selected architectures from NES on the D train + D val to make maximal use of the data available. We re-trained the ensembles constructed by NES and the best performing deep ensemble baseline, DeepEns (DARTS), on D train + D val , finding that both ensembles improve since the base learners improve due to more training data (see Figure 34). Note that, as shown in Appendix C.5, the performance of ensemble selection (which uses D val ) is relatively insensitive to the size of D val . Therefore, one way to bypass the additional cost of having to re-train the ensembles on D train + D val is to simply pick a very small D val , such that the performance of the models when trained on D train + D val is approximately the same as when trained on D train .

C.7 Averaging logits vs. averaging probabilities in an ensemble
In Figure 35, we explore whether ensembles should be constructed by averaging probabilities (i.e. post-softmax) as done in our work or by averaging the logits (i.e. pre-softmax). We find that while classification error of the resulting ensembles is similar, ensembles with averaged probabilities perform notably better in terms of uncertainty estimation (NLL) and are better calibrated as well (ECE).

C.8 Comparison of ensemble selection algorithms
In Figure 9a, we compare the performance of NES for different choices of the ensemble selection algorithm (ESA) used in stage 2 (ensemble selection from pool). In addition to our default ForwardSelect without replacement, we experiment with the following seven choices: • TopM: selects the top M (ensemble size) models by validation performance.
• QuickAndGreedy: Starting with the best network by validation performance, add the next best network to the ensemble only if it improves validation performance, iterating until the ensemble size is M or all models have been considered (returning an ensemble of size at most M). • ForwardSelect with replacement: ForwardSelect but with replacement which allows repetitions of base learners. • Stacking (Weighted): Linearly combines all base learners in the pool with learned stacking weights (which are positive and sum to 1). Then keeps only the base learners with the M largest stacking weights and weights them using the renormalized learned stacking weights. • Stacking (Unweighted): Linearly combines all base learners in the pool with learned stacking weights (which are positive and sum to 1). Then keeps only the base learners with the M largest stacking weights and combines them by a simple, unweighted average. • ForwardSelect with Bayesian model averaging by likelihood: Selects the base learners in the ensemble using ForwardSelect and then takes an average weighted by the (normalized) validation likelihoods of the base learners. • ForwardSelect with Bayesian model averaging by accuracy: Selects the base learners in the ensemble using ForwardSelect and then takes an average weighted by the (normalized) validation accuracies of the base learners.
In Figure 36 we provide additional results including the one in Figure 9a. From the plots we can observe that: • ForwardSelect performs better than or at par with all ESAs considered here.
• Stacking tends to perform competitively but still worse than ForwardSelect, and whether weighted averaging is used or not has a very minor impact on performance, since the learned weights tend to be close to uniform.
• Bayesian model averaging (BMA) also has a very minor impact (almost invisible in the plots in Figure 36b), because the likelihoods and accuracies of individual base learners are very similar (a consequence of multi-modal loss landscapes in neural networks with different models achieving similar losses) hence the BMA weights are close to uniform. The NLL achieved by the BMA ensemble only differed at the 4th or 5th decimal places compared to unweighted ensembles with the same base learners. Lastly, we assess the impact of explicitly regularizing for diversity during ensemble selection as follows: we use ForwardSelect as before but instead use it to minimize the objective "validation loss -diversity strength × diversity" where diversity is defined as the average (across base learners and data points) L 2 distance between a base learner's predicted class probabilities and the ensemble's predicted class probabilities. The plots below show the results for different choices of the "diversity strength" hyperparameter. In summary, if appropriately tuned, ForwardSelect with diversity performs slightly better than usual ForwardSelect for NES-RS, though for NES-RE, the diversity term seems to harm performance as shown in Figure 37-right.
C.9 Generating the architecture pool using weight-sharing NAS algorithms In order to accelerate the NES search phase, we generated the pool using the weight sharing schemes proposed by Random Search with Weight Sharing [38] and DARTS [40]. Specifically, we trained one-shot weight-sharing models using each of these two algorithms, then we sampled architectures from the weightshared models uniformly at random to build the pool. Next, we ran ForwardSelect as usual to select ensembles from these pools. Finally, we retrained the selected base learners from scratch using the same training pipeline as NES and the deep ensemble baselines. In Figure 38, we refer to the results of these methods as NES-RS-WS and NES-RS-DARTS, respectively. As we can see, their performance is worse than NES-RS and NES-RE, even though the computation cost is reduced significantly. This is likely due to the low correlation between the performance of the architectures when evaluated using the shared weights and the performance when re-trained in isolation. Prior work has shown this is caused by weight interference and co-adaptation which occurs during the weight-sharing model training [59,60]. Nonetheless, we believe this is a promising avenue for future research requiring more development, as it substantially reduces the cost of exploring the search space.