1 Foundations of least squares

1.1 Ordinary least squares (OLS) assumptions

Ordinary least squares is a parameter estimation strategy for linear regression that minimizes the sum of squared residuals. Its standard theoretical guarantees rely on assumptions about the error term: linearity in parameters, correct specification of the systematic part, and a particular stochastic structure for the residuals. In the classical setup, the errors have mean zero, are uncorrelated across observations, have constant variance (homoskedasticity), and are often taken to be normally distributed for exact finite-sample inference. When these conditions are violated—especially through correlated errors or non-constant variances—the OLS estimator can remain unbiased (under weaker conditions) but may lose efficiency and, under some forms of misspecification, can yield misleading inference.

1.2 Matrix formulation of linear models

In matrix form, a linear regression model is written as \[ y = X\beta + \varepsilon, \] where \(y\) is the vector of responses, \(X\) is the design matrix, \(\beta\) is the vector of unknown coefficients, and \(\varepsilon\) is the error vector. The OLS estimator is \[ \hat{\beta}_{\text{OLS}}=(X^\top X)^{-1}X^\top y, \] when \(X^\top X\) is invertible. The residual vector is \(r = y - X\hat{\beta}_{\text{OLS}}\). The covariance structure of \(\varepsilon\) is summarized by \[ \mathrm{Var}(\varepsilon)=\Sigma. \] Under the classical OLS assumptions, \(\Sigma\) is proportional to the identity matrix; departures from this structure motivate generalized least squares.

1.3 Motivation for generalized methods

Motivation for generalized least squares arises when the residuals cannot be treated as independent and identically distributed with equal variance. Common scenarios include time series, spatial measurements, repeated measurements, and clustered designs, where dependence is expected by construction. More broadly, if the error covariance \(\Sigma\) is known or can be modeled, incorporating it into the estimation criterion can yield estimators that better match the true noise mechanism. GLS is designed to produce estimates that are efficient for the specified covariance model, and to align inference with the true variability pattern.

2 Core idea of generalized least squares

2.1 Error covariance and why OLS may fail

OLS effectively assumes that each observation contributes equally and independently to the objective function. If \(\Sigma\) departs from a scaled identity matrix—through unequal variances or correlations—then the residuals from different observations do not carry equal information. In that case, minimizing unweighted squared errors can overemphasize observations with high noise or ignore correlation patterns, leading to a loss of statistical efficiency. While unbiasedness may still hold when \(E[\varepsilon]=0\) and the model for the mean is correct, the estimated standard errors and hypothesis tests can become inaccurate if they ignore the true error structure.

2.2 Weighted transformation approach

GLS can be understood through a transformation that converts the correlated-error problem into an equivalent one with uncorrelated, homoskedastic errors. If \(\Sigma\) is positive definite, one can use a matrix \(A\) such that \(A^\top A=\Sigma^{-1}\) (equivalently, \(A\) is a square-root factor of \(\Sigma^{-1}\)). Then premultiplying both sides of the model by \(A\) yields \[ Ay = AX\beta + A\varepsilon, \] with \[ \mathrm{Var}(A\varepsilon)=A\Sigma A^\top=I \] under the model. In this transformed space, ordinary least squares corresponds to a covariance-aware criterion in the original variables.

2.3 GLS estimator and its properties

The GLS estimator is \[ \hat{\beta}_{\text{GLS}}=(X^\top \Sigma^{-1}X)^{-1}X^\top \Sigma^{-1}y, \] assuming the inverse exists. Under the linear model with \(E[\varepsilon]=0\) and a correctly specified covariance matrix \(\Sigma\), GLS is unbiased. Under standard regularity conditions, it is also efficient in the sense that, among linear unbiased estimators, it has minimum variance. Its covariance matrix is \[ \mathrm{Var}(\hat{\beta}_{\text{GLS}})=(X^\top \Sigma^{-1}X)^{-1}. \] These properties motivate using GLS when \(\Sigma\) is known or can be reliably modeled.

2.4 Connection to efficient estimation

Efficiency in GLS stems from the way \(\Sigma^{-1}\) reweights the estimation criterion according to the noise level and dependence structure. Observations with higher variance or stronger correlation are downweighted relative to those carrying more precise information. This yields tighter coefficient estimates when the covariance model is accurate. In many applications, the practical challenge is not the algebra of GLS, but obtaining a stable and credible estimate of \(\Sigma\) or a parameterization that captures the dependence.

3 Assumptions and model structure

3.1 Linear model with correlated errors

GLS is formulated for the linear model \[ y = X\beta + \varepsilon, \] where correlations are permitted through \(\mathrm{Var}(\varepsilon)=\Sigma\). The matrix \(\Sigma\) typically must be symmetric and positive definite (or at least positive semidefinite with appropriate treatment). Correlation across observations means that residuals cannot be treated as independent random variables. Instead, their joint variability is modeled through \(\Sigma\), and GLS uses \(\Sigma^{-1}\) in the estimator.

3.2 Homoskedastic vs. heteroskedastic cases

In homoskedastic settings, all observations share the same marginal variance, and \(\Sigma\) may still contain off-diagonal terms if errors are correlated. In heteroskedastic settings, variances differ across observations, often producing a diagonal \(\Sigma\) with distinct entries. GLS handles both situations naturally: it uses the full covariance structure, not only variances. When \(\Sigma\) is diagonal with entries proportional to error variances, GLS reduces to a weighted least-squares form with inverse-variance weights.

3.3 Correlation across observations

Correlation may arise from temporal dependence (e.g., in time series), spatial proximity (e.g., in geostatistics), or shared random effects within clusters (e.g., repeated measures for the same individual). In each case, the covariance matrix encodes how errors co-move. Accurate modeling of \(\Sigma\) can be crucial because the estimator depends on \(\Sigma^{-1}\). Misspecification of dependence can degrade efficiency and distort inference, particularly in small samples or when correlation is strong.

3.4 Conditions for unbiasedness and consistency

Unbiasedness of GLS hinges on the mean model being correct: \(E[\varepsilon]=0\) and \(X\) being treated as fixed (or conditionally exogenous) in the relevant probability space. If \(\Sigma\) is known and correctly specified, GLS is unbiased under the standard linear model. Consistency generally requires that the model for the mean is properly specified and that covariance modeling does not introduce systematic bias in the estimator as sample size grows. When \(\Sigma\) is estimated (as in feasible GLS), additional assumptions about the quality and convergence of the covariance estimates are needed.

4 Computation of GLS

4.1 Using the covariance matrix Σ

When \(\Sigma\) is available, the computational task is to evaluate \(\Sigma^{-1}\) and solve \[ (X^\top \Sigma^{-1}X)\hat{\beta} = X^\top \Sigma^{-1}y. \] Direct inversion of \(\Sigma\) is often avoided for numerical reasons. Practical implementations typically compute matrix products involving \(\Sigma^{-1}\) without explicitly forming the inverse, using linear solvers or decompositions. Computational cost depends on the size of \(X\) and the dimensionality of \(\Sigma\), with large covariance matrices often requiring structured covariance models.

4.2 Cholesky and matrix square-root methods

A common method exploits the fact that \(\Sigma\) is symmetric positive definite. One may compute a Cholesky decomposition \[ \Sigma = LL^\top, \] then use it to apply \(\Sigma^{-1}\) through triangular solves rather than explicit inversion. This also provides a square-root factor for the whitening transformation. The approach improves stability because triangular solves are typically more reliable than naive inversion. When \(\Sigma\) is only positive semidefinite, more careful handling (such as generalized inverses or regularization) may be required.

4.3 Solving the GLS normal equations

After obtaining a way to compute \(\Sigma^{-1}X\) and \(\Sigma^{-1}y\), the GLS normal equations can be solved by linear algebra routines. If \(X^\top \Sigma^{-1}X\) is well-conditioned, direct solvers are efficient. In large problems, iterative methods may be used to solve the linear system while avoiding formation of dense intermediate matrices. The choice of solver and representation (dense versus sparse, structured versus unstructured) strongly affects runtime.

4.4 Numerical stability considerations

Numerical stability matters because GLS depends on \(\Sigma^{-1}\), which can amplify noise when \(\Sigma\) is ill-conditioned. Diagnostics such as condition numbers, checks for positive definiteness, and careful scaling of variables can mitigate issues. Regularization of covariance estimates or use of constrained covariance structures is sometimes adopted to prevent unstable inversion. Additionally, centering and scaling \(X\) may improve conditioning of the weighted normal equations.

5 Relationship to other estimation methods

5.1 GLS as weighted least squares

GLS is frequently presented as a weighted least-squares problem. In the special case where \(\Sigma\) is diagonal with entries \(\sigma_i^2\), the estimator becomes \[ \hat{\beta}=(X^\top W X)^{-1}X^\top W y \] with \(W=\mathrm{diag}(1/\sigma_i^2)\). Here each observation is weighted by the inverse of its error variance, aligning the objective with heteroskedasticity. If \(\Sigma\) has correlations, the weights are no longer simple scalar factors per observation, but the same principle—downweighting less reliable information—still underlies the computation through \(\Sigma^{-1}\).

5.2 GLS and maximum likelihood (Gaussian errors)

Under Gaussian errors, GLS connects directly to maximum likelihood estimation. If \(\varepsilon \sim \mathcal{N}(0,\Sigma)\) and the mean model is linear in parameters, the log-likelihood leads to the same estimator as GLS when \(\Sigma\) is treated as known. Thus, GLS can be seen as the maximum likelihood (or equivalently, least-squares) estimator for the regression coefficients in the Gaussian covariance model. If \(\Sigma\) is unknown, maximum likelihood may require joint estimation of covariance and regression parameters.

5.3 GLS vs. feasible GLS

In many real datasets, \(\Sigma\) is not known exactly. Feasible GLS replaces the unknown covariance matrix with an estimate \(\hat{\Sigma}\). This yields an estimator \[ \hat{\beta}_{\text{FGLS}}=(X^\top \hat{\Sigma}^{-1}X)^{-1}X^\top \hat{\Sigma}^{-1}y. \] FGLS retains much of GLS’s structure but inherits additional uncertainty from estimating \(\Sigma\). Under suitable regularity conditions, FGLS can be consistent and asymptotically efficient, but finite-sample performance depends on how accurately the covariance model is specified and estimated.

5.4 Special cases and reductions to OLS

When \(\Sigma\) is proportional to the identity matrix, \(\Sigma=\sigma^2 I\), GLS reduces to OLS because \(\Sigma^{-1}\) is a scalar multiple of \(I\). More generally, if \(\Sigma\) is correctly modeled as diagonal with equal entries, weighting has no effect on the coefficient estimates. These reductions provide a useful check: implementing GLS should reproduce OLS results in the homoskedastic independent-error limit.

6 Feasible GLS (FGLS)

6.1 Estimating the covariance structure

FGLS begins by specifying a form for \(\Sigma\), often via parameterized covariance models such as autoregressive structures for time series, distance-based correlation functions for spatial data, or variance components for clustered designs. Estimation may proceed by fitting residual-based methods from an initial model, estimating heteroskedastic variances, or jointly estimating covariance parameters using likelihood or method-of-moments techniques. The resulting \(\hat{\Sigma}\) is then used in the GLS formula.

6.2 Two-stage estimation workflows

A standard two-stage workflow uses OLS residuals to estimate covariance parameters. First, fit the regression by OLS to obtain residuals \(\hat{\varepsilon}\). Second, use these residuals to estimate \(\Sigma\) (or its parameters) and compute FGLS coefficients using \(\hat{\Sigma}\). This approach is conceptually straightforward but depends on the quality of the initial residuals; if the mean model is poor, covariance estimation can be contaminated.

6.3 When covariance estimates are poor

If \(\hat{\Sigma}\) is substantially misspecified—through incorrect functional form, unstable parameter estimates, or insufficient data—FGLS can perform worse than OLS. Because FGLS uses \(\hat{\Sigma}^{-1}\), errors in covariance estimation may lead to overconfident inference or biased standard errors. Practical remedies include using more robust covariance modeling, imposing structure consistent with the data-generating process, or switching to alternative inference methods when covariance modeling is uncertain.

6.4 Practical diagnostics for FGLS

Diagnostics often include examining residual autocorrelation (for time series), residual spatial dependence (for spatial models), or within-cluster correlation patterns. After estimating \(\hat{\Sigma}\), one may check whether standardized residuals behave closer to independence than under the OLS fit. Sensitivity analysis—recomputing results under alternative plausible covariance models—can help assess how dependent conclusions are on the covariance specification.

7 Inference in GLS frameworks

7.1 Variance of the GLS estimator

When \(\Sigma\) is known, inference is based on the known sampling variability of \(\hat{\beta}_{\text{GLS}}\): \[ \mathrm{Var}(\hat{\beta}_{\text{GLS}})=(X^\top \Sigma^{-1}X)^{-1}. \] Standard errors follow from the diagonal entries of this covariance matrix. In FGLS, the variance formula must account for estimation of \(\Sigma\); a naive substitution of \(\hat{\Sigma}\) can be inadequate in small samples, leading to underestimation of uncertainty.

7.2 Hypothesis testing (Wald-type concepts)

Hypothesis tests in GLS commonly use Wald-type statistics that compare linear restrictions on \(\beta\). For a restriction \(R\beta=r\), one can use \[ (R\hat{\beta}-r)^\top \left[R\,\widehat{\mathrm{Var}}(\hat{\beta})\,R^\top\right]^{-1}(R\hat{\beta}-r), \] with reference distributions depending on assumptions and sample size. These tests rely on correct covariance estimation and appropriate degrees-of-freedom accounting.

7.3 Confidence intervals for regression parameters

Confidence intervals for individual coefficients are constructed from estimated standard errors. Under the Gaussian linear model with known \(\Sigma\), intervals can be computed using normal or exact finite-sample distributions depending on variance knowledge. Under FGLS, intervals must reflect the extra variability from covariance estimation; bootstrap or sandwich-type adjustments are sometimes used as conceptual alternatives when covariance uncertainty is non-negligible.

7.4 Model comparison considerations

Comparing models in GLS contexts may involve likelihood-based criteria when Gaussian assumptions are adopted and covariance parameters are estimated by maximum likelihood. Alternatively, information criteria can compare different mean structures or covariance specifications. However, when covariance modeling is uncertain or misspecified, model comparison can become less reliable, and emphasis may shift to predictive performance or diagnostic-driven model selection.

8 Diagnostics and model checking

8.1 Assessing residual behavior

Residual diagnostics evaluate whether the fitted mean structure is adequate and whether residuals match the assumed covariance pattern. Useful checks include plotting residuals versus fitted values, examining residual distributions, and studying whether remaining dependence persists. In GLS workflows, residuals are often transformed (“whitened”) by the model’s covariance structure; approximate whiteness is an indicator that the covariance model is capturing dependence.

8.2 Checking covariance assumptions

Covariance assumptions can be assessed using residual-based methods consistent with the hypothesized form of \(\Sigma\). For time dependence, plots of residual autocorrelation can highlight mismatches. For spatial dependence, variograms or correlation decay curves can be inspected. For clustered designs, within-group residual correlations can be compared to what the model predicts. These checks help determine whether \(\Sigma\) needs refitting.

8.3 Sensitivity to covariance misspecification

Even when regression coefficients appear stable, misspecifying \(\Sigma\) can affect standard errors and test statistics. Sensitivity analysis can involve varying covariance model parameters, using alternative covariance structures, or comparing GLS results to robust methods. If conclusions change markedly across reasonable covariance specifications, that instability signals that inference is closely tied to covariance modeling assumptions.

8.4 Robust alternatives (conceptual overview)

When covariance modeling is uncertain, robust alternatives aim to protect inference against misspecification of \(\Sigma\). Conceptually, this includes using variance estimators that remain valid under weaker error structures, often producing “sandwich” covariance estimates. Another approach is to model dependence more flexibly (e.g., via structured semiparametric covariance models) or to base analysis on prediction error rather than strict covariance-dependent inference.

9 Applications in statistics and data analysis

9.1 Time series regression with correlated errors

Time series data frequently exhibit autocorrelated errors due to unobserved dynamics, inertia, or seasonality. GLS addresses this by using a covariance model aligned with temporal structure, such as autoregressive or moving-average–type covariance. By accounting for dependence across time points, GLS can improve coefficient estimation and produce more reliable standard errors than OLS.

9.2 Spatial modeling with dependence

In spatial regression, measurements from nearby locations tend to be more similar than those far apart, creating spatial correlation in errors. GLS can incorporate a spatial covariance model where covariance decreases with distance or depends on spatial directions. This yields coefficient estimates that properly reflect the effective information content of spatial observations.

9.3 Panel and clustered data

Panel data combine cross-sectional and time dimensions, often with correlations within individuals or groups. Clustered or hierarchical data have dependence induced by shared environments, shared measurement processes, or latent random effects. GLS can be applied using covariance structures that represent within-cluster correlation, such as compound symmetry or structured correlation by time within cluster.

9.4 Growth curve and repeated-measures settings

Repeated measurements on the same subject or experimental unit naturally create within-unit dependence. Growth curve models often describe trajectories over time with correlated residuals. GLS can be used when the covariance structure of residual variation across measurement occasions is modeled explicitly, enabling efficient estimation of trajectory parameters while accounting for the non-independence of repeated observations.

10 Extensions and advanced topics

Generalized linear models extend linear regression to non-Gaussian outcomes and nonlinear mean structures. While the classical GLS formulation is for linear models with correlated errors, analogous ideas appear in iterative estimation for correlated responses. In practice, methods that incorporate covariance structures into estimating equations can be viewed as “GLS-like” extensions when dependence and heteroskedasticity are present.

10.2 Nonlinear models and iterative GLS approaches

For nonlinear regression, the mean model is not linear in \(\beta\). Iterative procedures can linearize the model around current parameter estimates (e.g., via a Taylor expansion) and then apply GLS to the resulting approximate linear step. This yields iterative weighted estimation schemes that update both parameters and, sometimes, the covariance model.

10.3 Mixed models connection (random effects perspective)

Mixed-effects models represent dependence through random effects that induce correlation among observations. Under certain covariance structures, GLS estimation of fixed effects can correspond to marginal or conditional estimators derived from random-effects formulations. The connection is useful conceptually: GLS treats dependence as given through \(\Sigma\), while mixed models generate \(\Sigma\) from latent components.

10.4 GLS under missing or irregularly spaced data

Irregular spacing and missing observations complicate the definition of the covariance matrix and the alignment of covariance structure with observed times or locations. GLS can still be applied when the covariance model handles irregular observation patterns, often by computing covariances conditional on which time points or locations are present. In some settings, structured covariance models are adapted to work with incomplete schedules while preserving positive definiteness.