Climate change detection and attribution play a central role in establishing the causal influence of human activities on global warming. The dominant framework, optimal fingerprinting, is a linear errors-in-variables model in which each covariate is subject to measurement error with covariance proportional to that of the regression error. The reliability of such analyses depends critically on accurate inference of the regression coefficients. The optimal weight matrix for estimating these coefficients is the precision matrix of the regression error, which is typically unknown and must be estimated from climate model simulations. However, existing regularized optimal fingerprinting approaches often yield underestimated uncertainties and overly narrow confidence intervals that fail to attain nominal coverage, thereby compromising the reliability of analysis. In this paper, we first propose consistent variance estimators for the regression coefficients within the class of linear shrinkage weight matrices, addressing undercoverage in conventional methods. Building on this, we derive a linearly optimal weight matrix that directly minimizes the asymptotic variances of the estimated scaling factors. Numerical studies confirm improved empirical coverage and shorter interval lengths. When applied to annual mean temperature data, the proposed method produces narrower, more reliable intervals and provides new insights into detection and attribution across different regions. Successive assessments by the Intergovernmental Panel on Climate Change (IPCC) have firmly established that more than half of the observed increase in global average surface temperature in recent decades can be attributed to anthropogenic increases in greenhouse gas concentrations and other human-induced forcings (Hegerl et al. 2007, Bindoff et al. 2013, Eyring et al. 2021). Detection and attribution (D&A) analyses have played a central role in reaching these conclusions. In climate science, detection refers to the process of demonstrating that a climate variable has changed in a statistically significant manner, without necessarily identifying the cause of that change. Attribution, on the other hand, involves assessing the extent to which observed changes can be attributed to multiple external forcings, along with an assignment of statistical confidence (e.g., Hegerl & Zwiers 2011). By comparing climate model simulations with observed climate variables, detection and attribution analyses evaluate whether observed changes are statistically consistent with expected responses, also known as fingerprints or signals, to one or more external forcings in the climate system. Optimal fingerprinting (OF), the most widely used method in detection and attribution analyses, is a multiple linear regression framework in which observed climate variables are regressed onto the fingerprints of external forcings (e.g., Hegerl et al. 1996, Allen & Tett 1999, Allen & Stott 2003). The primary target of statistical inference in OF is the vector of regression coefficients, commonly referred to as scaling factors. These coefficients scale the fingerprints to match the observed climate changes best. An ideal point estimator of the scaling factors should be unbiased and exhibit minimal variance. Confidence intervals constructed around these estimates quantify uncertainty and form the basis for drawing detection and attribution conclusions. Specifically, if the confidence interval for a scaling factor lies significantly above zero, the effect of the corresponding external forcing is said to be "detected" in the observational data. If, in addition, the interval contains one, then there is sufficient statistical evidence to "attribute" the observed changes to that external forcing. Additionally, a proper confidence interval should have a coverage rate matching the nominal confidence level to ensure reliable and statistically robust conclusions. Historically, OF was deemed "optimal" in the context of generalized least squares (GLS), where the precision matrix of the regression error is used as a weight for prewhitening. This approach yields estimators of the scaling factors with minimum variance under idealized assumptions. However, it was later recognized that fingerprints are not directly observed but are instead estimated from climate model simulations, thereby introducing measurement errors and giving rise to an errors-in-variables (EIV) model. Under the standard assumption that internal climate variability in model simulations mirrors that in observations, the errors in the estimated fingerprints inherit the same covariance structure as the regression errors. As a result, the estimation framework shifted toward total least squares (TLS) (Allen & Stott 2003) with both the response and covariates are "prewhitened" using the covariance matrix  of the regression errors. In practice,  is typically unknown and must be estimated from ensembles of climate model simulations (Allen & Stott 2003, Ribes et al. 2013). This estimation step is generally handled separately from the regression analysis and treated as a preliminary step (e.g., Hannart et al. 2014). Estimating  poses significant challenges due to spatial and temporal dependencies and high dimensionality of the climate variables. The number of available control runs from climate model simulations is often insufficient to yield a reliable estimate of . In particular, when the number of control runs is smaller than the dimension, the sample covariance matrix becomes singular and cannot be directly inverted to construct a weight matrix. Early methods addressed this issue by projecting the data onto the leading empirical orthogonal functions (EOFs) of internal climate variability, as represented by the empirical covariance matrix of the control runs (Hegerl et al. 1996, Allen & Tett 1999). Later, Ribes et al. (2009) proposed using a linear shrinkage estimator of  developed by Ledoit & Wolf (2004), leading to the regularized optimal fingerprinting (ROF) method. Confidence intervals for the scaling factors in ROF can be constructed using normal approximation techniques (Ribes et al. 2013, DelSole et al. 2019, Li et al. 2021) or bootstrap methods (DelSole et al. 2019). Numerical studies by Ribes et al. (2013) demonstrated that ROF yields a more robust and accurate implementation of optimal fingerprinting than the EOF-based approach, primarily because the performance of EOF methods is highly sensitive to the selection of the number of retained EOFs. To address the limitations of EOF truncation, Katzfuss et al. (2017) introduced a Bayesian framework that treats the number of retained EOFs as a parameter, allowing robust inference through averaging over different truncation levels.
However, the Bayesian approach can be computationally intensive, particularly when the maximum number of truncations is large, and its performance can be sensitive to the choice of prior distributions. As an alternative, Hannart (2016) proposed an integrated likelihood formulation that derives a closed-form joint likelihood of the observational data and control runs by integrating out the unknown covariance matrix. This integrated approach is equivalent to a Bayesian model with an informative conjugate prior on the covariance matrix. It also enables shrinkage toward structured targets beyond the identity matrix, such as spatio-temporal covariance structures, as explored in the numerical studies of Hannart (2016). In practice, however, the prior structure on the covariance matrix may be uncertain or misspecified, potentially limiting the reliability of the resulting inference. The use of an estimated covariance matrix  has significant implications for optimal fingerprinting, especially when the estimation is based on relatively small samples. Recent studies have shown that the resulting scaling factor estimator under existing regularized optimal fingerprinting (ROF) is no longer optimal in terms of mean squared error (MSE) when  is estimated with substantial uncertainty (Li et al. 2023). In such settings, alternative weight matrices may yield more accurate estimators of the scaling factors in terms of MSE. Moreover, confidence intervals constructed using standard practices, such as asymptotic normal approximations based on two independent samples (Hegerl et al. 1996, Allen & Stott 2003) or bootstrap techniques (DelSole et al. 2019), often inadequately account for the uncertainty in the estimated . As a result, the coverage rates of these intervals are typically lower than the nominal level. To address this issue, Li et al. (2021) proposed a parametric bootstrap calibration method that inflates the confidence intervals to achieve nominal coverage. However, this method is computationally intensive and may perform poorly when the sample size for estimating  is limited. An alternative approach was developed by Ma et al. (2023), who proposed an estimation procedure based on estimating equations (EE). This method produces unbiased estimators that correct for the bias induced by EIV and constructs weights using the known structure of the error covariance matrix. The EE approach assumes temporal stationarity in climate variability, which implies a block Toeplitz structure for the covariance matrix, and uses a pseudobootstrap algorithm for constructing confidence intervals. This method improves both the efficiency of the estimator and the accuracy of the coverage rates. However, it may lose some efficiency in point estimation due to not directly estimating the temporal correlation and relying on a suboptimal weight matrix. In contrast, traditional fingerprinting methods make use of the full spatio-temporal covariance structure, which may maximize estimation efficiency but at the cost of large uncertainty in estimating the high dimensional covariance matrix, thereby compromising optimality. To avoid confusion, we henceforth distinguish regularized optimal fingerprinting (ROF) from regularized fingerprinting (RF) for clarity in the sequel. This distinction raises a natural and important question: Can optimality in terms of MSE be recovered for scaling factor estimators in regularized fingerprinting, at least within a suitably defined class of covariance estimators, when the error covariance matrix  is fully estimated? In this paper, we tackle the challenging problems in the regularized fingerprinting framework based on a class of linear shrinkage estimators towards an identity matrix of the form () = S + I, where S is the sample covariance matrix from controlled climate simulations and  > 0 is a tuning parameter. The contributions of our framework are twofold: First, we establish the asymptotic properties of the resulting RF estimators of the scaling factors when the weight matrix is constructed from (), under suitable regularity conditions. We then propose a data-driven procedure to consistently estimate the asymptotic covariance matrix of the scaling factor estimators, using only the estimated fingerprints and the sample covariance matrix. This directly addresses the undercoverage issue frequently encountered in regularized fingerprinting, where the uncertainty introduced by estimating the error covariance matrix is not fully accounted for, resulting in overly narrow confidence intervals. Second, we introduce a computational procedure for selecting the optimal regularization parameter  via grid search, aiming to minimize the total uncertainty in the scaling factor estimates. The method achieves asymptotic optimality within the class of linear shrinkage estimators, as both the sample size and the matrix dimension tend to infinity at a fixed ratio, potentially less than one. We refer to the proposed approach as regularized fingerprinting with linear optimality. Extensive simulation studies under realistic conditions demonstrate that our method provides accurate uncertainty quantification and leads to confidence intervals with empirical coverage rates close to the nominal level. Moreover, the optimally selected weight matrix produces substantially shorter intervals than existing approaches. These improvements have important practical implications, as precise uncertainty quantification is central to robust detection and attribution. In a real-world application to the detection and attribution of global near-surface air temperature over the period 1951-2020, our method consistently outperforms competing approaches in terms of interval length and inferential reliability, and provides new insights into detection and attribution conclusions across different regional scales. It is worth noting that, although developed in the context of climate detection and attribution, the proposed methodology is broadly applicable to multivariate errors-in-variables regression problems where the covariates are contaminated by structured noise and the error covariance must be estimated from limited auxiliary data. The rest of the paper is organized as follows. In Section 2, we briefly review the OF framework and the asymptotic properties of the class of linear shrinkage estimators of , then propose a consistent estimate of the asymptotic covariance and an optimal weight matrix that minimizes the total variances of the resulting scaling factors. A large-scale numerical study assessing the performance of the proposed method is reported in Section 3. In Section 4, we apply the proposed method to a detection and attribution analysis of changes in mean near-surface temperatures on continental and subcontinental scales. A discussion concludes in Section 5.
To improve readability, we relegate technical details, including proofs of the theoretical results, additional results from simulation studies, and details of the climate models, to the Supplementary Materials. Fingerprinting can be formulated as a linear regression problem with errors-in-variables (EIV): where Y  R N is the observed climate variable of interest, X i  R N is the true but unobserved fingerprint of the ith external forcing, and  i is the associated scaling factor. The error term  is assumed to follow a multivariate normal distribution with mean zero and covariance matrix . For each i, Xik is the kth simulation of the fingerprint X i from a climate model, contaminated by model-specific internal variability  ik . The number of ensemble members for the ith forcing is denoted by n i . Under the standard assumption that the internal variability in model simulations is consistent with that in observations, the errors  ik are modeled as independent Gaussian noise with covariance , i.e., Moreover,  is independent of  ik . Both the response Y and the covariates are assumed to be centered with respect to a common reference period, so no intercept is included in the regression model. Naturally, the ensemble mean Xik serves as an estimator for the true fingerprint X i . For notational convenience, define X = ( X1 , . . . , Xp ), X = (X 1 , . . . , X p ), D = diag(1/n 1 , . . . , 1/n p ). The primary objective is to estimate the scaling factors  = ( 1 , . . . ,  p )   R p and to construct reliable confidence intervals for inference, given the observed response vector Y , the ensemble-averaged fingerprints X, and the scaling matrix D. Historically, the term "optimal" in optimal fingerprinting originates from early formulations that relied on two idealized assumptions: (i) the error covariance matrix  is known, and (ii) the fingerprint matrix X is fully observed. Under these conditions, the frameworks of weighted linear regression and generalized least squares (GLS) motivate a prewhitening procedure, in which both Y and X are premultiplied by  1/2 . Then, the GLS estimator of  with weight matrix  is given by GLS The optimality of GLS follows from the fact that it is the best linear unbiased estimator (BLUE) of  under the classical Gauss-Markov assumptions (Chen et al. 2024). Subsequent developments recognized that the true fingerprints X are not directly observed; instead, only their estimates X are available. It has been shown that substituting X for X in the regression induces bias in the estimation of  (Allen & Stott 2003). Given knowledge of , and under the assumption that the simulation error  ik and the observational error  share the same covariance structure, both Y and X can be prewhitened accurately. In this setting, the method of total least squares (TLS) can be applied to the prewhitened data. Specifically, the TLS estimator of  with weight matrix  is given by TLS where    2 denotes the Euclidean norm. For further details, see Gleser (1981) and Allen & Stott (2003). Other approaches, such as the estimating equation estimator proposed by Ma et al. (2023), have also been developed, but are beyond the scope of the present work. Optimal fingerprinting assumes that the error covariance matrix  is known; however, in practice, it must be estimated. The standard approach proceeds in two steps. In the first step, preindustrial control runs from climate model simulations are used to estimate . These control runs are generated under the assumption that the models capture only the internal variability of the climate system, without the influence of external forcings. Specifically, suppose we have m independent, centered control runs Z 1 , . . . , Z m  R N . Under this modeling assumption, An estimator of  is then constructed from Z 1 , . . . , Z m . In the second step, the estimated covariance matrix is used to prewhiten both the response Y and the covariates X in the regression model. The TLS estimator of  is then computed from the prewhitened data. Nonetheless, consistently estimating  in the first step is a challenging task. The traditional estimator is the sample covariance matrix, Since  is an N xN matrix, it has N (N +1)/2 free parameters when no structural assumptions are imposed. This number is typically too large relative to the available sample size m, which in practice is often at most a few hundred. From a theoretical perspective, in the high-dimensional asymptotic regime where both N and m grow and N is comparable to or exceeds m, the sample covariance matrix S is known to be inconsistent (Silverstein & Bai 1995). More severely, when N > m, S is singular and cannot be directly used as a weight matrix in regression. As a consequence, the optimality of the original fingerprinting estimator, which relies on the true , is unlikely to hold when  is replaced by S. In such settings, the statistical properties of the resulting estimator of  are fundamentally governed by the behavior of the estimated weight matrix. To address the challenge of estimating  in high dimensions, Ribes et al. (2009) introduced the linear shrinkage framework developed by Ledoit & Wolf (2004) into the context of fingerprinting. Specifically, they considered an estimator of  of the form where the scalar tuning parameters  and  are chosen to minimize the expected mean squared error between the estimator and the true covariance matrix . Numerical studies by Ribes et al. (2009) demonstrated that this regularized estimator yields robust and stable results, particularly in settings where the sample size is limited relative to the dimension. The linear shrinkage family can be viewed as an application of classical ridge regularization in the context of fingerprinting. This approach has also been employed in high-dimensional linear hypothesis testing, as in Li, Aue & Paul (2020) and Li, Aue, Paul, Peng & Wang (2020). However, recent findings by Li et al.
(2023) have shown that the selection of  and  in the shrinkage estimator is not optimal with respect to MSE between TLS and the true scaling factor , particularly when the number of control runs m is relatively small. Moreover, confidence intervals constructed using either normal approximation or bootstrap techniques tend to exhibit low empirical coverage rates, often falling short of the nominal level. The current work aims to recover the optimality in terms of MSE of the TLS estimator TLS within the family of linearly regularized covariance estimators. In addition, we aim to develop reliable and computationally efficient procedures for constructing confidence intervals for the scaling factors, based on the optimally regularized estimator. In particular, our research is formulated as follows. Consider the family of linearly regularized estimators of  of the form where  > 0 is a regularization parameter. The total least squares (TLS) estimator of  based on () is defined as Notably, in contrast to Ribes et al. (2009), we fix  = 1 in the shrinkage formulation, as the TLS estimator TLS with weight matrix () is invariant to a positive scalar multiplication of (); that is, for any a > 0, replacing () with a () yields the same estimate. Our objective is to find the value of  that minimizes the asymptotic MSE of (), i.e. under the high-dimensional asymptotic regime where N and m grow and the ratio N/m  c for some c > 0. In addition, we aim to provide a reliable and computationally efficient estimator of the asymptotic covariance matrix of ( opt ). Based on this estimator, confidence intervals for the scaling factors  can be constructed in a principled and data-driven manner. The following assumptions are imposed for the analysis of the asymptotic properties of (). Assumption 1. There exist constants  and  such that 0 <  < lim inf N  l min ()  lim sup N  l max () < , where l min () and l max () are the smallest and largest eigenvalue of a matrix, respectively. Assumption 2. Assume that lim N  X T X/N exists and is a positive definite matrix. Assumption 3. For any fixed , lim N,m X T 1 ()X/N =  1 () exists, where  1 () is a nonsingular matrix. Assumption 4. For any fixed , lim N,m Theorem 1. Under Assumptions 1-6, as N, m   with N/m  c for some c > 0, where Note that the asymptotic MSE of () satisfies The optimal regularization parameter  opt is then the minimizer of tr(()). We now propose a consistent estimator of () that depends only on the simulated fingerprint matrix X and the sample covariance matrix S. where Further define Henceforth, o  (1) denotes a matrix A such that A F = o P (1), where    F denotes the Frobenius norm of a matrix. The following result is proved in Lemma 2 of Chen et al. (2011). Proposition 4. Under Assumptions 1-6, as N, m   with N/m  c for some c > 0, for any fixed . The estimators of the key parameters are then constructed as: Lemma 1. Under Assumptions 1-6, as N, m   with N/m  c for some c > 0, for any In practice, the empirical optimal regularization parameter is selected by where  and  are prespecified lower and upper search bounds, respectively. Given the eigen-decomposition of S, the computational complexity of evaluating 1 (), 2 (), and K() for each  is O(pN ). This allows the optimization over  to be performed efficiently via a grid search. In finite-sample settings, the quantities Q 1 () and Q 2 () may become unstable when N > n and  is close to zero, due to the near-singularity of (). To mitigate this issue, it is recommended that the lower bound  of the search range is not too small. A practical choice for the search interval is [0.01 , 10 ], where  = (1/N ) tr(S). Lastly, we propose marginal confidence intervals for each individual scaling factor  i , as well as a joint confidence region for the vector , at the asymptotic confidence level (1  ). Marginal confidence interval for  i : where z 1/2 denotes the upper (1  /2)-quantile of the standard normal distribution, and ii ( opt ) is the ith diagonal element of the estimated asymptotic covariance matrix. Joint confidence region for : where  2 p (1  ) is the upper (1  )-quantile of the  2 -distribution with p degrees of freedom. To evaluate the finite sample performance of the proposed method in comparison with existing practices in Regularized Fingerprinting, we conducted extensive simulation studies emulating realistic settings for detection and attribution analyses of global mean temperature changes, following the settings of Li et al. (2023). The climate variable of interest consisted of 11 decadal near-surface mean temperatures over 25 spatial grid boxes, resulting in a response vector of dimension N = 275. In each setting, we first set the true fingerprints X. Two external forcings were considered, anthropogenic (ANT) and natural (NAT) forcings. The expected fingerprints of these two forcings, denoted by X 1 and X 2 , respectively, were set to the average of all runs from the CNRM-CM5 model simulations, as in Ribes et al. (2013) and adopted in Li et al. (2023). To vary the strength of the signals, each X i was scaled by   1, 0.5, representing two signal-to-noise regimes. The case of  = 1 mimics a global scale study with strong signal strength, while  = 0.5 matches regional scale studies where the signals are weaker. The true scaling factors were set to For the true covariance matrix , two spatiotemporal structures were evaluated. The first  was set to be an unstructured matrix  U N obtained by manipulating the minimum variance estimate with a set of CNRM-CM5 model simulations considered, as used in Li et al. (2023). The resulting covariance structure resembles the pattern of an unstructured spatial-temporal covariance matrix with variance stationarity and weak dependence over the time dimension considered by (Hannart 2016).
The second structure, denoted as  ST , was set to be a separable spatiotemporal covariance matrix, where the diagonals were set to be the sample variances from the climate model simulations without imposing temporal stationarity, and the corresponding correlation matrix was set to be the Kronecker product of a spatial correlation matrix and a temporal correlation matrix, both with autoregressive of order 1 and coefficient 0.1. With the X i ,  i , , observed responses Y and noisy fingerprints ( X1 , X2 ) were generated from models (1) and ( 2). The regression errors  followed a multivariate normal distribution N (0, ). The distribution of the measurement error  i for Xi , i  {1, 2}, was N (0, n 1 i ) with (n 1 , n 2 ) = (35, 46), consistent with the number of simulation runs in a two-way detection and attribution analysis of the annual mean temperature conducted by Ribes et al. (2013). Control runs Z 1 , . . . , Z m were generated independently from N (0, ) with sample size m  {50, 100, 200, 400}. Here m = 50 is typical in OF studies, and m  200 is possible but not easily obtained unless runs from different climate models are pooled, ignoring the model structure differences. For each combination of ,  and m, we performed 1000 simulation replicates to evaluate the performance of the proposed method with optimally selected tuning parameter , denoted as "Optim", in comparison with two existing ROF methods based on TLS. The first competitor, denoted as "LS-CB", adopts the ROF method of Allen & Stott (2003) with a linear shrinkage estimator LS by Ledoit & Wolf (2004) for prewhitening (LS) and a calibration bootstrap (CB) for interval estimation (Li et al. 2021). The second, denoted as "MV-CB", applied the minimum variance estimator MV from Li et al. (2023) with the same CB adjustment for constructing confidence intervals. For the proposed method "Optim", we constructed the confidence interval from the newly proposed asymptotic results. Since LS falls into the class of linear shrinkage estimators defined in Section 2, our asymptotic results also apply to this setting. As an additional benchmark, we included an uncalibrated version of the ROF method with LS , denoted as "LS". The corresponding confidence intervals were constructed from the same normal approximation as our proposed method. Confidence intervals are essential in detection and attribution studies. Ideally, they should be as short as possible while maintaining empirical coverage rates close to the nominal level. To assess performance, we focus on two key metrics, empirical coverage rate and interval length. Here, we present in Figure 1 the empirical coverage rates and average lengths of the 95% confidence intervals derived from the four competing methods: Optim, LS, LS-CB, and MV-CB, for the ANT forcing, which is usually the main concern in detection and attribution analyses of climate changes. Full numerical results for both ANT and NAT forcings are provided in Table S1 of the Supplementary Material, which demonstrates that all methods yield unbiased point estimates of the scaling factors. For confidence intervals, our proposed method for estimating the asymptotic variance of the scaling factors estimators leads to empirical coverage rates consistently close to the nominal 95% level across nearly all settings. Even in the most challenging case with only m = 50 control runs, the coverage rate remains around 91%, supporting the accuracy of our estimated asymptotic variance. In contrast, both calibrated competitors (LS-CB and MV-CB) still suffer from undercoverage issues in proposed Optim method offers valid confidence intervals with near-nominal coverage and competitive or superior interval widths across a wide range of realistic scenarios. In addition, it also offers substantial computational advantages over calibration-based approaches, as it estimates the asymptotic variance without requiring any bootstrap procedure. To demonstrate the performance of proposed methods in real-world applications, we conducted a detection and attribution analysis of changes in the mean near-surface air temperature at global (GL), continental and subcontinent scales over the year period 1951-2020, utilizing the latest available climate observations and simulations. Following Zhang et al. (2006) and Li et al. (2023), we considered several regions: at the continental scale, Northern Hemisphere (NH), NH midlatitudes (NHM) between 30  and 70  , Eurasia (EA), and North America (NA); and at the subcontinental scale, Western North America (WNA), Central North America (CNA), and Eastern North America (ENA), where spatio-temporal correlation structures are more likely to hold. Detection and attribution analyses for two external forcings, anthropogenic (ANT) and natural (NAT) forcings, were conducted for each region. For each regional analysis, Models (1)-( 2) require three components: the observed mean temperature Y  R N , the estimated fingerprints XANT and XNAT , and independent control runs Z 1 , . . . , Z m for estimating the covariance matrix . We first obtained the observational vector Y from the latest HadCRUT5 dataset (Morice et al. 2021), which provides monthly anomalies of near-surface air temperature from January 1850 on 5  x 5  grid boxes relative to the 1961-1990 reference period. At each grid box, annual anomalies were computed from monthly values provided that at least nine months of data were available within a given year; otherwise, the annual mean was marked as missing. Nonoverlapping 5-year averages were subsequently calculated, requiring no more than two missing annual values within each 5-year period. After removing the 1961-1965 period due to centering, 13 values of 5-year averages were obtained per grid box. To reduce spatial dimensionality for the global and continental-scale analyses, available 5 In particular, the fingerprints XNAT for the NAT forcing were obtained directly by averaging over 40 available runs. For the ANT forcing, which is typically the primary focus in climate detection and attribution studies, direct simulation outputs were not available from CMIP6 models. Under the linear additivity assumption, X AN T = X GHG + X AER (Zhang et al. 2006), we constructed XANT by combining the greenhouse gas and aerosol fingerprints for each model.
During processing, the same missing data pattern observed in Y was imposed on the fingerprints, and the same averaging and spatial aggregation procedures were applied to reduce dimensionality. Since both the observational and model data were centered relative to the 1961-1990 period, the first 5-year block (1961)(1962)(1963)(1964)(1965) was excluded to maintain consistency. S3 of the Supplementary Material. confidence intervals from different methods. For weight matrix construction, "Optim" denotes the proposed linearly optimal method, "LS" denotes the linear shrinkage estimator and "MV-CB" denotes the minimum variance estimator with parametric calibration for confidence intervals. For the detection and attribution conclusions, results are largely consistent across the three methods for the ANT forcing, with detection achieved in all regional analyses, albeit with slight differences in magnitude. In contrast, results for the NAT forcing show greater variability across methods, owing to the weaker signal associated with natural influences on climate change (Gillett et al. 2021). In particular, for the Optim and LS methods lead to detection and attribution conclusions for the ANT forcing, while the MV-CB method supports detection only. For the NAT forcing, the MV-CB method suggests detection and attribution, whereas the Optim and LS methods do not, as their confidence intervals cover zero. In supercontinent scale analyses, such as GL and NH, MV-CB yields weaker detection statements compared to the other two methods, as the confidence intervals are close to zero. In other regions, the three methods produce consistent conclusions. Although it is not possible to definitively determine which method is more accurate for this single analysis, given the unknown true underlying scaling factors, we emphasize that the comparative performance of these methods has been thoroughly evaluated in the simulation studies. Overall, the proposed Optim method generally provides more reliable results, characterized by lower uncertainty and higher computational efficiency. Optimal fingerprinting, the principal methodological framework for detection and attribution studies in climate change research, has substantial influence on contemporary climate analysis. These analyses provide foundational support for observationally constrained climate projections and facilitate the estimation of critical climate system parameters, such as climate sensitivity. Nonetheless, the original optimality condition of optimal fingerprinting, which involves minimizing total uncertainty in scaling factor estimators, is compromised when estimated fingerprints of external forcings are subject to measurement errors and when the covariance matrix  of regression errors needs to be estimated rather than known as a priori. Furthermore, existing TLS approaches used to estimate scaling factors frequently exhibit undercoverage issues, primarily due to the underestimation of variances of the estimators. Our proposed methodology addresses these limitations and offers improvements in two key aspects. First, within the regularized fingerprinting framework utilizing a linear shrinkage weight matrix, we propose an efficient, data-driven procedure for consistently estimating the asymptotic covariance matrix of the scaling factor estimators. This procedure utilizes only the estimated fingerprints and the sample covariance matrix from control runs, enabling the effective construction of confidence intervals for scaling factors based on normal approximations. Within the class of linear shrinkage estimators, our approach yields valid confidence intervals with close-to-nominal coverage rates at substantially reduced computational cost compared to the existing calibration method (Li et al. 2021). Unlike methods dependent on explicit distributional assumptions (Hannart 2016, Katzfuss et al. 2017) or those requiring temporal stationarity to mitigate undercoverage at the expense of efficiency (Ma et al. 2023), our method imposes no additional assumptions and fully leverages the spatio-temporal covariance structure of , thus potentially achieving optimality. Our second contribution involves determining the optimal weight matrix within the linear shrinkage estimator class. By directly minimizing the asymptotic mean squared error (MSE) of the scaling factor estimator, we identify the optimal shrinkage parameter and construct the corresponding weight matrix by inverting the resulting optimal linear shrinkage estimator. This method is demonstrably more efficient than current regularized fingerprinting practices, as evidenced by both simulation studies and real-world applications. Consequently, the original optimality of fingerprinting is substantially restored in practical terms, reducing uncertainty in key quantities such as attributable warming and climate sensitivity. Overall, the proposed method represents a promising, easy-to-implement detection for practical applications, offering point estimates with lower MSE and confidence intervals with desirable coverage rates. By directly addressing the long-standing undercoverage issue, our methodology improves the reliability of inference. As demonstrated in our application, while qualitative conclusions regarding detection and attribution may remain largely consistent with those from existing methods, revisiting the main results that support attribution statements in IPCC Assessment Reports (Eyring et al. 2021) using our framework is both feasible and worthwhile. A publicly available software implementation further facilitates this task. The proposed method can be extended in several directions. It is particularly interesting to investigate how asymptotic results under conditions where N, m   and N/m  c might inform and refine regularized fingerprinting practices. Adjusting temporal and spatial resolution, which determines N , could significantly affect the efficiency of inference and thus alter the conclusions of detection and attribution analyses. Additionally, our current approach overlooks differences among climate models by treating the runs under each forcing as identical. In practice, these runs typically originate from different climate models. A more realistic modeling approach should explicitly account for the heterogeneity in variability among different models in estimating fingerprints under each external forcing. Furthermore, developing a goodness-of-fit procedure to verify the consistency of variability between observations and climate models would enhance the reliability of detection and attribution analyses. Finally, the proposed method is also applicable beyond climate science. In particular, multivariate errorsin-variables regression problems with high-dimensional structure and limited auxiliary data are common in fields such as genomics, econometrics, and neuroimaging. In these domains, where covariates are measured with structured noise and precise inference is essential, our regularized framework offers a principled approach to improving estimator efficiency and broader applications represent a natural direction for future methodological development and practical deployment.
any m  2, we have where C m is a constant only depending on m. We now consider a modified version of  in which the dependence on a specific Z k is removed, defined as The following result can be shown following similar arguments as those in Lemma 2.10 of Bai & Silverstein (1998). Lemma S4. For any matrix D and  > 0, We summarize several well-known results from Random Matrix Theory (RMT); for detailed derivations and background, see Bai & Silverstein (2010). Let l 1 , . . . , l N denote the eigenvalues of the population covariance matrix , and define its empirical spectral distribution as For any z  C + := {u + iv  C : v > 0}, there exists a unique solution s(z)  C + to the equation commonly referred to as the Marcenko-Pastur equation. Here, c = lim N  N/m  (0, ) denotes the limiting aspect ratio. The function s(z) admits a smooth extension to the negative real for any   R + , the limit exists and defines a real-valued function s() that is analytic on R  . Moreover, s() satisfies the transformed Marcenko-Pastur equation: . Also recall the definition of  1 () and  2 () as Under Assumptions 1-6, if N, m   with N/m  c > 0, then for any fixed  > 0, the quantities Q 1 () and Q 2 () satisfy where s() is the limiting Stieltjes transform defined in (A.3), and s  () denotes its derivative. Define Here, it is known that 1  c {1  s()} is bounded away from zero for all   R + . Then, it follows directly that  1 ()   1 () The following results are established in Lemma 2 of Chen et al. (2011). Note that the second result implies Proposition 4 of the manuscript, since  2 () =  2 () + o p (1). Following the strategy of Lemma 2 in Chen et al. (2011), we can also establish the following analogous result. In the following, we outline the key steps in the proof of Lemma S6, omitting detailed calculations as they closely follow the arguments used in the proof of Lemma 2 in Chen et al. (2011). Consider the identity Multiplying both sides by 2 () and taking the averaged trace, we have Using the Woodbury matrix identity, for appropriate residual terms R (j) , j = 1, 2, 3, 4. Following the strategy of the proof of Lemma 2 in Chen et al. (2011), we can use Lemma S2, Lemma S3, Lemma S4 to show that sup 1j4 R (j) = o p (1). The details are omitted, as they closely mirror those in the referenced proof. It follows then It completes the proof of Lemma S6. Next, we establish a series of analogous results concerning quadratic forms involving () and . Throughout the remainder of this section, let a  R N denote a sequence of deterministic vectors satisfying the normalization condition N 1 a 2 2 = 1. The following three lemmas are shown in El Karoui & Kosters (2011). Combining Lemma S7, Lemma S8, and Lemma S9, and using the relationship between  1 () and  2 (), we obtain The theorem is obtained by applying Theorem 2 of Li et al. (2023). Reparameterize the model as where X * = XD 1/2 and  * = D 1/2 straightforward to verify that all columns of X * have the same covariance matrix . Let  * () denote the TLS estimator of  * using the weight matrix (). Then, it holds that  * () = D 1/2 (). The conditions of Theorem 2 in Li et al. (2023) are satisfied under this reparameterized model. Applying the theorem yields The proof of Theorem 1 is complete by noting that D 1/2  * ()D 1/2 = (). Propositions 2 and 3 are proved by analyzing the bias introduced in the quadratic forms when X is replaced by X. Let W = (W 1 , . . . , W p ) be an N x p matrix with independent and identically distributed N (0, 1) entries, and assume that W is independent of . Then, the ensemble-averaged fingerprint matrix X has the same distribution as In the following analysis, we shall assume without loss X takes the form First, consider G 1 () = XT 1 () X/N . We decompose it as For i = 1, 2, . . . , p, using Lemma S2, Lemma S5 and Equation (A.4), Combining the results, we get Similarly, for any i = 1, 2, . . . , p, All together, we obtain Following analogous arguments as those in the analysis of G 1 (), for i = 1, 2, . . . , p, using Lemma S2, Lemma S5 and Equation (A.4), Combining the results, we get Similarly, for any i = 1, 2, . . . , p, All together, we obtain Combine the decomposition of G 1 () and G 2 () and Equation (A.5). It completes the proof of Proposition 3. The results of simulation studies in Section 3 of the main text are detailed in Table S1. The CMIP6 climate model simulations are detailed in Table S2 and Table S3. Table S1 Scaling factor estimates and 95% confidence interval. Two structures of considered:  UN for proposed shrinkage estimator from ensemble simulations and  ST for separable spatio-temporal covariance matrix. is the scale to control the signal-to-noise ratio. Error terms are generated from multivariate normal distribution. Number of ensembles for two forcings are n 1 = 35 and n 2 = 46. Four methods are compared, Optim, LS, LS-CB, and MV-CB, based on 1000 replicates. Bias and standard deviation (SD x 100) of scaling factors and average length of confidence intervals (CIL) and corresponding empirical rate (CR) from 1000 replicates are recorded.
Optim -0.00 3.2 0.11 90.6 0.00 8.5 0.33 92.6 LS -0.00 3.4 0.12 91.6 -0.00 8.6 0.35 94.2 LS-CB -0.00 3.2 0.11 91.8 0.00 8.7 0.39 96.7 MV-CB 0.00 3.1 0.11 91.0 0.00 8.5 0.35 95.1 100 Optim -0.00 2.1 0.08 95.2 -0.00 6.3 0.27 95.6 LS -0.00 2.2 0.09 94.6 -0.00 6.5 0.28 96.4 LS-CB -0.00 2.2 0.10 96.7 0.00 7.1 0.30 95.6 MV-CB -0.00 2.0 0.08 95.6 0.00 7.0 0.26 91.6 200 Optim -0.00 1.6 0.07 94.4 -0.00 5.1 0.22 94.0 LS -0.00 1.7 0.07 95.4 -0.00 5.2 0.23 95.4 LS-CB -0.00 1.6 0.08 97.3 0.00 5.9 0.25 96.5 MV-CB -0.00 1.6 0.06 95.4 0.00 5.6 0.21 93.7 400 Optim -0.00 1.3 0.05 96.6 -0.00 4.4 0.18 94.2 LS -0.00 1.4 0.06 96.2 -0.00 4.5 0.19 95.8 LS-CB -0.00 1.5 0.06 96.5 0.00 5.1 0.21 93.7 MV-CB -0.00 1.4 0.05 93.5 -0.00 4.8 0.18 90.7  Setting 2:  ST ; SNR  = 1 50 Optim -0.00 1.4 0.06 92.0 0.00 4.2 0.16 91.6 LS 0.00 1.5 0.06 93.0 0.00 4.3 0.16 93.0 LS-CB -0.00 1.6 0.06 90.5 0.00 4.4 0.18 94.8 MV-CB -0.00 1.6 0.05 89.6 0.00 4.5 0.16 91.6 100 Optim 0.00 1.1 0.04 94.6 0.00 3.5 0.13 94.4 LS 0.00 1.1 0.05 94.6 0.00 3.6 0.14 95.0 LS-CB -0.00 1.1 0.05 97.3 0.00 3.2 0.14 97.8 MV-CB -0.00 1.0 0.04 95.4 0.00 3.2 0.12 95.9 200 Optim 0.00 0.8 0.03 96.0 0.00 2.7 0.11 95.4 LS 0.00 0.8 0.03 95.4 0.00 2.8 0.11 95.0 LS-CB -0.00 0.8 0.04 97.0 0.00 2.8 0.12 96.7 MV-CB 0.00 0.7 0.03 95.4 0.00 2.7 0.10 94.6 400 Optim -0.00 0.7 0.03 95.2 0.00 2.3 0.09 96.8 LS -0.00 0.7 0.03 95.4 0.00 2.4 0.10 95.2 LS-CB 0.00 0.7 0.03 97.5 0.00 2.3 0.10 97.8 MV-CB 0.00 0.6 0.03 95.4 0.00 2.3 0.09 94.6 We collect several technical lemmas in this section. In what follows,    2 denotes the operator norm (spectral norm) of a matrix. for matrices A, U, C, V of conformable sizes, assuming all inverses exist and are well-defined. Lemma S2. Suppose that W  N (0, I N ), and let C be a symmetric N x N matrix with C 2  L. Then, for all 0 < t < L, The lemma is known in the literature. See, for example, Paul (2007). Lemma S3 (Lemma 2.7 of Bai & Silverstein (1998)). Let W = (w 1 , . . . , w p ) T , where w i 's are i.i.d. real r.v. 's with mean 0 and variance 1. Let B be a deterministic matrix. Then, for We collect several technical lemmas in this section. In what follows,    2 denotes the operator norm (spectral norm) of a matrix. for matrices A, U, C, V of conformable sizes, assuming all inverses exist and are well-defined. Lemma S2. Suppose that W  N (0, I N ), and let C be a symmetric N x N matrix with C 2  L. Then, for all 0 < t < L, The lemma is known in the literature. See, for example, Paul (2007). Lemma S3 (Lemma 2.7 of Bai & Silverstein (1998)). Let W = (w 1 , . . . , w p ) T , where w i 's are i.i.d. real r.v. 's with mean 0 and variance 1. Let B be a deterministic matrix. Then, for