Rust 製数理最適化ソルバー clarabel の使い方

clarabel とは

https://docs.rs/clarabel/latest/clarabel/
https://github.com/oxfordcontrol/Clarabel.rs

Apache License 2.0 のもと公開されている， Rust 🦀 で書かれた最適化ソルバー．
次の形の最適化問題を扱える：

\begin{align*} \min \quad & \frac{1}{2} \boldsymbol{x}^{\top} Q \boldsymbol{x} + \boldsymbol{q}^{\top} \boldsymbol{x} \\ \text{s.t.} \quad & A \boldsymbol{x} + \boldsymbol{s} = \boldsymbol{b} \\ & \boldsymbol{s} \in \mathcal{K} \end{align*}

ただし，Q \in \mathbb{R}^{n \times n} は半正定値対称行列，\boldsymbol{q} \in \mathbb{R}^{n}，A \in \mathbb{R}^{m \times n}，\boldsymbol{b} \in \mathbb{R}^{m} で，\mathcal{K} \subset \mathbb{R}^{m} は凸錐．

私が知る限り日本語での使い方解説は（2024/12/14 時点で）存在しないので，この記事で紹介する．
（この記事 の余談で触れられている以外，見つけられなかった）

具体例で学ぶ使い方

ここでは，次の二次計画問題を解く

\begin{align} \min \quad & \frac{1}{2} \boldsymbol{x}^{\top} Q \boldsymbol{x} + \boldsymbol{q}^{\top} \boldsymbol{x} \\ \text{s.t.} \quad & A \boldsymbol{x} \le \boldsymbol{b} \nonumber \\ \end{align}

ここで，

Q = \begin{pmatrix} 2 & 1 \\ 1 & 3 \end{pmatrix}, \quad \boldsymbol{q} = \begin{pmatrix} -3 \\ -4 \end{pmatrix}, \quad A = \begin{pmatrix} 2 & 2 \\ 1 & 0 \\ -1 & 0 \\ 0 & -1 \end{pmatrix}, \quad \boldsymbol{b} = \begin{pmatrix} 3 \\ 1 \\ 0 \\ 0 \\ \end{pmatrix}

と定義した．
この最適化問題は，以下の図の緑の領域で最小値を求める問題（desmos で描画）

Clarabel では 行列を CSC フォーマットで記述するので，Q, A を次のように書き直す．
CSC フォーマットについては，NVIDIA の解説ページが参考になる．

\begin{alignat*}{2} c_{Q} & = \begin{pmatrix} 0 & 2 & 4 \end{pmatrix} & \qquad & \text{(Column offsets)} \\ r_{Q} & = \begin{pmatrix} 0 & 1 & 0 & 1 \end{pmatrix} & \qquad & \text{(Row indices)} \\ v_{Q} & = \begin{pmatrix} 3 & -1 & -1 & 2 \end{pmatrix} & \qquad & \text{(Values)} \\ c_{A} & = \begin{pmatrix} 0 & 3 & 5 \end{pmatrix} & \qquad & \text{(Column offsets)} \\ r_{A} & = \begin{pmatrix} 0 & 1 & 2 & 0 & 3 \end{pmatrix} & \qquad & \text{(Row indices)} \\ v_{A} & = \begin{pmatrix} 2 & 1 & -1 & 2 & -1 \end{pmatrix} & \qquad & \text{(Values)} \end{alignat*}

1. プロジェクトを作成．ここでは，test-clarabel という名前にする．

   bash

   $ cargo new --bin test-clarabel
   $ cd test-clarabel
   $ ls
   Cargo.toml  src/
   

2. Cargo.toml に clarabel を追加．ここでは，バージョンは 0.9.0 を指定．

   Cargo.toml

   [package]
   name = "test-clarabel"
   version = "0.1.0"
   edition = "2021"
   
   [dependencies]
   + clarabel = { version = "0.9.0" }
   

3. main.rs に最適化問題を記述．
   数式の表記に合わせて大文字でプログラムを書きたいので，大文字の変数で警告を出さないようにする !#[allow(non_snake_case)] を書いている．

   main.rs

   #![allow(non_snake_case)]
   
   use clarabel::{
       algebra::*,
       solver::*,
   };
   
   fn main() {
       // 目的関数を記述
       let Q = CscMatrix::new(
           2usize, // 行列の行数
           2usize, // 行列の列数
           vec![0usize, 2usize, 4usize],         // column offsets
           vec![0usize, 1usize, 0usize, 1usize], // row offsets
           vec![2f64, -1f64, -1f64, 3f64],       // values
       );
       let q = vec![-3f64, -4f64];
   
   
       // 制約式を記述
       let A = CscMatrix::new(
           4usize, // 行列の行数
           2usize, // 行列の列数
           vec![0usize, 3usize, 5usize],                 // column offsets
           vec![0usize, 1usize, 2usize, 0usize, 3usize], // row offsets
           vec![2f64, 1f64, -1f64, 2f64, -1f64],         // values
       );
       // ≤ の制約が 4 つなので，`NonnegativeConeT(4)` と書く．
       // 統合製薬の場合は，`ZeroConeT(個数)` と書く．
       let sense = vec![NonnegativeConeT(4)];
       let b = vec![3f64, 1f64, 0f64, 0f64];
   
   
       // ソルバの設定情報を作る．
       let settings = DefaultSettingsBuilder::default()
           .build()
           .unwrap();
   
       // ソルバを初期化
       let mut solver = DefaultSolver::new(
           &Q,
           &q,
           &A,
           &b,
           &sense,
           settings,
       );
   
       // 最適化問題を解く
       solver.solve();
   
       // ソルバ実行後の状態を取得．解けていれば `Solved`
       let status = solver.solution.status;
       println!("status is {status:?}");
   
       // 最適値を取得
       let optval = solver.solution.obj_val;
       println!("optimal value is {optval}");
   
       // 最適解を取得
       let solution = &solver.solution.x[..];
       println!("optimal solution is {solution:?}");
   
       // 双対問題の最適解を取得
       let solution = &solver.solution.z[..];
       println!("optimal dual solution is {solution:?}");
   }
   

   これを実行すると，次のような出力を得る．

   実行

   $ cargo run
   -------------------------------------------------------------
              Clarabel.rs v0.9.0  -  Clever Acronym
                     *** debug build ***
                      (c) Paul Goulart
                   University of Oxford, 2022
   -------------------------------------------------------------
   
   problem:
     variables     = 2
     constraints   = 4
     nnz(P)        = 3
     nnz(A)        = 5
     cones (total) = 1
       : Nonnegative = 1,  numel = 4
   
   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  -4.0589e+00  -1.8313e+01  3.51e+00  2.06e-01  3.21e-01  1.00e+00  1.76e+00   ------
     1  -4.1033e+00  -5.2943e+00  2.90e-01  1.77e-02  3.01e-02  3.98e-02  2.29e-01  9.34e-01
     2  -4.4090e+00  -4.4473e+00  8.70e-03  3.99e-04  6.84e-04  1.53e-03  8.52e-03  9.77e-01
     3  -4.4106e+00  -4.4115e+00  1.84e-04  4.03e-06  6.89e-06  3.57e-05  1.96e-04  9.90e-01
     4  -4.4107e+00  -4.4107e+00  1.91e-06  4.03e-08  6.89e-08  3.73e-07  2.05e-06  9.90e-01
     5  -4.4107e+00  -4.4107e+00  1.92e-08  4.03e-10  6.89e-10  3.74e-09  2.05e-08  9.90e-01
     6  -4.4107e+00  -4.4107e+00  1.92e-10  4.03e-12  6.89e-12  3.74e-11  2.05e-10  9.90e-01
   ---------------------------------------------------------------------------------------------
   Terminated with status = Solved
   solve time = 918.625µs
   status is Solved
   optimal value is -4.410714285645115
   optimal solution is [0.7142857140248071, 0.7857142859458478]
   optimal dual solution is [1.1785714281169977, 1.9263888919739336e-9, 1.3093200867623522e-10, 5.11559596564861e-11]
   

   最後の 4 行は main.rs で書いた println! の出力．
   最適解は (x_{1}, x_{2}) = (0.7142857140248071, 0.7857142859458478) で，最適値が -4.410714285645115 であることがわかる．

先ほどの図にこの解を描画してみると，正しく解けていることが確認できる．

まとめ

clarabel で最適化問題を解く際は， CSC フォーマット問題を記述する必要がある．
公式ドキュメントは少ないものの，最適化問題の基本的な部分を知っていれば直感的に使える．

補遺

• log を消す
  clarabel では，デフォルトでログを出力する．これを無効にしたい場合，設定を次のように書く．

    let settings = DefaultSettingsBuilder::default()
        .verbose(false) // log を出力しないようにする．デフォルトでは true
  +     .build()
        .unwrap();
  

• 線形計画問題（LP）を解く
  上の例では二次計画問題（QP）を解いた．LP を解く場合は，Q に零行列を設定すればよい．
  幸運なことに，簡単な記法が用意されている．

  let n_variables = 100; // 変数の数．
  let zero_mat = CscMatrix::<f64>::zeros((n_variables, n_variables));
  

• 行生成法をやる場合
  現状，私が知る限りでは，何度も CscMatrix を再構築する羽目になる．
  双対問題を考え，Incremental にやるとよい．
• 他の例
  GitHub で公開されている examples を参照するとよい．