A nonconvex optimization approach to IMRT planning with dose-volume constraints

TL;DR

This paper introduces a novel relaxation method for nonconvex dose-volume constraints in IMRT planning, enabling more accurate optimization in radiation therapy with proven convergence and adaptability to complex clinical scenarios.

Contribution

It presents a new relaxation approach and efficient algorithms for nonconvex dose-volume constraints in IMRT, improving upon previous convex approximation methods.

Findings

Demonstrated effectiveness on the CORT dataset

Algorithms are provably convergent

Adaptable to multiple tumors and organs-at-risk

Abstract

Fluence map optimization for intensity-modulated radiation therapy planning can be formulated as a large-scale inverse problem with competing objectives and constraints associated with the tumors and organs-at-risk. Unfortunately, clinically relevant dose-volume constraints are nonconvex, so standard algorithms for convex problems cannot be directly applied. While prior work focused on convex approximations for these constraints, we propose a novel relaxation approach to handle nonconvex dose-volume constraints. We develop efficient, provably convergent algorithms based on partial minimization, and show how to adapt them to handle maximum-dose constraints and infeasible problems. We demonstrate our approach using the CORT dataset, and show that it is easily adaptable to radiation treatment planning with dose-volume constraints for multiple tumors and organs-at-risk.

Tables

Table 1. Table 1: Results for Section  4.1 .

	PTV	OAR
	D95 (Gy)	% $>$ 30 Gy	Time (s)
Initialization	79.65	64.14	0.18
Solution w/ Alg. 2	79.17	34.16	7.22
Solution w/ Alg. 3 (a)	79.03	29.61	33.61
Solution w/ Alg. 3 (b)	78.05	19.84	163.51

Table 2. Table 2: Results for Section  4.2 .

	PTV	PTV	OAR	OAR	OAR
	D95 (Gy)	% 81 Gy	% $>$ 20 Gy	% $>$ 40 Gy	% $>$ 60 Gy	Time (s)
Initialization	79.65	53.26	68.81	61.29	38.47	0.18
Solution w/ Alg. 2	80.26	24.27	61.95	42.90	21.84	8.89
Solution w/ Alg. 3	81.32	4.09	59.10	39.26	19.78	55.86

Table 3. Table 3: Results for Section  4.3 .

	PTV 1	PTV 2	OAR 1	OAR 2
	D95 (Gy)	D95 (Gy)	% $>$ 50 Gy	% $>$ 30 Gy	Time (s)
Initialization	76.96	58.53	82.47	91.87	0.37
Solution w/ Alg. 2	76.70	57.51	57.33	38.14	15.52
Solution w/ Alg. 3	76.48	57.31	47.53	29.93	84.83

Table 4. Table 4: Results for Section  4.4.1 with dose–volume objectives.

	PTV	PTV	OAR
	D95 (Gy)	% 81 Gy	% $>$ 30 Gy	Iter. ($k$)	Time (s)
Solution w/ $\left({𝒫}_2\right)$	79.65	55.32	38.35	5	1.83
Solution w/ (37)	78.74	44.33	56.55	173	5.20
Solution w/ (38)	77.89	43.74	20.63	5	198.70

Table 5. Table 5: Results for Section  4.4.1 with dose–volume constraints.

	PTV	PTV	OAR
	D95 (Gy)	% 81 Gy	% $>$ 30 Gy	Obj.	Time (s)
Solution w/ ${\left({𝒫}_2\right)}^a$	81.32	2.97	29.43	7.80	46.94
Solution w/ ${\left({𝒫}_2\right)}^b$	81.00	4.37	29.43	6.95	65.56
Solution w/ (27)	81.00	3.26	28.40	9.12	63.60
Solution w/ (39)	81.00	2.20	28.16	8.94	379.46

Table 6. Table 6: Results for Section  4.4.2 with dose–volume objectives.

	PTV 1	PTV 1	PTV 2	PTV 2	OAR 1	OAR 2
	D95 (Gy)	% 81 Gy	D95 (Gy)	% 60 Gy	% $>$ 50 Gy	% $>$ 30 Gy	Iter. ($k$)	Time (s)
Solution w/ $\left({𝒫}_2\right)$	76.95	57.86	57.51	52.95	64.58	49.47	3	2.80
Solution w/ (37)	76.57	42.07	56.87	41.18	77.48	44.45	195	12.98
Solution w/ (38)	74.18	45.73	57.34	42.25	25.68	27.55	3	2014.94

Table 7. Table 7: Results for Section  4.4.2 with dose–volume constraints.

	PTV 1	PTV 1	PTV 2	PTV 2	OAR 1	OAR 2
	D95 (Gy)	% 81 Gy	D95 (Gy)	% 60 Gy	% $>$ 50 Gy	% $>$ 30 Gy	Obj.	Time (s)
Solution w/ ${\left({𝒫}_2\right)}^a$	81.17	4.71	60.32	3.90	48.39	29.84	27.71	138.36
Solution w/ ${\left({𝒫}_2\right)}^b$	81.00	4.46	60.01	3.68	44.31	28.78	28.64	309.83
Solution w/ (39)	81.00	3.07	60.01	2.69	33.41	25.02	36.02	1118.13

Table 8. Table 8: PTV results for Section  4.4.3 with dose–volume objectives.

	PTV 1	PTV 1	PTV 2	PTV 2
	D95 (Gy)	% 81 Gy	D95 (Gy)	% 60 Gy	Iter. ($k$)
Solution w/ $\left({𝒫}_2\right)$	76.83	58.18	58.25	53.77	2
Solution w/ (37)	76.86	43.68	58.13	46.09	97
Solution w/ (38)	74.17	45.21	57.78	42.80	2

Table 9. Table 9: OAR results for Section  4.4.3 with dose–volume objectives.

	OAR 1	OAR 2	OAR 3	OAR 4
	% $>$ 50 Gy	% $>$ 40 Gy	% $>$ 30 Gy	% $>$ 30 Gy	Time (s)
Solution w/ $\left({𝒫}_2\right)$	65.61	52.89	2.37	2.38	2.16
Solution w/ (37)	73.89	53.02	3.88	3.63	7.38
Solution w/ (38)	21.06	19.37	1.65	1.79	29246.23

Table 10. Table 10: Results for Section  4.4.3 with dose–volume constraints.

	PTV 1	PTV 1	PTV 2	PTV 2	OAR 1	OAR 2
	D95 (Gy)	% 81 Gy	D95 (Gy)	% 60 Gy	% $>$ 50 Gy	% $>$ 40 Gy	Obj.	Time (s)
Solution w/ ${\left({𝒫}_2\right)}^a$	81.05	4.93	60.18	3.97	46.68	38.12	23.67	215.38
Solution w/ ${\left({𝒫}_2\right)}^b$	81.00	4.89	60.01	4.10	48.27	37.58	15.25	417.34
Solution w/ (27)	81.01	4.17	60.26	1.19	37.49	32.75	52.90	474.25

Equations

$$\frac{1}{2}\Vert A_{1}x-81{\Vert }_2^2+\frac{\lambda }{2}\Vert x{\Vert }_2^2$$

$$\Vert (A_{2}x-20{)}_{+}{\Vert }_0\le 1.$$

$$\Vert (A_{2}x-20{)}_{+}{\Vert }_0=\Vert \max \{0,A_{2}x-20\}{\Vert }_0=\text{number of voxels s.t. }{A}_{2}x>20,$$

$$\displaystyle\mathmakebox[\widthof{$\underset{\displaystyle x\in\mathcal{C}}{\mathrm{subject~{}to}}$}]{\underset{\displaystyle x\in\mathcal{C}}{\mathrm{minimize}}}\quad f(x)$$

$$\displaystyle\mathmakebox[\widthof{$\underset{\displaystyle\phantom{x\in\mathcal{C}}}{\mathrm{\,subject~{}to}}$}][c]{\mathmakebox[\widthof{$\mathrm{minimize}$}][l]{\mathrm{subject~{}to}}}\quad$$

$$({\mathcal{P}}_0)$$

$$\displaystyle\mathmakebox[\widthof{$\underset{\displaystyle\phantom{x\geq 0}}{\mathrm{\,subject~{}to}}$}][c]{\mathmakebox[\widthof{$\mathrm{minimize}$}][l]{\mathrm{subject~{}to}}}\quad$$

$$1^T{A}_{j}x\le n_j{d}_j^{\text{mean}},$$

$$\Vert (A_{j}x-d_j^{\text{dv}}{)}_{+}{\Vert }_0\le \frac{n_j{p}_j}{100},$$

$$({\mathcal{P}}_1)$$

$$\displaystyle\mathmakebox[\widthof{$\underset{\displaystyle\phantom{x\geq 0}}{\mathrm{\,subject~{}to}}$}][c]{\mathmakebox[\widthof{$\mathrm{minimize}$}][l]{\mathrm{subject~{}to}}}\quad$$

$$({\mathcal{P}}_2)$$

$$\text{mathmakebox}+\sum _{j\in {\mathcal{J}}^{\text{dv}}}\frac{{\alpha }_j}{2n_j}\Vert w_j-(A_{j}x-d_j^{\text{dv}}){\Vert }_2^2+\frac{\lambda }{2}\Vert x{\Vert }_2^2$$

$$\displaystyle\mathmakebox[\widthof{$\underset{\displaystyle\phantom{x\geq 0,w}}{\mathrm{\,subject~{}to}}$}][c]{\mathmakebox[\widthof{$\mathrm{minimize}$}][l]{\mathrm{subject~{}to}}}\quad$$

$$\text{epi}(f)=\{(x,\beta ):f(x)\le \beta \}.$$

$$f^{\alpha }{\to }_{\text{epi}}f\text{exactly when}\text{epi}(f^{\alpha })\to \text{epi}(f).$$

$$f^{{\alpha }_j}\le f^{{\alpha }_{j+1}}\text{when}{\alpha }_j\le {\alpha }_{j+1},$$

$$({\mathcal{P}}_3)$$

$$\text{mathmakebox}+\sum _{j\in {\mathcal{J}}^{\text{dv}}}\frac{{\alpha }_j}{2n_j}\Vert w_j-(A_{j}x-d_j^{\text{dv}}){\Vert }_2^2+\frac{\lambda }{2}\Vert x{\Vert }_2^2$$

$$\displaystyle\mathmakebox[\widthof{$\underset{\displaystyle\phantom{w}}{\mathrm{\,subject~{}to}}$}][c]{\mathmakebox[\widthof{$\mathrm{minimize}$}][l]{\mathrm{subject~{}to}}}\quad$$

$$w^{(k+1)}\in {\text{proj}}_{\Omega }\left[w^{(k)}-T\nabla g\left(w^{(k)}\right)\right],$$

$$g(w)=\underset{x\ge 0}{\min }\frac{{\alpha }_1}{2n_1}\Vert A_{1}x-d_1{\Vert }_2^2+\frac{{\alpha }_2}{2n_2}\Vert w-(A_{2}x-d_2^{\text{dv}}){\Vert }_2^2+\frac{\lambda }{2}\Vert x{\Vert }_2^2.$$

$${\delta }_{+}(x)=\{\begin{array}{ll}0& \text{if }x\ge 0,\\ +\infty & \text{otherwise,}\end{array}$$

$$g(w)$$

$$x(w)$$

$$\nabla g(w)=\frac{{\alpha }_2}{n_2}\left\{w-\left[A_{2}x(w)-d_2^{\text{dv}}\right]\right\}\text{and}\text{lip}(\nabla g)\le \frac{{\alpha }_2}{n_2}.$$

$${\Vert \nabla g(w)-\nabla g\left(w^{\prime }\right)\Vert }_2^2$$

$$\le \sum _{j\in {\mathcal{J}}^{\text{dv}}}{\left(\frac{{\alpha }_j}{n_j}\right)}^2{\Vert w_j-w_j^{\prime }\Vert }_2^2,$$

$$\le \underset{j\in {\mathcal{J}}^{\text{dv}}}{\max }{\left(\frac{{\alpha }_j}{n_j}\right)}^2\sum _{j\in {\mathcal{J}}^{\text{dv}}}{\Vert w_j-w_j^{\prime }\Vert }_2^2,$$

$$=\underset{j\in {\mathcal{J}}^{\text{dv}}}{\max }{\left(\frac{{\alpha }_j}{n_j}\right)}^2{\Vert w-w^{\prime }\Vert }_2^2,$$

$$0\in \nabla g(w)+N_{\Omega }(w),$$

$$\sum _{j\in {\mathcal{J}}^{\text{dv}}}\frac{1}{t_j^2}{\Vert w_j-{\bar{w}}_j\Vert }_2^2\le ϵ,$$

$$\frac{1}{t_j}\left(w_j-{\bar{w}}_j\right)\in \frac{\partial g}{\partial w_j}(w)+N_{{\Omega }_j}\left({\bar{w}}_j\right),$$

$$f(x,w^{(k)})=\frac{1}{2}\Vert A_{1}x-81{\Vert }_2^2+\frac{5}{2}\Vert w^{(k)}-(A_{2}x-20){\Vert }_2^2+\frac{\lambda }{2}\Vert x{\Vert }_2^2$$

$$g(w)=\underset{x\ge 0}{\min }\frac{1}{2}\Vert A_{1}x-81{\Vert }_2^2+\frac{5}{2}\Vert w-(A_{2}x-20){\Vert }_2^2+\frac{\lambda }{2}\Vert x{\Vert }_2^2$$

$$x^{(0)}={\mathrm{argmin}}_{x\ge 0}\sum _{i\in \mathcal{I}}\frac{{\alpha }_i}{2n_i}\Vert A_{i}x-d_i{\Vert }_2^2+\frac{\lambda }{2}\Vert x{\Vert }_2^2.$$

$${\mathrm{argmin}}_{w_j}f(x,w)=A_{j}x-d_j^{\text{dv}},$$

Full text

1 Introduction

Cancer is one of the leading causes of death worldwide, with approximately 38.4% of the population expected to develop cancer at some point in their lives [1]. Many cancer therapies target the weakened DNA damage response of cancer cells; in particular, radiotherapy uses high-energy ionizing radiation to damage cancer cell DNA [2, 3]. Unfortunately, this radiation also damages the healthy cells in its path, so the main goal of radiation treatment planning is to maximize the differential between damage to the tumor and to nearby healthy tissue. The most common form of radiotherapy is external beam radiation therapy (EBRT), where photon beams generated by a linear accelerator are directed at a patient from a rotating gantry [4]. By surrounding a tumor with beams from multiple angles, physicians can design treatment plans that target tumors with a variety of geometries while avoiding a large dose to critical structures.

One way to limit damage to healthy tissue is to require that treatment plans satisfy various constraints on the amount of radiation delivered to organs-at-risk (OARs). For example, it is often clinically relevant to impose maximum-dose constraints on a fraction of an organ, known as a dose–volume constraint. Unfortunately, these constraints are nonconvex, so convex formulations and algorithms cannot be directly applied to the treatment planning problem. We propose a novel approach to handle dose–volume constraints for intensity-modulated radiation therapy (IMRT) while preserving their nonconvex structure, as opposed to previous efforts that focused on convex approximations. The proposed method is amenable to efficient algorithms based on partial minimization and naturally adapts to handle maximum-dose constraints and cases of infeasiblity. We demonstrate our approach using the CORT dataset [5, 6], and we show that it easily adapts to radiation treatment planning with dose–volume constraints for multiple tumors and OARs. To start off, we give a brief overview of EBRT, including IMRT, and the specific challenge posed by dose–volume constraints.

1.1 External beam radiation therapy

Conventional EBRT began soon after the discovery of X-rays in 1895 and was significantly impacted by technological innovations such as advanced imaging and multileaf collimators (MLC) [7]. Specifically, the invention of computed tomography (CT) in the 1970s enabled accurate three-dimensional mapping of a patient’s internal anatomy, while the introduction of the MLC to radiotherapy in the 1990s provided more control over beam shape and fluence. For example, three-dimensional conformal radiation therapy (3D-CRT) uses the sliding leaves of the MLC to shape beam outlines to match tumor geometries, reducing radiation exposure to healthy tissue, while IMRT uses MLCs to not only shape beam outlines, but also modulate beam fluence, giving additional control over the patient dose distribution [8].

In the early history of radiotherapy, treatment plans were created using a trial-and-error approach called forward planning, where physicians would calculate dose distributions for different sets of beam parameters until an acceptable treatment plan was found. However, modulating beam fluence for IMRT involves discretizing each beam into subunits, called beamlets or bixels, which dramatically increases the number of variables under consideration, rendering forward planning infeasible [9, 10]. In the 1980s, the idea of inverse planning for IMRT was introduced by Brahme [11, 12], where mathematical optimization was used to calculate beam parameters to deliver a specified dose distribution.

Within radiotherapy inverse planning, there are a variety of different paradigms regarding dose objectives, mathematical models, and optimization methods. For instance, while our proposed model and many others seek to optimize dose distributions defined in terms of physical criteria such as dose and volume [9, 13], it is also possible to consider biological criteria such as tumor control probability (TCP) and normal-tissue complication probability (NTCP) [14], equivalent uniform dose (EUD) [15, 16, 17], and biologically effective dose (BED) [18, 19]. Additionally, many early approaches to inverse planning used analytical direct methods inspired by the similarities between IMRT planning and CT image reconstruction [20, 21, 22], but the trend has since moved towards approximate iterative methods [23, 24, 25].

Amongst iterative methods, the presence of nonconvexity within many IMRT models has fueled the debate between stochastic methods such as simulated annealing [26, 27, 28] and genetic algorithms [29, 30], which can escape local minima but tend to have higher computational costs, versus gradient-based methods using techniques such as Newton’s method [21, 31], L-BFGS [32, 33], and conjugate gradient methods [34, 35], which tend to be faster but can get trapped in local minima. Our method uses a projected gradient approach to handle issues related to nonconvexity, which we discuss in Section 1.3. For a more detailed history of radiotherapy and IMRT, see [7, 8, 9, 10, 13, 36].

1.2 Intensity-modulated radiation therapy

In current practice, the patient undergoes imaging (e.g., CT) before IMRT treatment, producing a three-dimensional model of their internal anatomy. This model is discretized into volume units (voxels) and labeled by structure. Next, the physician specifies the prescription, often in terms of a uniform target to the tumor with various constraints on critical structures. The goal of the treatment planner is then to determine the number of beams, beam angles, and beamlet intensities that will deliver the prescribed dose of radiation to the tumor while keeping doses to healthy tissue low. Inverse planning for IMRT typically involves three steps [10]:

Beam angles: Determine the number of beams and their orientations. 2. 2.

Fluence map: Calculate beamlet intensities to deliver the desired dose distribution to the patient. 3. 3.

Segmentation: Design MLC sequences that achieve the optimized fluence map.

These steps require solving optimization problems that can be approached separately or in combination [37, 38, 39, 40]. In this paper, we focus on the fluence map optimization (FMO) problem of assigning beamlet intensities for a given set of beam angles.

The FMO problem has been approached using a number of different mathematical frameworks including linear programming [41, 42, 43], piecewise linear–quadratic models [44, 45, 46], quadratic programming [31, 34, 47], and other nonlinear programming models [14, 15, 25, 48]. Using multi-objective optimization, many studies have generated solutions along the Pareto frontier to explore the trade-off between competing tumor and healthy tissue objectives and constraints [17, 32, 49]. Other approaches include formulating the FMO problem as a dynamical system [50] or eigenvalue problem [51], using a database of reference dose–volume histograms to guide optimization [52], and clustering beams and voxels to handle large-scale problems [53]. We develop a new approach that is well adapted to handle nonconvex dose–volume constraints, in addition to convex maximum-dose and mean-dose constraints. We first explain all three constraint types.

Radiation dose is measured in gray (Gy), where one Gy is defined as the absorption of one joule of energy per kilogram of matter. Although lower doses are always better for healthy tissue, tolerable doses for many tissues are known empirically and depend on the risk level that a patient and physician are willing to accept. In practice, tolerable doses for different tissues are formulated using the following three types of constraints:

• •

Maximum dose: No voxel in the organ receives more than $d^{\text{max}}$ Gy.

• •

Mean dose: Average radiation per voxel does not exceed $d^{\text{mean}}$ Gy.

• •

Dose–volume: At most $p\%$ of the organ volume receives more than $d^{\text{dv}}$ Gy.

Serial structures such as the spinal cord and brainstem lose functionality if any of their subvolumes are damaged, so maximum-dose constraints are a good indicator of tissue damage. Maximum-dose and mean-dose constraints define a convex feasible region for the decision variables. The associated optimization problems can be solved at scale; for example, inverse planning with linear inequality constraints can be formulated as a convex quadratic program, with state-of-the-art approaches and commercial implementations readily available (see e.g., [54]). On the other hand, parallel structures such as the liver and lungs are composed of semi-independent units. In this case, a fraction of the parallel organ can be damaged without loss of functionality, so dose–volume constraints are appropriate [34]. Assigning maximum-dose constraints on partial volumes to achieve the desired dose–volume criteria was first introduced by Langer and Leong in 1987 [41] and brought to the attention of the mathematical community by Shepard et al. in 1999 [36].

1.3 Challenge of dose–volume constraints

In contrast to maximum-dose and mean-dose constraints, dose–volume constraints are combinatorial, since the choice of the $p\%$ of organ voxels that may receive more than $d^{\text{dv}}$ Gy is left to the planner. In early IMRT research, using only a small number of beamlets and voxels, it was possible to solve the treatment planning problem for all possible voxel combinations satisfying a dose–volume constraint and then choose the best resulting plan [41]. This approach becomes infeasible as the problem size increases to clinically relevant situations: for an organ with $n$ voxels and a dose–volume constraint imposed on $m$ of these voxels, there are $\nicefrac{{n!}}{{(n-m)!m!}}$ possible voxel combinations that would satisfy the constraint [34].

To illustrate the nonconvex geometry of a feasible beamlet region induced by a dose–volume constraint, we recreate an example from Wu et al. [55] using the CORT dataset described in Section 4. We consider three body voxels of a patient, where voxel $u$ belongs to a prostate tumor and voxels $v_{1}$ and $v_{2}$ belong to the rectum (left panel of Figure 1). Using two beamlets $x_{1}$ and $x_{2}$, we aim to deliver a uniform dose of 81 Gy to the tumor while ensuring that no more than 50% of the rectum volume (one voxel in this case) exceeds 20 Gy.

We calculate the optimal beamlet intensities by minimizing the objective function

[TABLE]

subject to the dose–volume constraint

[TABLE]

The matrices $A_{1}\in\mathbb{R}^{1\times 2}$ and $A_{2}\in\mathbb{R}^{2\times 2}$ map beamlet intensities to prostate and rectum voxel doses, the vector $x\in\mathbb{R}^{2}$ contains the beamlet intensities, and the regularization term controlled by $\lambda=5\times 10^{-6}$ is added to stabilize beamlet solutions, amortizing the ill-conditioned beamlet-to-voxel map. The dose–volume constraint (2) forces the number of rectum voxels receiving a dose greater than 20 Gy to be less than or equal to 1:

[TABLE]

where the $\ell_{0}$-norm returns the number of nonzero elements in a vector (e.g., [56]). In Figure 1 (right), we plot the feasible beamlet region along with the contours of the objective function. Since negative beamlet intensities are impossible, we restrict the beamlet values to the first quadrant.

There are two possibilities for satisfying the dose–volume constraint: either the first rectum voxel is allowed to exceed 20 Gy (Figure 1, $v_{2}\leq 20$ region), or the second rectum voxel is allowed to exceed 20 Gy (Figure 1, $v_{1}\leq 20$ region). The union of these two regions forms an L-shaped nonconvex region with two local minima at points A and B. Depending on the beamlet initialization and algorithm used to minimize the objective function, it is possible to converge to either of these two local minima, and there is no guarantee that the global minimum at B will be reached. In general, the feasible beamlet region imposed by a dose–volume constraint comprises the nonconvex union of regions satisfying maximum-dose constraints for different choices of voxels [57].

The FMO problem with nonconvex dose–volume constraints is related to the general convex-cardinality problem

[TABLE]

where $f(x)$ is a convex function and $\mathcal{C}$ is a convex set. The $\ell_{0}$-norm constraint defines a highly nonconvex set: the union of all subspaces of dimension at most $k$. Problems with these cardinality (or sparsity) constraints have become widely used within the past three decades [58], with applications such as portfolio optimization [59], signal processing [60], and compressive sensing [61]. Due to the combinatorial nature of these constraints, the related feasibility problem is NP-complete [62], and the problem itself is known to be NP-hard in general [62, 63], except for in special cases [64]. More generally, it is an NP-hard problem just to certify whether a point of a nonlinear function is a local minimum (let alone a global one) [65]. For a more detailed overview of cardinality problems, see [58, 66, 67, 68].

Many of the standard approaches for solving the convex-cardinality problem have been applied to the FMO problem with nonconvex dose–volume constraints, most relying upon convex relaxations or approximations. One technique is to reformulate both problems as a mixed-integer programming (MIP) problem, which can be solved using branch-and-bound or branch-and-cut methods [37, 62, 69], simplex methods using big $M$ variables [70], or by replacing the binary constraint $x\in\{0,1\}^{n}$ with the linear constraint $x\in[0,1]^{n}$ [49, 71, 72]. Another popular technique for solving the convex-cardinality problem is to relax the nonconvex $\ell_{0}$-norm constraint (4b) using the convex $\ell_{1}$-norm as either a constraint or regularizer [58, 73, 74]. Convex approximations to the nonconvex dose–volume constraint $\|(Ax-d^{\text{dv}})_{+}\|_{0}\leq k$ include the convex constraint $\|(1+\lambda(Ax-d^{\text{dv}}))_{+}\|_{1}\leq k$ [46] and conditional value-at-risk (CVaR) constraints [75] that limit the average dose to the tail of the dose distribution [38, 43, 45, 76]. One common dose–volume penalty includes the convex objective term $\lambda\|(Ax-d^{\text{dv}})_{+}\|_{2}^{2}$ [10, 21, 31, 34], though the $\ell_{2}$-norm is less effective at promoting sparsity than the $\ell_{1}$-norm [74]. When replacing constraints with objectives, weight factors can be applied to individual vector elements, and iterative re-weighting schemes can then be used to improve convergence [47, 49, 61, 77].

The complication of nonconvexity arises from the combinatorial nature of the cardinality constraint (4b); once a sparsity pattern for $x$ has been chosen, the problem becomes convex. This observation is the motivation behind the multi-stage approach of 1) solving an approximation of the original problem, 2) determining a sparsity pattern for $x$ based on the resulting solution, and 3) polishing the solution from Step 1 by solving a convex problem using the sparsity pattern determined in Step 2 [68, 73, 78, 79]. In the FMO problem with nonconvex dose–volume constraints, a variety of heuristics have been used to choose the $p\%$ of organ voxels that may receive more than $d^{\text{dv}}$ Gy in Step 3. For example, Morrill et al. [80, 81] divide organ volumes a priori into high-dose and low-dose regions based on tumor proximity, Saberian et al. [18, 19] and Hou et al. [50] solve the FMO problem without dose–volume constraints to determine sparsity patterns, while Fu et al. [46] and Mukherjee et al. [72] use convex relaxations of the original problem to determine which voxels may receive more than $d$ Gy.

The great watershed in optimization is between convex and nonconvex problems [82]. When minimizing a convex objective function over a convex set, all local minima are also global. Conversely, optimizing any function over a nonconvex set leads to a nonconvex problem with potentially many local minima. Even in the small 2-D example in Figure 1, the dose–volume constraint creates a nonconvex feasible set where each choice of voxel corresponds to a different local minimum. As mentioned above, early IMRT treatment planning relied on stochastic methods such as simulated annealing and genetic algorithms to avoid converging to local minima, but these methods suffered from slow convergence rates [8, 34]. Fortunately, many investigations into the properties of local minima induced by nonconvex dose–volume constraints have concluded that, given a good initialization, it is possible to find local minima whose dose distributions are close to the global minimum using fast gradient methods [9, 55, 83], so the use of stochastic methods may not be warranted [13, 14, 84]. In this case, the emphasis in IMRT is often on finding a good local minimum that satisfies clinical objectives rather than finding the global minimum [13, 85].

1.4 Contributions

In this paper, we develop a novel approach to handle dose–volume constraints in the radiation treatment planning problem. In contrast to prior efforts, which focused on convex approximations, our formulation is nonconvex, yet can be efficiently solved using continuous optimization techniques. We start with the idealized problem, formulated using objectives for the tumor and constraints for the OARs, and extend it to a flexible formulation that allows re-weighting strategies to balance infeasible solutions (with respect to the OARs) against delivering a sufficient amount of radiation to the tumor. To solve this formulation, we develop a customized algorithm that combines quadratic programming with fast nonlinear operators, and show that it is guaranteed to converge to a stationary point of the proposed model. The resulting approach is scalable, handles multiple objectives on tumors and OARs, and returns actionable treatment plans on patient-scale datasets.

The paper proceeds as follows. In Section 2 we formulate radiation treatment planning as a nonconvex optimization problem over the fluence map. In Section 3 we develop a custom optimization algorithm and provide a convergence analysis. In Section 4 we present multiple radiation treatment use cases with data from an anonymized cancer patient dataset, demonstrating the feasibility of our approach. Finally, in Section 5 we summarize the results and discuss future avenues for research.

2 Problem Formulation

Given a fixed set of beam angles, the dose delivered to the patient is approximated with a linear mapping from beamlet intensity to voxel dose. The FMO step is to calculate beamlet intensities that achieve the desired dose distribution to planning target volumes (PTVs) indexed by $i\in\mathcal{I}$, while respecting constraints on OARs indexed by $j\in\mathcal{J}$. Let $d_{i}$ denote target delivery doses (often in Gy) to PTV $i\in\mathcal{I}$, and $d_{j}^{\max}$, $d_{j}^{\text{mean}}$, and $d_{j}^{\text{dv}}$ denote the maximum, mean, and dose–volume tolerance values for OAR $j\in\mathcal{J}$, where the size of dose parameters (i.e., scalar or vector) can be determined by context. The sets of OARs that have maximum-dose, mean-dose, and dose–volume constraints are denoted by $\mathcal{J}^{\mbox{\scriptsize{max}}},\mathcal{J}^{\mbox{\scriptsize{mean}}}$, and $\mathcal{J}^{\mbox{\scriptsize{dv}}}$. Let $A_{i}$ and $A_{j}$ be the linear beamlet-to-voxel maps for the PTV $i$ and OAR $j$. The nonnegative variable $x$ encodes the unknown vector of beamlet intensities. The idealized FMO problem with a uniform dose to the PTVs and the three constraints described in Section 1 on the OARs is given by

[TABLE]

where $n_{i}$ and $n_{j}$ indicate the number of voxels in the $i$th PTV and $j$th OAR, and the percentage of voxels in the $j$th OAR that may receive a dose exceeding $d_{j}^{\text{dv}}$ is specified by $p_{j}$. The least-squares term used for the uniform target dose is widely used [10, 21, 31, 34, 47]. In principle, Problem $(\mathcal{P}_{0})$ can be replaced by a feasibility problem with lower and upper dose–volume constraints also placed on the PTV [44, 56, 86]. Here we focus on $(\mathcal{P}_{0})$, but the techniques we propose extend to the feasibility problem.

Problem $(\mathcal{P}_{0})$ is difficult for several reasons. First, it is high-dimensional (beamlets can be on the order of $10^{3}-10^{5}$) [55]. Second, beamlet-to-voxel maps are generally ill-conditioned [8, 25, 30, 51]. Third, dose–volume constraints are combinatorial in nature, as constraint (5d) requires making a choice of voxels allowed to exceed $d_{j}^{\text{dv}}$, so finding the global solution is NP-hard [87]. Finally, while Problem $(\mathcal{P}_{0})$ is always feasible (i.e., $x=0$ is a solution), it may not yield a clinically useful radiation distribution. Feasibility-type reformulations of $(\mathcal{P}_{0})$ often fail to have any meaningful solutions because the target doses almost always compete with OAR doses, as it is common for tumors to be close to or within OARs. Practical radiation treatment planning converts constraints into objectives, and then looks for the best trade-off among competing objectives.

In the next section, we formulate the idealized problem using constraints, and then develop a relaxation that always has a feasible solution. This allows significant modeling and computational flexibility, yielding efficiently computable and clinically useful solutions that balance the need to deliver enough radiation to the tumor while approximately satisfying OAR constraints.

2.1 New formulation

To simplify the exposition, we explain the new formulation using only dose–volume constraints on the OARs:

[TABLE]

Solutions to $(\mathcal{P}_{1})$ may not be clinically useful if the dose–volume constraints prevent sufficient dose delivered to the tumor. We develop an extended, nonconvex formulation that relaxes the dose–volume constraints. The new formulation includes additional variables $w_{j}\approx A_{j}x-d_{j}^{\text{dv}}$, and penalizes deviations of the delivered OAR doses from $w_{j}$:

[TABLE]

In formulation $(\mathcal{P}_{2})$, the auxiliary variables $w_{j}$ are forced to satisfy the dose–volume constraints, while the residuals $A_{j}x-d_{j}^{\text{dv}}$ may not necessarily do so. The weights $\alpha_{j}$ control how closely $w_{j}$ must approximate $A_{j}x-d_{j}^{\text{dv}}$, and as $\alpha_{j}\uparrow\infty$, problem $(\mathcal{P}_{2})$ converges to the idealized problem $(\mathcal{P}_{1})$ in an epigraphical sense.

The epigraph of a function is defined by

[TABLE]

Epigraphs make it possible to go between sets and functions, and the notion of epigraphical convergence [88] is exactly that of set convergence for epigraphs, i.e.,

[TABLE]

We know that our objective functions (7a) converge epigraphically to the original problem by [88, Theorem 7.4(d)] because of monotonicity,

[TABLE]

which we always have when using quadratic penalties for prox-bounded functions (see discussion below Theorem 7.4(d)).

Formulation $(\mathcal{P}_{2})$ has several desirable features. First, while most heuristic approaches described in the introduction rely on convex approximations, $(\mathcal{P}_{2})$ fully captures the nonconvex structure of dose–volume constraints, while giving the modeler flexibility to match them approximately using weights $\alpha_{j}$. The weights $\alpha_{i}$ and $\alpha_{j}$ can be tuned to balance the trade-off between covering the tumor and approximately meeting OAR constraints. Second, from a computational perspective, we can design efficient and provably convergent algorithms that aggregate information from multiple tumors and OARs using partial minimization, as described in the next section. Furthermore, approximate solutions can be computed efficiently, then refined using techniques such as iterative re-weighting schemes or polishing.

3 Algorithmic Approach

Our main strategy is to partially minimize formulation $(\mathcal{P}_{2})$ with respect to $x$, viewing $(\mathcal{P}_{2})$ as a value-function optimization problem

[TABLE]

The algorithm is centered on optimizing the value function $g(w)$ in (11a), using a projected gradient approach to take care of the nonconvex constraints encoded by $\Omega_{j}$. The high-level iteration is simply

[TABLE]

where $w=(w_{j_{1}},w_{j_{2}},\dots)$, $T=\text{diag}(t_{j_{1}},\dots,t_{j_{2}},\dots)$, and $\Omega=\Omega_{j_{1}}\times\Omega_{j_{2}}\times\dots$ for $j_{i}\in\mathcal{J}^{\text{dv}}$. Evaluation of $g$ and $\nabla g$ is done through computing the partial minimum over $x$, as detailed in Algorithm 1. Projection onto $\Omega_{j}$ can be implemented efficiently by ordering the elements of $w_{j}$ and setting the lowest $\left\lceil\nicefrac{{(100-p_{j})n_{j}}}{{100}}\right\rceil$ entries of $w_{j}$ to $\min\{0,w_{j}\}$ [56].

The differentiability of $g(w)$, formula for the derivative, and its Lipschitz constant follow from [89, Theorem 1]. In particular, for one PTV and one OAR, define

[TABLE]

Next, taking $h(x)=\delta_{+}(x)+\frac{1}{2}\|A_{1}x-d_{1}\|_{2}^{2}+\frac{\lambda}{2}\|x\|_{2}^{2}$, where $\delta_{+}(x)$ is the indicator function of the nonnegative orthant,

[TABLE]

we have

[TABLE]

and from [89, Theorem 1], we get

[TABLE]

Similarly, for multiple OARs we have $\nicefrac{{\partial g}}{{\partial w_{j}}}(w)=\nicefrac{{\alpha_{j}}}{{n_{j}}}\left\{w_{j}-\left[A_{j}x(w)-d_{j}^{\text{dv}}\right]\right\}$ and $\text{lip}\left(\nicefrac{{\partial g}}{{\partial w_{j}}}\right)\leq\nicefrac{{\alpha_{j}}}{{n_{j}}}$ for all $j\in\mathcal{J}^{\text{dv}}$, and for all $w,w^{\prime}\in\mathbb{R}^{\sum_{j\in\mathcal{J}^{\text{dv}}}n_{j}}$ we have

[TABLE]

so $\text{lip}\left(\nabla g\right)\leq\max_{j\in\mathcal{J}^{\text{dv}}}\left(\nicefrac{{\alpha_{j}}}{{n_{j}}}\right)$.

Algorithm 1 is projected gradient descent for the value function (11a), and is equivalent to the simple iteration (12). Under fairy weak assumptions111Specifically, the sequence $\{w^{(k)}\}_{k\in\mathbb{N}}$ generated by (12) must be bounded, and $g+\delta_{\Omega}$ must be a Kurdyka-Łojasiewicz function (which covers a wide range of functions, including the class of real semi-algebraic functions)., the nonconvex gradient projection method converges for a fixed step size $t<\nicefrac{{1}}{{\text{lip}(\nabla g)}}$ [90, Theorem 5.3]. On the other hand, Algorithm 1 allows the flexibility to choose different step sizes $t_{j}\leq\nicefrac{{n_{j}}}{{\alpha_{j}}}$ for each $j\in\mathcal{J}^{\text{dv}}$. We prove the convergence of this variant of projected gradient descent, starting with some definitions to quantify stationarity and optimality.

Definition 1** (Stationary Point).**

A vector $w\in\Omega$ is called a stationary point for $(\mathcal{P}_{3})$ if

[TABLE]

where $N_{\Omega}$ is the normal cone to $\Omega$ at $w$.

This definition follows from the basic first-order conditions for optimality [88, Theorem 6.12], where (22) is a necessary condition for $w$ to be a local minimizer of $g$ on $\Omega$. This definition of stationarity motivates the following measure of optimality.

Definition 2** ($\epsilon$-accurate Solution).**

Given fixed step sizes $t_{j}>0$ for all $j\in\mathcal{J}^{\text{dv}}$, we call a vector $w\in\Omega$ an $\epsilon$-accurate solution for $(\mathcal{P}_{3})$ if

[TABLE]

where $\bar{w}_{j}\in\text{proj}_{\Omega_{j}}\left[w_{j}-t_{j}\nicefrac{{\partial g}}{{\partial w_{j}}}(w)\right]$.

By the first-order optimality conditions of the problem defining the projection in (12), for each $j\in\mathcal{J}^{\text{dv}}$ the vectors $w_{j}$ and $\bar{w}_{j}$ satisfy

[TABLE]

so if $\sum_{j\in\mathcal{J}^{\text{dv}}}\nicefrac{{1}}{{t_{j}^{2}}}\left\|w_{j}-\bar{w}_{j}\right\|_{2}^{2}=0$, then $0\in\nicefrac{{\partial g}}{{\partial w_{j}}}(w)+N_{\Omega_{j}}(w_{j})$ for all $j\in\mathcal{J}^{\text{dv}}$. Furthermore, by [88, Proposition 6.41], we have $0\in\nabla g(w)+N_{\Omega}(w)$, so $w$ is a stationary point for problem $(\mathcal{P}_{3})$.

We now establish the convergence of Algorithm 1, with proofs provided in Appendix B.

Theorem 1**.**

For fixed step sizes $t_{j}\in\left(0,\nicefrac{{n_{j}}}{{\alpha_{j}}}\right]$ for all $j\in\mathcal{J}^{\text{dv}}$, Algorithm 1 generates a sequence of iterates $\left\{w^{(k)}\right\}_{k\in\mathbb{N}}$ such that $g\left(w^{(k)}\right)\downarrow g^{*}\geq 0$, in particular $g^{*}\leq g\left(w^{(k)}\right)$ for all $k\in\mathbb{N}$. Furthermore, for any convergent subsequence $w^{(p_{k})}\to w_{p}^{*}$, the vector $w_{p}^{*}$ is a local minimizer of $g$ on $\Omega$, and $g\left(w_{p}^{*}\right)=g^{*}$.

We note that the convergence of our sequence of objective values $\left\{g\left(w^{(k)}\right)\right\}_{k\in\mathbb{N}}$ only implies subsequence convergence for our iterates $\left\{w^{(k)}\right\}_{k\in\mathbb{N}}$. In what follows, we prove the convergence of our sequence of iterates to a local minimizer of $g$ on $\Omega$ for the case where $t_{j}<\nicefrac{{n_{j}}}{{\alpha_{j}}}$ for all $j\in\mathcal{J}^{\text{dv}}$.

Theorem 2**.**

For fixed step sizes $t_{j}\in\left(0,\nicefrac{{n_{j}}}{{\alpha_{j}}}\right)$ for all $j\in\mathcal{J}^{\text{dv}}$, the sequence of iterates $\left\{w^{(k)}\right\}_{k\in\mathbb{N}}$ generated by Algorithm 1 converges to a local minimizer of $g$ on $\Omega$. Furthermore, we reach an $\epsilon$-accurate solution in $\mathcal{O}(\nicefrac{{1}}{{\epsilon}}$) iterations.

If we choose $t_{j}=\nicefrac{{n_{j}}}{{\alpha_{j}}}$ for each $j\in\mathcal{J}^{\text{dv}}$, Algorithm 1 reduces to block coordinate descent, detailed in Algorithm 2, and its convergence to a function value corresponding to a local minimizer follows immediately from Theorem 1.

3.1 Simple example

To demonstrate how the relaxed formulation $(\mathcal{P}_{2})$ and Algorithm 2 behave compared to the idealized problem $(\mathcal{P}_{1})$, we revisit our example from the introduction. For iterations $k=0,3,32$ of Algorithm 2 to solve $(\mathcal{P}_{2})$, we plot the contour lines of the objective function

[TABLE]

in terms of the beamlet intensity variable $x$ for fixed $w^{(k)}$. Here $A_{1}\in\mathbb{R}^{1\times 2}$ and $A_{2}\in\mathbb{R}^{2\times 2}$ are the beamlet-to-voxel maps for the prostate tumor and rectum, $x\in\mathbb{R}^{2}$ contains the beamlet intensities, $w^{(k)}\in\mathbb{R}^{2}$ approximates the deviation between the dose received by the rectum and the dose–volume constraint, and $\lambda$ is set at $5\times 10^{-6}$. In each figure, $x^{(k+1)}$ corresponds to the global minimum of the relaxed problem at iteration $k$, A and B correspond to the local and global minima of the idealized problem, and the nonconvex feasible region induced by the dose–volume constraint that no more than 50% of the rectum volume exceeds 20 Gy is given by the union of regions $v_{1}\leq 20$ and $v_{2}\leq 20$.

In Figure 2 (left), we initialize the problem with $x^{(0)}$ as the zero vector. The contours of $(\mathcal{P}_{2})$ are more circular than those of $(\mathcal{P}_{1})$ in Figure 1, indicating that the relaxation has improved the problem conditioning. Additionally, as $w^{(k)}$ changes with each iteration, the global minimum of the objective function in terms of $x$ at point $x^{(k+1)}$ moves closer to the global minimum of $(\mathcal{P}_{1})$ at point B. For stopping tolerance $\epsilon=10^{-3}$, our algorithm converges after 32 iterations, and in Figure 2 (right) we see that the global minimum of our relaxed problem, corresponding to our final beamlet intensities, is just shy of the global minimum of our original problem. The dose–volume constraint is met approximately, and we may choose larger weight coefficient $\alpha_{2}$ or decrease the dose level $d_{2}^{\text{dv}}$ to move closer to the feasible region if necessary.

Next, we plot our iterates along with the contour lines of the objective function

[TABLE]

in terms of the variable $w$. In Figure 3 (left), points A and B again correspond to the local and global minima of the idealized problem, and the nonconvex feasible region $\Omega=\{w\in\mathbb{R}^{2}\,:\,\|(w)_{+}\|_{0}\leq 1\}$ is given by the union of regions $w_{1}\leq 0$ and $w_{2}\leq 0$. The result of the gradient step is plotted as squares, and the projection onto $\Omega$ is plotted as circles. Just as we saw in $x$-space, our algorithm in $w$-space converges to a point near the global minimum at B.

3.2 Initialization

Because our formulation is nonconvex, the initialization can affect both the solution and the time to convergence. For example, using the initialization $x^{(0)}=(0,2\times 10^{3})$ for our simple example, we converge near the local minimum at A after 37 iterations (middle panel of Figure 3). Fortunately, as demonstrated in other papers (see e.g., [83, 84]), local minima do not have a large effect on the solutions to these least-squares formulations, especially when the convex tumor terms have weights relative in size to the OAR term weights [91]. Given enough time, different initializations tend to converge to similar objective values, fluence maps, and dose distributions.

For the remainder of our examples, we initialize the beamlet intensity vector as the solution to $(\mathcal{P}_{2})$ without any dose–volume terms, i.e.,

[TABLE]

In Figure 3 (right), we use this initialization for our simple example, converging near the global minimum at B after 26 iterations. For this initialization, we see that $A_{2}x^{(0)}-20$, the value of $w^{(0)}$ before projection onto $\Omega$, corresponds to the global minimum of $g(w)$. This is true for the general problem $(\mathcal{P}_{2})$ as well. Specifically, for all $j\in\mathcal{J}^{\text{dv}}$ we have

[TABLE]

which means that

[TABLE]

Therefore letting $x^{(0)}=\mathop{\rm arg\,min}_{x\geq 0}\left(\min_{w}f(x,w)\right)$, we are guaranteed that

[TABLE]

for all $j\in\mathcal{J}^{\text{dv}}$ because

[TABLE]

Since we are ultimately interested in solving $\min_{w\in\Omega}g(w)$, using an initialization for $x^{(0)}$ that corresponds to $w^{(0)}=\text{proj}_{\Omega}(\mathop{\rm arg\,min}_{w}g(w))$ is a natural choice.

4 Numerical Results

Our examples are performed on data from an anonymized prostate cancer patient included in the CORT dataset [5], available in the GigaScience repository, GigaDB [6]. The dataset includes beamlet-to-voxel maps for beam angles ranging from $0^{\circ}$ to $358^{\circ}$ in increments of $2^{\circ}$, PTVs in the prostate and lymph nodes, and OARs including the rectum, bladder, femoral heads, and unspecified tissue (Figure 4). In the CORT dataset, beamlet intensities are measured in monitor units (MU), where 100 MU delivers a dose of 1 Gy in 10 cm depth in water in the center of a 10 cm $\times$ 10 cm radiation field [5].

For a fixed set of beams, we consider the problem of treating tumors with dose–volume constraints of increasing difficulty on various combinations of critical structures. In Section 4.1, we look at an example with one PTV and one OAR, and show how the proposed approach works for a single dose–volume constraint on the OAR. In Section 4.2, we increase the complexity of the previous example by imposing multiple dose–volume constraints on the same PTV and OAR. In Section 4.3, we show how the approach handles multiple PTVs and OARs. Finally, in Section 4.4, we compare our approach to other published methods. For all of our examples, we use seven equally spaced beams with angles ranging from $0^{\circ}$ to $312^{\circ}$ in increments of $52^{\circ}$. We set our regularization parameter at $\lambda=10^{-8}$, our stopping tolerance at $\epsilon=10^{-3}$, and for simplicity we let $\alpha_{i}=\alpha_{j}=1$ for all $i\in\mathcal{I}$ and all $j\in\mathcal{J}^{\text{dv}}$, unless otherwise stated. Whenever there is overlap between PTVs and OARs, the voxels in common are assigned to the PTV. However, this choice is not required for our model, and a single voxel may be assigned to multiple structures if desired. All examples were run on a computer with a 2.9 GHz dual-core Intel Core i5 processor with 8 GB RAM.

We use the cumulative dose–volume histogram as a way of evaluating the quality of the resulting treatment plans. In a dose–volume histogram, each point represents the percent of the organ volume that receives at least a particular dose. For example, we aim to deliver a uniform dose of 81 Gy to the prostate tumor, corresponding to the dose–volume histogram on the left panel of Figure 5. In this case, 100% of the organ volume receives exactly 81 Gy. As we cannot deliver a perfectly uniform dose, our goal is to find a treatment plan that is as close as possible to this idealized dose profile. For instance, in the right panel of Figure 5 we see the dose–volume histogram of the solution to (27) for one prostate tumor with regularization parameter $\lambda=10^{-8}$. We use the beamlet intensities from this solution as our initialization for our examples in Sections 4.1, 4.2, and 4.4.1.

Cumulative dose–volume histograms can be used to determine whether or not an upper dose–volume constraint on an OAR has been met by a particular treatment plan. For example, we consider the dose–volume constraint that no more than 50% of the rectum volume exceeds 50 Gy. In Figure 6, the rectangular region within the dotted lines contains all points where at most 50% of the rectum volume receives more than 50 Gy. Therefore any dose curve that lies outside of the box does not meet the constraint, while any curve that travels inside the box meets the constraint. In the left panel of Figure 6, approximately 57% of the rectum volume receives more than 50 Gy, so the dose–volume constraint is not met. However, on the right panel of Figure 6 approximately 42% of the rectum volume receives more than 50 Gy, so the dose–volume constraint is met.

4.1 One PTV and one OAR with one dose–volume constraint

For our first example, we consider the problem of delivering a uniform dose of 81 Gy to a prostate tumor, while satisfying the dose–volume constraint that no more than 30% of the rectum volume receives more than 30 Gy. In this case, the PTV has 6770 voxels, the OAR has 1764 voxels, and there are 986 beamlets. In Figure 7, we see a representative axial slice of the calculated dose along with the intensities of four of the beams (due to symmetry, the remaining three beams have similar intensity patterns). In this slice, a nearly uniform dose of 81 Gy is delivered to the tumor, while only a small percentage of the rectum volume receives more than 30 Gy.

Figure 8 shows the behavior of our model using Algorithm 2. The objective $g(w)$ decreases monotonically (top left), and Algorithm 2 converges after 42 iterations. The trade-off between the competing PTV and OAR objectives is shown in the remaining left panels of Figure 8. Because our initial beamlet intensities correspond to a solution with no dose–volume constraints, at initialization the tumor term achieves its lowest value and the rectum term achieves its highest value. Once we begin iterating, the inclusion of the OAR term causes the tumor term value to increase and the rectum term to decrease. On the right panel of Figure 8, we show the behavior of the auxiliary variable $w$ (top right) and our approximation of the dose–volume constraint (bottom right). The percent of rectum voxels exceeding 30 Gy approaches the dose–volume constraint; however, because of the relaxation, we do not meet the dose–volume constraint exactly. Nevertheless, we improve upon the initialization value according to the trade-off implicitly specified by our weights $\alpha_{1}$ and $\alpha_{2}$.

To demonstrate how the weight factors and dose–volume parameters can influence solutions, we implement a variation of the penalty decomposition method for the convex-cardinality problem proposed in [92]. Specifically, we use the iterative re-weighting scheme given in Algorithm 3, where ($\mathcal{P}_{2}$) is solved multiple times with increasing OAR weights, decreasing dose–volume parameters, and decreasing stopping tolerance values. In this example, we use initial OAR weight $\alpha_{2}^{(0)}=1$, dose parameter $d_{2}^{\text{dv}(0)}=30$, volume parameter $p_{2}^{(0)}=30$, and stopping tolerance $\epsilon^{(0)}=10^{-3}$, with weight update parameter $\sigma=0.01$ and stopping tolerance update parameter $\gamma=0.99$.

We compare the approximation of the dose–volume constraint in our solutions with and without re-weighting in Table 1 and Figure 9. To illustrate the trade-off between the competing PTV and OAR objectives, we include results from Algorithm 3 using different convergence criteria: (a) stopping once the dose–volume constraint is met, and (b) stopping once the PTV D95 (i.e., the minimum dose delivered to 95% of the PTV) is less than 98% of its initial value, where weight factors and dose–volume parameters are updated regardless of whether or not the dose–volume constraint has been satisfied. There is little difference between the initial tumor dose and the final tumor dose for each solution, both in terms of the PTV D95 and the dose–volume histogram. At the same time, all solutions have reduced the dose to the rectum to approximately satisfy the dose–volume constraint, where re-weighting produces a lower dose at the expense of increased computation time. Without re-weighting, the constraint is roughly met, with only 34.16% of the volume receiving more than 30 Gy, in contrast to the initialization, where 64.14% of the volume exceeds 30 Gy. With re-weighting, the constraint is satisfied, with 29.61% and 19.84% of the volume receiving more then 30 Gy for stopping conditions (a) and (b), respectively.

In general, as the dose to the rectum decreases, the PTV D95 decreases, while the dose to the tumor deviates more from the target uniform dose. If the dose–volume constraint makes the problem completely infeasible, the approach still returns some solution, but it need not meet the dose–volume constraint nor the target tumor dose. In this case, the planner can influence the trade-off between treatment goals by modifying weighting factors, dose–volume parameters, and convergence criteria.

4.2 One PTV and one OAR with multiple dose–volume constraints

Next we consider the problem of delivering a uniform dose of 81 Gy to the prostate tumor, with multiple dose–volume constraints on the tumor and rectum:

• •

No more than 5% of the tumor volume receives less than 81 Gy,

• •

no more than 0% of the tumor volume receives more than 85 Gy,

• •

no more than 20% of the rectum volume receives more than 60 Gy,

• •

no more than 40% of the rectum volume receives more than 40 Gy, and

• •

no more than 60% of the rectum volume receives more than 20 Gy.

The first constraint, a lower dose–volume constraint on the PTV, is implemented in the same way as the upper dose–volume constraint on the OAR, with the exception that the corresponding term in the objective function is changed to

[TABLE]

and the projection onto $\Omega_{i}$ is computed by ordering the elements of $w_{i}$ and setting the highest $\lceil\nicefrac{{(100-p_{i})n_{i}}}{{100}}\rceil$ entries of $w_{i}$ to $\max\{0,w_{i}\}$. The second constraint corresponds to a maximum-dose constraint, but expressing it as an upper dose–volume constraint allows us to apply our relaxation rather than enforcing a hard constraint.

As in the previous example, the PTV has 6770 voxels, the OAR has 1764 voxels, and there are 986 beamlets. One downside to adding multiple dose–volume constraints on a particular structure is the inclusion of more terms in the objective function, increasing the problem dimension. However, we can reduce the number of objective terms by combining each of the $K_{j}$ dose–volume constraint terms corresponding to organ $j$,

[TABLE]

into one term,

[TABLE]

where $\tilde{\Omega}_{j}=\tilde{\Omega}_{j_{1}}\cap\cdots\cap\tilde{\Omega}_{j_{K_{j}}}$ represents the set of vectors that satisfy all of the organ’s $K_{j}$ dose–volume constraints, with each constraint defined as

[TABLE]

Projections onto $\tilde{\Omega}_{j}$ can be computed efficiently by ordering the elements of $y_{j}$, setting the lowest $(100-p_{j_{k}})$% of $y_{j}$ to $\min\{d_{j_{k}},y_{j}\}$ for upper dose–volume constraints, and setting the highest $(100-p_{j_{k}})$% of $y_{j}$ to $\max\{d_{j_{k}},y_{j}\}$ for lower dose–volume constraints, for all $k=1,\dots,K_{j}$.

Using Algorithm 2 with $\alpha_{1}=10$ for the PTV dose–volume constraint term, we converge after 21 iterations. In Table 2 and Figure 10, we see that the dose–volume constraints have been approximately met. Specifically, none of the PTV volume receives more than 85 Gy, 24.47% of the PTV volume receives less than 81 Gy (53.26% initially), and the PTV D95 increases from 79.65 Gy to 80.26 Gy. Furthermore, 61.95% of the OAR volume exceeds 20 Gy (68.81% initially), 42.9% exceeds 40 Gy (61.29% initially), and 21.84% exceeds 60 Gy (38.47% initially). Importantly, it is possible to specify the shape of the PTV and OAR dose–profiles in terms of multiple dose–volume constraints without any significant increase in computation.

Next, we use Algorithm 3 with weight factors and dose–volume parameters initialized with the values used in Algorithm 2, initial stopping tolerance $\epsilon^{(0)}=10^{-3}$, weight update parameter $\sigma=0.01$, and stopping tolerance update parameter $\gamma=0.99$. For the lower dose–volume constraint on the PTV, the dose parameter increases each iteration,

[TABLE]

rather than decreasing like the upper dose–volume constraints. Stopping once all dose–volume constraints are met, none of the PTV volume receives more than 85 Gy, 4.09% of the PTV volume receives less than 81 Gy, and the PTV D95 is 81.32 Gy. In addition, only 59.1% of the OAR volume exceeds 20 Gy, 39.26% exceeds 40 Gy, and 19.78% exceeds 60 Gy. In cases of infeasibility due to competing PTV and OAR dose–volume constraints, additional convergence criteria for Algorithm 3 can be specified, e.g., setting a maximum number of iterations.

4.3 Multiple PTVs and OARs

Now we consider the problem of delivering a uniform dose of 81 Gy to the prostate tumor and a uniform dose of 60 Gy to the lymph nodes, while satisfying the dose–volume constraints that no more than 50% of the rectum receives more than 50 Gy and no more than 30% of the bladder receives more than 30 Gy. In this case the first PTV has 6770 voxels, the second PTV has 8483 voxels, the first OAR has 1643 voxels, the second OAR has 9708 voxels, and there are 986 beamlets. We converge after 48 iterations of Algorithm 2.

Due to the proximity of the two PTVs and the difference between their target doses, it is much more difficult to deliver a uniform dose to both targets. The increase in geometric constraints and competing objectives results in more complicated dose distributions and beamlet intensity patterns, illustrated in Figure 11. More of the prostate volume receives less than 81 Gy than in previous examples, while more of the lymph node volume exceeds 60 Gy than we would like (Figure 12, top row). In this case, it may be important to add lower and upper dose–volume constraints on the PTVs to ensure that appropriate doses are delivered. The dose–volume constraints on the rectum and bladder are both approximately met, with 57.33% of the rectum exceeding 50 Gy (82.47% initially) and 38.14% of the bladder exceeding 30 Gy (91.87% initially) (Figure 12, bottom row).

Using Algorithm 3 with weight factors and dose–volume parameters initialized with the values used in Algorithm 2, initial stopping tolerance $\epsilon^{(0)}=10^{-3}$, weight update parameter $\sigma=0.01$, and stopping tolerance update parameter $\gamma=0.99$, we are able to meet both of the dose–volume constraints on the OARs. Specifically, only 47.53% of the rectum volume receives more than 50 Gy and 29.93% of the bladder volume exceeds 30 Gy, with little difference in the PTV doses.

4.4 Comparisons with other methods

Finally, we compare results from Algorithm 2, Algorithm 3, and four other methods for solving the FMO problem with dose–volume constraints. We start with two methods that, like our approach, approximate dose–volume constraints with objective terms, and thus provide flexibility for cases of infeasiblity. In particular, as mentioned in Section 1.3, many least-squares models include a penalty term for OAR voxels that exceed a specified dose (see e.g., [10, 21, 31, 34, 48, 84, 85, 93]). This penalty may apply to the entire OAR, or in the case of dose–volume objectives, to a fraction of the OAR volume. Following the approach used in [34, 48, 84], we replace our dose–volume terms in (7a) with

[TABLE]

where the vector $\tilde{d}_{j}^{\text{dv}}$ changes each iteration. Specifically, after each iteration we sort the voxel dose vector $A_{j}x$, assigning the highest $p_{j}$% to the dose value $\tilde{d}_{j}^{\text{dv}}=10^{6}$ and the remaining $(100-p_{j})$% to $\tilde{d}_{j}^{\text{dv}}=d_{j}^{\text{dv}}$. The beamlet vector $x$ is updated each iteration with a step of gradient descent, where only the the lowest $(100-p_{j})$% of voxels that exceed $\tilde{d}_{j}^{\text{dv}}$ influence the direction of the gradient, followed by a projection onto the nonnegative orthant. If the vector $\tilde{d}_{j}^{\text{dv}}$ converges, then the problem becomes convex [84].

Alternatively, we can penalize voxels that exceed a given dose by adding nonnegative slack variables to the OAR terms. For example, following the model used in the Sensitivity-Driven Greedy Algorithm [85, 93], we replace our dose–volume terms in (7a) with

[TABLE]

which is similar to our modified OAR terms in (34). The alternating projections approach of the Sensitivity-Driven Greedy Algorithm is equivalent to the block coordinate descent approach of Algorithm 2, with two important differences. First, the nonnegative least-squares step to compute $x^{(k+1)}$ also includes the slack variables $s_{j}^{(k+1)}$, so the number of variables in this subproblem increases by $\sum_{j\in\mathcal{J}^{\text{dv}}}n_{j}$. Second, to ensure the convergence of the dose variables, the $y^{(k+1)}$ update includes a projection onto the modified set $\tilde{\Omega}_{j}\cap\{y_{j}\in\mathbb{R}^{n_{j}}\,:\,y_{j}\geq y_{j}^{(k)}\}$. Using the initialization $y_{j}^{(0)}=d_{j}^{\text{dv}}$, this approach can be interpreted as iteratively selecting the voxels that may exceed $d_{j}^{\text{dv}}$, and once $\lfloor\nicefrac{{n_{j}p_{j}}}{{100}}\rfloor$ elements of vector $y_{j}^{(k)}$ surpass $d_{j}^{\text{dv}}$, the problem becomes convex.

Next we look at two methods that solve a convex problem by imposing maximum-dose constraints on a chosen subvolume of each OAR. As mentioned in Section 1.3, many methods take a multi-stage approach of 1) solving an approximation of the original problem, 2) determining which voxels may exceed $d_{j}^{\text{dv}}$ based on the resulting solution, and 3) polishing the solution from Step 1 by solving a convex problem using the subvolumes determined in Step 2. For example, one heuristic to choose the $p_{j}\%$ of the voxels that may receive more than $d_{j}^{\text{dv}}$ Gy is to solve the FMO problem without dose–volume constraints [18, 19, 50]. Following this approach, we first solve $(\mathcal{P}_{1})$ without any dose–volume constraints (which is equivalent to our initialization in (27)), and then re-solve with maximum-dose constraints applied to the $(100-p_{j})$% of the voxels that received the lowest dose in the previous solution.

Another heuristic is to solve an approximation of $(\mathcal{P}_{1})$ where the dose–volume constraints in (6b) are replaced by the convex constraint

[TABLE]

where $\nu_{j}$ is optimized along with the beamlet vector $x$ [46]. If this problem is feasible, the solution is guaranteed to satisfy the dose–volume constraints, but it may not be optimal for $(\mathcal{P}_{1})$. Therefore, the problem is polished by re-solving with maximum-dose constraints applied to the $(100-p_{j})$% of the voxels that received the lowest dose in the approximate solution. We compare results from these two convex methods with solutions from our re-weighting scheme in Algorithm 3 and solutions from Algorithm 2 which are then polished to satisfy the dose–volume constraints.

4.4.1 One PTV and one OAR.

First we consider the problem from Section 4.1 with stopping tolerance $\epsilon=10^{-2}$, where an additional lower dose–volume constraint that no more than 5% of the tumor volume receives less than 81 Gy has been applied to the PTV. In Table 4 and Figure 13 (found in Appendix C), we see the results of using approximate dose–volume objectives. Because none of the solutions come close to satisfying the lower dose–volume constraint on the PTV, and only one satisfies the upper dose–volume constraint on the OAR, we look for the best candidate method for re-weighting or polishing. Compared to the solution with (37), our solution achieves a higher PTV D95 and a lower OAR dose in less time. Because the method with (37) does not solve a subproblem to update the beamlet vector, it unsurprisingly requires less time per iteration than the other two methods, though it does require more iterations overall. The method with (38) requires significantly more time than the other two approaches due to the additional slack variables, but it produces the only solution that meets the dose–volume constraint on the OAR. However, its lower OAR dose comes at the expense of a lower PTV D95, and in less time we can use re-weighting or polishing with $(\mathcal{P}_{2})$ to meet all of the dose–volume constraints (first two rows of Table 5).

Next we examine results computed through either re-weighting or polishing. All four solutions in Table 5 and Figure 14 satisfy the dose–volume constraints on the PTV and OAR, and are thus feasible solutions of $(\mathcal{P}_{1})$. However, the lowest objective value for (6a) is achieved by solving $(\mathcal{P}_{2})$ with Algorithm 2, then re-solving with maximum-dose constraints applied to subvolumes of the PTV and OAR. The solution with (27) has the highest objective value, perhaps because it is the only method that does not include any approximations of the dose–volume constraints in Step 1, while the solution with (39) requires the most time. The solution computed with Algorithm 3, on the other hand, requires the least amount of time, though it is not as good with respect to the objective value as the solution computed with Algorithm 2 and polishing.

4.4.2 Two PTVs and two OARs.

Next we consider the problem from Section 4.3 with stopping tolerance $\epsilon=10^{-2}$, with additional lower dose–volume constraints that no more than 5% of the tumor volumes receives less than 81 Gy and 60 Gy for the prostate and lymph nodes, respectively. In Table 6 and Figure 15, we see the results of using approximate dose–volume objectives, with patterns similar to those in Table 4. Specifically, none of the solutions come close to meeting the lower dose–volume constraints on the PTVs, our solution from Algorithm 2 achieves the highest PTV D95 values and requires the least amount of time, and the solution with (38) is the only method to satisfy the upper dose–volume constraints on the OARs, at the expense of significantly more time. While the method with (37) produces a lower dose to the bladder, our method is still the best candidate for re-weighting or polishing, where it achieves solutions that meet all dose–volume constraints in less time than it takes the method with (38) to compute its approximate solution (first two rows of Table 7).

We see the results for this example computed through either re-weighting or polishing in Table 7 and Figure 16. This time, the PTV and OAR subvolumes determined through the solution computed with (27) does not lead to a feasible solution for $(\mathcal{P}_{1})$, though the subvolumes chosen by the other three methods do. Additionally, the method with (39) delivers the lowest dose to the OARs, at the expense of a higher objective value and computation time, while the solution found using $(\mathcal{P}_{2})$ and re-weighting produces the lowest objective value in the least amount of time.

4.4.3 Two PTVs and four OARs.

For our last example, we use the same PTVs, OARs, and prescriptions as the previous section, with the exception of the bladder, where the upper dose–volume constraint is changed so that at most 40% of the organ volume may receive more than 40 Gy. Moreover, we add the left and right femoral heads. For illustration purposes, these new OARs have maximum-dose constraints of 30 Gy, expressed as upper dose–volume constraints that no more than 0% of the OAR volumes receive more than 30 Gy. The left and right femoral heads have 5957 and 5974 voxels, respectively.

Tables 8 and 9 and Figure 17 contain the results computed using approximate dose–volume objectives. While the solution from (38) achieves the lowest OAR doses, meeting the dose–volume constraints for both the rectum and bladder, the method scales poorly. In particular, the addition of the femoral heads leads to an unreasonably long computation time. The other two methods produce solutions with similar PTV D95 doses, though neither do well at approximating the lower dose–volume constraints on the PTVs. On the other hand, the solution from Algorithm 2 requires less time and delivers a lower dose to the OARs than the method with (37).

Finally, Table 10 and Figure 18 include the results computed using either re-weighting or polishing. In this case, the solver for the convex approximation of the dose–volume constraints in (39) fails to find a feasible solution. The other three methods lead to beamlets that satisfy all of the dose–volume constraints. The subvolumes chosen by the initialization in (27) produces the solution with the lowest OAR doses, at the expense of a higher objective value and computation time, while the solutions using $(\mathcal{P}_{2})$ with re-weighting and polishing require the least amount of time and achieve the lowest objective value, respectively. Overall, these results demonstrate that, compared to other convex and nonconvex methods, our approach can efficiently deliver similar or better solutions to the FMO problem in the presence of competing PTV and OAR dose–volume constraints.

5 Conclusions and Future Research

The main challenge of radiation treatment planning is to deliver a prescribed dose to the tumor while sparing surrounding healthy tissues. The fluence map optimization problem for intensity-modulated radiation therapy can be formulated as a large-scale inverse problem with competing objectives and constraints on the planning target volumes and organs-at-risk. Unfortunately, the clinically relevant dose–volume constraints are nonconvex, so solving treatment planning problems exactly with dose–volume constraints is NP-hard. We proposed a new approach that is well adapted to handle nonconvex dose–volume constraints without relying on convex relaxations. To solve this formulation, we developed a customized algorithm that is guaranteed to converge to a stationary point of the proposed model. The overall approach is scalable, handles multiple objectives on PTVs and OARs, and returns actionable treatment plans as demonstrated with examples on the CORT dataset.

While this paper focuses primarily on the methodology and feasibility of our new approach, future work on this project will be directed towards practical considerations. We will test our model on additional datasets including different tumor types, patient geometries, and prescriptions. The comparison of the results using our approach and commercial treatment planning systems is left for future work.

Acknowledgements

This work is in part supported by NSF Grant DGE-1258485 AM002 and the Washington Research Foundation Data Science Professorship. The first author acknowledges the generous support of the University of Washington Applied Mathematics Boeing Fellowship.

Appendix A Availability of Materials

The figures in this article, as well as the plotting scripts necessary to reproduce them, are available openly under the CC-BY license [94]. All examples are performed on the CORT prostate dataset [5, 6]. Solutions to nonnegative least-squares problems are solved with the function minConf_TMP from the minConf package [95, 96], and solutions for (39) in Section 4.4 were solved with CVX, a package for specifying and solving convex programs [97, 98].

Appendix B Proofs for Section 3

In this section, we provide proofs for Theorems 1 and 2, which follow from several auxiliary lemmas. For simplicity, we let $m$ be the number of beamlets and $N:=\sum_{j\in\mathcal{J}^{\text{dv}}}n_{j}$ be the number of dose–volume constrained voxels, so that $x\in\mathbb{R}^{m}$ and $w\in\mathbb{R}^{N}$. We begin by proving that all fixed points of Algorithm 1 are local minimizers of $g$ on $\Omega$.

Lemma 1**.**

The objective $f(x,w)$ is strongly convex on $\mathbb{R}^{m}\times\mathbb{R}^{N}$.

Proof.

For simplicity, we consider the case with one PTV and one OAR, but the arguments used generalize to include additional terms. Rearranging variables, the objective $f(x,w)$ can be expressed using a single quadratic term,

[TABLE]

with symmetric Hessian matrix

[TABLE]

Consider the matrix $\left(\nicefrac{{\alpha_{2}}}{{n_{2}}}\right)I$ and the Schur complement of $\tilde{A}^{T}\tilde{A}$ in $\left(\nicefrac{{\alpha_{2}}}{{n_{2}}}\right)I$,

[TABLE]

Because the coefficients $\nicefrac{{\alpha_{1}}}{{n_{1}}}$, $\nicefrac{{\alpha_{2}}}{{n_{2}}}$, and $\lambda$ are positive, the matrices $\left(\nicefrac{{\alpha_{2}}}{{n_{2}}}\right)I$ and $\tilde{A}^{T}\tilde{A}\,/\left[\left(\nicefrac{{\alpha_{2}}}{{n_{2}}}\right)I\right]$ are positive definite (note that this does not depend on the rank of matrices $A_{1}$ or $A_{2}$), which implies that the Hessian matrix $\tilde{A}^{T}\tilde{A}$ is also positive definite [73]. Therefore the objective $f(x,w)$ with one PTV and one OAR is strongly convex on $\mathbb{R}^{m}\times\mathbb{R}^{n_{2}}$.

For multiple PTVs and OARs, the lower-right block of the Hessian matrix becomes a block diagonal matrix with matrices $\left(\nicefrac{{\alpha_{j}}}{{n_{j}}}\right)I$ for $j\in\mathcal{J}^{\text{dv}}$ along the diagonal, and the Schur complement becomes the matrix $\sum_{i\in\mathcal{I}}A_{i}^{T}A_{i}+\lambda I$. Because the coefficients $\nicefrac{{\alpha_{i}}}{{n_{i}}}$ for $i\in\mathcal{I}$, $\nicefrac{{\alpha_{j}}}{{n_{j}}}$ for $j\in\mathcal{J}^{\text{dv}}$, and $\lambda$ are positive, the arguments above can again be used to show that the objective $f(x,w)$ with multiple PTVs and OARs is strongly convex on $\mathbb{R}^{m}\times\mathbb{R}^{N}$. ∎

Lemma 2**.**

The objective $g(w)$ is strongly convex on $\mathbb{R}^{N}$.

Proof.

Because the objective $f(x,w)$ is strongly convex on $\mathbb{R}^{m}\times\mathbb{R}^{N}$ [Lemma 1], there exists a constant $\mu>0$ such that for all vector pairs $(x,w),\left(x^{\prime},w^{\prime}\right)\in\mathbb{R}^{m}\times\mathbb{R}^{N}$ and constants $t\in[0,1]$, we have

[TABLE]

Let $x(w)\in\mathop{\rm arg\,min}_{x\geq 0}f(x,w)$. Since $f(x,w)$ is a strongly convex function and the nonnegative orthant is a nonempty closed convex set, the vector $x(w)\in\mathbb{R}_{\geq 0}^{m}$ exists and is uniquely defined for all vectors $w\in\mathbb{R}^{N}$ [73], and $g(w)=\min_{x\geq 0}f(x,w)=f(x(w),w)$. Furthermore, for all vectors $w,w^{\prime}\in\mathbb{R}^{N}$ and constants $t\in[0,1]$, we have

[TABLE]

where line (50) follows from the fact that the nonnegative orthant is a convex set and $x(w)\geq 0$ for all vectors $w\in\mathbb{R}^{N}$, so $tx(w)+(1-t)x\left(w^{\prime}\right)\geq 0$ for all vectors $w,w^{\prime}\in\mathbb{R}^{N}$ and constants $t\in[0,1]$. Therefore $g(w)$ is strongly convex on $\mathbb{R}^{N}$. ∎

Theorem 3**.**

The fixed points of Algorithm 1 are local minimizers of $g$ on $\Omega$.

Proof.

Fixed points $w^{*}=(w_{j_{1}}^{*},w_{j_{2}}^{*},\dots)$ for $j_{i}\in\mathcal{J}^{\text{dv}}$ of Algorithm 1 satisfy

[TABLE]

Based on the form of the projection onto $\Omega_{j}$, for each element $i$ of a stationary vector $w_{j}^{*}$, we know that

if $(w_{j}^{*})_{i}<0$, then $[\nicefrac{{\partial g}}{{\partial w_{j}}}\left(w^{*}\right)]_{i}=(w_{j}^{*})_{i}-(A_{j}x(w^{*})-d_{j}^{\text{dv}})_{i}=0$, and the dose to voxel $i$ is less than $d_{j}^{\text{dv}}$, 2. 2.

if $(w_{j}^{*})_{i}>0$, then $[\nicefrac{{\partial g}}{{\partial w_{j}}}(w^{*})]_{i}=0$, $(w_{j}^{*})_{i}$ is one of the $\lfloor\nicefrac{{n_{j}p_{j}}}{{100}}\rfloor$ largest positive elements of vector $w_{j}^{*}$, and the dose to voxel $i$ is greater than $d_{j}^{\text{dv}}$, and 3. 3.

if $(w_{j}^{*})_{i}=0$, then either $[\nicefrac{{\partial g}}{{\partial w_{j}}}(w^{*})]_{i}=0$, and the dose to voxel $i$ is exactly $d_{j}^{\text{dv}}$, or $-\nicefrac{{1}}{{t_{j}}}[\nicefrac{{\partial g}}{{\partial w_{j}}}(w^{*})]_{i}>0$ is less than the $\lfloor\nicefrac{{n_{j}p_{j}}}{{100}}\rfloor$ largest positive elements of $w_{j}^{*}$, and the dose to voxel $i$ is greater than $d_{j}^{\text{dv}}$ (i.e., voxel $i$ violates the dose–volume constraint on OAR $j$).

Recall that for each $j\in\mathcal{J}^{\text{dv}}$, the set $\Omega_{j}=\left\{w\in\mathbb{R}^{n_{j}}\,:\,\left\|w_{+}\right\|_{0}\leq\nicefrac{{n_{j}p_{j}}}{{100}}\right\}$ is the union of finitely many convex sets, each representing a different combination of the $q_{j}:=\lfloor\nicefrac{{n_{j}p_{j}}}{{100}}\rfloor$ elements of the vector $w_{j}$ allowed to be positive. Specifically, for $\ell=1,\dots,\nicefrac{{n_{j}!}}{{\left(n_{j}-q_{j}\right)!q_{j}!}}$, let $M_{j}^{\ell}$ be the $n_{j}\times n_{j}$ identity matrix with $q_{j}$ of its diagonal elements set to zero. Then the convex set $\Omega_{j}^{\ell}:=\left\{w\in\mathbb{R}^{n_{j}}\,:\,M_{j}^{\ell}w\leq 0\right\}$ represents the set of vectors that satisfy combination $\ell$, and $\Omega_{j}=\cup_{\ell=1,\dots,\nicefrac{{n_{j}!}}{{\left(n_{j}-q_{j}\right)!q_{j}!}}}\Omega_{j}^{\ell}$ is the union of all such convex sets. Consequently, the set $\Omega=\prod_{j\in\mathcal{J}^{\text{dv}}}\Omega_{j}$, the Cartesian product of the feasible sets for vectors $w_{j}$, is also the union of finitely many convex sets, each representing a different combination of the $Q:=\sum_{j\in\mathcal{J}^{\text{dv}}}q_{j}$ elements of $w$ allowed to be positive. Because the objective $g$ is strongly convex on any convex subset of its domain $\mathbb{R}^{N}$ [Lemma 2], there is a unique global minimum of $g$ relative to any of these convex subsets of $\Omega$ [73], each corresponding to a local minimum of $g$ on $\Omega$. As a result, there are finitely many local minimizers of $g$ on $\Omega$.

Suppose $w^{*}$ is a fixed point of Algorithm 1. If $w^{*}\leq 0$, then the vector satisfies all feasible combinations of the $Q$ elements of $w$ allowed to be positive, and thus lies within the intersection of all corresponding convex subsets of $\Omega$. Furthermore, from observations 1 and 3 above, we have $\nabla g\left(w^{*}\right)=0$, and using the strong convexity of $g$, we can conclude that $w^{*}$ is the unique global minimimizer of $g$ on $\Omega$ [88, Theorem 6.12]. On the other hand, suppose that at least one element of $w^{*}$ is positive, and let $M^{*}$ be the $\mathbb{R}^{N}\times\mathbb{R}^{N}$ identity matrix where the diagonal elements corresponding to the at most $Q$ positive elements of $w^{*}$ are set to zero. Then $\Omega^{*}:=\left\{w\in\mathbb{R}^{N}\,:\,M^{*}w\leq 0\right\}$ is a convex subset of $\Omega$ containing $w^{*}$. Again, from observations 1–3 above, we have $[\nabla g\left(w^{*}\right)]_{i}=0$ where $\left(w^{*}\right)_{i}\neq 0$ and $[\nabla g\left(w^{*}\right)]_{i}\leq 0$ where $\left(w^{*}\right)_{i}=0$, so for all $w\in\Omega^{*}$ we have

[TABLE]

because positive elements of $w$ correspond to zero elements of the gradient. Therefore $-\nabla g\left(w^{*}\right)^{T}\in N_{\Omega^{*}}\left(w^{*}\right)$ [88, Theorem 6.9], and by the convexity of $g$ and $\Omega^{*}$, we know that $w^{*}$ is the unique global minimizer of $g$ on $\Omega^{*}$ [88, Theorem 6.12], and hence a local minimizer of $g$ on $\Omega$. ∎

We now establish the convergence of Algorithm 1 to a local minimum of $g$ on $\Omega$.

Lemma 3**.**

For all vectors $w,w^{\prime}\in\mathbb{R}^{N}$, we have

[TABLE]

Proof.

Following the proof from [99, Lemma 1.2.3], for all vectors $w,w^{\prime}\in\mathbb{R}^{N}$, we have

[TABLE]

Recall that for OARs $j\in\mathcal{J}^{\text{dv}}$, we have $\text{lip}\left(\nicefrac{{\partial g}}{{\partial w_{j}}}\right)\leq\frac{\alpha_{j}}{n_{j}}$. Therefore

[TABLE]

which implies that

[TABLE]

for all $w,w^{\prime}\in\mathbb{R}^{N}$. ∎

Lemma 4**.**

For fixed step sizes $t_{j}\in\left(0,\nicefrac{{n_{j}}}{{\alpha_{j}}}\right]$ for all $j\in\mathcal{J}^{\text{dv}}$, the sequence of objective values $\left\{g\left(w^{(k)}\right)\right\}_{k\in\mathbb{N}}$ generated by Algorithm 1 is nonincreasing.

Proof.

For each $j\in\mathcal{J}^{\text{dv}}$, let vectors $w_{j}^{(k+1)},w_{j}^{(k)}\in\Omega_{j}$ satisfy $w_{j}^{(k+1)}\in\text{proj}_{\Omega_{j}}\left[w_{j}^{(k)}-t_{j}\nicefrac{{\partial g}}{{\partial w_{j}}}\left(w^{(k)}\right)\right]$. Then from the definition of the projection operator, we have

[TABLE]

which simplifies to

[TABLE]

Noting that $\Omega\subset\mathbb{R}^{N}$, we can combine inequality (70) with Lemma 56 to get

[TABLE]

Therefore we have $g\left(w^{(k+1)}\right)\leq g\left(w^{(k)}\right)$ for all vectors $w_{j}^{(k+1)},w_{j}^{(k)}$ satisfying

$w_{j}^{(k+1)}\in\text{proj}_{\Omega_{j}}\left[w_{j}^{(k)}-t_{j}\frac{\partial g}{\partial w_{j}}\left(w^{(k)}\right)\right]$ with $t\in\left(0,\nicefrac{{n_{j}}}{{\alpha_{j}}}\right]$ for all $j\in\mathcal{J}^{\text{dv}}$. ∎

Proof.

(Proof of Theorem 1.) From Lemma 4, we know that if $t_{j}\in\left(0,\nicefrac{{n_{j}}}{{\alpha_{j}}}\right]$ for all $j\in\mathcal{J}^{\text{dv}}$, then the sequence of objective values $\left\{g\left(w^{(k)}\right)\right\}_{k\in\mathbb{N}}$ is nonincreasing, and $g(w)\geq 0$ for all $w\in\mathbb{R}^{N}$ by definition, so the sequence converges to some constant $g^{*}\geq 0$ where $g^{*}\leq g\left(w^{(k)}\right)$ for all $k\in\mathbb{N}$. Furthermore, for all convergent subsequences $w^{(p_{k})}\to w_{p}^{*}$ as $k\to\infty$, we have $g\left(w_{p}^{*}\right)=g^{*}$. Because $w_{p}^{*}$ is a fixed point of Algorithm 1, it must also be a local minimizer of $g$ on $\Omega$ [Theorem 3], so $g^{*}$ is a local minimum of $g$ on $\Omega$. ∎

Finally, we prove the convergence of the iterates of Algorithm 1 to a local minimizer of $g$ on $\Omega$ for the case when $t_{j}<\nicefrac{{n_{j}}}{{\alpha_{j}}}$ for all $j\in\mathcal{J}^{\text{dv}}$.

Proof.

(Proof of Theorem 2.) From Theorem 1 we have $g\left(w^{(k)}\right)\downarrow g^{*}$ as $k\to\infty$, so for all $\epsilon>0$, there exists some positive integer $K$ such that for all $k\geq K$, we have $g\left(w^{(k)}\right)-g^{*}<\epsilon$. Furthermore, from Lemma 2 we know that $g$ is strongly convex. In particular, there exists some scalar $\mu>0$ such that

[TABLE]

for all vectors $w,w^{\prime}\in\Omega$. Let $w^{*}$ be a fixed point of Algorithm 1 that satisfies $g(w^{*})=g^{*}$. Then for any $k\geq K$, we have

[TABLE]

where $\nabla g(w^{*})^{T}w^{*}=0$ because nonzero elements of $w^{*}$ correspond to zero elements of $\nabla g(w^{*})$ (see observations 1 and 2 from Theorem 3). We consider two cases based on the gradient of $g$ at $w^{*}$.

Case 1. Suppose $\nabla g(w^{*})=0$, and set $\delta=\sqrt{\nicefrac{{2\epsilon}}{{\mu}}}$. Then from Theorem 1 and inequality (78), we know that for all $\delta>0$, there exists some positive integer $K$ such that for all $k\geq K$, we have

[TABLE]

Therefore $w^{(k)}\to w^{*}$ as $k\to\infty$, and the sequence of iterates $\left\{w^{(k)}\right\}_{k\in\mathbb{N}}$ converges to a fixed point of Algorithm 1. In fact, because $g$ is strongly convex and $\nabla g(w^{*})=0$, the fixed point $w^{*}$ is the unique global minimizer of $g$ on $\Omega$, and $g(w^{*})=g^{*}$ is the global minimum.

Case 2. Suppose $[\nabla g(w^{*})]_{i}\neq 0$ for at least one element $i$. From observation 3 of Theorem 3, this will occur exactly when $-\nicefrac{{1}}{{t_{j}}}[\nicefrac{{\partial g}}{{\partial w_{j}}}(w^{*})]_{i}>0$ is less than the $q_{j}$ largest positive elements of $w_{j}^{*}$, which means that the subvector $w_{j}^{*}$ has reached its maximum number of positive entries. In other words, all elements of $\nabla g(w^{*})$ are zero, except for possibly in entries belonging to subvectors $w_{j}^{*}$ that have reached their maximum number of positive elements; in this case, we will have $w_{i}^{*}=0$ and $[\nabla g(w^{*})]_{i}<0$.

From Theorem 1, we know that all convergent subsequences of $\left\{w^{(k)}\right\}_{k\in\mathbb{N}}$ converge to vectors $w^{*}$ such that $g\left(w^{*}\right)=g^{*}$. Consider the convergent subsequence $w^{(p_{k})}\to w^{*}$ as $k\to\infty$. For all $\delta>0$, there exists some positive integer $L$ such that for all $k\geq L$, we have

[TABLE]

Let $\omega$ denote the smallest positive element of $w^{*}$, and let $\delta<\nicefrac{{\omega}}{{2}}$. Then for all $k\geq L$, the iterates $w^{(p_{k})}$ must be positive everywhere that $w^{*}$ is positive, including all subvectors $w_{j}^{*}$ that have reached their maximum number of positive elements. In particular, this means that $w^{(p_{k})}$ cannot be positive where $[\nabla g(w^{*})]_{i}<0$, which implies that $\nabla g(w^{*})^{T}w^{(p_{k})}\geq 0$.

Suppose that for some $k\geq L$ and $n>0$, we also have

[TABLE]

Then the iterate $w^{(p_{k}+n)}$ is positive everywhere that $w^{*}$ is positive, which means that $\nabla g(w^{*})^{T}w^{(p_{k}+n)}\geq 0$. If we can guarantee that the next iterate $w^{(p_{k}+n+1)}$ also satisfies $\nabla g(w^{*})^{T}w^{(p_{k}+n+1)}\geq 0$, and we let both $\sqrt{\nicefrac{{2\epsilon}}{{\mu}}}<\delta$ and $p_{k}+n+1\geq K$, then from inequality $\eqref{eq:wDIffBound}$ we will have

[TABLE]

The only way that we could have $\nabla g(w^{*})^{T}w^{(p_{k}+n+1)}<0$ would be if there were some element $i$ such that $w_{i}^{(p_{k}+n+1)}>0$ where $w_{i}^{*}=0$ and $[\nabla g(w^{*})]_{i}<0$, which would need to occur on a subvector $w_{j}^{*}$ that has reached its maximum number of positive elements. Because $w^{(p_{k}+n)}$ is positive everywhere that $w^{*}$ is positive, the subvector $w_{j}^{(p_{k}+n)}$ would also have reached its maximum number of positive elements. In order for $w_{i}^{(p_{k}+n+1)}>0$ where $w_{i}^{*}=0$ and $w_{i}^{(p_{k}+n)}\leq 0$, we would require $w_{\ell}^{(p_{k}+n+1)}\leq 0$ in an entry $\ell$ where $w_{\ell}^{(p_{k}+n)}>0$, that is, we would need to swap the location of a positive element in the subvector. However, since we assumed that

[TABLE]

this would require

[TABLE]

Now suppose $t_{j}<\nicefrac{{n_{j}}}{{\alpha_{j}}}$ for all $j\in\mathcal{J}^{\text{dv}}$, and define $\beta:=\min_{j\in\mathcal{J}^{\text{dv}}}(\nicefrac{{1}}{{2t_{j}}}-\nicefrac{{\alpha_{j}}}{{2n_{j}}})>0$. Then from inequality (74), we have for all $k\geq K$

[TABLE]

which means that

[TABLE]

In particular, if we add the assumption that $\sqrt{\nicefrac{{\epsilon}}{{\beta}}}<\nicefrac{{\omega}}{{2}}$, then

[TABLE]

and it is impossible to swap the location of the positive elements in the iterates $w^{(p_{k}+n)}$ and $w^{(p_{k}+n+1)}$. Therefore $\nabla g(w^{*})^{T}w^{(p_{k}+n+1)}\geq 0$, and we have

[TABLE]

By induction, if we choose $\delta<\nicefrac{{\omega}}{{2}}$ and $\epsilon<\min\left\{\nicefrac{{\mu\delta^{2}}}{{2}},\nicefrac{{\beta\omega^{2}}}{{4}}\right\}$, then for all integers $k\geq\max\left\{p_{L},K\right\}$, we have

[TABLE]

and the sequence of iterates $\left\{w^{(k)}\right\}_{k\in\mathbb{N}}$ converges to a fixed point $w^{*}$ of Algorithm 1 such that $g(w^{*})=g^{*}$, where $g^{*}$ is a local minimum of $g$ on $\Omega$.

To establish the convergence rate of Algorithm 1, we rearrange the terms in (86) and sum over $k$ to get

[TABLE]

Define $\zeta:=\min_{j\in\mathcal{J}^{\text{dv}}}\left(t_{j}\right)>0$. For any $\epsilon>0$, we let

[TABLE]

and solving for $k$, we find that we reach an $\epsilon$-accurate solution after $\mathcal{O}\left(\nicefrac{{1}}{{\epsilon}}\right)$ iterations. ∎

Appendix C Figures for Section 4.4