Projections onto the canonical simplex with additional linear inequalities

TL;DR

This paper develops efficient algorithms for projecting onto the canonical simplex with extra linear constraints, crucial for distributionally robust optimization, by transforming the problem into finding zeros of well-behaved functions.

Contribution

It introduces a novel approach to compute projections with additional linear inequalities using zero-finding of convex or monotonic functions, with guaranteed convergence and near-linear complexity.

Findings

Algorithms have guaranteed convergence.

Observed complexity is nearly linear.

Applicable to various distance functions in distributionally robust optimization.

Abstract

We consider the distributionally robust optimization and show that computing the distributional worst-case is equivalent to computing the projection onto the canonical simplex with additional linear inequality. We consider several distance functions to measure the distance of distributions. We write the projections as optimization problems and show that they are equivalent to finding a zero of real-valued functions. We prove that these functions possess nice properties such as monotonicity or convexity. We design optimization methods with guaranteed convergence and derive their theoretical complexity. We demonstrate that our methods have (almost) linear observed complexity.

Tables

Table 1. Table 1: Examples of the ϕ italic-ϕ \phi -divergences d 𝑑 d and generating functions ϕ italic-ϕ \phi . Note that ϕ italic-ϕ \phi is convex with ϕ ​ ( 1 ) = 0 italic-ϕ 1 0 \phi(1)=0 .

Name	Generating function $\varphi$	Formula $d$
Kullback-Leibler divergence	${\varphi }_1\left(t\right)=t\log t$	$d_1(p,q)=p\log \left(\frac{p}{q}\right)$
Burg entropy	${\varphi }_2\left(t\right)=-\log t$	$d_2(p,q)=q\log \left(\frac{q}{p}\right)$
Hellinger distance	${\varphi }_3\left(t\right)={\left(\sqrt{t}-1\right)}^2$	$d_3(p,q)={\left(\sqrt{p}-\sqrt{q}\right)}^2$
${\chi }^2$-distance	${\varphi }_4\left(t\right)=\frac{1}{t}{\left(t-1\right)}^2$	$d_4(p,q)=\frac{{\left(p-q\right)}^2}{p}$
Modified ${\chi }^2$-distance	${\varphi }_5\left(t\right)={\left(t-1\right)}^2$	$d_5(p,q)=\frac{{\left(p-q\right)}^2}{q}$

Table 2. Table 2: Table showing which problems have an exact algorithm and for which problems, the bisection and Newton’s methods are convergent. Note that all problems have at least one convergent algorithm.

		Guaranteed convergence
	Exact algorithm	Bisection	Newton
(DRO1) with Kullback-Leibler divergence	✗	✓	✗
(DRO1) with Burg entropy	✗	✓	✗
(DRO1) with Hellinger distance	✗	✓	✗
(DRO1) with ${\chi }^2$-distance	✗	✓	✓
(DRO1) with Modified ${\chi }^2$-distance	✗	✓	✓
(DRO2) with $l_1$ norm	✓	$\cdot$	$\cdot$
(DRO2) with $l_2$ norm	✗	✓	✓
(DRO2) with $l_{\infty }$ norm	✓	$\cdot$	$\cdot$
(SIMPLEX)	✓	$\cdot$	$\cdot$

Table 3. Table 3: Properties for the bisection and Newton’s method. The table shows which equation needs to be solved for ( 7 ), the bounds and whether the function h ℎ h is convex.

	Equation	Bounds	Convex
(DRO1) with Kullback-Leibler divergence	$h_1\left(\mu \right)=0$		✗
(DRO1) with Burg entropy	$h_2\left(\lambda \right)=0$		✗
(DRO1) with Hellinger distance	$h_3\left(\lambda \right)=0$		✗
(DRO1) with ${\chi }^2$-distance	$h_4\left(\lambda \right)=0$		✓
(DRO1) with Modified ${\chi }^2$-distance	$h_5\left(\lambda \right)=0$		✓
(DRO2) with $l_2$ norm	$h_6\left(\mu \right)=0$		✓

Table 4. Table 4: Computational complexity evaluating function h ℎ h and the total complexity for solving problems ( DRO1 ), ( DRO2 ) and ( SIMPLEX ). Here, n h subscript 𝑛 ℎ {n_{h}} refers to the number of evaluation of the function h ℎ h which is guaranteed to be constant for the bisection method whenever its bracketing interval from Table 3 is uniformly bounded. The observed complexity is shown in Table 5 later.

	Evaluation of $h$	Total
(DRO1)	$O\left(n\right)$	$O\left(n_{h}n\right)$
(DRO2) with $l_1$ norm	$-$	$O\left(n\log n\right)$
(DRO2) with $l_2$ norm	$O\left(n\right)$	$O\left(n_{h}n\right)$
(DRO2) with $l_{\infty }$ norm	$-$	$O\left(n\log n\right)$
(SIMPLEX)	$-$	$O\left(n\log n\right)$

Table 5. Table 5: Comparison of the observed and theoretical complexity for our methods. In most cases our methods exhibit the complexity of O ​ ( n ) 𝑂 𝑛 O(n) or O ​ ( n ​ log ⁡ n ) 𝑂 𝑛 𝑛 O(n\log n) which concurs with Table 4 .

	Observed complexity	Theoretical complexity
(DRO1) with Kullback-Leibler divergence	$3.328\cdot {10}^{-7}{n}^{1.102}$	$O\left(n{n}_h\right)$
(DRO2) with $l_1$ norm	$2.794\cdot {10}^{-8}{n}^{1.125}$	$O\left(n\log n\right)$
(DRO2) with $l_2$ norm	$3.759\cdot {10}^{-7}{n}^{1.056}$	$O\left(n{n}_h\right)$
(DRO2) with $l_{\infty }$ norm	$2.802\cdot {10}^{-7}{n}^{1.042}$	$O\left(n\log n\right)$
(SIMPLEX)	$4.911\cdot {10}^{-7}{n}^{1.039}$	$O\left(n\log n\right)$

Equations

$$\underset{\mathbf{p}}{\mathrm{minimize}}$$

$$subject\; to$$

$$0\le p_i,\forall i=1,\dots ,n.$$

$$\sum _{i=1}^n\max \{q_i-\mu ,0\}=1$$

$$p_i=\max \{q_i-\mu ,0\}.$$

$$\underset{\mathbf{x}}{\mathrm{minimize}}\underset{\mathbf{\xi }\in \Xi }{\mathrm{maximize}}f(\mathbf{x},\mathbf{\xi }).$$

$$\underset{\mathbf{x}}{\mathrm{minimize}}\underset{P\in \mathcal{P}}{\mathrm{maximize}}{\mathbb{E}}_{P}f(\mathbf{x},\mathbf{\xi }).$$

$$\underset{\mathbf{p}\in \mathcal{P}}{\mathrm{maximize}}{\mathbb{E}}_{P}f(\mathbf{x},\mathbf{\xi })=\underset{\mathbf{p}\in \mathcal{P}}{\mathrm{maximize}}\sum _{i=1}^n{p}_{i}f(\mathbf{x},{\mathbf{\xi }}_i)=\underset{\mathbf{p}\in \mathcal{P}}{\mathrm{maximize}}{\mathbf{c}}^{\top }\mathbf{p},$$

$$\underset{\mathbf{p}}{\mathrm{maximize}}$$

$$subject\; to$$

$$0\le p_i,\forall i=1,\dots ,n,$$

$$\hat{d}(\mathbf{p},\mathbf{q})\le \epsilon ,$$

$$\underset{\mathbf{p}}{\mathrm{minimize}}$$

$$subject\; to$$

$$0\le p_i,\forall i=1,\dots ,n,$$

$${\mathbf{c}}^{\top }\mathbf{p}\ge \delta .$$

$$c_{\max }$$

$$c_{\min }$$

$$I$$

$$\hat{d}(\bm{p},\bm{q})=\sum_{i=1}^{n}d(p_{i},q_{i})=\sum_{i=1}^{n}q_{i}\phi\mathopen{}\mathclose{{}\left(\frac{p_{i}}{q_{i}}}\right),$$

$$\underset{\mathbf{p}}{\mathrm{maximize}}$$

$$subject\; to$$

$$0\le p_i,\forall i=1,\dots ,n,$$

$$\displaystyle\sum_{i=1}^{n}q_{i}\phi\mathopen{}\mathclose{{}\left(\frac{p_{i}}{q_{i}}}\right)\leq{\varepsilon},$$

$${\hat{p}}_i=\{\begin{array}{ll}\frac{1}{{\sum }_{j\in I}{q}_j}{q}_i& \text{if }i\in I,\\ 0& \text{otherwise.}\end{array}$$

$$\sum_{i=1}^{n}q_{i}\phi\mathopen{}\mathclose{{}\left(\frac{\hat{p}_{i}}{q_{i}}}\right)\leq{\varepsilon}$$

$$h_{1}(\mu):=\sum_{i=1}^{n}q_{i}\exp\mathopen{}\mathclose{{}\left(\frac{c_{i}}{\mu}}\right)\mathopen{}\mathclose{{}\left(\frac{c_{i}}{\mu}-\log\mathopen{}\mathclose{{}\left(\sum_{j=1}^{n}q_{j}\exp\mathopen{}\mathclose{{}\left(\frac{c_{j}}{\mu}}\right)}\right)-{\varepsilon}}\right)=0.$$

$$\hat{p}_{i}=q_{i}\exp\mathopen{}\mathclose{{}\left(\frac{c_{i}}{\mu}}\right),$$

$$h_{2}(\lambda):=\sum_{i=1}^{n}q_{i}\log\mathopen{}\mathclose{{}\left(\lambda-c_{i}}\right)+\log\mathopen{}\mathclose{{}\left(\sum_{i=1}^{n}\frac{q_{i}}{\lambda-c_{i}}}\right)-{\varepsilon}=0.$$

$${\hat{p}}_i=q_i\frac{1}{\lambda -c_i},$$

$$h_3(\lambda ):=2\sum _{i=1}^n\frac{q_i}{\lambda -c_i}-(2-\epsilon )\sqrt{\sum _{i=1}^n\frac{q_i}{(\lambda -c_i{)}^2}}=0.$$

$${\hat{p}}_i=q_i\frac{1}{(\lambda -c_i{)}^2},$$

$$h_{4}(\lambda):=\mathopen{}\mathclose{{}\left(\sum_{i=1}^{n}q_{i}\sqrt{\lambda-c_{i}}}\right)\mathopen{}\mathclose{{}\left(\sum_{i=1}^{n}q_{i}\frac{1}{\sqrt{\lambda-c_{i}}}}\right)-1-{\varepsilon}=0.$$

$${\hat{p}}_i=q_i\frac{1}{\sqrt{\lambda -c_i}},$$

$$h_{5}(\lambda):=\sum_{i=1}^{n}q_{i}{\max}^{2}\mathopen{}\mathclose{{}\left(c_{i}+\lambda,0}\right)-(1+{\varepsilon})\mathopen{}\mathclose{{}\left(\sum_{i=1}^{n}q_{i}\max\mathopen{}\mathclose{{}\left(c_{i}+\lambda,0}\right)}\right)^{2}=0.$$

$$\hat{p}_{i}=q_{i}\max\mathopen{}\mathclose{{}\left(c_{i}+\lambda,0}\right),$$

$$\underset{\mathbf{p}}{\mathrm{maximize}}$$

Full text

Keywords: projection; simplex; distributionally robust optimization; linear obseved complexity.

1 Introduction

The projection of a vector onto the unit simplex appears in various fields such as portfolio optimization [18], multi-phase physics [4], mathematical optimization [20], knapsack problem [12] or machine learning applications [5]. Given a vector $\bm{q}\in{\mathbb{R}}^{n}$, this projection amounts to solving

[TABLE]

The Karush-Kuhn-Tucker optimality conditions imply that if one solves

[TABLE]

for $\mu$, then one recovers the optimal solution of (1) by thresholding

[TABLE]

This was discovered for the first time in [10] and then rediscovered many times later [17]. The simplest way to solve (2) is to sort $\bm{q}$, derive an iterative procedure computing the whole function on the left-hand side of (2) and then find when its value equals to $1$. Since the second part can be done in $O(n)$, the whole algorithm has complexity $O(n\log n)$ due to the sorting.

This procedure was improved in numerous papers. [24] observed that only those $q_{i}$ above $\mu$ need to be sorted in (2). Using a partially sorted structure called heap, they managed to reduce the complexity to $O(n+k\log n)$, where $k$ is the number of $q_{i}$ above the optimal $\mu$. [13] realized that many operations in quicksort may be ignored when it is used to solve (2) and reached complexity $O(n)$. [19] proposed a simple method based on the fixed-point theorem with observed complexity $O(n)$. [6] provided an excellent overview, pointed to some errors in previous papers and designed an improved algorithm.

Besides the standard projection (1), multiple versions appear in the literature. [2] considered an infinite-dimensional optimization problem with partial differential equations in the constraints. To get the independence of the number of iterations on the mesh size, they derived a path-following algorithm. [16] considered sparse learning problems containing a modified simplex with two vectors of variables whose sum had to be equal. They derived an improved bisection method and a fast-converging subgradient algorithm. A similar simplex appeared in [23] for ranking labels based on feedback and in [15] for a special binary classification problem. [14] considered SVM with top-${\rm k}$ error instead of the standard top-$1$ error. The resulting modified simplex contained a variable upper bound which was not fixed as in all previous cases. They penalized one constraint and computed an approximate projection. Note that these problems are difficult as observed in [1]. [21] considered a maximization of a linear function on reduced simplex. Their application came from financial stochastic dual dynamic programming.

In this paper, we also consider projections onto a modified simplex. Our motivation stems from the field of distributionally robust optimization [7] where one hedges against uncertainty by estimating a distribution $\bm{q}$ and considering the worse outcome when the true distribution $\bm{p}$ is not far away from $\bm{q}$. Since $\bm{p}$ is a distribution and we need to keep close to $\bm{q}$, this may be equivalently written as a projection with additional linear inequality. The projection is taken with respect to the distance between distributions and the additional linear inequality comes from considering the worst case. There are several possibilities of the considered distance function. Probably the most commonly used are the $\phi$-divergence [3] and the Wasserstein distance [9]. The authors in [11, 21] proposed an algorithm with quadratic complexity for the $l_{2}$ norm distance. In [22] the authors provided a closed-form algorithm for $l_{1}$ distance.

The resulting problem is convex and depending on the used distance function, it may be even quadratic or linear. This suggests using general-purpose solvers such as CPLEX which is guaranteed to converge. In our paper, we perform a comparison with CPLEX and IPOPT and show that our algorithm exhibit the observed (almost) linear complexity and outperform the above general-purpose solvers.

The paper is organized as follows. In Section 2 we give a brief introduction into the field of distributionally robust optimization. We derive the problems of interest and focus on measuring the distance by various $\phi$-divergences and $l_{p}$ norms. In Section 3 we present the main results. Instead of penalizing one constraint as in [14], we write the full KKT system and simplify it into one equation in one variable similar to (2). We derive the thresholding operator similar to (3). For readability, we postpone all proofs to the Appendix. In Section 4 we consider numerical properties and show that the derived equations have nice properties such as monotonicity or convexity. Finally, in Section 5 we focus on numerical results. All codes are available online.111https://github.com/VaclavMacha/Projections

2 Motivation from distributionally robust optimization

In this section, we provide motivation for the problems from the field of distributionally robust optimization. In the classical robust optimization, one minimizes a random function $f(\bm{x},\bm{\xi})$ with respect to a decision variable $\bm{x}$ while considering the worse possible random outcome of a random variable $\bm{\xi}$ which is bound to lie in $\Xi$. This leads to the problem

[TABLE]

Since (4) considers the worst possible scenario, this approach is usually too conservative. One way to alleviate it, is to consider the distributionally robust optimization, where one takes the worst outcome with respect to all probability distributions and not to all scenarios. Denoting the probability distribution by $P$, expectation with respect to $P$ by ${\mathbb{E}}_{P}$ and the set of all admissible probability distributions by ${\mathcal{P}}$, this results in

[TABLE]

Note that if ${\mathcal{P}}$ consists of all Dirac measures concentrated at $\Xi$, then (4) and (5) coincide.

The simplest case appears when we know possible realization $\bm{\xi}_{i}$ for the random variable and each may happen with the probability $p_{i}$. Then the inner maximization problem in (5) reduces to

[TABLE]

where we set $c_{i}:=f(\bm{x},\bm{\xi}_{i})$. However, the probability distribution $\bm{p}$ is often not known. Then we may assume that $\bm{p}$ is close to some known estimate $\bm{q}$ and we may want to hedge against the worst possible small deviation from $\bm{q}$. Measuring this deviation bu $\hat{d}$, the inner problem in (5) then takes form

[TABLE]

The first two constraints prescribe that $\bm{p}$ is a probability distribution while the last one determines that $\bm{p}$ is not far from $\bm{q}$.

2.1 Connection to projection onto the canonical simplex

Since convex programming allows to switch the objective and constraints, for each ${\varepsilon}$ there is some $\delta$ such that problem (7) is equivalent to

[TABLE]

This implies that problem (7) is equivalent to a projection (with respect to distance $\hat{d}$) onto the canonical simplex restricted by additional linear constraint. In Section 3.3 we will consider a similar problem where the canonical simplex is restricted by upper bounds.

2.2 Notation and assumptions

In the manuscript, we employ the following notation and assumption:

[TABLE]

Assumption 2.1**.**

We consider the following assumptions:

1. (A1)

Vector $\bm{q}$ has positive components which sum to one. 2. (A2)

Vector $\bm{c}$ is not a constant vector.

Assumption (A1) is natural since we want to measure the distance of $\bm{p}$ from $\bm{q}$ which is a distribution. Assumption (A2) is technical only. If $\bm{c}$ is a constant vector, then the objective $\bm{c}^{\top}\bm{p}={c_{\rm max}}$ is constant and the optimization is trivial.

3 Reduction of projections onto modified simplex to one equation

In the previous section, we mentioned how a projection onto the unit simplex with an additional constraint arises in the field of distributionally robust optimization. In the introduction, we recalled a way of solving the projection onto the unit simplex (2). Namely, one needs to solve the equation (2) and then apply the thresholding operator (3) to obtain the solution.

In this section, we follow a similar approach to solve problem (7) with two types of the distance function $\hat{d}$. Moreover, we use the same technique to derive an algorithm for projection onto the canonical simplex with additional upper bounds. These results allow us to propose numerical methods with linear complexity.

3.1 Distributionally robust optimization with $\phi$-divergences

In this section we consider the distributionally robust optimization where the distance function $\hat{d}$ is a $\phi$-divergence. In this case

[TABLE]

where $d$ is a $\phi$-divergence and $\phi$ is the convex generating function with $\phi(1)=0$. Then the distributionally robust problem (7) takes form

[TABLE]

Some examples of $\phi$-divergences are listed in Table 1. We also consider the variation distance ($l_{1}$ norm). Even though it is a $\phi$-divergence, we handle it in Section 3.2 as it is a norm as well. We start with the following result which states that if ${\varepsilon}$ is large enough, then the solution to (DRO1) is trivial.

Theorem 3.1**.**

Let Assumption 2.1 hold true and define vector $\hat{\bm{p}}$ with components

[TABLE]

If this solution satisfies

[TABLE]

then $\hat{\bm{p}}$ is the optimal solution of (DRO1).

Theorem 3.1 provides the optimal solution for the case of large ${\varepsilon}$. In the opposite case, we provide the solution in the following list of theorems. Each of them handles one $\phi$-divergence from Table 1.

Theorem 3.2** (Kullback-Leibler divergence).**

Assume that Assumption 2.1 holds true, that (11) is violated and that {\varepsilon}<-\log\mathopen{}\mathclose{{}\left(\sum_{i\in I}q_{i}}\right). Then there exists some $\mu\in(0,\frac{{c_{\rm max}}-{c_{\rm min}}}{{\varepsilon}}]$ which solves

[TABLE]

Moreover, defining the non-normalized weights

[TABLE]

the optimal solution of (DRO1) with the Kullback-Leibler divergence $\phi=\phi_{1}$ equals to $p_{i}=\frac{\hat{p}_{i}}{\sum\hat{p}_{j}}$.

Theorem 3.3** (Burg entropy).**

Assume that Assumption 2.1 holds true and that (11) is violated. Then there exists some $\lambda\in({c_{\rm max}},{c_{\rm max}}+\frac{{c_{\rm max}}-{c_{\rm min}}}{{\varepsilon}}]$ which solves

[TABLE]

Moreover, defining the non-normalized weights

[TABLE]

the optimal solution of (DRO1) with the Burg entropy distance $\phi=\phi_{2}$ equals to $p_{i}=\frac{\hat{p}_{i}}{\sum\hat{p}_{j}}$.

Theorem 3.4** (Hellinger distance).**

Assume that Assumption 2.1 holds true, that (11) is violated and that ${\varepsilon}<2-2\sqrt{\sum_{i\in I}q_{i}}$. Then there exists some $\lambda\in({c_{\rm max}},{c_{\rm max}}+\frac{(2-{\varepsilon})({c_{\rm max}}-{c_{\rm min}})}{{\varepsilon}}]$ which solves

[TABLE]

Moreover, defining the non-normalized weights

[TABLE]

the optimal solution of (DRO1) with the Hellinger distance $\phi=\phi_{3}$ equals to $p_{i}=\frac{\hat{p}_{i}}{\sum\hat{p}_{j}}$.

Theorem 3.5** ($\chi^{2}$-distance).**

Assume that Assumption 2.1 holds true and that (11) is violated. Then there exists some $\lambda>{c_{\rm max}}$ which solves

[TABLE]

Moreover, defining the non-normalized weights

[TABLE]

the optimal solution of (DRO1) with the $\chi^{2}$-distance $\phi=\phi_{4}$ equals to $p_{i}=\frac{\hat{p}_{i}}{\sum\hat{p}_{j}}$.

Theorem 3.6** (Modified $\chi^{2}$-distance).**

Assume that Assumption 2.1 holds true and that (11) is violated. Then there exists some $\lambda>-{c_{\rm max}}$ which solves

[TABLE]

Moreover, defining the non-normalized weights

[TABLE]

the optimal solution of (DRO1) with the modified $\chi^{2}$-distance $\phi=\phi_{5}$ equals to $p_{i}=\frac{\hat{p}_{i}}{\sum\hat{p}_{j}}$.

3.2 Distributionally robust optimization with norms

Anothor possibility to measure the distance between $\bm{p}$ and $\bm{q}$ is to use $l_{p}$ norms. This results in the following problem

[TABLE]

Note that the $l_{1}$ norm also generates a $\phi$-divergence but we handle it here. The simple algorithms for solving (DRO2) with the $l_{1}$ and $l_{\infty}$ norms are presented in Appendix A. Here, we show the results for the $l_{2}$ norm.

Theorem 3.7** ($l_{2}$ norm).**

Let Assumption 2.1 hold true and define vector $\hat{\bm{p}}$ with components

[TABLE]

If this solution satisfies

[TABLE]

then $\hat{\bm{p}}$ is the optimal solution of (DRO2). If $\hat{\bm{p}}$ violates (23), then there exists some $\mu>0$ and $\lambda$ which solve

[TABLE]

Moreover, the optimal solution of (DRO2) with the $l_{2}$ norm equals to

[TABLE]

Unlike in the previous cases, system (24) consists of two equations. To reduce them to one, we define first a function $g_{6}(\lambda;\mu)$ of $\lambda$ with fixed parameter $\mu$ by

[TABLE]

Lemma B.1 states that for each $\mu>0$, there is unique $\lambda$ solving $g_{6}(\lambda;\mu)=0$. We stress this dependence of $\lambda$ on $\mu$ by writing $\lambda(\mu)$. Moreover, the same lemma states that $\lambda(\mu)$ is a continuous function. Defining the continuous function

[TABLE]

we observe that solving system (24) can be reduced to solving the single equation $h_{6}(\mu)=0$.

3.3 Projection onto simplex with additional linear equality

In Section 2.1 we argued that (DRO2) with $p=2$ is equivalent to projecting onto the canonical simplex with additional linear inequality. In this section, we consider one more problem of the projection onto the simplex with additional upper bounds. This problem in a slightly more general form reads

[TABLE]

We obtain the reduced optimality conditions as follows, where $\operatorname{clip}$ is the projection operator.

Theorem 3.8** (Simplex with upper bounds).**

Assume that the feasible set of (SIMPLEX) is non-empty. Then there exists some $\lambda\in{\mathbb{R}}$ which solves

[TABLE]

Moreover, the optimal solution of (SIMPLEX) equals to

[TABLE]

Function $h_{7}$ is piecewise linear and non-increasing in $\lambda$. This allows us to find a simple algorithm to find the solution, we present it in Algorithm 1.

4 Numerical considerations

In the previous section, we derived theoretical results which will be the bases for numerical methods for solving (DRO1), (DRO2) and (SIMPLEX). In this section, we introduce these numerical methods and derive their complexity.

4.1 Computation of $\lambda$

For problems (DRO1) including the $\phi$-divergences and for (SIMPLEX), we reduced the optimality conditions into one equation in one variable. For (DRO2) with $l_{2}$ norm, we reduced it into two equations in two variables from which $\lambda$ is implicitly computed, further reducing the system into one equation in one variable. It is not difficult to show that this implicit equation (24b) is equivalent to

[TABLE]

Since (30) is identical to (2), there are algorithms in $O(n)$ which compute $\lambda$ for any fixed $\mu$. For simplicity, we implemented a simpler algorithm which first sorts $\mu\bm{q}+\bm{c}$ and then finds $\lambda$ in one pass through the sorted array.

For the analysis of (SIMPLEX) we realize that $h_{7}$ is a piecewise linear function which is non-increasing in $\lambda$. Since this problem differs from (2) only by the upper bound, we conjecture that there is an algorithm solving it in $O(n)$. Here, we present an algorithm with complexity $O(n\log n)$. We have

[TABLE]

The inequality holds due to the assumption that the feasible set of (SIMPLEX) is non-empty.

Denote $\bm{s}$ the sorted version of $\bm{q}-\bm{l}$ and $\bm{r}$ the sorted version of $\bm{q}-\bm{u}$. The main idea of Algorithm 1 is to utilize the piecewise linearity of $h_{7}$ with kinks at $s_{i}$ and $r_{j}$ by tracking the current slope $\hat{a}$. Algorithm 1 is an iterative procedure where at every iteration, we know the values of $h_{7}(s_{i-1})$ and $h_{7}(r_{j-1})$ and we want to evaluate $h_{7}$ at the next point. If $r_{j}\leq s_{i}$, then we consider $\lambda=r_{j}$ and increase $j$ by one. Since a new point enters the active set which contributes to the slope, we increase the slope $\hat{a}$ by $1$. If $r_{j}>s_{i}$, then we consider $\lambda=s_{i}$ and increase $i$ by one. Since one point leaves the active set, we decrease the slope $\hat{a}$ by $1$. In both cases, $g$ is decreased by $\hat{a}$ times the difference between the old value and the new values of $\lambda$. Once $g$ decreases below [math], we stop the algorithm and linearly interpolate between the last two values. To prevent an overflow, we set $r_{n+1}=\infty$. Concerning the initial values, since $r_{1}<s_{1}$, we set $i=1$ and $j=2$.

4.2 Numerical methods

From the proofs of Theorems 3.2-3.4 we obtain

[TABLE]

This implies that the bisection method is convergent for solving $h_{1}(\mu)=0$, $h_{2}(\lambda)=0$ and $h_{3}(\lambda)=0$ when starting from the bounds suggested by (31).

To solve $h_{4}(\lambda)=0$, $h_{5}(\lambda)=0$ and $h_{6}(\mu)=0$ we first observe that convexity is present due to the following result.

Proposition 4.1**.**

We have the following:

• •

If the assumptions of Theorem 3.5 are satisfied, then $h_{4}$ is decreasing and convex on $({c_{\rm max}},\infty)$.

• •

If the assumptions of Theorem 3.6 are satisfied, then $h_{5}$ is positive on $(-{c_{\rm max}},\lambda_{0})$ and decreasing and concave on $(\lambda_{0},\infty)$ for some $\lambda_{0}$.

• •

If the assumptions of Theorem 3.7 are satisfied, then $h_{6}$ is positive on $(0,\mu_{0})$ and decreasing and concave on $(\mu_{0},\infty)$ for some $\mu_{0}$.

Lemma C.1 states that if we start with a point with $h_{4}(\lambda)>0$, $h_{5}(\lambda)<0$ or $h_{6}(\mu)<0$, respectively, the Newton’s method is convergent. We summarize this discussion in Tables 2 and 3. The former shows which problem has an exact algorithm and which needs to be solved via an iterative method. Note that convergence is guaranteed for each problem. The latter comments more on the iterative methods and summarizes the equations needed to solve (DRO1) and (DRO2) with the $l_{2}$ norm. Moreover, it provides the bounds within which the solution lies and shows whether the function in question possesses convexity.

We summarize the whole procedure in Algorithm 2. The bisection method may be initialized based on the bounds from Table 3 while the Newton’s method must be initialized based on the paragraph following Proposition 4.1.

4.3 Complexity

In Table 4 we show the complexity of the evaluation of $h_{1},\dots,h_{6}$ and the total complexity of the algorithm. For (DRO1), the evaluation of $h_{1},\dots,h_{5}$ has the complexity of $O(n)$. Similarly for (DRO2) with the $l_{2}$ norm, for every $\mu$, the computation of $\lambda(\mu)$ can be done in $O(n)$ as shown is Section 4.1. Thus, the evaluation of $h_{6}$ consumes $O(n)$ as well. We get the total complexity by multiplying this by the number of evaluations ${n_{h}}$ of $h$. In order to have a good performance, the number of evaluations ${n_{h}}$ needs to stay constant. This happened in our numerical experiments as the second row of Figure 2 shows. Moreover, ${n_{h}}$ is guaranteed to be constant for the bisection method whenever the bracketing interval stays constant. In the table we also included problem (DRO2) with $l_{1}$ and $l_{\infty}$ norms and SIMPLEX. We provide a comparison of the theoretical and the observed complexity in Table 5 later.

5 Numerical results

In this section, we present the numerical results. We recall that our codes are available online.${}^{\ref{footnote:codes}}$ In Section 3, we derived the monotonicity and convexity of functions $h_{1},\dots,h_{7}$ corresponding to problems (DRO1), (DRO2) and (SIMPLEX) and in Section 4, we argued that finding a zero of these functions should be easy. This is confirmed in Figure 1. We see that $h_{1}$, which corresponds to (DRO1) with Kullback-Leibler divergence, is decreasing and seems to be convex. Similarly, $h_{6}$ corresponding to (DRO2) with $l_{2}$ norm is first increasing and after approximately $\mu=14$ decreasing and concave. The convexity for $h_{1}$ was not proven while the concavity for $h_{6}$ follows from Proposition 4.1.

For numerical comparison, we randomly generated the initial data $\bm{q}$ and $\bm{c}$ and solved problems (DRO1), (DRO2) and (SIMPLEX). This was repeated hundred times and the results were averaged to remove random bias. The main comparison is presented in Figure 2. The left column corresponds to problem (DRO1) while the right column to problems (DRO2) and (SIMPLEX). The $x$ axis always depicts the number of input data $n$ chosen in the range $n\in[10^{3},10^{6}]$. The first row depicts the computational time in seconds. The second row depicts the number of evaluations of $h_{1},\dots,h_{6}$.

We observe that the number of evaluations of $h$ in the second row stays relatively constant. Coming back to Table 4, this implies that ${n_{h}}$ is constant and the total complexity should be $O(n)$ or $O(n\log n)$. This is confirmed in the first row of Figure 2 where we see the (approximately) linear dependence of time on the data size $n$. To give a more quantitative result, we have interpolated the measured times with function $t(n)=an^{b}$ for the best possible parameters $a$ and $b$. We show this interpolation in Table 5. We see that the interpolation is close to linear. Note that the larger power of $n$ may hide the logarithm as the domain for $n$ is bounded.

In Figure 3 we compare our results to other solvers. We used the package JuMP for Julia [8]. It employed IPOPT for (DRO1) and CPLEX for (DRO2) and (SIMPLEX). Moreover, for (DRO2) with $l_{2}$ norm we compared ourselves to the algorithm from [21]. To keep the computation possible, we had to reduce the number of points from $n=10^{6}$ to $n=10^{4}$. We can see that our algorithms perform significantly better. CPLEX and IPOPT seem to also possess linear complexity unlike the algorithm from [21] with quadratic complexity. Its complexity estimate from Table 5 was $5.138\cdot 10^{-7}n^{1.632}$.

Acknowledgements

This work was supported by National Natural Science Foundation of China (Grant No. 61850410534), the Program for Guangdong Introducing Innovative and Enterpreneurial Teams (Grant No. 2017ZT07X386), Shenzhen Peacock Plan (Grant No. KQTD2016112514355531) and the Grant Agency of the Czech Republic (18-21409S).

Appendix A Problem (DRO2) with $l_{1}$ and $l_{\infty}$ norm

Here we present algorithms for solving (DRO2) for the $l_{1}$ and $l_{\infty}$ norms. Since $\bm{q}$ is a probability distribution due to Assumption 2.1, if we increase some components of $\bm{q}$, we have to decrease some components of $\bm{q}$ by the same margin. The priority is on increasing coordinates of $\bm{q}$ with the lowest value of $\bm{c}$ while decreasing those with the largest value. We summarize this procedure in Algorithms 1 and 2. Sorting $\bm{c}$, the lowest values have the lowest indices and similarly for the largest values. Thus, we start with $i=1$ and $j=n$. Then we increase $q_{i}$ by a possible maximal margin $\delta_{1}$ and start decreasing $q_{j}$, $q_{j-1}$ and so on until the total reduction $\delta_{2}$ equals to $\delta_{1}$. After doing so, we increase $i$ by one and continue until $i=j$. Note that $\delta_{\rm dec}$ in Algorithm 1 measures the decrease of $p_{j}$ while $\delta_{\rm tot}$ in Algorithm 2 measures the total reduction of $p_{j},\dots,p_{n}$. The first one has to be bounded by ${\varepsilon}$ while the other one by $\frac{{\varepsilon}}{2}$.

Appendix B Proofs

In this section, we present all proofs. The proofs are divided into subsections as in the manuscript body. We omit the proof of Theorem 3.8 as it is similar to other proofs.

B.1 Optimality conditions for Theorems 3.1-3.6

We start with a general part which is common to all Theorems 3.1-3.6. Since $\phi$ is convex in $\phi$, since ${\varepsilon}>0$, since $d(q_{i},q_{i})=0$ and since $\bm{q}$ defines a -probability distribution, the Slater constraint qualification is satisfied at $\bm{q}$. Thus, problem (DRO1) is equivalent to its KKT optimality conditions. The Lagrangian for (DRO1) reads

[TABLE]

The minus in front of the first term needs be to present since (DRO1) is a maximization problem. The KKT conditions then amount to the optimality conditions

[TABLE]

Since $\mu\geq 0$, there are two possibilities.

Case 1 of $\mu=0$:

If $\mu=0$, then from (32a) we get $\lambda=\alpha_{i}+c_{i}$. This, together with $\alpha_{i}\geq 0$ and $\alpha_{i}p_{i}=0$ implies that $\alpha_{i}=0$ for $i\in I$ and that $p_{i}=0$ for $i\notin I$. Then the KKT system is satisfied if there exists a feasible $\bm{p}$ with $p_{i}=0$ for $i\notin I$ such that $\sum_{i=1}^{n}d(p_{i},q_{i})\leq{\varepsilon}$. This problem can be verified by solving the convex problem

[TABLE]

and checking whether its optimal value is smaller or equal than ${\varepsilon}-\sum_{i\notin I}d(0,q_{i})$.

Since for $\phi$-divergences, we have $d(p_{i},q_{i})=q_{i}\phi(\frac{p_{i}}{q_{i}})$, it is not difficult to verify that the ratio $\frac{p_{i}}{q_{i}}$ is constant. This implies that $\hat{\bm{p}}$ defined in (10) is the solution to (33). This finishes the proof of Theorem 3.1.

Case 2 of $\mu>0$:

In the opposite case we have $\mu>0$, which due to (32c) implies that the feasibility conditions change into

[TABLE]

We now split the proof into five parts for Theorem 3.2 -3.6.

Proof of Theorem 3.2.

The feasibility constraint may be due to Assumption 2.1 written as

[TABLE]

For now we assume that $p_{i}>0$ for all $i$ and remove this assumption later. This implies $\alpha_{i}=0$. Then the optimality condition (32a) reads

[TABLE]

from which we deduce

[TABLE]

Plugging (35) into the feasibility conditions (34) yields

[TABLE]

We can express $\mu$ from (36) via

[TABLE]

which together with (37) gives the final equation (12). The optimal probabilities (13) then follow from (35).

Now we need to remove the assumption of $p_{i}>0$. Recall that

[TABLE]

For its middle part we have

[TABLE]

which implies that

[TABLE]

We consider now the limit of $h_{1}(\mu)$ as $\mu\downarrow 0$. Due to the properties of the exponential function, for all $\alpha>0$, there is some $\mu_{0}$ such that for all $\mu\in(0,\mu_{0})$ we have

[TABLE]

This implies

[TABLE]

The right-most term is positive and independent of $\mu$ whenever $i\in I$ and $\alpha$ is sufficiently small due to the assumptions of Theorem 3.2. Moreover as

[TABLE]

for all $i\notin I$, we deduce that $h_{1}(\mu)\to\infty$ as $\mu\downarrow 0$. This combined with (38) and the continuity of $h_{1}$ implies that the equation $h_{1}(\mu)=0$ has a solution on $(0,\frac{{c_{\rm max}}-{c_{\rm min}}}{{\varepsilon}}]$. Since the optimality conditions are equivalent to problem (DRO1) due to convexity, the existence of solution also implies that the assumption of $p_{i}>0$ may be alleviated. ∎

Proof of Theorem 3.3.

The form of $\phi_{2}$ implies $p_{i}>0$ for all $i$, which further means $\alpha_{i}=0$. Then the optimality condition (32a) reads

[TABLE]

from which we deduce

[TABLE]

Plugging (39) into the feasibility conditions (34) yields

[TABLE]

We can express $\mu$ from (40) via

[TABLE]

which together with (41) gives the final equation (14). The optimal probabilities (15) then follow from (39). The constraint $\mu>0$ transfers to $\lambda>{c_{\rm max}}$ due to (39).

Now we are interested in the limits. Recall that

[TABLE]

Then we have

[TABLE]

Thus

[TABLE]

We consider now the limit of $h_{2}(\lambda)$ as $\lambda\downarrow{c_{\rm max}}$. Denoting ${c_{\rm max2}}$ the second largest distinct component value of $\bm{c}$, we have

[TABLE]

Due to the assumptions of Theorem 3.6, we obtain $h_{2}(\lambda)\to\infty$ as $\lambda\downarrow{c_{\rm max}}$. This combined with (42) and the continuity of $h_{2}$ implies that the equation $h_{2}(\lambda)=0$ has a solution on $({c_{\rm max}},{c_{\rm max}}+\frac{{c_{\rm max}}-{c_{\rm min}}}{{\varepsilon}}]$. ∎

Proof of Theorem 3.4.

The feasibility constraint may due to Assumption 2.1 be written as

[TABLE]

For now we assume that $p_{i}>0$ for all $i$ and remove this assumption later. This implies $\alpha_{i}=0$. Then the optimality condition (32a) reads

[TABLE]

from which we deduce

[TABLE]

Plugging (44) into the feasibility conditions (34) yields

[TABLE]

We can express $\mu$ from (45) via

[TABLE]

which together with (46) gives the final equation (16). The optimal probabilities (17) then follow from (44). The constraint $\mu>0$ transfers to $\lambda>{c_{\rm max}}$ due to (43), which also allows us to remove the absolute value from (46).

Now we need to remove the assumption of $p_{i}>0$. Recall that

[TABLE]

Then we have

[TABLE]

Thus

[TABLE]

We consider now the limit of $h_{3}(\lambda)$ as $\lambda\downarrow{c_{\rm max}}$. Denoting ${c_{\rm max2}}$ the second largest distinct component value of $\bm{c}$, we have

[TABLE]

Due to the assumptions of Theorem 3.4, we obtain $h_{3}(\lambda)\to-\infty$ as $\lambda\downarrow{c_{\rm max}}$. This combined with (47) and the continuity of $h_{3}$ implies that the equation $h_{3}(\lambda)=0$ has a solution on $({c_{\rm max}},{c_{\rm max}}+\frac{(2-{\varepsilon})({c_{\rm max}}-{c_{\rm min}})}{{\varepsilon}}]$. Since the optimality conditions are equivalent to problem (DRO1) due to convexity, the existence of solution also implies that the assumption of $p_{i}>0$ may be alleviated. ∎

Proof of Theorem 3.5.

The form of $\phi_{4}$ implies $p_{i}>0$ for all $i$, which further means $\alpha_{i}=0$. Due to Assumption 2.1, the feasibility constraints on distance amounts to

[TABLE]

Then the optimality condition (32a) reads

[TABLE]

from which we deduce

[TABLE]

Plugging (48) into the feasibility conditions (34) yields

[TABLE]

We can express $\mu$ from (50) via

[TABLE]

which together with (49) gives the final equation (18). The optimal probabilities (19) then follow from (48). The constraint $\mu>0$ transfers to $\lambda>{c_{\rm max}}$ due to (48). ∎

Proof of Theorem 3.6.

The feasibility constraint may due to Assumption 2.1 be written as

[TABLE]

Then the optimality condition (32a) reads (for the uniformity of results, we flipped the sign of $\lambda$)

[TABLE]

from which we due to the complementarity conditions (32b) deduce

[TABLE]

Plugging (51) into the feasibility conditions (34) yields

[TABLE]

We can express $\mu$ from (52) via

[TABLE]

which together with (53) gives the final equation (20). The optimal probabilities (21) then follow from (51). The constraint $\mu>0$ transfers to $\lambda>-{c_{\rm max}}$ due to (51). ∎

B.2 Optimality conditions for Theorem 3.7

To prove Theorem 3.7, we first realize that the first part of Section B.1 including the first paragraph of “Case 1 of $\mu=0$” holds true with $d(p_{i},q_{i})=\frac{1}{2}(p_{i}-q_{i})^{2}$. Moreover, problem (33) takes the form

[TABLE]

Since $q_{i}>0$ and $\sum_{i\in I}q_{i}<1$ due to Assumption 2.1, it is not difficult to verify that $\hat{\bm{p}}$ defined in (22) is the optimal solution of (54). This proves the first part of Theorem 3.7.

For the second part, we realize that the feasibility constraint may be written as

[TABLE]

Then the optimality condition (32a) reads

[TABLE]

from which we due to the primal feasibility condition $p_{i}\geq 0$, the dual feasibility conditions $\alpha_{i}\geq 0$ and the complementarity condition (32b) deduce

[TABLE]

Plugging (55) into the feasibility conditions (34) yields

[TABLE]

The final equations (24) are obtained from (56) and (57) by simple calculus, the formula $\max(-x,-y)=-\min(x,y)$ and Assumption 2.1.

B.3 Convexity for Proposition 4.1

We first show an auxiliary result which states the continuity on $h_{6}$.

Lemma B.1**.**

For each $\mu>0$ there is a unique $\lambda(\mu)$ which solves $g_{6}(\lambda;\mu)$. Moreover, function $h_{6}(\mu)$ is continuous on $(0,\infty)$.

Proof.

Fix any $\mu>0$ and recall that

[TABLE]

Since $\mu>0$ for at least one $i$ we have $\lambda(\mu)-c_{i}<\mu q_{i}$ which implies that $g_{6}$ is strictly increasing in $\lambda$ around $\lambda(\mu)$. Thus, the solution $\lambda(\mu)$ is unique.

Consider now any $\mu_{k}\to\mu>0$. Since the corresponding $\lambda(\mu_{k})$ are bounded, we may select a converging subsequence, say $\lambda(\mu_{k})\to\lambda^{*}$. Then

[TABLE]

which implies that $\lambda^{*}=\lambda(\mu)$. Thus, $\lambda$ is a continuous function of $\mu$. This further implies that $h_{6}$ is a continuous function of $\mu$. ∎

Concerning the proof of Proposition 4.1, we divide it into three parts.

Proof for $h_{4}$:

After rearranging of terms, we get

[TABLE]

where for $i,j=1,\dots,n$ we define

[TABLE]

Computing the first derivative

[TABLE]

and the second derivative

[TABLE]

we realize that

[TABLE]

This together with (58) implies that $h_{4}$ is convex.

Proof for $h_{5}$:

Recall that

[TABLE]

For simplicity assume that $-c_{1}<\dots<-c_{n}$. Then on $\lambda\in(-c_{j},-c_{j+1})$ we have

[TABLE]

The derivative at this interval equals to

[TABLE]

Due to the assumption of $-c_{1}<\dots<-c_{n}$, the violation of (11) implies $(1+{\varepsilon})q_{1}<1$. This due to (60) means that $h_{4}^{\prime}(\lambda)>0$ on $\lambda\in(-c_{1},-c_{2})$. Since $h_{5}(-c_{1})=0$ due to (59), this implies that $h_{5}$ is positive on $(-c_{1},-c_{2})$. Moreover, $h_{5}^{\prime}$ equals to a product of two terms, the first one is always positive while the second one is piecewise constant with decreasing values on individual pieces. This implies that there is some $\lambda_{0}$ such that $h_{5}$ is non-decreasing on $(-{c_{\rm max}},\lambda_{0})$ and decreasing on $(\lambda_{0},\infty)$. By the same arguments we have that $h_{5}^{\prime}$ is decreasing on $(\lambda_{0},\infty)$, which implies that $h_{5}$ is concave on this interval. This finishes the proof of the second part.

If $-c_{1}\leq\dots\leq-c_{n}$ instead of the assumed $-c_{1}<\dots<-c_{n}$, then the proof can be performed in exactly the same way but $(1+{\varepsilon})q_{1}<1$ changes to $(1+{\varepsilon})\sum_{i=1}^{j}q_{i}<1$, where $j$ is the cardinaly of $I$. Then we can get the same estimates of (60) as in the previous paragraphs.

Proof for $h_{6}$:

Define

[TABLE]

and consider any $0<\mu_{1}<\mu_{2}$. Since $\lambda(\mu_{1})$ solves (24b) for $\mu=\mu_{1}$ and $\lambda(\mu_{2})$ solves the same equation for $\mu=\mu_{2}$, we obtain $\lambda(\mu_{1})\geq\lambda(\mu_{2})$. This due to (61) implies $I(\mu_{1})\subset I(\mu_{2})$. Thus, $I(\mu)$ is a non-decreasing function (with respect to set inclusion) of $\mu$.

From (24b) we have

[TABLE]

from which we deduce

[TABLE]

Since $I(\mu)$ is a non-decreasing function of $\mu$, this implies that $\lambda(\mu)$ is a piecewise linear function with a finite number of pieces. On each of these pieces, we have

[TABLE]

and consequently

[TABLE]

Differentiating this relation yields

[TABLE]

where in the third equality we used (62) and in the last one (63).

Denote ${c_{\rm max2}}$ the second largest distinct component value of $\bm{c}$ and define $\hat{\lambda}=\frac{1}{2}({c_{\rm max}}+{c_{\rm max2}})$. Fix any $\mu$ sufficiently small but positive. Then we have

[TABLE]

whenever $\mu$ is sufficiently small. Since $g_{6}(\lambda(\mu);\mu)=0$ due to definition and since $g_{6}$ is non-decreasing in $\lambda$, this means that $\lambda(\mu)\geq\hat{\lambda}=\frac{1}{2}({c_{\rm max}}+{c_{\rm max2}})$. This due to (61) implies $I(\mu)=I$ whenever $\mu$ is sufficiently small. It is possible to show that the violation of (23) is equivalent to

[TABLE]

Combining this with (64) and $I(\mu)=I$, this implies that $h_{6}$ is strictly increasing on some interval $(0,\mu_{1})$. Moreover, since $\lambda(0)={c_{\rm max}}$, we have $I(0)=\emptyset$ and thus $h_{6}(0)=0$. This implies that $h_{6}$ is positive on $(0,\mu_{1})$.

Moreover, $h_{6}^{\prime}$ equals to a product of two terms, the first one is always positive while the second one is piecewise constant with decreasing values on individual pieces. This implies that there is some $\mu_{0}$ such that $h_{6}$ is non-decreasing on $(0,\mu_{0})$ and decreasing on $(\mu_{0},\infty)$. By the same arguments, we have that $h_{6}^{\prime}$ is decreasing on $(\mu_{0},\infty)$, which implies that $h_{6}$ is concave on this interval. This finishes the proof of the second part.

Appendix C Convergence of Newton’s method

We show the following theorem only for differentiable functions. However, it holds for any concave function by replacing the derivative by its concave superdifferential.

Lemma C.1**.**

Consider a continuous concave function $h:[a,b]\to{\mathbb{R}}$ with $h(a)>0$ and $h(b)<0$. Then the Newton’s method

[TABLE]

started from any point $\lambda^{0}$ with $h(\lambda^{0})<0$ gives a decreasing sequence which converges to some $\bar{\lambda}\in(a,b)$ with $h(\bar{\lambda})=0$.

Proof.

Due to concavity of $h$ and $h(a)>0$ with $h(b)<0$, there exists unique $\lambda^{*}\in(a,b)$ with $h(\lambda^{*})=0$. Moreover, due to concavity again we have $h^{\prime}(\lambda^{*})<0$. Since $h$ is concave and since $h(a)>0$ and $h(\lambda^{0})<0$, we obtain $\lambda^{0}>\lambda^{*}$. Then the Newton’s method forms a decreasing sequence $\{\lambda^{k}\}$ bounded below by $\lambda^{*}$, which is therefore convergent. At the same time, $h^{\prime}(\lambda^{k})$ is uniformly bounded above by zero as $h^{\prime}(\lambda^{k})<h^{\prime}(\lambda^{*})<0$. This implies

[TABLE]

which due to the uniform boundedness of $h^{\prime}(\lambda^{k})$ from zero implies $h(\lambda^{k})\to 0$. But this due to the continuity of $h$ means that $\lambda^{k}\to\bar{\lambda}$. ∎