1 Introduction↩︎

Random feature method (RFM) collocation converts a linear partial differential equation (PDE) into a linear least-squares problem once the trial functions have been drawn. The fixed trial space makes the coefficient solve explicit, but it also exposes a discretization question: which residual equations represent the action of the differential operator on that space? Geometric coverage of the domain can miss weakly observed residual directions and produce a poorly conditioned empirical norm. We address this problem by constructing the collocation measure, equation weights, and coefficient coordinates from the geometry of the operator-applied trial functions.

Let \(\Omega\) be the computational domain, \(L\) the interior differential operator, \(B\) the boundary operator, and \(f\) and \(g\) the corresponding data in \(Lu=f\) in \(\Omega\) and \(Bu=g\) on \(\partial\Omega\). For fixed random trial functions \(\{\varphi_m\}_{m=1}^M\), write \[u_M(\mathbf{x})=\sum_{m=1}^M c_m\varphi_m(\mathbf{x}),\qquad \psi(\mathbf{x})=\big(L\varphi_1(\mathbf{x}),\ldots,L\varphi_M(\mathbf{x})\big)^\top,\] where \(M\) is the number of features, \(c_m\) are the unknown coefficients, \(\mathbf{x}\in\Omega\), and \({}^\top\) denotes transpose. An interior point \(\mathbf{x}_i\) contributes the row \(\psi(\mathbf{x}_i)^\top\); boundary rows are obtained by applying \(B\) at points on \(\partial\Omega\). The Gram matrix generated by \(\psi\) describes the residual directions that the collocation equations must resolve.

RFM collocation has been developed for stationary and time-dependent PDEs [1], [2]. Related fixed-feature solvers include physics-informed extreme learning machines (ELMs) [3], local and high-dimensional ELMs [4], [5], and local randomized neural-network discontinuous Galerkin methods [6]. Strong-form RFM also shares fixed trial functions and collocated operator evaluations with Kansa-type and radial-basis-function methods [7]–[9]. Across these settings, stable coefficient recovery requires the discrete residual norm to capture the relevant trial-space directions.

More broadly, operator- and regime-adapted classical solvers show how numerical performance can improve when the discretization reflects the PDE structure, as in uniformly accurate treatments of highly oscillatory Dirac equations and nonuniform-FFT methods for nonlocal Schrödinger interactions [10]–[12]. Our construction brings the same design principle to the selection and scaling of residual equations in fixed random-feature spaces.

Existing point-selection methods encode different information. Geometric and low-discrepancy rules cover the physical domain. Residual-adaptive physics-informed neural network (PINN) methods refine regions where the current approximation violates the equation [13]–[15]. PDE-greedy kernel methods select operator functionals through power-function or current-error criteria [16]. Optimal weighted least squares and Christoffel sampling provide trial-space-aware designs for regression [17]–[20]. The Christoffel Sampling for Machine Learning (CS4ML) framework further permits general linear observations, including differential observations, and develops Christoffel adaptive sampling for PINNs [21]. Christoffel adaptive sampling has also been paired with sparse random-feature expansions for multivariate function approximation, where the sampling density is formed from the feature values and the features are subsequently sparsified [22].

The fixed linear structure of RFM collocation permits a direct operator-residual construction. In the PINN realization of CS4ML, the computable surrogate used for interior sampling omits the differential operator [21]. Here every feature \(L\varphi_m\) is explicit, so a reference quadrature determines the operator-residual Gram, its retained eigenspace, and the associated leverage scores before the coefficient solve. The same Gram defines a residual-Christoffel sampling density and a coefficient whitening map. This coupling targets the reference residual norm over all retained directions and yields a sampling measure that can be reused for right-hand sides sharing the operator, domain, and fixed trial space. Randomized least-squares preconditioners, including count-sketch preconditioning of the random-feature system itself, instead act on an assembled row set with a fixed point set [23]–[25], while feature filtering changes the retained columns and hence the approximation space [26].

For a reference probability measure \(\mu\) on \(\Omega\), define \[G=\mathbb{E}_{\mu}[\psi(\mathbf{x})\psi(\mathbf{x})^\top],\qquad K_M^L(\mathbf{x})=\psi(\mathbf{x})^\top G_r^\dagger\psi(\mathbf{x}),\] where \(\mathbb{E}_\mu\) denotes expectation with respect to \(\mu\), \(G\) is the population residual Gram, \(r\) is the retained residual rank, \(G_r^\dagger\) is its rank-\(r\) truncated spectral inverse, and \(K_M^L\) is the residual-Christoffel function. The randomized method samples from the density \(K_M^L/r\) relative to \(\mu\) and uses inverse-density weights together with the whitening map derived from \(G\). For scalar residual rows, a deterministic companion selects the candidate giving the largest regularized log-determinant increment. The analysis establishes retained-space Gram concentration and conditioning. For uniformly analytic residual kernels, it further connects the effective residual dimension to spectral decay of the associated random-feature operator.

The main contributions are as follows.

1. We construct residual-Christoffel sampling, inverse-density weighting, and a matching whitening map from the fixed operator-residual features, together with deterministic leverage-greedy selection.

2. We prove retained-space Gram concentration and conditioning, and relate the sampling dimension to the spectrum and ridge effective dimension of the random-feature residual-kernel operator.

3. Numerical experiments across scalar and vector PDEs, varied geometries, and one to three dimensions demonstrate improved conditioning, faster iterative solution, and the lowest small-sample errors from greedy selection in the scalar tests.

The remainder of this paper is organized as follows. Section 2 formulates the residual-space discretization, Section 3 develops the sampling algorithms, and Section 4 gives the conditioning and spectral analysis. Section 5 presents the numerical experiments, followed by the conclusions and outlook in Section 6.

2 Problem formulation and residual-space discretization↩︎

We distinguish the reference residual Gram, sampled interior Gram, boundary-augmented least-squares system, and coefficient map because they govern sampling, assembly, and solver coordinates separately.

2.1 Model problem and fixed-feature trial space↩︎

Let \(\Omega\subset\mathbb{R}^d\) be a bounded computational domain, where \(\mathbb{R}\) denotes the real numbers and \(d\) is the spatial dimension, and let \(L\) and \(B\) be linear interior and boundary operators, respectively. We consider \[\label{eq:model-problem} Lu=f\quad\text{in }\Omega,\qquad Bu=g\quad\text{on }\partial\Omega,\tag{1}\] for the unknown solution \(u\), using scalar notation. Systems of equations are handled by stacking component equations and their coefficient blocks, so one physical point may contribute several equation rows, as in the elasticity experiments. Let \(\mu\) be a reference probability measure on \(\Omega\) and let \(\nu\) be a reference probability measure on \(\partial\Omega\). Unless stated otherwise, these are normalized volume and boundary reference measures. The positive integer \(M\) is the number of realized trial functions, and the analysis is conditional on the fixed realization \(\{\varphi_m\}_{m=1}^M\). We assume that \(L\varphi_m\in L^2(\mu)\) and \(B\varphi_m\in L^2(\nu)\), with pointwise representatives on the collocation sets, and that the data \(f\) and \(g\) can be evaluated at the corresponding points. Here \(L^2(\mu)\) and \(L^2(\nu)\) are the usual spaces of square-integrable functions with respect to the indicated measures.

Throughout, \(\|\cdot\|_2\) denotes the Euclidean norm for vectors and the spectral norm for matrices; \({}^\top\) denotes transpose, \(I_k\) denotes the \(k\times k\) identity matrix, and \(\mathbb{E}\) and \(\Pr\) denote expectation and probability.

In computations, \(\mu\) is represented by a probability-normalized quadrature with \(Q\) nodes \(\tilde{\mathbf{x}}_q\) and nonnegative weights \(w_q\), \(\mathcal{Q}=\{(\tilde{\mathbf{x}}_q,w_q)\}_{q=1}^Q\), where \(\sum_qw_q=1\), so that, for an integrable function \(h\), \(\mathbb{E}_\mu h(\mathbf{x})\) is approximated by \(\sum_qw_qh(\tilde{\mathbf{x}}_q)\). When a boundary Gram is required, an independent quadrature with \(Q_b\) boundary nodes \(\tilde{y}_q\) and weights \(w_q^b\) represents \(\nu\): \(\mathcal{Q}_b=\{(\tilde{y}_q,w_q^b)\}_{q=1}^{Q_b}\) with \(\sum_qw_q^b=1\). All expectations below are over spatial variables after the feature realization has been fixed.

The RFM trial space is \[\label{eq:rfm-trial-space} V_M=\mathrm{span}\{\varphi_1,\ldots,\varphi_M\},\qquad u_M(\mathbf{x})=\sum_{m=1}^M c_m\varphi_m(\mathbf{x}).\tag{2}\] For the fixed trial space 2 , define the operator-residual feature vector \[\label{eq:operator-residual-vector} \psi(\mathbf{x})=\big(L\varphi_1(\mathbf{x}),\ldots,L\varphi_M(\mathbf{x})\big)^\top .\tag{3}\]

2.2 Interior residual Gram and sampling objective↩︎

For symmetric matrices \(C\) and \(D\), write \(C\preceq D\) when \(D-C\) is positive semidefinite (PSD). The interior residual Gram associated with 3 is \[\label{eq:interior-residual-gram} G:=\mathbb{E}_{\mathbf{x}\sim\mu}[\psi(\mathbf{x})\psi(\mathbf{x})^\top]\succeq0.\tag{4}\] The Gram in 4 is the default object used for residual-space sampling. Boundary information forms a separate block, and the augmented Gram below provides full interior-plus-boundary conditioning. Assume \(G\ne0\), and let \((\lambda_j,v_j)\) denote its eigenpairs, ordered so that the eigenvectors are orthonormal and \(\lambda_1\ge\cdots\ge\lambda_M\ge0\). The zero-Gram case contains no interior residual information and is excluded.

The analysis allows any fixed retained rank \(r\) satisfying \(1\le r\le\operatorname{rank}(G)\). This rank is chosen after the feature realization and reference Gram are fixed but before the collocation rows are sampled. In the computations, the following gap-aware numerical rule is used. First set \[\label{eq:rank-selection-rule} \begin{align} r_0&=\#\{j:\lambda_j>10^{-10}\lambda_1\},\\ \mathcal{J}(r_0)&=\big\{j:\max(1,r_0-4)\le j\le\min(M-1,r_0+4),\;\lambda_j>0\big\},\\ r&=\begin{cases} \displaystyle\min\operatorname*{arg\,max}_{j\in\mathcal{J}(r_0)} \frac{\lambda_j}{\lambda_{j+1}}, &2\le r_0\le M-1,\\[2mm] r_0, &r_0\in\{1,M\}. \end{cases} \end{align}\tag{5}\] Here \(\#\) denotes the cardinality of a finite set. When \(\lambda_j>0=\lambda_{j+1}\), the ratio in 5 is interpreted as \(+\infty\). The minimum fixes a tie convention. The analysis requires a fixed \(r\) with \(\lambda_r>0\). For \(M=1\), the second branch gives \(r=1\). After the rank is chosen, define the truncated inverse and whitening map by \[\label{eq:truncated-whitening} \begin{align} G_r^\dagger&=\sum_{j=1}^r\lambda_j^{-1}v_jv_j^\top,\\ T_r&=[v_1/\sqrt{\lambda_1},\dots,v_r/\sqrt{\lambda_r}] \in\mathbb{R}^{M\times r},\\ \psi^\perp(\mathbf{x})&=T_r^\top\psi(\mathbf{x})\in\mathbb{R}^r . \end{align}\tag{6}\] Then \(G_r^\dagger=T_rT_r^\top\) and \(T_r^\top GT_r=I_r\), so \(\mathbb{E}_\mu[\psi^\perp\psi^{\perp\top}]=I_r\) in the retained eigencoordinates. The corresponding residual-Christoffel function is \[\label{eq:residual-christoffel-function} K_M^L(\mathbf{x})=\psi(\mathbf{x})^\top G_r^\dagger\psi(\mathbf{x})=\|\psi^\perp(\mathbf{x})\|_2^2, \qquad \mathbb{E}_\mu[K_M^L]=r .\tag{7}\] The randomized sampling rule uses 7 directly. The deterministic rule uses the same operator-residual rows but updates their leverage metric as points are selected. Points where \(K_M^L=0\) contain no information in the retained coordinates and receive zero probability in the randomized rule.

For a system with total coefficient dimension \(M_{\rm sys}\) and \(s\) equation rows per physical point, collect the residual rows in \(\Psi(\mathbf{x})\in\mathbb{R}^{s\times M_{\rm sys}}\). The same construction applies after replacing \(\psi\psi^\top\) by \(\Psi^\top\Psi\) in \(G\) and replacing \(K_M^L\) by \[K_{M,\mathrm{sys}}^L(\mathbf{x}) =\operatorname{tr}\!\big(\Psi(\mathbf{x})G_r^\dagger\Psi(\mathbf{x})^\top\big) =\|\Psi(\mathbf{x})T_r\|_F^2,\] where \(\|\cdot\|_F\) is the Frobenius norm. In particular, \[\mathbb{E}_\mu[K_{M,\mathrm{sys}}^L] =\operatorname{tr}(T_r^\top G T_r)=r.\] Sampling a point assigns the same inverse-Christoffel weight to all of its equation rows, and its empirical Gram contribution is \[\omega_iT_r^\top\Psi(\mathbf{x}_i)^\top\Psi(\mathbf{x}_i)T_r.\] The resulting point contribution is PSD and its spectral norm is bounded by its trace, so the concentration argument below retains the same \(r/N\) summand bound. This block formulation is used for the elasticity experiment, where \(M_{\rm sys}=2M\) and \(s=2\). The population and finite-candidate formulas below apply unchanged after replacing each scalar outer product by its block counterpart.

For selected interior points or point blocks \(\{(\mathbf{x}_i,\omega_i)\}_{i=1}^{N_i}\), let \(\omega_i>0\) denote the sampling weight before the common assembly normalization introduced below. The empirical whitened Gram is \[\label{eq:empirical-whitened-gram} \widehat G^\perp = \sum_{i=1}^{N_i}\omega_i\, \psi^\perp(\mathbf{x}_i)\psi^\perp(\mathbf{x}_i)^\top .\tag{8}\] The sampling target is therefore \(\widehat G^\perp\approx I_r\) on the retained operator-residual coordinates. Thus row selection constructs a weighted empirical approximation of the reference residual Gram rather than only a geometric covering of \(\Omega\). The immediate objective is algebraic conditioning on the retained interior block. Solution accuracy is evaluated separately. The concentration analysis uses the sampling weights \(\omega_i\); the normalized weights \(\bar\omega_i\) below are reserved for the assembled PDE objective.

2.3 Boundary block and assembled least-squares system↩︎

The actual PDE discretization also contains boundary equations. For interior points \(\{\mathbf{x}_i\}_{i=1}^{N_i}\) with \(N_i\ge1\), normalize the sampling weights by \[\bar\omega_i=\frac{\omega_i}{\sum_{s=1}^{N_i}\omega_s}, \qquad \sum_{i=1}^{N_i}\bar\omega_i=1.\] Let \(\{y_j\}_{j=1}^{N_b}\), with \(N_b\ge1\), be the boundary points and let \(\eta_j>0\) satisfy \(\sum_{j=1}^{N_b}\eta_j=1\). The assembled objective is \[\label{eq:weighted-ls-objective} \begin{align} \min_{c\in\mathbb{R}^M}\quad & \sum_{i=1}^{N_i}\bar\omega_i\,\big(Lu_M(\mathbf{x}_i)-f(\mathbf{x}_i)\big)^2 +\beta\sum_{j=1}^{N_b}\eta_j\,\big(Bu_M(y_j)-g(y_j)\big)^2,\\ & u_M=\sum_{m=1}^M c_m\varphi_m . \end{align}\tag{9}\] Here \(\beta>0\) is the total boundary weight relative to the unit-normalized interior block. Uniform and greedy rules use constant unnormalized interior weights, whereas the randomized rule uses inverse-Christoffel weights. The normalization from \(\omega_i\) to \(\bar\omega_i\) does not change the condition number of the interior block because it is a uniform scalar rescaling, but it fixes the interior scale relative to the boundary block.

Define \[\label{eq:assembled-blocks} \begin{array}{ll} A_i(i,m)=L\varphi_m(\mathbf{x}_i),\quad A_i\in\mathbb{R}^{N_i\times M}, & A_b(j,m)=B\varphi_m(y_j),\quad A_b\in\mathbb{R}^{N_b\times M}, \\[2mm] \bar\Omega_i=\mathrm{diag}(\bar\omega_1,\ldots,\bar\omega_{N_i}), & \Omega_b=\mathrm{diag}(\eta_1,\ldots,\eta_{N_b}). \end{array}\tag{10}\] Here \(\operatorname{diag}\) forms a diagonal matrix from the listed entries. Define the sampled data vectors \(f_i:=(f(\mathbf{x}_1),\ldots,f(\mathbf{x}_{N_i}))^\top\) and \(g_b:=(g(y_1),\ldots,g(y_{N_b}))^\top\). The corresponding base weighted least-squares problem is \[\label{eq:assembled-ls-system} \min_{c\in\mathbb{R}^M}\|\mathcal{A} c-b\|_2^2,\qquad \mathcal{A}:=\begin{bmatrix} \bar\Omega_i^{1/2}A_i\\ \sqrt{\beta}\,\Omega_b^{1/2}A_b \end{bmatrix},\qquad b:=\begin{bmatrix} \bar\Omega_i^{1/2}f_i\\ \sqrt{\beta}\,\Omega_b^{1/2}g_b \end{bmatrix},\tag{11}\] Here \(\mathcal{A}\in\mathbb{R}^{(N_i+N_b)\times M}\) and \(b\in\mathbb{R}^{N_i+N_b}\).

Whitening or preconditioning introduces a full-column-rank coefficient map \(P_{\rm coef}\in\mathbb{R}^{M\times m_{\rm tr}}\), where \(m_{\rm tr}\) is the dimension of the transformed coefficient space, and solves for the transformed coefficient vector \(\widetilde{c}\in\mathbb{R}^{m_{\rm tr}}\): \[\label{eq:transformed-ls-system} \min_{\widetilde{c}\in\mathbb{R}^{m_{\rm tr}}}\|\mathcal{A} P_{\rm coef}\widetilde{c}-b\|_2^2, \qquad c=P_{\rm coef}\widetilde{c}.\tag{12}\] For retained residual whitening, \(P_{\rm coef}=T_r\) and \(m_{\rm tr}=r\). Because \(T_r\) is rectangular when \(r<M\), this choice restricts the coefficient vector to \(\operatorname{range}(T_r)\) as well as whitening that subspace. An untransformed solve in the full coefficient space based on a singular value decomposition (SVD) uses \(P_{\rm coef}=I_M\), where \(I_M\) is the \(M\times M\) identity matrix; an invertible full-rank \(P_{\rm coef}\) changes only the coefficient coordinates. Full-rank ridge and boundary-aware transforms are identified explicitly where they are used. Reported condition numbers and LSQR iterations refer to the matrix \(\mathcal{A} P_{\rm coef}\) actually passed to the solver. Experiments that isolate point selection hold \(P_{\rm coef}\) fixed, whereas combined stabilization results are labeled as sampling plus whitening. Within each collocation-sampling comparison, the boundary points, \(\eta_j\), and \(\beta\) are shared.

Thus the reference Gram and sampling weights determine the interior discretization, whereas \(\beta\) and \(P_{\rm coef}\) determine boundary scaling and solver coordinates, respectively.

3 Residual-space sampling strategies↩︎

The objective is to approximate the weighted residual Gram on the retained eigencoordinates. Randomized residual-Christoffel sampling uses the fixed reference Gram \(G_r^\dagger\), whereas deterministic leverage-greedy selection updates the Gram of the rows already selected. Boundary rows are appended afterward as a fixed block. Throughout this section, \(N=N_i\) denotes the number of selected interior points or point blocks. A point contributes one residual row in the scalar formulation and may contribute several rows for a system. The notation \(K_M\) abbreviates \(K_M^L\) when the operator is clear.

3.1 Randomized residual-Christoffel sampling↩︎

The randomized rule aims to select a small set of weighted interior residual rows whose empirical whitened Gram 8 approximates \(I_r\) in the retained eigencoordinates. The population sampling rule draws interior points from the residual-Christoffel measure \[\label{eq:population-christoffel-design} dp=(K_M/r)\,d\mu,\qquad \omega_i=\frac{r}{N K_M(\mathbf{x}_i)} .\tag{14}\] Here \(p\) is the probability measure obtained by tilting \(\mu\) with the residual-Christoffel density \(K_M/r\) on the support where \(K_M>0\); it is normalized because \(\mathbb{E}_\mu K_M=r\). Conditional on the fixed features, Gram, and retained rank, the points \(\mathbf{x}_1,\ldots,\mathbf{x}_N\) are drawn independently with replacement from \(p\). Each sampled row contributes an unbiased \(1/N\) share of the identity in the retained eigencoordinates: \[\label{eq:population-unbiased-row} \mathbb{E}_p\!\left[\omega_i\,\psi^\perp(\mathbf{x}_i)\psi^\perp(\mathbf{x}_i)^\top\right] =\frac{1}{N} I_r .\tag{15}\] Moreover, every weighted whitened row contribution has the same trace, \[\label{eq:population-row-norm} \omega_i\|\psi^\perp(\mathbf{x}_i)\|_2^2=\frac{r}{N}.\tag{16}\] Equation 15 gives the mean interior Gram, while 16 bounds every PSD rank-one summand by \(r/N\) in spectral norm. Before boundary assembly, they are rescaled to \(\bar\omega_i\) as defined in Section 2; this common scalar rescaling leaves the interior condition number unchanged.

For the discrete sampler, a reference quadrature is first used to approximate the residual Gram and construct the whitening map. The leverage scores used for sampling are then evaluated on a finite candidate set \(\mathcal{C}_P=\{(z_q,\varpi_q)\}_{q=1}^{P}\) with \(\sum_q\varpi_q=1\), where \(\varpi_q\) is the probability mass assigned to candidate \(z_q\). Unless otherwise stated, this candidate set is drawn from the same normalized Lebesgue/reference measure \(\mu\); for a uniform finite set, \(\varpi_q=1/P\). Let \(K_q=K_M(z_q)\). The discrete sampling probability of candidate index \(q\) is \[\label{eq:finite-pool-probability} S_P=\sum_s \varpi_s K_s, \qquad \pi_q=\frac{\varpi_q K_q}{S_P},\tag{17}\] where we require \(S_P>0\). The index \(q_i\) is drawn independently from the discrete distribution in 17 , and the sampled row is assigned the sampling weight \[\label{eq:finite-pool-weight} \omega(q)= \frac{S_P}{N K_q} .\tag{18}\] Candidates with \(K_q=0\) have zero sampling probability in 17 and are omitted. Conditional on the candidate set and whitening map, the finite-set estimator is unbiased for the candidate-measure Gram \[\label{eq:finite-pool-unbiased-gram} \mathbb{E}\!\left[ \sum_{i=1}^N \omega(q_i)\, \psi^\perp(z_{q_i})\psi^\perp(z_{q_i})^\top \mathrel{\Big|}\mathcal{C}_P,T_r\right] =H_P:=\sum_{q=1}^P \varpi_q\, \psi^\perp(z_q)\psi^\perp(z_q)^\top .\tag{19}\] Thus 19 is the finite-set counterpart of 15 , but its target is \(H_P\), not necessarily \(I_r\). When the candidate measure approximates the reference residual Gram, \(H_P\approx I_r\) and \(S_P=\operatorname{tr}(H_P)\approx r\), where \(\operatorname{tr}\) denotes the matrix trace. Hence 18 approaches the population normalization \[\omega(q)\approx \frac{r}{N K_M(z_q)} .\] This finite-set rule is used in the experiments. Sampling is with replacement, so a selected point may appear more than once and is retained as a repeated weighted equation. Christoffel and greedy sampling rules use the same finite candidate set when they are compared directly. Uniform and residual-adaptive points are drawn independently from the reference distribution. The candidate mismatch \(\|H_P-I_r\|_2\) is isolated in Remark [rem:finite-candidate]. Unlike greedy selection, randomized sampling allows arbitrary \(N\) through sampling with replacement.

For experiments that avoid explicit rank truncation, we use \(\gamma>0\) and the ridge score \[\label{eq:ridge-christoffel-score} K_{M,\gamma}(\mathbf{x})=\psi(\mathbf{x})^\top(G+\gamma I_M)^{-1}\psi(\mathbf{x}).\tag{20}\] For the population measure, the normalizing constant is \[d_{\mathrm{eff}}^G(\gamma)=\mathrm{tr}\big(G(G+\gamma I_M)^{-1}\big).\] This quantity is the ridge effective dimension of the finite Gram \(G\) at ridge parameter \(\gamma\). For a finite candidate set, it is \[S_{P,\gamma}=\sum_q\varpi_qK_{M,\gamma}(z_q).\] Equations 17 –18 then apply with \(K_M\) and \(S_P\) replaced by \(K_{M,\gamma}\) and \(S_{P,\gamma}\).

The constructions in this subsection sample only the interior residual rows; boundary rows are appended after converting \(\omega_i\) to the normalized assembly weights \(\bar\omega_i\) in 9 . Algorithm [alg:christoffel] summarizes the resulting randomized finite-candidate construction.

3.2 Deterministic leverage-greedy selection↩︎

The deterministic rule uses leverage relative to the rows already selected, rather than the fixed reference score \(K_M\). This variant is useful when the candidate residual rows have already been evaluated and a nonrandom set is desired. Consider an unweighted candidate set \(\mathcal{C}_P=\{z_q\}_{q=1}^{P}\), choose \(N\le P\), and let \(\lambda_{\rm g}>0\). Weighted candidate measures can be incorporated into the corresponding weighted log-determinant design by replacing \(\psi(z_q)\) with \(\sqrt{\varpi_q}\psi(z_q)\). Starting from \(G_0=\lambda_{\rm g} I_M\), define the remaining set \(\mathcal{C}_P^{(k)}=\mathcal{C}_P\setminus\{\mathbf{x}_1,\ldots,\mathbf{x}_k\}\) and choose \[\label{eq:greedy-selection-rule} \mathbf{x}_{k+1}\in\arg\max_{\mathbf{x}\in\mathcal{C}_P^{(k)}} \psi(\mathbf{x})^\top G_k^{-1}\psi(\mathbf{x}),\qquad G_k=\lambda_{\rm g} I_M+\sum_{j=1}^{k}\psi(\mathbf{x}_j)\psi(\mathbf{x}_j)^\top .\tag{21}\] Here \(G_k\) is the regularized Gram of the residual rows already selected, and \(\psi(\mathbf{x})^\top G_k^{-1}\psi(\mathbf{x})\) is the leverage of a candidate row in the metric induced by that selected set. A large score identifies a direction that is weakly represented by the current rows. For a system of equations, the implementation uses the sum of the equation-row scores as a trace surrogate and incorporates a selected point by successive rank-one updates. If the rows at a point form \(A(\mathbf{x})\in\mathbb{R}^{s\times M_{\rm sys}}\), the exact block log-determinant increment is \(\log\det(I_s+A(\mathbf{x})G_k^{-1}A(\mathbf{x})^\top)\); it is not generally equal to the summed row scores. The parameter \(\lambda_{\rm g}\) ensures \(G_k\succ0\) from the first step and controls the scale assigned to directions not yet represented. The selection rule 21 has a local log-determinant characterization.

The selected greedy rows enter the common least-squares assembly with unit sampling weights before the same boundary and normalization convention is applied. Its accuracy as a function of sample size is evaluated numerically. For scalar equations, each accepted row gives the largest available local increase in the regularized log determinant. For systems, the summed row-leverage score is the trace surrogate described above. Algorithm [alg:greedy] summarizes this deterministic finite-candidate construction, and Figure 3 visualizes the interior collocation sets produced by the three row-selection rules.

3.3 Computational regimes and cost↩︎

Let \(Q\) be the number of reference quadrature points, \(P\) the candidate-pool size, and \(r\) the retained residual dimension. On the reference quadrature, residual-row evaluation costs \(O(QM c_L)\), where \(c_L\) is the per-feature cost of evaluating the action of \(L\). Forming \(G_Q=\Psi_Q^\top W\Psi_Q\) costs \(O(QM^2)\), where \(\Psi_Q\in\mathbb{R}^{Q\times M}\) has row \(q\) equal to \(\psi(\tilde{\mathbf{x}}_q)^\top\) and \(W=\operatorname{diag}(w_1,\ldots,w_Q)\), and diagonalizing this dense \(M\times M\) Gram costs \(O(M^3)\) with \(O(M^2)\) Gram storage.

Candidate-row evaluation costs \(O(PM c_L)\). Once the retained eigenspace is available, rank-\(r\) Christoffel scores cost \(O(PMr)\). Probability normalization costs \(O(P)\); an alias table gives \(O(P+N)\) total preprocessing and sampling work. Approximate ridge scores have the same projection cost when the inverse is represented by a rank-\(r\) or sketched factor. An exact full-dimensional ridge factor costs \(O(PM^2)\) to apply over the candidate set.

For the deterministic greedy rule, recomputing the inverse and all scores at every step costs \(O(N(M^3+PM^2))\). With the \(P\times M\) candidate-row matrix stored, Sherman–Morrison updates cost \(O(PM+M^2)\) per step, giving \(O(N(PM+M^2))\) work and \(O(PM+M^2)\) storage after candidate evaluation. For larger feature counts, a randomized range finder can be applied to the weighted residual matrix \(R_Q=W^{1/2}\Psi_Q\). With sketch size \(s_{\rm sk}\), each pass costs \(O(QMs_{\rm sk})\); a streamed or matrix-free implementation uses \(O((Q+M)s_{\rm sk})\) auxiliary memory in addition to the data needed to evaluate residual rows, whereas explicitly storing \(R_Q\) requires \(O(QM)\) memory. These estimates exclude boundary assembly and the final least-squares solve.

4 Conditioning theory↩︎

Write \(K_M:=K_M^L\). Conditional on the fixed features, reference Gram, whitening map, and rank \(r\), inverse-Christoffel sampling makes the whitened interior Gram concentrate about \(I_r\). For a finite candidate set, the target is \(H_P\), which equals \(I_r\) only when the candidate measure also defines the whitening Gram. These statements quantify interior conditioning. The finite-dimensional probability is over the sampled points or point blocks; the operator-limit probability below is over an independent random-feature realization.

4.1 Residual-Gram concentration and interior conditioning↩︎

4.2 Operator spectrum and ridge effective dimension↩︎

The retained-rank concentration bound is most informative when the effective residual dimension is small relative to the ambient feature count. We relate the normalized feature Gram \(G/M\) to a residual-kernel operator and then state the finite-matrix ridge sampling result separately. Multiplicative convergence of the ridge effective dimension additionally requires trace-resolvent control beyond this operator-norm approximation.

In this subsection, specialize to single-layer random features \(\varphi(\mathbf{x};\theta)=\sigma(w\!\cdot\!\mathbf{x}+b_\theta)\), where \(\sigma\) is the activation function, \(\theta=(w,b_\theta)\) consists of a feature weight \(w\in\mathbb{R}^d\) and bias \(b_\theta\in\mathbb{R}\), and \(\theta\sim\rho\) is drawn from the parameter law \(\rho\). Write \(\zeta_\theta(\mathbf{x})=(L\varphi)(\mathbf{x};\theta)\), so the residual features are independent and identically distributed (i.i.d.) random fields. Assume that \(\rho\) has bounded support and that the required derivatives of \(\sigma\) are bounded on the corresponding compact affine image of \(\Omega\). We also assume \(R_\zeta:=\sup_{\theta\in\operatorname{supp}(\rho)}\|\zeta_\theta\|_{L^2(\mu)}^2<\infty\); this holds, for example, for bounded-coefficient differential operators of finite order under the preceding feature assumptions. For \(\mathbf{x},\mathbf{x}'\in\Omega\), define the residual kernel and its operator on \(H=L^2(\mu)\) by \[\label{eq:residual-kernel-operator} \begin{align} k_L(\mathbf{x},\mathbf{x}')&=\mathbb{E}_{\theta}[\zeta_\theta(\mathbf{x})\zeta_\theta(\mathbf{x}')],\\ (\mathcal{T} v)(\mathbf{x})&=\int_\Omega k_L(\mathbf{x},\mathbf{x}')\,v(\mathbf{x}')\,d\mu(\mathbf{x}'), \qquad v\in H. \end{align}\tag{24}\] The operator \(\mathcal{T}\) is self-adjoint, positive semidefinite, and trace-class, with \[\mathrm{tr}(\mathcal{T})=\mathbb{E}_\theta\|\zeta_\theta\|_H^2\le R_\zeta.\] Let its eigenvalues be \(\tau_1\ge\tau_2\ge\cdots\ge0\). We exclude the trivial case \(\mathcal{T}=0\) and define \(\widetilde{d}=\mathrm{tr}(\mathcal{T})/\tau_1\), the dimensionless trace ratio of \(\mathcal{T}\). We use \(\|\cdot\|_2\) for the induced operator norm on \(H\) as well as for the matrix spectral norm. We reserve \(T_r\) for the finite-dimensional whitening matrix used in the algorithms. For a realized feature set, define \(\Phi_M:\mathbb{R}^M\to H\) by \(\Phi_Mc=\sum_{m=1}^M c_m\zeta_{\theta_m}\). With \(*\) denoting the Hilbert-space adjoint, we have \(G=\Phi_M^*\Phi_M\) and \(S_M=\Phi_M\Phi_M^*\), so their nonzero eigenvalues coincide. This is a sample-covariance-operator setting; sharper concentration results are available under additional distributional assumptions [29], while the elementary bound below is sufficient for the present normalization argument.

The proof is given in Appendix 7.

For a constant-coefficient differential operator of order \(m_L\), \[L=\sum_{|\alpha|\le m_L}a_\alpha\partial^\alpha,\] where \(a_\alpha\in\mathbb{R}\) are the operator coefficients. Here \(\alpha=(\alpha_1,\ldots,\alpha_d)\in\mathbb{N}_0^d\) is a multi-index, \(\mathbb{N}_0\) denotes the nonnegative integers, \(|\alpha|=\sum_{j=1}^d\alpha_j\), and \(\partial^\alpha=\partial_{x_1}^{\alpha_1}\cdots\partial_{x_d}^{\alpha_d}\). The residual feature has the general form \[\label{eq:constant-coefficient-residual-feature} \zeta_\theta(\mathbf{x}) =\sum_{|\alpha|\le m_L}a_\alpha w^\alpha \sigma^{(|\alpha|)}(w\!\cdot\!\mathbf{x}+b_\theta).\tag{25}\] Here \(w^\alpha=\prod_{j=1}^d w_j^{\alpha_j}\). Formula 25 is for scalar constant-coefficient operators; variable coefficients and vector-valued systems require the corresponding product and component structure. Thus operators containing several derivative orders, including Helmholtz and convection–diffusion operators, require a finite sum rather than a single derivative of \(\sigma\).

The proof is given in Appendix 7.

The stretched-exponential form is consistent with classical spectral results for analytic kernels [30] and related kernel-eigenfunction approximation theory [31]; the theorem above states the precise one-variable analyticity condition used here.

A sufficient feature-level condition for the hypothesis of Theorem 3 is that every derivative \(\sigma^{(q)}\), \(0\le q\le m_L\), occurring in 25 be analytic and uniformly bounded on the compact affine image \[\mathcal{Z}_\varrho :=\{w\!\cdot\!\xi+b_\theta:\xi\in\mathcal{E}_\varrho(D),\; (w,b_\theta)\in\mathrm{supp}(\rho)\}.\] Sine satisfies this condition on every fixed polyellipse. For tanh, the parameter support and polyellipse must be chosen so that \(\mathcal{Z}_\varrho\) has positive distance from the pole set \(\{i\pi(k+\tfrac12):k\in\mathbb{Z}\}\), where \(i\) is the imaginary unit and \(\mathbb{Z}\) is the set of integers. The analytic decay conclusion therefore requires this complex-domain separation from the pole set. For variable-coefficient and elasticity residual kernels, the same conclusion follows whenever the kernel satisfies the analytic continuation hypothesis of Theorem 3; Section 5 examines their spectra numerically.

The proof is given in Appendix 7.

We next return to the realized finite Gram. The following continuous-row statement is the finite-dimensional ridge-leverage guarantee used to interpret the sampling experiment.

For the relative finite-matrix ridges used in the experiments, let the dimensionless parameter \(\rho_{\mathrm{rel}}>0\) specify \(\gamma=\rho_{\mathrm{rel}}\lambda_1(G)\), and set \(\ell=\gamma/M\). On the event ?? , if \(\varepsilon_{\rm op}<1\), then \[\rho_{\mathrm{rel}}(1-\varepsilon_{\rm op})\tau_1\le\ell \le\rho_{\mathrm{rel}}(1+\varepsilon_{\rm op})\tau_1.\] The numerical section separately checks normalized-spectrum collapse, the empirical growth of \(d_{\mathrm{eff}}^G(\gamma)\), ridge-relative Gram concentration, and the retained-rank ratio. Over the reported range \(M=400,\ldots,6400\), the observed ratio \(r/M\) decreases from \(0.44\) to \(0.05\).

5 Numerical results↩︎

The numerical study examines effective residual dimension and conditioning. It then considers LSQR convergence and accuracy versus sample size. Controlled comparisons distinguish sampling from boundary treatment, preconditioning, feature pruning, trial-space enrichment, and direct solution. Section 5.1 gives the common setup.

5.1 Experimental setup↩︎

The compared interior sampling strategies are uniform sampling, randomized residual-Christoffel sampling (Algorithm [alg:christoffel]), and deterministic leverage-greedy sampling (Algorithm [alg:greedy]). A residual-adaptive control is included in the accuracy-versus-sample-size experiments: three enrichment stages sample from a density proportional to the squared current PDE residual. For elasticity, the score is the squared Euclidean norm of the two equation residuals. All methods use the same trial space and boundary rows within a test. Ranks, whitening maps, and Christoffel weights follow the constructions in Sections 2 and 3, with the default rank chosen by 5 . Throughout the numerical section, \(r\) denotes this retained residual rank, selected separately for each benchmark before the collocation-sampling rule is changed; ratios such as \(N/r\) therefore report selected interior points per retained residual direction. In scalar tests each point contributes one equation row, whereas an elasticity point contributes a two-row block. Boundary rows are fixed across samplers.

Unless otherwise stated, boundary points use equal weights normalized to total boundary weight \(\beta=1\), so \(\eta_j=1/N_b\) in scalar tests; for vector-valued systems, the same point weight is applied to each boundary component row. Boundary weighting and boundary-aware whitening are varied only in the ablation study.

Unless stated otherwise, the trial functions are \(\varphi_m(\mathbf{x})=\tanh(w_m\!\cdot\!\mathbf{x}+b_m)\) with independent \(w_m\sim\operatorname{Unif}([-a_{\rm rf},a_{\rm rf}]^d)\) and \(b_m\sim\operatorname{Unif}([-a_{\rm rf},a_{\rm rf}])\), where \(\operatorname{Unif}\) denotes the uniform distribution. The feature bandwidth is \(a_{\rm rf}=6\) in 1D and in the variable-coefficient, convection–diffusion, elasticity, and geometry tests, \(a_{\rm rf}=5\) in the smooth 2D/3D Poisson tests, and \(a_{\rm rf}=8\) in the high-frequency Poisson and \(k=6\) Helmholtz accuracy tests. The dense scaling and sketched-construction tests instead use \(\varphi_m(\mathbf{x})=\sin(w_m\!\cdot\!\mathbf{x}+b_m)\) with \(a_{\rm rf}=18\) and \(b_m\sim\operatorname{Unif}[-\pi,\pi]\); the \(k=16\) Helmholtz test uses the same sine law with \(a_{\rm rf}=28\).

Errors are reported as relative discrete \(L^2\) errors on held-out interior test sets. Except in solver comparisons, accuracy is computed by a singular value decomposition (SVD) in the full coefficient space with relative cutoff \(10^{-13}\); runs of the iterative least-squares solver LSQR use absolute and relative tolerances \(10^{-10}\) with the stated iteration caps. Randomized results are reported as medians over the repetitions stated in each caption. Spectrum experiments vary both the feature and reference-quadrature realizations. Accuracy-versus-sample-size experiments fix the feature map, reference quadrature, boundary set, and test set, and vary the point-selection rule and finite candidate set. The solver-scaling experiment varies the complete feature, quadrature, and sampling realization. Dotted large-sample uniform reference errors provide comparison levels. Wall-clock timings were measured on a server with an Intel Xeon Platinum 8476C central processing unit (CPU) using four numerical-library threads.

All scalar tests impose Dirichlet data from a manufactured solution, with the right-hand side obtained by applying the stated operator. The Poisson operator is \(L=-\Delta\). On unit boxes, the Poisson solutions are products of sine functions; the high-frequency Poisson accuracy and scaling tests use \(u=\sin(4\pi x)\sin(4\pi y)\) on \([0,1]^2\). Writing \(x\) and \(y\) for the Cartesian coordinates, the two Helmholtz tests on \([0,1]^2\) use wavenumber \(k\) in \(L=-\Delta-k^2I\), where \(I\) is the identity operator. They use \(k=6\) with \(u=\sin(2\pi x)\sin(2\pi y)\) and \(k=16\) with \(u=\sin(5\pi x)\sin(5\pi y)\). The remaining square-domain tests use \(-\nabla\cdot(a\nabla u)\) with diffusion coefficient \(a=1+0.5\sin(\pi x)\sin(\pi y)\), \(-0.1\Delta u+v\cdot\nabla u\) with convection velocity \(v=(2,1)\), and \(u=\sin(2\pi x)\sin(2\pi y)\) in both scalar cases.

The Navier–Cauchy operator is \[\mathcal{L}_{\rm el}\mathbf{u} =-\mu_{\mathrm{el}}\Delta\mathbf{u} -(\lambda_{\mathrm{el}}+\mu_{\mathrm{el}})\nabla(\nabla\!\cdot\!\mathbf{u}),\] with Lamé parameters \(\mu_{\mathrm{el}}=\lambda_{\mathrm{el}}=1\) and manufactured displacement \(\mathbf{u}=(\sin(\pi x)\sin(\pi y),\sin(2\pi x)\sin(2\pi y))^\top\). The subscripts distinguish these elasticity constants from the reference measure \(\mu\) and the Gram eigenvalues \(\lambda_j\); \(\nabla\) and \(\Delta\) denote the spatial gradient and Laplacian.

Geometry transfer is tested on the L-shaped domain \([-1,1]^2\setminus((0,1]\times[-1,0))\) and the annulus \(0.3\le\|\mathbf{x}\|_2\le1\).

5.2 Effective residual dimension and conditioning↩︎

We first examine the dimension of the resolved residual space and then test how its geometry controls the conditioning and iterative solution of the collocation system.

Figure 4 compares normalized residual-Gram spectra and ridge effective dimensions for smooth Poisson problems. The leading 2D spectra nearly coincide after normalization, in agreement with the operator limit in Theorem 2, while the resolved tails separate as \(M\) increases. The effective dimension remains below the ambient feature count in every dimension, but occupies a larger fraction of the feature space as the spatial dimension increases.

Figure 5 tests whether effective dimension organizes the sample size required for Gram approximation. For representative 1D, 2D, and 3D Poisson residual Grams, it plots \(\|(G+\gamma I_M)^{-1/2}(\widehat G-G)(G+\gamma I_M)^{-1/2}\|_2\) against the normalized sample size \(N/d_{\mathrm{eff}}^G(\gamma)\). The error decreases with normalized sample size in all three cases. The 1D curve crosses unit relative deviation, the 2D curve approaches it, and the 3D curve remains larger, revealing dimension-dependent finite-sample constants. Because the candidate pool grows with \(N\), the curves combine row-sampling concentration with increasingly accurate finite-pool resolution.

Table 1 reports the effect of sampling and whitening on condition numbers. Unless stated otherwise, \(\kappa_2\) refers to the weighted matrix actually solved, including the common boundary block. The untransformed uniform matrices are numerically rank deficient in all but the 3D case, whereas the sampled and whitened matrices are numerically full rank with substantially smaller condition estimates. The scalar Poisson rows use greedy selection and the remaining rows use randomized residual-Christoffel sampling; later controls separate point selection from whitening.

2pt

Figure 6 connects conditioning to iterative convergence on the 2D Poisson benchmark. As \(M\) increases, the untransformed uniform matrix becomes numerically rank deficient and LSQR reaches its iteration cap. Both residual-space constructions maintain numerical full rank and reduce the solve to tens or hundreds of iterations. The paired curves show that the conditioning improvement is directly reflected in LSQR convergence.

5.3 Accuracy versus number of collocation points↩︎

We next examine PDE accuracy for a fixed number of interior collocation points. Conditioning measures the stability of the residual system, whereas discretization error also depends on the approximation space; the experiment therefore evaluates both quantities.

Figure 7 compares the three sampling rules at \(N=\lfloor2.5r\rfloor\), where \(\lfloor\cdot\rfloor\) denotes the floor function. Uniform sampling is most accurate for the two smallest feature spaces. Residual-Christoffel and greedy sampling become more accurate from \(M=400\) onward, and the greedy curve reaches the numerical accuracy floor at the largest feature count.

Figure 8 examines accuracy as a function of \(N/r\) for fixed feature maps. At the smallest ratios, greedy selection has the lowest median error in all four scalar problems and improves the high-frequency Poisson error by more than an order of magnitude relative to both uniform and residual-adaptive sampling. The elasticity errors remain comparable across the four methods. Residual-Christoffel sampling supplies the randomized rule linked to Gram concentration, while the greedy rule is especially effective when only a small number of scalar equations can be assembled.

The high-frequency Helmholtz results also separate stability from approximation quality. Greedy selection approaches the large-sample accuracy reference with only a few equations per retained direction, while residual-Christoffel sampling with whitening produces the well-conditioned transformed system reported in Table 1. The two constructions therefore act on the same residual geometry but emphasize small-sample coverage and randomized norm preservation, respectively.

5.4 Geometry transfer and field-level validation↩︎

Table 2 examines smooth manufactured solutions on L-shaped and annular domains for four scalar operators. Every untransformed uniform system reaches the LSQR cap, whereas residual-Christoffel sampling with whitening produces convergent iterative solves and reduces the condition estimates by several orders of magnitude. The result shows that the algebraic stabilization persists across nonrectangular geometries and operator classes.

2.4pt

Figure 9 complements the table with one field-level Poisson realization on a square, an L-shaped domain, and an annulus. All three use residual-Christoffel rows and an SVD solve in the full coefficient space. The pointwise maps place the largest errors near the outer boundary of the square, the reentrant geometry and exterior edges of the L-shape, and the two annular boundary components.

5.5 Ablations and controls↩︎

This subsection presents controlled ablations that isolate five effects otherwise easily conflated with residual-space sampling: boundary weighting, boundary-aware whitening, solver preconditioning, column pruning after assembly, and trial-space deficiency. Except where explicitly varied, the boundary convention is the rebalanced system with total boundary weight \(\beta=1\).

The first controls concern the appended boundary block. Interior-only whitening stabilizes the sampled residual rows; the condition number of the full matrix additionally depends on the boundary block. Figure 10 (a) varies only the augmented whitening Gram \(G_\alpha=(1-\alpha)G+\alpha G_b\). Incorporating the boundary Gram reduces the full-system condition number by nearly three orders of magnitude while leaving the error unchanged. Boundary-aware whitening therefore stabilizes the appended block without altering the collocation points or boundary penalty.

The boundary penalty controls the relative strength of boundary enforcement. Figure 10 (b) and Table 3 show the same error-conditioning tradeoff across the scalar operators: stronger boundary enforcement improves accuracy while increasing the condition number. The greedy error varies by only about a factor of two over seven orders of magnitude in its regularization parameter, as shown in panel (c). Panel (d) isolates whitening itself; removing it returns the sampled system to numerical rank deficiency.

4pt

The next control fixes the collocation set and separates point selection, solver preconditioning, and column pruning. Direct SVD supplies the accuracy reference for each row set. Unpreconditioned LSQR, ridge-LSQR, and the least-squares solver with randomized normal projection (LSRN) [23], [24] act on the assembled system; rank-revealing QR factorization (RRQR), implemented by column-pivoted QR, changes the retained feature columns [26]. Table 4 reports the resulting accuracy, iterations, and algebraic solve times.

Unpreconditioned and ridge-transformed LSQR reach the iteration cap, whereas LSRN converges in about seventy iterations for both row sets. Thus solver preconditioning remains effective after either collocation rule has been chosen. RRQR pruning to \(k=r\) removes essential approximation directions; retaining larger column sets recovers accuracy, with the nonmonotone errors confirming that column identity matters as well as column count.

2.2pt

The RRQR rows show one side of the approximation-space issue: removing too many feature columns can destroy accuracy even when the collocation set is fixed. A complementary control examines whether point selection can compensate for a missing solution direction. On the singular L-shaped problem, a smooth RFM with \(8000\) interior points leaves a relative \(L^2\) error of \(1.4\times10^{-2}\), whereas appending the known corner-singularity function reaches the numerical precision of the solve. This enrichment separates trial-space deficiency from point placement: residual leverage stabilizes represented directions, whereas enrichment expands the approximation space.

5.6 Computational regime and approximate leverage construction↩︎

Figure 11 compares dense post-assembly solvers on the high-frequency Poisson benchmark. Direct SVD remains the fastest and most accurate method over the displayed feature range. Among the iterative methods, LSRN is fastest and most accurate. Unpreconditioned LSQR and ridge-LSQR are slower and less accurate. The residual-Christoffel curve reports ridge-LSQR on the selected rows.

Figure 12 compares dense Gram eigendecomposition with the randomized range-finder construction from Section 3.3. Beyond the observed crossover, the sketch reduces both construction time and dominant-array storage, while its smaller retained residual subspace produces higher PDE error. The crossover therefore exposes a direct cost-accuracy tradeoff in leverage construction.

6 Conclusion and outlook↩︎

This paper developed an operator-aware collocation strategy for linear PDEs in fixed random feature spaces. The operator-residual Gram defines the residual-Christoffel density, inverse-density weights, and coefficient whitening map; a deterministic scalar-row companion maximizes regularized log-determinant increments. Conditional on the trial space, the randomized construction approximates the whitened residual Gram with \(O(r\log(r/\delta)/\varepsilon^2)\) samples. For uniformly analytic residual kernels, spectral decay yields a polylogarithmic ridge effective dimension. Experiments across scalar and vector PDEs, varied geometries, and one to three dimensions show substantially improved conditioning and iterative convergence, together with the lowest small-sample scalar errors from the greedy construction. These results establish residual-space geometry as the organizing principle for stable strong-form random feature collocation.

Natural extensions combine localized or enriched trial spaces with operator-aware sampling for nonsmooth and multiscale solutions, use streamed residual-Gram factorizations for larger feature sets, and employ block leverage with augmented interior-boundary Grams for coupled or parameterized PDEs. Each direction preserves the principle of resolving the operator-induced geometry of the trial space.

7 Proofs of the spectral-decay results↩︎

This appendix collects the proofs of the kernel-spectral results of Section 4.2.

Acknowledgements↩︎

This work was carried out during Jiale Linghu’s visit to Professor Weizhu Bao at the National University of Singapore. Both authors thank Professor Bao for his valuable advice and helpful discussions.

Data availability↩︎

An initial code release is publicly available in the RC-RFM repository at https://github.com/linghujiale/RC-RFM. The complete experiment scripts and numerical records supporting this manuscript will be deposited in the same repository before publication and are available from the corresponding author on reasonable request.

Funding↩︎

The authors declare that no funds, grants, or other support were received during the preparation of this manuscript.

Competing interests↩︎

The authors have no relevant financial or non-financial interests to disclose.

Ethics approval↩︎

Not applicable. This article does not contain studies with human participants or animals performed by any of the authors.

Author contributions↩︎

Jiale Linghu and Yangshuai Wang contributed equally to the conceptualization, methodology, formal analysis, investigation, validation, numerical implementation, interpretation of results, and manuscript writing. Both authors reviewed, edited, and approved the final manuscript.

References↩︎