ISRN Applied Mathematics

Volume 2013 (2013), Article ID 635263, 12 pages

http://dx.doi.org/10.1155/2013/635263

Research Article

Iterative Scheme for Solving Optimal Transportation Problems Arising in Reflector Design

Tilmann Glimm¹ and Nick Henscheid²

¹Department of Mathematics, Western Washington University, Bellingham, WA 98225, USA

²Program in Applied Mathematics, University of Arizona, Tucson, AZ 85716, USA

Received 29 July 2013; Accepted 5 September 2013

Academic Editors: M.-H. Hsu, G. Mishuris, L. Rebollo-Neira, and Q. Song

Copyright © 2013 Tilmann Glimm and Nick Henscheid. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

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 previously that this problem is equivalent to an infinite-dimensional linear programming (LP) problem. Here we investigate techniques for constructing the two reflectors numerically by considering the finite-dimensional LP problems which arise as approximations to the infinite-dimensional problem. A straightforward discretization has the disadvantage that the number of constraints increases rapidly with the mesh size, so 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, where information from the previous iteration step is used to reduce the number of necessary constraints. As an illustration, we apply our proposed scheme to solve a problem with synthetic data, demonstrating that the method allows for much finer meshes than a simple discretization. We also give evidence that the scheme converges. There exists a growing literature for the application of optimal transportation theory to other beam shaping problems, and our proposed scheme is easy to adapt for these problems as well.

1. Introduction

The following beam shaping problem from geometric optics was described in [1]; see Figure 1. Suppose a point source emits a spherical wavefront with a given intensity distribution. The problem   we are concerned with 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 [1].

Figure 1: Sketch of the reflector problem. The point source is located at the origin of the coordinate system. The coordinate system in the lower left-hand corner explains our use of coordinates. The output beam propagates in the direction of the negative -axis. Points in the plane perpendicular to the -axis are denoted by the vector . Correspondingly, points in three-dimensional space are denoted by . See Section 2 for more. Reprinted from [2].

The first rigorous mathematical solution to the problem was provided in [2], using an approach based on the theory of optimal transportation [3–6]. See also the references [7–11] which deal with other beam shaping problems using related techniques.

The result of [2] is summarized in Section 2 below, with Theorem A.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 linear 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 Theorem 3 in Section 2 for the exact statement).

This reformulation of the problem not only is 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. This was already observed for similar beam shaping problems by Glimm and Oliker [7, 8] and independently by Wang [9].

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 with quadratic cost (also called optimal transport) [12–20], although some authors have treated costs proportional to the distance ( optimal transport) [21]. These methods are based on fluid mechanical approaches [13] or various finite difference approaches [17]. 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.

A linear programming approach similar to the one we investigate here for more general situations for optimal approach was investigated by Rüschendorf and Uckelmann in 2000 [22] and more recently in two papers by Canavesi et al. [23, 24] who also propose a new algorithm inspired by, but distinct from, a linear programming approach for the single reflector problem treated in [8, 9].

As noted in [22] and other works, an immediate obstacle for the linear programming approach is that the number of constraints becomes very large even for relatively coarse meshes—if there are sample points in the input aperture (denoted by in Figure 1) and points in the output aperture ( in Figure 1), then the number of constraints is . For instance, Rüschendorf and Uckelmann note in [22] that the LP solver they used, SOPLEX, was designed to handle up to 2 million variables. This corresponds to mesh sizes of approximately 1400 points on the input and output aperture in our problem if one employs a straightforward discretization scheme. We found that, in our implementation, the solver MOSEK (Version 7) could handle about 3 million constraints, corresponding to approximately 1750 sample points on the domains and (see Table 1 for details). Rüschendorf and Uckelmann noted in [22] that, for better results, one needs carefully designed programs. This is what we supply here for our problem: 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 drastically slows the growth of the number of constraints. Details are described in Section 3. With this step-wise mesh refinement scheme, as an illustration, we were able to obtain solutions on meshes with about 10,300 points on each aperture using MATLAB and MOSEK 7, a commercially available LP solver. The number of constraints needed for this was only a fraction of the size of the “full” system. For instance, the run displayed in Table 2 used only 0.33% of all possible constraints.

Table 1: Results of runs with the “Simple Discretization” Scheme 5 for several edge lengths . Memory was insufficient for the LP solver to handle the case , so the maximum number of mesh points that could be handled was about 1,750 on each aperture. The data used was the data from Section 4.1. Symbols used: : edge length; : number of mesh points in ; : number of mesh points in ; time (s): time in seconds; constr.: number of constraints.

Max error
error
Max error
error
Time (s)
Constr.

0.09
511
522

2.97016
266,742

0.08
660
663

4.67627
437,580

0.07
880
874

10.1342
769,120

0.06
1,214
1,214

17.4949
1,473,796

0.05
1,752
1,757

45.9136
3,078,264

Table 2: Results of an iterative run with the “Iterative Refinement” Scheme 6. Memory was insufficient to handle the sixth iteration due to matrix size limitations in MATLAB. Note the run times (last two rows). The bulk of the algorithm run times is taken in setting up the linear programming problem, that is, deciding which constraints to include. Solving the LP once it is set up is fairly fast (last column). The data used was that of Section 4.1, the number of nearest neighbors searched over , and the quasi-activity threshold was . : iteration number; : number of mesh points in ; : number of mesh points in ; constr.: number of constraints; con. dens.: constraint density, that is, used constraints as a percentage of all possible constraints ; time: computation time for iterative step in seconds, including generating meshes, building, and solving the LP problem; LP time: time in seconds to solve the LP problem.

Max error
error
Max error
error
Constr.
Con. dens.
Time (s)
LP time (s)

0

1,059
1,059

1,121,481
100%
15.962
12.3146

1

1,675
1,678

56,080
2.00%
208.533
0.938877

2

2,647
2,647

89,724
1.28%
392.982
1.59831

3

4,171
4,175

141,550
0.81%
969.239
2.54405

4

6,582
6,567

224,327
0.52%
2363.79
4.32524

5

10,313
10,318

351,514
0.33%
5909.25
9.21891

We note that our proposed algorithm does not make any a priori symmetry assumptions on the form the reflectors. It can also be adapted for the numerical solution of other beam shaping problems for which a formulation using optimal transportation theory has been found, for example those in [7, 8, 25].

This article is organized as follows: in Section 2, we recall the reformulation of the reflector construction problem as a linear programming problem as given in [2] and fix some notation. We then describe a basic discretization scheme for the numerical solution and propose an improved “mesh refinement” scheme in Section 3. The following Section 4 is devoted to numerical illustrations and 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. We conclude with a brief summary of the advantages and drawbacks of the method and propose future work.

2. Formulation as Linear Programming Problem

We now briefly review the notation for the problem posed in [1], as well as the result of [2], which reformulates the reflector design problem as a linear programming problem. More details are contained in Appendix A.

Consider the configuration shown in Figure 1. A point source at the origin generates a spherical wavefront over a given input aperture contained in the unit sphere . We use the bar in the notation for to stress that is a closed set. The input beam has a given intensity distribution. By means of two reflectors, this wavefront is to be transformed into a beam of parallel rays propagating in the direction of the negative -axis. This output beam is required 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 is denoted by . Again, we use the notation to stress that is closed. Certain regularity conditions apply to and ; see [2] for the technical details.

Denote points in space by pairs , where is the position vector in a plane perpendicular to the direction of propagation and is the coordinate in the (negative) direction of propagation. See again Figure 1 for our convention on the direction of the -axis. Points on the unit sphere will typically be denoted by ; their components are also written as with .

We fix the output aperture in the plane . We will seek to represent the first reflector as the graph of its polar radius (for and the second reflector as the graph of a function (for ). 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 to the plane ; this length is called the optical path length and will be denoted by . We define the reduced optical path length as .

Oliker [1] showed that local energy conservation translates into a complicated partial differential equation of Monge-Ampère type for . As noted in [1], the resulting equation is quite involved, and a rigorous analysis of this equation seems very difficult. See equation (59) in [1].

To amend this, the problem was reformulated in [2] as an infinite dimensional 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 , 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 1. Let be a continuous function defined on . Then define the function

Definition 2. Let be a continuous function defined on with . Then define the function

The following characterization of the solution of the reflector problem described above was given in [2].

Theorem 3. If the pair is a solution of the above reflector design problem, then the transformed pair minimizes the functional

on the function space   for all.

One can solve the reflector construction problem by solving the infinite-dimensional LP problem in the above theorem and then applying the inverse transformations of (3) and (4) to recover and . In the remainder of this paper, we will concentrate on numerical solutions for this problem.

It is also possible to characterize the ray-tracing map analytically.

Definition 4. Let be a solution of the reflector problem. Define its ray-tracing map as a set-valued map via

In [2], it is shown that is in fact single-valued for almost all . Therefore, may be regarded as a transformation . Indeed, one can show that a ray emitted from the origin in the direction will be reflected to a ray in the negative -direction labeled by ; see [2].

More details and rigorous statements are found in Appendix A.

3. Numerical Schemes

In the following, we describe two schemes for solving the minimization problem in Theorem 3 numerically. First, we list the given data on which the solution depends.

(i)Domains   , (The rigorous result in [2]   required   the additional technical assumption . For practical purposes, these assumptions can often be dropped.)(ii)Nonnegative integrable functions , , and , , with .(iii)A reduced optical length . (iv)The value for some fixed , where is the radial function of the first reflector. (Note that without this constraint, the reflectors are not uniquely determined. This can be seen in the expression for the functional in Theorem 3. Indeed, any transformation of the form , for constant will leave the objective functional unchanged.)

The first numerical scheme, Scheme 5 is a straightforward discretization of the minimization problem in Theorem 3 via meshes on the domains and .

As described in the introduction, this scheme is memory-intensive because every pair of points from and gives rise to a constraint. To address this issue, we propose a more sophisticated scheme, Scheme 6, 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.

3.1. “Simple Discretization” Scheme

The characterization of solutions to the reflector problem by means of the minimization of a linear functional in Theorem 3 immediately suggests a solution algorithm by means of a discretization. Essentially, the same algorithm was used for different problems in [22–24].

Scheme 5 (simple discretization). (i) Create a mesh in the input domain by choosing sets , where the interiors of any and are disjoint (for ) and . For instance, may be a triangulation of .

(ii) Similarly, create a mesh in the output domain by choosing sets , where the interiors of any and are disjoint (for ) and

(iii) Choose sample points (for ) and (for ). Here, we may assume that is the same point as given in the data in (iv) above.

 (iv) Find the solution of the following LP problem:

(v) Find the numbers and such that and . This is straightforward by taking the inverse of the transformations given in Definition 1 and Scheme 5.

Then is an approximation for the true value of the radial function of the first reflector, evaluated at the sample point for . Similarly, is an approximation for the true value of the function describing the second reflector, evaluated at the sample point for . We are using the superscript “(0)” here to distinguish this solution from additional iterative approximations we will obtain with the iterative scheme described below.

Note that (7) is a discretization of the sum of integrals in the functional . We have chosen the simplest discretization here. Depending on the geometry of the problem, other discretization schemes [26] may yield better results. In practice, we have chosen partitions and of approximately uniform size and diameter, but depending on the form of the intensity distributions, other partitions may be more appropriate.

It is important to point out that this discretization also yields a discretized version of the ray-tracing map (see Definition 4). Namely, for each index , there is at least one corresponding index , where the constraint corresponding to the pair is active; that is, where we have

This means that the point is approximately the image of the point under the ray-tracing map .

3.2. “Iterative Refinement” Scheme

As noted in the introduction, one of the drawbacks of Scheme 5 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 about 1750 mesh points on each aperture on a standard PC with 4 GB RAM ().

We addressed this problem by developing an iterative scheme. First, the problem is solved for a mesh with relatively few sample points (say ). 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 , which we call the “quasi-activity threshold,” and an integer , the nearest neighbor search depth, to be explained in detail below.

Scheme 6 (iterative refinement). Use Scheme 5 to find an initial solution for sample points on and sample points on , respectively.

Definition 7. Call a pair with and “quasi-active” if

Note that the pair is “active” in the standard sense of linear programming theory if the above term is equal to zero.

(i)Create a second, finer mesh on by choosing sets , where the interiors of any and are disjoint (for ) and Here, is chosen larger than , meaning that we have a finer mesh. (ii) Similarly, create a second mesh on by choosing sets , where the interiors of any and are disjoint (for ) and (iii) Choose sample points (for ) and (for ). Here, we may assume that is the same point as given in the data in (iv) above.(iv) For each pair with , and , find the nearest neighbors of from the previous mesh on and the nearest neighbors of from the previous mesh on .

Definition 8. Say that the pair is potentially active if at least one of the pairs consisting of a nearest neighbor of and a nearest neighbor of is quasi-active in the above sense.

(i) Find the solution of the following LP problem:

(ii) Find the numbers and such that and

The idea of Scheme 6 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. A key difference between Schemes 5 and 6 is that in the discretized LP problem, in Scheme 6, not all pairs of sample points from and , represented by pairs of indices , 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 the corresponding point on the second reflector is close to the image of the corresponding point on the first reflector under the ray-tracing map. In terms of the LP problem formulation, the corresponding pair of points is located close to an active pair from the previous iteration. We call such a constraint “potentially active.” Our heuristic is that a constraint should be considered potentially active if there is a pair of its respective nearest neighbors which corresponds to a quasi-active constraint form the previous iteration. Here, “quasi-active” constraints include active constraints (where the constraint is attained, i.e., the term is zero) and those which are “almost active” in the sense that the term is positive, but small. “Smallness” is encoded in the quasi-activity threshold . Clearly, increasing 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 and so that the approximations for the integrals and yield exactly the same result.

Note that the algorithm in Scheme 6 can be iterated again. We can use the first iterative solution and , in Scheme 6, to obtain a second iterative solution and , and so on.

Note that (7) is a discretization of the sum of integrals in the functional . We note that we have chosen again the simplest discretization here for the sum of integrals . As in Scheme 6, depending on the geometry of the problem, other discretization schemes [26] may yield better results.

3.3. Choice of Thresholds for Iterative Refinement Scheme 6

For the iteration of Scheme 6 described above, we have to choose a corresponding sequence of mesh sizes . We also need to choose a corresponding sequence of thresholds .

This last question is of a certain practical importance. Indeed, if the thresholds are “high” (e.g. if they decreases slowly), then “many” constraints are included in each iteration, meaning that the size of the LP problems increases quickly. If the sequence is “low”, then “few” constraints are included in each iteration. This is of course in principle desirable, but if the thresholds are too low, this could cause that some index (with ) may not even be included in any of the constraints with any of the indices (with or vice versa. This would cause the LP problem to be unbounded, and hence, there would be no solution.

We can make a very rough estimate for a good choice of the sequence assuming for simplicity that . Let us assume that each sample point in the input aperture is paired via the ray-tracing map with a unique point and vice versa. Thus, there are such pairs. Let and denote the solutions of the reflector problem. Then the set of all points in -space close to a given pair that satisfy is approximately an ellipsoid for small . This can be seen by using the Taylor expansion around the pair in -space. The volume of this ellipsoid is proportional to . So the proportion of pairs of points that satisfy the above constraint among all constraints is roughly proportional to . Since each pair that satisfies this inequality corresponds to a constraint included in the LP of the iteration, the total number of constraints should be roughly proportional to . So to keep the number of constraints growing linearly throughout the iterations, based on these heuristics, one may choose

where is a constant. These heuristics are very rough. In practice, we have found that when using the formula , there is a range of values of which all produce similar results. See Section 4.2.

4. 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.

4.1. A Data Set

We have constructed a particular data set using an explicit analytical solution described in more detail in Appendix B. The data set of this section is obtained from the more general solution of Appendix B by choosing and , , and . We consider the input aperture

which is a spherical cap centered at the point . The data do not represent a physically possible set of reflectors since there would be blockage. However, the example is useful as an illustration, and the numerical scheme itself is completely independent of the shape of the apertures or any a priori symmetries of the problem. Using the results from Appendix B, we can now find the output aperture in a straightforward manner as follows:

We choose a constant output intensity

Using the relation , this gives the input intensity

The two reflectors are given as follows:

The corresponding reduced optical path length is

The two reflectors are shown in Figure 2. is part of a spheroid, and is part of a paraboloid; see Appendix B for details.

(a)

(b)

Figure 2: Plots of the reflectors given by (18) (a) and (19) (b) in the explicit data set given in Section 4.1.

4.2. Results

Using the data set from Section 4.1, we conducted a series of numerical tests to establish the validity of the numerical Scheme 6. Since   there are   explicit analytical expressions for the two solutions and describing the two reflectors as given in (18) and (19), respectively, we were able to compute the error of approximation directly. See also Figure 2 for plots of these reflectors.

The Schemes 5 and 6 were implemented in MATLAB with the package distmesh [27] for mesh generation in each iteration step and the commercial linear programming solver MOSEK (version 7) [28] for solving the resulting LP problems. Computations were performed on an Intel Core i7-2620M/2.70 GHz with 4 GB RAM, running Windows 7 (64 bit).

The mesh generation algorithm is based on a relaxation scheme of forces in a truss structure [27]. For this, the input is the desired edge length , not the number of mesh points (for the domain ) or (for the domain ), respectively. However, it can be seen from geometric considerations that the relation between the average edge length and the number of mesh points on obeys

In the tables below, we show both the number of mesh points and the edge length , as we believe both are informative.

Our primary aim is to demonstrate the practicability of the proposed Scheme 6. For comparison, we first tested the performance of the “Simple Discretization” Scheme 5. The results are summarized in Table 1. The maximum mesh sizes attainable are about 1750 points, corresponding to approximately 3 million constraints. We found that this was the maximum size for the LP solver.

To test Scheme 6, we chose a sequence of desired edge lengths by reducing by in each step; that is, . Then, based on the heuristics given above and in Section 3.3, we determined the corresponding constraint thresholds with the following   formula:

where and are constants. The heuristic arguments in Section 3.3 yield , but we also tested other values. As indicated in Section 3.3, “large” values for mean that “many” constraints are included in the LP in each step, which potentially means good approximations, but also high memory usage. In contrast, “small” values for improve memory usage but may mean that the LP problem may become unbounded. This is the case if constraints that would be active in the LP problem comprising of all constraints are left out in an iterated LP problem. In our tests, the only time this happened was when we chose , that is, if the quasi-activity threshold is set to zero and only active constraints are considered “quasi-active.” We conclude that our results are very robust with respect to and , meaning for a broad range of values, we reached the terminal iteration with ; see Table 2. This is because the memory limiting issue in our tests was not the LP solver but the maximum size of arrays in MATLAB when computing the cost matrix; see paragraph below. We also conducted some informal tests of the algorithm dependence on the output and input intensities and , specifically to know whether the algorithm is stable under small perturbations of and . For this, we added small perturbation to and to , then renormalized so that the two integrals were the same. We found that the algorithm was indeed stable in the sense that the same maximum number of points was reached for different values of .

Results of a typical run are summarized in Table 2. 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 largest number of points were about 10,000 on each aperture. It is notable that what prevented our tests to include more plot points was not the LP solver but memory problems related to the maximum size of matrices in MATLAB during the set-up phase of the LP problem when computing the matrix of cost coefficients , which is of size . For the configuration we used, the maximum size of arrays in MATLAB is about 280 million variables, see [29].

These results give indication that the scheme converges and that the maximum error for the first reflector decreases proportional to , or equivalently proportional to , where , data from Table 2. The same formulation for scheme 5 yields an exponent of , unsurprisingly a slightly more rapid decrease, data from Table 1. The geometric distribution of errors is indicated in Figure 3.

(a)

(b)

Figure 3: Error surface plots for reflector 1 (a) and reflector 2 (b) for the data of Table 3 with . For better visibility, the plots were obtained by interpolating the error at each mesh point to a continuous function. Note that the error distribution is relatively uniform, although slight spikes of the error tend to occur near the boundary.

To investigate the role of the nearest neighbor search depth , we performed several tests with different values. The results are summarized in Table 3. The value of has very little influence on the error of approximation. This makes intuitive sense as the crucial issue in each iteration is that those constraints which are active in the “full” LP problem (i.e., the one that includes all possible constraints) are included. As long as all active constraints are included, one expects to get the same solution. Note, however, that the number of constraints changes quite drastically with , ranging from 96,256 ( of all possible constraints) for to 48,2058 ( of all possible constraints) for in this example. The time to solve the LP problems, thus, also varies from 2.82 seconds for to seconds for , although those gains are dwarfed by the overall run times of the algorithm which include assembling the LP problem.

Table 3: Results of the final iteration of the “Iterative Refinement” Scheme 6 with different values of the nearest neighbor search depth . All other data used was the same as in Table 2. Note that the nearest neighbor search depth has a large influence on the size of the LP problem but relatively little influence on the errors or run times. Column labels are identical to Table 2 except for the first column.

Max error
error
Max error
error
Constr.
Con. dens.
Time (s)
LP time (s)

2
0.02097152
10,313
10,318

96256
0.090%
5947.76
2.82425

4
0.02097152
10,313
10,318

226006
0.212%
5752.17
6.09849

6
0.02097152
10,313
10,318

351514
0.330%
5909.25
9.21891

8
0.02097152
10,313
10,318

482058
0.453%
6322.4
12.9115

5. Conclusion and Future Work

We investigated two numerical schemes for solving an infinite-dimensional optimal transportation problem arising in reflector design, a straightforward discretization and a novel iterative scheme, which uses knowledge of the previous solution in each step to reduce the number of constraints. The latter scheme is easily adapted to similar transportation problems arising in beam shaping problems, for example see [7–9].

As a proof of concept, we implemented both schemes using MATLAB and the commercial LP solver MOSEK and showed that the new scheme is easy to implement and makes it possible to solve the problem on much finer meshes. In particular, we demonstrated that the number of constraints are greatly reduced by the new scheme, to as little as 0.09% of the full number of constraints. Our aim was to show the validity of the scheme, and we contend that performance could be further improved, for instance, by using a compiled computer language or a more specialized LP solver for transportation problems.

One drawback of the scheme is that it is quite slow. See Table 2. The most time-consuming part of the algorithm was not the solution of the LP problem but assembling the problem, that is, deciding which constraints to be included. We believe that this step, in particular, could be improved by using a compiled computer language instead of MATLAB.

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 2 gives evidence for this. It would be valuable to have a rigorous proof of convergence. We also expect that PDE-based schemes for solving the transportation cost with quadratic costs [13, 17–19] could be adapted to the problem at hand, leading to faster algorithms.

Appendices

A. Details on Theorem 3

We give here some more details on Theorem 3 from Section 2. We use the same notation as in this section. This is a summary of definitions and results from [2].

In [2], the following notion of a reflector pair and its associated ray-tracing map is used.

Definition A.1. Let denote the set of all continuous functions on whose range is contained in the interval .   A pair is called a reflector pair if and

Here, notations (4) and (3) were used.

Definition A.2. Let be a reflector pair. Define its reflector map, or ray-tracing map, as a set-valued map via

In [2], it is shown that is, in fact, single-valued for almost all . Therefore, may be regarded as a transformation .

If is a “reflector pair” in the above sense, the following can be shown [2], 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 will be reflected to a ray in the negative direction labeled by .

With this, the reflector problem can be formulated rigorously.

Problem A.3 (reflector problem). For given input and output intensities , and , respectively, satisfying global energy conservation , find a pair that satisfies the following conditions:

(i) is a reflector pair in the sense of Definition A.1; (ii)the ray-tracing map satisfies

for any Borel set .

Here, denotes the standard area element on the sphere . The second condition is local energy conservation.

As we indicate below, the reflector problem can be reformulated in the following form.

Problem A.4 (minimize the functional).

on the set.

Problem A.4 is the dual problem to the problem of finding a measure-preserving map from to which minimizes the transportation functional . Problems A.3 and A.4 are equivalent, as expressed in the following theorem.

Theorem A.5 (see [2]). Let be a reflector pair. Then . The following statements are equivalent:

(i) solves the reflector Problem A.3. (ii) solves Problem A.4.

Thus, solving the reflector construction problem in Problem A.3 is equivalent to solving the infinite-dimensional LP problem in Problem A.4. Indeed, it can be shown that a solution as in Theorem A.5 exists. See Corollary 6.5 in [2].

B. An Explicit Analytic Solution to the Reflector Problem

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. This general family of solutions we obtained was used to generate the explicit synthetic data set in Section 4.1.

B.1. Construction of a Pair of Reflectors ,

Consider the following set of two reflectors and , sketched in Figure 4. Let be a given point in , and let and be two positive numbers with . Now, let the first reflector be the boundary of a prolate spheroid (ellipsoid of revolution) with foci at the origin and at and with major diameter . Thus, for each point on , the sum of the distances from each of the two foci equals . Let the second reflector be the boundary of a paraboloid whose main axis is the negative axis and whose focus is at with focal parameter . Thus, for each point on , the sum of the distance to the focus and the shifted -component equals .

Figure 4: Sketch of the reflectors and . is the boundary of an ellipsoid whose foci are the points and . is the boundary of a paraboloid whose focus is . The sketch shows a two-dimensional cross section. A ray given by the direction will be reflected by and then by to a ray in the direction of the negative axis.

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 -axis. Indeed, a ray emitted from O will be reflected off toward the focus . This ray will then be reflected in the direction of the negative -axis by . See again Figure 4.

It is not hard to find explicit expressions for the two reflectors. If we write and , then we have the following expressions:

where we used the shifted coordinates .

B.2. Ray-Tracing Map

One can now determine the ray-tracing map corresponding to the reflector pair , . See again Figure 4.

Let be a given direction. Thus, with and . To determine the image , consider the ray emitted from the origin in the direction . The ray will encounter the first reflector at the point and will then be reflected towards the second focus . Thus, the reflected ray can be parameterized by , where is a parameter. With (B.2), this yields that the point where the reflected ray hits the second reflector is , where

The projection of this point to the plane is the value of the ray-tracing map. Thus, we have the following explicit formula for the ray-tracing map:

B.3. Solution of the Forward Problem

We can now solve the forward problem. Given the reflector pair as defined above, an input aperture and an intensity distribution , , find the output aperture and the output intensity , .

The output aperture is simply the image of under the ray-tracing map : . To find the induced intensity , use the defining property

for all Borel sets . This is an energy balance equation.

The above integral equation allows us to find an explicit expression for given or vice versa. We consider, for simplicity, the case that is contained in the left hemisphere . We can then use coordinates:

In these coordinates, the standard measure on is given by , where . Using this, (B.5) 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 .

Acknowledgment

Tilmann Glimm gratefully acknowledges the contribution of V. Oliker for the unpublished collaborative work on an algorithm similar to the one described in this paper for a different reflector design problem. This work was done during 2002 to 2004 at Emory University.

References

V. Oliker, “Mathematical aspects of design of beam shaping surfaces in geometrical optics,” in Trends in Nonlinear Analysis, M. Kirkilionis, S. Krömker, R. Rannacher, and F. Tomi, Eds., pp. 191–222, Springer, Berlin, Germany, 2002.

T. Glimm, “A rigorous analysis using optimal transport theory for a two-reflector design problem with a point source,” Inverse Problems, vol. 26, no. 4, Article ID 045001, p. 16, 2010.

Y. Brenier, “Polar factorization and monotone rearrangement of vector-valued functions,” Communications on Pure and Applied Mathematics, vol. 44, no. 4, pp. 375–417, 1991.

L. A. Caffarelli, “Allocation maps with general cost functions,” in Partial Differential Equations and Applications, vol. 177 of Lecture Notes in Pure and Applied Mathematics, pp. 29–35, Dekker, New York, NY, USA, 1996.

W. Gangbo and R. J. McCann, “Optimal maps in Monge's mass transport problem,” Comptes Rendus de l'Académie des Sciences, vol. 321, no. 12, pp. 1653–1658, 1995.

C. Villani, Optimal Transport: Old and New, vol. 338 of Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences], Springer, Berlin, Germany, 2009.

T. Glimm and V. Oliker, “Optical design of two-reflector systems, the Monge-Kantorovich mass transfer problem and Fermat's principle,” Indiana University Mathematics Journal, vol. 53, no. 5, pp. 1255–1277, 2004.

T. Glimm and V. Oliker, “Optical design of single reflector systems and the Monge-Kantorovich mass transfer problem,” Journal of Mathematical Sciences, vol. 117, no. 3, pp. 4096–4108, 2003.

X.-J. Wang, “On the design of a reflector antenna. II,” Calculus of Variations and Partial Differential Equations, vol. 20, no. 3, pp. 329–341, 2004.

T. Graf, On the near-field reector problem and optimal transport [Ph.D. thesis], Emory University, 2010.

V. Oliker, “Designing freeform lenses for intensity and phase control of coherent light with help from geometry and mass transport,” Archive for Rational Mechanics and Analysis, vol. 201, no. 3, pp. 1013–1045, 2011.

E. Haber, T. Rehman, and A. Tannenbaum, “An efficient numerical method for the solution of the optimal mass transfer problem,” SIAM Journal on Scientific Computing, vol. 32, no. 1, pp. 197–211, 2010.

J.-D. Benamou and Y. Brenier, “A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem,” Numerische Mathematik, vol. 84, no. 3, pp. 375–393, 2000.

G. L. Delzanno, L. Chacón, J. M. Finn, Y. Chung, and G. Lapenta, “An optimal robust equidistribution method for two-dimensional grid adaptation based on Monge-Kantorovich optimization,” Journal of Computational Physics, vol. 227, no. 23, pp. 9841–9864, 2008.

G. Awanou, “Spline element method for the Monge-Ampere equation,” December 2010

L. Chacón, G. L. Delzanno, and J. M. Finn, “Robust, multidimensional mesh-motion based on Monge-Kantorovich equidistribution,” Journal of Computational Physics, vol. 230, no. 1, pp. 87–103, 2011.

J.-D. Benamou, B. D. Froese, and A. M. Oberman, “Two numerical methods for the elliptic Monge-Ampère equation,” Mathematical Modelling and Numerical Analysis, vol. 44, no. 4, pp. 737–758, 2010.

B. D. Froese and A. M. Oberman, “Fast finite difference solvers for singular solutions of the elliptic Monge-Ampère equation,” Journal of Computational Physics, vol. 230, no. 3, pp. 818–834, 2011.

V. Zheligovsky, O. Podvigina, and U. Frisch, “The Monge-Ampère equation: various forms and numerical solution,” Journal of Computational Physics, vol. 229, no. 13, pp. 5043–5061, 2010.

V. I. Oliker and L. D. Prussner, “On the numerical solution of the equation and its discretizations. I,” Numerische Mathematik, vol. 54, no. 3, pp. 271–293, 1988.

J. W. Barrett and L. Prigozhin, “A mixed formulation of the Monge-Kantorovich equations,” Mathematical Modelling and Numerical Analysis, vol. 41, no. 6, pp. 1041–1060, 2007.

L. Rüschendorf and L. Uckelmann, “Numerical and analytical results for the transportation problem of Monge-Kantorovich,” Metrika, vol. 51, no. 3, pp. 245–258, 2000.

C. Canavesi, W. J. Cassarly, and J. P. Rolland, “Observations on the linear programming formulation of the single reflector design problem,” Optics Express, vol. 20, no. 4, pp. 4050–4055, 2012.

C. Canavesi, W. J. Cassarly, and J. P. Rolland, “Implementation of the linear programming algorithm for freeform reflector design,” p. 84850E-84850E-6, 2012

X.-J. Wang, “On the design of a reflector antenna,” Inverse Problem, vol. 12, no. 3, pp. 351–375, 1996.

R. Cools, “An encyclopaedia of cubature formulas,” Journal of Complexity, vol. 19, no. 3, pp. 445–453, 2003.

P.-O. Persson and G. Strang, “A simple mesh generator in Matlab,” SIAM Review, vol. 46, no. 2, pp. 329–345, 2004.

“The MOSEK optimization software,” http://www.mosek.com/

I. The MathWorks. Support pages for matlab, 2012, http://www.mathworks.com/support/solutions/en/data/1-IHYHFZ/index.html