You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
#### Note: The package is at its early stages of development. At this point, we're only including our recent work [1] where we used wavelets to induce fluence smoothness. This work was originally written in Matlab.
Compress radiation therapy planning (CompressRTP) is a python package that aims to leverage and re-purpose the dimensionality reduction tools to improve the speed and quality of the radiation treatment planning. The dimensionality reduction has a rich history in statistics and has gained attention recently and re-emerged as a powerful tool to deal with the increasingly dimensional problems arising in the fields of big-data and machine learning. They also have made profound impacts in imaging science where the high-dimensional images are often lying in a low-dimensional subspace.
10
8
11
-
The optimization problems arising in radiation therapy also suffer from the curse of dimensionality with so many radiotherapy machine parameters to be optimized (e.g., beamlets, MLC leaf positions, beam angles) and the discretization of the patient body to many small 3-dimensional cubes (known as voxels). This is an on-going project and at this point we only included the work based on our recent publication [1] which uses low-frequency wavelets to represent and characterize the beamlet intensities in order to induce fluence smoothness in intensity modulated radiation therapy (IMRT).
9
+
Radiotherapy is used to treat more than half of all cancer patients, either on its own or alongside other treatments like surgery, chemotherapy, or immunotherapy. It works by directing high-energy radiation beams at the patient's body to destroy cancer cells. Because every patient's anatomy is unique, radiotherapy must be personalized. This means customizing the radiation beams to effectively target the tumor while minimizing harm to nearby healthy tissue.
12
10
13
-
The key dimensionality reduction idea in the wavelet-induced smoothness project is that the beamlet intensity maps are lying in a low-dimensional subspace of all possible intensity maps since they need to be smooth (to reduce the plan complexity, an unnecessary modulation, and improve the plan delivery efficiency). We propose representing the beamlet intensities using an incomplete wavelet basis of low-frequency wavelets that explicitly excludes fluctuating intensity maps from the decision space (explicit hard constraint). This technique provides a built-in wavelet-induced smoothness and excludes complex and clinically irrelevant radiation plans from the search space that improves both dosimetric plan quality and delivery efficiency.
11
+
# The Challenge
14
12
15
-
From the implementation perspective, the current tool is an add-on for the PortPy project. However, the proposed technique can be easily integrated with any optimization approach by simply adding a set of linear constraints in the form of (_x=Wy_), where _x_ represents the beamlet intensity, _W_ is the matrix including low-frequency wavelets (can be generated using method `get_low_dim_basis` in the code), and _y_ is a free variable.
13
+
Personalizing radiotherapy involves solving large and complex optimization problems. These problems need to be solved quickly due to the limited time available in clinical settings. Currently, they are often solved using rough approximations, which can lead to less effective treatments. This might result in the tumor not receiving enough radiation or healthy tissues being exposed to too much.
14
+
15
+
# The CompressRTP Project
16
+
17
+
The CompressRTP project aims to solve these optimization problems both rapidly and accurately. This ongoing project currently includes tools introduced in our latest three publications [1, 2, 3]. These tools are available as three Jupyter Notebooks in this repository:
The optimization problems in radiotherapy are highly complex due to the "curse of dimensionality" since they involve many beams, beamlets (small segments of beams), and voxels (3D pixels representing volume). However, much of this data is redundant because it comes from discretizing a system that is inherently continuous. For example, radiation doses delivered from adjacent beamlets are highly correlated, and radiation doses delivered to neighboring voxels are very similar. This redundancy means that large-scale radiotherapy optimization problems are highly compressible, which is the foundation of **CompressRTP**.
27
+
28
+
Dimensionality reduction and compression have a rich history in statistics and engineering. Recently, these techniques have re-emerged as powerful tools for addressing increasingly high-dimensional problems in fields like big data and machine learning. Our goal is to adapt and adopt these versatile methods to embed high-dimensional radiotherapy optimization problems into lower-dimensional spaces so they can be solved efficiently. A general radiotherapy optimization problem can be formulated as:
29
+
30
+
$Minimize f(Ax,x)$
31
+
32
+
Subject to $g(Ax,x)\leq 0,x\geq 0$
33
+
34
+
35
+
**CompressRTP** currently addresses the following two issues with this problem:
36
+
37
+
1.**Matrix Compression to Address the Computational Intractability of $A$:**
38
+
39
+
-**The Challenge:** The matrix $𝐴$ is large and dense (approximately 100,000–500,000 rows and 5,000–20,000 columns) and is the main source of computational difficulty in solving radiotherapy optimization problems.
40
+
-**Traditional Approach:** This matrix is often sparsified in practice by simply ignoring small elements (e.g., zeroing out elements less than 1% of the maximum value in $𝐴$, which can potentially lead to sub-optimal treatment plans.
41
+
-**CompressRTP Solutions:** We provide a compressed and accurate representation of matrix $𝐴$ using two different techniques:
42
+
-**(1.1) Sparse-Only Compression:** This technique sparsifies $𝐴$ using advanced tools from probability and randomized linear algebra.
43
+
-**(1.2) Sparse-Plus-Low-Rank Compression:** This method decomposes $𝐴$ into a sum of a sparse matrix and a low-rank matrix.
44
+
2.**Fluence Compression to Enforce Smoothness on $𝑥$:**
45
+
-**The Need for Smoothness:** The beamlet intensities $𝑥$ need to be smooth for efficient and accurate delivery of radiation. Smoothness refers to small variations in the intensity of neighboring beamlets in two dimensions.
46
+
-**Traditional Approach:** Smoothness is often achieved implicitly by adding regularization terms to the objective function that discourage variations between neighboring beamlets.
47
+
-**CompressRTP Solution:** We enforce smoothness explicitly by representing the beamlet intensities using low-frequency wavelets, resulting in built-in wavelet-induced smoothness. This can be easily integrated into the optimization problem by adding a set of linear constraints.
48
+
49
+
50
+
# 1) Matrix Compression to Address the Computational Intractability of $𝐴$
51
+
52
+
## 1.1) Sparse-Only Matrix Compression
53
+
Matrix sparsification has been extensively studied in the machine learning community for applications such as low-rank approximation and Principal Component Analysis (PCA). This technique is also a key part of an emerging field known as randomized linear algebra. The main idea is to carefully sample and scale elements from the original dense matrix $𝐴$ to create a sparse "sketch" matrix $𝑆$ that closely resembles the characteristics of $𝐴$ (for example, ensuring that
54
+
$||𝐴-S||_2$ is small).
55
+
56
+
In radiotherapy optimization, we can replace the original dense matrix $𝐴$ with this sparse matrix $𝑆$ and solve the following surrogate optimization problem:
57
+
58
+
$Minimize f(Sx,x)$
59
+
60
+
Subject to $g(Sx,x)\leq 0,x\geq 0$
61
+
$(S≈A,S$ is sparse, $A$ is dense)
62
+
63
+
64
+
We have developed a simple yet effective matrix sparsification technique with desirable mathematical properties (see **Algorithm 1** below). The main idea is to **retain the large elements of the matrix deterministically** and **handle the smaller elements probabilistically** (see the [paper](./images/RMR_NeurIPS_Paper.pdf) for a detailed explanation).
Our **Randomized Minor Rectification (RMR)** algorithm transforms a dense matrix $𝐴$ into a sparse matrix $𝑆$ that typically contains only 2–4% of the non-zero elements of $𝐴$, while maintaining negligible accuracy loss (i.e., small feasibility and sub-optimality gaps).
71
+
72
+
This approach ensures that an optimal solution of the surrogate problem is a near-optimal solution of the original problem, as stated in the following two theorems:
73
+
74
+
**Small Feasibility Gap (Theorem 3.6 in the Paper)**:
75
+
An optimal point of the surrogate problem violates each constraint of the original problem by no more than $(19 + 5 log m)\epsilon ∥x∥_2$ with a probability of at least 95%.
76
+
77
+
**Small Sub-Optimality Gap (Theorem 3.9 in the Paper)**: An optimal point of the surrogate problem, $x_A$, is a near-optimal solution to the original problem with a probability of at least 95%, and the sub-optimality gap of O(e), where $e = (19 + 5 logm) \epsilon max (∥x_A∥_2, ∥x_S∥_2)$ ($x_A$ is an optimal solution of the original problem).
**Figure Explanation:** The figure above compares the proposed RMR algorithm with four existing sparsification algorithms in terms of feasibility and optimality gaps. These gaps were calculated by solving both the original and surrogate optimization problems for 10 lung cancer patients, whose data is publicly available on PortPy. The results demonstrate that the RMR algorithm outperforms the existing methods.
**Figure Explanation:** The figure above illustrates the discrepancies in Dose Volume Histogram (DVH) plots between the actual dose ($𝐴𝑥$, shown as a solid line) and the approximated dose ($𝑆𝑥$, shown as a dotted line), where
90
+
$𝑥$ is the optimal solution of the surrogate optimization problem. A smaller gap between the dotted and solid lines indicates a more accurate dose approximation. **Left figure** demonstrates a significant dose discrepancy when the matrix
91
+
$𝐴$ is sparsified by simply zeroing out small elements—a technique commonly used in practice. **Right figure** shows a minimal dose discrepancy when the matrix $𝐴$ is sparsified using the RMR algorithm. Importantly, in both cases, the sparsified matrix contained only 2% non-zero elements.
92
+
93
+
94
+
**Implementation in PortPy:**
95
+
96
+
If you are using PortPy for your radiotherapy research, you can apply RMR sparsification by simply adding the following lines of code:
97
+
98
+
```python
99
+
from compress_rtp.utils.get_sparse_only import get_sparse_only
100
+
101
+
# Apply RMR sparsification to the matrix A
102
+
S = get_sparse_only(matrix=A, threshold_perc=10, compression='rmr')
103
+
104
+
# Replace the original matrix A with the sparsified matrix S
105
+
inf_matrix.A = S
106
+
```
107
+
108
+
## 1.2) Sparse-Plus-Low-Rank Matrix Compression
109
+
The rows of matrix $𝐴$ correspond to the patient's voxels, and the similarity of doses delivered to neighboring voxels makes these rows highly correlated (the same argument applies to the columns). This high correlation means that matrix $𝐴$
**Figure Explanation:** The low-rank nature of matrix $𝐴$ can be verified by observing the exponential decay of its singular values, as shown by the blue line in the left figure. If we decompose matrix
117
+
$𝐴$ into $𝐴=𝑆+𝐿$, where $𝑆$ is a sparse matrix containing large-magnitude elements (e.g., elements greater than 1% of the maximum value of $𝐴$), and $𝐿$ includes smaller elements mainly representing scattering doses, then the singular values of the scattering matrix $𝐿$ reveal an even sharper exponential decay (depicted by the red line). This suggests the use of “sparse-plus-low-rank” compression, $𝐴≈𝑆+𝐻𝑊$, as schematically shown in the right figure.
118
+
119
+
120
+
The matrix $𝑆$ is sparse, $𝐻$ is a “tall skinny matrix” with only a few columns, and $𝑊$ is a “wide short matrix” with only a few rows. Therefore, $𝐴≈𝑆+𝐻𝑊$ provides a compressed representation of the data. This allows us to solve the following surrogate problem instead of the original problem
**Figure Explanation:** The figure above demonstrates improvements in plan quality across multiple clinical criteria when using the sparse-plus-low-rank compression method instead of naively sparsifying the matrix
131
+
$𝐴$ by simply removing small elements. This comparison was performed on data from 10 lung cancer patients. Additionally, using compression reduced the average optimization time by approximately 20% compared to the naïve sparsification. More aggressive compression can be applied to gain even greater acceleration.
132
+
133
+
Decomposing a matrix into the sum of a sparse matrix and a low-rank matrix has found numerous applications in fields such as computer vision, medical imaging, and statistics. Historically, this structure has been employed as a form of prior knowledge to recover objects of interest that manifest themselves in either the sparse or low-rank components.
134
+
135
+
However, the application presented here represents a novel departure from conventional uses of sparse-plus-low-rank decomposition. Unlike traditional settings where specific components (sparse or low-rank) hold intrinsic importance, our primary goal is not to isolate or interpret these structures. Instead, we leverage them for computationally efficient matrix representation. In this case, the structure serves purely as a tool for optimizing computational efficiency while maintaining data integrity.
136
+
137
+
**Note:** Both sparse-only and sparse-plus-low-rank compression techniques serve the same purpose. We are currently investigating the strengths and weaknesses of each technique and their potential combination. Stay tuned for more results.
138
+
139
+
**Implementation in PortPy:**
140
+
141
+
In PortPy, you can apply the sparse-plus-low-rank compression using the following lines of code. Unlike the sparse-only compression using RMR, which did not require any changes other than replacing $𝐴x$ with $𝑆x$ in your optimization formulation and code, this compression requires adding a linear constraint $x=𝑊𝑦$ and replacing $Ax$ with $𝑆x+Hy$. These changes can be easily implemented using CVXPy (see the [Sparse-Plus-Low-Rank Matrix Compression](https://github.com/PortPy-Project/CompressRTP/blob/main/examples/matrix_sparse_plus_low_rank.ipynb) for details).
142
+
143
+
```python
144
+
from compress_rtp.utils.get_sparse_plus_low_rank import get_sparse_plus_low_rank
145
+
146
+
S, H, W = get_sparse_plus_low_rank(A=A, threshold_perc=1, rank=5)
147
+
```
148
+
149
+
150
+
## 2) Fluence Compression to Enforce Smoothness on $𝑥$
151
+
152
+
The fluence smoothness required for efficient and accurate plan delivery is typically achieved by adding an additional "regularization" term to the objective function. This term measures local variations in adjacent beamlets to discourage fluctuating beamlet intensities. However, a significant limitation of this method is its focus on **local complexity** within each beam—it assesses variations between adjacent beamlets but overlooks the **global complexity** of the entire plan. Another challenge is that achieving an optimal balance between plan complexity and dosimetric quality requires careful fine-tuning of the importance weight associated with the smoothness term in the objective function.
153
+
154
+
To address these challenges, we treat the intensity map of each beam as a **2D image** and represent it using wavelets corresponding to **low-frequency changes**. The compressed representation of fluence using low-frequency wavelets induces built-in local and global smoothness that can be achieved **without any hyperparameter fine-tuning**. This approach can be easily incorporated into the optimization by adding a set of linear constraints in the form of $𝑥=𝑊𝑦$, where $W$ is the matrix including low-frequency wavelets.
**Figure Explanation:** As illustrated in the figure above, the treatment plan achieved using wavelet compression is not only more conformal to the tumor but also less complex. This is evidenced by a smaller duty cycle compared to the plan achieved by adding only a regularization term to the objective function.
162
+
163
+
164
+
**Implementation in PortPy:**
165
+
166
+
In **PortPy**, you can incorporate wavelet smoothness by adding the following lines of code:
167
+
168
+
```python
169
+
from compress_rtp.utils.get_low_dim_basis import get_low_dim_basis
170
+
171
+
# Generate the matrix W of low-frequency wavelets
172
+
W = get_low_dim_basis(inf_matrix=inf_matrix, compression='wavelet')
173
+
174
+
# Add the variable y
175
+
y = cp.Variable(W.shape[1])
176
+
177
+
# Add the constraint Wy = x
178
+
opt.constraints += [W @ y == opt.vars['x']]
179
+
```
16
180
17
181
## License
18
182
CompressRTP code is distributed under **Apache License 2.0 with Commons Clause**, and is available for non-commercial academic purposes.
19
183
20
-
## Team
21
-
1.[Mojtaba Tefagh](https://mtefagh.github.io/) ([Sharif University of Technology](https://en.sharif.edu/))
22
-
2.[Masoud Zarepisheh](https://masoudzp.github.io/) ([Memorial Sloan Kettering Cancer Center](https://www.mskcc.org/))
23
-
3.[Gourav Jhanwar](https://github.com/gourav3017) ([Memorial Sloan Kettering Cancer Center](https://www.mskcc.org/))
24
-
25
184
## Reference
26
185
If you find our work useful in your research or if you use parts of this code, please cite the following paper:
27
186
```
187
+
188
+
@article{Adeli2024Randomized,
189
+
title={Randomized Sparse Matrix Compression for
190
+
Large-Scale Constrained Optimization in Cancer
191
+
Radiotherapy},
192
+
author={Adeli, Shima and Tefagh, Mojtaba and Jhanwar, Gourav and Zarepisheh, Masoud},
193
+
journal={accepted at NeurIPS},
194
+
year={2024}
195
+
}
196
+
197
+
28
198
@article{tefagh2023built,
29
199
title={Built-in wavelet-induced smoothness to reduce plan complexity in intensity modulated radiation therapy (IMRT)},
30
200
author={Tefagh, Mojtaba and Zarepisheh, Masoud},
@@ -35,5 +205,14 @@ If you find our work useful in your research or if you use parts of this code, p
35
205
year={2023},
36
206
publisher={IOP Publishing}
37
207
}
208
+
209
+
@article{tefagh2024compressed,
210
+
title={Compressed radiotherapy treatment planning (CompressRTP): A new paradigm for rapid and high-quality treatment planning optimization},
211
+
author={Tefagh, Mojtaba and Jhanwar, Gourav and Zarepisheh, Masoud},
0 commit comments