GPU implementation of Clarabel solver for Julia

Features • Installation • License • Documentation

CuClarabel.jl is the GPU implementation of the Clarabel solver, which can solve conic problems of the following form:

$$ \begin{array}{r} \text{minimize} & \frac{1}{2}x^T P x + q^T x\\[2ex] \text{subject to} & Ax + s = b \\[1ex] & s \in \mathcal{K} \end{array} $$

with decision variables $x \in \mathbb{R}^n$, $s \in \mathbb{R}^m$ and data matrices $P=P^\top \succeq 0$, $q \in \mathbb{R}^n$, $A \in \mathbb{R}^{m \times n}$, and $b \in \mathbb{R}^m$. The set $\mathcal{K}$ is a composition of convex cones; we support zero cones (linear equality constraints), nonnegative cones (linear inequality constraints), second-order cones, exponential cone, power cones and semidefinite cones of the same dimensionality. Our package relies on the external package CUDSS.jl for the linear system solver CUDSS. We also support linear system solves in lower (mixed) precision.

Installation

The project has been merged back to the original repo of Clarabel.jl under the branch CuClarabel. The repo CuClarabel.jl will no longer be updated.

Please go to https://github.com/oxfordcontrol/Clarabel.jl/tree/CuClarabel for the latest update of the GPU solver.

Tutorial

Use in Julia

Modeling a conic optimization problem is the same as in the original Clarabel solver, except with the additional parameter direct_solve_method. This can be set to :cudss or :cudssmixed. Here is a portfolio optimization problem modelled via JuMP: ``` using LinearAlgebra, SparseArrays, Random, JuMP using Clarabel

generate the data

rng = Random.MersenneTwister(1) k = 5; # number of factors n = k * 10; # number of assets D = spdiagm(0 => rand(rng, n) .* sqrt(k)) F = sprandn(rng, n, k, 0.5); # factor loading matrix μ = (3 .+ 9. * rand(rng, n)) / 100. # expected returns between 3% - 12% γ = 1.0; # risk aversion parameter d = 1 # we are starting from all cash x0 = zeros(n);

a = 1e-3 b = 1e-1 γ = 1.0;

model = JuMP.Model(Clarabel.Optimizer) setoptimizerattribute(model, "directsolvemethod", :cudss)

@variable(model, x[1:n]) @variable(model, y[1:k])
@variable(model, s[1:n]) @variable(model, t[1:n]) @objective(model, Min, x' * D * x + y' * y - 1/γ * μ' * x); @constraint(model, y .== F' * x); @constraint(model, x .>= 0);

transaction costs

@constraint(model, sum(x) + a * sum(s) == d + sum(x0) ); @constraint(model, [i = 1:n], x0[i] - x[i] == t[i]) @constraint(model, [i = 1:n], [s[i], t[i]] in MOI.SecondOrderCone(2)); JuMP.optimize!(model) ```

Use in Python

We can call julia code within a python file by using JuliaCall package. We can download the package by pip install juliacall Then, we load the package as a single variable jl which represents the Main module in Julia, and we can write Julia code and call it via jl.seval() in Python. ``` from juliacall import Main as jl import numpy as np import cupy as cp from cupyx.scipy.sparse import csr_matrix

Load Clarabel in Julia

jl.seval('using Clarabel, LinearAlgebra, SparseArrays') jl.seval('using CUDA, CUDA.CUSPARSE') Here we build up a simple optimization problem with a second-order cone, which is fully written by Julia. jl.seval(''' P = CuSparseMatrixCSR(sparse([2.0 1.0 0.0; 1.0 2.0 0.0; 0.0 0.0 2.0])) q = CuVector([0, -1., -1]) A = CuSparseMatrixCSR(SparseMatrixCSC([1. 0 0; -1 0 0; 0 -1 0; 0 0 -1])) b = CuVector([1, 0., 0., 0.])

# 0-cone dimension 1, one second-order-cone of dimension 3
cones = Dict("f" => 1, "q"=> [3])

settings = Clarabel.Settings(direct_solve_method = :cudss)

solver   = Clarabel.Solver(P,q,A,b,cones, settings)
Clarabel.solve!(solver)

# Extract solution
x = solver.solution

''') It is also possible to call the julia functions directly via JuliaCall. For example, if we want to reuse the solver object and update only coefficients in the problem, we can call the following blocks,

Update b vector

bpy = cp.array([2.0, 1.0, 1.0, 1.0], dtype=cp.float64) bjl = jl.Clarabel.cupytocuvector(jl.Float64, int(bpy.data.ptr), bpy.size)

"_b" is the replacement of "!" in julia function

jl.Clarabel.updatebb(jl.solver,bjl) #Clarabel.update_b!()

Update P matrix

Define a new CSR sparse matrix on GPU

Ppy = csr_matrix(cp.array([ [3.0, 0.5, 0.0], [0.5, 2.0, 0.0], [0.0, 0.0, 1.0] ], dtype=cp.float64))

Extract the pointers (as integers)

dataptr = int(Ppy.data.data.ptr) indicesptr = int(Ppy.indices.data.ptr) indptr_ptr = int(Ppy.indptr.data.ptr)

Get matrix shape and non-zero count

nrows, ncols = Ppy.shape nnz = Ppy.nnz

jl.Pjl = jl.Clarabel.cupytocucsrmat(jl.Float64, dataptr, indicesptr, indptrptr, nrows, n_cols, nnz)

jl.Clarabel.updatePb(jl.solver, jl.Pjl) #Clarabel.update_P!()

Solve the new problem without creating memory

jl.Clarabel.solveb(jl.solver) #Clarabel.solve!() `` where we update the linear costbby values in a cupy vectorbpyand the quadratic costPby values in a cupy csr matrixPpy. Note that we need to replace!in a julia function withb. Reversely, we can also extract value from a Julia object back to Python, ``

Retrieve the solution from Julia to Python

solution = np.array(jl.solver.solution.x) print("Solution:", solution) `` The example file can be found under thepython` folder.

Performance tips

Due to the just-in-time (JIT) compilation in Julia, the first call of CuClarabel will also be slow in python and it is recommended to solve a mini problem first to trigger the JIT-compilation and get full performance on the subsequent solve of the actual problem.

Citing

@misc{CuClarabel, title={CuClarabel: GPU Acceleration for a Conic Optimization Solver}, author={Yuwen Chen and Danny Tse and Parth Nobel and Paul Goulart and Stephen Boyd}, year={2024}, eprint={2412.19027}, archivePrefix={arXiv}, primaryClass={math.OC}, url={https://arxiv.org/abs/2412.19027}, }

License 🔍

This project is licensed under the Apache License 2.0 - see the LICENSE.md file for details.