diff --git a/src/nns/__init__.py b/src/nns/__init__.py
index ba5b401..1d8f8ce 100644
--- a/src/nns/__init__.py
+++ b/src/nns/__init__.py
@@ -26,6 +26,7 @@
     "d_lpm": ("nns.co_moments", "d_lpm"),
     "dpm_nd": ("nns.dependence", "dpm_nd"),
     "dy_d": ("nns.diff", "dy_d"),
+    "dy_d_best": ("nns.diff", "dy_d_best"),
     "dy_dx": ("nns.diff", "dy_dx"),
     "d_upm": ("nns.co_moments", "d_upm"),
     "ecdf_pm": ("nns.classical", "ecdf_pm"),
diff --git a/src/nns/diff.py b/src/nns/diff.py
index fdf397f..781612f 100644
--- a/src/nns/diff.py
+++ b/src/nns/diff.py
@@ -262,6 +262,299 @@ def _combine_dy_d_outputs(
     }
 
 
+def dy_d_best(
+    x: NDArray[Any],
+    y: NDArray[Any],
+    wrt: int | NDArray[np.int64],
+    eval_points: str | float | NDArray[np.float64] = "obs",
+    *,
+    mixed: bool = False,
+    messages: bool = False,
+    factor_levels: Sequence[Sequence[Any] | None] | None = None,
+) -> dict[str, NDArray[np.float64]]:
+    """Reconciled ``dy.d_`` partial-derivative estimator.
+
+    This is the reconciliation of the original NNS 0.5.7 ``dy.d_`` (Vinod & Viole
+    2020, SSRN 3681436) with the current regression engine. It keeps the v0.5.7
+    finite-difference scaffolding but makes two changes that make the estimates
+    uniform across identically-distributed regressors and independent of the old
+    data.table machinery:
+
+    * **Step (``h_step``) shares the ``dy.dx`` logic** - a locally-adaptive step
+      centred on the evaluation point's percentile,
+      ``h_step = VaR(p + H, 1, x) - VaR(p - H, 1, x)`` with
+      ``p = LPM.ratio(1, eval, x)`` - instead of a global quantile spacing with a
+      cumulative window. This is what removes the cross-regressor scatter.
+    * **Estimates come from** ``nns_stack`` on the equal-weight synthetic regressor
+      ``X*`` via the increased-dimension trick ``cbind(X*, X*)``, with
+      ``method=(1, 2)``, ``dim_red_method="equal"``, ``order="max"`` and
+      ``folds=5``. The stack's cross-validated ``n.best`` regularises the (sharper)
+      current engine back toward the paper's regime.
+
+    Bandwidths follow v0.5.7: ``h_s = 1/log(size(x), [2, 10])`` extended by
+    ``10 * h_s`` and doubled when ``nns_dep(x[:, wrt], y) < 0.5``. First and second
+    derivatives use the central-difference forms
+    ``First = (upper - lower) / (2 h_step)`` and
+    ``Second = (upper - 2 f(x) + lower) / h_step ** 2`` (matching ``dy.dx``),
+    blended across bandwidths with a plain ``nanmean``.
+
+    ``wrt`` uses R-style 1-based indexing into the factor-expanded predictor
+    matrix. A scalar/1-D ``eval_points`` evaluates only the selected regressor
+    (averaged over the distribution of the other regressors); a 2-D array
+    evaluates complete predictor tuples.
+    """
+    raw_x = np.asarray(x)
+    if raw_x.ndim != 2:
+        raise ValueError("Please ensure (x) is a matrix or data.frame type object.")
+    if raw_x.shape[1] < 2:
+        raise ValueError("Please use NNS::dy.dx(...) for univariate partial derivatives.")
+
+    raw_y = np.asarray(y).reshape(-1)
+    if raw_y.size != raw_x.shape[0]:
+        raise ValueError("x and y must have compatible row counts.")
+    if _dy_d_has_missing(raw_x) or _dy_d_has_missing(raw_y):
+        raise ValueError("You have some missing values, please address.")
+
+    x_values = _dy_d_expand_predictors(x, factor_levels=factor_levels)
+    y_values = np.asarray(raw_y, dtype=np.float64)
+    if _dy_d_has_missing(x_values) or _dy_d_has_missing(y_values):
+        raise ValueError("You have some missing values, please address.")
+
+    wrt_values = _dy_d_validate_wrt(wrt, x_values.shape[1])
+    outputs = [
+        _dy_d_best_scalar(
+            x_values,
+            y_values,
+            int(wrt_value) - 1,
+            eval_points,
+            mixed=bool(mixed),
+            messages=bool(messages),
+            wrt_label=int(wrt_value),
+        )
+        for wrt_value in wrt_values
+    ]
+    if len(outputs) == 1:
+        return outputs[0]
+    return _combine_dy_d_outputs(outputs)
+
+
+def _dy_d_best_reg_estimates(
+    x: NDArray[np.float64],
+    y: NDArray[np.float64],
+    test_points: NDArray[np.float64],
+) -> NDArray[np.float64]:
+    """f(x +/- h) via NNS.stack on the equal-weight synthetic regressor X*.
+
+    Reduces both the training design and the test points to the equal-weight
+    synthetic regressor X* = rowMeans(.), then estimates with the NNS
+    increased-dimension trick cbind(X*, X*) under
+    ``method=(1, 2), dim_red_method="equal", order="max", folds=5``. The
+    cross-validated ``n.best`` is what regularises the current engine.
+    """
+    from nns.stack import nns_stack
+
+    x_star = np.asarray(x, dtype=np.float64).mean(axis=1)
+    test_star = np.asarray(test_points, dtype=np.float64).mean(axis=1)
+    result = nns_stack(
+        ivs_train=np.column_stack([x_star, x_star]),
+        dv_train=y,
+        ivs_test=np.column_stack([test_star, test_star]),
+        method=(1, 2),
+        dim_red_method="equal",
+        status=False,
+        order="max",
+        folds=5,
+        ncores=1,
+        dist=None,
+    )
+    return np.asarray(result["stack"], dtype=np.float64).reshape(-1)
+
+
+def _r_seq(start: float, stop: float, by: float) -> NDArray[np.float64]:
+    """Reproduce R's seq(start, stop, by)."""
+    if not np.isfinite(by) or by == 0.0:
+        return np.asarray([start], dtype=np.float64)
+    count = int(np.floor((stop - start) / by + 1e-9))
+    return start + by * np.arange(0, count + 1, dtype=np.float64)
+
+
+def _v057_bandwidths(x_size: int, wrt_dependence: float) -> NDArray[np.float64]:
+    """h_s = 1/log(length(x), c(2, 10)); c(h_s, 10*h_s); doubled if dependence < .5."""
+    base = np.log(np.asarray([2.0, 10.0])) / np.log(float(x_size))
+    h_s = np.concatenate([base, 10.0 * base])
+    if wrt_dependence < 0.5:
+        h_s = 2.0 * h_s
+    return h_s
+
+
+def _dydx_step(column: NDArray[np.float64], eval_value: float, band: float) -> float:
+    """Locally-adaptive dy.dx step: VaR(p + H, 1, x) - VaR(p - H, 1, x), p = CDF(eval)."""
+    from nns.core import lpm_ratio
+    from nns.var import lpm_var
+
+    p = float(lpm_ratio(1.0, float(eval_value), column))
+    upper = lpm_var(min(1.0, p + band), 1.0, column)
+    lower = lpm_var(max(0.0, p - band), 1.0, column)
+    return float(upper - lower)
+
+
+def _dy_d_best_scalar(
+    x_values: NDArray[np.float64],
+    y_values: NDArray[np.float64],
+    wrt_index: int,
+    eval_points: str | float | NDArray[np.float64],
+    *,
+    mixed: bool,
+    messages: bool,
+    wrt_label: int,
+) -> dict[str, NDArray[np.float64]]:
+    from nns.dependence import nns_dep
+    from nns.var import lpm_var
+
+    _n_rows, n_predictors = x_values.shape
+    if wrt_index < 0 or wrt_index >= n_predictors:
+        raise ValueError("`wrt` must select exactly one column of the expanded predictor matrix.")
+    if n_predictors != 2:
+        mixed = False
+
+    if messages:
+        print(
+            "Currently generating NNS.reg finite difference estimates...Regressor "
+            f"{wrt_label}\r"
+        )
+
+    eval_values, vector_branch = _dy_d_eval_points(x_values, wrt_index, eval_points)
+
+    column = x_values[:, wrt_index]
+    dependence = float(nns_dep(column, y_values)["Dependence"])
+    h_s = _v057_bandwidths(x_values.size, dependence)
+
+    firsts: list[NDArray[np.float64]] = []
+    seconds: list[NDArray[np.float64]] = []
+    mixeds: list[NDArray[np.float64]] = []
+
+    if vector_branch:
+        eval_vec = np.asarray(eval_values, dtype=np.float64).reshape(-1)
+        grid = np.column_stack(
+            [
+                np.asarray(
+                    [lpm_var(float(p), 0.0, x_values[:, col]) for p in _r_seq(0.0, 1.0, 0.05)],
+                    dtype=np.float64,
+                )
+                for col in range(n_predictors)
+            ]
+        )
+        sample_size = grid.shape[0]
+        k = eval_vec.size
+        for band in h_s:
+            steps = np.array([_dydx_step(column, ev, float(band)) for ev in eval_vec])
+            blocks: list[NDArray[np.float64]] = []
+            position: list[str] = []
+            ids: list[int] = []
+            for g in range(k):
+                lower = grid.copy()
+                middle = grid.copy()
+                upper = grid.copy()
+                lower[:, wrt_index] = eval_vec[g] - steps[g]
+                middle[:, wrt_index] = eval_vec[g]
+                upper[:, wrt_index] = eval_vec[g] + steps[g]
+                blocks.extend((lower, middle, upper))
+                position.extend(["l"] * sample_size + ["m"] * sample_size + ["u"] * sample_size)
+                ids.extend([g] * (3 * sample_size))
+            estimates = _dy_d_best_reg_estimates(x_values, y_values, np.vstack(blocks))
+            pos = np.asarray(position, dtype=object)
+            idx = np.asarray(ids)
+            band_first = np.empty(k)
+            band_second = np.empty(k)
+            for g in range(k):
+                lo = np.mean(estimates[(pos == "l") & (idx == g)])
+                mid = np.mean(estimates[(pos == "m") & (idx == g)])
+                up = np.mean(estimates[(pos == "u") & (idx == g)])
+                h = steps[g]
+                if np.isfinite(h) and h != 0.0:
+                    band_first[g] = (up - lo) / (2.0 * h)
+                    band_second[g] = (up - 2.0 * mid + lo) / (h**2)
+                else:
+                    band_first[g] = np.nan
+                    band_second[g] = np.nan
+            firsts.append(band_first)
+            seconds.append(band_second)
+        mixed_eval: NDArray[np.float64] | None = None
+    else:
+        eval_mat = _as_eval_matrix(eval_values, n_predictors)
+        n_eval = eval_mat.shape[0]
+        for band in h_s:
+            steps = np.array(
+                [_dydx_step(column, eval_mat[i, wrt_index], float(band)) for i in range(n_eval)]
+            )
+            finite = np.isfinite(steps) & (steps != 0.0)
+            lower = eval_mat.copy()
+            upper = eval_mat.copy()
+            lower[:, wrt_index] = eval_mat[:, wrt_index] - steps
+            upper[:, wrt_index] = eval_mat[:, wrt_index] + steps
+            if messages:
+                print(
+                    "Currently generating NNS.reg finite difference estimates...bandwidth\r"
+                )
+            estimates = _dy_d_best_reg_estimates(
+                x_values, y_values, np.vstack((lower, eval_mat, upper))
+            )
+            lo = estimates[:n_eval]
+            mid = estimates[n_eval : 2 * n_eval]
+            up = estimates[2 * n_eval :]
+            with np.errstate(invalid="ignore", divide="ignore"):
+                first = (up - lo) / (2.0 * steps)
+                second = (up - 2.0 * mid + lo) / (steps**2)
+            first[~finite] = np.nan
+            second[~finite] = np.nan
+            firsts.append(first)
+            seconds.append(second)
+        mixed_eval = eval_mat
+
+    def _row_nanmean(bands: list[NDArray[np.float64]]) -> NDArray[np.float64]:
+        matrix = np.column_stack([np.asarray(b, dtype=np.float64).reshape(-1) for b in bands])
+        with np.errstate(invalid="ignore"):
+            return np.asarray(np.nanmean(matrix, axis=1), dtype=np.float64)
+
+    output = {"First": _row_nanmean(firsts), "Second": _row_nanmean(seconds)}
+
+    if mixed:
+        if vector_branch:
+            tuple_eval = np.asarray(eval_values, dtype=np.float64).reshape(-1)
+            if tuple_eval.size != 2:
+                raise ValueError("Mixed Derivatives are only for 2 IV")
+            mixed_points = tuple_eval.reshape(1, 2)
+        else:
+            assert mixed_eval is not None
+            if mixed_eval.shape[1] != 2:
+                raise ValueError("Mixed Derivatives are only for 2 IV")
+            mixed_points = mixed_eval
+        for band in h_s:
+            band_vals: list[float] = []
+            for point in mixed_points:
+                s1 = _dydx_step(x_values[:, 0], point[0], float(band))
+                s2 = _dydx_step(x_values[:, 1], point[1], float(band))
+                if not (np.isfinite(s1) and np.isfinite(s2) and s1 != 0.0 and s2 != 0.0):
+                    band_vals.append(np.nan)
+                    continue
+                corners = np.array(
+                    [
+                        [point[0] + s1, point[1] + s2],
+                        [point[0] - s1, point[1] + s2],
+                        [point[0] + s1, point[1] - s2],
+                        [point[0] - s1, point[1] - s2],
+                    ]
+                )
+                z = _dy_d_best_reg_estimates(x_values, y_values, corners)
+                band_vals.append((z[0] + z[3] - z[1] - z[2]) / (4.0 * s1 * s2))
+            mixeds.append(np.asarray(band_vals, dtype=np.float64))
+        output["Mixed"] = _row_nanmean(mixeds)
+
+    if messages:
+        print("\r")
+    return output
+
+
 def _dy_d_scalar(
     x_values: NDArray[np.float64],
     y_values: NDArray[np.float64],
diff --git a/tests/invariants/test_diff.py b/tests/invariants/test_diff.py
index 1c74520..81a589f 100644
--- a/tests/invariants/test_diff.py
+++ b/tests/invariants/test_diff.py
@@ -115,3 +115,80 @@ def test_dy_d_point_modes_preserve_linear_slope_direction() -> None:
     for eval_points in ("mean", "median", "last"):
         assert dy_d(x, positive, wrt=1, eval_points=eval_points)["First"][0] > 0.0
         assert dy_d(x, negative, wrt=1, eval_points=eval_points)["First"][0] < 0.0
+
+
+def test_dy_d_best_is_exported() -> None:
+    import nns
+
+    assert "dy_d_best" in nns.__all__
+    assert callable(nns.dy_d_best)
+
+
+def test_dy_d_best_vectorized_wrt_tuple_returns_first_second() -> None:
+    from nns import dy_d_best
+
+    x = np.random.RandomState(0).randn(60, 3)
+    y = x[:, 0] + 2.0 * x[:, 1] - x[:, 2]
+
+    result = dy_d_best(x, y, wrt=np.array([1, 2, 3]), eval_points=np.zeros((1, 3)))
+
+    assert result.keys() == {"First", "Second"}
+    assert result["First"].shape == (1, 3)
+    assert result["Second"].shape == (1, 3)
+    assert np.all(np.isfinite(result["First"]))
+
+
+def test_dy_d_best_single_wrt_returns_one_dimensional_first() -> None:
+    from nns import dy_d_best
+
+    x = np.random.RandomState(0).randn(60, 3)
+    y = x[:, 0] + 2.0 * x[:, 1] - x[:, 2]
+
+    # A single `wrt` returns 1-D arrays (one value per evaluation point).
+    result = dy_d_best(x, y, wrt=1, eval_points=np.array([[0.0, 0.0, 0.0]]))
+
+    assert result["First"].shape == (1,)
+    assert np.all(np.isfinite(result["First"]))
+
+
+def test_dy_d_best_two_column_mixed_returns_mixed() -> None:
+    from nns import dy_d_best
+
+    x = np.random.RandomState(1).randn(60, 2)
+    y = x[:, 0] * x[:, 1]
+
+    result = dy_d_best(x, y, wrt=np.array([1, 2]), eval_points=np.array([[0.0, 0.0]]), mixed=True)
+
+    assert result.keys() == {"First", "Second", "Mixed"}
+    assert result["First"].shape == (1, 2)
+    assert result["Mixed"].shape == (1, 2)
+
+
+def test_dy_d_best_is_deterministic() -> None:
+    from nns import dy_d_best
+
+    x = np.random.RandomState(0).randn(60, 3)
+    y = x[:, 0] + 2.0 * x[:, 1] - x[:, 2]
+
+    a = dy_d_best(x, y, wrt=np.array([1, 2, 3]), eval_points=np.zeros((1, 3)))
+    b = dy_d_best(x, y, wrt=np.array([1, 2, 3]), eval_points=np.zeros((1, 3)))
+
+    assert np.allclose(a["First"], b["First"])
+    assert np.allclose(a["Second"], b["Second"])
+
+
+def test_dy_d_best_matches_pinned_values() -> None:
+    # Golden values guard the v0.5.7 finite-difference design against regressions.
+    from nns import dy_d_best
+
+    x = np.random.RandomState(0).randn(60, 3)
+    y = x[:, 0] + 2.0 * x[:, 1] - x[:, 2]
+
+    result = dy_d_best(x, y, wrt=np.array([1, 2, 3]), eval_points=np.zeros((1, 3)))
+
+    np.testing.assert_allclose(
+        np.ravel(result["First"]),
+        np.array([0.754963, 0.792164, 0.699243]),
+        rtol=0.0,
+        atol=1e-4,
+    )