diff --git a/README.md b/README.md
index cad0353..a25cf15 100644
--- a/README.md
+++ b/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,
diff --git a/distributed/distributed_utils.cu b/distributed/distributed_utils.cu
index adc088b..ba86354 100644
--- a/distributed/distributed_utils.cu
+++ b/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);
diff --git a/docs/C_API.md b/docs/C_API.md
index c82a1f3..4199b48 100644
--- a/docs/C_API.md
+++ b/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
diff --git a/docs/c/functions.md b/docs/c/functions.md
index 3491cf4..00d3751 100644
--- a/docs/c/functions.md
+++ b/docs/c/functions.md
@@ -7,6 +7,7 @@ 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,
@@ -14,7 +15,8 @@ qp_problem_t *create_qp_problem(
 );
 ```
 
-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:**
 
@@ -22,7 +24,8 @@ Creates a QP problem from matrix descriptors.
 |-----------|-------------|
 | `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) |
diff --git a/docs/c/overview.md b/docs/c/overview.md
index 09766a9..6764266 100644
--- a/docs/c/overview.md
+++ b/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,6 +95,7 @@ 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,
@@ -102,7 +103,9 @@ qp_problem_t *create_qp_problem(
 );
 ```
 
-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
 
diff --git a/docs/c/types.md b/docs/c/types.md
index d83b9c0..b986048 100644
--- a/docs/c/types.md
+++ b/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;
 
diff --git a/docs/examples.md b/docs/examples.md
index 5262f6e..3bc1cb2 100644
--- a/docs/examples.md
+++ b/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
diff --git a/docs/index.md b/docs/index.md
index 48de299..759b918 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -4,11 +4,11 @@ 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}
@@ -16,8 +16,9 @@ $$
 
 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)
diff --git a/docs/python/model.md b/docs/python/model.md
index ef30c6c..a10f289 100644
--- a/docs/python/model.md
+++ b/docs/python/model.md
@@ -13,6 +13,7 @@
         - setObjectiveConstant
         - setObjectiveMatrix
         - setObjectiveMatrixLowRank
+        - setObjectiveMatrixLowRankMiddle
         - setConstraintMatrix
         - setConstraintLowerBound
         - setConstraintUpperBound
diff --git a/docs/python/quickstart.md b/docs/python/quickstart.md
index 014652b..2d9097b 100644
--- a/docs/python/quickstart.md
+++ b/docs/python/quickstart.md
@@ -9,8 +9,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, i.e. 0.5 * x'(Q + R^T R)x + c'x.)
 
 # 1. Define Standard QP terms
 Q = sp.csc_matrix([[1.0, -1.0], [-1.0, 2.0]])
@@ -20,16 +21,20 @@ c = np.array([-2.0, -6.0])
 # This adds 0.5 * ||Rx||^2 to the objective
 R = sp.csc_matrix([[1.0, 0.0]])
 
-# 3. Define Constraints
+# 3. (Optional) Middle matrix D for R^T D R; 1-D = diag, 2-D = dense, sparse OK.
+# D = np.array([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)
+# 5. Create QP model with Low-Rank term (R) and optional middle D
 m = Model(objective_matrix=Q,
           objective_matrix_low_rank=R,
+          # objective_matrix_low_rank_middle=D,
           objective_vector=c,
           constraint_matrix=A,
           constraint_lower_bound=l,
@@ -56,7 +61,7 @@ The `Model` class is the core interface for defining QP problems. The problem fo
 
 $$
 \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}
@@ -69,7 +74,8 @@ $$
 ### Optional Parameters
 
 - `objective_matrix` ($Q$): Sparse quadratic coefficients
-- `objective_matrix_low_rank` ($R$): Low-rank quadratic component (stores $R$, objective gets $R^\top R$)
+- `objective_matrix_low_rank` ($R$): Low-rank quadratic factor of shape $(k, n)$
+- `objective_matrix_low_rank_middle` ($D$, $k\times k$): Middle matrix in $R^\top D R$. 1-D array → diagonal $D$; 2-D array → dense symmetric $D$; scipy sparse → sparse $D$. May be indefinite. Defaults to identity
 - `constraint_matrix` ($A$): Linear constraint matrix
 - `constraint_lower_bound` ($\ell_c$): Constraint lower bounds
 - `constraint_upper_bound` ($u_c$): Constraint upper bounds
diff --git a/include/pdhcg.h b/include/pdhcg.h
index 7e14626..5b096d7 100644
--- a/include/pdhcg.h
+++ b/include/pdhcg.h
@@ -24,10 +24,12 @@ extern "C"
 {
 #endif
 
-    // create an qp_problem_t from a matrix descriptor
+    // create an qp_problem_t from matrix descriptors.
+    // pass NULL for the default D = I.
     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,
diff --git a/include/pdhcg_types.h b/include/pdhcg_types.h
index 38559de..1ac91ad 100644
--- a/include/pdhcg_types.h
+++ b/include/pdhcg_types.h
@@ -79,6 +79,9 @@ extern "C"
         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;
 
diff --git a/include/presolve_wrapper.h b/include/presolve_wrapper.h
index 112b730..4a34b82 100644
--- a/include/presolve_wrapper.h
+++ b/include/presolve_wrapper.h
@@ -80,7 +80,9 @@ extern "C"
                                            const int *Pp,
                                            size_t Pnnz);
 
-    /* Presolve QP in QR format: P = Q + RR^T */
+    /* Presolve QP in QR format: P = Q + R^T R.
+     * Note: PSQP does not support the optional middle matrix D from
+     * Q + R^T D R; the solver auto-disables presolve when D != I. */
     PDHCG_PresolvedData *pdhcg_presolve_qr(const double *Ax,
                                            const int *Ai,
                                            const int *Ap,
diff --git a/internal/internal_types.h b/internal/internal_types.h
index 28c1811..ccf97c8 100644
--- a/internal/internal_types.h
+++ b/internal/internal_types.h
@@ -60,6 +60,11 @@ typedef struct
     cusparseDnVecDescr_t vec_RRx_prod;
     int num_rank_lowrank_obj;
 
+    int lowrank_middle_type;
+    double *d_middle_diag;
+    double *d_middle_dense;
+    double *Rx_buffer;
+
     // Buffer for Distributed Version
     double *global_primal_obj_product;
     cusparseDnVecDescr_t vec_global_primal_obj_prod;
@@ -194,6 +199,13 @@ typedef struct
     grid_context_t *grid_context;
 } pdhg_solver_state_t;
 
+typedef enum
+{
+    PDHCG_D_NONE = 0,
+    PDHCG_D_DIAG = 1,
+    PDHCG_D_DENSE = 2
+} lowrank_middle_kind_t;
+
 typedef struct
 {
     int num_variables;
@@ -213,6 +225,10 @@ typedef struct
     CsrComponent *objective_lowrank_matrix;
     int objective_lowrank_matrix_num_nonzeros;
 
+    lowrank_middle_kind_t objective_lowrank_middle_kind;
+    double *objective_lowrank_middle_diag;
+    double *objective_lowrank_middle_dense;
+
     double *diagonal_quad_objective;
 
     double *constraint_lower_bound;
diff --git a/internal/pdhcg_kernels.cuh b/internal/pdhcg_kernels.cuh
index 787bea6..64592ff 100644
--- a/internal/pdhcg_kernels.cuh
+++ b/internal/pdhcg_kernels.cuh
@@ -32,6 +32,8 @@ extern "C"
     __global__ void
     element_wise_mul_kernel(const double *__restrict__ A, const double *__restrict__ B, double *__restrict__ C, int n);
 
+    __global__ void element_wise_mul_inplace_kernel(double *__restrict__ x, const double *__restrict__ d, int n);
+
     // ======================================================================
     // Advanced Metrics & Reduced Costs
     // ======================================================================
@@ -179,6 +181,17 @@ extern "C"
 
     __global__ void compute_csr_row_sq_norm_kernel(const int *row_ptr, const double *val, double *out, int num_rows);
 
+    __global__ void compute_csr_row_sq_norm_weighted_kernel(
+        const int *row_ptr, const int *col_ind, const double *val, const double *weights, double *out, int num_rows);
+
+    __global__ void compute_csr_row_quad_form_dense_kernel(const int *row_ptr,
+                                                           const int *col_ind,
+                                                           const double *val,
+                                                           const double *D_dense,
+                                                           int rank,
+                                                           double *out,
+                                                           int num_rows);
+
     __global__ void refresh_inner_precond_kernel(
         const double *diag_h_static, double inv_tau, double *m_diag, double *m_inv, int n_vars);
 
diff --git a/pdhcg/_pdhcg_core.pyi b/pdhcg/_pdhcg_core.pyi
index aab82df..fee1491 100644
--- a/pdhcg/_pdhcg_core.pyi
+++ b/pdhcg/_pdhcg_core.pyi
@@ -8,5 +8,5 @@ def get_default_params() -> dict:
     """
     Return default PDHG parameters as a dict
     """
-def solve_once(Q: typing.Any, R: typing.Any, A: typing.Any, objective_vector: typing.Any, objective_constant: typing.Any = None, variable_lower_bound: typing.Any = None, variable_upper_bound: typing.Any = None, constraint_lower_bound: typing.Any = None, constraint_upper_bound: typing.Any = None, zero_tolerance: typing.SupportsFloat | typing.SupportsIndex = 0.0, params: typing.Any = None, primal_start: typing.Any = None, dual_start: typing.Any = None) -> dict:
+def solve_once(Q: typing.Any, R: typing.Any, A: typing.Any, objective_vector: typing.Any, objective_constant: typing.Any = None, variable_lower_bound: typing.Any = None, variable_upper_bound: typing.Any = None, constraint_lower_bound: typing.Any = None, constraint_upper_bound: typing.Any = None, zero_tolerance: typing.SupportsFloat | typing.SupportsIndex = 0.0, params: typing.Any = None, primal_start: typing.Any = None, dual_start: typing.Any = None, D: typing.Any = None) -> dict:
     ...
diff --git a/python/README.md b/python/README.md
index cf94550..b7c2924 100644
--- a/python/README.md
+++ b/python/README.md
@@ -60,8 +60,9 @@ import numpy as np
 import scipy.sparse as sp
 from pdhcg import Model, PDHCG
 
-# Example: minimize 0.5 * x'Qx + 0.5 * ||Rx||^2 + 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 0.5 * x'Qx + 0.5 * ||Rx||^2.)
 
 # 1. Define Standard QP terms
 Q = sp.csc_matrix([[1.0, -1.0], [-1.0, 2.0]])
@@ -72,16 +73,21 @@ c = np.array([-2.0, -6.0])
 # This effectively adds 0.5 * (x1)^2 to the objective
 R = sp.csc_matrix([[1.0, 0.0]])
 
-# 3. Define Constraints
+# 3. (Optional) Middle matrix D. Pass a 1-D numpy array for a diagonal D,
+# a 2-D array for dense D, or a scipy.sparse matrix. Omit to use D = I.
+# D = np.array([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
+# 5. Create QP model with Low-Rank term
 m = Model(objective_matrix=Q,
-          objective_matrix_low_rank=R,  # <--- Pass R here
+          objective_matrix_low_rank=R,                # <--- Pass R here
+          # objective_matrix_low_rank_middle=D,      # <--- and optionally D
           objective_vector=c,
           constraint_matrix=A,
           constraint_lower_bound=l,
@@ -112,7 +118,7 @@ print("Dual solution:", m.Pi)
 The `Model` class represents a quadratic programming problem of the form:
 
 $$
-\min \frac{1}{2} x^\top (Q + R^\top R) x + c^\top x + c_0 \quad
+\min \frac{1}{2} x^\top (Q + R^\top D R) x + c^\top x + c_0 \quad
 \text{s.t.} \; \ell \le A x \le u, \quad
 \text{lb} \le x \le \text{ub}.
 $$
@@ -121,6 +127,7 @@ $$
 
 - **objective_matrix** (`Q`, optional): Quadratic part of the objective function. Both dense (`numpy.ndarray`) and sparse (`scipy.sparse.csr_matrix`) inputs are supported.
 - **objective_matrix_low_rank** (`R`, optional): Low-rank factor matrix in the quadratic objective term. Both dense (`numpy.ndarray`) and sparse (`scipy.sparse.csr_matrix`) inputs are supported.
+- **objective_matrix_low_rank_middle** (`D`, optional): Middle matrix in $R^\top D R$, $\text{rank}\times\text{rank}$. Accepts a 1-D numpy array (diagonal $D$), a 2-D numpy array (dense symmetric $D$), or a scipy sparse matrix. May be indefinite — useful for quasi-Newton updates, kernel-based formulations, and weighted least squares. Defaults to identity, recovering $R^\top R$.
 - **objective_vector** (`c`): Linear part of the objective function.
 - **constraint_matrix** (`A`): Coefficient matrix for the constraints. Both dense (`numpy.ndarray`) and sparse (`scipy.sparse.csr_matrix`) inputs are supported.
 - **constraint_lower_bound** (`l`): Lower bounds for each constraint. Use `-np.inf` or `None` for no lower bound.
diff --git a/python/pdhcg/model.py b/python/pdhcg/model.py
index 901c139..87da0e5 100644
--- a/python/pdhcg/model.py
+++ b/python/pdhcg/model.py
@@ -87,11 +87,14 @@ class Model:
     The quadratic programming problem is defined as:
 
     ```
-    minimize      1/2 x^T (Q + R^T R) x + c^T x
+    minimize      1/2 x^T (Q + R^T D R) x + c^T x
     subject to    l_c <= A x <= u_c
                   l_v <= x <= u_v
     ```
 
+    where D is an optional rank-by-rank middle matrix; if omitted, D defaults to
+    the identity, recovering the standard Q + R^T R formulation.
+
     Solver Parameters:
         Parameters can be set via the `Params` attribute or `setParam()` method.
         Common parameters include:
@@ -115,6 +118,7 @@ def __init__(
         constraint_upper_bound: Optional[ArrayLike] = None,
         objective_matrix: Optional[Union[np.ndarray, sp.spmatrix]] = None,
         objective_matrix_low_rank: Optional[Union[np.ndarray, sp.spmatrix]] = None,
+        objective_matrix_low_rank_middle: Optional[ArrayLike] = None,
         variable_lower_bound: Optional[ArrayLike] = None,
         variable_upper_bound: Optional[ArrayLike] = None,
         objective_constant: float = 0.0,
@@ -129,6 +133,9 @@ def __init__(
             constraint_upper_bound: Upper bounds for the linear constraints.
             objective_matrix: Quadratic coefficients of the objective function (Q).
             objective_matrix_low_rank: Low-rank quadratic coefficients of the objective (R).
+            objective_matrix_low_rank_middle: Optional middle matrix D in Q + R^T D R.
+                Accepts a 1-D array of length `rank` (treated as diag(D)) or a
+                2-D `rank` x `rank` symmetric array. Defaults to identity.
             variable_lower_bound: Lower bounds for the decision variables.
             variable_upper_bound: Upper bounds for the decision variables.
             objective_constant: Constant term in the objective function.
@@ -197,6 +204,7 @@ def __init__(
         self.setObjectiveConstant(objective_constant)
         self.setObjectiveMatrix(objective_matrix)
         self.setObjectiveMatrixLowRank(objective_matrix_low_rank)
+        self.setObjectiveMatrixLowRankMiddle(objective_matrix_low_rank_middle)
         self.setConstraintMatrix(constraint_matrix)
         self.setConstraintLowerBound(constraint_lower_bound)
         self.setConstraintUpperBound(constraint_upper_bound)
@@ -300,6 +308,56 @@ def setObjectiveMatrixLowRank(self, R_like: ArrayLike) -> None:
 
         self._clear_solution_cache()
 
+    def setObjectiveMatrixLowRankMiddle(
+        self, D_like: Optional[Union[np.ndarray, sp.spmatrix, ArrayLike]]
+    ) -> None:
+        """
+        Overwrite the middle matrix D in Q + R^T D R.
+
+        Accepts any of:
+          - None -> D = I (identity, original behavior).
+          - 1-D array of length `rank` -> diagonal D.
+          - 2-D `rank` x `rank` numpy array -> dense (symmetric) D.
+          - scipy.sparse `rank` x `rank` matrix -> sparse D.
+
+        The backend's `preprocess_qp_problem` inspects the nonzero pattern and
+        picks either a diagonal or a dense runtime representation.
+        """
+        if D_like is None:
+            self.D = None
+            self._clear_solution_cache()
+            return
+        if getattr(self, "R", None) is None:
+            raise ValueError(
+                "setObjectiveMatrixLowRankMiddle: D is only meaningful when R is set; "
+                "call setObjectiveMatrixLowRank first."
+            )
+        rank = int(self.R.shape[0])
+        if sp.issparse(D_like):
+            if D_like.shape != (rank, rank):
+                raise ValueError(
+                    f"setObjectiveMatrixLowRankMiddle: sparse D shape {D_like.shape} must be ({rank}, {rank})"
+                )
+            self.D = _as_csr_f64_i32(D_like)
+        else:
+            d_arr = _as_dense_f64_c(D_like)
+            if d_arr.ndim == 1:
+                if d_arr.size != rank:
+                    raise ValueError(
+                        f"setObjectiveMatrixLowRankMiddle: diag D length {d_arr.size} must equal rank {rank}"
+                    )
+            elif d_arr.ndim == 2:
+                if d_arr.shape != (rank, rank):
+                    raise ValueError(
+                        f"setObjectiveMatrixLowRankMiddle: dense D shape {d_arr.shape} must be ({rank}, {rank})"
+                    )
+            else:
+                raise ValueError(
+                    f"setObjectiveMatrixLowRankMiddle: D must be 1D, 2D, or scipy.sparse; got ndim={d_arr.ndim}"
+                )
+            self.D = d_arr
+        self._clear_solution_cache()
+
     def setConstraintMatrix(self, A_like: ArrayLike) -> None:
         """
         Overwrite the linear constraint matrix A.
@@ -518,6 +576,7 @@ def optimize(self):
             params=self._params,
             primal_start=self._primal_start,
             dual_start=self._dual_start,
+            D=getattr(self, "D", None),
         )
         # solutions
         self._x = np.asarray(info.get("X")) if info.get("X") is not None else None
diff --git a/python_bindings/_core_bindings.cpp b/python_bindings/_core_bindings.cpp
index 1919a42..db81950 100644
--- a/python_bindings/_core_bindings.cpp
+++ b/python_bindings/_core_bindings.cpp
@@ -493,7 +493,8 @@ static py::dict solve_once(py::object Q,
                            double zero_tolerance = 0.0,
                            py::object params = py::none(),
                            py::object primal_start = py::none(),
-                           py::object dual_start = py::none())
+                           py::object dual_start = py::none(),
+                           py::object D = py::none())
 {
     static std::once_flag cuda_init_flag;
     std::call_once(cuda_init_flag, []() { cudaFree(0); });
@@ -554,8 +555,75 @@ static py::dict solve_once(py::object Q,
     const matrix_desc_t *r_desc_ptr = R.is_none() ? nullptr : &view_r.desc;
     const matrix_desc_t *a_desc_ptr = A.is_none() ? nullptr : &view_a.desc;
 
+    PyMatrixView view_d;
+    std::vector<int32_t> d_diag_rp, d_diag_ci;
+    std::vector<double> d_diag_vv;
+    const matrix_desc_t *d_desc_ptr = nullptr;
+    if (D && !D.is_none())
+    {
+        if (R.is_none())
+        {
+            throw std::invalid_argument("D was provided but R is None; D is only meaningful with a low-rank R.");
+        }
+        int rank = view_r.desc.m;
+        bool is_1d = false;
+        if (py::isinstance<py::array>(D))
+        {
+            py::array d_arr = py::cast<py::array>(D);
+            if (d_arr.ndim() == 1)
+                is_1d = true;
+        }
+        if (is_1d)
+        {
+            /* Build a CSR diag(d) directly. */
+            py::array_t<double, py::array::c_style | py::array::forcecast> d64(py::cast<py::array>(D));
+            if ((int)d64.size() != rank)
+            {
+                throw std::invalid_argument("D (diag) length " + std::to_string((int)d64.size()) + " must equal rank " +
+                                            std::to_string(rank));
+            }
+            const double *p = d64.data();
+            d_diag_rp.reserve(rank + 1);
+            d_diag_ci.reserve(rank);
+            d_diag_vv.reserve(rank);
+            int nz = 0;
+            d_diag_rp.push_back(0);
+            for (int i = 0; i < rank; ++i)
+            {
+                if (p[i] != 0.0)
+                {
+                    d_diag_ci.push_back(i);
+                    d_diag_vv.push_back(p[i]);
+                    ++nz;
+                }
+                d_diag_rp.push_back(nz);
+            }
+            view_d.desc.m = rank;
+            view_d.desc.n = rank;
+            view_d.desc.fmt = matrix_csr;
+            view_d.desc.zero_tolerance = 0.0;
+            view_d.desc.data.csr.nnz = (int)d_diag_vv.size();
+            view_d.desc.data.csr.row_ptr = d_diag_rp.data();
+            view_d.desc.data.csr.col_ind = d_diag_ci.data();
+            view_d.desc.data.csr.vals = d_diag_vv.data();
+            d_desc_ptr = &view_d.desc;
+        }
+        else
+        {
+            view_d = get_matrix_from_python(D, 0.0);
+            if (view_d.desc.m != rank || view_d.desc.n != rank)
+            {
+                throw std::invalid_argument("D shape (" + std::to_string(view_d.desc.m) + ", " +
+                                            std::to_string(view_d.desc.n) + ") must be (" + std::to_string(rank) +
+                                            ", " + std::to_string(rank) + ")");
+            }
+            view_a.keep.owners.insert(view_a.keep.owners.end(), view_d.keep.owners.begin(), view_d.keep.owners.end());
+            d_desc_ptr = &view_d.desc;
+        }
+    }
+
     qp_problem_t *prob =
-        create_qp_problem(c_ptr, q_desc_ptr, r_desc_ptr, a_desc_ptr, l_ptr, u_ptr, lb_ptr, ub_ptr, c0_ptr);
+        create_qp_problem(c_ptr, q_desc_ptr, r_desc_ptr, d_desc_ptr, a_desc_ptr, l_ptr, u_ptr, lb_ptr, ub_ptr, c0_ptr);
     if (!prob)
     {
         throw std::runtime_error("create_qp_problem failed.");
@@ -668,5 +736,6 @@ PYBIND11_MODULE(_pdhcg_core, m)
           py::arg("zero_tolerance") = 0.0,
           py::arg("params") = py::none(),
           py::arg("primal_start") = py::none(),
-          py::arg("dual_start") = py::none());
+          py::arg("dual_start") = py::none(),
+          py::arg("D") = py::none());
 }
diff --git a/src/pdhcg.c b/src/pdhcg.c
index 674f9b6..2299e69 100644
--- a/src/pdhcg.c
+++ b/src/pdhcg.c
@@ -34,6 +34,7 @@ volatile sig_atomic_t g_pdhcg_cancel_request = 0;
 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,
@@ -313,6 +314,68 @@ qp_problem_t *create_qp_problem(const double *objective_c,
         prob->objective_lowrank_matrix_num_nonzeros = 0;
     }
 
+    prob->objective_lowrank_middle_matrix = NULL;
+    prob->objective_lowrank_middle_matrix_num_nonzeros = 0;
+    if (D_desc)
+    {
+        int k = prob->num_rank_lowrank_obj;
+        if (k <= 0)
+        {
+            fprintf(stderr, "[interface] D matrix ignored: problem has no low-rank component.\n");
+        }
+        else if (D_desc->m != k || D_desc->n != k)
+        {
+            fprintf(stderr, "[interface] D matrix shape (%d, %d) must be (%d, %d).\n", D_desc->m, D_desc->n, k, k);
+            qp_problem_free(prob);
+            return NULL;
+        }
+        else
+        {
+            prob->objective_lowrank_middle_matrix = (CsrComponent *)safe_calloc(1, sizeof(CsrComponent));
+            int *rp = NULL, *ci = NULL;
+            double *vv = NULL;
+            int nnz = 0;
+            int rc = 0;
+            switch (D_desc->fmt)
+            {
+                case matrix_dense:
+                    rc = dense_to_csr(D_desc, &rp, &ci, &vv, &nnz);
+                    break;
+                case matrix_csc:
+                    rc = csc_to_csr(D_desc, &rp, &ci, &vv, &nnz);
+                    break;
+                case matrix_coo:
+                    rc = coo_to_csr(D_desc, &rp, &ci, &vv, &nnz);
+                    break;
+                case matrix_csr:
+                    nnz = D_desc->data.csr.nnz;
+                    rp = (int *)safe_malloc((size_t)(k + 1) * sizeof(int));
+                    ci = (int *)safe_malloc((size_t)nnz * sizeof(int));
+                    vv = (double *)safe_malloc((size_t)nnz * sizeof(double));
+                    memcpy(rp, D_desc->data.csr.row_ptr, (size_t)(k + 1) * sizeof(int));
+                    memcpy(ci, D_desc->data.csr.col_ind, (size_t)nnz * sizeof(int));
+                    memcpy(vv, D_desc->data.csr.vals, (size_t)nnz * sizeof(double));
+                    break;
+                default:
+                    rc = -1;
+                    fprintf(stderr, "[interface] D matrix: unsupported format %d.\n", D_desc->fmt);
+                    break;
+            }
+            if (rc != 0)
+            {
+                free(rp);
+                free(ci);
+                free(vv);
+                qp_problem_free(prob);
+                return NULL;
+            }
+            prob->objective_lowrank_middle_matrix->row_ptr = rp;
+            prob->objective_lowrank_middle_matrix->col_ind = ci;
+            prob->objective_lowrank_middle_matrix->val = vv;
+            prob->objective_lowrank_middle_matrix_num_nonzeros = nnz;
+        }
+    }
+
     // default fill values
     prob->objective_constant = objective_constant ? *objective_constant : 0.0;
     fill_or_copy(&prob->objective_vector, prob->num_variables, objective_c, 0.0);
@@ -358,6 +421,8 @@ void qp_problem_free(qp_problem_t *prob)
     free(prob->constraint_upper_bound);
     free(prob->primal_start);
     free(prob->dual_start);
+    csr_component_free(prob->objective_lowrank_middle_matrix);
+    free(prob->objective_lowrank_middle_matrix);
     memset(prob, 0, sizeof(*prob));
     free(prob);
 }
diff --git a/src/pdhcg_kernels.cu b/src/pdhcg_kernels.cu
index ad35f14..a555ed5 100644
--- a/src/pdhcg_kernels.cu
+++ b/src/pdhcg_kernels.cu
@@ -322,6 +322,61 @@ __global__ void compute_csr_row_sq_norm_kernel(const int *row_ptr, const double
     }
 }
 
+__global__ void element_wise_mul_inplace_kernel(double *__restrict__ x, const double *__restrict__ d, int n)
+{
+    int i = blockIdx.x * blockDim.x + threadIdx.x;
+    if (i < n)
+    {
+        x[i] *= d[i];
+    }
+}
+
+__global__ void compute_csr_row_sq_norm_weighted_kernel(
+    const int *row_ptr, const int *col_ind, const double *val, const double *weights, double *out, int num_rows)
+{
+    int i = blockIdx.x * blockDim.x + threadIdx.x;
+    if (i < num_rows)
+    {
+        double sum = 0.0;
+        int start = row_ptr[i];
+        int end = row_ptr[i + 1];
+        for (int k = start; k < end; ++k)
+        {
+            double v = val[k];
+            sum += weights[col_ind[k]] * v * v;
+        }
+        out[i] = sum;
+    }
+}
+
+__global__ void compute_csr_row_quad_form_dense_kernel(const int *row_ptr,
+                                                       const int *col_ind,
+                                                       const double *val,
+                                                       const double *D_dense,
+                                                       int rank,
+                                                       double *out,
+                                                       int num_rows)
+{
+    int i = blockIdx.x * blockDim.x + threadIdx.x;
+    if (i < num_rows)
+    {
+        int start = row_ptr[i];
+        int end = row_ptr[i + 1];
+        double sum = 0.0;
+        for (int ka = start; ka < end; ++ka)
+        {
+            int ra = col_ind[ka];
+            double va = val[ka];
+            const double *Drow = D_dense + (size_t)ra * (size_t)rank;
+            for (int kb = start; kb < end; ++kb)
+            {
+                sum += va * val[kb] * Drow[col_ind[kb]];
+            }
+        }
+        out[i] = sum;
+    }
+}
+
 __global__ void
 refresh_inner_precond_kernel(const double *diag_h_static, double inv_tau, double *m_diag, double *m_inv, int n_vars)
 {
diff --git a/src/pdhg_core_op.cu b/src/pdhg_core_op.cu
index 3145b11..abb532f 100644
--- a/src/pdhg_core_op.cu
+++ b/src/pdhg_core_op.cu
@@ -37,6 +37,39 @@ limitations under the License.
 #include "distributed_types.h"
 #endif
 
+static void apply_lowrank_middle(pdhg_solver_state_t *state)
+{
+    quadratic_objective_term_t *qot = state->quadratic_objective_term;
+    int rank = qot->num_rank_lowrank_obj;
+    if (qot->lowrank_middle_type == 0 || rank <= 0)
+        return;
+
+    if (qot->lowrank_middle_type == 1)
+    {
+        int nb = (rank + THREADS_PER_BLOCK - 1) / THREADS_PER_BLOCK;
+        element_wise_mul_inplace_kernel<<<nb, THREADS_PER_BLOCK>>>(qot->Rx_product, qot->d_middle_diag, rank);
+        return;
+    }
+
+    cublasPointerMode_t prev_mode;
+    CUBLAS_CHECK(cublasGetPointerMode(state->blas_handle, &prev_mode));
+    CUBLAS_CHECK(cublasSetPointerMode(state->blas_handle, CUBLAS_POINTER_MODE_HOST));
+    CUBLAS_CHECK(cublasDsymv(state->blas_handle,
+                             CUBLAS_FILL_MODE_LOWER,
+                             rank,
+                             &HOST_ONE,
+                             qot->d_middle_dense,
+                             rank,
+                             qot->Rx_product,
+                             1,
+                             &HOST_ZERO,
+                             qot->Rx_buffer,
+                             1));
+    CUBLAS_CHECK(cublasSetPointerMode(state->blas_handle, prev_mode));
+    CUDA_CHECK(
+        cudaMemcpyAsync(qot->Rx_product, qot->Rx_buffer, (size_t)rank * sizeof(double), cudaMemcpyDeviceToDevice));
+}
+
 void update_obj_product(pdhg_solver_state_t *state, double *primal_solution)
 {
     switch (state->quadratic_objective_term->quad_obj_type)
@@ -83,6 +116,8 @@ void update_obj_product(pdhg_solver_state_t *state, double *primal_solution)
                                    PDHCG_SCOPE_ROW,
                                    0);
 
+            apply_lowrank_middle(state);
+
             pdhcg_spmv_execute(state->sparse_handle,
                                state->quadratic_objective_term->spmv_ctx_Rt,
                                &HOST_ONE,
@@ -120,6 +155,8 @@ void update_obj_product(pdhg_solver_state_t *state, double *primal_solution)
                                    PDHCG_SCOPE_ROW,
                                    0);
 
+            apply_lowrank_middle(state);
+
             pdhcg_spmv_execute(state->sparse_handle,
                                state->quadratic_objective_term->spmv_ctx_Rt,
                                &HOST_ONE,
diff --git a/src/preconditioner.c b/src/preconditioner.c
index 2266115..035be99 100644
--- a/src/preconditioner.c
+++ b/src/preconditioner.c
@@ -62,6 +62,12 @@ qp_problem_t *deepcopy_problem(const qp_problem_t *prob)
     new_prob->objective_lowrank_matrix = deepcopy_csr_component(
         prob->objective_lowrank_matrix, prob->num_rank_lowrank_obj, prob->objective_lowrank_matrix_num_nonzeros);
 
+    new_prob->objective_lowrank_middle_matrix =
+        deepcopy_csr_component(prob->objective_lowrank_middle_matrix,
+                               prob->num_rank_lowrank_obj,
+                               prob->objective_lowrank_middle_matrix_num_nonzeros);
+    new_prob->objective_lowrank_middle_matrix_num_nonzeros = prob->objective_lowrank_middle_matrix_num_nonzeros;
+
     if (prob->primal_start)
     {
         new_prob->primal_start = safe_malloc(var_bytes);
@@ -394,6 +400,57 @@ processed_qp_problem_t *preprocess_qp_problem(const qp_problem_t *raw_problem)
     processed->objective_sparse_matrix = raw_problem->objective_sparse_matrix;
     processed->objective_lowrank_matrix = raw_problem->objective_lowrank_matrix;
 
+    processed->objective_lowrank_middle_kind = PDHCG_D_NONE;
+    processed->objective_lowrank_middle_diag = NULL;
+    processed->objective_lowrank_middle_dense = NULL;
+    {
+        CsrComponent *D_csr = raw_problem->objective_lowrank_middle_matrix;
+        int D_nnz = raw_problem->objective_lowrank_middle_matrix_num_nonzeros;
+        int k = raw_problem->num_rank_lowrank_obj;
+        if (D_csr && D_csr->row_ptr && D_nnz > 0 && k > 0)
+        {
+            bool diag_only = true;
+            for (int r = 0; r < k && diag_only; ++r)
+            {
+                for (int j = D_csr->row_ptr[r]; j < D_csr->row_ptr[r + 1]; ++j)
+                {
+                    if (D_csr->col_ind[j] != r)
+                    {
+                        diag_only = false;
+                        break;
+                    }
+                }
+            }
+
+            if (diag_only)
+            {
+                processed->objective_lowrank_middle_kind = PDHCG_D_DIAG;
+                processed->objective_lowrank_middle_diag = (double *)safe_calloc(k, sizeof(double));
+                for (int r = 0; r < k; ++r)
+                {
+                    for (int j = D_csr->row_ptr[r]; j < D_csr->row_ptr[r + 1]; ++j)
+                    {
+                        processed->objective_lowrank_middle_diag[r] = D_csr->val[j];
+                    }
+                }
+            }
+            else
+            {
+                processed->objective_lowrank_middle_kind = PDHCG_D_DENSE;
+                size_t rank2 = (size_t)k * (size_t)k;
+                processed->objective_lowrank_middle_dense = (double *)safe_calloc(rank2, sizeof(double));
+                for (int r = 0; r < k; ++r)
+                {
+                    for (int j = D_csr->row_ptr[r]; j < D_csr->row_ptr[r + 1]; ++j)
+                    {
+                        int c = D_csr->col_ind[j];
+                        processed->objective_lowrank_middle_dense[(size_t)r * (size_t)k + (size_t)c] = D_csr->val[j];
+                    }
+                }
+            }
+        }
+    }
+
     if ((!raw_problem->objective_sparse_matrix || raw_problem->objective_sparse_matrix_num_nonzeros == 0) &&
         (!raw_problem->objective_lowrank_matrix || raw_problem->objective_lowrank_matrix_num_nonzeros == 0))
     {
@@ -444,5 +501,10 @@ void free_processed_qp_problem(processed_qp_problem_t *processed)
         processed->diagonal_quad_objective = NULL;
     }
 
+    free(processed->objective_lowrank_middle_diag);
+    processed->objective_lowrank_middle_diag = NULL;
+    free(processed->objective_lowrank_middle_dense);
+    processed->objective_lowrank_middle_dense = NULL;
+
     free(processed);
 }
diff --git a/src/solver_state.cu b/src/solver_state.cu
index 26e54d5..2829323 100644
--- a/src/solver_state.cu
+++ b/src/solver_state.cu
@@ -192,6 +192,25 @@ static void initialize_lowrank_component_obj(pdhg_solver_state_t *state, const p
                               state->quadratic_objective_term->objective_lowrank_matrix_t->val,
                               state->quadratic_objective_term->vec_Rx_prod,
                               state->quadratic_objective_term->vec_primal_obj_prod);
+
+    state->quadratic_objective_term->lowrank_middle_type = (int)problem->objective_lowrank_middle_kind;
+    state->quadratic_objective_term->d_middle_diag = NULL;
+    state->quadratic_objective_term->d_middle_dense = NULL;
+    state->quadratic_objective_term->Rx_buffer = NULL;
+
+    int rank = problem->num_rank_lowrank_obj;
+    if (problem->objective_lowrank_middle_kind == PDHCG_D_DIAG && rank > 0)
+    {
+        ALLOC_AND_COPY(state->quadratic_objective_term->d_middle_diag,
+                       problem->objective_lowrank_middle_diag,
+                       (size_t)rank * sizeof(double));
+    }
+    else if (problem->objective_lowrank_middle_kind == PDHCG_D_DENSE && rank > 0)
+    {
+        size_t bytes = (size_t)rank * (size_t)rank * sizeof(double);
+        ALLOC_AND_COPY(state->quadratic_objective_term->d_middle_dense, problem->objective_lowrank_middle_dense, bytes);
+        ALLOC_ZERO(state->quadratic_objective_term->Rx_buffer, (size_t)rank * sizeof(double));
+    }
 }
 
 static void initialize_quadratic_obj_term(pdhg_solver_state_t *state, const processed_qp_problem_t *problem)
@@ -352,8 +371,29 @@ static void initialize_inner_solver(pdhg_solver_state_t *state, const pdhg_param
                 {
                     cu_sparse_matrix_csr_t *Rt = state->quadratic_objective_term->objective_lowrank_matrix_t;
                     double *out = state->inner_solver->bb_step_size->Ms_buffer;
-                    compute_csr_row_sq_norm_kernel<<<state->num_blocks_primal, THREADS_PER_BLOCK>>>(
-                        Rt->row_ptr, Rt->val, out, n);
+                    int mtype = state->quadratic_objective_term->lowrank_middle_type;
+                    if (mtype == 1)
+                    {
+                        compute_csr_row_sq_norm_weighted_kernel<<<state->num_blocks_primal, THREADS_PER_BLOCK>>>(
+                            Rt->row_ptr, Rt->col_ind, Rt->val, state->quadratic_objective_term->d_middle_diag, out, n);
+                    }
+                    else if (mtype == 2)
+                    {
+                        int rank = state->quadratic_objective_term->num_rank_lowrank_obj;
+                        compute_csr_row_quad_form_dense_kernel<<<state->num_blocks_primal, THREADS_PER_BLOCK>>>(
+                            Rt->row_ptr,
+                            Rt->col_ind,
+                            Rt->val,
+                            state->quadratic_objective_term->d_middle_dense,
+                            rank,
+                            out,
+                            n);
+                    }
+                    else
+                    {
+                        compute_csr_row_sq_norm_kernel<<<state->num_blocks_primal, THREADS_PER_BLOCK>>>(
+                            Rt->row_ptr, Rt->val, out, n);
+                    }
                     CUDA_CHECK(cudaGetLastError());
                     const double one = 1.0;
                     CUBLAS_CHECK(cublasDaxpy(
@@ -863,6 +903,13 @@ void pdhg_solver_state_free(pdhg_solver_state_t *state)
         if (state->quadratic_objective_term->spmv_ctx_Rt)
             pdhcg_spmv_ctx_destroy(state->quadratic_objective_term->spmv_ctx_Rt);
 
+        if (state->quadratic_objective_term->d_middle_diag)
+            CUDA_CHECK(cudaFree(state->quadratic_objective_term->d_middle_diag));
+        if (state->quadratic_objective_term->d_middle_dense)
+            CUDA_CHECK(cudaFree(state->quadratic_objective_term->d_middle_dense));
+        if (state->quadratic_objective_term->Rx_buffer)
+            CUDA_CHECK(cudaFree(state->quadratic_objective_term->Rx_buffer));
+
         free(state->quadratic_objective_term);
     }
 
diff --git a/test/test_lowrank_middle_D.py b/test/test_lowrank_middle_D.py
new file mode 100644
index 0000000..3574238
--- /dev/null
+++ b/test/test_lowrank_middle_D.py
@@ -0,0 +1,188 @@
+"""
+Smoke test for the optional middle matrix D in Q + R^T D R.
+
+For each variant (D = identity, diag, dense PSD), we solve the same QP two ways:
+  (a) Using the new D parameter on a sparse R.
+  (b) Folding D into R as R' = sqrt(D) R (only valid when D is PSD), so the
+      "no D" path solves the same effective QP.
+The two solvers should return primal solutions that agree to a few digits.
+
+We also check the indefinite-D case (where folding into R via real sqrt does
+not exist) by comparing the explicit Q-only formulation: build Q_eff = R^T D R
+and feed it as the sparse Q, then compare to the (R, D) formulation.
+"""
+import numpy as np
+import scipy.sparse as sp
+
+from pdhcg import Model, PDHCG
+
+
+def _make_problem(n=40, rank=6, n_cons=20, seed=0):
+    rng = np.random.default_rng(seed)
+    R_dense = rng.standard_normal((rank, n)) * 0.5
+    # Make R modestly sparse but valid CSR
+    R_dense[np.abs(R_dense) < 0.25] = 0.0
+    R = sp.csr_matrix(R_dense)
+
+    A = sp.csr_matrix(rng.standard_normal((n_cons, n)))
+    c = rng.standard_normal(n)
+    b = rng.standard_normal(n_cons)
+    return n, rank, R, A, c, b
+
+
+def _solve(R, A, c, b, *, D=None, Q=None, iter_limit=5000):
+    m = Model(
+        objective_vector=c,
+        constraint_matrix=A,
+        constraint_lower_bound=b,
+        constraint_upper_bound=b,
+        objective_matrix=Q,
+        objective_matrix_low_rank=R,
+        objective_matrix_low_rank_middle=D,
+        variable_lower_bound=-np.full(c.size, 10.0),
+        variable_upper_bound=np.full(c.size, 10.0),
+    )
+    m.ModelSense = PDHCG.MINIMIZE
+    m.setParams(
+        LogLevel=0,
+        OptimalityTol=1e-4,
+        FeasibilityTol=1e-4,
+        Presolve=False,
+        IterationLimit=iter_limit,
+        TimeLimit=30,
+    )
+    m.optimize()
+    if m.Status not in ("OPTIMAL", "ITERATION_LIMIT"):
+        raise AssertionError(f"unexpected status {m.Status}")
+    return m.X, m.ObjVal
+
+
+def test_identity_matches_no_D():
+    n, rank, R, A, c, b = _make_problem()
+    D = np.ones(rank)  # diag identity
+    x_ref, obj_ref = _solve(R, A, c, b)
+    x_d, obj_d = _solve(R, A, c, b, D=D)
+    assert np.allclose(x_ref, x_d, atol=1e-2, rtol=1e-2), (
+        f"identity D mismatch, max diff {np.max(np.abs(x_ref - x_d))}"
+    )
+    print(f"[identity]    obj_ref={obj_ref:.6e}  obj_D={obj_d:.6e}")
+
+
+def test_diag_D_matches_folded_R():
+    n, rank, R, A, c, b = _make_problem(seed=1)
+    d = np.array([0.3, 1.7, 0.9, 2.4, 0.5, 1.1])
+    assert d.size == rank
+
+    # Fold into R: R' = sqrt(D) R  ->  R'^T R' = R^T D R, valid since d > 0
+    R_folded = sp.diags(np.sqrt(d)) @ R
+
+    x_ref, obj_ref = _solve(R_folded, A, c, b)
+    x_d, obj_d = _solve(R, A, c, b, D=d)
+    assert np.allclose(x_ref, x_d, atol=1e-2, rtol=1e-2), (
+        f"diag D mismatch, max diff {np.max(np.abs(x_ref - x_d))}"
+    )
+    print(f"[diag PSD]    obj_ref={obj_ref:.6e}  obj_D={obj_d:.6e}")
+
+
+def test_dense_D_matches_folded_R():
+    n, rank, R, A, c, b = _make_problem(seed=2)
+    # Build a random symmetric PSD D = M^T M + alpha I
+    rng = np.random.default_rng(3)
+    M = rng.standard_normal((rank, rank))
+    D = M.T @ M + 0.2 * np.eye(rank)
+    D = 0.5 * (D + D.T)  # enforce symmetry
+
+    # Fold into R via Cholesky: R' = L^T R where D = L L^T
+    L = np.linalg.cholesky(D)
+    R_folded = sp.csr_matrix(L.T @ R.toarray())
+
+    x_ref, obj_ref = _solve(R_folded, A, c, b)
+    x_d, obj_d = _solve(R, A, c, b, D=D)
+    assert np.allclose(x_ref, x_d, atol=1e-2, rtol=1e-2), (
+        f"dense PSD D mismatch, max diff {np.max(np.abs(x_ref - x_d))}"
+    )
+    print(f"[dense PSD]   obj_ref={obj_ref:.6e}  obj_D={obj_d:.6e}")
+
+
+def test_indefinite_dense_D_matches_explicit_Q():
+    """D with negative eigenvalues cannot be folded into R via real sqrt;
+    compare against an explicit Q built from Q = R^T D R + tiny shift."""
+    n, rank, R, A, c, b = _make_problem(seed=3)
+    rng = np.random.default_rng(4)
+    M = rng.standard_normal((rank, rank))
+    D = M.T @ M
+    # Inject a negative eigenvalue
+    eigvals, eigvecs = np.linalg.eigh(D)
+    eigvals[0] = -0.3
+    D = eigvecs @ np.diag(eigvals) @ eigvecs.T
+    D = 0.5 * (D + D.T)
+    # Verify indefinite
+    assert np.min(np.linalg.eigvalsh(D)) < 0
+
+    # Build effective Q = R^T D R + small PSD shift so the resulting Q is convex enough
+    R_dense = R.toarray()
+    Q_eff_dense = R_dense.T @ D @ R_dense
+    # Symmetrize numerically
+    Q_eff_dense = 0.5 * (Q_eff_dense + Q_eff_dense.T)
+
+    # Solve via explicit Q (no R, no D); nonconvex Q so run more iterations.
+    Q_eff_sparse = sp.csr_matrix(Q_eff_dense)
+    x_ref, obj_ref = _solve(None, A, c, b, Q=Q_eff_sparse, iter_limit=20000)
+
+    # Solve via (R, D=indefinite)
+    x_d, obj_d = _solve(R, A, c, b, D=D, iter_limit=20000)
+
+    # Objective is the most stable thing to compare on nonconvex problems.
+    assert abs(obj_ref - obj_d) / (abs(obj_ref) + 1.0) < 1e-3, (
+        f"indefinite D objective mismatch: ref={obj_ref:.6e} D={obj_d:.6e}"
+    )
+    print(f"[indef dense] obj_ref={obj_ref:.6e}  obj_D={obj_d:.6e}")
+
+
+def test_sparse_csr_D_matches_dense_D():
+    """Pass D as a scipy CSR matrix; should match feeding the equivalent dense D."""
+    n, rank, R, A, c, b = _make_problem(seed=4)
+    rng = np.random.default_rng(5)
+    M = rng.standard_normal((rank, rank))
+    D_dense = M.T @ M + 0.2 * np.eye(rank)
+    D_dense = 0.5 * (D_dense + D_dense.T)
+    # Sparsify the dense D so the CSR input genuinely has fewer than rank^2 nnz
+    mask = np.abs(D_dense) > 0.3
+    mask |= mask.T
+    np.fill_diagonal(mask, True)
+    D_dense_sparsified = np.where(mask, D_dense, 0.0)
+
+    D_sparse = sp.csr_matrix(D_dense_sparsified)
+
+    x_dense, obj_dense = _solve(R, A, c, b, D=D_dense_sparsified)
+    x_sparse, obj_sparse = _solve(R, A, c, b, D=D_sparse)
+    assert np.allclose(x_dense, x_sparse, atol=1e-2, rtol=1e-2), (
+        f"sparse vs dense D mismatch, max diff {np.max(np.abs(x_dense - x_sparse))}"
+    )
+    print(f"[sparse CSR D] obj_dense={obj_dense:.6e}  obj_sparse={obj_sparse:.6e}")
+
+
+def test_diagonal_sparse_D_detected_as_diag():
+    """A sparse D with only diagonal entries should auto-detect to the DIAG kind
+    and give the same answer as a 1-D diag array."""
+    n, rank, R, A, c, b = _make_problem(seed=6)
+    d = np.array([0.3, 1.7, 0.9, 2.4, 0.5, 1.1])
+    assert d.size == rank
+    D_sparse_diag = sp.diags(d).tocsr()
+
+    x_1d, obj_1d = _solve(R, A, c, b, D=d)
+    x_csr, obj_csr = _solve(R, A, c, b, D=D_sparse_diag)
+    assert np.allclose(x_1d, x_csr, atol=1e-2, rtol=1e-2), (
+        f"1-D vs sparse-diag mismatch, max diff {np.max(np.abs(x_1d - x_csr))}"
+    )
+    print(f"[diag via CSR] obj_1d={obj_1d:.6e}  obj_csr={obj_csr:.6e}")
+
+
+if __name__ == "__main__":
+    test_identity_matches_no_D()
+    test_diag_D_matches_folded_R()
+    test_dense_D_matches_folded_R()
+    test_indefinite_dense_D_matches_explicit_Q()
+    test_sparse_csr_D_matches_dense_D()
+    test_diagonal_sparse_D_detected_as_diag()
+    print("\nAll D-middle tests passed.")