On the Minimization of Convex Functionals of Probability Distributions Under Band Constraints

TL;DR

This paper develops theoretical optimality conditions and proposes two numerical algorithms for minimizing convex functionals of probability distributions with bounded densities, ensuring convergence and efficiency.

Contribution

It introduces new first-order optimality conditions and two algorithms for constrained convex functional minimization of probability distributions.

Findings

Algorithms efficiently solve the problem with high accuracy.

The second algorithm guarantees convergence under mild assumptions.

Numerical examples demonstrate theoretical and practical effectiveness.

Abstract

The problem of minimizing convex functionals of probability distributions is solved under the assumption that the density of every distribution is bounded from above and below. A system of sufficient and necessary first-order optimality conditions as well as a bound on the optimality gap of feasible candidate solutions are derived. Based on these results, two numerical algorithms are proposed that iteratively solve the system of optimality conditions on a grid of discrete points. Both algorithms use a block coordinate descent strategy and terminate once the optimality gap falls below the desired tolerance. While the first algorithm is conceptually simpler and more efficient, it is not guaranteed to converge for objective functions that are not strictly convex. This shortcoming is overcome in the second algorithm, which uses an additional outer proximal iteration, and, which is proven to converge under mild assumptions. Two examples are given to demonstrate the theoretical usefulness of the optimality conditions as well as the high efficiency and accuracy of the proposed numerical algorithms.

Figures

Equations

$$\int_{\Omega}f\biggl{(}\omega,\frac{dP_{1}}{dP}(\omega),\ldots,\frac{dP_{N}}{dP}(\omega)\biggr{)}\,\mathrm{d}P(\omega),$$

$${\partial }_{x_n}f(\mathbf{x}):=\left\{\gamma \in \mathbb{R}:\frac{f(\mathbf{y})-f(\mathbf{x})}{y_n-x_n}\le \gamma \forall \mathbf{y}\in \mathcal{X}\setminus \{\mathbf{x}\}\right\}.$$

$$f^{-1}(c)=\inf \left\{x\in \mathcal{X}:f(x)\ge c\right\}.$$

$$f^{-1}(c)$$

$${\mathcal{P}}^{=}=\{p\in {\mathcal{L}}_{\mu }^1:p^{\prime }\le p\le p^{\prime \prime }\},$$

$$0\le p^{\prime }\le p^{\prime \prime }\le \infty ,{\int }_{\Omega }{p}^{\prime }\mathrm{d}\mu \le 1,{\int }_{\Omega }{p}^{\prime \prime }\mathrm{d}\mu \ge 1.$$

$$I_f(p_1,\dots ,p_N):={\int }_{\Omega }f(\omega ,p_1(\omega ),\dots ,p_N(\omega ))\mathrm{d}\mu (\omega ),$$

$$f:\Omega \times [0,\infty {)}^N$$

$$(\omega ,x_1,\dots ,x_N)$$

$$f(\mathbf{x}):=f(x_1,\dots ,x_N):=f(\omega ,x_1,\dots ,x_N).$$

$$\underset{\{p_n\in {\mathcal{P}}_n^{=}{\}}_{n=1}^N}{\min }{I}_f(p_1,\dots ,p_N),$$

$$\underset{\{p_n\in {\mathcal{L}}_{\mu }^1{\}}_{n=1}^N}{\min }{\int }_{\Omega }f(\omega ,p_1(\omega ),\dots ,p_N(\omega ))\mathrm{d}\mu (\omega )$$

$$\text{s.t.}{p}_n^{\prime }\le p_n\le p_n^{\prime \prime },{\int }_{\Omega }{p}_n\mathrm{d}\mu =1,n=1,\dots ,N.$$

$$q_n=\{\begin{array}{ll}{p}_n^{\prime \prime },& f_n(\mathbf{q})<c_n\\ f_n^{-1}({\mathbf{q}}_{[n]},c_n),& f_n(\mathbf{q})=c_n\\ p_n^{\prime },& f_n(\mathbf{q})>c_n\end{array}$$

$$f_n(\mathbf{x})\in {\partial }_{x_n}f(\mathbf{x}),$$

$$q_n=\min \{p_n^{\prime \prime },\max \{f_n^{-1}({\mathbf{q}}_{[n]},c_n),p_n^{\prime }\}\}$$

$$I_f(\mathbf{p})-I_f(\mathbf{q})\le ⟨\mathbf{e},1⟩$$

$$e_n^{\prime \prime }$$

$$e_n^{\prime }$$

$$g(\omega )\approx g_{\mathbf{a}}(\omega ):=\sum _{k=1}^K{a}_k{\psi }_k(\omega -{\omega }_k),$$

$${\mu }_k:={\int }_{\Omega }{\psi }_k(\omega )\mathrm{d}\mu (\omega )<\infty ,k=1,\dots ,K.$$

$${\psi }_k({\omega }_l-{\omega }_k)=\{\begin{array}{ll}1,& l=k\\ 0,& l\ne k\end{array}$$

$$g_{\mathbf{a}}({\omega }_k)=a_k,k=1,\dots ,K.$$

$$p_n(\omega )\approx p_{{\mathbf{a}}_n}(\omega ):=\sum _{k=1}^K{a}_{n,k}{\psi }_k(\omega -{\omega }_k),n=1,\dots ,N,$$

$$\underset{\mathbf{A}\in {\mathbb{R}}^{N\times K}}{\min }{\int }_{\Omega }f(\omega ,p_{{\mathbf{a}}_1}(\omega ),\dots ,p_{{\mathbf{a}}_N}(\omega ))\mathrm{d}\mu$$

$$\text{s.t.}{p}_n^{\prime }(\mathbf{\omega })\le {\mathbf{a}}_n\le p_n^{\prime \prime }(\mathbf{\omega }),⟨{\mathbf{a}}_n,\mathbf{\mu }⟩=1,n=1,\dots ,N.$$

$${\mathbf{a}}_n=\min \{p_n^{\prime \prime }(\mathbf{\omega }),\max \{f_n^{-1}({\mathbf{A}}_{[n]},c_n),p_n^{\prime }(\mathbf{\omega })\}\},$$

$$f_n^{-1}({\mathbf{A}}_{[n]},c_n)=\inf \{{\mathbf{a}}_n\in {\mathbb{R}}^K:f_n(\mathbf{A})\ge c_{n}1\},$$

$${\tilde{e}}_n^{\prime \prime }$$

$${\tilde{e}}_n^{\prime }$$

$$n^{\ast }\in \underset{n=1,\dots ,N}{\mathrm{argmax}}\{{\tilde{e}}_n\}.$$

$$⟨\min \{p_{n^{\ast }}^{\prime \prime }(\mathbf{\omega }),\max \{f_{n^{\ast }}^{-1}({\mathbf{A}}_{[n^{\ast }]},c),p_{n^{\ast }}^{\prime }(\mathbf{\omega })\}\},\mathbf{\mu }⟩=1.$$

Full text

I Introduction

Functionals of probability distributions play a central role in probability theory and statistics. To clarify, a functional is a mapping from a vector space to the real line, i.e., a function which maps an element of a possibly high- or infinite-dimensional space to a scalar value. Omnipresent examples are the moments of a real-valued random variable, which map a distribution to a real number.

Convex functionals of probability distributions occur naturally in problems of statistical inference and decision making. In general, the expected cost of any inference procedure with a convex loss function can be shown to be a convex functional of probability distributions [1]. Consequently, examples of convex functionals can be found in detection [2], estimation [3], and joint detection and estimation theory, as well as in Bayesian inference [4].

In practice, it is often the case that the distributions of random variables that describe a random phenomenon are not known exactly, but are subject to uncertainty. This uncertainty can, for example, be caused by a lack of information about the random phenomenon or the absence of an appropriate model for its mathematical description. Uncertainty can also be introduced intentionally in order to safeguard an inference procedure against deviations from the underlying assumptions. The type and degree of uncertainty is usually specified by means of an uncertainty set that contains all feasible distributions. Given such a set, a question that naturally arises is which distributions are most favorable and which are least favorable. In order to answer this question, a given cost function has to be minimized or maximized over the uncertainty set. This is the problem addressed in this paper.

The uncertainty model that is assumed to hold throughout the paper restricts the densities of feasible distributions to lie within a band that is defined by two non-intersecting functions, which bound the density from above an below. It can be thought of as a confidence interval for the true density function. This model is known as the density band model and is commonly used in robust statistics [5]. In the context of robust hypothesis testing it was first studied in [6] and was recently revisited in [7]. It is discussed in more detail in Section II.

The problem that motivated the work in this paper relates to the design of minimax optimal sequential tests for multiple hypotheses. The least favorable distributions for this type of test can be shown to be minimizers of functionals of the form111 For binary tests, this result can be found in [8]. The corresponding results on tests for multiple hypotheses are to be presented in a forthcoming publication.

[TABLE]

where $P,P_{1},\ldots,P_{N}$ denote distributions that are subject to uncertainty and the function $f$ is jointly convex in the likelihood ratios $dP_{n}/dP$. A particular difficulty that arises in the design of minimax sequential tests is that $f$ is itself the solution of an optimization problem so that its value and its derivatives can only be evaluated numerically. Minimizing (1) over $P,P_{1},\ldots,P_{N}$, and under band constraints, is a challenging task because it involves two nested optimization problems. Deriving analytic solutions, or determining approximate solutions is usually not possible. Therefore, a suitable numerical algorithm is required that is:

• •

accurate enough to closely approximate continuous density functions,

• •

efficient enough to handle multiple distributions and reasonably fine grids for their discretization;

• •

parallelizable, in order to leverage modern hardware;

• •

and robust against mild numerical noise in the evaluation of $f$ and its derivatives.

Looking into existing convex optimization frameworks and algorithms, it was found that most off-the-shelf methods did not satisfy these requirements. On the one hand, high performance solvers such as Gurobi, MOSEK, or CPLEX were found to be too restrictive in terms of feasible objective functions. On the other hand, commonly used generic convex optimization algorithms, such as interior point, steepest descent or conjugate gradient methods, turned out to be too inefficient to be useful in practice. In addition, all algorithms considered in the survey suffered from severe accuracy issues, especially in the tails of the optimal distributions.

In this work, an approach for the minimization of convex functionals of probability distributions, under density band uncertainty, is detailed, which is efficient, reliable, offers control over the achieved accuracy, and is applicable beyond the particular use case of robust sequential hypothesis testing. The proposed algorithms offer a good trade-off between being generic and specific, in the sense that they heavily exploit assumptions about the structure of the objective function and the constraints, while at the same time providing enough flexibility to be applicable to a large class of problems in statistical signal processing and robust statistics in particular.

The existing literature on convex functionals is large and dispersed. While in the functional analysis literature the term convex functionals prevails [9], in statistics, signal processing, and information theory, similar classes of functions often go by the names divergence, distance, dissimilarity, or disparity. Early results on the subject are due to, among others, Pearson [10], Mahalanobis [11], Shannon [12], and Kullback [13]—see [14, 15] and references therein for a detailed treatment. The minimization of convex functionals, and their relation to robust decision making, has been addressed by Huber [16], Poor [17], Kassam [6], and Guntuboyina [18] to name just a few. An approach similar to the one presented here was used in [19] to solve a constrained minimization problem with a separable cost function. However, neither is the problem in [19] a special case of the problem investigated in this paper, nor vice versa.

The paper is organized as follows: the density band uncertainty model is briefly reviewed in Section II. In Section III, the functional minimizaton problem is stated in a formal manner. Necessary and sufficient optimality conditions as well a bound on the optimality gap of a feasible candidate solution are stated in Section IV. Both are useful results in their own right and constitute a main contribution of the paper. In Section V, the functional optimization problem is discretized and an algorithm is proposed that iteratively solves the optimality conditions using a block coordinate descent (BCD) strategy. Since the latter is not guaranteed to converge to a global minimum for general convex objective functions, a second algorithm is introduced that augments the objective function with an additional proximal term. Guaranteed convergence is shown for the proximal algorithm. In Section VI, two examples are provided to illustrate how the optimality conditions as well as the proposed numerical algorithms can be used in practice. Section VII concludes the paper.

Notation: Probability distributions are denoted by upper case letters, their densities by the corresponding lower case letters. Boldface lower case letters $\bm{x}$ are used to indicate row vectors and boldface upper case letters $\bm{X}$ to indicate matrices. The notations $\bm{x}_{[n]}$ and $\bm{X}_{[n]}$ are used to denote a vector whose $n$th element has been removed and matrices whose $n$th row has been removed, respectively. The inner product of two vectors $\bm{x}$ and $\bm{y}$ is denoted by $\langle\bm{x},\bm{y}\rangle$, the element-wise product by $\bm{xy}$. All comparisons between vectors are defined element-wise. The all-ones vector is denoted by $\bm{1}$. In the pseudocode of the algorithms, an R-style arrow notation $x\leftarrow y$ is used to assign a value $y$ to a variable $x$. $\mathcal{L}_{\mu}^{p}$ denotes the space of functions whose $p$th power of the absolute value is integrable with respect to the measure $\mu$. The operators $(f)^{+}$ and $(f)^{-}$ denote the positive and negative parts of a function $f$, i.e., $(f(\omega))^{+}=\max\{f(\omega),0\}$ and $(f(\omega))^{-}=\min\{f(\omega),0\}$. The notation $f(\bm{\omega})$ is used as shorthand for the vector $(f(\omega_{1}),\ldots,f(\omega_{N}))$. For functions that are defined directly on the sample space, the explicit argument $\omega$ is often omitted for the sake of a more compact notation, in particular if the function is integrated or is itself an argument of a higher-order function. Comparisons between functions are defined point-wise. The notation $\partial_{x_{n}}f(\bm{x})$ is used for the partial subdifferential of a convex function $f\colon\mathcal{X}\subset\mathbb{R}^{N}\to\mathbb{R}$ with respect to $x_{n}$ at $\bm{x}$, i.e.,

[TABLE]

Finally, the generalized inverse [20] of a nondecreasing function $f\colon\mathcal{X}\subset\mathbb{R}\to\mathbb{R}$ is defined as

[TABLE]

Note that this implies that $f^{-1}$ is nondecreasing in $c$ and that for $\mathcal{X}=\mathbb{R}$

[TABLE]

II The Density Band Uncertainty Model

Let $(\Omega,\mathcal{F})$ be a measurable space and let $\mu$ be an absolutely continuous $\sigma$-finite measure on this space. The density-band uncertainty model specifies sets of the form

[TABLE]

where

[TABLE]

Thus, all feasible densities are upper bounded by $p^{\prime\prime}$ and lower bounded by $p^{\prime}$. The case where the upper bound is infinity and the lower bound is zero is the unconstrained case.

The density band uncertainty model is useful for several reasons. First, it provides a great amount of flexibility to the designer of an inference procedure, as it allows for varying local degrees of uncertainty on different regions of the sample space. Depending on the application, it can be constructed by hand, based on expert knowledge, or statistically, via confidence interval estimators. Second, in contrast to many parametric uncertainty models, it provides clear visualization and easy interpretation. Third, from a theoretical point of view, the density band model is of interest because it generalizes several popular uncertainty models such as the $\varepsilon$-contamination model [21] and the bounded distribution function model [22]. A more detailed discussion of the band model and its properties can be found in [6] and [7].

III Problem Formulation

Let $N$ be a positive integer. The functionals considered in this work are of the form

[TABLE]

where $p_{1},\ldots,p_{N}$ are probability densities on $(\Omega,\mathcal{F})$ and

[TABLE]

is a function that is convex with respect to $(x_{1},\ldots,x_{N})$. In order to guarantee that (5) is well defined, it is further assumed that $f(\omega,p_{1}(\omega),\ldots,p_{N}(\omega))$ is $\mu$-measurable for all feasible densities $p_{1},\ldots,p_{N}$. A proof, and a more in-depth analysis, of the existence and well-definedness of (5) is detailed in [9]. Also note that $f$ being convex with respect to $(x_{1},\ldots,x_{N})$ does not imply that $f(\omega,p_{1}(\omega),\ldots,p_{N}(\omega))$ is convex with respect to $\omega$. To facilitate compact notation, the arguments of $f$ are occasionally written in vector notation and the direct dependence on $\omega$ is omitted, i.e.,

[TABLE]

It is important to note that $I_{f}(p_{1},\ldots,p_{N})$ is used instead of $I_{f}(P_{1},\ldots,P_{N})$. The latter notation is commonly used in the context of distance measures between distributions, such as $f$-divergences [23] and $f$-dissimilarities [24], in order to emphasize that the distance does not depend on how the reference measure $\mu$ is chosen. However, since $f$ in (6) is not assumed to be homogeneous and is allowed to directly depend on $\omega$, this independence does not hold in general.

The optimization problem considered in this paper is

[TABLE]

where all $\mathcal{P}^{=}_{n}$ are of the form (3) and $\{p_{n}\in\mathcal{P}^{=}_{n}\}_{n=1}^{N}$ is used as a shorthand notation for $p_{n}\in\mathcal{P}^{=}_{n}$ for all $n=1,\ldots,N$. Expressing the objective function and constraints explicitly, (7) becomes

[TABLE]

In the next section, a sufficient and necessary condition for a density vector $\bm{q}=(q_{1},\ldots,q_{N})$ to be a solution of (8) is given as well as a bound on the optimality gap $I_{f}(\bm{p})-I_{f}(\bm{q})$ for a feasible density vector $\bm{p}$.

IV Optimality Conditions and Optimality Gap

The optimality conditions for (8) are given in the following theorem.

Theorem 1** (Optimality Conditions)**

A sufficient and necessary condition for the densities $\bm{q}=(q_{1},\ldots,q_{N})$ to be a solution of (8) is that they satisfy

[TABLE]

for some $c_{1},\ldots,c_{N}\in\mathbb{R}$ and all $n=1,\ldots,N$. Here $f_{n}$ denotes a partial subderivative of $f$ with respect to $x_{n}$, i.e.,

[TABLE]

and $f_{n}^{-1}(\bm{x}_{[n]},c)$ denotes the generalized inverse of $f_{n}(\bm{x})$ with respect to $x_{n}$ in the sense of (2).

A proof for Theorem 1 is detailed in Appendix A.

Making use of the fact that $f_{n}$ is a subderivative of a convex function and, hence, is nondecreasing, Theorem 1 can further be written in an alternative, more expressive form that eliminates the need for the explicit case-by-case definition.

Corollary 1** (Optimality Conditions)**

A sufficient and necessary condition for the densities $\bm{q}=(q_{1},\ldots,q_{N})$ to be a solution of (8) is that they satisfy

[TABLE]

for some $c_{1},\ldots,c_{N}\in\mathbb{R}$ and all $n=1,\ldots,N$.

Corollary 1 is proved in Appendix B. It gives an expression for $q_{n}$ solely in terms of the remaining optimal densities $\bm{q}_{[n]}$ and the scalar $c_{n}$. That is, knowing $\bm{q}_{[n]}$, the missing density $q_{n}$ can be found via a search over $c_{n}$. The iterative algorithms presented in the next section are based on this idea. It can further be seen from (9) that $q_{n}$ is a projection of $f_{n}^{-1}(\bm{q}_{[n]},c_{n})$ onto the band of feasible densities $\mathcal{P}^{=}_{n}$. In the limit, i.e., $p^{\prime\prime}_{n}\to\infty$ and $p^{\prime}_{n}\to 0$, it follows that $q_{n}=f_{n}^{-1}(\bm{q}_{[n]},c_{n})$. The example in Section VI-A demonstrates this.

In practice, it might not be possible to find an exact analytic solution to the optimality conditions in Theorem 1 or Corollary 1. Typically, in such cases, a sequence of approximations is constructed until one of the candidate solutions is sufficiently close to the true optimum. In order to quantify the deviation from the optimum, it is useful to have an upper bound on the optimality gap, i.e., a bound on the difference between the value of the objective function at the candidate solution and at the exact solution. Such a bound is given in the next theorem.

Theorem 2** (Optimality Gap)**

Let $\bm{q}=(q_{1},\ldots,q_{N})$ be a solution of (8). For every feasible vector of densities $\bm{p}=(p_{1},\ldots,p_{N})$ and every $\bm{c}=(c_{1},\ldots,c_{N})\in\mathbb{R}^{N}$ it holds that

[TABLE]

where $\bm{e}=(e_{1},\ldots,e_{N})$, with $e_{n}=e^{\prime\prime}_{n}+e^{\prime}_{n}$ and

[TABLE]

A proof of Theorem 2 is detailed in Appendix C. The bound in Corollary 2 can be interpreted as a sum of $N$ residuals, where each residual $e_{n}$ is a measure for how much the density $p_{n}$ violates the complementary slackness constraint (39). In turn, $\langle\bm{e},\bm{1}\rangle=0$ implies that the optimality conditions in Theorem 1 are satisfied and that $\bm{p}$ solves (8).

In the next section, two numerical algorithms are presented that use the optimality conditions and the bound on the optimality gap to construct a sequence of increasingly accurate approximations to the solution of (8).

V Numerical Calculation of the Optimal Densities

For functions $f$ with analytically invertible subderivatives and simple bounds $p^{\prime\prime}_{n}$, $p^{\prime}_{n}$, it can be possible to solve the optimality conditions in Theorem 1 and Corollary 1 analytically. However, for more complicated problems, the optimality conditions need to be solved numerically. For this case, two efficient and numerically stable algorithms to iteratively approximate the solution of (8) are detailed. Both algorithms are based on fixed-point iterations that are solved via a block coordinate descent strategy. The proposed coordinate selection rule is based on the residuals in Theorem 2. While the first algorithm is more efficient and conceptually simpler, it is not guaranteed to converge for every feasible choice of $f$. This shortcoming is overcome in the second algorithm, where an additional outer proximal iteration is introduced to ensure convergence. Both algorithms use the bound on the optimality gap in Corollary 2 to track how closely the current iterate approximates the true minimum so that the iteration can be terminated when the desired level of accuracy is achieved.

V-A Discretization of the Continuous Problem

In order to make (8) tractable for numerical optimization techniques, it first needs to be reduced to a finite-dimensional problem. A common approach to represent a continuous function $g\colon\Omega\to\mathbb{R}$ by a finite-dimensional vector $\bm{a}\in\mathbb{R}^{K}$, $K\geq 1$, is to express it in terms of a linear combination of basis functions [25], i.e.,

[TABLE]

where $\psi_{1},\ldots\psi_{K}\colon\Omega\to\mathbb{R}$ denote $K$ basis functions centered at grid points $\bm{\omega}=(\omega_{1},\ldots,\omega_{K})$ and $\bm{a}=(a_{1},\ldots,a_{K})$ denote the combination weights.

In an optimization context, the basis functions are typically chosen a priori, while the weights are subject to the optimization. For the problem considered in this paper, all basis functions need to be chosen such that they are $\mu$-integrable, i.e.,

[TABLE]

The masses of the basis functions are collected in a vector $\bm{\mu}=(\mu_{1},\ldots,\mu_{K})$. Moreover, it is useful to require the basis functions to satisfy

[TABLE]

for all $k=1,\ldots,K$. This assumption is relatively mild and is satisfied by many common sets of basis functions, in particular by appropriately chosen M-splines and squared sinc functions. It guarantees that on the grid points the approximated function $g_{\bm{a}}$ evaluates to the corresponding combination weight, i.e.,

[TABLE]

This property decouples the combination weights and enables a parallel implementation of the algorithms presented in the following sections. Note that (13) does not imply orthogonality.

V-B Discrete Optimality Conditions

In order to obtain a discrete version of the optimality conditions in Theorem 1 and Corollary 1, all densities $p_{1},\ldots,p_{N}$ are again expressed as combinations of a fixed set of basis functions, i.e.,

[TABLE]

where $\bm{a}_{n}\in\mathbb{R}^{K}$. To facilitate compact notation, the row vectors $\bm{a}_{n}$ are stacked into a matrix $\bm{A}\in\mathbb{R}^{N\times K}$. The minimization in (7) can then be written as the finite-dimensional convex optimization problem

[TABLE]

Since the degrees of freedom in (15) are reduced to $NK$, the constraints on each density $p_{\bm{a}_{n}}$ can only be satisfied at the $K$ grid points. Accordingly, the optimality conditions for problem (15) are obtained from the general conditions in Theorem 1 or Corollary 1 by evaluating them on the grid $\bm{\omega}$. In particular, the optimality condition in Corollary 1 becomes

[TABLE]

for all $n=1,\ldots,N$. In analogy to $f_{n}^{-1}(\bm{q}_{[n]},c_{n})$, the inverse function $f_{n}^{-1}(\bm{A}_{[n]},c_{n})$ in (16) is defined as

[TABLE]

where the infimum is taken element-wise and $f_{n}(\bm{A})$ is shorthand for the vector $f_{n}(\bm{a}_{1},\ldots,\bm{a}_{N})$.

V-C Discrete Optimality Gap

Given a feasible vector of densities $\bm{p}$, the integrals in (10) and (11) can in principle be evaluated using standard numerical integration techniques. However, this might be problematic in practice. First, depending on the dimensions of $\Omega$ and the chosen technique, numerical integration can be computationally costly. Second, when solving the discrete problem (15), the set of feasible solutions is restricted to densities of the form $p_{\bm{a}_{n}}$ in (14), which will typically not solve the continuous problem (8) exactly. Hence, even if the matrix $\bm{A}$ satisfies the discrete optimality conditions in (16), the optimality gap can be nonzero. This property is undesirable since in order to provide a useful termination criterion the optimality gap should converge to zero. That is, it should reflect how accurately a candidate solution solves the discrete problem (15), instead of the continuous version (8).

In order to avoid these two problems, the following discrete versions of the residuals in (10) and (11) are proposed:

[TABLE]

The expressions in (18) and (19) are obtained from (10) and (11) by approximating all continuous functions of $\omega$, i.e., $p^{\prime}_{n}(\omega),p^{\prime\prime}_{n}(\omega)$, and $f_{n}(\bm{p}(\omega))$, via linear combinations of the basis functions $\psi_{1},\ldots,\psi_{K}$ and assuming all $p_{n}$ to be of the form (14). Hence, (18) and (19) are approximations of the true residuals whose accuracy depends on how well the chosen basis functions represent the true functions. Being an inner product of two vectors, (18) and (19) can be evaluated efficiently, even for large grid sizes $K$. Moreover, by construction, the discrete residuals become zero whenever a density satisfies the discrete optimality conditions in (16).

V-D Block Coordinate Descent

The first proposed algorithm for numerically approximating the optimal densities is detailed below. It is based on the observation that if $N-1$ optimal densities are given, then the remaining one can be determined by a simple line search over the scalar $c_{n}$.

Algorithm 1 implements a block coordinate descent (BCD) [26, 27] in the vectors $\bm{a}_{1},\ldots,\bm{a}_{N}$ with a custom coordinate selection rule. It is composed of two main steps, namely:

the coordinate selection step: calculate the residuals $\tilde{e}_{1},\ldots,\tilde{e}_{N}$ and pick the index $n^{*}$ corresponding to the largest residual; 2. 2.

the coordinate descent step: determine the optimal vector $\bm{a}_{n^{*}}$ while keeping the remaining rows of $\bm{A}$ fixed.

The idea behind the coordinate selection step is to optimize the density that has the largest residual and hence the largest contribution to the optimality gap. Although this selection scheme is a heuristic, it provides a well-motivated, systematic rule to choose the next coordinate. See the example in Section VI for a comparison of the proposed rule to a cyclic and a random selection rule.

In the coordinate descend step of Algorithm 1, the weight vector $\bm{a}_{n}$ is chosen such that $p_{\bm{a}_{n}}$ satisfies the discrete optimality condition in (16) and, therefore, is a coordinate-wise minimum of (15). The central task in the coordinate descend step is the evaluation of $f_{n}^{-1}(\bm{A}_{[n]},c)$. In case the inverse of $f_{n}$ has an analytical from, the corresponding explicit function can simply be substituted for $f_{n}^{-1}$. An example for this case is given in Section VI-A. Otherwise, $f_{n}^{-1}(\bm{A}_{[n]},c)$ needs to be evaluated numerically by solving (17) for $\bm{a}_{n}$. Owing to property (13) of the basis functions, the elements of the vector $\bm{a}_{n}$ in (17) can be determined by solving

[TABLE]

for each $k=1,\ldots,K$ individually. This decoupling makes it possible to evaluate $f_{n}^{-1}(\bm{A}_{[n]},c)$ in a highly parallel manner, using up to $K$ compute cores simultaneously. In addition, the memory requirements per core are minimal since apart from $f_{n}$ and $c$ only the vector $(a_{1,k},\ldots,a_{N,k})$ needs to be stored, which is of dimension $N$, irrespective of how $K$ is chosen.

Since $f_{n}$ is monotonically nondecreasing for all $n$, (20) can be solved via standard one-dimensional root-finding techniques. In order to obtain $\bm{a}_{n}$, the vector $f_{n}^{-1}(\bm{A}_{[n]},c)$ is then projected onto the density band $\mathcal{P}^{=}_{n}$. Consequently, no root-finding needs to be performed if the root is guaranteed to lie outside the feasible interval $[p^{\prime}_{n}(\omega_{k}),p^{\prime\prime}_{n}(\omega_{k})]$. Whether or not this is the case can be determined by simply evaluating $f_{n}$ at the endpoints of the interval. From

[TABLE]

it follows that $a_{n,k}<p^{\prime}_{n}(\omega_{k})$ so that, after the projection onto the feasible band, $a_{n,k}=p^{\prime}_{n}(\omega_{k})$. Analogously, it follows from

[TABLE]

that $a_{n,k}=p^{\prime\prime}_{n}(\omega_{k})$. These additional checks are consistent with the definition of the generalized inverse in (2) and simplify the evaluation of $f_{n}^{-1}(\bm{A}_{[n]},c)$ significantly.

Algorithm 1 terminates if the optimality gap becomes smaller than the tolerance $\varepsilon$. It should be emphazised again that $\tilde{\bm{e}}$ in line 8 corresponds to the optimality gap of the discrete problem (15) so that $\langle\tilde{\bm{e}},\bm{1}\rangle\leq\varepsilon$ does not guarantee $\langle\bm{e},\bm{1}\rangle\leq\varepsilon$. However, for reasonably small grid sizes both values are sufficiently close. Moreover, in case of strict optimality requirements, the equations in Theorem 2 can be used to obtain a tighter bound on the true optimality gap.

Algorithm 1 is a straightforward attempt at solving the system of optimality conditions in Theorem 1. By construction, its limit points satisfy (16) and, hence, are global minimizers of (15). However, for general convex functions $f$, its convergence cannot be guaranteed since the block coordinate descent strategy implemented in Algorithm 1 is not guaranteed to converge for objective functions that are not strictly convex [28]. This problem is intrinsic to coordinate-wise minimization techniques and exists irrespective of the coordinate selection rule and the termination criterion. A proximal algorithm that addresses these shortcomings at the cost of a reduced efficiency is detailed in the next section.

V-E Proximal BCD With Guaranteed Convergence

The algorithm presented in this section solves the system of optimality conditions by means of a proximal iteration instead of a regular fixed-point iteration. It extends Algorithm 1 to cover cases where it might otherwise fail to converge.

Proximal algorithms are well studied and offer a reliable tool for iteratively solving optimization problems that are not strictly convex. The underlying idea is to construct a sequence of strictly convex objective functions whose minimizers converge to the minimizer of the original problem. To this end, the original objective function is augmented with a strictly convex term such that the solution of the augmented problem is guaranteed to be unique. More precisely, in each iteration, a proximal algorithm seeks to minimize a weighted sum of the objective function and a term that penalizes some distance between the current and the previous iterate. The additional distance term ensures that the problem is strictly convex and automatically vanishes when a minimum of the orginal objective function is approached. A comprehensive introduction to proximal algorithms can be found in [29].

For the function $I_{f}$ in (5), a proximal operator based on the $\mathcal{L}_{\mu}^{2}$-norm can be defined as

[TABLE]

where $h_{1},\ldots,h_{N}\colon\Omega\to\mathbb{R}$ are $\mu$-integrable functions. The proximal algorithm iteratively approximates the solution of the fixed-point equation

[TABLE]

In order to implement an iteration of the proximal algorithm, the inner problem in (21) needs to be solved. By inspection, this problem is equivalent to the original problem (8) with $f$ replaced by

[TABLE]

which is strictly convex in $\bm{x}$ for every convex function $f$ and admits the partial subderivatives

[TABLE]

The proximal block coordinate descent algorithm is denoted Algorithm 2 and is specified below.

Proposition 1** (Convergence of Algorithm 2)**

Let the basis functions in (12) be chosen such that $\bm{\mu}>0$. For all convex functions $f$ and all bands $\mathcal{P}^{=}$ that satisfy

[TABLE]

Algorithm 2 is guaranteed to converge to a global minimizer of (15).

A proof for Proposition 1 can be found in Appendix D. The two conditions on the density bands are vectorized versions of the integral inequalities in (4) and ensure that in each iteration a feasible density of the form (12) can be constructed.

The price for the guaranteed convergence of the proximal version of the block coordinate descent algorithm is reduced efficiency. Since Algorithm 2 repeatedly calls Algorithm 1 to solve the minimization in (21), it requires as least as many iterations as Algorithm 1, assuming that the latter converges. The total number of coordinate descent steps of Algorithm 2 is roughly given by the product of the number of iterations required by both algorithms. Moreover, in some cases an analytic expression for the inverse of $f_{n}$ exists, but not for the inverse of $\tilde{f}_{n}$. In general, the use of Algorithm 2 is recommended only when Algorithm 1 indeed fails to converge.

V-F Remarks

There are several options to improve the performance of the presented algorithms. First, if evaluating the residuals is significantly more expensive than the coordinate descent step, the coordinate selection rule can be modified to only evaluate the residuals every $M$th $(M>1)$ iteration or a simple cyclic or random selection rule can be used [30, 31]. Second, in Algorithm 2, the weight of the penalty term can be reduced in order to reduce the number of outer iterations. In practice, this means trading off speed of convergence for numerical stability.

Apart from the two algorithms presented in this section, the finite-dimensional problem (15) can, in principle, be solved by many off-the-shelf convex optimization programs. The proposed algorithms are often preferable over generic solvers for the following reasons:

First, the approach presented in this paper often allows an analytical general form for the optimal densities to be obtained. The numerical part of the solution process then reduces to determining the scalars $c_{1},\ldots,c_{N}$, which usually is a much simpler problem. Moreover, having a parametric expression for the least favorable densities can facilitate the derivation of analytical bounds or approximations to the exact solution. An example where the system of optimality conditions leads to a useful lower bound is given in the next section.

Second, the accuracy of generic convex optimization algorithms can be hard to control for certain density bands. If, for example, in some region of the sample space the width of the band is of the same order of magnitude as the tolerance of the solver, all points within the band become equivalent. While this effect is negligible in terms of the shape of the optimal densities themselves, it can be critical when calculating functions that involve products or ratios of optimal densities. Although the amount of numerical noise in the solution can be reduced by applying suitable variable transformations and carefully tuning the absolute and relative tolerances of the solvers, the general issue of having to deal with, potentially, very badly scaled problems is an important consideration.

The presented algorithms avoid scaling problems by reducing the entire optimization process to repeated searches for the root of a nondecreasing real-valued function. In each iteration, both $c_{n}$ and, if necessary, $\bm{a}_{n}$ are determined via one-dimensional root-finding. The latter is a standard problem in numerical mathematics and can be solved fast and reliably—see, for example, [32, Chapter 8.3] for an overview of suitable algorithms. In cases where $f_{n}$ can only be evaluated with (numerical) noise, stochastic root-finding methods, which handle the additional uncertainty in a systematic manner [33], can be applied. Owing to the explicit projection on $\mathcal{P}^{=}_{n}$, the optimal densities are exact in regions of $\Omega$ where the band constraints are active.

The third advantage of the proposed algorithms is that they can be significantly faster than generic solvers. The most substantial performance gains can be achieved with Algorithm 1 if an analytic expression for the inverse of $f_{n}$ exists; see the example in the next section. However, even if $f_{n}$ has to be inverted numerically, this can be done in a highly parallel manner that facilitates the use of high performance graphics processing units and distributed computing systems. Note that neither of the presented algorithms needs to perform matrix operation (multiplication, inversion) which makes them highly scalable in terms of the number of densities $N$ and the number of grid points $K$. This is illustrated with an example in the next section.

In summary, the proposed approach to the minimization of convex functionals offers useful theoretical insights and at the same time provides a fast and reliable way to obtain accurate numerical results.

VI Examples

In this section, the usefulness of the results presented in the previous sections is demonstrated by means of two examples. First, the convex functional is chosen as a weighted sum of Kullback–Leibler divergences. This example shows how analytical results can be obtained by means of Theorem 1 and it demonstrates the increased efficiency and accuracy of the proposed algorithms in comparison to a state-of-the-art generic solver. In the second example, the least favorable distributions for a binary decision making problem with an observation dependent cost function are derived. This example illustrates how the concept of minimizing convex functionals is applicable beyond its traditional context of statistical distance measures.

In order to simplify the presentation, the sample space is chosen to be the real line, i.e., $\Omega=\mathbb{R}$. For the finite-dimensional representation of the densities, a regular grid with step size $\omega_{k}-\omega_{k-1}=\Delta\omega$ is used in combination with a linear interpolation scheme, i.e.,

[TABLE]

The masses of the basis functions calculate to $\mu_{k}=\Delta\omega$ for all $k=1,\ldots,K$ and a simple bisection algorithm was used to perform the root finding.

VI-A Weighted Sum of Kullback–Leibler Divergences

Sums of $f$-divergences, and Kullback–Leibler divergences in particular, have been shown to have applications in a variety of fields, including minimax robust statistics [18], geoscience [34] and biology [35]. In this example, $I_{f}$ is a weighted sum of Kullback–Leibler divergences with respect to a common reference distribution $P_{N}$, i.e.,

[TABLE]

where $D_{\text{KL}}(\cdot\|\cdot)$ denotes the Kullback-Leibler divergence and $\alpha_{1},\ldots,\alpha_{N-1}$ are convex combination weights, i.e., they satisfy

[TABLE]

The corresponding function $f$ is given by

[TABLE]

Its partial derivatives are

[TABLE]

for $n=1,\ldots,N-1$ and

[TABLE]

The inverse functions are obtained by solving $f_{n}=c_{n}$ for $x_{n}$ and are given by

[TABLE]

for $n=1,\ldots,N-1$ and

[TABLE]

Here $b_{1},\ldots,b_{N-1}\in\mathbb{R}$ and $b_{N}>0$ are introduced for the sake of a more compact notation. From Corollary 1 it follows that the optimal densities are of the form

[TABLE]

for $n=1,\ldots,N-1$ and

[TABLE]

That is, $q_{1},\ldots,q_{N-1}$ are the projections of $q_{N}$ onto the bands $\mathcal{P}^{=}_{1},\ldots,\mathcal{P}^{=}_{N-1}$, respectively, and $q_{N}$ is the projection of the weighted geometric mean of $q_{1},\ldots,q_{N-1}$ onto the band $\mathcal{P}^{=}_{N}$.

Before presenting numerical results, it is shown that the above expressions for the optimal densities can be used to derive a tight lower bound on (23) for the special case where $p_{1},\ldots,p_{N-1}$ are given and the optimization is performed only over $p_{N}$. Problems of this kind are, for example, considered in [18] in order to derive lower bounds on the minimax risk of a decision making procedure. Making use of (25), it can be shown that

[TABLE]

where

[TABLE]

is a generalized version of the Bhattacharyya coefficient [36]. To the best of our knowledge, this bound has not been stated in the literature so far.

For the numerical minimization of (23) consider $N=3$ densities and uncertainty bands defined by scaling and shifting Gaussian densities according to

[TABLE]

where $p_{\mathcal{N}}(m,\sigma^{2})$ denotes the density function of a Gaussian distribution with mean $m$ and variance $\sigma^{2}$.

Three triplets of optimal densities, for different weights $\alpha_{1}$ and $\alpha_{2}$, are depicted in Fig. 1. As can be seen, $q_{1}$ and $q_{2}$ are independent of the weights in this particular example, but $q_{3}$ changes significantly and different combinations push it either towards $q_{1}$ or $q_{2}$.

The densities in Fig. 1 were calculated using Algorithm 1 on the interval $[-5,5]$ with step size $\Delta\omega=0.01$ and initial densities $q_{1}=p_{\mathcal{N}}(-0.5,1)$, $q_{2}=p_{\mathcal{N}}(0.5,1)$, $q_{3}=p_{\mathcal{N}}(0,1)$. The tolerance for the optimality gap was set to $\varepsilon=10^{-7}$. The number of iterations required to reach convergence is shown in Table I for different coordinate selection rules. As can be seen, a selection based on the largest residual is preferable for asymmetric weights $\alpha_{1},\alpha_{2}$, while for symmetric weights the cyclic rule performs best. As expected [31], random coordinate selection, i.e., drawing the next coordinate from the set $\{1,\ldots,N\}\setminus\{n^{*}\}$, performs uniformly worst. Although these results are not conclusive, they indicate that a selection based on the largest residual is less sensitive to the shape of the objective function compared to a cyclic rule.

In this example, the inverse functions $f_{n}^{-1}$ do not need to be evaluated numerically and Algorithm 1 is highly efficient. Our Matlab© implementation terminated in well under a second on a regular desktop computer.222 All simulations were performed on an Intel® Core™ i5-760@2.80GHz using Matlab© 2016a.

For comparison, the minimization problem was also solved using version 2.0.4 of the ECOS solver [37], which is a state-of-the-art software package for solving conic optimization problems and is written in C. The reason for choosing ECOS over other options is that it is one of the few high performance solvers that support the exponential cone and, hence, logarithmic objective functions. The average run-times of Algorithm 1, using the largest residual coordinate selection rule, and the ECOS solver are given in Table II. The results were obtained on the same machine with absolute tolerances set to $10^{-7}$ for all algorithms. The ECOS solver was called via its CVX interface [38], but only the time spent in the C routine was used for the benchmark. As can be seen, Algorithm 1 is consistently faster than the ECOS solver, despite the additional handicap of being written in an interpreted language. Especially for large numbers of grid points $K$, the advantage of the proposed approach becomes obvious.

Table I and Table II also include the number of iterations and run-times of Algorithm 2. It can be shown that for $f$ of the form (24) the inverses of the partial derivatives in (22) calculate to

[TABLE]

for $n\in\{1,\ldots,N\}$ and

[TABLE]

Here, $W$ denotes the Lambert W function. Table I shows the number of inner iterations for Algorithm 2, i.e., the total number of coordinate descent steps. The results were obtained by calling Algorithm 1 with the specified selection rule in each iteration of Algorithm 2. The number of outer iterations, i.e., the number of times Algorithm 1 was called from within Algorithm 2, is (approximately) independent of the coordinate selection rule. In this example, Algorithm 2 required roughly $40$, $50$, and $110$ outer iterations for $\alpha_{1}=0.5$, $\alpha_{1}=0.7$, and $\alpha_{1}=0.1$, respectively. Owing to the increased number of coordinate descent steps, and the computationally costly evaluation of the Lambert W function, Algorithm 2 is also significantly slower than Algorithm 1. However, for medium to large problem sizes its performance is comparable to, or even better than, that of the ECOS solver.

These results arise from the chosen example and it is noted that raw execution time is not a reliable performance metric. Nevertheless, the example is non-trivial and the fact that the proposed algorithms are able to outperform an optimized software package is indicative of a good performance in general.

Another important aspect of the proposed algorithms is that they lead to results with high accuracy. This can be seen by inspection of the ratios of the optimal densities, which are of particular interest in detection problems, where they determine the optimal test statistic. For $\alpha_{1}=0.7$ and $\alpha_{2}=0.3$, the log-likelihood ratios

[TABLE]

as calculated by the proposed algorithms and the ECOS solver, are depicted in Fig. 2. Since the objective function in (23) is strictly convex, Algorithm 1 and Algorithm 2 converge to the same solution and only the result of Algorithm 1 is depicted in Fig. 2. Note that the interval of the sample space is increased to $[-10,10]$. It can clearly be seen how the ECOS solver produces artifacts at the tails of the densities, where their values fall below the tolerance of $10^{-7}$. Algorithm 1, in contrast, correctly identifies the tails as regions where the band constraints are active and is, therefore, able to calculate the likelihood ratios exactly. If at all possible, obtaining results of comparable quality with generic solvers requires careful parameter tuning and tolerances close to the machine precision.

VI-B Minimax Detection with Observation Dependent Cost

The purpose of the second example is to show that the proposed algorithms handle cases where $f$ directly depends on $\omega$ and that this dependence widens the scope of problems which can be solved within the presented framework.

Consider a binary decision making problem with a cost function

[TABLE]

where $E_{P}$ denotes expectation with respect to a distribution $P$, $\delta\colon\Omega\to[0,1]$ denotes a decision rule, and $r_{0},r_{1}\colon\Omega\to[0,\infty)$ denote observation dependent costs for each decision. Cost functions of this form occur, for example, in detection problems, where the cost for an incorrect decision depends on the true state of the system. In a collision avoidance system for vehicles, for instance, the cost for not detecting an obstacle becomes higher when the obstacle is close to the vehicle or the vehicle is moving at a high speed.

The optimal decision rule for the cost function (27) can be shown to be given by

[TABLE]

Without loss of generality, it is assumed that $\kappa=1$. Using this decision rule, the expected cost is given by

[TABLE]

It is assumed that $p_{1}$ and $p_{2}$ are subject to uncertainties of the density band type. In order to design a minimax detector, i.e., a detector that minimizes the worst-case cost, (28) needs to be maximized with respect to the densities $p_{1}$ and $p_{2}$. This problem is of the form (7) with

[TABLE]

and

[TABLE]

For illustration purposes, the cost functions are chosen as

[TABLE]

Their graphs are shown in Fig. 3. The same density bands as in (26) are used to constrain $p_{1}$ and $p_{2}$. The grid for the discrete representation is constructed on $[-5,5]$ with grid size $\Delta\omega=0.01$. Since $f$ in (29) is not strictly convex, Algorithm 2 was used with initial densities $q_{1}=p_{\mathcal{N}}(-0.5,1)$, $q_{2}=p_{\mathcal{N}}(0.5,1)$. It reached the required tolerance of $\varepsilon=10^{-7}$ after $57$ iterations. The resulting least favorable densities are depicted in Fig. 4. Their effect on the cost function can be seen in Fig. 5, where $r^{*}(p_{0},p_{1})$ is plotted for the least favorable and two Gaussian densities. Interestingly, the shape of the cost function $r^{*}$ is preserved, despite the “ragged” shape of the least favorable densities. Again, the proposed algorithm is able to accurately identify abrupt changes as well as smooth variations in the optimal densities.

VII Conclusion

The problem of minimizing convex functionals of probability distributions under density band constraints was analyzed. A system of sufficient and necessary first-order optimality conditions was derived as well as a bound on the gap between the exact minimum and the value of the objective function at a candidate solution. The optimality conditions were used to characterize global minimizers of the constrained optimization problem as solutions of a nonlinear fixed-point equation. Two algorithms were proposed that iteratively solve this equation by following a block coordinate descent strategy until the optimality gap falls below a given tolerance. While the first algorithm proved to be efficient in practice, it is not guaranteed to converge for objective functions that are not strictly convex. This problem was overcome by introducing an additional outer proximal iteration. The modified algorithm was then shown to admit guaranteed convergence properties for all band constraints and all convex objective functions, at the cost of a reduced efficiency. Two examples were given to illustrate how the optimality conditions can be used in practice to derive analytical results and to demonstrate the high efficiency and accuracy of the proposed numerical algorithms.

Appendix A Proof of Theorem 1

Theorem 1 is proven by showing that if the densities $(q_{1},\ldots,q_{N})$ satisfy the conditions in Theorem 1, they also satisfy the Karush–Kuhn–Tucker (KKT) conditions of (8). Since $I_{f}$ is convex in $(p_{1},\ldots,p_{N})$ by assumption and the constraints in (8) are linear, the KKT conditions are necessary and sufficient for $(q_{1},\ldots,q_{N})$ to be a global minimizer [39].

The proof makes use of some basic results of infinite-dimensional optimization theory, more precisely, the theory of Lagrange multipliers on Banach spaces. An introduction to the topic is beyond the scope of this paper. A comprehensive treatment can be found, for example, in [40]. An elegant standalone proof for the method of Lagrange multipliers and the sufficiency of the KKT conditions is given in [41] and [42], respectively. In brief, the method of Lagrange multipliers can be applied to convex functions on Banach spaces by using Fréchet subderivatives instead of subgradients and elements of the dual space instead of scalar- or vector-valued Lagrange multipliers.

Let

[TABLE]

denote the space of all essentially bounded functions on $\Omega$, which is the dual space of $\mathcal{L}_{\mu}^{1}$. The Lagrange function $L\colon(\mathcal{L}_{\mu}^{1}\times\mathcal{L}_{\mu}^{\infty}\times\mathcal{L}_{\mu}^{\infty}\times\mathbb{R})^{N}\to\mathbb{R}$ of (8) is given by

[TABLE]

where

[TABLE]

and $u_{n},v_{n}\in\mathcal{L}_{\mu}^{\infty}$, $c_{n}\in\mathbb{R}$ denote the Lagrange multipliers corresponding to the constraints $p_{n}-p^{\prime\prime}\leq 0$, $p_{n}-p^{\prime}\geq 0$ and $\int_{\Omega}p_{n}\,\mathrm{d}\mu=1$, respectively. The dual problem is given by

[TABLE]

The partial Fréchet-subdifferential of $L$ with respect to $p_{n}$ can be shown to be

[TABLE]

The KKT conditions for the optimal densities require that for all $n=1,\ldots,N$

[TABLE]

Let all $q_{n}$ and $c_{n}$ be chosen such that they comply with the conditions in Theorem 1. By construction, this implies that $q_{n}$ satisfies the primal feasibility constraints. Since $f$ is convex and $\mu$ is continuous, it follows from Rademacher’s theorem [43, Chapter 9.J.] that the partial subderivatives of $f$ are finite $\mu$-almost everywhere, i.e., $f_{n}\in\mathcal{L}^{\infty}$ for all $n$. Without violating dual feasibility, the functions $v_{n}$ and $u_{n}$ can be chosen as

[TABLE]

so that

[TABLE]

Inserting (42) back into the stationarity condition yields

[TABLE]

for all $n=1,\ldots,N$. The last step in the proof is to show that these choices for $q_{n}$, $u_{n}$, and $v_{n}$ also satisfy the complementary slackness constraints, i.e.,

[TABLE]

for all $\omega\in\Omega$. By construction of $u_{n}$ and $v_{n}$, $u_{n}(\omega)>0$ implies $f_{n}(\bm{q}(\omega))<c_{n}$, which in turn implies $q_{n}=p^{\prime\prime}_{n}$. Analogously, $v_{n}(\omega)>0$ implies $f_{n}(\bm{q}(\omega))>c_{n}$ and in turn $q_{n}(\omega)=p^{\prime}_{n}(\omega)$. $\hfill\square$

Appendix B Proof of Corollary 1

Corollary 1 is a consequence of the fact that $f_{n}$, being a subderivative of a convex function, is nondecreasing in $x_{n}$. For the three cases in Theorem 1, it then follows that

[TABLE]

In words, the three equations state that:

• •

$q_{n}$ equals its upper bound $p_{n}^{\prime\prime}$, if $f_{n}^{-1}(\bm{q}_{[n]},c_{n})$ is larger than $p_{n}^{\prime\prime}$;

• •

$q_{n}$ equals its lower bound $p_{n}^{\prime}$, if $f_{n}^{-1}(\bm{q}_{[n]},c_{n})$ is smaller than $p_{n}^{\prime}$;

• •

$q_{n}$ equals $f_{n}^{-1}(\bm{q}_{[n]},c_{n})$ otherwise.

The expression for $q_{n}$ given in Corollary 1 is merely a more compact way of writing this case-by-case definition. $\hfill\square$

Appendix C Proof of Corollary 2

The bound on the optimality gap in Corollary 2 can be obtained from the Lagrange dual in (35). By construction, it holds that

[TABLE]

for all feasible Lagrange multipliers $\bm{u}$, $\bm{v}$, and $\bm{c}$. This minimization is hard to solve in general, but a lower bound on $I_{f}(\bm{q})$ can be obtained as follows: instead of optimizing over $\bm{p}$, the Lagrange multipliers $\bm{u}$, $\bm{v}$ and $\bm{c}$ can be chosen such that a given $\bm{p}$ satisfies the stationarity conditions in (36) and, consequently, solves the minimization in (43). Let $\bm{u}^{*}$ and $\bm{v}^{*}$ denote multipliers that satisfy (40) and (41), respectively. It is shown in Appendix A that this choice fulfills the stationarity conditions. It hence holds that

[TABLE]

$L(\bm{p},\bm{u}^{*},\bm{v}^{*},\bm{c})$ can be shown to evaluate to

[TABLE]

with $\bm{e}$ defined in Theorem 2. Note that the sum over $c_{n}$ in (31) cancels with the $p_{n}c_{n}$ terms in (32) since all $p_{n}$ are assumed to be valid densities. From (44), the bound on the optimality gap follows:

[TABLE]

$\hfill\square$

Appendix D Proof of Convergence of Algorithm 2

The convergence of proximal iterations is a well-established result in the convex optimization literature. It follows directly from the contractive property of the proximal operator and can be applied to Algorithm 2 in a straightforward manner. For a selection of convergence proofs, see, for example, [29, Chapter 2.3], [44, 45, 46, 47], and the references therein.

In order to prove that Algorithm 2 converges, it hence suffices to show that Algorithm 1 indeed solves the inner minimization in (21). This is guaranteed if, first, the function $f$ is strictly convex and, second, if, in every iteration, the equation in line 6 of Algorithm 1 has a solution. The first condition ensures that the algorithm is able to find a block-coordinate-wise minimum in each iteration; the second condition ensures that the block-coordinate-wise descent indeed converges to a global minimum.

The first condition is fulfilled by construction. Since $\tilde{f}$ is chosen such that it is strictly convex in $(x_{1},\ldots,x_{n})$, $I_{f}(p_{\bm{a}_{1}},\ldots,p_{\bm{a}_{N}})$ is strictly convex in $(\bm{a}_{1},\ldots,\bm{a}_{N})$.

The second condition can be shown to be fulfilled for arbitrary convex functions $f$ of the form (6). By definition of the inverse function in (2), it holds that for every $n=1,\ldots,N$

[TABLE]

for all $\bm{x}\in[0,\infty)^{N}$. Therefore,

[TABLE]

and

[TABLE]

where the last inequalities hold by assumption. Consequently, some $c\in\mathbb{R}$ is guaranteed to exist that solves the equation in line 6 of Algorithm 1. $\hfill\square$

Acknowledgment

The authors would like to thank the anonymous reviewers for their helpful and constructive comments that greatly improved the quality of the paper.