diff --git a/packages/tabarena/src/tabarena/simulation/ensemble/__init__.py b/packages/tabarena/src/tabarena/simulation/ensemble/__init__.py
index 0b468281e..14c7666bd 100644
--- a/packages/tabarena/src/tabarena/simulation/ensemble/__init__.py
+++ b/packages/tabarena/src/tabarena/simulation/ensemble/__init__.py
@@ -15,6 +15,7 @@
     TopKAverageEnsembler,
 )
 from tabarena.simulation.ensemble.greedy_ensembler import GreedyEnsembler
+from tabarena.simulation.ensemble.hill_climbing_ensembler import HillClimbingEnsembler
 from tabarena.simulation.ensemble.stacking_ensembler import StackingEnsembler
 
 __all__ = [
@@ -23,6 +24,7 @@
     "AutoGluonStackerRegressor",
     "FixedWeightsEnsembler",
     "GreedyEnsembler",
+    "HillClimbingEnsembler",
     "LegacyEnsemblerAdapter",
     "SingleBestEnsembler",
     "StackingEnsembler",
diff --git a/packages/tabarena/src/tabarena/simulation/ensemble/hill_climbing_ensembler.py b/packages/tabarena/src/tabarena/simulation/ensemble/hill_climbing_ensembler.py
new file mode 100644
index 000000000..bb6ba30f9
--- /dev/null
+++ b/packages/tabarena/src/tabarena/simulation/ensemble/hill_climbing_ensembler.py
@@ -0,0 +1,242 @@
+"""Kaggle-style hill-climbing ensemble (Matt-OP / Deotte family).
+
+Same approach family as :class:`GreedyEnsembler` (Caruana GES / AG ``EnsembleSelection``):
+greedy local search over linear combinations of OOF predictions. Differences are the step
+parameterization (continuous ``w`` grid vs discrete multiset votes), not a different problem.
+
+Each step evaluates all ``(model, w)`` candidates and accepts the single best improvement.
+Optionally includes the Caruana step weight ``w = 1/(n_support+1)`` on the grid.
+
+This class is optional research / post-hoc; TabArena's default remains
+:class:`GreedyEnsembler`.
+"""
+
+from __future__ import annotations
+
+import time
+from typing import TYPE_CHECKING
+
+import numpy as np
+
+from tabarena.simulation.ensemble.abstract_ensembler import WeightedEnsembler
+
+if TYPE_CHECKING:
+    from autogluon.core.metrics import Scorer
+
+
+class HillClimbingEnsembler(WeightedEnsembler):
+    """Iterative convex blending of base-model predictions (hill climbing).
+
+    Starts from the best single model (or optional warm-start weights), then repeatedly
+    blends ``(1 - w) * ensemble + w * model`` and keeps the best improving step.
+
+    Parameters
+    ----------
+    precision : float, default 0.01
+        Weight step on ``(0, 1]``. Smaller is slower and can overfit small OOF sets.
+    max_rounds : int, default 100
+        Maximum accepted improvement steps. Stops early when a step finds no improvement.
+    include_caruana_step : bool, default True
+        Also try ``w = 1/(n_support + 1)`` each step (GES-style blend weight when treating
+        the current ensemble as one unit).
+    allow_negative_weights : bool, default False
+        If True, also tries ``w`` on ``[-1, 0)``. Safer off on small OOF.
+    max_models : int | None, default None
+        Optional cap on how many models may receive non-zero weight. ``None`` means no cap.
+    initial_weights : array-like | None, default None
+        Optional warm-start weights (one entry per model). When set, skip best-single
+        initialization. Must be finite; non-negative unless ``allow_negative_weights``;
+        renormalized to sum to 1.
+    random_state : int | np.random.RandomState | None, default None
+        Tie-breaking when selecting the initial best single model.
+    """
+
+    def __init__(
+        self,
+        *,
+        problem_type: str,
+        metric: Scorer,
+        precision: float = 0.01,
+        max_rounds: int = 100,
+        include_caruana_step: bool = True,
+        allow_negative_weights: bool = False,
+        max_models: int | None = None,
+        initial_weights: np.ndarray | list[float] | None = None,
+        random_state: int | np.random.RandomState | None = None,
+    ):
+        super().__init__(problem_type=problem_type, metric=metric)
+        if precision <= 0 or precision > 1:
+            raise ValueError(f"precision must be in (0, 1], got {precision}")
+        if max_rounds < 1:
+            raise ValueError(f"max_rounds must be >= 1, got {max_rounds}")
+        if max_models is not None and max_models < 1:
+            raise ValueError(f"max_models must be >= 1 or None, got {max_models}")
+        self.precision = float(precision)
+        self.max_rounds = int(max_rounds)
+        self.include_caruana_step = bool(include_caruana_step)
+        self.allow_negative_weights = bool(allow_negative_weights)
+        self.max_models = max_models
+        self.initial_weights = None if initial_weights is None else np.asarray(initial_weights, dtype=np.float64)
+        if isinstance(random_state, np.random.RandomState):
+            self.random_state = random_state
+        else:
+            self.random_state = np.random.RandomState(0 if random_state is None else int(random_state))
+
+        self.trajectory_: list[float] = []
+        self.n_rounds_: int = 0
+        self.init_val_error_: float | None = None
+
+    def _weight_grid(self, *, n_support: int | None = None) -> np.ndarray:
+        pos = np.arange(self.precision, 1.0 + self.precision * 0.5, self.precision)
+        pos = np.clip(pos, self.precision, 1.0)
+        pos = np.unique(np.round(pos / self.precision) * self.precision)
+        if self.include_caruana_step and n_support is not None and n_support >= 1:
+            w_c = 1.0 / float(n_support + 1)
+            if 0 < w_c <= 1:
+                pos = np.unique(np.concatenate([pos, np.asarray([w_c], dtype=np.float64)]))
+        if self.allow_negative_weights:
+            return np.concatenate([-pos[::-1], pos])
+        return pos
+
+    def _allowed_model(self, j: int, w: float, weights: np.ndarray) -> bool:
+        if self.max_models is None or weights[j] != 0 or w == 0:
+            return True
+        return self._n_nonzero(weights) < self.max_models
+
+    def _try_blend(
+        self,
+        *,
+        ensemble_pred: np.ndarray,
+        pred_j: np.ndarray,
+        w: float,
+        labels: np.ndarray,
+        best_error: float,
+    ) -> tuple[float, float, np.ndarray] | None:
+        trial = (1.0 - w) * ensemble_pred + w * pred_j
+        trial = self._renormalize_proba(trial)
+        err = self._score_error(labels, trial)
+        if err < best_error - 1e-15:
+            return err, float(w), trial
+        return None
+
+    def _combine(self, predictions: list[np.ndarray], weights: np.ndarray) -> np.ndarray:
+        preds_norm = [pred * w for pred, w in zip(predictions, weights, strict=True) if w != 0]
+        if not preds_norm:
+            return np.mean(predictions, axis=0)
+        return np.sum(preds_norm, axis=0)
+
+    def _n_nonzero(self, weights: np.ndarray) -> int:
+        return int(np.sum(np.abs(weights) > 0))
+
+    def _renormalize_proba(self, trial: np.ndarray) -> np.ndarray:
+        """Keep multiclass / softclass rows on the probability simplex."""
+        if trial.ndim == 2 and self.problem_type in ("multiclass", "softclass"):
+            row_sum = trial.sum(axis=1, keepdims=True)
+            row_sum = np.where(row_sum == 0, 1.0, row_sum)
+            trial = trial / row_sum
+            np.maximum(trial, 0.0, out=trial)
+            row_sum = trial.sum(axis=1, keepdims=True)
+            row_sum = np.where(row_sum == 0, 1.0, row_sum)
+            trial = trial / row_sum
+        return trial
+
+    def _fit(self, *, predictions: np.ndarray, labels: np.ndarray, time_limit: float | None = None) -> None:
+        start = time.time()
+        # float64 for stable blends; renorm multiclass so float32 OOF stays on the simplex.
+        predictions = [self._renormalize_proba(np.asarray(p, dtype=np.float64)) for p in predictions]
+        n_models = len(predictions)
+        if n_models == 0:
+            raise ValueError("HillClimbingEnsembler requires at least one model")
+
+        if self.initial_weights is not None:
+            weights = np.asarray(self.initial_weights, dtype=np.float64).copy()
+            if weights.shape != (n_models,):
+                raise ValueError(f"initial_weights length {weights.shape} != n_models={n_models}")
+            if not np.isfinite(weights).all():
+                raise ValueError("initial_weights must be finite")
+            if not self.allow_negative_weights:
+                weights = np.maximum(weights, 0.0)
+            total0 = float(np.sum(weights))
+            if abs(total0) <= 1e-12:
+                raise ValueError("initial_weights sum to ~0; cannot warm-start")
+            weights = weights / total0
+            ensemble_pred = self._renormalize_proba(self._combine(predictions, weights))
+            best_error = float(self._score_error(labels, ensemble_pred))
+            start_idx = int(np.argmax(np.abs(weights)))
+        else:
+            single_errors = np.array([self._score_error(labels, pred) for pred in predictions], dtype=np.float64)
+            best_error = float(np.nanmin(single_errors))
+            candidates = np.flatnonzero(np.isclose(single_errors, best_error, atol=0, rtol=1e-12))
+            start_idx = int(self.random_state.choice(candidates))
+            weights = np.zeros(n_models, dtype=np.float64)
+            weights[start_idx] = 1.0
+            ensemble_pred = predictions[start_idx].copy()
+
+        self.init_val_error_ = best_error
+        self.trajectory_ = [best_error]
+
+        for round_i in range(self.max_rounds):
+            if time_limit is not None and (time.time() - start) >= time_limit:
+                break
+
+            n_support = max(self._n_nonzero(weights), 1)
+            weight_grid = self._weight_grid(n_support=n_support)
+            best_local_error = best_error
+            best_j = None
+            best_local_w = None
+            best_local_pred = None
+            for j in range(n_models):
+                if time_limit is not None and (time.time() - start) >= time_limit:
+                    break
+                pred_j = predictions[j]
+                for w in weight_grid:
+                    if not self._allowed_model(j, float(w), weights):
+                        continue
+                    hit = self._try_blend(
+                        ensemble_pred=ensemble_pred,
+                        pred_j=pred_j,
+                        w=float(w),
+                        labels=labels,
+                        best_error=best_local_error,
+                    )
+                    if hit is not None:
+                        best_local_error, best_local_w, best_local_pred = hit
+                        best_j = j
+
+            self.n_rounds_ = round_i + 1
+            if best_j is None or best_local_w is None:
+                break
+
+            weights *= 1.0 - best_local_w
+            weights[best_j] += best_local_w
+            ensemble_pred = best_local_pred
+            best_error = best_local_error
+            self.trajectory_.append(best_error)
+
+        weights[np.abs(weights) < 1e-12] = 0.0
+        if not self.allow_negative_weights:
+            weights = np.maximum(weights, 0.0)
+        total = float(np.sum(weights))
+        if abs(total) > 1e-12:
+            weights = weights / total
+        else:
+            weights = np.zeros(n_models, dtype=np.float64)
+            weights[start_idx] = 1.0
+
+        self.weights_ = weights
+
+    def predict_proba(self, predictions: np.ndarray) -> np.ndarray:
+        out = super().predict_proba(predictions)
+        return self._renormalize_proba(np.asarray(out, dtype=np.float64))
+
+    def info(self) -> dict:
+        return {
+            "n_rounds": self.n_rounds_,
+            "trajectory_len": len(self.trajectory_),
+            "init_val_error": self.init_val_error_,
+            "final_val_error": self.trajectory_[-1] if self.trajectory_ else None,
+            "warm_started": self.initial_weights is not None,
+            "include_caruana_step": self.include_caruana_step,
+            "precision": self.precision,
+            "allow_negative_weights": self.allow_negative_weights,
+        }
diff --git a/packages/tabarena/src/tabarena/utils/rank_utils.py b/packages/tabarena/src/tabarena/utils/rank_utils.py
index e2216ac86..21cf1ed32 100644
--- a/packages/tabarena/src/tabarena/utils/rank_utils.py
+++ b/packages/tabarena/src/tabarena/utils/rank_utils.py
@@ -91,10 +91,17 @@ def __init__(
         self.pct = pct
         self.include_partial = include_partial
         df_pivot = df_results.pivot_table(values=metric_error_col, index=task_col, columns=framework_col)
-        df_pivot.values.sort(axis=1)  # NOTE: The framework columns are now no longer correct. Do not use them.
-
-        # tolist to drop the framework col name, since it is no longer ordered.
-        self.error_dict = {dataset: df_pivot.loc[dataset].dropna().tolist() for dataset in tasks}
+        # Copy: recent pandas/numpy expose read-only views for .values, so in-place sort fails.
+        sorted_errors = np.ascontiguousarray(df_pivot.to_numpy(dtype=np.float64, copy=True))
+        sorted_errors.flags.writeable = True
+        sorted_errors.sort(axis=1)
+        # NOTE: Framework columns are no longer meaningful after row-wise sort.
+        index_list = list(df_pivot.index)
+        self.error_dict = {}
+        for task in tasks:
+            i = index_list.index(task)
+            row = sorted_errors[i]
+            self.error_dict[task] = row[~np.isnan(row)].tolist()
 
     def rank(self, task: str, error: float) -> float:
         """Get the rank of a result on a dataset given an error."""
diff --git a/tests/tabarena/simulation/test_ensembler.py b/tests/tabarena/simulation/test_ensembler.py
index 09e03e54f..ab6ca2f10 100644
--- a/tests/tabarena/simulation/test_ensembler.py
+++ b/tests/tabarena/simulation/test_ensembler.py
@@ -545,3 +545,200 @@ def test_autogluon_stacker_leaves_no_cwd_artifacts(monkeypatch, tmp_path):
 
     assert np.isfinite(regressor.predict(preds.T)).all()
     assert list(tmp_path.iterdir()) == []
+
+
+# -------------------------
+# Hill climbing (autogluon/autogluon#4505)
+# -------------------------
+def test_hill_climbing_beats_or_matches_single_best_on_synthetic():
+    """HC starts from the best single model and only accepts improving blends, so
+    validation error must be <= single-best error on the fit split.
+    """
+    from tabarena.simulation.ensemble import HillClimbingEnsembler, SingleBestEnsembler
+
+    y, preds = _make_binary_task(n_models=12, n_samples=800, seed=7)
+    metric = get_metric(metric="roc_auc", problem_type="binary")
+
+    single = SingleBestEnsembler(problem_type="binary", metric=metric)
+    single.fit(predictions=preds, labels=y)
+    single_err = metric.error(y, single.predict_proba(preds))
+
+    hc = HillClimbingEnsembler(
+        problem_type="binary",
+        metric=metric,
+        precision=0.05,
+        max_rounds=20,
+        random_state=0,
+    )
+    assert hc.include_caruana_step is True
+    hc.fit(predictions=preds, labels=y)
+    hc_err = metric.error(y, hc.predict_proba(preds))
+
+    assert hc_err <= single_err + 1e-12
+    assert hc.model_weights() is not None
+    assert np.isclose(hc.model_weights().sum(), 1.0)
+    assert hc.models_used().any()
+
+
+def test_hill_climbing_regression_and_weights_sum():
+    from tabarena.simulation.ensemble import HillClimbingEnsembler
+
+    y, preds = _make_regression_task(n_models=10, n_samples=600, seed=3)
+    metric = get_metric(metric="rmse", problem_type="regression")
+    hc = HillClimbingEnsembler(
+        problem_type="regression",
+        metric=metric,
+        precision=0.05,
+        max_rounds=15,
+        random_state=1,
+    )
+    hc.fit(predictions=preds, labels=y)
+    w = hc.model_weights()
+    assert w is not None
+    assert np.isclose(w.sum(), 1.0)
+    assert (w >= -1e-12).all()
+    # Prediction is a weighted sum of base preds
+    combined = hc.predict_proba(preds)
+    expected = sum(p * wi for p, wi in zip(preds, w, strict=True) if wi != 0)
+    np.testing.assert_allclose(combined, expected, rtol=1e-5, atol=1e-5)
+
+
+def test_hill_climbing_task_evaluator_runs():
+    """Smoke: HillClimbingEnsembler plugs into TaskEvaluator like other ensemblers."""
+    from tabarena.simulation.ensemble import HillClimbingEnsembler
+
+    y, preds = _make_binary_task(seed=11)
+    metric = get_metric(metric="roc_auc", problem_type="binary")
+    results, ensemble = _run_task_evaluator(
+        HillClimbingEnsembler,
+        {"precision": 0.05, "max_rounds": 10, "random_state": 0},
+        problem_type="binary",
+        eval_metric=metric,
+        fit_eval_metric=metric,
+        y=y,
+        preds=preds,
+    )
+    assert "metric_error" in results
+    assert "ensemble_weights" in results
+    assert results["ensemble_weights"] is not None
+    assert ensemble is not None
+
+
+def test_hill_climbing_invariants_and_edge_cases():
+    """Trajectory monotonicity, sparsity cap, single-model, determinism, hyperparam guards."""
+    from tabarena.simulation.ensemble import HillClimbingEnsembler
+
+    y, preds = _make_binary_task(n_models=12, n_samples=600, seed=13)
+    metric = get_metric(metric="roc_auc", problem_type="binary")
+
+    hc = HillClimbingEnsembler(problem_type="binary", metric=metric, precision=0.05, max_rounds=25, random_state=0)
+    hc.fit(predictions=preds, labels=y)
+    traj = np.asarray(hc.trajectory_)
+    assert traj.ndim == 1 and len(traj) >= 1
+    assert np.all(np.diff(traj) <= 1e-15)
+    single_best_err = min(metric.error(y, p) for p in preds)
+    assert np.isclose(traj[0], single_best_err, rtol=1e-10, atol=1e-12)
+    assert traj[-1] <= single_best_err + 1e-12
+    assert np.isclose(hc.model_weights().sum(), 1.0)
+    assert (hc.model_weights() >= -1e-12).all()
+    assert hc.info()["n_rounds"] >= 1
+
+    # max_models sparsity
+    hc_sparse = HillClimbingEnsembler(
+        problem_type="binary",
+        metric=metric,
+        precision=0.05,
+        max_rounds=25,
+        max_models=3,
+        random_state=0,
+    )
+    hc_sparse.fit(predictions=preds, labels=y)
+    assert int(hc_sparse.models_used().sum()) <= 3
+    assert metric.error(y, hc_sparse.predict_proba(preds)) <= single_best_err + 1e-12
+
+    # single model pool
+    hc_one = HillClimbingEnsembler(problem_type="binary", metric=metric, precision=0.1, max_rounds=5, random_state=0)
+    hc_one.fit(predictions=preds[:1], labels=y)
+    np.testing.assert_allclose(hc_one.model_weights(), [1.0])
+
+    # determinism
+    a = HillClimbingEnsembler(problem_type="binary", metric=metric, precision=0.05, max_rounds=15, random_state=42)
+    b = HillClimbingEnsembler(problem_type="binary", metric=metric, precision=0.05, max_rounds=15, random_state=42)
+    a.fit(predictions=preds, labels=y)
+    b.fit(predictions=preds, labels=y)
+    np.testing.assert_allclose(a.model_weights(), b.model_weights())
+
+    # does not mutate caller arrays
+    preds_copy = preds.copy()
+    HillClimbingEnsembler(problem_type="binary", metric=metric, precision=0.1, max_rounds=5, random_state=0).fit(
+        predictions=preds, labels=y
+    )
+    assert np.array_equal(preds, preds_copy)
+
+    # hyperparam validation
+    with pytest.raises(ValueError, match="precision"):
+        HillClimbingEnsembler(problem_type="binary", metric=metric, precision=0)
+    with pytest.raises(ValueError, match="max_rounds"):
+        HillClimbingEnsembler(problem_type="binary", metric=metric, max_rounds=0)
+    with pytest.raises(ValueError, match="max_models"):
+        HillClimbingEnsembler(problem_type="binary", metric=metric, max_models=0)
+    with pytest.raises(ValueError, match="at least one"):
+        HillClimbingEnsembler(problem_type="binary", metric=metric).fit(predictions=np.zeros((0, 10)), labels=y[:10])
+
+
+def test_hill_climbing_multiclass_simplex():
+    """Multiclass blends stay on the probability simplex; fit error ≤ single-best."""
+    from tabarena.simulation.ensemble import HillClimbingEnsembler
+
+    y, preds = _make_multiclass_task(n_models=8, n_samples=500, n_classes=4, seed=14)
+    metric = get_metric(metric="log_loss", problem_type="multiclass")
+    hc = HillClimbingEnsembler(problem_type="multiclass", metric=metric, precision=0.05, max_rounds=20, random_state=0)
+    hc.fit(predictions=preds, labels=y)
+    out = hc.predict_proba(preds)
+    assert out.shape == (len(y), preds.shape[2])
+    np.testing.assert_allclose(out.sum(axis=1), 1.0, atol=1e-5)
+    assert (out >= -1e-8).all()
+    single_best_err = min(metric.error(y, p) for p in preds)
+    assert metric.error(y, out) <= single_best_err + 1e-10
+
+
+def test_hill_climbing_time_limit_returns_valid_weights():
+    """time_limit may stop early but must still yield a valid weight vector."""
+    from tabarena.simulation.ensemble import HillClimbingEnsembler
+
+    y, preds = _make_binary_task(n_models=30, n_samples=2000, seed=15)
+    metric = get_metric(metric="roc_auc", problem_type="binary")
+    hc = HillClimbingEnsembler(problem_type="binary", metric=metric, precision=0.01, max_rounds=200, random_state=0)
+    hc.fit(predictions=preds, labels=y, time_limit=0.02)
+    w = hc.model_weights()
+    assert w is not None
+    assert np.isclose(w.sum(), 1.0)
+    assert np.isfinite(w).all()
+
+
+def test_hill_climbing_warm_start_from_greedy_weights():
+    """Warm-start HC from Greedy weights: refined val error must be <= Greedy val error."""
+    from tabarena.simulation.ensemble import GreedyEnsembler, HillClimbingEnsembler
+
+    y, preds = _make_binary_task(n_models=12, n_samples=800, seed=21)
+    metric = get_metric(metric="roc_auc", problem_type="binary")
+
+    greedy = GreedyEnsembler(
+        problem_type="binary", metric=metric, ensemble_size=30, random_state=np.random.RandomState(0)
+    )
+    greedy.fit(predictions=preds, labels=y)
+    g_err = metric.error(y, greedy.predict_proba(preds))
+
+    hc = HillClimbingEnsembler(
+        problem_type="binary",
+        metric=metric,
+        precision=0.05,
+        max_rounds=20,
+        initial_weights=greedy.model_weights(),
+        random_state=0,
+    )
+    hc.fit(predictions=preds, labels=y)
+    hc_err = metric.error(y, hc.predict_proba(preds))
+    assert hc_err <= g_err + 1e-12
+    assert hc.info()["warm_started"] is True
+    assert hc.info()["init_val_error"] is not None