cupdlpxのソリューション(2vector)をCOPTのソリューション(4vector)にセットする方法

cupdlpxで計算させると時々crossoverさせたくなります。そこで、COPTで計算させることにします。しかし、折角近くまで来ているのにまた最初から計算させるのは、勿体ないです。（COPTは、CrossoverAPIを備えています。）そこで、cuplpxの２vector solutionをCOPTの4vector solutionに変換する方法をcopilotに聞きました。

C++コードまで示してくれたのですが、Desnse行列になっているので、Sparse行列用に手直しする必要があります。また、COPTのslacks=b-Axとなっていますが、実際にCOPTが返してくる値をみると、slacks=Axとなっていました。マニュアルを見てもこの辺のことは書いていないのですが、これも修正する必要があります。(それはともかく、CopilotにOSI実装を相談した覚えはないのですが..)

---

🔎 背景整理

• COPT の LP/MIP インターフェース（例えば OSI や Python API）では、setSolution() を呼ぶ際に 4種類のベクトルを渡す必要があります。典型的には：

  1. primal solution（変数値）
  2. dual solution（制約の双対値）
  3. reduced costs（変数の双対値）
  4. slacks（制約の残差）

• 一方、cuPDLPX（GPUベースの PDLP 系ソルバ）は、主に primal と dual の2種類しか返しません。これはアルゴリズムの性質上、reduced cost や slack を直接保持していないためです。

---

🛠 残り2つを生成する方法

cuPDLPX が返す primal (x) と dual (y) から、COPT が要求する残りのベクトルを自前で計算する必要があります。

• slacks
  制約式 (Ax \leq b) に対して
  [ s = b - Ax ]
  を計算すれば得られます。

• reduced costs
  目的関数 (c^{T x) に対して
  [ r = c - A}T y ]
  を計算すれば得られます。

---

📌 実務的な流れ

1. cuPDLPX から primal 解 (x) と dual 解 (y) を取得する。
2. 問題データ（行列 (A)、ベクトル (b)、目的係数 (c)）を使って上記の式で slack と reduced cost を計算する。
3. 4つのベクトル（primal, dual, reduced cost, slack）を揃えて、COPT の setSolution() に渡す。

---

💡 まとめ

• cuPDLPX は 2種類しか返さない → 残りは 線形代数的に計算可能。
• slack = b - Ax
• reduced cost = c - Aᵀy
• これで COPT の要求する4つの Solution を満たせます。

---

孝幸さんは OSI 実装を触っているので、もし COPT OSI に cuPDLPX を接続するなら、この「補完計算」をラッパー側で実装するのが自然です。

👉 実際のコード例（C++/Python）で計算式を埋め込んだ形を見てみますか？

#include <vector>

#include <iostream>

// A: 制約行列 (m×n)、b: 右辺ベクトル (m)、c: 目的関数係数 (n)

// x: primal 解 (n)、y: dual 解 (m)

// slacks: b - A*x

// reducedCosts: c - A^T*y

void computeSlackAndReducedCosts(

    const std::vector<std::vector<double>>& A, // m×n

    const std::vector<double>& b,              // m

    const std::vector<double>& c,              // n

    const std::vector<double>& x,              // n (primal)

    const std::vector<double>& y,              // m (dual)

    std::vector<double>& slacks,               // m 出力

    std::vector<double>& reducedCosts          // n 出力

) {

    int m = A.size();

    int n = A[0].size();

    // slacks = b - A*x

    slacks.assign(m, 0.0);

    for (int i = 0; i < m; ++i) {

        double Ax_i = 0.0;

        for (int j = 0; j < n; ++j) {

            Ax_i += A[i][j] * x[j];

        }

        slacks[i] = b[i] - Ax_i;

    }

    // reducedCosts = c - A^T*y

    reducedCosts.assign(n, 0.0);

    for (int j = 0; j < n; ++j) {

        double ATy_j = 0.0;

        for (int i = 0; i < m; ++i) {

            ATy_j += A[i][j] * y[i];

        }

        reducedCosts[j] = c[j] - ATy_j;

    }

}

int main() {

    // 例: min c^T x, s.t. Ax <= b

    std::vector<std::vector<double>> A = {

        {1.0, 2.0},

        {3.0, 1.0}

    };

    std::vector<double> b = {4.0, 6.0};

    std::vector<double> c = {1.0, 1.0};

    // cuPDLPXから得られた解 (仮の値)

    std::vector<double> x = {1.0, 1.0}; // primal

    std::vector<double> y = {0.5, 0.2}; // dual

    std::vector<double> slacks, reducedCosts;

    computeSlackAndReducedCosts(A, b, c, x, y, slacks, reducedCosts);

    // 出力確認

    std::cout << "Slacks: ";

    for (double s : slacks) std::cout << s << " ";

    std::cout << "\nReduced Costs: ";

    for (double r : reducedCosts) std::cout << r << " ";

    std::cout << std::endl;

    return 0;

}