Guaranteed ε-optimal solutions with the linear optimizer ART3+O
aa r X i v : . [ phy s i c s . m e d - ph ] J a n Guaranteed ε -optimal solutions with the linearoptimizer ART3+O Bram L. Gorissen
Abstract
The linear optimization algorithm ART3+O introduced by Chen et al. (2010) canefficiently solve large scale inverse planning problems encountered in radiation therapyby iterative projection. Its major weakness is that it cannot guarantee ε -optimality ofthe final solution due to an arbitrary stopping criterion. We propose an improvementto ART3+O where the stopping criterion is based on Farkas’ lemma. The same theorycan be used to detect inconsistency in other projection methods as well. The proposedalgorithm guarantees to find an ε -optimal solution in finite time. The algorithm isdemonstrated on numerical examples in radiation therapy. This is an author-created, un-copyedited version of an article published in Physics in Medicine and BiologyDOI:10.1088/1361-6560/ab0a2e.
A linear optimization problem is a linear objective subject to linear constraints:max x { c T x : Ax ≤ b, x ≥ } . (1)It is easy to show that any linear optimization problem can be written in this form, evenif it is originally a minimization problem, has equality constraints or allows variablesto be negative.One application of linear optimization is inverse planning in radiation therapy,where the goal is to customize a treatment plan to the patient’s anatomy. A 3d grid issuperimposed onto a patient for discretization purposes. The cubes of this grid (voxels)that lie in or near the tumor (the planning target volume, PTV) get prescribed a certaindose level, while the other voxels, especially those in critical organs, should get as littledose as possible. In many treatment modalities, the dose in a voxel is linear in thefluence vector x ∈ R n (Bortfeld et al., 1990; Lessard and Pouliot, 2001; Lomax, 1999): d i = X j D ij x j . (2)The coefficients D ij depend on the patient’s anatomy and can be computed analyticallyor with monte carlo simulation. Surprisingly many dose objectives and constraintscan be expressed in the framework of linear optimization. Mean and maximum dose(Bahr et al., 1968), mean dose above or under a given threshold dose (Shepard et al., A , a vector b , and iteratively changes an arbitrarystarting point x until Ax ≤ b . In iteration k , it considers the i th row of A , where i = 1 , , . . . , m, , , . . . , m, , , etc. If the i th inequality is violated, x is reflected intothe hyperplane A i x = b . Figure 1 depicts three iterations of ART3, and shows thata feasible point is found after two reflections. If the interior of the feasible set is fulldimensional, ART3 finds a feasible point in a finite number of steps (vs in the limit)(Herman, 1975). The key advantages of ART3 are its low memory requirements, fastpractical running time, and simplicity.The performance of ART3 has been improved by temporarily removing nonviolatedconstraints from the inner loop (Alg. 2). This improved algorithm, ART3+, benefitsfrom not having to check for a violation for those constraints that are not likely tobe violated, which yields a significant speed-up (Herman and Chen, 2008). Similarto ART3, ART3+ terminates in a finite number of steps if the feasible set is fulldimensional.The algorithm ART3+O is an optimization algorithm that uses ART3+ as itsworkhorse. Given a linear optimization problem (1) and an interval [ L, U ] that isknown to contain the optimal objective value, ART3+O uses bisection search to shrinkthe interval until it finds a ε -optimal value of (1). A description is given in Alg. 3. Ineach iteration, ART3+O determines the midpoint of the interval M = ( L + U ) / M by calling ART3+. If ART3+ finds such a feasible point, the search continues on theinterval [ c T x, U ], otherwise it continues on [ L, M ], and repeats until the width of theinterval is less than ε . The x ART O that is found, is at most ε from the optimum,meaning that c T x ART O ≥ c T x − ε for all x ≥ Ax ≤ b .A practical problem with ART3+O is that ART3+ runs indefinitely when there isno x that satisfies Ax ≤ b . Therefore, the algorithm needs to be stopped prematurelywhen it becomes unlikely that a feasible point will be found, e.g., after a large number ata: A ∈ R m × n , b ∈ R m Result: x such that Ax ≤ b Let the starting point x ∈ R n be an arbitrary vector;0 → k ; while Ax (cid:2) b do k mod m ) → i ; A i x − b i → v i ; if v i > then x − v i || A i || A Ti → x end k + 1 → k ; endAlgorithm 1: The image reconstruction algorithm ART3 (Herman, 1975).
Data: A ∈ R m × n , b ∈ R m Result: x such that Ax ≤ b Let the starting point x ∈ R n be an arbitrary vector; while Ax (cid:2) b do { , , . . . , m } → S ; while S = ∅ do Let i be the next element in S ; A i x − b i → v i ; if v i > then x − v i || A i || A Ti → x else S \{ i } → S endendend Algorithm 2: The algorithm ART3+ (Herman and Chen, 2008).3 ata: A ∈ R m × n , b ∈ R m , c ∈ R n , L ∈ R , U ∈ R Result: x that solves (1) while U − L > ε do ( L + U ) / → M ;Call ART3+ to find x such that − c T A − I x ≤ − Mb ; if ART3+ found such an x then c T x → L else M → U endend Algorithm 3: The algorithm ART3+O (Chen et al., 2010). x x x Figure 1: Three iterations of ART3. The feasible region is shaded while the constraints aredashed. The starting point x is reflected into the first constraint to get to x , then into thesecond constraint to get to x . Since x is in the feasible region, the algorithm terminates.4 f iterations, after a certain runtime, or when the iterates show a certain behavior(De Pierro and Iusem, 1985; Censor and Tom, 2003). This carries the risk of incor-rectly concluding that an objective value of M is unachievable. An ε -optimal solutioncan therefore not be guaranteed. In addition, the user is faced with a difficult trade-off:using a conservative stopping criterion means many cpu cycles are wasted on findinga solution that does not exist whereas a more aggressive criterion means that the finalsolution can be far from ε -optimal. What makes matters worse is that there is noguidance to what constitutes a “conservative” stopping criterion. As of today there isno known upper bound on the number of iterations ART3+ may take to find a feasiblepoint.This paper proposes an improvement to ART3+O so that it can accurately detect ifan objective value of M is attainable. It eliminates the need for an arbitrary stoppingcriterion, and guarantees an ε -optimal solution. With this improvement, the algorithmstill terminates in a finite number of steps. The cornerstone of the proposed method is Farkas’ lemma. Although this lemma canbe stated in many different ways, our results can be derived from the following form:
Theorem 1 (Farkas’ lemma) Exactly one of the following statements is true (Farkas,1902):1. There exists an x ∈ R n such that F x = b and x ≥ .2. There exists a y ∈ R m such that F T y ≥ and b T y ≤ − . This lemma is a statement of alternatives. If the first statement is true, the secondstatement is not true and vice versa.The algorithm ART3+O calls ART3+ to find a point x in the set { x : c T x ≥ M, Ax ≤ b, x ≥ } . (3)Such a point exists if and only if there exists s ∈ R and s ∈ R m such that: (cid:18) − c T OA I (cid:19) xs s = (cid:18) − Mb (cid:19) , x ≥ , s ≥ , s ≥ . We use these equations for the first statement of Farkas’ lemma. The equations in thesecond statement are: − cy + A T y ≥ , − M y + b T y ≤ − , y ≥ , y ≥ . (4)Therefore, the set (3) is empty if and only if there is a y that satisfies (4).Currently, ART3+O cannot conclude with certainty that the set (3) is empty. Itpresumes that the set is empty when ART3+ cannot find a point in that set withinthe allocated runtime or number of iterations, but that presumption may be wrong. ata: A ∈ R m × n , b ∈ R m , c ∈ R n , L ∈ R , U ∈ R Result: x that solves (1) while U − L > ε do ( L + U ) / → M ;Call two instances of ART3+ in parallel:1. to find x such that − c T A − I x ≤ − Mb
2. to find y such that c − A T − M b T − O − I (cid:18) y y (cid:19) ≤ − ;When one instance of ART3+ finds a solution, terminate the other instance; if a feasible x was found then c T x → L else M → U endend Algorithm 4: The suggested improved algorithm.
We propose to run ART3+ both on the set (3) and on the inequalities (4) inparallel, until it finds a point in either set. When it finds a point that satisfies (4),that is mathematical proof that the set (3) is empty and that an objective value of M is unachievable. Therefore, bisection search can continue on the interval [ L, M ]. Thisalgorithm is outlined in Alg. 4.ART3+O requires that the rows of the constraint matrix A can be accessed effi-ciently, while the improved algorithm also requires efficient access to the columns of A .If only a single copy of A is kept in memory, one of the two instances of ART3+ cannothave sequential data access, causing unnecessary slowness. In radiation therapy, thisis an issue because the matrix D in (2) is often stored in sparse format, meaning thatonly the indices and values of nonzero elements get stored. For sparse formats, eitherthe rows or the columns are efficiently accessible, but not both. Therefore, the D ma-trix needs to be transposed. Consequently, the memory requirements of the improvedalgorithm are roughly twice as high as those of ART3+O. The difference in runtimeis harder to predict. Although Alg. 4 runs two threads in parallel and therefore usestwice as much cpu power when a feasible point x is found, it is possible that the secondthread terminates early enough when no feasible x exists to make up for this.Since ART3+ terminates in a finite number of steps if the feasible set is full di-mensional (Herman and Chen, 2008), it trivially follows that Alg. 4 also terminatesin a finite number of steps. The requirement of full dimensionality is not an actualrestriction: techniques from the ellipsoid method can be used to slightly perturb the nequalities in a way that the feasible region becomes full dimensional if it is not empty,while it remains empty if the original feasible region is empty (Bland et al., 1981). Nev-ertheless, we shall validate the full dimensionality assumption in our numerical resultsby subtracting a small number (10 − ) from the right-hand side of the inequalities givento ART3+ by Alg. 4 and show that a point satisfying those inequalities exists.We have taken the implementation of ART3+O from the in-house treatment plan-ning system Astroid and modified it to make it accurately correspond to Alg. 3, witha stopping criterion of 2 · iterations and ε = 0 . D matrix has 288,610 nonzero elements ( ≈ M to avoid unnecessary overhead. So, the optimization problem is max x ≥ ,t ≥ { t : D x ≥ t, D x ≤ e } , but the inequalities given to ART3+ in Alg. 4 are { D x ≥ M e, D x ≤ e, x ≥ } and {− D y + D y ≤ , M e T y − e T y ≤ − , y ≥ } .The method is also tested on a larger case with 5384 pencil beams and approxi-mately 2 · constraints that limit the mean, minimum or maximum dose, and anobjective to minimize the maximum dose in a structure with 3 · voxels. The objec-tive was treated in the same way as the small case.The method is tested on another 77 clinical test cases spanning 11 patients withdifferent tumor sites. Each patient has a common set of constraints, but differentobjectives (at most 10 per patient). In addition to minimum and maximum doseobjectives and constraints, some cases have mean dose constraints or a mean underdoseobjective. The D matrices vary in size between 12 MB and 2.5 GB. These cases will berun on the ERISone cluster (Partners Healthcare, Boston, MA, USA), which uses IntelE5-2690 v3 and Intel X5670 CPUs, with a time limit of 1 day of cpu-time (around 12hours of wall-clock time).The first two cases are run on a machine with two Intel E5-2687W v3 cpus. CPLEX12.8, a state of the art solver for linear optimization, is used to obtain the optimalsolution for both cases and to validate the full dimensionality assumption. Althoughthe second case is large and CPLEX runs out of memory if it attempts to solve it, it cansolve an equivalent reformulation. The case has one D matrix per beam, D and D ,with corresponding x and x , and constraints are formulated on D x , on D x , andon D x + D x . By adding auxiliary variables d = D x and d = D x , we avoidthat rows of D k get duplicated across constraints, which resolves the memory issue. Toresolve the memory issue when trying to validate the full dimensionality assumptionfor the y -space, we take the feasibility problem for x and turn it into an optimizationproblem by adding the objective of maximizing − P x i . The dual variables of theinequality constraints then satisfy cy − A T y ≤ − b T y <
0, and y ≥ L U
Farkas Time (s) Iterations1 0.00 50.00 1 0.1 62,3562 25.01 50.00 2 1.9 24,948,3273 25.01 37.51 1 0.0 332,4984 31.26 37.51 2 11.1 147,789,6185 31.26 34.38 2 3160.4 39,374,210,6166 31.26 32.82 1 0.0 53,6627 32.38 32.82 1 0.1 2,477,7878 32.60 32.82 1 1345.9 29,435,448,9469 32.71 32.82 2 4618.7 57,541,614,242Table 1: Results of Alg. 4 on the small radiation therapy case. [
L, U ] is the search intervalof bisection search, Farkas indicates whether the first or second statement in Farkas’ lemmawas found to be true, and the last two columns show the runtime statistics of the instanceof ART3+ that found a solution.
We then slightly increase each element of y such that y > CPLEX solves the small case in 0.8 seconds and shows that the optimal value is 32.71Gy. We confirmed that the feasible set is full dimensional in the first instance ofART3+ as long as M is less than the optimal value, while it is full dimensional in thesecond instance of ART3+ when M is larger than the optimal value. ART3+O takes12.5 seconds to obtain a solution with value 32.57. Table 1 shows the steps taken byAlg. 4 during its runtime of around 2.5 hours. In steps 2, 4, 5 and 9, it was proven thatthe objective values 37.51, 34.38, 32.82 and 32.77 are unachievable, while in steps 1,3, 6, 7 and 8 the objective value was gradually improved to 32.71. Although this finalsolution is guaranteed to be within 0.1 Gy of the optimum, it is in fact only 0.0015Gy from the optimum. Looking at the number of iterations, there are two remarkableobservations. The first is that it takes a considerable number of iterations to concludethat a value of M is unattainable. Even for M = 37 . · iterations are necessary. The second is that close to optimality,even the regular ART3+ may take a large number of iterations. In Step 8, 29 · iterations are needed to obtain a solution for M = 32 .
71. This is orders of magnitudehigher than the iteration limit of 2 · used by Chen et al. (2010).CPLEX (barrier) solves the large case in just under 3.5 hours with an optimal valueof 50.8 Gy. ART3+O finds a value of 51.6 Gy in 847.4 seconds. Alg. 4 takes 8 secondsto transpose the D matrix. In the two steps of bisection search, the second instance ofART3+ proves that M = 26 .
47 and M = 39 .
70 are unattainable in approximately 6minutes and 6 hours, respectively. We terminated bisection search in the third step ofbisection search, because even after 100 hours and 218 · iterations, ART3+ could notprove that M = 46 .
32 is unattainable, despite that the full dimensionality assumption olds for this value of M .For the test set of 77 cases, ART3+O does not find a feasible point within theiteration limit for half the patients, affecting 43 of the optimization runs, and incor-rectly concludes that the problem is infeasible. We therefore re-ran ART3+O with aniteration limit of 10 instead of 2 · , affecting the runtimes by a median factor of 26.Still, 2 out of 10 patients are incorrectly classified as infeasible, affecting 13 optimiza-tion runs. The remaining 64 runs took between 8.5 seconds and 3:05 hours (median:35 minutes). The suboptimality is between 0 and 3.6 Gy (median: 0.2 Gy), and is lessthan 0.1 Gy in 13 out of 64 cases. Alg. 4 finds a feasible point in all cases. For the 64cases that ART3+O could solve, Alg. 4 finishes before the time limit in 7 cases, butis always slower than ART3+O. The suboptimality ranges between 0 Gy and 12.5 Gy(median: 1.1 Gy). Alg. 4 outperforms ART3+O by at least 0.1 Gy in 14 cases, whilethe opposite is true in 41 cases. Alg. 4 solves the 13 cases that ART3+O incorrectlydeclares infeasible, with the suboptimality ranging between 0 and 5.7 Gy (median 0.2Gy). The results presented by Chen et al. (2010) show that ART3+O can be used to quicklyoptimize a treatment plan. Its main weakness is the arbitrary iteration limit. In ournumerical examples, we found out that this weakness is not purely theoretical, butthat it has practical impact. The iteration limit suggested by Chen et al. (2010) is toolow, and even using a considerably higher number still leads to suboptimalities of upto 3.6 Gy. The theoretical advantage of the proposed alternative does not translateinto a practical advantage. Despite the significantly higher cpu time (24 hours versus3:05 hours), the final objective value is often worse.The dual problem of (1) is min { b T y : A T y ≥ c, y ≥ } . Alternative solvers for linearoptimization, notably the simplex method and an interior point method, generate notjust an optimal solution x ∗ to (1), but also an optimal dual solution y ∗ as a byproduct.Due to strong duality, c T x ∗ = b T y ∗ if the problem is feasible and bounded. For M ≤ c T x ∗ , x ∗ is a solution for the first instance of ART3+, while for M > c T x ∗ , y = 1 / ( M − b T y ∗ ) and y = y y ∗ is a solution for the second instance of ART3+ inAlg. 4. That means that alternative solvers solve the same two problems as Alg. 4.It is an open question as to why ART3+ often finds near optimal solutions in amoderate number of iterations, while it takes significantly more iterations to prove thata certain objective value is unattainable. A possible explanation is that although theprimal can have a large number of constraints, many are correlated or redundant, suchas constraints on the maximum dose in two neighboring voxels. Correlated constraintsdo probably not occur in the second instance of ART3+ in Alg. 4, because that wouldrequire two different x j to deliver almost the same dose. The phenomenon that it iseasy to find an x and hard to find a y is problem specific. Outside of the domainof radiation therapy there are problems for which the opposite is true. An artificialexample can be constructed by taking the dual of (1) for a radiation therapy problem.Our improvement to ART3+O can be applied to ART3+ as well. ART3+ wasoriginally proposed for image reconstruction by Herman (1975) where it approximately olves a large system Ax = b by finding a solution to b − ε ≤ Ax ≤ b + ε. (5)Herman (1975) suggested to select ε based on prior knowledge or experimentally, butthere is a delicate balance. When ε is chosen too small, there may be no solution,while a large value allows for solutions that do not accurately reconstruct the originaldistribution. It is therefore desirable to determine the smallest ε for which a solutionexists. By Farkas’ lemma, there is either an x that satisfies (5), or a y that satisfies A T y − A T y ≥
0, ( b + ε ) T y − ( b − ε ) T y ≤ − y ≥ y ≥
0. One could thereforerun ART3+ in parallel on both sets of inequalities to establish if there is a solution to(5) for a given value of ε .It is well known that the order in which the projections are performed has a signif-icant impact on the performance of a projection algorithm. Ideally, consecutive con-straints are as orthogonal as possible (Guan and Gordon, 1994). To test the impacton our results, we have reordered the constraints in a way that consecutive constraintsare as orthogonal as possible. We did this greedily, starting with the nonnegativityconstraints (which are fully orthogonal), and then adding the constraints one by onesuch that the added constraint was as orthogonal as possible to the previous one. Dueto the computational cost of this ordering heuristic, we only tested it on the small case.To our surprise, the number of iterations in creased for the second instance of ART3+,making the overall performance worse. An alternative ordering based on the goldenratio (Kohler, 2004) did not improve the performance either.The lack of a good stopping criterion has also been noted for projection methodsin convex optimization (Gibali et al., 2014). The conic variant of Farkas’ lemma couldbe used for such problems. Let f i be closed and convex functions and let f ∗ i denotetheir convex conjugate functions. Exactly one of the following sets is not empty:1 . { x ∈ R n : f i ( x ) < i = 1 , . . . , m } (6)2 . { λ ∈ R m + , v i ∈ R n : m X i =1 v i = 0 , X i λ i f ∗ (cid:18) v i λ i (cid:19) ≤ − } . (7)The derivation is given in Appendix A. A necessary step toward full dimensionality issubstituting out one of the vectors v i . By using the two sets (6) and (7), the methodsin this paper extend to convex optimization. ART3+O cannot guarantee ε -optimality of the final solution due to an arbitrary stop-ping criterion. An improvement that provides such a guarantee was suggested, butdespite its theoretical basis, its runtime is too high, especially in comparison with al-ternative solvers. If an optimality guarantee is not required, ART3+O remains a viablealgorithm, but their users should be aware that a solution may be suboptimal. Extension to convex optimization
Theorem 2 (Conic variant of Farkas’ lemma) Let K be a closed convex cone withdual cone K ∗ and let F ∈ R m × n . Exactly one of the following statements is true(Craven and Koliha, 1977):1. There exists an x ∈ R n such that F x = b and x ∈ K .2. There exists a y ∈ R m such that F T y ∈ K ∗ and b T y ≤ − . The closure of the set (6) can be denoted as { ( x, t ) : t = 1 , ( x, t ) ∈ K = ∩ mi =1 K i } ,with K i = { ( x, t ) : tf i ( x/t ) ≤ , t ≥ } and 0 f i ( x/
0) = lim t ↓ tf i ( x/t ). Note that K is indeed a closed convex cone. The x in Farkas’ lemma is actually the vector( x, t ) ∈ R n +1 , F = (0 , . . . , , ∈ R × ( n +1) , and the second set in Farkas’ lemma is { y ∈ R : (0 , . . . , , y ) ∈ K ∗ , y ≤ − } . The condition (0 , . . . , , y ) ∈ K ∗ is by definitionequivalent to ty ≥ ∀ ( x, t ) ∈ ∩ mi =1 K i . To eliminate the “ ∀ ” quantifier, we rewrite thisdefinition as min x,t ≥ { ty : tf i ( x/t ) ≤ i = 1 , . . . , m } ≥
0, and then apply Lagrangeduality to the minimization problem using the method by Roos et al. (2018). Theobjective function of this problem is g ( x, t ) = ty and the constraint functions are g i ( x, t ) = tf i ( x/t ). The corresponding conjugate functions are g ∗ ( v , u ) = 0 if v = 0and u = y ( ∞ otherwise) and g ∗ i ( v i , u i ) = 0 if f ∗ ( v i )+ u i ≤ ∞ otherwise). Therefore,by Roos et al. (2018), the Lagrange dual is max λ ≥ ,v,u { P mi =1 v i = 0 , y + P mi =1 u i =0 , λ i f ∗ ( v i /λ i ) + u i ≤ } . The second statement in Farkas lemma reduces to: thereexists a y ∈ R , λ ∈ R m + , v i ∈ R n and u i ∈ R such that P mi =1 v i = 0, y + P mi =1 u i = 0, λ i f ∗ ( v i /λ i ) + u i ≤ y ≤ −
1. Substituting out y and P i u i , we get the set (7). Acknowledgment
Supported in part by NIH U19 Grant 5U19CA021239-38.
References
G. Bahr, J. Kereiakes, H. Horwitz, R. Finney, J. Galvin, and K. Goode. The methodof linear programming applied to radiation treatment planning.
Radiology , 91(4):686–693, 1968.H. H. Bauschke and J. M. Borwein. On projection algorithms for solving convexfeasibility problems.
SIAM review , 38(3):367–426, 1996.R. G. Bland, D. Goldfarb, and M. J. Todd. The ellipsoid method: A survey.
OperationsResearch , 29(6):1039–1091, 1981.T. Bortfeld, J. B¨urkelbach, R. Boesecke, and W. Schlegel. Methods of image reconstruc-tion from projections applied to conformation radiotherapy.
Physics in Medicine &Biology , 35(10):1423, 1990. . Breedveld, B. van den Berg, and B. Heijmen. An interior-point implementationdeveloped and tuned for radiation therapy treatment planning. Computational Op-timization and Applications , 68(2):209–242, 2017.Y. Censor and E. Tom. Convergence of string-averaging projection schemes for in-consistent convex feasibility problems.
Optimization Methods and Software , 18(5):543–554, 2003.W. Chen, D. Craft, T. M. Madden, K. Zhang, H. M. Kooy, and G. T. Herman. Afast optimization algorithm for multicriteria intensity modulated proton therapyplanning.
Medical physics , 37(9):4938–4945, 2010.W. Chen, J. Unkelbach, A. Trofimov, T. Madden, H. Kooy, T. Bortfeld, and D. Craft.Including robustness in multi-criteria optimization for intensity-modulated protontherapy.
Physics in Medicine & Biology , 57(3):591–608, 2012.B. D. Craven and J. J. Koliha. Generalizations of Farkas’ theorem.
SIAM Journal onMathematical Analysis , 8(6):983–997, 1977.A. R. De Pierro and A. N. Iusem. A simultaneous projections method for linearinequalities.
Linear Algebra and its applications , 64:243–253, 1985.J. Farkas. ¨Uber die Theorie der Einfachen Ungleichungen.
Journal f¨ur die Reine undAngewandte Mathematik , 124:1–27, 1902.A. Gibali, K.-H. K¨ufer, D. Reem, and P. S¨uss. A generalized projection-based schemefor solving convex constrained optimization problems.
Computational Optimizationand Applications , 1–26, 2014.H. Guan and R. Gordon. A projection access order for speedy convergence of ART(algebraic reconstruction technique): a multilevel scheme for computed tomography.
Physics in Medicine & Biology , 39(11):2005, 1994.G. T. Herman. A relaxation method for reconstructing objects from noisy x-rays.
Mathematical Programming , 8(1):1–19, 1975.G. T. Herman and W. Chen. A fast algorithm for solving a linear feasibility problemwith application to intensity-modulated radiation therapy.
Linear algebra and itsapplications , 428(5-6):1207–1217, 2008.T. Kohler. A projection access scheme for iterative reconstruction based on the goldensection. In
Nuclear Science Symposium Conference Record, 2004 IEEE , volume 6,3961–3965, 2004.E. Lessard and J. Pouliot. Inverse planning anatomy-based dose optimization forHDR-brachytherapy of the prostate using fast simulated annealing algorithm anddedicated objective function.
Medical Physics , 28(5):773–779, 2001.A. Lomax. Intensity modulation methods for proton radiotherapy.
Physics in Medicine& Biology , 44(1):185, 1999. . E. Romeijn, R. K. Ahuja, J. F. Dempsey, A. Kumar, and J. G. Li. A novel linearprogramming approach to fluence map optimization for intensity modulated radia-tion therapy treatment planning. Physics in Medicine & Biology , 48(21):3521–3542,2003.C. Roos, M. Balvert, B. L. Gorissen, and D. den Hertog. A universal and structuredway to derive dual optimization problem formulations.
Optimization Online , 2018.D. M. Shepard, M. C. Ferris, G. H. Olivera, and T. R. Mackie. Optimizing the deliveryof radiation therapy to cancer patients.
Siam Review , 41(4):721–744, 1999.A. Trofimov, J. Unkelbach, and D. Craft. Treatment-planning optimization. In H. Pa-ganetti, editor,
Proton therapy physics , 461–484. CRC Press/Taylor & Francis, BocaRaton, FL, USA, 2012., 461–484. CRC Press/Taylor & Francis, BocaRaton, FL, USA, 2012.