Iterative scheme for solving optimal transportation problems arising in reflector design

We consider the geometric optics problem of finding a system of two reflectors that transform a spherical wavefront into a beam of parallel rays with prescribed intensity distribution. Using techniques from optimal transportation theory, it has been shown before that this problem is equivalent to an infinite dimensional linear programming (LP) problem. We investigate techniques for constructing the two reflectors numerically. A straightforward discretization of this problem has the disadvantage that the number of constraints increases rapidly with the mesh size. So with this technique only very coarse meshes are practical. To address this well-known issue we propose an iterative solution scheme. In each step an LP problem is solved. Information from the previous iteration step is used to reduce the number of constraints necessary. As a proof of concept we apply our proposed scheme to solve a problem with synthetic data. We give evidence that the scheme converges. We also show that it allows for much finer meshes than a simple discretization scheme. There exists a growing literature for the application of optimal transportation theory to other beam shaping problems. Our proposed scheme is easy to adapt for these problems as well.


Introduction
The following beam shaping problem from geometric optics was described in [15]; see Figure 1: Suppose a point source emits a spherical wavefront with a given intensity distribution. The problem at hand consists of transforming this input beam into an output beam of parallel light rays with a prescribed intensity distribution. This transformation is to be achieved with a system of two reflectors. The problem has some practical importance in engineering, see further literature cited in [15].
The first rigorous mathematical solution to the problem was provided in [10], using an approach based on the theory of optimal transportation [5,6,8,19]. See also the references [11,9,21,12,16] which deal with other beam shaping problems using related techniques.
The result of [10] is summarized in section 2 below, with Theorem 2.5 stating the main result. The central feature is that the original reflector design problem is reformulated as an infinite dimensional constrained optimization problem, namely the problem of minimizing a certain functional on a function space. It is the dual problem for the problem of finding a map from the input aperture to the output aperture which minimizes a certain cost functional. (See Problem II in section 2 for the exact statement.) This reformulation of the problem is not only of theoretical value for questions of existence and uniqueness of solutions, but it also translates into a practical method for finding the solution. In fact, the discretization of the infinite dimensional constrained optimization problem is a standard linear programming problem and can be solved numerically. In this paper, we describe a numerical scheme for solving these problems. An immediate obstacle is that the number of constraints becomes very large even for relatively coarse meshes -if there are M sample points in the input aperture (denoted byD in Figure 1) and N points in the output aperture (T in Figure 1), then the number of constraints is M · N . We found that with a simple straightforward discretization method, standard LP solvers (on PCs with up to 4 GB RAM) could not handle more than approximately 500 sample points on the domainsT andD -a number which is arguably too small for many applications. For this reason, we devised a more elaborate iterative method, where discretized systems are solved on finer and finer grids, in each step using information from the previous solution to choose only a subset of all possible constraints. This slows the growth of the number of constraints. (Details are described in Section 3.) With this step-wise mesh refinement scheme, as a proof of concept, we were able to obtain solutions on meshes with about 4500 points on each aperture using MATLAB and a standard LP solver. Using compiled computer language and more specialized solvers for optimal transportation problems, we expect that this number could still be substantially improved.
There has been increased interest in numerical methods for optimal transportation problems in the last 10 years. Most work has concentrated on the Monge-Ampère equation, which arises in the special case of optimal transport in R n with quadratic cost (also called L 2 optimal transport) [14,2,3,7,13,23], although some authors have treated costs proportional to the distance (L 1 optimal transport) [1]. These methods are based on fluid mechanical approaches [2] or various finite difference approaches [3]. They are generally faster and allow for much larger mesh sizes than methods based on a discretization of the linear programming problem, but since they use the special structure of the Monge-Ampère equation, they are not directly applicable to more general cost functions or more general manifolds.
One of the few papers that numerically addressed more general situations for optimal transport is [18]. A similar approach to the solution was taken there as we do in this paper, namely a discretization of the linear programming problem. (The authors in [18] used the primal formulation of the optimal transportation problem, whereas our approach is equivalent to using the dual formulation.) The authors noted that the number of variables is growing very fast with mesh sizes, making it impossible to solve the problem numerically even for moderate mesh sizes due to memory limitations. In fact, they noted that the software they used, the package Soplex [22], was designed to handle up to 2 million variables. This corresponds to mesh sizes of approximately 700 points on the input and output aperture in our problem if one employs a straightforward discretization scheme. The authors noted that for better results, one needs carefully designed programs. This is what we supply here for our problem: As explained above, our "mesh refinement" approach addresses this problem and makes it feasible to solve the problem for finer meshes.
We note that our proposed algorithm does not make any a priori symmetry assumptions on the form the reflectors. It can also very easily be adapted for the numerical solution of other beam shaping problems for which a formulation using optimal transportation theory have been found, for example those in [9,11,20].
This article is organized as follows: In section 2, we recall the reformulation of the reflector construction problem as a linear programming transportation as given in [10] and fix some notation. We then describe a basic discretization scheme for the numerical solution and propose Figure 1: Sketch of the reflector problem. The point source is located at the origin 0 of the coordinate system. Note the conventions on coordinates as illustrated by the coordinate system in the lower left hand corner: The output beam propagates in the direction of the negative z−axis, and points in the plane perpendicular to the z−axis are denoted by the vector x ∈ R 2 . Thus points in three dimensional space are denoted by (x, z). See Section 2 for more. (Adapted from [10].) an improved "mesh refinement" scheme in section 3. The following section 4 is devoted to numerical tests of this scheme. We first derive an explicit analytical solution to be used as synthetic data for our numerical work. Then we compute the solution numerically and analyze the error of approximation. Our scheme requires a sequence ε (1) , ε (2) , . . . of "constraint inclusion thresholds." We analyze different choices of this sequence. We conclude with a brief summary and propose future work.

Formulation as Linear Programming Problem
We now briefly review the notation for the problem posed in [15], as well as the result of [10], which reformulates the reflector design problem as a linear programming problem.
Consider the configuration show in Figure 1. A point source at the origin O = (0, 0, 0) generates a spherical wave front over a given input apertureD contained in the unit sphere S 2 . The input beam has a given intensity distribution. By means of two reflectors, this wave front is to be transformed into a beam of parallel rays propagating in the direction of the negative z−axis. This output beam has to have a prescribed intensity distribution. The cross section of the output beam in a plane perpendicular to the direction of propagation is called the output aperture, and denoted byT . (Certain regularity conditions apply toD andT , see [10] for the technical details.) Denote points in space R 3 by pairs (x, z), where x ∈ R 2 is the position vector in a plane perpendicular to the direction of propagation and z ∈ R is the coordinate in the (negative) direction of propagation. See again Figure 1 for our convention on the direction of the z−axis. Points on the unit sphere S 2 will typically be denoted by m ∈ S 2 ; their components are also written as m = (m x , m z ) with |m x | 2 + m 2 z = 1. We fix the output aperture in the plane z = −d. We will seek to represent the first reflector as the graph of its polar radius ρ(m) (for m ∈D), and the second reflector as the graph of a function z(x) (for z ∈T ). (See Figure 1.) That is The geometrical optics approximation is assumed. It follows from general principles of geometric optics that all rays will have equal length from (0, 0, 0) to the plane z = −d; this length is called the optical path length and will be denoted by L. We define the reduced optical path length as = L − d.
Oliker [15] showed that local energy conservation translates into a complicated partial differential equation of Monge-Ampère type for ρ(m). As noted in [15], the resulting equation is quite involved and a rigorous analysis of this equation seems very difficult. (See equation (59) in [15].) To amend this, the problem was reformulated in [10] as a linear programming problem, which makes a complete analysis possible, both concerning theoretical results on existence and uniqueness, and gives a method for practical computations.
For this, the following function K(m, x), called the cost function in analogy with the theory of optimal transportation, plays an important role: In further preparation, the following two transformations are needed: Definition 2.1. Let z = z(x) be a continuous function defined onT ⊆ R 2 . Then define the functionz Definition 2.2. Let ρ = ρ(m) be a continuous function defined onD ⊆ S 2 with ρ > 0. Then define the functionρ With these preparations, in [10], the following notion of a quasi-reflector pair and its associated ray tracing map is used. (The term used in [10] is "reflector pair", but we use "quasi-reflector pair" here for clarity.) Definition 2.4. Let (ρ, z) ∈ C >0 (D) × C(T ) be a quasi-reflector pair. Define its reflector map, or ray tracing map, as a set-valued map γ :D → {subsets ofT } via In [10], it is shown that γ(m) is in fact single-valued for almost all m ∈D. Therefore γ may be regarded as a transformation γ :D →T .
If (ρ, z) is a "quasi-reflector pair" in the above sense,the following can be shown [10], justifying the choice of nomenclature: If physical copies of the corresponding surfaces (1) are made from a reflective material, then a ray emitted from the origin in the direction m ∈D will be reflected to a ray in the negative z−direction labeled by x = γ(m) ∈T .
With this, the reflector problem can be formulated rigorously: Problem I. (Reflector Problem) For given input and output intensities I(m), m ∈D, and L(x), x ∈T , respectively, satisfying global energy conservation D I(m) dσ = T L(x)dx, find a pair (ρ, z) ∈ C >0 (D) × C(T ) that satisfies the following conditions: (i) (ρ, z) is a quasi-reflector pair in the sense of Definition 2.3.
(ii) The ray tracing map γ :D →T satisfies for any Borel set τ ⊆T .
Here dσ denotes the standard area element on the sphere S 2 . The second condition is local energy conservation.
As we indicate below, the Reflector Problem can be reformulated in the following form: Problem II. Minimize the functional Problem II is the dual problem to the problem of finding a map P fromD ⊆ S 2 tō T ⊆ R 2 which minimizes the transportation functional P → D K(m, P (m))dσ. Problem I and Problem II are equivalent, as expressed in the following theorem.
be a quasi-reflector pair. Then (logρ, logz) ∈ Adm(D,T). The following statements are equivalent: Thus solving the reflector construction problem in Problem I is equivalent to solving the infinite dimensional LP problem in Problem II. Indeed, it can be shown that a solution as in Theorem 2.5 exists. (See Corollary 6.5 in [10].) In the remainder of this paper, we will concentrate on numerical solutions for this problem.

Numerical Schemes
In the following, we describe two scheme for solving Problem II numerically. First, let us list explicitly the given data on which the solution depends: The first numerical scheme, Scheme 3.1 is a very straightforward discretization of Problem II via meshes on the domainsD andT . Because of its simplicity and straightforwardness, we can informally call this a "brute force" technique. The problem with this scheme is that each pair of points fromD andT gives rise to a constraint, and so the technique requires a lot of memory, so that only coarse meshes can be used. (This problem is very well known, see for example [18].) To address this issue, we propose a more sophisticated scheme, Scheme 3.2, which is an iterative scheme using finer and finer meshes, where in each iteration, information from the previous solution is utilized to reduce the number of constraints.
The technique depends on choosing a parameter ε, the constraint threshold, in each iteration. We then discuss the problem of choosing an appropriate sequence of such ε.

"Simple Discretization" Scheme
• Similarly, create a mesh in the output domainT by choosing sets t where the interiors of any t Here we may assume that m (0) 1 = m 1 is the same point as given in the data in (iv) above.
• Find the solution {r of the following LP problem: Minimize 1 = logρ 1 (whereρ 1 is as given in the data in (iv) above), and r for j = 1, . . . , N (0) . (We are using the superscript "(0)" here to distinguish this solution from additional iterative approximations we will obtain with the iterative scheme described below.) It is important to point out that this discretization also yields a discretized version of the ray tracing map γ (see Definition 2.4). Namely, for each index i = 1, . . . , M (0) , there is at least one corresponding index j * , 1 ≤ j * ≤ N (0) where the constraint corresponding to the pair (i, j * ) is active, that is, where we have This means that the point m

"Iterative Refinement" Scheme
As noted in the introduction, one of the drawbacks of scheme 3.1 is that the constraint set in the linear programming problem in the penultimate step becomes large very fast with finer meshes, and the corresponding problem becomes too large to handle for standard LP solvers. We were not able to solve problems with more than very roughly 1000 mesh points on a standard PC with 4 GB RAM (M (0) ≈ 500, N (0) ≈ 500), which is arguably too small for many applications.
We addressed this problem by developing an iterative scheme. First the problem is solved for a mesh with relativey few sample points (say M (0) = N (0) = 250). Then a finer mesh is chosen, and the previous solution is used to reduce the number of constraints.
Specifically, we have the following scheme, depending on a number ε (1) > 0, which we call the "inclusion threshold," to be explained in detail below.
• Create a mesh onD by choosing sets d (1) is chosen larger than M (0) , meaning that we have a finer mesh.
• Similarly, create a mesh onT by choosing sets t (1) (1) ). Here we may assume that m (1) 1 = m 1 is the same point as given in the data in (iv) above.
• Interpolate the initial solution {r i=1 to find an approximation for the first reflector. That is, by some interpolation method, find a continuous function r(m), m ∈D, with j=1 to find an approximation for the second reflector. That is, by some interpolation method, find a continuous function

• Find the solution {r
(1) of the following LP problem: • Find the numbers ρ The idea of the second scheme 3.2 is to solve the discretized LP on a coarse mesh first and then use this information to reduce the number of constraints needed for the LP on a finer mesh. To wit, the difference between schemes 3.1 and 3.2 is that in the discretized LP problem, in scheme 3.2 not all pairs of sample points fromD andT , represented by pairs of indices (i, j), are included in the list of constraints. In fact, we would in principle only need to include those where the constraint is active, that is, where the corresponding inequality holds with equality. Of course, this information is not available a priori. Instead, we use the following heuristic: A constraint should only be included in the LP problem if it is "almost" active when we use an interpolation of the solution on the coarse mesh as an approximate solution on the fine mesh. This is where the "inclusion threshold" ε (1) come into play. Namely, we only include those constraints (i, j) for which the difference r(m is less than the inclusion threshold ε (1) . (Here r(m) and ζ(x) are interpolations of the solution of the LP problem for the coarse mesh.) Clearly, increasing ε (1) means more constraints are included, but also potentially the solution of the LP represents a better approximation of the exact reflector pair.
In our numerical tests, we found it advantageous to scale the functions I(m) and L(x) so that the approximations for the integrals D Idσ and T bLdx yield exactly the same result.
Note that the algorithm in scheme 3.2 can be iterated again. We can use the first iterative solution ρ  This is a question of great practical importance. To wit, if the sequence {ε (k) } k is constant or decreases too slowly, then "many" constraints are included in each iteration, meaning that the LP problems will become large fast and we will run out of memory quickly after a few iterations. If the sequence {ε (k) } k decreases too fast, then "few" constraints are included in each iteration. This could cause for instance that some index i (with 1 ≤ i ≤ M (k) ) may not even be included in any of the constraints with any of the indices j (with 1 ≤ j ≤ N (k) ), or vice versa. This would cause the LP problem to be unbounded and hence there would be no solution. (This behavior is illustrated in the practical tests summarized in Table 1; see also the discussion in section 4.) We can make a very rough estimate for a good choice of the sequence {ε (k) } k , assuming for simplicity that M (k) ≈ N (k) . Let us assume that each sample point m . Then the set of all points in (m, x)−space that satisfy is approximately an ellipsoid for small ε (k) . This can be seen by using the Taylor expansion around the pair (m  i , x (k) j ) that satisfies this inequality corresponds to a constraint included in the LP of the k th iteration, the total number of constraints should be roughly proportional to M (k) · (ε (k) ) 2 . So to keep the number of constraints approximately constant throughout the iterations, based on these heuristics, one may choose where C is a constant. These heuristics are very rough and in practice the number of constraints is still increasing with M (k) , albeit at a much slower rate than the case of constant ε (k) . In section 4.3, we describe an implementation of this idea.

Numerical Tests
To test the validity of the numerical scheme described above, we used it on a case where the solution is known in analytic form. In the next subsection, we first describe this special analytic solution, and then discuss our results.

An Analytic Solution
In order to obtain a data set where an explicit analytic solution is known, we consider a special given configuration of reflectors, and then solve the "forward" problem of determining the output intensity this system produces for a given input intensity.

Construction of a pair of reflectors R 1 , R 2
Consider the following set of two reflectors R 1 and R 2 , sketched in Figure 2: Let a = (a, b, c) be a given point in R 3 , and let R > 0 and α > 0 be two positive numbers with R > |a| = √ a 2 + b 2 + c 2 . Now let the first reflector R 1 be the boundary of a spheroid with foci at the origin O and at a and with major diameter R/2. (Thus for each point P on R 1 , the sum of the distances from each of the two foci PO + Pa equals R.) Let the second reflector R 2 be the boundary of a paraboloid whose main axis is the negative z−axis and whose focus is at a with focal parameter 2α. (Thus for each point p = (x, y, z) on R 2 , the sum of the distance to the focus and the shifted z−component Pa + (z − c) equals 2α.) Note that by definition of the two reflectors, any cone of light rays emitted from the origin O will be transformed into a beam of parallel rays traveling in the direction of the negative z−axis. Indeed, a ray emitted from O will be reflected off R 1 towards the focus a. This ray will then be reflected in the direction of the negative z−axis by R 2 . See again Figure 2.
It is not hard to find explicit expressions for the two reflectors. If we write R 1 = {ρ(m) · m m ∈ S 2 } and R 2 = {(x, z(x)) x = (x, y) ∈ R 2 }, then we have the following expressions: and where we used the shifted coordinates (p, q) Figure 2: Sketch of the reflectors R 1 and R 2 . R 1 is the boundary of an ellipsoid whose foci are the points O and a. R 2 is the boundary of a paraboloid whose focus is a. The sketch shows a two-dimensional cross section. A ray given by the direction m ∈ S 2 will be reflected by R 1 and then by R 2 to a ray in the direction of the negative z−axis.

Ray Tracing Map γ
One can now determine the ray tracing map γ corresponding to the reflector pair R 1 , R 2 . See again Figure 2.
Let m ∈ S 2 be a given direction. Thus m = (m x , m z ) with m x ∈ R 2 and |m x | 2 + m 2 z = 1. To determine the image γ(m), consider the ray emitted from the origin in the direction m. The ray will encounter the first reflector at the point ρ(m) · m and will then be reflected towards the second focus a. Thus the reflected ray can be parameterized by a + λ · (a − ρ(m) · m), where λ is a parameter. With equation (7), this yields that the point where the reflected ray hits the second reflector is a + λ * · (a − ρ · m) where .
The projection of this point to the xy−plane is the value of the ray tracing map. Thus we have the following explicit formula for the ray tracing map:

Solution of the Forward Problem
We can now solve the forward problem: Given the reflector pair R 1 , R 2 as defined above, an input apertureΩ ⊆ S 2 and an intensity distribution I(m), m ∈Ω, find the output aperturē T ⊆ R 2 and the output intensity L(x), x ∈T . The output aperture is simply the image ofΩ under the ray tracing map γ:T = γ(Ω). To find the induces intensity L(x), use the defining property ω I(m)dσ = γ(ω) L(x)dx (9) for all Borel sets ω ⊆Ω. This is an energy balance equation. The above integral equation allows us to find an explicit expression for L(x) given I(m), or vice versa. We consider for simplicity the case thatΩ is contained in the left hemisphere S 2 − = S 2 ∩ {(x, y, z) ∈ R 3 z < 0}. We can then use coordinates In these coordinates, the standard measure on S 2 is given by dσ = dm x dm y |m z | , where Using this, (9) is then equivalent to the equation Here denotes the Jacobian of the map γ • τ , With the help of a computer algebra system like Mathematica, this can be evaluated explicitly as where again m z = − 1 − m 2 x − m 2 y .

A data set
We can now construct a particular data set using the results of the previous sections. We pick a = b = 0, and c = −0.4, α = 1, R = 1.3, and consider the input aperturē which is a spherical cap centered at the point (0, 0, −1) 3 . Using the results from the previous subsections, we can now find the output aperture in a straightforward manner: We choose a constant output intensity Using the relation (12), this gives the input intensity The two reflectors are given as follows: radial radius ρ for R 1 : ρ(m) = 0.765 and equation for R 2 : z = −0.25 x 2 + y 2 + 0.6.
The corresponding reduced optical path length is This follows from the construction of the two reflectors and the properties of ellipsoids and paraboloids. The two reflectors are shown in Figure 3.

Results
Using the data set from subsection 4.2, we conducted a series of numerical tests to establish the validity of the numerical scheme 3.2. Since we have explicit analytical expressions for the two solutions ρ(m) and z(x) describing the two reflectors as given in (14) and (15), respectively, we were able to compute the error of approximation directly. (See also Figure 3 for plots of these reflectors.) The iterative scheme 3.2 was implemented in MATLAB with the package distmesh [17] for mesh generation in each iteration step and lp solve [4] for solving the resulting sparse LP. Computations were performed on an Intel I5 Dual Core/2.53GHz with 4GB RAM. Our primary aim is to demonstrate the practicability of the proposed scheme. While there are several ways in which the implementation could be made more efficient, for instance by using a compiled computer language or a more specialized LP solver for transportation problems, we believe that our results show that the proposed scheme is easy to implement and makes it possible to use much finer meshes than a simple discretization scheme.
The mesh generation algorithm employed by the software we used is based on a relaxation scheme of forces in a truss structure [17]. For this, user input is the desired edge length h, not the number of mesh points M (for the domainΩ) or N (for the domainT ), respectively. However, it is not hard to see that the relation between the average edge length h and the number of mesh points M onD obeys Indeed, the total area A is proportional to the number of mesh points M times the area of one triangle of a Delaunay triangulation. The area of a triangle is proportional to the square of the average edge length h. The same relation holds of course for the number of mesh points N onT . (In Tables 1 and 2 in this section, we show the number of mesh points M, N instead of the edge length h, as we believe the former are more informative.) We chose a sequence of desired edge lengths h (0) , h (1) , . . .. Then, based on the heuristics given above and in section 3.3, we determined the corresponding constraint thresholds with the formula where C and a are constants. We tested several values for the proportionality constant C and the exponent a. (The heuristic arguments in section 3.3 yields a = 1, but we also tested other values.) Some results are summarized in Table 1. As indicated in section 3.3, "large" values for C mean that "many" constraints are included in the LP in each step, which potentially means good approximations, but also high memory usage, and thus we may run out of memory after a few iterations. In contrast, "small" values for C improve memory usage, but may mean that the numerical results may be less accurate approximations of the exact results, or indeed the problem may become unbounded, as typically seems to happen in practice. (See Table 1.) A similar logic holds for the exponent a: For larger a, the threshold ε (k) decreases faster, leading to less memory usage, but the danger of the problem becoming unbounded. Table 1 confirms and illustrates these principles. At each iteration, the LP is either unbounded or there is a solution. Interestingly, if there is a solution, it seems largely independent of the number of constraints, as indicated by the fact that the error of approximation varies very little. (See fourth column in Table 1.) This makes some intuitive sense as the most important issue in each iteration is really that we have included the active constraints. As long as all active constraints are included, we get the same solution. If however we did not include any of the constraints for a certain mesh point, the problem becomes unbounded.
In fact, Table 1 illustrates the following feature of the algorithm, which is intuitively quite obvious, although we do not give a formal proof: At each iteration step, there is a critical value ε crit such that if ε < ε crit , then the LP is unbounded because no constraints involving a certain point are included. If ε > ε crit , the LP has a solution. As indicated above, it also appears that if there is a solution, it is likely typically close to the solution of the "full" problem (that is, the problem including all possible constraints). This leads to a possible alternative algorithm for choosing the thresholds ε (1) , ε (2) , . . . : Instead of choosing each ε (k) with a pre-determined formula, use the corresponding critical value ε crit in each iteration step. This critical value can be found by simple trial and error. While this idea is promising, it is possible that such a scheme may lead to a much slower decay of the maximum error, or even an increase of the error. We did not test this idea in this article, but it would be interesting to analyze it further.
To find out the maximum mesh size on which a solution could be obtained without running out of memory or arriving at an unbounded problem, we "calibrated" the values of C and a. The process is summarized in Table 2.
We obtained best results (that is, largest mesh sizes) for the run with C = 1.7 and a = 1. Table 3 now summarizes the results for this specific run in more detail. The table shows the errors between the numerical computed approximations and the exact formulas as given in (14) and (15). Figure 4 shows plots of the approximation errors for the reflectors as functions on the input and the output apertures.
These results are a good indication that the scheme converges, and that the error decreases proportional to M −α , or equivalently proportional to h 2α , where roughly 0.5 α 1. (The results from Figure 1 appear to indicate α 1, but in other data, we also encountered results where α seemed to be closer to 0.5.) The results also illustrate the practicability of the scheme. Note that the number of mesh points was increased by a factor of about 16 from iteration 1 to iteration 6. In fact the problem in iteration 6 would be impossible to solve with a simple discretization scheme due to the size of the problem. In iteration 6, we were able to reduce the number of constraints used to only about 15% of the number of constraints used for the "naive" simple discretization scheme 3.1.   (17). The third column shows the sequence of mesh points M (k) forD. (So for instance, in the first row, in the initial discretization, there were 284 points. In the first iteration, there were 455 mesh points. The second iteration had 724 points, but the problem was unbounded.) The third column notes if the problem became unbounded before reaching 1147 mesh points. If the total run was completed up to 1147 mesh points, the third column gives the maximum norm ρ num − ρ ana ∞ of the error between the resulting numerical solution ρ num and the analytic solution ρ ana . The fifth column shows the number of constraints for the case of 1147 mesh points, if the run reached this maximum number. The last column indicates the corresponding percentage of number of constraints relative to the total number of constraints in the "brute force" scheme 3.   (17) for maximum mesh sizes. Each row represents a different run of the iterative numerical scheme. The third column is the number k of iterations in the run and the fourth is the number M (k) of points in the mesh forD in the last iteration. The fifth column notes the result of the run -if a solution was found, the maximum error is given. If not, the linear program either became unbounded, or the computer ran out of memory. The last two columns are the number of constraints used for the last iteration and the percentage of total possible constraints (compared to the "brute force" scheme 3.1).    Table 3 for reflector 1 (left) and reflector 2 (right). For better visibility, the plots were obtained by interpolating the error at each mesh point to a continuous function. Note that the error distribution appears to be relatively uniform, although the "spikes" are noticably arranged in triangular patterns and in a fashion somewhat reminiscent of spokes on a wheel. This may be an artifact of an iterated error in conjuncture with triangular meshes.

Conclusion and Future Work
We proposed two numerical schemes for solving an infinite dimensional optimal transportation problem arising in reflector design: A straightforward discretization and an improved iterative scheme, which uses knowledge of the previous solution in each step to reduce the number of constraints. This scheme is easily adapted to similar transportation problems arising in beam shaping problems, e.g. in [9] and [11]. We showed that this new scheme is easy to implement and makes it possible to solve the problem on much finer meshes. There are a number of possible directions for future research. We did not rigorously prove that the scheme converges, although we strongly expect that it does. In fact the decreasing error shown in Table 3 gives evidence for this, at least in our example. It would be valuable to have a rigorous proof of convergence.
Another direction for further investigation is to change the selection algorithm for the threshold values ε (1) , ε (2) , . . . as suggested in section 4.3: It is a straightforward conjecture that in each step, there is a critical threshold inclusion number such that the problem is unbounded for values of ε below that number. One way to select a value of ε is to pick a value for ε (k) close to the critical value ε crit . It would be interesting to flesh out this ideas and analyze the results more closely.
Another avenue of research is to see whether our proposed solution algorithm can be made parallizable.
It would also be worthwhile to investigate whether the existing PDE-based schemes for solving the transportation cost with quadratic costs [2,3,7,23] can be adapted to the problem at hand.