problemequation problem

#  Benchopt: Reproducible, efficient and collaborative
optimization benchmarks

Thomas Moreau  Université Paris-Saclay, Inria, CEA, 91120 Palaiseau, France
corresponding authors<EMAIL_ADDRESS>Mathurin Massias  Univ Lyon, Inria,
CNRS, ENS de Lyon, UCB Lyon 1, LIP UMR 5668, F-69342, Lyon, France
corresponding authors<EMAIL_ADDRESS>Alexandre Gramfort  Université
Paris-Saclay, Inria, CEA, 91120 Palaiseau, France  corresponding authors,
<EMAIL_ADDRESS>Pierre Ablin  Université Paris-Dauphine, PSL University,
CNRS, 75016, Paris, France  Pierre-Antoine Bannier  Benjamin Charlier  IMAG,
Univ Montpellier, CNRS, Montpellier, France  Mathieu Dagréou  Université
Paris-Saclay, Inria, CEA, 91120 Palaiseau, France  Tom Dupré la Tour
University of California, Berkeley, CA 94720, USA  Ghislain Durif  IMAG, Univ
Montpellier, CNRS, Montpellier, France  Cassio F. Dantas  IMAG, Univ
Montpellier, CNRS, Montpellier, France  Quentin Klopfenstein  University of
Luxembourg, LCSB, Esch-sur-Alzette, Luxembourg  Johan Larsson  The Department
of Statistics, Lund University  En Lai  Université Paris-Saclay, Inria, CEA,
91120 Palaiseau, France  Tanguy Lefort  IMAG, Univ Montpellier, CNRS,
Montpellier, France  Benoit Malézieux  Université Paris-Saclay, Inria, CEA,
91120 Palaiseau, France  Badr Moufad  Univ Lyon, Inria, CNRS, ENS de Lyon, UCB
Lyon 1, LIP UMR 5668, F-69342, Lyon, France  Binh T. Nguyen  LTCI, Télécom
Paris, 91120 Palaiseau, France  Alain Rakotomamonjy  Criteo AI Lab, Paris,
France  Zaccharie Ramzi  ENS Ulm, CNRS, UMR 8553, Paris, France  Joseph Salmon
IMAG, Univ Montpellier, CNRS, Montpellier, France  Institut Universitaire de
France (IUF)  Samuel Vaiter  CNRS & Université Côte d’Azur, Laboratoire J.A.
Dieudonné, CNRS, Nice, France

###### Abstract

Numerical validation is at the core of machine learning research as it allows
to assess the actual impact of new methods, and to confirm the agreement
between theory and practice. Yet, the rapid development of the field poses
several challenges: researchers are confronted with a profusion of methods to
compare, limited transparency and consensus on best practices, as well as
tedious re-implementation work. As a result, validation is often very partial,
which can lead to wrong conclusions that slow down the progress of research.
We propose Benchopt, a collaborative framework to automate, reproduce and
publish optimization benchmarks in machine learning across programming
languages and hardware architectures. Benchopt simplifies benchmarking for the
community by providing an off-the-shelf tool for running, sharing and
extending experiments. To demonstrate its broad usability, we showcase
benchmarks on three standard learning tasks: $\ell_{2}$-regularized logistic
regression, Lasso, and ResNet18 training for image classification. These
benchmarks highlight key practical findings that give a more nuanced view of
the state-of-the-art for these problems, showing that for practical
evaluation, the devil is in the details. We hope that Benchopt will foster
collaborative work in the community hence improving the reproducibility of
research findings.

##  1 Introduction

Numerical experiments have become an essential part of statistics and machine
learning (ML). It is now commonly accepted that every new method needs to be
validated through comparisons with existing approaches on standard problems.
Such validation provides insight into the method’s benefits and limitations
and thus adds depth to the results. While research aims at advancing knowledge
and not just improving the state of the art, experiments ensure that results
are reliable and support theoretical claims [120]. Practical validation also
helps the ever-increasing number of ML users in applied sciences to choose the
right method for their task. Performing rigorous and extensive experiments is,
however, time-consuming [111], particularly because comparisons against
existing methods in new settings often requires reimplementing baseline
methods from the literature. In addition, ingredients necessary for a proper
reimplementation may be missing, such as important algorithmic details,
hyperparameter choices, and preprocessing steps [105].

In the past years, the ML community has actively sought to overcome this
“reproducibility crisis” [66] through collaborative initiatives such as open
datasets (OpenML, [130]), standardized code sharing [52], benchmarks (MLPerf,
[93]), the NeurIPS and ICLR reproducibility challenges [105, 104] and new
journals (e.g.,  [115]). As useful as these endeavors may be, they do not,
however, fully address the problems in optimization for ML since, in this
area, there are no clear community guidelines on how to perform, share, and
publish benchmarks.

Optimization algorithms pervade almost every area of ML, from empirical risk
minimization, variational inference to reinforcement learning [123]. It is
thus crucial to know which methods to use depending on the task and setting
[6]. While some papers in optimization for ML provide extensive validations
[89], many others fall short in this regard due to lack of time and resources,
and in turn feature results that are hard to reproduce by other researchers.
In addition, both performance and hardware evolve over time, which eventually
makes static benchmarks obsolete. An illustration of this is the recent work
by [119], which extensively evaluates the performances of 15 optimizers across
8 deep-learning tasks. While their benchmark gives an overall assessment of
the considered solvers, this assessment is bound to become out-of-date if it
is not updated with new solvers and new architectures. Moreover, the benchmark
does not reproduce state-of-the-art results on the different datasets,
potentially indicating that the considered architectures and optimizers could
be improved.

We firmly believe that this critical task of maintaining an up-to-date
benchmark in a field cannot be solved without a collective effort. We want to
empower the community to take up this challenge and build a living,
reproducible and standardized state of the art that can serve as a foundation
for future research.

Figure 1: A visual summary of Benchopt. Each Solver is run (in parallel) on
each Dataset and each variant of the Objective. Results are exported as a CSV
file that is easily shared and can be automatically plotted as interactive
html visualizations or PDF figures.

Benchopt provides the tools to structure the optimization for machine learning
(Opt-ML) community around standardized benchmarks, and to aggregate individual
efforts for reproducibility and results sharing. Benchopt can handle
algorithms written in Python, R, Julia or C/C++ via binaries. It offers built-
in functionalities to ease the execution of benchmarks: parallel runs,
caching, and automatical results archiving. Benchmarks are meant to evolve
over time, which is why Benchopt offers a modular structure through which a
benchmark can be easily extended with new objective functions, datasets, and
solvers by the addition of a single file of code.

The paper is organized as follows. We first detail the design and usage of
Benchopt, before presenting results on three canonical problems:

  * •

$\ell_{2}$-regularized logistic regression: a convex and smooth problem which
is central to the evaluation of many algorithms in the Opt-ML community, and
remains of high relevance for practitioners;

  * •

the Lasso: the prototypical example of non-smooth convex problem in ML;

  * •

training of ResNet18 architecture for image classification: a large scale non-
convex deep learning problem central in the field of computer vision.

The reported benchmarks, involving dozens of implementations and datasets,
shed light on the current state-of-the-art solvers for each problem, across
various settings, highlighting that the best algorithm largely depends on the
dataset properties (e.g., size, sparsity), the hyperparameters, as well as
hardware. A variety of other benchmarks (e.g., MCP, TV1D, etc.) are also
presented in Appendix, with the goal to facilitate contributions from the
community.

By the open source and collaborative design of Benchopt (BSD 3-clause
license), we aim to open the way towards community-endorsed and peer-reviewed
benchmarks that will improve the tracking of progress in optimization for ML.

##  2 The Benchopt library

The Benchopt library aims to provide a test bed to benchmark algorithms for
optimization in ML, where the problems depend on some input dataset
$\mathcal{D}$. The considered problems are of the form

$\theta^{*}\in\mathop{\mathrm{arg\,min}}_{\theta\in\Theta}f(\theta;\mathcal{D},\Lambda)\enspace,$
(1)

where $f$ is the objective function, $\Lambda$ are its hyperparameters, and
$\Theta$ is the feasible set for $\theta$. The scope of the library is to
evaluate optimization methods in their wide sense by considering the sequence
$\\{\theta^{t}\\}_{t}$ produced to approximate $\theta^{*}$. To provide a
flexible and extendable test bed, benchmarks are defined as the association of
three types of object classes:

.1 benchmark/. .2 objective.py. .2 datasets/. .3 dataset1.py. .3 dataset2.py.
.2 solvers/. .3 solver1.py. .3 solver2.py.

Figure 2: Standard benchmark structure

  * Objective: It defines the function $f$ to be minimized as well as the hyperparameters $\Lambda$ or the set $\Theta$, and the metrics to track along the iterations (e.g., objective value, gradient norm for smooth problems, or validation loss). Multiple metrics can be registered for each $\theta^{t}$.

  * Datasets: The Dataset objects provide the data $\mathcal{D}$ to be passed to the Objective class. They control how data is loaded and preprocessed. Datasets are separated from the Objective, making it easy to add new ones, provided they are coherent with the Objective.

  * Solvers: The Solver objects define how to run the algorithm. They are provided with the Objective and Dataset objects and output a sequence $\\{\theta^{t}\\}_{t}$. This sequence can be obtained using a single run of the method, or with multiple runs in case the method only returns its final iterate.

Each of these objects can have parameters that change their behavior, e.g.,
the regularization parameters for the Objective, the choice of preprocessing
for the Datasets, or the step size for the Solvers. By exposing these
parameters in the different objects, Benchopt can evaluate their influence on
the benchmark results. The Benchopt library defines an application programming
interface (API) for each of these concepts and provides a command line
interface (CLI) to make them work together. A benchmark is defined as a folder
that contains an Objective as well as subfolders containing the Solvers and
Datasets. Code examples of the three classes can be found in Appendix C.
Benchopt is written in Python, but Solvers run with implementations in
different languages (e.g., R and Julia, as in Section 4) and frameworks (e.g.,
PyTorch and TensorFlow, as in Section 5).

For each Dataset and Solver, and for each set of parameters, Benchopt
retrieves a sequence $\\{\theta^{t}\\}_{t}$ and evaluates the metrics defined
in the Objective for each $\theta^{t}$. To ensure fair evaluation, the
computation of these metrics is done off-line. The metrics are gathered in a
CSV file that can be used to display the benchmark results, either locally or
as HTML files published on a website that reference the benchmarks run with
Benchopt. This workflow is described in Figure 1.

This modular and standardized organization for benchmarks empowers the
optimization community by making numerical experiments easily reproducible,
shareable, flexible and extendable. The benchmark can be shared as a git
repository or a folder containing the different definitions for the Objective,
Datasets and Solvers and it can be run with the Benchopt CLI, hence becoming a
convenient reference for future comparisons. This also ensures fair evaluation
of baselines in follow-up experiments, as implementations validated by the
community are available.

Finally, by supporting multiple metrics the Objective class limits the risk of
mismatch in the evaluation. Moreover, benchmarks can be run in parallel with
checkpointing of the results, enabling large scale evaluations on many CPU or
GPU nodes. Benchmark results are automatically exported as interactive
visualizations, helping with the exploration of the many different settings.

Benchmarks All presented benchmarks are run on 10 cores of an Intel Xeon Gold
6248 CPUs @ 2.50GHz and NVIDIA V100 GPUs (16GB). The results’ interactive
plots and data are available at
https://benchopt.github.io/results/preprint_results.html.

##  3 First example: $\ell_{2}$-regularized logistic regression

Logistic regression is a very popular method for binary classification. From a
design matrix $X\in\mathbb{R}^{n\times p}$ with rows $X_{i}$ and a vector of
labels $y\in\\{-1,1\\}^{n}$ with corresponding element $y_{i}$,
$\ell_{2}$-regularized logistic regression provides a generalized linear model
indexed by $\theta^{*}\in\mathbb{R}^{p}$ to discriminate the classes by
solving

$\theta^{*}=\mathop{\mathrm{arg\,min}}_{\theta\in\mathbb{R}^{p}}\sum_{i=1}^{n}\log\big{(}1+\exp(-y_{i}X_{i}^{\top}\theta)\big{)}+\frac{\lambda}{2}\|\theta\|_{2}^{2}\enspace,$
(1)

where $\lambda>0$ is the regularization hyperparameter. Thanks to the
regularization part, Section 3 is strongly convex with a Lipschitz gradient,
and thus its solution can be estimated efficiently using many iterative
optimization schemes.

The most classical methods to solve this problem take inspiration from
Newton’s method [133]. On the one hand, quasi-Newton methods aim at
approximating the Hessian of the cost function with cheap to compute
operators. Among these methods, L-BFGS [85] stands out for its small memory
footprint, its robustness and fast convergence in a variety of settings. On
the other hand, truncated Newton methods [41] try to directly approximate
Newton’s direction by using e.g.,  the conjugate gradient method [51] and
Hessian-vector products to solve the associated linear system. Yet, these
methods suffer when $n$ is large: each iteration requires a pass on the whole
dataset.

In this context, methods based on stochastic estimates of the gradient have
become standard [18], with Stochastic Gradient Descent (SGD) as a main
instance. The core idea is to use cheap and noisy estimates of the gradient
[113, 73]. While SGD generally converges either slowly due to decreasing step
sizes, or to a neighborhood of the solution for constant step sizes, variance-
reduced adaptations such as SAG [118], SAGA [40] and SVRG [69] make it
possible to solve the problem more efficiently and are often considered to be
state-of-the-art for large scale problems.

Finally, methods based on coordinate descent [11] have also been proposed to
solve Section 3. While these methods are usually less popular, they can be
efficient in the context of sparse datasets, where only few samples have non-
zero values for a given feature, or when accelerated on distributed systems or
GPU [43].

Figure 3:  Benchmark for the $\ell_{2}$-regularized logistic regression, on 13
solvers, 4 datasets (_rows_), and 3 variants of the Objective (_columns_) with
$\lambda=1$. The curves display the suboptimality of the iterates,
$f(\theta^{t})-f(\theta^{*})$, as a function of time. The first column
corresponds to the objective function detailed in Section 3. In the second
column, datasets were preprocessed by normalizing each feature to unit
standard deviation. The third column is for an objective function which
includes an unregularized intercept.

The code for the benchmark is available at
https://github.com/benchopt/benchmark_logreg_l2/. To reflect the diversity of
solvers available, we showcase a Benchopt benchmark with 3 datasets, 10
optimization strategies implemented in 5 packages, leveraging GPU hardware
when possible. We also consider different scenarios for the objective
function: (i) scaling (or not) the features, a recommended data preprocessing
step, crucial in practice to have comparable regularization strength on all
variables; (ii) fitting (or not) an unregularized intercept term, important in
practice and making optimization harder when omitted from the regularization
term [75]; (iii) working (or not) with sparse features, which prevent explicit
centering during preprocessing to keep memory usage limited. Details on
packages, datasets and additional scenarios are available in Appendix D.

Results Figure D.1 presents the results of the benchmarks, in terms of
suboptimality of the iterates, $f(\theta^{t})-f(\theta^{*})$, for three
datasets and three scenarios. Here, because the problem is convex,
$\theta^{*}$ is approximated by the best iterate across all runs (see Section
C.5). Overall, the benchmark shows the benefit of using GPU solvers (cuML and
snapML), even for small scale tasks such as _ijcnn1_. It can also be noted
that these two accelerated solvers converge to a higher suboptimality level
compared to other solvers, due to operating with 32-bit float precision.
Another observation is that data scaling can drastically change the picture.
In the case of _madelon_ , most solvers have a hard time converging for the
scaled data. For the solvers that converge, we note that the convergence time
is one order of magnitude smaller with the scaled dataset compared to the raw
one. This stems from the fact that in this case, the scaling improves the
conditioning of the dataset.111The condition number of the dataset is divided
by 5.9 after scaling. For _news20.binary_ , the stochastic solvers such as SAG
and SAGA have degraded performances on scaled data. Here, the scaling makes
the problem harder.222The condition number is multiplied by 407 after scaling.

On CPU, quasi-Newton solvers are often the most efficient ones, and provide a
reasonable choice in most situations. For large scale _news20.binary_ ,
stochastic solvers such as SAG, SAGA or SVRG –that are often considered as
state of the art for such problem– have worst performances for the presented
datasets. While this dataset is often used as a test bed for benchmarking new
stochastic solvers, we fail to see an improvement over non-stochastic ones for
this experimental setup. In contrast, the last row in Figure D.1 displays an
experiment with the larger scale _criteo_ dataset, which demonstrates a regime
where variance-reduced stochastic gradient methods outperform quasi-Newton
methods. For future benchmarking of stochastic solvers, we therefore recommend
using such a large dataset.

Finally, the third column in Figure D.1 illustrates a classical problem when
benchmarking different solvers: their specific (and incompatible) definition
and resolution of the corresponding optimization problem. Here, the objective
function is modified to account for an intercept (bias) in the linear model.
In most situations, this intercept is not regularized when it is fitted.
However, snapML and liblinear solvers do regularize it, leading to
incomparable losses.

##  4 Second example: The Lasso

The Lasso, [126, 30], is an archetype of non-smooth ML problems, whose impact
on ML, statistics and signal processing in the last three decades has been
considerable [26, 61]. It consists of solving

$\theta^{*}\in\mathop{\mathrm{arg\,min}}_{\theta\in\mathbb{R}^{p}}\tfrac{1}{2}\left\|y-X\theta\right\|^{2}+\lambda\left\|\theta\right\|_{1}\enspace,$
(1)

where $X\in\mathbb{R}^{n\times p}$ is a design matrix containing $p$ features
as columns, $y\in\mathbb{R}^{n}$ is the target vector, and $\lambda>0$ is a
regularization hyperparameter. The Lasso estimator was popularized for
variable selection: when $\lambda$ is high enough, many entries in
$\theta^{*}$ are exactly equal to 0. This leads to more interpretable models
and reduces overfitting compared to the least-squares estimator.

Solvers for Section 4 have evolved since its introduction by [126]. After
generic quadratic program solvers, new dedicated solvers were proposed based
on iterative reweighted least-squares (IRLS) [58], followed by LARS [44], a
homotopy method computing the full Lasso path333The Lasso path is the set of
solutions of Section 4 as $\lambda$ varies in $(0,\infty)$.. The LARS solver
helped popularize the Lasso, yet the algorithm suffers from stability issues
and can be very slow for worst case situations [90]. General purpose solvers
became popular for Lasso-type problems with the introduction of the iterative
soft thresholding algorithm (ISTA, [39]), an instance of forward-backward
splitting [34]. These algorithms became standard in signal and image
processing, especially when accelerated (FISTA, [7]).

In parallel, proximal coordinate descent has proven particularly relevant for
the Lasso in statistics. Early theoretical results were proved by [129, 117],
before it became the standard solver of the widely distributed packages glmnet
in R and scikit-learn in Python. For further improvements, some solvers
exploit the sparsity of $\theta^{*}$, trying to identify its support to reduce
the problem size. Best performing variants of this scheme are screening rules
(e.g.,  [45, 17, 98]) and working/active sets (e.g.,  [70, 92]), including
strong rules [125].

While reviews of Lasso solvers have already been performed [2, Sec. 8.1], they
are limited to certain implementation and design choices, but also naturally
lack comparisons with more recent solvers and modern hardware, hence drawing
biased conclusions.

The code for the benchmark is available at
https://github.com/benchopt/benchmark_lasso/. Results obtained on 4 datasets,
with 9 standard packages and some custom reimplementations, possibly
leveraging GPU hardware, and 17 different solvers written in
Python/numba/Cython, R, Julia or C++ (Table E.1) are presented in Figure 4.
All solvers use efficient numerical implementations, possibly leveraging calls
to BLAS, precompiled code in Cython or just-in-time compilation with numba.

The different parameters influencing the setup are

  * •

the regularization strength $\lambda$, controlling the sparsity of the
solution, parameterized as a fraction of
$\lambda_{\max}=\left\|X^{\top}y\right\|_{\infty}$ (the minimal hyperparameter
such that $\theta^{*}=0$),

  * •

the dataset dimensions: _MEG_ has small $n$ and medium $p$; _rcv1.binary_ has
medium $n$ and $p$; _news20.binary_ has medium $n$ and very large $p$ while
_MillionSong_ has very large $n$ and small $p$ (Table E.2).

Figure 4:  Benchmark for the Lasso, on 17 solvers, 4 datasets (_rows_), and 3
variants of the Objective (_columns_) with decreasing regularization
$\lambda$. The curves display the suboptimality of the objective function,
$f(\theta^{t})-f(\theta^{*})$, as a function of time.

Results Figure 4 presents the result of the benchmark on the Lasso, in terms
of objective suboptimality $f(\theta^{t})-f(\theta^{*})$ as a function of
time.

Similarly to Section 3, the GPU solvers obtain good performances in most
settings, but their advantage is less clear. A consistent finding across all
settings is that coordinate descent-based methods outperform full gradient
ones (ISTA and FISTA, even restarted), and are improved by the use of working
set strategies (blitz, celer, skglm, glmnet). This observation is even more
pronounced when the regularization parameter is large, as the solution is
sparser.

When observing the influence of the dataset dimensions, we observe 3 regimes.
When $n$ is small (_MEG_), the support of the solution is small and coordinate
descent, LARS and noncvx-pro perform the best. When $n$ is much larger than
$p$ (_MillionSong_), noncvx-pro clearly outperforms other solvers, and working
set methods prove useless. Finally, when $n$ and $p$ are large (_rcv1.binary_
, _news20.binary_), CD and working sets vastly outperforms the rest while
noncvx-pro fails, as it requires solving a linear system of size $\min(n,p)$.
We note that this setting was not tested in the original experiment of [106],
which highlights the need for extensive, standard experimental setups.

When the support of the solution is small (either small $\lambda$, either
small $n$ since the Lasso solution has at most $n$ nonzero coefficients), LARS
is a competitive algorithm. We expect this to degrade when $n$ increases, but
as the LARS solver in scikit-learn does not support sparse design matrices we
could not include it for _news20.binary_ and _rcv1.binary_.

This benchmark is the first to evaluate solvers across languages, showing the
competitive behavior of lasso.jl and glmnet compared to Python solvers. We
observe that both solvers have a large initialization time, and then converge
very fast. To ensure that the benchmark is fair, even though the Benchopt
library is implemented in Python, we made sure to track conversion overhead,
as well as just-in-time compilation cost, and ignore them. We also checked the
timing’s consistency with native calls to the libraries.

Since the Lasso is massively used for it feature selection properties, the
speed at which the solvers identify the support of the solution is also an
important performance measure. Monitoring this with Benchopt is
straightforward, and a figure reporting this benchmark is in Appendix E.

##  5 Third example: How standard is a benchmark on ResNet18?

As early successes of deep learning have been focused on computer vision tasks
[80], image classification has become a _de facto_ standard to validate novel
methods in the field. Among the different network architectures, ResNets [63]
are extensively used in the community as they provide strong and versatile
baselines [134, 124, 42, 25, 86]. While many papers present results with such
model on classical datasets, with sometimes extensive ablation studies [64,
132, 8, 119], the lack of standardized codebase and missing implementation
details makes it hard to replicate their results.

The code for the benchmark is available at
https://github.com/benchopt/benchmark_resnet_classif/. We provide a cross-
dataset – _SVHN_ , [100]; _MNIST_ , [82] and _CIFAR-10_ , [79]– and cross-
framework –TensorFlow/Keras, [1, 32]; PyTorch, [101]– evaluation of the
training strategies for image classification with ResNet18 (see Appendix F for
details on architecture and datasets). We train the network by minimizing the
cross entropy loss relatively to the weights $\theta$ of the model. Contrary
to logistic regression and the Lasso, this problem is non-convex due to the
non-linearity of the model $f_{\theta}$. Another notable difference is that we
report the evolution of the test error rather than the training loss.

The Solvers are defined as the combination of an optimization algorithm, its
hyperparameters, the learning rate (LR) and weight decay schedules, and the
data augmentation strategy. We focus on 2 standard methods: stochastic
gradient descent (SGD) with momentum and Adam [74], as well as a more recently
published one: Lookahead [136]. The LR schedules are chosen among fixed LR,
step LR444decreasing the learning rate by a factor 10 at mid-training, and
again at 3/4 of the training, and cosine annealing [87]. We also consider
decoupled weight decay for Adam [88], and coupled weight decay (i.e.,
$\ell_{2}$-regularization) for SGD. Regarding data augmentation, we use random
cropping for all datasets and add horizontal flipping only for _CIFAR-10_ , as
the digits datasets do not exhibit a mirror symmetry. We detail the remaining
hyperparameters in Table F.2.

Aligning cross-framework implementations Due to some design choices,
components with the same name in the different frameworks do not have the same
behavior. For instance, when it comes to applying weight decay, PyTorch’s SGD
uses coupled weight decay, while in TensorFlow/Keras weight decay always
refers to decoupled weight decay. These two methods lead to significantly
different performance and it is not straightforward to apply coupled weight
decay in a post-hoc manner in TensorFlow/Keras (see further details in Section
F.2). We conducted an extensive effort to align them using unit testing to
make the conclusions of our benchmarks independent of the chosen framework. We
found additional significant differences (reported in Table F.3) in the
initialization, the batch normalization, the convolutional layers and the
weight decay scaling.

Results

Figure 5: ResNet18 image classification benchmark with PyTorch Solvers. The
best SGD configuration features data augmentation, momentum, cosine learning
rate schedule and weight decay. In dashed black is the state of the art for
the corresponding datasets with a ResNet18 measured by [136] for _CIFAR-10_ ,
by [138] for _SVHN_ with a PreAct ResNet18, by PapersWithCode for _MNIST_ with
all networks considered. Off-the-shelf ResNet implementations in
TensorFlow/Keras do not support images smaller than $32\times 32$ and is hence
not shown for _MNIST_. Curves are exponentially smoothed.

The results of the benchmark are reported in Figure 5. Each graph reports the
test error relative to time, with an ablation study on the solvers parameters.
Note that we only report selected settings for clarity but that we run every
possible combinations. The results are available in the supplementary
materials as a user-friendly interactive html file.

Firstly, we see that reaching the state of the art for a vanilla ResNet18 is
not straightforward. On the popular website Papers with code it has been so
far underestimated. It can achieve 4.45% and 2.65% test error rates on
_CIFAR-10_ and _SVHN_ respectively (compared to 4.73% and 2.95% – for a PreAct
ResNet18 – before that). What our ablation study shows, is that a variety of
techniques is required to reach it. The most significant one is an appropriate
data augmentation strategy, which lowers the error rate on _CIFAR-10_ from
about 18% to about 8%. The second most important one is weight decay, but it
has to be used in combination with a proper LR schedule, as well as momentum.
While these techniques are not novel, they are regularly overlooked in
baselines, resulting in underestimation of their performance level.

This reproducible benchmark not only allows a researcher to get a clear
understanding of how to achieve the best performances for this model and
datasets, but also provides a way to reproduce and extend these performances.
In particular, we also include in this benchmark the original implementation
of Lookahead [136]. We confirm that it slightly accelerates the convergence of
the Best SGD, even with a cosine LR schedule – a setting that had not been
studied in the original paper.

Our benchmark also evaluates the relative computational performances of the
different frameworks. We observe that PyTorch-Lightning is significantly
slower than the other frameworks we tested, in large part due to their
callbacks API. We also notice that our TensorFlow/Keras implementation is
significantly slower ($\approx$ 28%) than the PyTorch ones, despite following
the best practices and our profiling efforts. Note that we do not imply that
TensorFlow is intrinsically slower than PyTorch, but a community effort is
needed to ensure that the benchmark performances are framework-agnostic.

Figure 6: ResNet18 image classification benchmark with a validation split. In
dashed black is the state of the art (see caption of Figure 5 for more
details). In addition, we show in colored horizontal dashed lines, the test
results for early stopping on the validation and on the test set for the
different solvers, the square mark indicating the moment this stopping would
happen. The curves for the train-val splits show the exponentially smoothed
median results for five different random seeds.

A recurrent criticism of such benchmarks is that only the best test error is
reported. In Figure 6, we measure the effect of using a train-validation-test
split, by keeping a fraction of the training set as a validation set. The
splits we use are detailed in Table F.1. Our finding is that the results of
the ablation study do not change significantly when using such procedure, even
though their validity is reinforced by the use of multiple trainings. Yet, a
possible limitation of our findings is that some of the hyperparameters we
used for our study, coming from the PyTorch-CIFAR GitHub repository, may have
been tuned while looking at the test set.

##  6 Conclusion and future work

We have introduced Benchopt, a library that makes it easy to collaboratively
develop fair and extensive benchmarks of optimization algorithms, which can
then be seamlessly published, reproduced, and extended. In the future, we plan
on supporting the creation of new benchmarks, that could become the standards
the community builds on. This work is part of a wider effort to improve
reproducibility of machine learning results. It aims to contribute to raising
the standard of numerical validation for optimization, which is pervasive in
the statistics and ML community as well as for the experimental sciences that
rely more and more on these tools for research.

## References

  * [1] Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Soprano Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, D. Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu and Xiaoqiang Zheng “TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems” Software available from tensorflow.org, 2015
  * [2] Francis Bach, Rodolphe Jenatton, Julien Mairal and Guillaume Obozinski “Optimization with sparsity-inducing penalties” In _Foundations and Trends in Machine Learning_ 4.1, 2012, pp. 1–106
  * [3] E. Bacry, M. Bompaire, S. Gaïffas and S. Poulsen “Tick: A Python Library for Statistical Learning, with a Particular Emphasis on Time-Dependent Modeling” In _ArXiv e-prints_ , 2017 eprint: 1707.03003
  * [4] Alvero Barbero and Suvrit Sra “Modular Proximal Optimization for Multidimensional Total-Variation Regularization” In _Journal of Machine Learning Research_ 19.56, 2018, pp. 1–82 URL: http://jmlr.org/papers/v19/13-538.html
  * [5] R.. Barlow and H.. Brunk “The Isotonic Regression Problem and Its Dual” In _Journal of the American Statistical Association_ 67.337 [American Statistical Association, Taylor & Francis, Ltd.], 1972, pp. 140–147
  * [6] T. Bartz-Beielstein, C. Doerr, D. Berg, J. Bossek, S. Chandrasekaran, T. Eftimov, A. Fischbach, P. Kerschke, W. La Cava and M. Lopez-Ibanez “Benchmarking in optimization: Best practice and open issues” In _arXiv preprint arXiv:2007.03488_ , 2020
  * [7] A. Beck and M. Teboulle “A fast iterative shrinkage-thresholding algorithm for linear inverse problems” In _SIAM J. Imaging Sci._ 2.1, 2009, pp. 183–202
  * [8] I. Bello, W. Fedus, X. Du, E.. Cubuk, A. Srinivas, T.-Y. Lin, J. Shlens and B. Zoph “Revisiting ResNets: Improved Training and Scaling Strategies” In _Advances in Neural Information Processing Systems_ , 2021 arXiv:2103.07579v1
  * [9] Thierry Bertin-Mahieux, Daniel P.W. Ellis, Brian Whitman and Paul Lamere “The Million Song Dataset” In _Proceedings of the 12th International Conference on Music Information Retrieval (ISMIR 2011)_ , 2011
  * [10] Q. Bertrand, Q. Klopfenstein, P.-A. Bannier, G. Gidel and M. Massias “Beyond L1: Faster and Better Sparse Models with skglm” In _arXiv preprint arXiv:2204.07826_ , 2022
  * [11] D.. Bertsekas “Nonlinear programming” Athena Scientific, 1999
  * [12] Jeff Bezanson, Alan Edelman, Stefan Karpinski and Viral B Shah “Julia: A fresh approach to numerical computing” In _SIAM Review_ 59.1 SIAM, 2017, pp. 65–98 DOI: 10.1137/141000671
  * [13] Kevin Bleakley and Jean-Philippe Vert “The group fused Lasso for multiple change-point detection” arXiv, 2011
  * [14] Mathieu Blondel and Fabian Pedregosa “Lightning: Large-Scale Linear Classification, Regression and Ranking in Python”, 2016
  * [15] A. Boisbunon, R. Flamary and A. Rakotomamonjy “Active set strategy for high-dimensional non-convex sparse optimization problems” In _ICASSP_ , 2014, pp. 1517–1521 IEEE
  * [16] Jérôme Bolte, Shoham Sabach and Marc Teboulle “Proximal alternating linearized minimization for nonconvex and nonsmooth problems” In _Mathematical Programming_ 146.1 Springer, 2014, pp. 459–494
  * [17] A. Bonnefoy, V. Emiya, L. Ralaivola and R. Gribonval “Dynamic screening: accelerating first-order algorithms for the Lasso and Group-Lasso” In _IEEE Trans. Signal Process._ 63.19, 2015, pp. 20
  * [18] Léon Bottou “Large-Scale Machine Learning with Stochastic Gradient Descent” In _COMPSTAT_ Physica-Verlag, 2010, pp. 177–186
  * [19] Stephen Boyd, Neal Parikh, Eric Chu, Borja Peleato and Jonathan Eckstein “Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers” In _Foundations and Trends in Machine Learning_ 3.1, 2011
  * [20] Stephen P. Boyd and Lieven Vandenberghe “Convex Optimization” Cambridge, UK ; New York: Cambridge University Press, 2004
  * [21] Y. Boykov, O. Veksler and R. Zabih “Fast approximate energy minimization via graph cuts” In _IEEE Transactions on Pattern Analysis and Machine Intelligence_ 23.11, 2001, pp. 1222–1239 DOI: 10.1109/34.969114
  * [22] Georg Brandl “Sphinx documentation” In _URL http://sphinx-doc. org/sphinx. pdf_ , 2010
  * [23] Kristian Bredies, Karl Kunisch and Thomas Pock “Total generalized variation” In _SIAM J. Imaging Sci._ 3.3 SIAM, 2010, pp. 492–526
  * [24] P. Breheny and J. Huang “Coordinate descent algorithms for nonconvex penalized regression, with applications to biological feature selection” In _Ann. Appl. Stat._ 5.1 NIH Public Access, 2011, pp. 232
  * [25] A. Brock, S. De, S.. Smith and K. Simonyan “High-performance large-scale image recognition without normalization” In _ICML_ , 2021, pp. 1059–1071 PMLR
  * [26] P. Bühlmann and S. van de Geer “Statistics for high-dimensional data” Methods, theory and applications, Springer Series in Statistics Heidelberg: Springer, 2011
  * [27] E.. Candès, M.. Wakin and S.. Boyd “Enhancing Sparsity by Reweighted $l_{1}$ Minimization” In _J. Fourier Anal. Applicat._ 14.5-6, 2008, pp. 877–905
  * [28] Antonin Chambolle and Pierre-Louis Lions “Image recovery via total variation minimization and related problems” In _Numerische Mathematik_ 76.2 Springer, 1997, pp. 167–188
  * [29] Antonin Chambolle and Thomas Pock “A first-order primal-dual algorithm for convex problems with applications to imaging” In _Journal of Mathematical Imaging and Vision_ 40.1, 2011
  * [30] S.. Chen, D.. Donoho and M.. Saunders “Atomic decomposition by basis pursuit” In _SIAM J. Sci. Comput._ 20.1, 1998, pp. 33–61
  * [31] Hamza Cherkaoui, Thomas Moreau, Abderrahim Halimi and Philippe Ciuciu “Sparsity-Based Semi-Blind Deconvolution of Neural Activation Signal in fMRI” In _IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2019
  * [32] Francois Chollet “Keras” GitHub, 2015
  * [33] Alex Clark “Pillow (PIL Fork) Documentation” readthedocs, 2015 URL: https://buildmedia.readthedocs.org/media/pdf/pillow/latest/pillow.pdf
  * [34] P.. Combettes and V. Wajs “Signal recovery by proximal forward-backward splitting” In _Multiscale modeling & simulation_ 4.4 SIAM, 2005, pp. 1168–1200
  * [35] Patrick L. Combettes and Lilian E. Glaudin “Solving Composite Fixed Point Problems with Block Updates” In _Advances in Nonlinear Analysis_ , 2021
  * [36] Laurent Condat “A Direct Algorithm for 1-D Total Variation Denoising” In _IEEE SIGNAL PROCESSING LETTERS_ 20.12, 2013
  * [37] Laurent Condat “A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms” In _Journal of Optimization Theory and Applications, Springer Verlag_ , 2013
  * [38] Criteo-Labs “Criteo releases industry’s largest-ever dataset for machine learning to academic community”, 2015 URL: https://www.criteo.com/news/press-releases/2015/07/criteo-releases-industrys-largest-ever-dataset
  * [39] I. Daubechies, M. Defrise and C. De Mol “An iterative thresholding algorithm for linear inverse problems with a sparsity constraint” In _Commun. Pure Appl. Math._ 57.11 Wiley Online Library, 2004, pp. 1413–1457
  * [40] Aaron Defazio, Francis Bach and Simon Lacoste-Julien “SAGA: A Fast Incremental Gradient Method With Support for Non-Strongly Convex Composite Objectives” In _Advances in Neural Information Processing Systems_ 28, 2014, pp. 1646–1654
  * [41] Ron S Dembo, Stanley C Eisenstat and Trond Steihaug “Inexact Newton Methods” In _SIAM J. Numer. Anal._ 19.2 SIAM, 1982, pp. 400–408
  * [42] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, Jakob Uszkoreit and Neil Houlsby “An Image is Worth 16x16 Words: Transformers for Image Recognition at Scale” In _ICLR_ , 2021
  * [43] C. Dünner, T. Parnell, D. Sarigiannis, N. Ioannou, A. Anghel, G. Ravi, M. Kandasamy and H. Pozidis “Snap ML: A hierarchical framework for machine learning” In _Advances in Neural Information Processing Systems_ 31, 2018
  * [44] B. Efron, T.. Hastie, I.. Johnstone and R. Tibshirani “Least angle regression” With discussion, and a rejoinder by the authors In _Ann. Statist._ 32.2, 2004, pp. 407–499
  * [45] L. El Ghaoui, V. Viallon and T. Rabbani “Safe feature elimination in sparse supervised learning” In _J. Pacific Optim._ 8.4, 2012, pp. 667–698
  * [46] Michael Elad, Peyman Milanfar and Ron Rubinstein “Analysis versus synthesis in signal priors” In _2006 14th European Signal Processing Conference_ , 2006
  * [47] William Falcon, Jirka Borovec, Adrian Wälchli, Nic Eggert, Justus Schock, Jeremy Jordan, Nicki Skafte, Ir1dXD, Vadim Bereznyuk, Ethan Harris, Tullie Murrell, Peter Yu, Sebastian Præsius, Travis Addair, Jacob Zhong, Dmitry Lipin, So Uchida, Shreyas Bapat, Hendrik Schröter, Boris Dayma, Alexey Karnachev, Akshay Kulkarni, Shunta Komatsu, Martin.B, Jean-Baptiste SCHIRATTI, Hadrien Mary, Donal Byrne, Cristobal Eyzaguirre, cinjon and Anton Bakhtin “PyTorchLightning/pytorch-lightning: 0.7.6 release” Zenodo, 2020
  * [48] J. Fan and R. Li “Variable selection via nonconcave penalized likelihood and its oracle properties” In _J. Amer. Statist. Assoc._ 96.456, 2001, pp. 1348–1360
  * [49] Rong-En Fan, Kai-Wei Chang, Cho-Jui Hsieh, Xiang-Rui Wang and Chih-Jen Lin “LIBLINEAR: A Library for Large Linear Classification” In _J. Mach. Learn. Res._ 9 Cambridge, MA, USA: MIT Press, 2008
  * [50] S Farrens, A Grigis, L El Gueddari, Z Ramzi, GR Chaithya, S Starck, B Sarthou, H Cherkaoui, P Ciuciu and J-L Starck “PySAP: Python Sparse Data Analysis Package for multidisciplinary image processing” In _Astronomy and Computing_ 32 Elsevier, 2020, pp. 100402
  * [51] R. Fletcher and C.. Reeves “Function Minimization by Conjugate Gradients” In _The Computer Journal_ 7.2, 1964, pp. 149–154 eprint: https://academic.oup.com/comjnl/article-pdf/7/2/149/959725/070149.pdf
  * [52] J. Forde, T. Head, C. Holdgraf, Y. Panda, G. Nalvarete, B. Ragan-Kelley and E. Sundell “Reproducible research environments with repo2docker”, 2018
  * [53] J. Friedman, T. Hastie and R. Tibshirani “Regularization paths for generalized linear models via coordinate descent” In _J. Stat. Softw._ 33.1 NIH Public Access, 2010, pp. 1–22
  * [54] X. Glorot and Y. Bengio “Understanding the difficulty of training deep feedforward neural networks” In _AISTATS_ 9, 2010, pp. 249–256
  * [55] Todd R Golub, Donna K Slonim, Pablo Tamayo, Christine Huard, Michelle Gaasenbeek, Jill P Mesirov, Hilary Coller, Mignon L Loh, James R Downing and Mark A Caligiuri “Molecular classification of cancer: class discovery and class prediction by gene expression monitoring” In _science_ 286.5439 American Association for the Advancement of Science, 1999, pp. 531–537
  * [56] P. Gong, C. Zhang, Z. Lu, J. Huang and J. Ye “A general iterative shrinkage and thresholding algorithm for non-convex regularized optimization problems” In _ICML_ , 2013, pp. 37–45
  * [57] A. Gramfort, M. Luessi, E. Larson, D.. Engemann, D. Strohmeier, C. Brodbeck, L. Parkkonen and M.. Hämäläinen “MNE software for processing MEG and EEG data” In _NeuroImage_ 86, 2014, pp. 446–460
  * [58] Y. Grandvalet “Least absolute shrinkage is equivalent to quadratic penalization” In _International Conference on Artificial Neural Networks_ , 1998, pp. 201–206 Springer
  * [59] Isabelle Guyon, Steve Gunn, Asa Ben-Hur and Gideon Dror “Result analysis of the nips 2003 feature selection challenge” In _Advances in neural information processing systems_ 17, 2004
  * [60] Charles R Harris, K Jarrod Millman, Stéfan J Van Der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg and Nathaniel J Smith “Array programming with NumPy” In _Nature_ 585.7825 Nature Publishing Group, 2020, pp. 357–362
  * [61] T.. Hastie, R. Tibshirani and M. Wainwright “Statistical Learning with Sparsity: The Lasso and Generalizations” CRC Press, 2015
  * [62] Kaiming He, Xiangyu Zhang, Shaoqing Ren and Jian Sun “Delving deep into rectifiers: Surpassing human-level performance on imagenet classification” In _CVPR_ , 2015, pp. 1026–1034
  * [63] Kaiming He, Xiangyu Zhang, Shaoqing Ren and Jian Sun “Deep residual learning for image recognition” In _CVPR_ , 2016
  * [64] Tong He, Zhi Zhang, Hang Zhang, Zhongyue Zhang, Junyuan Xie and Mu Li “Bag of tricks for image classification with convolutional neural networks” In _CVPR_ 2019-June, 2019, pp. 558–567 arXiv:1812.01187
  * [65] J.. Hunter “Matplotlib: A 2D graphics environment” In _Computing in Science & Engineering_ 9.3 IEEE COMPUTER SOC, 2007, pp. 90–95 DOI: 10.1109/MCSE.2007.55
  * [66] M. Hutson “Artificial intelligence faces reproducibility crisis” In _Science_ 359.6377 American Association for the Advancement of Science, 2018, pp. 725–726
  * [67] Plotly Technologies Inc. “Collaborative data science” Montreal, QC: Plotly Technologies Inc., 2015 URL: https://plot.ly
  * [68] S. Ioffe and C. Szegedy “Batch normalization: Accelerating deep network training by reducing internal covariate shift” In _ICML_ , 2015, pp. 448–456
  * [69] Rie Johnson and Tong Zhang “Accelerating Stochastic Gradient Descent Using Predictive Variance Reduction” In _Advances in Neural Information Processing Systems_ 26, 2013
  * [70] T.. Johnson and C. Guestrin “Blitz: A Principled Meta-Algorithm for Scaling Sparse Optimization” In _ICML_ 37, 2015, pp. 1171–1179
  * [71] Fikret Işık Karahanoğlu, César Caballero-Gaudes, François Lazeyras and Dimitri Van De Ville “Total Activation: fMRI Deconvolution through Spatio-Temporal Regularization” In _NeuroImage_ 73, 2013, pp. 121–134
  * [72] S Sathiya Keerthi, Dennis DeCoste and Thorsten Joachims “A modified finite Newton method for fast solution of large scale linear SVMs.” In _Journal of Machine Learning Research_ 6.3, 2005
  * [73] Jack Kiefer and Jacob Wolfowitz “Stochastic estimation of the maximum of a regression function” In _Ann. Math. Stat._ , 1952, pp. 462–466
  * [74] Diederik P. Kingma and Jimmy Ba “Adam: A Method for Stochastic Optimization” In _ICLR_ , 2015, pp. 1–10
  * [75] K. Koh, J. Kim and S. Boyd “An interior-point method for large-scale l1-regularized logistic regression.” In _J. Mach. Learn. Res._ 8.8, 2007, pp. 1519–1555
  * [76] V. Kolmogorov and R. Zabin “What energy functions can be minimized via graph cuts?” In _IEEE Transactions on Pattern Analysis and Machine Intelligence_ 26.2, 2004, pp. 147–159 DOI: 10.1109/TPAMI.2004.1262177
  * [77] Nikos Komodakis and Jean-Christophe Pesquet “Playing with Duality: An overview of recent primal-dual approaches for solving large-scale optimization problems” In _IEEE Signal Processing Magazine, Institute of Electrical and Electronics Engineers_ , 2015
  * [78] Simon Kornblith “Lasso.jl”, 2021 JuliaStats URL: https://github.com/JuliaStats/Lasso.jl
  * [79] Alex Krizhevsky “Learning multiple layers of features from tiny images”, 2009
  * [80] Alex Krizhevsky, Ilya Sutskever and Geoffrey E. Hinton “Imagenet Classification with Deep Convolutional Neural Networks” In _Advances in Neural Information Processing Systems_ , 2012, pp. 1097–1105
  * [81] Clement Lalanne, Maxence Rateaux, Laurent Oudre, Matthieu P. Robert and Thomas Moreau “Extraction of Nystagmus Patterns from Eye-Tracker Data with Convolutional Sparse Coding” In _Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC)_ Montreal, QC, Canada: IEEE, 2020, pp. 928–931
  * [82] Yann LeCun, Corinna Cortes and CJ Burges “MNIST handwritten digit database” In _ATT Labs [Online]. Available: http://yann.lecun.com/exdb/mnist_ 2, 2010
  * [83] David D Lewis, Yiming Yang, Tony Russell-Rose and Fan Li “Rcv1: A new benchmark collection for text categorization research” In _Journal of machine learning research_ 5.Apr Goldsmiths, University of London, 2004, pp. 361–397
  * [84] Jingwei Liang, Tao Luo and Carola-Bibiane Schöenlieb “Improving “Fast Iterative Shrinkage-Thresholding Algorithm”: Faster, Smarter, and Greedier” In _SIAM J. Sci. Comput._ 44.3 SIAM, 2022, pp. A1069–A1091
  * [85] Dong C. Liu and Jorge Nocedal “On the Limited Memory BFGS Method for Large Scale Optimization” In _Math. Program._ 45.1-3, 1989, pp. 503–528
  * [86] Zhuang Liu, Hanzi Mao, Chao-Yuan Wu, Christoph Feichtenhofer, Trevor Darrell and Saining Xie “A ConvNet for the 2020s” In _CVPR_ , 2022
  * [87] Ilya Loshchilov and Frank Hutter “SGDR: Stochastic Gradient Descent with Warm Restarts” In _ICLR_ , 2017
  * [88] Ilya Loshchilov and Frank Hutter “Decoupled Weight Decay Regularization” In _ICLR_ , 2019
  * [89] Jan-Matthis Lueckmann, Jan Boelts, David S Greenberg, Pedro J Gonçalves and Jakob H Macke “Benchmarking Simulation-Based Inference” In _AISTATS_ 130, PMLR, 2021, pp. 343–351
  * [90] J. Mairal and B. Yu “Complexity analysis of the Lasso regularization path” In _ICML_ , 2012, pp. 353–360
  * [91] Julien Mairal “Cyanure: An open-source toolbox for empirical risk minimization for python, c++, and soon more” In _arXiv preprint arXiv:1912.08165_ , 2019
  * [92] M. Massias, A. Gramfort and J. Salmon “Celer: a fast solver for the lasso with dual extrapolation” In _ICML_ , 2018, pp. 3315–3324
  * [93] P. Mattson, V. Reddi, C. Cheng, C. Coleman, G. Diamos, D. Kanter, P. Micikevicius, D. Patterson, G. Schmuelling and H. Tang “MLPerf: An industry standard benchmark suite for machine learning performance” In _IEEE Micro_ 40.2 IEEE, 2020, pp. 8–16
  * [94] R. Mazumder, J.. Friedman and T. Hastie “Sparsenet: Coordinate descent with nonconvex penalties” In _J. Amer. Statist. Assoc._ 106.495 Taylor & Francis, 2011, pp. 1125–1138
  * [95] Wes McKinney “Data structures for statistical computing in python” In _Proceedings of the 9th Python in Science Conference_ 445, 2010, pp. 51–56 Austin, TX
  * [96] Thomas Moreau, Pierre Glaser, Roman Yurchak and Olivier Grisel “Loky”, 2017 DOI: 10.5281/zenodo.5770406
  * [97] Óscar Nájera, Eric Larson, Loïc Estève, Gael Varoquaux, Lucy Liu, Jaques Grobler, Elliott Sales Andrade, Chris Holdgraf, Alexandre Gramfort, Mainak Jas, Joel Nothman, Olivier Grisel, Nelle Varoquaux, Emmanuelle Gouillart, Martin Luessi, Antony Lee, Jake Vanderplas, Tim Hoffmann, Thomas A Caswell, Bane Sullivan, Alyssa Batula, jaeilepp, Thomas Robitaille, Stefan Appelhoff, Patrick Kunzmann, Matthias Geier, Lars, Kyle Sunden, Dominik Stańczak and Albert Y. Shih “sphinx-gallery/sphinx-gallery: Release v0.7.0” Zenodo, 2020 DOI: 10.5281/zenodo.3838216
  * [98] E. Ndiaye, O. Fercoq, A. Gramfort and J. Salmon “Gap Safe screening rules for sparsity enforcing penalties” In _J. Mach. Learn. Res._ 18.128, 2017, pp. 1–33
  * [99] Yurii E Nesterov “A method for solving the convex programming problem with convergence rate $O(1/k^{2})$” In _Dokl. akad. nauk Sssr_ 269, 1983, pp. 543–547
  * [100] Yuval Netzer, Tao Wang, Adam Coates, Alessandro Bissacco, Bo Wu and Andrew Y Ng “Reading Digits in Natural Images with Unsupervised Feature Learning” In _Advances in Neural Information Processing Systems_ , 2011
  * [101] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai and Soumith Chintala “PyTorch: An Imperative Style, High-Performance Deep Learning Library” In _Advances in Neural Information Processing Systems_ , 2019, pp. 8024–8035
  * [102] Fabian Pedregosa, Geoffrey Negiar and Gideon Dresdner “copt: composite optimization in Python”, 2020
  * [103] Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss and Vincent Dubourg “Scikit-learn: Machine learning in Python” In _J. Mach. Learn. Res._ 12 JMLR. org, 2011, pp. 2825–2830
  * [104] J. Pineau, P. Vincent-Lamarre, K. Sinha, V. Larivière, A. Beygelzimer, F. d’Alché-Buc, E. Fox and H. Larochelle “Improving reproducibility in machine learning research: a report from the NeurIPS 2019 reproducibility program” In _J. Mach. Learn. Res._ 22 Microtome Publishing, 2021
  * [105] Joelle Pineau, Koustuv Sinha, Genevieve Fried, Rosemary Nan Ke and Hugo Larochelle “ICLR reproducibility challenge 2019” In _ReScience C_ 5.2, 2019, pp. 5
  * [106] Clarice Poon and Gabriel Peyré “Smooth Bilevel Programming for Sparse Regularization” In _Advances in Neural Information Processing Systems_ 34, 2021, pp. 1543–1555
  * [107] William H Press, Saul A Teukolsky, William T Vetterling and Brian P Flannery “Numerical recipes 3rd edition: The art of scientific computing” Cambridge university press, 2007
  * [108] Danil Prokhorov “IJCNN 2001 neural network competition” In _Slide presentation in IJCNN_ 1.97, 2001, pp. 38
  * [109] Emmanuel Quemener and Marianne Corvellec “SIDUS—the solution for extreme deduplication of an operating system” In _Linux Journal_ 2013.235 Belltown Media Houston, TX, 2013, pp. 3
  * [110] R Core Team “R: A Language and Environment for Statistical Computing”, 2017 R Foundation for Statistical Computing URL: https://www.R-project.org/
  * [111] E. Raff “A step toward quantifying independently reproducible machine learning research” In _Advances in Neural Information Processing Systems_ 32, 2019, pp. 5486–5496
  * [112] Sebastian Raschka, Joshua Patterson and Corey Nolet “Machine Learning in Python: Main Developments and Technology Trends in Data Science, Machine Learning, and Artificial Intelligence” In _Information-an International Interdisciplinary Journal_ 11.4 Multidisciplinary Digital Publishing Institute, 2020, pp. 193
  * [113] Herbert Robbins and Sutton Monro “A stochastic approximation method” In _Ann. Math. Stat._ , 1951, pp. 400–407
  * [114] Giampaolo Rodola “Psutil package: a cross-platform library for retrieving information on running processes and system utilization” In _Google Scholar_ , 2016
  * [115] N.. Rougier and K. Hinsen “ReScience C: a journal for reproducible replications in computational science” In _International Workshop on Reproducible Research in Pattern Recognition_ , 2018, pp. 150–156 Springer
  * [116] Leonid I. Rudin, Stanley Osher and Emad Fatemi “Nonlinear total variation based noise removal algorithms” In _Physica D: Nonlinear Phenomena_ 60.1-4 Elsevier, 1992, pp. 259–268
  * [117] S. Sardy, A.. Bruce and P. Tseng “Block coordinate relaxation methods for nonparametric wavelet denoising” In _J. Comput. Graph. Stat._ 9.2 Taylor & Francis, 2000, pp. 361–379
  * [118] Mark Schmidt, Nicolas Le Roux and Francis Bach “Minimizing Finite Sums with the Stochastic Average Gradient” In _Math. Program._ 162.arXiv:1309.2388, 2017, pp. 83–112 arXiv:1309.2388 [cs, math, stat]
  * [119] Robin M. Schmidt, Frank Schneider and Philipp Hennig “Descending through a Crowded Valley - Benchmarking Deep Learning Optimizers” In _ICML_ 139 PMLR, 2021, pp. 9367–9376
  * [120] David Sculley, Jasper Snoek, Alex Wiltschko and Ali Rahimi “Winner’s curse? On pace, progress, and empirical rigor”, 2018
  * [121] Thalles Santos Silva “How to Add Regularization to Keras Pre-trained Models the Right Way” In _https://sthalles.github.io_ , 2019
  * [122] Karen Simonyan and Andrew Zisserman “Very Deep Convolutional Networks for Large-Scale Image Recognition” In _ICLR_ , 2015
  * [123] S. Sra, S. Nowozin and S.. Wright “Optimization for machine learning” MIT Press, 2012
  * [124] Mingxing Tan and Quoc Le “Efficientnet: Rethinking model scaling for convolutional neural networks” In _ICML_ , 2019, pp. 6105–6114
  * [125] R. Tibshirani, J. Bien, J. Friedman, T. Hastie, N. Simon, J. Taylor and R.. Tibshirani “Strong rules for discarding predictors in lasso-type problems” In _J. R. Stat. Soc. Ser. B Stat. Methodol._ 74.2 Wiley Online Library, 2012, pp. 245–266
  * [126] Robert Tibshirani “Regression shrinkage and selection via the lasso” In _J. R. Stat. Soc. Ser. B Stat. Methodol._ 58.1 Wiley Online Library, 1996, pp. 267–288
  * [127] Ryan J. Tibshirani “Adaptive Piecewise Polynomial Estimation via Trend Filtering” In _The Annals of Statistics_ 42.1, 2014
  * [128] Ryan J. Tibshirani and Jonathan Taylor “The solution path of the generalized lasso” In _Ann. Statist._ 39.3 Institute of Mathematical Statistics, 2011, pp. 1335–1371
  * [129] P. Tseng “Dual coordinate ascent methods for non-strictly convex minimization” In _Math. Program._ 59.1 Springer, 1993, pp. 231–247
  * [130] J. Vanschoren, J. Rijn, B. Bischl and L. Torgo “OpenML: networked science in machine learning” In _SIGKDD Explorations_ 15.2 ACM, 2013, pp. 49–60
  * [131] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Millman, Nikolay Mayorov, Andrew R.. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E.. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt and SciPy 1.0 Contributors “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python” In _Nature Methods_ 17, 2020, pp. 261–272
  * [132] R. Wightman, H. Touvron and H. Jégou “ResNet strikes back: An improved training procedure in timm”, 2021, pp. 1–22 arXiv:2110.00476
  * [133] S. Wright and J. Nocedal “Numerical Optimization” Science Springer, 1999
  * [134] Saining Xie, Ross Girshick, Piotr Dollár, Zhuowen Tu and Kaiming He “Aggregated residual transformations for deep neural networks” In _CVPR_ , 2017, pp. 1492–1500
  * [135] . Zhang “Nearly unbiased variable selection under minimax concave penalty” In _Ann. Statist._ 38.2, 2010, pp. 894–942
  * [136] Michael Zhang, James Lucas, Jimmy Ba and Geoffrey E Hinton “Lookahead optimizer: k steps forward, 1 step back” In _Advances in Neural Information Processing Systems_ 32, 2019
  * [137] T. Zhang “Analysis of multi-stage convex relaxation for sparse regularization” In _J. Mach. Learn. Res._ 11.Mar, 2010, pp. 1081–1107
  * [138] Yaowei Zheng, Richong Zhang and Yongyi Mao “Regularizing neural networks via adversarial model perturbation” In _CVPR_ , 2021, pp. 8156–8165

##  Appendix A Contribution from the different authors

Here, we list the different contributions of the authors to the paper:

The original design and ideas where proposed and discussed by TM, AG, JS and
TDLT. The lead developers of the library are TM, MM and AG. Key features have
been added by TDLT and TL. Other contributors includes BaM, BC, BeM, JS, GD,
JL and MD.

The $\ell_{2}$ logistic regression benchmark was coordinated by TM with
contributions from AG, BN, MM, TDLT, TL. The Lasso benchmark was coordinated
by MM, with contributions from AG, BN, CFD, JS, PAB, QK, TM, TL, ZR, and JL.
The ResNet benchmark was coordinated by TM and ZR and most of the
implementation was done by ZR with contributions from PA. The $\ell_{1}$
logistic regression benchmark was coordinated by MM, with contributions from
AG, TL and TM. The MCP benchmark was coordinated by MM, with contributions
from AR, BaM, JL, and JS. The TV-1D benchmark was coordinated by TM, with
contributions from AG, EL and SV. The Ridge regression benchmark was
coordinated by TM and most of the implementation was done by TL with
contributions from CFD and MM.

All authors were involved in the writing of the paper.

##  Appendix B Acknowledgements

It can not be stressed enough how much the Benchopt library relies on
contributions from the community and in particular the Python open source
ecosystem. The pipeline presented in Figure 1 could not exist without the
libraries mentioned hereafter.

The command line interface and API use the
click555https://github.com/pallets/click,
pyyaml666https://github.com/yaml/pyyaml and psutil [114] libraries.

Numerical computations involve numpy [60] and scipy [131]. For cross-language
processing, we used rpy2777https://rpy2.github.io for calling R [110]
libraries and PyJulia888https://pyjulia.readthedocs.io for interfacing with
Julia [12]. The benchmark runs extensively use joblib
999https://github.com/joblib/joblib, loky [96] and
submitit101010https://github.com/facebookincubator/submitit for
parallelization.

The results are stored and processed for visualizations using pandas [95],
matplotlib [65] for static rendering, mako
111111https://github.com/sqlalchemy/mako/ and plotly [67] for interactive
webpages. The participative results website relies partially on
pygithub121212https://github.com/PyGithub.

Our documentation is generated by multiple sphinx-based [22] libraries
(sphinx-bootstrap-theme131313https://github.com/ryan-roemer/sphinx-bootstrap-
theme, sphinx-click141414https://github.com/click-contrib/sphinx-click,
sphinx-gallery [97] and sphinx-prompt151515https://github.com/sbrunner/sphinx-
prompt), and also the numpydoc161616https://github.com/numpy/numpydoc and
pillow [33] libraries.

MM gratefully acknowledges the support of the Centre Blaise Pascal’s IT test
platform at ENS de Lyon (Lyon, France) for Machine Learning facilities. The
platform operates the SIDUS solution [109]. This work was also granted access
to the HPC resources of IDRIS under the allocation 2022-AD011011172R2 and
2022-AD011013570 made by GENCI.

TL, CFD and JS contributions were supported by the Chaire IA CaMeLOt
(ANR-20-CHIA-0001-01).

AG, EL and TM contributions were supported by the Chaire IA ANR BrAIN
(ANR-20-CHIA-0016).

BM contributions were supported by a grant from Digiteo France.

MD contributions were supported by DATAIA convergence institute as part of the
“Programme d’Investissement d’Avenir”, (ANR-17-CONV-0003) operated by Inria.

BN work was supported by the Télécom Paris’s Chaire DSAIDIS (Data Science &
Artificial Intelligence for Digitalized Industry Services).

##  Appendix C A complete Benchmark example: Objective, Dataset and Solver
classes for Ridge regression

Here, we provide code examples for a simple benchmark on Ridge regression. The
Ridge regression – also called $\ell_{2}$-regularized least-squares or
Tikhonov regression – is a popular method to solve least-square problems in
the presence of noisy observations or correlated features. The problem reads:

$\min_{\theta}\frac{1}{2}\|y-X\theta\|_{2}^{2}+\frac{\lambda}{2}\|\theta\|_{2}^{2}\enspace,$
(1)

where $X\in\mathbb{R}^{n\times p}$ is a design matrix, $y\in\mathbb{R}^{n}$ is
the target vector and $\lambda$ is the regularization parameter. This problem
is strongly convex and many methods can be used to solve it. Direct
computation of the close form solution $\theta^{*}=(X^{\top}X+\lambda
Id)^{-1}X^{\top}y$ can be obtained using matrix factorization methods such as
Cholesky decomposition or the SVD [107] or iterative linear solver such as
Conjugate-Gradient [85]. One can also resort on first order methods such as
gradient descent, coordinate descent (known as the Gauss-Seidel method in this
context), or their stochastic variant.

The code for the benchmark is available at
https://github.com/benchopt/benchmark_ridge/. The following code snippets are
provided in the documentation as a template for new benchmarks.

###  C.1 Objective class

The Objective class is the central part of the benchmark, defining the
objective function. This class allows us to monitor the quantities of interest
along the iterations of the solvers, amongst which the objective function
value. An Objective class should define 3 methods:

  * •

set_data(**data): it allows specifying the nature of the data used in the
benchmark. The data is passed as a dictionary of Python variables, so no
constraint is enforced to what can be passed here.

  * •

compute(theta): it allows evaluating the objective function for a given value
of the iterate, here called $\theta$. This method should take only one
parameter, the output returned by the Solver. All other parameters should be
stored in the class with the set_data method. The compute function should
return a float (understood as the objective value) or a dictionary. If a
dictionary is returned it should contain a key called value (the objective
value) and all other keys should correspond to float values allowing tracking
more than one quantity of interest (e.g. train and test errors).

  * •

to_dict(): a method that returns a dictionary to be passed to the
set_objective() method of a Solver.

An Objective class needs to inherit from a base class, benchopt.BaseObjective.
Below is the implementation of the Ridge regression Objective class.

from benchopt import BaseObjective

class Objective(BaseObjective): name = "Ridge regression" parameters = "reg":
[0.1, 1, 10]

def __init__(self, reg=1): self.reg = reg

def set_data(self, X, y): self.X, self.y = X, y

def compute(self, theta): res = self.y - self.X @ theta return .5 * res @ res
+ 0.5 * self.reg * theta @ theta

def to_dict(self): return dict(X=self.X, y=self.y, reg=self.reg)

###  C.2 Dataset class

A Dataset class defines data to be passed to the Objective. More specifically,
a Dataset class should implement one method:

  * •

get_data(): A method outputting a dictionary that can be passed as keyword
arguments **data to the set_data method of the Objective.

A Dataset class also needs to inherit from a base class, benchopt.BaseDataset.

If a Dataset requires some packages to function, Benchopt allows listing some
requirements. The necessary packages should be available via conda or pip.

Below is an example of a Dataset definition using the libsvmdata library,
which exposes datasets from libsvm, such as _leukemia_ , _bodyfat_ and
_gisette_ – described in Table C.1.

from benchopt import BaseDataset from benchopt import safe_import_context

# This context allow to manipulate the Dataset object even if # libsvmdata is
not installed. It is used in ‘benchopt install‘. with safe_import_context() as
import_ctx: from libsvmdata import fetch_libsvm

class Dataset(BaseDataset): name = "libsvm" install_cmd = "conda" requirements
= ["libsvmdata"] parameters = "dataset": ["bodyfat", "leukemia", "gisette"]

def __init__(self, dataset="bodyfat"): self.dataset = dataset

def get_data(self): X, y = fetch_libsvm(self.dataset) return dict(X=self.X,
y=self.y)

###  C.3 Solver class

A Solver class must define three methods:

  * •

set_objective(**objective_dict): This method will be called with the
dictionary objective_dict returned by the method to_dict from the Objective.
The goal of this method is to provide all necessary information to the Solver
so it can optimize the objective function.

  * •

run(stop_value): This method takes only one parameter that controls the
stopping condition of the Solver. Typically this is either a number of
iterations n_iter or a tolerance parameter tol. Alternatively, a callback
function that will be called at each iteration can be passed. The callback
should return False once the computation should stop. The parameter stop_value
is controlled by the stopping_strategy, see below for details.

  * •

get_result(): This method returns a variable that can be passed to the compute
method from the Objective. This is the output of the Solver.

If a Python Solver requires some packages such as scikit-learn, Benchopt
allows listing some requirements. The necessary packages must be available via
conda or pip.

Below is a simple Solver example using scikit-learn implementation of Ridge
regression with different optimization algorithms.

from benchopt import BaseSolver from benchopt import safe_import_context

# This context allow to manipulate the Solver object even if # scikit-learn is
not installed. It is used in ‘benchopt install‘. with safe_import_context() as
import_ctx: from sklearn.linear_model import Ridge

class Solver(BaseSolver): name = "scikit-learn" install_cmd = "conda"
requirements = ["scikit-learn"] parameters = "alg": ["svd", "cholesky",
"lsqr", "sparse_cg", "saga"],

def __init__(self, alg="svd"): self.alg = alg

def set_objective(self, X, y, reg=1): self.X, self.y = X, y self.clf = Ridge(
fit_intercept=False, alpha=reg, solver=self.alg, tol=1e-10 )

def run(self, n_iter): self.clf.max_iter = n_iter + 1 self.clf.fit(self.X,
self.y)

def get_result(self): return self.clf.coef_

###  C.4 Results from the benchmark

Descriptions of datasets Table C.1 describes the datasets used in this
benchmarks.

Table C.1: List of the datasets used in Ridge regression in Appendix C Datasets | References | Samples ($\mathbf{n}$) | Features ($\mathbf{p}$)
---|---|---|---
_leukemia_ | [55] | $38$ | $7129$
_bodyfat_ | [59] | $252$ | $8$
_gisette_ | [59] | $6000$ | $5000$

We also run the solvers on the simulated data described bellow.

Generation process for simulated dataset We generate a linear regression
scenario with decaying correlation for the design matrix, i.e., the ground-
truth covariance matrix is a Toeplitz matrix, with each element
$\Sigma_{ij}=\rho^{|i-j|}$. As a consequence, the generated features have 0
mean, a variance of 1, and the correlation structure as:

$\mathbb{E}[X_{i}]=0~{},\quad\mathbb{E}[X_{i}^{2}]=1\quad\text{and}\quad\mathbb{E}[X_{i}X_{j}]=\rho^{|i-j|}\
.$ (2)

Our simulation scheme also includes the parameter $\texttt{density}=0.2$ that
controls the proportion of non-zero elements in $\theta^{*}$. The target
vector is generated according to linear relationship with Gaussian noise:

$y=X\theta^{*}+\varepsilon\ ,$

such that the signal-to-noise ratio is
$\texttt{snr}=\frac{\|X\theta^{*}\|_{2}}{\|\varepsilon\|_{2}}$.

We use a signal-to-noise ratio $\texttt{snr}=3$, a correlation $\rho$ of $0$
or $0.6$ with $n=500$ samples and $p=1000$ features.

Description of the solvers Table C.2 describes the different solvers compared
in this benchmark.

Table C.2: List of solvers used in the Ridge benchmark in Appendix C Solver | References |  Short Description | Language
---|---|---|---
GD | [20] |  Gradient Descent | Python
Accelerated GD | [99] |  Gradient Descent + acceleration | Python
scikit-learn[svd] | [103] |  SVD (Singular Value Decomposition) |  Python (Cython)
scikit-learn[cholesky] | [103] |  Cholesky decomposition |  Python (Cython)
scikit-learn[lsqr] | [103] |  regularized least-squares |  Python (Cython)
scikit-learn[saga] | [103] |  SAGA (Varianced reduced stochastic method) |  Python (Cython)
scikit-learn[cg] | [103] |  Conjugate Gradient |  Python (Cython)
CD | [11] |  Cyclic Coordinate Descent |  Python (Numba)
lightning[cd] | [14] |  Cyclic Coordinate Descent |  Python (Cython)
snapML[cpu] | [43] |  CD |  Python, C++
snapML[gpu] | [43] |  CD + GPU |  Python, C++

Results Figure D.1 presents the performance of the different methods for
different values of the regularization parameter in the benchmark. The
algorithms based on the direct computation of the closed-form solution
outperform iterative ones in a majority of presented datasets. Among closed-
form algorithms, the Cholesky solver converges faster.

Figure C.1:  Benchmark for the Ridge regression, on 10 solvers, 5 datasets
(_rows_), and 3 variants of the Objective (_columns_) each with a different
regularization value $\lambda\in\\{0.01,0.1,1\\}$. The curves display the
suboptimality of the iterates, $f(\theta^{t})-f(\theta^{*})$, as a function of
time.

###  C.5 Estimating $\theta^{*}$ for convex problems

When the problem is convex, many solvers are guaranteed to converge to a
global minimizer $\theta^{*}$ of the objective function $f$. To estimate
$\theta^{*}$ and $f(\theta^{*})$, Benchopt approximates $\theta^{*}$ by the
iterate achieving the lowest objective_value among all solvers for a given
Dataset and Objective. This means that the sub-optimality plot proposed by
Benchopt are only valid if at least one solver has converged to the optimal
solution. Else, the curves are a lower bound estimate of the sub-optimality.
In practice, for most considered convex problems, running the Solver for long
enough ensures that $f(\theta^{*})$ is correctly estimated.

##  Appendix D $\ell_{2}$-regularized logistic regression

###  D.1 List of solvers and datasets used in the benchmark in Section 3

Table D.1 and Table D.2 respectively present the Solvers and Datasets used in
this benchmark.

Table D.1: List of solvers used in the $\ell_{2}$-regularized logistic regression benchmark in Section 3 Solver |  References |  Short Description | Language
---|---|---|---
lightning[sag] |  [14] |  SAG |  Python (Cython)
lightning[saga] |  [14] |  SAGA |  Python (Cython)
lightning[cd] |  [14] |  Cyclic Coordinate Descent |  Python (Cython)
Tick[svrg] |  [3] |  Stochastic Variance Reduced Gradient |  Python, C++
scikit-learn[sgd] |  [103] |  Stochastic Gradient Descent |  Python (Cython)
scikit-learn[sag] |  [103] |  SAG |  Python (Cython)
scikit-learn[saga] |  [103] |  SAGA |  Python (Cython)
scikit-learn[liblinear] |  [103], [49] |  Truncated Newton Conjugate-Gradient |  Python (Cython)
scikit-learn[lbfgs] |  [103], [131] |  L-BFGS (Quasi-Newton Method) |  Python (Cython)
scikit-learn[newton-cg] |  [103], [131] |  Truncated Newton Conjugate-Gradient |  Python (Cython)
snapml[cpu] |  [43] |  CD |  Python, C++
snapml[gpu] |  [43] |  CD + GPU |  Python, C++
cuML[gpu] |  [112] |  L-BFGS + GPU |  Python, C++
Table D.2: List of the datasets used in $\ell_{2}$-regularized logistic regression in Section 3 Datasets | References | Samples ($\mathbf{n}$) | Features ($\mathbf{p}$) | Density
---|---|---|---|---
_ijcnn1_ | [108] | $49\,990$ | $22$ | $0.045$
_madelon_ | [59] | $2000$ | $500$ | $0.0020$
_news20.binary_ | [72] | $19\,996$ | $1\,355\,191$ | $0.000\,34$
_criteo_ | [38] | $45\,840\,617$ | $1\,000\,000$ | $0.000\,039$

###  D.2 Results

Figure D.1 presents the performance results for the different solvers on the
different datasets using various regularization parameter values, on unscaled
raw data. We observe that when the regularization parameter $\lambda$
increases, the problem tends to become easier and faster to solve for most
methods. Also, the relative order of the method does not change significantly
for the considered range of regularization.

Figure D.1:  Additional benchmark for the $\ell_{2}$-regularized logistic
regression on variants of the Objective (_columns_) with fit_intercept=False.
The curves display the suboptimality of the iterates,
$f(\theta^{t})-f(\theta^{*})$, as a function of time. The columns correspond
to the objective detailed in Section 3 with different value of $\lambda$:
(_first_) $\lambda=0.1$, (_second_) $\lambda=1$ and (_third_) $\lambda=10$.

##  Appendix E Lasso

###  E.1 List of solvers and datasets used in the Lasso benchmark in Section 4

Table E.1 and Table E.2 respectively present the Solvers and Datasets used in
this benchmark.

Table E.1: List of solvers used in the Lasso benchmark in Section 4 Solver |  References |  Short Description | Language
---|---|---|---
blitz |  [70] |  CD + working set |  Python, C++
coordinate descent |  [53] |  (Cyclic) Minimization along coordinates |  Python (Numba)
celer |  [92] |  CD + working set + dual extrapolation |  Python (Cython)
cuML[cd] |  [112] |  (Cyclic) Minimization along coordinates |  Python, C++
cuML[qn] |  [112] |  Orthant-Wise Limited Memory Quasi-Newton (OWL-QN) |  Python, C++
FISTA |  [7] |  ISTA + acceleration | Python
glmnet |  [53] |  CD + working set + strong rules |  R, C++
ISTA |  [39] |  ISTA (Proximal GD) | Python
LARS |  [44] |  Least-Angle Regression algorithm (LARS) |  Python (Cython)
FISTA[adaptive-1] |  [84, Algo 4 ], [50] |  FISTA + adaptive restart | Python
FISTA[greedy] |  [84, Algo 5 ], [50] |  FISTA + greedy restart | Python
noncvx-pro |  [106] |  Bilevel optim + L-BFGS |  Python (Cython)
skglm |  [10] |  CD + working set + primal extrapolation |  Python (Numba)
scikit-learn |  [103] |  CD |  Python (Cython)
snapML[gpu] |  [43] |  CD + GPU |  Python, C++
snapML[cpu] |  [43] |  CD |  Python, C++
lasso.jl |  [78] |  CD | Julia
Table E.2: List of datasets used in the Lasso benchmark in Section 4 Dataset | References | Samples ($\mathbf{n}$) |  Features ($\mathbf{p}$) | Density
---|---|---|---|---
_MEG_ | [57] | $305$ | $7498$ | $1.0$
_news20_ | [72] | $19\,996$ | $1\,355\,191$ | $0.000\,34$
_rcv1_ | [83] | $20\,242$ | $47\,236$ | $0.0036$
_MillionSong_ | [9] | $463\,715$ | $90$ | $1$

###  E.2 Support identification speed benchmark

Since the Lasso is massively used for its feature selection properties, the
speed at which the solvers identify the support of the solution is also an
important performance measure. To evaluate the behavior of solvers in this
task, it is sufficient to add a single new variable in the Objective, namely
the $\ell_{0}$ pseudonorm of the iterate, allowing to produce Figure E.1 in
addition to Figure 4.

Figure E.1: Additional benchmark for the Lasso on variants of the Objective
(_columns_). The curves display the fraction of non-zero coefficients in
iterates $\theta_{t}$ ($\|\theta_{t}\|_{0}/p$), as a function of time.

##  Appendix F ResNet18

###  F.1 Description of the benchmark

Setting up the benchmark The three currently supported frameworks are
TensorFlow/Keras [1, 32], PyTorch [101] and PyTorch Lightning [47]. We report
here results for TensorFlow/Keras and PyTorch. To guarantee that the model
behaves consistently across the different considered frameworks, we
implemented several consistency unit tests. We followed the best practice of
each framework to make sure to achieve the optimal computational efficiency.
In particular, we tried as much as possible to use official code from the
frameworks, and not third-party code. We also optimized and profiled the data
pipelines to make sure that our training was not IO-bound. Our benchmarks were
run using TensorFlow version 2.8 and PyTorch version 1.10.

Descriptions of the datasets

Table F.1: Description of the datasets used in the ResNet18 image classification benchmark Dataset | Content | References | Classes | Train Size | Val. Size | Test Size | Image Size | RGB
---|---|---|---|---|---|---|---|---
_CIFAR-10_ | natural images | [79] | 10 |  $40$k |  $10$k |  $10$k | 32 | ✓
_SVHN_ | digits in natural images | [100] | 10 |  $58.6$k |  $14.6$k |  $26$k | 32 | ✓
_MNIST_ | handwritten digits | [82] | 10 |  $50$k |  $10$k |  $10$k | 28 | ✗

In Table F.1, we present some characteristics of the different datasets used
for the ResNet18 benchmark. In particular, we specify the size of each splits
when using the train-validation-test split strategy. The test split is always
fixed, and is the official one.

While the datasets are downloaded and preprocessed using the official
implementations of the frameworks, we made sure to test that they matched
using a unit test.

ResNet The ResNet18 is the smallest variant of the architecture introduced by
[63]. It consists in 3 stages:

  1. 1.

A feature extension convolution that goes from 3 channels (RGB, or a repeated
grayscale channel in the _MNIST_ case) to 64, followed by a batch
normalization and a ReLU.

  2. 2.

A series of residual blocks. Residual blocks are grouped by scale, and each
individual group starts with a strided convolution to reduce the image scale
(except the first one). As the scale increases, so does the number of features
(64, 128, 256, 512). In the ResNet18 case, each scale group has two individual
residual blocks and there are four scales. A residual block is comprised of
three convolution operations, all followed by a batch normalization layer, and
the first two also followed by a ReLU. The input is then added to the output
of the third batch normalization layer before being fed to a ReLU.

  3. 3.

A classification head that performs global average pooling, before applying a
fully connected (i.e. dense) layer to obtain logits.

Training hyperparameters

Table F.2: Hyperparameters used for each solver. If a hyperparameter’s value is not specified in the table, it was set as the default of the implementation (checked to be consistent across frameworks). Hyperparameter | SGD | Adam
---|---|---
Learning Rate | 0.1 | 0.001
Momentum | 0.9 | N/A
Weight Decay | $5\times 10^{-4}$ | 0.02
Batch Size | 128 | 128

In Table F.2, we specify the hyperparameters we used for the benchmark. For
the SGD, the values were taken from the pytorch-cifar GitHub repository, while
for Adam we took the most relevant ones from the work of [132].

###  F.2 Aligning TensorFlow and PyTorch ResNet18 training

Table F.3:  Differences in off-the-shelf implementations of various components when training ResNet18 for image classification in PyTorch and TensorFlow. The selected versions are put in bold font for components that we were able to reconcile. This highlights the numerous details to consider when comparing experimental results. Component | PyTorch | TensorFlow/Keras
---|---|---
Bias in convolutions | ✗ | ✓
Decoupled weight decay scaling | Multiplied by learning rate | Completely decoupled
Batch normalization momentum | 0.9 | 0.99
Conv2D weights init. | Fan out, normal | Fan average, uniform
Classification head init. (weights) | Fan in, uniform | Fan average, uniform
Classification head init. (bias) | Fan in, uniform | Zeros
Striding in convolutions | Starts one off | Ends one off
Variance estimation in batch norm | unbiased (eval)/biased (training) | biased

We summarized in Table F.3 the different elements that have to be considered
to align the training of a ResNet18 in PyTorch and TensorFlow. Let us detail
here some lines of this table:

  * •

Bias in convolutions: It can be seen in TensorFlow/Keras official
implementation, that convolutions operations use a bias. This is in contrast
to PyTorch’s official implementation in torchvision which does not. Since the
convolutions are followed by batch normalization layers, with a mean removal,
the convolutions’ bias is a spurious parameter, as was noted by [68]. We
therefore chose to use unbiased convolutions.

  * •

Decoupled weight decay scaling: this led us to scale manually the weight decay
used in TensorFlow by the learning rate when setting it. Moreover, because the
weight decay is completely decoupled from the learning rate, it is important
to update it accordingly when using a learning rate schedule, as noted in the
TensorFlow documentation.

  * •

Batch normalization momentum: an important note here is that the convention
used to implement the batch normalization momentum is not the same in the 2
frameworks. Indeed we have the relationship
$\text{{momentum}}_{\text{TF}}=1-\text{{momentum}}_{\text{PT}}$.

  * •

Conv2D weights intialization: TensorFlow/Keras uses the default intialization
which is a Glorot uniform intialization [54]. PyTorch uses a He normal
initialization [62]. We used TensorFlow’s Variance Scaling framework to
differentiate the 2.

  * •

Striding in convolutions: when using a stride of 2 in convolutions on an even-
size image, one needs to specify where to start the convolution in order to
know which lines (one in every two) in the image will be removed. The decision
is different between TensorFlow and PyTorch. This is not expected to have an
effect on the final performance, but it makes it more difficult to compare the
architectures when unit testing. We therefore decided to align the models on
this aspect as well.

  * •

Variance estimation in batch normalization: in order to estimate the batch
variance during training for batch normalization layers, it is possible to
chose between the unbiased and the biased variance estimator. The unbiased
variance estimator applies a Bessel correction to the biased variance
estimator, namely a multiplication by a factor $\frac{m}{m-1}$, where $m$ is
the number of samples used to estimate. It is to be noted that PyTorch does
uses the biased estimator in training, but stores the unbiased estimator for
use during inference. TensorFlow does not allow for such a behaviour, and the
2 are therefore not reconcilable171717It is possible to use the unbiased
estimator in TensorFlow for the batch normalization, even if not documented,
but its application is consistent between training and inference unlike
PyTorch.. Arguably this inconsistency should not play a big role with large
batch sizes, but can be significant for smaller batches, especially in deeper
layers where the feature map size (and therefore the number of samples used to
compute the estimates) is reduced.

Adapting official ResNet implementations to small images In addition to these
elements, it is important to adapt the reference implementations of both
frameworks to the small image case. Indeed, for the case of ImageNet, the
ResNet applies two downsampling operations (a stride-2 convolution and a max
pooling) at the very beginning of the network to make the feature maps size
more manageable. In the case of smaller images, it is necessary to do without
these downsampling operations (i.e. perform the convolution with stride 1 and
get rid of the max pooling).

Coupled weight decay in TensorFlow In TensorFlow, the SGD implementation does
not allow the setting of coupled weight decay. Rather, one has to rely on the
equivalence (up to a scale factor of 2) between coupled weight decay and L2
regularization. However, in TensorFlow/Keras, adding L2 regularization on an
already built model (which is the case for the official ResNet
implementation), is not straightforward and we relied on the workaround of
[121].

###  F.3 VGG benchmark on CIFAR-10

In order to show how flexible Benchopt is, we also ran a smaller version of
the ResNet benchmark using a VGG16 [122] network instead of a ResNet18. In the
Benchopt framework, this amounts to specifying a different model in the
objective, while all the other pieces of code in the benchmark remain
unchanged. Note that the VGG official implementations also need to be adapted
to the CIFAR-10 case by changing the classification head. This was not
specified in the original paper, where no experiment was conducted on small-
scale datasets, and we relied on available open source implementations
(cifar10-vgg16 and cifar-vgg) to make this decision. Importantly, these
implementations use batch normalization to make the training of VGG more
robust to initialization, which is not the case in the official framework
implementations.

Figure F.1: VGG16 image classification benchmark with PyTorch solvers. The
best SGD configuration features data augmentation, momentum, step learning
rate schedule and weight decay.

In Figure F.1, we see that for the case of VGG, the application of weight
decay is so important that without it, in cases with momentum, the model does
not converge.

##  Appendix G $\ell_{1}$-regularized logistic regression

This additional benchmark is dedicated to $\ell_{1}$-regularized logistic
regression, in the same setting as Section 3 but this time with an
$\ell_{1}$-regularization for the parameters of the model:

$\theta^{*}=\mathop{\mathrm{arg\,min}}_{\theta\in\mathbb{R}^{p}}\sum_{i=1}^{n}\log\big{(}1+\exp(-y_{i}X_{i}^{\top}\theta)\big{)}+\lambda\|\theta\|_{1}\enspace.$
(2)

###  G.1 List of solvers and datasets used in the $\ell_{1}$-regularized
logistic regression benchmark

The code for the benchmark is available at
https://github.com/benchopt/benchmark_logreg_l1/. Table G.1 and Table G.2
present the solvers and datasets used in this benchmark.

Table G.1: List of solvers used in the $\ell_{1}$-regularized logistic regression benchmark Solver |  References |  Short Description | Language
---|---|---|---
blitz |  [70] |  CD + working set |  Python, C++
coordinate descent |  [53] |  (Cyclic) Minimization along coordinates |  Python (Numba)
coordinate descent (Newton) |  [53] |  CD + Newton |  Python (Numba)
celer |  [92] |  CD + working set + dual extrapolation |  Python (Cython)
copt[FISTA line search] |  [102], [7] |  FISTA (ISTA + acceleration) + line search |  Python (Cython)
copt[PGD] |  [102, 34] |  Proximal Gradient Descent |  Python (Cython)
copt[PGD linesearch] |  [102, 34] |  Proximal Gradient Descent + linesearch |  Python (Cython)
copt[saga] |  [102] |  SAGA (Variance reduced stochastic method) |  Python (Cython)
copt[svrg] |  [102] |  SVRG (Variance reduced stochastic method) |  Python (Cython)
cuML[gpu] |  [112] |  L-BFGS + GPU |  Python, C++
cuML[qn] |  [112] |  Orthant-Wise Limited Memory Quasi-Newton (OWL-QN) |  Python, C++
cyanure |  [91] |  Proximal Minimization by Incremental Surrogate Optimization (MISO) |  C++ \+ Python wrapper
lightning |  [14] |  (Cyclic) Coordinate Descent |  Python (Cython)
scikit-learn[liblinear] |  [103, 49] |  Truncated Newton Conjugate-Gradient |  Python (Cython)
scikit-learn[lbfgs] |  [103, 131] |  L-BFGS (Quasi-Newton Method) |  Python (Cython)
scikit-learn[newton-cg] |  [103, 131] |  Truncated Newton Conjugate-Gradient |  Python (Cython)
snapml[gpu=True] |  [43] |  CD + GPU |  Python, C++
snapml[gpu=False] |  [43] |  CD |  Python, C++
Table G.2: List of the datasets used in the $\ell_{1}$-regularized logistic regression benchmark Datasets | References | Samples ($\mathbf{n}$) | Features ($\mathbf{p}$) | Density
---|---|---|---|---
_gisette_ | [59] | $6000$ | $5000$ | $9.9\text{\times}{10}^{-1}$
_colon-cancer_ | [59] | $62$ | $2000$ | $1.0$
_news20.binary_ | [72] | $19\,996$ | $1\,355\,191$ | $0.000\,34$
_rcv1.binary_ | [59] | $20\,242$ | $19\,959$ | $3.6\text{\times}{10}^{-3}$

###  G.2 Results

The results of the $\ell_{1}$-regularized logistic regression benchmark are in
Figure G.1.

Figure G.1:  Benchmark for the $\ell_{1}$-regularized logistic regression on
variants of the Objective (_columns_). The curves display the suboptimality of
the iterates, $f(\theta^{t})-f(\theta^{*})$, as a function of time. The first
column corresponds to the objective detailed in Appendix G with
$\lambda=0.1\|X^{\top}y\|_{\infty}/2$, the second one with
$\lambda=0.01\|X^{\top}y\|_{\infty}/2$ and the third column with
$\lambda=0.001\|X^{\top}y\|_{\infty}/2$.

##  Appendix H Unidimensional total variation

The use of 1D Total Variation regularization takes its root in the taut-string
algorithm [5] and can be seen as a special case of either the Rudin-Osher-
Fatemi model [116] or the Generalized Lasso [128] for a quadratic data fit
term. It reads

$\theta^{*}\in\mathop{\mathrm{arg\,min}}_{\theta\in\mathbb{R}^{p}}F(y,X\theta)+\lambda\left\|D\theta\right\|_{1}\enspace,$
(2)

where $F$ is a data fidelity term, $X\in\mathbb{R}^{n\times p}$ is a design
matrix with $n$ samples and $p$ features, $y\in\mathbb{R}^{n}$ is the target
vector, $\lambda>0$ is a regularization hyperparameter, and
$D\in\mathbb{R}^{(p-1)\times p}$ is a finite difference operator defined by
$(D\theta)_{k}=\theta_{k+1}-\theta_{k}$ for all $1\leq k\leq p-1$ (it is also
possible to use cyclic differences).

Most often, the data fidelity term is the $\ell^{2}$-loss
$F(y,z)=\tfrac{1}{2}\left\|y-z\right\|_{2}^{2}$, following an additive
Gaussian noise hypothesis. But the data fit term can also account for other
types of noises, such as noises with heavy tails using the Huber loss
$F(y,z)=|y-z|_{\mu}$ where $|\cdot|_{\mu}$ is defined coordinate-wise by

$|t|_{\mu}=\begin{cases}\frac{1}{2}t^{2}&\textrm{if }|t|\leq\mu\\\
\mu|t|-\frac{\mu^{2}}{2}&\textrm{otherwise.}\end{cases}$

Appendix H promotes piecewise-constant solutions – alternatively said,
solutions such that their gradients is sparse – and was proved to be useful in
several applications, in particular for change point detection [13, 127], for
BOLD signal deconvolution in functional MRI [71, 31] or for detrending in
oculomotor recordings [81].

The penalty $\theta\mapsto\left\|D\theta\right\|_{1}$ is convex but non-
smooth, and its proximity operator has no closed form. Yet as demonstrated by
[36], the taut-string algorithm allows to compute this proximity operator in
$O(p^{2})$ operations in the worst case, but it enjoys a $O(p)$ complexity in
most cases. Other methods do not rely on this proximity operator and directly
solve Appendix H, using either primal-dual approaches [29, 37], or solving the
dual problem [77]. Finally, for 1-dimensional TV regularization, one can also
use the synthesis formulation [46] to solve the problem. By setting
$z=D\theta$ and $\theta=Lz+\rho$ where $L\in\mathbb{R}^{p\times p-1}$ is a
lower trianglar matrix representing an integral operator (cumulative sum), the
problem is equivalent to a Lasso problem, and $\rho^{*}$ has a closed-form
expression (see e.g., [13] for a proof). As a consequence, any lasso solver
can be used to obtain the solution of the Lasso problem $z^{*}$ and the
solution of the original Appendix H $u^{*}$ is retrieved as
$u^{*}=Lz^{*}+\rho^{*}$.

The code for the benchmark is available at
https://github.com/benchopt/benchmark_tv_1d/ and Table H.1 details the
different algorithms used in this benchmark.

Table H.1: List of solvers used in the 1D Total Variation benchmarks Solver | References | Formulation |  Short Description
---|---|---|---
ADMM | [19] | Analysis |  Primal-Dual Augmented Lagrangian
ChambollePock | [29] | Analysis |  Primal-Dual Hybrid Gradient
CondatVu | [37] | Analysis |  Primal-Dual Hybrid Gradient
DPGD | [77] | Analysis |  Dual proximal GD (+ acceleration)
PGD | [36] | Analysis |  Proximal GD + taut–string
| [4] |  |  ProxTV (+ acceleration)
celer | [92] | Synthesis |  CD + working set (lasso)
|  |  |  _only for $\ell_{2}$ data-fit_
FP | [35] | Synthesis |  Fixed point with block updates
ISTA | [39] | Synthesis |  Proximal GD (+ acceleration)
skglm | [10] | Synthesis |  CD + working set

Figure H.1:  Benchmark for the $TV$-regularized regression regression, on 13
solvers, 4 variants of the Objective (_rows_), and 3 configurations for a
simulated dataset (_columns_). The curves display the suboptimality of the
iterates, $f(\theta^{t})-f(\theta^{*})$, as a function of time. The solvers in
this benchmark showcase the three resolution approaches with the Analysis (A),
Dual (D) and Synthesis (S) formulations.

Simulated dataset We use here simulated data, as applications based on fMRI
and EOG signals require access to open and preprocessed data that we will make
available on OpenML [130] in the future. The data are generated as follows: a
block signal $\bar{\theta}\in\mathbb{R}^{p}$ is generated by sampling first a
sparse random vector $z\in\mathbb{R}^{p}$ with $K$ non-zero coefficients
positioned randomly, and taking random values following a $\mathcal{N}(0,1)$
distribution. Finally, $\bar{\theta}$ is obtained by discrete integration as
$\bar{\theta}_{i}=\sum_{k=1}^{i}z_{k}$ for $1\leq i\leq p$. The design matrix
$X\in\mathbb{R}^{n\times p}$ is a Gaussian random design with
$X_{ij}\sim\mathcal{N}(0,1)$. The observations $y$ are obtained as
$y=X\bar{\theta}+\epsilon$, with $\epsilon\sim\mathcal{N}(0,0.01)$ a Gaussian
white noise. For all experiments, we used $p=500$ and vary the number of non-
zero coefficient $K$, and the number of rows $n$ of the matrix $X$.

Results Figure H.1 shows that the solvers using the synthesis formulation and
coordinate descent-based solvers for the Lasso ($\ell_{2}$ data fit term) work
best on this type of problem. For the Huber data fit term, the solver using
the analysis formulation and the taut-string algorithm for the proximal
operator are faster. An interesting observation from this benchmark is the
behavior of the solvers based on primal-dual splitting or dual formulation. We
observe that for all these solvers, the objective starts by increasing. This
is probably due to a sub-optimal initialization of the dual variable compared
to the primal one. While this initialization is seldom described in the
literature, it seems to have a large impact on the speed of these algorithms.
This shows how Benchopt allows to reveal such behavior, and could lead to
practical guidelines on how to select this dual initialization.

Extensions We plan to extend this benchmark in the future to consider higher
dimensional problems – e.g.,  2D TV problems for images – or higher order TV
regularization, such as Total Generalized Variation [23] or inf-convolution of
TV functionals [28] – used of instance for change point detection [127]. Yet,
for 2D or higher dimensional problems, we can no longer use the synthesis
formulation. It is however possible to apply the taut-string method of [36]
and graph-cut methods of [21] and [76] for anisotropic TV, and dual or primal-
dual methods for isotropic, such as Primal-Dual Hybrid Gradient algorithm
[29].

##  Appendix I Linear regression with minimax concave penalty (MCP)

The Lasso problem [126] is a least-squares regression problem with a convex
non-smooth penalty that induces sparsity in its solution. However, despite its
success and large adoption by the machine learning and signal processing
communities, it is plagued with some statistical drawbacks, such as bias for
large coefficients. To overcome these issues, the standard approach is to
consider non-convex sparsity-inducing penalties. Several penalties have been
proposed: _Smoothly Clipped Absolute Deviation_ (SCAD, [48]), the _Log Sum
penalty_ [27], the _capped- $\ell_{1}$ penalty_ [137] or the _Minimax Concave
Penalty_ (MCP, [135]).

This benchmark is devoted to least-squares regression with the latter, namely
the problem:

$\theta^{*}\in\mathop{\mathrm{arg\,min}}_{\theta\in\mathbb{R}^{p}}\tfrac{1}{2n}\left\|y-X\theta\right\|^{2}+\sum_{j=1}^{p}\rho_{\lambda,\gamma}(\theta_{j})\enspace,$
(2)

where $X\in\mathbb{R}^{n\times p}$ is a design matrix containing $p$ features
as columns, $y\in\mathbb{R}^{n}$ is the target vector, and
$\rho_{\lambda,\gamma}$ the penalty function that reads as:

$\rho_{\lambda,\gamma}(t)=\begin{cases}\lambda|t|-\frac{t^{2}}{2\gamma}\enspace,\quad&\text{if
}|t|\leq\gamma\lambda\enspace,\\\
\frac{\lambda^{2}\gamma}{2}\enspace,\quad&\text{if
}|t|>\gamma\lambda\enspace.\end{cases}$

Similarly to the Lasso, Appendix I promotes sparse solutions but the
optimization problem raises some difficulties due to the non-convexity and
non-smoothness of the penalty. Nonetheless, several efficient algorithms have
been derived for solving it. The ones we use in the benchmark are listed in
Table I.1.

Table I.1: List of solvers used in the MCP benchmark Solver |  References | Short Description
---|---|---
CD |  [24], [94] | Proximal coordinate descent
PGD |  [16] | Proximal operator + Gradient descent
GIST |  [56] | Proximal gradient + Barzilai-Borwein rule
WorkSet CD |  [15] | Coordinate descent + working set
skglm |  [10] | Accelerated coordinate descent + Working set

The code for the benchmark is available at
https://github.com/benchopt/benchmark_mcp/. For this benchmark, we run the
solvers on the _colon-cancer_ dataset and on the simulated dataset described
in Table C.1. We use a signal-to-noise ratio $\texttt{snr}=3$, a correlation
$\rho=0.6$ with $n=500$ observations and $p=2000$ features.

###  I.1 Results

The result of the benchmark is presented in Figure I.1 The problem is non-
convex and solvers are only guaranteed to converge to local minima; hence in
Figure I.1 we monitor the distance of the negative gradient to the Fréchet
subdifferential of the MCP, representing the violation of the first order
optimality condition. Other metrics, such as objective of iterates sparsity,
are monitored in the full benchmark, allowing to compare the different limit
points obtained by the solvers.

Figure I.1:  Benchmark for the MCP regression on variants of the Objective
(_columns_). The curves display the violation of optimality conditions,
$\operatorname{dist}(-X^{\top}(X\theta_{t}-y)/n,\partial\rho_{\lambda,\gamma}(\theta_{t}))$,
as a function of time. $\gamma$ is set to 3, and $\lambda$ is parameterized as
a fraction of the Lasso’s $\lambda_{\max}$,
$\left\|X^{\top}y\right\|_{\infty}/n$.