diff --git a/.github/workflows/unit_tests.yml b/.github/workflows/unit_tests.yml
index 0da55d8..2963a91 100644
--- a/.github/workflows/unit_tests.yml
+++ b/.github/workflows/unit_tests.yml
@@ -50,8 +50,8 @@ jobs:
             c_compiler: clang
             cpp_compiler: clang++
           - os: windows-latest
-            c_compiler: MSVC
-            cpp_compiler: MSVC
+            c_compiler: cl
+            cpp_compiler: cl
           - precision: long_double
             long_double_flag: -DQOCO_LONG_DOUBLE_PRECISION:BOOL=True
           - precision: double
diff --git a/algebra/builtin/qdldl_backend.c b/algebra/builtin/qdldl_backend.c
index 07a5737..f04b216 100644
--- a/algebra/builtin/qdldl_backend.c
+++ b/algebra/builtin/qdldl_backend.c
@@ -224,7 +224,9 @@ static QOCOFloat compute_linsys_residual(LinSysData* linsys_data,
                                          QOCOWorkspace* work, QOCOFloat* b,
                                          QOCOFloat* x_scratch)
 {
-  QOCOFloat* Wfull = get_data_vectorf(work->Wfull);
+  QOCOFloat* nt_scaling = get_data_vectorf(work->nt_scaling);
+  QOCOInt* nt_scaling_soc_idx = get_data_vectori(work->nt_scaling_soc_idx);
+  QOCOInt* soc_idx = get_data_vectori(work->soc_idx);
   QOCOFloat* xbuff = get_data_vectorf(work->xbuff);
   QOCOFloat* ubuff1 = get_data_vectorf(work->ubuff1);
   QOCOFloat* ubuff2 = get_data_vectorf(work->ubuff2);
@@ -239,8 +241,8 @@ static QOCOFloat compute_linsys_residual(LinSysData* linsys_data,
   // Compute K_true * x_scratch -> xyzbuff2 against the unregularized matrix.
   // data->P stores the regularized P (P + eps_P * I), so subtract the P
   // regularization contribution from the x block to recover the true product.
-  kkt_multiply(x_scratch, linsys_data->xyzbuff2, work->data, Wfull, NULL, NULL,
-               xbuff, ubuff1, ubuff2);
+  kkt_multiply(x_scratch, linsys_data->xyzbuff2, work->data, nt_scaling,
+               nt_scaling_soc_idx, soc_idx, xbuff, ubuff1, ubuff2);
   for (QOCOInt k = 0; k < n; ++k) {
     linsys_data->xyzbuff2[k] -= linsys_data->kkt_static_reg_P * x_scratch[k];
   }
diff --git a/algebra/cuda/cudss_backend.cu b/algebra/cuda/cudss_backend.cu
index 00067c6..ea0df3c 100644
--- a/algebra/cuda/cudss_backend.cu
+++ b/algebra/cuda/cudss_backend.cu
@@ -715,8 +715,8 @@ static QOCOFloat compute_linsys_residual(LinSysData* linsys_data,
                                          const QOCOFloat* x,
                                          QOCOFloat* residual_scratch)
 {
-  QOCOFloat* Wfull = get_data_vectorf(work->Wfull);
-  QOCOInt* Wsoc_idx = get_data_vectori(work->Wsoc_idx);
+  QOCOFloat* nt_scaling = get_data_vectorf(work->nt_scaling);
+  QOCOInt* nt_scaling_soc_idx = get_data_vectori(work->nt_scaling_soc_idx);
   QOCOInt* soc_idx = get_data_vectori(work->soc_idx);
   QOCOFloat* xbuff = get_data_vectorf(work->xbuff);
   QOCOFloat* ubuff1 = get_data_vectorf(work->ubuff1);
@@ -726,8 +726,8 @@ static QOCOFloat compute_linsys_residual(LinSysData* linsys_data,
 
   // d_rhs_matrix_data is scratch here; cudss_solve_system overwrites it before
   // every cuDSS solve.
-  kkt_multiply((QOCOFloat*)x, linsys_data->d_rhs_matrix_data, work->data, Wfull,
-               Wsoc_idx, soc_idx, xbuff, ubuff1, ubuff2);
+  kkt_multiply((QOCOFloat*)x, linsys_data->d_rhs_matrix_data, work->data,
+               nt_scaling, nt_scaling_soc_idx, soc_idx, xbuff, ubuff1, ubuff2);
 
   // data->P stores P + eps_P * I, so remove the P regularization from the
   // product before measuring the true KKT residual.
diff --git a/include/cone.h b/include/cone.h
index ec16a38..f34c6ae 100644
--- a/include/cone.h
+++ b/include/cone.h
@@ -23,16 +23,17 @@ extern "C" {
 #include "structs.h"
 
 /**
- * @brief Sets Wfull to I.
+ * @brief Sets NT scaling data to identity.
  *
- * @param Wfull Full NT scaling matrix.
- * @param Wnnzfull Number of elements in Wfull.
- * @param Wsoc_idx Vector pointing to the start of each SOC block in Wfull.
+ * @param nt_scaling NT scaling data.
+ * @param nt_scaling_nnz Number of elements in nt_scaling.
+ * @param nt_scaling_soc_idx Vector pointing to the start of each SOC block.
  *                 Only used in the GPU (cone.cu) implementation.
  * @param data Pointer to problem data.
  */
-void set_Wfull_identity(QOCOVectorf* Wfull, QOCOInt Wnnzfull,
-                        QOCOVectori* Wsoc_idx, QOCOProblemData* data);
+void set_nt_scaling_identity(QOCOVectorf* nt_scaling, QOCOInt nt_scaling_nnz,
+                             QOCOVectori* nt_scaling_soc_idx,
+                             QOCOProblemData* data);
 
 /**
  * @brief Computes cone product u * v = p with respect to C.
@@ -79,13 +80,13 @@ void cone_division(const QOCOFloat* lambda, const QOCOFloat* v, QOCOFloat* d,
 void bring2cone(QOCOFloat* u, QOCOInt* soc_idx, QOCOProblemData* data);
 
 /**
- * @brief Computes z = W * x where W is a full Nesterov-Todd scaling matrix.
- * The NT scaling array for the LP cones are stored first, then the NT
- * scalings for the second-order cones are stored in column major order.
+ * @brief Computes z = W * x using Nesterov-Todd scaling data.
+ * The LP cone scalings are stored first. Each SOC block stores the fast
+ * scaling parameters [eta, w0, w1...] in the builtin backend.
  *
- * @param W Nesterov Todd scaling matrix.
- * @param Wsoc_idx Vector pointing to the start of each SOC block in W.
- *                 Only used in the GPU (cone.cu) implementation.
+ * @param W Nesterov Todd scaling data.
+ * @param nt_scaling_soc_idx Vector pointing to the start of each SOC block in
+ * W. Only used in the GPU (cone.cu) implementation.
  * @param soc_idx Array pointing to the start of each SOC block in x and z.
  *                Only used in the GPU (cone.cu) implementation.
  * @param x Input vector.
@@ -95,10 +96,18 @@ void bring2cone(QOCOFloat* u, QOCOInt* soc_idx, QOCOProblemData* data);
  * @param nsoc Number of second-order cones in C.
  * @param q Array of second-order cone dimensions.
  */
-void nt_multiply(QOCOFloat* W, QOCOInt* Wsoc_idx, QOCOInt* soc_idx,
+void nt_multiply(QOCOFloat* W, QOCOInt* nt_scaling_soc_idx, QOCOInt* soc_idx,
                  QOCOFloat* x, QOCOFloat* z, QOCOInt l, QOCOInt m, QOCOInt nsoc,
                  QOCOInt* q);
 
+/**
+ * @brief Computes z = W^{-1} * x using the same NT scaling data as
+ * nt_multiply().
+ */
+void nt_multiply_inv(QOCOFloat* W, QOCOInt* nt_scaling_soc_idx,
+                     QOCOInt* soc_idx, QOCOFloat* x, QOCOFloat* z, QOCOInt l,
+                     QOCOInt m, QOCOInt nsoc, QOCOInt* q);
+
 /**
  * @brief Compute Nesterov-Todd scalings and scaled variables.
  *
diff --git a/include/kkt.h b/include/kkt.h
index c96049f..9c34966 100644
--- a/include/kkt.h
+++ b/include/kkt.h
@@ -152,14 +152,16 @@ void predictor_corrector(QOCOSolver* solver);
  * @param x Pointer to input vector.
  * @param y Pointer to output vector.
  * @param data Pointer to problem data.
- * @param Wfull Pointer to full NT scaling matrix W.
- * @param Wsoc_idx Vector pointing to the start of each SOC block in Wfull.
+ * @param nt_scaling Pointer to NT scaling data.
+ * @param nt_scaling_soc_idx Vector pointing to the start of each SOC block in
+ * nt_scaling.
  * @param soc_idx Array pointing to the start of each SOC block in x and y.
  * @param nbuff Temporary buffer of length n.
  * @param mbuff1 Temporary buffer of length m.
  * @param mbuff2 Temporary buffer of length m.
  */
 void kkt_multiply(QOCOFloat* x, QOCOFloat* y, QOCOProblemData* data,
-                  QOCOFloat* Wfull, QOCOInt* Wsoc_idx, QOCOInt* soc_idx,
-                  QOCOFloat* nbuff, QOCOFloat* mbuff1, QOCOFloat* mbuff2);
+                  QOCOFloat* nt_scaling, QOCOInt* nt_scaling_soc_idx,
+                  QOCOInt* soc_idx, QOCOFloat* nbuff, QOCOFloat* mbuff1,
+                  QOCOFloat* mbuff2);
 #endif /* #ifndef QOCO_KKT_H */
\ No newline at end of file
diff --git a/include/structs.h b/include/structs.h
index ad8f133..8e0eda4 100644
--- a/include/structs.h
+++ b/include/structs.h
@@ -224,26 +224,26 @@ typedef struct {
   /** Number of nonzeros in upper triangular part of Nesterov-Todd Scaling. */
   QOCOInt Wnnz;
 
-  /** Number of nonzeros in full Nesterov-Todd Scaling. */
-  QOCOInt Wnnzfull;
+  /** Number of entries in nt_scaling. */
+  QOCOInt nt_scaling_nnz;
 
   /** Upper triangular part of Nesterov-Todd Scaling */
   QOCOVectorf* W;
 
-  /** Full Nesterov-Todd Scaling */
-  QOCOVectorf* Wfull;
+  /** NT scaling data used by nt_multiply(). Shared layout across backends:
+   *   LP entries: scalar scales sqrt(s_i / z_i), length l.
+   *   SOC i block: [eta, w0, w1[0], ..., w1[q_i - 2]], length q_i + 1.
+   */
+  QOCOVectorf* nt_scaling;
 
   /** Upper triangular part of inverse of Nesterov-Todd Scaling */
   QOCOVectorf* Winv;
 
-  /** Full inverse of Nesterov-Todd Scaling */
-  QOCOVectorf* Winvfull;
-
   /** Nesterov-Todd Scaling squared */
   QOCOVectorf* WtW;
 
-  /** Vector which points to the start of the ith soc block in Wfull */
-  QOCOVectori* Wsoc_idx;
+  /** Vector which points to the start of the ith SOC block in nt_scaling. */
+  QOCOVectori* nt_scaling_soc_idx;
 
   /** Vector which points to the start of the start of the ith soc variable
    * block */
diff --git a/src/cone.c b/src/cone.c
index 833f97b..0d8dfcb 100644
--- a/src/cone.c
+++ b/src/cone.c
@@ -11,24 +11,24 @@
 #include "cone.h"
 #include "qoco_utils.h"
 
-void set_Wfull_identity(QOCOVectorf* Wfull, QOCOInt Wnnzfull,
-                        QOCOVectori* Wsoc_idx, QOCOProblemData* data)
+void set_nt_scaling_identity(QOCOVectorf* nt_scaling, QOCOInt nt_scaling_nnz,
+                             QOCOVectori* nt_scaling_soc_idx,
+                             QOCOProblemData* data)
 {
-  (void)Wsoc_idx;
-  QOCOFloat* Wfull_data = get_data_vectorf(Wfull);
-  for (QOCOInt i = 0; i < Wnnzfull; ++i) {
-    Wfull_data[i] = 0.0;
+  (void)nt_scaling_soc_idx;
+  QOCOFloat* nt_scaling_data = get_data_vectorf(nt_scaling);
+  for (QOCOInt i = 0; i < nt_scaling_nnz; ++i) {
+    nt_scaling_data[i] = 0.0;
   }
   for (QOCOInt i = 0; i < data->l; ++i) {
-    Wfull_data[i] = 1.0;
+    nt_scaling_data[i] = 1.0;
   }
   QOCOInt idx = data->l;
   for (QOCOInt i = 0; i < data->nsoc; ++i) {
     QOCOInt qi = get_element_vectori(data->q, i);
-    for (QOCOInt k = 0; k < qi; ++k) {
-      Wfull_data[idx + k * qi + k] = 1.0;
-    }
-    idx += qi * qi;
+    nt_scaling_data[idx] = 1.0;
+    nt_scaling_data[idx + 1] = 1.0;
+    idx += qi + 1;
   }
 }
 
@@ -212,44 +212,71 @@ void bring2cone(QOCOFloat* u, QOCOInt* soc_idx, QOCOProblemData* data)
   }
 }
 
-void nt_multiply(QOCOFloat* W, QOCOInt* Wsoc_idx, QOCOInt* soc_idx,
-                 QOCOFloat* x, QOCOFloat* z, QOCOInt l, QOCOInt m, QOCOInt nsoc,
-                 QOCOInt* q)
+static void nt_multiply_impl(QOCOFloat* W, QOCOInt* nt_scaling_soc_idx,
+                             QOCOInt* soc_idx, QOCOFloat* x, QOCOFloat* z,
+                             QOCOInt l, QOCOInt m, QOCOInt nsoc, QOCOInt* q,
+                             QOCOInt inverse)
 {
-  (void)Wsoc_idx;
+  (void)nt_scaling_soc_idx;
   (void)soc_idx;
+  (void)m;
   // Compute product for LP cone part of W.
   for (QOCOInt i = 0; i < l; ++i) {
-    z[i] = (W[i] * x[i]);
+    z[i] = inverse ? (safe_div(1.0, W[i]) * x[i]) : (W[i] * x[i]);
   }
 
-  // Compute product for second-order cones.
+  // Compute product for second-order cones using fast O(m) operations
+  // from equations (14) and (15) in the ECOS paper.
   QOCOInt nt_idx = l;
   QOCOInt idx = l;
 
-  // Zero out second-order cone block of result z.
-  for (QOCOInt i = l; i < m; ++i) {
-    z[i] = 0;
-  }
-
-  // Loop over all second-order cones.
   for (QOCOInt i = 0; i < nsoc; ++i) {
-    // Loop over elements within a second-order cone.
-    for (QOCOInt j = 0; j < q[i]; ++j) {
-      z[idx + j] += qoco_dot(&W[nt_idx + j * q[i]], &x[idx], q[i]);
+    QOCOFloat scale = inverse ? safe_div(1.0, W[nt_idx]) : W[nt_idx];
+    QOCOFloat w0 = W[nt_idx + 1];
+    QOCOFloat* w1 = &W[nt_idx + 2];
+    QOCOFloat x0 = x[idx];
+    QOCOFloat zeta = qoco_dot(w1, &x[idx + 1], q[i] - 1);
+    QOCOFloat w0p1_inv = safe_div(1.0, 1.0 + w0);
+
+    if (inverse) {
+      z[idx] = scale * (w0 * x0 - zeta);
+      QOCOFloat coeff = -x0 + zeta * w0p1_inv;
+      for (QOCOInt j = 1; j < q[i]; ++j) {
+        z[idx + j] = scale * (x[idx + j] + coeff * w1[j - 1]);
+      }
+    }
+    else {
+      z[idx] = scale * (w0 * x0 + zeta);
+      QOCOFloat coeff = x0 + zeta * w0p1_inv;
+      for (QOCOInt j = 1; j < q[i]; ++j) {
+        z[idx + j] = scale * (x[idx + j] + coeff * w1[j - 1]);
+      }
     }
     idx += q[i];
-    nt_idx += q[i] * q[i];
+    nt_idx += q[i] + 1;
   }
 }
 
+void nt_multiply(QOCOFloat* W, QOCOInt* nt_scaling_soc_idx, QOCOInt* soc_idx,
+                 QOCOFloat* x, QOCOFloat* z, QOCOInt l, QOCOInt m, QOCOInt nsoc,
+                 QOCOInt* q)
+{
+  nt_multiply_impl(W, nt_scaling_soc_idx, soc_idx, x, z, l, m, nsoc, q, 0);
+}
+
+void nt_multiply_inv(QOCOFloat* W, QOCOInt* nt_scaling_soc_idx,
+                     QOCOInt* soc_idx, QOCOFloat* x, QOCOFloat* z, QOCOInt l,
+                     QOCOInt m, QOCOInt nsoc, QOCOInt* q)
+{
+  nt_multiply_impl(W, nt_scaling_soc_idx, soc_idx, x, z, l, m, nsoc, q, 1);
+}
+
 void compute_nt_scaling(QOCOWorkspace* work)
 {
   QOCOFloat* W = get_data_vectorf(work->W);
   QOCOFloat* WtW = get_data_vectorf(work->WtW);
-  QOCOFloat* Wfull = get_data_vectorf(work->Wfull);
+  QOCOFloat* nt_scaling = get_data_vectorf(work->nt_scaling);
   QOCOFloat* Winv = get_data_vectorf(work->Winv);
-  QOCOFloat* Winvfull = get_data_vectorf(work->Winvfull);
   QOCOFloat* sbar = get_data_vectorf(work->sbar);
   QOCOFloat* zbar = get_data_vectorf(work->zbar);
   QOCOFloat* lambda = get_data_vectorf(work->lambda);
@@ -259,17 +286,16 @@ void compute_nt_scaling(QOCOWorkspace* work)
     WtW[idx] = safe_div(get_element_vectorf(work->s, idx),
                         get_element_vectorf(work->z, idx));
     W[idx] = qoco_sqrt(WtW[idx]);
-    Wfull[idx] = W[idx];
+    nt_scaling[idx] = W[idx];
     Winv[idx] = safe_div(1.0, W[idx]);
-    Winvfull[idx] = Winv[idx];
   }
 
   // Compute Nesterov-Todd scaling for second-order cones.
   QOCOInt nt_idx = idx;
-  QOCOInt nt_idx_full = idx;
+  QOCOInt nt_idx_fast = idx;
   for (QOCOInt i = 0; i < work->data->nsoc; ++i) {
     QOCOInt qi = get_element_vectori(work->data->q, i);
-    // Compute normalized vectors.
+    // Compute normalized vectors. Bottom of page 8 in coneprog.
     QOCOFloat s_scal = soc_residual2(get_pointer_vectorf(work->s, idx), qi);
     s_scal = qoco_sqrt(s_scal);
     QOCOFloat f = safe_div(1.0, s_scal);
@@ -284,27 +310,47 @@ void compute_nt_scaling(QOCOWorkspace* work)
 
     f = safe_div(1.0, (2 * gamma));
 
-    // Overwrite sbar with wbar.
+    // Overwrite sbar with wbar. Eq (15) in coneprog.
     sbar[0] = f * (sbar[0] + zbar[0]);
     for (QOCOInt j = 1; j < qi; ++j) {
       sbar[j] = f * (sbar[j] - zbar[j]);
     }
 
-    // Overwrite zbar with v.
-    f = safe_div(1.0, qoco_sqrt(2 * (sbar[0] + 1)));
-    zbar[0] = f * (sbar[0] + 1.0);
+    // eta = sqrt(s_scal / z_scal). Eq (7) in ecos.
+    f = qoco_sqrt(safe_div(s_scal, z_scal));
+
+    // Store fast scaling parameters: [eta, w0, w1[0], ..., w1[qi-2]]
+    nt_scaling[nt_idx_fast] = f;
+    nt_scaling[nt_idx_fast + 1] = sbar[0];
+    for (QOCOInt j = 1; j < qi; ++j) {
+      nt_scaling[nt_idx_fast + 1 + j] = sbar[j];
+    }
+
+    // Compute W (upper triangular) and Winv for the sparse KKT update.
+    // Also compute WtW = eta^2 * (2*w*w' - J) in upper triangular storage.
+    QOCOFloat finv = safe_div(1.0, f);
+    QOCOFloat eta2 = f * f;
+
+    // W is given by Eq (7) in ecos, but in terms of wbar that has a block
+    // structure (scalar, first row/column, and bottom block) with differing
+    // coefficients, so it cannot be built with a single uniform outer product.
+    // We instead introduce v so that W = eta * (2 * v * v' - J), which has the
+    // same rank-1-minus-J form as WtW = eta^2 * (2 * wbar * wbar' - J). This
+    // lets W, Winv, and WtW share one outer-product loop, and makes Winv just a
+    // sign flip of the cross terms of W. v is defined as
+    //   v0 = (wbar0 + 1) / sqrt(2 * (wbar0 + 1)),
+    //   vj = wbar_j   / sqrt(2 * (wbar0 + 1)) for j >= 1,
+    // i.e. wbar with its leading entry shifted by 1 and the whole vector scaled
+    // by 1 / sqrt(2 * (wbar0 + 1))..
+    QOCOFloat fv = safe_div(1.0, qoco_sqrt(2 * (sbar[0] + 1)));
+    zbar[0] = fv * (sbar[0] + 1.0);
     for (QOCOInt j = 1; j < qi; ++j) {
-      zbar[j] = f * sbar[j];
+      zbar[j] = fv * sbar[j];
     }
 
-    // Compute W for second-order cones.
     QOCOInt shift = 0;
-    f = qoco_sqrt(safe_div(s_scal, z_scal));
-    QOCOFloat finv = safe_div(1.0, f);
     for (QOCOInt j = 0; j < qi; ++j) {
       for (QOCOInt k = 0; k <= j; ++k) {
-        QOCOInt full_idx1 = nt_idx_full + j * qi + k;
-        QOCOInt full_idx2 = nt_idx_full + k * qi + j;
         W[nt_idx + shift] = 2 * (zbar[k] * zbar[j]);
         if (j != 0 && k == 0) {
           Winv[nt_idx + shift] = -W[nt_idx + shift];
@@ -322,31 +368,26 @@ void compute_nt_scaling(QOCOWorkspace* work)
         }
         W[nt_idx + shift] *= f;
         Winv[nt_idx + shift] *= finv;
-        Wfull[full_idx1] = W[nt_idx + shift];
-        Wfull[full_idx2] = W[nt_idx + shift];
-        Winvfull[full_idx1] = Winv[nt_idx + shift];
-        Winvfull[full_idx2] = Winv[nt_idx + shift];
-        shift += 1;
-      }
-    }
 
-    // Compute WtW for second-order cones.
-    shift = 0;
-    for (QOCOInt j = 0; j < qi; ++j) {
-      for (QOCOInt k = 0; k <= j; ++k) {
-        WtW[nt_idx + shift] = qoco_dot(&Wfull[nt_idx_full + j * qi],
-                                       &Wfull[nt_idx_full + k * qi], qi);
+        QOCOFloat val = eta2 * 2.0 * sbar[j] * sbar[k];
+        if (j == k && j == 0) {
+          val -= eta2;
+        }
+        else if (j == k) {
+          val += eta2;
+        }
+        WtW[nt_idx + shift] = val;
         shift += 1;
       }
     }
 
     idx += qi;
     nt_idx += (qi * qi + qi) / 2;
-    nt_idx_full += qi * qi;
+    nt_idx_fast += qi + 1;
   }
 
   // Compute scaled variable lambda. lambda = W * z.
-  nt_multiply(Wfull, NULL, NULL, get_pointer_vectorf(work->z, 0), lambda,
+  nt_multiply(nt_scaling, NULL, NULL, get_pointer_vectorf(work->z, 0), lambda,
               work->data->l, work->data->m, work->data->nsoc,
               get_data_vectori(work->data->q));
 }
@@ -557,4 +598,4 @@ void add_e(QOCOFloat* x, QOCOFloat a, QOCOInt l, QOCOInt nsoc, QOCOVectori* q)
     x[idx] -= a;
     idx += get_element_vectori(q, i);
   }
-}
\ No newline at end of file
+}
diff --git a/src/cone.cu b/src/cone.cu
index 1dfda40..5b9a07b 100644
--- a/src/cone.cu
+++ b/src/cone.cu
@@ -85,11 +85,12 @@ __device__ void soc_division(const QOCOFloat* lam, const QOCOFloat* v,
   }
 }
 
-__global__ void set_Wfull_linear(QOCOFloat* W, QOCOInt Wnnzfull, QOCOInt l)
+__global__ void set_nt_scaling_linear(QOCOFloat* W, QOCOInt nt_scaling_nnz,
+                                      QOCOInt l)
 {
   QOCOInt i = blockIdx.x * blockDim.x + threadIdx.x;
 
-  if (i < Wnnzfull) {
+  if (i < nt_scaling_nnz) {
     W[i] = 0.0;
   }
   if (i < l) {
@@ -97,21 +98,24 @@ __global__ void set_Wfull_linear(QOCOFloat* W, QOCOInt Wnnzfull, QOCOInt l)
   }
 }
 
-__global__ void set_Wfull_soc(QOCOFloat* W, QOCOInt* q, QOCOInt* Wsoc_idx,
-                              QOCOInt nsoc, QOCOInt l)
+__global__ void set_nt_scaling_soc(QOCOFloat* W, QOCOInt* q,
+                                   QOCOInt* nt_scaling_soc_idx, QOCOInt nsoc,
+                                   QOCOInt l)
 {
+  (void)q;
+  (void)l;
   QOCOInt soc = blockIdx.x;
   if (soc >= nsoc)
     return;
 
-  QOCOInt dim = q[soc];
-  QOCOInt k = threadIdx.x;
-
-  if (k >= dim)
+  // Only one thread per SOC writes the compact identity scaling
+  // [eta = 1, w0 = 1, w1 = 0, ...].
+  if (threadIdx.x != 0)
     return;
 
-  // diagonal element
-  W[Wsoc_idx[soc] + k * dim + k] = 1.0;
+  QOCOInt base = nt_scaling_soc_idx[soc];
+  W[base] = 1.0;     // eta
+  W[base + 1] = 1.0; // w0
 }
 
 __global__ void cone_residual_stage1(const QOCOFloat* u, QOCOInt l,
@@ -218,8 +222,8 @@ __global__ void bring2cone_kernel(QOCOFloat* u, QOCOInt* q, QOCOInt l,
 }
 
 __global__ void compute_nt_scaling_kernel(QOCOFloat* W, QOCOFloat* WtW,
-                                          QOCOFloat* Wfull, QOCOFloat* Winv,
-                                          QOCOFloat* Winvfull, QOCOFloat* s,
+                                          QOCOFloat* nt_scaling,
+                                          QOCOFloat* Winv, QOCOFloat* s,
                                           QOCOFloat* z, QOCOFloat* sbar,
                                           QOCOFloat* zbar, QOCOInt l,
                                           QOCOInt nsoc, const QOCOInt* q)
@@ -233,11 +237,10 @@ __global__ void compute_nt_scaling_kernel(QOCOFloat* W, QOCOFloat* WtW,
     QOCOFloat w = qoco_sqrt(val);
 
     W[tid] = w;
-    Wfull[tid] = w;
+    nt_scaling[tid] = w;
 
     QOCOFloat winv = safe_div((QOCOFloat)1.0, w);
     Winv[tid] = winv;
-    Winvfull[tid] = winv;
     return;
   }
 
@@ -249,12 +252,12 @@ __global__ void compute_nt_scaling_kernel(QOCOFloat* W, QOCOFloat* WtW,
   /* ---- compute SOC offsets ---- */
   QOCOInt idx = l;
   QOCOInt nt_idx = l;
-  QOCOInt nt_idx_full = l;
+  QOCOInt nt_idx_fast = l;
 
   for (QOCOInt k = 0; k < soc_id; ++k) {
     idx += q[k];
     nt_idx += (q[k] * q[k] + q[k]) / 2;
-    nt_idx_full += q[k] * q[k];
+    nt_idx_fast += q[k] + 1;
   }
 
   QOCOInt qi = q[soc_id];
@@ -284,7 +287,19 @@ __global__ void compute_nt_scaling_kernel(QOCOFloat* W, QOCOFloat* WtW,
   for (QOCOInt j = 1; j < qi; ++j)
     sbar[idx + j] = f * (sbar[idx + j] - zbar[idx + j]);
 
-  /* overwrite zbar with v */
+  /* eta = sqrt(s_scal / z_scal) */
+  QOCOFloat eta = qoco_sqrt(safe_div(s_scal, z_scal));
+  QOCOFloat finv = safe_div((QOCOFloat)1.0, eta);
+  QOCOFloat eta2 = eta * eta;
+
+  /* Store compact fast scaling parameters [eta, w0, w1[0], ..., w1[qi-2]].
+   * sbar currently holds the wbar vector. */
+  nt_scaling[nt_idx_fast] = eta;
+  nt_scaling[nt_idx_fast + 1] = sbar[idx + 0];
+  for (QOCOInt j = 1; j < qi; ++j)
+    nt_scaling[nt_idx_fast + 1 + j] = sbar[idx + j];
+
+  /* overwrite zbar with v (needed for the sparse W / Winv blocks) */
   f = safe_div((QOCOFloat)1.0,
                qoco_sqrt((QOCOFloat)2.0 * (sbar[idx + 0] + (QOCOFloat)1.0)));
 
@@ -292,10 +307,8 @@ __global__ void compute_nt_scaling_kernel(QOCOFloat* W, QOCOFloat* WtW,
   for (QOCOInt j = 1; j < qi; ++j)
     zbar[idx + j] = f * sbar[idx + j];
 
-  /* --- build W and Winv --- */
-  QOCOFloat fwd = qoco_sqrt(safe_div(s_scal, z_scal));
-  QOCOFloat finv = 1.0 / fwd;
-
+  /* --- build W, Winv (sparse upper triangular) and WtW = eta^2 (2 w w' - J)
+   * --- */
   QOCOInt shift = 0;
   for (QOCOInt j = 0; j < qi; ++j) {
     for (QOCOInt k = 0; k <= j; ++k) {
@@ -315,56 +328,68 @@ __global__ void compute_nt_scaling_kernel(QOCOFloat* W, QOCOFloat* WtW,
         winv_val += (QOCOFloat)1.0;
       }
 
-      val *= fwd;
+      val *= eta;
       winv_val *= finv;
 
-      QOCOInt full1 = nt_idx_full + j * qi + k;
-      QOCOInt full2 = nt_idx_full + k * qi + j;
-
       W[nt_idx + shift] = val;
       Winv[nt_idx + shift] = winv_val;
 
-      Wfull[full1] = val;
-      Wfull[full2] = val;
-      Winvfull[full1] = winv_val;
-      Winvfull[full2] = winv_val;
+      QOCOFloat wtw = eta2 * (QOCOFloat)2.0 * sbar[idx + j] * sbar[idx + k];
+      if (j == k && j == 0)
+        wtw -= eta2;
+      else if (j == k)
+        wtw += eta2;
+      WtW[nt_idx + shift] = wtw;
 
       shift++;
     }
   }
-
-  /* --- compute WtW --- */
-  shift = 0;
-  for (QOCOInt j = 0; j < qi; ++j) {
-    for (QOCOInt k = 0; k <= j; ++k) {
-      WtW[nt_idx + shift] = qoco_dot_dev(&Wfull[nt_idx_full + j * qi],
-                                         &Wfull[nt_idx_full + k * qi], qi);
-      shift++;
-    }
-  }
 }
 
-__global__ void nt_multiply_kernel(const QOCOFloat* W, const QOCOInt* Wsoc_idx,
+__global__ void nt_multiply_kernel(const QOCOFloat* W,
+                                   const QOCOInt* nt_scaling_soc_idx,
                                    const QOCOInt* soc_idx, const QOCOFloat* x,
                                    QOCOFloat* z, QOCOInt l, QOCOInt m,
-                                   QOCOInt nsoc, const QOCOInt* q)
+                                   QOCOInt nsoc, const QOCOInt* q,
+                                   QOCOInt inverse)
 {
+  (void)m;
   QOCOInt i = blockIdx.x * blockDim.x + threadIdx.x;
   if (i >= l + nsoc)
     return;
 
   /* ================= LP cone ================= */
   if (i < l) {
-    z[i] = W[i] * x[i];
+    z[i] = inverse ? (safe_div((QOCOFloat)1.0, W[i]) * x[i]) : (W[i] * x[i]);
     return;
   }
 
-  /* ================= SOC cones ================= */
+  /* ================= SOC cones =================
+   * Fast O(q) product using the compact scaling [eta, w0, w1...] and
+   * equations (14) and (15) in the ECOS paper. */
   QOCOInt soc = i - l;
   QOCOInt qi = q[soc];
-  for (QOCOInt k = 0; k < qi; ++k) {
-    z[soc_idx[soc] + k] =
-        qoco_dot_dev(&W[Wsoc_idx[soc] + k * qi], &x[soc_idx[soc]], qi);
+  QOCOInt nt_idx = nt_scaling_soc_idx[soc];
+  QOCOInt xi = soc_idx[soc];
+
+  QOCOFloat scale = inverse ? safe_div((QOCOFloat)1.0, W[nt_idx]) : W[nt_idx];
+  QOCOFloat w0 = W[nt_idx + 1];
+  const QOCOFloat* w1 = &W[nt_idx + 2];
+  QOCOFloat x0 = x[xi];
+  QOCOFloat zeta = qoco_dot_dev(w1, &x[xi + 1], qi - 1);
+  QOCOFloat w0p1_inv = safe_div((QOCOFloat)1.0, (QOCOFloat)1.0 + w0);
+
+  if (inverse) {
+    z[xi] = scale * (w0 * x0 - zeta);
+    QOCOFloat coeff = -x0 + zeta * w0p1_inv;
+    for (QOCOInt j = 1; j < qi; ++j)
+      z[xi + j] = scale * (x[xi + j] + coeff * w1[j - 1]);
+  }
+  else {
+    z[xi] = scale * (w0 * x0 + zeta);
+    QOCOFloat coeff = x0 + zeta * w0p1_inv;
+    for (QOCOInt j = 1; j < qi; ++j)
+      z[xi + j] = scale * (x[xi + j] + coeff * w1[j - 1]);
   }
 }
 
@@ -439,27 +464,28 @@ __global__ void add_e_kernel(QOCOFloat* x, QOCOFloat a, QOCOInt l, QOCOInt nsoc,
   x[idx] -= a;
 }
 
-void set_Wfull_identity(QOCOVectorf* Wfull, QOCOInt Wnnzfull,
-                        QOCOVectori* Wsoc_idx, QOCOProblemData* data)
+void set_nt_scaling_identity(QOCOVectorf* nt_scaling, QOCOInt nt_scaling_nnz,
+                             QOCOVectori* nt_scaling_soc_idx,
+                             QOCOProblemData* data)
 {
   CUDA_CHECK(cudaGetLastError());
-  QOCOFloat* W = get_data_vectorf(Wfull);
+  QOCOFloat* W = get_data_vectorf(nt_scaling);
 
   const int threads = 256;
-  const int blocks = (Wnnzfull + threads - 1) / threads;
+  const int blocks = (nt_scaling_nnz + threads - 1) / threads;
 
   // kernel 1: zero + linear cone
   if (data->l > 0) {
-    set_Wfull_linear<<<blocks, threads>>>(W, Wnnzfull, data->l);
+    set_nt_scaling_linear<<<blocks, threads>>>(W, nt_scaling_nnz, data->l);
     CUDA_CHECK(cudaGetLastError());
   }
 
   // kernel 2: SOC blocks
   const int blocks2 = data->nsoc;
   if (data->nsoc > 0) {
-    set_Wfull_soc<<<blocks2, 256>>>(W, get_data_vectori(data->q),
-                                    get_data_vectori(Wsoc_idx), data->nsoc,
-                                    data->l);
+    set_nt_scaling_soc<<<blocks2, 256>>>(W, get_data_vectori(data->q),
+                                         get_data_vectori(nt_scaling_soc_idx),
+                                         data->nsoc, data->l);
     CUDA_CHECK(cudaGetLastError());
   }
 }
@@ -545,15 +571,28 @@ void bring2cone(QOCOFloat* u, QOCOInt* soc_idx, QOCOProblemData* data)
   }
 }
 
-void nt_multiply(QOCOFloat* W, QOCOInt* Wsoc_idx, QOCOInt* soc_idx,
+void nt_multiply(QOCOFloat* W, QOCOInt* nt_scaling_soc_idx, QOCOInt* soc_idx,
                  QOCOFloat* x, QOCOFloat* z, QOCOInt l, QOCOInt m, QOCOInt nsoc,
                  QOCOInt* q)
 {
   int threads = 256;
   int blocks = (l + nsoc + threads - 1) / threads;
   if (m > 0) {
-    nt_multiply_kernel<<<blocks, threads>>>(W, Wsoc_idx, soc_idx, x, z, l, m,
-                                            nsoc, q);
+    nt_multiply_kernel<<<blocks, threads>>>(W, nt_scaling_soc_idx, soc_idx, x,
+                                            z, l, m, nsoc, q, 0);
+  }
+  CUDA_CHECK(cudaGetLastError());
+}
+
+void nt_multiply_inv(QOCOFloat* W, QOCOInt* nt_scaling_soc_idx,
+                     QOCOInt* soc_idx, QOCOFloat* x, QOCOFloat* z, QOCOInt l,
+                     QOCOInt m, QOCOInt nsoc, QOCOInt* q)
+{
+  int threads = 256;
+  int blocks = (l + nsoc + threads - 1) / threads;
+  if (m > 0) {
+    nt_multiply_kernel<<<blocks, threads>>>(W, nt_scaling_soc_idx, soc_idx, x,
+                                            z, l, m, nsoc, q, 1);
   }
   CUDA_CHECK(cudaGetLastError());
 }
@@ -562,10 +601,9 @@ void compute_nt_scaling(QOCOWorkspace* work)
 {
   QOCOFloat* W = get_data_vectorf(work->W);
   QOCOFloat* WtW = get_data_vectorf(work->WtW);
-  QOCOFloat* Wfull = get_data_vectorf(work->Wfull);
+  QOCOFloat* nt_scaling = get_data_vectorf(work->nt_scaling);
   QOCOFloat* Winv = get_data_vectorf(work->Winv);
-  QOCOFloat* Winvfull = get_data_vectorf(work->Winvfull);
-  QOCOInt* Wsoc_idx = get_data_vectori(work->Wsoc_idx);
+  QOCOInt* nt_scaling_soc_idx = get_data_vectori(work->nt_scaling_soc_idx);
   QOCOInt* soc_idx = get_data_vectori(work->soc_idx);
   QOCOFloat* s = get_data_vectorf(work->s);
   QOCOFloat* z = get_data_vectorf(work->z);
@@ -581,14 +619,14 @@ void compute_nt_scaling(QOCOWorkspace* work)
   QOCOInt grid = (total_threads + block - 1) / block;
 
   if (work->data->m > 0) {
-    compute_nt_scaling_kernel<<<grid, block>>>(W, WtW, Wfull, Winv, Winvfull, s,
-                                               z, sbar, zbar, l, nsoc, q);
+    compute_nt_scaling_kernel<<<grid, block>>>(W, WtW, nt_scaling, Winv, s, z,
+                                               sbar, zbar, l, nsoc, q);
   }
   CUDA_CHECK(cudaGetLastError());
 
   /* ================= lambda = W * z ================= */
-  nt_multiply(Wfull, Wsoc_idx, soc_idx, z, lambda, work->data->l, work->data->m,
-              work->data->nsoc, q);
+  nt_multiply(nt_scaling, nt_scaling_soc_idx, soc_idx, z, lambda, work->data->l,
+              work->data->m, work->data->nsoc, q);
 }
 
 void compute_centering(QOCOSolver* solver)
diff --git a/src/kkt.c b/src/kkt.c
index 0d35629..d119b4c 100644
--- a/src/kkt.c
+++ b/src/kkt.c
@@ -172,9 +172,10 @@ void initialize_ipm(QOCOSolver* solver)
   QOCOWorkspace* work = solver->work;
   QOCOProblemData* data = solver->work->data;
 
-  // Set Nesterov-Todd block in Wfull to -I (need for kkt_multiply in iterative
-  // refinement).
-  set_Wfull_identity(work->Wfull, work->Wnnzfull, work->Wsoc_idx, data);
+  // Set Nesterov-Todd scaling data to identity (needed by kkt_multiply in
+  // iterative refinement).
+  set_nt_scaling_identity(work->nt_scaling, work->nt_scaling_nnz,
+                          work->nt_scaling_soc_idx, data);
   solver->linsys->linsys_set_nt_identity(solver->linsys_data, data->m);
 
   // Needs to be set to 1.0 not 0.0 due to low tolerance stopping criteria
@@ -309,8 +310,8 @@ void construct_kkt_aff_rhs(QOCOWorkspace* work)
 {
   QOCOFloat* rhs = get_data_vectorf(work->rhs);
   QOCOFloat* kktres = get_data_vectorf(work->kktres);
-  QOCOFloat* Wfull = get_data_vectorf(work->Wfull);
-  QOCOInt* Wsoc_idx = get_data_vectori(work->Wsoc_idx);
+  QOCOFloat* nt_scaling = get_data_vectorf(work->nt_scaling);
+  QOCOInt* nt_scaling_soc_idx = get_data_vectori(work->nt_scaling_soc_idx);
   QOCOInt* soc_idx = get_data_vectori(work->soc_idx);
   QOCOFloat* lambda = get_data_vectorf(work->lambda);
   QOCOFloat* ubuff1 = get_data_vectorf(work->ubuff1);
@@ -321,8 +322,8 @@ void construct_kkt_aff_rhs(QOCOWorkspace* work)
                          work->data->n + work->data->p + work->data->m);
 
   // Compute W*lambda
-  nt_multiply(Wfull, Wsoc_idx, soc_idx, lambda, ubuff1, work->data->l,
-              work->data->m, work->data->nsoc, q);
+  nt_multiply(nt_scaling, nt_scaling_soc_idx, soc_idx, lambda, ubuff1,
+              work->data->l, work->data->m, work->data->nsoc, q);
 
   // Add W*lambda to z portion of rhs.
   qoco_axpy(ubuff1, &rhs[work->data->n + work->data->p],
@@ -334,9 +335,8 @@ void construct_kkt_comb_rhs(QOCOWorkspace* work)
   QOCOFloat* rhs = get_data_vectorf(work->rhs);
   QOCOFloat* xyz = get_data_vectorf(work->xyz);
   QOCOFloat* kktres = get_data_vectorf(work->kktres);
-  QOCOFloat* Wfull = get_data_vectorf(work->Wfull);
-  QOCOFloat* Winvfull = get_data_vectorf(work->Winvfull);
-  QOCOInt* Wsoc_idx = get_data_vectori(work->Wsoc_idx);
+  QOCOFloat* nt_scaling = get_data_vectorf(work->nt_scaling);
+  QOCOInt* nt_scaling_soc_idx = get_data_vectori(work->nt_scaling_soc_idx);
   QOCOInt* soc_idx = get_data_vectori(work->soc_idx);
   QOCOFloat* Ds = get_data_vectorf(work->Ds);
   QOCOFloat* lambda = get_data_vectorf(work->lambda);
@@ -353,13 +353,13 @@ void construct_kkt_comb_rhs(QOCOWorkspace* work)
   /// cone_product((W' \ Dsaff), (W * Dzaff), pdata) + sigma * mu * e.
 
   // ubuff1 = Winv * Dsaff.
-  nt_multiply(Winvfull, Wsoc_idx, soc_idx, Ds, ubuff1, work->data->l,
-              work->data->m, work->data->nsoc, q);
+  nt_multiply_inv(nt_scaling, nt_scaling_soc_idx, soc_idx, Ds, ubuff1,
+                  work->data->l, work->data->m, work->data->nsoc, q);
 
   // ubuff2 = W * Dzaff.
   QOCOFloat* Dzaff = &xyz[work->data->n + work->data->p];
-  nt_multiply(Wfull, Wsoc_idx, soc_idx, Dzaff, ubuff2, work->data->l,
-              work->data->m, work->data->nsoc, q);
+  nt_multiply(nt_scaling, nt_scaling_soc_idx, soc_idx, Dzaff, ubuff2,
+              work->data->l, work->data->m, work->data->nsoc, q);
 
   // ubuff3 = cone_product((W' \ Dsaff), (W * Dzaff), pdata).
   cone_product(ubuff1, ubuff2, ubuff3, work->data->l, work->data->nsoc, q,
@@ -384,8 +384,8 @@ void construct_kkt_comb_rhs(QOCOWorkspace* work)
                 soc_idx);
 
   // ubuff1 = W * cone_division(lambda, ds).
-  nt_multiply(Wfull, Wsoc_idx, soc_idx, ubuff2, ubuff1, work->data->l,
-              work->data->m, work->data->nsoc, q);
+  nt_multiply(nt_scaling, nt_scaling_soc_idx, soc_idx, ubuff2, ubuff1,
+              work->data->l, work->data->m, work->data->nsoc, q);
 
   // rhs = [dx;dy;dz-W'*cone_division(lambda, ds, pdata)];
   qoco_axpy(ubuff1, &rhs[work->data->n + work->data->p],
@@ -397,8 +397,8 @@ void predictor_corrector(QOCOSolver* solver)
   QOCOWorkspace* work = solver->work;
   QOCOProblemData* data = solver->work->data;
 
-  QOCOFloat* Wfull = get_data_vectorf(work->Wfull);
-  QOCOInt* Wsoc_idx = get_data_vectori(work->Wsoc_idx);
+  QOCOFloat* nt_scaling = get_data_vectorf(work->nt_scaling);
+  QOCOInt* nt_scaling_soc_idx = get_data_vectori(work->nt_scaling_soc_idx);
   QOCOInt* soc_idx = get_data_vectori(work->soc_idx);
   QOCOFloat* lambda = get_data_vectorf(work->lambda);
   QOCOFloat* Ds = get_data_vectorf(work->Ds);
@@ -421,12 +421,12 @@ void predictor_corrector(QOCOSolver* solver)
   // Compute Dsaff. Dsaff = W' * (-lambda - W * Dzaff).
   QOCOFloat* xyz = get_data_vectorf(work->xyz);
   QOCOFloat* Dzaff = &xyz[data->n + data->p];
-  nt_multiply(Wfull, Wsoc_idx, soc_idx, Dzaff, ubuff1, data->l, data->m,
-              data->nsoc, q);
+  nt_multiply(nt_scaling, nt_scaling_soc_idx, soc_idx, Dzaff, ubuff1, data->l,
+              data->m, data->nsoc, q);
   copy_and_negate_arrayf(ubuff1, ubuff1, data->m);
   qoco_axpy(lambda, ubuff1, ubuff1, -1.0, data->m);
-  nt_multiply(Wfull, Wsoc_idx, soc_idx, ubuff1, Ds, data->l, data->m,
-              data->nsoc, q);
+  nt_multiply(nt_scaling, nt_scaling_soc_idx, soc_idx, ubuff1, Ds, data->l,
+              data->m, data->nsoc, q);
 
   // Compute centering parameter.
   compute_centering(solver);
@@ -450,12 +450,12 @@ void predictor_corrector(QOCOSolver* solver)
   // computed in construct_kkt_comb_rhs() and stored in work->Ds.
   QOCOFloat* Dz = &xyz[data->n + data->p];
   cone_division(lambda, Ds, ubuff1, data->l, data->nsoc, q, soc_idx);
-  nt_multiply(Wfull, Wsoc_idx, soc_idx, Dz, ubuff2, data->l, data->m,
-              data->nsoc, q);
+  nt_multiply(nt_scaling, nt_scaling_soc_idx, soc_idx, Dz, ubuff2, data->l,
+              data->m, data->nsoc, q);
 
   qoco_axpy(ubuff2, ubuff1, ubuff3, -1.0, data->m);
-  nt_multiply(Wfull, Wsoc_idx, soc_idx, ubuff3, Ds, data->l, data->m,
-              data->nsoc, q);
+  nt_multiply(nt_scaling, nt_scaling_soc_idx, soc_idx, ubuff3, Ds, data->l,
+              data->m, data->nsoc, q);
 
   // Compute step-size.
   QOCOFloat a =
@@ -480,8 +480,9 @@ void predictor_corrector(QOCOSolver* solver)
 }
 
 void kkt_multiply(QOCOFloat* x, QOCOFloat* y, QOCOProblemData* data,
-                  QOCOFloat* Wfull, QOCOInt* Wsoc_idx, QOCOInt* soc_idx,
-                  QOCOFloat* nbuff, QOCOFloat* mbuff1, QOCOFloat* mbuff2)
+                  QOCOFloat* nt_scaling, QOCOInt* nt_scaling_soc_idx,
+                  QOCOInt* soc_idx, QOCOFloat* nbuff, QOCOFloat* mbuff1,
+                  QOCOFloat* mbuff2)
 {
 
   // Compute y[1:n] = P * x[1:n] + A^T * x[n+1:n+p] + G^T * x[n+p+1:n+p+m].
@@ -506,11 +507,12 @@ void kkt_multiply(QOCOFloat* x, QOCOFloat* y, QOCOProblemData* data,
     SpMv(data->G, x, &y[data->n + data->p]);
   }
 
-  if (Wfull) {
-    nt_multiply(Wfull, Wsoc_idx, soc_idx, &x[data->n + data->p], mbuff1,
+  if (nt_scaling) {
+    nt_multiply(nt_scaling, nt_scaling_soc_idx, soc_idx, &x[data->n + data->p],
+                mbuff1, data->l, data->m, data->nsoc,
+                get_data_vectori(data->q));
+    nt_multiply(nt_scaling, nt_scaling_soc_idx, soc_idx, mbuff1, mbuff2,
                 data->l, data->m, data->nsoc, get_data_vectori(data->q));
-    nt_multiply(Wfull, Wsoc_idx, soc_idx, mbuff1, mbuff2, data->l, data->m,
-                data->nsoc, get_data_vectori(data->q));
     qoco_axpy(mbuff2, &y[data->n + data->p], &y[data->n + data->p], -1.0,
               data->m);
   }
diff --git a/src/qoco_api.c b/src/qoco_api.c
index e3c2fe0..d6e6ceb 100644
--- a/src/qoco_api.c
+++ b/src/qoco_api.c
@@ -129,21 +129,21 @@ QOCOInt qoco_setup(QOCOSolver* solver, QOCOInt n, QOCOInt m, QOCOInt p,
   QOCOInt Wnnz = m + Wsoc_nnz;
   work->Wnnz = Wnnz;
 
-  QOCOInt* Wsoc_idx = NULL;
+  QOCOInt* nt_scaling_soc_idx = NULL;
   QOCOInt* soc_idx = NULL;
   if (nsoc > 0) {
-    Wsoc_idx = (QOCOInt*)qoco_malloc(nsoc * sizeof(QOCOInt));
+    nt_scaling_soc_idx = (QOCOInt*)qoco_malloc(nsoc * sizeof(QOCOInt));
     soc_idx = (QOCOInt*)qoco_malloc(nsoc * sizeof(QOCOInt));
-    Wsoc_idx[0] = l;
+    nt_scaling_soc_idx[0] = l;
     soc_idx[0] = l;
     for (QOCOInt i = 1; i < nsoc; ++i) {
-      Wsoc_idx[i] = Wsoc_idx[i - 1] + q[i - 1] * q[i - 1];
+      nt_scaling_soc_idx[i] = nt_scaling_soc_idx[i - 1] + q[i - 1] + 1;
       soc_idx[i] = soc_idx[i - 1] + q[i - 1];
     }
   }
-  work->Wsoc_idx = new_qoco_vectori(Wsoc_idx, nsoc);
+  work->nt_scaling_soc_idx = new_qoco_vectori(nt_scaling_soc_idx, nsoc);
   work->soc_idx = new_qoco_vectori(soc_idx, nsoc);
-  qoco_free(Wsoc_idx);
+  qoco_free(nt_scaling_soc_idx);
   qoco_free(soc_idx);
 
   solver->linsys = &backend;
@@ -180,19 +180,17 @@ QOCOInt qoco_setup(QOCOSolver* solver, QOCOInt n, QOCOInt m, QOCOInt p,
   work->best_valid = 0;
 
   // Allocate Nesterov-Todd scalings and scaled variables.
-  QOCOInt Wnnzfull = data->l;
+  QOCOInt nt_scaling_nnz = data->l;
   set_cpu_mode(1);
   for (QOCOInt i = 0; i < data->nsoc; ++i) {
-    Wnnzfull +=
-        get_element_vectori(data->q, i) * get_element_vectori(data->q, i);
+    nt_scaling_nnz += get_element_vectori(data->q, i) + 1;
   }
   set_cpu_mode(0);
 
   work->W = new_qoco_vectorf(NULL, work->Wnnz);
-  work->Wfull = new_qoco_vectorf(NULL, Wnnzfull);
-  work->Wnnzfull = Wnnzfull;
+  work->nt_scaling = new_qoco_vectorf(NULL, nt_scaling_nnz);
+  work->nt_scaling_nnz = nt_scaling_nnz;
   work->Winv = new_qoco_vectorf(NULL, work->Wnnz);
-  work->Winvfull = new_qoco_vectorf(NULL, Wnnzfull);
   work->WtW = new_qoco_vectorf(NULL, work->Wnnz);
   work->lambda = new_qoco_vectorf(NULL, m);
 
@@ -591,11 +589,10 @@ QOCOInt qoco_cleanup(QOCOSolver* solver)
 
   // Free Nesterov-Todd scalings and scaled variables.
   free_qoco_vectorf(solver->work->W);
-  free_qoco_vectorf(solver->work->Wfull);
+  free_qoco_vectorf(solver->work->nt_scaling);
   free_qoco_vectorf(solver->work->Winv);
-  free_qoco_vectorf(solver->work->Winvfull);
   free_qoco_vectorf(solver->work->WtW);
-  free_qoco_vectori(solver->work->Wsoc_idx);
+  free_qoco_vectori(solver->work->nt_scaling_soc_idx);
   free_qoco_vectori(solver->work->soc_idx);
   free_qoco_vectorf(solver->work->lambda);
   free_qoco_vectorf(solver->work->sbar);