diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml new file mode 100644 index 0000000..da49e42 --- /dev/null +++ b/.github/workflows/docs.yml @@ -0,0 +1,70 @@ +name: Documentation + +on: + push: + branches: [main] + paths: + - "docs/**" + - "zensical.toml" + - ".github/workflows/docs.yml" + pull_request: + branches: [main] + paths: + - "docs/**" + - "zensical.toml" + - ".github/workflows/docs.yml" + workflow_dispatch: + +permissions: + contents: read + +concurrency: + group: pages-${{ github.ref }} + cancel-in-progress: true + +jobs: + build: + name: Build site + runs-on: ubuntu-latest + steps: + - uses: actions/configure-pages@v6 + if: github.event_name != 'pull_request' + - uses: actions/checkout@v7 + - uses: actions/setup-python@v7 + with: + python-version: "3.14" + - name: Install Zensical + run: python -m pip install zensical==0.0.56 + - name: Build documentation + run: | + zensical build --clean --strict + find site -type f -name '.DS_Store' -delete + - name: Upload PR preview artifact + if: github.event_name == 'pull_request' + uses: actions/upload-artifact@v7 + with: + name: docs-preview + path: site + if-no-files-found: error + retention-days: 7 + - name: Upload Pages artifact + if: github.event_name != 'pull_request' + uses: actions/upload-pages-artifact@v5 + with: + path: site + + deploy: + name: Deploy to GitHub Pages + if: github.event_name != 'pull_request' + needs: build + permissions: + pages: write + id-token: write + environment: + name: github-pages + url: ${{ steps.deployment.outputs.page_url }} + runs-on: ubuntu-latest + steps: + - name: Deploy + id: deployment + uses: actions/deploy-pages@v5 diff --git a/README.md b/README.md index b14f5dc..edb7f45 100644 --- a/README.md +++ b/README.md @@ -3,10 +3,11 @@ [![License](https://img.shields.io/badge/License-Apache%202.0-blue.svg)](LICENSE) [![GitHub release](https://img.shields.io/github/release/MIT-Lu-Lab/cuPDLPx.svg)](https://github.com/MIT-Lu-Lab/cuPDLPx/releases) [![PyPI version](https://badge.fury.io/py/cupdlpx.svg)](https://pypi.org/project/cupdlpx/) +[![Documentation](https://img.shields.io/badge/docs-latest-blue.svg)](https://mit-lu-lab.github.io/cuPDLPx/) [![arXiv](https://img.shields.io/badge/arXiv-2407.16144-B31B1B.svg)](https://arxiv.org/abs/2407.16144) [![arXiv](https://img.shields.io/badge/arXiv-2507.14051-B31B1B.svg)](https://arxiv.org/abs/2507.14051) -**cuPDLPx** is a GPU-accelerated linear programming solver based on a restarted Halpern PDHG method specifically tailored for GPU architectures. It incorporates a Halpern update scheme, an adaptive restart scheme, and a PID-controlled primal weight, resulting in substantial empirical improvements over its predecessor, **[cuPDLP](https://github.com/jinwen-yang/cuPDLP.jl)**, on standard LP benchmark suites. +**cuPDLPx** is a GPU-accelerated linear programming solver based on a restarted Halpern PDHG method specifically tailored for GPU architectures. It incorporates a Halpern update scheme, an adaptive restart scheme, and a PID-controlled primal weight, resulting in substantial empirical improvements over its predecessor, **[cuPDLP](https://github.com/jinwen-yang/cuPDLP.jl)**, on standard LP benchmark suites. cuPDLPx solves linear programs of the form ```math @@ -23,137 +24,18 @@ Our work is presented in two papers: * **Theoretical Paper:** [Restarted Halpern PDHG for Linear Programming](https://arxiv.org/pdf/2407.16144) provides the mathematical foundation for our method. -## Installation +For installation instructions, examples, solver parameters, and algorithm details, see the [cuPDLPx documentation](https://mit-lu-lab.github.io/cuPDLPx/). -### Requirements -* **GPU:** NVIDIA GPU with CUDA 12.4+, or AMD GPU with ROCm 7.2+. -* **Build Tools:** CMake (≥ 3.20), GCC, and NVCC (CUDA) or hipcc (ROCm). +## Interfaces -> **SpMV backend** is selected automatically at compile time based on cuSPARSE version: -> - `cusparseSpMV` — CUDA 12.4 – 13.2 (cuSPARSE < 12.8.1) -> - `cusparseSpMVOp` — CUDA 13.3+ (cuSPARSE ≥ 12.8.1) -> -> On AMD GPUs the solver uses the `hipsparseSpMV` backend via hipSPARSE. +| Interface | Description | +| --- | --- | +| [Command line](https://mit-lu-lab.github.io/cuPDLPx/guides/command-line/) | Solve MPS files from a shell. | +| [Python](https://mit-lu-lab.github.io/cuPDLPx/guides/python/) | Build and solve LPs with [NumPy](https://numpy.org/doc/stable/) and [SciPy](https://docs.scipy.org/doc/scipy/). | +| [Julia](https://mit-lu-lab.github.io/cuPDLPx/guides/julia/) | Use cuPDLPx through [JuMP](https://jump.dev/JuMP.jl/stable/) and [MathOptInterface](https://jump.dev/MathOptInterface.jl/stable/). | +| [C](https://mit-lu-lab.github.io/cuPDLPx/guides/c-api/) | Embed cuPDLPx in native applications. | -### Build from Source -Clone the repository and compile the project using CMake. -```bash -git clone git@github.com:MIT-Lu-Lab/cuPDLPx.git -cd cuPDLPx -cmake -B build -cmake --build build --clean-first -``` -This will create the solver binary at `./build/cupdlpx`. - -#### Building for AMD GPUs (ROCm/HIP) -To target AMD GPUs, configure with `-DUSE_HIP=ON` and select the GPU -architecture with `-DCMAKE_HIP_ARCHITECTURES`. The CUDA sources are compiled -as HIP and the cuBLAS/cuSPARSE/CUB calls are mapped to hipBLAS/hipSPARSE/hipCUB. -```bash -cmake -B build -DUSE_HIP=ON -DCMAKE_HIP_ARCHITECTURES=gfx90a -DCMAKE_PREFIX_PATH=/opt/rocm -cmake --build build --clean-first -``` -Set `CMAKE_HIP_ARCHITECTURES` to match your GPU (for example `gfx90a` for -MI200, `gfx1100` for RDNA3 desktop, or `gfx1201` for RDNA4). If the ROCm -install is not on CMake's default search path, point `-DCMAKE_PREFIX_PATH` at -it (e.g. `/opt/rocm`) so `find_package` can locate hip, hipBLAS, hipSPARSE, -and hipCUB. The resulting `./build/cupdlpx` binary is used exactly as in the -CUDA build. - -#### Verifying the Installation -Run a small test problem to confirm that the solver was built correctly. -```bash -# 1. Download a test instance from the MIPLIB library -wget -P test/ https://miplib.zib.de/WebData/instances/2club200v15p5scn.mps.gz - -# 2. Solve the problem and write output to the current directory (.) -./build/cupdlpx test/2club200v15p5scn.mps.gz test/ -``` -If the solver runs and creates output files, your installation is successful. - -### Python Package Installation -To use cuPDLPx in Python, you can install the pre-built package `cupdlpx` directly from PyPI: -```bash -pip install cupdlpx -``` -Or build from source: -``` -git clone https://github.com/MIT-Lu-Lab/cuPDLPx.git -cd cuPDLPx -pip install . -``` - -## Usage & Interfaces -### Command-line Interface - -After building the project, the `./build/cupdlpx` binary can be invoked from the command line as follows: - -```bash -./build/cupdlpx [OPTIONS] -``` - -#### Arguments -- ``: The path to the input linear programming problem. Both plain (`.mps`) and gzipped (`.mps.gz`) files are supported. -- ``: The directory where the output files will be saved. - -#### Solver Options - -| Option | Type | Description | Default | -| :--- | :--- | :--- | :--- | -| `-h`, `--help` | `flag` | Display the help message. | N/A | -| `-v`, `--verbose` | `flag` | Verbose logging (enabled by default). | `true` | -| `-q`, `--quiet` | `flag` | Disable verbose logging. | `false` | -| `--debug` | `flag` | Developer diagnostics (implies verbose). | `false` | -| `--time_limit` | `double` | Time limit in seconds. | `3600.0` | -| `--iter_limit` | `int` | Iteration limit. | `2147483647` | -| `--opt_norm` | `string` | Norm for optimality criteria: `l2` or `linf` | `l2` | -| `--eps_opt` | `double` | Relative optimality tolerance. | `1e-4` | -| `--eps_feas` | `double` | Relative feasibility tolerance. | `1e-4` | -| `--geo_mean_iter` | `int` | Iterations of geometric-mean scaling | `12` | -| `--l_inf_ruiz_iter` | `int` | Iterations of L-inf Ruiz rescaling| `10` | -| `--no_pock_chambolle` | `flag` | Disable Pock-Chambolle rescaling | `enabled` | -| `--pock_chambolle_alpha` | `float` | Value for Pock-Chambolle alpha | `1.0` | -| `--no_bound_obj_rescaling` | `flag` | Disable bound objective rescaling | `enabled` | -| `--eval_freq` | `int` | Termination evaluation frequency | `200` | -| `--sv_max_iter` | `int` | Max iterations for singular value estimation | `5000` | -| `--sv_tol` | `float` | Tolerance for singular value estimation | `1e-4` | -| `--no_presolve` | `flag` | Disable presolve | `enabled` | -| `-f`,`--feasibility_polishing` |`flag` | Run the polishing loop | `false` | -| `--eps_feas_polish` | `double` | Relative tolerance for polishing | `1e-6` | -| `--no_active_set_boost` | `flag` | Disable the active-set stepsize boost | `enabled` | -| `--asb_activation_tol` | `double` | Residual threshold at which the boost activates | `1e-4` | -| `--asb_window_iter` | `int` | Number of recent iterations used to identify the active set | `10000` | -| `--asb_safety_factor` | `double` | Boosted step = factor / estimated singular value | `0.9` | -| `--asb_max_reverts` | `int` | Divergences tolerated before the boost turns off; a diverged step is always reverted | `2` | -| `--asb_min_raise_ratio` | `double` | Minimum ratio for a step increase | `1.1` | -| `--asb_reestimate_change_ratio` | `double` | Fraction of the active set that must change before re-estimating | `0.01` | -| `--asb_constraint_tol` | `double` | Tolerance for treating a constraint as binding | `1e-8` | -| `--asb_variable_tol` | `double` | Tolerance for treating a variable as at its bound | `1e-8` | -| `--asb_divergence_ceiling_ratio` | `double` | Step ceiling after a revert, relative to the diverged step | `0.7` | -| `--asb_divergence_margin` | `double` | Allowed fixed-point error increase before a revert | `0.05` | - -The active-set boost enlarges the stepsize late in the solve, based on the constraints and bounds identified as active. - -#### Output Files -The solver generates three text files in the specified . The filenames are derived from the input file's basename. For an input `INSTANCE.mps.gz`, the output will be: -``` -/ -├── INSTANCE_summary.txt # Statistics, timings, and termination status -├── INSTANCE_primal_solution.txt # Primal solution vector -└── INSTANCE_dual_solution.txt # Dual solution vector -``` - -### Python Interface -The `cupdlpx` Python package supports building and solving LPs directly with `NumPy` and `SciPy`. -Documentation and examples are available in the [Python API Guide](python/README.md). - -### Julia Interface -`CuPDLPx.jl` provides a `JuMP`/`MathOptInterface` wrapper for cuPDLPx, enabling cuPDLPx to be used as a solver backend for `JuMP` models. Documentation and examples are available at: [CuPDLPx.jl](https://github.com/MIT-Lu-Lab/CuPDLPx.jl). - -### C Interface -The public C API is defined in header file [`include/cupdlpx.h`](include/cupdlpx.h). A detailed description with usage examples can be found in the [C API Guide](docs/C_API.md). - -## Reference +## References If you use cuPDLPx or the ideas in your work, please cite the source below. ```bibtex diff --git a/docs/C_API.md b/docs/C_API.md deleted file mode 100644 index 4ae6950..0000000 --- a/docs/C_API.md +++ /dev/null @@ -1,92 +0,0 @@ - -### C Interface - -cuPDLPx provides a C API for directly solving LPs in memory, defined in header file [`include/cupdlpx.h`](../include/cupdlpx.h). This is useful when integrating cuPDLPx into other C/C++ projects or when generating problems programmatically. - -#### Functions and Parameters - -The C API involves two main functions: - -```c -lp_problem_t *create_lp_problem( - const double *objective_c, // objective vector (length n) - const matrix_desc_t *A_desc, // constraint matrix (m×n) - const double *con_lb, // constraint lower bounds (length m) - const double *con_ub, // constraint upper bounds (length m) - const double *var_lb, // variable lower bounds (length n) - const double *var_ub, // variable upper bounds (length n) - const double *objective_constant, // scalar objective offset - const objective_sense_t *objective_sense // objective sense (NULL → minimize) -); - -cupdlpx_result_t* solve_lp_problem( - lp_problem_t* prob, - const pdhg_parameters_t* params // NULL → use default parameters -); -``` - -`create_lp_problem` parameters: -- `objective_c`: Objective vector. If `NULL`, defaults to all zeros. -- `A_desc`: Matrix descriptor. Supports `matrix_dense`, `matrix_csr`, `matrix_csc`, `matrix_coo`. -- `con_lb`: Constraint lower bounds. If `NULL`, defaults to all `-INFINITY`. -- `con_ub`: Constraint upper bounds. If `NULL`, defaults to all `+INFINITY`. -- `var_lb`: Variable lower bounds. If `NULL`, defaults to all `-INFINITY`. -- `var_ub`: Variable upper bounds. If `NULL`, defaults to all `+INFINITY`. -- `objective_constant`: Scalar constant term added to the objective value. If `NULL`, defaults to `0.0`. -- `objective_sense`: Objective sense, `OBJECTIVE_SENSE_MINIMIZE` or `OBJECTIVE_SENSE_MAXIMIZE`. If `NULL`, defaults to minimize. - - -`solve_lp_problem` parameters: -- `prob`: An LP problem built with `create_LP_problem`. The solver may clean up the matrix (e.g., drop near-zero entries), so the struct must be mutable. -- `params`: Solver parameters. If `NULL`, the solver will use default parameters. - -#### Example: Solving a Small LP -```c -#include "cupdlpx.h" -#include -#include - -int main() { - int m = 3; // number of constraints - int n = 2; // number of variables - - // Dense matrix A - double A[3][2] = { - {1.0, 2.0}, - {0.0, 1.0}, - {3.0, 2.0} - }; - - // Describe A - matrix_desc_t A_desc; - A_desc.m = m; A_desc.n = n; - A_desc.fmt = matrix_dense; - A_desc.data.dense.A = &A[0][0]; - - // Objective coefficients - double c[2] = {1.0, 1.0}; - - // Constraint bounds: l <= A x <= u - double l[3] = {5.0, -INFINITY, -INFINITY}; - double u[3] = {5.0, 2.0, 8.0}; - - // Build the problem - lp_problem_t* prob = create_lp_problem( - c, &A_desc, l, u, NULL, NULL, NULL, NULL); - - // Solve (NULL → use default parameters) - cupdlpx_result_t* res = solve_lp_problem(prob, NULL); - - printf("Termination reason: %d\n", res->termination_reason); - printf("Primal objective: %.6f\n", res->primal_objective_value); - printf("Dual objective: %.6f\n", res->dual_objective_value); - for (int j = 0; j < res->num_variables; ++j) { - printf("x[%d] = %.6f\n", j, res->primal_solution[j]); - } - - lp_problem_free(prob); - cupdlpx_result_free(res); - - return 0; -} -``` \ No newline at end of file diff --git a/docs/algorithm/base-algorithm.md b/docs/algorithm/base-algorithm.md new file mode 100644 index 0000000..088cd04 --- /dev/null +++ b/docs/algorithm/base-algorithm.md @@ -0,0 +1,165 @@ +--- +description: PDHG, Halpern anchoring, and reflection in the base iteration of cuPDLPx. +--- + +# Base algorithm + +cuPDLPx combines PDHG, Halpern anchoring, and reflection. Within each restart +epoch, let $(x^k,y^k)$ be the current primal–dual iterate and +$(x^0,y^0)$ the anchor. The index $k$ counts iterations within the epoch; +we omit the epoch index used in +[Algorithm 1](index.md#restarted-reflected-halpern-pdhg). + +## PDHG + +For the [bound-form LP](index.md#lp-formulation), let +$\mathcal X=[\ell_v,u_v]$ and $\mathcal S=[\ell_c,u_c]$. +With primal and dual [step sizes](step-size.md) $\tau$ and $\sigma$, a +primal–dual hybrid gradient (PDHG) step[\[1\]](#ref-chambolle-pock) is + +$$ +\begin{aligned} +\widehat x^{k+1} +&=\operatorname{proj}_{\mathcal X} +\left(x^k-\tau(c-A^\top y^k)\right),\\[0.4em] +\widehat y^{k+1} +&=y^k-\sigma A(2\widehat x^{k+1}-x^k) +-\sigma\operatorname{proj}_{-\mathcal S} +\left(\sigma^{-1}y^k-A(2\widehat x^{k+1}-x^k)\right). +\end{aligned} +$$ + +The primal step projects a gradient update onto the variable bounds. The +dual step uses the extrapolated primal iterate $2\widehat x^{k+1}-x^k$ and +a projection onto $-\mathcal S=[-u_c,-\ell_c]$. Both projections clip each +coordinate to its allowed interval. + +Standard PDHG sets +$(x^{k+1},y^{k+1})=(\widehat x^{k+1},\widehat y^{k+1})$. +The fixed points of this update are primal–dual solutions of the LP. + +## Halpern iteration + +The Halpern iteration[\[2\]](#ref-halpern) takes a convex combination +of the PDHG update and a fixed anchor: + +$$ +\begin{aligned} +x^{k+1} +&=\frac{k+1}{k+2}\widehat x^{k+1} ++\frac{1}{k+2}x^0,\\[0.4em] +y^{k+1} +&=\frac{k+1}{k+2}\widehat y^{k+1} ++\frac{1}{k+2}y^0. +\end{aligned} +$$ + +The anchor remains fixed within each epoch, while its weight $1/(k+2)$ +decreases with each iteration. + +## Restart + +Restart is an important technique for accelerating the convergence of +Halpern PDHG. The anchor draws the iterates toward the epoch's starting +point. When this point is far from a solution, anchoring can slow further +progress. + +cuPDLPx therefore uses an [adaptive restart scheme](restart.md) based on the +fixed-point error and the number of iterations since the last restart. +At each restart, the anchor is replaced with the latest PDHG iterate and +$k$ is reset to zero. + +## Reflection + +Reflection (a special case of +**over-relaxation**[\[4\]](#ref-eckstein-bertsekas)) doubles the displacement +from the current iterate to its PDHG update: + +$$ +\begin{aligned} +\bar x^{k+1}&=2\widehat x^{k+1}-x^k,\\[0.4em] +\bar y^{k+1}&=2\widehat y^{k+1}-y^k. +\end{aligned} +$$ + +The reflected operator has the same fixed points as PDHG and, under the +step-size condition discussed in [Theoretical insights](#theoretical-insights), +is nonexpansive in the +[PDHG metric](restart.md#fixed-point-error): applying the operator to two +points does not increase the distance between them in that metric. +The Halpern iteration can therefore be applied to this +operator[\[3\]](#ref-lu-yang),[\[5\]](#ref-hpr-lp). + +With reflection coefficient $\gamma$, the combined update is + +$$ +\begin{aligned} +x^{k+1} +&=\frac{k+1}{k+2} +\left[\gamma\bar x^{k+1}+(1-\gamma)x^k\right] ++\frac{1}{k+2}x^0,\\[0.4em] +y^{k+1} +&=\frac{k+1}{k+2} +\left[\gamma\bar y^{k+1}+(1-\gamma)y^k\right] ++\frac{1}{k+2}y^0. +\end{aligned} +$$ + +Here $\gamma$ is +[`reflection_coefficient`](../reference/parameters.md#step-size-and-reflection). +The default $\gamma=1$ uses full reflection; $\gamma=\tfrac12$ recovers +the unreflected Halpern PDHG update. The primal reflection +$2\widehat x^{k+1}-x^k$ is also the extrapolated point used in the PDHG dual +step. + +## Theoretical insights + +For fixed step sizes satisfying $\tau\sigma\lVert A\rVert_2^2<1$, the PDHG +operator $\mathcal T$ is **firmly nonexpansive** in the +[PDHG metric](restart.md#fixed-point-error). Its reflection $2\mathcal T-I$ +is therefore **nonexpansive**, which is sufficient for applying Halpern +iteration. + +**Halpern iteration** uses anchoring to accelerate fixed-point iterations. +When the LP has a primal–dual solution, both Halpern PDHG and its reflected +variant have an $O(1/k)$ bound on the +[fixed-point error](restart.md#fixed-point-error) +(Lemmas 3 and 11)[\[3\]](#ref-lu-yang). +Reflection thus allows a larger step while retaining this convergence +guarantee. + +**Restart** improves this sublinear rate by exploiting the **sharpness of +LP**. For feasible and bounded LPs, sharpness bounds the distance to the +solution set by a constant multiple of the fixed-point error on bounded +regions. Combining this relationship with Halpern's $O(1/k)$ bound, suitable +restart rules yield accelerated linear convergence by refreshing the +anchor after sufficient progress. + +For the precise assumptions, convergence rates, and proofs, see Sections 3 +and 6 of the theoretical paper, +[*Restarted Halpern PDHG for Linear Programming*](https://arxiv.org/abs/2407.16144)[\[3\]](#ref-lu-yang). + +## References + +\[1\] Antonin Chambolle and Thomas Pock. +[*A First-Order Primal-Dual Algorithm for Convex Problems with Applications +to Imaging*](https://doi.org/10.1007/s10851-010-0251-1). +*Journal of Mathematical Imaging and Vision*, 40:120–145, 2011. + +\[2\] Benjamin Halpern. +[*Fixed Points of Nonexpanding Maps*](https://doi.org/10.1090/S0002-9904-1967-11864-0). +*Bulletin of the American Mathematical Society*, 73:957–961, 1967. + +\[3\] Haihao Lu and Jinwen Yang. +[*Restarted Halpern PDHG for Linear Programming*](https://arxiv.org/abs/2407.16144), +2024. + +\[4\] Jonathan Eckstein and Dimitri P. Bertsekas. +[*On the Douglas–Rachford Splitting Method and the Proximal Point Algorithm +for Maximal Monotone Operators*](https://doi.org/10.1007/BF01581204). +*Mathematical Programming*, 55:293–318, 1992. + +\[5\] Kaihuang Chen, Defeng Sun, Yancheng Yuan, Guojun Zhang, +and Xinyuan Zhao. +[*HPR-LP: An Implementation of an HPR Method for Solving Linear +Programming*](https://arxiv.org/abs/2408.12179), 2024. diff --git a/docs/algorithm/feasibility-polishing.md b/docs/algorithm/feasibility-polishing.md new file mode 100644 index 0000000..0c985ce --- /dev/null +++ b/docs/algorithm/feasibility-polishing.md @@ -0,0 +1,80 @@ +--- +description: Primal and dual feasibility polishing after the main solve in cuPDLPx. +--- + +# Feasibility polishing + +Feasibility polishing solves separate primal and dual feasibility problems +after the main solve. It can reduce feasibility residuals, but may increase +the primal–dual gap. Polishing is disabled by default. + +!!! note "Scheduling" + + cuPDLPx uses the feasibility problems from PDLP + (Algorithm 4)[\[1\]](#ref-pdlp), but runs polishing only after the + main solve. PDLP interleaves polishing with the main iterations. + +## Activation + +When enabled, polishing runs unless + +- the main solve ended with primal or dual infeasibility, or +- both relative feasibility residuals already meet the polishing tolerance, + `1e-6` by default. + +The primal phase runs first, followed by the dual phase. + +## Primal polishing + +$$ +\begin{aligned} +\operatorname*{minimize}_{x\in\mathbb R^n}\quad + & 0 \\ +\text{subject to}\quad + & \ell_c \le Ax \le u_c, \\ + & \ell_v \le x \le u_v. +\end{aligned} +$$ + +The primal phase initializes $x$ from the main solve and sets $y=0$. +It applies restarted reflected Halpern PDHG, using the relative primal +residual as the convergence criterion. + +## Dual polishing + +$$ +\begin{aligned} +\operatorname*{maximize}_{y\in\mathbb R^m,\,r\in\mathbb R^n}\quad + & 0 \\ +\text{subject to}\quad + & c-A^\top y=r, \\ + & y\in\mathcal Y, \\ + & r\in\mathcal R. +\end{aligned} +$$ + +Here $r$ is the dual-slack vector, and $\mathcal Y$ and $\mathcal R$ impose +the [dual sign constraints](index.md#dual-problem). The dual phase +initializes $y$ from the main solve and sets $x=0$. It uses the relative +dual residual as the convergence criterion. + +## Result + +Each phase replaces its corresponding primal or dual values only if it +reaches the polishing tolerance. Otherwise, the values from the main solve +are retained. The objective values and primal–dual gap reflect the accepted +updates, but the main solve's termination status is unchanged. The final +gap can therefore exceed the optimality tolerance even when the status is +`OPTIMAL`. Check the final residuals and gap before using the result; see +[Results and status](../getting-started/results.md). + +## Parameters + +See [Parameters](../reference/parameters.md#feasibility-polishing) to enable +polishing and set its tolerance. + +## References + +\[1\] David Applegate et al. +[*PDLP: A Practical First-Order Method for Large-Scale Linear +Programming*](https://arxiv.org/abs/2501.07018), 2025. diff --git a/docs/algorithm/index.md b/docs/algorithm/index.md new file mode 100644 index 0000000..752b893 --- /dev/null +++ b/docs/algorithm/index.md @@ -0,0 +1,254 @@ +--- +description: LP formulation and the restarted reflected Halpern PDHG algorithm implemented by cuPDLPx. +--- + +# Algorithm + +cuPDLPx solves linear programs using **restarted reflected Halpern +primal–dual hybrid gradient (r2HPDHG)**[\[1\]](#ref-cupdlpx). +Each iteration updates the primal and dual iterates using sparse +matrix–vector products and projections onto the bound constraints. + +## LP formulation + +### Primal problem + +cuPDLPx solves linear programs in bound form: + +$$ +\begin{aligned} +\operatorname*{minimize}_{x\in\mathbb R^n}\quad + & c^\top x + c_0 \\ +\text{subject to}\quad + & \ell_c \le Ax \le u_c, \\ + & \ell_v \le x \le u_v. +\end{aligned} +$$ + +Here $A\in\mathbb R^{m\times n}$ is the constraint matrix, +$c\in\mathbb R^n$ is the objective vector, and $c_0$ is the objective +constant. The vectors $\ell_c,u_c$ and $\ell_v,u_v$ are the constraint +and variable bounds, respectively. Bounds may be infinite. + +Define the variable and constraint sets + +$$ +\mathcal X=\{x:\ell_v\le x\le u_v\},\qquad +\mathcal S=\{s:\ell_c\le s\le u_c\}. +$$ + +### Dual problem + +$$ +\begin{aligned} +\operatorname*{maximize}_{y\in\mathcal Y,\,r\in\mathcal R}\quad + & -p(-y;\ell_c,u_c)-p(-r;\ell_v,u_v)+c_0 \\ +\text{subject to}\quad + & c-A^\top y=r. +\end{aligned} +$$ + +The support function of the interval $[\ell,u]$ is + +$$ +p(z;\ell,u)=u^\top z^+-\ell^\top z^-, +\qquad +z^+=\max(z,0),\quad z^-=\max(-z,0). +$$ + +Let $\mathbb R_+=[0,+\infty)$ and $\mathbb R_-=(-\infty,0]$. The dual and +dual-slack domains are the Cartesian products +$\mathcal Y=\prod_{i=1}^m\mathcal Y_i$ and +$\mathcal R=\prod_{j=1}^n\mathcal R_j$, where + +$$ +\mathcal Y_i= +\begin{cases} +\{0\}, & \ell_{c,i}=-\infty,\ u_{c,i}=+\infty,\\ +\mathbb R_-, & \ell_{c,i}=-\infty,\ u_{c,i}\in\mathbb R,\\ +\mathbb R_+, & \ell_{c,i}\in\mathbb R,\ u_{c,i}=+\infty,\\ +\mathbb R, & \text{otherwise}, +\end{cases} +\qquad +\mathcal R_j= +\begin{cases} +\{0\}, & \ell_{v,j}=-\infty,\ u_{v,j}=+\infty,\\ +\mathbb R_-, & \ell_{v,j}=-\infty,\ u_{v,j}\in\mathbb R,\\ +\mathbb R_+, & \ell_{v,j}\in\mathbb R,\ u_{v,j}=+\infty,\\ +\mathbb R, & \text{otherwise}. +\end{cases} +$$ + +Here $r$ is the dual-slack variable. + +### Saddle-point formulation + +The LP is equivalent to + +$$ +\min_{x\in\mathcal X}\max_{y\in\mathcal Y} +L(x,y) +=c^\top x+c_0-y^\top Ax-p(y;-u_c,-\ell_c). +$$ + +The objective constant $c_0$ shifts both objective values without changing +the feasible set or the PDHG iterates. The +paper[\[1\]](#ref-cupdlpx) uses $c_0=0$. + +## Restarted reflected Halpern PDHG + +The algorithm applies reflection and Halpern anchoring to the PDHG update, +with adaptive restarts[\[1\]](#ref-cupdlpx). + +
+
+ Algorithm 1. Restarted reflected Halpern PDHG +
+
+ Input: An LP, an initial point + $(x^{0,0},y^{0,0})$, and a termination tolerance $\varepsilon$. +
+
    +
  1. + Precondition. Apply geometric mean scaling, Ruiz equilibration, + Pock–Chambolle scaling, and objective and bound scaling. +
    2. +
  2. + Initialize. Set + $$ + \eta=\frac{0.998}{\lVert A\rVert_2},\qquad + \omega^0=1,\qquad + \tau=\frac{\eta}{\omega^0},\qquad + \sigma=\eta\omega^0. + $$ +
    3. +
  3. + for + $n=0,1,\ldots$ do +
    4. +
  4. + Set $k\leftarrow 0$. +
    5. +
  5. + repeat +
    6. +
  6. +
    + PDHG update. + $$ + \begin{aligned} + \widehat x^{n,k+1} + &=\operatorname{proj}_{\mathcal X} + \left(x^{n,k}-\tau(c-A^\top y^{n,k})\right),\\[0.4em] + \widehat y^{n,k+1} + &=y^{n,k} + -\sigma A(2\widehat x^{n,k+1}-x^{n,k}) + -\sigma\operatorname{proj}_{[-u_c,-\ell_c]} + \left(\sigma^{-1}y^{n,k} + -A(2\widehat x^{n,k+1}-x^{n,k})\right). + \end{aligned} + $$ +
    +
    7. +
  7. +
    + Termination check. +
    + if + $\operatorname{KKT}(\widehat x^{n,k+1},\widehat y^{n,k+1}) + \le\varepsilon$ + then +
    +
    + return + $(\widehat x^{n,k+1},\widehat y^{n,k+1})$. +
    +
    + end if +
    +
    +
    8. +
  8. +
    + Reflection and Halpern update. + $$ + \begin{aligned} + x^{n,k+1} + &=\frac{k+1}{k+2} + \left(\gamma\bigl(2\widehat x^{n,k+1}-x^{n,k}\bigr) + +(1-\gamma)x^{n,k}\right) + +\frac{1}{k+2}x^{n,0},\\ + y^{n,k+1} + &=\frac{k+1}{k+2} + \left(\gamma\bigl(2\widehat y^{n,k+1}-y^{n,k}\bigr) + +(1-\gamma)y^{n,k}\right) + +\frac{1}{k+2}y^{n,0}. + \end{aligned} + $$ +
    +
    9. +
  9. + Set $k\leftarrow k+1$. +
    10. +
  10. + until a restart condition + holds. +
    11. +
  11. +
    + Restart. +
    Set the next anchors:
    + $$ + (x^{n+1,0},y^{n+1,0}) + =(\widehat x^{n,k},\widehat y^{n,k}). + $$ +
    + Update $\omega^{n+1}$ with the PID controller, then update $\tau$ + and $\sigma$. +
    +
    +
    12. +
  12. + end for +
    13. +
+
+ +The algorithm uses the following notation: + +- $\operatorname{proj}_{\mathcal X}$: projection onto the variable bounds, + computed by clipping each coordinate to its interval; +- $\eta$: the step size, based on the scaled matrix's spectral norm; + the [active-set step size boost](step-size.md#active-set-step-size-boost) may increase it at a restart; +- $\omega^n$: the primal weight for epoch $n$; the initial value + $\omega^0=1$ assumes the default objective and bound scaling described + under [Preconditioning](preconditioning.md); +- $\tau$ and $\sigma$: the primal and dual step sizes, respectively; +- $\gamma$: the reflection coefficient + ([`reflection_coefficient`](../reference/parameters.md#step-size-and-reflection), default `1`), + which weights the reflected point $2\widehat x^{n,k+1}-x^{n,k}$ against + the current iterate $x^{n,k}$; +- $n$ and $k$: the restart epoch and the iteration within that epoch + (here $n$ is an epoch index, rather than the variable count in the LP); +- $(x^{n,0},y^{n,0})$: the primal and dual anchors for epoch $n$; +- $(\widehat x^{n,k+1},\widehat y^{n,k+1})$: the PDHG iterate before the + reflection and Halpern update, at which the KKT residual is evaluated. + +The **KKT** (Karush–Kuhn–Tucker) check tests primal feasibility, dual +feasibility, and the primal–dual gap. See [Termination criteria](termination.md) +for the residuals and tolerances. + +The following pages describe the components: +[Base algorithm](base-algorithm.md), +[Presolve and postsolve](presolve.md), [Preconditioning](preconditioning.md), +[Step size](step-size.md), +[Adaptive restart](restart.md), [Primal weight](primal-weight.md), and +[Termination criteria](termination.md). Optional +[feasibility polishing](feasibility-polishing.md) can reduce primal and dual +feasibility residuals after the main solve. + +## References + +\[1\] Haihao Lu, Zedong Peng, and Jinwen Yang. +[*cuPDLPx: A Further Enhanced GPU-Based First-Order Solver for Linear +Programming*](https://arxiv.org/abs/2507.14051), 2025. diff --git a/docs/algorithm/preconditioning.md b/docs/algorithm/preconditioning.md new file mode 100644 index 0000000..5cd5b24 --- /dev/null +++ b/docs/algorithm/preconditioning.md @@ -0,0 +1,125 @@ +--- +description: Row, column, objective, and bound scaling in cuPDLPx. +--- + +# Preconditioning + +cuPDLPx rescales the constraint matrix, objective, and bounds to improve +numerical conditioning and accelerate convergence. + +Positive diagonal matrices $D_1$ and $D_2$ scale the rows and columns of $A$, +and two positive scalars $\theta_b$ and $\theta_c$ scale the bounds and the +objective. The scaled data are + +$$ +\begin{aligned} +\widetilde A&=D_1AD_2,\\ +\widetilde c&=\theta_c\,D_2c,\\ +(\widetilde\ell_c,\widetilde u_c)&=\theta_b\,D_1(\ell_c,u_c),\\ +(\widetilde\ell_v,\widetilde u_v)&=\theta_b\,D_2^{-1}(\ell_v,u_v). +\end{aligned} +$$ + +To compute $D_1$ and $D_2$, cuPDLPx applies geometric mean scaling as used +in Gurobi[\[1\]](#ref-gurobi-scaleflag), followed by Ruiz +equilibration and Pock–Chambolle scaling as used in +PDLP[\[2\]](#ref-pdlp). Each stage rescales $A$ and accumulates +its row and column factors in $D_1$ and $D_2$. Empty rows and columns have +scaling factor $1$. + +The final [objective and bound scaling](#objective-and-bound-scaling) stage +follows HPR-LP[\[3\]](#ref-hpr-lp) and determines $\theta_b$ and +$\theta_c$. Both are $1$ when this stage is disabled. + +## Geometric mean scaling + +One geometric mean scaling pass[\[1\]](#ref-gurobi-scaleflag) first scales the rows, +then the columns. For each row, the factor is the reciprocal geometric mean +of its smallest and largest nonzero entry magnitudes: + +$$ +d_{1,i}= +\left(\min_{j:\,a_{ij}\ne 0}|a_{ij}| +\;\max_{j:\,a_{ij}\ne 0}|a_{ij}|\right)^{-1/2}. +$$ + +After applying the row factors, compute and apply the column factors from the +row-scaled matrix: + +$$ +d_{2,j}= +\left(\min_{i:\,a_{ij}\ne 0}|a_{ij}| +\;\max_{i:\,a_{ij}\ne 0}|a_{ij}|\right)^{-1/2}. +$$ + +The number of row–column passes defaults to `12`; setting it to `0` +disables this stage. + +## Ruiz equilibration scaling + +One Ruiz iteration[\[2\]](#ref-pdlp) computes the row and column factors +from the same matrix and applies them together: + +$$ +d_{1,i}=\left(\max_{j}|a_{ij}|\right)^{-1/2}, +\qquad +d_{2,j}=\left(\max_{i}|a_{ij}|\right)^{-1/2}. +$$ + +Repeating the iteration balances the row and column $\ell_\infty$ norms. +The number of iterations defaults to `10`. + +## Pock–Chambolle scaling + +The Pock–Chambolle factors[\[2\]](#ref-pdlp) are computed from the +same matrix and applied together: + +$$ +d_{1,i}=\left(\sum_{j}|a_{ij}|^\alpha\right)^{-1/2}, +\qquad +d_{2,j}=\left(\sum_{i}|a_{ij}|^{2-\alpha}\right)^{-1/2}. +$$ + +At the default $\alpha=1$, they are the reciprocal square roots of the +corresponding row and column $\ell_1$ norms. This stage is enabled by +default. + +## Objective and bound scaling + +The final stage runs after the diagonal scaling and multiplies the scaled +data by two scalars[\[3\]](#ref-hpr-lp): + +$$ +\theta_b=\frac{1}{1+\|b\|_2}, +\qquad +\theta_c=\frac{1}{1+\|c\|_2}. +$$ + +Here $b$ collects the finite entries of the scaled $\ell_c$ and $u_c$, +counting each equality bound once, and $c$ is the scaled +objective vector. The constraint and variable bounds are multiplied by +$\theta_b$, and the objective vector by $\theta_c$. +This stage is enabled by default. Whether it runs also decides the +[initial primal weight](primal-weight.md#initial-weight). + +## Parameters + +See [Parameters](../reference/parameters.md#scaling-and-preprocessing) for +scaling settings. + +## References + +\[1\] Gurobi Optimization. +[*ScaleFlag*](https://docs.gurobi.com/projects/optimizer/en/current/reference/parameters.html#scaleflag). +Gurobi Optimizer Reference Manual. + +\[2\] David Applegate, Mateo Díaz, Oliver Hinder, Haihao Lu, +Miles Lubin, Brendan O'Donoghue, and Warren Schudy. +[*Practical Large-Scale Linear Programming Using Primal-Dual Hybrid +Gradient*](https://proceedings.neurips.cc/paper/2021/hash/a8fbbd3b11424ce032ba813493d95ad7-Abstract.html). +*NeurIPS*, 2021. + +\[3\] Kaihuang Chen, Defeng Sun, Yancheng Yuan, Guojun Zhang, +and Xinyuan Zhao. +[*HPR-LP: An Implementation of an HPR Method for Solving Linear +Programming*](https://arxiv.org/abs/2408.12179), 2024. diff --git a/docs/algorithm/presolve.md b/docs/algorithm/presolve.md new file mode 100644 index 0000000..51126f0 --- /dev/null +++ b/docs/algorithm/presolve.md @@ -0,0 +1,47 @@ +--- +description: LP reductions with PSLP and recovery of the original model's solution. +--- + +# Presolve and postsolve + +cuPDLPx uses [PSLP](https://github.com/dance858/PSLP)[\[1\]](#ref-pslp) +to reduce the LP before applying r2HPDHG. Postsolve recovers the +primal and dual variables for the original LP. + +## Presolve + +Presolve is enabled by default and applies reductions such as: + +| Reduction | What it does | +| --- | --- | +| Fixed-variable elimination | Substitute variables with equal lower and upper bounds into the constraints and objective. | +| Redundant-constraint removal | Remove constraints already implied by the remaining model. | +| Bound tightening | Use constraints and existing bounds to derive tighter variable bounds. | + +Reducing the number of rows, columns, and nonzeros can lower the cost of +the sparse matrix–vector products in each PDHG iteration. PSLP records +these transformations for postsolve. + +cuPDLPx skips the main iteration if presolve solves the LP or returns +an infeasible or infeasible-or-unbounded status. + +!!! note + + When presolve is enabled, [termination checks](termination.md) are performed + on the presolved model, not the original model. + +## Postsolve + +After the reduced LP is solved and any feasibility polishing is complete, +PSLP recovers the primal variables, dual variables, and dual slacks for the +original LP. cuPDLPx returns these recovered vectors. + +!!! warning "Presolve and warm starts" + + User-provided warm starts require presolve to be disabled. + +## References + +\[1\] Daniel Cederberg and Stephen Boyd. +[*Presolving for GPU-Accelerated First-Order LP +Solvers*](https://arxiv.org/abs/2604.23951), 2026. diff --git a/docs/algorithm/primal-weight.md b/docs/algorithm/primal-weight.md new file mode 100644 index 0000000..54d1320 --- /dev/null +++ b/docs/algorithm/primal-weight.md @@ -0,0 +1,119 @@ +--- +description: Primal and dual step sizes and the adaptive primal-weight update in cuPDLPx. +--- + +# Primal weight + +The primal weight $\omega>0$ determines the relative sizes of the primal +and dual steps[\[1\]](#ref-pdlp): + +$$ +\tau=\frac{\eta}{\omega},\qquad +\sigma=\eta\omega. +$$ + +Increasing $\omega$ decreases the primal step size and increases the dual +step size. cuPDLPx updates $\omega$ at each restart to balance the scaled +primal and dual displacements. + +## Initial weight + +With +[objective and bound scaling](preconditioning.md#objective-and-bound-scaling) +enabled, the initial weight is $\omega^0=1$. Otherwise, + +$$ +\omega^0=\frac{1+\lVert c\rVert}{1+\lVert b\rVert}, +$$ + +where $c$ is the objective vector and $b$ collects the finite constraint +bounds before preconditioning, counting each equality bound once. The norm +is selected by +[`optimality_norm`](../reference/parameters.md#accuracy-and-termination). + +## PID controller + +cuPDLPx uses a PID controller with a discounted integral +term[\[2\]](#ref-cupdlpx). + +For epoch $n$, let $(x^{n,0},y^{n,0})$ be the anchor and +$(\widehat x^{n,k},\widehat y^{n,k})$ the PDHG iterate at the restart. +The imbalance between the scaled displacements is + +$$ +e^n=\log\left( +\frac{\sqrt{\omega^n}\lVert \widehat x^{n,k}-x^{n,0}\rVert_2} +{(1/\sqrt{\omega^n})\lVert \widehat y^{n,k}-y^{n,0}\rVert_2} +\right). +$$ + +Positive $e^n$ indicates a larger scaled primal displacement; negative +$e^n$ indicates a larger scaled dual displacement. The update is + +$$ +\log\omega^{n+1}=\log\omega^n- +\left[ +K_Pe^n ++K_I\sum_{i=1}^{n}\rho^{n-i}e^i ++K_D(e^n-e^{n-1}) +\right]. +$$ + +The proportional, integral, and derivative terms use the current imbalance, +its discounted history, and its change since the previous epoch. +Updating $\log\omega$ preserves $\omega>0$. + +cuPDLPx uses + +$$ +K_P=0.99,\qquad K_I=0.01,\qquad K_D=0,\qquad \rho=0.3 +$$ + +by default. The weight and step sizes remain constant within each epoch. + +## Safeguard + +cuPDLPx resets $\omega$ and skips the PID update if either of the following +conditions holds[\[3\]](#ref-hpr-lp): + +- either displacement norm, $\lVert\widehat x^{n,k}-x^{n,0}\rVert_2$ or + $\lVert\widehat y^{n,k}-y^{n,0}\rVert_2$, lies outside + $[10^{-16},10^{12}]$; +- the ratio of the relative dual residual to the relative primal residual + lies outside $[10^{-8},10^{8}]$. + +The reset restores $\omega_{\mathrm{best}}$ and clears the controller's +integral and derivative state. This weight was produced at the restart +where the relative primal and dual residuals, $\delta_p^{\,n}$ and +$\delta_d^{\,n}$, were closest on a logarithmic scale: + +$$ +\omega_{\mathrm{best}}=\omega^{m+1}, +\qquad +m=\operatorname*{arg\,min}_{n} +\left|\log_{10}\frac{\delta_d^{\,n}}{\delta_p^{\,n}}\right|. +$$ + +The stored weight is initialized to $\omega^0$ and updated when a restart +improves this measure. + +## Parameters + +See [Parameters](../reference/parameters.md#adaptive-restart) for the PID gains. + +## References + +\[1\] David Applegate, Mateo Díaz, Oliver Hinder, Haihao Lu, +Miles Lubin, Brendan O'Donoghue, and Warren Schudy. +[*Practical Large-Scale Linear Programming Using Primal-Dual Hybrid +Gradient*](https://proceedings.neurips.cc/paper/2021/hash/a8fbbd3b11424ce032ba813493d95ad7-Abstract.html). +*NeurIPS*, 2021. + +\[2\] Haihao Lu, Zedong Peng, and Jinwen Yang. +[*cuPDLPx: A Further Enhanced GPU-Based First-Order Solver for Linear +Programming*](https://arxiv.org/abs/2507.14051), 2025. + +\[3\] Kaihuang Chen, Defeng Sun, Yancheng Yuan, Guojun Zhang, +and Xinyuan Zhao. +[*HPR-LP: An Implementation of an HPR Method for Solving Linear +Programming*](https://arxiv.org/abs/2408.12179), 2024. diff --git a/docs/algorithm/restart.md b/docs/algorithm/restart.md new file mode 100644 index 0000000..ee1c5a9 --- /dev/null +++ b/docs/algorithm/restart.md @@ -0,0 +1,94 @@ +--- +description: Fixed-point-error restart criteria and epoch transitions in cuPDLPx. +--- + +# Adaptive restart + +cuPDLPx restarts the Halpern iteration based on the fixed-point error and +the number of iterations since the last restart[\[1\]](#ref-lu-yang),[\[2\]](#ref-cupdlpx). +The iterations between restarts form an epoch. + +## Fixed-point error + +For the PDHG operator $\mathcal T$, define the fixed-point error of a point +$(x,y)$ as + +$$ +r(x,y)=\left\lVert (x,y)-\mathcal T(x,y)\right\rVert_P, +$$ + +where $P$ is the canonical PDHG metric: + +$$ +P= +\begin{bmatrix} +\dfrac{\omega}{\eta}I & A^\top \\ +A & \dfrac{1}{\eta\omega}I +\end{bmatrix}\succ0. +$$ + +PDLP uses the normalized duality gap for restarting, while the earlier +cuPDLP uses the KKT error[\[2\]](#ref-cupdlpx). + +Let $(x^{n,0},y^{n,0})$ be the anchor of epoch $n$, +$(x^{n,k},y^{n,k})$ its current iterate, and $N$ the total iteration count. A +restart occurs when any of the following conditions is satisfied. + +## Sufficient reduction + +$$ +r(x^{n,k},y^{n,k})\le +\beta_{\mathrm{sufficient}}r(x^{n,0},y^{n,0}). +$$ + +This condition requires the fixed-point error to fall to a fraction +$\beta_{\mathrm{sufficient}}$ of its value at the anchor. The default is `0.2`. + +## Necessary reduction and local increase + +$$ +\begin{aligned} +r(x^{n,k},y^{n,k})&\le +\beta_{\mathrm{necessary}}r(x^{n,0},y^{n,0}),\\ +r(x^{n,k},y^{n,k})&>r(x^{n,k'},y^{n,k'}). +\end{aligned} +$$ + +Here $(x^{n,k'},y^{n,k'})$ is the iterate at the previous termination check, +[`termination_evaluation_frequency`](../reference/parameters.md#limits-and-logging) +iterations earlier. The error must be sufficiently below its value at the +anchor but larger than at the previous check. The default +$\beta_{\mathrm{necessary}}$ is `0.5`. + +## Artificial epoch limit + +$$ +k\ge\beta_{\mathrm{artificial}}N. +$$ + +This condition limits the epoch length relative to the total iteration count. +The default $\beta_{\mathrm{artificial}}$ is `0.36`. + +## Restart operation + +At a restart, cuPDLPx: + +1. replaces the anchor with the latest PDHG iterate; +2. resets the local Halpern iteration count; +3. updates the [primal weight](primal-weight.md). + +The restart criteria are evaluated at every termination check. + +## Parameters + +See [Parameters](../reference/parameters.md#adaptive-restart) for restart settings. + +## References + +\[1\] Haihao Lu and Jinwen Yang. +[*Restarted Halpern PDHG for Linear Programming*](https://arxiv.org/abs/2407.16144), +2024. + +\[2\] Haihao Lu, Zedong Peng, and Jinwen Yang. +[*cuPDLPx: A Further Enhanced GPU-Based First-Order Solver for Linear +Programming*](https://arxiv.org/abs/2507.14051), 2025. diff --git a/docs/algorithm/step-size.md b/docs/algorithm/step-size.md new file mode 100644 index 0000000..e045d64 --- /dev/null +++ b/docs/algorithm/step-size.md @@ -0,0 +1,198 @@ +--- +description: Spectral-norm-based step size and active-set step size boost in cuPDLPx. +--- + +# Step size + +cuPDLPx initializes the step size from the spectral norm of the scaled +constraint matrix. The [active-set step size boost](#active-set-step-size-boost) +can increase it near a solution. + +## Initial constant step size + +For the scaled constraint matrix $A$, the initial step size +is[\[1\]](#ref-cupdlpx) + +$$ +\eta=\frac{0.998}{\lVert A\rVert_2}. +$$ + +The spectral norm is estimated by power iteration on $AA^\top$ after +[preconditioning](preconditioning.md). +The primal and dual step sizes are + +$$ +\tau=\frac{\eta}{\omega},\qquad +\sigma=\eta\omega. +$$ + +The factor `0.998` provides a safety margin relative to the estimated norm. +The canonical PDHG metric is positive definite when +$\eta<1/\lVert A\rVert_2$. At each restart, the +[primal-weight controller](primal-weight.md) can change $\tau$ and $\sigma$ +while preserving $\tau\sigma=\eta^2$ for a fixed $\eta$. + +## Active-set step size boost + +The initial step size $\eta^0=0.998/\lVert A\rVert_2$ uses the full constraint +matrix. Near a solution, inactive constraints and variables that remain at +their bounds can permit a larger step size. The active-set step size boost +(ASB) retains potentially binding constraints in $B$ and variables not +identified as persistently clamped at a bound in $F$. The tests below +define these sets. The resulting submatrix $A_{B,F}$ satisfies + +$$ +\lVert A_{B,F}\rVert_2\le\lVert A\rVert_2. +$$ + +ASB estimates $\lVert A_{B,F}\rVert_2$ and uses it to propose a larger step +size at an adaptive restart. A divergence check can trigger rollback and +start a new epoch. ASB is enabled by default. + +### Active set + +The sets $B$ and $F$ retain constraints and variables that do not meet the +exclusion tests over the [trailing window](#trailing-window). +These tests use the PDHG update from $(x,y)$ to $(\widehat x,\widehat y)$: + +$$ +\begin{aligned} +\widehat x +&=\operatorname{proj}_{[\ell_v,u_v]} +\left(x-\tau(c-A^\top y)\right),\\[0.4em] +\widehat y +&=y-\sigma A(2\widehat x-x) +-\sigma\operatorname{proj}_{[-u_c,-\ell_c]} +\left(\sigma^{-1}y-A(2\widehat x-x)\right). +\end{aligned} +$$ + +**Variables.** The primal projection gives the dual-slack estimate + +$$ +\tilde r=c-A^\top y+\tau^{-1}(\widehat x-x). +$$ + +A positive $\tilde r_i$ indicates projection onto the lower bound; a +negative value indicates projection onto the upper bound. Variable $i$ +meets the exclusion test when $\ell_{v,i}=u_{v,i}$ or when + +$$ +\begin{aligned} +\tilde r_i&>\varepsilon_r\max\{1,|c_i|\} +&&\text{if } \ell_{v,i}>-\infty, +\qquad\text{or}\\ +\tilde r_i&<-\varepsilon_r\max\{1,|c_i|\} +&&\text{if } u_{v,i}<\infty. +\end{aligned} +$$ + +**Constraints.** The dual update projects + +$$ +\tilde s=\sigma^{-1}y-A(2\widehat x-x) +$$ + +onto $[-u_c,-\ell_c]$. It gives $\widehat y_j=0$ when $\tilde s_j$ lies +in this interval. Constraint $j$ meets the exclusion test when +$\ell_{c,j}\varepsilon_c\max\{1,|\tilde s_j|,|u_{c,j}|\} +&&\text{if } u_{c,j}<\infty,\\ +-\ell_{c,j}-\tilde s_j&>\varepsilon_c\max\{1,|\tilde s_j|,|\ell_{c,j}|\} +&&\text{if } \ell_{c,j}>-\infty. +\end{aligned} +$$ + +Both tolerances default to `1e-8`. The submatrix $A_{B,F}$ contains the rows in +$B$ and the columns in $F$. + +### Activation + +#### Trailing window + +The exclusion tests run at every termination check. An index is removed +from $B$ or $F$ only if it meets the exclusion test at every check over the +last `10000` iterations by default. It is restored as soon as it fails the test. + +#### Activation tolerance + +ASB activates when the relative primal residual, relative dual residual, +and relative primal–dual gap are all below the activation tolerance +($10^{-4}$ by default), and at least one index has been excluded. +See [Termination criteria](termination.md#optimality-conditions) +for the residual definitions. + +### Step-size increase + +After activation, ASB uses adaptive restarts to estimate + +$$ +\hat\sigma\approx\lVert A_{B,F}\rVert_2 +$$ + +by power iteration, setting entries outside $B$ and $F$ to zero after each +product and warm-starting from the previous eigenvector. + +ASB re-estimates the norm when the accumulated additions and removals +reach `1%` of the current $|B|+|F|$ by default. It skips the estimate if the +step-size limit rules out an increase, and stops power iteration early if +the running estimate does so. These decisions are retained until the same +change threshold is reached. + +For a positive norm estimate, the proposed step size and acceptance condition are + +$$ +\eta^{+} +=\max\left\{\eta^0,\ +\min\left\{\frac{\alpha}{\hat\sigma},\ \eta_{\mathrm{ceil}}\right\}\right\}, +\qquad +\eta\leftarrow\eta^{+} +\ \text{ if }\ \eta^{+}>\eta +\ \text{ and }\ \eta^{+}\ge\rho\,\eta. +$$ + +Here $\alpha$ is a safety factor (default `0.9`), $\rho$ is the minimum +increase ratio (default `1.1`), and the upper limit +$\eta_{\mathrm{ceil}}$ is initially $\infty$. +When an increase is accepted, ASB saves the restart iterate, primal weight, +and controller state for rollback. + +### Divergence protection + +While $\eta>\eta^0$, ASB monitors the +[fixed-point error](restart.md#fixed-point-error) $r$, with the metric $P$ +evaluated at the initial step size. The divergence check triggers when the +error or residuals are nonfinite, or when + +$$ +r(z^k)>(1+\delta)\,r(z^0). +$$ + +Here $z^0$ is the first iterate of the current epoch, and $\delta$ defaults +to `0.05`. A rollback + +- restores the saved iterate, primal weight, and controller state; +- resets the step size to $\eta^0$ and discards $\hat\sigma$; +- resets the local iteration count; +- sets $\eta_{\mathrm{ceil}}$ to `0.7` times the rejected step size by default. + +By default, ASB is disabled for the rest of the solve after two rollbacks. +The C result and command-line output report step-size increases, rollbacks, +and power iteration counts. + +## Parameters + +See [Step size and reflection](../reference/parameters.md#step-size-and-reflection) +for power iteration settings and +[Active-set step size boost](../reference/parameters.md#active-set-step-size-boost) for +ASB settings. + +## References + +\[1\] Haihao Lu, Zedong Peng, and Jinwen Yang. +[*cuPDLPx: A Further Enhanced GPU-Based First-Order Solver for Linear +Programming*](https://arxiv.org/abs/2507.14051), 2025. diff --git a/docs/algorithm/termination.md b/docs/algorithm/termination.md new file mode 100644 index 0000000..b518b66 --- /dev/null +++ b/docs/algorithm/termination.md @@ -0,0 +1,130 @@ +--- +description: Primal feasibility, dual feasibility, and objective-gap termination checks used by cuPDLPx. +--- + +# Termination criteria + +cuPDLPx checks optimality using relative primal and dual residuals and the +primal–dual gap. These quantities are evaluated on the unscaled LP, after +presolve if enabled. The formulas below use the default $\ell_2$ norm. +Selecting $\ell_\infty$ replaces the norms in both residual numerators +and denominators. + +## Optimality conditions + +The optimality test requires all three quantities to satisfy the following +bounds[\[1\]](#ref-cupdlpx). Here $x$ denotes the primal variables, +$y$ the dual multipliers, $r$ the dual slacks, and $p$ the +[support function](index.md#dual-problem). The vector $b$ collects the finite +constraint bounds, counting each equality bound once. + +- **Primal feasibility** + + $$ + \frac{ + \left\lVert Ax-\operatorname{proj}_{[\ell_c,u_c]}(Ax)\right\rVert_2 + }{ + 1+\lVert b\rVert_2 + } + <\varepsilon_{\mathrm{feas}}. + $$ + +- **Dual feasibility** + + $$ + \frac{\left\lVert c-A^\top y-r\right\rVert_2} + {1+\lVert c\rVert_2} + <\varepsilon_{\mathrm{feas}}. + $$ + +- **Relative objective gap** + + $$ + \frac{ + \left|c^\top x+p(-y;\ell_c,u_c)+p(-r;\ell_v,u_v)\right| + }{ + 1+\left|c^\top x+c_0\right| + +\left|-p(-y;\ell_c,u_c)-p(-r;\ell_v,u_v)+c_0\right| + } + <\varepsilon_{\mathrm{opt}}. + $$ + +Variable-bound feasibility and the dual-slack sign conditions are satisfied +by the primal projection and Moreau recovery, respectively. Both tolerances +default to `1e-4`. + +## Dual-slack recovery + +At each termination check, cuPDLPx recovers the dual-slack vector $r$ +using the Moreau identity: + +$$ +\widehat r^{\,k+1} +=\tau^{-1}\left(\widehat x^{k+1}-x^{k}\right)+(c-A^\top y^k). +$$ + +The recovered vector is used to evaluate dual feasibility and the +primal–dual gap at the PDHG iterate $\widehat x^{k+1}$. The returned dual +slacks are computed separately from $c-A^\top y$, with each component +restricted to the sign permitted by its variable bounds. + +## Unscaling + +The iterations and dual-slack recovery use the +[scaled LP](preconditioning.md), with $\widetilde A=D_1AD_2$ and +$\widetilde c=\theta_cD_2c$. The variables map to the unscaled LP as + +$$ +x=\theta_b^{-1}D_2\widetilde x, +\qquad +y=\theta_c^{-1}D_1\widetilde y, +\qquad +r=\theta_c^{-1}D_2^{-1}\widetilde r. +$$ + +Residuals are unscaled before taking norms: + +$$ +\begin{aligned} +Ax-\operatorname{proj}_{[\ell_c,u_c]}(Ax) +&=\theta_b^{-1}D_1^{-1} +\Bigl(\widetilde A\widetilde x +-\operatorname{proj}_{[\widetilde\ell_c,\widetilde u_c]} +(\widetilde A\widetilde x)\Bigr),\\[0.4em] +c-A^\top y-r +&=\theta_c^{-1}D_2^{-1} +\bigl(\widetilde c-\widetilde A^\top\widetilde y-\widetilde r\bigr). +\end{aligned} +$$ + +The objective values are recovered from +$c^\top x=(\theta_b\theta_c)^{-1}\,\widetilde c^{\,\top}\widetilde x$ and its +dual counterpart. The denominators $1+\lVert b\rVert_2$ and +$1+\lVert c\rVert_2$ use the unscaled data. When presolve is enabled, +[postsolve](presolve.md#postsolve) then recovers the original model's +solution vectors. + +## Check frequency + +cuPDLPx checks optimality, infeasibility, and restart criteria every `200` +iterations by default. Time and iteration limits are checked at the same +interval, so the iteration count can exceed its limit by less than one +check interval. + +An iterate returned after a time or iteration limit may not satisfy the +optimality conditions. See [Results and status](../getting-started/results.md) +for interpretation and [CUDA Graphs](../implementation/cuda-graphs.md) +for execution between checks. + +## Parameters + +See [Parameters](../reference/parameters.md#accuracy-and-termination) for +the tolerances and norm, and +[Limits and logging](../reference/parameters.md#limits-and-logging) for the +check interval. + +## References + +\[1\] Haihao Lu, Zedong Peng, and Jinwen Yang. +[*cuPDLPx: A Further Enhanced GPU-Based First-Order Solver for Linear +Programming*](https://arxiv.org/abs/2507.14051), 2025. diff --git a/docs/assets/cuda-graphs-launch-overhead-dark.svg b/docs/assets/cuda-graphs-launch-overhead-dark.svg new file mode 100644 index 0000000..b307b01 --- /dev/null +++ b/docs/assets/cuda-graphs-launch-overhead-dark.svg @@ -0,0 +1,80 @@ + + CUDA Graphs reduce repeated CPU kernel-launch overhead + Separate CPU launches leave latency gaps between short GPU kernels, while one graph launch submits the recorded kernels together. + + + + + + + + CPU Time + GPU Time + + + + + + + + Launch A + Launch B + Launch C + Launch D + Launch E + + + + + + + + + + + + + A + B + C + D + E + launch latency + + + + time + + + + Build + graph + + + Launch graph + + + + + + + + A + B + C + D + E + + + time saved + + diff --git a/docs/assets/cuda-graphs-launch-overhead.svg b/docs/assets/cuda-graphs-launch-overhead.svg new file mode 100644 index 0000000..6b27a66 --- /dev/null +++ b/docs/assets/cuda-graphs-launch-overhead.svg @@ -0,0 +1,80 @@ + + CUDA Graphs reduce repeated CPU kernel-launch overhead + Separate CPU launches leave latency gaps between short GPU kernels, while one graph launch submits the recorded kernels together. + + + + + + + + CPU Time + GPU Time + + + + + + + + Launch A + Launch B + Launch C + Launch D + Launch E + + + + + + + + + + + + + A + B + C + D + E + launch latency + + + + time + + + + Build + graph + + + Launch graph + + + + + + + + A + B + C + D + E + + + time saved + + diff --git a/docs/assets/mit-favicon.svg b/docs/assets/mit-favicon.svg new file mode 100644 index 0000000..451d380 --- /dev/null +++ b/docs/assets/mit-favicon.svg @@ -0,0 +1,6 @@ + diff --git a/docs/citation.md b/docs/citation.md new file mode 100644 index 0000000..5f5cbb9 --- /dev/null +++ b/docs/citation.md @@ -0,0 +1,127 @@ +--- +description: Recommended papers, BibTeX entries, and related projects for cuPDLPx. +--- + +# Citation + +If you use cuPDLPx in published work, please cite the computational paper. If +your work also relies on the restarted Halpern PDHG method or its convergence +theory, please cite the theoretical paper as well. + +## Computational paper + +Algorithmic enhancements, GPU implementation, and numerical results are +described in this paper: + +```bibtex +@article{lu2025cupdlpx, + title = {{cuPDLPx}: A Further Enhanced GPU-Based First-Order Solver + for Linear Programming}, + author = {Lu, Haihao and Peng, Zedong and Yang, Jinwen}, + journal = {arXiv preprint arXiv:2507.14051}, + year = {2025}, + url = {https://arxiv.org/abs/2507.14051} +} +``` + +## Theoretical paper + +Restarted Halpern PDHG and its reflected variant are developed in this paper: + +```bibtex +@article{lu2024restarted, + title = {Restarted Halpern PDHG for Linear Programming}, + author = {Lu, Haihao and Yang, Jinwen}, + journal = {arXiv preprint arXiv:2407.16144}, + year = {2024}, + url = {https://arxiv.org/abs/2407.16144} +} +``` + +## Background + +### Papers + +- David Applegate et al., [*Practical Large-Scale Linear Programming Using + Primal-Dual Hybrid + Gradient*](https://proceedings.neurips.cc/paper/2021/hash/a8fbbd3b11424ce032ba813493d95ad7-Abstract.html), + NeurIPS 2021. +- David Applegate et al., [*PDLP: A Practical First-Order Method for + Large-Scale Linear Programming*](https://arxiv.org/abs/2501.07018), 2025. +- Haihao Lu and Jinwen Yang, [*cuPDLP.jl: A GPU Implementation of Restarted + Primal-Dual Hybrid Gradient for Linear Programming in + Julia*](https://arxiv.org/abs/2311.12180), 2023. +- Haihao Lu et al., [*cuPDLP-C: A Strengthened Implementation of cuPDLP for + Linear Programming by C language*](https://arxiv.org/abs/2312.14832), 2023. +- Haihao Lu and Jinwen Yang, [*An Overview of GPU-based First-Order Methods + for Linear Programming and Extensions*](https://arxiv.org/abs/2506.02174), + 2025. +- Daniel Cederberg and Stephen Boyd, [*Presolving for GPU-Accelerated + First-Order LP Solvers*](https://arxiv.org/abs/2604.23951), 2026. + +### Blogs + +- Cara Touretzky, Robert Luce, and David Torres Sanchez, [*Using GPUs to Solve + LPs: What's in It for + Me?*](https://www.gurobi.com/resources/blog/using-gpus-to-solve-lps-what-s-in-it-for-me), + Gurobi blog, 2025. +- Cara Touretzky and Robert Luce, [*Introducing Gurobi's First GPU-Accelerated + Solver: Test a Beta Version of Gurobi's PDHG Implementation on NVIDIA's GPU + Hardware*](https://www.gurobi.com/resources/blog/introducing-gurobi-s-first-gpu-accelerated-solver-test-a-beta-version-of-gurobi-s-pdhg-implementation-on-nvidia-s-gpu-hardware), + Gurobi blog, 2025. +- Imre Pólik, [*GPU Acceleration of the Hybrid Gradient Algorithm in FICO + Xpress*](https://www.fico.com/blogs/gpu-acceleration-hybrid-gradient-algorithm-fico-xpress), + FICO blog, 2025. +- Nicolas Blin, [*Accelerate Large Linear Programming Problems with NVIDIA + cuOpt*](https://developer.nvidia.com/blog/accelerate-large-linear-programming-problems-with-nvidia-cuopt/), + NVIDIA Technical Blog, 2024. +- Artelys, [*Artelys Knitro 15.0: New Tools for Your Large-Scale + Models*](https://www.artelys.com/news/knitro-15-0-new-tools-for-your-large-scale-models/), + Artelys news, 2025. +- [Scaling up linear programming with PDLP](https://research.google/blog/scaling-up-linear-programming-with-pdlp/), + a Google Research article by Haihao Lu and David Applegate. +- [Mathematical background for PDLP](https://developers.google.com/optimization/lp/pdlp_math), + the OR-Tools reference for PDLP formulations, residuals, rescaling, and + infeasibility certificates. + +## Related projects + +### Open-source solvers + +- [cuPDLP.jl](https://github.com/jinwen-yang/cuPDLP.jl), the earlier Julia GPU solver. +- [cuPDLP-C](https://github.com/COPT-Public/cuPDLP-C), the C implementation of + cuPDLP. +- [PDQP.jl](https://github.com/jinwen-yang/PDQP.jl), a Julia first-order solver + for convex quadratic programming on CPUs and NVIDIA GPUs. +- [PDHCG](https://github.com/Lhongpei/PDHCG), a GPU-accelerated first-order + solver for convex quadratic and conic quadratic programming. +- [HPR-LP-C](https://github.com/PolyU-IOR/HPR-LP-C), a C GPU solver for LP + based on the Halpern Peaceman–Rachford method. +- [cuOpt](https://github.com/NVIDIA/cuopt), NVIDIA's open-source GPU-accelerated + optimization library for LP, MIP, and vehicle routing. +- [HiGHS](https://github.com/ERGO-Code/HiGHS), an open-source solver for LP, + MIP, and QP that includes a PDLP-based first-order LP solver. +- [D-PDLP](https://github.com/Lhongpei/D-PDLP), a distributed LP solver built + on cuPDLPx for execution across multiple GPUs. +- [CoolPDLP.jl](https://github.com/JuliaDecisionFocusedLearning/CoolPDLP.jl), + a Julia implementation of PDLP and its variants with support for CPUs, + multiple GPU architectures, and batched solves. + +### Commercial solvers (implementing PDLP) + +- [Gurobi](https://www.gurobi.com/) +- [FICO Xpress](https://www.fico.com/en/products/fico-xpress-optimization) +- [Knitro](https://www.artelys.com/solvers/knitro/) +- [COPT](https://www.shanshu.ai/copt) + +### Benchmarks + +- [MIPLIB 2017](https://miplib.zib.de/), a source of LP relaxation benchmark instances. +- [Mittelmann benchmarks](https://plato.asu.edu/bench.html), independent + benchmarks and test sets for optimization software. + +## Acknowledgements + +The development of cuPDLPx is partially supported by AFOSR Grant No. +FA9550-24-1-0051, ONR Grant No. N000142412735, and the NVIDIA Academic Grant +Program. diff --git a/docs/getting-started/index.md b/docs/getting-started/index.md new file mode 100644 index 0000000..bd87dca --- /dev/null +++ b/docs/getting-started/index.md @@ -0,0 +1,26 @@ +--- +description: Choose a Python, Julia, command-line, or C interface for cuPDLPx. +--- + +# Get started + +cuPDLPx provides Python, Julia, command-line, and C interfaces. Choose an +interface for installation instructions and examples. + +## Hardware requirements + +cuPDLPx requires a supported GPU and toolkit: + +| GPU | Toolkit | +| --- | --- | +| NVIDIA GPU | CUDA 12.4+ | +| AMD GPU | ROCm 7.2+ | + +## Interfaces + +| Interface | Best for | +| --- | --- | +| [Python](../guides/python.md) | Applications using NumPy or SciPy data | +| [Julia](../guides/julia.md) | JuMP and MathOptInterface workflows | +| [Command line](../guides/command-line.md) | Solving MPS files from a shell | +| [C](../guides/c-api.md) | Embedding cuPDLPx in native applications | diff --git a/docs/getting-started/log-interpretation.md b/docs/getting-started/log-interpretation.md new file mode 100644 index 0000000..b1ed692 --- /dev/null +++ b/docs/getting-started/log-interpretation.md @@ -0,0 +1,235 @@ +--- +description: Interpret cuPDLPx progress tables, residuals, timings, and feasibility-polishing logs. +--- + +# Log interpretation + +All interfaces use the same solver log. Logging is controlled by +`OutputFlag` in Python, `verbose` in Julia and C, and `--verbose` or +`--quiet` on the command line. + +The example uses CUDA with feasibility polishing disabled. Timings depend +on the problem, GPU, and system configuration. + +??? example "Complete solver log" + + ```text + --------------------------------------------------------------------------------------- + cuPDLPx v0.3.0 + A GPU-Accelerated First-Order LP Solver + (c) Haihao Lu, Massachusetts Institute of Technology, 2025 + --------------------------------------------------------------------------------------- + Problem: 17013 rows, 200 columns, 104811 nonzeros + Settings: + iter_limit : 2147483647 + time_limit : 3600.00 sec + eps_opt : 1.0e-04 + eps_feas : 1.0e-04 + spmv_backend : cusparseSpMVOp (auto) + + Running presolver (PSLP v0.0.8)... + status : REDUCED + presolve time : 0.00853 sec + reduced problem : 16997 rows, 200 columns, 104747 nonzeros + + Preconditioning + Ruiz scaling (10 iterations) + Pock-Chambolle scaling (alpha=1.0000) + Bound-objective scaling + + --------------------------------------------------------------------------------------- + runtime | objective | absolute residuals | relative residuals + iter time | pr obj du obj | pr res du res gap | pr res du res gap + --------------------------------------------------------------------------------------- + 0 0.0e+00 | 0.0e+00 0.0e+00 | 0.0e+00 0.0e+00 0.0e+00 | 0.0e+00 0.0e+00 0.0e+00 + 200 0.0e+00 | -1.3e+02 -1.2e+02 | 1.1e+00 2.5e-01 9.7e+00 | 8.5e-03 1.7e-02 3.9e-02 + 400 8.7e-03 | -1.2e+02 -1.2e+02 | 2.3e-01 7.6e-02 9.2e-01 | 1.7e-03 5.0e-03 3.8e-03 + 600 1.2e-02 | -1.2e+02 -1.2e+02 | 4.6e-02 8.5e-02 5.2e-03 | 3.5e-04 5.6e-03 2.1e-05 + 800 1.5e-02 | -1.2e+02 -1.2e+02 | 2.5e-02 4.7e-02 3.3e-02 | 1.9e-04 3.1e-03 1.4e-04 + 1000 1.8e-02 | -1.2e+02 -1.2e+02 | 1.2e-02 3.4e-02 7.3e-03 | 8.8e-05 2.2e-03 3.0e-05 + 1200 2.1e-02 | -1.2e+02 -1.2e+02 | 4.3e-03 2.5e-02 7.6e-04 | 3.2e-05 1.6e-03 3.1e-06 + 1400 2.4e-02 | -1.2e+02 -1.2e+02 | 3.1e-03 1.9e-02 1.4e-03 | 2.4e-05 1.3e-03 5.9e-06 + 1600 2.7e-02 | -1.2e+02 -1.2e+02 | 3.2e-04 1.8e-02 2.6e-03 | 2.4e-06 1.2e-03 1.1e-05 + 1800 3.0e-02 | -1.2e+02 -1.2e+02 | 2.8e-04 1.2e-02 6.5e-03 | 2.1e-06 8.0e-04 2.7e-05 + 2000 3.3e-02 | -1.2e+02 -1.2e+02 | 2.1e-04 1.3e-02 3.5e-03 | 1.6e-06 8.3e-04 1.4e-05 + 2200 3.6e-02 | -1.2e+02 -1.2e+02 | 3.0e-03 1.1e-02 5.8e-03 | 2.3e-05 7.2e-04 2.4e-05 + 2400 4.0e-02 | -1.2e+02 -1.2e+02 | 3.1e-03 5.9e-03 1.9e-03 | 2.4e-05 3.9e-04 7.8e-06 + 2600 4.3e-02 | -1.2e+02 -1.2e+02 | 2.2e-03 3.1e-03 1.4e-03 | 1.7e-05 2.1e-04 5.7e-06 + 2800 4.6e-02 | -1.2e+02 -1.2e+02 | 1.0e-03 2.1e-03 1.7e-03 | 7.6e-06 1.4e-04 6.8e-06 + 3000 4.8e-02 | -1.2e+02 -1.2e+02 | 6.4e-04 3.6e-04 4.6e-04 | 4.9e-06 2.4e-05 1.9e-06 + --------------------------------------------------------------------------------------- + Solution Summary + Status : OPTIMAL + Presolve time : 0.00853 sec + Precondition time : 0.003735 sec + Solve time : 0.0517 sec + Iterations : 3000 + Primal objective : -121.2216698 + Dual objective : -121.2221271 + Objective gap : 1.879e-06 + Primal infeas : 4.889e-06 + Dual infeas : 2.399e-05 + ``` + +## Problem and settings + +
+ +```text +--------------------------------------------------------------------------------------- + cuPDLPx v0.3.0 + A GPU-Accelerated First-Order LP Solver + (c) Haihao Lu, Massachusetts Institute of Technology, 2025 +--------------------------------------------------------------------------------------- +Problem: 17013 rows, 200 columns, 104811 nonzeros +Settings: + iter_limit : 2147483647 + time_limit : 3600.00 sec + eps_opt : 1.0e-04 + eps_feas : 1.0e-04 + spmv_backend : cusparseSpMVOp (auto) +``` + +
+ +The problem dimensions, four core settings, and selected sparse matrix–vector +backend are always shown. Additional settings appear when they differ from +their defaults. See [Parameters](../reference/parameters.md) for the +corresponding interface names. + +## Presolve and preconditioning + +
+ +```text +Running presolver (PSLP v0.0.8)... + status : REDUCED + presolve time : 0.00853 sec + reduced problem : 16997 rows, 200 columns, 104747 nonzeros + +Preconditioning + Ruiz scaling (10 iterations) + Pock-Chambolle scaling (alpha=1.0000) + Bound-objective scaling +``` + +
+ +`REDUCED` means that presolve produced an equivalent smaller LP. Postsolve +recovers the solution vectors for the original LP. The reported residuals +and gap are computed on the presolved LP. + +Preconditioning rescales the LP before the main iteration. Its work is not +included in the PDHG iteration count. + +## Progress table + +For this example, the iteration log is: + +
+ +```text +--------------------------------------------------------------------------------------- + runtime | objective | absolute residuals | relative residuals + iter time | pr obj du obj | pr res du res gap | pr res du res gap +--------------------------------------------------------------------------------------- + 0 0.0e+00 | 0.0e+00 0.0e+00 | 0.0e+00 0.0e+00 0.0e+00 | 0.0e+00 0.0e+00 0.0e+00 + 200 0.0e+00 | -1.3e+02 -1.2e+02 | 1.1e+00 2.5e-01 9.7e+00 | 8.5e-03 1.7e-02 3.9e-02 + 400 8.7e-03 | -1.2e+02 -1.2e+02 | 2.3e-01 7.6e-02 9.2e-01 | 1.7e-03 5.0e-03 3.8e-03 + 600 1.2e-02 | -1.2e+02 -1.2e+02 | 4.6e-02 8.5e-02 5.2e-03 | 3.5e-04 5.6e-03 2.1e-05 + 800 1.5e-02 | -1.2e+02 -1.2e+02 | 2.5e-02 4.7e-02 3.3e-02 | 1.9e-04 3.1e-03 1.4e-04 + 1000 1.8e-02 | -1.2e+02 -1.2e+02 | 1.2e-02 3.4e-02 7.3e-03 | 8.8e-05 2.2e-03 3.0e-05 + 1200 2.1e-02 | -1.2e+02 -1.2e+02 | 4.3e-03 2.5e-02 7.6e-04 | 3.2e-05 1.6e-03 3.1e-06 + 1400 2.4e-02 | -1.2e+02 -1.2e+02 | 3.1e-03 1.9e-02 1.4e-03 | 2.4e-05 1.3e-03 5.9e-06 + 1600 2.7e-02 | -1.2e+02 -1.2e+02 | 3.2e-04 1.8e-02 2.6e-03 | 2.4e-06 1.2e-03 1.1e-05 + 1800 3.0e-02 | -1.2e+02 -1.2e+02 | 2.8e-04 1.2e-02 6.5e-03 | 2.1e-06 8.0e-04 2.7e-05 + 2000 3.3e-02 | -1.2e+02 -1.2e+02 | 2.1e-04 1.3e-02 3.5e-03 | 1.6e-06 8.3e-04 1.4e-05 + 2200 3.6e-02 | -1.2e+02 -1.2e+02 | 3.0e-03 1.1e-02 5.8e-03 | 2.3e-05 7.2e-04 2.4e-05 + 2400 4.0e-02 | -1.2e+02 -1.2e+02 | 3.1e-03 5.9e-03 1.9e-03 | 2.4e-05 3.9e-04 7.8e-06 + 2600 4.3e-02 | -1.2e+02 -1.2e+02 | 2.2e-03 3.1e-03 1.4e-03 | 1.7e-05 2.1e-04 5.7e-06 + 2800 4.6e-02 | -1.2e+02 -1.2e+02 | 1.0e-03 2.1e-03 1.7e-03 | 7.6e-06 1.4e-04 6.8e-06 + 3000 4.8e-02 | -1.2e+02 -1.2e+02 | 6.4e-04 3.6e-04 4.6e-04 | 4.9e-06 2.4e-05 1.9e-06 +--------------------------------------------------------------------------------------- +``` + +
+ +| Column | Meaning | +| --- | --- | +| `iter` | Total PDHG iterations completed. | +| `time` | Elapsed main-solve time in seconds. | +| `pr obj` | Primal objective value. | +| `du obj` | Dual objective value. | +| `pr res` | Primal feasibility residual. | +| `du res` | Dual feasibility residual. | +| `gap` | Primal–dual objective gap. | + +The log reports absolute and relative primal residuals, dual residuals, and +gaps. The optimality test requires the relative residuals to be below +`FeasibilityTol` and the relative gap to be below `OptimalityTol`. +See [Termination criteria](../algorithm/termination.md#optimality-conditions) +for the formulas. + +`TermCheckFreq` sets the termination-check interval, `200` in this example. +Progress is logged at these checks, with fewer rows printed as the iteration +count grows. Restarts are checked at the same interval and have no separate +log marker. Residuals and objective values need not decrease monotonically. + +## Feasibility polishing + +When feasibility polishing runs, the log shows separate primal and dual +tables with objective values and absolute and relative residuals. The +polishing summary reports each phase's status, iteration count, and time. + +Polishing can reduce feasibility residuals while increasing the primal–dual +gap. The main solve's termination status is unchanged. See +[Feasibility polishing](../algorithm/feasibility-polishing.md). + +## Solution summary + +Read the final summary in this order: + +1. **Status:** check why the solver stopped. `OPTIMAL` means it reached the + requested accuracy; a time or iteration limit does not establish that. +2. **Primal and dual residuals:** compare `Primal infeas` and `Dual infeas` + with `FeasibilityTol`. +3. **Objective gap:** compare `Objective gap` with `OptimalityTol`. + +For example: + +
+ +```text +Solution Summary + Status : OPTIMAL + Presolve time : 0.00853 sec + Precondition time : 0.003735 sec + Solve time : 0.0517 sec + Iterations : 3000 + Primal objective : -121.2216698 + Dual objective : -121.2221271 + Objective gap : 1.879e-06 + Primal infeas : 4.889e-06 + Dual infeas : 2.399e-05 +``` + +
+ +Both residuals and the relative gap are below the requested `1e-4` tolerances. + +The fields have the following meanings: + +| Field | Meaning | +| --- | --- | +| `Status` | Termination reason, such as `OPTIMAL`, `TIME_LIMIT`, or `ITERATION_LIMIT`. | +| `Presolve time` | Time spent in PSLP, shown when presolve is enabled. | +| `Precondition time` | Time spent scaling the LP and preparing the scaled problem. | +| `Solve time` | Time spent in the main iteration. | +| `Iterations` | Total main-solve iteration count. | +| `Primal objective`, `Dual objective` | Final objective values, including the objective constant and original objective sense. | +| `Objective gap` | Final relative primal–dual gap on the presolved LP when presolve is enabled. | +| `Primal infeas`, `Dual infeas` | Final relative residuals on the presolved LP when presolve is enabled. | + +`Primal infeas` and `Dual infeas` are residual magnitudes. Use `Status` to +determine whether the solver detected infeasibility. See +[Results and status](results.md) for termination reasons and result fields. diff --git a/docs/getting-started/results.md b/docs/getting-started/results.md new file mode 100644 index 0000000..b93e8ef --- /dev/null +++ b/docs/getting-started/results.md @@ -0,0 +1,98 @@ +--- +description: Solution fields, termination statuses, and infeasibility information returned by cuPDLPx. +--- + +# Results and status + +cuPDLPx reports a termination status, primal and dual vectors, residuals, +objective values, and solve statistics. Check the status before using the +returned vectors. Field names differ across interfaces. + +## Termination status + +
+ +| Status | Meaning | Next step | +| --- | --- | --- | +| Optimal | Presolve solved the LP, or the main iteration met the feasibility and optimality tolerances. | Check the final residuals and gap, especially if polishing is enabled. | +| Primal infeasible | The solver detected primal infeasibility. | Review the constraints and bounds, along with any returned infeasibility information. | +| Dual infeasible | The solver detected dual infeasibility, which can indicate an unbounded primal problem. | Review any returned infeasibility information and check for unbounded directions. | +| Infeasible or unbounded | Presolve or the solver could not distinguish infeasibility from unboundedness. | Inspect the model and log for more detail. | +| Time limit | The wall-clock limit was reached. | Check residuals and the gap; increase the time limit if needed. | +| Iteration limit | The iteration limit was reached. | Check residuals and the gap; increase the iteration limit if needed. | +| Feasibility polishing succeeded | A feasibility-polishing phase reached the polishing tolerance. | Check the final residuals and gap; polishing alone does not establish optimality. | +| Unspecified | No more specific termination reason is available. | Inspect the log and input data. | + +
+ +[Feasibility polishing](../algorithm/feasibility-polishing.md) preserves the +main solve's termination status. The final gap can exceed the optimality +tolerance even when the status is `OPTIMAL`. + +!!! warning "Always inspect the status" + + Reaching a time or iteration limit does not establish feasibility or + optimality. Check the residuals and primal–dual gap before using the + returned iterate. + +## Solution and quality measures + +The result contains: + +| Result | Description | +| --- | --- | +| Primal solution $x$ | One value for each variable. | +| Dual solution $y$ | One multiplier for each constraint. | +| Dual slacks $r$ | One dual slack for each variable. | +| Primal objective | $c^\top x+c_0$ in the original objective sense. | +| Dual objective | Reported dual objective value; a valid bound requires dual feasibility. | +| Objective gap | Absolute and relative primal–dual gaps. | +| Primal residual | Absolute and relative violation of the primal constraints. | +| Dual residual | Absolute and relative violation of $c-A^\top y-r=0$. | +| Infeasibility information | Primal- and dual-ray quality measures when applicable. | +| Work statistics | Iteration counts and phase timings. | + +With presolve enabled, the reported residuals and objective gap refer to +the [presolved model](../algorithm/presolve.md#presolve), while the solution +vectors are recovered for the original model. + +See [Python results](../guides/python.md#solve-and-inspect-results), the [Julia +interface](../guides/julia.md), [C result fields](../guides/c-api.md#result-fields), +and [command-line output files](../guides/command-line.md#output-files) +for field names and access methods. + +## Status constants + +Use symbolic constants to compare statuses; Python and C use different +integer values. Julia maps the native termination reason to a +MathOptInterface status; see the [Julia interface](../guides/julia.md). + +
+ +=== "Python" + + | Status | Constant | + | --- | --- | + | Optimal | `PDLP.OPTIMAL` | + | Primal infeasible | `PDLP.PRIMAL_INFEASIBLE` | + | Dual infeasible | `PDLP.DUAL_INFEASIBLE` | + | Infeasible or unbounded | `PDLP.INFEASIBLE_OR_UNBOUNDED` | + | Time limit | `PDLP.TIME_LIMIT` | + | Iteration limit | `PDLP.ITERATION_LIMIT` | + | Feasibility polishing succeeded | `PDLP.FEAS_POLISH_SUCCESS` | + | Unspecified | `PDLP.UNSPECIFIED` | + +=== "C" + + | Status | Constant | + | --- | --- | + | Optimal | `TERMINATION_REASON_OPTIMAL` | + | Primal infeasible | `TERMINATION_REASON_PRIMAL_INFEASIBLE` | + | Dual infeasible | `TERMINATION_REASON_DUAL_INFEASIBLE` | + | Infeasible or unbounded | `TERMINATION_REASON_INFEASIBLE_OR_UNBOUNDED` | + | Time limit | `TERMINATION_REASON_TIME_LIMIT` | + | Iteration limit | `TERMINATION_REASON_ITERATION_LIMIT` | + | Feasibility polishing succeeded | `TERMINATION_REASON_FEAS_POLISH_SUCCESS` | + | Unspecified | `TERMINATION_REASON_UNSPECIFIED` | + +
diff --git a/docs/guides/c-api.md b/docs/guides/c-api.md new file mode 100644 index 0000000..583b89a --- /dev/null +++ b/docs/guides/c-api.md @@ -0,0 +1,188 @@ +--- +description: Native C data structures, solver lifecycle, parameters, and result handling for cuPDLPx. +--- + +# C interface + +The C API solves LPs from arrays in host memory. Functions are declared in +[`include/cupdlpx.h`](https://github.com/MIT-Lu-Lab/cuPDLPx/blob/main/include/cupdlpx.h), +with public structures and enums in `include/cupdlpx_types.h`. + +## Installation + +The C library and command-line executable are built from the same source tree. +See [hardware requirements](../getting-started/index.md#hardware-requirements) +for supported GPUs and required CUDA or ROCm versions. + +Follow the [native installation instructions](command-line.md#installation) +to build the library. + +## Create a problem + +```c +lp_problem_t *create_lp_problem( + const double *objective_c, + const matrix_desc_t *A_desc, + const double *con_lb, + const double *con_ub, + const double *var_lb, + const double *var_ub, + const double *objective_constant, + const objective_sense_t *objective_sense +); +``` + +`create_lp_problem` copies the input arrays into a new `lp_problem_t`. +The caller retains ownership of the inputs. Free the problem with +`lp_problem_free` and the result of `solve_lp_problem` with +`cupdlpx_result_free`. + +Only `A_desc` is required. Passing `NULL` for another argument selects its +default: + +| Argument | Length | `NULL` default | +| --- | ---: | --- | +| `objective_c` | $n$ | all zeros | +| `con_lb` | $m$ | all $-\infty$ | +| `con_ub` | $m$ | all $+\infty$ | +| `var_lb` | $n$ | all $-\infty$ | +| `var_ub` | $n$ | all $+\infty$ | +| `objective_constant` | 1 | `0.0` | +| `objective_sense` | 1 | `OBJECTIVE_SENSE_MINIMIZE` | + +## Matrix descriptors + +`matrix_desc_t` accepts four host-memory layouts: + +| Format | Enum | Required arrays | +| --- | --- | --- | +| Row-major dense | `matrix_dense` | `A` with $m n$ values | +| CSR | `matrix_csr` | `row_ptr`, `col_ind`, `vals`, `nnz` | +| CSC | `matrix_csc` | `col_ptr`, `row_ind`, `vals`, `nnz` | +| COO | `matrix_coo` | `row_ind`, `col_ind`, `vals`, `nnz` | + +Indices are zero-based `int` values and numeric data uses `double`. + +See the +[cuSPARSE matrix formats documentation](https://docs.nvidia.com/cuda/cusparse/index.html#matrix-formats). + +## Solve a small LP + +```c +#include "cupdlpx.h" +#include +#include + +int main(void) { + double A[3][2] = { + {1.0, 2.0}, + {0.0, 1.0}, + {3.0, 2.0} + }; + double c[2] = {1.0, 1.0}; + double l[3] = {5.0, -INFINITY, -INFINITY}; + double u[3] = {5.0, 2.0, 8.0}; + + matrix_desc_t matrix = { + .m = 3, + .n = 2, + .fmt = matrix_dense, + .data.dense = {.A = &A[0][0]}, + }; + + lp_problem_t *problem = create_lp_problem( + c, &matrix, l, u, NULL, NULL, NULL, NULL); + if (problem == NULL) { + return 1; + } + + pdhg_parameters_t params; + set_default_parameters(¶ms); + params.verbose = false; + params.termination_criteria.eps_optimal_relative = 1e-6; + params.termination_criteria.eps_feasible_relative = 1e-6; + + cupdlpx_result_t *result = solve_lp_problem(problem, ¶ms); + if (result == NULL) { + lp_problem_free(problem); + return 1; + } + + printf("termination reason: %d\n", result->termination_reason); + printf("objective: %.6f\n", result->primal_objective_value); + for (int j = 0; j < result->num_variables; ++j) { + printf("x[%d] = %.6f\n", j, result->primal_solution[j]); + } + + cupdlpx_result_free(result); + lp_problem_free(problem); + return 0; +} +``` + +## Warm starts + +To warm-start the solver, disable presolve and provide a primal vector, +a dual vector, or both: + +```c +double x0[2] = {1.0, 2.0}; +double y0[3] = {1.0, -1.0, 0.0}; + +params.presolve = false; +set_start_values(problem, x0, y0); +``` + +The function copies the supplied arrays. Passing `NULL` clears the +corresponding starting vector. + +## Default parameters + +Call `set_default_parameters` before changing fields: + +```c +pdhg_parameters_t params; +set_default_parameters(¶ms); + +params.termination_criteria.time_sec_limit = 300.0; +params.feasibility_polishing = true; +``` + +Passing `NULL` as the second argument of `solve_lp_problem` also selects all +defaults: + +```c +cupdlpx_result_t *result = solve_lp_problem(problem, NULL); +``` + +See the [parameter reference](../reference/parameters.md#usage) for +the main fields and defaults. + +## Result fields + +`solve_lp_problem` returns a `cupdlpx_result_t *`. Its main fields are: + +| Result | `cupdlpx_result_t` field | +| --- | --- | +| Termination status | `termination_reason` | +| Original dimensions | `num_variables`, `num_constraints`, `num_nonzeros` | +| Reduced dimensions | `num_reduced_variables`, `num_reduced_constraints`, `num_reduced_nonzeros` | +| Primal solution, dual solution, dual slacks | `primal_solution`, `dual_solution`, `reduced_cost` | +| Objectives | `primal_objective_value`, `dual_objective_value` | +| Gaps | `objective_gap`, `relative_objective_gap` | +| Primal residuals | `absolute_primal_residual`, `relative_primal_residual` | +| Dual residuals | `absolute_dual_residual`, `relative_dual_residual` | +| Iterations | `total_count`, `feasibility_iteration` | +| Timings | `cumulative_time_sec`, `rescaling_time_sec`, `presolve_time`, `feasibility_polishing_time` | +| Ray measures | `max_primal_ray_infeasibility`, `max_dual_ray_infeasibility`, `primal_ray_linear_objective`, `dual_ray_objective` | + +The arrays remain valid until `cupdlpx_result_free(result)` is called. +Compare `termination_reason` with symbolic members of `termination_reason_t` +rather than their integer values. See [results and +status](../getting-started/results.md) for interpretation. + +## Error handling + +Creation and solve functions return `NULL` on failure. Check each returned +pointer before dereferencing it and release any objects already created on the +error path. diff --git a/docs/guides/command-line.md b/docs/guides/command-line.md new file mode 100644 index 0000000..3c49e27 --- /dev/null +++ b/docs/guides/command-line.md @@ -0,0 +1,151 @@ +--- +description: Build and use the cuPDLPx command-line interface with CUDA or ROCm. +--- + +# Command-line interface + +The `cupdlpx` executable solves LPs in `.mps` or `.mps.gz` format and writes +solution vectors and a solve summary to the output directory. +Build it from source with the C libraries. + +## Build requirements + +| Component | CUDA build | ROCm build | +| --- | --- | --- | +| Compiler | GCC and NVCC | GCC and `hipcc` | +| Build system | CMake 3.20+ | CMake 3.20+ | + +See [hardware requirements](../getting-started/index.md#hardware-requirements) +for supported GPUs and required CUDA or ROCm versions. + +## Installation + +Clone the repository: + +```bash +git clone https://github.com/MIT-Lu-Lab/cuPDLPx.git +cd cuPDLPx +``` + +### CUDA + +Configure CMake and build: + +```bash +cmake -B build -DCMAKE_BUILD_TYPE=Release +cmake --build build --clean-first --parallel +``` + +The build selects the sparse matrix–vector backend from the cuSPARSE +version: + +| CUDA toolkit | Backend | +| --- | --- | +| CUDA 12.4–13.2 | `cusparseSpMV` | +| CUDA 13.3+ | `cusparseSpMVOp` | + +To target a specific CUDA architecture, pass +`-DCMAKE_CUDA_ARCHITECTURES=` when configuring. + +### ROCm + +Enable HIP and provide the architecture of the target AMD GPU: + +```bash +cmake -B build \ + -DUSE_HIP=ON \ + -DCMAKE_BUILD_TYPE=Release \ + -DCMAKE_HIP_ARCHITECTURES=gfx90a \ + -DCMAKE_PREFIX_PATH=/opt/rocm +cmake --build build --clean-first --parallel +``` + +Common architecture values include `gfx90a` for MI200, `gfx1100` for RDNA3, +and `gfx1201` for RDNA4. Replace the value with the architecture supported by +the installed ROCm toolchain. + +On HIP builds, cuBLAS, cuSPARSE, and CUB calls are mapped to hipBLAS, +hipSPARSE, and hipCUB. + +### Build outputs + +Both configurations create: + +- `build/cupdlpx`, the command-line solver; +- a static core library; and +- a shared `cupdlpx` library. + +## Usage + +```text +cupdlpx [OPTIONS] +``` + +For a source build, the executable is normally `./build/cupdlpx`: + +```bash +mkdir -p results +./build/cupdlpx problem.mps.gz results +``` + +The output directory must already exist. + +## Solver options + +Pass solver options before the input and output paths. See +[Parameters](../reference/parameters.md) for defaults and allowed values, +or run `./build/cupdlpx --help`. + +```bash +./build/cupdlpx \ + --time_limit 600 \ + --eps_opt 1e-6 \ + --eps_feas 1e-6 \ + --opt_norm linf \ + problem.mps.gz results +``` + +## Output files + +For `problem.mps.gz`, cuPDLPx creates: + +```text +results/ +├── problem_summary.txt +├── problem_primal_solution.txt +└── problem_dual_solution.txt +``` + +The primal and dual files contain one value per line. The summary records the +termination reason, model dimensions, objective values, residuals, iteration +count, and phase timings. Dual slacks are available through the in-memory +Python, Julia, and C interfaces rather than as a separate CLI output file. + +!!! tip "In-memory results" + + Use the [Python interface](python.md) or [C interface](c-api.md) when a + program needs results in memory instead of text files. + +## Advanced build options + +The native build exposes the following CMake options: + +| Option | Default | Purpose | +| --- | --- | --- | +| `CUPDLPX_BUILD_STATIC_LIB` | `ON` | Build the static core library. | +| `CUPDLPX_BUILD_SHARED_LIB` | `ON` | Build the shared library. | +| `CUPDLPX_BUILD_CLI` | `ON` | Build the `cupdlpx` executable. | +| `CUPDLPX_BUILD_PYTHON` | `OFF` | Build the Python extension. | +| `CUPDLPX_BUILD_TESTS` | `OFF` | Build the native test suite. | + +For example, build only the shared library: + +```bash +cmake -B build \ + -DCUPDLPX_BUILD_STATIC_LIB=OFF \ + -DCUPDLPX_BUILD_CLI=OFF +cmake --build build --parallel +``` + +See [GPU backends](../implementation/gpu-backends.md) for architecture +defaults and platform requirements. diff --git a/docs/guides/julia.md b/docs/guides/julia.md new file mode 100644 index 0000000..a24aa8a --- /dev/null +++ b/docs/guides/julia.md @@ -0,0 +1,30 @@ +--- +description: Requirements and installation for using cuPDLPx from Julia through CuPDLPx.jl. +--- + +# Julia interface + +CuPDLPx.jl provides a Julia wrapper for cuPDLPx and an interface to JuMP. + +## Software requirements + +The Julia interface requires Julia 1.10+. + +See [hardware requirements](../getting-started/index.md#hardware-requirements) +for supported GPUs and required CUDA or ROCm versions. + +## Installation + +Install CuPDLPx.jl from the Julia package manager: + +```julia +import Pkg +Pkg.add("CuPDLPx") +``` + +The package installs the cuPDLPx binaries automatically; a separate source +build is not normally required. + +CuPDLPx.jl uses the standard JuMP and MathOptInterface APIs for model +construction, parameters, and results. See the [CuPDLPx.jl +documentation](https://jump.dev/JuMP.jl/stable/packages/CuPDLPx/){ target="_blank" rel="noopener noreferrer" }. diff --git a/docs/guides/python.md b/docs/guides/python.md new file mode 100644 index 0000000..8f20bd0 --- /dev/null +++ b/docs/guides/python.md @@ -0,0 +1,328 @@ +--- +description: Install and use the cuPDLPx Python interface, including its model API, parameters, warm starts, and results. +--- + +# Python interface + +The `cupdlpx` package solves LPs from [NumPy](https://numpy.org/doc/stable/) +arrays or [SciPy sparse](https://docs.scipy.org/doc/scipy/reference/sparse.html) +matrices. The `Model` class stores problem data, solver settings, and results. +Solution vectors are returned as read-only NumPy arrays. + +To solve an MPS file, see [Read MPS data](#read-mps-data). +To build an LP from arrays, see the [Quickstart](#quickstart). + +## Software requirements + +| Component | Requirement | +| --- | --- | +| Python | 3.8+ | +| NumPy | 1.21+ | +| SciPy | 1.8+ | + +NumPy and SciPy are installed automatically with the package. + +See [hardware requirements](../getting-started/index.md#hardware-requirements) +for supported GPUs and required CUDA or ROCm versions. + +## Installation + +Install from PyPI: + +```bash +python -m pip install cupdlpx +``` + +Confirm that the package imports: + +```bash +python -c "import cupdlpx; print(cupdlpx.__version__)" +``` + +To install the current source tree instead: + +```bash +git clone https://github.com/MIT-Lu-Lab/cuPDLPx.git +cd cuPDLPx +python -m pip install . +``` + +Building from source requires CMake 3.21+ and a +[CUDA or ROCm toolchain](command-line.md#build-requirements). Pip may build +the package locally if no compatible wheel is available. + +## Read MPS data + +To solve an existing LP, replace `instance.mps` with the path to your file: + +```python +import cupdlpx + +model = cupdlpx.read("instance.mps") +model.optimize() + +if model.Status == cupdlpx.PDLP.OPTIMAL: + print("objective:", model.ObjVal) + print("x:", model.X) +else: + print("solver stopped with", model.StatusName) +``` + +`cupdlpx.read(filename)` accepts any path-like object and returns a new `Model`; +both `.mps` and `.mps.gz` files are supported. A missing path raises +`FileNotFoundError`; a file rejected by the native parser raises `RuntimeError`. + +See [Solve and inspect results](#solve-and-inspect-results) for residuals, +timings, and other solution fields. + +## Quickstart + +This example solves + +$$ +\begin{aligned} +\operatorname*{minimize}_{x_1,x_2}\quad & x_1+x_2 \\ +\text{subject to}\quad +& x_1 + 2x_2 = 5, \\ +& x_2 \le 2, \\ +& 3x_1 + 2x_2 \le 8, \\ +& x_1,x_2 \ge 0. +\end{aligned} +$$ + +Every constraint is represented by a lower and an upper bound: + +```python +import numpy as np +from cupdlpx import Model, PDLP + +c = np.array([1.0, 1.0]) +A = np.array([ + [1.0, 2.0], + [0.0, 1.0], + [3.0, 2.0], +]) + +model = Model( + objective_vector=c, + constraint_matrix=A, + constraint_lower_bound=np.array([5.0, -np.inf, -np.inf]), + constraint_upper_bound=np.array([5.0, 2.0, 8.0]), + variable_lower_bound=np.zeros(2), + variable_upper_bound=None, +) + +model.Params.TimeLimit = 60 +model.setParams( + OptimalityTol=1e-6, + FeasibilityTol=1e-6, + OutputFlag=False, +) +model.optimize() + +if model.Status == PDLP.OPTIMAL: + print(f"objective: {model.ObjVal:.4f}") + print("x:", np.round(model.X, 4)) +else: + print("solver stopped with", model.StatusName) +``` + +Output: + +```text +objective: 3.0000 +x: [1. 2.] +``` + +Passing `None` omits the corresponding lower or upper bounds. +Set `OutputFlag=True` to print the solver progress table; its fields are explained under [Log +interpretation](../getting-started/log-interpretation.md). + +## Create a model + +```python +from cupdlpx import Model + +model = Model( + objective_vector=c, + constraint_matrix=A, + constraint_lower_bound=l, + constraint_upper_bound=u, + variable_lower_bound=lb, + variable_upper_bound=ub, + objective_constant=0.0, +) +``` + +`constraint_matrix` accepts a two-dimensional `numpy.ndarray` or a SciPy +sparse matrix. + +
+ +| Argument | Type | Shape | Default | Description | +| --- | --- | --- | --- | --- | +| `objective_vector` | array-like | `(n,)` | required | Finite objective coefficients $c$. | +| `constraint_matrix` | NumPy array or SciPy sparse matrix | `(m,n)` | required | Finite constraint matrix $A$. | +| `constraint_lower_bound` | array-like or `None` | `(m,)` | required | Lower bounds $\ell_c$; `None` means all $-\infty$. | +| `constraint_upper_bound` | array-like or `None` | `(m,)` | required | Upper bounds $u_c$; `None` means all $+\infty$. | +| `variable_lower_bound` | array-like or `None` | `(n,)` | `None` | Lower bounds $\ell_v$; `None` means all $-\infty$. | +| `variable_upper_bound` | array-like or `None` | `(n,)` | `None` | Upper bounds $u_v$; `None` means all $+\infty$. | +| `objective_constant` | real number | scalar | `0.0` | Finite objective offset $c_0$. | + +
+ +Dense vectors and matrices are copied into C-contiguous `float64` arrays. +Sparse matrices are copied to CSR with `float64` values and 32-bit indices. +The stored arrays are read-only. + +Construction raises `TypeError` for unsupported input types, `ValueError` for +invalid dimensions, bounds, or nonfinite data, and `OverflowError` when a +sparse matrix cannot use 32-bit CSR indices. + +Models minimize by default. To maximize: + +```python +from cupdlpx import PDLP + +model.ModelSense = PDLP.MAXIMIZE +``` + +## Solve and inspect results + +Call `optimize()` to solve the model, then read its status and solution attributes: + +```python +model.optimize() +status = model.StatusName +``` + +| Property | Type after a solve | Description | +| --- | --- | --- | +| `Status` | `int` | Termination code; compare with a constant in `cupdlpx.PDLP`. | +| `StatusName` | `str` | Symbolic termination name. | +| `X` | `numpy.ndarray`
shape `(n,)` | Primal solution. | +| `Pi` | `numpy.ndarray`
shape `(m,)` | Dual solution. | +| `RC` | `numpy.ndarray`
shape `(n,)` | Dual slacks. | +| `ObjVal` | `float` | Primal objective value. | +| `DualObj` | `float` | Dual objective value. | +| `Gap` | `float` | Absolute primal–dual objective gap. | +| `RelGap` | `float` | Relative primal–dual objective gap. | +| `RelPrimalResidual` | `float` | Relative primal feasibility residual. | +| `RelDualResidual` | `float` | Relative dual feasibility residual. | +| `IterCount` | `int` | Total iteration count. | +| `Runtime` | `float` | Solve time in seconds. | +| `RescalingTime` | `float` | Preconditioning time in seconds. | +| `MaxPrimalRayInfeas` | `float` | Maximum primal-ray infeasibility. | +| `MaxDualRayInfeas` | `float` | Maximum dual-ray infeasibility. | +| `PrimalRayLinObj` | `float` | Primal-ray linear objective. | +| `DualRayObj` | `float` | Dual-ray objective. | + +Solution attributes are `None` before a solve and are cleared whenever model +data changes. `PrimalInfeas` and `DualInfeas` are aliases for +`RelPrimalResidual` and `RelDualResidual`. Invalid model bounds or objective +sense raise `ValueError`; native solver failures raise `RuntimeError`. See +[results and status](../getting-started/results.md) for status interpretation. + +## Set parameters + +```python +model.setParam(name, value) +model.getParam(name) +model.setParams(**kwargs) +model.resetParams() +``` + +Each of these forms sets the same time limit: + +```python +model.Params.TimeLimit = 120 +model.setParam("TimeLimit", 120) +model.setParams(TimeLimit=120) +``` + +To set several parameters at once, use +`model.setParams(TimeLimit=120, OptimalityTol=1e-6)`. + +C parameter names are also accepted: + +```python +model.setParam("time_sec_limit", 120) +``` + +Read or reset parameters with: + +```python +print(model.getParam("TimeLimit")) +print(dict(model.Params.items())) +model.resetParams() +``` + +`setParams` validates every update before applying any of them. `model.Params` +is a live parameter view that supports attribute and item access, iteration, +and `keys()`, `values()`, and `items()`. + +An unknown name raises `KeyError`, a value of the wrong type raises +`TypeError`, and a value outside its allowed range raises `ValueError`. See the +[parameter reference](../reference/parameters.md). + +## Update model data + +Update model data by assigning properties or calling setter methods. +Arrays returned by model properties are read-only; assign a new array +instead of modifying individual entries in place: + +```python +model.c = new_objective +model.A = new_constraint_matrix +model.lb = new_variable_lower_bounds +model.ub = new_variable_upper_bounds + +model.setObjectiveConstant(2.0) +model.setConstraintLowerBound(new_l) +model.setConstraintUpperBound(new_u) +``` + +| Property | Type | Validating setter | +| --- | --- | --- | +| `c` | `numpy.ndarray`
shape `(n,)` | `setObjectiveVector(c)` | +| `c0` | `float` | `setObjectiveConstant(c0)` | +| `A` | NumPy array or CSR matrix
shape `(m,n)` | `setConstraintMatrix(A)` | +| `constr_lb` | `numpy.ndarray` or `None`
shape `(m,)` | `setConstraintLowerBound(lower)` | +| `constr_ub` | `numpy.ndarray` or `None`
shape `(m,)` | `setConstraintUpperBound(upper)` | +| `lb` | `numpy.ndarray` or `None`
shape `(n,)` | `setVariableLowerBound(lower)` | +| `ub` | `numpy.ndarray` or `None`
shape `(n,)` | `setVariableUpperBound(upper)` | +| `ModelSense` | `int` | Assign `PDLP.MINIMIZE` or `PDLP.MAXIMIZE`. | +| `num_vars` | `int` | — | +| `num_constrs` | `int` | — | + +Changing model data clears the cached solution. The number of variables is +fixed at construction. To change the number of constraints, clear the +constraint bounds, replace `A`, and set new bounds with matching lengths. + +## Warm start + +```python +model.setWarmStart(primal=..., dual=...) +model.clearWarmStart() +``` + +To warm-start the solver, disable presolve and provide a primal vector, +a dual vector, or both: + +```python +model.setParam("Presolve", False) +model.setWarmStart(primal=x0, dual=y0) +model.optimize() +``` + +Omitting an argument preserves its previous starting vector. Passing `None` +clears it: + +```python +model.setWarmStart(primal=None) +model.clearWarmStart() +``` + +The primal vector must have length `model.num_vars`; the dual vector must have +length `model.num_constrs`. A length mismatch or nonfinite value raises +`ValueError`. diff --git a/docs/implementation/cuda-graphs.md b/docs/implementation/cuda-graphs.md new file mode 100644 index 0000000..58400a7 --- /dev/null +++ b/docs/implementation/cuda-graphs.md @@ -0,0 +1,55 @@ +--- +description: How cuPDLPx uses CUDA Graphs to reduce launch overhead and reuse GPU work across restarts. +--- + +# CUDA Graphs + +cuPDLPx uses [CUDA Graphs](https://docs.nvidia.com/cuda/cuda-programming-guide/04-special-topics/cuda-graphs.html) +to reduce kernel launch overhead. A graph records GPU operations and their +dependencies, allowing the sequence to be replayed with a single CPU launch. +[Kernel fusion](kernel-fusion.md) also reduces memory traffic within each +iteration. + +
+ ![Individual kernel launches compared with a single CUDA Graph launch](../assets/cuda-graphs-launch-overhead.svg#only-light) + ![Individual kernel launches compared with a single CUDA Graph launch](../assets/cuda-graphs-launch-overhead-dark.svg#only-dark) +
+ Launch overhead can leave gaps between short GPU kernels; + a graph launch submits the recorded sequence together. Adapted from + NVIDIA, Effortless CUDA Graphs, GTC 2021. +
+
+ +## Graph capture and replay + +cuPDLPx evaluates termination and restart criteria every +[`termination_evaluation_frequency`](../reference/parameters.md#limits-and-logging) iterations +(`200` by default). The iterations between consecutive checks form a +window, executed as follows: + +1. Execute the first iteration of each window outside the graph. After a + restart, its PDHG iterate gives the initial fixed-point error of the new + epoch, the baseline for the [restart criteria](../algorithm/restart.md) + and the [active-set step size boost](../algorithm/step-size.md#divergence-protection). +2. **Launch a graph** covering iterations 2 through + `termination_evaluation_frequency` within the window. Before its first + launch, the graph is captured and instantiated; subsequent windows reuse + it. The final iteration also saves the intermediate values needed for + the checks. +3. Compute residual vectors and perform reductions on the GPU. The CPU uses + the resulting scalar measures to evaluate termination and restart criteria. + +## Graph reuse across restarts + +The iterate arrays and work buffers keep the same device addresses across +restarts. cuPDLPx updates the anchor, iteration counter, and step sizes in +GPU memory, so the same graph can be reused without recapture. + +!!! note "Check frequency" + + Increasing `termination_evaluation_frequency` reduces graph launches, + reductions, and host synchronization, but delays convergence checks, + restarts, and detection of time or iteration limits. + +ROCm builds use the corresponding HIP Graph APIs through the +[GPU backend compatibility layer](gpu-backends.md). diff --git a/docs/implementation/gpu-backends.md b/docs/implementation/gpu-backends.md new file mode 100644 index 0000000..72d457f --- /dev/null +++ b/docs/implementation/gpu-backends.md @@ -0,0 +1,74 @@ +--- +description: CUDA and ROCm build requirements, architecture selection, and backend capabilities for cuPDLPx. +--- + +# GPU backends + +cuPDLPx supports NVIDIA GPUs through CUDA and AMD GPUs through ROCm/HIP. +Both backends use the same algorithm and Python, command-line, and C APIs. +The compatibility layer selects the GPU libraries. + +## Backend matrix + +| Target | Compiler | Dense algebra | Sparse algebra | Parallel primitives | +| --- | --- | --- | --- | --- | +| NVIDIA CUDA | NVCC | cuBLAS | cuSPARSE | CUB | +| AMD ROCm | `hipcc` | hipBLAS | hipSPARSE | hipCUB / rocPRIM | + +With `USE_HIP=ON`, compatibility headers map CUDA calls to HIP so that the +same `.cu` sources compile for AMD GPUs. + +See [SpMV](sparse-matrix-vector-products.md) for sparse matrix storage, +backend selection, and workspace reuse. + +## Select an architecture + +### CUDA + +CMake selects default architectures based on the CUDA version. To reduce +binary size, specify the target architecture: + +```bash +cmake -B build -DCMAKE_CUDA_ARCHITECTURES=90 +cmake --build build --parallel +``` + +Use an architecture supported by the installed CUDA toolkit and the deployment +GPU. + +### ROCm + +HIP builds require a target architecture. The default is `gfx90a`; override it +for other devices: + +```bash +cmake -B build \ + -DUSE_HIP=ON \ + -DCMAKE_HIP_ARCHITECTURES=gfx1100 \ + -DCMAKE_PREFIX_PATH=/opt/rocm +``` + +## Platform notes + +- CI builds both the CUDA and ROCm configurations on Linux. +- CUDA builds are also exercised on Windows. +- The native CLI is disabled automatically for MSVC builds because its current + argument parser depends on POSIX headers; the libraries and Python binding are + separate build targets. +- ROCm packages installed outside CMake's search path may require an explicit + `CMAKE_PREFIX_PATH`. + +## Vendor documentation + +- NVIDIA: [CUDA Toolkit](https://docs.nvidia.com/cuda/), + [cuBLAS](https://docs.nvidia.com/cuda/cublas/), + [cuSPARSE](https://docs.nvidia.com/cuda/cusparse/), and + [CUB](https://nvidia.github.io/cccl/unstable/cub/index.html). +- AMD: [ROCm](https://rocm.docs.amd.com/), + [hipBLAS](https://rocm.docs.amd.com/projects/hipBLAS/en/latest/), + [hipSPARSE](https://rocm.docs.amd.com/projects/hipSPARSE/en/latest/), and + [rocPRIM](https://rocm.docs.amd.com/projects/rocPRIM/en/latest/). +- The backend dispatch described here is implemented in + [`internal/cusparse_compat.h`](https://github.com/MIT-Lu-Lab/cuPDLPx/blob/main/internal/cusparse_compat.h) + and + [`src/spmv_backend.cu`](https://github.com/MIT-Lu-Lab/cuPDLPx/blob/main/src/spmv_backend.cu). diff --git a/docs/implementation/index.md b/docs/implementation/index.md new file mode 100644 index 0000000..422b638 --- /dev/null +++ b/docs/implementation/index.md @@ -0,0 +1,35 @@ +--- +description: Overview of the cuPDLPx solve pipeline and its GPU-oriented implementation. +--- + +# Implementation + +cuPDLPx preprocesses the LP on the CPU and runs the main iteration on the GPU. +Sparse matrices and iterate vectors remain in GPU memory until solution +recovery. + +## The solve pipeline + +| Stage | Main job | Where it runs | +| --- | --- | --- | +| Preprocess | Remove small matrix entries; map sufficiently negative lower bounds to $-\infty$ and sufficiently positive upper bounds to $+\infty$ ([thresholds](../reference/parameters.md#scaling-and-preprocessing)) | CPU | +| [Presolve](../algorithm/presolve.md) | Reduce the numbers of variables, constraints, and matrix nonzeros with PSLP | CPU | +| GPU setup | Transfer data, [preconditioning](../algorithm/preconditioning.md), allocate workspaces, estimate $\lVert A\rVert_2$ | CPU + GPU | +| Main iteration | Apply reflected Halpern PDHG between termination checks | GPU | +| Termination and restart | Form residuals, check termination and restart criteria, update primal weight | GPU + CPU | +| Optional polishing | Solve primal and dual feasibility problems | GPU | +| Recovery | Undo scaling and apply postsolve, then assemble the result | CPU + GPU | + +## Implementation details + +- [GPU backends](gpu-backends.md) covers CUDA, ROCm/HIP, and build targets. +- [SpMV](sparse-matrix-vector-products.md) explains why + both $A$ and $A^\top$ are stored and how the sparse backend is selected. +- [Kernel fusion](kernel-fusion.md) explains the fused primal and dual updates. +- [CUDA Graphs](cuda-graphs.md) explains capture, replay, and synchronization. +- [Results and status](../getting-started/results.md) documents the values + returned after recovery. +- [Log interpretation](../getting-started/log-interpretation.md) explains the + progress table and final solve summary. + +See [Algorithm](../algorithm/index.md) for the mathematical formulation. diff --git a/docs/implementation/kernel-fusion.md b/docs/implementation/kernel-fusion.md new file mode 100644 index 0000000..50ef0e3 --- /dev/null +++ b/docs/implementation/kernel-fusion.md @@ -0,0 +1,84 @@ +--- +description: How cuPDLPx fuses memory-bound PDHG vector operations into primal and dual GPU kernels. +--- + +# Kernel fusion + +cuPDLPx fuses the affine updates, projections, reflection, and Halpern +updates into two vector kernels. These operations have low arithmetic +intensity; fusion reduces global memory traffic and kernel launches. +See the [Roofline model](https://modal.com/gpu-glossary/perf/roofline-model) +for the relationship between arithmetic intensity and memory bandwidth. + +## Fused operations + +The fused updates combine PDHG, reflection, and Halpern anchoring, with +vector operations highlighted below. Within each epoch, $(x^0,y^0)$ is the +fixed anchor and $k$ is the local iteration index. + +$$ +\begin{aligned} +x^{k+1} +&=\frac{k+1}{k+2} +\left[ +\gamma\,\underbrace{\left( +2\htmlClass{vector-operation-highlight}{\operatorname{proj}_{\mathcal X}} +\left( +x^k\htmlClass{vector-operation-highlight}{\boldsymbol{-}} +\tau\left(c\htmlClass{vector-operation-highlight}{\boldsymbol{-}} +A^\top y^k +\right) +\right) +\htmlClass{vector-operation-highlight}{\boldsymbol{-}}x^k +\right)}_{\bar x^{k+1}} +\htmlClass{vector-operation-highlight}{\boldsymbol{+}} +(1-\gamma)x^k +\right] +\htmlClass{vector-operation-highlight}{\boldsymbol{+}} +\frac{1}{k+2}x^0,\\[0.6em] +y^{k+1} +&=\frac{k+1}{k+2} +\left[ +y^k +\htmlClass{vector-operation-highlight}{\boldsymbol{-}} +2\gamma\sigma +\left( +A\bar x^{k+1} +\htmlClass{vector-operation-highlight}{\boldsymbol{+}} +\htmlClass{vector-operation-highlight}{\operatorname{proj}_{-\mathcal S}} +\left( +\sigma^{-1}y^k +\htmlClass{vector-operation-highlight}{\boldsymbol{-}} +A\bar x^{k+1} +\right) +\right) +\right] +\htmlClass{vector-operation-highlight}{\boldsymbol{+}} +\frac{1}{k+2}y^0. +\end{aligned} +$$ + +The primal kernel computes $\bar x^{k+1}$ once and reuses it for both the +primal reflection and the PDHG extrapolation in the dual update. + +## Kernel execution + +Each iteration uses two SpMV calls and two fused vector kernels: + +1. **SpMV:** compute $A^\top y^k$ for the primal update. +2. **Fused primal kernel:** perform the primal PDHG, reflection, and Halpern updates. +3. **SpMV:** compute $A\bar x^{k+1}$ for the dual update. +4. **Fused dual kernel:** perform the dual PDHG, reflection, and Halpern updates. + +Each GPU thread computes all three updates for one coordinate. +The primal kernel writes $x^{k+1}$ and +$\bar x^{k+1}$ to global memory, and the dual kernel writes $y^{k+1}$. +The second SpMV uses the stored vector $\bar x^{k+1}$. + +!!! note "Major iterations" + + cuPDLPx evaluates termination and restart criteria at regular intervals; + these iterations are called **major iterations**. For these checks, the + fused kernels also save the PDHG iterates and reflected dual values to global + memory. These values otherwise remain in registers. The primal kernel also performs + [dual-slack recovery](../algorithm/termination.md#dual-slack-recovery). diff --git a/docs/implementation/sparse-matrix-vector-products.md b/docs/implementation/sparse-matrix-vector-products.md new file mode 100644 index 0000000..14229b8 --- /dev/null +++ b/docs/implementation/sparse-matrix-vector-products.md @@ -0,0 +1,67 @@ +--- +description: Sparse matrix storage and CUDA and ROCm sparse backends in cuPDLPx. +--- + +# SpMV + +The two sparse matrix–vector products (SpMV) in the PDHG update are +highlighted below: + +$$ +\begin{aligned} +\widehat x^{k+1} +&=\operatorname{proj}_{\mathcal X} +\left( +x^k-\tau\left(c- +\htmlClass{spmv-highlight}{\boldsymbol{A^\top y^k}} +\right) +\right),\\[0.4em] +\widehat y^{k+1} +&=y^k +-\sigma\htmlClass{spmv-highlight}{\boldsymbol{A(2\widehat x^{k+1}-x^k)}} +-\sigma\operatorname{proj}_{-\mathcal S} +\left( +\sigma^{-1}y^k +-\htmlClass{spmv-highlight}{\boldsymbol{A(2\widehat x^{k+1}-x^k)}} +\right). +\end{aligned} +$$ + +cuPDLPx computes $A(2\widehat x^{k+1}-x^k)$ once and reuses the result in both +terms of the dual update. + +cuPDLPx stores both $A$ and $A^\top$ in compressed sparse row (CSR) form. +This increases matrix storage, but lets both products use optimized +non-transposed SpMV kernels throughout the solve. + +## CUDA sparse backends + +For CUDA versions before 13.3, cuPDLPx uses `cusparseSpMV` with +`CUSPARSE_SPMV_CSR_ALG2`. This path is deterministic across runs. The solver +also calls `cusparseSpMV_preprocess` so the one-time analysis cost is amortized +over the repeated products. + +Starting with CUDA 13.3, cuPDLPx uses `cusparseSpMVOp`. On the LP workloads +reported in the cuPDLPx paper, this backend is often faster, but it may request +a larger temporary buffer. + +The sparse descriptors, preprocessing state, and temporary buffers are created +during initialization and reused. No sparse workspace allocation occurs in the +main iteration loop. + +See the [cuSPARSE documentation for CUDA +13.3](https://docs.nvidia.com/cuda/archive/13.3.0/cusparse/contents.html) +for `cusparseSpMV`, `cusparseSpMV_preprocess`, and `cusparseSpMVOp`. + +## ROCm/HIP + +The HIP build uses `hipsparseSpMV` for the same two products. `cusparseSpMVOp` +is CUDA-specific, so it is not selected on the HIP path. The higher-level +iteration and solver interfaces are shared between both backends. + +Build requirements and architecture selection are listed under +[GPU backends](gpu-backends.md). + +The AMD [hipSPARSE generic API +reference](https://rocm.docs.amd.com/projects/hipSPARSE/en/latest/reference/generic.html) +documents `hipsparseSpMV`, preprocessing, and workspace requirements. diff --git a/docs/index.md b/docs/index.md new file mode 100644 index 0000000..98db30c --- /dev/null +++ b/docs/index.md @@ -0,0 +1,29 @@ +--- +title: Home +description: cuPDLPx is an open-source, GPU-accelerated first-order solver for large-scale linear programming. +hide: + - navigation + - toc +--- + +
+ +
+ +
+ +# cuPDLPx: Solve linear programs at GPU scale. + +An open‑source, GPU‑accelerated first‑order solver for large‑scale linear programming.
+Supports NVIDIA GPUs via CUDA and AMD GPUs via ROCm. + +
+ Get started + Algorithm overview +
+ +
+ +
+ +
diff --git a/docs/javascripts/katex.js b/docs/javascripts/katex.js new file mode 100644 index 0000000..4439062 --- /dev/null +++ b/docs/javascripts/katex.js @@ -0,0 +1,11 @@ +document$.subscribe(({ body }) => { + renderMathInElement(body, { + delimiters: [ + { left: "$$", right: "$$", display: true }, + { left: "$", right: "$", display: false }, + { left: "\\(", right: "\\)", display: false }, + { left: "\\[", right: "\\]", display: true }, + ], + trust: (context) => context.command === "\\htmlClass", + }); +}); diff --git a/docs/reference/parameters.md b/docs/reference/parameters.md new file mode 100644 index 0000000..b20358f --- /dev/null +++ b/docs/reference/parameters.md @@ -0,0 +1,332 @@ +--- +description: cuPDLPx parameter names, defaults, ranges, and usage across Python, Julia, C, and the command line. +hide: + - navigation +--- + +# Parameters + +Solver settings are grouped by function. Each table lists the parameter +name, default, type, and allowed values for the selected interface. + +## Usage + +=== "Python" + + Set parameters with `setParams`, `setParam`, or `Params`: + + ```python + model.setParams( + TimeLimit=300, + OptimalityTol=1e-6, + FeasibilityTol=1e-6, + OutputFlag=True, + ) + ``` + + Python also accepts the leaf C field names, such as `time_sec_limit` and + `eps_optimal_relative`. + + Read and reset values with: + + ```python + time_limit = model.getParam("TimeLimit") + same_value = model.Params.TimeLimit + + for key, value in model.Params.items(): + print(key, value) + + model.resetParams() + ``` + + `setParams(...)` validates all updates before applying them. If any name + or value is invalid, no parameters are changed. + +=== "Julia" + + CuPDLPx.jl exposes supported solver settings as JuMP attributes using the + C leaf field names: + + ```julia + using JuMP, CuPDLPx + + model = Model(CuPDLPx.Optimizer) + set_attribute(model, "time_sec_limit", 300.0) + set_attribute(model, "eps_optimal_relative", 1e-6) + set_attribute(model, "eps_feasible_relative", 1e-6) + set_attribute(model, "verbose", true) + ``` + + See the [CuPDLPx.jl documentation](https://jump.dev/JuMP.jl/stable/packages/CuPDLPx/#Supported-parameters) + for the attributes exposed by the current Julia wrapper. + +=== "C" + + Initialize `pdhg_parameters_t` before overriding individual fields: + + ```c + pdhg_parameters_t params; + set_default_parameters(¶ms); + + params.termination_criteria.time_sec_limit = 300.0; + params.termination_criteria.eps_optimal_relative = 1e-6; + params.termination_criteria.eps_feasible_relative = 1e-6; + params.verbose = false; + ``` + + Use values in the documented ranges; assigning C fields does not perform + the Python interface's runtime parameter validation. + +=== "Command line" + + Pass CLI settings before the input and output paths: + + ```bash + cupdlpx \ + --time_limit 300 \ + --eps_opt 1e-6 \ + --eps_feas 1e-6 \ + --quiet \ + model.mps.gz output + ``` + + Run `cupdlpx --help` for the full list of flags. + +
+ +## Limits and logging + +=== "Python" + + | Parameter | Default | Type and range | Description | + | --- | --- | --- | --- | + | `TimeLimit` | `3600.0` | float
$[0,\infty)$ | Wall-clock limit in seconds. | + | `IterationLimit` | $2^{31}-1$ | integer
$[0,2^{31}-1]$ | Maximum number of iterations. | + | `OutputFlag`, `LogToConsole` | `true` | boolean
`false`, `true` | Print progress and the solve summary. | + | `TermCheckFreq` | `200` | integer
$[3,2^{31}-1]$ | Iterations between termination and restart checks. | + | `Debug` | `false` | boolean
`false`, `true` | Print developer diagnostics: restart reasons, active-set step size boost events, and a primal-weight log column. | + +=== "C and Julia" + + | Parameter | Default | Type and range | Description | + | --- | --- | --- | --- | + | `time_sec_limit` | `3600.0` | float
$[0,\infty)$ | Wall-clock limit in seconds. | + | `iteration_limit` | $2^{31}-1$ | integer
$[0,2^{31}-1]$ | Maximum number of iterations. | + | `verbose` | `true` | boolean
`false`, `true` | Print progress and the solve summary. | + | `termination_evaluation_frequency` | `200` | integer
$[3,2^{31}-1]$ | Iterations between termination and restart checks. | + | `debug` | `false` | boolean
`false`, `true` | Print developer diagnostics: restart reasons, active-set step size boost events, and a primal-weight log column. | + +=== "Command line" + + | Parameter | Default | Type and range | Description | + | --- | --- | --- | --- | + | `--time_limit` | `3600.0` | float
$[0,\infty)$ | Wall-clock limit in seconds. | + | `--iter_limit` | $2^{31}-1$ | integer
$[0,2^{31}-1]$ | Maximum number of iterations. | + | `--verbose`, `--quiet` | `true` | boolean
`false`, `true` | Print progress and the solve summary. | + | `--eval_freq` | `200` | integer
$[3,2^{31}-1]$ | Iterations between termination and restart checks. | + | `--debug` | `false` | boolean
`false`, `true` | Print developer diagnostics: restart reasons, active-set step size boost events, and a primal-weight log column. | + +## Accuracy and termination + +=== "Python" + + | Parameter | Default | Type and range | Description | + | --- | --- | --- | --- | + | `OptimalityTol` | `1e-4` | float
$(0,\infty)$ | Relative objective-gap tolerance. | + | `FeasibilityTol` | `1e-4` | float
$(0,\infty)$ | Relative primal and dual feasibility tolerance. | + | `OptimalityNorm` | `l2` | string
`l2`, `linf` | Norm used by the termination criteria. | + +=== "C and Julia" + + | Parameter | Default | Type and range | Description | + | --- | --- | --- | --- | + | `eps_optimal_relative` | `1e-4` | float
$(0,\infty)$ | Relative objective-gap tolerance. | + | `eps_feasible_relative` | `1e-4` | float
$(0,\infty)$ | Relative primal and dual feasibility tolerance. | + | `optimality_norm` | `NORM_TYPE_L2` | `norm_type_t`
`NORM_TYPE_L2`, `NORM_TYPE_L_INF` | Norm used by the termination criteria. | + +=== "Command line" + + | Parameter | Default | Type and range | Description | + | --- | --- | --- | --- | + | `--eps_opt` | `1e-4` | float
$(0,\infty)$ | Relative objective-gap tolerance. | + | `--eps_feas` | `1e-4` | float
$(0,\infty)$ | Relative primal and dual feasibility tolerance. | + | `--opt_norm` | `l2` | string
`l2`, `linf` | Norm used by the termination criteria. | + +Residuals are evaluated in the unscaled model; when +[presolve](../algorithm/presolve.md) is enabled, this is the presolved model. + +## Scaling and preprocessing + +=== "Python" + + | Parameter | Default | Type and range | Description | + | --- | --- | --- | --- | + | `GeoMeanIters` | `12` | integer
$[0,2^{31}-1]$ | Number of [geometric mean scaling](../algorithm/preconditioning.md#geometric-mean-scaling) passes; `0` disables this stage. | + | `RuizIters` | `10` | integer
$[0,2^{31}-1]$ | Number of $\ell_\infty$ [Ruiz equilibration scaling](../algorithm/preconditioning.md#ruiz-equilibration-scaling) passes. | + | `UsePCAlpha` | `true` | boolean
`false`, `true` | Enable Pock–Chambolle scaling. | + | `PCAlpha` | `1.0` | float
$\mathbb R$ | Exponent used by Pock–Chambolle scaling. | + | `BoundObjRescaling` | `true` | boolean
`false`, `true` | Enable objective and bound scaling. | + | `Presolve` | `true` | boolean
`false`, `true` | Enable presolve. | + | `MatrixZeroTol` | `1e-9` | float
$[0,\infty)$ | Remove entries of $A$ at or below this absolute magnitude. | + | `InfiniteBound` | `1e20` | float
$(0,\infty)$ | Replace lower bounds $\le -t$ and upper bounds $\ge t$ with infinities, where $t$ is this threshold. | + +=== "C and Julia" + + | Parameter | Default | Type and range | Description | + | --- | --- | --- | --- | + | `geometric_mean_iterations` | `12` | integer
$[0,2^{31}-1]$ | Number of [geometric mean scaling](../algorithm/preconditioning.md#geometric-mean-scaling) passes; `0` disables this stage. | + | `l_inf_ruiz_iterations` | `10` | integer
$[0,2^{31}-1]$ | Number of $\ell_\infty$ [Ruiz equilibration scaling](../algorithm/preconditioning.md#ruiz-equilibration-scaling) passes. | + | `has_pock_chambolle_alpha` | `true` | boolean
`false`, `true` | Enable Pock–Chambolle scaling. | + | `pock_chambolle_alpha` | `1.0` | float
$\mathbb R$ | Exponent used by Pock–Chambolle scaling. | + | `bound_objective_rescaling` | `true` | boolean
`false`, `true` | Enable objective and bound scaling. | + | `presolve` | `true` | boolean
`false`, `true` | Enable presolve. | + | `matrix_zero_tol` | `1e-9` | float
$[0,\infty)$ | Remove entries of $A$ at or below this absolute magnitude. | + | `infinite_bound` | `1e20` | float
$(0,\infty)$ | Replace lower bounds $\le -t$ and upper bounds $\ge t$ with infinities, where $t$ is this threshold. | + +=== "Command line" + + | Parameter | Default | Type and range | Description | + | --- | --- | --- | --- | + | `--geo_mean_iter` | `12` | integer
$[0,2^{31}-1]$ | Number of [geometric mean scaling](../algorithm/preconditioning.md#geometric-mean-scaling) passes; `0` disables this stage. | + | `--l_inf_ruiz_iter` | `10` | integer
$[0,2^{31}-1]$ | Number of $\ell_\infty$ [Ruiz equilibration scaling](../algorithm/preconditioning.md#ruiz-equilibration-scaling) passes. | + | `--no_pock_chambolle` | Not passed | flag | Disable Pock–Chambolle scaling (enabled by default). | + | `--pock_chambolle_alpha` | `1.0` | float
$\mathbb R$ | Exponent used by Pock–Chambolle scaling. | + | `--no_bound_obj_rescaling` | Not passed | flag | Disable objective and bound scaling (enabled by default). | + | `--no_presolve` | Not passed | flag | Disable presolve (enabled by default). | + | `--matrix_zero_tol` | `1e-9` | float
$[0,\infty)$ | Remove entries of $A$ at or below this absolute magnitude. | + | `--infinite_bound` | `1e20` | float
$(0,\infty)$ | Replace lower bounds $\le -t$ and upper bounds $\ge t$ with infinities, where $t$ is this threshold. | + +## Step size and reflection + +=== "Python" + + | Parameter | Default | Type and range | Description | + | --- | --- | --- | --- | + | `SVMaxIter` | `5000` | integer
$[1,2^{31}-1]$ | Maximum iterations for estimating $\lVert A\rVert_2$. | + | `SVTol` | `1e-4` | float
$(0,\infty)$ | Stopping tolerance for the norm estimate. | + | `ReflectionCoeff` | `1.0` | float
$\mathbb R$ | Weight $\gamma$ on the reflected point in the Halpern update; `1.0` uses the fully reflected point $2\widehat x-x$. | + +=== "C and Julia" + + | Parameter | Default | Type and range | Description | + | --- | --- | --- | --- | + | `sv_max_iter` | `5000` | integer
$[1,2^{31}-1]$ | Maximum iterations for estimating $\lVert A\rVert_2$. | + | `sv_tol` | `1e-4` | float
$(0,\infty)$ | Stopping tolerance for the norm estimate. | + | `reflection_coefficient` | `1.0` | float
$\mathbb R$ | Weight $\gamma$ on the reflected point in the Halpern update; `1.0` uses the fully reflected point $2\widehat x-x$. | + +=== "Command line" + + | Parameter | Default | Type and range | Description | + | --- | --- | --- | --- | + | `--sv_max_iter` | `5000` | integer
$[1,2^{31}-1]$ | Maximum iterations for estimating $\lVert A\rVert_2$. | + | `--sv_tol` | `1e-4` | float
$(0,\infty)$ | Stopping tolerance for the norm estimate. | + + The command-line interface does not expose the reflection coefficient. + +## Active-set step size boost + +=== "Python" + + | Parameter | Default | Type and range | Description | + | --- | --- | --- | --- | + | `ActiveSetBoost` | `true` | boolean
`false`, `true` | Enable the [active-set step size boost](../algorithm/step-size.md#active-set-step-size-boost). | + | `ASBActivationTol` | `1e-4` | float
$[0,\infty)$ | Threshold for the maximum relative primal residual, dual residual, and primal–dual gap; see [activation](../algorithm/step-size.md#activation). | + | `ASBWindowIter` | `10000` | integer
$[1,2^{31}-1]$ | Length of the [trailing window](../algorithm/step-size.md#trailing-window) in iterations. | + | `ASBVariableTol` | `1e-8` | float
$[0,\infty)$ | Dual-slack tolerance for identifying a variable at a bound. | + | `ASBConstraintTol` | `1e-8` | float
$[0,\infty)$ | Margin for identifying an inactive constraint. | + | `ASBReestimateChangeRatio` | `0.01` | float
$[0,\infty)$ | Threshold for accumulated additions and removals, as a fraction of the current active-set size, before norm re-estimation. | + | `ASBSafetyFactor` | `0.9` | float
$(0,\infty)$ | Safety factor $\alpha$ in the target step $\alpha/\hat\sigma$. | + | `ASBMinRaiseRatio` | `1.1` | float
$[1,\infty)$ | Minimum ratio of the target step to the current step for a step-size increase. | + | `ASBDivergenceMargin` | `0.05` | float
$[0,\infty)$ | Allowed relative increase in fixed-point error within an epoch before rollback. | + | `ASBDivergenceCeilingRatio` | `0.7` | float
$(0,1]$ | Limit on subsequent step-size increases, as a fraction of the rejected step size; the step never falls below its initial value. | + | `ASBMaxReverts` | `2` | integer
$[0,2^{31}-1]$ | Number of rollbacks before ASB is disabled. | + +=== "C and Julia" + + | Parameter | Default | Type and range | Description | + | --- | --- | --- | --- | + | `active_set_boost` | `true` | boolean
`false`, `true` | Enable the [active-set step size boost](../algorithm/step-size.md#active-set-step-size-boost). | + | `asb_activation_tol` | `1e-4` | float
$[0,\infty)$ | Threshold for the maximum relative primal residual, dual residual, and primal–dual gap; see [activation](../algorithm/step-size.md#activation). | + | `asb_window_iter` | `10000` | integer
$[1,2^{31}-1]$ | Length of the [trailing window](../algorithm/step-size.md#trailing-window) in iterations. | + | `asb_variable_tol` | `1e-8` | float
$[0,\infty)$ | Dual-slack tolerance for identifying a variable at a bound. | + | `asb_constraint_tol` | `1e-8` | float
$[0,\infty)$ | Margin for identifying an inactive constraint. | + | `asb_reestimate_change_ratio` | `0.01` | float
$[0,\infty)$ | Threshold for accumulated additions and removals, as a fraction of the current active-set size, before norm re-estimation. | + | `asb_safety_factor` | `0.9` | float
$(0,\infty)$ | Safety factor $\alpha$ in the target step $\alpha/\hat\sigma$. | + | `asb_min_raise_ratio` | `1.1` | float
$[1,\infty)$ | Minimum ratio of the target step to the current step for a step-size increase. | + | `asb_divergence_margin` | `0.05` | float
$[0,\infty)$ | Allowed relative increase in fixed-point error within an epoch before rollback. | + | `asb_divergence_ceiling_ratio` | `0.7` | float
$(0,1]$ | Limit on subsequent step-size increases, as a fraction of the rejected step size; the step never falls below its initial value. | + | `asb_max_reverts` | `2` | integer
$[0,2^{31}-1]$ | Number of rollbacks before ASB is disabled. | + +=== "Command line" + + | Parameter | Default | Type and range | Description | + | --- | --- | --- | --- | + | `--no_active_set_boost` | Not passed | flag | Disable the [active-set step size boost](../algorithm/step-size.md#active-set-step-size-boost) (enabled by default). | + | `--asb_activation_tol` | `1e-4` | float
$[0,\infty)$ | Threshold for the maximum relative primal residual, dual residual, and primal–dual gap; see [activation](../algorithm/step-size.md#activation). | + | `--asb_window_iter` | `10000` | integer
$[1,2^{31}-1]$ | Length of the [trailing window](../algorithm/step-size.md#trailing-window) in iterations. | + | `--asb_variable_tol` | `1e-8` | float
$[0,\infty)$ | Dual-slack tolerance for identifying a variable at a bound. | + | `--asb_constraint_tol` | `1e-8` | float
$[0,\infty)$ | Margin for identifying an inactive constraint. | + | `--asb_reestimate_change_ratio` | `0.01` | float
$[0,\infty)$ | Threshold for accumulated additions and removals, as a fraction of the current active-set size, before norm re-estimation. | + | `--asb_safety_factor` | `0.9` | float
$(0,\infty)$ | Safety factor $\alpha$ in the target step $\alpha/\hat\sigma$. | + | `--asb_min_raise_ratio` | `1.1` | float
$[1,\infty)$ | Minimum ratio of the target step to the current step for a step-size increase. | + | `--asb_divergence_margin` | `0.05` | float
$[0,\infty)$ | Allowed relative increase in fixed-point error within an epoch before rollback. | + | `--asb_divergence_ceiling_ratio` | `0.7` | float
$(0,1]$ | Limit on subsequent step-size increases, as a fraction of the rejected step size; the step never falls below its initial value. | + | `--asb_max_reverts` | `2` | integer
$[0,2^{31}-1]$ | Number of rollbacks before ASB is disabled. | + +## Adaptive restart + +=== "Python" + + | Parameter | Default | Type and range | Description | + | --- | --- | --- | --- | + | `RestartArtificialThresh` | `0.36` | float
$\mathbb R$ | Maximum epoch length relative to the total iteration count. | + | `RestartSufficientReduction` | `0.2` | float
$\mathbb R$ | Fixed-point-error ratio for a restart triggered by sufficient reduction. | + | `RestartNecessaryReduction` | `0.5` | float
$\mathbb R$ | Required reduction before a local-increase restart. | + | `RestartKp` | `0.99` | float
$\mathbb R$ | Proportional gain in the primal-weight controller. | + + The Python interface does not expose the integral gain, derivative gain, + or integral decay. Their defaults are listed in the C and Julia tab. + +=== "C and Julia" + + | Parameter | Default | Type and range | Description | + | --- | --- | --- | --- | + | `artificial_restart_threshold` | `0.36` | float
$\mathbb R$ | Maximum epoch length relative to the total iteration count. | + | `sufficient_reduction_for_restart` | `0.2` | float
$\mathbb R$ | Fixed-point-error ratio for a restart triggered by sufficient reduction. | + | `necessary_reduction_for_restart` | `0.5` | float
$\mathbb R$ | Required reduction before a local-increase restart. | + | `k_p` | `0.99` | float
$\mathbb R$ | Proportional gain in the primal-weight controller. | + | `k_i` | `0.01` | float
$\mathbb R$ | Integral gain in the primal-weight controller. | + | `k_d` | `0.0` | float
$\mathbb R$ | Derivative gain in the primal-weight controller. | + | `i_smooth` | `0.3` | float
$\mathbb R$ | Decay factor for the accumulated integral error. | + +=== "Command line" + + These settings are not exposed by the command-line interface. + +See [Adaptive restart](../algorithm/restart.md) for the restart conditions and +[Primal weight](../algorithm/primal-weight.md) for the controller update. + +## Feasibility polishing + +=== "Python" + + | Parameter | Default | Type and range | Description | + | --- | --- | --- | --- | + | `FeasibilityPolishing` | `false` | boolean
`false`, `true` | Run the feasibility-polishing phases after the main solve. | + | `FeasibilityPolishingTol` | `1e-6` | float
$(0,\infty)$ | Target relative feasibility residual. | + +=== "C and Julia" + + | Parameter | Default | Type and range | Description | + | --- | --- | --- | --- | + | `feasibility_polishing` | `false` | boolean
`false`, `true` | Run the feasibility-polishing phases after the main solve. | + | `eps_feas_polish_relative` | `1e-6` | float
$(0,\infty)$ | Target relative feasibility residual. | + +=== "Command line" + + | Parameter | Default | Type and range | Description | + | --- | --- | --- | --- | + | `-f`, `--feasibility_polishing` | `false` | boolean
`false`, `true` | Run the feasibility-polishing phases after the main solve. | + | `--eps_feas_polish` | `1e-6` | float
$(0,\infty)$ | Target relative feasibility residual. | + +
diff --git a/docs/stylesheets/extra.css b/docs/stylesheets/extra.css new file mode 100644 index 0000000..05fcc59 --- /dev/null +++ b/docs/stylesheets/extra.css @@ -0,0 +1,540 @@ +/* Small cuPDLPx-specific adjustments to the Modern theme. */ +.md-header__button.md-logo { + display: none; +} + +/* Documentation-first landing page, inspired by Zensical's home layout. */ +body:has(.cupdlpx-home) { + --cupdlpx-home-chrome-height: 4.8rem; + --cupdlpx-ink-strong: hsl(218 40% 12%); + --cupdlpx-ink: hsl(218 28% 24%); + --cupdlpx-ink-muted: hsl(218 18% 40%); + --cupdlpx-space-xs: 0.25rem; + --cupdlpx-space-sm: 0.5rem; + --cupdlpx-space-md: 1rem; + --cupdlpx-space-lg: 1.5rem; + --cupdlpx-space-xl: 2rem; + --cupdlpx-space-2xl: 3rem; + --cupdlpx-space-3xl: 4rem; + --md-typeset-a-color: #165bb8; + background: var(--md-default-bg-color); +} + +body[data-md-color-scheme="slate"]:has(.cupdlpx-home) { + --cupdlpx-ink-strong: hsl(215 25% 94%); + --cupdlpx-ink: hsl(215 20% 84%); + --cupdlpx-ink-muted: hsl(215 14% 66%); + --md-typeset-a-color: #82b1ff; +} + +body:has(.cupdlpx-home) .md-header:not(.md-header--shadow), +body:has(.cupdlpx-home) .md-tabs { + background-color: transparent; + box-shadow: none; + color: var(--cupdlpx-ink-strong); +} + +body:has(.cupdlpx-home) .md-tabs__link { + color: var(--cupdlpx-ink-muted); +} + +body:has(.cupdlpx-home) .md-tabs__item--active .md-tabs__link, +body:has(.cupdlpx-home) .md-tabs__link:hover { + color: var(--cupdlpx-ink-strong); +} + +body:has(.cupdlpx-home) .md-main__inner { + margin: 0; + max-width: none; + padding: 0; +} + +body:has(.cupdlpx-home) .md-content { + max-width: none; +} + +body:has(.cupdlpx-home) .md-content__inner { + margin: 0; + max-width: none; + padding: 0; +} + +body:has(.cupdlpx-home) .md-content__button { + display: none; +} + +/* The homepage has its own calls to action, so omit the next-page prompt. */ +body:has(.cupdlpx-home) .md-footer__inner { + display: none; +} + +.md-typeset .cupdlpx-home { + color: var(--cupdlpx-ink); +} + +.md-typeset .cupdlpx-home a { + color: var(--md-typeset-a-color); +} + +.md-typeset .cupdlpx-home__hero { + align-items: center; + box-sizing: border-box; + display: flex; + margin: 0; + min-height: calc(80svh - var(--cupdlpx-home-chrome-height)); + padding: clamp( + var(--cupdlpx-space-2xl), + 5vh, + var(--cupdlpx-space-3xl) + ) + 0 var(--cupdlpx-space-3xl); +} + +.md-typeset .cupdlpx-home__hero-inner { + margin: 0 auto; + max-width: 61rem; + width: calc(100% - 3rem); +} + +.md-typeset .cupdlpx-home__hero h1 { + color: var(--cupdlpx-ink-strong); + font-size: clamp(2.25rem, 3.5vw, 3.2rem); + font-weight: 700; + letter-spacing: -0.045em; + line-height: 1.12; + margin: 0; + max-width: none; + text-wrap: balance; + white-space: normal; +} + +.md-typeset .cupdlpx-home .headerlink { + display: none; +} + +.md-typeset .cupdlpx-home__hero-inner > p { + color: var(--cupdlpx-ink-muted); + font-size: clamp(0.95rem, 1.1vw, 1.05rem); + line-height: 1.55; + margin: var(--cupdlpx-space-lg) 0 0; + max-width: none; +} + +.md-typeset .cupdlpx-home__tagline { + white-space: nowrap; +} + +.md-typeset .cupdlpx-home__actions { + align-items: center; + display: flex; + flex-wrap: wrap; + gap: var(--cupdlpx-space-md); + margin: var(--cupdlpx-space-xl) 0 0; +} + +.md-typeset .cupdlpx-home__actions .md-button { + border: 0; + border-radius: 0.25rem; + color: var(--cupdlpx-ink); + margin: 0; + padding: var(--cupdlpx-space-sm) var(--cupdlpx-space-md); + text-decoration: none; +} + +.md-typeset .cupdlpx-home__actions .md-button--primary { + background: var(--cupdlpx-ink-strong); + color: var(--md-default-bg-color); +} + +.md-typeset .cupdlpx-home__actions .md-button--primary:hover { + background: var(--cupdlpx-ink); + color: var(--md-default-bg-color); +} + +.md-typeset .cupdlpx-home__actions .md-button:focus-visible { + outline: 0.15rem solid var(--md-typeset-a-color); + outline-offset: var(--cupdlpx-space-xs); +} + +.md-typeset .cupdlpx-home__actions .md-button:not(.md-button--primary) { + background: transparent; + padding-inline: var(--cupdlpx-space-sm); +} + +.md-typeset + .cupdlpx-home__actions + .md-button:not(.md-button--primary)::after { + content: " →"; + display: inline-block; + margin-left: var(--cupdlpx-space-xs); + transition: transform 150ms ease; +} + +.md-typeset + .cupdlpx-home__actions + .md-button:not(.md-button--primary):hover { + color: var(--md-typeset-a-color); +} + +.md-typeset + .cupdlpx-home__actions + .md-button:not(.md-button--primary):hover::after { + transform: translateX(0.15rem); +} + +@media screen and (max-width: 76.2344em) { + body:has(.cupdlpx-home) { + --cupdlpx-home-chrome-height: 2.4rem; + } +} + +@media screen and (max-width: 44.9844em) { + .md-typeset .cupdlpx-home__tagline { + white-space: normal; + } + + .md-typeset .cupdlpx-home__hero { + padding-block: var(--cupdlpx-space-2xl); + } +} + +/* KaTeX auto-render wraps display math in a span; avoid extra line boxes. */ +.md-typeset div.arithmatex > span { + display: block; +} + +.md-typeset div.arithmatex .katex-display { + margin: 0.4em 0; + overflow-x: auto; + overflow-y: hidden; + padding-block: 0.05em; +} + +.md-typeset .spmv-highlight, +.md-typeset .vector-operation-highlight { + color: #c62828; + font-weight: 700; +} + +[data-md-color-scheme="slate"] .md-typeset .spmv-highlight, +[data-md-color-scheme="slate"] .md-typeset .vector-operation-highlight { + color: #ff6b6b; +} + +/* Keep model argument names readable as single identifiers. */ +.md-typeset .model-arguments td:first-child code { + white-space: nowrap; +} + +.md-typeset .parameter-reference td:first-child code { + background-color: rgb(21 101 192 / 8%); + color: #1565c0; + font-weight: 600; +} + +[data-md-color-scheme="slate"] + .md-typeset .parameter-reference td:first-child code { + background-color: rgb(100 181 246 / 12%); + color: #64b5f6; +} + +.md-typeset .parameter-reference th, +.md-typeset .parameter-reference td { + padding: 0.75em 0.9em; + vertical-align: top; +} + +.md-typeset .parameter-reference .md-typeset__table { + display: block; +} + +.md-typeset .parameter-reference table { + display: table; + width: 100%; +} + +.md-typeset .parameter-reference th:first-child, +.md-typeset .parameter-reference td:first-child { + width: 28%; +} + +.md-typeset .parameter-reference th:nth-child(2), +.md-typeset .parameter-reference td:nth-child(2) { + width: 10%; +} + +.md-typeset .parameter-reference th:nth-child(3), +.md-typeset .parameter-reference td:nth-child(3) { + width: 16%; +} + +.md-typeset .parameter-reference th:last-child, +.md-typeset .parameter-reference td:last-child { + width: 46%; +} + +.md-typeset .parameter-reference th, +.md-typeset .parameter-reference td:nth-child(2), +.md-typeset .parameter-reference td code { + white-space: nowrap; +} + +@media screen and (max-width: 59.9844em) { + .md-typeset :is(.parameter-reference, .status-guide) .md-typeset__table { + display: block; + width: 100%; + } + + .md-typeset .status-guide table, + .md-typeset .status-guide tbody, + .md-typeset .parameter-reference table, + .md-typeset .parameter-reference tbody { + display: block; + width: 100%; + } + + .md-typeset :is(.parameter-reference, .status-guide) thead { + clip: rect(0 0 0 0); + clip-path: inset(50%); + display: block; + height: 1px; + overflow: hidden; + position: absolute; + white-space: nowrap; + width: 1px; + } + + .md-typeset :is(.parameter-reference, .status-guide) tbody tr { + border-top: 0.05rem solid var(--md-typeset-table-color); + display: grid; + gap: 0.75rem; + grid-template-columns: minmax(0, 1fr) minmax(0, 1fr); + padding: 0.75rem; + } + + .md-typeset .status-guide tbody tr { + grid-template-columns: minmax(0, 1fr); + } + + .md-typeset :is(.parameter-reference, .status-guide) tbody tr:first-child { + border-top: 0; + } + + .md-typeset .status-guide table td, + .md-typeset .parameter-reference td, + .md-typeset .parameter-reference td:first-child, + .md-typeset .parameter-reference td:nth-child(2), + .md-typeset .parameter-reference td:nth-child(3), + .md-typeset .parameter-reference td:last-child { + border-top: 0; + display: block; + padding: 0; + white-space: normal; + width: auto; + } + + .md-typeset .parameter-reference td code { + overflow-wrap: anywhere; + white-space: normal; + } + + .md-typeset .parameter-reference td:first-child, + .md-typeset .parameter-reference td:last-child { + grid-column: 1 / -1; + } + + .md-typeset .status-guide td:first-child { + font-size: 1.1em; + font-weight: 700; + } + + .md-typeset :is(.parameter-reference, .status-guide) td::before { + color: var(--md-default-fg-color--light); + display: block; + font-size: 0.75em; + font-weight: 700; + letter-spacing: 0.025em; + line-height: 1.4; + margin-bottom: 0.2rem; + text-transform: uppercase; + } + + .md-typeset .parameter-reference td:first-child::before { + content: "Parameter"; + } + + .md-typeset .parameter-reference td:nth-child(2)::before { + content: "Default"; + } + + .md-typeset .parameter-reference td:nth-child(3)::before { + content: "Type and range"; + } + + .md-typeset .parameter-reference td:last-child::before { + content: "Description"; + } + + .md-typeset .status-guide td:nth-child(2)::before { + content: "Meaning"; + } + + .md-typeset .status-guide td:last-child::before { + content: "Next step"; + } +} + +.md-typeset .status-reference td code { + overflow-wrap: anywhere; + white-space: normal; +} + +/* Paper-style pseudocode with live, searchable KaTeX formulas. */ +.md-typeset .latex-algorithm { + border-block: 0.12rem solid var(--md-default-fg-color); + margin: 1.25rem 0; +} + +.md-typeset .latex-algorithm__caption { + border-bottom: 0.05rem solid var(--md-default-fg-color--lightest); + padding: 0.45rem 0.7rem; +} + +.md-typeset .latex-algorithm__input { + border-bottom: 0.05rem solid var(--md-default-fg-color--lightest); + padding: 0.45rem 0.7rem; +} + +.md-typeset .latex-algorithm__lines { + counter-reset: algorithm-line; + list-style: none; + margin: 0; + padding: 0.5rem 0; +} + +.md-typeset .latex-algorithm__lines > li { + counter-increment: algorithm-line; + line-height: 1.55; + margin: 0; + padding: 0.18rem 0.8rem 0.18rem 2.8rem; + position: relative; +} + +.md-typeset .latex-algorithm__lines > li::before { + color: var(--md-default-fg-color--light); + content: counter(algorithm-line); + font-variant-numeric: tabular-nums; + left: 0.45rem; + position: absolute; + text-align: right; + width: 1.45rem; +} + +.md-typeset .latex-algorithm__lines > li.latex-algorithm__indent-1 { + padding-left: 4rem; +} + +.md-typeset .latex-algorithm__lines > li.latex-algorithm__indent-1::after, +.md-typeset .latex-algorithm__lines > li.latex-algorithm__indent-2::after { + border-left: 0.05rem solid var(--md-default-fg-color--lightest); + bottom: 0; + content: ""; + left: 2.75rem; + position: absolute; + top: 0; +} + +.md-typeset .latex-algorithm__lines > li.latex-algorithm__indent-2 { + padding-left: 5.2rem; +} + +.md-typeset .latex-algorithm__lines > li.latex-algorithm__indent-2::after { + border-right: 0.05rem solid var(--md-default-fg-color--lightest); + width: 1.25rem; +} + +.md-typeset .latex-algorithm__keyword { + font-weight: 700; +} + +.md-typeset .latex-algorithm__operation { + background: color-mix( + in srgb, + var(--md-accent-fg-color) 7%, + var(--md-default-bg-color) + ); + border-left: 0.16rem solid var(--md-accent-fg-color); + border-radius: 0.2rem; + padding: 0.4rem 0.65rem; +} + +.md-typeset .latex-algorithm__operation--check { + background: color-mix(in srgb, #2e7d32 7%, var(--md-default-bg-color)); + border-left-color: #2e7d32; +} + +.md-typeset .latex-algorithm__operation--restart { + background: color-mix(in srgb, #ed6c02 7%, var(--md-default-bg-color)); + border-left-color: #ed6c02; +} + +.md-typeset .latex-algorithm__statement { + margin-top: 0.12rem; +} + +.md-typeset .latex-algorithm__statement--indent { + padding-left: 1.4rem; +} + +[data-md-color-scheme="slate"] + .md-typeset + .latex-algorithm__operation--check { + background: color-mix(in srgb, #66bb6a 10%, var(--md-default-bg-color)); + border-left-color: #66bb6a; +} + +[data-md-color-scheme="slate"] + .md-typeset + .latex-algorithm__operation--restart { + background: color-mix(in srgb, #ffb74d 10%, var(--md-default-bg-color)); + border-left-color: #ffb74d; +} + +.md-typeset .latex-algorithm .katex-display { + margin: 0.3rem 0; + overflow-x: auto; + overflow-y: hidden; + padding-block: 0.05rem; +} + +.md-typeset .latex-algorithm .katex-display > .katex { + font-size: 1em; +} + +@media screen and (max-width: 44.9844em) { + .md-typeset .latex-algorithm__lines > li { + padding-left: 2.35rem; + padding-right: 0.35rem; + } + + .md-typeset .latex-algorithm__lines > li.latex-algorithm__indent-1 { + padding-left: 3rem; + } + + .md-typeset .latex-algorithm__lines > li.latex-algorithm__indent-1::after, + .md-typeset .latex-algorithm__lines > li.latex-algorithm__indent-2::after { + left: 2.35rem; + } + + .md-typeset .latex-algorithm__lines > li.latex-algorithm__indent-2 { + padding-left: 3.75rem; + } + + .md-typeset .latex-algorithm__lines > li.latex-algorithm__indent-2::after { + width: 0.75rem; + } + + .md-typeset .latex-algorithm .katex-display > .katex { + font-size: 0.88em; + } +} diff --git a/pyproject.toml b/pyproject.toml index 62b0535..f60c034 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -5,7 +5,7 @@ build-backend = "scikit_build_core.build" [project] name = "cupdlpx" version = "0.3.0" -description = "Python bindings for cuPDLPx (GPU-accelerated first-order LP solver)" +description = "Python bindings for cuPDLPx, a GPU-accelerated first-order solver for large-scale linear programming" readme = "README.md" license = { text = "Apache-2.0" } requires-python = ">=3.8" diff --git a/python/README.md b/python/README.md deleted file mode 100644 index 7486b2b..0000000 --- a/python/README.md +++ /dev/null @@ -1,300 +0,0 @@ -# **Python Interface for cuPDLPx** - -[![License](https://img.shields.io/badge/License-Apache%202.0-blue.svg)](LICENSE) -[![GitHub release](https://img.shields.io/github/release/MIT-Lu-Lab/cuPDLPx.svg)](https://github.com/MIT-Lu-Lab/cuPDLPx/releases) -[![PyPI version](https://badge.fury.io/py/cupdlpx.svg)](https://pypi.org/project/cupdlpx/) -[![arXiv](https://img.shields.io/badge/arXiv-2407.16144-B31B1B.svg)](https://arxiv.org/abs/2407.16144) -[![arXiv](https://img.shields.io/badge/arXiv-2507.14051-B31B1B.svg)](https://arxiv.org/abs/2507.14051) - -There is the Python interface to **[`cuPDLPx`](../README.md)**, a GPU-accelerated first-order solver for large-scale linear programming (LP). -It provides a high-level, Pythonic API for constructing, modifying, and solving LPs using NumPy and SciPy data structures. - -## Installation - -### Requirements -- Python ≥ 3.9 -- NumPy ≥ 1.21 -- SciPy ≥ 1.8 -- An NVIDIA GPU with CUDA support (≥12.4 required) -- A C/C++ toolchain with GCC and NVCC - -> **SpMV backend** is selected automatically at compile time based on cuSPARSE version: -> - `cusparseSpMV` — CUDA 12.4 – 13.1 (cuSPARSE < 12.7.3) -> - `cusparseSpMVOp` — CUDA 13.1 Update 1+ (cuSPARSE ≥ 12.7.3) - - -### Install -Install from PyPI: - -```bash -pip install cupdlpx -``` - -Or build from source: - -```bash -git clone https://github.com/MIT-Lu-Lab/cuPDLPx.git -cd cuPDLPx -pip install . -``` - -## Quick Start - -```python -import numpy as np -from cupdlpx import Model, PDLP - -# Example: minimize c^T x -# subject to l <= A x <= u, lb <= x <= ub -c = np.array([1.0, 1.0]) -A = np.array([[1.0, 2.0], - [0.0, 1.0], - [3.0, 2.0]]) -l = np.array([5.0, -np.inf, -np.inf]) -u = np.array([5.0, 2.0, 8.0]) -lb = np.zeros(2) # x >= 0 -ub = None # no upper bound - -# Create LP model -m = Model(objective_vector=c, - constraint_matrix=A, - constraint_lower_bound=l, - constraint_upper_bound=u, - variable_lower_bound=lb, - variable_upper_bound=ub) - -# Set model sense -m.ModelSense = PDLP.MAXIMIZE - -# Parameters can be set in multiple ways -m.Params.TimeLimit = 60 # attribute style -m.setParam("FeasibilityTol", 1e-6) -m.setParams(OutputFlag=True, OptimalityTol=1e-8) - -# Solve -m.optimize() - -# Retrieve results -print("Status:", m.StatusName) -print("Objective:", m.ObjVal) -print("Primal solution:", m.X) -print("Dual solution:", m.Pi) -print("Reduced cost:", m.RC) -``` - -## Modeling - -The `Model` class represents a linear programming problem of the form: - -$$ -\min c^\top x + c_0 \quad -\text{s.t.} \; \ell \le A x \le u, \quad -\text{lb} \le x \le \text{ub}. -$$ - -### Arguments - -- **objective_vector** (`c`): Coefficients of the objective function. -- **constraint_matrix** (`A`): Coefficient matrix for the constraints. Both dense (`numpy.ndarray`) and sparse (`scipy.sparse.csr_matrix`) inputs are supported. Internally stored in double precision (`float64`). -- **constraint_lower_bound** (`l`): Lower bounds for each constraint. Use `-np.inf` or `None` for no lower bound. -- **constraint_upper_bound** (`u`): Upper bounds for each constraint. Use `+np.inf` or `None` for no upper bound. -- **variable_lower_bound** (`lb`, optional): Lower bounds for the decision variables. Defaults to `-np.inf` for all variables if not provided. -- **variable_upper_bound** (`ub`, optional): Upper bounds for the decision variables. Defaults to `+np.inf` for all variables if not provided. -- **objective_constant** (`c0`, optional): Constant offset in the objective function. Defaults to `0.0`. - -To initialize a linear programming problem, you need to provide the objective vector, constraint matrix, and bounds on both constraints and variables. - -```python -c = np.array([1.0, 1.0]) -A = np.array([[1.0, 2.0], - [0.0, 1.0], - [3.0, 2.0]]) -l = np.array([5.0, -np.inf, -np.inf]) -u = np.array([5.0, 2.0, 8.0]) -lb = np.zeros(2) # x >= 0 -ub = None # no upper bound - -# Create LP model -m = Model(objective_vector=c, - constraint_matrix=A, - constraint_lower_bound=l, - constraint_upper_bound=u, - variable_lower_bound=lb, - variable_upper_bound=ub) -``` - -### Reading from MPS Files - -A `Model` can also be created directly from an MPS file (plain or gzip-compressed) with `cupdlpx.read`: - -```python -import cupdlpx - -m = cupdlpx.read("problem.mps") # or "problem.mps.gz" -m.optimize() -``` - -The objective sense declared in the file (`OBJSENSE` section) is applied to `ModelSense` automatically. - - -## Model Sense - -By default, `cupdlpx` solves **minimization problems**. - -To switch between minimization and maximization, set the attribute `ModelSense`: - -```python -# Set model sense -m.ModelSense = PDLP.MAXIMIZE -``` - -## Parameters - -Solver parameters control termination criteria, logging, scaling, and restart behavior. - -Below is a list of commonly used parameters, their internal keys, and descriptions. - -| Alias | Internal Key | Type | Default | Description | -|---|---|---|---|---| -| `TimeLimit` | `time_sec_limit` | float | `3600.0` | Maximum wall-clock time in seconds. The solver terminates if the limit is reached. | -| `IterationLimit` | `iteration_limit` | int | `2147483647` | Maximum number of iterations. | -| `OutputFlag`, `LogToConsole` | `verbose` | bool | `False` | Enable (`True`) or disable (`False`) console logging output. | -| `TermCheckFreq` | `termination_evaluation_frequency` | int | `200` | Frequency (in iterations) at which termination conditions are evaluated. | -| `OptimalityNorm` | `optimality_norm` | string | `"l2"` | Norm for optimality criteria. Use `"l2"` for L2 norm or `"linf"` for infinity norm. | -| `OptimalityTol` | `eps_optimal_relative` | float | `1e-4` | Relative tolerance for optimality gap. Solver stops if the relative primal-dual gap ≤ this value. | -| `FeasibilityTol` | `eps_feasible_relative` | float | `1e-4` | Relative feasibility tolerance for primal/dual residuals. | -| `InfeasibleTol` | `eps_infeasible_relative` | float | `1e-14` | Relative tolerance on the ray certificate for declaring primal/dual infeasibility. | -| `GeoMeanIters` | `geometric_mean_iterations` | int | `12` | Number of iterations of geometric-mean scaling. Improves numerical conditioning. | -| `RuizIters` | `l_inf_ruiz_iterations` | int | `10` | Number of iterations of L∞ Ruiz scaling. Improves numerical conditioning. | -| `UsePCAlpha` | `has_pock_chambolle_alpha` | bool | `True` | Whether to use the Pock–Chambolle α step size adjustment. | -| `PCAlpha` | `pock_chambolle_alpha` | float | `1.0` | Value of the Pock–Chambolle α parameter. | -| `BoundObjRescaling` | `bound_objective_rescaling` | bool | `True` | Whether to rescale the objective vector during preprocessing. | -| `RestartArtificialThresh` | `artificial_restart_threshold` | float | `0.36` | Threshold for artificial restart. | -| `RestartSufficientReduction` | `sufficient_reduction_for_restart` | float | `0.2` | Sufficient reduction factor to justify a restart. | -| `RestartNecessaryReduction` | `necessary_reduction_for_restart` | float | `0.5` | Necessary reduction factor required for a restart. | -| `RestartKp` | `k_p` | float | `0.99` | Proportional coefficient for PID-controlled primal weight updates. | -| `ReflectionCoeff` | `reflection_coefficient` | float | `1.0` | Reflection coefficient. | -| `SVMaxIter` | `sv_max_iter` | int | 5000 | Maximum number of iterations for the power method | -| `SVTol`| `sv_tol` | float | `1e-4` | Termination tolerance for the power method | -| `Presolve`| `presolve` | bool | `True` | Whether to use presolve. | -| `FeasibilityPolishing` | `feasibility_polishing` | bool | `False` | Run feasibility polishing process.| -| `FeasibilityPolishingTol` | `eps_feas_polish_relative` | float | `1e-6` | Relative tolerance for primal/dual residual. | -| `Debug` | `debug` | bool | `False` | Developer diagnostics (implies `OutputFlag`). | -| `ActiveSetBoost` | `active_set_boost` | bool | `True` | Enable the active-set stepsize boost. | -| `ASBActivationTol` | `asb_activation_tol` | float | `1e-4` | Residual threshold at which the boost activates. | -| `ASBWindowIter` | `asb_window_iter` | int | `10000` | Number of recent iterations used to identify the active set. | -| `ASBSafetyFactor` | `asb_safety_factor` | float | `0.9` | Boosted step = factor / estimated singular value. | -| `ASBMaxReverts` | `asb_max_reverts` | int | `2` | Divergences tolerated before the boost turns off; a diverged step is always reverted. | -| `ASBMinRaiseRatio` | `asb_min_raise_ratio` | float | `1.1` | Minimum ratio for a step increase. | -| `ASBReestimateChangeRatio` | `asb_reestimate_change_ratio` | float | `0.01` | Fraction of the active set that must change before re-estimating. | -| `ASBConstraintTol` | `asb_constraint_tol` | float | `1e-8` | Tolerance for treating a constraint as binding. | -| `ASBVariableTol` | `asb_variable_tol` | float | `1e-8` | Tolerance for treating a variable as at its bound. | -| `ASBDivergenceCeilingRatio` | `asb_divergence_ceiling_ratio` | float | `0.7` | Step ceiling after a revert, relative to the diverged step. | -| `ASBDivergenceMargin` | `asb_divergence_margin` | float | `0.05` | Allowed fixed-point error increase before a revert. | - -They can be set in multiple ways: - -```python -# Method 1: single parameter -m.setParam("TimeLimit", 300) -m.setParam("FeasibilityTol", 1e-6) - -# Method 2: multiple parameters -m.setParams(TimeLimit=300, FeasibilityTol=1e-6) - -# Method 3: attribute-style access -m.Params.TimeLimit = 300 -m.Params.FeasibilityTol = 1e-6 - -# Reset all parameters to backend defaults -m.resetParams() -``` - -## Solution Attributes - -After calling `m.optimize()`, the solver stores results in a set of read-only attributes. `optimize()` returns the model itself, so chained access like `m.optimize().Status` is also supported. These attributes provide access to primal/dual solutions, objective values, residuals, and runtime statistics. - -### Attribute Reference - -| Attribute | Type | Description | -|---|---|---| -| `Status` | int | Integer termination status code; compare against `cupdlpx.PDLP` constants: `OPTIMAL=0`, `PRIMAL_INFEASIBLE=1`, `DUAL_INFEASIBLE=2`, `TIME_LIMIT=3`, `ITERATION_LIMIT=4`, `INFEASIBLE_OR_UNBOUNDED=5`, `FEAS_POLISH_SUCCESS=6`, `UNSPECIFIED=-1`. | -| `StatusName` | str | Human-readable status name, e.g. `"OPTIMAL"`, `"PRIMAL_INFEASIBLE"`. | -| `ObjVal` | float | Primal objective value at termination (sign-adjusted according to `ModelSense`). | -| `DualObj` | float | Dual objective value at termination. | -| `Gap` | float | Absolute primal-dual gap. | -| `RelGap` | float | Relative primal-dual gap. | -| `X` | numpy.ndarray | Primal solution vector \(x\). May be `None` if no feasible solution was found. | -| `Pi` | numpy.ndarray | Dual solution vector (Lagrange multipliers). | -| `RC` | numpy.ndarray | Reduced Cost vector. | -| `IterCount` | int | Number of iterations performed. | -| `Runtime` | float | Total wall-clock runtime in seconds. | -| `RescalingTime` | float | Time spent on preprocessing and rescaling (seconds). | -| `RelPrimalResidual` | float | Relative primal residual. | -| `RelDualResidual` | float | Relative dual residual. | -| `MaxPrimalRayInfeas` | float | Maximum primal ray infeasibility (indicator for infeasibility). | -| `MaxDualRayInfeas` | float | Maximum dual ray infeasibility. | -| `PrimalRayLinObj` | float | Linear objective value along a primal ray (used in infeasibility detection). | -| `DualRayObj` | float | Objective value along a dual ray (used in unboundedness detection). | - -All solution-related information can then be queried directly from the `Model` object: - -```python -m.optimize() - -print("Status:", m.StatusName, "(code:", m.Status, ")") -print("Primal objective:", m.ObjVal) -print("Dual objective:", m.DualObj) -print("Relative gap:", m.RelGap) -print("Iterations:", m.IterCount, " Runtime (s):", m.Runtime) - -# Access solutions -print("Primal solution:", m.X) -print("Dual solution:", m.Pi) -print("Reduced cost:", m.RC) - -# Check residuals -print("Primal residual:", m.RelPrimalResidual) -print("Dual residual:", m.RelDualResidual) -``` - -## Warm Start - -`cupdlpx` supports warm starting from user-provided primal and/or dual solutions. -This allows resuming from a previous iterate or reusing solutions from a related instance, often reducing iterations needed to reach optimality. - -```python -# Warm starting solution -x_init = [1.0, 2.0] -pi_init = [1.0, -1.0, 0.0] - -# Set warm start -m.setWarmStart(primal=x_init, dual=pi_init) - -# Solve -m.optimize() -``` - -Both primal and dual arguments are optional. You may specify only one of them if desired: - -```python -# Only provide primal start -m.setWarmStart(primal=x_init) - -# Only provide dual start -m.setWarmStart(dual=pi_init) -``` - -If the warm-start vectors have incorrect dimensions, `setWarmStart` raises a `ValueError`. Omitting an argument leaves that side unchanged; passing `None` clears it. - -To clear existing warm-start values: - -```python -m.clearWarmStart() -``` - -or - -```python -m.setWarmStart(primal=None, dual=None) -``` diff --git a/python/cupdlpx/__init__.py b/python/cupdlpx/__init__.py index fed1f83..874845a 100644 --- a/python/cupdlpx/__init__.py +++ b/python/cupdlpx/__init__.py @@ -12,7 +12,7 @@ # See the License for the specific language governing permissions and # limitations under the License. -"""cuPDLPx: Python bindings for the GPU-accelerated first-order LP solver.""" +"""Python bindings for cuPDLPx, a GPU-accelerated first-order solver for large-scale linear programming.""" import os import platform diff --git a/zensical.toml b/zensical.toml new file mode 100644 index 0000000..a4946de --- /dev/null +++ b/zensical.toml @@ -0,0 +1,143 @@ +[project] +site_name = "cuPDLPx" +site_url = "https://mit-lu-lab.github.io/cuPDLPx/" +site_description = "Documentation for cuPDLPx, an open-source, GPU-accelerated first-order solver for large-scale linear programming." +site_author = "Haihao Lu, Zedong Peng, Jinwen Yang, and contributors" +copyright = '© 2025–2026 MIT Lu Lab · Apache License 2.0
Website made with Zensical' +repo_url = "https://github.com/MIT-Lu-Lab/cuPDLPx" +repo_name = "MIT-Lu-Lab/cuPDLPx" +edit_uri = "edit/main/docs/" +docs_dir = "docs" +site_dir = "site" +extra_css = [ + "https://unpkg.com/katex@0.18.4/dist/katex.min.css", + "stylesheets/extra.css", +] +extra_javascript = [ + "javascripts/katex.js", + "https://unpkg.com/katex@0.18.4/dist/katex.min.js", + "https://unpkg.com/katex@0.18.4/dist/contrib/auto-render.min.js", +] + +nav = [ + { "Home" = "index.md" }, + { "Get started" = [ + { "Overview" = "getting-started/index.md" }, + { "Python interface" = "guides/python.md" }, + { "Julia interface" = "guides/julia.md" }, + { "Command-line interface" = "guides/command-line.md" }, + { "C interface" = "guides/c-api.md" }, + { "Log interpretation" = "getting-started/log-interpretation.md" }, + { "Results and status" = "getting-started/results.md" }, + ] }, + { "Algorithm" = [ + { "Overview" = "algorithm/index.md" }, + { "Base algorithm" = "algorithm/base-algorithm.md" }, + { "Presolve and postsolve" = "algorithm/presolve.md" }, + { "Preconditioning" = "algorithm/preconditioning.md" }, + { "Step size" = "algorithm/step-size.md" }, + { "Adaptive restart" = "algorithm/restart.md" }, + { "Primal weight" = "algorithm/primal-weight.md" }, + { "Termination criteria" = "algorithm/termination.md" }, + { "Feasibility polishing" = "algorithm/feasibility-polishing.md" }, + ] }, + { "Implementation" = [ + { "Overview" = "implementation/index.md" }, + { "GPU backends" = "implementation/gpu-backends.md" }, + { "SpMV" = "implementation/sparse-matrix-vector-products.md" }, + { "Kernel fusion" = "implementation/kernel-fusion.md" }, + { "CUDA Graphs" = "implementation/cuda-graphs.md" }, + ] }, + { "Parameters" = "reference/parameters.md" }, + { "Citation" = "citation.md" }, +] + +[project.theme] +font = false +variant = "modern" +language = "en" +favicon = "assets/mit-favicon.svg" +features = [ + "content.action.edit", + "content.action.view", + "content.code.copy", + "content.code.select", + "content.tabs.link", + "content.tooltips", + "navigation.footer", + "navigation.indexes", + "navigation.instant", + "navigation.instant.prefetch", + "navigation.path", + "navigation.sections", + "navigation.tabs", + "navigation.tabs.sticky", + "navigation.top", + "navigation.tracking", + "search.highlight", + "toc.follow", +] + +[project.theme.icon] +repo = "fontawesome/brands/github" +edit = "material/pencil" +view = "material/eye" + +[project.extra] +generator = false + +[[project.theme.palette]] +media = "(prefers-color-scheme)" +toggle.icon = "lucide/sun-moon" +toggle.name = "Use light theme" + +[[project.theme.palette]] +media = "(prefers-color-scheme: light)" +scheme = "default" +toggle.icon = "lucide/sun" +toggle.name = "Use dark theme" + +[[project.theme.palette]] +media = "(prefers-color-scheme: dark)" +scheme = "slate" +toggle.icon = "lucide/moon" +toggle.name = "Use system theme" + +[[project.extra.social]] +icon = "fontawesome/brands/github" +link = "https://github.com/MIT-Lu-Lab/cuPDLPx" +name = "cuPDLPx on GitHub" + +[[project.extra.social]] +icon = "fontawesome/solid/cube" +link = "https://pypi.org/project/cupdlpx/" +name = "cuPDLPx on PyPI" + +[project.validation] +invalid_links = true +invalid_link_anchors = true +unresolved_references = true +unresolved_footnotes = true + +[project.markdown_extensions] +abbr = {} +admonition = {} +attr_list = {} +def_list = {} +footnotes = {} +md_in_html = {} +toc.permalink = true +pymdownx.arithmatex.generic = true +pymdownx.details = {} +pymdownx.highlight.anchor_linenums = true +pymdownx.highlight.line_spans = "__span" +pymdownx.highlight.pygments_lang_class = true +pymdownx.inlinehilite = {} +pymdownx.superfences.custom_fences = [ + { name = "mermaid", class = "mermaid", format = "pymdownx.superfences.fence_code_format" }, +] +pymdownx.tabbed.alternate_style = true +pymdownx.tabbed.combine_header_slug = true +pymdownx.tasklist.custom_checkbox = true + +[project.plugins.offline]