random forest, ensemble-size selection, plateau search, stationary distribution, birth–death Markov chain

1 Introduction↩︎

The setting is supervised learning on tabular data with \(n\) observations and \(p\) input features. After standard preprocessing of categorical variables, such data can be represented by a design matrix \(X\in\mathbb{R}^{n\times p}\), whose rows correspond to observations and columns to features, together with a target vector \(y\). Tree-based ensembles remain state-of-the-art for tabular data, consistently outperforming deep learning in many benchmark studies [1]–[3]. Among them, Gradient Boosting [4], [5] often delivers superior predictive accuracy, while Random Forest [6], [7] offers greater stability. This stability stems from averaging over many randomized trees: the Monte Carlo component of the ensemble prediction decreases as the number of trees grows. The independence of trees also enables parallel training and the use of out-of-bag (OOB) scores, which provide an internal performance estimate without cross-validation.

Beyond prediction, Random Forests provide variable importance measures (VIMs), such as the Mean Decrease in Impurity (MDI), which naturally arise from the tree-building procedure [6], [8]. These VIMs are widely used for feature selection, network inference, and scientific discovery [9]–[11]. However, stabilizing VIMs in high-dimensional settings with correlated features may require substantially more trees than stabilizing the predictive score itself [12], [13]. A well-calibrated predictive score is therefore a necessary, though not sufficient, condition for trustworthy VIMs. This observation reinforces the motivation of the present line of work: before one can reliably assess variable importance, one must first obtain a Random Forest with stable and sufficiently accurate predictive performance.

The number of trees \(T\) is therefore a central computational parameter of Random Forests. Increasing \(T\) reduces finite-ensemble variability, but it also increases training and prediction cost. Standard hyperparameter optimization (HPO) methods, such as TPE [14] or Hyperband [15], require the user to specify a search range \([T_{\min},T_{\max}]\). Because adding trees does not induce the usual overfitting behavior of many other hyperparameters, the selected value of \(T\) tends to be driven toward the upper boundary \(T_{\max}\). Raising \(T_{\max}\) shifts the selected value further toward the boundary, offering no guarantee that the chosen bound is either sufficient or computationally efficient. Early-stopping heuristics avoid an explicit upper bound by monitoring incremental score improvements, but they may stop prematurely when OOB score fluctuations make the observed improvement appear small.

In a recent paper, [16] introduced a triplet-based plateau search procedure that adapts the number of trees without requiring a fixed \(T_{\max}\). At each HPO trial, the non-\(T\) hyperparameters are sampled in the usual way, while the ensemble size is represented by a geometric triplet \(L=B/\mathrm{sf}\), \(B\), and \(R=B\cdot\mathrm{sf}\), with a fixed scale factor \(\mathrm{sf}>1\). The geometric construction keeps the relative separation between neighboring ensemble sizes approximately constant, since \(R-B=(\mathrm{sf}-1)B\) and hence \((R-B)/B=\mathrm{sf}-1\). This is important because, for a fixed additive step \(R=B+\Delta\), the relative separation satisfies \((R-B)/B=\Delta/B\to0\) as \(B\to\infty\), making neighboring forests increasingly difficult to distinguish by their OOB scores. Fixed increments of this type, for example \(\Delta=10\), were used in earlier ensemble-size studies [12], [17].

Thus, the scale factor \(\mathrm{sf}\) acts as a resolution parameter for the plateau comparison. The forest is trained sequentially up to \(R\) trees, and the OOB scores at \(L\), \(B\), and \(R\) are recorded. The relative score gaps \[\label{eq:rel95diff95intro} d_L = \frac{|S_B-S_L|}{|S_B|}, \qquad d_R = \frac{|S_R-S_B|}{|S_B|}\tag{1}\] indicate how close the central ensemble size \(B\) is to the plateau region. They are compared with a user-specified tolerance \(\varepsilon\), typically of the order \(10^{-3}\).

If both inequalities \(d_L\le\varepsilon\) and \(d_R\le\varepsilon\) hold, the ensemble is considered stable but potentially unnecessarily large, and the triplet is shifted left for the next trial. In the original rule, if \(d_R>\varepsilon\), the improvement from \(B\) to \(R\) is still above the tolerance, and the current ensemble is treated as insufficient; the triplet is therefore shifted right. In the remaining original-rule case, \(d_L>\varepsilon\) and \(d_R\le\varepsilon\), the triplet remains at the current level. A symmetric modified rule is also analyzed below; in this variant, the mixed case \(d_L\le\varepsilon<d_R\) is assigned to staying rather than to a right shift. Importantly, trials that trigger a right shift are not considered later for selecting the trial with the highest score, because they are deemed too unstable to be trusted. This adaptive mechanism eliminates the need for an arbitrary upper bound \(T_{\max}\) and avoids the systematic underestimation bias of one-shot early stopping. Moreover, it jointly optimizes the ensemble size with other hyperparameters such as tree depth and \(m_{\mathrm{try}}\), recognizing that the required number of trees interacts with the remaining hyperparameters [18], [19].

However, this adaptivity changes the nature of the problem. The central triplet point \(B_t\) evolves across HPO trials according to random OOB score comparisons. This randomness has several sources. First, the OOB scores are computed from finite data and therefore inherit sampling variability. Second, even for fixed hyperparameters and fixed data, Random Forest training involves algorithmic randomness from bootstrap sampling and feature subsampling, which produces finite-ensemble fluctuations of the OOB score. In addition, the non-\(T\) hyperparameters are sampled by the HPO procedure, so changes in \(B_t\) are coupled with the stochastic exploration of the remaining search space.

The triplet-based plateau algorithm partially regularizes this randomness. Unlike one-shot early stopping, it updates the central ensemble size across trials, so left and right corrections can compensate over time. Nevertheless, the process remains random. Consequently, the selected number of trees exhibits substantial run-to-run variability across independent optimization runs with different random seeds. This variability limits the reliability of a single-point estimate of \(T\): a practitioner who runs the same HPO procedure twice may obtain noticeably different recommendations for the ensemble size. This variability is not merely a nuisance, but a fundamental property of the adaptive tuning process. The selected number of trees should therefore be viewed not as a deterministic quantity, but as a random variable induced by the optimization dynamics.

Figure 1 illustrates this behavior. After an initial transient phase, the central triplet point \(B_t\) does not converge to a fixed deterministic value, but fluctuates within a problem-dependent range that is naturally interpreted as a stationary regime. The visibly non-negligible amplitude of these fluctuations motivates the variance analysis below, which shows that the stationary standard deviation on the ensemble-size scale can be a sizeable fraction of the stationary mean. The tolerance \(\varepsilon\) controls how small the relative OOB score gaps in 1 must be before the plateau tests are regarded as passed. The equilibrium center of this regime is denoted by \(B_*\). As \(\varepsilon\) decreases, the plateau condition becomes more stringent, so the equilibrium level shifts toward larger ensemble sizes. This shift also increases the absolute size of the fluctuations when measured in the number of trees. At the same time, the OOB score-difference noise in 1 becomes smaller at larger ensemble sizes because the finite-ensemble variance decays with the number of trees.

This paper develops a theoretical description of this stationary regime. The analysis assumes that the non-\(T\) hyperparameters have effectively stabilized and isolates the stochastic evolution of the central ensemble size \(B_t\). Under this reduced description, the plateau procedure induces a birth–death Markov chain on the geometric grid \(B_j=T_0\cdot\mathrm{sf}^j\). The goal is not to construct a practical estimator from a finite HPO trajectory, but to characterize the population-level stationary quantities of this process: its stationary distribution, equilibrium center, and local spread. This separates the theoretical question studied here from the separate practical problem of estimating these quantities from finite trajectories.

The analysis proceeds in three steps. First, the transition probabilities induced by the original plateau update rule and by a symmetric modified variant are derived. These probabilities are expressed through plateau-pass probabilities for the left and right OOB score gaps. Second, the finite-ensemble variance scaling of the signed relative gaps is used to approximate these pass probabilities by a leading centered folded-normal model. This yields explicit equilibrium equations and shows that the stationary center satisfies \(B_*=O(\varepsilon^{-2})\) as \(\varepsilon\downarrow0\). Third, the local shape of the stationary distribution is analyzed. A product-form local-balance argument gives a local Gaussian approximation, and an equivalent Fokker–Planck interpretation yields explicit variance constants for both update rules. After converting from the grid-level scale to the ensemble-size scale, the stationary standard deviation also satisfies \(\sigma_{B,*}=O(\varepsilon^{-2})\), while the corresponding variance is \(O(\varepsilon^{-4})\).

The main contributions of this paper are as follows. First, triplet-based plateau search is formulated as a birth–death Markov chain on a geometric ensemble-size grid. Second, the product-form stationary distribution is derived and conditions under which it is normalizable are established. Third, under a leading folded-normal approximation, explicit equilibrium equations are obtained for the original and modified update rules, yielding the scaling law \(B_*=O(\varepsilon^{-2})\). Fourth, local stationary variance constants are derived, interpreted through a drift–diffusion representation, and used to show that the ensemble-size standard deviation satisfies \(\sigma_{B,*}=O(\varepsilon^{-2})\). This appears to be the first work to explicitly model the Random Forest ensemble size selected by an adaptive HPO procedure as a stochastic process and to characterize its stationary behavior. Together, these results provide a theoretical basis for interpreting plateau-based Random Forest tuning as a stationary stochastic process rather than as a deterministic stopping rule.

The present paper is an analytical study. It characterizes the population-level stationary quantities of the plateau process and does not attempt to construct a finite-trajectory estimator. Empirical calibration, synthetic Markov-chain simulations, and moment-based estimators from finite HPO trajectories are natural follow-up problems and are discussed in Section 5. The remainder of the paper is organized as follows. Section 2 reviews related work on selecting the number of trees, plateau and stopping rules, and noisy hyperparameter optimization. Section 3 reformulates triplet-based plateau search as a stochastic update process and derives the signed-gap variance asymptotics used later. Section 4 develops the birth–death Markov model, characterizes its stationary distribution, and derives the equilibrium and variance scaling laws. Section 5 discusses the implications and limitations of the theory, including the role of the geometric grid and directions for empirical follow-up. Appendix 6 contains the proofs.

2 Related Work↩︎

2.1 Selecting the Number of Trees in Random Forest↩︎

The challenge of selecting the number of trees \(T\) has been studied extensively. Early empirical work by [20] and [21] confirmed that the predictive score plateaus as \(T\) increases, with [20] relating the required \(T\) to data set density \(\log_p n\). [22] and [23] proposed heuristic mappings for \(T\) based on information-theoretic measures and data set size, respectively. Statistical stopping criteria were developed by [17], who used a McNemar test, and by [24], who derived an asymptotic estimate for the ensemble size needed to match an infinite forest. [25] analyzed the algorithmic variance of randomized ensembles, showing that the conditional error variance decays as \(O(T^{-1})\) and discussing relative stopping conditions. More recently, [12] modeled the relationship between \(T\) and the stability of variable importance measures using a logistic function, extrapolating to very large ensemble sizes.

These works provide valuable guidelines for determining a sufficient number of trees. Most of them, however, treat \(T\) as a deterministic target for a fixed learning problem or for a fixed set of hyperparameters. In contrast, the plateau search procedure studied here adapts \(T\) across HPO trials, and the central ensemble size evolves as a stochastic process driven by OOB score comparisons. The focus here is therefore not a single stopping time or a single sufficient value of \(T\), but the stationary distribution induced by the adaptive triplet update.

2.2 Using the HPO Trace↩︎

A growing body of research argues that discarding all but the single best HPO trial is wasteful. [26] introduced fANOVA, a post-hoc method that uses the entire history of HPO runs to decompose the variance of the observed performance, identifying which hyperparameters are important and how they interact. [27] demonstrated that constructing an ensemble from a library of models generated during a single HPO run can outperform selecting the single best model. Both approaches exploit information contained in the optimization trace rather than relying only on the best observed configuration.

The goal of the present work is different. The present work does not aggregate several fitted models, nor does it estimate global hyperparameter importances. Instead, the trajectory of a single hyperparameter, the number of trees, is used as the object of analysis. The central point \(B_t\) of the plateau triplet is treated as a stochastic process, and the HPO trace motivates studying its stationary center and spread. This view is especially relevant for Random Forests because the effect of \(T\) is primarily computational and variance-reducing rather than a conventional overfitting tradeoff.

2.3 Noisy and Robust Hyperparameter Optimization↩︎

Noisy observations are a central concern in Bayesian optimization (BO). When function evaluations are noisy, standard BO algorithms can be misled by random fluctuations. [28] proposed a modification of the expected improvement acquisition function that integrates over the posterior distribution of noisy observations. Such methods are designed to improve decision making in the presence of observation noise, often by modifying the acquisition function or the identification step.

In the Random Forest plateau setting, the noise has a specific structure. As shown in [16], the variance of relative OOB score differences 1 decays as \(O(1/T)\), or equivalently their standard deviation scales as \(O(T^{-1/2})\). Thus, the uncertainty is heteroscedastic and directly tied to the ensemble size. A generic Gaussian process (GP) surrogate over \(T\) would need to account for this structure, while also separating true saturation of the forest from finite-ensemble OOB fluctuations. More importantly, the objective of noisy BO methods is typically to identify a single best configuration under uncertainty. The objective here is instead to characterize the distribution of the ensemble size produced by an adaptive plateau process.

Related distinctions arise in robust and risk-seeking Bayesian optimization. [29] formulated a risk-seeking BO problem, aiming to find the single best possible reward under uncontrollable environmental noise. For a monotone Random Forest score, such a risk-seeking strategy would tend to favor larger values of \(T\) and does not by itself define a sufficient ensemble size. [30] considered robust optimization under uncertainty by building a GP surrogate for an integral of the objective over uncontrollable variables. These methods address robustness or risk preferences in optimization, whereas the present paper studies the stationary stochastic dynamics generated by a specific adaptive tuning rule.

2.4 Position of the Present Work↩︎

The present paper is closest in spirit to work on stopping rules and sufficient ensemble size, but it changes the object of analysis. Rather than asking when a Random Forest should stop growing, the question considered here is where the balance between upward and downward ensemble-size moves is established on the geometric grid, and how broadly the process fluctuates around this balance point. This perspective is motivated by the empirical observation that the central triplet point need not settle at a single deterministic value. It may instead fluctuate around a problem-dependent region after an initial transient phase.

The contribution is theoretical. The plateau update is formulated as a birth–death Markov chain, its product-form stationary distribution is derived, and conditions for existence and uniqueness of the stationary distribution are established. Under a leading folded-normal approximation for the OOB score-gap probabilities, explicit balance equations are obtained for the original update rule and for a symmetric modified variant. These equations yield the scaling law \(B_*=O(\varepsilon^{-2})\) for the stationary center. Local variance constants are further derived, showing that the stationary standard deviation on the ensemble-size scale also satisfies \(\sigma_{B,*}=O(\varepsilon^{-2})\).

Thus, for a user-specified tolerance such as \(\varepsilon=10^{-3}\), decreasing the tolerance increases both the stationary mean ensemble size and the absolute magnitude of its fluctuations. At the same time, their leading ratio remains controlled by the scale factor and the update rule, rather than by \(\varepsilon\). For commonly used scale factors, this relative spread can be substantial: for example, it is approximately one half of the stationary mean for \(\mathrm{sf}=1.5\). This shows that the variability of the selected number of trees is not only a finite-sample nuisance, but a structural feature of the stationary plateau process.

This distinguishes the paper from post-hoc heuristics, robust BO methods, and model-aggregation approaches. Those methods aim to improve selection, prediction, or robustness in an HPO workflow. Here, the aim is to understand the stochastic process induced by the plateau rule itself. The derived relation between the stationary center and the stationary spread can also serve as a theoretical basis for future trajectory-based estimators of the mean number of trees, including moment-based approaches such as the generalized method of moments. The construction and empirical benchmarking of such estimators are left for separate work.

3 Plateau Search as a Stochastic Process↩︎

The triplet-based plateau search was introduced in [16] as an adaptive mechanism for selecting the number of trees without specifying an explicit upper bound \(T_{\max}\). The full algorithmic description was summarized in the Introduction. This section extracts the part of the construction that is needed for the stochastic analysis: after the non-\(T\) hyperparameters have effectively stabilized, the central triplet point \(B_t\) evolves across HPO trials according to random threshold comparisons of OOB scores.

3.1 Stochastic Plateau Updates↩︎

The full plateau-based HPO procedure jointly samples the non-\(T\) hyperparameters by TPE and updates the ensemble size through a triplet rule. Since the present analysis focuses on the stationary regime after the non-\(T\) hyperparameters have effectively stabilized, only the ensemble-size update component is isolated. At trial \(t\), given the current central point \(B_t\), the triplet points are \(L_t=B_t/\mathrm{sf}\), \(B_t\), and \(R_t=B_t\cdot\mathrm{sf}\). For the sampled non-\(T\) hyperparameter configuration, a nested forest is trained up to \(R_t\) trees, and the OOB scores \(S_{L_t}\), \(S_{B_t}\), and \(S_{R_t}\) are recorded.

The plateau rule is based on the two absolute relative gaps \[d_{L,t} = \left| \frac{S_{B_t}-S_{L_t}}{S_{B_t}} \right|, \qquad d_{R,t} = \left| \frac{S_{R_t}-S_{B_t}}{S_{B_t}} \right|.\] The next central point is then updated according to the four threshold cases \[B_{t+1} = \begin{cases} B_t\cdot\mathrm{sf}, & d_{L,t}>\varepsilon,\quad d_{R,t}>\varepsilon, \\[1mm] B_t, & d_{L,t}>\varepsilon,\quad d_{R,t}\le\varepsilon, \\[1mm] B_t/\mathrm{sf}, & d_{L,t}\le\varepsilon,\quad d_{R,t}\le\varepsilon, \\[1mm] B_t\cdot\mathrm{sf}\;\text{in the original rule, or } B_t\;\text{in the modified rule}, & d_{L,t}\le\varepsilon,\quad d_{R,t}>\varepsilon . \end{cases}\] Thus, the original algorithm and the modified variant considered below differ only in the mixed case where the left plateau test passes while the right plateau test fails. The modified rule assigns this case to a stay decision, which will lead to a more symmetric transition structure in the Markov model.

For the analysis below, it is useful to view the absolute gaps in 1 as absolute values of their signed counterparts, \[\label{eq:signed95relative95gaps} \frac{S_{B_t}-S_{L_t}}{S_{B_t}}, \qquad \frac{S_{R_t}-S_{B_t}}{S_{B_t}} .\tag{2}\] The first quantity measures the local score change from the left neighbor \(L_t\) to the central point \(B_t\), while the second measures the change from \(B_t\) to the right neighbor \(R_t\).

Because the OOB scores are random, the quantities \(d_{L,t}\) and \(d_{R,t}\) are random even when the data and the non-\(T\) hyperparameters are fixed. Consequently, the update of \(B_t\) is stochastic. The unambiguous cases correspond to shifting right, staying, or shifting left according to the plateau logic. The remaining mixed case can be assigned either to a right shift or to a stay decision. Both variants are handled in the Markov formulation below; the distinction affects only the transition probabilities, not the asymptotic score model used to derive them.

3.2 Signed-Gap Variance Asymptotics↩︎

The finite-ensemble asymptotic model developed in [16] is used. Conditionally on the training data \(D\), let \(\mu_T=\mathbb{E}[S_T\mid D]\) denote the expected OOB score of a Random Forest with \(T\) trees. The conditional mean score is assumed to converge to an infinite-forest limit according to \[\label{eq:score95convergence} \mu_T = S_\infty + cT^{-\gamma} + o(T^{-\gamma}), \qquad T\to\infty,\tag{3}\] where \(S_\infty\ne 0\), \(c\ne 0\), and \(\gamma>0\). The stronger condition \(\gamma>1/2\) will be used below when deriving the centered folded-normal approximation for plateau-pass probabilities.

The variance calculation follows the same finite-ensemble Gaussian approximation used in [16]. For the variance asymptotics themselves, however, only its second-order consequences are needed: the \(O(T^{-1})\) decay of algorithmic variance and the covariance scaling for nested warm-start forests, formalized in 4 . The variance component in 4 is motivated by finite-ensemble variance results for randomized ensembles, in particular the \(O(T^{-1})\) decay analyzed by [25]. The covariance component in 4 is an additional nested-forest approximation tailored to the warm-start construction used by the plateau algorithm: when \(T_1<T_2\), the larger forest contains the smaller forest as a prefix, leading to the leading-order covariance \(v/T_2\).

The resulting right-gap variance in 6 was derived in [16]. The corresponding left-gap variance in 5 is derived here by the same nested-forest covariance argument; the proof is given in Appendix 6. Both expressions are needed because the Markov model below uses both the left and right plateau tests. The factor difference between 5 and 6 reflects only the fact that the left comparison involves the smaller ensemble \(L=B/\mathrm{sf}\), whose finite-forest variance is larger. It does not assume any scaling relation between the scores \(S_L\), \(S_B\), and \(S_R\) themselves. In particular, \((\mathrm{sf}-1)/(1-\mathrm{sf}^{-1})=\mathrm{sf}\), so the left signed relative gap has asymptotically \(\mathrm{sf}\) times larger variance than the right signed relative gap under the nested covariance approximation.

3.3 Absolute Gaps and Plateau-Pass Probabilities↩︎

The actual plateau rule uses the absolute values of the signed relative gaps in 2 . For a generic central value \(B\), define the left and right plateau-pass probabilities by \[\label{eq:plateau95pass95probabilities95B} \alpha_L(B;\varepsilon) = \mathbb{P}\left[ \left| \frac{S_B-S_{B/\mathrm{sf}}}{S_B} \right| \le \varepsilon \;\middle|\;D \right], \qquad \alpha_R(B;\varepsilon) = \mathbb{P}\left[ \left| \frac{S_{B\cdot\mathrm{sf}}-S_B}{S_B} \right| \le \varepsilon \;\middle|\;D \right].\tag{7}\]

Figure 2 gives a schematic interpretation of \(\alpha_L(B;\varepsilon)\) and \(\alpha_R(B;\varepsilon)\). The shaded regions correspond to the events that the corresponding signed relative gaps fall inside the tolerance interval \([-\varepsilon,\varepsilon]\). This notation is intentionally distribution-agnostic: the probabilities in 7 may be evaluated using a folded-normal approximation, an empirical or bootstrap approximation, or another model for the OOB score fluctuations.

To obtain explicit analytical formulas, a stronger approximation than the variance calculation above is now imposed. Proposition [prop:signed95gap95variance95asymptotics] uses only the second-order asymptotics of the signed gaps in 2 , whereas the next result assumes an approximate Gaussian law for these transformed signed relative gaps themselves.

The proof is given in Appendix 6. Substituting the variance scales from 5 –6 into 8 –9 gives the leading centered folded-normal approximations \[\label{eq:alpha95leading} \alpha_L(B;\varepsilon) \approx 2\Phi\left( \varepsilon \sqrt{ \frac{S_\infty^2 B}{v(\mathrm{sf}-1)} } \right)-1,\qquad \alpha_R(B;\varepsilon) \approx 2\Phi\left( \varepsilon \sqrt{ \frac{S_\infty^2 B}{v(1-\mathrm{sf}^{-1})} } \right)-1 .\tag{10}\]

For a signed gap \(X_B\) with conditional mean \(m(B)\) and standard deviation \(s(B)\), the non-centered folded-normal approximation is \[\label{eq:folded95normal95noncentered} \mathbb{P}\left[ |X_B|\le\varepsilon \;\middle|\;D \right] \approx \Phi\left(\frac{\varepsilon-m(B)}{s(B)}\right) - \Phi\left(\frac{-\varepsilon-m(B)}{s(B)}\right).\tag{11}\] The asymptotics in 8 –9 follow by expanding 11 in powers of the normalized mean \(m(B)/s(B)\). The first-order terms cancel because the two Gaussian CDF terms enter with opposite signs, so the first nonzero correction is quadratic, that is, of order \(O(m^2(B)/s^2(B))\). The behavior \(m(B)=O(B^{-\beta})\), where \(\beta=\min\{\gamma,1\}>1/2\), follows from a second-order delta-method expansion. Together with \(s(B)\asymp B^{-1/2}\), this implies that \(m(B)/s(B)\to0\), and hence the correction due to the nonzero conditional mean vanishes asymptotically. Thus, the approximations in 10 are centered leading-order folded-normal approximations.

Note that the leading terms in 10 do not involve \(\gamma\); they are determined only by the variance scales in 5 –6 . The condition \(\gamma>1/2\) ensures the vanishing of \(m_L(B)/s_L(B)\) and \(m_R(B)/s_R(B)\). Moreover, the approximation does not require any specific limiting behavior of the ratios \(\varepsilon/s_L(B)\) and \(\varepsilon/s_R(B)\). The tolerance \(\varepsilon\) is a user-controlled parameter, and decreasing it shifts the stationary regime toward larger ensembles, where the noise scales \(s_L(B)\) and \(s_R(B)\) are smaller. However, the correction due to the nonzero conditional means is bounded by \(O(B^{1-2\beta})\) uniformly in \(\varepsilon/s(B)\).

The asymptotic forms in 8 –9 , and their explicit leading versions in 10 , provide the probabilistic input for the Markov model. As \(B\) increases, the standard deviations \(s_L(B)\) and \(s_R(B)\) decrease as \(B^{-1/2}\), so the leading plateau-pass probabilities in 10 increase toward one. This creates an inward drift: small ensembles tend to shift right because the plateau tests often fail, whereas sufficiently large ensembles tend to shift left because both tests are likely to pass. In the next section, this intuition is formalized as a birth–death Markov chain and its stationary distribution is analyzed.

4 Markov Model and Stationary Distribution↩︎

The stochastic plateau updates are now converted into a birth–death Markov chain. The central point \(B_t\) is assumed to evolve on the geometric grid \(T_j=T_0 \cdot \mathrm{sf}^j\), \(j=0,1,\ldots\), where \(T_0>0\) is the initial ensemble size. Thus, for an interior level, the event \(B_t=T_j\) means that at trial \(t\) the central triplet point is located at level \(j\) of the fixed geometric grid, with \(L_t=T_{j-1}=B_t/\mathrm{sf}\) and \(R_t=T_{j+1}=B_t\cdot\mathrm{sf}\), up to rounding to integer tree counts. When no ambiguity is possible, the notation \(B_j=T_j\) is also used for the central ensemble size associated with grid level \(j\). Although the implementation may occasionally move to negative grid indices \(j=-1,-2,\ldots\), this can only occur while the rounded tree counts remain positive and distinct; such lower-boundary effects are rare and do not affect the large-\(B\) stationary analysis below.

On this grid, for simplicity, write \(\alpha_{L,j}=\alpha_L(T_j;\varepsilon)\) and \(\alpha_{R,j}=\alpha_R(T_j;\varepsilon)\), where \(\alpha_L\) and \(\alpha_R\) were defined in 7 . Thus, \(\alpha_{L,j}\) is the probability that the left plateau test passes when the central point is \(B_t=T_j\), and \(\alpha_{R,j}\) is the corresponding probability for the right plateau test. When there is no ambiguity, the argument \(\varepsilon\) is suppressed.

4.1 Plateau Cases and Update Variants↩︎

The update rule is determined by the outcomes of the left and right plateau tests. Here “pass” means that the corresponding absolute relative gap is at most \(\varepsilon\), and “fail” means that it exceeds \(\varepsilon\). The four possible cases are shown in Table 1.

The only difference between the two variants is the mixed case in which the left test passes but the right test fails. In the original algorithm of [16], this case is assigned to a right shift. In the modified variant, it is assigned to staying at the current level. The modified rule is more symmetric: a left shift occurs only when both tests pass, a right shift occurs only when both tests fail, and mixed evidence leads to no move.

4.2 Transition Probabilities↩︎

The transition probabilities of the Markov chain are introduced as \[p_j^-=\mathbb{P}(B_{t+1}=T_{j-1}\mid B_t=T_j), \qquad p_j^+=\mathbb{P}(B_{t+1}=T_{j+1}\mid B_t=T_j),\] and \(p_j^0=1-p_j^- - p_j^+\). Here \(p_j^-\) is the probability of a left shift, \(p_j^+\) is the probability of a right shift, and \(p_j^0\) is the probability of staying at the same level.

The exact transition probabilities are determined by the joint distribution of the left and right plateau tests. For explicit analytical formulas, the factorized approximation is used, in which the two tests are treated as conditionally independent at the same level. Under this approximation, the original update rule gives \[\label{eq:transition95original} p_j^- \approx \alpha_{L,j}\alpha_{R,j}, \qquad p_j^+ \approx 1-\alpha_{R,j}, \qquad p_j^0 \approx \alpha_{R,j}(1-\alpha_{L,j}).\tag{12}\] Indeed, a left shift occurs when both tests pass, while a right shift occurs whenever the right test fails.

For the modified update rule, the transition probabilities become \[\label{eq:transition95modified} p_j^- \approx \alpha_{L,j}\alpha_{R,j}, \qquad p_j^+ \approx (1-\alpha_{L,j})(1-\alpha_{R,j}), \qquad p_j^0 \approx \alpha_{R,j}(1-\alpha_{L,j}) + \alpha_{L,j}(1-\alpha_{R,j}).\tag{13}\] In this case, the chain moves right only when both tests fail, moves left only when both tests pass, and remains at the current level in the two mixed cases. If the dependence between the two plateau tests is modeled explicitly, the products in 12 –13 should be replaced by the corresponding joint probabilities. The birth–death analysis below remains unchanged after this replacement.

In the remainder of the paper, the auxiliary birth–death chains defined by the factorized transition probabilities in 12 and 13 are analyzed. Within these auxiliary chains, the displayed transition probabilities are treated as exact once the plateau-pass probabilities \(\alpha_{L,j}\) and \(\alpha_{R,j}\) are specified. In contrast, the leading centered folded-normal relation 10 is an approximation linking these probabilities to the OOB score-gap model. For the equilibrium analysis, this approximation is used to obtain leading asymptotic balance equations. For the local variance analysis, the corresponding auxiliary chain is introduced by treating 10 as exact on the grid \(B_j=T_0\cdot\mathrm{sf}^j\), which allows the local derivatives of the transition probabilities to be evaluated explicitly.

4.3 Stationary Distribution↩︎

The transition probabilities in 12 or 13 define a birth–death Markov chain on the levels \(j=0,1,\ldots\) of the geometric grid. If the limiting distribution exists, it is denoted by \(\pi_j=\lim_{t\to\infty}\mathbb{P}(B_t=T_j)\). Any such limiting distribution must satisfy the stationary Kolmogorov–Chapman equations. The existence and normalizability of this distribution are established in the next subsection.

The lower-boundary convention is \(p_0^-=0\). The stationary Kolmogorov–Chapman equations are \[\begin{align} \pi_0 &= \pi_0p_0^0 + \pi_1p_1^-, \\ \pi_j &= \pi_{j-1}p_{j-1}^+ + \pi_jp_j^0 + \pi_{j+1}p_{j+1}^-, \qquad j > 0. \end{align}\] These equations follow from the law of total probability: to be at level \(j\) after one step, the chain must have shifted right from \(j-1\), stayed at \(j\), or shifted left from \(j+1\). Since \(p_0^0=1-p_0^+\) and \(p_j^0=1-p_j^- -p_j^+\) for \(j\ge1\), these equations are equivalently written as \[\begin{align} \pi_0p_0^+ &= \pi_1p_1^-, \\ \pi_j(p_j^-+p_j^+) &= \pi_{j-1}p_{j-1}^+ + \pi_{j+1}p_{j+1}^-, \qquad j > 0. \end{align}\] Substituting the first equation into the second one for \(j=1\) yields \(\pi_1p_1^+=\pi_2p_2^-\). Repeating the same argument recursively gives the local balance relations \[\label{eq:detailed95balance} \pi_jp_j^+ = \pi_{j+1}p_{j+1}^-, \qquad j=0,1,\ldots.\tag{14}\] Equation 14 gives the recursion \(\pi_{j+1}=\pi_j p_j^+/p_{j+1}^-\), from which the stationary probabilities have the product form \[\label{eq:stationary95product} \pi_j = \pi_0 \prod_{k=0}^{j-1} \frac{p_k^+}{p_{k+1}^-}, \qquad j > 0.\tag{15}\]

4.4 Existence and Uniqueness of a Stationary Distribution↩︎

The proof is given in Appendix 6. Theorem [thm:stationary95existence] verifies the standard product-form normalization condition for birth–death chains. In the notation of 15 , this condition is the convergence of the series \[\sum_{j=1}^{\infty} \prod_{k=0}^{j-1} \frac{p_k^+}{p_{k+1}^-} <\infty .\] This is the discrete birth–death analogue of the classical Karlin–McGregor normalization condition [31]. The proof also uses Proposition [prop:signed95gap95variance95asymptotics] to show that the noise scales \(s_L(T_j)\) and \(s_R(T_j)\) vanish as \(T_j\to\infty\), and Proposition [prop:plateau95pass95probability95asymptotics] to convert this variance decay into \(\alpha_{L,j}\to1\) and \(\alpha_{R,j}\to1\). Consequently, at large ensemble sizes the chain shifts left with probability tending to one and shifts right with probability tending to zero.

Under the folded-normal transition model, the neighboring transition probabilities are strictly positive, so the birth–death chain is irreducible on the geometric grid. Irreducibility implies that the stationary distribution, once it exists, is unique. Moreover, for an irreducible countable-state Markov chain, the existence of a stationary distribution is equivalent to positive recurrence [32]. In practical terms, positive recurrence means that the chain returns to its recurrent states in finite expected time; hence fluctuations away from the typical stationary region are not transient excursions to infinity. Thus, the stationary regime is not only formally normalizable, but also dynamically stable.

4.5 Equilibrium Level↩︎

The stationary distribution is expected to concentrate near the grid levels where the local drift changes sign. Figure 3 illustrates the transition from a realized trajectory of the grid-level process to its stationary distribution. The equilibrium center \(j_*\) need not coincide with an integer grid level; it is interpreted as the center of the local Gaussian approximation to the stationary masses \(\pi_j\).

In the continuous approximation, the equilibrium level \(j_*\) is determined by equality of the right- and left-shift probabilities, \[\label{eq:equilibrium95general} p_{j_*}^+ \approx p_{j_*}^- .\tag{16}\] Equation 16 is the zero-drift condition that identifies the center of the stationary regime in the continuous approximation.

For compactness, write \(F_0(y)=2\Phi(y)-1\) and \(y_j=\varepsilon/s_R(T_j)\). It follows from 5 –6 that \(s_L(T_j)\sim\sqrt{\mathrm{sf}}\,s_R(T_j)\). Therefore, using the leading centered folded-normal approximation 10 , \[\alpha_{L,j} \approx F_0\left(\frac{y_j}{\sqrt{\mathrm{sf}}}\right), \qquad \alpha_{R,j} \approx F_0(y_j).\]

For the original update rule, substituting 12 into 16 gives \(1-\alpha_{R,j_*}=\alpha_{L,j_*}\alpha_{R,j_*}\). Hence \(y_*^{\mathrm{orig}}=y_{j_*}\) satisfies \[\label{eq:balance95original} F_0(y_*^{\mathrm{orig}}) \left[ 1+ F_0\left(\frac{y_*^{\mathrm{orig}}}{\sqrt{\mathrm{sf}}}\right) \right] = 1 .\tag{17}\] For the modified update rule, substituting 13 into 16 gives \[\left[ 1-F_0\left(\frac{y_*^{\mathrm{mod}}}{\sqrt{\mathrm{sf}}}\right) \right] \left[ 1-F_0(y_*^{\mathrm{mod}}) \right] = F_0\left(\frac{y_*^{\mathrm{mod}}}{\sqrt{\mathrm{sf}}}\right) F_0(y_*^{\mathrm{mod}}).\] Equivalently, this reduces to the symmetric balance equation \[\label{eq:balance95modified95symmetric} F_0\left(\frac{y_*^{\mathrm{mod}}}{\sqrt{\mathrm{sf}}}\right) + F_0(y_*^{\mathrm{mod}}) = 1 .\tag{18}\] This simplification is a direct consequence of assigning the mixed case, in which one plateau test passes and the other fails, to the stay decision. Since \(F_0\) is strictly increasing from \(0\) to \(1\), both 17 and 18 have unique positive solutions for each fixed \(\mathrm{sf}>1\).

The balance equation 18 has a simple interpretation. Since mixed outcomes are assigned to the stay decision, only the two joint events “both tests fail” and “both tests pass” contribute to the drift. Thus, equilibrium requires these two probabilities to be equal: \((1-\alpha_L)(1-\alpha_R)=\alpha_L\alpha_R\), which reduces to \(\alpha_L+\alpha_R=1\). Because the left gap has the larger variance, typically \(\alpha_L<\alpha_R\). Together with \(\alpha_L+\alpha_R=1\), this gives \(\alpha_L<1/2<\alpha_R\).

Once \(y_*^{\mathrm{orig}}\) or \(y_*^{\mathrm{mod}}\) is found from 17 or 18 , respectively, the corresponding equilibrium ensemble size follows from \(y_*=\varepsilon/s_R(B_*)\), where \(B_*=T_0\cdot\mathrm{sf}^{j_*}\) denotes the tree count at the equilibrium level. Using 6 , this gives \[\label{eq:B95star95variants} B_*^{\mathrm{orig}} \approx \frac{v(1-\mathrm{sf}^{-1})}{S_\infty^2} \frac{\left[y_*^{\mathrm{orig}}(\mathrm{sf})\right]^2}{\varepsilon^2}, \qquad B_*^{\mathrm{mod}} \approx \frac{v(1-\mathrm{sf}^{-1})}{S_\infty^2} \frac{\left[y_*^{\mathrm{mod}}(\mathrm{sf})\right]^2}{\varepsilon^2}.\tag{19}\]

Since the \(y_*\)-constants are independent of \(\varepsilon\), the preceding derivation yields the following scaling result.

Proposition [prop:equilibrium95scaling] has a direct practical interpretation. The tolerance \(\varepsilon\) controls the equilibrium ensemble size quadratically: decreasing \(\varepsilon\) by a factor of two increases the equilibrium number of trees by approximately a factor of four, all else being equal.

Moreover, since the constants \(y_*^{\mathrm{orig}}(\mathrm{sf})\) and \(y_*^{\mathrm{mod}}(\mathrm{sf})\) depend only on the scale factor and the update rule, they can be precomputed once for commonly used values such as \(\mathrm{sf}=1.5\) or \(\mathrm{sf}=2.0\). The formulas 19 still contain the unknown problem-dependent factor \(v/S_\infty^2\). This does not prevent their use: the result identifies the functional dependence on \(\varepsilon\) and \(\mathrm{sf}\) up to a single scale coefficient. This coefficient may later be treated as a nuisance parameter in trajectory-based calibration, but such estimation is outside the scope of the present paper.

4.6 Original Versus Modified Equilibrium↩︎

The two update rules differ only in the mixed case \(d_{L,t}\le\varepsilon\), \(d_{R,t}>\varepsilon\). The original rule assigns this case to a right shift, whereas the modified rule assigns it to staying at the current level. Therefore, relative to the modified rule, the original rule has a stronger preference for right shifts and should be expected to place the stationary regime at a larger ensemble size.

This effect can be quantified directly from 19 . The unknown factor \(v(1-\mathrm{sf}^{-1})/(S_\infty^2\varepsilon^2)\) cancels in the ratio, giving \[\label{eq:orig95mod95ratio} \frac{B_*^{\mathrm{orig}}}{B_*^{\mathrm{mod}}} \approx \left( \frac{y_*^{\mathrm{orig}}(\mathrm{sf})}{y_*^{\mathrm{mod}}(\mathrm{sf})} \right)^2 .\tag{20}\] Thus, the relative inflation of the equilibrium ensemble size caused by the original mixed-case assignment depends only on the scale factor \(\mathrm{sf}\), and not on \(\varepsilon\) or on the problem-dependent factor \(v/S_\infty^2\).

The limiting case \(\mathrm{sf}\to1^+\) is available analytically. Since \(F_0(y/\sqrt{\mathrm{sf}})\to F_0(y)\), the modified balance equation 18 gives \(2F_0(y)=1\), hence \(y_*^{\mathrm{mod}}(1^+)=\Phi^{-1}(3/4)\). For the original rule, 17 gives \(F_0(y)(1+F_0(y))=1\), and therefore \(F_0(y)=(\sqrt{5}-1)/2\). Consequently, \[y_*^{\mathrm{orig}}(1^+) = \Phi^{-1}\left(\frac{1+\sqrt{5}}{4}\right), \qquad y_*^{\mathrm{mod}}(1^+) = \Phi^{-1}\left(\frac{3}{4}\right),\] and taking the ratio gives \[R_0=\lim_{\mathrm{sf}\to1^+} \frac{B_*^{\mathrm{orig}}}{B_*^{\mathrm{mod}}} = \left[ \frac{ \Phi^{-1}\left((1+\sqrt{5})/4\right) }{ \Phi^{-1}(3/4) } \right]^2 \approx 1.68 .\] Moreover, expanding the implicit equations around \(\mathrm{sf}=1\) gives \[\frac{B_*^{\mathrm{orig}}}{B_*^{\mathrm{mod}}} = R_0 \left[ 1-\frac{\mathrm{sf}-1}{2\sqrt{5}} + O\left((\mathrm{sf}-1)^2\right) \right].\]

Thus, near the fine-grid limit \(\mathrm{sf}=1\), the original rule inflates the equilibrium ensemble size by about \(68\%\) relative to the symmetric modified rule, and this inflation decreases as \(\mathrm{sf}\) moves away from one. Figure 4 shows the numerical solution of 17 and 18 over \(\mathrm{sf}\in(1,2.5]\), with the commonly used values \(\mathrm{sf}=1.5\) and \(\mathrm{sf}=2.0\) highlighted. For \(\mathrm{sf}=1.5\), the ratio in 20 is approximately \(1.545\), meaning that the original update rule places the equilibrium ensemble size about \(54.5\%\) higher than the modified rule. For \(\mathrm{sf}=2.0\), the corresponding ratio is approximately \(1.467\), or about \(46.7\%\) higher.

From the perspective of HPO efficiency, this creates a trade-off. The previous experiments in [16] showed that using \(\mathrm{sf}=1.5\) can be statistically significantly faster than the classical doubling strategy of [20], \(\mathrm{sf}=2.0\), because the procedure spends less time training forests corresponding to the expensive right endpoints \(R_t=B_t\cdot\mathrm{sf}\). At the same time, the relative inflation induced by the conservative original rule is larger for \(\mathrm{sf}=1.5\) than for \(\mathrm{sf}=2.0\). Thus, smaller scale factors lower the cost of individual right probes, but the conservative mixed-case assignment raises the equilibrium level relative to the symmetric modified benchmark.

This comparison also clarifies the role of the triplet-based plateau mechanism relative to conventional early-stopping heuristics. Early-stopping methods typically scan ensemble sizes from left to right and stop once a local improvement criterion falls below a tolerance, for example when the left-gap condition \(d_{L,t}\le\varepsilon\) is first satisfied. Such one-sided stopping can lead to systematic underestimation of the sufficient number of trees. In contrast, the plateau algorithm performs stochastic updates around an equilibrium region: erroneous left and right corrections can compensate over trials. The original rule was designed as a more conservative variant by assigning the mixed case \(d_{L,t}\le\varepsilon<d_{R,t}\) to a right shift. The analysis above shows that this conservatism has a measurable cost: relative to the symmetric modified rule, it shifts the stationary ensemble size upward by a factor depending only on \(\mathrm{sf}\).

4.7 Local Gaussian Approximation of the Stationary Distribution↩︎

The product-form representation 15 gives the stationary distribution exactly, up to normalization. However, for interpretation and estimation it is useful to approximate this distribution locally near its equilibrium center. This subsection derives such an approximation directly from the discrete local-balance relation.

Let \(J_t\) denote the grid-level process, so that \(B_t=T_{J_t}=T_0\cdot\mathrm{sf}^{J_t}\). From 14 , the adjacent stationary probabilities satisfy \[\label{eq:local95log95ratio} \log\frac{\pi_{j+1}}{\pi_j} = g_j, \qquad g_j = \log\frac{p_j^+}{p_{j+1}^-}.\tag{21}\] Thus, the local shape of the stationary distribution is controlled by the log-ratio \(g_j\). If \(g_j>0\), the stationary mass increases from level \(j\) to level \(j+1\); if \(g_j<0\), it decreases. Therefore, the center of the stationary region is located near the point where this log-ratio changes sign.

For a local continuous approximation, regard \(g_j\) as the value of a smooth function \(g(u)\) at the midpoint \(u=j+1/2\) between the adjacent grid levels \(j\) and \(j+1\). This midpoint convention is natural because \(\log(\pi_{j+1}/\pi_j)\) compares the masses on the two sides of the midpoint \(j+1/2\). Let \(j_*\) be the local equilibrium center, defined by \(g(j_*)=0\). Around this point, \(\pi_j\) is expected to have a local maximum: to the left of \(j_*\) the mass increases with \(j\), while to the right of \(j_*\) it decreases. Equivalently, \(g(u)\) crosses zero with negative slope. Thus, near \(j_*\), write \(g(u)\approx g'(j_*)(u-j_*)\), where \(g'(j_*)<0\) under inward drift.

Now compare this with a Gaussian approximation on the grid-level scale, \(\pi_j\propto\exp\{-(j-j_*)^2/(2\sigma_{J,*}^2)\}\). For such a Gaussian sequence, \[\log\frac{\pi_{j+1}}{\pi_j} = -\frac{j+1/2-j_*}{\sigma_{J,*}^2}.\] Matching this expression with the linear expansion of \(g(u)\) at \(u=j+1/2\) yields \[\label{eq:sigma95J95local95balance} \sigma_{J,*}^2 \approx -\frac{1}{g'(j_*)}.\tag{22}\] The sign condition \(g'(j_*)<0\) guarantees that the variance in 22 is positive.

The approximation in 22 is discrete in origin. It comes from the exact local-balance relation 21 ; the only approximation is the local linearization of \(g(u)\), or equivalently the quadratic approximation of \(\log\pi_j\) near its maximum. In practice, \(g'(j_*)\) can be evaluated from the transition probabilities 12 or 13 , using the plateau-pass probabilities \(\alpha_{L,j}\) and \(\alpha_{R,j}\). Equivalently, it can be estimated numerically from a local fit of \(g_j\) around the equilibrium region.

This also gives a regression interpretation of the Gaussian approximation. Over a window of grid levels near the stationary center, one may fit \(g_j\approx\beta_0+\beta_1(j+1/2)\). Comparing this with \(g_j\approx-(j+1/2-j_*)/\sigma_{J,*}^2\) gives \[\sigma_{J,*}^2\approx-\frac{1}{\beta_1}, \qquad j_*\approx-\frac{\beta_0}{\beta_1}.\] This is a least-squares approximation to the local-balance equations, rather than a likelihood model for the observed HPO trajectory.

Finally, the grid-level approximation can be translated to the original ensemble-size scale. Since \(B=T_0\cdot\mathrm{sf}^J\), the derivative with respect to the grid level is \(dB/dJ=B\log(\mathrm{sf})\). At the equilibrium point \(B_*=T_0\cdot\mathrm{sf}^{j_*}\), the first-order delta-method conversion gives \[\label{eq:sigma95B95delta} \sigma_{B,*} \approx \left. \frac{d}{dJ} \left(T_0\cdot\mathrm{sf}^{J}\right) \right|_{J=j_*} \sigma_{J,*} \approx B_*\log(\mathrm{sf})\,\sigma_{J,*}.\tag{23}\] Thus, even a moderate spread on the grid-level scale can correspond to a large absolute spread in the number of trees when \(B_*\) is large. This \(B\)-scale fluctuation should not be confused with the relative OOB score-gap fluctuations used in the plateau tests; the latter enter the pass probabilities through the variance asymptotics of Proposition [prop:signed95gap95variance95asymptotics].

4.8 Fokker–Planck Interpretation↩︎

The local variance formula 22 also has a diffusion interpretation. For one step of the grid-level chain, let \(\Delta J_t=J_{t+1}-J_t\). Then \[\mathbb{E}[\Delta J_t\mid J_t=j]=p_j^+-p_j^-, \qquad \mathbb{E}[(\Delta J_t)^2\mid J_t=j]=p_j^++p_j^- .\] Thus, in a continuous approximation, \(a(j)=p_j^+-p_j^-\) plays the role of the drift coefficient, while \(b(j)=p_j^++p_j^-\) plays the role of the one-step second-moment coefficient. The condition \(p_{j_*}^+=p_{j_*}^-\) is the zero-drift condition \(\mathbb{E}[\Delta J_t\mid J_t=j_*]=a(j_*)=0\), which defines the equilibrium point \(j_*\). The conditional variance of one step is \(\operatorname{Var}(\Delta J_t\mid J_t=j)=b(j)-a^2(j)\). At \(j_*\), however, \(a(j_*)=0\), so the second moment \(b(j_*)\) and the conditional variance coincide.

The corresponding stationary Fokker–Planck equation has the formal zero-flux form \[a(j)\pi(j) - \frac{1}{2}\frac{d}{dj}\{b(j)\pi(j)\} = 0.\] Linearizing near the equilibrium point, where \(a(j_*)=0\), gives \(a(j)\approx a'(j_*)(j-j_*)\) and \(b(j)\approx b(j_*)\). This is a local Ornstein–Uhlenbeck-type approximation: the restoring drift is linearized, while the local step activity is treated as constant over the neighborhood in which the stationary mass is concentrated. Substituting these approximations into the zero-flux equation yields \[a'(j_*)(j-j_*)\pi(j) - \frac{b(j_*)}{2}\pi'(j) = 0.\] Hence \(\pi'(j)/\pi(j)=2a'(j_*)(j-j_*)/b(j_*)\), and integration gives \[\pi(j) \propto \exp\left\{ \frac{a'(j_*)}{b(j_*)}(j-j_*)^2 \right\}.\] Since inward drift implies \(a'(j_*)<0\), this is a Gaussian density with local variance \[\label{eq:fokker95planck95variance} \sigma_{J,*}^2 \approx \frac{b(j_*)}{-2a'(j_*)}.\tag{24}\]

To relate this expression to 22 , consider a local continuous analogue of 21 . Instead of the discrete neighboring ratio \(p_j^+/p_{j+1}^-\), write \(h(j)=\log[p^+(j)/p^-(j)]\), where \(p^+(j)\) and \(p^-(j)\) are smooth local versions of the right- and left-shift probabilities. At the equilibrium point \(j_*\), the zero-drift condition gives \(p^+(j_*)=p^-(j_*)=:p_*\). Therefore, \[h'(j_*) = \frac{(p^+)'(j_*)}{p^+(j_*)} - \frac{(p^-)'(j_*)}{p^-(j_*)} = \frac{(p^+)'(j_*)-(p^-)'(j_*)}{p_*} = \frac{a'(j_*)}{p_*}.\] Since \(b(j_*)=p^+(j_*)+p^-(j_*)=2p_*\), it follows that \(h'(j_*)=2a'(j_*)/b(j_*)\). Thus the local-balance variance \(\sigma_{J,*}^2\approx-1/h'(j_*)\) reduces to the Fokker–Planck expression 24 . The only difference is that the exact discrete local-balance relation in 21 uses the neighboring denominator \(p_{j+1}^-\), whereas the Fokker–Planck argument replaces it by the local smooth value \(p^-(j)\).

The Fokker–Planck argument is therefore best viewed as an interpretation of the discrete result rather than as its replacement. It introduces an additional continuum approximation, whereas 22 is obtained directly from the product-form local balance of the birth–death chain. In the next subsection, 24 is used rather than 22 , since it separates the local drift term \(a(j)\) from the one-step fluctuation term \(b(j)\).

4.9 Evaluating the Stationary Variance↩︎

The local Gaussian approximation in 22 reduces the problem of estimating the stationary spread to the evaluation of the local slope \(g'(j_*)\). This slope can be evaluated analytically under the leading folded-normal transition model, or estimated numerically from the local behavior of the log-ratio \(g_j\) around the stationary region.

The exact discrete expression in 21 uses the neighboring denominator \(p_{j+1}^-\). For the local calculation below, this neighboring term is replaced by its local smooth value, and \(h(j)=\log[p^+(j)/p^-(j)]\) is used. This replacement is consistent with the zero-drift approximation \(p^+(j_*)\approx p^-(j_*)\). Evaluating \(h'(j)\) requires differentiating the pass probabilities with respect to the grid level. To make this operation well defined, the auxiliary birth–death chain is considered in which the leading approximation 10 is treated as an equality, rather than as an asymptotic approximation. For this auxiliary chain, with \(B_j=T_0\cdot\mathrm{sf}^j\), the transition probabilities are smooth functions of the continuous variable \(j\).

Under this leading model, \(y_j=\varepsilon/s_R(B_j)\) has the form \[y_j = C\cdot\mathrm{sf}^{j/2}, \qquad C = \varepsilon \sqrt{\frac{ S_\infty^2 T_0}{v(1-\mathrm{sf}^{-1})}}.\] Hence \(dy_j/dj=\log(\mathrm{sf})\,y_j/2\), so the constant \(C\), and therefore the factors \(\varepsilon\), \(S_\infty\), \(v\), and \(T_0\), drop out of the local slope calculation. Consequently, the stationary variance constants below depend only on \(\mathrm{sf}\) and on the update rule.

The proof is given in Appendix 6. The result uses the transition probabilities 12 and 13 , together with the leading plateau-pass approximation 10 . The proof differentiates the drift \(a(j)=p^+(j)-p^-(j)\), uses \(dy_j/dj=\log(\mathrm{sf})\,y_j/2\), and then substitutes the balance equations 17 and 18 into the Fokker–Planck variance formula 24 . Proposition [prop:stationary95variance95constants] shows that, in the leading model, the grid-level stationary variance constants in 25 and 26 are independent of the tolerance \(\varepsilon\) and depend only on \(\mathrm{sf}\) and on the update rule. Therefore, similarly to the equilibrium constants in 19 , they can be precomputed for commonly used scale factors. Changing \(\varepsilon\) shifts the equilibrium location but does not change the local variance constant on the \(J\)-scale.

For the modified rule, 26 has a useful physical interpretation. It is derived from the Fokker–Planck variance formula 24 . At the modified equilibrium, 18 gives \(\alpha_L(j_*)+\alpha_R(j_*)=1\). Together with the modified transition probabilities in 13 , this implies \[p_* = p_{j_*}^+ = p_{j_*}^- = \alpha_L(j_*)\alpha_R(j_*) = [1-\alpha_L(j_*)][1-\alpha_R(j_*)].\] Thus \(p_*\) is the common probability of a non-stay move to the right or to the left at equilibrium. Moreover, 13 gives \(a(j)=p^+(j)-p^-(j)=1-\alpha_L(j)-\alpha_R(j)\). Therefore, \(-a'(j_*)=\left.d(\alpha_L+\alpha_R)/dj\right|_{j=j_*}\). Since the one-step second-moment coefficient at equilibrium is \(b(j_*)=p_{j_*}^+ + p_{j_*}^- = 2p_*\), the Fokker–Planck expression 24 becomes \[\sigma_{J,*}^{2,\mathrm{mod}} \approx \frac{p_*}{ \left. \dfrac{d}{dj}\bigl(\alpha_L(j)+\alpha_R(j)\bigr) \right|_{j=j_*} }.\] Hence the stationary variance is proportional to the local probability of actual shifts, rather than stay decisions, and inversely proportional to the combined sensitivity of the left and right plateau-pass probabilities to perturbations of the equilibrium level. In the local-balance view, this is the same mechanism as \(\sigma_{J,*}^2\approx-1/g'(j_*)\): a sharper peak of the stationary distribution corresponds to a smaller variance.

The grid-level variance can be converted to the ensemble-size scale by 23 . For the original and modified update rules, \[\label{eq:sigma95B95variants} \sigma_{B,*}^{\mathrm{orig}} \approx B_*^{\mathrm{orig}}\log(\mathrm{sf})\,\sigma_{J,*}^{\mathrm{orig}}, \qquad \sigma_{B,*}^{\mathrm{mod}} \approx B_*^{\mathrm{mod}}\log(\mathrm{sf})\,\sigma_{J,*}^{\mathrm{mod}}.\tag{27}\] Combining the delta-method conversion 27 with the equilibrium-size formulas 19 gives the following scaling result.

The proof is given in Appendix 6. Although the explicit forms of the grid-level variance constants 25 –26 are not substituted into 28 –29 , Proposition [prop:stationary95variance95constants] is used through the fact that \(\sigma_{J,*}^{\mathrm{orig}}(\mathrm{sf})\) and \(\sigma_{J,*}^{\mathrm{mod}}(\mathrm{sf})\) are independent of \(\varepsilon\).

Thus, \(B_*\) and \(\sigma_{B,*}\) have the same order in \(\varepsilon\), whereas the variance on the \(B\)-scale has the squared order. This is different from a Poisson-type scaling, where the variance is of the same order as the mean. The difference is caused by the geometric grid: a one-level fluctuation around \(B_*\) corresponds to an absolute change of order \(B_*\) in the number of trees.

The relative leading spread, i.e. the coefficient-of-variation-like quantity on the ensemble-size scale, is \[\frac{\sigma_{B,*}}{B_*} \approx \log(\mathrm{sf})\,\sigma_{J,*},\] which is independent of \(\varepsilon\) in the leading model. Therefore, decreasing \(\varepsilon\) increases both the equilibrium ensemble size and the absolute spread around it, while preserving the leading relative spread.

Table 2 provides a numerical evaluation of the constants that enter the stationary-center and stationary-spread formulas. For each fixed value of \(\mathrm{sf}\), the one-dimensional nonlinear balance equations defining \(y_*^{\mathrm{orig}}\) and \(y_*^{\mathrm{mod}}\) are first solved, and the resulting roots are then substituted into 25 –26 . Thus, the table is not an empirical validation on real data sets, but a numerical consequence of the analytical stationary model. It quantifies the magnitude of the constants that determine the relative spread on the ensemble-size scale.

The value \(\mathrm{sf}=\sqrt{2}\) is included because it corresponds to splitting the classical doubling step \(\mathrm{sf}=2\) into two equal multiplicative substeps. Although \(\sigma_{J,*}\) decreases as \(\mathrm{sf}\) increases, the relative spread coefficient \(\log(\mathrm{sf})\sigma_{J,*}\) increases. This is consistent with the geometric parameterization \(B_j=T_0\cdot\mathrm{sf}^j\): larger scale factors make the stationary distribution narrower in grid levels, but each one-level fluctuation corresponds to a larger multiplicative change in the ensemble size.

The last two columns also clarify the equilibrium activity of the two update rules. For the modified rule, \(p_*^{\mathrm{mod}}\) remains close to \(1/4\) over the reported range of scale factors. Since \(p_*^{\mathrm{mod}}\) is the common left- and right-shift probability at equilibrium, the total probability of a non-stay move is \(2p_*^{\mathrm{mod}}\approx 1/2\). Thus, in the symmetric modified rule, the process spends roughly half of the stationary steps moving to a neighboring grid level and roughly half staying at the current level.

For the original rule, the common shift probability \(p_*^{\mathrm{orig}}\) is slightly above \(1/3\) for the practically relevant scale factors \(\mathrm{sf}=\sqrt{2}\) and \(\mathrm{sf}=1.5\). Consequently, the left- and right-shift probabilities are each a little larger than one third, while the stay probability is a little smaller than one third. This is qualitatively consistent with Figure 1, where the central triplet point continues to move frequently in the stationary regime rather than remaining fixed for long stretches.

This also reframes classical stopping criteria based on comparing two neighboring ensemble sizes, for example \(T\) and \(2T\). In the present notation, when the central point is \(B=2T\) and \(\mathrm{sf}=2\), such a rule corresponds to a single left-gap test \(d_L\le\varepsilon\). Equivalently, it can be viewed as a right-gap test at the previous grid level, but the left-gap formulation is the one naturally associated with a left-to-right early-stopping scan. However, Table 2 shows that, for \(\mathrm{sf}=2\), the grid-level stationary standard deviations \(\sigma_{J,*}^{\mathrm{orig}}\) and \(\sigma_{J,*}^{\mathrm{mod}}\) are both already close to one grid step. For \(\mathrm{sf}=\sqrt{2}\) and \(\mathrm{sf}=1.5\), both standard deviations are even larger than one grid step, although still below one and a half steps. Thus, the two ensemble sizes involved in an adjacent comparison are separated by a distance comparable to the intrinsic stationary spread of the plateau process. This places the classical, still widely used adjacent-score stopping criterion under a specific criticism: it should be interpreted as a noisy local diagnostic of plateau behavior rather than as a deterministic stopping certificate.

4.10 Implications for Stationary Estimation↩︎

The preceding results identify \(B_*\) and \(\sigma_{B,*}\) as population-level stationary quantities of the plateau process. This suggests that practical estimators should target the stationary regime rather than a first hitting time or a single best trial. The construction and benchmarking of such estimators are left to separate work.

5 Discussion↩︎

The analysis in this paper shifts the interpretation of triplet-based plateau search from a stopping mechanism to a stationary stochastic process. After the non-\(T\) hyperparameters have effectively stabilized, the central triplet point \(B_t\) is not expected to converge to a deterministic ensemble size. Instead, it fluctuates around a stationary region whose center and spread are determined by the tolerance, the scale factor, and the stochastic behavior of the OOB score gaps.

A first consequence is that the tolerance \(\varepsilon\) has a quadratic effect on the stationary center. Under the leading centered folded-normal model, both the original and modified update rules satisfy \(B_*=O(\varepsilon^{-2})\). Thus, decreasing \(\varepsilon\) by a factor of two shifts the equilibrium ensemble size by approximately a factor of four. The same order also appears for the stationary standard deviation on the ensemble-size scale, \(\sigma_{B,*}=O(\varepsilon^{-2})\). Consequently, decreasing \(\varepsilon\) increases both the typical number of trees and the absolute magnitude of the stationary fluctuations around it.

This behavior should not be interpreted as a weakness of the plateau procedure. It reflects the fact that the process evolves on a geometric grid. A one-level fluctuation near \(B_*\) corresponds to a multiplicative change in the ensemble size, and hence to an absolute change of order \(B_*\). For this reason, the relative leading spread \(\sigma_{B,*}/B_*\) is independent of \(\varepsilon\) to first order and is controlled mainly by \(\mathrm{sf}\) and by the update rule.

The numerical constants in Table 2 make this effect concrete. For the scale factor \(\mathrm{sf}=1.5\), used in the plateau experiments illustrated in Figure 1, the conservative original rule gives \(\log(\mathrm{sf})\sigma_{J,*}\approx0.526\), while the symmetric modified rule gives \(\log(\mathrm{sf})\sigma_{J,*}\approx0.486\). Thus, the leading standard deviation on the ensemble-size scale is about one half of the stationary mean in both cases. This order of magnitude is qualitatively consistent with Figure 1, where the stationary fluctuations of the central triplet point are visibly substantial rather than negligible. A similar one-half order is obtained for the half-doubling grid \(\mathrm{sf}=\sqrt{2}\), where the corresponding constants are \(0.488\) and \(0.449\). This is a substantial spread, and it becomes even larger under the classical doubling grid \(\mathrm{sf}=2\), where the constants increase to \(0.677\) and \(0.633\), respectively. This behavior is an intrinsic consequence of using stochastic plateau decisions on a geometric grid rather than a deterministic one-shot stopping threshold.

The comparison between the original and modified update rules also clarifies the role of the mixed case \(d_L\le\varepsilon<d_R\). The original rule assigns this case to a right shift, making it more conservative and moving the stationary center to larger ensemble sizes. The modified rule assigns the same case to staying, producing a more symmetric transition structure. The theory shows that this design choice changes both the equilibrium location and the relative spread by constants depending only on \(\mathrm{sf}\). Therefore, the choice between the two rules is not merely an implementation detail: it controls a tradeoff between conservative overestimation and a more symmetric stationary regime.

The scale factor \(\mathrm{sf}\) plays a second, distinct role. Smaller scale factors make each right probe less expensive because the right endpoint \(R_t=B_t\cdot\mathrm{sf}\) is closer to the current central point. They also reduce the relative spread on the ensemble-size scale. However, smaller scale factors imply a finer grid and can require more update steps to move across a large range of tree counts. If \(\mathrm{sf}\) is chosen too close to one, the neighboring ensembles become nearly indistinguishable, and the relative score gaps may become negligible compared with the tolerance. In that case, both plateau tests tend to pass, and the update rule is biased toward left shifts. Thus, \(\mathrm{sf}\) should be viewed as a resolution parameter of the stationary search process, rather than only as a multiplicative expansion factor.

The geometric grid also has a scaling justification. For a general grid with local increment \(\Delta(B)\), the relevant separation between neighboring ensembles is the relative increment \(\Delta(B)/B\). If \(\Delta(B)=o(B)\), then this relative separation vanishes as \(B\to\infty\), and the neighboring forests become asymptotically indistinguishable on the scale of the plateau tests. For any fixed tolerance, such a regime makes the signed score gaps increasingly likely to fall inside the plateau interval, producing an artificial tendency toward left shifts. Thus, a non-degenerate plateau comparison requires \(\Delta(B)\) to be of order \(B\), at least asymptotically. The geometric grid \(B_j=T_0\cdot\mathrm{sf}^j\) is the simplest construction with this property, since \(B_{j+1}-B_j=(\mathrm{sf}-1)B_j\).

The theory is derived under several idealizations. It isolates the ensemble-size dynamics after the remaining hyperparameters have stabilized, approximates the signed OOB score gaps by a centered Gaussian model whose absolute values determine the plateau-pass probabilities, and uses an auxiliary birth–death chain in which the leading folded-normal transition probabilities are treated as exact. These assumptions are useful because they expose the dominant scaling laws and yield closed-form constants, but they should not be interpreted as claiming that finite HPO trajectories follow the limiting model exactly.

These limitations also indicate natural directions for empirical follow-up. One can test the normal approximation for the signed relative gaps directly, for example using diagnostic plots or normality tests such as Shapiro–Wilk on suitably standardized local gap samples. One can also examine whether the empirical relative spread \(\sigma_{B,*}/B_*\) remains approximately stable as \(\varepsilon\) varies, as predicted by the leading theory. Such validation requires repeated trajectories or sufficiently long stationary segments, and is therefore better suited to a separate empirical study focused on finite-sample estimation.

The present work deliberately stops before constructing a full trajectory-based estimator. The results identify the stationary center \(B_*\) and the stationary spread \(\sigma_{B,*}\) as population-level quantities of the plateau process. A practical estimator could be based on empirical transition probabilities, local regression of the log-balance ratio, or moment equations derived from the stationary conditions. The construction and benchmarking of such estimators are left to separate work, where they can be evaluated against early-stopping heuristics and fixed-grid HPO baselines without overloading the theoretical development presented here.

6 Proofs↩︎

References↩︎