Gram–Schmidt Forward Regression (GSFR)

• Gram–Schmidt Forward Regression (GSFR) is a technique that uses sequential orthogonalization to robustly estimate coefficients and reduce multicollinearity.
• It adapts to contexts such as temporally ordered regressions, high-dimensional variable selection, and exact least squares computation while ensuring causal interpretability.
• Empirical studies show that GSFR yields unbiased estimates with lower standard errors and improved computational efficiency compared to traditional OLS.

Searching arXiv for papers on Gram-Schmidt Forward Regression and related formulations. Gram–Schmidt Forward Regression (GSFR) is a label used in the arXiv literature for regression procedures that apply Gram–Schmidt orthogonalization in a forward order and then use the resulting orthogonal directions for estimation, causal interpretation, or variable selection. In the papers considered here, the term covers three related but non-identical constructions: a causal-inference estimator for temporally ordered regressors, a forward selection procedure for ultra-high dimensional linear regression, and a computational route to exact ordinary least squares without forming a pseudo-inverse matrix (Cross et al., 2024, Chen et al., 7 Jul 2025, Christopoulos, 2013, Madar et al., 2023). This suggests a common computational motif—sequential orthogonalization—but not a single invariant target parameter.

1. Sequential orthogonalization and the forward-regression idea

In the causal formulation of Cross and Buccola, let $X=[x_1,\dots,x_K]$ be a full-rank $N\times K$ matrix of centered regressors and let $y$ be the $N$-vector of outcomes. GSFR orthogonalizes the columns of $X$ in their given order by a Modified Gram–Schmidt recursion:

$z_1 \leftarrow x_1,$

and for $j=2,\dots,K$,

$$z_j \leftarrow x_j-\sum_{i=1}^{j-1}\frac{x_j' z_i}{z_i' z_i}\,z_i.$$

The resulting matrix $Z$ has pairwise orthogonal columns, so $Z'Z$ is diagonal. The least-squares fit of $N\times K$0 on $N\times K$1 therefore has the closed-form “one-step” estimator

$N\times K$2

In that formulation, GSFR is literally “forward” Gram–Schmidt followed by “one-step” regression (Cross et al., 2024).

The QR-based formulations use the same forward orthogonalization idea for exact least squares computation. In one version, GSFR replaces the direct inversion of $N\times K$3 by a thin QR factorization $N\times K$4, with $N\times K$5 and $N\times K$6 upper triangular, and then computes

$N\times K$7

A non-normalized variant constructs orthogonal vectors $N\times K$8 by

$N\times K$9

forms the upper-triangular matrix $y$0 with entries $y$1, and solves $y$2 by backward substitution, where $y$3 (Christopoulos, 2013, Madar et al., 2023).

2. Finite-sample properties in ordered linear systems

For the causal GSFR estimator, finite-sample theory is established under the usual Gauss–Markov assumptions: linear data-generating recursion or LSEM in the same order as $y$4, $y$5, homoskedasticity, and no autocorrelation. Under those conditions, several properties are stated explicitly (Cross et al., 2024).

• Orthogonality: $y$6 for all $y$7.
• Unbiasedness: $y$8 equals the corresponding reduced-form total-effect parameter $y$9 in the recursive system.
• Variance formula:

$N$0

• Variance comparison with OLS: by the Frisch–Waugh–Lovell theorem and properties of the Schur complement, one proves $N$1 is at least as large as the OLS denominator, hence

$N$2

• Stability:

$N$3

for all $N$4, so there is no coefficient-covariance.

• Information preservation: the $N$5 of the $N$6 regression equals the $N$7 of $N$8.
• Omitted-variable behavior: dropping any later regressor $N$9 with $X$0 does not change $X$1, and including irrelevant later regressors does not inflate the variance of $X$2.

These properties distinguish the ordered-regressor GSFR estimator from ordinary least squares in settings with temporally ordered and collinear regressors. In the language of the paper, coefficients are unbiased and stable with lower standard errors than those from Ordinary Least Squares (Cross et al., 2024).

3. Temporal order, recursive structure, and treatment-effect interpretation

The causal interpretation of GSFR enters when the regressors $X$3 arise in time order and the structural system is recursive, so that each $X$4 has a direct-plus-indirect total effect on $X$5 through later regressors. In that setting, $X$6 is an unbiased estimate of the total derivative $X$7, understood as the sum of direct and downstream indirect effects (Cross et al., 2024).

A central distinction is between “late” and “early” treatments. When $X$8 is a late treatment, meaning it is applied after all other covariates are fixed, Frisch–Waugh–Lovell implies that $X$9. The coefficient then admits the usual Angrist–Imbens “convex-weight” interpretation or ATTT/ATTU decomposition of a heterogeneous treatment effect. When $z_1 \leftarrow x_1,$0 is an early treatment, meaning it precedes and causally drives later regressors, the paper shows that

$z_1 \leftarrow x_1,$1

Accordingly, under zero-conditional-mean and linear-probability ignorability assumptions, $z_1 \leftarrow x_1,$2 is directly interpretable as the Average Total Treatment Effect on the Treated (ATTT) (Cross et al., 2024).

This formulation is designed for settings in which multicollinearity is entangled with causal ordering. The paper’s claim is not merely computational. Rather, chronological ordering changes coefficient interpretation: GSFR is intended to recover direct-plus-indirect causal total effects when a sensible temporal ordering of covariates is available (Cross et al., 2024).

4. Extension to ordered blocks and simultaneous regressors

The same paper extends GSFR to data in which some regressors are simultaneous within blocks while blocks themselves are ordered in time. The design matrix is partitioned into $z_1 \leftarrow x_1,$3 ordered blocks $z_1 \leftarrow x_1,$4, where within each block regressors may be simultaneous and mutual partial regressions are therefore omitted (Cross et al., 2024).

The block procedure is defined as follows. Block 1 is kept as is. For block 2, each column of $z_1 \leftarrow x_1,$5 is regressed on all of $z_1 \leftarrow x_1,$6, and the residuals are stored as $z_1 \leftarrow x_1,$7. More generally, for block $z_1 \leftarrow x_1,$8, $z_1 \leftarrow x_1,$9 is regressed on $j=2,\dots,K$0, and the residuals are stored as $j=2,\dots,K$1. The final regression is then run on

$j=2,\dots,K$2

The coefficient on each $j=2,\dots,K$3 is interpreted as the partial derivative of $j=2,\dots,K$4 along the path that enters through block $j=2,\dots,K$5’s $j=2,\dots,K$6-th regressor, preserving all indirects via later blocks but excluding unidentifiable feedback within a block (Cross et al., 2024).

The paper states that the finite-sample properties established for scalar ordering carry through mutatis mutandis in the block setting: orthogonality across blocks, unbiasedness for block-total effects, variance less than or equal to OLS, and zero omitted-block bias. This generalization places GSFR between fully recursive systems and designs with contemporaneous simultaneity, rather than forcing a single-variable time order on all regressors (Cross et al., 2024).

5. Reanalyses of the NSW program and NLSY reading scores

Cross and Buccola illustrate GSFR by reanalyzing two studies that controlled for temporally ordered and collinear characteristics, including race, education, and income. In both applications, the paper states that GSFR removes the collinearity among regressors so that each coefficient is estimated in isolation, yields strictly lower standard errors than OLS, recovers direct-plus-indirect causal total effects when regressors are chronologically ordered, and generalizes to mixed recursive/simultaneous designs by orthogonalizing across blocks rather than within them (Cross et al., 2024).

In the National Supported Work (NSW) program application, $j=2,\dots,K$7 included race (black/not), age, education, degree-status, marriage, and final treatment program-participation. OLS reports a direct “black” effect of $j=2,\dots,K$8, whereas GSFR reports a total effect of $j=2,\dots,K$9. The extra $$z_j \leftarrow x_j-\sum_{i=1}^{j-1}\frac{x_j' z_i}{z_i' z_i}\,z_i.$$0 points is described as the systemic component via lower education, marriage, and related downstream channels. The program itself was a late treatment, so its GSFR estimate equals OLS at $$z_j \leftarrow x_j-\sum_{i=1}^{j-1}\frac{x_j' z_i}{z_i' z_i}\,z_i.$$1. Race, treated as an early treatment, admits $$z_j \leftarrow x_j-\sum_{i=1}^{j-1}\frac{x_j' z_i}{z_i' z_i}\,z_i.$$2 with $$z_j \leftarrow x_j-\sum_{i=1}^{j-1}\frac{x_j' z_i}{z_i' z_i}\,z_i.$$3, and OLS mixing-weight decompositions confirm that GSFR recovers the ATTT exactly (Cross et al., 2024).

In the NLSY reading-scores application, $$z_j \leftarrow x_j-\sum_{i=1}^{j-1}\frac{x_j' z_i}{z_i' z_i}\,z_i.$$4 included mother’s race, age, education, and AFQT; spouse presence and education; child’s age and gender; family size and log-income; and year dummies. OLS yields a direct “nonwhite” effect of $$z_j \leftarrow x_j-\sum_{i=1}^{j-1}\frac{x_j' z_i}{z_i' z_i}\,z_i.$$5, which is insignificant, whereas GSFR yields a total effect of $$z_j \leftarrow x_j-\sum_{i=1}^{j-1}\frac{x_j' z_i}{z_i' z_i}\,z_i.$$6 points. The paper attributes the shift to mother’s lower school-grade completion, reported as $$z_j \leftarrow x_j-\sum_{i=1}^{j-1}\frac{x_j' z_i}{z_i' z_i}\,z_i.$$7 grades, multiplied by a child-grade effect of $$z_j \leftarrow x_j-\sum_{i=1}^{j-1}\frac{x_j' z_i}{z_i' z_i}\,z_i.$$8 points, plus AFQT chains and related indirect paths. Family log-income is a late treatment: OLS gives $$z_j \leftarrow x_j-\sum_{i=1}^{j-1}\frac{x_j' z_i}{z_i' z_i}\,z_i.$$9 and GSFR gives the same $Z$0 with no inflation of standard error, together with a decomposed ATTT of $Z$1 in a larger heterogeneous-effect decomposition. The authors summarize the method as expanding Bohren et al.’s decomposition of systemic discrimination into channel-specific effects and improving significance levels (Cross et al., 2024).

6. GSFR for ultra-high dimensional forward variable selection

A different use of the name GSFR appears in the ultra-high dimensional linear-regression literature. Here the model is centered,

$Z$2

with $Z$3, $Z$4, and mean-zero errors independent of $Z$5. At stage $Z$6, after selecting indices $Z$7, the method maintains for each candidate variable its Gram–Schmidt residual $Z$8, namely the part of $Z$9 orthogonal to the span of the selected variables, together with the current fitted value of $Z'Z$0 on the selected set. GSFR then selects the next variable by the largest absolute unique-contribution correlation, defined in population form as

$Z'Z$1

The sample algorithm initializes with $Z'Z$2, $Z'Z$3, and $Z'Z$4, iterates up to a pre-specified maximum $Z'Z$5 such as $Z'Z$6, updates the fitted value by regressing the residual on the newly selected orthogonal direction, and updates all remaining candidates by rank-one Gram–Schmidt residualization (Chen et al., 7 Jul 2025).

This GSFR is explicitly compared with Forward Regression (FR) and the Orthogonal Greedy Algorithm (OGA). FR adds the variable that gives the greatest decrease in residual sum of squares and requires refitting OLS on increasingly large subsets via matrix projections. OGA selects by the absolute marginal correlation of the raw residual with each original $Z'Z$7, then re-orthogonalizes the design. GSFR differs from OGA only in the denominator: OGA normalizes by the total variance of $Z'Z$8, whereas GSFR normalizes by the variance of the unique part of $Z'Z$9 uncorrelated with selected variables. The paper states that this correction avoids over-valuing collinear noise variables. It also states that GSFR and FR produce the same selection path up to stopping, but GSFR replaces full-matrix projections by rank-one Gram–Schmidt updates, reducing per-step cost from $N\times K$00 to $N\times K$01 for $N\times K$02 (Chen et al., 7 Jul 2025).

The method introduces a model-size selection rule based on ratio drops in marginal-variance reductions along the GSFR path. With nested sets $N\times K$03, the ideal stopping index is

$N\times K$04

or $N\times K$05 if none. The final model size is selected by

$N\times K$06

where $N\times K$07 is the sample ratio-drop statistic defined in the paper. Under assumptions including $N\times K$08, $N\times K$09, light tails for $N\times K$10 and $N\times K$11, weak sparsity $N\times K$12, and $N\times K$13, the paper proves a convergence rate

$N\times K$14

sure screening under stronger sparsity and invertibility conditions, and stopping consistency in the sense that $N\times K$15 (Chen et al., 7 Jul 2025).

The empirical study compares GSFR with OGA+HDBIC, FR+BIC, and GSFR with full $N\times K$16 steps in simulated models with $N\times K$17 and $N\times K$18. Metrics include coverage probability, false-negative and false-positive rates, best model size, selected model size, running time, in-sample RSS, and out-of-sample MSPE. The reported findings are that GSFR achieved comparable or higher coverage than FR, much lower false positives and smaller models, and dramatically faster computation than FR; against OGA, GSFR had higher coverage and lower MSPE, especially when predictors were strongly correlated; and stopping at $N\times K$19 gave almost identical performance to full-step GSFR with approximately $N\times K$20 speedup. In the riboflavin gene-expression example $N\times K$21, GSFR selected approximately $N\times K$22–$N\times K$23 genes, ran in approximately $N\times K$24, and achieved the lowest MSPE, approximately $N\times K$25, compared with OGA at approximately $N\times K$26 and FR at approximately $N\times K$27 (Chen et al., 7 Jul 2025).

7. Exact OLS computation, numerical issues, and scope of the term

In another strand of work, GSFR denotes an exact method for solving least squares without computing a pseudo-inverse matrix. One formulation observes the standard linear model

$N\times K$28

and contrasts GSFR with the ordinary least squares expression

$N\times K$29

The stated motivation is to avoid explicit computation of $N\times K$30, which costs $N\times K$31 and can be numerically unstable if $N\times K$32 is ill-conditioned. Using classical Gram–Schmidt to build $N\times K$33 and $N\times K$34, GSFR computes the fitted values by orthogonal projection and recovers the coefficients by solving the upper-triangular system. The paper reports factorization cost $N\times K$35, projection cost $N\times K$36, back-substitution cost $N\times K$37, and total cost $N\times K$38, compared with $N\times K$39 for the pseudo-inverse route. It also states exact equivalence to OLS in exact arithmetic and recommends Modified Gram–Schmidt or reorthogonalization for numerical reliability (Christopoulos, 2013).

The non-normalized formulation dispenses with square-root normalizations. It constructs orthogonal vectors $N\times K$40, forms the upper-triangular matrix $N\times K$41 with diagonal entries $N\times K$42, forms $N\times K$43, and solves

$N\times K$44

by backward substitution. Its stated total cost is approximately $N\times K$45, versus $N\times K$46 for direct computation through $N\times K$47. The paper notes that non-normalized classical Gram–Schmidt can lose orthogonality in floating-point arithmetic when columns are nearly collinear, and suggests re-orthogonalization, Modified Gram–Schmidt ordering, pivoting, or Householder QR for greater stability. It also states that the same code generalizes to weighted least squares by replacing the inner product $N\times K$48 with $N\times K$49 (Madar et al., 2023).

Across these papers, GSFR therefore names related but distinct procedures. In the causal setting, it changes coefficient interpretation by imposing temporal order and targeting total effects. In the ultra-high dimensional selection setting, it is theoretically equivalent to FR except for the stopping rule and uses unique-part correlations to construct the selection path. In the QR-based least-squares setting, it is a computational mechanism for obtaining the exact OLS minimizer without forming or inverting the normal equations. A common source of confusion is to treat these as interchangeable; the literature surveyed here instead indicates that GSFR is best understood as a Gram–Schmidt-based forward framework whose statistical meaning depends on the surrounding model class and inferential objective.