cpcert_bindings.cpp

// --------------------------------------------------------------------------
// pybind11 bindings for RankTools::CPCert and supporting types
// --------------------------------------------------------------------------
//
// The C++ CPCert class stores *references* to the constraint matrices
// (A) and RHS vector (b).  To make that safe from Python we introduce a small
// wrapper (PyCPCert) that *owns* copies of A and b and forwards every
// public method to the real object.
// --------------------------------------------------------------------------

#include <pybind11/eigen.h>  // automatic Eigen <-> numpy conversion
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>  // std::vector / std::pair conversion

#include "cpcert.hpp"
#include "interior_point_sdp.hpp"
#include "max_clique_certifier.hpp"
#include "rank_reduction.hpp"

namespace py = pybind11;
using namespace RankTools;

// ---------------------------------------------------------------------------
// Wrapper: owns copies of A and b so Python GC can't pull the rug out.
// ---------------------------------------------------------------------------
class PyCPCert {
 public:
  PyCPCert(const Matrix& C, double rho,
                   std::vector<Eigen::SparseMatrix<double>> A,
                   std::vector<double> b, CPCertParams params)
      : A_(std::move(A)),
        b_(std::move(b)),
        cpcert_(C, rho, A_, b_, std::move(params)) {}

  // ---- forwarded public API ----

  int dim() const { return cpcert_.dim; }
  int m() const { return cpcert_.m; }
  CPCertParams& params() { return cpcert_.params_; }
  const CPCertParams& params() const { return cpcert_.params_; }

  Vector eval_constraints(const Matrix& X) const {
    return cpcert_.eval_constraints(X);
  }

  SDPResult solve_sdp_mosek() const { return cpcert_.solve_sdp_mosek(); }

  CPCertResult certify(const Matrix& Y_0,
                               const Matrix* perturb = nullptr) const {
    return cpcert_.certify(Y_0, perturb);
  }

  std::pair<Matrix, Vector> get_central_path_point(
      const Matrix& Y_0, const Matrix* perturb = nullptr) const {
    return cpcert_.get_central_path_point(Y_0, perturb);
  }

  Matrix build_adjoint(const Vector& coeffs) const {
    return cpcert_.build_adjoint(coeffs);
  }

  Matrix build_certificate_from_dual(const Vector& multipliers) const {
    return build_adjoint(multipliers);
  }

  void export_problem(const std::string& file_path,
                      const std::string& problem_name,
                      const Matrix& solution) const {
    cpcert_.export_problem(file_path, problem_name, solution);
  }

  std::pair<double, double> eval_certificate(const Matrix& H,
                                             const Matrix& Y) const {
    return cpcert_.eval_certificate(H, Y);
  }

  std::pair<double, double> check_certificate(const Matrix& H,
                                              const Matrix& Y) const {
    return cpcert_.check_certificate(H, Y);
  }

 private:
  std::vector<Eigen::SparseMatrix<double>> A_;
  std::vector<double> b_;
  CPCert cpcert_;
};

// ---------------------------------------------------------------------------
// Wrapper: MaxCliqueCertifier builds (and owns copies of) its own constraint
// matrices from the cost matrix M, so unlike PyCPCert this wrapper does
// not need to hold A and b itself. It simply forwards the inherited public API.
// ---------------------------------------------------------------------------
class PyMaxCliqueCertifier {
 public:
  PyMaxCliqueCertifier(const Matrix& M, double rho, CPCertParams params)
      : mcc_(M, rho, std::move(params)) {}

  // ---- forwarded public API ----

  int dim() const { return mcc_.dim; }
  int m() const { return mcc_.m; }
  CPCertParams& params() { return mcc_.params_; }
  const CPCertParams& params() const { return mcc_.params_; }

  Vector eval_constraints(const Matrix& X) const {
    return mcc_.eval_constraints(X);
  }

  SDPResult solve_sdp_mosek() const { return mcc_.solve_sdp_mosek(); }

  CPCertResult certify(const Matrix& Y_0,
                               const Matrix* perturb = nullptr) const {
    return mcc_.certify(Y_0, perturb);
  }

  std::pair<Matrix, Vector> get_central_path_point(
      const Matrix& Y_0, const Matrix* perturb = nullptr) const {
    return mcc_.get_central_path_point(Y_0, perturb);
  }

  Matrix build_adjoint(const Vector& coeffs) const {
    return mcc_.build_adjoint(coeffs);
  }

  Matrix build_certificate_from_dual(const Vector& multipliers) const {
    return build_adjoint(multipliers);
  }

  void export_problem(const std::string& file_path,
                      const std::string& problem_name,
                      const Matrix& solution) const {
    mcc_.export_problem(file_path, problem_name, solution);
  }

  std::pair<double, double> eval_certificate(const Matrix& H,
                                             const Matrix& Y) const {
    return mcc_.eval_certificate(H, Y);
  }

  std::pair<double, double> check_certificate(const Matrix& H,
                                              const Matrix& Y) const {
    return mcc_.check_certificate(H, Y);
  }

 private:
  MaxCliqueCertifier mcc_;
};

// ---------------------------------------------------------------------------
// Module definition
// ---------------------------------------------------------------------------
PYBIND11_MODULE(ranktools, m) {
  m.doc() =
      "Python bindings for the RankTools library (CPCert, "
      "rank_reduction, solve_sdp_mosek)";

  // ---- LinearSolverType ----
  py::enum_<LinearSolverType>(
      m, "LinearSolverType", "Linear solver type for the CPCert step.")
      .value("LDLT", LinearSolverType::LDLT, "Cholesky-based LDLT solver")
      .value("CG", LinearSolverType::CG, "Conjugate gradient solver")
      .value(
          "MFCG_DP", LinearSolverType::MFCG_DP,
          "Matrix-free conjugate gradient solver with diagonal preconditioner")
      .value(
          "MFCG_LRP", LinearSolverType::MFCG_LRP,
          "Matrix-free conjugate gradient solver with low-rank preconditioner")
      .export_values();

  // ---- LowRankPrecondMethod ----
  py::enum_<LowRankPrecondMethod>(
      m, "LowRankPrecondMethod",
      "Method for building the low-rank preconditioner.")
      .value("DenseLDLT", LowRankPrecondMethod::DenseLDLT,
             "Dense LDLT factorization approach")
      .value("SparseLDLT", LowRankPrecondMethod::SparseLDLT,
             "Sparse LDLT factorization approach using alternate top-right "
             "formulation")
      .value("SparseLDLT_ZL", LowRankPrecondMethod::SparseLDLT_ZL,
             "Sparse LDLT factorization approach using the Zhang-Lavaei (2017) "
             "top-right formulation")
      .value("DenseQR", LowRankPrecondMethod::DenseQR,
             "Dense QR factorization approach")
      .value("SparseQR", LowRankPrecondMethod::SparseQR,
             "Sparse QR factorization approach")
      .value("DenseLU", LowRankPrecondMethod::DenseLU,
             "Dense LU factorization approach")
      .value("DirectInverse", LowRankPrecondMethod::DirectInverse,
             "Direct inverse approach")
      .export_values();

  // ---- LowRankPrecondParams ----
  py::class_<LowRankPrecondParams>(m, "LowRankPrecondParams")
      .def(py::init<>())
      .def_readwrite("tau", &LowRankPrecondParams::tau)
      .def_readwrite("method", &LowRankPrecondParams::method)
      .def_readwrite("use_approx", &LowRankPrecondParams::use_approx)
      .def_readwrite("ldlt_zero_thresh",
                     &LowRankPrecondParams::ldlt_zero_thresh)
      .def("__repr__", [](const LowRankPrecondParams& p) {
        return "<LowRankPrecondParams tau=" + std::to_string(p.tau) +
               " method=" + std::to_string(static_cast<int>(p.method)) + ">";
      });

  // ---- CPCertParams ----
  py::class_<CPCertParams>(m, "CPCertParams")
      .def(py::init<>())
      // General
      .def_readwrite("verbose", &CPCertParams::verbose)
      .def_readwrite("tol_rank_sol", &CPCertParams::tol_rank_sol)
      .def_readwrite("tol_step_norm", &CPCertParams::tol_step_norm)
      .def_readwrite("max_iter", &CPCertParams::max_iter)
      .def_readwrite("rescale_lin_sys", &CPCertParams::rescale_lin_sys)
      .def_readwrite("rescaling_factor",
                     &CPCertParams::rescaling_factor)
      .def_readwrite("reuse_multipliers",
                     &CPCertParams::reuse_multipliers)
      .def_readwrite("check_indep_constr",
                     &CPCertParams::check_indep_constr)
      .def_readwrite("tol_indep_constr",
                     &CPCertParams::tol_indep_constr)
      .def_readwrite("delta", &CPCertParams::delta)
      .def_readwrite("perturb_constraints",
                     &CPCertParams::perturb_constraints)
      .def_readwrite("perturb_cost", &CPCertParams::perturb_cost)
      .def_readwrite("eps_cost", &CPCertParams::eps_cost)
      .def_readwrite("eps_constr", &CPCertParams::eps_constr)
      .def_readwrite("adaptive_perturb",
                     &CPCertParams::adaptive_perturb)
      .def_readwrite("eps_mult_min", &CPCertParams::eps_mult_min)
      .def_readwrite("eps_inc_step_thresh",
                     &CPCertParams::eps_inc_step_thresh)
      .def_readwrite("eps_inc", &CPCertParams::eps_inc)
      .def_readwrite("eps_dec_step_thresh",
                     &CPCertParams::eps_dec_step_thresh)
      .def_readwrite("eps_dec", &CPCertParams::eps_dec)
      .def_readwrite("lin_solver", &CPCertParams::lin_solver)
      // Iterative linear solve
      .def_readwrite("lin_solve_max_iter",
                     &CPCertParams::lin_solve_max_iter)
      .def_readwrite("lin_solve_tol", &CPCertParams::lin_solve_tol)
      .def_property(
          "lrp_params",
          [](CPCertParams& p) -> LowRankPrecondParams& {
            return p.lrp_params;
          },
          [](CPCertParams& p, const LowRankPrecondParams& lrp) {
            p.lrp_params = lrp;
          },
          py::return_value_policy::reference_internal)
      // Backward-compatible alias for older scripts.
      .def_property(
          "tau_lrp",
          [](const CPCertParams& p) { return p.lrp_params.tau; },
          [](CPCertParams& p, double tau) { p.lrp_params.tau = tau; })
      // Line search
      .def_readwrite("enable_line_search",
                     &CPCertParams::enable_line_search)
      .def_readwrite("ln_search_red_factor",
                     &CPCertParams::ln_search_red_factor)
      .def_readwrite("alpha_init", &CPCertParams::alpha_init)
      .def_readwrite("alpha_min", &CPCertParams::alpha_min)
      // Certificate
      .def_readwrite("early_stop_cert", &CPCertParams::early_stop_cert)
      .def_readwrite("tol_cert_psd", &CPCertParams::tol_cert_psd)
      .def_readwrite("tol_cert_complementarity",
                     &CPCertParams::tol_cert_complementarity)
      .def_readwrite("tol_cert_primal_feas",
                     &CPCertParams::tol_cert_primal_feas)
      .def_readwrite("early_stop_angle",
                     &CPCertParams::early_stop_angle)
      .def_readwrite("max_angle", &CPCertParams::max_angle)
      .def_readwrite("use_cert_centrality_metric",
                     &CPCertParams::use_cert_centrality_metric)
      .def_readwrite("tol_cert_centrality",
                     &CPCertParams::tol_cert_centrality)
      .def("__repr__", [](const CPCertParams& p) {
        return "<CPCertParams max_iter=" + std::to_string(p.max_iter) +
               " verbose=" + (p.verbose ? "True" : "False") + ">";
      });

  // ---- CPCertResult ----
  py::class_<CPCertResult>(m, "CPCertResult")
      .def_readonly("X", &CPCertResult::X)
      .def_readonly("H", &CPCertResult::H)
      .def_readonly("multipliers", &CPCertResult::multipliers)
      .def_readonly("violation", &CPCertResult::violation)
      .def_readonly("certified", &CPCertResult::certified)
      .def_readonly("min_eig", &CPCertResult::min_eig)
      .def_readonly("complementarity", &CPCertResult::complementarity)
      .def_readonly("solver_time", &CPCertResult::solver_time)
      .def_readonly("num_iterations", &CPCertResult::num_iterations)
      .def("__repr__", [](const CPCertResult& r) {
        return "<CPCertResult certified=" +
               std::string(r.certified ? "True" : "False") +
               " min_eig=" + std::to_string(r.min_eig) + ">";
      });

  // ---- CPCert (via PyCPCert wrapper) ----
  py::class_<PyCPCert>(m, "CPCert")
      .def(py::init<const Matrix&, double,
                    std::vector<Eigen::SparseMatrix<double>>,
                    std::vector<double>, CPCertParams>(),
           py::arg("C"), py::arg("rho"), py::arg("A"), py::arg("b"),
           py::arg("params") = CPCertParams(),
           R"pbdoc(
Construct an CPCert problem.

Parameters
----------
C : numpy.ndarray (n, n)
    Cost matrix.
rho : float
    Optimal cost value (scalar offset).
A : list of scipy.sparse.csc_matrix or numpy.ndarray
    Constraint matrices (each n×n, upper-triangular storage).
b : list of float
    Right-hand side values for trace(A_i X) = b_i.
params : CPCertParams, optional
    Algorithm parameters.
)pbdoc")
      .def_property_readonly("dim", &PyCPCert::dim)
      .def_property_readonly("m", &PyCPCert::m)
      .def_property(
          "params", py::overload_cast<>(&PyCPCert::params),
          [](PyCPCert& self, const CPCertParams& p) {
            self.params() = p;
          },
          py::return_value_policy::reference_internal)
      .def("eval_constraints", &PyCPCert::eval_constraints,
           py::arg("X"), "Evaluate constraint violations at X.")
      .def("solve_sdp_mosek", &PyCPCert::solve_sdp_mosek,
           R"pbdoc(
Solve the problem's SDP using MOSEK.

Builds the primal SDP from this object's cost matrix and constraints and
solves it with MOSEK.

Returns
-------
SDPResult
)pbdoc")
      .def(
          "certify",
          [](const PyCPCert& self, const Matrix& Y_0,
             const py::object& perturb = py::none()) {
            if (perturb.is_none()) {
              return self.certify(Y_0, nullptr);
            } else {
              auto perturb_mat = perturb.cast<Matrix>();
              return self.certify(Y_0, &perturb_mat);
            }
          },
          py::arg("Y_0"), py::arg("perturb") = py::none(),
          R"pbdoc(
Run CPCert to certify the local solution Y_0.

If perturb is provided, it is used as the initial perturbation matrix.
Otherwise, params.delta * Identity is used as the fallback.

Parameters
----------
Y_0 : numpy.ndarray (n, r)
    Initial low-rank factor.
perturb : numpy.ndarray (n, n), optional
    Initial perturbation matrix. If not provided, uses params.delta * Identity.

Returns
-------
CPCertResult
)pbdoc")
      .def(
          "get_central_path_point",
          [](const PyCPCert& self, const Matrix& Y_0,
             const py::object& perturb = py::none()) {
            if (perturb.is_none()) {
              return self.get_central_path_point(Y_0, nullptr);
            } else {
              auto perturb_mat = perturb.cast<Matrix>();
              return self.get_central_path_point(Y_0, &perturb_mat);
            }
          },
          py::arg("Y_0"), py::arg("perturb") = py::none(),
          R"pbdoc(
Compute the central path point starting from Y_0.

If perturb is provided, it is used as the initial perturbation matrix.
Otherwise, params.delta * Identity is used as the fallback.

Parameters
----------
Y_0 : numpy.ndarray (n, r)
    Initial point (low-rank factor; X_0 = Y_0 @ Y_0.T).
perturb : numpy.ndarray (n, n), optional
    Initial perturbation matrix. If not provided, uses params.delta * Identity.

Returns
-------
tuple(X, multipliers)
    X : numpy.ndarray (n, n)  — centered primal solution.
    multipliers : numpy.ndarray (m,) — optimal dual multipliers.
)pbdoc")
      .def("build_adjoint", &PyCPCert::build_adjoint, py::arg("coeffs"),
           "Build the adjoint matrix sum_i coeffs[i] * A_i + coeffs[-1] * C.")
      .def("build_certificate_from_dual",
           &PyCPCert::build_certificate_from_dual,
           py::arg("multipliers"),
           "Backward-compatible alias for build_adjoint(multipliers).")
      .def("export_problem", &PyCPCert::export_problem,
           py::arg("file_path"), py::arg("problem_name"), py::arg("solution"),
           R"pbdoc(
Export the current problem to a text file.

The file format matches `load_problem_from_file` in the C++ test helpers.

Parameters
----------
file_path : str
    Destination file path.
problem_name : str
    Value written to the `name` field.
solution : numpy.ndarray
    Solution matrix written in the `soln` block.
)pbdoc")
      .def("eval_certificate", &PyCPCert::eval_certificate,
           py::arg("H"), py::arg("Y"),
           R"pbdoc(
Evaluate the optimality certificate.

Returns the minimum eigenvalue of the certificate matrix and the
evaluation of the certificate matrix at the solution (first order
condition).

Parameters
----------
H : numpy.ndarray (n, n)
    Certificate matrix.
Y : numpy.ndarray (n, r)
    Low-rank factor of the solution.

Returns
-------
tuple(min_eig, first_order_cond)
)pbdoc")
      .def("check_certificate", &PyCPCert::check_certificate,
           py::arg("H"), py::arg("Y"),
           R"pbdoc(
Check global optimality of a solution.

Returns whether the certificate matrix is PSD and the complementarity
of the provided solution.

Returns
-------
tuple(min_eig, complementarity)
)pbdoc");

  // ---- MaxCliqueCertifier (via PyMaxCliqueCertifier wrapper) ----
  py::class_<PyMaxCliqueCertifier>(m, "MaxCliqueCertifier")
      .def(py::init<const Matrix&, double, CPCertParams>(),
           py::arg("M"), py::arg("rho"),
           py::arg("params") = CPCertParams(),
           R"pbdoc(
Construct a MaxCliqueCertifier for the maximum-clique / Lovasz-theta SDP.

Behaves exactly like CPCert, except that the constraint matrices A and
right-hand side b are not supplied by the caller. They are constructed from the
cost matrix M: there is one constraint per non-edge enforcing that the
corresponding off-diagonal entry of the solution is zero, plus a trace
constraint fixing the trace to one. The non-edges are taken to be the
off-diagonal (i, j) indices of M whose entries are exactly equal to zero.

Parameters
----------
M : numpy.ndarray (n, n)
    Cost matrix. Its off-diagonal zero entries define the non-edge constraints.
rho : float
    Optimal cost value (scalar offset).
params : CPCertParams, optional
    Algorithm parameters.
)pbdoc")
      .def_property_readonly("dim", &PyMaxCliqueCertifier::dim)
      .def_property_readonly("m", &PyMaxCliqueCertifier::m)
      .def_property(
          "params", py::overload_cast<>(&PyMaxCliqueCertifier::params),
          [](PyMaxCliqueCertifier& self, const CPCertParams& p) {
            self.params() = p;
          },
          py::return_value_policy::reference_internal)
      .def("eval_constraints", &PyMaxCliqueCertifier::eval_constraints,
           py::arg("X"), "Evaluate constraint violations at X.")
      .def("solve_sdp_mosek", &PyMaxCliqueCertifier::solve_sdp_mosek,
           R"pbdoc(
Solve the problem's SDP using MOSEK.

Builds the primal SDP from this object's cost matrix and constraints and
solves it with MOSEK.

Returns
-------
SDPResult
)pbdoc")
      .def(
          "certify",
          [](const PyMaxCliqueCertifier& self, const Matrix& Y_0,
             const py::object& perturb = py::none()) {
            if (perturb.is_none()) {
              return self.certify(Y_0, nullptr);
            } else {
              auto perturb_mat = perturb.cast<Matrix>();
              return self.certify(Y_0, &perturb_mat);
            }
          },
          py::arg("Y_0"), py::arg("perturb") = py::none(),
          R"pbdoc(
Run CPCert to certify the local solution Y_0.

If perturb is provided, it is used as the initial perturbation matrix.
Otherwise, params.delta * Identity is used as the fallback.

Parameters
----------
Y_0 : numpy.ndarray (n, r)
    Initial low-rank factor.
perturb : numpy.ndarray (n, n), optional
    Initial perturbation matrix. If not provided, uses params.delta * Identity.

Returns
-------
CPCertResult
)pbdoc")
      .def(
          "get_central_path_point",
          [](const PyMaxCliqueCertifier& self, const Matrix& Y_0,
             const py::object& perturb = py::none()) {
            if (perturb.is_none()) {
              return self.get_central_path_point(Y_0, nullptr);
            } else {
              auto perturb_mat = perturb.cast<Matrix>();
              return self.get_central_path_point(Y_0, &perturb_mat);
            }
          },
          py::arg("Y_0"), py::arg("perturb") = py::none(),
          R"pbdoc(
Compute the central path point starting from Y_0.

If perturb is provided, it is used as the initial perturbation matrix.
Otherwise, params.delta * Identity is used as the fallback.

Parameters
----------
Y_0 : numpy.ndarray (n, r)
    Initial point (low-rank factor; X_0 = Y_0 @ Y_0.T).
perturb : numpy.ndarray (n, n), optional
    Initial perturbation matrix. If not provided, uses params.delta * Identity.

Returns
-------
tuple(X, multipliers)
    X : numpy.ndarray (n, n)  — centered primal solution.
    multipliers : numpy.ndarray (m,) — optimal dual multipliers.
)pbdoc")
      .def("build_adjoint", &PyMaxCliqueCertifier::build_adjoint,
           py::arg("coeffs"),
           "Build the adjoint matrix sum_i coeffs[i] * A_i + coeffs[-1] * C.")
      .def("build_certificate_from_dual",
           &PyMaxCliqueCertifier::build_certificate_from_dual,
           py::arg("multipliers"),
           "Backward-compatible alias for build_adjoint(multipliers).")
      .def("export_problem", &PyMaxCliqueCertifier::export_problem,
           py::arg("file_path"), py::arg("problem_name"), py::arg("solution"),
           R"pbdoc(
Export the current problem to a text file.

The file format matches `load_problem_from_file` in the C++ test helpers.

Parameters
----------
file_path : str
    Destination file path.
problem_name : str
    Value written to the `name` field.
solution : numpy.ndarray
    Solution matrix written in the `soln` block.
)pbdoc")
      .def("eval_certificate", &PyMaxCliqueCertifier::eval_certificate,
           py::arg("H"), py::arg("Y"),
           R"pbdoc(
Evaluate the optimality certificate.

Returns the minimum eigenvalue of the certificate matrix and the
evaluation of the certificate matrix at the solution (first order
condition).

Parameters
----------
H : numpy.ndarray (n, n)
    Certificate matrix.
Y : numpy.ndarray (n, r)
    Low-rank factor of the solution.

Returns
-------
tuple(min_eig, first_order_cond)
)pbdoc")
      .def("check_certificate", &PyMaxCliqueCertifier::check_certificate,
           py::arg("H"), py::arg("Y"),
           R"pbdoc(
Check global optimality of a solution.

Returns whether the certificate matrix is PSD and the complementarity
of the provided solution.

Returns
-------
tuple(min_eig, complementarity)
)pbdoc");

  // ---- SDPResult ----
  py::class_<SDPResult>(m, "SDPResult")
      .def_readonly("X", &SDPResult::X)
      .def_readonly("y", &SDPResult::y)
      .def_readonly("S", &SDPResult::S)
      .def_readonly("obj_value", &SDPResult::obj_value)
      .def_readonly("solver_time", &SDPResult::solver_time)
      .def("__repr__", [](const SDPResult& r) {
        return "<SDPResult obj_value=" + std::to_string(r.obj_value) + ">";
      });

  // ---- solve_sdp_mosek ----
  m.def(
      "solve_sdp_mosek",
      [](const Eigen::MatrixXd& C, std::vector<Eigen::SparseMatrix<double>> As,
         std::vector<double> b,
         bool verbose) { return solve_sdp_mosek(C, As, b, verbose); },
      py::arg("C"), py::arg("As"), py::arg("b"), py::arg("verbose") = true,
      R"pbdoc(
Solve a semidefinite program using MOSEK.

Solves the primal SDP:
    min   trace(C @ X)
    s.t.  trace(A_i @ X) = b_i  for i = 0 ... m-1
          X >= 0  (positive semidefinite)

Parameters
----------
C : numpy.ndarray (n, n)
    Cost matrix.
As : list of scipy.sparse matrices
    Constraint matrices (each n x n).
b : list of float
    Right-hand side values.
verbose : bool, optional
    Enable MOSEK solver output (default True).

Returns
-------
SDPResult
    Struct with fields:
      X         : numpy.ndarray (n, n) — primal solution.
      y         : numpy.ndarray (m,)   — dual multipliers for equality constraints.
      S         : numpy.ndarray (n, n) — dual PSD matrix C - sum_i y_i A_i.
      obj_value : float                — objective value trace(C @ X).
      solver_time : float                — time taken by the solver (seconds).
)pbdoc");

  // ---- RankReductionParams ----
  py::class_<RankReductionParams>(m, "RankReductionParams")
      .def(py::init<>())
      .def_readwrite("targ_rank", &RankReductionParams::targ_rank)
      .def_readwrite("null_tol", &RankReductionParams::null_tol)
      .def_readwrite("eig_tol", &RankReductionParams::eig_tol)
      .def_readwrite("max_iter", &RankReductionParams::max_iter)
      .def_readwrite("verbose", &RankReductionParams::verbose)
      .def("__repr__", [](const RankReductionParams& p) {
        return "<RankReductionParams targ_rank=" + std::to_string(p.targ_rank) +
               " verbose=" + (p.verbose ? "True" : "False") + ">";
      });

  // ---- rank_reduction ----
  m.def(
      "rank_reduction",
      [](std::vector<Eigen::SparseMatrix<double>> As, const Matrix& V_init,
         RankReductionParams params) {
        return rank_reduction(As, V_init, std::move(params));
      },
      py::arg("As"), py::arg("V_init"),
      py::arg("params") = RankReductionParams(),
      R"pbdoc(
Reduce the rank of an SDP solution.

Implements the rank reduction algorithm from:
  Lemon, So & Ye, "Low-rank semidefinite programming: Theory and
  applications", Foundations and Trends in Optimization, 2016.

Parameters
----------
As : list of scipy.sparse matrices
    Constraint matrices (each n×n, upper-triangular storage).
V_init : numpy.ndarray (n, r)
    Initial low-rank factor; the SDP solution is X = V_init @ V_init.T.
params : RankReductionParams, optional
    Algorithm parameters.

Returns
-------
numpy.ndarray (n, r')
    Reduced low-rank factor V such that V @ V.T satisfies all constraints
    with rank r' <= r.
)pbdoc");
}