AN ENTROPIC OPTIMAL TRANSPORT NUMERICAL APPROACH TO THE REFLECTOR PROBLEM

A BSTRACT . The point source far ﬁeld reﬂector design problem is one of the main classic optimal transport problems with a non-euclidean displacement cost [Wang, 2004] [Glimm and Oliker, 2003]. This work describes the use of Entropic Optimal Transport and the associated Sinkhorn algorithm [Cuturi, 2013] to solve it numerically. As the reﬂector modelling is based on the Kantorovich potentials, several questions arise. First, on the convergence of the discrete entropic approximation and here we follow the recent work of [Berman, 2017] and in particular the imposed discretization requirements therein. Secondly, the correction of the Entropic bias induced by the Entropic OT, as discussed in particular in [Ramdas et al., 2017] [Genevay et al., 2018] [Feydy et al., 2018], is another important tool to achieve reasonable results. The paper reviews the necessary mathematical and numerical tools needed to produce and discuss the obtained numerical results. We ﬁnd that Sinkhorn algorithm may be adapted, at least in simple academic cases, to the resolution of the far ﬁeld reﬂector problem. Sinkhorn canonical extension to continuous potentials is needed to generate continuous reﬂector approximations. The use of Sinkhorn divergences [Feydy et al., 2018] is useful to mitigate the entropic bias.

it can be regarded as a point in space.It can therefore be modelled as a probability distribution on the sphere, it will be denoted µ in this paper.The light hits a perfect mirror and we are also given a desired target light distribution, the "illumination" in the far field.From the far field the reflecting surface can be regarded as a point and the illumination again modelled as a probability distribution, denoted ν, on the sphere.Total light conservation is assumed.The design problem is to determine the shape of the mirror which produces the specular reflection from the source to the target distribution.This can be interpreted as the inverse problem of generating some illumination given an illuminance and a reflector (see figure 1.1).
This problem has an elegant mathematical modelization and solution based on the optimal transportation (OT) theory due to [Glimm and Oliker, 2003] and [Wang, 2004].We briefly recall the main result as presented in [Wang, 2004].In its Kantorovich primal and dual form (see [Villani, 2008]) : Theorem 1 (Kantorovich duality).Given two compact manifold X and Y endowed with a continuous, bounded from below cost function c : X × X → R and two borel probability measures (µ, ν) ∈ P (X) × P (Y).Then, Kantorovich problem in primal and dual forms (1.1) has solutions. (1.1) OT(µ, ν) := min with respectively primal : and dual : constraints sets The notation f , α Ω stands for the duality product Ω f dα between bounded continuous functions f ∈ C(Ω) and probability measures α ∈ P (Ω), { f ⊕ g}(x, y) = f (x) + g(y) is the direct sum and µ ⊗ ν ∈ P (X × Y) the tensor product.Finally 1 Ω is the characteristic function, i.e. a constant 1 on Ω.
Under suitable hypothesis c (as they are technical and satisfied for the costs in this paper, we skip this part), the OT problem is well posed and the the optimal transference plan γ is concentrated on a graph of the OT map y = T(x) implicitely defined by the saturation of the dual constraint : (1.2) f (x) + g(T(x)) = c(x, T(x)), µ a.e.
The pair ( f , g) are called the Kantorovich potentials and is unique up to a constant .
By construction T is a measure preserving map characterizing the transport.The measure preserving property is usually denoted ν = T#µ (T pushes forward µ to ν).The pushforward of µ is the measure defined as ) for all ν measurable subset A Remark 1 (L p Wasserstein metric).For complete separable metric space X and L p costs c := 1/pd p (x, y), this OT problems defines a separable metric on the set of probability measures with finite second moments: the "Wasserstein" distance, which is given by W p p (µ, ν) := OT(µ, ν).This metric metrizes weak convergence of measures is a fundamential tool in image processing (see [Peyré and Cuturi, 2018]).
In [Wang, 2004], Wang showed that the point source reflector model can be translated to an OT problem.More precisely, he proved the following theorem : Theorem 2. Let S 0 ∈ S d−1 and S ∞ ∈ S d−1 be connected domains in northern and southern hemispheres respectively, µ and ν which represent the given illuminance and illumination probability distributions.Then theorem 1 applies to the cost function A transport map T satisfying (1.2) exists and the solution of the corresponding OT problem can be used to build the desired reflector.
The construction of the reflector can be summarized as follows : Taking the exponential of the dual constraints and the saturation property (1.2) we get We now define in R d a family of parabolic reflectors with axis y ∈ S ∞ : x ∈ S 0 → P y (x) := e −g(y) 1 − x • y .
And directly infer that the reflector shape parameterized over the directions in S 0 and given as : Under this choice the map x → T(x) can be interpreted as the specular reflection of an optical ray at R(x) onto a parabola of axis T(x) while the illumination and illuminance constraints are enforced by (1.3).

Numerical resolution of (1).
A good overview of the available numerical methods applying to the general problem (1) are described in [Peyré and Cuturi, 2018].We summarize the general classes of methods which can be used, in particular, to tackle the reflector cost (1.4).
Linear Programming.The first and obvious one, as already suggested in [Wang, 2004] is the linear programming approach.This assumes that the data in discrete form is given by Of course the number of discrete points for µ and ν may differ, we keep N for both to simplify the presentation.
(1.8) is a discrete linear programming problem which can be solved numerically using standard linear programming solvers.The main drawback of this method is it's high dimensionality.It is a linear problem with N × N unknowns and 2 N constraints.Numerical resolution with linear solvers which have cubic complexity in the number of unknowns is therefore out of reach for reasonable discretizations (typically N > 100).For more details see [Peyré and Cuturi, 2018].
A Partial Differential Approach.At least formally, taking the gradient of (1.2) gives the OT map as a function of f : (1.10) Assuming sufficient smoothness to interpret (1.3) in a pointwise sense and plugging this expression gives a Monge-Ampère PDE for L 2 cost.This formulation has been studied in detail by Urbas [Urbas, 1997].For the reflector cost (1.4) the PDE (see equation (1.5) in [Wang, 2004]) is even more non-linear and there is no available numerical approximation theory.A B-spline collocation method was nevertheless proposed and implemented in [Brix et al., 2014], [Brix et al., 2015] with convincing numerical results.An alternative solution method of the Monge-Ampere problem based on a least square approach is presented [Romijn et al., 2020] [Romijn et al., 2019].Wu [Wu et al., 2013] derives the Monge-Ampere equation for a lens surface and solve the equations using standard finite differences and Newton iteration.
Semi-discrete OT.Semi-discrete OT deals with the special case where µ has continuous density and the other ν N is discrete as in (1.7).Then the dual formulation in 1) can be reduced to discrete optimisation problem through the elimination of the constraints replacing f by its c-transform : and optimizing over the vector {g(y j )} j=1..N .For the L 2 cost, It is now well understood that a Newton method can be used [Mérigot, 2011].The implementation relies on fast (linear time) computations of a tesselation of the target domain called Laguerre cells : . The method was applied to the reflector cost (1.4) in [Machado Manhães De Castro et al., 2016] but the Laguerre cells can be only computed in R d and not on S d−1 thus lifting by one the dimension of the problem.Solvers in dimension 2 and 3 are available [Leclerc, 2019] [Lévy, 2015].See also the pioneering work of [Kochengin and Oliker, 2003] using this methods on the reflector problem.1.3.Our contribution.Our contribution is based on a strict convexification of the OT problem called "entropic regularization".In a nutshell : It has been introduced for OT computations in [Cuturi, 2013] (see [Peyré and Cuturi, 2018] for a comprehensive review).This topic has been very active since for two reasons : first, the method comes with a simple computational method call Sinkhorn algorithm for which there are now countless variants and acceleration techniques.It can be applied to discretization on the order of N = 10 6 .Secondly, it can be applied to a wide range of costs but does not seem to have been investigated yet for the reflector cost.The "entropic regularization" of the Kantorovich problem (1) is based on the following KullBack-Leibler divergence or "relative entropy" (KL) penalization : (1.12) where ε > 0 is a small "temperature" parameter (see [Léonard, 2013] for a Statistical Physics interpretation of this problem due to Schroedinger) and ) dγ if γ is absolutely continuous w.r.t to µ ⊗ ν and +∞ else.
The primal-dual optimality condition is given by (1.13) The optimal entropic plan is therefore the scaling by the Kantorovich potentials of a fixed Kernel exp(− 1 ε c).It is diffuse, i.e. not concentrated on a map, ε can be interpreted as a bandwidth under which transport is blurred.
Numerical solutions are produced under the discretization (1.7) of this problem, i.e. replacing (X, Y , c , µ, ν) by (X N , Y N , c N , µ N ν N ) in the above formula.The resulting discrete non-linear system of optimality is solvable in particular with Sinkhorn algorithm.
In section 2 we recall basic convergence properties of this approximation when ε → 0 and how N and ε are connected.We also discuss the resulting bias when ε > 0 and how the notion of Sinkhorn Divergence introduced in [Ramdas et al., 2017] [Genevay et al., 2018] [Feydy et al., 2018] can produce a correction.
We then adapt this framework and Sinkhorn algorithm in section 3 to the reflector cost OT problem and discuss specifics that need to be addressed for this transition, such as de-biasing for reflector cost, choice of discretization and interpolation.We present several example of re-simulations using ray tracing.
Wasserstein distances is a natural distance to measure discrepancies between measures.In section 4 we therefore propose to evaluate it between the illumination ν and a sampling of the computed approximation.The sampling is obtained via ray tracing reflections on a continuous interpolation of the reflector as a residual error.We use this setting to study the convergence with N and ε.We also show how the continuous approximation of the reflector behaves when the number of rays used for the resimulation increase.

ENTROPIC OT
2.1.Entropic OT and Sinkhorn algorithm.We will restrict from now on to the discrete case (see (1.7)) and the dual optimal transport problem : where we use the same notation ( f ε , g ε ) for discrete vectors in R N .We solve (2.1) with Sinkhorn algorithm.It corresponds to a block coordinate ( f ε and g ε ) ascent : Initialize with g 0 ε = 0 Y and then iterate (in k) : (2.2) Sinkhorn linear convergence for fixed ε has been largely studied (see [Peyré and Cuturi, 2018]) again).Its numerical stability depends on typical transport scale of the data τ = sup x∈X c(x, T(x)) and ε.The computer hits overflows or underflows when the ratio τ/ε is too large.A good review of existing hacks and methods to mitigate this problem can be found in [Schmitzer, 2016].
The potentials computed using the iterations (2.2) are defined on the discrete sets X N and Y N but they admit a canonical extention by replacing c N (x i , y j ) respectively by c(x, y j ), x ∈ X and c(x i , y), y ∈ Y. Omitting the iteration index : As our goal is to form the reflector (1.6) from the discrete numerical potential f ε , we are interested in the convergence of f ε as ε → 0 and N → ∞.To the best of our knowledge the joint convergence in N and ε has only be studied in [Berman, 2017], we reproduce partially his results : Theorem 3 (Berman joint convergence -corollary 1.3 [Berman, 2017] ).We assume µ and ν are in C 2,α and positive, and that N and ε are dependent parameters : N = (1/ε) d where d is the dimension of the problem.A technical condition on the sequence of discretization (X N , Y N , c N , µ N ν N ) called "density propery" (see remark 3 below) is also necessary.Then there exists a positive constant A 0 such that for any A > A 0 the folowing holds : setting for some constant C (depending on A) and f is an optimal potential for (1).
The following remarks can also be extracted from [Berman, 2017] : Remark 2 (Domains and Cost ).Theorem 3 holds for compact manifold X and Y and in particular the Euclidean torus endowed with the periodic L 2 Euclidean costs.But it also on the sphere for the reflector cost (see section 6.3.3 [Berman, 2017] ).
Remark 3 ( Density property Lemma 3.1 [Berman, 2017] ).For any given open set U intersecting the support X of µ (same for Y and ν) For the flat space X ⊂ R d , this condition is enough.For curved surfaces, a technical generalization is required.But in both cases, this density property ensures the discretization of X and µ (1.7) is such that, for U the sequence of approximations µ N (U) never converges faster to 0 than ε (remember that N = (1/ε) d ).
For X = R d , uniform grids are fine.See section 3.2 for the reflector case.
Remark 4 (Domains and Cost ).Theorem 3 holds for compa particular for the torus and the L 2 Euclidean costs but also on the sphere for the reflector cost (see section 6.3.3 [Berman, 2017] ).
Remark 5 (Transportation plan convergence).The support of the entropic transportation plan γ ε (see (1.13)), built trough the interpolation procedure in theorem (3), converges exponentially fast to the graph of the non entropic OT map when ε → 0.
Remark 6 (Gradient convergence).The convergence rate in (2.4) is sublinear in ε.As ε is the "edge" length of the discretization, discrete approximations of the derivative of f may be non convergent.

Entropic bias and Sinkhorn Divergence.
From now on, the potentials computed using Sinkhorn iterations (2.2) are denoted by f OT ε and g OT ε .
Entropic OT is popular to approximate transport distances between probability measures.The Entropic cost OT ε (µ, ν) is know to converge exponentially fast to OT(µ, ν) = OT 0 (µ, ν) with ε [Cominetti and Martin, 1994].Regarding the reflector problem, we are not interested in the transport cost but in f OT ε itself, the approximation of f .As can be seen formally on the dual formulations of ( 1) and (1.12) the gradient of the transport and entropic transport costs as a function of measures µ and ν are given by the Kantorovich potentials : Proposition 2 [Feydy et al., 2018]).The convergence rate in theorem 3 with an infinite slope at ε = 0 and the numerical stability limit imposed when decreasing ε is a difficulty.The pollution of the potential by the entropy induced bias is well known issue.
It arises for instance when using the L 2 Wasserstein distance as a loss over the space of probability measures.Using the easier to compute f ε as a gradient proxy for f to minimize µ → OT(µ, ν) will therefore be biased by the entropy.This problem is discussed in depth in [Feydy et al., 2018] where it is proposed, in order to correct the bias, to add "diagonal terms" to correct the entropic cost : (2.5) Quite remarkably, the authors show that this quantity, called Sinkhorn divergence, remains positive and is convex.It also obviously vanishes for µ = ν wich is not the case for OT ε .Thanks to the symmetry, there is only one dual potentials for each of diagonal problems.We denote them f µ ε and f ν ε .They can be computed using independent Sinkhorn iterations : (2.6) The µ gradient of S ε , denoted f S ε may be formed by a simple substraction.An open question is whether this correction : (2.7) is a better approximation to f than f OT ε .Numerical simulations of gradient flows in [Feydy et al., 2018] do indicate this is the case.
Remark 7 (Asymptotics of OT ε ).The entropic bias may be related formally to asymptotic results on the difference between OT and OT ε , see [Conforti andTamanini, 2019] [Pal, 2019] for recent publications on this subject.For the L 2 cost (Theorem 1 [Pal, 2019] treats more costs in the form c := g(x − y), g convex which includes the 1D reflector cost), it gives the following assymptotic behavior for small ε : This result is for continuous and compactly supported densities µ and ν with respect to L(R d ), the Lebesgue measure.Formally taking the gradient in µ of (2.8) we get for small ε : In a smooth continuous setting and after convergence of Sinkhorn (2.6) we get ) has maximum at x ′ = x.This is expected to be asymptotically true as ε → 0 (remark that T(x) = x in (1.2) for µ = ν).Then, the Laplace method (also used in [Berman, 2017] and [Pal, 2019]) gives for ε → 0 for some constant C depending on x and the hessian of f µ OT ε .It gives : which shows that the Sinkhorn divergence de-biasing (2.7) is consistent with (2.9).This correction can also be observed on the numerical simulation presented in section 3.1.Remark finally that the density property discussed in remark 3 is requiring for the first order bias terms to converge to 0 when discrete mesures µ N converge to original measure µ.

APPLICATION OF ENTROPIC REGULARIZATION TO THE REFLECTOR PROBLEM
3.1.Entropic Bias and de-biasing.The Entropic bias may be observed on numerical solutions.In the simplest identity reflector case for example, where each ray is reflected in its opposite direction, i.e. when ν(x) = µ(x + π), x ∈ S, the exact potentials f and g are constant and the reflector R = {x e f (x) |x ∈ X} is a portion of circle.It is represented, for the plane problem (d = 2), in dashed line in figure 3.1.The solid line is the entropic reflector R ε = {x e f OTε (x) |x ∈ X} computed with Sinkhorn algorithm.The discretization is N = 128 and ε = 0.1 is chosen large enough to see the bias.Remember that the potentials are defined up to a constant, we therefore adjust the plot such that the reflectors superimpose at x = 0.The illuminance and illumination µ and ν are perfect centered Gaussians and represented respectively in blue and red.This example shows the effect of the entropic bias (2.9).We also trace rays : the bottom and middle rays are reflected if the exact opposite "identity" direction.The top rays, reflecting on the biased entropic reflector will send light strictly inside the support of ν.This shrinking effect of the re-simulation offers a former explanation of the third plot of figure (5.1)where the density is a characteristic function which decreases to 0 at the boundary.
In order to deal with this bias, the notion of Sinkhorn divergences (2.5) is required.When the spaces X and Y are different and a for a general cost c(x, y), we suggest the following extension : where µ ′ = argminOT(µ, •) and ν ′ = argminOT (•, ν).This is consistent with the distance costs c(x, y) = d p p (x, y) as µ ′ = µ and ν ′ = ν in this case.Simillarly, in the reflector cost case , µ ′ and ν ′ correspond to the reflections of respectively µ and ν on the circular reflector mentioned above.Indeed one has arg min y∈Y −log(1− < x, y >) = −x.
The de-biased Kantorovich potentials for the Sinkhorn divergence follow the same formula as for L 2 -distance cost : 3.2.Choice of Discretization.In theorem 3, the requirement that ǫ is of order (1/N) 1 d and that µ N , ν N satisfy the density condition (remark 3) are closely related.Intuitively this condition means that while choosing the discretizations µ N , ν N with N points, it is important to make sure that they approximate corresponding distributions µ, ν with integration error of order (1/N) 1 d for functions which do not oscillate on finer scale then (1/N) 1 d .This density property is satisfied in particular by Quasi Monte-Carlo point clouds ([Berman, 2017]), defined by bounding the worst case error of integration over a desired function space W : Definition 1.Given a Riemmanian manifold X of dimension d, a corresponding normalized volume form dV, and a discretization X N = x i{i≤N} , the worst case error of X N with respect to function space W is: Then, X N = x i{i≤N} is a Quasi Monte-Carlo discretization of X on lengthscale (1/N) Where W s p is a Sobolev space of functions f such that all derivatives of order s are in L p (X) and C s,p is a uniform constant depending only on s and p. [Berman, 2017] shows that for a density µ which is absolutely continuous with respect to the volume form dV with density ρ µ , the discrete approximation µ N which satisfies the density requirements of theorem 3 can be constructed using Quasi Monte-Carlo discretization X N of X as the empirical measure : When the space X is linear, a standard uniform square grids with step-size h is a Quasi Monte-Carlo systems on lengthscale h.But for curved spaces, such as a sphere in the case of the reflector problem, constructing a Quasi Monte-carlo systems is not straightforwad and usually they do not have such simple structure.
One way of constructing such a discretization on the sphere, which was used for the simulations in this work, is given in [Womersley, 2017].

Ray tracing, Re-Simulation and Interpolation.
The reflector problem is the inverse problem associated to the following forward problem : Given the source illumination µ as a distribution on the unit sphere S d and the reflecting surface R as a radial graph over S d , what is the distribution of illuminance after µ is reflected from R? The geometric optics answer to this problem is the simulation of the light propagation known as Ray-Tracing.It can be simply described as follows : The distribution µ over the sphere is sampled with an empirical measure µ M := ∑ i=1..M p i δ x i .Then M rays are shot according to their direction {x i }s and reflected in the directions {y i = x i − 2 (x i • n i )n i }s according to Snells's law (n i is the normal vector to R at R(x i )).The resulting directions provide a sample ν M := ∑ i=1..M p i δ y i of the illumination distribution in the far field ν = T # µ.In the optics community, ray tracing is used to evaluate visually the quality of the computed reflector.This procedure is called Re-Simulation.The quality of course also depends on the sampling method and the number of samples M.
Remark 8 (Ray tracing and Push-Forward).Since by construction, the exact continuous reflector reflects the ray x in the direction T(x) (the associated OT map), the above described procedure is the same as computing the empirical push-forward ν M = T # µ M and this is a sampling of ν = T # µ.
In order to use this procedure with our numerical approximation of the reflector, we need a continuous extension of the reflector from its discrete version computed on the N discretization (3.5).For standard cartesian or structured grids this is not a problem, since there are various wellstudied interpolation techniques.Our solutions on the contrary are computed on an unstructured Quasi Monte-Carlo grid.We use the Optimal transport and entropic OT canonical continuous extension procedures based on duality and using the g potentials (respectively (1.11) and (2.3)).
Definition 2 (c-concave interpolation).For a discrete pair of Kantorovich potentials ( f N , g N ), defined on the discrete grids X N and Y N respectively, the c-convex interpolation of the first potential f N is (3.6) f (x) := min As discussed in section 1.1 (see (1.5)), the continuous reflector R = {xe f (x) |x ∈ X}. is the envelope of N paraboloids with focal directions {y j }.This means that this approach is able to generate a continuous reflector with discontinuous derivatives.This may be an advantage when the target distribution is discontinuous itself.The downside is that in order to achieve continuous effect for human eye, a very fine discretization (a large N) is needed.
The second approach uses either the smoothed transform (2.3) for OT ǫ or (2.10) for S ǫ Definition 3 (Entropic interpolation).For a discrete pair of Kantorovich potentials ( f OT ε ,N , g OT ε ,N ), defined respectively on the discrete grids X N and Y N and solutions of the Sinkhorn algorithm applied to the OT ǫ problem, the extension of the first potential f OT ε ,N is Note that (3.7) is the "softmin" version of the mininization (3.6).
Finally using the "'identity" potentials ( ) computed when applying Sinkhorn algorithm applied to S ǫ , the smooth interpolation is : This approach would result in a smoother reflector than in the previous case but will also blur any sharp details in the target image.
3.4.Summary of the numerical procedure.Our input are analytical descriptions of the illumination/illuminance µ and ν described in the test cases section below.All test cases presented in this paper will have the same source and target domains X and Y.The source domain X ⊂ S 2 will be the inverse stereograpic projection in the northern hemisphere of the square domain centered at the origin {(x 1 , x 2 ) ∈ R 2 | − 0.6 ≤ x 1 ≤ 0.6, −0.6 ≤ x 2 ≤ 0.6}.Similarly, Y ⊂ S 2 will be the inverse stereographic projection in the southern hemisphere of same domain.The outputs are resimulations/RT computed according to the following procedure : (1) Computation of f N : The potentials, optimal maps and reflectors obtained from the minimization of the regularized functional (2.1) will be denoted using the subscript • OT ε , and those obtained from the minimization of the de-biased functional (3.1) by subscript • S ε .The discrete Kantorovich potentials are computed for the discretizations (µ N , ν N ) induced by the Quasi Monte-Carlo grids X N and Y N respectively according to the formula (3.5).
(2) Interpolation/ray tracing: The choice of interpolation when building the continuous reflector for RT will be denoted with the upperscript • when using the c-concave interpolation (3.6), and • when using the Entropic interpolation (3.7,3.8).The ray tracing/re-simulation will be shot using a Quasi Monte-Carlo sampling of the source domain of size M.
(3) Domain of computation: The distribution of reflected directions (obtained as distribution of directions in R 3 ) is projected using the stereographic projection from the south pole to the equator plane: Definition 4 ( Parameters and notation for the forward map).In order to discuss the numerics we need to introduce a notation for the forward map induced by the above procedure.This is cumbersome as many parameters are involved : ε the entropic regularization, N the discretization of the distributions µ and ν used in the Sinkhorn algorithm.The choice of using or not the Sinkhorn divergence correction.The choice and notation of the interpolation for the reflector has discussed above and in section 3.3.We will use the notation RN OT ε , RN OT ε , RN S ε , RN S ε for the reflectors defined using (1.6) and respectively f We will denote ν M = RN OT ε [µ M ] the ray traced sampling with M rays (and a similar notation for the other versions of the reflector).Of course ν M also depends on the choice of the method and other parameters.Finally the N and M discretizations of the domain are QMC discretizations as explained in section 3.2.
3.5.Numerical illustration.We give here several re-simulations as an illustration of the different approximation of the reflectors described in section 3.3.See section 3.4 and 4 for a precise description of the numerical procedure.
The test cases we used are described below and are set on bounded subsets of the sphere.We are not exactly in the setting of theorem 3 and remark (2).
Test Case 1: Square To Circle.The source distribution µ will be the uniform distribution over the set with a square stereographic projection The target distribution ν in the first example will be an uniform distribution over the set with a disck (circle) stereographic projection Even though the densities are constant, mapping from a non smooth support geometry of the square to the smooth geometry of a circle is not a trivial task.
We show in figure 5.1, from left to right and from top to bottom the re-simulations using the reflectors ROT ε , ROT ε , RS ε , RS ε and finally also the QMC discretization used for ν N , N = 64 * 64 points.The regularization parameter was taken to be ε = 1 2 * 64 .Figures 5.1 (a) and (c) correspond to the c-concave interpolations.This builds the minimal envelope of the family of parabolae with focal axis given by the ν discretization.All rays hitting one parabola end up at the same position.So up to numerical errors we indeed recover (e), that is ν N even though M >> N. Also the S ε solution performs better at the boundary.
Figures 5.1 (b) and (d) correspond to the entropic interpolations.This is a smoothing of the the c-concave interpolation.The M rays are therefore distributed more evenly over the support of the target.We observe, for (b), at the boundary of the support, an important shrinking effect.It can be explained by the OT ε entropic bias as in section 3.1, which smoothing the density at the boundary of the support would produce, according to formula (2.9) and the effect observed in figure 3.1.As expected, the S ε (d) solution is effective to de-bias the solution.Except, at the corners of the square where the the map is singular, we see extra artefacts probably induced by the smooth interpolation.
Test Case 2 : Square to Gaussian.The target distribution ν is a gaussian distribution with density on the projected domain ρ(x 1 , x 2 ) = e − x 2 1 +x 2 2 2 over whole target domain Y, and the source distribution will be the same as in Test Case 1.
Test Case 3 : Square to Square.The source and target distribution is a uniform distribution over the set with stereographic projection (square) Test Case 4 : Circle To Steep Gaussian.The source distribution µ is uniform distribution over the set with stereographic projection (circle) The target distribution ν is gaussian distribution with density on the projected domain ρ(x 1 , x 2 ) = e −16 * (x 2 1 +x 2 2 ) over whole target domain Y. Test Case 5 : Circle To Two Steep Gaussians.The source distribution µ is the same as in the fourth.The target distribution ν is a gaussian distribution with density on the projected domain ρ(x 1 , x 2 ) = e −16 * ((x 1 −0.25) 2 +(x 2 −0.25) 2 ) + e −16 * ((x 1 −0.25) 2 +(x 2 +0.25) 2 ) over whole target domain Y.

ERROR ESTIMATIONS AND CONVERGENCE
4.1.Wasserstein metrics as error estimator.From the OT point of view, a straightforward evaluation of our numerical approach would be to build an approximate transport map T app obtained by plugging the different approximations of the potential f 3) and compare it against the exact solution T. For the reflector problem, however, only trivial analytical solutions such as the circle (identity) or a parabola (dirac target) are known.On the other hand and from the optics application point of view, the main concern is not the shape of reflector itself (except for designing constraints ignored here), but rather the quality of the re-simulation.As mentioned in remark 8, mathematically this re-simulation is a sampling of the push-forward T app # µ.Therefore, in order to build an error estimator we will consider the difference between the push-forward T app # µ and the desired distribution ν.For this we will use the standard L 2 Wasserstein distance W 2 (T app # µ, ν), or its entropic counterparts presented in section 2. They are known to provide a smooth estimators between on empirical distributions [Peyré and Cuturi, 2018].In order to build a numerically computable error estimate we will use the ray traced sample µ M of the source measure µ.
First, using the triangle inequality for the Wasserstein metrics one gets : The first term on the right hand side (W 2 (T app # µ, T app # µ M )) can be estimated using Lemma 4 and inequality (5.3) in the appendix 5 giving : where K depends on the data of the problem.The convergence in W 2 norm of such sampling has been studied in [Fournier and Guillin, 2015] and is known to behave asymptotically as M − 1 2 in expectation.
The second term on the right hand side of (4.1) can be approximated using νM = RN OT ε [µ M ] the sample obtained using the different combinations of ray tracing and interpolation (we only show one of these combinations there see definition 4 and remark 8) : The continuous densities µ and ν still appear in estimates (4.2-4.3).The W 2 distance can be computed either using semi-discrete OT (section 1.2) which relies on a P1 discretization of the continuous densities or again using Sinkhorn divergence but this time for the L 2 cost.This is the choice we made in this paper and µ and ν are discretized on a finer M ∞ = 512 2 grid.
The error estimates are therefore computed using the Sinkhorn divergence approximation of the L 2 Wasserstein distance on the projection plane with a smaller ε : W 2 (., .)≃ S ε=1e−06 (., .),but we will keep the W 2 notation below.

Numerical convergence Study in N.
In this section M = 128 2 the number of rays for the resimulation is fixed.The densities µ and ν are discretized with a finer M ∞ = 512 2 points orthogonal grid on the plane.
In (4.2), for µ as in Test Cases 1,2 and 3 , we obtained with the above parameters : For Test cases 4 and 5, we obtained 2), we plot for different tests cases described below the convergence curves in N for the different reflector approximations methods explained in definition 4 and two values of ε : 1/2 √ N, 1/8 √ N. See figures 5.7-5.11for the 5 test cases.We observe : (1) Convergence to error levels comparable to (4.4-4.5) which makes sense as ν M is at best a QMC sampling of ν.
(3) Sinkhorn divergence S ε de-besiasing is effective, both in the sense of obtaining lower errors, and being less dependent on the choice of ε. (4) As the target is smooth,the entropic interpolation of S ǫ solutions are less dependent on the choice of the discretization, and moderate values of N are enough to achieve the same error as when using the highest value of N used here.

Numerical convergence Study in M.
In this section, we use the Entropic interpolation method and generate the potential with the Sinkhorn divergence method as they seemed to performed the best.We then fix N = 128 2 and ǫ = 128/8 and study the dependence of the error term W In figures 5.12-5.13we plot the error curves for Test Cases 1 and 5 in original and log scales.The curves demonstrates that the computed continuous numerical approximation of the reflector preserves after ray tracing the quality of the illumination/source ray sampling.In particular : (1) The convergence curves closely are similar to the convergence curves in N obtained using the Sinkhorn divergence method with ǫ = 1/(8 √ N) and the c-concave interpolation.Indeed, with this interpolation method and a good approximation of the potential, the reflector will send all the rays onto ν N wich is also discretized using a QMC system (see also section 3.3).So increasing N there and increasing M here results in the same empirical measure ν (2) In logarithmic scales the curves agree with the M − 1 2 rate predicted by the theory [Fournier and Guillin, 2015].

CONCLUSION
We have studied numerically the use of entropic regularization to solve the OT/Reflector problem.Based on theorem 3 [Berman, 2017] we checked the convergence on the method the entropic regularization and space discretization parameter.Error quantities on the reflector shape itself are not available.We therefore defined and used an error criterium based on re-simulations/ray tracing and the Wasserstein 2 distance, a standard tool to compare images and empirical distribution.We also tested a novel approach [Feydy et al., 2018] called Sinkhorn Divergence to correct the bias induced by the entropic regularization.
Our experiments indicate first that the method converges.It also shows that the combination of a QMC discretization satisfying the density property (remark 3) the entropic interpolation (2.1) and correcting the Entropic bias with the Sinkhorn divergence technique produce the best results.W p (T # µ 1 , T # µ 2 ) ≤ Lip(T)W p (µ 1 , µ 2 ) Proof.Let ( f , g) be optimal pair of Kantorovich potentials for W p (T # µ 1 , T # µ 2 ).Then, for all x, x ′ ∈ X When the target support Y is convex and the densities of µ and ν smooth enough (C 1,α is enough) and bounded below and above by positive constants, a classic regularity result by Caffarelli [Caffarelli, 1992] gives sufficient regularity to get (5.3)Lip(T) ≤ K where the constant K only depends on the dimension d and the the data µ and ν.See [Philippis and Figalli, 2014] for further refinements and discussions.
FIGURE 1.1.Reflector problem from Point source O to Far Field.

FIGURE 3 . 1 .
FIGURE 3.1.Polar plot of a 1D computation.In blue and red are µ and ν.Solid and dashed black curves are R ε and R. Also shown are rays traced.The bottom and middle rays are reflected if the exaxt opposite "identity" directon.The top rays, reflecting on the biased entropic reflector will send light strictly inside the support of ν.

Figures 5
Figures 5.2-5.6 show the re-simulations for the 5 test cases.

FIGURES
FIGURES FIGURE 5.1.Test Case 1 : From left to right and from top to bottom the resimulations of µ M (M=128*128) using the reflectors ROT ε , ROT ε , RS ε , RS ε (check definition 4 for explanation of the notations) and finally also the QMC discretization used for ν N , N = 64 * 64 points.The regularization parameter was taken to be ε = 1 2 * 64 for all four solutions.