Support general low-rank object format RDR^T and add middle D test

README.md

@@ -12,16 +12,20 @@ For a detailed explanation of the methodology, please refer to our papers: [A Re
 
 ## Problem Formulation
 
-PDHCG solves convex quadratic programs in the following form, which allows a flexibile input of convex quadratic objective matrix, a sparse component and a low-rank component:
+PDHCG solves convex quadratic programs in the following form, which allows a flexible input of the quadratic objective matrix as a sparse component plus a structured low-rank component:
 
 ```math
 \begin{aligned}
-\min_{x} \quad & \frac{1}{2}x^\top (Q + R^\top R) x + c^\top x \\
+\min_{x} \quad & \frac{1}{2}x^\top (Q + R^\top D R) x + c^\top x \\
 \text{s.t.} \quad & \ell_c \le Ax \le u_c, \\
                   & \ell_v \le x \le u_v.
 \end{aligned}
 ```
 
+- $Q$ is the sparse symmetric quadratic component (optional).
+- $R \in \mathbb{R}^{k\times n}$ is a tall low-rank factor (optional, $k$ = rank).
+- $D \in \mathbb{R}^{k\times k}$ is an optional middle matrix that scales / weights / signs the low-rank term. When omitted it defaults to the identity, recovering the standard $Q + R^\top R$ formulation. $D$ may be **diagonal, sparse, dense, or indefinite** — the backend auto-detects the cheapest runtime representation.
+
 
 ## Installation (C++ Executable)
 
@@ -139,8 +143,9 @@ import numpy as np
 import scipy.sparse as sp
 from pdhcg import Model
 
-# Example: minimize 0.5 * x'(Q + R'R)x + c'x
+# Example: minimize 0.5 * x'(Q + R^T D R)x + c'x
 # subject to l <= A x <= u,  lb <= x <= ub
+# (D defaults to identity, recovering the classic Q + R^T R form.)
 
 # 1. Define Standard QP terms
 Q = sp.csc_matrix([[1.0, -1.0], [-1.0, 2.0]])
@@ -151,16 +156,22 @@ c = np.array([-2.0, -6.0])
 # This adds 0.5 * (x1)^2 to the objective
 R = sp.csc_matrix([[1.0, 0.0]])
 
-# 3. Define Constraints
+# 3. (Optional) Middle matrix D in 0.5 * x^T R^T D R x. Pass a 1-D array
+# for a diagonal D, a 2-D array for dense D, or a scipy.sparse matrix.
+# Omit entirely (or pass None) to use D = identity.
+# D = np.array([2.5])  # e.g., weight the low-rank term by 2.5
+
+# 4. Define Constraints
 A = sp.csc_matrix([[1.0, 1.0], [-1.0, 2.0], [2.0, 1.0]])
 l = np.array([-np.inf, -np.inf, -np.inf])
 u = np.array([2.0, 2.0, 3.0])
 lb = np.zeros(2)
 ub = np.array([np.inf, np.inf])
 
-# 4. Create QP model with Low-Rank term (R), where Q and R are both optional.
+# 5. Create QP model. Q, R and D are all optional.
 m = Model(objective_matrix=Q,
           objective_matrix_low_rank=R,
+          # objective_matrix_low_rank_middle=D,  # uncomment to use D
           objective_vector=c,
           constraint_matrix=A,
           constraint_lower_bound=l,

distributed/distributed_utils.cu

@@ -328,6 +328,29 @@ qp_problem_t *partition_qp_problem(const qp_problem_t *global_qp,
                                                           global_qp->objective_lowrank_matrix,
                                                           &loc->objective_lowrank_matrix_num_nonzeros);
 
+    loc->objective_lowrank_middle_matrix = NULL;
+    loc->objective_lowrank_middle_matrix_num_nonzeros = 0;
+    if (global_qp->objective_lowrank_middle_matrix && loc->num_rank_lowrank_obj > 0 &&
+        global_qp->objective_lowrank_middle_matrix_num_nonzeros > 0)
+    {
+        int k = loc->num_rank_lowrank_obj;
+        int nnz = global_qp->objective_lowrank_middle_matrix_num_nonzeros;
+        loc->objective_lowrank_middle_matrix_num_nonzeros = nnz;
+        loc->objective_lowrank_middle_matrix = (CsrComponent *)malloc(sizeof(CsrComponent));
+        loc->objective_lowrank_middle_matrix->row_ptr = (int *)malloc((size_t)(k + 1) * sizeof(int));
+        loc->objective_lowrank_middle_matrix->col_ind = (int *)malloc((size_t)nnz * sizeof(int));
+        loc->objective_lowrank_middle_matrix->val = (double *)malloc((size_t)nnz * sizeof(double));
+        memcpy(loc->objective_lowrank_middle_matrix->row_ptr,
+               global_qp->objective_lowrank_middle_matrix->row_ptr,
+               (size_t)(k + 1) * sizeof(int));
+        memcpy(loc->objective_lowrank_middle_matrix->col_ind,
+               global_qp->objective_lowrank_middle_matrix->col_ind,
+               (size_t)nnz * sizeof(int));
+        memcpy(loc->objective_lowrank_middle_matrix->val,
+               global_qp->objective_lowrank_middle_matrix->val,
+               (size_t)nnz * sizeof(double));
+    }
+
     loc->objective_vector = copy_slice(global_qp->objective_vector, n_start, loc->num_variables);
     loc->variable_lower_bound = copy_slice(global_qp->variable_lower_bound, n_start, loc->num_variables);
     loc->variable_upper_bound = copy_slice(global_qp->variable_upper_bound, n_start, loc->num_variables);
@@ -391,6 +414,24 @@ rescale_info_t *partition_rescale_info(rescale_info_t *global_info,
     loc_processed->objective_lowrank_matrix = loc_lp->objective_lowrank_matrix;
     loc_processed->quad_type = global_processed->quad_type;
 
+    loc_processed->objective_lowrank_middle_kind = global_processed->objective_lowrank_middle_kind;
+    loc_processed->objective_lowrank_middle_diag = NULL;
+    loc_processed->objective_lowrank_middle_dense = NULL;
+    if (global_processed->objective_lowrank_middle_kind == PDHCG_D_DIAG &&
+        global_processed->objective_lowrank_middle_diag && loc_lp->num_rank_lowrank_obj > 0)
+    {
+        size_t bytes = (size_t)loc_lp->num_rank_lowrank_obj * sizeof(double);
+        loc_processed->objective_lowrank_middle_diag = (double *)malloc(bytes);
+        memcpy(loc_processed->objective_lowrank_middle_diag, global_processed->objective_lowrank_middle_diag, bytes);
+    }
+    else if (global_processed->objective_lowrank_middle_kind == PDHCG_D_DENSE &&
+             global_processed->objective_lowrank_middle_dense && loc_lp->num_rank_lowrank_obj > 0)
+    {
+        size_t bytes = (size_t)loc_lp->num_rank_lowrank_obj * (size_t)loc_lp->num_rank_lowrank_obj * sizeof(double);
+        loc_processed->objective_lowrank_middle_dense = (double *)malloc(bytes);
+        memcpy(loc_processed->objective_lowrank_middle_dense, global_processed->objective_lowrank_middle_dense, bytes);
+    }
+
     if (global_processed->quad_type == PDHCG_DIAG_Q && global_processed->diagonal_quad_objective != NULL)
     {
         loc_processed->diagonal_quad_objective =
@@ -433,6 +474,11 @@ size_t get_qp_problem_size(const qp_problem_t *qp)
     ADD_CSR_SIZE(qp->objective_sparse_matrix, qp->num_variables, qp->objective_sparse_matrix_num_nonzeros);
     ADD_CSR_SIZE(qp->objective_lowrank_matrix, qp->num_rank_lowrank_obj, qp->objective_lowrank_matrix_num_nonzeros);
 
+    size += sizeof(int);
+    ADD_CSR_SIZE(qp->objective_lowrank_middle_matrix,
+                 qp->num_rank_lowrank_obj,
+                 qp->objective_lowrank_middle_matrix_num_nonzeros);
+
     size += sizeof(int) * 2;
     if (qp->primal_start)
         size += sizeof(double) * qp->num_variables;
@@ -487,6 +533,11 @@ void serialize_qp_problem_to_ptr(const qp_problem_t *qp, char **ptr_ref)
     S_CSR(qp->objective_sparse_matrix, qp->num_variables, qp->objective_sparse_matrix_num_nonzeros);
     S_CSR(qp->objective_lowrank_matrix, qp->num_rank_lowrank_obj, qp->objective_lowrank_matrix_num_nonzeros);
 
+    S_COPY(qp->objective_lowrank_middle_matrix_num_nonzeros, int);
+    S_CSR(qp->objective_lowrank_middle_matrix,
+          qp->num_rank_lowrank_obj,
+          qp->objective_lowrank_middle_matrix_num_nonzeros);
+
     int has_primal = (qp->primal_start != NULL);
     int has_dual = (qp->dual_start != NULL);
     S_COPY(has_primal, int);
@@ -547,6 +598,12 @@ qp_problem_t *deserialize_qp_problem_from_ptr(const char **ptr_ref)
     D_CSR(qp->objective_sparse_matrix, qp->num_variables, qp->objective_sparse_matrix_num_nonzeros);
     D_CSR(qp->objective_lowrank_matrix, qp->num_rank_lowrank_obj, qp->objective_lowrank_matrix_num_nonzeros);
 
+    D_VAL(qp->objective_lowrank_middle_matrix_num_nonzeros, int);
+    qp->objective_lowrank_middle_matrix = NULL;
+    D_CSR(qp->objective_lowrank_middle_matrix,
+          qp->num_rank_lowrank_obj,
+          qp->objective_lowrank_middle_matrix_num_nonzeros);
+
     int has_primal, has_dual;
     D_VAL(has_primal, int);
     D_VAL(has_dual, int);

docs/C_API.md

@@ -10,9 +10,10 @@ The C API involves two main functions:
 ```c
 qp_problem_t *create_qp_problem(
     const double *objective_c,           // objective vector (length n)
-    const matrix_desc_t *Q_desc,         // quadratic sparse matrix (n×n)
-    const matrix_desc_t *R_desc,         // quadratic low-rank matrix (n×m)
-    const matrix_desc_t *A_desc,         // constraint matrix (m×n)
+    const matrix_desc_t *Q_desc,         // sparse quadratic matrix (n x n)
+    const matrix_desc_t *R_desc,         // low-rank factor (k x n)
+    const matrix_desc_t *D_desc,         // middle matrix in R^T D R (k x k)
+    const matrix_desc_t *A_desc,         // constraint matrix (m x n)
     const double *con_lb,                // constraint lower bounds (length m)
     const double *con_ub,                // constraint upper bounds (length m)
     const double *var_lb,                // variable lower bounds (length n)
@@ -26,10 +27,13 @@ pdhcg_result_t* solve_qp_problem(
 );
 ```
 
+The objective minimized is `0.5 * x^T (Q + R^T D R) x + c^T x + c0`. `Q`, `R`, and `D` are all optional; `D` defaults to identity, recovering the standard `Q + R^T R` form.
+
 `create_qp_problem` parameters:
 - `objective_c`: Objective vector. If `NULL`, defaults to all zeros.
-- `Q_desc`: Matrix descriptor. Supports `matrix_dense`, `matrix_csr`, `matrix_csc`, `matrix_coo`.
-- `R_desc`: Matrix descriptor. Supports `matrix_dense`, `matrix_csr`, `matrix_csc`, `matrix_coo`.
+- `Q_desc`: Matrix descriptor. Supports `matrix_dense`, `matrix_csr`, `matrix_csc`, `matrix_coo`. Pass `NULL` to omit.
+- `R_desc`: Matrix descriptor for the low-rank factor (shape `k x n`). Same supported formats. Pass `NULL` to omit.
+- `D_desc`: Matrix descriptor for the middle matrix in `R^T D R` (shape `k x k`). Same supported formats. May be diagonal, sparse, dense, or indefinite — the runtime auto-detects the cheapest representation in `preprocess_qp_problem`. Pass `NULL` for `D = I`.
 - `A_desc`: Matrix descriptor. Supports `matrix_dense`, `matrix_csr`, `matrix_csc`, `matrix_coo`.
 - `con_lb`: Constraint lower bounds. If `NULL`, defaults to all `-INFINITY`.
 - `con_ub`: Constraint upper bounds. If `NULL`, defaults to all `+INFINITY`.
@@ -103,11 +107,13 @@ int main() {
     double ub[2] = {INFINITY, INFINITY};
 
     // 6. Build the QP problem
-    // Note: We pass NULL for R_desc (low-rank factor) and objective_constant
+    // Note: We pass NULL for R_desc (low-rank factor), D_desc (middle matrix),
+    // and objective_constant.
     qp_problem_t* prob = create_qp_problem(
         c,          // objective_c
         &Q_desc,    // Q_desc
         NULL,       // R_desc
+        NULL,       // D_desc (NULL -> D = I)
         &A_desc,    // A_desc
         l,          // con_lb
         u,          // con_ub

docs/c/functions.md

@@ -7,22 +7,25 @@ qp_problem_t *create_qp_problem(
     const double *objective_c,
     const matrix_desc_t *Q_desc,
     const matrix_desc_t *R_desc,
+    const matrix_desc_t *D_desc,
     const matrix_desc_t *A_desc,
     const double *con_lb, const double *con_ub,
     const double *var_lb, const double *var_ub,
     const double *objective_constant
 );
 ```
 
-Creates a QP problem from matrix descriptors.
+Creates a QP problem of the form
+`min 0.5 * x^T (Q + R^T D R) x + c^T x` subject to `con_lb <= A x <= con_ub` and `var_lb <= x <= var_ub`.
 
 **Parameters:**
 
 | Parameter | Description |
 |-----------|-------------|
 | `objective_c` | Linear objective coefficients (size n) |
 | `Q_desc` | Sparse quadratic matrix descriptor (can be NULL) |
-| `R_desc` | Low-rank quadratic matrix descriptor (can be NULL) |
+| `R_desc` | Low-rank factor descriptor, shape `k x n` (can be NULL) |
+| `D_desc` | Middle matrix in `R^T D R`, shape `k x k` (can be NULL). |
 | `A_desc` | Constraint matrix descriptor |
 | `con_lb` | Constraint lower bounds (size m) |
 | `con_ub` | Constraint upper bounds (size m) |

docs/c/overview.md

@@ -46,9 +46,9 @@ int main() {
     double var_lb[] = {0.0, 0.0};
     double var_ub[] = {1e30, 1e30};
 
-    // Create problem
+    // Create problem (NULL for Q, R, and D -> linear problem)
     qp_problem_t *prob = create_qp_problem(
-        c, NULL, NULL, &A_desc,
+        c, NULL, NULL, NULL, &A_desc,
         con_lb, con_ub, var_lb, var_ub, NULL
     );
 
@@ -95,14 +95,17 @@ qp_problem_t *create_qp_problem(
     const double *objective_c,
     const matrix_desc_t *Q_desc,
     const matrix_desc_t *R_desc,
+    const matrix_desc_t *D_desc,
     const matrix_desc_t *A_desc,
     const double *con_lb, const double *con_ub,
     const double *var_lb, const double *var_ub,
     const double *objective_constant
 );
 ```
 
-Creates a QP problem from matrix descriptors. The `Q_desc` (sparse quadratic) and `R_desc` (low-rank quadratic) are optional (pass `NULL` if not needed).
+Creates a QP problem of the form
+`min 0.5 * x^T (Q + R^T D R) x + c^T x  s.t.  con_lb <= A x <= con_ub,  var_lb <= x <= var_ub`
+from matrix descriptors. `Q_desc` (sparse quadratic), `R_desc` (low-rank factor, shape `k x n`), and `D_desc` (rank-by-rank middle matrix in `R^T D R`) are all optional — pass `NULL` to omit any of them. `D_desc` defaults to identity, recovering the standard `Q + R^T R` formulation; it may be diagonal, sparse, dense, or indefinite, and the runtime auto-detects the cheapest representation.
 
 ### Setting Start Values
 

docs/c/types.md

@@ -129,6 +129,9 @@ typedef struct {
   CsrComponent *objective_lowrank_matrix;
   int objective_lowrank_matrix_num_nonzeros;
 
+  CsrComponent *objective_lowrank_middle_matrix;
+  int objective_lowrank_middle_matrix_num_nonzeros;
+
   double *constraint_lower_bound;
   double *constraint_upper_bound;
 

docs/examples.md

@@ -59,6 +59,45 @@ m = Model(
 m.optimize()
 ```
 
+### Low-Rank with a Middle Matrix D
+
+```python
+import numpy as np
+import scipy.sparse as sp
+from pdhcg import Model
+
+# Minimize 0.5 * x^T (Q + R^T D R) x + c^T x
+# D may be diagonal (1-D), dense (2-D), or scipy.sparse;
+# may also be indefinite. Defaults to identity if omitted.
+
+n = 1000
+r = 10  # rank
+Q = None
+R = np.random.randn(r, n)
+c = np.random.randn(n)
+
+# (a) Weighted least-squares-style: D = diag(w)
+D_diag = np.array([0.5, 1.0, 1.5, 2.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])
+m = Model(
+    objective_matrix=Q,
+    objective_matrix_low_rank=R,
+    objective_matrix_low_rank_middle=D_diag,
+    objective_vector=c,
+)
+m.optimize()
+
+# (b) Dense (possibly indefinite) D, e.g., from a quasi-Newton compact form
+M = np.random.randn(r, r)
+D_dense = 0.5 * (M + M.T)  # symmetric, no PSD requirement
+m = Model(
+    objective_matrix=Q,
+    objective_matrix_low_rank=R,
+    objective_matrix_low_rank_middle=D_dense,
+    objective_vector=c,
+)
+m.optimize()
+```
+
 ### Warm Starting
 
 ```python

docs/index.md

@@ -4,20 +4,21 @@ PDHCG-II is a high-performance, GPU-accelerated implementation of the Primal-Dua
 
 ## Problem Formulation
 
-PDHCG solves convex quadratic programs in the following form:
+PDHCG solves quadratic programs in the following form:
 
 $$
 \begin{aligned}
-\min_{x} \quad & \frac{1}{2}x^\top (Q + R^\top R) x + c^\top x \\
+\min_{x} \quad & \frac{1}{2}x^\top (Q + R^\top D R) x + c^\top x \\
 \text{s.t.} \quad & \ell_c \le Ax \le u_c, \\
                   & \ell_v \le x \le u_v.
 \end{aligned}
 $$
 
 Where:
 
-- $Q$ is a sparse positive semi-definite matrix (optional)
-- $R$ is a low-rank matrix such that $R^\top R$ represents a low-rank component (optional)
+- $Q$ is a sparse symmetric matrix (optional)
+- $R \in \mathbb{R}^{k\times n}$ is a low-rank factor of rank $k$ (optional)
+- $D \in \mathbb{R}^{k\times k}$ is an optional middle matrix; defaults to the identity, recovering the standard $Q + R^\top R$ form. May be diagonal, sparse, dense, or indefinite — the backend auto-detects the cheapest representation
 - $A$ is the constraint matrix
 - $c$ is the linear objective vector
 - $\ell_c, u_c$ are constraint bounds
@@ -26,7 +27,7 @@ Where:
 ## Key Features
 
 - **GPU Acceleration**: Fully leverages NVIDIA CUDA for extreme-scale QP problems
-- **Flexible Problem Structure**: Supports sparse quadratic terms, low-rank quadratic terms, or both
+- **Flexible Problem Structure**: Supports sparse, low-rank, and middle-weighted low-rank ($R^\top D R$) quadratic terms — alone or combined
 - **High Performance**: Competitive with commercial solvers on large-scale problems
 - **SpMVOp Auto-Detection**: Automatically uses cuSPARSE SpMVOp on CUDA 13+ while falling back to standard SpMV on CUDA 12.x
 - **Multi-GPU Distributed Solving**: Supports parallel solving across multiple GPUs via MPI and NCCL (optional, enabled at compile time)

docs/python/model.md

@@ -13,6 +13,7 @@
         - setObjectiveConstant
         - setObjectiveMatrix
         - setObjectiveMatrixLowRank
+        - setObjectiveMatrixLowRankMiddle
         - setConstraintMatrix
         - setConstraintLowerBound
         - setConstraintUpperBound