GLM Plug-in Method

• GLM Plug-in Method is a robust approach that reformulates the GLS estimator as an augmented saddle-point system solved via iterative PCG methods.
• It leverages the structure of model matrices and employs OLS plug-in updates to ensure unbiased intermediate iterates with finite convergence.
• The method adapts efficiently to constrained and multivariate GLMs, providing scalability and improved performance in high-dimensional settings.

A General Linear Model (GLM) Plug-in Method refers to a numerically efficient and theoretically robust approach for estimating GLM parameters—particularly the Generalized Least Squares (GLS) estimator—by iteratively combining direct Ordinary Least Squares (OLS) computations with iterative Preconditioned Conjugate Gradient (PCG) solvers. The PCG-based GLM Plug-in Method reformulates the GLS estimator as the solution to an augmented saddle-point linear system, which is then solved using PCG with an indefinite preconditioner. This approach exploits the structure of the model matrices, guarantees unbiasedness at every iteration, and achieves convergence to the full GLS estimator in a finite number of steps, with the ability to extend efficiently to constrained and structured multivariate GLMs (Foschi, 16 Oct 2025).

1. Conceptual Foundations and Motivation

The GLM Plug-in Method originates from the need to compute the GLS estimator efficiently without direct inversion of large or structured matrices, especially when the covariance Σ is nontrivial or when additional constraints and multivariate structures are present. The GLS estimator in a general linear model (y = Xβ + ε, ε ∼ N(0, Σ)) is defined as:

$b_{GLS} = (X^\top Σ^{-1} X)^{-1} X^\top Σ^{-1} y$

Direct calculation can be prohibitive for large-scale or highly-structured problems due to computational and memory costs associated with forming and inverting $X^\top Σ^{-1} X$. The plug-in approach resolves this by expressing the GLS problem as a solution to an augmented system:

$$\begin{bmatrix} Σ & X \ X^\top & 0 \end{bmatrix} \begin{bmatrix} w \ b_{Aug} \end{bmatrix} = \begin{bmatrix} y \ 0 \end{bmatrix}$$

The PCG Plug-in Method operates by iteratively solving this system, where at every iteration, an OLS-like computation is "plugged in" within the PCG update process.

2. Mathematical Structure and Algorithmic Framework

The algorithm centers on the augmented system for estimating β:

$$G = \begin{bmatrix} Σ & X \ X^\top & 0 \end{bmatrix}, \quad x = \begin{bmatrix} w \ z \end{bmatrix}, \quad h = \begin{bmatrix} y \ 0 \end{bmatrix},$$

with $z$ iteratively approximating $b_{Aug}$ (the estimator for β). The system is solved using the PCG method, equipped with an indefinite preconditioner:

$$K = \begin{bmatrix} D & X \ X^\top & 0 \end{bmatrix}^{-1}$$

where $D$ closely approximates $Σ$ but exploits structural features (block diagonal, banded, Kronecker, or sparse) for computational tractability. The explicit action of $K$ involves forms:

• $X^\top Σ^{-1} X$0
• $X^\top Σ^{-1} X$1

Given $X^\top Σ^{-1} X$2, the preconditioned system $X^\top Σ^{-1} X$3 would become the identity and converge in a single step, but in practice $X^\top Σ^{-1} X$4 is chosen as a convenient surrogate to retain computational efficiency.

Each PCG iteration involves the following steps:

• Compute residual based on current solution
• Update direction vectors using preconditioned OLS computations
• Update $X^\top Σ^{-1} X$5 (the plug-in β estimator) until convergence

The iterative process continues until the projected residual norm (measured with respect to $X^\top Σ^{-1} X$6) falls below a threshold.

3. Statistical and Numerical Properties

A distinctive property of this PCG Plug-in Method is the unbiasedness of the intermediate iterates. Theorem 3.7 in (Foschi, 16 Oct 2025) establishes that if the initial $X^\top Σ^{-1} X$7 is set in the null space of $X^\top Σ^{-1} X$8 (e.g. $X^\top Σ^{-1} X$9), then at every iteration $$\begin{bmatrix} Σ & X \ X^\top & 0 \end{bmatrix} \begin{bmatrix} w \ b_{Aug} \end{bmatrix} = \begin{bmatrix} y \ 0 \end{bmatrix}$$0, i.e., no bias is introduced at any stage before final convergence. This property is not shared by PCG methods applied directly to the normal equations.

Convergence is finite in exact arithmetic: the spectral properties of $$\begin{bmatrix} Σ & X \ X^\top & 0 \end{bmatrix} \begin{bmatrix} w \ b_{Aug} \end{bmatrix} = \begin{bmatrix} y \ 0 \end{bmatrix}$$1 imply that the number of iterations is generally determined by the number of distinct eigenvalues (often $$\begin{bmatrix} Σ & X \ X^\top & 0 \end{bmatrix} \begin{bmatrix} w \ b_{Aug} \end{bmatrix} = \begin{bmatrix} y \ 0 \end{bmatrix}$$2), with a large block of unity eigenvalues corresponding to directions already resolved by OLS updates.

Error bounds are given by quadratic forms—e.g., $$\begin{bmatrix} Σ & X \ X^\top & 0 \end{bmatrix} \begin{bmatrix} w \ b_{Aug} \end{bmatrix} = \begin{bmatrix} y \ 0 \end{bmatrix}$$3 non-increasing—and similar forms for $$\begin{bmatrix} Σ & X \ X^\top & 0 \end{bmatrix} \begin{bmatrix} w \ b_{Aug} \end{bmatrix} = \begin{bmatrix} y \ 0 \end{bmatrix}$$4. The error decreases predictably according to a sequence $$\begin{bmatrix} Σ & X \ X^\top & 0 \end{bmatrix} \begin{bmatrix} w \ b_{Aug} \end{bmatrix} = \begin{bmatrix} y \ 0 \end{bmatrix}$$5 governed by the spectral gap.

Numerical experiments demonstrate that the augmented approach (PCG-Aug) outperforms PCG applied to the normal equations (PCG-NE) in terms of robustness (insensitivity to ill-conditioning), statistical efficiency, and convergence rate—particularly in worst-case tail scenarios.

4. Adaptation to Constrained and Structured GLMs

The plug-in principle extends naturally to estimation in constrained and structured GLMs:

• Linear Restrictions: For models where β is partially fixed (e.g., $$\begin{bmatrix} Σ & X \ X^\top & 0 \end{bmatrix} \begin{bmatrix} w \ b_{Aug} \end{bmatrix} = \begin{bmatrix} y \ 0 \end{bmatrix}$$6), the augmented system incorporates the constraints as additional rows, embedding them directly in the iterative scheme. The preconditioner $$\begin{bmatrix} Σ & X \ X^\top & 0 \end{bmatrix} \begin{bmatrix} w \ b_{Aug} \end{bmatrix} = \begin{bmatrix} y \ 0 \end{bmatrix}$$7 is augmented analogously, and the residuals and projections adapt accordingly.
• Multivariate Structures: For systems such as Seemingly Unrelated Regressions (SUR) or Vector Autoregressive (VAR) models, the design matrix $$\begin{bmatrix} Σ & X \ X^\top & 0 \end{bmatrix} \begin{bmatrix} w \ b_{Aug} \end{bmatrix} = \begin{bmatrix} y \ 0 \end{bmatrix}$$8 and covariance $$\begin{bmatrix} Σ & X \ X^\top & 0 \end{bmatrix} \begin{bmatrix} w \ b_{Aug} \end{bmatrix} = \begin{bmatrix} y \ 0 \end{bmatrix}$$9 often possess block diagonal or Kronecker structures. PCG iterates exploit these by efficient OLS solves for each block, leaving only the "unexploited" global structure for the PCG to resolve iteratively.
• Efficient Implementation: Precomputing QR decompositions, block-wise OLS, and using structure in $$G = \begin{bmatrix} Σ & X \ X^\top & 0 \end{bmatrix}, \quad x = \begin{bmatrix} w \ z \end{bmatrix}, \quad h = \begin{bmatrix} y \ 0 \end{bmatrix},$$0 yields rapid per-iteration complexity while minimizing computational overhead. For restricted models, specialized QR updating algorithms (as in Wilkinson and Reinsch) further speed up the required computations.

5. Performance and Practical Implications

Empirical evaluations confirm that, for large or structured design matrices, the PCG-Aug method achieves comparable or superior accuracy to direct solvers at a fraction of the computational time. The plug-in approach allows intermediate iterates to be used as statistically meaningful estimates (since they are unbiased and their error can be bounded), providing flexibility for early stopping or real-time estimation scenarios.

The PCG-Aug framework is particularly valuable for:

• Large-scale models: Where forming or inverting $$G = \begin{bmatrix} Σ & X \ X^\top & 0 \end{bmatrix}, \quad x = \begin{bmatrix} w \ z \end{bmatrix}, \quad h = \begin{bmatrix} y \ 0 \end{bmatrix},$$1 is computationally infeasible.
• Econometric and multivariate time-series models: Where constraints and structured covariance are ubiquitous.
• Situations demanding robust statistical properties and numerical stability: The plug-in iterates maintain desirable expectations, ensuring reliable inference even at intermediate stages.

A plausible implication is that this approach enables scalable and robust GLM estimation in contemporary high-dimensional and structured-data applications.

6. Extensions and Research Directions

The plug-in methodology, as demonstrated, advocates a hybrid approach blending direct and iterative methods. Future directions suggested in (Foschi, 16 Oct 2025) include:

• Exploration of alternative preconditioners $$G = \begin{bmatrix} Σ & X \ X^\top & 0 \end{bmatrix}, \quad x = \begin{bmatrix} w \ z \end{bmatrix}, \quad h = \begin{bmatrix} y \ 0 \end{bmatrix},$$2 tailored to further exploit model or covariance structure, potentially accelerating convergence.
• Extension to models with even more intricate structures (panel data, spatial models), or where the covariance is only approximately sparse or block structured.
• Development of hybrid direct-iterative strategies that dynamically select between OLS, PCG, and plug-in estimators based on problem scale and computational resources.

This trajectory points toward a broad class of structure-exploiting, statistically robust GLM plug-in methods for scalable statistical computing.

---

In summary, the PCG Plug-in Method for GLM estimation encapsulates an iterative, preconditioned approach that leverages OLS-type computations at each step to efficiently and reliably converge to the GLS estimator. It generalizes flexibly to constrained and multivariate settings, maintains unbiasedness at all stages, and achieves computational gains—positioning it as an advanced tool for large-scale and structured GLM analysis (Foschi, 16 Oct 2025).