HPR-QP: A dual Halpern Peaceman–Rachford (HPR) method for solving large-scale convex composite quadratic programming (CCQP).
This version represents a major architectural overhaul with significant improvements:
- Single codebase for all problem types (standard QP, QAP, LASSO)
- Modular Q operator system for extensibility: easily add custom problem types
- Full CPU implementation in addition to GPU acceleration
- Automatic device selection via
use_gpuparameter
- Native JuMP/MathOptInterface (MOI) support for easy modeling
- Use HPR-QP directly as a JuMP optimizer
- Initialize via
initial_xandinitial_yparameters - Resume optimization from previous (auto-saved) solutions
- Automatically save best solution during optimization (
auto_save=true) - Resume from saved states for long-running problems
-
$Q$ is a positive semidefinite self-adjoint linear operator; -
$Q$ 's matrix representation may not be computable in large-scale instances, such as QAP relaxation and LASSO problems; -
$\phi$ is a proper closed convex function.
Before using HPR-QP, make sure the following dependencies are installed:
- Julia (Recommended version:
1.10.4) - CUDA (Required for GPU acceleration; install the appropriate version for your GPU and Julia)
- Required Julia packages
To install the required Julia packages and build the HPR-QP environment, run:
julia --project -e 'using Pkg; Pkg.instantiate()'To verify that CUDA is properly installed and working with Julia, run:
using CUDA
CUDA.versioninfo()Before running the scripts, please modify
run_single_file.jlorrun_dataset.jlin the demo directory to specify the data path and result path according to your setup.
To test the script on a single instance:
julia --project demo/run_single_file.jlTo process all files in a directory:
julia --project demo/run_dataset.jlQAP Instances:
The.matfile for QAP should include the matrices$A$ ,$B$ ,$S$ , and$T$ .
For details, refer to Section 4.5 of the paper.
SeeHPR-QP_QAP_LASSO/demo/demo_QAP.jlfor an example of generating such files.LASSO Instances:
The.matfile for LASSO should contain the matrix$A$ , vector$b$ .
This example demonstrates how to construct a CQP model using the JuMP modeling language in Julia and export it to MPS format for use with the HPR-QP solver.
julia --project demo/demo_JuMP.jlThe script:
- Builds a CQP model.
- Uses HPR-QP to solve the CQP instance.
Remark: If the model may be infeasible or unbounded, you can use HiGHS to check it.
using JuMP, HiGHS
## read a model from file (or create in other ways)
mps_file_path = "xxx" # your file path
model = read_from_file(mps_file_path)
## set HiGHS as the optimizer
set_optimizer(model, HiGHS.Optimizer)
## solve it
optimize!(model)This example demonstrates how to construct and solve a CQP problem directly in Julia without relying on JuMP.
julia --project demo/demo_QAbc.jlThis example demonstrates how to construct and solve a random LASSO instance.
julia --project demo/demo_LASSO.jlYou may notice that solving a single instance — or the first instance in a dataset — appears slow. This is due to Julia’s Just-In-Time (JIT) compilation, which compiles code on first execution.
💡 Tip for Better Performance:
To reduce repeated compilation overhead, it’s recommended to run scripts from an IDE like VS Code or the Julia REPL in the terminal.
julia --projectThen, at the Julia REPL, run demo/demo_QAbc.jl (or other scripts):
include("demo/demo_QAbc.jl")CAUTION:
If you encounter the error message:
Error: Error during loading of extension AtomixCUDAExt of Atomix, use Base.retry_load_extensions() to retry.Don’t panic — this is usually a transient issue. Simply wait a few moments; the extension typically loads successfully on its own.
Below is a list of the parameters in HPR-QP along with their default values and usage:
| Parameter | Default Value | Description |
|---|---|---|
stoptol | 1e-6 | Stopping tolerance for convergence checks. |
sigma | -1 (auto) | Initial value of the σ parameter used in the algorithm. |
max_iter | typemax(Int32) | Maximum number of iterations allowed. |
sigma_fixed | false | Whether σ is fixed throughout the optimization process. |
time_limit | 3600.0 | Maximum allowed runtime (seconds) for the algorithm. |
eig_factor | 1.05 | Factor used to scale the maximum eigenvalue estimation. |
check_iter | 100 | Frequency (in iterations) to check for convergence or perform other checks. |
warm_up | false | Determines if a warm-up phase is performed before main execution. |
spmv_mode_Q | "auto" | Mode for Q matrix-vector multiplication (e.g., "auto", "CUSPARSE", "customized", "operator"). |
spmv_mode_A | "auto" | Mode for A matrix-vector multiplication (e.g., "auto", "CUSPARSE", "customized"). |
print_frequency | -1 (auto) | Frequency (in iterations) for printing progress or logging information. |
device_number | 0 | GPU device number (e.g., 0, 1, 2, 3). |
use_Ruiz_scaling | true | Whether to apply Ruiz scaling to the problem data. |
use_bc_scaling | false | Whether to apply bc scaling. (For QAP and LASSO, only this scaling is applicable) |
use_l2_scaling | false | Whether to apply L2-norm based scaling. |
use_Pock_Chambolle_scaling | true | Whether to apply Pock-Chambolle scaling to the problem data. |
problem_type | "QP" | Type of problem being solved (e.g., "QP", "LASSO", "QAP"). |
lambda | 0.0 | Regularization parameter for LASSO problems. |
initial_x | nothing | Initial primal solution for warm-start. |
initial_y | nothing | Initial dual solution for warm-start. |
auto_save | false | Automatically save best x, y, z, w, and sigma during optimization. |
save_filename | "hprqp_autosave.h5" | Filename for auto-save HDF5 file. |
verbose | true | Enable verbose output during optimization. |
use_gpu | true | Whether to use GPU acceleration (requires CUDA). |
After solving an instance, you can access the result variables as shown below:
# Example from /demo/demo_QAbc.jl
println("Objective value: ", result.primal_obj)
println("x1 = ", result.x[1])
println("x2 = ", result.x[2])| Category | Variable | Description |
|---|---|---|
| Iteration Counts | iter | Total number of iterations performed by the algorithm. |
iter_4 | Number of iterations required to achieve an accuracy of 1e-4. | |
iter_6 | Number of iterations required to achieve an accuracy of 1e-6. | |
| Time Metrics | time | Total time in seconds taken by the algorithm. |
time_4 | Time in seconds taken to achieve an accuracy of 1e-4. | |
time_6 | Time in seconds taken to achieve an accuracy of 1e-6. | |
power_time | Time in seconds used by the power method. | |
| Objective Values | primal_obj | The primal objective value obtained. |
gap | The gap between the primal and dual objective values. | |
| Residuals | residuals | Relative residuals of the primal feasibility, dual feasibility, and duality gap. |
| Algorithm Status | status | The final status of the algorithm: - OPTIMAL: Found optimal solution- MAX_ITER: Max iterations reached- TIME_LIMIT: Time limit reached |
| Solution Vectors | x | The final solution vector x. |
y | The final solution vector y. | |
z | The final solution vector z. | |
w | The final solution vector w. |
@article{chen2025hpr,
title={HPR-QP: A dual Halpern Peaceman-Rachford method for solving large-scale convex composite quadratic programming},
author={Chen, Kaihuang and Sun, Defeng and Yuan, Yancheng and Zhang, Guojun and Zhao, Xinyuan},
journal={arXiv preprint arXiv:2507.02470},
year={2025}
}