# Iterative Scheme for Solving Optimal Transportation Problems Arising in Reflector Design

Authors: Tilmann Glimm; Nick Henscheid
Journal: ISRN Applied Mathematics
 (2013)

---

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


---

## Body

## 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 origin0 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 negativez-axis. Points in the plane perpendicular to the z-axis are denoted by the vector x∈ℝ2. Correspondingly, points in three-dimensional space are denoted by (x,z). 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ℝn with quadratic cost (also called L2 optimal transport) [12–20], although some authors have treated costs proportional to the distance (L1 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 M sample points in the input aperture (denoted by D- in Figure 1) and N points in the output aperture (T- in Figure 1), then the number of constraints is M·N. 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 T- and D- (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” Scheme5 for several edge lengths h. Memory was insufficient for the LP solver to handle the case h=0.04, 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: h: edge length; N: number of mesh points in D-; M: number of mesh points in T-; time (s): time in seconds; constr.: number of constraints.
h
M
N
Max errorR1
L
2 error R1
Max errorR2
L
2 error R2
Time (s)
Constr.
0.09
511
522
3.68393
E
-
3
1.25443
E
-
3
4.58391
E
-
3
1.91131
E
-
3
2.97016
266,742
0.08
660
663
2.96216
E
-
3
7.30738
E
-
4
3.39579
E
-
3
1.08343
E
-
3
4.67627
437,580
0.07
880
874
2.33561
E
-
3
7.43916
E
-
4
2.83982
E
-
3
1.15916
E
-
3
10.1342
769,120
0.06
1,214
1,214
1.3051
E
-
3
7.12914
E
-
4
2.10579
E
-
3
8.28341
E
-
4
17.4949
1,473,796
0.05
1,752
1,757
9.11319
E
-
4
2.66152
E
-
4
1.20934
E
-
3
4.61313
E
-
4
45.9136
3,078,264Table 2
Results of an iterative run with the “Iterative Refinement” Scheme6. 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 κnn=6, and the quasi-activity threshold was ε=0.01h2. k: iteration number; M: number of mesh points in D-; N: number of mesh points in T-; constr.: number of constraints; con. dens.: constraint density, that is, used constraints as a percentage of all possible constraints M·N; 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.
k
h
M
N
Max errorR1
L
2 error R1
Max errorR2
L
2 error R2
Constr.
Con. dens.
Time (s)
LP time (s)
0
0.064
1,059
1,059
1.2884
E
-
3
3.77819
E
-
4
1.88783
E
-
3
4.09084
E
-
4
1,121,481
100%
15.962
12.3146
1
0.0512
1,675
1,678
9.40425
E
-
4
3.36918
E
-
4
1.26944
E
-
3
3.63526
E
-
4
56,080
2.00%
208.533
0.938877
2
0.04096
2,647
2,647
7.56437
E
-
4
1.67174
E
-
4
7.53459
E
-
4
1.85922
E
-
4
89,724
1.28%
392.982
1.59831
3
0.032768
4,171
4,175
3.64902
E
-
4
9.63559
E
-
5
5.25941
E
-
4
1.21839
E
-
4
141,550
0.81%
969.239
2.54405
4
0.0262144
6,582
6,567
3.03932
E
-
4
8.8103
E
-
5
5.30343
E
-
4
1.41908
E
-
4
224,327
0.52%
2363.79
4.32524
5
0.02097152
10,313
10,318
2.49854
E
-
4
9.9045
E
-
5
3.9372
E
-
4
1.51851
E
-
4
351,514
0.33%
5909.25
9.21891We note that our proposed algorithm does not make anya 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 Section2, 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 Figure1. A point source at the origin O=(0,0,0) generates a spherical wavefront over a given input aperture D- contained in the unit sphere S2. We use the bar in the notation for D- to stress that D- 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 z-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 T-. Again, we use the notation T- to stress that T- is closed. Certain regularity conditions apply to D- and T-; see [2] for the technical details.Denote points in spaceℝ3 by pairs (x,z), where x∈ℝ2 is the position vector in a plane perpendicular to the direction of propagation and z∈ℝ 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 S2 will typically be denoted by m∈S2; their components are also written as m=(mx,mz) with |mx|2+mz2=1.We fix the output aperture in the planez=-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
(1)Reflector1:R1={ρ(m)·m∣m∈D-},Reflector2:R2={(x,z(x))∣x∈T-}.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 [1] showed that local energy conservation translates into a complicated partial differential equation of Monge-Ampère type for ρ(m). 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 functionK(m,x), called the cost function in analogy with the theory of optimal transportation, plays an important role:
(2)K(m,x)=ℓ-〈mx,x〉2ℓ(ℓ2-|x|2)(1+mz)-14ℓ2111111111form=(mx,mz)∈D-,x∈T-.In further preparation, the following two transformations are needed.Definition 1.
Letz=z(x) be a continuous function defined on T-⊆ℝ2. Then define the function
(3)z~(x)=12ℓ-z(x)ℓ2-|x|2forx∈T-.Definition 2.
Letρ=ρ(m) be a continuous function defined on D-⊆S2 with ρ>0. Then define the function
(4)ρ^(m)=-12ℓ+12ρ(m)·(mz+1)form∈D-.The following characterization of the solution(ρ,z) of the reflector problem described above was given in [2].Theorem 3.
If the pair(ρ,z) is a solution of the above reflector design problem, then the transformed pair (logρ^,logz~) minimizes the functional
(5)ℱ(r,ζ)=∫D-‍r(m)I(m)dσ+∫T-‍ζ(x)L(x)dx
on the function space Adm(D-,T-)={(r,ζ)∈C(D-)×C(T-)∣r(m)+ζ(x)≥logK(m,x)  for allm∈D-,x∈T-}.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 ρ(m) and z(x). 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(ρ,z)∈C(D-)×C(T-) be a solution of the reflector problem. Define its ray-tracing map as a set-valued map γ:D-→{subsets ofT-} via
(6)γ(m)={x∈T-∣ρ^(m)z~(x)=K(m,x)}form∈D-.In [2], it is shown that γ(m) is in fact single-valued for almost all m∈D-. Therefore, γ may be regarded as a transformation γ:D-→T-. Indeed, one can show that 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-; see [2].More details and rigorous statements are found in AppendixA.


## 3. Numerical Schemes

In the following, we describe two schemes for solving the minimization problem in Theorem3 numerically. First, we list the given data on which the solution depends.(i)
DomainsD-⊆S2, T-⊆ℝ2 (The rigorous result in [2]   required   the additional technical assumption (0,0,-1),(0,0,1)∉D-. For practical purposes, these assumptions can often be dropped.)(ii)
Nonnegative integrable functionsI(m), m∈D-, and L(x), x∈T-, with ∫D-‍Idσ=∫T-‍Ldx.(iii)
A reduced optical lengthℓ.(iv)
The valueρ(m1)=ρ1 for some fixed m1∈D-, 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 r↦r+c, ζ↦z-c for constant c will leave the objective functional unchanged.)The first numerical scheme, Scheme5 is a straightforward discretization of the minimization problem in Theorem 3 via meshes on the domains D- and T-.As described in the introduction, this scheme is memory-intensive because every pair of points fromD- and T- 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 Theorem3 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 domainD- by choosing sets d1(0),d2(0),…,dM(0)(0)⊆D-, where the interiors of any di(0) and dj(0) are disjoint (for i≠j) and D-=⋃i=1M(0)‍di(0). For instance, d1(0),d2(0),…,dM(0) may be a triangulation of D-.
(ii) Similarly, create a mesh in the output domainT- by choosing sets t1(0),t2(0),…,tN(0)(0)⊆D-, where the interiors of any ti(0) and tj(0) are disjoint (for i≠j) and T-=⋃j=1N(0)tj(0).
(iii) Choose sample pointsmi(0)∈di(0) (for i=1,…,M(0)) and xj(0)∈tj(0) (for j=1,…,N(0)). Here, we may assume that m1(0)=m1 is the same point as given in the data in (iv) above.
(iv) Find the solution({ri(0)}i=1M(0),{ζj(0)}j=1N(0)) of the following LP problem:
(7)Minimize∑i=1M(0)‍ri(0)I(mi(0))area(di(0))1111111+∑j=1N(0)‍ζj(0)L(xj(0))area(tj(0))(8)subjecttor1(0)=logρ^1,whereρ^1isasgiveninthedatain(iv)above,and(9)ri(0)+ζj(0)≥logK(mi(0),xj(0))fori=1,…,M(0),j=1,…,N(0).(v)
Find the numbersρi(0),i=1,…,M(0) and zj(0),j=1,…,N(0) such that ri(0)=logρ^i(0) and ζj(0)=logz~j(0). This is straightforward by taking the inverse of the transformations given in Definition 1 and Scheme 5.Thenρi(0) is an approximation for the true value ρ(mi(0)) of the radial function of the first reflector, evaluated at the sample point mi(0) for i=1,…,M(0). Similarly, zj(0) is an approximation for the true value z(xj(0)) of the function describing the second reflector, evaluated at the sample point xj(0) 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.Note that (7) is a discretization of the sum of integrals in the functional ℱ(r,ζ)=∫D-‍r(m)I(m)dσ+∫T-‍ζ(x)L(x)dx. 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 di and tj 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 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
(10)ri(0)+ζj*(0)=logK(mi(0),xj*(0)).
This means that the point mj*(0) is approximately the image of the point xi(0) under the ray-tracing map γ.


### 3.2. “Iterative Refinement” Scheme

As noted in the introduction, one of the drawbacks of Scheme5 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 (M(0)≈1750,N(0)≈1750).We addressed this problem by developing an iterative scheme. First, the problem is solved for a mesh with relatively few sample points (sayM(0)=N(0)=1000). 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 “quasi-activity threshold,” and an integer κnn, the nearest neighbor search depth, to be explained in detail below.Scheme 6 (iterative refinement).
Use Scheme5 to find an initial solution ({ri(0)}i=1M(0),{ζj(0)}j=1N(0)) for M(0) sample points on D- and N(0) sample points on T-, respectively.Definition 7.
Call a pair(i,j) with i∈{1,…,M(0)} and j∈{1,…,N(0)} “quasi-active” if
(11)r(mj(0))+ζ(xj(0))-logK(mi(0),xj(0))≤ε(0).
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 onD- by choosing sets d1(1),d2(1),…,dM(1)(1)⊆D-, where the interiors of any di(1) and dj(1) are disjoint (for i≠j) and D-=⋃i=1M(1)di(1). Here, M(1) is chosen larger than M(0), meaning that we have a finer mesh.(ii)
Similarly, create a second mesh onT- by choosing sets t1(1),t2(1),…,tN(1)(1)⊆D-, where the interiors of any ti(1) and tj(1) are disjoint (for i≠j) and T-=⋃j=1N(1)tj(1).(iii)
Choose sample pointsmi(1)∈di(1) (for i=1,…,M(1)) and xi(1)∈tj(1) (for j=1,…,N(1)). Here, we may assume that m1(1)=m1 is the same point as given in the data in (iv) above.(iv)
For each pair(i,j) with i=1,…,M(1), and j=1,…,M(1), find the κnn nearest neighbors of mi(1) from the previous mesh on D- and the κnn nearest neighbors of xj(1) from the previous mesh on T-.Definition 8.
Say that the pair(i,j) is potentially active if at least one of the κnn2 pairs consisting of a nearest neighbor of mi(1) and a nearest neighbor of xj(1) is quasi-active in the above sense. (i)
Find the solution({ri}i=1M(1),{ζj(1)}j=1N(1)) of the following LP problem:
(12)Minimize∑i=1M(1)‍ri(1)I(mi(1))area(di(1))1111111111111+∑j=1N(1)‍ζi(1)L(xj(1))area(tj(1))subject  tor1(1)=logρ^1,1111111111ri(1)+ζj(1)≥logK(mi(1),xj(1))11111111111111for  any  potentially  active  pair11111111111111111111(in  the  above  sense).(ii)
Find the numbersρi(1),i=1,…,M(1) and zj(1),j=1,…,N(1) such that ri(1)=logρ^i(1) and ζj(1)=logz~j(1).The idea of Scheme6 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 D- and T-, 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 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 κnn 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 r(mj(0))+ζ(xj(0))-logK(mi(0),xj(0)) is zero) and those which are “almost active” in the sense that the term r(mj(0))+ζ(xj(0))-logK(mi(0),xj(0)) is positive, but small. “Smallness” is encoded in the quasi-activity threshold ε(1). 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 functionsI(m) and L(x) so that the approximations for the integrals ∫D-‍Idσ and ∫T-‍Ldx yield exactly the same result.Note that the algorithm in Scheme6 can be iterated again. We can use the first iterative solution ρi(1),i=1,…,M(1) and zj(1),j=1,…,N(1), in Scheme 6, to obtain a second iterative solution ρi(2),i=1,…,M(2) and zj(2),j=1,…,N(2), and so on.Note that (7) is a discretization of the sum of integrals in the functional ℱ(r,ζ)=∫D-‍r(m)I(m)dσ+∫T-‍ζ(x)L(x)dx. We note that we have chosen again the simplest discretization here for the sum of integrals ∫D-‍r(m)I(m)dσ+∫T-‍ζ(x)L(x)dx. 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 Scheme6

For the iteration of Scheme6 described above, we have to choose a corresponding sequence of mesh sizes {(M(0),N(0)),(M(1),N(1)),(M(2),N(2)),…}. We also need to choose a corresponding sequence of thresholds ε(1),ε(2),ε(3),….This last question is of a certain practical importance. Indeed, if the thresholds{ε(k)}k 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 {ε(k)}k 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 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.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 mi(k) in the input aperture D- is paired via the ray-tracing map γ with a unique point xj(i)(k)≈γ(mi(k)) and vice versa. Thus, there are M(k) such pairs. Let r(m),m∈D- and ζ(x),x∈T- denote the solutions of the reflector problem. Then the set of all points in (m,x)-space close to a given pair (mi(k),xj(i)(k)) that satisfy r(x)-z(x)-logK(m,x)<ε(k) is approximately an ellipsoid for small ε(k). This can be seen by using the Taylor expansion around the pair (mi(k),xj(i)(k)) in (m,x)-space. The volume of this ellipsoid is proportional to (ε(k))2. So the proportion of pairs of points that satisfy the above constraint among all constraints is roughly proportional to (ε(k))2. Since each pair (mi(k),xj(k)) that satisfies this inequality corresponds to a constraint included in the LP of the kth iteration, the total number of constraints should be roughly proportional to (M(k))3·(ε(k))2. So to keep the number of constraints growing linearly throughout the iterations, based on these heuristics, one may choose
(13)ε(k)=CM(k),
where C is a constant. These heuristics are very rough. In practice, we have found that when using the formula ε(k)=C(M(k))α, there is a range of values of α which all produce similar results. See Section 4.2.


## 3.1. “Simple Discretization” Scheme

The characterization of solutions to the reflector problem by means of the minimization of a linear functional in Theorem3 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 domainD- by choosing sets d1(0),d2(0),…,dM(0)(0)⊆D-, where the interiors of any di(0) and dj(0) are disjoint (for i≠j) and D-=⋃i=1M(0)‍di(0). For instance, d1(0),d2(0),…,dM(0) may be a triangulation of D-.
(ii) Similarly, create a mesh in the output domainT- by choosing sets t1(0),t2(0),…,tN(0)(0)⊆D-, where the interiors of any ti(0) and tj(0) are disjoint (for i≠j) and T-=⋃j=1N(0)tj(0).
(iii) Choose sample pointsmi(0)∈di(0) (for i=1,…,M(0)) and xj(0)∈tj(0) (for j=1,…,N(0)). Here, we may assume that m1(0)=m1 is the same point as given in the data in (iv) above.
(iv) Find the solution({ri(0)}i=1M(0),{ζj(0)}j=1N(0)) of the following LP problem:
(7)Minimize∑i=1M(0)‍ri(0)I(mi(0))area(di(0))1111111+∑j=1N(0)‍ζj(0)L(xj(0))area(tj(0))(8)subjecttor1(0)=logρ^1,whereρ^1isasgiveninthedatain(iv)above,and(9)ri(0)+ζj(0)≥logK(mi(0),xj(0))fori=1,…,M(0),j=1,…,N(0).(v)
Find the numbersρi(0),i=1,…,M(0) and zj(0),j=1,…,N(0) such that ri(0)=logρ^i(0) and ζj(0)=logz~j(0). This is straightforward by taking the inverse of the transformations given in Definition 1 and Scheme 5.Thenρi(0) is an approximation for the true value ρ(mi(0)) of the radial function of the first reflector, evaluated at the sample point mi(0) for i=1,…,M(0). Similarly, zj(0) is an approximation for the true value z(xj(0)) of the function describing the second reflector, evaluated at the sample point xj(0) 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.Note that (7) is a discretization of the sum of integrals in the functional ℱ(r,ζ)=∫D-‍r(m)I(m)dσ+∫T-‍ζ(x)L(x)dx. 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 di and tj 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 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
(10)ri(0)+ζj*(0)=logK(mi(0),xj*(0)).
This means that the point mj*(0) is approximately the image of the point xi(0) under the ray-tracing map γ.


## 3.2. “Iterative Refinement” Scheme

As noted in the introduction, one of the drawbacks of Scheme5 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 (M(0)≈1750,N(0)≈1750).We addressed this problem by developing an iterative scheme. First, the problem is solved for a mesh with relatively few sample points (sayM(0)=N(0)=1000). 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 “quasi-activity threshold,” and an integer κnn, the nearest neighbor search depth, to be explained in detail below.Scheme 6 (iterative refinement).
Use Scheme5 to find an initial solution ({ri(0)}i=1M(0),{ζj(0)}j=1N(0)) for M(0) sample points on D- and N(0) sample points on T-, respectively.Definition 7.
Call a pair(i,j) with i∈{1,…,M(0)} and j∈{1,…,N(0)} “quasi-active” if
(11)r(mj(0))+ζ(xj(0))-logK(mi(0),xj(0))≤ε(0).
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 onD- by choosing sets d1(1),d2(1),…,dM(1)(1)⊆D-, where the interiors of any di(1) and dj(1) are disjoint (for i≠j) and D-=⋃i=1M(1)di(1). Here, M(1) is chosen larger than M(0), meaning that we have a finer mesh.(ii)
Similarly, create a second mesh onT- by choosing sets t1(1),t2(1),…,tN(1)(1)⊆D-, where the interiors of any ti(1) and tj(1) are disjoint (for i≠j) and T-=⋃j=1N(1)tj(1).(iii)
Choose sample pointsmi(1)∈di(1) (for i=1,…,M(1)) and xi(1)∈tj(1) (for j=1,…,N(1)). Here, we may assume that m1(1)=m1 is the same point as given in the data in (iv) above.(iv)
For each pair(i,j) with i=1,…,M(1), and j=1,…,M(1), find the κnn nearest neighbors of mi(1) from the previous mesh on D- and the κnn nearest neighbors of xj(1) from the previous mesh on T-.Definition 8.
Say that the pair(i,j) is potentially active if at least one of the κnn2 pairs consisting of a nearest neighbor of mi(1) and a nearest neighbor of xj(1) is quasi-active in the above sense. (i)
Find the solution({ri}i=1M(1),{ζj(1)}j=1N(1)) of the following LP problem:
(12)Minimize∑i=1M(1)‍ri(1)I(mi(1))area(di(1))1111111111111+∑j=1N(1)‍ζi(1)L(xj(1))area(tj(1))subject  tor1(1)=logρ^1,1111111111ri(1)+ζj(1)≥logK(mi(1),xj(1))11111111111111for  any  potentially  active  pair11111111111111111111(in  the  above  sense).(ii)
Find the numbersρi(1),i=1,…,M(1) and zj(1),j=1,…,N(1) such that ri(1)=logρ^i(1) and ζj(1)=logz~j(1).The idea of Scheme6 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 D- and T-, 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 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 κnn 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 r(mj(0))+ζ(xj(0))-logK(mi(0),xj(0)) is zero) and those which are “almost active” in the sense that the term r(mj(0))+ζ(xj(0))-logK(mi(0),xj(0)) is positive, but small. “Smallness” is encoded in the quasi-activity threshold ε(1). 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 functionsI(m) and L(x) so that the approximations for the integrals ∫D-‍Idσ and ∫T-‍Ldx yield exactly the same result.Note that the algorithm in Scheme6 can be iterated again. We can use the first iterative solution ρi(1),i=1,…,M(1) and zj(1),j=1,…,N(1), in Scheme 6, to obtain a second iterative solution ρi(2),i=1,…,M(2) and zj(2),j=1,…,N(2), and so on.Note that (7) is a discretization of the sum of integrals in the functional ℱ(r,ζ)=∫D-‍r(m)I(m)dσ+∫T-‍ζ(x)L(x)dx. We note that we have chosen again the simplest discretization here for the sum of integrals ∫D-‍r(m)I(m)dσ+∫T-‍ζ(x)L(x)dx. 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 Scheme6

For the iteration of Scheme6 described above, we have to choose a corresponding sequence of mesh sizes {(M(0),N(0)),(M(1),N(1)),(M(2),N(2)),…}. We also need to choose a corresponding sequence of thresholds ε(1),ε(2),ε(3),….This last question is of a certain practical importance. Indeed, if the thresholds{ε(k)}k 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 {ε(k)}k 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 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.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 mi(k) in the input aperture D- is paired via the ray-tracing map γ with a unique point xj(i)(k)≈γ(mi(k)) and vice versa. Thus, there are M(k) such pairs. Let r(m),m∈D- and ζ(x),x∈T- denote the solutions of the reflector problem. Then the set of all points in (m,x)-space close to a given pair (mi(k),xj(i)(k)) that satisfy r(x)-z(x)-logK(m,x)<ε(k) is approximately an ellipsoid for small ε(k). This can be seen by using the Taylor expansion around the pair (mi(k),xj(i)(k)) in (m,x)-space. The volume of this ellipsoid is proportional to (ε(k))2. So the proportion of pairs of points that satisfy the above constraint among all constraints is roughly proportional to (ε(k))2. Since each pair (mi(k),xj(k)) that satisfies this inequality corresponds to a constraint included in the LP of the kth iteration, the total number of constraints should be roughly proportional to (M(k))3·(ε(k))2. So to keep the number of constraints growing linearly throughout the iterations, based on these heuristics, one may choose
(13)ε(k)=CM(k),
where C is a constant. These heuristics are very rough. In practice, we have found that when using the formula ε(k)=C(M(k))α, 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 AppendixB. The data set of this section is obtained from the more general solution of Appendix B by choosing a=b=0 and c=-0.4, α=1, and R=1.3. We consider the input aperture
(14)Ω-={(mx,my,mz)⊆S2∣mx2+my2≤0.8,mz<0},
which is a spherical cap centered at the point (0,0,-1). 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:
(15)T-={(x,y)∈ℝ2∣x2+y2≤179}.
We choose a constant output intensity
(16)L(x)=1forx∈T-.
Using the relation (B.7), this gives the input intensity
(17)I(m)=14.2716049383(1-mz)2form=(mx,my,mz)∈Ω-.
The two reflectors are given as follows:
(18)radial radiusρforR1:ρ(m)=0.7651.3+0.4mz,(19)equationforR2:z=-0.25(x2+y2)+0.6.
The corresponding reduced optical path length is
(20)ℓ=R+2α=2.9.
The two reflectors are shown in Figure 2. R1 is part of a spheroid, and R2 is part of a paraboloid; see Appendix B for details.Plots of the reflectors given by (18) (a) and (19) (b) in the explicit data set given in Section 4.1.
(a)
(b)


### 4.2. Results

Using the data set from Section4.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 ρ(m) and z(x) 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 Schemes5 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 h, not the number of mesh points M (for the domain Ω-) or N (for the domain T-), respectively. However, it can be seen from geometric considerations that the relation between the average edge length h and the number of mesh points M on D- obeys
(21)h~1M.
In the tables below, we show both the number of mesh points M,N and the edge length h, as we believe both are informative.Our primary aim is to demonstrate the practicability of the proposed Scheme6. 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 Scheme6, we chose a sequence of desired edge lengths h(0),h(1),… by reducing h by 20% in each step; that is, h(k)=0.8kh(0). Then, based on the heuristics given above and in Section 3.3, we determined the corresponding constraint thresholds with the following   formula:
(22)ε(k)=C·(h(k))μ,(k=1,2,3,…),
where C and μ are constants. The heuristic arguments in Section 3.3 yield μ=2, but we also tested other values. 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. In contrast, “small” values for C 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 C=0, 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 C and μ, meaning for a broad range of values, we reached the terminal iteration with h=0.02097152; 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 I(m) and L(x), specifically to know whether the algorithm is stable under small perturbations of I(m) and L(x). For this, we added small perturbation δcos(my) to I(m) and δcos(x) to L(x), then renormalized L(x) 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 Table2. 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 Kij, which is of size M(k)×N(k). 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 toMσ, or equivalently proportional to h2σ, where σ=-0.773, data from Table 2. The same formulation for scheme 5 yields an exponent of σ=-1.181, unsurprisingly a slightly more rapid decrease, data from Table 1. The geometric distribution of errors is indicated in Figure 3.Error surface plots for reflector 1 (a) and reflector 2 (b) for the data of Table3 with κnn=2. 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.
(a)
(b)To investigate the role of the nearest neighbor search depthκnn, we performed several tests with different values. The results are summarized in Table 3. The value of κnn 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 κnn, ranging from 96,256 (0.09% of all possible constraints) for κnn=2 to 48,2058 (0.45% of all possible constraints) for κnn=8 in this example. The time to solve the LP problems, thus, also varies from 2.82 seconds for κnn=2 to 12.91 seconds for κnn=8, 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” Scheme6 with different values of the nearest neighbor search depth κnn. All other data used was the same as in Table 2. Note that the nearest neighbor search depth κnn 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.
κ
nn
h
M
N
Max errorR1
L
2 error R1
Max errorR2
L
2 error R2
Constr.
Con. dens.
Time (s)
LP time (s)
2
0.02097152
10,313
10,318
2.71973
E
-
4
9.04017
E
-
5
4.38374
E
-
4
1.37601
E
-
4
96256
0.090%
5947.76
2.82425
4
0.02097152
10,313
10,318
2.48492
E
-
4
9.77635
E
-
5
3.92226
E
-
4
1.50005
E
-
4
226006
0.212%
5752.17
6.09849
6
0.02097152
10,313
10,318
2.49854
E
-
4
9.9045
E
-
5
3.9372
E
-
4
1.51851
E
-
4
351514
0.330%
5909.25
9.21891
8
0.02097152
10,313
10,318
2.49816
E
-
4
9.90729
E
-
5
3.94393
E
-
4
1.5188
E
-
4
482058
0.453%
6322.4
12.9115


## 4.1. A Data Set

We have constructed a particular data set using an explicit analytical solution described in more detail in AppendixB. The data set of this section is obtained from the more general solution of Appendix B by choosing a=b=0 and c=-0.4, α=1, and R=1.3. We consider the input aperture
(14)Ω-={(mx,my,mz)⊆S2∣mx2+my2≤0.8,mz<0},
which is a spherical cap centered at the point (0,0,-1). 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:
(15)T-={(x,y)∈ℝ2∣x2+y2≤179}.
We choose a constant output intensity
(16)L(x)=1forx∈T-.
Using the relation (B.7), this gives the input intensity
(17)I(m)=14.2716049383(1-mz)2form=(mx,my,mz)∈Ω-.
The two reflectors are given as follows:
(18)radial radiusρforR1:ρ(m)=0.7651.3+0.4mz,(19)equationforR2:z=-0.25(x2+y2)+0.6.
The corresponding reduced optical path length is
(20)ℓ=R+2α=2.9.
The two reflectors are shown in Figure 2. R1 is part of a spheroid, and R2 is part of a paraboloid; see Appendix B for details.Plots of the reflectors given by (18) (a) and (19) (b) in the explicit data set given in Section 4.1.
(a)
(b)


## 4.2. Results

Using the data set from Section4.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 ρ(m) and z(x) 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 Schemes5 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 h, not the number of mesh points M (for the domain Ω-) or N (for the domain T-), respectively. However, it can be seen from geometric considerations that the relation between the average edge length h and the number of mesh points M on D- obeys
(21)h~1M.
In the tables below, we show both the number of mesh points M,N and the edge length h, as we believe both are informative.Our primary aim is to demonstrate the practicability of the proposed Scheme6. 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 Scheme6, we chose a sequence of desired edge lengths h(0),h(1),… by reducing h by 20% in each step; that is, h(k)=0.8kh(0). Then, based on the heuristics given above and in Section 3.3, we determined the corresponding constraint thresholds with the following   formula:
(22)ε(k)=C·(h(k))μ,(k=1,2,3,…),
where C and μ are constants. The heuristic arguments in Section 3.3 yield μ=2, but we also tested other values. 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. In contrast, “small” values for C 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 C=0, 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 C and μ, meaning for a broad range of values, we reached the terminal iteration with h=0.02097152; 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 I(m) and L(x), specifically to know whether the algorithm is stable under small perturbations of I(m) and L(x). For this, we added small perturbation δcos(my) to I(m) and δcos(x) to L(x), then renormalized L(x) 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 Table2. 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 Kij, which is of size M(k)×N(k). 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 toMσ, or equivalently proportional to h2σ, where σ=-0.773, data from Table 2. The same formulation for scheme 5 yields an exponent of σ=-1.181, unsurprisingly a slightly more rapid decrease, data from Table 1. The geometric distribution of errors is indicated in Figure 3.Error surface plots for reflector 1 (a) and reflector 2 (b) for the data of Table3 with κnn=2. 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.
(a)
(b)To investigate the role of the nearest neighbor search depthκnn, we performed several tests with different values. The results are summarized in Table 3. The value of κnn 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 κnn, ranging from 96,256 (0.09% of all possible constraints) for κnn=2 to 48,2058 (0.45% of all possible constraints) for κnn=8 in this example. The time to solve the LP problems, thus, also varies from 2.82 seconds for κnn=2 to 12.91 seconds for κnn=8, 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” Scheme6 with different values of the nearest neighbor search depth κnn. All other data used was the same as in Table 2. Note that the nearest neighbor search depth κnn 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.
κ
nn
h
M
N
Max errorR1
L
2 error R1
Max errorR2
L
2 error R2
Constr.
Con. dens.
Time (s)
LP time (s)
2
0.02097152
10,313
10,318
2.71973
E
-
4
9.04017
E
-
5
4.38374
E
-
4
1.37601
E
-
4
96256
0.090%
5947.76
2.82425
4
0.02097152
10,313
10,318
2.48492
E
-
4
9.77635
E
-
5
3.92226
E
-
4
1.50005
E
-
4
226006
0.212%
5752.17
6.09849
6
0.02097152
10,313
10,318
2.49854
E
-
4
9.9045
E
-
5
3.9372
E
-
4
1.51851
E
-
4
351514
0.330%
5909.25
9.21891
8
0.02097152
10,313
10,318
2.49816
E
-
4
9.90729
E
-
5
3.94393
E
-
4
1.5188
E
-
4
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 Table2. 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 Table2 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.



---

*Source: 635263-2013-11-24.xml*