Reference MATLAB implementation of the R-Trust / R3-Trust algorithms from:

M. Tefagh, G. Jhanwar, and M. Zarepisheh, "Sparse plus low-rank matrix embedding with applications in cancer radiotherapy optimization."

License: Apache License 2.0 with the Commons Clause — free for non-commercial, academic use (see LICENSE).

---

Overview

Given a large dense matrix A (m × n), SLME builds a computationally efficient surrogate

A ≈ S + H*W

where S is sparse, and H (m × r) and W (r × n) form a low-rank factor with r ≪ min(m, n). The point of the surrogate is speed: a matrix–vector product

A*x  ≈  S*x + H*(W*x)

costs nnz(S) + r*(m+n) operations instead of the m*n operations of a dense A*x.

Unlike Robust PCA and related recovery models, SLME does not try to recover interpretable latent components. It seeks the best accuracy-vs-cost trade-off, formulated as a bi-objective, nonconvex problem that trades approximation error ‖A − (S + H*W)‖_F² against representation cost η·(m+n)·rank(L) + ‖S‖₀. The proposed trust-region algorithms approximate the entire Pareto frontier of this trade-off in a single, parameter-free run.

The motivating application is cancer radiotherapy treatment planning, where a large dense dose-influence matrix is the main computational bottleneck. SLME also works well as a matrix-recovery method on synthetic sparse-plus-low-rank data (see SLME_Synthetic.m).

---

The algorithms

Algorithm	File	Description
R-Trust	R_Trust.m	Exact recursive trust-region baseline. Each iteration compares one exact rank-1 SVD update against the exact sparse update and keeps whichever removes more residual energy.
R2-Trust	—	General framework (in the paper) that allows approximate sparse/low-rank projection subroutines with a rectification step.
R3-Trust	R3_Trust.m	Fast, scalable instantiation of R2-Trust using randomized SVD and a sampled quantile threshold, plus batched low-rank updates. Recommended for large matrices.

Both R_Trust and R3_Trust are parameter-free: a single run returns the whole error-vs-cost surface. R3_Trust's batch_size (default 10) trades speed for granularity, like a step size — it does not need tuning to produce the surface.

---

Requirements

• MATLAB (developed and tested on R2023a). No additional toolboxes are required for the core R-Trust / R3-Trust algorithms.
• An internet connection only if you want to reproduce the radiotherapy experiment, which downloads patient data from the PortPy dataset on Hugging Face.

---

Quick start

Decompose any dense matrix:

% Build a sparse-plus-low-rank test matrix  A = S + H*W
m = 2000; n = 500; r = 20;                       % dimensions and rank of the low-rank part
Htrue = randn(m, r);  Wtrue = randn(r, n);       % random low-rank factors (Htrue*Wtrue is rank r)
density = 0.05;                                  % fraction of nonzero entries in the sparse part
Strue = (rand(m, n) < density) .* (20*rand(m, n) - 10);   % random sparse matrix, entries in [-10, 10]
A = Strue + Htrue*Wtrue;                          % the dense matrix to embed

out = R3_Trust(A);               % single, parameter-free run

% Final embedding  A ≈ H*W + S
% (the recovered factors need not match Strue/Htrue/Wtrue — SLME targets the
%  accuracy-vs-cost trade-off, not recovery of the underlying components)
H = out(end).H;                  % m × r  (low-rank left factor)
W = out(end).W;                  % r × n  (low-rank right factor, orthonormal rows)
S = out(end).S;                  % m × n  sparse matrix

% Fast matrix-vector product
x  = randn(size(A,2), 1);
Ax = S*x + H*(W*x);              % ~ nnz(S) + r*(m+n) flops, vs m*n for A*x

Each entry of the returned struct array out is one point on the accuracy-vs-cost surface:

Field	Meaning
out(k).relative_error	‖A − (H*W + S)‖_F / ‖A‖_F at iteration k
out(k).relative_NNZ	(nnz(S) + (m+n)*rank) / (m*n) — the relative storage/flop cost
out(k).time	cumulative wall-clock time
out(end).H/.W/.S	the finalized factors

Common options (pass as a struct):

opts = struct('relative_error_bound', 0.01, ...  % stop at 1% relative error
              'batch_size',           10,   ...  % R3-Trust only
              'verbose',              1);
out  = R3_Trust(A, opts);

Plot the full trade-off surface:

plot_pareto({out}, {'R3-Trust'}, 'YScale', 'log');

Use R_Trust in place of R3_Trust when you want the exact reference behavior on smaller matrices.

---

Reproducing the paper experiments

1. Radiotherapy (matrix embedding vs. matrix recovery)

SLME_Radiotherapy

At the top of the script, set the patient_id and beam_ids. If the corresponding data/<patient_id>.mat file is not present, the script calls Radiotherapy_Data_Gen.m to download the beams from the PortPy Hugging Face dataset and assemble the dose-influence matrix automatically (this can take a while and a fair amount of disk on the first run).

As of August 2026 the PortPy dataset provides roughly 200 lung cases (Lung_Patient_1 … Lung_Patient_202) and about 100 prostate cases (Prostate_Patient_1 … Prostate_Patient_99), plus a lung phantom (Lung_Phantom_Patient_1); see the dataset page for the current list.

The script compares R-Trust / R3-Trust against the recovery baselines (GoDec, VB, ADMM) and against sparse-only and low-rank-only approximations, and displays the Pareto plots.

2. Synthetic sparse-plus-low-rank recovery

SLME_Synthetic

Generates synthetic A = L + S + E instances following the protocol of Zhou et al. (2014) and compares R3-Trust with GoDec, VB, and ADMM across sparsity and noise regimes.

---

Repository contents

File	Purpose
R_Trust.m	Exact R-Trust algorithm
R3_Trust.m	Randomized R3-Trust algorithm (recommended)
SLME_Radiotherapy.m	Radiotherapy experiment driver
SLME_Synthetic.m	Synthetic-recovery experiment driver
Radiotherapy_Data_Gen.m	Builds a dose-influence matrix from the PortPy dataset
plot_pareto.m	Accuracy-vs-cost / accuracy-vs-time plotting
sparse_only_pareto.m, lowrank_only_pareto.m	Sparse-only and low-rank-only baselines
GoDec.m	GoDec baseline (Zhou & Tao, ICML 2011)
VBRPCA.m	Variational Bayesian RPCA baseline (Babacan et al., 2012)
ADMM_PCP.m, ADMM_SPCP.m	ADMM baselines (based on Boyd et al.)
denoise_decomposition.m	Numerical-rank / round-off cleanup helper
LICENSE	Apache 2.0 with Commons Clause

The baseline solvers (GoDec, VBRPCA, the ADMM_* files) are adapted from their original authors' implementations and retain their own attributions; they are included so the paper's comparisons can be reproduced.

---

Citation

If you use this code, please cite:

@article{tefagh2026slme,
  title   = {Sparse plus low-rank matrix embedding with applications in cancer radiotherapy optimization},
  author  = {Tefagh, Mojtaba and Jhanwar, Gourav and Zarepisheh, Masoud},
  year    = {2026},
  note    = {Preprint; bibliographic details will be updated upon publication.}
}

---

References

The baseline solvers and datasets included here are due to their original authors:

• GoDec — T. Zhou and D. Tao, "GoDec: Randomized low-rank & sparse matrix decomposition in noisy case," in Proc. 28th International Conference on Machine Learning (ICML), 2011.
• Variational Bayesian RPCA (VB) — S. D. Babacan, M. Luessi, R. Molina, and A. K. Katsaggelos, "Sparse Bayesian methods for low-rank matrix estimation," IEEE Transactions on Signal Processing, vol. 60, no. 8, pp. 3964–3977, 2012.
• ADMM — S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, "Distributed optimization and statistical learning via the alternating direction method of multipliers," Foundations and Trends in Machine Learning, vol. 3, no. 1, pp. 1–122, 2011. The matrix-decomposition solver is adapted from N. Parikh and S. Boyd, "Proximal algorithms," Foundations and Trends in Optimization, vol. 1, no. 3, pp. 127–239, 2014.
• Stable Principal Component Pursuit (the model solved by ADMM_SPCP.m) — Z. Zhou, X. Li, J. Wright, E. Candès, and Y. Ma, "Stable principal component pursuit," in Proc. IEEE International Symposium on Information Theory (ISIT), 2010.
• Robust PCA / Principal Component Pursuit — E. J. Candès, X. Li, Y. Ma, and J. Wright, "Robust principal component analysis?," Journal of the ACM, vol. 58, no. 3, pp. 1–37, 2011.
• Synthetic-data protocol — X. Zhou, C. Yang, H. Zhao, and W. Yu, "Low-rank modeling and its applications in image analysis," ACM Computing Surveys, vol. 47, no. 2, article 36, 2014.
• PortPy dataset — PortPy: an open-source platform for cancer radiotherapy treatment planning. Project: https://github.com/PortPy-Project/PortPy; dataset: https://huggingface.co/datasets/PortPy-Project/PortPy_Dataset.

---

License

Licensed under the Apache License 2.0 with the Commons Clause — you may use, modify, and redistribute this software for non-commercial, academic purposes, but you may not sell it. See LICENSE for the full terms.

---

Authors

• Mojtaba Tefagh — University of Edinburgh (Email: m.tefagh@ed.ac.uk)
• Gourav Jhanwar — Memorial Sloan Kettering Cancer Center (Email: JhanwarG@mskcc.org)
• Masoud Zarepisheh — Memorial Sloan Kettering Cancer Center (Email: zarepism@mskcc.org)