oxfordcontrol / clarabel.rs Goto Github PK

Clarabel.rs: Interior-point solver for convex conic optimisation problems in Rust.

License: Apache License 2.0

Rust 97.72% HTML 0.10% Julia 2.16% Python 0.02%

conic-optimization conic-programs convex-optimization interior-point-method linear-programming optimization optimization-algorithms quadratic-programming rust-lang semidefinite-programming

clarabel.rs's Introduction

Interior Point Conic Optimization for Rust and Python

Features • Installation • License • Documentation

Clarabel.rs is a Rust implementation of an interior point numerical solver for convex optimization problems using a novel homogeneous embedding. Clarabel.rs solves the following problem:

$$ \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 convex set $\mathcal{K}$ is a composition of convex cones.

For more information see the Clarabel Documentation (stable | dev).

Clarabel is also available in a Julia implementation. See here.

Features

Versatile: Clarabel.rs solves linear programs (LPs), quadratic programs (QPs), second-order cone programs (SOCPs) and semidefinite programs (SDPs). It also solves problems with exponential, power cone and generalized power cone constraints.
Quadratic objectives: Unlike interior point solvers based on the standard homogeneous self-dual embedding (HSDE), Clarabel.rs handles quadratic objectives without requiring any epigraphical reformulation of the objective. It can therefore be significantly faster than other HSDE-based solvers for problems with quadratic objective functions.
Infeasibility detection: Infeasible problems are detected using a homogeneous embedding technique.
Open Source: Our code is available on GitHub and distributed under the Apache 2.0 License

Installation

Clarabel can be imported to Cargo based Rust projects by adding

[dependencies]
clarabel = "0"

to the project's Cargo.toml file. To install from source, see the Rust Installation Documentation.

To use the Python interface to the solver:

pip install clarabel

To install the Python interface from source, see the Python Installation Documentation.

Citing

@misc{Clarabel_2024,
      title={Clarabel: An interior-point solver for conic programs with quadratic objectives}, 
      author={Paul J. Goulart and Yuwen Chen},
      year={2024},
      eprint={2405.12762},
      archivePrefix={arXiv},
      primaryClass={math.OC}
}

License 🔍

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

clarabel.rs's People

Contributors

Stargazers

Watchers

clarabel.rs's Issues

Clarabel 7.1 errors with insufficient progress on a seemingly simple SOCP

While upgrading Drake to clarabel 7.1 (RobotLocomotion/drake#21053), one of our (seemingly simple) unit tests started falling down with "insufficient progress". Mosek solves it fine.

I've carved out the following reproduction

import clarabel
import numpy as np
from scipy import sparse

q = [0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1]
b = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, -1, 1, -0, -0, -0, 1, -0, -0, -0, 1, 1, -0, 1, -0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
P = sparse.csc_matrix((18, 18))

data = [-1, 1, 1, -1, 1, -1, -1, -1, -1, 1, -0.5, -0.5, 1, -0.7071067811865476, 1, -0.7071067811865476, 1, -1, 0.7071067811865475, -1.0536712127723509e-08, 1, -1, 0.7071067811865475, -1.0536712127723509e-08, -0.5, 0.5, -1, -1, 1, -1, 1, -1, -1, -1, -1, 1, -0.5, -0.5, 1, -0.7071067811865476, 1, -0.7071067811865476, 1, -1, 0.7071067811865475, -1.0536712127723509e-08, 1, -1, 0.7071067811865475, -1.0536712127723509e-08, -0.5, 0.5, -1, 1, -2, 1, -1, 1, 1, -1, 1, -0.5, -0.5, 1, 1, -0.7071067811865476, 1, 1, -0.7071067811865476, 1, 0.7071067811865475, -1.0536712127723509e-08, 1, 0.7071067811865475, -1.0536712127723509e-08, -0.5, 0.5]
rows = [2, 3, 12, 13, 20, 21, 22, 23, 30, 31, 34, 35, 0, 36, 1, 37, 2, 14, 36, 38, 3, 15, 37, 39, 34, 35, 6, 7, 12, 16, 24, 25, 26, 27, 30, 31, 40, 41, 4, 42, 5, 43, 6, 17, 42, 44, 7, 18, 43, 45, 40, 41, 8, 9, 10, 13, 19, 28, 29, 32, 33, 46, 47, 8, 14, 48, 9, 15, 49, 10, 48, 50, 11, 49, 51, 46, 47]
cols = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 8, 8, 9, 9, 9, 9, 10, 10, 10, 10, 11, 11, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 13, 13, 13, 14, 14, 14, 15, 15, 15, 16, 16, 16, 17, 17]
A = sparse.csc_matrix((data, (rows, cols)))
cones = [
  clarabel.ZeroConeT(20),
  clarabel.NonnegativeConeT(1),
  clarabel.NonnegativeConeT(1),
  clarabel.NonnegativeConeT(1),
  clarabel.NonnegativeConeT(1),
  clarabel.NonnegativeConeT(1),
  clarabel.NonnegativeConeT(1),
  clarabel.NonnegativeConeT(1),
  clarabel.NonnegativeConeT(1),
  clarabel.NonnegativeConeT(1),
  clarabel.ZeroConeT(1),
  clarabel.NonnegativeConeT(4),
  clarabel.SecondOrderConeT(6),
  clarabel.SecondOrderConeT(6),
  clarabel.SecondOrderConeT(6),
]

settings = clarabel.DefaultSettings()
solver = clarabel.DefaultSolver(P, q, A, b, cones, settings)
solver.solve()

Results in

-------------------------------------------------------------
           Clarabel.rs v0.7.1  -  Clever Acronym              

                   (c) Paul Goulart                          
                University of Oxford, 2022                   
-------------------------------------------------------------

problem:
  variables     = 18
  constraints   = 52
  nnz(P)        = 0
  nnz(A)        = 77
  cones (total) = 15
    :        Zero = 2,  numel = (20,1)
    : Nonnegative = 10,  numel = (1,1,1,1,...,4)
    : SecondOrder = 3,  numel = (6,6,6)

settings:
  linear algebra: direct / qdldl, precision: 64 bit
  max iter = 200, time limit = Inf,  max step = 0.990
  tol_feas = 1.0e-8, tol_gap_abs = 1.0e-8, tol_gap_rel = 1.0e-8,
  static reg : on, ϵ1 = 1.0e-8, ϵ2 = 4.9e-32
  dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-7
  iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
               max iter = 10, stop ratio = 5.0
  equilibrate: on, min_scale = 1.0e-5, max_scale = 1.0e5
               max iter = 10

iter    pcost        dcost       gap       pres      dres      k/t        μ       step      
---------------------------------------------------------------------------------------------
  0  +0.0000e+00  -5.7518e+00  5.75e+00  9.93e-01  6.89e+03  1.00e+00  1.33e+03   ------   
  1  +4.4363e-02  +7.9443e+01  7.94e+01  9.93e-01  3.89e-01  8.61e+01  2.96e+02  7.90e-01  
  2  +4.4016e-01  +4.4876e+02  4.48e+02  9.92e-01  4.14e-01  4.55e+02  2.93e+02  8.17e-02  
  3  +2.1750e+00  +1.5761e+03  7.24e+02  9.91e-01  1.30e+00  1.58e+03  2.79e+02  3.00e-01  
  4  +5.5621e+00  +2.2199e+03  3.98e+02  9.69e-01  1.81e-01  2.22e+03  4.48e+01  8.55e-01  
  5  +6.7417e+00  +1.0322e+02  1.43e+01  4.21e-01  4.90e-03  9.66e+01  1.22e+00  9.83e-01  
  6  +5.0782e+00  +6.1775e+00  2.16e-01  4.76e-02  5.32e-04  1.10e+00  1.31e-01  9.05e-01  
  7  +5.0782e+00  +6.1775e+00  2.16e-01  4.76e-02  5.32e-04  1.10e+00  1.31e-01  0.00e+00  
---------------------------------------------------------------------------------------------
Terminated with status = InsufficientProgress
solve time = 286.812µs

Admittedly there are some ugly small coefficients in there... this was generated through our higher-level GCS abstraction. But I thought it was best to capture the instance and report.

Julia vs Rust performance?

Hello and thank you for sharing this software with the world.

This is the only project I've enountered so far written in both Julia and Rust.

I'm curious whether you've compared the performance between the two implementations. Theoretically, they should be quite similar after Julia's precompilation, but I'm curious to hear what you've observed in practice.

Thanks,
Oliver

Normed constraint, stupid question

First of all many thanks for the outstanding library.

I have a short mathematical question concerning the mathematical problem with soc.

I would like to understand the mathematical translation in terms of second order cone for the following problem

$$ \inf x^\top P x + q^\top x $$

under the constraint

$$ ||x||_1 =\sum |x_k| \leq 1 $$

It is clear that it can be transformed with $n$ second order constraints $||x_k||_2 \leq t_k$ and $\sum t_k \leq 1$. However I don't get how to define this auxiliary variable in terms of the matrix $A$ and vector $b$ as well as the vector of cones.

Many thanks in advance for answering this trivial question.

Support devcontainers

Could help to simplify dive into Rust

Assertion failed

I started using clarabel in a project and managed to make the following assertion fail:

Clarabel.rs/src/qdldl/qdldl.rs

Line 741 in d1ca25d

debug_assert!(*p.iter().max().unwrap_or(&0) < x.len());

Unfortunately I don't have a test case to reproduce it, as I am generating the LP problems dynamically and they are massive.

Question: Will Clarabel support box cone?

Hello there,

I am starting to benchmark Clarabel, and would like to know if support for box cone (to manage box constraints on variables more "natively" than with the positive cone) is something planned?

(Like SCS -> https://www.cvxgrp.org/scs/api/cones.html)

Cheers,

Roman

textbook sos example that worked on v0.7.1 and fails on v0.9.0

We recently upgraded Drake from v0.7.1 to v0.9.0, and I experienced a regression in one of my "simple" examples using sums-of-squares for Lyapunov analysis. Mosek, SCS, and Clarabel v0.7.1 solved it fine.

For context, you can find the original notebook here. Unfortunately, the failure happens inside some alternations... specifically in the OptimizeLyapunov method at the bottom of the notebook, and if OptimizeMultipliers uses Mosek, and only OptimizeLyapunov uses Clarabel, then everything works.

Nevertheless, here is the Clarabel-only reproduction of the final solve which now results in Terminated with status = InsufficientProgress.

import clarabel
import numpy as np
from scipy import sparse

q = [1.8435820570082377, 1.8435820570082377, 2.7653730855123553, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
b = [-1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
P = sparse.csc_matrix((81, 81))

data = [-4.6713756474581505, -7.366046824401246e-16, -0.33482275518489235, 1.248525969168963e-19, 2.6328516974490616e-09, -3.9681215422316167e-16, -0.9587196352099683, -2.8174889909481824e-16, 2.6184498937672173e-07, -0.10625452869742728, 1.4101154172740503e-16, 0.5484683093767915, -1, 2.031287962092983e-21, -0.2742382564473529, -2.5439681835905437e-10, -1, -4.6713756474581505, -7.366046824401246e-16, -0.33482275518489235, 1.248525969168963e-19, 2.6328516974490616e-09, -4.6713756474581505, -1.1334168366632864e-15, -2.293542390394861, -2.8162404649790135e-16, 2.644778410741708e-07, 4.6713756474581505, 3.397925282169629e-16, 0.2698485912774968, -1.408622099643301e-16, 0.5484685685889291, 3.9681215422316167e-16, 2.523840753970691, 1.1593671545503101e-15, 0.6090525462693416, -1.248525969168963e-19, -2.6328516974490616e-09, 0.10625452869742728, 2.5580264378371873e-16, -0.8639869308685729, 2.8174889909481824e-16, -2.6184498937672173e-07, -2.031287962092983e-21, 0.38049278489038335, -1.4101154172740503e-16, -0.5484683093767915, 2.5439681835905437e-10, -2.031287962092983e-21, 0.2742382564473529, 2.5439681835905437e-10, -1.4142135623730951, -4.6713756474581505, -7.366046824401246e-16, -1.3348227551848924, 1.248525969168963e-19, 2.6328516974490616e-09, 4.6713756474581505, 3.397925282169629e-16, 0.37610311997492407, -2.8187375169173514e-16, 2.5921213767927265e-07, 3.9681215422316167e-16, 3.523840753970692, 1.159365123262348e-15, 0.8832908027166945, -1.248525969168963e-19, -2.6328516974490616e-09, 0.10625452869742728, 2.5580264378371873e-16, -0.8639869306141761, 2.8174889909481824e-16, -2.6184498937672173e-07, -2.031287962092983e-21, 0.38049278489038335, -1.4101154172740503e-16, -0.5484683093767915, 2.5439681835905437e-10, -2.031287962092983e-21, 0.2742382564473529, 2.5439681835905437e-10, -1, 1, -1, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 1, -1, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 1, -1, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 1, -1, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 1, -1, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 1, -1, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 1, -1, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 1, -1, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 1, -1, 2, -1.4142135623730951, 2, -1.4142135623730951, 2, -1.4142135623730951, 1, -1, 2, -1.4142135623730951, 2, -1.4142135623730951, 1, -1, 2, -1.4142135623730951, 1, -1]
rows = [5, 6, 7, 8, 9, 12, 13, 14, 15, 19, 20, 21, 24, 25, 26, 30, 115, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 11, 12, 13, 14, 15, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 116, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 12, 13, 14, 15, 16, 17, 19, 20, 21, 22, 23, 25, 26, 27, 28, 30, 31, 32, 35, 117, 17, 37, 10, 38, 23, 40, 16, 43, 9, 47, 28, 52, 22, 58, 15, 65, 8, 73, 27, 82, 21, 92, 14, 103, 4, 39, 16, 41, 9, 44, 3, 48, 22, 53, 15, 59, 8, 66, 2, 74, 21, 83, 14, 93, 7, 104, 28, 42, 22, 45, 15, 49, 32, 54, 27, 60, 21, 67, 14, 75, 31, 84, 26, 94, 20, 105, 15, 46, 8, 50, 27, 55, 21, 61, 14, 68, 7, 76, 26, 85, 20, 95, 13, 106, 2, 51, 21, 56, 14, 62, 7, 69, 1, 77, 20, 86, 13, 96, 6, 107, 35, 57, 31, 63, 26, 70, 20, 78, 34, 87, 30, 97, 25, 108, 26, 64, 20, 71, 13, 79, 30, 88, 25, 98, 19, 109, 13, 72, 6, 80, 25, 89, 19, 99, 12, 110, 0, 81, 19, 90, 12, 100, 5, 111, 33, 91, 29, 101, 24, 112, 24, 102, 18, 113, 11, 114]
cols = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15, 16, 16, 17, 17, 18, 18, 19, 19, 20, 20, 21, 21, 22, 22, 23, 23, 24, 24, 25, 25, 26, 26, 27, 27, 28, 28, 29, 29, 30, 30, 31, 31, 32, 32, 33, 33, 34, 34, 35, 35, 36, 36, 37, 37, 38, 38, 39, 39, 40, 40, 41, 41, 42, 42, 43, 43, 44, 44, 45, 45, 46, 46, 47, 47, 48, 48, 49, 49, 50, 50, 51, 51, 52, 52, 53, 53, 54, 54, 55, 55, 56, 56, 57, 57, 58, 58, 59, 59, 60, 60, 61, 61, 62, 62, 63, 63, 64, 64, 65, 65, 66, 66, 67, 67, 68, 68, 69, 69, 70, 70, 71, 71, 72, 72, 73, 73, 74, 74, 75, 75, 76, 76, 77, 77, 78, 78, 79, 79, 80, 80]
A = sparse.csc_matrix((data, (rows, cols)))
cones = [
  clarabel.ZeroConeT(37),
  clarabel.NonnegativeConeT(0),
  clarabel.PSDTriangleConeT(12),
  clarabel.PSDTriangleConeT(2),
]

settings = clarabel.DefaultSettings()
solver = clarabel.DefaultSolver(P, q, A, b, cones, settings)
solver.solve()

I know at least one relatively simple thing I can do to improve the numerics for this problem, but thought I should report the regression before I do.

Clarabel solver is occupying all cores

Hello,

Thanks for sharing this beautiful software.

I've encountered that for large problems, Clarabel tends to occupy all threads/cores. How can we control the paralellization of the solver? Are there some environment variables which we can set to control this behaviour?

I observed this happens especially when the P matrix has a big block of dense data (filled in matrix); in our case, the P matrix looks like:
P = [ [diag, O, O O], [O, S, O, O], [O, O, diag, O], [O, O, O, diag]], where diag are diagonal matrices, O is matrix of zeros, and S is a matrix filled with non-zero elements, each around 200x200 dimensions.

Many thanks in advance!

Python interface does not print output in colab

The solver does not produce output to the terminal when used via the python interface within google colab.

Example here: https://colab.research.google.com/drive/1HgEPqJAk__I6QKRkLz7jynTy07YEyKmD#scrollTo=KsnL1mMgOBIf

Checking results status in Python

Thanks for releasing Clarabel! As I was trying it out in Python, I noticed something confusing about the status of a result: currently it can be both == and != to Solved at the same time.

Here is a small example, following: https://oxfordcontrol.github.io/ClarabelDocs/stable/python/getting_started_py/

    solver = clarabel.DefaultSolver(P, q, A, b, cones, settings)
    result = solver.solve()
    print(type(result.status))
    print(result.status == clarabel.SolverStatus.Solved)
    print(result.status != clarabel.SolverStatus.Solved)
    print(result.status is clarabel.SolverStatus.Solved)  # to check the underlying references

The output on my machine is:

<class 'builtins.SolverStatus'>
True
True
False

Clarabel was installed from PyPI. Version:

>>> import clarabel
>>> clarabel.__version__
<built-in function __version__>  # a bit confusing, __version__ is usually a string in Python
>>> clarabel.__version__()
'0.3.0'

ClarabelRs with `faer` result in `NumericalError`. Saving to `json` and reloading fixes issue?

Hi there,

I am slightly puzzled by this. I am reading a model from a cbf using JuMP and setting the direct method to :faer. Unfortunately something happens and the solver bails with a NumericalError.

using JuMP, ClarabelRs
path_to_file = joinpath(@__DIR__,path_to_files, file_name)
model_faer = read_from_file(path_to_file)
set_optimizer(model_faer, ClarabelRs.Optimizer)
set_attribute(model_faer, "direct_solve_method", :faer)
optimize!(model_faer)

-------------------------------------------------------------
           Clarabel.rs v0.9.0  -  Clever Acronym

                   (c) Paul Goulart
                University of Oxford, 2022
-------------------------------------------------------------

problem:
  variables     = 86778
  constraints   = 112351
  nnz(P)        = 0
  nnz(A)        = 580925
  cones (total) = 8251
    :        Zero = 1,  numel = 58738
    : Nonnegative = 2,  numel = (1,4124)
    : PSDTriangle = 8248,  numel = (6,6,6,6,...,6)

settings:
  linear algebra: direct / faer, precision: 64 bit
  max iter = 200, time limit = Inf,  max step = 0.990
  tol_feas = 1.0e-8, tol_gap_abs = 1.0e-8, tol_gap_rel = 1.0e-8,
  static reg : on, ϵ1 = 1.0e-8, ϵ2 = 4.9e-32
  dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-7
  iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
               max iter = 10, stop ratio = 5.0
  equilibrate: on, min_scale = 1.0e-4, max_scale = 1.0e4
               max iter = 10

iter    pcost        dcost       gap       pres      dres      k/t        μ       step
---------------------------------------------------------------------------------------------
  0  +9.2749e-01  +9.2749e-01  1.99e-13  1.34e-01  2.21e-02  1.00e+00  5.72e+01   ------
  1  -2.6265e+01  -2.0492e+01  2.82e-01  6.95e-02  1.10e-02  6.28e+00  3.43e+01  5.22e-01
  2  -5.2928e+04  -5.2918e+04  1.87e-04  1.62e-02  3.85e-03  1.05e+01  3.42e+01  5.27e-02
  3  -2.0302e+05  -2.0300e+05  9.90e-05  5.32e-03  1.29e-03  2.09e+01  2.91e+01  2.45e-01
  4  -3.0882e+05  -3.0879e+05  9.61e-05  3.69e-03  8.97e-04  3.05e+01  2.55e+01  2.20e-01
  5  -5.1792e+05  -5.1788e+05  8.41e-05  2.22e-03  5.39e-04  4.44e+01  1.98e+01  3.30e-01
  6  -7.7519e+05  -7.7514e+05  7.13e-05  1.37e-03  3.32e-04  5.60e+01  1.45e+01  3.72e-01
  7  -8.6242e+05  -8.6237e+05  6.44e-05  1.10e-03  2.66e-04  5.62e+01  1.21e+01  2.38e-01
  8  -8.7572e+05  -8.7566e+05  5.90e-05  9.33e-04  2.26e-04  5.23e+01  1.03e+01  2.35e-01
  9  -8.7381e+05  -8.7377e+05  5.12e-05  7.88e-04  1.91e-04  4.52e+01  8.64e+00  2.13e-01
 10  -8.0331e+05  -8.0328e+05  3.94e-05  5.53e-04  1.34e-04  3.20e+01  5.77e+00  4.97e-01
 11  -7.6529e+05  -7.6526e+05  3.36e-05  4.20e-04  1.01e-04  2.59e+01  4.25e+00  4.44e-01
 12  -7.4962e+05  -7.4961e+05  2.49e-05  2.95e-04  7.11e-05  1.88e+01  2.94e+00  4.23e-01
 13  -7.5090e+05  -7.5089e+05  2.17e-05  2.05e-04  4.95e-05  1.64e+01  2.03e+00  6.92e-01
 14  -7.4220e+05  -7.4219e+05  1.34e-05  1.17e-04  2.83e-05  9.98e+00  1.15e+00  6.79e-01
 15  -7.3116e+05  -7.3116e+05  7.48e-06  6.46e-05  1.55e-05  5.50e+00  6.25e-01  5.42e-01
 16  -7.2435e+05  -7.2434e+05  3.51e-06  2.90e-05  6.97e-06  2.55e+00  2.78e-01  7.46e-01
 17  -7.2435e+05  -7.2434e+05  3.51e-06  2.90e-05  6.97e-06  2.55e+00  2.78e-01  0.00e+00
---------------------------------------------------------------------------------------------
Terminated with status = NumericalError
solve time = 4.263004331s

What I find puzzling is that if I save the model to a json file and then load the model from said json file it seems as no NumericalError is encountered even though the initial iterations are exactly the same?

json_file = joinpath(@__DIR__,path_to_files, file_name)[1:end-4] * ".json"
ClarabelRs.write_to_file(unsafe_backend(model_faer).solver,json_file)
json_solver = ClarabelRs.read_from_file(json_file)
ClarabelRs.solve!(json_solver)

-------------------------------------------------------------
           Clarabel.rs v0.9.0  -  Clever Acronym

                   (c) Paul Goulart
                University of Oxford, 2022
-------------------------------------------------------------

problem:
  variables     = 86778
  constraints   = 112351
  nnz(P)        = 0
  nnz(A)        = 580925
  cones (total) = 8251
    :        Zero = 1,  numel = 58738
    : Nonnegative = 2,  numel = (1,4124)
    : PSDTriangle = 8248,  numel = (6,6,6,6,...,6)

settings:
  linear algebra: direct / faer, precision: 64 bit
  max iter = 200, time limit = Inf,  max step = 0.990
  tol_feas = 1.0e-8, tol_gap_abs = 1.0e-8, tol_gap_rel = 1.0e-8,
  static reg : on, ϵ1 = 1.0e-8, ϵ2 = 4.9e-32
  dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-7
  iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
               max iter = 10, stop ratio = 5.0
  equilibrate: on, min_scale = 1.0e-4, max_scale = 1.0e4
               max iter = 10

iter    pcost        dcost       gap       pres      dres      k/t        μ       step
---------------------------------------------------------------------------------------------
  0  +9.2749e-01  +9.2749e-01  2.02e-13  1.34e-01  2.21e-02  1.00e+00  5.72e+01   ------
  1  -2.6265e+01  -2.0492e+01  2.82e-01  6.95e-02  1.10e-02  6.28e+00  3.43e+01  5.22e-01
  2  -5.2928e+04  -5.2918e+04  1.87e-04  1.62e-02  3.85e-03  1.05e+01  3.42e+01  5.27e-02
  3  -2.0302e+05  -2.0300e+05  9.90e-05  5.32e-03  1.29e-03  2.09e+01  2.91e+01  2.45e-01
  4  -3.0882e+05  -3.0879e+05  9.61e-05  3.69e-03  8.97e-04  3.05e+01  2.55e+01  2.20e-01
  5  -5.1792e+05  -5.1788e+05  8.41e-05  2.22e-03  5.39e-04  4.44e+01  1.98e+01  3.30e-01
  6  -7.7519e+05  -7.7514e+05  7.13e-05  1.37e-03  3.32e-04  5.60e+01  1.45e+01  3.72e-01
  7  -8.6242e+05  -8.6237e+05  6.44e-05  1.10e-03  2.66e-04  5.62e+01  1.21e+01  2.38e-01
  8  -8.7572e+05  -8.7566e+05  5.90e-05  9.33e-04  2.26e-04  5.23e+01  1.03e+01  2.35e-01
  9  -8.7381e+05  -8.7377e+05  5.12e-05  7.88e-04  1.91e-04  4.52e+01  8.64e+00  2.13e-01
 10  -8.0331e+05  -8.0328e+05  3.94e-05  5.53e-04  1.34e-04  3.20e+01  5.77e+00  4.97e-01
 11  -7.6529e+05  -7.6526e+05  3.36e-05  4.20e-04  1.01e-04  2.59e+01  4.25e+00  4.44e-01
 12  -7.4962e+05  -7.4961e+05  2.49e-05  2.95e-04  7.11e-05  1.88e+01  2.94e+00  4.23e-01
 13  -7.5090e+05  -7.5089e+05  2.17e-05  2.05e-04  4.95e-05  1.64e+01  2.03e+00  6.92e-01
 14  -7.4220e+05  -7.4219e+05  1.34e-05  1.17e-04  2.83e-05  9.98e+00  1.15e+00  6.79e-01
 15  -7.3116e+05  -7.3116e+05  7.48e-06  6.46e-05  1.55e-05  5.50e+00  6.25e-01  5.42e-01
 16  -7.2435e+05  -7.2434e+05  3.51e-06  2.90e-05  6.97e-06  2.55e+00  2.78e-01  7.46e-01
 17  -7.1920e+05  -7.1920e+05  1.74e-06  1.42e-05  3.41e-06  1.26e+00  1.35e-01  6.34e-01
 18  -7.1909e+05  -7.1909e+05  1.63e-06  1.33e-05  3.19e-06  1.18e+00  1.26e-01  2.01e-01
 19  -7.1601e+05  -7.1600e+05  7.13e-07  5.74e-06  1.38e-06  5.14e-01  5.44e-02  7.54e-01
 20  -7.1424e+05  -7.1424e+05  2.61e-07  2.09e-06  5.01e-07  1.87e-01  1.97e-02  8.37e-01
 21  -7.1373e+05  -7.1373e+05  1.37e-07  1.10e-06  2.64e-07  9.86e-02  1.04e-02  8.37e-01
 22  -7.1341e+05  -7.1341e+05  6.61e-08  5.29e-07  1.27e-07  4.74e-02  4.98e-03  6.06e-01
 23  -7.1320e+05  -7.1320e+05  2.16e-08  1.73e-07  4.14e-08  1.55e-02  1.62e-03  8.38e-01
 24  -7.1315e+05  -7.1315e+05  9.32e-09  7.45e-08  1.79e-08  6.69e-03  7.01e-04  9.90e-01
 25  -7.1312e+05  -7.1312e+05  3.67e-09  2.93e-08  7.04e-09  2.63e-03  2.76e-04  8.70e-01
 26  -7.1311e+05  -7.1311e+05  1.52e-09  1.22e-08  2.92e-09  1.09e-03  1.14e-04  7.66e-01
 27  -7.1311e+05  -7.1311e+05  8.49e-10  6.79e-09  1.63e-09  6.09e-04  6.37e-05  7.22e-01
---------------------------------------------------------------------------------------------
Terminated with status = Solved
solve time = 6.56012246s

Increasing float precision

Hey,
I was wondering if there is a possibility to recompile Clarabel with increased floating point precision for increased overall precision for special problems that have such requirements. If I were to do it myself would there be anything I need to know about?
Thanks~
Matz

Convergence Difference Between Sedumi and Clarabel for Linear Problems

Im running sedumi/clarabel for a linear optimization problem P = 0. Settings are as in the example i.e. equilibrate is enabled

Sedumi:

K.l = size(A,2); % Linear Cone equivalent to Clarabel NonNegative
pars.eps = 1e-8;
pars.maxiter = 1000;
pars.theta = 0.25;
pars.beta = 0.5;
[X,Y,info] = sedumi(A,B,C,K,pars);



SeDuMi 1.3.7 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
eqs m = 2328, order n = 3721, dim = 3721, blocks = 1
nnz(A) = 8060 + 0, nnz(ADA) = 13799, nnz(L) = 8064
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            5.16E+03 0.000
  1 :  -9.23E+03 2.83E+03 0.000 0.5492 0.9000 0.9000   4.33  1  1  4.1E+01
  2 :  -8.81E+02 8.64E+02 0.000 0.3051 0.9000 0.9000   2.57  1  1  8.2E+00
  3 :   1.11E+02 7.55E+01 0.000 0.0874 0.9900 0.9900   1.29  1  1  8.5E-01
  4 :  -4.18E-01 5.54E+00 0.000 0.0733 0.9900 0.9900   1.02  1  1  4.0E-01
  5 :  -1.50E+01 5.29E-01 0.000 0.0956 0.9900 0.9900   1.00  1  1  5.4E-02
  6 :  -1.55E+01 1.18E-05 0.000 0.0000 0.9999 0.9999   1.00  1  1
iter seconds digits       c*x               b*y
  6      0.0  14.0 -1.5500000000e+01 -1.5500000000e+01
|Ax-b| =   1.0e-11, [Ay-c]_+ =   1.3E-16, |x|=  3.9e+01, |y|=  1.7e+04

Detailed timing (sec)
   Pre          IPM          Post
1.422E-01    4.195E-02    1.509E-03
Max-norms: ||b||=1, ||c|| = 100,
Cholesky |add|=0, |skip| = 3, ||L.L|| = 2.53566.

Clarabel:

-------------------------------------------------------------
           Clarabel.rs v0.5.0  -  Clever Acronym              

                   (c) Paul Goulart                          
                University of Oxford, 2022                   
-------------------------------------------------------------

problem:
  variables     = 3720
  constraints   = 2328
  nnz(P)        = 0
  nnz(A)        = 8660160
  cones (total) = 1
    : Nonnegative = 1,  numel = 2328

settings:
  linear algebra: direct / qdldl, precision: 64 bit
  max iter = 50, time limit = Inf,  max step = 0.990
  tol_feas = 1.0e-8, tol_gap_abs = 1.0e-8, tol_gap_rel = 1.0e-8,
  static reg : on, ϵ1 = 1.0e-8, ϵ2 = 4.9e-32
  dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-7
  iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
               max iter = 10, stop ratio = 5.0
  equilibrate: on, min_scale = 1.0e-4, max_scale = 1.0e4
               max iter = 10

iter    pcost        dcost       gap       pres      dres      k/t        μ       step      
---------------------------------------------------------------------------------------------
  0  +1.5933e+04  +1.3028e+06  8.08e+01  6.27e-01  2.29e+00  1.00e+00  7.14e+02   ------   
  1  +1.5803e+04  +1.3028e+06  8.14e+01  8.54e-02  2.16e+00  1.30e+02  7.06e+02  1.91e-01  
  2  +1.5748e+04  +1.3028e+06  8.17e+01  3.81e-02  2.03e+00  1.85e+02  7.01e+02  9.74e-02  
  3  +1.5650e+04  +1.3028e+06  8.22e+01  8.29e-02  1.98e+00  2.83e+02  6.45e+02  2.19e-01  
  4  +1.5390e+04  +1.3028e+06  8.37e+01  2.02e-02  1.81e+00  5.43e+02  5.40e+02  5.01e-01  
  5  +1.4692e+04  +1.3028e+06  8.77e+01  9.36e-03  2.29e+00  1.24e+03  2.90e+02  7.72e-01  
  6  +8.8915e+03  +1.3028e+06  1.46e+02  4.08e-03  1.68e+00  7.04e+03  5.50e+01  9.90e-01  
  7  -5.8897e+05  +1.3028e+06  3.21e+00  4.67e-05  7.57e-02  6.05e+05  6.40e-01  9.89e-01  
  8  -6.0478e+07  +1.3028e+06  4.74e+01  4.67e-07  7.58e-04  6.05e+07  6.40e-03  9.90e-01  
  9  -6.0494e+09  +1.3028e+06  4.64e+03  4.67e-09  7.58e-06  6.05e+09  6.40e-05  9.90e-01  
---------------------------------------------------------------------------------------------
Terminated with status = PrimalInfeasible
solve time = 151.522276563s
The application panicked (crashed).
Message:  capacity overflow
Location: library/alloc/src/raw_vec.rs:525

The only potential issue I might have is converting sparse nalgebra matricies to clarabel. But currently there is no way to probe the CSC matrices. Where A is a DMatrix and B and C are DVectors. Changing type to f32 doesnt help getting Nans

For reference:

    let B_vec = B.iter().map(|v| *v).collect::<Vec<_>>();
    let C_vec = C.iter().map(|v| *v).collect::<Vec<_>>();

    let mut A_coo_na = nalgebra_sparse::CooMatrix::<Float>::new(A.nrows(),A.ncols());
    for r in 0..A.nrows() {
        for c in 0..A.ncols() {
            A_coo_na.push(r,c,A[(r,c)]);
        }
    }
    let A_csc_na = nalgebra_sparse::CscMatrix::<Float>::from(&A_coo_na);

    let A_scs_clarabel = clarabel::algebra::CscMatrix::<Float>::new(
        A.nrows(),
        A.ncols(),
        A_csc_na.col_offsets().to_vec(),
        A_csc_na.row_indices().to_vec(),
        A_csc_na.values().to_vec()
    );

    let cones = [clarabel::solver::NonnegativeConeT::<Float>(A.nrows())];
    let P_cl = clarabel::algebra::CscMatrix::<Float>::spalloc(C.nrows(), C.nrows(), 0);

    (P_cl, A_scs_clarabel,B_vec,C_vec,cones)

Extend python interface to support data updating

Python wrapper does not yet support data updating features implemented in Rust in #63.

ImportError: DLL load failed while importing clarabel: The specified procedure could not be found.

I am using Python 3.8 on a conda venv , I installed clarabel using pip and I am facing errors in importing clarabel

Python 3.8.0 (default, Nov  6 2019, 16:00:02) [MSC v.1916 64 bit (AMD64)] :: Anaconda, Inc. on win32 
>>> from clarabel import *
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "D:\Anaconda3\envs\<env>\lib\site-packages\clarabel\__init__.py", line 1, in <module>
    from .clarabel import *
ImportError: DLL load failed while importing clarabel: The specified procedure could not be found.

Resolve new clippy warnings in rustc >=1.75

Rustc v1.75 is issuing new lint warnings, all of which relate to unused imports. This does not effect solver functionality, but putting a marker down here to fix so I don't forget.

What about https://github.com/oxfordcontrol/Clarabel.rs/blob/main/src/python/.github/workflows/CI.yml

I don't understand the exact purpose of this file. When would it be triggered?
I should also be taken care of by the dependabot as some of the GitHub actions are rather dated by now

Clarabel Panics rather than indicating insufficient progress with mixed cones.

When Clarabel is solving problems with mixed PSD and exponential cone constraints, Clarabel panics when it cannot make sufficient progress rather than setting it's status to Insufficient Progress.

An example is the following program:

In it's default settings, Clarabel solves the program with no issues.
However, by setting max_step_fraction = 1e-10, we can get Clarabel to throw the stack trace at the bottom of this thread. If the the log det Z >= 0 constraint is removed and the small step-size is used then Clarabel accurately throws Insufficient Progress.

If the max_step_fraction = 1e-3 then Clarabel also panics on the original program, but the program without the log det Z >= 0 constraint reports MaxIterations

Please note that reducing the max_step_fraction to be so small is only done for illustrative purposes. Such small step sizes are currently causing issues in the Drake CI: RobotLocomotion/drake#20799 (comment). A discussion is available in the Drake Developer Slack:
https://drakedevelopers.slack.com/archives/C3L92BM2Q/p1705935880972559

Stack Trace:

-------------------------------------------------------------
           Clarabel.rs v0.6.0  -  Clever Acronym              

                   (c) Paul Goulart                          
                University of Oxford, 2022                   
-------------------------------------------------------------

problem:
  variables     = 20
  constraints   = 32
  nnz(P)        = 0
  nnz(A)        = 39
  cones (total) = 6
    :        Zero = 1,  numel = 13
    : Nonnegative = 1,  numel = 0
    : Exponential = 2,  numel = (3,3)
    : PSDTriangle = 2,  numel = (3,10)

settings:
  linear algebra: direct / qdldl, precision: 64 bit
  max iter = 200, time limit = Inf,  max step = 0.001
  tol_feas = 1.0e-8, tol_gap_abs = 1.0e-8, tol_gap_rel = 1.0e-8,
  static reg : on, ϵ1 = 1.0e-8, ϵ2 = 4.9e-32
  dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-7
  iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
               max iter = 10, stop ratio = 5.0
  equilibrate: on, min_scale = 1.0e-4, max_scale = 1.0e4
               max iter = 10

iter    pcost        dcost       gap       pres      dres      k/t        μ       step      
---------------------------------------------------------------------------------------------
  0  +0.0000e+00  -1.1128e+00  1.11e+00  7.77e-01  6.58e-01  1.00e+00  1.00e+00   ------   
  1  +0.0000e+00  -1.1128e+00  1.11e+00  7.77e-01  6.58e-01  1.00e+00  1.00e+00  0.00e+00  
thread '<unnamed>' panicked at external/crate__clarabel-0.6.0/src/solver/core/cones/psdtrianglecone.rs:292:9:
not implemented: Mixed PSD and Exponential/Power cones are not yet supported
stack backtrace:
   0: rust_begin_unwind
             at /rustc/82e1608dfa6e0b5569232559e3d385fea5a93112/library/std/src/panicking.rs:645:5
   1: core::panicking::panic_fmt
             at /rustc/82e1608dfa6e0b5569232559e3d385fea5a93112/library/core/src/panicking.rs:72:14
   2: <clarabel::solver::core::cones::psdtrianglecone::PSDTriangleCone<T> as clarabel::solver::core::cones::Cone<T>>::compute_barrier
             at external/crate__clarabel-0.6.0/src/solver/core/cones/psdtrianglecone.rs:292:9
   3: <clarabel::solver::core::cones::supportedcone::SupportedCone<T> as clarabel::solver::core::cones::Cone<T>>::compute_barrier
             at external/crate__clarabel-0.6.0/src/solver/core/cones/mod.rs:42:1
   4: <clarabel::solver::core::cones::compositecone::CompositeCone<T> as clarabel::solver::core::cones::Cone<T>>::compute_barrier
             at external/crate__clarabel-0.6.0/src/solver/core/cones/compositecone.rs:363:24
   5: <clarabel::solver::implementations::default::variables::DefaultVariables<T> as clarabel::solver::core::traits::Variables<T>>::barrier
             at external/crate__clarabel-0.6.0/src/solver/implementations/default/variables.rs:213:20
   6: <clarabel::solver::core::solver::Solver<D,V,R,K,C,I,SO,SE> as clarabel::solver::core::solver::internal::IPSolverInternals<T,D,V,R,K,C,I,SO,SE>>::backtrack_step_to_barrier
             at external/crate__clarabel-0.6.0/src/solver/core/solver.rs:483:31
   7: <clarabel::solver::core::solver::Solver<D,V,R,K,C,I,SO,SE> as clarabel::solver::core::solver::internal::IPSolverInternals<T,D,V,R,K,C,I,SO,SE>>::get_step_length
             at external/crate__clarabel-0.6.0/src/solver/core/solver.rs:473:22
   8: <clarabel::solver::core::solver::Solver<D,V,R,K,C,I,SO,SE> as clarabel::solver::core::solver::IPSolver<T,D,V,R,K,C,I,SO,SE>>::solve
             at external/crate__clarabel-0.6.0/src/solver/core/solver.rs:330:18
   9: clarabel_cpp_rust_wrapper::solver::implementations::default::solver::_internal_DefaultSolver_solve
             at external/clarabel_cpp_internal/rust_wrapper/src/solver/implementations/default/solver.rs:109:5
  10: clarabel_DefaultSolver_f64_solve
             at external/clarabel_cpp_internal/rust_wrapper/src/solver/implementations/default/solver.rs:117:5

Large Scale Sparse Linear Programs

I am attempting to use Clarabel to solve a large scale, sparse linear program, as accessed through the CVXPY interface to Clarabel. However, I'm having some trouble, in that it appears the linear system solving step takes prohibitively long. For example, the following instance with about 5 million constraints and 400 thousand variables hangs for over 30 minutes with the following output.

-------------------------------------------------------------
           Clarabel.rs v0.7.1  -  Clever Acronym

                   (c) Paul Goulart                                                                                                                                                                           University of Oxford, 2022
-------------------------------------------------------------

problem:
  variables     = 412686
  constraints   = 4749546
  nnz(P)        = 0
  nnz(A)        = 12581604
  cones (total) = 2
    :        Zero = 1,  numel = 8145
    : Nonnegative = 1,  numel = 4741401

settings:
  linear algebra: direct / qdldl, precision: 64 bit
  max iter = 200, time limit = Inf,  max step = 0.990
  tol_feas = 1.0e-8, tol_gap_abs = 1.0e-8, tol_gap_rel = 1.0e-8,
  static reg : on, ϵ1 = 1.0e-8, ϵ2 = 4.9e-32
  dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-7
  iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
               max iter = 10, stop ratio = 5.0
  equilibrate: on, min_scale = 1.0e-5, max_scale = 1.0e5
               max iter = 10

iter    pcost        dcost       gap       pres      dres      k/t        μ       step
---------------------------------------------------------------------------------------------

In attempting to understand this, I have the following three questions:

Is Clarabel appropriate for large scale, sparse linear programming problems?
Could this be an issue with the linear system solver QDLDL?
On your benchmarks in the manuscript, Clarabel outperforms HiGHS, however, my tests show that Clarabel is massively outperformed by HiGHS for the large-scale, sparse linear programs I am studying. Do you have any thoughts on this?

Thank you for your assistance.

SDP option to use blas-sys and lapack-sys

Here's a feature request ...

Is your feature request related to a problem? Please describe.

In RobotLocomotion/drake#19449, I'm working to add Clarabel integration into Drake. Most users of Drake use our pre-compiled binary releases instead of building from source.

As such, it is not practical for us to force a choice of BLAS/LAPACK at compile time. Instead, we dynamically link to the libblas.so and liblapack.so (libblas.dylib and liblapack.dylib) from the user's operating system, so that their choice of the best BLAS for their own computer applications is respected. It allows us to ship one binary release of Drake that can perform well on many different computers.

At the moment, it does not appear that Clarabel.rs or Clarabel.cpp offer a feature to use the system BLAS, which prevents us from extending our BLAS convention to Clarabel as well.

Describe the solution you'd like

See https://github.com/blas-lapack-rs/blas-lapack-rs.github.io/wiki for background.

The Cargo.toml as of 0.6.0 currently offers these features:

# enables blas/lapack for SDP support, with blas/lapack src unspecified
sdp = ["blas","lapack"]

# explicit configuration options for different blas flavours
sdp-accelerate = ["sdp", "blas-src/accelerate", "lapack-src/accelerate"]
sdp-netlib     = ["sdp", "blas-src/netlib", "lapack-src/netlib"]
sdp-openblas   = ["sdp", "blas-src/openblas", "lapack-src/openblas"]
sdp-mkl        = ["sdp", "blas-src/intel-mkl", "lapack-src/intel-mkl"]
sdp-r          = ["sdp", "blas-src/r", "lapack-src/r"]

I'd like to request a feature option to use the system blas, something like:

sdp-sys        = ["sdp", "blas-sys", "lapack-sys"]

I believe it also might require changes to src/algebra/dense/blas.rs etc to spell things a little differently.

Describe alternatives you've considered

I've been trying to prepare my own patches to accomplish this, but doing it with downstream patching it fairly awkward, and my unfamiliarity with Rust is slowing me down a bit. Since I would plan to upstream the patches anyway, I thought I'd ask here to see if we could find a way to ship 0.7.0 with blas-sys out of the box. Y'all might be able to do this more quickly than myself.

Regression from ECOS causing CVXPY test failures

We were investigating some of the test case failures as we migrate the CVXPY test suite from ECOS to Clarabel in cvxpy/cvxpy#2479

It seems at least some of the cases are related to when ECOS correctly reports PRIMAL INFEASIBLE, but Clarabel reports InsufficientProgress instead.

This breaks our disciplined quasiconvex programing work and is consistent with issues that the EE 364a course staff had with Clarabel when assigning bisection methods in the homework.

Frankly, I'm a little worried about what to do here; the 1.5.x CVXPY release apparently regresses the DQCP support compared to prior releases. I'm not sure if this is just an issue of different default feastols or if this is from a new infeasibility detection approach.

Here's a minimized example:

import scipy.sparse as sp
import numpy as np

def make_array(d):
    out = sp.dok_matrix((8, 3))
    for (i, j), v in d.items():
        out[i, j] = v
    return out.tocsc()

G = make_array({
    (0, 0): -1.0,
    (2, 0): 1.0,
    (4, 0): -1.0,
    (1, 1): -1.0,
    (3, 1): 1.0,
    (5, 1): -1.0,
    (6, 1): -1.0,
    (4, 2): 72.001953125,
    (5, 2): -1.0,
    (6, 2): 1.0})
h = np.array([ 0.,  0., 12.,  6.,  0.,  0.,  0.,  2.])
c = np.array([0., 0., 0.])
nonneg = 5
soc = [3]

import ecos

solution = ecos.solve(c, G, h, {'l': 5, 'q': [3]})


import clarabel

P = sp.dok_matrix((3, 3)).tocsc()

cones = [clarabel.NonnegativeConeT(5), clarabel.SecondOrderConeT(3)]
settings = clarabel.DefaultSettings()
solver = clarabel.DefaultSolver(P,c,G,h,cones,settings)
solution = solver.solve()

Custom linear solver

Hey!

I was trying to use MKL as a custom direct solver via the Python setting direct_solve_method, however, I noticed in the Rust source that so far only qdldl is supported. Is there any plan to have support for MKL or other linear solvers?

return cost_dual in the solution

Thank you for the project. It is both fast and stable for our usage.

Could we add the cost_dual in the solution? (and the Python API )

Clarabel.rs/src/solver/implementations/default/solution.rs

Line 61 in b6f7bf9

self.obj_val = info.cost_primal;

Clarabel.rs/src/python/impl_default_py.rs

Line 68 in b6f7bf9

obj_val: result.obj_val,

Use case:
Sometimes, we got AlmostSolved. Currently, we print a warning message for ourselves and recheck manually later.

However, from the the value returned in solution, we don't know the tol_gap, so we would need to regenerate the whole data and re-run the solver with verbose=True to get the that information.

Ideally, we would like to investigate cases that has wider gap, e.g. we would like to check cases with 5e-5 gap instead of 1.1e-8 case. But currently, we would have to run the solver again first to know that.

Another benefit is that, we probably can downgrade the 1e-7 case from warning to info. (We still want to keep the 1e-8 for tol_gap so the solver tries harder instead of declare optimal at 1e-7)

crates.io tags

Currently these labels are defined in the cargo.toml: keywords = ["convex", "optimization", "QP", "SOCP", "SDP"]

There is only one package with the tags SDP and QP. The Tag SDP seems to be used mostly for a "Session Description Protocol".

I would propose to replace one of these labels with linear programming. The packages using it are mostly mathematical optimization packages: https://crates.io/search?q=linear%20programming
This would it make easier to find this amaizing package.

Recover s (slack?) after optimisation?

Is it possible to recover the s vector after the optimization has completed? I cant seem to find any methods that will return this vector.

Best

Infinite values in nonnegative cones?

Would it be possible for Clarabel to handle $\pm\infty$ values in non-negative cones? (Alternatively the solver could provide a "big float" threshold above which it considers values are infinite.)

At present having an infinity values (e.g. np.inf in Python which maps to an IEEE 754 infinity) results in a numerical error.

Context: this issue arose from testing Clarabel on the QPTEST problem from the Maros-Meszaros test set.

Reproduction steps

import clarabel
import numpy as np
import scipy.sparse as spa

# This is not the full QPTEST problem, just the gist of the issue
P = spa.csc_matrix(np.array([[8.0, 2.0], [2.0, 10.0]]))
q = np.array([1.5, -2.0])
A = spa.eye(2, format="csc")
b = np.array([20.0, np.inf])

cones = [clarabel.NonnegativeConeT(2)]
settings = clarabel.DefaultSettings()
settings.verbose = True
solver = clarabel.DefaultSolver(P, q, A, b, cones, settings)
result = solver.solve()
print(f"{result.status}")

Output

-------------------------------------------------------------
           Clarabel.rs v0.4.0  -  Clever Acronym              

                   (c) Paul Goulart                          
                University of Oxford, 2022                   
-------------------------------------------------------------

problem:
  variables     = 2
  constraints   = 2
  nnz(P)        = 3
  nnz(A)        = 2
  cones (total) = 1
    : Nonnegative = 1,  numel = 2

settings:
  linear algebra: direct / qdldl, precision: 64 bit
  max iter = 200, time limit = Inf,  max step = 0.990
  tol_feas = 1.0e-8, tol_gap_abs = 1.0e-8, tol_gap_rel = 1.0e-8,
  static reg : on, ϵ1 = 1.0e-8, ϵ2 = 4.9e-32
  dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-7
  iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
               max iter = 10, stop ratio = 5.0
  equilibrate: on, min_scale = 1.0e-4, max_scale = 1.0e4
               max iter = 10

iter    pcost        dcost       gap       pres      dres      k/t        μ       step      
---------------------------------------------------------------------------------------------
  0       NaN       NaN     NaN     NaN     NaN  1.00e+00     NaN   ------   
  1       NaN       NaN     NaN     NaN     NaN  1.00e+00     NaN  0.00e+00  
---------------------------------------------------------------------------------------------
Terminated with status = NumericalError
solve time = 62.359µs
NumericalError

Clarabel 0.4.1 v.s. 0.5.0 - Is a (relatively) high difference in dual variables expected?

Hello there,

I am trying to migrate to version 0.5.0 of Clarabel in Python.

On my test set of problems, I noticed neglictible differences in primal variables (let's say, about the order of the default eps 1e-8), but I also noticed non-neglictible differences in dual variables of ~5e-3.

Attached a reproducer, with output.

=> My question would be: are these differences expected?

Cheers,

Roman
clarabel.txt

Solver should detect and report nonconvexity

PanicException: assertion `left == right` failed

Hi and thanks again for maintaining this open source solver 😃 I noticed the following exception thrown by Clarabel on the QUADCMPC1 problem of the free-for-all QP test set:

$ python free_for_all.py run --solver clarabel --problem QUADCMPC1
results_file='src/free_for_all_qpbenchmark/results/free_for_all.csv'
[2024-07-31 14:30:39,925] [info] Loading existing results from src/free_for_all_qpbenchmark/results/free_for_all.csv (results.py:78)
[2024-07-31 14:30:39,953] [info] Solving QUADCMPC1 by clarabel with default settings... (run.py:108)
thread '<unnamed>' panicked at src/solver/implementations/default/solver.rs:56:9:
assertion `left == right` failed
  left: 1536
 right: 1408
note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace
Traceback (most recent call last):
  File "src/free_for_all_qpbenchmark/free_for_all.py", line 231, in <module>
    main(test_set_path=os.path.abspath(__file__))
  File ".micromamba/envs/qpbenchmark/lib/python3.12/site-packages/qpbenchmark/benchmark.py", line 268, in main
    run(
  File ".micromamba/envs/qpbenchmark/lib/python3.12/site-packages/qpbenchmark/run.py", line 113, in run
    solution, runtime = time_solve_problem(
                        ^^^^^^^^^^^^^^^^^^^
  File ".micromamba/envs/qpbenchmark/lib/python3.12/site-packages/qpbenchmark/utils.py", line 129, in time_solve_problem
    solution = qpsolvers.solve_problem(problem, solver=solver, **kwargs)
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File ".micromamba/envs/qpbenchmark/lib/python3.12/site-packages/qpsolvers/solve_problem.py", line 81, in solve_problem
    return solve_function[solver](problem, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File ".micromamba/envs/qpbenchmark/lib/python3.12/site-packages/qpsolvers/solvers/clarabel_.py", line 118, in clarabel_solve_problem
    solver = clarabel.DefaultSolver(P, q, A_stack, b_stack, cones, settings)
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
pyo3_runtime.PanicException: assertion `left == right` failed

I checked quickly that A_stack.shape=(2432, 768) and b_stack.shape=(2432,) (could be "left" and "right"?). Other solvers like OSQP solve the problem successfully.

Sorry I don't have the time right now to make a reproducer that does not depend on the benchmarking library, but the exception should be reproducible (after conda/PyPI installing qpbenchmark) by:

$ python free_for_all.py run --solver clarabel --problem QUADCMPC1

Finally I noticed that Clarabel 0.7.0 does not raise the exception, while Clarabel 0.8.0 does.

Hoping this helps!

Rust implementation significantly slower than Julia implementation

As noted in #36 the performance of the Julia and Rust implementation should be fairly similar, however for the models that I am currently solving the performance of the Rust implementation is significantly slower (e.g. around 18 times slower).

I am a complete beginner at Rust, so I initially thought that the issue could be due to the way that I set up the problem. I therefore tried setting up the same model using the C++ interface which ended up with similar same results. For me this indicates that it is not the model formulation is not the issue.

That being said the C++ and Rust models have something in common, namely that the model matrices A, b, P and q as well as the constraint type information are exported from the Julia model into .mtx files which are then loaded in both Rust and C++. Could there be an issue hiding somewhere here? If not I am quite lost as to where my issue might be.

Below is the results from one the models that I solved using Rust, the C++ interface, and Julia. The number of iterations are similar, so it is seemingly just that each iterations takes significantly longer for the Rust/C++ models than the Julia model.

-------------------------------------------------------------
           Clarabel.rs v0.7.1  -  Clever Acronym              

                   (c) Paul Goulart                          
                University of Oxford, 2022                   
-------------------------------------------------------------

problem:
  variables     = 12776
  constraints   = 19626
  nnz(P)        = 0
  nnz(A)        = 60848
  cones (total) = 1155
    :        Zero = 1,  numel = 8105
    : Nonnegative = 2,  numel = (1,8064)
    : SecondOrder = 1152,  numel = (3,3,3,3,...,3)

settings:
  linear algebra: direct / qdldl, precision: 64 bit
  max iter = 200, time limit = Inf,  max step = 0.990
  tol_feas = 1.0e-8, tol_gap_abs = 1.0e-8, tol_gap_rel = 1.0e-8,
  static reg : on, ϵ1 = 1.0e-8, ϵ2 = 4.9e-32
  dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-7
  iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
               max iter = 10, stop ratio = 5.0
  equilibrate: on, min_scale = 1.0e-5, max_scale = 1.0e5
               max iter = 10

iter    pcost        dcost       gap       pres      dres      k/t        μ       step      
---------------------------------------------------------------------------------------------
  0  -1.8583e-01  -4.2644e+04  4.26e+04  3.83e-01  2.45e-01  1.00e+00  9.47e+00   ------   
  1  -7.5069e+00  -2.2558e+04  3.00e+03  2.54e-01  7.43e-02  2.14e+00  4.96e+00  6.51e-01  
  2  -7.5394e+04  -2.3308e+05  2.09e+00  7.78e-03  6.01e-03  2.95e+01  4.95e+00  6.87e-02  
  3  -7.5199e+04  -1.7746e+05  1.36e+00  5.06e-03  3.87e-03  6.20e+01  3.41e+00  4.70e-01  
  4  -1.0995e+05  -1.6523e+05  5.03e-01  1.87e-03  1.43e-03  5.72e+01  1.58e+00  5.80e-01  
  5  -8.7170e+04  -1.1454e+05  3.14e-01  1.17e-03  8.74e-04  2.50e+01  1.00e+00  7.78e-01  
  6  -9.2001e+04  -1.1352e+05  2.34e-01  8.71e-04  6.48e-04  1.96e+01  7.95e-01  3.61e-01  
  7  -9.7567e+04  -1.0801e+05  1.07e-01  3.99e-04  2.94e-04  9.66e+00  4.15e-01  9.90e-01  
  8  -1.0425e+05  -1.0817e+05  3.76e-02  1.40e-04  1.02e-04  2.63e+00  1.77e-01  6.97e-01  
  9  -1.0142e+05  -1.0354e+05  2.09e-02  7.80e-05  5.65e-05  1.74e+00  1.09e-01  6.00e-01  
 10  -1.0096e+05  -1.0148e+05  5.06e-03  1.89e-05  1.36e-05  4.10e-01  2.93e-02  7.74e-01  
 11  -1.0085e+05  -1.0097e+05  1.15e-03  4.30e-06  3.08e-06  9.64e-02  6.91e-03  8.40e-01  
 12  -1.0079e+05  -1.0082e+05  2.53e-04  9.44e-07  6.76e-07  2.15e-02  1.53e-03  9.45e-01  
 13  -1.0078e+05  -1.0080e+05  1.68e-04  6.24e-07  4.46e-07  1.41e-02  1.01e-03  5.76e-01  
 14  -1.0076e+05  -1.0077e+05  5.50e-05  2.05e-07  1.46e-07  4.63e-03  3.33e-04  7.09e-01  
 15  -1.0076e+05  -1.0076e+05  3.93e-05  1.47e-07  1.05e-07  3.30e-03  2.38e-04  4.54e-01  
 16  -1.0076e+05  -1.0076e+05  2.09e-05  7.79e-08  5.56e-08  1.75e-03  1.27e-04  7.13e-01  
 17  -1.0075e+05  -1.0075e+05  7.31e-06  2.72e-08  1.94e-08  6.05e-04  4.43e-05  9.90e-01  
 18  -1.0075e+05  -1.0075e+05  3.07e-06  1.14e-08  8.16e-09  2.54e-04  1.86e-05  5.95e-01  
 19  -1.0075e+05  -1.0075e+05  6.27e-07  2.34e-09  1.67e-09  5.19e-05  3.80e-06  9.22e-01  
 20  -1.0075e+05  -1.0075e+05  7.33e-08  2.73e-10  1.95e-10  6.03e-06  4.43e-07  9.90e-01  
 21  -1.0075e+05  -1.0075e+05  1.01e-08  3.77e-11  2.69e-11  8.29e-07  6.11e-08  8.89e-01  
 22  -1.0075e+05  -1.0075e+05  2.87e-09  1.07e-11  7.64e-12  2.34e-07  1.73e-08  8.64e-01  
---------------------------------------------------------------------------------------------
Terminated with status = Solved
solve time = 5.490045665s

-------------------------------------------------
           Clarabel.rs v0.7.0  -  Clever Acronym              

                   (c) Paul Goulart                          
                University of Oxford, 2022                   
-------------------------------------------------------------

problem:
  variables     = 12776
  constraints   = 19626
  nnz(P)        = 0
  nnz(A)        = 60848
  cones (total) = 1155
    :        Zero = 1,  numel = 8105
    : Nonnegative = 2,  numel = (1,8064)
    : SecondOrder = 1152,  numel = (3,3,3,3,...,3)

settings:
  linear algebra: direct / qdldl, precision: 64 bit
  max iter = 200, time limit = Inf,  max step = 0.990
  tol_feas = 1.0e-8, tol_gap_abs = 1.0e-8, tol_gap_rel = 1.0e-8,
  static reg : on, ϵ1 = 1.0e-8, ϵ2 = 4.9e-32
  dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-7
  iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
               max iter = 10, stop ratio = 5.0
  equilibrate: on, min_scale = 1.0e-5, max_scale = 1.0e5
               max iter = 10

iter    pcost        dcost       gap       pres      dres      k/t        μ       step      
---------------------------------------------------------------------------------------------
  0  -1.8583e-01  -4.2644e+04  4.26e+04  3.83e-01  2.45e-01  1.00e+00  9.47e+00   ------   
  1  -7.5069e+00  -2.2558e+04  3.00e+03  2.54e-01  7.43e-02  2.14e+00  4.96e+00  6.51e-01  
  2  -7.5394e+04  -2.3308e+05  2.09e+00  7.78e-03  6.01e-03  2.95e+01  4.95e+00  6.87e-02  
  3  -7.5199e+04  -1.7746e+05  1.36e+00  5.06e-03  3.87e-03  6.20e+01  3.41e+00  4.70e-01  
  4  -1.0995e+05  -1.6523e+05  5.03e-01  1.87e-03  1.43e-03  5.72e+01  1.58e+00  5.80e-01  
  5  -8.7170e+04  -1.1454e+05  3.14e-01  1.17e-03  8.74e-04  2.50e+01  1.00e+00  7.78e-01  
  6  -9.2001e+04  -1.1352e+05  2.34e-01  8.71e-04  6.48e-04  1.96e+01  7.95e-01  3.61e-01  
  7  -9.7567e+04  -1.0801e+05  1.07e-01  3.99e-04  2.94e-04  9.66e+00  4.15e-01  9.90e-01  
  8  -1.0425e+05  -1.0817e+05  3.76e-02  1.40e-04  1.02e-04  2.63e+00  1.77e-01  6.97e-01  
  9  -1.0142e+05  -1.0354e+05  2.09e-02  7.80e-05  5.65e-05  1.74e+00  1.09e-01  6.00e-01  
 10  -1.0096e+05  -1.0148e+05  5.06e-03  1.89e-05  1.36e-05  4.10e-01  2.93e-02  7.74e-01  
 11  -1.0085e+05  -1.0097e+05  1.15e-03  4.30e-06  3.08e-06  9.64e-02  6.91e-03  8.40e-01  
 12  -1.0079e+05  -1.0082e+05  2.53e-04  9.44e-07  6.76e-07  2.15e-02  1.53e-03  9.45e-01  
 13  -1.0078e+05  -1.0080e+05  1.68e-04  6.24e-07  4.46e-07  1.41e-02  1.01e-03  5.76e-01  
 14  -1.0076e+05  -1.0077e+05  5.50e-05  2.05e-07  1.46e-07  4.63e-03  3.33e-04  7.09e-01  
 15  -1.0076e+05  -1.0076e+05  3.93e-05  1.47e-07  1.05e-07  3.30e-03  2.38e-04  4.54e-01  
 16  -1.0076e+05  -1.0076e+05  2.09e-05  7.79e-08  5.56e-08  1.75e-03  1.27e-04  7.13e-01  
 17  -1.0075e+05  -1.0075e+05  7.31e-06  2.72e-08  1.94e-08  6.05e-04  4.43e-05  9.90e-01  
 18  -1.0075e+05  -1.0075e+05  3.07e-06  1.14e-08  8.16e-09  2.54e-04  1.86e-05  5.95e-01  
 19  -1.0075e+05  -1.0075e+05  6.27e-07  2.34e-09  1.67e-09  5.19e-05  3.80e-06  9.22e-01  
 20  -1.0075e+05  -1.0075e+05  7.33e-08  2.73e-10  1.95e-10  6.03e-06  4.43e-07  9.90e-01  
 21  -1.0075e+05  -1.0075e+05  1.01e-08  3.77e-11  2.69e-11  8.29e-07  6.11e-08  8.89e-01  
 22  -1.0075e+05  -1.0075e+05  2.87e-09  1.07e-11  7.64e-12  2.34e-07  1.73e-08  8.64e-01  
---------------------------------------------------------------------------------------------
Terminated with status = Solved
solve time = 5.331283753s

-------------------------------------------------------------
           Clarabel.jl v0.7.1  -  Clever Acronym              
                   (c) Paul Goulart                          
                University of Oxford, 2022                   
-------------------------------------------------------------

problem:
  variables     = 12776
  constraints   = 19626
  nnz(P)        = 0
  nnz(A)        = 60848
  cones (total) = 1155
    : Zero        = 1,  numel = 8105
    : Nonnegative = 2,  numel = (1,8064)
    : SecondOrder = 1152,  numel = (3,3,3,3,...,3)

settings:
  linear algebra: direct / qdldl, precision: Float64
  max iter = 200, time limit = Inf,  max step = 0.990
  tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
  static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
  dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
  iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12, 
               max iter = 10, stop ratio = 5.0
  equilibrate: on, min_scale = 1.0e-05, max_scale = 1.0e+05
               max iter = 10

iter    pcost        dcost       gap       pres      dres      k/t        μ       step      
---------------------------------------------------------------------------------------------
  0  -1.8432e-01  -4.2644e+04  4.26e+04  3.38e-01  3.24e-01  1.00e+00  9.37e+00   ------   
  1  -8.4105e+00  -2.2688e+04  2.70e+03  2.38e-01  2.34e-01  2.85e+00  4.63e+00  7.43e-01  
  2  -6.3248e+04  -2.0439e+05  2.23e+00  9.73e-03  1.36e-02  3.41e+01  4.62e+00  6.27e-02  
  3  -5.2950e+04  -1.4098e+05  1.66e+00  7.25e-03  1.01e-02  7.96e+01  3.42e+00  5.97e-01  
  4  -1.1244e+05  -1.5073e+05  3.41e-01  1.49e-03  2.09e-03  8.53e+01  1.04e+00  7.13e-01  
  5  -9.7052e+04  -1.2673e+05  3.06e-01  1.33e-03  1.88e-03  4.78e+01  9.40e-01  3.04e-01  
  6  -1.1303e+05  -1.2533e+05  1.09e-01  4.75e-04  6.67e-04  1.30e+01  3.96e-01  7.03e-01  
  7  -1.0243e+05  -1.1024e+05  7.63e-02  3.33e-04  4.66e-04  8.02e+00  3.11e-01  5.05e-01  
  8  -1.0229e+05  -1.0417e+05  1.84e-02  8.03e-05  1.12e-04  1.26e+00  9.27e-02  8.16e-01  
  9  -1.0172e+05  -1.0275e+05  1.01e-02  4.41e-05  6.12e-05  8.10e-01  5.43e-02  5.91e-01  
 10  -1.0111e+05  -1.0135e+05  2.38e-03  1.04e-05  1.44e-05  2.31e-01  1.36e-02  9.66e-01  
 11  -1.0085e+05  -1.0090e+05  5.48e-04  2.39e-06  3.31e-06  5.12e-02  3.21e-03  8.71e-01  
 12  -1.0083e+05  -1.0087e+05  4.15e-04  1.81e-06  2.50e-06  3.77e-02  2.43e-03  4.68e-01  
 13  -1.0077e+05  -1.0078e+05  8.93e-05  3.90e-07  5.39e-07  7.97e-03  5.25e-04  8.40e-01  
 14  -1.0076e+05  -1.0077e+05  5.70e-05  2.49e-07  3.44e-07  5.00e-03  3.36e-04  5.56e-01  
 15  -1.0076e+05  -1.0076e+05  3.09e-05  1.35e-07  1.86e-07  2.59e-03  1.82e-04  9.90e-01  
 16  -1.0076e+05  -1.0076e+05  1.50e-05  6.54e-08  9.04e-08  1.25e-03  8.82e-05  6.16e-01  
 17  -1.0075e+05  -1.0075e+05  2.44e-06  1.06e-08  1.47e-08  2.04e-04  1.44e-05  9.04e-01  
 18  -1.0075e+05  -1.0075e+05  8.96e-07  3.91e-09  5.40e-09  7.48e-05  5.27e-06  8.33e-01  
 19  -1.0075e+05  -1.0075e+05  1.59e-07  6.94e-10  9.59e-10  1.33e-05  9.36e-07  8.42e-01  
 20  -1.0075e+05  -1.0075e+05  2.73e-08  1.19e-10  1.64e-10  2.27e-06  1.60e-07  9.19e-01  
 21  -1.0075e+05  -1.0075e+05  6.27e-09  2.73e-11  3.78e-11  5.20e-07  3.69e-08  7.90e-01  
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time =  314ms

feature request - API to solve multiple times with different data

I'm interested in solving a sequence of problems with the same structure for online control, and I wasn't able to find documentation on this. Is it intended that solve can be called multiple times on the same solver, while updating problem data between each call?

I'm interested:

If it's possible to avoid allocation after the solver is created
If two calls to solve will give slightly different results due to re-use of x or s
If any of {P, q, A, b} are more expensive or cheap to update, due to avoiding re-factorization

Infeasibility depends on objective

I'm getting a behavior I don't understand where the solver sometimes fails with PrimalInfeasible depending on the objective function. My naive assumption is that feasibility of the primal is independent of the objective. As shown in the snippet, if I alter the objective the solver finishes without issue. Altering q does not appear to make a difference (aside from getting AlmostPrimalInfeasible instead of PrimalInfeasible). I'm using clarabel 0.6.0 (from conda forge).

Is this expected behavior, and is there a work-around? Thanks in advance.

import numpy as np
from clarabel import NonnegativeConeT, ZeroConeT, DefaultSolver, DefaultSettings

P=np.array([[0.,         0.,         0.        ],
            [0.,         0.,         0.        ],
            [0.,         0.,         0.683013]])

# With these entries result is Solved. Without the result is PrimalInfeasible
#P[0,0]=0.01
#P[1,1]=0.01

q=np.array([0.0,0.0, -7.47697849e-06])
A=np.array([[ 1.        ,  0.        ,  0.        ],
            [ 0.        ,  1.        ,  0.        ],
            [-1.        ,  0.        ,  0.        ],
            [ 0.        ,  0.        ,  1.        ],
            [ 0.        ,  0.        , -1.        ],
            [ 0.        , -1.        ,  0.        ],
            [ 1.        ,  1.        ,  1.        ],
            [ 0.02479339,  0.02479339, -0.05785124]] )

b=np.array([2000.,         2000.,         -400.        , 2000.,         -230.,
            -200.,         2200.,          -77.68595041,])
S=[NonnegativeConeT(1), NonnegativeConeT(1), NonnegativeConeT(1), NonnegativeConeT(1),
   NonnegativeConeT(1), NonnegativeConeT(1), NonnegativeConeT(1), ZeroConeT(1)]

solver = DefaultSolver(sparse.csc_matrix(P), q,
                       sparse.csc_matrix(A), b, S,
                       DefaultSettings())
soln=solver.solve()

InsufficientProgress status

I raised the following issue with CVXPY two weeks ago, but since it's related to both CVXPY and Clarabel, I would also like to raise it here.

For infeasible problems, I noticed that Clarabel sometimes returns an InsufficientProgress status. What is the interpretation of this?

I checked the docs, but this didn't help me much. I'm asking because the InsufficientProgress is mapped to a SolverError in CVXPY, which is creating some issues in my workflow where I catch infeasible solutions. Would it be appropriate to map this to a different status in CVXPY?

Thanks a lot!

Custom Newton's Step Solver: Exploiting Problem Structure

Is it possible to pass in a custom Newton's step solver which is able to exploit problem specific structure? In particular, one may like to pass in a custom linear solver to solve the linear system (eq. 10/11) for the update $(\Delta x, \Delta z, \Delta \tau)$. Such a feature would allow one to exploit problem specific structure, while still fitting within your efficient interior point method framework, which has extremely impressive performance compared to state-of-the-art solvers.

For example, in an application of mine to inferring the evolutionary history of tumors, I am able to solve the linear system (eq. 10 in preprint) in $\mathcal{O}(n + m)$ time using a combinatorial algorithm ($P = 0$ in my application and the cone $\mathcal{K}$ is the non-negative orthant), where $m$ and $n$ are the dimensions of the constraint matrix $A$. Incorporating such a procedure into your general framework would be extremely useful, as current state-of-the-art linear programming solvers such as Gurobi/CPLEX do not allow problem specific customization.

Issues when `tril` is chosen for a new ldlsolver

Hi Paul,

I am trying to add PanuaPardiso as a possible solver for the KKT system. Preferably I would like to use tril as in Julia (PanuaPardiso wants the upper part in a CSR format, so if I get tril in CSC I can simply swap the row/cols in order to get triu in CSR format). However, I've encountered two issues so far:

Adding the following to a new PanuaPardisoDirectLDLSolver result in an initial KKT system for which the diagonal is not completely filled (for a triu format the diagonal is filled - as is required by PanuaPardiso).

fn required_matrix_shape() -> MatrixTriangle {
      MatrixTriangle::Tril
  }

Later the code also panics (unsurprisingly) when it hits the following line

Clarabel.rs/src/algebra/csc/core.rs

Line 302 in d1ca25d

debug_assert!(self.is_triu());

I've tried to look into the code, but seems a little complicated, so my hope is that you might know an easy fix.

Best,
Mikkel