Efficient Cross-Validation in Julia

This archive is distributed in association with the INFORMS Journal on Computing under the MIT License.

The software and data in this repository are a snapshot of the software and data that were used in the research reported on in the paper Efficient Cross-Validation for Sparse Linear Regression by R. Cory-Wright and A. Gomez.

Cite

To cite the contents of this repository, please cite both the paper and this repo, using their respective DOIs.

https://doi.org/10.1287/ijoc.2025.1336

https://doi.org/10.1287/ijoc.2025.1336.cd

Below is the BibTex for citing this snapshot of the repository.

@misc{
    EfficientCrossValidation,
    author =        {R. Cory-Wright and A. G{\'o}mez},
    publisher =     {INFORMS Journal on Computing},
    title =         {{Efficient Cross-Validation for Sparse Linear Regression}},
    year =          {2026},
    doi =           {https://doi.org/10.1287/ijoc.2025.1336.cd},
    url =           {https://github.com/INFORMSJoC/2025.1336},
    note =          {Available for download at https://github.com/INFORMSJoC/2025.1336},
}

The code is written primarily in Julia, with R called from Julia for several benchmark methods.

Introduction

The software in this package is designed to minimize the problem

$$h(\gamma,\tau) = \frac{1}{n}\sum_{j=1}^k\sum_{i\in \mathcal{N}_j} (y_i -x_i^\top \beta^{(\mathcal{N}_j)}(\gamma,\tau))^2$$

where $\beta^{(\mathcal{N}_j)}(\gamma,\tau)$ denotes an optimal solution to the following lower-level problem for any $\mathcal{N}_j$ such that

$$\beta^{(\mathcal{N}_j)}(\gamma,\tau) \in \arg\min_{\beta \in \mathbb{R}^p} \ \frac{\gamma}{2}\Vert \beta\Vert_2^2 +\Vert y^{(\mathcal{N}_j)}-X^{(\mathcal{N}_j)}\beta\Vert_2^2\quad \text{s.t.}\quad \Vert \beta\Vert_0 \leq \tau$$

using the methods described in the paper "Optimal Cross-Validation for Sparse Linear Regression" by Cory-Wright and Gomez. The rest of this readme file explains how to install/setup the code, and how to reproduce each section of the numerical results in the paper.

Installation and set up

In order to run this software, you must install a recent version of Julia from http://julialang.org/downloads/, and a recent version of the Mosek solver (academic licenses are freely available at https://www.mosek.com/products/academic-licenses/). To run parts of the code you will also need to install the Gurobi solver (academic licenses are freely available at https://www.gurobi.com/academia/academic-program-and-licenses/) Finally, you will need to have an R installation and several R packages installed, including L0Learn and GLMNet among others. We recommend following the RCall installation instructions here https://juliainterop.github.io/RCall.jl/stable/installation/

You will then need to install a number of Julia packages. Among others, you will need to install the following packages: Gurobi, StatsBase, CSV, DataFrames, JuMP, LinearAlgebra, Distributions, Random, Suppressor, Mosek, MosekTools, Polynomials, ForwardDiff, Roots, LinearAlgebra, RCall, MathOptInterface, Plots, JLD2

At this point, the "real_data_benchmark_all_IJOC.jl" files should run successfully, which reproduces Tables 3 and 4 in the paper (see also the Excel sheet "real_data_results_analysis_forIJOC_final.xlsm" which provides the raw output and the summarized output from running this script).

Reproducing the numerical experiments

Section 4.1: exact k-fold optimization and node-count reductions

The main code for this section is available in the "Crossval" folder. More information on the code can be found in the readme.md file therein.

Section 4.2: statistical results with real data

The main code for this section (run via the script file entitled "real_data_benchmark_all_IJOC") implements an optimization-based approach for selecting hyperparameters in sparse ridge regression. For a response vector Y and feature matrix X, the method searches over:

• a sparsity budget, called k in the code and corresponding to tau in the paper; and
• a ridge parameter, called gamma or γ in the code.

For each fold of a cross-validation split, the code solves or approximates a sparse ridge regression problem. The paper's main computational idea is to avoid solving all fold/sparsity mixed-integer optimization problems. Instead, it first solves perspective relaxations, derives lower and upper bounds on cross-validation loss, and solves exact mixed-integer problems only where the bounds indicate they are needed.

The real-data benchmark then compares the MIO-based approach with:

• Elastic Net via GLMNet.jl;
• SCAD and MCP via the R package ncvreg; and
• L0-regularized regression via the R package L0Learn.

To reproduce all rows in Tables 3 and 4, run the default datasets, then separately run the Housing and alcohol second-order interaction datasets as described above. Depending on solver licenses, hardware, and package versions, exact runtimes and MIP gaps may vary.

Section 4.2 Outputs and checkpoints

The benchmark appends rows to:

real_data_results_forIJOC_withgrid.csv

The columns are:

DataName, n, repNum,
kcv_MIO, k_MIO, t_MIO, mse_MIO,
kcv_glm, k_glm, t_glm, mse_glm,
kcv_mcp, k_mcp, t_mcp, mse_mcp,
kcv_scd, k_scd, t_scd, mse_scd,
kcv_l0l, k_l0l, t_l0l, mse_l0l

Here kcv_* entries are cross-validation sums of squared errors on the training/validation split, while mse_* entries are average test-set mean squared errors on the held-out 20% split. The k_* entries are estimated sparsities.

Checkpoints

Each MIO cross-validation run writes:

checkpoints/<DataName><repNum>_.jld2

for example:

checkpoints/winedata1_.jld2

The checkpoint stores at least:

k_current
gamma_current
k_best_global
kcv_best_global
iterNum
gamma_best_global
runtime_sofar

These are used to resume interrupted runs.

Thank you

Thank you for your interest in the paper. Please let us know if you encounter any issues using this code, or have comments or questions. Feel free to email us anytime.

Ryan Cory-Wright r.cory-wright@imperial.ac.uk

Andres Gomez gomezand@usc.edu