The reference is GENERALIZED RANDOM FORESTS by ATHEY, TIBSHIRANI and WAGER (2019). They construct a general algorithm to grow trees and forest for estimation of target parameters that are conditional on a covariate vector $X$.
Notation and Setup
Observable data $O_{i}=(Y_{i},W_{i})$, where $Y$ defines the outcome, $W\in\{0,1\}$ defines the binary treatment, and $X\in\mathbb{R}^{p}$ represents a covariates matrix. In matrix form
$$W=\left(\begin{array}{c} W_{1}\\ \vdots\\ W_{N} \end{array}\right)_{N\times1}Y=\left(\begin{array}{c} Y_{1}\\ \vdots\\ Y_{N} \end{array}\right)_{N\times1}X=\left(\begin{array}{ccc} X_{11} & \cdots & X_{1p}\\ \vdots & \ddots\\ X_{N1} & & X_{Np} \end{array}\right)_{N\times p} $$
The paper defines the target parameter $\theta(x)$ (a possible nuisance parameters $\nu(x)$) as solutions to local moment conditions $$\mathbb{E}\left[\psi_{\theta,\nu}\left(O\right)\mid X=x\right]=0.$$
In the example of conditional treatment effect (CATE) estimation we assume $Y=W_ib_i+e_i$ and conditional independence, and the CATE is defined as $\beta(x)=\mathbb{E}[b_i\mid X_i=x]$. Then a for a given contrast $\xi$, the target parameter $\theta(x)=\xi\beta(x)$ is identified via the score function \begin{align*} \psi_{\beta,c}(Y_{i},W_{i})_{\mid X_{i}=x} & =\left[\begin{array}{c} Y_{i}-\beta(x)W_{i}-c(x)\\ \left(Y_{i}-\beta(x)W_{i}-c(x)\right)W_{i} \end{array}\right]_{2\times1} \end{align*} which defines the two moment conditions in experiments, \begin{align*} \mathbb{E}\left[\psi_{\beta,c}\left(Y,W\right)\mid X_{i}=x\right]=0\leftrightarrow\left[\begin{array}{c} \mathbb{E}\left[Y-\beta(x)W-c(x)\mid X_{i}=x\right]\\ \mathbb{E}\left[\left(Y-\beta(x)W-c(x)\right)W\mid X_{i}=x\right] \end{array}\right]_{2\times1}=\left[\begin{array}{c} 0\\ 0 \end{array}\right]_{2\times1} \end{align*} for all $x\in\mathcal{X}$.
A gradient of $\mathbb{E}[\psi\mid X_i=x]$
The generalized random forest performs CART splits on psuedo-outcomes, calculated as follows. Given a parent node $P$, we define the pseudo-outcome $\rho$ by
$$\rho_{i}=-\xi^{T}A_{p}^{-1}\psi_{\hat{\theta}_{P},\hat{\nu}_{P}}(O_{i}),$$ where $A_{P}$ is a consistent estimate for the gradient of $\mathbb{E}\left[\psi_{b,c}\left(Y,W\right)\mid X_{i}\in P\right]$.
In the CATE estimation this amounts to \begin{align} \rho_i &= A_P^{-1}\left(W_i-\bar{W}_P\right)\left(Y_i-\bar{Y}_P-\left(W_i-\bar{W}_P\right)\hat{\beta}_P\right) \\ A_P &= \frac{1}{\left| \{i:X_i\in P\} \right|} \left[ \sum_{\{i:X_i\in P\}} \left(W_i-\bar{W}_P\right)\left(W_i-\bar{W}_P\right)^T \right] \end{align}
My question is: why? Why is $A_p$ the gradient of the local expectation score function $\mathbb{E}[\psi\mid X_i=x]$? How can $A_p$ be defined as a gradient, when the $\mathbb{E}[\psi\mid X_i=x]$ is in fact a matrix of vector?