Generalized Random Forests Using Fixed-Point Trees

... explain why the derivative matrix can be replaced by an arbitrary scalar

A1: In order to show that choosing a different scalar $\eta$ does not change the stability or outcome of the splitting process, let us compare two types of fixed-point pseudo-outcomes. Consider first $\rho\_{i}=-\eta\_{0}\xi^{\top} n_{C_j}^{-1}\sum_{i\in C_j} \psi_{\hat{\theta}\_{P},\hat{\nu}\_{P}}(O_i)$ with an arbitrary scalar $\eta=\eta\_{0}$. Next, consider $\delta\_{i}=-\xi^{\top} n_{C_j}^{-1}\sum_{i \in C_j}\psi_{\hat{\theta}\_{P},\hat{\nu}\_{P}}(O_i)$, where we choose $\eta=1$. Clearly, we have the relationship: $\rho\_{i}=\eta\_{0}\delta\_{i}$. Using these pseudo-outcomes, the fixed-point approximations for the child node parameters $\hat{\theta}\_{C\_{j}}$ become: $\tilde{\theta}\_{C\_{j}}:=\hat{\theta}\_{P}-n_{C_j}^{-1}\sum\_{i:X\_{i}\in C\_{j}}\rho\_{i}$ for $\rho\_{i}$, and $\bar{\theta}\_{C\_{j}}:=\hat{\theta}\_{P}-n_{C_j}^{-1}\sum\_{i:X\_{i}\in C\_{j}}\delta\_{i}$ for $\delta\_{i}$. Then, the difference between child nodes satisfies a simple scaling relationship

$$||\tilde{\theta}\_{C\_{1}}-\tilde{\theta}\_{C\_{2}}||^{2}=\eta\_{0}^{2}||\bar{\theta}\_{C\_{1}}-\bar{\theta}\_{C\_{2}}||$$

this is simply because

$$||\frac{1}{n\_{C\_{1}}}\sum\_{i:X\_{i}\in C\_{1}}\rho\_{i}-\frac{1}{n\_{C\_{2}}}\sum\_{i:X\_{i}\in C\_{2}}\rho\_{i}||^{2}=\eta\_{0}^{2}||\frac{1}{n\_{C\_{1}}}\sum\_{i:X\_{i}\in C\_{1}}\delta\_{i}-\frac{1}{n\_{C\_{2}}}\sum\_{i:X\_{i}\in C\_{2}}\delta\_{i}||^{2}$$

Based on eq (17), the CART splitting criterion is $\tilde{\Delta}(C\_{1},C\_{2})=\frac{n\_{C\_{1}}n\_{C\_{2}}}{n\_{P}^{2}}||\tilde{\theta}\_{C\_{1}}-\tilde{\theta}\_{C\_{2}}||^{2}$ for $\eta=\eta\_{0}$, $\bar{\Delta}(C\_{1},C\_{2})=\frac{n\_{C\_{1}}n\_{C\_{2}}}{n\_{P}^{2}}||\bar{\theta}\_{C\_{1}}-\bar{\theta}\_{C\_{2}}||^{2}$ for $\eta=1$. Again, these criteria satisfy a clear scaling relationship:

$$\tilde{\Delta}(C\_{1},C\_{2})=\eta\_{0}^{2}\bar{\Delta}(C\_{1},C\_{2})$$

This scaling shows clearly that multiplying the splitting criterion by any positive scalar $\eta\_{0}^{2}$ does not affect the CART splits. Since CART only ranks splits based on the criterion--it selects the partition $(C\_{1},C\_{2})$ that maximizes $\Delta(C\_{1},C\_{2})$:

$$\arg\max\_{(C\_{1},C\_{2})} \tilde{\Delta}(C\_{1},C\_{2})=\arg\max\_{(C\_{1},C\_{2})} \eta\_{0}^{2}\bar{\Delta}(C\_{1},C\_{2})=\arg\max\_{(C\_{1},C\_{2})} \bar{\Delta}(C\_{1},C\_{2}).$$

Therefore, choosing a different scalar $\eta$ does not change the stability or outcome of the splitting process. The absolute scale of $\Delta(C\_{1},C\_{2})$ does not matter; only its ranking across candidate splits determines the final partition.

... Including visualizations ... to contrast the gradient-based and fixed-point splitting .... could enhance comprehension

A2: Thank you for this valuable suggestion. Below, we provide a simple, clear visualization that illustrates why the fixed-point approximation is effective for tree splitting, comparing it directly with the gradient-based method.

We simulate data from a simple VCM with 1D covariates $X\_{i}\sim U([0,1])$ and 2D regressors $W\_{i}\sim N\_{2}(0,\mathbb{I})$. Outcomes are generated from: $Y\_{i} = W\_{i}^{\top}\theta^{\star}(X\_{i})+\epsilon\_{i}$ where $\theta^{\star}(x):=(\sin(2\pi x),x)$, $\epsilon\_{i}\sim N(0,1)$ with $n=5000$.

To simplify, we focus clearly on just a single candidate split at a threshold $t$. We define candidate partitions $C\_{1}(t):=\\{X\_{i}:X\_{i}\leq t\\}$ and $C\_{2}(t):=\\{X\_{i}:X\_{i}>t\\}$ and denote the true splitting criterion by $\Delta(t)$. We compare $\Delta(t)$ with three approximations:

1. Gradient based: $\tilde{\Delta}^{grad}(t)$

2. Fixed point with $\eta=1$: $\tilde{\Delta}\_{1}^{FPT}(t)$

3. Fixed point with $\eta=\sqrt{2}$: $\tilde{\Delta}\_{2}^{FPT}(t)$

We plot these criteria against all possible splits $t$ from 0 to 1.

https://tinyurl.com/bdkd9tpb

The visualization clearly shows:

• The curves for the true criterion $\Delta(t)$, gradient-based $\tilde{\Delta}^{grad}(t)$, and fixed-point with $\eta=1$, i.e. $\tilde{\Delta}\_{1}^{FPT}(t)$, are all very close.

• Even the fixed-point criterion with $\eta=\sqrt{2}$, i.e. $\tilde{\Delta}\_{2}^{FPT}(t)$, although scaled differently, identifies the same optimal split.

• A vertical dashed line shows that the optimal split chosen by four methods are the same:

$$\arg\max\_{t}\Delta(t)=\arg\max\_{t}\tilde{\Delta}^{grad}(t)=\arg\max\_{t}\tilde{\Delta}\_{1}^{FPT}(t)=\arg\max\_{t}\tilde{\Delta}\_{2}^{FPT}(t)$$

Since CART chooses splits by ranking candidate splits $t$, the absolute scale of each splitting criterion does not matter -- only their rankings determine the chosen split. Thus, our fixed-point method is effective and stable.

... should demonstrate metrics beyond ... a smaller sample size.

A3: Please refer to the comments under A5 for Reviewer 1 (Gxac) regarding the larger-scale experiments.

Why is the saving not $O(K^3)$

A1: In the general case, removing the inverse $A_P^{-1}$ should indeed provide a computational saving of order $O(K^3)$. However, for the specific applications discussed in our paper—such as varying coefficient models (VCM) and heterogeneous treatment effect (HTE) models—the product $\xi^\top A_P^{-1}\psi_{\hat\theta_P,\hat\nu_P}(O_i)$ (eq 20) has a simple analytical form without the need for a full matrix inversion, making the computational saving scale instead as $O(K^2)$, due to the required matrix-vector multiplications. In general settings without an analytical form for $A_P^{-1}$, our method would indeed yield $O(K^3)$ savings.

In (27) an approximation to $\hat\theta_P$, how would this affect the theoretical results in Prop. 4.1?

A2: The approximation does not affect our theoretical results. Here we provide a brief and direct argument. Consider pseudo-outcomes defined with the one-step approximation $\tilde\theta_P$:

$$\phi_i^{FPT} := -(W_i - \bar W_P)(Y_i - \bar Y_P - (W_i - \bar W_P)^\top \tilde\theta_P).$$

Let $\bar\theta_{C_j}$ be the child node estimator derived from these approximate pseudo-outcomes. A direct calculation shows that the difference between our original estimator $\tilde\theta_{C_j}$ and this approximate estimator $\bar\theta_{C_j}$ satisfies

$$||\tilde\theta_{C_j} - \bar\theta_{C_j}|| = ||n_{C_j}^{-1}\sum_{i \in C_j} (W_i - \bar W_P)(W_i - \bar W_P)^\top (\hat\theta_P - \tilde\theta_P)|| = ||A_{C_j} (\hat\theta_P - \tilde\theta_P)||\leq ||A_{C_j}|| || \hat\theta_P - \tilde\theta_P||,$$

where $A_{C_j} := n_{C_j}^{-1}\sum_{i \in C_j} (W_i - \bar W_P)(W_i - \bar W_P)^\top$. Under GRF's conditions, in a limit where $n_{C_j}\to\infty$, $A_{C_j}\stackrel{p}{\to}Q$ for a PSD matrix $Q$, and hence $||A_{C_j}|| = O_P(1)$. Let $W_P$ and $Y_P$ denote the centered observations with $i$-th rows $W_i - \bar W_P$ and $Y_i - \bar Y_P$, respectively. Then

$$||\hat\theta_P - \tilde\theta_P|| = ||(A_P^{-1} - \gamma\mathbb I) n_P^{-1} W_P^\top Y_P|| \leq ||A_P^{-1} -\gamma\mathbb I|||| n_P^{-1} W_P^\top Y_P||$$

Under GRF's conditions, in a limit where $r := \sup_{i:X_i\in P}||X_i -\bar X_P|| \to 0$ (the radius of the parent node) and $n_{C_j}\to\infty$, we have $||n_P^{-1} W_P^\top Y_P|| = o_P(r, 1/\sqrt{n_P})$, and $||A_P^{-1} - \gamma\mathbb I|| = O_P(1)$. Therefore $||\hat\theta_P - \tilde\theta_P|| = o_P(r,1/\sqrt{n_P})$ and so

$$||\tilde\theta_{C_j} - \bar\theta_{C_j}|| \leq ||A_{C_j}|| || \hat\theta_P - \tilde\theta_P|| = o_P(r,1/\sqrt{n_P})$$

This means the approximation does not alter the asymptotic behavior of our estimator. The theoretical results in Prop. 4.1 remain fully valid and unaffected.

How to generalize approximation (27) beyond linear regression?

A3: The closed form approximation (27) specifically applies to linear-like settings--such as VCM and HTE models--where $\hat\theta_P$ can be replaced by the one-step exact-line-search approximation $\tilde\theta_P$ (27).

In a more general setting, such approximation is not available, but one can still use the general fix-point approximation, which would still lead to computational gain.

additional settings. 1) Small sample size

A4: To address this, we have run an experiment for small n with data generated according to VCM Setting 3, tested on a separate set of 1,000 observations over 50 replications of each model using a forest of 50 trees

In all cases, both the exact and approximate FPT methods are still observe a computational gain relative to GRF-grad, with no material impact on MSE

additional settings. 2) Correlated W

A5: The case of correlated W has been already considered for HTE models in App C.4, Table 2.

additional settings. 3) Nonlinear models beyond VCM and HTE

A6: Due to space constraints, we will demonstrate other nonlinear models in the final manuscript.

Fig 4 and 18, it seems that the speed gain for VCM is higher than HTE. Explain?

A7: The computational gain is higher for VCM than HTE because VCM uses continuous regressors $W_i$, while HTE uses categorical ones. Continuous regressors provide a larger set of candidate splits. Evaluating these extra splits makes the gradient-based method more computationally demanding, so the simpler FPT approach achieves greater relative savings.

comment on the negative computational gain

A8: It appears to be caused by random fluctuations in computation rather than by the FPT algorithm itself. To clarify this, we repeated the experiment five times under the same conditions. The median computational gain is clearly positive, as shown below

These stabilized results are consistent with the gains observed in other HTE experiments

The proofs from theoretical analysis are majorly motivated by Athey et al. (2019). What are the differences and the challenges?

A1: The main difference is our new splitting criterion based on a fixed-point approximation instead of a gradient-based approximation. Proposition 4.1 in our paper, which proves the asymptotic equivalence of our criterion, differs clearly from Proposition 2 in Athey et al. (2019). Our proof directly addresses the unique challenges posed by removing the Jacobian term. These challenges include establishing asymptotic equivalence and ensuring statistical accuracy without relying explicitly on the gradient information.

... lacks a rigorous convergence analysis for fixed-point iteration... does not discuss whether the proposed approach satisfies ... the contraction mapping principle, which are essential for ensuring convergence.

A2: We would like to clarify a misunderstanding regarding the role of fixed-point method in our approach.

• Numerical convergence: Our method does not perform iterative fixed-point optimization. Instead, we take inspiration from fixed-point iteration but use only a single step of it to simplify the tree-splitting criterion. This single step gives us simplified pseudo-outcomes. Then, these pseudo-outcomes go directly into a standard CART algorithm, which decides how to split the data. Because we rely on CART, our method inherits CART’s proven numerical stability. The stopping rules for CART are clear and well-established (such as minimum node size: fewer than 5 observations or less than 5% of the parent samples). Thus, numerical convergence does not depend on the contraction mapping and is always guaranteed.

• Statistical convergence: From a statistical viewpoint, Prop 4.1 shows clearly that our simplified fixed-point criterion $\tilde\Delta^{FPT}$ is asymptotically equivalent to the original gradient-based criterion. Hence, all statistical convergence guarantees from the original GRF method still hold.

In Section 3.4, the choice of step size $\eta$ is arbitrary, with the paper setting $\eta = 1$ without theoretical justification. This assumption may not always hold and could potentially impact the stability of the splitting criterion.

A3: Please refer to our response to Reviewer 4 (DWaw) (A1 and A2) for the justifications supporting the choice of $\eta = 1$, as well as an explanation of why this choice does not impact the stability of the splitting criterion.

The paper primarily compares FPT with gradient-based approaches but does not benchmark it against other potential improvements in recursive partitioning...

A4: Thank you for raising this point. We would like to clarify a misunderstanding. Our method (FPT) only simplifies the criterion used for tree-splitting in generalized random forests. It does not attempt to improve the recursive partitioning algorithm itself. Both our approach (GRF-FPT) and the gradient-based method (GRF-grad) use exactly the same standard CART algorithm for recursive partitioning.

Moreover, to our knowledge, GRF-FPT and GRF-grad are the only existing methods for approximating the splitting criterion in generalized random forests. Thus, no other comparable methods exist for benchmarking.

It seems that the contributions are somewhat limited.

A5: We respectfully disagree. Our method applies broadly to a wide range of important statistical models—local linear regression (Friedberg et al., 2020), survival analysis, missing data problems (Yifan et al., 2023), nonparametric quantile regression, heterogeneous treatment effect estimation, and nonlinear instrumental variables regression (Athey et al., 2019).

More generally, our approach works for any model defined by nonparametric estimating equations of the form in eq (1). It is especially useful in high-dimensional settings, where the computational cost of GRF becomes a bottleneck. Our fixed-point method offers a practical and scalable alternative while preserving statistical guarantees. We believe this broad applicability and efficiency represent a clear and substantial contribution.

... are there specific scenarios where the gradient-based approach still offers advantages?

A6: No. In all our extensive studies, we have not found any clear scenario where GRF-grad outperforms GRF-FPT, either statistically or computationally.

In Fig 10 of App F, e.g. the MSE of the FPT method is higher than that of the gradient-based method under the conditions of VCM Setting 2. The FPT method is not always better than grad. Please elaborate on these results.

A7: Our simulations show that sometimes GRF-grad slightly outperforms GRF-FPT in accuracy, while at other times GRF-FPT slightly outperforms GRF-grad. The differences are typically small (about 3-5%). The key benefit we emphasize is that GRF-FPT achieves nearly identical statistical accuracy with a substantial improvement in computational speed.

...Additional real-world applications in different domains could strengthen the practical relevance of the method.

A1: We agree with the reviewer. Space limits our discussion here. In the camera-ready version, we will add real-world applications from different domains.

...Some additional experiments designed to test robustness to multicollinearity would be valuable.

A2: We conducted an additional experiment using the VCM with highly correlated features. We modified our original VCM Setting 3 by generating correlated features $X_i \sim N(0,\Sigma)$, where $\Sigma_{i,j} = \omega^{|i-j|}$ for $\omega = 0, 0.5, 0.9$.

Below is a clear summary of computational performance (speedup factor of GRF-FPT over GRF-grad) and statistical accuracy (MSE, multiplied by 100 for readability). Results are averages over 50 replications with $n=10,000$, tested on a separate set of $5,000$ samples using a forest of 10 trees:

These results demonstrate clearly that GRF-FPT remains robust, stable, and computationally efficient, even under high multicollinearity. We will include this analysis in the revised manuscript.

The paper could benefit from more discussion of the hyperparameter sensitivity of the approach, particularly regarding the choice of subsample size in the honest splitting procedure.

A3: Because our method changes only the splitting criterion for the recursive partitioning, both GRF-FPT and GRF-grad share exactly the same subsample size $s$. We do not expect subsample size to change the relative computational advantage or statistical accuracy of our method.

To demonstrate, we ran an experiment varying the subsample ratio $s/n$ from 0.25 to 0.75 using VCM Setting 3. The table below summarizes our results (averaged over 50 replications, with a test set of 5,000 samples):

These results show clearly that the statistical accuracy (MSE) of GRF-FPT relative to GRF-grad does not depend strongly on the subsample ratio.

...A deeper discussion of the substantive insights gained from the spatially varying coefficient estimates would enhance the applied contribution.

A4: Our model (eq 36) clearly shows how the influence of housing features changes across different locations in California. Our results from the housing data uncover interesting geographic patterns.

Fig 5 clearly illustrates this point. In major urban centers such as LA, San Francisco, and Sacramento, housing prices decrease as the number of households increases. This likely reflects overcrowding in densely populated areas. In contrast, rural regions show the opposite trend: prices rise slightly as the number of households grows. This suggests that, in sparsely populated rural areas, a modest increase in households makes these places more attractive and livable.

Population size, however, consistently shows a negative effect on housing prices across nearly all of California, with less geographic variation. Higher populations generally mean lower housing prices, highlighting broader statewide pressures on housing affordability.

Have you explored how the computational gains scale with very large datasets? ...

A5: We ran additional experiments to clearly show how our method scales for very large datasets. Using a forest of 10 trees, we tested our method on VCM Setting 3 with $K = 256$, dim($X$) = 2, and data size $n$ up to $500,000$. The results appear below:

These results clearly demonstrate that, even as the dataset grows very large, our method consistently remains faster than GRF-grad. While the relative speedup slightly decreases at first, it stabilizes at a consistent advantage as $n$ grows very large. Thus, our method has no clear computational bottleneck and maintains a robust advantage at scale.

...Are there any practical upper limits to K where memory or other constraints might become problematic?

A6: As the dimension $K$ grows large, memory usage increases for both methods, potentially hitting hardware constraints. Yet our method (GRF-FPT) is clearly less memory-intensive compared to GRF-grad.

Specifically, GRF-grad scales in memory as $O$(number.of.nodes * $K^2$) for each tree. In contrast, our GRF-FPT scales as $O$(number.of.nodes * $K$). Thus, GRF-FPT remains practical even when $K$ becomes large, and memory constraints become significant.