Adaptive, anisotropic and hierarchical cones of discrete convex functions

Abstract

We introduce a new class of adaptive methods for optimization problems posed on the cone of convex functions. Among the various mathematical problems which possess such a formulation, the Monopolist problem (Rochet and Choné, Econometrica 66:783–826, 1998; Ekeland and Moreno-Bromberg, Numer Math 115:45–69, 2010) arising in economics is our main motivation. Consider a two dimensional domain $\Omega$, sampled on a grid $X$ of N points. We show that the cone ${{\mathrm{Conv}}}(X)$ of restrictions to $X$ of convex functions on $\Omega$ is typically characterized by $\approx N^2$ linear inequalities; a direct computational use of this description therefore has a prohibitive complexity. We thus introduce a hierarchy of sub-cones ${{\mathrm{Conv}}}(\mathcal{V})$ of ${{\mathrm{Conv}}}(X)$, associated to stencils $\mathcal{V}$ which can be adaptively, locally, and anisotropically refined. We show, using the arithmetic structure of the grid, that the trace $U_{|X}$ of any convex function U on $\Omega$ is contained in a cone ${{\mathrm{Conv}}}(\mathcal{V})$ defined by only $\mathcal{O}( N \ln ^2 N)$ linear constraints, in average over grid orientations. Numerical experiments for the Monopolist problem, based on adaptive stencil refinement strategies, show that the proposed method offers an unrivaled accuracy/complexity trade-off in comparison with existing methods. We also obtain, as a side product of our theory, a new average complexity result on edge flipping based mesh generation.

Appendices

Appendix A: Convergence of the discretization

For completeness, and at the request of a referee, we reproduce (with minor adaptations) arguments of [7] which establish the convergence of our discrete numerical approximations of the principal agent problem towards the continuous solution, as the grid scale tends to 0.

Let $\Omega \subseteq {\mathbb R}^2$ be an open bounded convex polygonal domain, and let $(\mathcal{T}_h)_{h >0}$ be a family of regular triangulations of $\Omega$: there exists a constant C such that, for all $h >0$, $T \in \mathcal{T}_h$, one has

$$\begin{aligned} \text {Correct scaling: } {{\mathrm{diam}}}(T)&\le h.&\text {Shape regularity: }{{\mathrm{diam}}}(T) \le&C \rho (T). \end{aligned}$$

We denoted by ${{\mathrm{diam}}}(T)$ the largest distance between two points of a triangle T, and by $\rho (T)$ its inner circle radius. Let $X_h$ be the collection of vertices of $\mathcal{T}_h$, and let ${{\mathrm{\mathrm{I}}}}_h : \mathcal{F}(X_h) \rightarrow C^0(\Omega )$ be the piecewise linear interpolation operator associated to $\mathcal{T}_h$.

The classical principal agent problem, and its finite element discretization, are as follows:

$$\begin{aligned} E_0&:= \inf \left\{ \mathcal{E}(u); \, u \!\in \! {{\mathrm{Conv}}}(\Omega ), \, u \!\ge \! 0\right\} , \text { with } \mathcal{E}(u)\!:=\! \int _\Omega \left( \frac{1}{2} \Vert \nabla u\Vert ^2 -\langle \nabla u, z\rangle \!+\! u\right) .\\ E_h&:= \inf \{ \mathcal{E}({{\mathrm{\mathrm{I}}}}_h v); \, v \in {{\mathrm{Conv}}}(\Omega ), \, v \ge 0\} = \inf \{ \mathcal{E}({{\mathrm{\mathrm{I}}}}_h v); \, v \in {{\mathrm{Conv}}}(X_h), \, v \ge 0\}. \end{aligned}$$

As proposed in [7], the triangulations $\mathcal{T}_h$ are fixed, non-adaptive, and solely used for constructing the discrete energies $u \mapsto \mathcal{E}({{\mathrm{\mathrm{I}}}}_h u)$. In particular ${{\mathrm{\mathrm{I}}}}_h u$ need not be convex; instead we use the rigidity result Proposition 8.1 to establish the convexity of the weak limit of the minimizers to the discretized problems. Our numerical experiments all involve rectangular domains $\Omega$ and discretization points $X_h$ distributed on a cartesian grid, for which $\mathcal{T}_h$ is the standard structured mesh built of half squares.

Proposition 8.1

For all $h>0$ let $v_h \in {{\mathrm{Conv}}}(\Omega )$ and $u_h := {{\mathrm{\mathrm{I}}}}_h v_h$. Assume that $u_h \rightarrow u_0$ as $h \rightarrow 0$ in the sense of distributions. Then $u_0$ is convex, and both $u_h$ and $v_h$ converge uniformly to $u_0$ on compact subsets of $\Omega$.

Proof

We fix a radius $r>0$, and limit our attention to grid scales $0<h\le r$.

Non-uniform bound on $u_h$. There exists M, independent of $0<h \le r$, such that for any ball $B(x,r) \subseteq \Omega$, one does not have $u_h \ge M$ (resp. $u_h \le -M$) uniformly on B(x, r). Indeed otherwise there exists a sequence $h_n \rightarrow 0$, $M_n\rightarrow \infty$, $B(x_n,r) \subseteq \Omega$, such that $u_n :=u_{h_n} \ge M_n$ (resp. $\le -M_n$) on $B(x_n,r)$. By compactness we can assume that $x_n \rightarrow x_* \in \Omega$. Considering a smooth bump function $\varphi$ supported on $B(x_*,r)$ we obtain $\int u_n \varphi \rightarrow + \infty$ (resp. $-\infty$) which contradicts the distributional convergence of $u_h$ as $h\rightarrow 0$.

Non-uniform bound on $v_h$. If for contradiction $|v_h| \ge M$ uniformly on a ball $B(z,2 r)\subseteq \Omega$, then by continuity $v_h \ge M$ (resp. $v_h \le -M$) on B(z, 2r). Thus by interpolation $u_h = {{\mathrm{\mathrm{I}}}}_h v_h \ge M$ (resp. $u_h \le -M$) on B(z, r), which was excluded.

Upper bound on $v_h$. Consider a ball $B(z,4r) \subseteq \Omega$, a subgradient $g \in \partial _z v_h$, and note that by convexity $v_h(w) \ge v_h(z)+\langle w-z,g\rangle \ge v_h(z)$ for all $w \in B(z+2r g / \Vert g\Vert , 2r) \subseteq B(z,4r) \subseteq \Omega$. From the non-uniform bound on $v_h$, we conclude that $v_h(z) \le M$.

Lower bound on $v_h$. Consider a ball $B(x,10 r) \subseteq \Omega$. By the non-uniform bound on $v_h$, there exists $y \in B(x,2r)$ such that $v_h(y) \ge -M$. Denoting $z := y+2(y-x)$ one has by convexity $v_h(z) \ge v_h(y)+2 (v_h(y)-v_h(x))$. In addition $B(z,4 r) \subseteq B(y,8 r) \subseteq B(x,10 r) \subseteq \Omega$, thus $M \ge v_h(z) \ge v_h(y)+2 (v_h(y)-v_h(x)) \ge -3 M - 2 v_h(x)$, so that $v_h(x) \ge -2 M$.

Lipschitz bound on $v_h$. Consider a ball $B(x,10r) \subseteq \Omega$, and a sub-gradient element $g \in \partial _x v_h$. Let $y := x+r g/\Vert g\Vert$. Then by convexity $r \Vert g\Vert = \langle g, y-x\rangle \le v_h(y)-v_h(x) \le M-(-2M)=3M$, thus $\Vert g\Vert \le C_L := 3M/r$. In particular, if $B(x,11 r) \subseteq \Omega$ and $x \in T \in \mathcal{T}_h$, then $v_h$ is $C_L$-Lipschitz on T, so that by interpolation $|u_h(x) - v_h(x)| \le C_L h$.

Conclusion. Consider the convex compact $K := \{x \in \Omega ; \, B(x,11r)\}$. The maps $(v_h)_{0<h \le r}$ are uniformly bounded and Lipschitz on K, hence they form by Arzela–Ascoli’s theorem a pre-compact family: from any sub-sequence, one can extract a sub-sub-sequence which converges uniformly on K to a limit $v_*$ as $h \rightarrow 0$. Since $\Vert u_h-v_h\Vert _{L^\infty (K)} \rightarrow 0$, the sub-sub-sequence limit can only be $v_* = u_0$. This implies that $u_0$ is convex and that $v_h \rightarrow u_0$ uniformly on K as $h \rightarrow 0$. Thus also $u_h \rightarrow u_0$ uniformly on K, which concludes the proof. $\square$

The next steps involve the precise form of the optimized energy, which for non-negative $u : \Omega \rightarrow {\mathbb R}$ reads

$$\begin{aligned} \mathcal{E}(u) = \Vert u\Vert _{L^1(\Omega )} + \frac{1}{2} \Vert \nabla u-z\Vert _{L^2(\Omega )} - \Vert z\Vert _{L^2(\Omega )} \end{aligned}$$

(38)

Lemma 8.2

(Upper estimate) One has $\limsup _{h \rightarrow 0} E_h \le E_0$.

Proof

Consider $u \in {{\mathrm{Conv}}}(\Omega )$, non-negative, and such that $\mathcal{E}(u) < \infty$. In particular, $u \in H^1(\Omega )$, by (38). Let $x_*$ be a fixed point in the interior of $\Omega$, let $\alpha >0$, and let $u_\alpha := u \circ D_\alpha$, where $D_\alpha$ is the dilation of center $x_*$ and factor $1/(1+\alpha )$. By construction $u_\alpha$ is convex, non-negative, and defined on a neighborhood of $\Omega$. Furthermore $\mathcal{E}(u_\alpha ) \rightarrow \mathcal{E}(u)$ as $\alpha \rightarrow 0$, because $u_\alpha \rightarrow u$ strongly in $H^1(\Omega )$.

Fix $\alpha$. Let $\beta >0$, and let $u_\alpha ^\beta := u_\alpha * K_\beta$, where $K_\beta (x) := \frac{1}{\beta ^2} K( \frac{x}{\beta })$, and K is a mollifier: a smooth non-negative function of unit integral and supported on the unit ball. Since $u_\alpha$ is defined on a neighborhood of $\Omega$, the convolution $u_\alpha * K_\beta$ is well defined on $\Omega$ for sufficiently small $\beta$. By construction $u_\alpha ^\beta$ is convex, non-negative, and smooth. Furthermore $\mathcal{E}(u_\alpha ^\beta ) \rightarrow \mathcal{E}(u_\alpha )$ as $\beta \rightarrow 0$, because $u_\alpha ^\beta \rightarrow u_\alpha$ strongly in $H^1(\Omega )$.

Fix $\alpha$ and $\beta$. Since $u_\alpha ^\beta$ is smooth one has ${{\mathrm{\mathrm{I}}}}_h u_\alpha ^\beta \rightarrow u_\alpha ^\beta$ strongly in $H^1(\Omega )$ as $h \rightarrow 0$, hence $\mathcal{E}({{\mathrm{\mathrm{I}}}}_h u_\alpha ^\beta ) \rightarrow \mathcal{E}(u_\alpha ^\beta )$. Thus $L := \limsup _{h \rightarrow 0} E_h \le \mathcal{E}(u_\alpha ^\beta )$. Letting $\beta \rightarrow 0$ yields $L\le \mathcal{E}(u_\alpha )$. Letting $\alpha \rightarrow 0$ yields $L \le \mathcal{E}(u)$. Optimizing over u yields the announced result. $\square$

Proposition 8.3

(Convergence of the discretization) One has $\lim _{h\rightarrow 0} E_h = E_0$. For each $h>0$ consider a minimizer $v_h \in {{\mathrm{Conv}}}(X_h)$ for $E_h$, and denote $u_h = {{\mathrm{\mathrm{I}}}}_h v_h$. Then $(u_h)_{h > 0}$ is bounded in $H^1(\Omega )$, and the limit of any weakly converging subsequence is a minimizer for $E_0$.

Proof

In view of (38), the minimizers $v_h$ exist by a finite dimensional compactness argument. Using again (38), and noting that $\mathcal{E}(u_h) \le \mathcal{E}(0) < \infty$, we obtain as announced that $u_h$ is bounded in $H^1(\Omega )$. Since $H^1(\Omega )$ is a Hilbert space, its bounded sets are pre-compact for the weak topology. Assuming, up to extraction, that $u_h \rightarrow u_0$ weakly in $H^1(\Omega )$ as $h\rightarrow 0$ we obtain

$$\begin{aligned}&\lim _{h \rightarrow 0} \int _{\Omega } \langle \nabla u_h,z\rangle = \int _{\Omega } \langle \nabla u_0, z\rangle , \quad \lim _{h \rightarrow 0} \Vert u_h\Vert _{L^1(\Omega )} = \Vert u_0\Vert _{L^1(\Omega )},\\&\limsup _{h \rightarrow 0} \Vert \nabla u_h\Vert _{L^2(\Omega )} \ge \Vert \nabla u_0\Vert _{L^2(\Omega )}. \end{aligned}$$

We used weak $H^1(\Omega )$ convergence for the first identity, the strong $L^1(\Omega )$ convergence for the second identity (due to the compact inclusion $H^1(\Omega ) \rightarrow L^1(\Omega )$), and the weak lower semi-continuity of the $H^1(\Omega )$ semi-norm for the third inequality. It follows that $\mathcal{E}(u_0) \le \limsup _{h \rightarrow 0} \mathcal{E}(u_h) = \limsup _{h \rightarrow 0} E_h \le E_0$, using Lemma 8.2. The limit $u_0$ non-negative, and convex by Proposition 8.1, hence it is a minimizer for $E_0$ which concludes the proof. Note that by Proposition 8.1 the sub-sequence’s convergence is in fact locally uniform. $\square$

Minor modifications of the proof allow to address other optimization problems.

• Projections in the $H^1(\Omega )$ norm Sect. 7.3: the proof applies with almost no change.

• Projections in the $L^2(\Omega )$ norm Sect. 7.3: the Hilbert space $H^1(\Omega )$ needs to be replaced with $L^2(\Omega )$ in all steps of the proof, which does not introduce new difficulties.

• Product bundles and lottery tickets Sect. 7.2. This model involves the Lipschitz regularity constraint ${\mathbb P}(u) : \forall x \in \Omega , \nabla u(x) \in [0,1]^2$. We make the additional assumption that each triangle $T \in \mathcal{T}_h$, $h>0$, has two edges aligned with the grid axes, so as to preserve this constraint: ${\mathbb P}(u) \Rightarrow {\mathbb P}({{\mathrm{\mathrm{I}}}}_h u)$. Lemma 8.2 is then easily adapted. In Proposition 8.3, replace Hilbert weak compactness with Arzela–Ascoli compactness, and rewrite the optimized functional without gradient terms:

  $$\begin{aligned} \int _\Omega (u - \langle \nabla u,z\rangle ) = \int _\Omega (u - {{\mathrm{div}}}(u z)+u{{\mathrm{div}}}z) = 3 \int _\Omega u - \int _{\partial \Omega } u \langle z,n\rangle . \end{aligned}$$

The convergence of higher order quantities, such as $\nabla u$ or $\det (\nabla ^2 u)$, is empirically observed but not claimed. The more complex functional (37) is not easily expressed using norms of classical functional spaces, and its analysis is cumbersome [6]; the convergence of the numerical method is thus left as an open question for this specific model.

Appendix B: Directional convexity

We briefly explore a weaker notion of discrete convexity [22, 23] based solely on the directional constraints $S_z^e$, and omitting the triangular constraints $T_z^e$.

Definition 9.1

We denote by ${{\mathrm{DConv}}}(X)$ the collection of elements in $\mathcal{F}(X)$ on which all the linear forms $S_x^e$, $x \in X$, $e\in {\mathbb Z}^2$ irreducible, supported on X, take non-negative values.

Some elements of ${{\mathrm{DConv}}}(X)$ cannot be extended into global convex maps on ${{\mathrm{Hull}}}(X)$. Their existence follows from the second point of Theorem 2.5 (minimality of the collection of constraints $S_x^e$, $T_x^e$), but for completeness we give (without proof) a concrete example.

Proposition 9.2

Let $u \in \mathcal{F}({\mathbb Z}^2)$ be defined by $u(1,1) = 1$, $u(-1,0)=u(0,-1) = -1$, and $u(x) := 2\Vert x\Vert ^2$ for other $x \in {\mathbb Z}^2$. Then for all $x \in {\mathbb Z}^2$, and all irreducible $e \in {\mathbb Z}^2$, one has $S_x^e (u) \ge 1$, and $T_x^e(u) \ge 2$ if $\Vert e\Vert >1$, with the exception $T_0^{(1,1)}(u) = -1$.

Elements of ${{\mathrm{DConv}}}(X)$ are nevertheless “almost” convex, in the sense that their restriction to a coarsened grid is convex. This might be a starting point for adapting the convergence analysis of Appendix A, which so far is limited to ${{\mathrm{Conv}}}(X)$.

Proposition 9.3

If $u \in {{\mathrm{DConv}}}(X)$, then $u_{|X'} \in {{\mathrm{Conv}}}(X')$, with $X' : =X \cap 2 {\mathbb Z}^2$.

Proof

Let $x \in X$, and let $e \in {\mathbb Z}^2$, $\Vert e\Vert >1$, be irreducible and of parents f, g. Assuming that $x+2 e, x-2 f, x-2 g \in X$, and observing that $x \pm e \in X$ by convexity, we obtain

$$\begin{aligned} u(x+2 e)+u(x-2 f) +u(x- 2 g) - 3 u(x) = 2 S_x^e(u) + S_{x+e}^e(u) + S_{x-e}^{f-g}(u) \ge 0. \end{aligned}$$

Likewise $u(x+2 e)- 2 u(x) + u(x-2 e) = S_{x-e}^e(u)+2 S_x^e(u)+S_{x+e}^e(u) \ge 0$. $\square$

The cone ${{\mathrm{DConv}}}(X)$ of directionally convex functions admits, just like ${{\mathrm{Conv}}}(X)$, a hierarchy of sub-cones ${{\mathrm{DConv}}}(\mathcal{V})$ associated to stencils. It is used in [3] to construct a monotone discretization of the Monge-Ampere operator.

Definition 9.4

Let $\mathcal{V}$ be a family of stencils on X, and let $u \in \mathcal{F}(X)$. The cone ${{\mathrm{DConv}}}(\mathcal{V})$ is defined by the non-negativity of the following linear forms: for all $x \in X$

• For all $e \in \mathcal{V}(x)$, the linear form $S_x^e$, if supported on X.

• For all $e \in \hat{\mathcal{V}}(x)$, the linear form $H_x^e := P_x^e + P_x^{-e}$, if supported on X.

Proposition 9.5

• For any stencils $\mathcal{V}, \mathcal{V}'$, one has ${{\mathrm{DConv}}}(\mathcal{V}) \subseteq {{\mathrm{DConv}}}(X)$, ${{\mathrm{DConv}}}(\mathcal{V}) \cap {{\mathrm{DConv}}}(\mathcal{V}') = {{\mathrm{DConv}}}(\mathcal{V}\cap \mathcal{V}')$, and ${{\mathrm{DConv}}}(\mathcal{V}) \cup {{\mathrm{DConv}}}(\mathcal{V}') \subseteq {{\mathrm{DConv}}}(\mathcal{V}\cup \mathcal{V}')$.

• For any stencils $\mathcal{V}$, one has

  $$\begin{aligned} {{\mathrm{DConv}}}(\mathcal{V}) \!=\! \{u \in {{\mathrm{DConv}}}(X); \, H_x^e(u) \!\ge \! 0 \text { for all } x \in X, \, e \in \mathcal{V}_{\max }(x) {\setminus } \mathcal{V}(x)\}.\nonumber \\ \end{aligned}$$

  (39)

Proof

As observed in Sect. 4.2, the second point of this proposition implies the first one. We denote by ${\mathbb P}(\mathcal{V})$ the identity (39), and prove it by decreasing induction over $\#(\mathcal{V})$. Since ${\mathbb P}(\mathcal{V}_{\max })$ clearly holds, we consider stencils $\mathcal{V}\subsetneq \mathcal{V}_{\max }$.

Let $x \in X$, $e \in \mathcal{V}_{\max }(x) {\setminus } \mathcal{V}(x)$, be such that $\Vert e\Vert$ is minimal. Similarly to Proposition 4.9 we find that e belongs to the set $\hat{\mathcal{V}}(x)$ of candidates for refinement at x, and define stencils $\mathcal{V}'$ by $\mathcal{V}(x) := \mathcal{V}(x) \cup \{e\}$, and $\mathcal{V}'(y) := \mathcal{V}(y)$ for $y \ne x$. The cones ${{\mathrm{DConv}}}(\mathcal{V})$ and ${{\mathrm{DConv}}}(\mathcal{V}')$ are defined by a common collection of constraints, with the addition respectively of $H_x^e$ for ${{\mathrm{DConv}}}(\mathcal{V})$, and $S_x^e$, $H_x^{e+f}$, $H_x^{e+g}$ for ${{\mathrm{DConv}}}(\mathcal{V}')$. Expressing the latter linear forms as combinations of those defining ${{\mathrm{DConv}}}(\mathcal{V})$

$$\begin{aligned} S_x^e = H_x^e + S_x^f + S_x^g,\quad H_x^{e+f} = H_x^e + S_{x+e}^f + S_{x-e}^f, \quad H_x^{e+g} = H_x^e + S_{x+e}^g + S_{x-e}^g, \end{aligned}$$

and observing that ${\mathbb P}(\mathcal{V}')$ holds by induction, we conclude the proof of ${\mathbb P}(\mathcal{V})$. $\square$

Appendix C: Numerical results

See Figs. 10, 11, 12, 13, 14.

Fig. 10

Results for the classical principal agent model, with a uniform density of customers on the domain obtained by rotating the square $[1,2]^2$ around its center (3 / 2, 3 / 2) by the indicated angle $\theta$. From top to bottom level sets of U (vanishing set indicated in white), level sets of $\det (\nabla ^2U)$ (allows to discriminate between the three categories of customers), density of products bought (which explodes on some curves), margin of the monopolist

Fig. 11

Three dimensional plot of the optimal U for the product-bundles variant of the monopolist problem, with respect to various customer distributions. Left uniform distribution on $[0,1]^2$. Center, and right uniform distribution on the illustrated black polygon

Fig. 12

Optimal pricing of risky assets, with the parameters $\alpha =1$, $\xi =1$, and a uniform customer density on $[1,2]^2$. Top left level sets of the estimated solution U of (37), with exclusion region $U=0$ in white. Top right level sets of $\det (\nabla ^2U)$, the darkest one is the estimated bunching region $\det (\nabla ^2U)=0$. Bottom left and right (detail): optimal product line, colored with the monopolist margin (color figure online)

Fig. 13

Top the projected function $u_0$ is the sum of three gaussian bumps arranged at the vertices of a triangle. The $L^2$ projection u inherits this triangular structure: curvature is concentrated at the triangle vertices, and the convexity is degenerate along the triangle edges as revealed by the u-Delaunay triangulation. Middle the projected function is a centered gaussian bump. The $L^2$ projection is non-smooth along the four diagonals, in contrast with the $H^1$ projection. Bottom $L^2$ projection onto the convex functions is used as a parameterless denoising algorithm, so as to recover the orientation and axes-ratio of a convex quadratic function. Experiments on $50\times 50$ grids, except for the $20\times 20$ grid top right

Fig. 14

Comparison of different numerical methods for the classical Monopolist problem