diff --git a/aeon/base/_estimators/interval_based/base_interval_forest.py b/aeon/base/_estimators/interval_based/base_interval_forest.py
index a1555e0010..cd01081886 100644
--- a/aeon/base/_estimators/interval_based/base_interval_forest.py
+++ b/aeon/base/_estimators/interval_based/base_interval_forest.py
@@ -36,7 +36,13 @@ class BaseIntervalForest(ABC):
     ----------
     base_estimator : BaseEstimator or None, default=None
         scikit-learn BaseEstimator used to build the interval ensemble. If None, use a
-        simple decision tree.
+        simple decision tree, i.e. ``DecisionTreeClassifier(criterion="entropy")`` for
+        classification and ``DecisionTreeRegressor(criterion="squared_error")`` for
+        regression. Versions prior to and including 1.5.0 used
+        ``DecisionTreeRegressor(criterion="absolute_error")`` for regression by
+        default, which scales poorly with the number of cases and can be orders of
+        magnitude slower for larger datasets. Pass such an estimator explicitly to
+        restore the old behaviour.
     n_estimators : int, default=200
         Number of estimators to build for the ensemble.
     interval_selection_method : "random", "supervised" or "random-supervised",
@@ -225,21 +231,34 @@ def _fit(self, X, y):
 
         return self
 
+    def _eval_estimators(self, tasks):
+        """Run a list of joblib ``delayed`` tasks over the estimators.
+
+        Runs sequentially when ``n_jobs == 1`` to avoid the joblib dispatch
+        overhead (a ``Parallel`` object plus per-task wrapping) for what is the
+        default and dominant case, otherwise runs them in parallel. The tasks
+        must be built as a list so that any random draws in their arguments
+        happen in the same order joblib would consume a generator, keeping the
+        result identical to the parallel path.
+        """
+        if self._n_jobs == 1:
+            return [func(*args, **kwargs) for func, args, kwargs in tasks]
+        return Parallel(n_jobs=self._n_jobs, backend=self.parallel_backend)(tasks)
+
     def _predict(self, X):
         if is_regressor(self):
             Xt = self._predict_setup(X)
 
-            y_preds = Parallel(
-                n_jobs=self._n_jobs,
-                backend=self.parallel_backend,
-            )(
-                delayed(self._predict_for_estimator)(
-                    Xt,
-                    self.estimators_[i],
-                    self.intervals_[i],
-                    predict_proba=False,
-                )
-                for i in range(self._n_estimators)
+            y_preds = self._eval_estimators(
+                [
+                    delayed(self._predict_for_estimator)(
+                        Xt,
+                        self.estimators_[i],
+                        self.intervals_[i],
+                        predict_proba=False,
+                    )
+                    for i in range(self._n_estimators)
+                ]
             )
 
             return np.mean(y_preds, axis=0)
@@ -251,14 +270,16 @@ def _predict(self, X):
     def _predict_proba(self, X):
         Xt = self._predict_setup(X)
 
-        y_probas = Parallel(n_jobs=self._n_jobs, backend=self.parallel_backend)(
-            delayed(self._predict_for_estimator)(
-                Xt,
-                self.estimators_[i],
-                self.intervals_[i],
-                predict_proba=True,
-            )
-            for i in range(self._n_estimators)
+        y_probas = self._eval_estimators(
+            [
+                delayed(self._predict_for_estimator)(
+                    Xt,
+                    self.estimators_[i],
+                    self.intervals_[i],
+                    predict_proba=True,
+                )
+                for i in range(self._n_estimators)
+            ]
         )
 
         output = np.sum(y_probas, axis=0) / (
@@ -272,14 +293,16 @@ def _fit_predict(self, X, y) -> np.ndarray:
         if is_regressor(self):
             Xt = self._fit_forest(X, y, save_transformed_data=True)
 
-            p = Parallel(n_jobs=self._n_jobs, backend=self.parallel_backend)(
-                delayed(self._train_estimate_for_estimator)(
-                    Xt,
-                    y,
-                    i,
-                    check_random_state(rng.randint(np.iinfo(np.int32).max)),
-                )
-                for i in range(self._n_estimators)
+            p = self._eval_estimators(
+                [
+                    delayed(self._train_estimate_for_estimator)(
+                        Xt,
+                        y,
+                        i,
+                        check_random_state(rng.randint(np.iinfo(np.int32).max)),
+                    )
+                    for i in range(self._n_estimators)
+                ]
             )
             y_preds, oobs = zip(*p)
 
@@ -314,15 +337,17 @@ def _fit_predict_proba(self, X, y) -> np.ndarray:
 
         rng = check_random_state(self.random_state)
 
-        p = Parallel(n_jobs=self._n_jobs, backend=self.parallel_backend)(
-            delayed(self._train_estimate_for_estimator)(
-                Xt,
-                y,
-                i,
-                check_random_state(rng.randint(np.iinfo(np.int32).max)),
-                probas=True,
-            )
-            for i in range(self._n_estimators)
+        p = self._eval_estimators(
+            [
+                delayed(self._train_estimate_for_estimator)(
+                    Xt,
+                    y,
+                    i,
+                    check_random_state(rng.randint(np.iinfo(np.int32).max)),
+                    probas=True,
+                )
+                for i in range(self._n_estimators)
+            ]
         )
         y_probas, oobs = zip(*p)
 
@@ -352,7 +377,7 @@ def _fit_forest(self, X, y, save_transformed_data=False):
             if is_classifier(self):
                 self._base_estimator = DecisionTreeClassifier(criterion="entropy")
             elif is_regressor(self):
-                self._base_estimator = DecisionTreeRegressor(criterion="absolute_error")
+                self._base_estimator = DecisionTreeRegressor(criterion="squared_error")
             else:
                 raise ValueError(
                     f"{self} must be a scikit-learn compatible classifier or "
@@ -520,7 +545,7 @@ def _fit_forest(self, X, y, save_transformed_data=False):
             or self.max_interval_length == np.inf
         ):
             self._max_interval_length = [self.max_interval_length] * len(Xt)
-        # max_interval_length must be at less than one if it is a float  (proportion of
+        # max_interval_length must be less than one if it is a float (proportion
         # of the series length)
         elif (
             isinstance(self.max_interval_length, float)
@@ -610,7 +635,7 @@ def _fit_forest(self, X, y, save_transformed_data=False):
                                 self._interval_function[i] = True
                             else:
                                 raise ValueError(
-                                    "Individual items in a interval_features list or "
+                                    "Individual items in an interval_features list or "
                                     "tuple must be a transformer or function. Input "
                                     f"{feature} does not contain only transformers and "
                                     f"functions."
@@ -625,7 +650,7 @@ def _fit_forest(self, X, y, save_transformed_data=False):
                         self._interval_features.append([feature])
                     else:
                         raise ValueError(
-                            "Individual items in a interval_features list or tuple "
+                            "Individual items in an interval_features list or tuple "
                             f"must be a transformer or function. Found {feature}"
                         )
         # use basic summary stats by default if None
@@ -646,8 +671,8 @@ def _fit_forest(self, X, y, save_transformed_data=False):
                 )
 
             self._att_subsample_size = [self.att_subsample_size] * len(Xt)
-        # att_subsample_size must be at less than one if it is a float (proportion of
-        # total attributed to subsample)
+        # att_subsample_size must be less than one if it is a float (proportion of
+        # total attributes to subsample)
         elif isinstance(self.att_subsample_size, float):
             if self.att_subsample_size > 1 or self.att_subsample_size <= 0:
                 raise ValueError(
@@ -804,17 +829,16 @@ def _fit_forest(self, X, y, save_transformed_data=False):
                 train_time < time_limit
                 and self._n_estimators < self.contract_max_n_estimators
             ):
-                fit = Parallel(
-                    n_jobs=self._n_jobs,
-                    backend=self.parallel_backend,
-                )(
-                    delayed(self._fit_estimator)(
-                        Xt,
-                        y,
-                        rng.randint(np.iinfo(np.int32).max),
-                        save_transformed_data=save_transformed_data,
-                    )
-                    for _ in range(self._n_jobs)
+                fit = self._eval_estimators(
+                    [
+                        delayed(self._fit_estimator)(
+                            Xt,
+                            y,
+                            rng.randint(np.iinfo(np.int32).max),
+                            save_transformed_data=save_transformed_data,
+                        )
+                        for _ in range(self._n_jobs)
+                    ]
                 )
 
                 (
@@ -832,17 +856,16 @@ def _fit_forest(self, X, y, save_transformed_data=False):
         else:
             self._n_estimators = self.n_estimators
 
-            fit = Parallel(
-                n_jobs=self._n_jobs,
-                backend=self.parallel_backend,
-            )(
-                delayed(self._fit_estimator)(
-                    Xt,
-                    y,
-                    rng.randint(np.iinfo(np.int32).max),
-                    save_transformed_data=save_transformed_data,
-                )
-                for _ in range(self._n_estimators)
+            fit = self._eval_estimators(
+                [
+                    delayed(self._fit_estimator)(
+                        Xt,
+                        y,
+                        rng.randint(np.iinfo(np.int32).max),
+                        save_transformed_data=save_transformed_data,
+                    )
+                    for _ in range(self._n_estimators)
+                ]
             )
 
             (
@@ -1063,11 +1086,9 @@ def _predict_setup(self, X):
         return Xt
 
     def _predict_for_estimator(self, Xt, estimator, intervals, predict_proba=False):
-        interval_features = np.empty((Xt[0].shape[0], 0))
-
-        for r in range(len(Xt)):
-            f = intervals[r].transform(Xt[r])
-            interval_features = np.hstack((interval_features, f))
+        interval_features = np.hstack(
+            [intervals[r].transform(Xt[r]) for r in range(len(Xt))]
+        )
 
         if isinstance(self.replace_nan, str) and self.replace_nan.lower() == "nan":
             interval_features = np.nan_to_num(
diff --git a/aeon/classification/interval_based/_drcif.py b/aeon/classification/interval_based/_drcif.py
index f31ce4a722..314c7b82ae 100644
--- a/aeon/classification/interval_based/_drcif.py
+++ b/aeon/classification/interval_based/_drcif.py
@@ -270,7 +270,7 @@ def _get_test_params(cls, parameter_set="default"):
                     previously generated results where the default set of parameters
                     cannot produce suitable probability estimates
                 "contracting" - used in classifiers that set the
-                    "capability:contractable" tag to True to test contacting
+                    "capability:contractable" tag to True to test contracting
                     functionality
                 "train_estimate" - used in some classifiers that set the
                     "capability:train_estimate" tag to True to allow for more efficient
diff --git a/aeon/regression/interval_based/_cif.py b/aeon/regression/interval_based/_cif.py
index c26deafbeb..f2690f0da2 100644
--- a/aeon/regression/interval_based/_cif.py
+++ b/aeon/regression/interval_based/_cif.py
@@ -31,7 +31,11 @@ class CanonicalIntervalForestRegressor(BaseIntervalForest, BaseRegressor):
     ----------
     base_estimator : BaseEstimator or None, default=None
         scikit-learn BaseEstimator used to build the interval ensemble. If None, use a
-        simple decision tree.
+        simple decision tree, i.e.
+        ``DecisionTreeRegressor(criterion="squared_error")``. Versions prior to and
+        including 1.5.0 used ``DecisionTreeRegressor(criterion="absolute_error")``
+        by default, which scales poorly with the number of cases. Pass such an
+        estimator explicitly to restore the old behaviour.
     n_estimators : int, default=200
         Number of estimators to build for the ensemble.
     n_intervals : int, str, list or tuple, default="sqrt"
diff --git a/aeon/regression/interval_based/_drcif.py b/aeon/regression/interval_based/_drcif.py
index a8ef6baf55..9f123ed89a 100644
--- a/aeon/regression/interval_based/_drcif.py
+++ b/aeon/regression/interval_based/_drcif.py
@@ -45,7 +45,11 @@ class DrCIFRegressor(BaseIntervalForest, BaseRegressor):
     ----------
     base_estimator : BaseEstimator or None, default=None
         scikit-learn BaseEstimator used to build the interval ensemble. If None, use a
-        simple decision tree.
+        simple decision tree, i.e.
+        ``DecisionTreeRegressor(criterion="squared_error")``. Versions prior to and
+        including 1.5.0 used ``DecisionTreeRegressor(criterion="absolute_error")``
+        by default, which scales poorly with the number of cases. Pass such an
+        estimator explicitly to restore the old behaviour.
     n_estimators : int, default=200
         Number of estimators to build for the ensemble.
     n_intervals : int, str, list or tuple, default="sqrt"
@@ -240,7 +244,7 @@ def _get_test_params(cls, parameter_set="default"):
                     previously generated results where the default set of parameters
                     cannot produce suitable probability estimates
                 "contracting" - used in classifiers that set the
-                    "capability:contractable" tag to True to test contacting
+                    "capability:contractable" tag to True to test contracting
                     functionality
 
         Returns
diff --git a/aeon/regression/interval_based/_interval_forest.py b/aeon/regression/interval_based/_interval_forest.py
index 568a7f5d58..a7d82421e0 100644
--- a/aeon/regression/interval_based/_interval_forest.py
+++ b/aeon/regression/interval_based/_interval_forest.py
@@ -25,7 +25,11 @@ class IntervalForestRegressor(BaseIntervalForest, BaseRegressor):
     ----------
     base_estimator : BaseEstimator or None, default=None
         scikit-learn BaseEstimator used to build the interval ensemble. If None, use a
-        simple decision tree.
+        simple decision tree, i.e.
+        ``DecisionTreeRegressor(criterion="squared_error")``. Versions prior to and
+        including 1.5.0 used ``DecisionTreeRegressor(criterion="absolute_error")``
+        by default, which scales poorly with the number of cases. Pass such an
+        estimator explicitly to restore the old behaviour.
     n_estimators : int, default=200
         Number of estimators to build for the ensemble.
     interval_selection_method : "random", default="random"
diff --git a/aeon/regression/interval_based/_rise.py b/aeon/regression/interval_based/_rise.py
index 9f5f0f8422..8a8320ce0a 100644
--- a/aeon/regression/interval_based/_rise.py
+++ b/aeon/regression/interval_based/_rise.py
@@ -27,7 +27,11 @@ class RandomIntervalSpectralEnsembleRegressor(BaseIntervalForest, BaseRegressor)
     ----------
     base_estimator : BaseEstimator or None, default=None
         scikit-learn BaseEstimator used to build the interval ensemble. If None, use a
-        simple decision tree.
+        simple decision tree, i.e.
+        ``DecisionTreeRegressor(criterion="squared_error")``. Versions prior to and
+        including 1.5.0 used ``DecisionTreeRegressor(criterion="absolute_error")``
+        by default, which scales poorly with the number of cases. Pass such an
+        estimator explicitly to restore the old behaviour.
     n_estimators : int, default=200
         Number of estimators to build for the ensemble.
     min_interval_length : int, float, list, or tuple, default=3
diff --git a/aeon/regression/interval_based/_tsf.py b/aeon/regression/interval_based/_tsf.py
index df7bc3b704..695e27bd27 100644
--- a/aeon/regression/interval_based/_tsf.py
+++ b/aeon/regression/interval_based/_tsf.py
@@ -32,7 +32,11 @@ class TimeSeriesForestRegressor(BaseIntervalForest, BaseRegressor):
     ----------
     base_estimator : BaseEstimator or None, default=None
         scikit-learn BaseEstimator used to build the interval ensemble. If None, use a
-        simple decision tree.
+        simple decision tree, i.e.
+        ``DecisionTreeRegressor(criterion="squared_error")``. Versions prior to and
+        including 1.5.0 used ``DecisionTreeRegressor(criterion="absolute_error")``
+        by default, which scales poorly with the number of cases. Pass such an
+        estimator explicitly to restore the old behaviour.
     n_estimators : int, default=200
         Number of estimators to build for the ensemble.
     n_intervals : int, str, list or tuple, default="sqrt"
diff --git a/aeon/testing/expected_results/expected_regressor_results.py b/aeon/testing/expected_results/expected_regressor_results.py
index 8254e480dc..f9bdb421b1 100644
--- a/aeon/testing/expected_results/expected_regressor_results.py
+++ b/aeon/testing/expected_results/expected_regressor_results.py
@@ -6,16 +6,16 @@
 
 univariate_expected_results = {
     "CanonicalIntervalForestRegressor": [
+        0.0302,
+        0.0301,
         0.0277,
-        0.027,
-        0.0243,
-        0.0187,
-        0.0186,
-        0.0109,
-        0.0224,
-        0.0173,
-        0.0195,
-        0.0215,
+        0.0276,
+        0.0197,
+        0.0309,
+        0.0279,
+        0.0203,
+        0.0208,
+        0.0329,
     ],
     "Catch22Regressor": [
         0.0164,
@@ -54,16 +54,16 @@
         0.0243,
     ],
     "IntervalForestRegressor": [
-        0.0309,
-        0.0317,
-        0.0141,
-        0.0139,
-        0.0198,
-        0.0133,
-        0.026,
-        0.0151,
-        0.0251,
-        0.008,
+        0.0209,
+        0.0367,
+        0.0188,
+        0.0165,
+        0.0276,
+        0.0159,
+        0.0298,
+        0.0186,
+        0.0323,
+        0.0132,
     ],
     "KNeighborsTimeSeriesRegressor": [
         0.0069,
@@ -162,16 +162,16 @@
         0.0234,
     ],
     "RandomIntervalSpectralEnsembleRegressor": [
-        0.0105,
-        0.0194,
-        0.0206,
-        0.0203,
-        0.0149,
-        0.0173,
-        0.0174,
-        0.0144,
-        0.0129,
-        0.0258,
+        0.0231,
+        0.0246,
+        0.0283,
+        0.0323,
+        0.0232,
+        0.0224,
+        0.0233,
+        0.0327,
+        0.0246,
+        0.0265,
     ],
     "RocketRegressor": [
         0.026,
@@ -222,30 +222,30 @@
         0.0074,
     ],
     "TimeSeriesForestRegressor": [
-        0.0428,
-        0.0233,
         0.017,
+        0.0292,
+        0.0216,
         0.018,
-        0.018,
+        0.0438,
         0.018,
         0.0363,
-        0.0144,
-        0.0214,
-        0.0059,
+        0.0203,
+        0.0273,
+        0.0069,
     ],
 }
 multivariate_expected_results = {
     "CanonicalIntervalForestRegressor": [
-        0.0261,
-        0.0202,
-        0.1841,
-        0.0258,
-        0.0604,
-        0.1304,
-        0.0959,
-        0.0973,
-        0.0848,
-        -0.1254,
+        0.0798,
+        0.0755,
+        0.1795,
+        0.1158,
+        0.0884,
+        0.1452,
+        0.1303,
+        0.1665,
+        0.0785,
+        0.005,
     ],
     "Catch22Regressor": [
         -0.0209,
@@ -284,16 +284,16 @@
         -0.0986,
     ],
     "IntervalForestRegressor": [
-        0.1287,
-        -0.0063,
-        -0.1061,
-        0.0114,
-        -0.0838,
-        0.1019,
-        0.0477,
-        -0.0433,
-        0.2159,
-        -0.0108,
+        0.1347,
+        0.0636,
+        -0.1039,
+        0.0073,
+        -0.0003,
+        0.1399,
+        0.0834,
+        -0.0232,
+        0.1065,
+        -0.022,
     ],
     "KNeighborsTimeSeriesRegressor": [
         0.1805,
@@ -392,16 +392,16 @@
         -0.124,
     ],
     "RandomIntervalSpectralEnsembleRegressor": [
-        0.0941,
-        0.1009,
-        -0.0619,
-        0.0621,
-        0.06,
-        0.1229,
-        0.0327,
-        -0.031,
-        -0.0026,
-        -0.0074,
+        0.0545,
+        0.1026,
+        0.1199,
+        -0.0065,
+        0.1737,
+        0.102,
+        0.0164,
+        -0.0182,
+        -0.019,
+        0.1003,
     ],
     "RocketRegressor": [
         0.0869,
@@ -440,15 +440,15 @@
         -0.1682,
     ],
     "TimeSeriesForestRegressor": [
-        0.0848,
-        -0.0445,
-        0.1356,
+        0.0904,
+        -0.1176,
+        0.2084,
         -0.0132,
-        0.0005,
-        0.3315,
-        -0.247,
-        -0.089,
-        0.3315,
+        0.0732,
+        -0.021,
+        0.0566,
+        0.124,
+        0.3893,
         -0.1123,
     ],
 }
diff --git a/aeon/transformations/collection/feature_based/_catch22.py b/aeon/transformations/collection/feature_based/_catch22.py
index 7baac1b3b1..67d91f39a3 100644
--- a/aeon/transformations/collection/feature_based/_catch22.py
+++ b/aeon/transformations/collection/feature_based/_catch22.py
@@ -287,28 +287,82 @@ def _transform(self, X, y=None):
                     stacklevel=2,
                 )
 
-        c22_list = Parallel(
-            n_jobs=n_jobs, backend=self.parallel_backend, prefer="threads"
-        )(
-            delayed(
-                self._transform_case_pycatch22
-                if use_pycatch22_transform
-                else self._transform_case
-            )(
-                X[i],
-                f_idx,
-                features,
-            )
-            for i in range(n_cases)
-        )
+        if (
+            not use_pycatch22_transform
+            and n_jobs == 1
+            and isinstance(X, np.ndarray)
+            and X.ndim == 3
+        ):
+            # transform the whole collection in a single numba call, avoiding
+            # Python-level per-case overhead. This is the hot path for
+            # interval-based forests, which call transform once per interval with
+            # a single thread. Produces the same output as _transform_case.
+            c22_list = self._transform_batch_equal_length(X, f_idx)
+        else:
+            if use_pycatch22_transform:
+                func = self._transform_case_pycatch22
+                case_args = [(X[i], f_idx, features) for i in range(n_cases)]
+            else:
+                # The two welch power-spectrum features (indices 15 and 20) each
+                # need np.fft.fft of the mean-centred series. np.fft.fft has a
+                # high fixed per-call cost on short interval series, so compute it
+                # once for the whole batch (all cases share a length) instead of
+                # once per case. The result is bit-identical: a 1D FFT equals the
+                # matching row of the batched FFT, and the subtracted mean uses
+                # the same numba mean().
+                fft_cache = self._welch_fft_cache(X, f_idx, n_cases)
+                # The autocorrelation features (2, 3, 8, 10, 12) share the FFT
+                # twiddle table, which depends only on the (common) series length
+                # - compute it once and reuse instead of rebuilding it in every
+                # _compute_autocorrelations call.
+                ac_tw, ac_nfft = self._ac_twiddle_cache(X, f_idx)
+                func = self._transform_case
+                case_args = [
+                    (
+                        X[i],
+                        f_idx,
+                        features,
+                        None if fft_cache is None else fft_cache[i],
+                        ac_tw,
+                        ac_nfft,
+                    )
+                    for i in range(n_cases)
+                ]
+
+            # Run cases sequentially when not parallelising: a joblib Parallel
+            # still wraps every task in delayed() and copies it, which is pure
+            # overhead here.
+            if n_jobs == 1:
+                c22_list = [func(*args) for args in case_args]
+            else:
+                c22_list = Parallel(
+                    n_jobs=n_jobs, backend=self.parallel_backend, prefer="threads"
+                )(delayed(func)(*args) for args in case_args)
 
-        c22_array = np.array(c22_list)
+        c22_array = np.asarray(c22_list)
         if self.replace_nans:
             c22_array = np.nan_to_num(c22_array, False, 0, 0, 0)
 
         return c22_array
 
-    def _transform_case(self, X, f_idx, features):
+    def _ac_twiddle_cache(self, X, f_idx):
+        """Precompute the FFT twiddle table shared by the autocorrelation features.
+
+        Features 2, 3, 8, 10 and 12 all need the autocorrelation, whose FFT twiddle
+        factors depend only on the (common) series length. Build them once here and
+        reuse across every case instead of rebuilding per autocorrelation call.
+        Returns (None, 0) for np-list input or when no autocorrelation feature runs.
+        """
+        if not (isinstance(X, np.ndarray) and X.ndim == 3):
+            return None, 0
+        if not any(fi in (2, 3, 8, 10, 12) for fi in f_idx):
+            return None, 0
+        nfft = _ac_nfft(X.shape[2])
+        return _ac_twiddles(nfft), nfft
+
+    def _transform_case(
+        self, X, f_idx, features, fft_cache=None, ac_tw=None, ac_nfft=0
+    ):
         c22 = np.zeros(len(f_idx) * len(X))
 
         if hasattr(self, "_transform_features") and len(
@@ -359,21 +413,35 @@ def _transform_case(self, X, f_idx, features):
                     if smean is None:
                         smean = mean(series)
                     if fft is None:
-                        nfft = int(
-                            np.power(2, np.ceil(np.log(len(series)) / np.log(2)))
-                        )
-                        fft = np.fft.fft(series - smean, n=nfft)
+                        if fft_cache is not None:
+                            fft = fft_cache[i]
+                        else:
+                            nfft = int(
+                                np.power(2, np.ceil(np.log(len(series)) / np.log(2)))
+                            )
+                            fft = np.fft.fft(series - smean, n=nfft)
                     args = [series, fft]
                 elif feature == 2 or feature == 3:
                     if ac is None:
-                        ac = _compute_autocorrelations(series)
+                        ac = (
+                            _autocorrelations_with_tw(series, ac_tw, ac_nfft)
+                            if ac_tw is not None
+                            else _compute_autocorrelations(series)
+                        )
                     args = [ac, len(series)]
                 elif feature == 12 or feature == 10 or feature == 8:
                     if ac is None:
-                        ac = _compute_autocorrelations(series)
+                        ac = (
+                            _autocorrelations_with_tw(series, ac_tw, ac_nfft)
+                            if ac_tw is not None
+                            else _compute_autocorrelations(series)
+                        )
                     if acfz is None:
                         acfz = _ac_first_zero(ac)
-                    args = [series, acfz]
+                    if feature == 12:
+                        args = [series, acfz, ac_tw, ac_nfft]
+                    else:
+                        args = [series, acfz]
                 if feature == 22:
                     c22[dim + n] = smean
                 elif feature == 23:
@@ -383,6 +451,48 @@ def _transform_case(self, X, f_idx, features):
 
         return c22
 
+    def _welch_fft_cache(self, X, f_idx, n_cases):
+        """Batch the welch-feature FFT across all cases when it will be used.
+
+        Returns an array of shape (n_cases, n_channels, nfft) with the FFT of each
+        mean-centred series, or None if the welch features (15/20) are absent, the
+        input is not an equal-length 3D array, or attribute-skipping means neither
+        welch feature will be computed. Matches the per-case computation exactly.
+        """
+        if (15 not in f_idx and 20 not in f_idx) or not (
+            isinstance(X, np.ndarray) and X.ndim == 3
+        ):
+            return None
+
+        n_channels = X.shape[1]
+
+        # Respect the transform_features skip mask set for efficient predictions:
+        # if no welch output is going to be produced, don't build the cache.
+        tf = getattr(self, "_transform_features", None)
+        if tf is not None and len(tf) == len(f_idx) * n_channels:
+            will_run = False
+            for c in range(n_channels):
+                base = c * len(f_idx)
+                for n, feat in enumerate(f_idx):
+                    if (feat == 15 or feat == 20) and tf[base + n]:
+                        will_run = True
+                        break
+                if will_run:
+                    break
+            if not will_run:
+                return None
+
+        m = X.shape[2]
+        nfft = int(np.power(2, np.ceil(np.log(m) / np.log(2))))
+        flat = X.reshape(n_cases * n_channels, m)
+        # Use the same numba mean() as _transform_case so the centred series, and
+        # therefore the FFT, are identical to the per-case path.
+        fmeans = np.empty(flat.shape[0])
+        for i in range(flat.shape[0]):
+            fmeans[i] = mean(flat[i])
+        fft = np.fft.fft(flat - fmeans[:, np.newaxis], n=nfft, axis=1)
+        return fft.reshape(n_cases, n_channels, nfft)
+
     def _transform_case_pycatch22(self, X, f_idx, features):
         c22 = np.zeros(len(f_idx) * len(X))
 
@@ -417,6 +527,65 @@ def _transform_case_pycatch22(self, X, f_idx, features):
 
         return c22
 
+    def _transform_batch_equal_length(self, X, f_idx):
+        """Transform a whole equal length collection in a single numba call.
+
+        Produces the same output as calling _transform_case for each case, without
+        the Python-level per-case overhead. The FFT used by the two SP_Summaries
+        features is precomputed with numpy for all cases, as np.fft is not
+        supported in numba.
+
+        Note that as the feature functions are compiled with fastmath enabled,
+        calling them from another numba function rather than from Python can
+        produce differences in the last few float digits.
+        """
+        n_cases, n_channels, n_timepoints = X.shape
+        n_out = len(f_idx) * n_channels
+
+        transform_features = getattr(self, "_transform_features", None)
+        if transform_features is not None and len(transform_features) == n_out:
+            transform_features = np.asarray(transform_features, dtype=np.bool_)
+        else:
+            transform_features = np.ones(n_out, dtype=np.bool_)
+
+        if 15 in f_idx or 20 in f_idx:
+            # the series mean must be calculated with a Python call to mean() to
+            # exactly match the per-case path here. The FFT DC term is a
+            # catastrophic cancellation of sum(series - mean), so even a one ULP
+            # difference in the mean can significantly change the value of the
+            # SP_Summaries features for short series.
+            means = np.empty((n_cases, n_channels))
+            for c in range(n_cases):
+                for i in range(n_channels):
+                    means[c, i] = mean(X[c, i])
+            nfft = int(np.power(2, np.ceil(np.log(n_timepoints) / np.log(2))))
+            ffts = np.fft.fft(X - means[:, :, np.newaxis], n=nfft, axis=2)
+        else:
+            ffts = np.zeros((n_cases, n_channels, 1), dtype=np.complex128)
+
+        c22 = _transform_collection_3d(
+            X,
+            ffts,
+            np.asarray(f_idx, dtype=np.int64),
+            transform_features,
+            self.outlier_norm,
+        )
+
+        # mean and standard deviation for catch24 are calculated here, matching
+        # the functions used in _transform_case
+        for n, feature in enumerate(f_idx):
+            if feature == 22 or feature == 23:
+                for i in range(n_channels):
+                    col = i * len(f_idx) + n
+                    if not transform_features[col]:
+                        continue
+                    for c in range(n_cases):
+                        c22[c, col] = (
+                            mean(X[c, i]) if feature == 22 else np.std(X[c, i])
+                        )
+
+        return c22
+
     @property
     def get_features_arguments(self):
         """Return feature names for the estimators features argument."""
@@ -610,49 +779,25 @@ def _SB_BinaryStats_mean_longstretch1(X, smean):
     @staticmethod
     @njit(fastmath=True, cache=True)
     def _SB_MotifThree_quantile_hh(X):
+        # Entropy of adjacent-symbol transitions after 3-symbol coarse-graining.
+        # The original built per-symbol position lists (r1) and per-transition
+        # lists (r2) and trimmed the final position; that is exactly a 3x3 count
+        # of adjacent pairs (yt[p], yt[p+1]) for p in 0..len-2 (the last position
+        # has no successor, which the range already excludes).
         alphabet_size = 3
-        yt = np.zeros(len(X), dtype=np.int32)
+        n = len(X)
+        yt = np.zeros(n, dtype=np.int32)
         _sb_coarsegrain(X, 3, yt)
-        r1 = [np.zeros(len(X), np.int32) for i in range(alphabet_size)]
-        sizes_r1 = np.zeros(alphabet_size, np.int32)
-        for i in range(alphabet_size):
-            r_idx = 0
-            sizes_r1[i] = 0
-            for j in range(len(X)):
-                if yt[j] == i + 1:
-                    r1[i][r_idx] = j
-                    r_idx += 1
-                    sizes_r1[i] += 1
 
-        for i in range(alphabet_size):
-            if sizes_r1[i] != 0 and r1[i][sizes_r1[i] - 1] == len(X) - 1:
-                tmp_ar = np.zeros(sizes_r1[i], np.int32)
-                # isn't this doing the same thing?
-                for x in range(sizes_r1[i]):
-                    tmp_ar[x] = r1[i][x]
-                for y in range(sizes_r1[i] - 1):
-                    r1[i][y] = tmp_ar[y]
-                sizes_r1[i] -= 1
-
-        r2 = [
-            [np.zeros(len(X), np.int32) for j in range(alphabet_size)]
-            for i in range(alphabet_size)
-        ]
-        sizes_r2 = [np.zeros(alphabet_size, np.int32) for i in range(alphabet_size)]
-        out2 = [np.zeros(alphabet_size, np.float64) for i in range(alphabet_size)]
+        counts = np.zeros((alphabet_size, alphabet_size), dtype=np.float64)
+        for p in range(n - 1):
+            counts[yt[p] - 1][yt[p + 1] - 1] += 1.0
 
+        out2 = np.zeros((alphabet_size, alphabet_size), dtype=np.float64)
         for i in range(alphabet_size):
             for j in range(alphabet_size):
-                sizes_r2[i][j] = 0
-                dynamic_idx = 0
-                for k in range(sizes_r1[i]):
-                    tmp_idx = yt[r1[i][k] + 1]
-                    if tmp_idx == j + 1:
-                        r2[i][j][dynamic_idx] = r1[i][k]
-                        dynamic_idx += 1
-                        sizes_r2[i][j] += 1
-                tmp = np.float64(sizes_r2[i][j]) / (np.float64(len(X)) - 1.0)
-                out2[i][j] = tmp
+                out2[i][j] = counts[i][j] / (np.float64(n) - 1.0)
+
         hh = 0.0
         for i in range(alphabet_size):
             f = 0.0
@@ -663,13 +808,20 @@ def _SB_MotifThree_quantile_hh(X):
         return hh
 
     @staticmethod
-    def _FC_LocalSimple_mean1_tauresrat(X, acfz):
+    @njit(fastmath=True, cache=True)
+    def _FC_LocalSimple_mean1_tauresrat(X, acfz, ac_tw, ac_nfft):
         # Change in correlation length after iterative differencing.
         if len(X) < 2:
             return 0
         res = _local_simple_mean(X, 1)
 
-        ac = _compute_autocorrelations(res)
+        # Reuse the transform's precomputed twiddle table when the differenced
+        # series lands on the same FFT size (it does except for length 2**k + 1);
+        # otherwise fall back to computing them. Bit-identical either way.
+        if ac_tw is not None and _ac_nfft(len(res)) == ac_nfft:
+            ac = _autocorrelations_with_tw(res, ac_tw, ac_nfft)
+        else:
+            ac = _compute_autocorrelations(res)
         return _ac_first_zero(ac) / acfz
 
     @staticmethod
@@ -906,8 +1058,31 @@ def _long_stretch(X_binary, val):
     return max_stretch
 
 
+@njit(fastmath=True, cache=True)
+def _bit_update(bit, i, n):
+    # Fenwick tree point increment at 1-indexed position i (size n).
+    while i <= n:
+        bit[i] += 1
+        i += i & (-i)
+
+
+@njit(fastmath=True, cache=True)
+def _bit_select(bit, k, n, logn):
+    # Return the 1-indexed position of the k-th smallest present element.
+    pos = 0
+    pw = 1 << logn
+    while pw > 0:
+        nxt = pos + pw
+        if nxt <= n and bit[nxt] < k:
+            pos = nxt
+            k -= bit[nxt]
+        pw >>= 1
+    return pos + 1
+
+
 @njit(fastmath=True, cache=True)
 def _outlier_include(X):
+    n = len(X)
     total = 0
     threshold = 0
     for v in X:
@@ -923,26 +1098,64 @@ def _outlier_include(X):
     means = np.zeros(num_thresholds)
     dists = np.zeros(num_thresholds)
     medians = np.zeros(num_thresholds)
-    for i in range(num_thresholds):
-        d = i * 0.01
 
-        count = 0
-        r = np.zeros(len(X))
-        for n in range(len(X)):
-            if X[n] >= d:
-                r[count] = n + 1
-                count += 1
+    # For each position, membership "X[pos] >= i * 0.01" holds for a contiguous
+    # range of thresholds i = 0 .. k. Bucket each position by its highest such k
+    # (found with the same comparison the original loop used) so thresholds can
+    # be swept without rescanning X. Positions sharing a k are chained via nxt.
+    head = np.full(num_thresholds, -1, dtype=np.int64)
+    nxt = np.full(n, -1, dtype=np.int64)
+    for pos in range(n):
+        x = X[pos]
+        if x < 0:
+            continue
+        k = int(x / 0.01)
+        while k * 0.01 > x:
+            k -= 1
+        while (k + 1) * 0.01 <= x:
+            k += 1
+        if k > num_thresholds - 1:
+            k = num_thresholds - 1
+        nxt[pos] = head[k]
+        head[k] = pos
+
+    logn = 0
+    while (1 << (logn + 1)) <= n:
+        logn += 1
+    bit = np.zeros(n + 1, dtype=np.int64)
+
+    # Sweep thresholds from high to low, inserting positions as they enter the
+    # set. count/min/max update incrementally; the median is read from the
+    # Fenwick tree. This reproduces the per-threshold r-array statistics exactly.
+    count = 0
+    cur_min = n + 1
+    cur_max = -1
+    for i in range(num_thresholds - 1, -1, -1):
+        pos = head[i]
+        while pos != -1:
+            rv = pos + 1
+            _bit_update(bit, rv, n)
+            if rv < cur_min:
+                cur_min = rv
+            if rv > cur_max:
+                cur_max = rv
+            count += 1
+            pos = nxt[pos]
 
         if count == 0:
             continue
 
-        diff = np.zeros(count - 1)
-        for n in range(len(diff)):
-            diff[n] = r[n + 1] - r[n]
+        # mean of consecutive gaps telescopes to (last - first) / (count - 1).
+        means[i] = (cur_max - cur_min) / (count - 1) if count > 1 else 9999999999
+        dists[i] = (count - 1) * 100 / total
 
-        means[i] = np.mean(diff) if len(diff) > 0 else 9999999999
-        dists[i] = len(diff) * 100 / total
-        medians[i] = np.median(r[:count]) / (len(X) / 2) - 1
+        if count % 2 == 1:
+            median_val = _bit_select(bit, (count + 1) // 2, n, logn)
+        else:
+            a = _bit_select(bit, count // 2, n, logn)
+            b = _bit_select(bit, count // 2 + 1, n, logn)
+            median_val = (a + b) / 2.0
+        medians[i] = median_val / (n / 2) - 1
 
     mj = 0
     fbi = num_thresholds - 1
@@ -1074,21 +1287,46 @@ def _fluct_prop(X, og_length, dfa):
                 buffer[n][j] = X[count]
                 count += 1
 
-        d = np.zeros(tau)
-        for n in range(tau):
-            d[n] = n + 1
+        # Bespoke least squares over the fixed grid d = [1..tau]: its moments are
+        # closed-form and computed once per tau instead of re-summed per window.
+        # For DFA the fluctuation is the residual sum of squares, obtained from
+        # sums (SSR = sumy2 - c1*sumxy - c2*sumy) without detrending in place.
+        sumx = tau * (tau + 1) / 2.0
+        sumx2 = tau * (tau + 1) * (2 * tau + 1) / 6.0
+        denom = tau * sumx2 - sumx * sumx
 
         for n in range(buff_size):
-            c1, c2 = _linear_regression(d, buffer[n], tau, 0)
-
+            sumy = 0.0
+            sumxy = 0.0
+            sumy2 = 0.0
             for j in range(tau):
-                buffer[n][j] = buffer[n][j] - (c1 * (j + 1) + c2)
+                v = buffer[n][j]
+                sumy += v
+                sumxy += (j + 1) * v
+                sumy2 += v * v
+
+            if denom == 0:
+                c1 = 0.0
+                c2 = 0.0
+            else:
+                c1 = (tau * sumxy - sumx * sumy) / denom
+                c2 = (sumy * sumx2 - sumx * sumxy) / denom
 
             if dfa:
-                for j in range(tau):
-                    f[i] += buffer[n][j] * buffer[n][j]
+                ssr = sumy2 - c1 * sumxy - c2 * sumy
+                if ssr < 0.0:
+                    ssr = 0.0
+                f[i] += ssr
             else:
-                f[i] += np.power(np.max(buffer[n]) - np.min(buffer[n]), 2)
+                rmax = -np.inf
+                rmin = np.inf
+                for j in range(tau):
+                    resid = buffer[n][j] - (c1 * (j + 1) + c2)
+                    if resid > rmax:
+                        rmax = resid
+                    if resid < rmin:
+                        rmin = resid
+                f[i] += np.power(rmax - rmin, 2)
 
         if dfa:
             f[i] = np.sqrt(f[i] / (buff_size * tau))
@@ -1396,24 +1634,35 @@ def _verify_features(features, catch24):
 
 
 @njit(fastmath=True, cache=True)
-def _compute_autocorrelations(X):
-    mean = np.mean(X)
-    nFFT = int(np.log2(len(X)))
-    if 2**nFFT == len(X):
-        nFFT = len(X) * 2
+def _ac_nfft(n):
+    nFFT = int(np.log2(n))
+    if 2**nFFT == n:
+        nFFT = n * 2
     else:
         nFFT = (2 ** (nFFT + 1)) * 2
-    F = np.zeros(nFFT * 2, dtype=np.complex128)
-    for i in range(len(X)):
-        F[i] = complex(X[i] - mean, 0.0)
-    for i in range(len(X), nFFT):
-        F[i] = complex(0.0, 0.0)
+    return nFFT
+
+
+@njit(fastmath=True, cache=True)
+def _ac_twiddles(nFFT):
+    # FFT twiddle factors; depend only on nFFT so can be reused across the many
+    # autocorrelation calls of a transform (all series share a length).
     tw = np.zeros(nFFT * 2, dtype=np.complex128)
-    # twiddles
     PI = np.pi
     for i in range(nFFT):
         tmp = 0.0 - PI * i / nFFT * 1j
         tw[i] = np.exp(tmp)
+    return tw
+
+
+@njit(fastmath=True, cache=True)
+def _autocorrelations_with_tw(X, tw, nFFT):
+    mean = np.mean(X)
+    F = np.zeros(nFFT * 2, dtype=np.complex128)
+    for i in range(len(X)):
+        F[i] = complex(X[i] - mean, 0.0)
+    for i in range(len(X), nFFT):
+        F[i] = complex(0.0, 0.0)
     F = _fft(F, tw)
     # dot multiply
     F = np.multiply(F, np.conj(F))
@@ -1426,6 +1675,12 @@ def _compute_autocorrelations(X):
     return out
 
 
+@njit(fastmath=True, cache=True)
+def _compute_autocorrelations(X):
+    nFFT = _ac_nfft(len(X))
+    return _autocorrelations_with_tw(X, _ac_twiddles(nFFT), nFFT)
+
+
 @njit(fastmath=True, cache=True)
 def _fft(a, tw):
     n = a.shape[0]
@@ -1471,30 +1726,197 @@ def _sb_coarsegrain(y, num_groups, labels):
     for i in range(num_groups + 1):
         ls[i] = start
         start += step_size
+    # Sort once and reuse for every quantile: _quantile would otherwise re-sort
+    # the same array num_groups + 1 times.
+    tmp = np.sort(y)
     for i in range(num_groups + 1):
-        th[i] = _quantile(y, ls[i])
+        th[i] = _quantile_sorted(tmp, ls[i])
     th[0] -= 1
-    for i in range(num_groups):
-        for j in range(len(y)):
+    # Single pass over y: the group ranges (th[i], th[i+1]] are disjoint, so the
+    # first match is the only match (equivalent to the original num_groups passes).
+    for j in range(len(y)):
+        for i in range(num_groups):
             if y[j] > th[i] and y[j] <= th[i + 1]:
                 labels[j] = i + 1
+                break
 
 
 @njit(fastmath=True, cache=True)
-def _quantile(X, quant):
-    tmp = np.sort(X)
-    q = 0.5 / len(X)
+def _quantile_sorted(tmp, quant):
+    # Linear-interpolation quantile of an already-sorted array.
+    n = len(tmp)
+    q = 0.5 / n
     if quant < q:
-        value = tmp[0]
-        return value
+        return tmp[0]
     elif quant > (1 - q):
-        value = tmp[len(X) - 1]
-        return value
+        return tmp[n - 1]
 
-    quant_idx = len(X) * quant - 0.5
+    quant_idx = n * quant - 0.5
     idx_left = int(np.floor(quant_idx))
     idx_right = int(np.ceil(quant_idx))
-    value = tmp[idx_left] + (quant_idx - idx_left) * (
-        tmp[idx_right] - tmp[idx_left]
-    ) / (idx_right - idx_left)
-    return value
+    return tmp[idx_left] + (quant_idx - idx_left) * (tmp[idx_right] - tmp[idx_left]) / (
+        idx_right - idx_left
+    )
+
+
+@njit(fastmath=True, cache=True)
+def _quantile(X, quant):
+    return _quantile_sorted(np.sort(X), quant)
+
+
+# module level references to the feature staticmethods, required to call them
+# from the numba function below
+_co_f1ecac = Catch22._CO_f1ecac
+_co_first_min_ac = Catch22._CO_FirstMin_ac
+_co_histogram_ami_even_2_5 = Catch22._CO_HistogramAMI_even_2_5
+_co_trev_1_num = Catch22._CO_trev_1_num
+_md_hrv_classic_pnn40 = Catch22._MD_hrv_classic_pnn40
+_sb_binary_stats_mean_longstretch1 = Catch22._SB_BinaryStats_mean_longstretch1
+_sb_transition_matrix_3ac_sumdiagcov = Catch22._SB_TransitionMatrix_3ac_sumdiagcov
+_pd_periodicity_wang_th0_01 = Catch22._PD_PeriodicityWang_th0_01
+_co_embed2_dist_tau_d_expfit_meandiff = Catch22._CO_Embed2_Dist_tau_d_expfit_meandiff
+_in_auto_mutual_info_stats_40_gaussian_fmmi = (
+    Catch22._IN_AutoMutualInfoStats_40_gaussian_fmmi
+)
+_dn_outlier_include_n_001_mdrmd = Catch22._DN_OutlierInclude_n_001_mdrmd
+_sb_binary_stats_diff_longstretch0 = Catch22._SB_BinaryStats_diff_longstretch0
+_sb_motif_three_quantile_hh = Catch22._SB_MotifThree_quantile_hh
+_sc_fluct_anal_2_rsrangefit_50_1_logi_prop_r1 = (
+    Catch22._SC_FluctAnal_2_rsrangefit_50_1_logi_prop_r1
+)
+_sc_fluct_anal_2_dfa_50_1_2_logi_prop_r1 = (
+    Catch22._SC_FluctAnal_2_dfa_50_1_2_logi_prop_r1
+)
+_fc_local_simple_mean3_stderr = Catch22._FC_LocalSimple_mean3_stderr
+
+
+@njit(fastmath=True, cache=True)
+def _transform_collection_3d(X, ffts, f_idx, transform_features, outlier_norm):
+    # numba version of the per-case loop over _transform_case for equal length
+    # collections. The branching and helper functions used mirror _transform_case
+    # exactly. As the helper functions use fastmath, output can differ from the
+    # per-case path in the last few float digits. The catch24 mean and standard
+    # deviation features (22 and 23) are left as zero here and filled in by the
+    # caller.
+    n_cases, n_channels, n_timepoints = X.shape
+    nf = len(f_idx)
+    c22 = np.zeros((n_cases, nf * n_channels))
+
+    for c in range(n_cases):
+        f_count = -1
+        for i in range(n_channels):
+            series = X[c, i]
+            dim = i * nf
+
+            smin = 0.0
+            smax = 0.0
+            has_min_max = False
+            smean = 0.0
+            has_mean = False
+            outlier_series = series
+            has_outlier_series = False
+            ac = np.zeros(0)
+            has_ac = False
+            acfz = 0
+            has_acfz = False
+
+            for n in range(nf):
+                f_count += 1
+                if not transform_features[f_count]:
+                    continue
+
+                feature = f_idx[n]
+
+                # lazily compute the shared precomputed values used by the
+                # current feature
+                if feature == 0 or feature == 1 or feature == 4:
+                    if not has_min_max:
+                        smin = numba_min(series)
+                        smax = numba_max(series)
+                        has_min_max = True
+                elif feature == 7:
+                    if not has_mean:
+                        smean = mean(series)
+                        has_mean = True
+                elif (feature == 13 or feature == 14) and outlier_norm:
+                    if not has_mean:
+                        smean = mean(series)
+                        has_mean = True
+                    if not has_outlier_series:
+                        outlier_series = z_normalise_series_with_mean(series, smean)
+                        has_outlier_series = True
+                elif feature == 2 or feature == 3:
+                    if not has_ac:
+                        ac = _compute_autocorrelations(series)
+                        has_ac = True
+                elif feature == 12 or feature == 10 or feature == 8:
+                    if not has_ac:
+                        ac = _compute_autocorrelations(series)
+                        has_ac = True
+                    if not has_acfz:
+                        acfz = _ac_first_zero(ac)
+                        has_acfz = True
+
+                if feature == 0:
+                    c22[c, dim + n] = _histogram_mode(series, 5, smin, smax)
+                elif feature == 1:
+                    c22[c, dim + n] = _histogram_mode(series, 10, smin, smax)
+                elif feature == 2:
+                    c22[c, dim + n] = _co_f1ecac(ac, n_timepoints)
+                elif feature == 3:
+                    c22[c, dim + n] = _co_first_min_ac(ac, n_timepoints)
+                elif feature == 4:
+                    c22[c, dim + n] = _co_histogram_ami_even_2_5(series, smin, smax)
+                elif feature == 5:
+                    c22[c, dim + n] = _co_trev_1_num(series)
+                elif feature == 6:
+                    c22[c, dim + n] = _md_hrv_classic_pnn40(series)
+                elif feature == 7:
+                    c22[c, dim + n] = _sb_binary_stats_mean_longstretch1(series, smean)
+                elif feature == 8:
+                    c22[c, dim + n] = _sb_transition_matrix_3ac_sumdiagcov(series, acfz)
+                elif feature == 9:
+                    c22[c, dim + n] = _pd_periodicity_wang_th0_01(series)
+                elif feature == 10:
+                    c22[c, dim + n] = _co_embed2_dist_tau_d_expfit_meandiff(
+                        series, acfz
+                    )
+                elif feature == 11:
+                    c22[c, dim + n] = _in_auto_mutual_info_stats_40_gaussian_fmmi(
+                        series
+                    )
+                elif feature == 12:
+                    # mirrors _FC_LocalSimple_mean1_tauresrat, which is not a
+                    # numba function
+                    if n_timepoints < 2:
+                        c22[c, dim + n] = 0
+                    else:
+                        res = _local_simple_mean(series, 1)
+                        res_ac = _compute_autocorrelations(res)
+                        c22[c, dim + n] = _ac_first_zero(res_ac) / acfz
+                elif feature == 13:
+                    c22[c, dim + n] = _outlier_include(
+                        outlier_series if outlier_norm else series
+                    )
+                elif feature == 14:
+                    c22[c, dim + n] = _dn_outlier_include_n_001_mdrmd(
+                        outlier_series if outlier_norm else series
+                    )
+                elif feature == 15:
+                    c22[c, dim + n] = _summaries_welch_rect(series, False, ffts[c, i])
+                elif feature == 16:
+                    c22[c, dim + n] = _sb_binary_stats_diff_longstretch0(series)
+                elif feature == 17:
+                    c22[c, dim + n] = _sb_motif_three_quantile_hh(series)
+                elif feature == 18:
+                    c22[c, dim + n] = _sc_fluct_anal_2_rsrangefit_50_1_logi_prop_r1(
+                        series
+                    )
+                elif feature == 19:
+                    c22[c, dim + n] = _sc_fluct_anal_2_dfa_50_1_2_logi_prop_r1(series)
+                elif feature == 20:
+                    c22[c, dim + n] = _summaries_welch_rect(series, True, ffts[c, i])
+                elif feature == 21:
+                    c22[c, dim + n] = _fc_local_simple_mean3_stderr(series)
+
+    return c22
diff --git a/aeon/transformations/collection/interval_based/_random_intervals.py b/aeon/transformations/collection/interval_based/_random_intervals.py
index f85538b890..f214f45bbb 100644
--- a/aeon/transformations/collection/interval_based/_random_intervals.py
+++ b/aeon/transformations/collection/interval_based/_random_intervals.py
@@ -139,19 +139,34 @@ def __init__(
 
     transformer_feature_skip = ["transform_features_", "_transform_features"]
 
+    def _eval_intervals(self, tasks):
+        """Run a list of joblib ``delayed`` interval tasks.
+
+        Runs sequentially when ``n_jobs == 1`` to skip the joblib dispatch
+        overhead (this transform is invoked once per representation per tree
+        inside the interval forests, always with n_jobs=1), otherwise in
+        parallel. Tasks are pre-built as a list so any random draws in their
+        arguments occur in the same order joblib would consume a generator.
+        """
+        if self._n_jobs == 1:
+            return [func(*args, **kwargs) for func, args, kwargs in tasks]
+        return Parallel(
+            n_jobs=self._n_jobs, backend=self.parallel_backend, prefer="threads"
+        )(tasks)
+
     def _fit_transform(self, X, y=None):
         X, rng = self._fit_setup(X)
 
-        fit = Parallel(
-            n_jobs=self._n_jobs, backend=self.parallel_backend, prefer="threads"
-        )(
-            delayed(self._generate_interval)(
-                X,
-                y,
-                rng.randint(np.iinfo(np.int32).max),
-                True,
-            )
-            for _ in range(self.n_intervals)
+        fit = self._eval_intervals(
+            [
+                delayed(self._generate_interval)(
+                    X,
+                    y,
+                    rng.randint(np.iinfo(np.int32).max),
+                    True,
+                )
+                for _ in range(self.n_intervals)
+            ]
         )
 
         (
@@ -160,7 +175,7 @@ def _fit_transform(self, X, y=None):
         ) = zip(*fit)
 
         current = []
-        removed_idx = []
+        kept_parts = []
         self.n_intervals_ = 0
         for i, interval in enumerate(intervals):
             new_interval = (
@@ -173,29 +188,24 @@ def _fit_transform(self, X, y=None):
                 current.append(new_interval)
                 self.intervals_.extend(interval)
                 self.n_intervals_ += 1
-            else:
-                removed_idx.append(i)
+                kept_parts.append(transformed_intervals[i])
 
-        Xt = transformed_intervals[0]
-        for i in range(1, self.n_intervals):
-            if i not in removed_idx:
-                Xt = np.hstack((Xt, transformed_intervals[i]))
-
-        return Xt
+        # Concatenate once rather than growing Xt with a fresh copy per interval.
+        return np.hstack(kept_parts)
 
     def _fit(self, X, y=None):
         X, rng = self._fit_setup(X)
 
-        fit = Parallel(
-            n_jobs=self._n_jobs, backend=self.parallel_backend, prefer="threads"
-        )(
-            delayed(self._generate_interval)(
-                X,
-                y,
-                rng.randint(np.iinfo(np.int32).max),
-                False,
-            )
-            for _ in range(self.n_intervals)
+        fit = self._eval_intervals(
+            [
+                delayed(self._generate_interval)(
+                    X,
+                    y,
+                    rng.randint(np.iinfo(np.int32).max),
+                    False,
+                )
+                for _ in range(self.n_intervals)
+            ]
         )
 
         (
@@ -232,22 +242,19 @@ def _transform(self, X, y=None):
                         transform_features.append(self._transform_features[count])
                         count += 1
 
-        transform = Parallel(
-            n_jobs=self._n_jobs, backend=self.parallel_backend, prefer="threads"
-        )(
-            delayed(self._transform_interval)(
-                X,
-                i,
-                transform_features[i],
-            )
-            for i in range(len(self.intervals_))
+        transform = self._eval_intervals(
+            [
+                delayed(self._transform_interval)(
+                    X,
+                    i,
+                    transform_features[i],
+                )
+                for i in range(len(self.intervals_))
+            ]
         )
 
-        Xt = transform[0]
-        for i in range(1, len(self.intervals_)):
-            Xt = np.hstack((Xt, transform[i]))
-
-        return Xt
+        # Concatenate once rather than growing Xt with a fresh copy per interval.
+        return np.hstack(transform)
 
     def _fit_setup(self, X):
         self.intervals_ = []
@@ -381,10 +388,9 @@ def _generate_interval(self, X, y, seed, transform):
                         y,
                     )
             elif transform:
-                t = [
-                    [f]
-                    for f in feature(X[:, dim, interval_start:interval_end:dilation])
-                ]
+                t = np.asarray(
+                    feature(X[:, dim, interval_start:interval_end:dilation])
+                ).reshape(-1, 1)
                 Xt = np.hstack((Xt, t))
 
             intervals.append((interval_start, interval_end, dim, feature, dilation))
@@ -401,7 +407,7 @@ def _transform_interval(self, X, idx, keep_transform):
                         setattr(feature, n, keep_transform)
                         break
             elif not keep_transform:
-                return [[0] for _ in range(X.shape[0])]
+                return np.zeros((X.shape[0], 1))
 
         if isinstance(feature, BaseTransformer):
             Xt = feature.transform(
@@ -411,7 +417,9 @@ def _transform_interval(self, X, idx, keep_transform):
             if Xt.ndim == 3:
                 Xt = Xt.reshape((Xt.shape[0], Xt.shape[2]))
         else:
-            Xt = [[f] for f in feature(X[:, dim, interval_start:interval_end:dilation])]
+            Xt = np.asarray(
+                feature(X[:, dim, interval_start:interval_end:dilation])
+            ).reshape(-1, 1)
 
         return Xt