Is it possible to shape a piece of glass so that it refracts and concentrates sunlight in order to produce a given image? The modelling of this kind of problem leads to nonlinear second-order partial differential equations, which belong to the family of Monge–Ampère equations. We will see how semi-discrete methods, that can be traced back to Minkowski’s works, allow us to numerically solve such equations.
Figure 1. Mirror transforming a parallel, uniform light source into the shape of a train 🅭🅯 CC BY 4.0

1 Anidolic optics and Monge Ampère type equations

In anidolic optics, or non-imaging optics, one studies the design of devices that transfer light energy between a source and a target. The general problem is to design the shape of a mirror (or a lens) that reflects (or refracts) the light emitted from a given source towards a target whose geometry and intensity distributions are prescribed (see Figures 1 and 2). Applications of anidolic optics include the design of solar ovens, public lighting, car headlights, and more generally the optimization of the use of light energy and the reduction of light pollution.

Figure 2. Lens transforming a parallel light source into a hikari 🅭🅯 CC BY 4.0

Near-field and far-field light sources

There exists many different problems in anidolic optics, depending for instance on the geometry of the light source, the type of optical component, and the target to be illuminated. These problems are distinguished in particular by the spatial position of the target illumination. A problem is called near-field when the target is located within a finite distance, i.e., when one wishes to illuminate an area of space such as a screen. In Figures 1 and 2, the target illumination is on a wall, making the problem near-field. Most of the illustrations in this article correspond near-field targets. However, we will first consider the far-field case, which is mathematically simpler, and in which the target lies “at infinity”, in the space of directions. In practice, this means that when light is reflected or refracted from a point of the optical component, we can forget the spatial position of the reflected or refracted ray, keeping only its direction, which can then be encoded by a unit vector. Note that if the near-field target illumination is far away from the optical component, each point of the target almost corresponds to a direction, so that the far-field problem is a good approximation of the near-field one in this situation. We will see in Section 4.2 that one can solve a problem involving a near-field target by iteratively solving problems with far-field targets.

We will first present two far-field mirror problems in their continuous form, as illustrated in Figure 1. Then, we will explain how these continuous problems can be approached by discrete problems, following the so-called supporting quadric method introduced by Luis Caffarelli and Vladimir Oliker. This method can be traced back to work of Hermann Minkowski and Aleksandr Aleksandrov in convex geometry.

Figure 3. Mirror transforming light from a point source (left) or collimated light (right) 🅭🅯 CC BY 4.0

1.1 Mirror for a point light source

In this first problem, light is emitted from a point OO, which we assume to be located at the origin of the space R3\mathbb{R}^{3}. The intensity of the light source is modeled by a probability density μ\mu on the sphere of directions S2\mathbb{S}^{2}. Let XS2X\subset\mathbb{S}^{2} denote the support of the measure μ\mu. For example, if the light is emitted in a solid cone, then XX is a disc on the sphere. The quantity of light emanating from a measurable set of directions AXA\subset X is given by μ(A)\mu(A). For the far-field problem, the target is described by a probability measure ν\nu on the sphere of directions S2\mathbb{S}^{2}, which then represents the directions “at infinity”, i.e., after reflection. Let YS2Y\subset\mathbb{S}^{2} denote the support of the measure ν\nu.

The inverse problem considered here consists in constructing the surface R\mathcal{R} of a mirror which will transport the intensity μ\mu of the light source to the desired light distribution ν\nu at infinity using Snell’s law of reflection. For example, if the target measure is a Dirac mass δy\delta_{y}, meaning that we want to reflect all the light in a single direction yy, then the shape of the mirror is given by a paraboloid of revolution.

Let \langle\cdot\,|\,\cdot\rangle denote the Euclidean scalar product on R3\mathbb{R}^{3}. An incident ray xS2x\in\mathbb{S}^{2} is reflected by a surface R\mathcal{R} in the direction R(x)=xxn(x)n(x)R(x)=x-\langle x\,|\,n(x)\rangle n(x), where n(x)n(x) is the unit vector normal to the surface R\mathcal{R} at the point touched by the direction xx and oriented so that xn(x)0\langle x\,|\,n(x)\rangle\leq 0. The surface R\mathcal{R} solves the inverse mirror problem if RR transports the source measure μ\mu to the target measure ν\nu, in the sense that for any measurable subset BB of the sphere one has

BS2,ν(B)=μ(R1(B)).\forall B\subseteq\mathbb{S}^{2},\quad\nu(B)=\mu\bigl(R^{-1}(B)\bigr).

Note that the preservation of overall light quantity was already ensured by having chosen probability measures, i.e., μ(S2)=ν(S2)=1\mu(\mathbb{S}^{2})=\nu(\mathbb{S}^{2})=1. Now assume that μ\mu and ν\nu are absolutely continuous measures with respect to the area measure on the sphere. Let μ(x)=ρ(x)dx\mu(x)=\rho(x)\mathrm{d}x and ν(x)=σ(x)dx\nu(x)=\sigma(x)\mathrm{d}x, where ρ\rho and σ\sigma are the densities of μ\mu and ν\nu respectively. The previous equation then reads

BS2,Bσ(x)dx=R1(B)ρ(x)dx.\forall B\subseteq\mathbb{S}^{2},\quad\int_{B}\sigma(x)\mathrm{d}x=\int_{R^{-1}(B)}\rho(x)\mathrm{d}x.

Suppose furthermore that the densities ρ\rho et σ\sigma are continuous and that RR is a diffeomorphism from XX to YY. By the change of variable y=R(x)y=R(x), the last equation is then equivalent to σ(R(x))det(DR(x))=ρ(x)\sigma(R(x))\det(DR(x))=\rho(x) for any xXx\in X.

Since the mirror reflects rays emitted from the origin, we will assume that the surface R\mathcal{R} is radially parametrized by xS2u(x)xx\in\mathbb{S}^{2}\mapsto u(x)x, where u:S2R+u:\mathbb{S}^{2}\to\mathbb{R}^{+} is a positive function that must be determined. The unit normal to the surface R\mathcal{R} at the point xu(x)xu(x) and the direction of the reflected ray can both be expressed as a function of xx and of the gradient u(x)TxS2\nabla u(x)\in\mathrm{T}_{x}\mathbb{S}^{2}:

Ru(x)=xxnu(x)nu(x)andnu(x)=u(x)u(x)xu(x)2+u(x)2.R_{u}(x)=x-\langle x\,|\,n_{u}(x)\rangle n_{u}(x)\quad\text{and}\quad n_{u}(x)=\frac{\nabla u(x)-u(x)x}{\sqrt{\left\|\nabla u(x)\right\|^{2}+u(x)^{2}}}.

This allows us to formulate the problem as a system of partial differential equations, i.e., the problem of finding a positive function u:S2R+u:\mathbb{S}^{2}\to\mathbb{R}^{+} of class C2\mathcal{C}^{2} which satisfies

{σ(Ru(x))det(DRu(x))=ρ(x),Ru be a diffeomorphism from X to Y.\begin{cases}\sigma\bigl(R_{u}(x)\bigr)\det\bigl(DR_{u}(x)\bigr)=\rho(x),\\ \text{$R_{u}$ be a diffeomorphism from $X$ to $Y$.}\end{cases}

The first line of equation (Mir-Ponc-C) involves the determinant of a quantity which depends on the second derivatives of uu. This equation belongs to the family of Monge–Ampère equations. Note that the requirement that RuR_{u} is a diffeomorphism is non-standard and difficult to handle. In practice, it is replaced by a condition on uu which is akin to convexity, and by the so-called second boundary condition Ru(X)=YR_{u}(X)=Y. These two conditions ensure the ellipticity of the problem. Caffarelli and Oliker proved in 1994 [1 L. Caffarelli and V. Oliker, Weak solutions of one inverse problem in geometric optics. Journal of Mathematical Sciences154, 39–49 (2008) ] the existence of weak solutions to this equation, i.e., the existence of a locally Lipschitz function uu such that the application RR defined by the last two lines of (Mir-Ponc-C) satisfies (1). The existence of regular solutions to the problem (Mir-Ponc-C) is due to Wang and Guan [6 X.-J. Wang, On the design of a reflector antenna. Inverse Problems12, 351–375 (1996) , 2 P. Guan and X.-J. Wang, On a monge-ampere equation arising in geometric optics. J. Differential Geom48, 205–223 (1998) ].

1.2 Mirror for a collimated light source

We now present a second inverse problem arising in anidolic optics. This time the light source is collimated, which means that all the rays of light emitted by the source are parallel. We furthermore assume that they are positively collinear to the vertical vector ez=(0,0,1)e_{z}=(0,0,1) and emitted from a domain of the horizontal plane XR2×{0}X\subset\mathbb{R}^{2}\times\{0\}. For convenience, we will identify R2\mathbb{R}^{2} and R2×{0}\mathbb{R}^{2}\times\{0\}. We assume that the surface of the optical component is smooth and parametrized by a height function u:XRu:X\to\mathbb{R}. The intensity of the light source is modeled by a probability measure μ\mu on XX. As in the previous case, the intensity of the target illumination is modeled by a probability measure ν\nu on the sphere of directions at infinity. At each point (x,u(x))(x,u(x)) of the optical component, the gradient u(x)\nabla u(x) encodes the direction of the normal to the surface and we denote by F(u(x))S2F(\nabla u(x))\in\mathbb{S}^{2} the direction of the ray reflected by Snell’s law. The reflector defined by uu solves the inverse mirror problem between μ\mu and ν\nu if for any measurable set BS2B\subset\mathbb{S}^{2} one has

ν(B)=μ((Fu)1(B)).\nu(B)=\mu\bigl((F\circ\nabla u)^{-1}(B)\bigr).

Let us introduce the measure ν~\tilde{\nu} defined by ν~(B)=ν(F(B))\tilde{\nu}(B)=\nu(F(B)), which is supported on R2\mathbb{R}^{2}. We assume that μ\mu and ν~\tilde{\nu} are absolutely continuous with respect to the Lebesgue measure, with continuous densities ρ\rho and σ\sigma, and that xu(x)x\mapsto\nabla u(x) is a diffeomorphism on its image. Then, with the change of variable y=u(x)y=\nabla u(x), the inverse mirror problem for a collimated source can also be phrased as a partial differential equation:

{σ(u(x))det(D2u(x))=ρ(x),Fu is a diffeomorphism from X to Y.\begin{cases}\sigma\bigl(\nabla u(x)\bigr)\det\bigl(D^{2}u(x)\bigr)=\rho(x),\\ \text{$F\circ\nabla u$ is a diffeomorphism from $X$ to $Y$.}\end{cases}

We finally note that if uu is smooth and strongly convex (or strongly concave), the application u\nabla u is a diffeomorphism on its image.

1.3 Lenses

The construction of lenses that transform a light source into a target illumination prescribed at infinity are similar and are also described by Monge–Ampère type equations. As with mirrors, when the light is emitted from a point, the equation to be solved is on the sphere and when the light source is collimated, it is on the plane. In these problems, the light source passes through the surface of one side of the lens, either flat or spherical, and the aim is to construct the surface of the other side of the lens such that it refracts the light onto a prescribed target illumination at infinity. We do not detail the modelling of these problems here, but we will show results with lenses at the end.

2 Geometric discretization of Monge–Ampère equations

The two inverse problems in anidolic optics described in the previous section each involve two sets XX and YY, on which we have two probability measures μ\mu and ν\nu respectively, representing the light source and the desired target illumination. We saw that when these measures are absolutely continuous, the problems of construction of optical components correspond to partial differential equations of Monge–Ampère type.

The most direct method to solve a partial differential equation numerically is to approximate the domain XX with a discrete grid and to replace the partial derivatives with differences of the values of the function at the points of the grid divided by the grid step. In the case of Monge–Ampère equations, the application of these methods is made difficult by the non-linearity of the Monge–Ampère operator and by the diffeomorphism condition. We refer to the work of Adam Oberman, Brittany Froese, Jean-David Benamou and Jean-Marie Mirebeau for this line of research.

In recent years, alternative methods, called semi-discrete methods, have been used to discretize and numerically solve Monge–Ampère type equations arising from optimal transport. In order to apply this method, one assumes that one of the two measures μ\mu or ν\nu is absolutely continuous, while the other is finitely supported. Here we assume that μ\mu has a density μ(x)=ρ(x)dx\mu(x)=\rho(x)\mathrm{d}x on the space XX and that ν\nu is a discrete measure on the space Y={y1,,yN}Y=\{y_{1},\ldots,y_{N}\}, i.e., ν=1iNδyiνi\nu=\sum_{1\leq i\leq N}\delta_{y_{i}}\nu_{i} where δyi\delta_{y_{i}} is the Dirac mass in yiy_{i}.

In this section, we describe the semi-discrete variant of the two far-field mirror problems seen in the previous section, leaving aside the problem of convergence of the solutions to the discrete problems towards those to the continuous problems. These constructions give rise to equations which can naturally be seen as discrete Monge–Ampère equations. We also propose an economic interpretation by addressing the bakeries problem.

2.1 Mirror for a point light source

Let us go back to the problem of constructing mirrors that transform the light emitted by a point light source (see Section 1.1). As in the previous section, the light source is modeled by a continuous probability density ρ\rho on the sphere of directions S2\mathbb{S}^{2}, whose support Xρ:={xS2, ρ(x)>0}X_{\rho}:=\{x\in\mathbb{S}^{2},\ \rho(x)>0\} corresponds to the set of directions in which light is emitted. This time we assume that the desired target illumination is described by a probability measure ν=1iNδyiνi\nu=\sum_{1\leq i\leq N}\delta_{y_{i}}\nu_{i} supported on a set Y={y1,,yN}S2Y=\{y_{1},\ldots,y_{N}\}\subset\mathbb{S}^{2} of distinct directions. The problem is still to find the mirror surface R\mathcal{R} that will reflect the measure μ\mu onto the measure ν\nu under Snell’s law, but this time the target measure ν\nu is discrete.

Mirror composed of paraboloid pieces

We use the method of supporting paraboloids proposed by Caffarelli and Oliker in 1994 [1 L. Caffarelli and V. Oliker, Weak solutions of one inverse problem in geometric optics. Journal of Mathematical Sciences154, 39–49 (2008) ], which was originally developed to show the existence of weak solutions in the case where both measures are absolutely continuous. Caffarelli and Oliker’s idea is based on a well-known property of paraboloids of revolution: a paraboloid of revolution with focal point OO and direction yy reflects any ray coming from point OO to the direction yy. It is thus natural to seek to construct a mirror whose surface is composed of pieces of paraboloids, each paraboloid illuminating a direction yiy_{i}.

More precisely, we take ψ=(ψ1,,ψN)RN\psi=(\psi_{1},\ldots,\psi_{N})\in\mathbb{R}^{N} and denote by P(yi,ψi)P(y_{i},\psi_{i}) the solid (i.e. filled) paraboloid of direction yiy_{i}, with focal point at the origin OO and focal distance ψi\psi_{i}. This means that 12ψi\frac{1}{2}\psi_{i} is the distance between OO and the paraboloid’s closest point to OO. We define by Rψ\mathcal{R}_{\psi} the surface bordering the intersection of the solid paraboloids P(yi,ψi)P(y_{i},\psi_{i}):

Rψ=(1iNP(yi,ψi)).\mathcal{R}_{\psi}=\partial\left(\bigcap_{1\leq i\leq N}P(y_{i},\psi_{i})\right).

For each i{1,,N}i\in\{1,\ldots,N\} we denote by Vi(ψ)\mathrm{V}_{i}(\psi) the set of rays xS2x\in\mathbb{S}^{2} emitted by the light source and reflected by Snell’s law in the direction yiy_{i}. This set is called the ii-th visibility cell of the mirror Rψ\mathcal{R}_{\psi}. By construction, it corresponds to the radial projection of RψP(yi,ψi)\mathcal{R}_{\psi}\cap\partial P(y_{i},\psi_{i}) onto the sphere (see Figure 2.1).

A simple calculation shows that the intersection of two confocal paraboloids P(yi,ψi)\partial P(y_{i},\psi_{i}) and P(yj,ψj)\partial P(y_{j},\psi_{j}) is included in a plane curve. Projecting radially onto the unit sphere, this implies that the intersection of two visibility cells Vi(ψ)Vj(ψ)\mathrm{V}_{i}(\psi)\cap\mathrm{V}_{j}(\psi) is included in a curve on the sphere. We deduce that the set of visibility cells forms a partition of the sphere S2\mathbb{S}^{2}, up to a set of measure zero.

The paraboloid of revolution P(yk,ψk)\partial P(y_{k},\psi_{k}) can be parametrized radially by the function xS2xρk(x)x\in\mathbb{S}^{2}\mapsto x\rho_{k}(x), where ρk(x)=ψk/(1xyi)R\rho_{k}(x)=\psi_{k}/(1-\langle x\,|\,y_{i}\rangle)\in\mathbb{R}. We deduce that xx belongs to the visibility cell Vi(ψ)\mathrm{V}_{i}(\psi) if and only if the distance ρi(x)\rho_{i}(x) is smaller than the distances ρj(x)\rho_{j}(x) for j{1,,N}j\in\{1,\ldots,N\}. Composing with the logarithm to linearize the expression in ψ\psi, we obtain an explicit expression for the visibility cells

Vi(ψ)={xS2 j, c(x,yi)+ln(ψi)c(x,yj)+ln(ψj)},\mathrm{V}_{i}(\psi)=\left\{x\in\mathbb{S}^{2}\mid~{}\forall j,\ c(x,y_{i})+\ln(\psi_{i})\leq c(x,y_{j})+\ln(\psi_{j})\right\},

where c(x,y)=ln(1xy)c(x,y)=-\ln(1-\langle x\,|\,y\rangle).

By construction, each ray emitted by the point source and belonging to the cell Vi(ψ)\mathrm{V}_{i}(\psi) hits the mirror Rψ\mathcal{R}_{\psi} at a point which belongs to the paraboloid P(yi,ψi)\partial P(y_{i},\psi_{i}) and which is reflected in the direction yiy_{i}. The quantity of light received in the direction yiy_{i} is therefore exactly the quantity of light emanating from the visibility cell Vi(ψ)\mathrm{V}_{i}(\psi), i.e., μ(Vi(ψ))\mu(\mathrm{V}_{i}(\psi)). The desired quantity of light in the direction yiy_{i} is νi\nu_{i}. The equation to be solved is therefore μ(Vi(ψ))=νi\mu(\mathrm{V}_{i}(\psi))=\nu_{i} for any i{1,,N}i\in\{1,\ldots,N\}. Moreover, note that a paraboloid of revolution is only determined by its focal point, its direction and its focal distance. The free parameter remaining for each paraboloid P(yi,ψi)\partial P(y_{i},\psi_{i}) is the focal distance ψi\psi_{i}.

Figure 4. Mirror composed of three pieces of paraboloides reflecting in three directions 🅭🅯 CC BY 4.0

Formulation of the problem

The semi-discrete far-field mirror problem for a point source can be formulated as the problem of finding focal distances ψ=(ψ1,,ψN)RN\psi=(\psi_{1},\ldots,\psi_{N})\in\mathbb{R}^{N} that satisfy

i,  μ(Vi(ψ))=νi\forall i,~{}~{}\mu\bigl(\mathrm{V}_{i}(\psi)\bigr)=\nu_{i}

where c(x,y)=ln(1xy)c(x,y)=-\ln(1-\langle x\,|\,y\rangle) and where

Vi(ψ)={x j, c(x,yi)+ln(ψi)c(x,yj)+ln(ψj)}.\mathrm{V}_{i}(\psi)=\left\{x\mid~{}\forall j,\ c(x,y_{i})+\ln(\psi_{i})\leq c(x,y_{j})+\ln(\psi_{j})\right\}.

We will see in Section 3 how to solve such systems of equations. Note that if ψRN\psi\in\mathbb{R}^{N} is a vector of focal distances solving the mirror problem for a point-like source, then the surface of the mirror is parametrized by

Rψ: xS2  miniψi1xyi x.\mathcal{R}_{\psi}:\ x\in\mathbb{S}^{2}\ \mapsto\ \min_{i}\frac{\psi_{i}}{1-\langle x\,|\,y_{i}\rangle}\ x.

In the numerical experiments, we assume that the target illumination ν\nu is included in the half-sphere S2:={xS2, xez0}\mathbb{S}^{2}_{-}:=\{x\in\mathbb{S}^{2},\ \langle x\,|\,e_{z}\rangle\leq 0\}, that the support XρX_{\rho} of ρ\rho is included in the half-sphere S+2:={xS2, xez0}\mathbb{S}^{2}_{+}:=\{x\in\mathbb{S}^{2},\ \langle x\,|\,e_{z}\rangle\geq 0\}, and that the mirror is parametrized above the domain XρX_{\rho}.

Remark 2.1. The mirror surface is by construction the boundary of a convex set, i.e., the intersection of the solid paraboloids P(y1,ψ1),,P(yN,ψN)P(y_{1},\psi_{1}),\ldots,P(y_{N},\psi_{N}). It is also possible to construct a mirror contained in the boundary of the union of solid paraboloids rather than an intersection. This produces mirrors that are somewhat less interesting in practice, as they are neither convex nor concave.

2.2 Mirror for a collimated light source

Let us now consider the mirror problem for a collimated light source, already seen in Section 1.2. As before, the probability measure modelling the light source has a density ρ\rho with respect to the Lebesgue measure on the plane. However the probability measure modelling the target illumination intensity is discrete ν=1iNδyiνi\nu=\sum_{1\leq i\leq N}\delta_{y_{i}}\nu_{i}, supported on a finite set Y={y1,,yN}S2Y=\{y_{1},\ldots,y_{N}\}\subset\mathbb{S}^{2} of distinct directions. The problem is, again, to find the surface R\mathcal{R} of a mirror which reflects the measure μ\mu to the measure ν\nu.

Mirror with planar faces

We choose to construct the mirror surface R\mathcal{R} as the graph of affine height functions of the form xR2maxixpiψix\in\mathbb{R}^{2}\mapsto\max_{i}\langle x\,|\,p_{i}\rangle-\psi_{i} (see Figure 2.2). The vector pip_{i} is chosen so that the plane Pi={(x,xpi)xR2}P_{i}=\{(x,\langle x\,|\,p_{i}\rangle)\mid x\in\mathbb{R}^{2}\} reflects vertical rays, i.e., with direction eze_{z}, into the direction yiS2y_{i}\in\mathbb{S}^{2}. We need to determine the heights ψi\psi_{i} of those planes. Given a family of heights ψRN\psi\in\mathbb{R}^{N}, we define the ii-th visibility cell as

Vi(ψ)={xR2×{0} j, xpi+ψixpj+ψj}.\mathrm{V}_{i}(\psi)=\{x\in\mathbb{R}^{2}\times\{0\}\mid~{}\forall j,\ -\langle x\,|\,p_{i}\rangle+\psi_{i}\leq-\langle x\,|\,p_{j}\rangle+\psi_{j}\}.
Figure 5. Convex mirror for a collimated light source 🅭🅯 CC BY 4.0

By construction, for each i{1,,N}i\in\{1,\ldots,N\}, any vertical ray emitted from a point xVi(ψ)x\in V_{i}(\psi) hits the mirror R\mathcal{R} at a height xpiψi\langle x\,|\,p_{i}\rangle-\psi_{i} and is reflected in the direction yiy_{i}. Thus, the amount of light reflected in the direction yiy_{i} is equal to μ(Vi(ψ))\mu(\mathrm{V}_{i}(\psi)).

Formulation of the problem

Solving the semi-discrete far-field mirror problem for a collimated light source amounts to finding the heights ψRN\psi\in\mathbb{R}^{N} that satisfies

i,  μ(Vi(ψ))=νi\forall i,~{}~{}\mu\bigl(\mathrm{V}_{i}(\psi)\bigr)=\nu_{i}

where c(x,y)=xyc(x,y)=-\langle x\,|\,y\rangle and

Vi(ψ)={xj, c(x,yi)+ψic(x,yj)+ψj}.\mathrm{V}_{i}(\psi)=\left\{x\mid\forall j,\ c(x,y_{i})+\psi_{i}\leq c(x,y_{j})+\psi_{j}\right\}.

A solution of the equation (Mir-Colli-SD) induces a parametrization of the convex mirror R\mathcal{R} that reflects μ\mu onto ν\nu:

Rψ: xR2  (x,maxixpiψi)R3.\mathcal{R}_{\psi}:\ x\in\mathbb{R}^{2}\ \mapsto\ \bigl(x,\max_{i}\langle x\,|\,p_{i}\rangle-\psi_{i}\bigr)\in\mathbb{R}^{3}.

In practice, we only consider the part of the mirror located above the domain Xρ:={xR2×{0}, ρ(x)0}X_{\rho}:=\{x\in\mathbb{R}^{2}\times\{0\},~{}\rho(x)\neq 0\}.

Remark 2.2. The function Rψ\mathcal{R}_{\psi} being the maximum of affine functions, it is convex. The optical component which is parametrized by the graph of this application is also convex. Note that one could have the same construction by replacing the max\max in the formula by a min\min. This would result in a concave function Rψ\mathcal{R}_{\psi} and a concave mirror.

Remark 2.3. Problem (Mir-Colli-SD) is very similar (but not equivalent) to Minkowski’s problem in convex geometry which is also an inverse problem. Given a set of unit vectors yiy_{i} and real numbers νi>0\nu_{i}>0, this problem consists in building a convex polyhedron whose ii-th facet has normal yiy_{i} and area νi\nu_{i} – which is possible only under some assumptions on the directions and areas. We also note that Oliker, who was the first to introduce semi-discrete methods for the numerical resolution of Monge–Ampère equations, was a doctoral student of the famous geometer Aleksandrov who is known (among other) for introducing and studying the “continuous” formulation of Minkowski’s problem.

Figure 6. Bakeries: The city XX with its boundary drawn in blue is endowed with a probability density pictured in grayscale representing the population density. The set YY (in red) represents the location of bakeries. Here, X,YR2X,Y\subseteq\mathbb{R}^{2} and c(x,y)=xy2c(x,y)=|x-y|^{2}. We see the Voronoi tessellation of the city (in the middle, uniform price) as well as its Laguerre tessellation (on the right, only the bread price ψ1\psi_{1} has increased). 🅭🅯 CC BY 4.0

2.3 The bakeries problem

We now present an economic analogy which leads to an equation having the same structure as in the two optical problems presented above. We assume that XX represents a city whose population density is described by a probability density ρ\rho, that the finite set Y={y1,,yN}Y=\{y_{1},\ldots,y_{N}\} represents the locations of the city’s bakeries and that νi\nu_{i} represents the quantity of bread available in bakery yiy_{i}. Customers living at a location xx in XX naturally will look for the bakery minimizing the cost of walking from xx to yiy_{i}, denoted c(x,yi)c(x,y_{i}). This leads to a decomposition of the city space into Voronoi cells,

Vori:={xΩXj, c(x,yi)c(x,yj)}.\operatorname{Vor}_{i}:=\bigl\{x\in\Omega_{X}\mid\forall j,\ c(x,y_{i})\leq c(x,y_{j})\bigr\}.

The number of customers going to a bakery yiy_{i} is equal to the integral of the density ρ\rho over Vori\operatorname{Vor}_{i}. Suppose that a bakery yiy_{i} receives too many customers in comparison to its bread’s production capacity νi\nu_{i} – this could be the case in Figure 6 for the downtown bakery y1y_{1} where the population density is high. This means we have μ(Vor1)>ν1\mu(\operatorname{Vor}_{1})>\nu_{1}, where we denote μ(x)=ρ(x)dx\mu(x)=\rho(x)\mathrm{d}x. The baker y1y_{1} will then seek to increase the price of his bread. This will reduce the number of potential customers, but will increase the baker’s profit as long as he manages to sell all his stock. We write νi0\nu_{i}\geq 0 for the proportion of the population that the bakery yiy_{i} is able to serve, and ψi\psi_{i} the price of the bread in the bakery yiy_{i}. If we assume that the customers living at point xx make a compromise between walking cost and price of bread by minimizing the sum (c(x,yi)+ψic(x,y_{i})+\psi_{i}), the city is then decomposed into Laguerre cells

Lagi(ψ)={xX j, c(x,yi)+ψic(x,yj)+ψj}.\operatorname{Lag}_{i}(\psi)=\left\{x\in X\mid~{}\forall j,\ c(x,y_{i})+\psi_{i}\leq c(x,y_{j})+\psi_{j}\right\}.

Note that we do not necessarily have yiLagi(ψ)y_{i}\in\operatorname{Lag}_{i}(\psi), and that it is even possible to have Lagi(ψ)=\operatorname{Lag}_{i}(\psi)=\emptyset: indeed, if the bread is very expensive in a certain bakery, even people living next door may prefer going to a more distant one.

Problem formulation

The bakeries problem therefore boils down to finding a price vector ψRN\psi\in\mathbb{R}^{N} such that each bakery sells all its stock of bread νi\nu_{i}. This is described by the system of equations

μ(Lagi(ψ))=νii{1,,N},\mu\bigl(\operatorname{Lag}_{i}(\psi)\bigr)=\nu_{i}\quad\forall i\in\{1,\ldots,N\},

This equation has exactly the same structure as (Mir-Ponc-SD) and (Mir-Colli-SD). We will see in the next section how to solve this class of equations.

3 Numerical resolution

The discrete problems mentioned in the previous section all show the same structure; our focus will now be on their numerical resolution. We start by introducing the semi-discrete Monge–Ampère equation, and show that its solution is equivalent to finding the maximum of a concave function. Subsequently, we present a Newton method that allows us to solve these equations efficiently.

3.1 Semi-discrete Monge–Ampère equation

Let XX be a compact subset of the space R2\mathbb{R}^{2} or of the sphere S2\mathbb{S}^{2}, let Y={y1,,yN}Y=\{y_{1},\ldots,y_{N}\}, and let cC1(X×Y)c\in\mathcal{C}^{1}(X\times Y) be a cost function. The Laguerre cell (which corresponds to a visibility cell in optics) associated with a family of real numbers ψ=(ψ1,,ψN)RN\psi=(\psi_{1},\ldots,\psi_{N})\in\mathbb{R}^{N} is given by

Lagi(ψ)={xX j, c(x,yi)+ψic(x,yj)+ψj}.\operatorname{Lag}_{i}(\psi)=\left\{x\in X\mid~{}\forall j,\ c(x,y_{i})+\psi_{i}\leq c(x,y_{j})+\psi_{j}\right\}.

Suppose that the cost function satisfies the Twist condition

xX,yxc(x,y) is injective,\forall x\in X,\quad y\mapsto\nabla_{x}c(x,y)\ \textrm{is injective},

which ensures that the Laguerre cells form a partition of the domain XX up to a negligible set.

Semi-discrete Monge–Ampère equation

Let μ\mu be a probability measure on XX with density ρ\rho with respect to the area measure, and let ν=iνiδyi\nu=\sum_{i}\nu_{i}\delta_{y_{i}} be a probability measure on YY. In the following equation, the discrete probability measure ν\nu is conflated with the vector ν=(νi)1iN\nu=(\nu_{i})_{1\leq i\leq N}. We are seeking ψRN\psi\in\mathbb{R}^{N} satisfying


where the function G:RNRNG:\mathbb{R}^{N}\to\mathbb{R}^{N} is defined by

G(ψ)=(G1(ψ),,GN(ψ))andGi(ψ)=μ(Lagi(ψ)).G(\psi)=\bigl(G_{1}(\psi),\ldots,G_{N}(\psi)\bigr)\quad\text{and}\quad G_{i}(\psi)=\mu\bigl(\operatorname{Lag}_{i}(\psi)\bigr).

Remark 3.1. The visibility cells used in optics in (Mir-Ponc-SD) and (Mir-Colli-SD) are Laguerre cells, with

c(x,y)=log(1xy)andc(x,y)=xyc(x,y)=-\log\bigl(1-\langle x\,|\,y\rangle\bigr)\quad\text{and}\quad c(x,y)=-\langle x\,|\,y\rangle

respectively. Equation (MA) is a reformulation of equations (Mir-Ponc-SD) and (Mir-Colli-SD). Note that the Laguerre cells are invariant under addition of a constant to ψ\psi, and that the solution of (MA) is therefore defined up to an additive constant. Optical problems have a similar invariance: for example, if a surface R\mathcal{R} is a solution of the mirror problem for a point source, then so is λR\lambda\mathcal{R} for all λ>0\lambda>0.

3.2 Variational formulation

The following theorem shows that the function GG in the semi-discrete Monge–Ampère equation is the gradient of a concave function.

Theorem 3.1.

We assume that the cost function cc satisfies (Twist). Then the function K:RNR\mathcal{K}:\mathbb{R}^{N}\to\mathbb{R} defined by

K(ψ)=1iNLagi(ψ)(c(x,y)+ψi)ρ(x)dx1iNψiνi\mathcal{K}(\psi)=\sum_{1\leq i\leq N}\int_{\operatorname{Lag}_{i}(\psi)}\bigl(c(x,y)+\psi_{i}\bigr)\rho(x)\mathrm{d}x-\sum_{1\leq i\leq N}\psi_{i}\nu_{i}

is concave, of class C1\mathcal{C}^{1} and with gradient

K(ψ)=G(ψ)ν=(μ(Lagi(ψ))νi)1iN.\nabla\mathcal{K}(\psi)=G(\psi)-\nu=\bigl(\mu\bigl(\operatorname{Lag}_{i}(\psi)\bigr)-\nu_{i}\bigr)_{1\leq i\leq N}.

As we will see in the next paragraph, the function K\mathcal{K} is related to the Kantorovitch duality in optimal transport theory, and we will therefore call it the Kantorovitch functional. Moreover, since a concave function of class C1\mathcal{C}^{1} reaches its maximum exactly at its critical points, we obtain the following corollary:

Corollary 3.2. Under the assumptions of Theorem 3.1, a vector ψRN\psi\in\mathbb{R}^{N} is a solution to equation (MA) if and only if ψ\psi is a maximizer of K\mathcal{K}.

Since the function K\mathcal{K} is invariant under addition of a constant, one can choose to work on the set M0\mathcal{M}_{0} of vectors whose coordinates sum to zero. It can be shown that the function K\mathcal{K} is proper on M0\mathcal{M}_{0}, i.e., limψ+,ψM0K(ψ)=\lim_{\left\|\psi\right\|\to+\infty,\psi\in\mathcal{M}_{0}}\mathcal{K}(\psi)=-\infty, which ensures that it reaches its maximum: the problem (MA) thus admits a solution.

3.3 Relation to optimal transport

The variational formulation of the Monge–Ampère equation, i.e., the search for a maximizer of the Kantorovitch functional, corresponds in fact to the dual of the Monge–Kantorovitch problem in optimal transport theory. We discuss this link in detail below in the semi-discrete case. The reader interested in the proofs may refer for instance to the book chapter [4 Q. Mérigot and B. Thibert, Optimal transport: Discretization and algorithms. In Handbook of Numerical Analysis, Geometric PDEs 22, Elsevier, 621–679 (2020) ].

Monge’s problem

The image of a probability measure μ\mu on XX under a measurable application T:XYT:X\to Y is the measure T#μT_{\#}\mu on YY defined by T#μ(B)=μ(T1(B))T_{\#}\mu(B)=\mu(T^{-1}(B)). If T#μ=νT_{\#}\mu=\nu, we say that TT transports μ\mu to ν\nu. Since the set YY is finite, we have T#μ=1iNμ(T1(yi))δyiT_{\#}\mu=\sum_{1\leq i\leq N}\mu(T^{-1}(y_{i}))\delta_{y_{i}}. Monge’s optimal transport problem consists in finding a transport map TT that transports μ\mu to ν\nu and that minimizes the total cost Xc(x,T(x))dμ(x)\int_{X}c(x,T(x))d\mu(x). If the cost function cc satisfies the Twist condition, Brenier and Gangbo–McCann, relying on Kantorovich duality, proved the existence of a minimizer for this problem when the source μ\mu is absolutely continuous. For example, one can state the following:

Theorem 3.3 (Kantorovitch duality).

Suppose that cc satisfies the condition (Twist) and that μ\mu is absolutely continuous. Then

minT:XYTμ=νXc(x,T(x))dμ(x)=maxψRNK(ψ).\min_{{T:X\to Y}\atop{T_{\sharp}\mu=\nu}}\int_{X}c\bigl(x,T(x)\bigr)d\mu(x)=\max_{\psi\in\mathbb{R}^{N}}\mathcal{K}(\psi).

If moreover ψ\psi is a maximizer of K\mathcal{K}, then the function Tψ:XYT_{\psi}:X\to Y defined μ\mu-a.e. by TψLagy(ψ)=yT_{\psi}|_{\operatorname{Lag}_{y}(\psi)}=y realizes the minimum in Monge’s problem.

Remark 3.2. Not all Monge–Ampère equations derive from an optimal transport problem and not all of them admit a variational formulations. These two strong properties come in fact from the very particular structure of Laguerre cells, which follows from the functions ψc(x,y)+ψ(y)\psi\mapsto c(x,y)+\psi(y) being affine.

We saw that the far-field optics problems presented in Section 2 possess this structure. On the other hand, if we consider the mirror construction problems for a target illumination in the near-field (i.e., we are illuminating points in R3\mathbb{R}^{3} and not directions at infinity), we still have semi-discrete Monge–Ampère equations to solve, but the Laguerre cells are of the form

Lagi(ψ)={xXj, G(x,yi,ψi)G(x,yj,ψj)},\operatorname{Lag}_{i}(\psi)=\bigl\{x\in X\mid\forall j,\ G(x,y_{i},\psi_{i})\leq G(x,y_{j},\psi_{j})\bigr\},

where the function GG is nonlinear in ψ\psi. These equations do not derive from the optimal transport problem and in fact do not admit a variational formulation. They are called prescribed Jacobian equations by Trudinger, and are the subject of recent research both in analysis and in more applied fields (optics, economics).

3.4 Laguerre cells and derivatives

Before applying Newton’s method to solve the equation G(ψ)=νG(\psi)=\nu, we need to show that the function GG is of class C1\mathcal{C}^{1} (or equivalently that K\mathcal{K} is of class C2\mathcal{C}^{2}), calculate its partial derivatives and study the (strict) concavity of K\mathcal{K}. To do this, we need a genericity assumption which is somewhat technical, but which is natural and not restrictive in practice. In the optical cases mentioned in this paper, this assumption is satisfied if the intersection of three distinct Laguerre cells is finite and if the intersection of two Laguerre cells with the boundary of XX is also finite. For more details, the reader may refer to the book chapter [4 Q. Mérigot and B. Thibert, Optimal transport: Discretization and algorithms. In Handbook of Numerical Analysis, Geometric PDEs 22, Elsevier, 621–679 (2020) ].

Theorem 3.4 (Differental of GG).

Suppose that the cost satisfies (Twist), that YY is generic (see above), and ρ\rho is continuous. Then the application G:RNRNG:\mathbb{R}^{N}\to\mathbb{R}^{N} is of class C1\mathcal{C}^{1} and

ji,  Giψj(ψ)=Lagij(ψ)ρ(x)xc(x,yi)xc(x,yj)dx,\displaystyle\forall j\neq i,~{}~{}\frac{\partial G_{i}}{\partial\psi_{j}}(\psi)=\int_{\operatorname{Lag}_{ij}(\psi)}\frac{\rho(x)}{\left\|\nabla_{x}c(x,y_{i})-\nabla_{x}c(x,y_{j})\right\|}\mathrm{d}x,
i,  Giψi(ψ)=jiGiψj(ψ),\displaystyle\forall i,~{}~{}\frac{\partial G_{i}}{\partial\psi_{i}}(\psi)=-\sum_{j\neq i}\frac{\partial G_{i}}{\partial\psi_{j}}(\psi),

where Lagij(ψ)=Lagi(ψ)Lagj(ψ)\operatorname{Lag}_{ij}(\psi)=\operatorname{Lag}_{i}(\psi)\cap\operatorname{Lag}_{j}(\psi).

The formula for the partial derivatives of GG has a geometric interpretation. In the following two figures, which are obtained for the cost c(x,y)=xy2c(x,y)=\left\|x-y\right\|^{2} on R2\mathbb{R}^{2}, we explain why the formula for partial derivatives involves integrals over the interfaces between Laguerre cells and how the singularities of DGDG may occur depending on the geometry of the points yiy_{i}.

Figure 3.4 illustrates that the partial derivative Gi/ψj(ψ)\partial G_{i}/\partial\psi_{j}(\psi) is an integral over the interface Lagij(ψ)\operatorname{Lag}_{ij}(\psi): the value Gi(ψ)G_{i}(\psi) is an integral over the Laguerre cell Lagi(ψ)\operatorname{Lag}_{i}(\psi) (in grey on the left); we increase the value ψj\psi_{j} by ε>0\varepsilon>0 considering ψ+εej\psi+\varepsilon e_{j}; the rate of increase (Gi(ψ)Gi(ψ+εej))/ε(G_{i}(\psi)-G_{i}(\psi+\varepsilon e_{j}))/\varepsilon is proportional to an integral over the symmetric difference between two Laguerre cells (in grey in the middle); passing to the limit we obtain an integral over the green segment Lagij(ψ)\operatorname{Lag}_{ij}(\psi). The signs that occurs in the formula for the partial derivatives can also be interpreted with the bakeries metaphor: when the price of bread ψi\psi_{i} increases, the number of customers of the bakery yiy_{i} decreases (i.e., the Laguerre cell Lagi(ψ)\operatorname{Lag}_{i}(\psi) shrinks) and the number of customers for the other bakeries increases, so that Gi/ψi(ψ)0\partial G_{i}/\partial\psi_{i}(\psi)\leq 0 and Gi/ψj(ψ)0\partial G_{i}/\partial\psi_{j}(\psi)\geq 0 for jij\neq i.

In Figure 3.4, the genericity condition is not satisfied because y1y_{1}, y2y_{2} and y3y_{3} are aligned, and there exists ψRN\psi\in\mathbb{R}^{N} for which Lag1(ψ)Lag2(ψ)Lag3(ψ)\operatorname{Lag}_{1}(\psi)\cap\operatorname{Lag}_{2}(\psi)\cap\operatorname{Lag}_{3}(\psi) is a line segment. The partial derivative G2/ψ3(ψ)\partial G_{2}/\partial\psi_{3}(\psi) is an integral on the (green) segment Lag23(ψ)\operatorname{Lag}_{23}(\psi). If we simultaneously decrease ψ1\psi_{1} and ψ2\psi_{2} by the same amount, we can see that the segment Lag23(ψ)\operatorname{Lag}_{23}(\psi) varies continuously and then suddenly becomes empty when the cell Lag2(ψ)\operatorname{Lag}_{2}(\psi) gets empty (bottom right of Figure 3.4). Thus, G2/ψ3(ψ)\partial G_{2}/\partial\psi_{3}(\psi) is not continuous. Newton’s method requires a certain regularity, and we will see below that it converges under the above genericity assumptions.

Figure 7. The partial derivatives are boundary integrals 🅭🅯 CC BY 4.0

Figure 8. Non-continuous partial derivative: G2/ψ3\partial G_{2}/\partial\psi_{3} is an integral on the green segment Lag23\operatorname{Lag}_{23} which is discontinuous. 🅭🅯 CC BY 4.0

To establish the convergence of Newton’s method, we also need to study the concavity of the Kantorovitch functional K\mathcal{K}, or equivalently the monotonicity of its gradient K=Gν\nabla\mathcal{K}=G-\nu. The functions K\mathcal{K} and GG are invariant by addition of a constant vector (i.e., K(ψ+C(1,,1))=K(ψ)\mathcal{K}(\psi+C(1,\ldots,1))=\mathcal{K}(\psi)), which can be seen in the definition of Laguerre cells. Thus, we can only hope to establish strong concavity of K\mathcal{K} in the directions orthogonal to constant vectors, i.e., belonging to the set

M0:={vRN1iNvi=0}.\mathcal{M}_{0}:=\Bigl\{v\in\mathbb{R}^{N}\mid\sum_{1\leq i\leq N}v_{i}=0\Bigr\}.

Another reason for the lack of strong concavity of K\mathcal{K} is that if ψi\psi_{i} is very large, then Lagi(ψ)\operatorname{Lag}_{i}(\psi) is empty and remains empty in a neighborhood of ψ\psi. In this case, Gi(ψ)G_{i}(\psi) is constant equal to zero, and the Hessian matrix D2G(ψ)=DGD^{2}G(\psi)=DG has a row of zeros. We can therefore hope to establish strong concavity only if ψ\psi belongs to the set

C+:={ψRNi, Gi(ψ)>0}.\mathcal{C}_{+}:=\bigl\{\psi\in\mathbb{R}^{N}\mid\forall i,\ G_{i}(\psi)>0\bigr\}.

The next theorem shows, in a nutshell, that these are the only two obstructions to strong concavity.

Theorem 3.5 (Strict concavity).

We assume the hypotheses of the previous theorem hold. If the set {ρ>0}\{\rho>0\} is connected, the function K\mathcal{K} is locally strongly concave on C+\mathcal{C}_{+} in the direction M0\mathcal{M}_{0}:

ψC+, vM0{0},DG(ψ)vv<0.\forall\psi\in\mathcal{C}_{+},~{}\forall v\in\mathcal{M}_{0}\setminus\{0\},\quad\langle DG(\psi)v\,|\,v\rangle<0.

Remark 3.3 (Uniqueness). We saw above that the function K\mathcal{K} has a maximum, and thus equation (MA) has a solution. The previous theorem implies that this maximum is unique if we impose that ψM0\psi\in\mathcal{M}_{0}, i.e., ψ\psi has zero average, since a strictly concave function admits at most one local maximum.

We will see in the next paragraph how these results of regularity and monotonicity allow us to iteratively construct a sequence (ψ(k))k0(\psi^{(k)})_{k\geq 0} converging to the unique zero-average ψ\psi^{*} satisfying G(ψ)=νG(\psi^{*})=\nu.

3.5 Newton’s method

Newton’s method in 1D

We begin by recalling Newton’s method for solving the equation g(x)=0g(x)=0, where g:RRg:\mathbb{R}\to\mathbb{R} is a real function. Newton’s method starts from x0Rx_{0}\in\mathbb{R} and constructs the sequence xk+1=xkg(xk)/g(xk)x^{k+1}=x^{k}-g(x^{k})/g^{\prime}(x^{k}) by induction. If we assume that gg is of class C1\mathcal{C}^{1} and that there exists aRa\in\mathbb{R} such that g(a)=0g(a)=0 and g(a)0g^{\prime}(a)\neq 0, then one can show, using Taylor–Lagrange formulas, that for x0x^{0} sufficiently close to aa, the sequence (xk)k0(x^{k})_{k\geq 0} converges to aa. The convergence is then said to be local. Thus, under a regularity hypothesis (gC1g\in\mathcal{C}^{1}) and monotonicity (gg^{\prime} has constant sign in a neighborhood of aa), Newton’s method converges locally.

Newton’s method (local)

Assume that we are given a zero-average vector ψ0M0\psi^{0}\in\mathcal{M}_{0} such that the mass of all Laguerre cells is strictly positive:

ε0:=12min[minyYGi(ψ0), min1iNνyi]>0.\varepsilon_{0}:=\frac{1}{2}\min\left[\min_{y\in Y}G_{i}(\psi^{0}),~{}\min_{1\leq i\leq N}\nu_{y_{i}}\right]>0.

We define ψk+1\psi^{k+1} in the following way: we start by calculating the Newton direction dkd^{k}, i.e., the vector dkd^{k} satisfying

DG(ψk)dk=(G(ψk)ν)anddikM0,DG(\psi^{k})d^{k}=-(G(\psi^{k})-\nu)\quad\text{and}\quad d_{i}^{k}\in\mathcal{M}_{0},

which exists and is unique by according to Theorem 3.5. The second equation enables us to overcome the invariance of GG and thus the non-invertibility of DG(ψk)DG(\psi^{k}). We then define ψk+1=ψk+dk\psi^{k+1}=\psi^{k}+d^{k}. As in the 1D case, it can be shown that the method converges locally: if ψ0\psi^{0} is chosen close enough to the ψ\psi^{*} solution, then the sequence (ψk)(\psi^{k}) converges to ψ\psi^{*}.

Globally convergent Newton’s method

However, the condition ψ0\psi^{0} is close to the solution ψ\psi^{*} is impossible to fulfill in practice. Fortunately, a very simple modification of the method allows to ensure a global convergence, allowing us to drop this closeness assumption. To do this, one must construct ψk+1\psi^{k+1} in such a way that the kernel of the Jacobian DG(ψk+1)DG(\psi^{k+1}) remains equal to constant vectors, so that the system defining the direction dk+1d^{k+1} admits a unique solution. For this purpose, we define the stepτk\tau^{k} as the largest real of the form 22^{-\ell} (with N\ell\in\mathbb{N}) such that ψk,:=ψk+2dk\psi^{k,\ell}:=\psi^{k}+2^{-\ell}d^{k} satisfies

{i{1,,N},Gi(ψk,)ε0,G(ψk,)ν(12(+1))G(ψk)ν.\begin{cases}\forall i\in\{1,\ldots,N\},\quad G_{i}(\psi^{k,\ell})\geq\varepsilon_{0},\\ \left\|G(\psi^{k,\ell})-\nu\right\|\leq(1-2^{-(\ell+1)})\left\|G(\psi^{k})-\nu\right\|.\end{cases}

Finally, we define ψk+1=ψk+τkdk\psi^{k+1}=\psi^{k}+\tau^{k}d^{k}.

By using the regularity and concavity results on K\mathcal{K}, the step τk\tau^{k} can be bounded from below, thus ensuring the convergence of the sequence constructed above to a solution of the optimal transport problem [4 Q. Mérigot and B. Thibert, Optimal transport: Discretization and algorithms. In Handbook of Numerical Analysis, Geometric PDEs 22, Elsevier, 621–679 (2020) ]:

Theorem 3.6.

Under the assumptions of Theorem 3.5, there exists τ>0\tau^{*}>0 such that


In particular, the sequence (ψk)k0(\psi^{k})_{k\geq 0} converges to the unique solution ψ\psi^{*} of (MA) satisfying iψi=0\sum_{i}\psi_{i}^{*}=0.

Remark 3.4 (Quadratic convergence). The above theorem shows that the convergence of Newton’s method is globally exponential. This convergence is actually called linear convergence in optimization. When the cost cc satisfies the Ma–Trudinger–Wang (MTW) condition that appears in the theory of optimal transport regularity, and the density ρ\rho is Lipschitz, then the convergence is even locally quadratic [3 J. Kitagawa, Q. Mérigot and B. Thibert, Convergence of a Newton algorithm for semi-discrete optimal transport. Journal of the European Mathematical Society (2019) ]: for sufficiently large kk, we have


In practice, the convergence is very fast and the basin where quadratic convergence occurs seems to be quite large. This last observation is empirical, and not mathematically explained yet. In Figure 9, X=[0,1]2X=[0,1]^{2} is the large white square and YY is a set of points in the lower left corner and c(x,y)=xy2c(x,y)=\left\|x-y\right\|^{2}. With N=100N=100 points, after three iterations the error G(ψ3)ν1\|G(\psi^{3})-\nu\|_{1} is already of order 10910^{-9}. Even difficult examples of size N=107N=10^{7} in dimension d=3d=3 can be solved to high numerical precision with less than 20 iterations!

4 Applications to anidolic optics

In this section we present the adaptation of semi-discrete methods to the practical resolution of inverse problems in optics. These results were obtained in the PhD thesis of Jocelyn Meyron and the images are taken from the article [5 J. Meyron, Q. Mérigot and B. Thibert, Light in power: A general and parameter-free algorithm for caustic design. In SIGGRAPH Asia 2018 Technical Papers, ACM Transaction on Graphics, p. 224 (2018) ].

Figure 9. Convergence of the sequence (ψk)(\psi^{k}). On images 2, 3 and 4 we see the Laguerre cells Lagi(ψk)\operatorname{Lag}_{i}(\psi^{k}) for k=0,1,3k{=}0{,}1{,}3.

4.1 Far-field problems

We saw in Section 2 that in several far-field problems, i.e., when the target illumination is at infinity, solving the Monge–Ampère equation (MA) allows us to construct an optical component. This involves modelling mirrors or lenses, with a point or collimated light source, and in each case there are two components that may be produced (one of which is convex), so that in all we have formulated eight different near-field optical problems.

The main difficulty in implementing Newton’s algorithm to solve (MA) lies in the evaluation of the function GG and its differential DGDG at point ψk\psi^{k}, and more precisely in the calculation of the set of Laguerre cells Lagi(ψk)\operatorname{Lag}_{i}(\psi^{k}). For cells from non-imaging optics problems, also called visibility cells, it is possible to perform this calculation in almost linear time in the number NN of Dirac masses. Take for example the mirror problem for a point source. The visibility cells are obtained by projecting radially onto the sphere an intersection of “solid” confocal paraboloids, and we have already seen that the intersection of two confocal paraboloids is included in a plane. Another simple calculation shows that the radial projection of such an intersection is also included in a (different) plane. This shows that the visibility cells are separated by hyperplanes. In fact, it can be shown that there exists a partition of R3\mathbb{R}^{3} into convex polyhedra P1,,PNP_{1},\ldots,P_{N} – called a power diagram in computational geometry – such that each visibility cell is of the form Vi(ψ)=S2PiV_{i}(\psi)=\mathbb{S}^{2}\cap P_{i} (Figure 4.1). A similar property holds for each of the eight problems. The point of this reformulation is that there are powerful libraries – for example Cgal or Geogram – that allow us to compute power diagrams in dimensions 22 and 33, and thus also the Laguerre cells associated with the optics problems. It is therefore possible to implement the damped Newton algorithm, and to use it to construct – numerically and even physically – mirrors and lenses for far-field targets in anidolic optics.

Figure 10. Visibility cell structure 🅭🅯 CC BY 4.0

4.2 Near-field problems

It is also possible to deal with more realistic target illuminations in the near-field – i.e., when illuminating points at a finite distance rather than directions – with an iterative method that solves a far-field solution at each step [5 J. Meyron, Q. Mérigot and B. Thibert, Light in power: A general and parameter-free algorithm for caustic design. In SIGGRAPH Asia 2018 Technical Papers, ACM Transaction on Graphics, p. 224 (2018) ]. The convergence is very fast, requiring only a few iterations, as illustrated in Figure 4.2.

In all the experiments presented below, the light source is assumed to be uniform, so that the light source μ\mu has a constant density on its support. The reflection or refraction of this light on a wall is simulated in the computer by the physically realistic rendering software LuxRender.

Generic method

The different problems of anidolic optics having the same structure (point or collimated light sources, mirrors or lenses, convex or concave components, near-field or far-field), it is possible to solve them in a unified, precise and automatic manner with the same

Figure 11. Convergence of far-field mirrors to near-field mirrors 🅭🅯 CC BY 4.0

Figure 12. Mirrors for collimated (top) and point (bottom) light; visibility cells (left), component mesh (middle) and rendering with LuxRender (right) 🅭🅯 CC BY 4.0