Thomas Moreau

Mathurin Massias

Alexandre Gramfort

Pierre Ablin

Pierre-Antoine Bannier

Benjamin Charlier

Mathieu Dagréou

Tom Dupré La Tour

Ghislain Durif

Cassio F Dantas

Quentin Klopfenstein

Johan Larsson

En Lai

Tanguy Lefort

Benoit Malézieux

Badr Moufad

Binh T Nguyen

Alain Rakotomamonjy

Zaccharie Ramzi

Joseph Salmon

Samuel Vaiter

Benchopt: Reproducible, efficient and collaborative optimization benchmarks

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: ℓ 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.

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 [START_REF] Sculley | Winner's curse? On pace, progress, and empirical rigor[END_REF]. 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 [START_REF] Raff | A step toward quantifying independently reproducible machine learning research[END_REF], particularly because comparisons 36th Conference on Neural Information Processing Systems (NeurIPS 2022). 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 [START_REF] Pineau | ICLR reproducibility challenge 2019[END_REF].

In the past years, the ML community has actively sought to overcome this "reproducibility crisis" [START_REF] Hutson | Artificial intelligence faces reproducibility crisis[END_REF] through collaborative initiatives such as open datasets (OpenML, Vanschoren et al. 2013), standardized code sharing [START_REF] Forde | Reproducible research environments with repo2docker[END_REF], benchmarks (MLPerf, [START_REF] Mattson | MLPerf: An industry standard benchmark suite for machine learning performance[END_REF], the NeurIPS and ICLR reproducibility challenges [START_REF] Pineau | ICLR reproducibility challenge 2019[END_REF][START_REF] Pineau | Improving reproducibility in machine learning research: a report from the NeurIPS 2019 reproducibility program[END_REF] and new journals (e.g., [START_REF] Rougier | ReScience C: a journal for reproducible replications in computational science[END_REF]. 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 [START_REF] Sra | Optimization for machine learning[END_REF]. It is thus crucial to know which methods to use depending on the task and setting [START_REF] Bartz-Beielstein | Benchmarking in optimization: Best practice and open issues[END_REF]. While some papers in optimization for ML provide extensive validations [START_REF] Lueckmann | Benchmarking Simulation-Based Inference[END_REF], 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 [START_REF] Schmidt | Descending through a Crowded Valley -Benchmarking Deep Learning Optimizers[END_REF], 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. 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:

• ℓ 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.

The Benchopt library

The Benchopt library aims to provide a standard toolset and structure to implement benchmarks for optimization in ML, where the problems depend on some input dataset D. The considered problems are of the form θ * ∈ arg min

θ∈Θ f (θ; D, Λ) , ( 1
)
where f is the objective function, Λ are its hyperparameters, and Θ is the feasible set for θ. The scope of the library is to evaluate optimization methods in their wide sense by considering the sequence {θ t } t produced to approximate θ * . We emphasize than Benchopt does not provide a fixed set of benchmarks, but a framework to create, extend and share benchmarks on any problem of the form (1). To provide a flexible and extendable coding standard, benchmarks are defined as the association of three types of object classes: Objective: It defines the function f to be minimized as well as the hyperparameters Λ or the set Θ, 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 θ t .

Datasets:

The Dataset objects provide the data 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 {θ 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. Appendix B presents a concrete example on Ridge regression of how to construct a Benchopt benchmark while additional design design choices of Benchopt are discussed in Appendix C.

For each Dataset and Solver, and for each set of parameters, Benchopt retrieves a sequence {θ t } t and evaluates the metrics defined in the Objective for each θ 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 ensures fair evaluation of baselines in follow-up experiments, as implementations validated by the community are available. Moreover, benchmarks can be extended easily as one can add a Dataset or a Solver to the comparison by adding a single file. Finally, by supporting multiple metrics -e.g., training and testing losses, error on parameter estimates, sparsity of the estimate -the Objective class offers the flexibility to define the concurrent evaluation, which can be extended to track extra metrics on a per benchmark basis, depending on the problem at hand.

As one of the goal of Benchopt is to make benchmarks as simple as possible, it also provides a set of features to make them easy to develop and run. 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). Moreover, benchmarks can be run in parallel with checkpointing of the results, enabling large scale evaluations on many CPU or GPU nodes. Benchopt also makes it possible to run solvers with many different hyperparameters' values , allowing to assess their sensitivity on the method performance. Benchmark results are also 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: ℓ 2 -regularized logistic regression

Logistic regression is a very popular method for binary classification. From a design matrix X ∈ R n×p with rows X i and a vector of labels y ∈ {-1, 1} n with corresponding element y i , ℓ 2 -regularized logistic regression provides a generalized linear model indexed by θ * ∈ R p to discriminate the classes by solving

θ * = arg min θ∈R p n i=1 log 1 + exp(-y i X ⊤ i θ) + λ 2 ∥θ∥ 2 2 , (2)
where λ > 0 is the regularization hyperparameter. Thanks to the regularization part, Problem (2) 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 [START_REF] Wright | Numerical Optimization[END_REF]. 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 [START_REF] Liu | On the Limited Memory BFGS Method for Large Scale Optimization[END_REF] stands out for its small memory footprint, its robustness and fast convergence in a variety of settings. On the other hand, truncated Newton methods [START_REF] Dembo | Inexact Newton Methods[END_REF] try to directly approximate Newton's direction by using e.g., the conjugate gradient method [START_REF] Fletcher | Function Minimization by Conjugate Gradients[END_REF] 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 [START_REF] Bottou | Large-Scale Machine Learning with Stochastic Gradient Descent[END_REF], with Stochastic Gradient Descent (SGD) as a main instance. The core idea is to use cheap and noisy estimates of the gradient [START_REF] Robbins | A stochastic approximation method[END_REF][START_REF] Kiefer | Stochastic estimation of the maximum of a regression function[END_REF]. 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 [START_REF] Schmidt | Minimizing Finite Sums with the Stochastic Average Gradient[END_REF], SAGA [START_REF] Defazio | SAGA: A Fast Incremental Gradient Method With Support for Non-Strongly Convex Composite Objectives[END_REF] and SVRG [START_REF] Johnson | Accelerating Stochastic Gradient Descent Using Predictive Variance Reduction[END_REF]) 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 [START_REF] Bertsekas | Nonlinear programming[END_REF] have also been proposed to solve Problem (2). 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 [START_REF] Dünner | Snap ML: A hierarchical framework for machine learning[END_REF].

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 [START_REF] Koh | An interior-point method for large-scale l1-regularized logistic regression[END_REF] 10 -8 madelon 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 -3 10 -2 10 -1 10 0 10 1 10 2 10

Time (s)

Figure 3: Benchmark for the ℓ 2 -regularized logistic regression, on 13 solvers, 4 datasets (rows), and 3 variants of the Objective (columns) with λ = 1. The curves display the suboptimality of the iterates, f (θ t )f (θ * ), as a function of time. The first column corresponds to the objective function detailed in Problem (2). 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.

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 (θ t )f (θ * ), for three datasets and three scenarios. Here, because the problem is convex, θ * is approximated by the best iterate across all runs (see Section C.1). Overall, the benchmark shows the benefit of using GPU solvers (cuML and snapML), even for small scale tasks such as ijcnn1. Note 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. 1 For news20.binary, the stochastic solvers such as SAG and SAGA have degraded performances on scaled data. Here, the scaling makes the problem harder. 2

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, [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF][START_REF] Chen | Atomic decomposition by basis pursuit[END_REF], is an archetype of non-smooth ML problems, whose impact on ML, statistics and signal processing in the last three decades has been considerable [START_REF] Bühlmann | Statistics for high-dimensional data[END_REF][START_REF] Hastie | Statistical Learning with Sparsity: The Lasso and Generalizations[END_REF]. It consists of solving

θ * ∈ arg min θ∈R p 1 2 ∥y -Xθ∥ 2 + λ ∥θ∥ 1 , (3)
where X ∈ R n×p is a design matrix containing p features as columns, y ∈ R n is the target vector, and λ > 0 is a regularization hyperparameter. The Lasso estimator was popularized for variable selection: when λ is high enough, many entries in θ * are exactly equal to 0. This leads to more interpretable models and reduces overfitting compared to the least-squares estimator.

Solvers for Problem (3) have evolved since its introduction by [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. After generic quadratic program solvers, new dedicated solvers were proposed based on iterative reweighted leastsquares (IRLS) [START_REF] Grandvalet | Least absolute shrinkage is equivalent to quadratic penalization[END_REF], followed by LARS [START_REF] Efron | Least angle regression[END_REF], a homotopy method computing the full Lasso path 3 . The LARS solver helped popularize the Lasso, yet the algorithm suffers from stability issues and can be very slow for worst case situations [START_REF] Mairal | Complexity analysis of the Lasso regularization path[END_REF]. General purpose solvers became popular for Lasso-type problems with the introduction of the iterative soft thresholding algorithm (ISTA, [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF]), an instance of forward-backward splitting [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]. These algorithms became standard in signal and image processing, especially when accelerated (FISTA, [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF].

In parallel, proximal coordinate descent has proven particularly relevant for the Lasso in statistics.

Early theoretical results were proved by [START_REF] Tseng | Dual coordinate ascent methods for non-strictly convex minimization[END_REF] and [START_REF] Sardy | Block coordinate relaxation methods for nonparametric wavelet denoising[END_REF], 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 θ * , trying to identify its support to reduce the problem size. Best performing variants of this scheme are screening rules (e.g., El [START_REF] El Ghaoui | Safe feature elimination in sparse supervised learning[END_REF][START_REF] Bonnefoy | Dynamic screening: accelerating first-order algorithms for the Lasso and Group-Lasso[END_REF][START_REF] Ndiaye | Gap Safe screening rules for sparsity enforcing penalties[END_REF] and working/active sets (e.g., [START_REF] Johnson | Blitz: A Principled Meta-Algorithm for Scaling Sparse Optimization[END_REF][START_REF] Massias | Celer: a fast solver for the lasso with dual extrapolation[END_REF], including strong rules [START_REF] Tibshirani | Strong rules for discarding predictors in lasso-type problems[END_REF].

While reviews of Lasso solvers have already been performed (Bach et al., 2012, 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 λ, controlling the sparsity of the solution, parameterized as a fraction of λ max = X ⊤ y ∞ (the minimal hyperparameter such that θ * = 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 10 -2 10 -1 10 0 10 1 10 2 10 3 rcv1.binary 10 -2 10 -1 10 0 10 1 10 2 10 -2 10 -1 10 0 10 1 10 2 10 -1 10 0 10 1 10 2 10 3 Time (s)

10 3 10 0 10 -3 10 -6
news20.binary 10 -1 10 0 10 1 10 2 10 3

Time (s)

10 -1 10 0 10 1 10 2 10 3

Time (s)

Figure 4: Benchmark for the Lasso, on 17 solvers, 4 datasets (rows), and 3 variants of the Objective (columns) with decreasing regularization λ. The curves display the suboptimality of the objective function, f (θ t )f (θ * ), as a function of time.

Results Figure 4 presents the result of the benchmark on the Lasso, in terms of objective suboptimality f (θ t )f (θ * ) 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 [START_REF] Poon | Smooth Bilevel Programming for Sparse Regularization[END_REF], which highlights the need for extensive, standard experimental setups.

When the support of the solution is small (either small λ, 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. 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 ignore conversion overhead, as well as just-in-time compilation cost. 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 [START_REF] Krizhevsky | Imagenet Classification with Deep Convolutional Neural Networks[END_REF], image classification has become a de facto standard to validate novel methods in the field. Among the different network architectures, ResNets [START_REF] He | Deep residual learning for image recognition[END_REF] are extensively used in the community as they provide strong and versatile baselines [START_REF] Xie | Aggregated residual transformations for deep neural networks[END_REF][START_REF] Tan | Efficientnet: Rethinking model scaling for convolutional neural networks[END_REF][START_REF] Dosovitskiy | An Image is Worth 16x16 Words: Transformers for Image Recognition at Scale[END_REF][START_REF] Brock | High-performance large-scale image recognition without normalization[END_REF][START_REF] Liu | A ConvNet for the 2020s[END_REF]. While many papers present results with such model on classical datasets, with sometimes extensive ablation studies [START_REF] He | Bag of tricks for image classification with convolutional neural networks[END_REF][START_REF] Wightman | ResNet strikes back: An improved training procedure in timm[END_REF][START_REF] Bello | Revisiting ResNets: Improved Training and Scaling Strategies[END_REF][START_REF] Schmidt | Descending through a Crowded Valley -Benchmarking Deep Learning Optimizers[END_REF], 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, Netzer et al. ( 2011 Another notable difference is that we report the evolution of the test error rather than the training loss.

Because we chose to monitor the test 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. This is in contrast to a case where we would monitor the train loss, and therefore make the LR and weight decay schedules, as well as the data augmentation policy, part of the objective. We focus on 2 standard methods: stochastic gradient descent (SGD) with momentum and Adam [START_REF] Kingma | Adam: A Method for Stochastic Optimization[END_REF], as well as a more recently published one: Lookahead [START_REF] Zhang | Lookahead optimizer: k steps forward, 1 step back[END_REF]. The LR schedules are chosen among fixed LR, step LR4 , and cosine annealing [START_REF] Loshchilov | SGDR: Stochastic Gradient Descent with Warm Restarts[END_REF]. We also consider decoupled weight decay for Adam [START_REF] Loshchilov | Decoupled Weight Decay Regularization[END_REF], and coupled weight decay (i.e., ℓ 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, and discuss their selection as well as their sensitivity in Appendix F.

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.3). We conducted an extensive effort to align the networks implementation in different frameworks 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

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. 5 .

Firstly, 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 ResNet18before that). Our ablation study shows 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 [START_REF] Zhang | Lookahead optimizer: k steps forward, 1 step back[END_REF]. 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 (≈ 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.

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. 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.

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.

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. In particular, it could not exist without the libraries mentioned in Appendix A.

This work was granted access to the HPC resources of IDRIS under the allocation 2022-AD011011172R2 and 2022-AD011013570 made by GENCI, which was used to run all the benchmarks. MM also 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 [START_REF] Quemener | SIDUS-the solution for extreme deduplication of an operating system[END_REF]. A Software ecosystem acknowledgement

The command line interface and API use the click, pyyaml and psutil [START_REF] Rodola | Psutil package: a cross-platform library for retrieving information on running processes and system utilization[END_REF] libraries.

Numerical computations involve numpy [START_REF] Harris | Array programming with NumPy[END_REF] and scipy [START_REF] Virtanen | SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python[END_REF]. For cross-language processing, we used rpy2 for calling R (R Core Team, 2017) libraries and PyJulia for interfacing with Julia [START_REF] Bezanson | Julia: A fresh approach to numerical computing[END_REF]. The benchmark runs extensively use joblib, loky [START_REF] Moreau | [END_REF] and submitit for parallelization.

The results are stored and processed for visualizations using pandas [START_REF] Mckinney | Data structures for statistical computing in python[END_REF], matplotlib [START_REF] Hunter | Matplotlib: A 2D graphics environment[END_REF] for static rendering, mako and plotly [START_REF] Inc | Collaborative data science[END_REF] for interactive webpages. The participative results website relies partially on pygithub.

Our documentation is generated by multiple sphinx-based [START_REF] Brandl | Sphinx documentation[END_REF] libraries (sphinx-bootstrap-theme, sphinx-click, sphinx-gallery [START_REF] Nájera | sphinx-gallery/sphinx-gallery: Release v0[END_REF] and sphinx-prompt), and also the numpydoc and pillow [START_REF] Clark | Pillow (PIL Fork) Documentation[END_REF] libraries.

B 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 ℓ 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 θ 1 2 ∥y -Xθ∥ 2 2 + λ 2 ∥θ∥ 2 2 , (4)
where X ∈ R n×p is a design matrix, y ∈ R n is the target vector and λ 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 θ * = (X ⊤ X + λId) -1 X ⊤ y can be obtained using matrix factorization methods such as Cholesky decomposition or the SVD [START_REF] Press | Numerical recipes 3rd edition: The art of scientific computing[END_REF] or iterative linear solver such as Conjugate-Gradient [START_REF] Liu | On the Limited Memory BFGS Method for Large Scale Optimization[END_REF]. 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.

B.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 θ. 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. 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 )

An

B.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 B

B.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. 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 Σ ij = ρ |i-j| . As a consequence, the generated features have 0 mean, a variance of 1, and the correlation structure as:

E[X i ] = 0 , E[X 2 i ] = 1 and E[X i X j ] = ρ |i-j| . (5)
Our simulation scheme also includes the parameter density = 0.2 that controls the proportion of non-zero elements in θ * . The target vector is generated according to linear relationship with Gaussian noise:

y = Xθ * + ε ,
such that the signal-to-noise ratio is snr = ∥Xθ * ∥2 ∥ε∥2 .

We use a signal-to-noise ratio snr = 3, a correlation ρ of 0 or 0.6 with n = 500 samples and p = 1000 features. 10 -8 leukemia 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 4 10 0 10 -4 10 -8 gisette 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 4 10 0 10 -4 10 -8 Simulated 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 Time (s)

Description of the solvers

10 4 10 0 10 -4 10 -8
Simulated[ρ = 0.6] 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1

Time (s)

10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 Time (s)

C Design choices

Benchopt has made some design choices, while trying as much as possible to leave users free of customizing the behavior on each benchmark. We detail the most important ones in this section.

C.1 Estimating θ * for convex problems

When the problem is convex, many solvers are guaranteed to converge to a global minimizer θ * of the objective function f . To estimate θ * and f (θ * ), Benchopt approximates θ * 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 (θ * ) is correctly estimated.

C.2 Stopping solvers

Benchopt offers many ways to stop running a solver. The most common is to stop the solver when the objective value does not decrease significantly between iterations. For some convex problems, we also propose to track the duality gap (which upper bounds the suboptimality), as is done for the Lasso. For non convex problems, criteria such as gradient norm or violation of first order conditions can be used, as users do in practice. These criteria can easily be customized.

C.3 Wall-clock time versus number of iterations

Measuring time or iteration are two alternatives that make sense in their respective contexts. Practitioners mostly care about the time it takes to solve their problem, while researchers in mathematical optimization may want to abstract away the implementation and hardware details and only consider iteration. The benchmarks we have presented showcase efficient implementations and are also interested in hardware and implementation differences (e.g. CPU vs GPU solvers for ??, torch versus tensorflow for Section F.4), hence our focus on time. However, Benchopt does not impose a choice between the two measures: it is perfectly possible to create plots as a function of the number of iterations as evidenced for example in ??.

D ℓ 2 -regularized logistic regression D.1 List of solvers and datasets used in the benchmark in Section 3

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 λ 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.

E Lasso

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

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 ℓ 0 pseudonorm of the iterate, allowing to produce Figure 10 -2 10 -1 10 0 10 1 10 2 10 3 10 0 10 -1 10 -2 10 -3 MEG λ = 0.1λ max 10 -2 10 -1 10 0 10 1 10 2 10 3 λ = 0.01λ max 10 -2 10 -1 10 0 10 1 10 2 10 3 λ = 0.001λ max 10 -2 10 -1 10 0 10 1 10 2 10 0 10 -1 10 -2 10 -3 rcv1.binary 10 -2 10 -1 10 0 10 1 10 2 10 -2 10 -1 10 0 10 1 10 2 10 -1 10 0 10 1 10 2 10 3 10 0 10 -1 10 -2 10 -3 10 -4 10 -5 news20.binary 10 -1 10 0 10 1 10 2 10 3 10 -1 10 0 10 1 10 2 10 3 10 -2 10 -1 10 0 10 1 Time (s)

10 0 10 -1 10 -2
MillionSong 10 -2 10 -1 10 0 10 1

Time (s)

10 -2 10 -1 10 0 10 1

Time (s)

Figure E.1: Additional benchmark for Lasso on variants of the Objective (columns). The curves display the fraction of non-zero coefficients in iterates θ t (∥θ t ∥ 0 /p), as a function of time.

While practitioners are mainly concerned with the time it takes to solve their optimization problem, one may also be interested in the convergence as a function of the number of iterations. This is particularly relevant to compare theoretical convergence rates with experiments. Benchopt natively supports such functionality. Yet, this makes sense only if one iteration of each algorithm costs the same. Figure E.2 presents such a case on the leukemia dataset, using algorithms for which one iteration costs n × p. One can observe that cyclic coordinate descent as implemented in Cython in scikit-learn or in Numba lead to identical results, while they outperform proximal gradient methods.

F ResNet18

F.1 Description of the benchmark Setting up the benchmark The three currently supported frameworks are TensorFlow/Keras [START_REF] Abadi | TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems[END_REF][START_REF] Chollet | [END_REF], PyTorch [START_REF] Paszke | PyTorch: An Imperative Style, High-Performance Deep Learning Library[END_REF] and PyTorch Lightning [START_REF] Falcon | PyTorchLightning/pytorchlightning: 0.7.6 release[END_REF]. 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. 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.

Descriptions of the datasets

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.

The ResNet18 is the smallest variant of the architecture introduced by [START_REF] He | Deep residual learning for image recognition[END_REF]. It consists in 3 stages:

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. 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. A classification head that performs global average pooling, before applying a fully connected (i.e. dense) layer to obtain logits.

Training's hyperparameters F.2. In dashed black is the state of the art for CIFAR-10 with a ResNet18 measured by [START_REF] Zhang | Lookahead optimizer: k steps forward, 1 step back[END_REF]. Curves are exponentially smoothed.

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 [START_REF] Wightman | ResNet strikes back: An improved training procedure in timm[END_REF].

F.2 Hyperparameter sensitivity

In the benchmark presented in Section 5, we consider fixed hyperparameters chosen from common practices to train ResNet18 models for an image classification task. However, in practice, these hyperparameters must be carefully set, either via a grid search, or via more adapted algorithms such as random search (bergstra2012random) or bayesian optimization [START_REF] Paszke | PyTorch: An Imperative Style, High-Performance Deep Learning Library[END_REF]. It is therefore important to evaluate how sensitive an optimizer is to choosing the right parameters, as more sensitive methods will require more exhaustive hyperparameters search. In Figure F.1, we study this issue using Benchopt for image classification on CIFAR-10. Despite achieving the best results in terms of accuracy, SGD is way more sensitive to the choice of hyperparameters than Adam. 6Another way to look at hyperparameter sensitivity is to evaluate how a given selection of hyperparameters performs for different tasks. Figure 5 shows that while SGD is sensitive to the choice of learning rate and weight decay, the selected values work very well across 3 different datasets.

F.3 Aligning TensorFlow and PyTorch ResNet18 training

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 [START_REF] Ioffe | Batch normalization: Accelerating deep network training by reducing internal covariate shift[END_REF]. 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 momentum TF = 1momentum PT . [START_REF] Glorot | Understanding the difficulty of training deep feedforward neural networks[END_REF]. PyTorch uses a He normal initialization [START_REF] He | Delving deep into rectifiers: Surpassing human-level performance on imagenet classification[END_REF]. 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 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 reconcilable7 . 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 Silva (2019).

F.4 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 [START_REF] Simonyan | Very Deep Convolutional Networks for Large-Scale Image Recognition[END_REF] 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.

In Figure F.2, 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.

G ℓ 1 -regularized logistic regression

This additional benchmark is dedicated to ℓ 1 -regularized logistic regression, in the same setting as Problem (2) but this time with an ℓ 1 -regularization for the parameters of the model:

θ * = arg min θ∈R p n i=1 log 1 + exp(-y i X ⊤ i θ) + λ∥θ∥ 1 . (6)
G.1 List of solvers and datasets used in the ℓ 1 -regularized logistic regression benchmark

The code for the benchmark is available at https://github.com/benchopt/benchmark_logreg_ l1/. colon-cancer 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 -2 10 -1 10 0 10 1 10 2 10 3 10 0 10 -3 10 -6

rcv1.binary 10 -2 10 -1 10 0 10 1 10 2 10 -2 10 -1 10 0 10 1 10 2 10 -1 10 0 10 1 10 2 10 3 Time (s) The first column corresponds to the objective detailed in Problem ( 6) with λ = 0.1∥X ⊤ y∥ ∞ /2, the second one with λ = 0.01∥X ⊤ y∥ ∞ /2 and the third column with λ = 0.001∥X ⊤ y∥ ∞ /2.

G.2 Results

The results of the ℓ 1 -regularized logistic regression benchmark are in Huber[µ = 0.9] λ = 0.5λ max 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 Time (s) Huber[µ = 0.9] λ = 0.1λ max 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 Time (s)

10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 Time (s)

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 Vanschoren et al., 2013 in the future. The data are generated as follows: a block signal θ ∈ R p is generated by sampling first a sparse random vector z ∈ R p with K non-zero coefficients positioned randomly, and taking random values following a N (0, 1) distribution. Finally, θ is obtained by discrete integration as θi = i k=1 z k for 1 ≤ i ≤ p. The design matrix X ∈ R n×p is a Gaussian random design with X ij ∼ N (0, 1). The observations y are obtained as y = X θ + ϵ, with ϵ ∼ 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 descentbased solvers for the Lasso (ℓ 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

I Linear regression with minimax concave penalty (MCP)

The Lasso problem [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] 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, [START_REF] Fan | Variable selection via nonconcave penalized likelihood and its oracle properties[END_REF], the Log Sum penalty [START_REF] Candès | Enhancing Sparsity by Reweighted l 1 Minimization[END_REF], the capped-ℓ 1 penalty [START_REF] Zhang | Analysis of multi-stage convex relaxation for sparse regularization[END_REF] or the Minimax Concave Penalty (MCP, Zhang 2010a). This benchmark is devoted to least-squares regression with the latter, namely the problem:

θ * ∈ arg min θ∈R p 1 2n ∥y -Xθ∥ 2 + p j=1 ρ λ,γ (θ j ) , (8)
where X ∈ R n×p is a design matrix containing p features as columns, y ∈ R n is the target vector, and ρ λ,γ the penalty function that reads as:

ρ λ,γ (t) = λ|t| -t 2 2γ , if |t| ≤ γλ , λ 2 γ 2 , if |t| > γλ .
Similarly to the Lasso, Problem (8) 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. 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 B.4. We use a signal-to-noise ratio snr = 3, a correlation ρ = 0.6 with n = 500 observations and p = 2000 features.

I.1 Results

The result of the benchmark is presented in can observe that the CMA and TwoPointsDE implementations from nevergrad consistently reaches the global minimum. In addition, the CMA implementation from optuna is a bit slower than the one from nevergrad. Also one can notice that random search offers reasonable results. The TPE method seems to suffer from the curse of dimensionality, as most kernel methods in non-parametric estimation. Finally regarding the scipy solvers, Powell can be competitive, while Nelder-Mead and BFGS suffer a lot from local minima.

Figure 1 :

 1 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.

Figure

  Figure 2: Standard benchmark structure

  ); MNIST, LeCun et al. (2010) and CIFAR-10, Krizhevsky (2009)-and cross-framework -TensorFlow/Keras, Abadi et al. (2015) and Chollet et al. (2015); PyTorch, Paszke et al. (2019)-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 θ of the model. Contrary to logistic regression and the Lasso, this problem is non-convex due to the non-linearity of the model f θ .

Checklist 1 .

 1 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). BMa contributions were supported by a grant from Digiteo France. MD contributions were supported by a public grant overseen by the French National Research Agency (ANR) through the program UDOPIA, project funded by the ANR-20-THIA-0013-01 and DATAIA convergence institute (ANR-17-CONV-0003). BN work was supported by the Télécom Paris's Chaire DSAIDIS (Data Science & Artificial Intelligence for Digitalized Industry Services). BMo contributions were supported by a grant from the Labex MILYON. For all authors... (a) Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? [Yes] (b) Did you describe the limitations of your work? [Yes] (c) Did you discuss any potential negative societal impacts of your work? [Yes] (d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes] 2. If you are including theoretical results... (a) Did you state the full set of assumptions of all theoretical results? [N/A] (b) Did you include complete proofs of all theoretical results? [N/A] 3. If you ran experiments... (a) Did you include the code, data, and instructions needed to reproduce the main experimental results (either in the supplemental material or as a URL)? [Yes] (b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes] These are specified in Appendix, on a per-benchmark basis. (c) Did you report error bars (e.g., with respect to the random seed after running experiments multiple times)? [No] Error bars are not reported for clarity, but Benchopt allows this in particular in html versions of the plots that can be found in https://benchopt.github.io/results/preprint_results.html. (d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [Yes] 4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets... (a) If your work uses existing assets, did you cite the creators? [Yes] (b) Did you mention the license of the assets? [Yes] In the introduction. (c) Did you include any new assets either in the supplemental material or as a URL? [Yes] (d) Did you discuss whether and how consent was obtained from people whose data you're using/curating? [N/A] (e) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A] 5. If you used crowdsourcing or conducted research with human subjects... (a) Did you include the full text of instructions given to participants and screenshots, if applicable? [N/A] (b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, if applicable? [N/A] (c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]

  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

  Results Figure B.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. 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 λ = 0.1 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2

Figure B. 1 :

 1 Figure B.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 λ ∈ {0.01, 0.1, 1}. The curves display the suboptimality of the iterates, f (θ t )f (θ * ), as a function of time.

Figure D. 1 :

 1 Figure D.1: Additional benchmark for the ℓ 2 -regularized logistic regression on variants of the Objective (columns) with fit_intercept=False. The curves display the suboptimality of the iterates, f (θ t )f (θ * ), as a function of time. The columns correspond to the objective detailed in Problem (2) with different value of λ: (first) λ = 0.1, (second) λ = 1 and (third) λ = 10.

Figure E. 2 :

 2 Figure E.2: Convergence speed with respect to the number of iterations for some solvers of the Lasso benchmark on the leukemia dataset.

Figure F. 1 :

 1 Figure F.1: ResNet18 image classification benchmark on CIFAR-10 for different values of learning rate and weight decay for SGD and Adam. The default values are that reported in TableF.2. In dashed black is the state of the art for CIFAR-10 with a ResNet18 measured by[START_REF] Zhang | Lookahead optimizer: k steps forward, 1 step back[END_REF]. Curves are exponentially smoothed.

Figure G. 1 :

 1 Figure G.1: Benchmark for the ℓ 1 -regularized logistic regression on variants of the Objective (columns). The curves display the suboptimality of the iterates, f (θ t )f (θ * ), as a function of time.The first column corresponds to the objective detailed in Problem (6) with λ = 0.1∥X ⊤ y∥ ∞ /2, the second one with λ = 0.01∥X ⊤ y∥ ∞ /2 and the third column with λ = 0.001∥X ⊤ y∥ ∞ /2.

Figure

  10 -2 10 -1 10 0 10 1 10 2 10 3 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3

Figure H. 1 :

 1 Figure H.1: Benchmark for the T V -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 (θ t )f (θ * ), 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.

Figure I. 1 Figure I. 1 :Figure J. 1 :

 111 Figure I.1: Benchmark for the MCP regression on variants of the Objective (columns). The curves display the violation of optimality conditions, dist(-X ⊤ (Xθ ty)/n, ∂ρ λ,γ (θ t )), as a function of time. γ is set to 3, and λ is parameterized as a fraction of the Lasso's λ max , X ⊤ y ∞ /n.

  .2).

	blitz	celer	coordinate descent	cuML[cd]
	cuML[qn]	FISTA	FISTA[adaptive-1]	FISTA[greedy]
	glmnet	ISTA	LARS	lasso jl
	noncvx-pro	scikit-learn	skglm	snapML[cpu]
	snapML[gpu]

  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 Zhang et al. (2019) for CIFAR-10, by Zheng et al. (2021) 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 × 32 and is hence not shown for MNIST.Curves are exponentially smoothed.

			Vanilla SGD			SGD + data aug. + momentum + cosine LR sched.	Best Adam
			SGD + data aug.			Best SGD						Lookahead
			SGD + data aug. + momentum	Best SGD (TF/Keras)
		20		CIFAR-10	8		SVHN		3		MNIST
	Test error (%)	10 15				4 6				1 2
		5	0	1000	2000	3000	0	2000	4000	0	0	1000	2000
				Time (s)			Time (s)				Time (s)
	Figure 5:

  ResNet18 image classification benchmark with a validation split. In dashed black is the state of the art (see caption of Figure

		Best SGD with step LR schedule (no val test)	Best SGD with cosine LR schedule (no val test)
		Best SGD with step LR schedule (val test)	Best SGD with cosine LR schedule (val test)
	Test/Val error (%)	4 6 8 10	1000	1500 CIFAR-10 2000	2500	4 6 8	1000	2000 SVHN	3000	0 1 2	1000	2000 MNIST
				Time (s)				Time (s)				Time (s)
	Figure 6:

  .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)

  Table B.2 describes the different solvers compared in this benchmark. Table B.2: List of solvers used in the Ridge benchmark in Appendix B

	Solver	References	Description		Language
	GD	Boyd and Vandenberghe (2004) Gradient Descent	Python
	Accelerated GD	Nesterov (1983)	Gradient Descent + ac-	Python
			celeration
	scikit-learn[svd]	Pedregosa et al. (2011)	SVD (Singular Value	Python (Cython)
			Decomposition)
	scikit-learn[cholesky] Pedregosa et al. (2011)	Cholesky decomposi-	Python (Cython)
			tion
	scikit-learn[lsqr]	Pedregosa et al. (2011)	regularized	least-	Python (Cython)
			squares
	scikit-learn[saga]	Pedregosa et al. (2011)	SAGA	(Varianced	Python (Cython)
			reduced	stochastic
			method)
	scikit-learn[cg]	Pedregosa et al. (2011)	Conjugate Gradient	Python (Cython)
	CD	Bertsekas (1999)	Cyclic Coordinate De-	Python (Numba)
			scent
	lightning[cd]	Blondel and Pedregosa (2016)	Cyclic Coordinate De-	Python (Cython)
			scent
	snapML[cpu]	Dünner et al. (2018)	CD		Python, C++
	snapML[gpu]

Table D .

 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 ℓ 2 -regularized logistic regression benchmark in Section 3

	Solver	References	Description	Language
	lightning[sag]	Blondel and Pedregosa (2016)	SAG	Python (Cython)
	lightning[saga]	Blondel and Pedregosa (2016)	SAGA	Python (Cython)
	lightning[cd]	Blondel and Pedregosa (2016)	Cyclic Coordi-	Python (Cython)
			nate Descent
	Tick[svrg]	Bacry et al. (2017)	Stochastic Vari-	Python, C++
			ance Reduced
			Gradient
	scikit-learn[sgd]	Pedregosa et al. (2011)	Stochastic Gra-	Python (Cython)
			dient Descent
	scikit-learn[sag]	Pedregosa et al. (2011)	SAG	Python (Cython)
	scikit-learn[saga]	Pedregosa et al. (2011)	SAGA	Python (Cython)
	scikit-learn[liblinear] Pedregosa et al. (2011), Fan et al.	Truncated New-	Python (Cython)
		(2008)	ton Conjugate-
			Gradient
	scikit-learn[lbfgs]	Pedregosa et al. (2011), Virtanen	L-BFGS	Python (Cython)
		et al. (2020)	(Quasi-Newton
			Method)
	scikit-learn[newton-cg] Pedregosa et al. (2011), Virtanen	Truncated New-	Python (Cython)
		et al. (2020)	ton Conjugate-
			Gradient
	snapml[cpu]	Dünner et al. (2018)	CD	Python, C++
	snapml[gpu]	Dünner et al. (2018)	CD + GPU	Python, C++
	cuML[gpu]	Raschka et al. (2020)	L-BFGS + GPU Python, C++

Table D .

 D 2: List of the datasets used in ℓ 2 -regularized logistic regression in Section 3

	Datasets	References	Samples (n) Features (p)	Density
	ijcnn1 madelon news20.binary Keerthi et al. (2005) Prokhorov (2001) Guyon et al. (2004) criteo Criteo-Labs (2015)	49 990 2000 19 996 45 840 617	22 500 1 355 191 1 000 000	4.5 × 10 -2 2.0 × 10 -3 3.4 × 10 -4 3.9 × 10 -5

Table E .

 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		Description	Language
	blitz		Johnson and Guestrin (2015)	CD + working set	Python, C++
	coordinate descent Friedman et al. (2010)		(Cyclic) Minimization	Python (Numba)
					along coordinates
	celer		Massias et al. (2018)		CD + working set + dual	Python (Cython)
					extrapolation
	cuML[cd]		Raschka et al. (2020)		(Cyclic) Minimization	Python, C++
					along coordinates
	cuML[qn]		Raschka et al. (2020)		Orthant-Wise Limited	Python, C++
					Memory Quasi-Newton
					(OWL-QN)
	FISTA		Beck and Teboulle (2009)		ISTA + acceleration	Python
	glmnet		Friedman et al. (2010)		CD + working set +	R, C++
					strong rules
	ISTA		Daubechies et al. (2004)		ISTA (Proximal GD)	Python
	LARS		Efron et al. (2004)		Least-Angle Regression	Python (Cython)
					algorithm (LARS)
	FISTA[adaptive-1]	Liang et al. (2022, Algo 4), Far-	FISTA + adaptive restart Python
			rens et al. (2020)
	FISTA[greedy]		Liang et al. (2022, Algo 5), Far-	FISTA + greedy restart Python
			rens et al. (2020)
	noncvx-pro		Poon and Peyré (2021)		Bilevel optim + L-	Python (Cython)
					BFGS
	skglm		Bertrand et al. (2022)		CD + working set + pri-	Python (Numba)
					mal extrapolation
	scikit-learn		Pedregosa et al. (2011)		CD	Python (Cython)
	snapML[gpu]		Dünner et al. (2018)		CD + GPU	Python, C++
	snapML[cpu]		Dünner et al. (2018)		CD	Python, C++
	lasso.jl		Kornblith (2021)		CD	Julia
	Table E.2: List of datasets used in the Lasso benchmark in Section 4
	Dataset	References	Samples (n) Features (p)	Density
	MEG	Gramfort et al. (2014)		305	7498	1.0
	news20 rcv1 MillionSong Bertin-Mahieux et al. (2011) Keerthi et al. (2005) Lewis et al. (2004)	19 996 20 242 463 715	1 355 191 47 236 90	3.4 × 10 -4 3.6 × 10 -3 1

  E.1 in addition to Figure 4.

	blitz	celer	coordinate descent	cuML[cd]
	cuML[qn]	FISTA	FISTA[adaptive-1]	FISTA[greedy]
	glmnet	ISTA	LARS	lasso.jl
	noncvx-pro	scikit-learn	skglm	snapML[cpu]
	snapML[gpu]
	E.3 Convergence in terms of iteration
	??

Table F .

 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	Krizhevsky (2009)	10	40k	10k	10k	32	✓
	SVHN	digits in natural images Netzer et al. (2011)	10	58.6k	14.6k	26k	32	✓
	MNIST	handwritten digits	LeCun et al. (2010)	10	50k	10k	10k	28	✗

Table F .

 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).

				SGD -Learning rate				SGD -Weight decay
		20	1.0e-03	1.0e-02	1.0e-01	5.0e-01	20	5.0e-05	5.0e-04	1.0e-03	5.0e-03
	Test error (%)	10 20	0 1.0e-04	1000 1.0e-03 Adam -Learning rate 2000 5.0e-03	1.0e-02	10 20	0 2.0e-03	1000 2.0e-02 Adam -Weight decay 2000 2.0e-01	5.0e-01
		10					10
			0	1000	2000	Time (s)	0	1000	2000
					Hyperparameter	SGD	Adam
					Learning Rate		0.1	0.001
					Momentum		0.9		N/A
					Weight Decay Batch Size		5 × 10 -4 128	0.02 128

Table F .

 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
	• Conv2D weights intialization: TensorFlow/Keras uses the default intialization which is a Glorot
	uniform intialization

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

		Best Adam			Vanilla SGD			SGD + data aug.
		SGD + data aug. + momentum		SGD + data aug. + momentum + step LR sched.	Best SGD
		100			CIFAR-10
	Test error (%)	25 50 75
		0	200	400	600	800	1000	1200
					Time (s)
	Figure F.2:

  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 ℓ 1 -regularized logistic regression benchmark

	Solver		References	Description		Language
	blitz		Johnson and Guestrin	CD + working set	Python, C++
			(2015)
	coordinate descent	Friedman et al. (2010)	(Cyclic) Minimization	Python (Numba)
				along coordinates
	coordinate	descent	Friedman et al. (2010)	CD + Newton	Python (Numba)
	(Newton)
	celer		Massias et al. (2018)	CD + working set + dual	Python (Cython)
				extrapolation
	copt[FISTA	line	Pedregosa et al. (2020),	FISTA (ISTA + acceler-	Python (Cython)
	search]		Beck and Teboulle (2009)	ation) + line search
	copt[PGD]		Pedregosa et al. (2020),	Proximal Gradient De-	Python (Cython)
			Combettes and Wajs	scent
			(2005)
	copt[PGD linesearch] Pedregosa et al. (2020),	Proximal Gradient De-	Python (Cython)
			Combettes and Wajs	scent + linesearch
			(2005)
	copt[saga]		Pedregosa et al. (2020)	SAGA	(Variance	Python (Cython)
				reduced	stochastic
				method)
	copt[svrg]		Pedregosa et al. (2020)	SVRG	(Variance	Python (Cython)
				reduced	stochastic
				method)
	cuML[gpu]		Raschka et al. (2020)	L-BFGS + GPU	Python, C++
	cuML[qn]		Raschka et al. (2020)	Orthant-Wise Limited	Python, C++
				Memory Quasi-Newton
				(OWL-QN)
	cyanure		Mairal (2019)	Proximal	Minimiza-	Python, C++
				tion by Incremental
				Surrogate Optimization
				(MISO)
	lightning		Blondel and Pedregosa	(Cyclic) Coordinate De-	Python (Cython)
			(2016)	scent
	scikit-learn[liblinear]Pedregosa et al. (2011),	Truncated	Newton	Python (Cython)
			Fan et al. (2008)	Conjugate-Gradient
	scikit-learn[lbfgs]	Pedregosa et al. (2011),	L-BFGS	(Quasi-	Python (Cython)
			Virtanen et al. (2020)	Newton Method)
	scikit-learn[newton-cg]Pedregosa et al. (2011),	Truncated	Newton	Python (Cython)
			Virtanen et al. (2020)	Conjugate-Gradient
	snapml[gpu=True]	Dünner et al. (2018)	CD + GPU		Python, C++
	snapml[gpu=False]	Dünner et al. (2018)	CD		Python, C++

Table G .

 G 2: List of the datasets used in the ℓ 1 -regularized logistic regression benchmark

	Datasets	References	Samples (n) Features (p)	Density
	gisette colon-cancer	Guyon et al. (2004) Guyon et al., 2004	6000 62	5000 2000	9.9 × 10 -1 1.0
	news20.binary Keerthi et al. (2005) rcv1.binary Guyon et al., 2004	19 996 20 242	1 355 191 19 959	3.4 × 10 -4 3.6 × 10 -3
	blitz		celer	coordinate descent	coordinate descent (Newton)
	copt (FISTA line search)	copt (FISTA)	copt (PGD line search)	copt (PGD)
	copt (SAGA)		copt (SVRG)	cuML[qn]	cyanure
	liblinear		lightning	snapML[cpu]	snapML[gpu]

Table I .

 I 1: List of solvers used in the MCP benchmark

	Solver	References	Short Description
	CD	Breheny and Huang (2011), Mazumder	Proximal coordinate descent
		et al. (2011)
	PGD	Bolte et al. (2014)	Proximal gradient descent
	GIST	Gong et al. (2013)	Proximal gradient + Barzilai-Borwein rule
	WorkSet CD Boisbunon et al. (2014)	Coordinate descent + working set
	skglm	Bertrand et al. (2022)	Accelerated coordinate descent + Working set

The condition number of the dataset is divided by 5.9 after scaling.

The condition number is multiplied by 407 after scaling.

The Lasso path is the set of solutions of Problem (3) as λ varies in (0, ∞).

decreasing the learning rate by a factor 10 at mid-training, and again at 3/4 of the training

The results are available online as a user-friendly interactive HTML file

We ran the same experiment on two other datasets obtaining similar figures.

It 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.

H Unidimensional total variation

The use of 1D Total Variation regularization takes its root in the taut-string algorithm [START_REF] Barlow | The Isotonic Regression Problem and Its Dual[END_REF] and can be seen as a special case of either the Rudin-Osher-Fatemi model [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] or the Generalized Lasso [START_REF] Tibshirani | The solution path of the generalized lasso[END_REF] for a quadratic data fit term. It reads

where F is a data fidelity term, X ∈ R n×p is a design matrix with n samples and p features, y ∈ R n is the target vector, λ > 0 is a regularization hyperparameter, and D ∈ R (p-1)×p is a finite difference operator defined by (Dθ) k = θ k+1θ k for all 1 ≤ k ≤ p -1 (it is also possible to use cyclic differences).

Most often, the data fidelity term is the ℓ 2 -loss F (y, z) = 1 2 ∥y -z∥ 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| µ where | • | µ is defined coordinate-wise by

Problem (7) 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 [START_REF] Bleakley | The group fused Lasso for multiple change-point detection[END_REF][START_REF] Tibshirani | Adaptive Piecewise Polynomial Estimation via Trend Filtering[END_REF], for BOLD signal deconvolution in functional MRI [START_REF] Karahanoglu | Total Activation: fMRI Deconvolution through Spatio-Temporal Regularization[END_REF][START_REF] Cherkaoui | Sparsity-Based Semi-Blind Deconvolution of Neural Activation Signal in fMRI[END_REF] or for detrending in oculomotor recordings [START_REF] Lalanne | Extraction of Nystagmus Patterns from Eye-Tracker Data with Convolutional Sparse Coding[END_REF].

The penalty θ → ∥Dθ∥ 1 is convex but non-smooth, and its proximity operator has no closed form. Yet as demonstrated by Condat (2013a), 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 Problem (7), using either primal-dual approaches [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Condat | A primal-dual splitting method for convex optimization involving Lipschitzian, proximable and linear composite terms[END_REF], or solving the dual problem [START_REF] Komodakis | Playing with Duality: An overview of recent primal-dual approaches for solving large-scale optimization problems[END_REF]. Finally, for 1-dimensional TV regularization, one can also use the synthesis formulation [START_REF] Elad | Analysis versus synthesis in signal priors[END_REF] to solve the problem. By setting z = Dθ and θ = Lz + ρ where L ∈ R p×p-1 is a lower trianglar matrix representing an integral operator (cumulative sum), the problem is equivalent to a Lasso problem, and ρ * has a closed-form expression (see e.g., Bleakley and Vert 2011 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 Problem (7) u * is retrieved as u * = Lz * + ρ * .

The code for the benchmark is available at https://github.com/benchopt/benchmark_tv_1d/ and Table H

J Zero-order optimization on standard functions

Zero-order optimization refers to scenarios where only calls to the function to minimize are possible. This is in contrast with first-order optimization where gradient information is available. Grid search, random search, evolution strategies (ES) or Bayesian optimization (BO) are popular methods to tackle such a problem and are most commonly employed for hyperparameter optimization. This setting is also known as black-box optimization.

This benchmark demonstrates the usability of Benchopt for zero-order optimization considering functions classically used in the literature [START_REF] Hansen | COCO: A Platform for Comparing Continuous Optimizers in a Black-Box Setting[END_REF]. The functions are available in the PyBenchFCN package https://github.com/Y1fanHE/PyBenchFCN/. In particular, among the 61 functions of interest we present here (see Figure J.1) a benchmark for three functions, defined for any

For each function, the domain is restricted to a box: ∥x∥ ∞ ≤ 32 for Ackley, ∥x∥ ∞ ≤ 30 for Rosenbrock and ∥x∥ ∞ ≤ 5.12 for Rastrigin. The algorithms considered in the benchmark are listed in Table J.1. As BFGS requires first-order information, gradients are approximated with finitedifferences. The code for the benchmark is available at https://github.com/benchopt/benchmark_zero_ order/.

J.1 Results

The results of the benchmark are presented in