Polynomial whitening for high-dimensional data

The inverse square root of a covariance matrix is often desirable for performing data whitening in the process of applying many common multivariate data analysis methods. Direct calculation of the inverse square root is not available when the covariance matrix is either singular or nearly singular, as often occurs in high dimensions. We develop new methods, which we broadly call polynomial whitening, to construct a low-degree polynomial in the empirical covariance matrix which has similar properties to the true inverse square root of the covariance matrix (should it exist). Our method does not suffer in singular or near-singular settings, and is computationally tractable in high dimensions. We demonstrate that our construction of low-degree polynomials provides a good substitute for high-dimensional inverse square root covariance matrices, in both d<N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d < N$$\end{document} and d≥N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d \ge N$$\end{document} cases. We offer examples on data whitening, outlier detection and principal component analysis to demonstrate the performance of the proposed method.


Introduction
Let X ∈ ℝ d×N be a matrix of data, with N observations in d dimensions. We denote the empirical d-dimensional mean vector and the empirical d × d covariance matrix of X by and respectively, and make no other assumptions about the generation or structure of the data. In this paper we consider transformations of X of the form X A = A(X − ), where A is a d × d matrix, with the aim of whitening the data X . Data whitening is a transformation of the data intended 1 3 to decorrelate and standardize the variables. Fully decorrelated data possesses a diagonal covariance matrix, and standardized data has unit variance for each variable (Hossain 2016). Applying a whitening transformation both decorrelates and standardizes the data, so that in the case of non-degenerate data, the covariance matrix of the whitened data X A will be the identity matrix. By removing the simple elliptic structure of the data through such whitening transformations, we can uncover more interesting and complex structures in the data that may have previously been hidden in correlations, such as clusters or outliers (Li and Zhang 1998). Furthermore, the orthogonality of whitened variables can improve computational time and performance of many statistical methods (Koivunen and Kostinski 1999;Huang et al. 2020;Zuber and Strimmer 2009).
Examples of existing whitening methods include the so-called Mahalanobis whitening defined by which is popularly used to whiten data before performing many classical methods of multivariate analyses (Kessy et al. 2018). The transformed data X −1∕2 has zerovalued mean and the d × d identity matrix I d as the covariance matrix. The success of the Mahalanobis whitening depends on the ability to compute −1∕2 in a way that is both accurate and stable.
It is common for big, high-dimensional data to be close to degeneracy/lowrank (Udell and Townsend 2019) yielding unstable computations of −1∕2 , with numerous examples of this problem observed in: recommender systems data (Zhou et al. 2008;Li et al. 2016); finance (Bai and Shi 2011); medicine (Schuler et al. 2016); genomics (Wu et al. 2015); and social networks (Liben-Nowell and Kleinberg 2007). These issues also arise in generalized mixture models (Xiao 2020); multiple regression (Healy 1968;Hoang and Baraille 2012); adaptive algorithms (Baktash et al. 2017); and linear discriminant analysis (Ye and Xiong 2006). This is because variables often possess (approximate) linear dependencies, resulting in a covariance matrix that is singular, or very close to singularity. As such, the inverse of the covariance matrix therefore does not exist or is at least unstable, and it becomes inadvisable or impossible to calculate −1∕2 . Consequently Mahalanobis whitening, and many other methods which directly rely on the inverse of the covariance matrix (such as those described in the survey of the recent paper (Kessy et al. 2018)), are not recommended.
Nevertheless, it has been demonstrated that applying Mahalanobis whitening prior to clustering or outlier detection (to give just two out of many possible examples) often results in better empirical results. This has been observed in several practical examples (Zafeiriou and Laskaris 2008;Shi et al. 2015). Theoretically, Mahalanobis whitening underpins weighted least squares (Seber and Lee 2012), PCA (Jolliffe 1986;Hyvärinen and Oja 2000), canonical correlation analysis (Härdle and Simar 2007) and most of the array of classic multivariate statistics methods (Li and Zhang 1998;Malsiner-Walli et al. 2016). Crucially, decorrelated and standardized data greatly simplifies both theoretical and practical multivariate data analysis (Agostinelli and Greco 2019; Anaya-Izquierdo et al. Polynomial whitening for high-dimensional data 2011; Martens et al. 2003;Thameri et al. 2011;Chen et al. 2015;Yang and Jin 2006).
The need to find stable ways of calculating (or approximating) −1 in settings when X is (close to) degeneracy is evidenced in other application domains. In high-dimensional examples one often would like to use the Mahalanobis distance (Mahalanobis 1936) to measure the 'scatter' or 'spread' of the data (Pronzato et al. 2017(Pronzato et al. , 2018, as a basis for proximity-dependent techniques such as clustering (Zuanetti et al. 2019) and Approximate Bayesian Computation (ABC) (Akeret et al. 2015). But again, for the reasons outlined earlier, the covariance matrix may be singular or ill-conditioned. ABC is used to find estimates of distribution parameters by simulating data over a parameter space (informed by some prior) and finding simulated data closest to the observed data. The ideal measure of closeness is to use the Mahalanobis distance, but the practical need to use ABC is often informed by degeneracy of the observed data, rendering the construction of a computationally tractable distance measure one of the fundamental problems of applying ABC (Wegmann et al. 2009;Beaumont 2019;Prangle 2017).
Recent literature has shown that whitening can also be used to improve the training of neural networks (Huang et al. 2018). Often, normalization is used in such training, rather than whitening, due to ill-conditioned problems (Luo 2017) and the great expense of computing a large inverse square root covariance matrix (Ioffe and Szegedy 2015), despite whitening being preferable if it is possible (Huang et al. 2020).
Naturally, many methods attempt to circumvent the aforementioned problems by use of the Moore-Penrose pseudo-inverse − to , and then to take the square root of − if needed. However, in the case of high-dimensional data, it has been shown that the Moore-Penrose pseudo-inverse is not always suitable (Hoyle 2011;Bodnar et al. 2016;Bickel and Levina 2004), particularly when there are small eigenvalues. This is a problem appearing in several branches of mathematics, and work on this topic can be found in the statistics literature (which tend to use shrinkage-type estimators (Ledoit and Wolf 2004;Fisher and Sun 2011;Ito et al. 2015)and/or assume sparsity of the covariance matrix (Cai et al. 2011(Cai et al. , 2016Janková and van de Geer 2017)) and in the linear algebra/numerical analysis literature. In this latter work, algorithms (alternating projections (Higham and Strabić 2016) or Newton-type (Qi and Sun 2011;Higham 2008)) to compute 'substitute' covariance matrices are developed. Challenges in the successful use of these algorithms include computational tractability, stability, and the ability to find good starting-points. Our work differs from this in that we provide explicit formulae for the construction of our substitute to −1∕2 and we are able to quickly generate a family of whitening matrices based on the order of our polynomials.
In this paper, we introduce polynomial whitening, and what we call the minimalvariance polynomial matrix, to be used in place of the square root of the inverse of the empirical covariance matrix. In view of the celebrated Cayley-Hamilton theorem (Cayley 1858; Hamilton 1853) the true inverse of a full-rank d × d matrix can be calculated through a d − 1 degree polynomial in . In Gillard et al. (2022), it is shown that an alternative to the inverse of a matrix can be found using low degree polynomials. Our work follows on from this, as we consider polynomials 1 3 of low degree in to now provide an alternative to the square root of its inverse. Our method is applicable in cases where the true inverse square root of the covariance matrix does not exist, making it a viable alternative for degenerate and close-to degenerate datasets. Parameter options also allow for a trade-off between data whitening accuracy and time complexity.
The main practical focus of this paper is in data whitening, but in view of the discussion above, we envisage other settings where our work may be useful. The structure of the paper is as follows. Section 2 introduces the form of our matrix polynomial, and the optimization problem we solve in order to obtain an alternative to the inverse square root of the covariance matrix. The main theorem of this paper which is studied in later examples is given in Sect. 2.3. We address different parameter choices in Sect. 2.4, and in Sect. 2.6 we discuss using our procedure in conjunction with random projection methods, which can be useful when dealing with very high-dimensional data. Examples applying our method to data whitening, outlier detection and dimension reduction are given in Sect. 3, before we conclude the paper in Sect. 4.

Covariance matrix of transformed data
The mean vector E(X A ) and covariance matrix D(X A ) of X A = A(X − ) are respectively: where 0 d is the d-dimensional vector of zeroes (Mathai and Provost 1992). Data transformed by Mahalanobis whitening, X −1∕2 = −1∕2 (X − ) , has covariance matrix where I d is the d × d identity matrix. The total variation in X A is given by trace D(X A ) = trace(A A ⊤ ), and for the Mahalanobis whitening the total variation in X −1∕2 is given by trace D X −1∕2 = trace(I d ) = d.

The minimal-variance polynomial alternative to ˙− 1∕2
We define A k to be a (k − 1)-degree matrix polynomial in , of the form: For a chosen integer k such that k − 1 < d , our objective is to find the k coefficients of the matrix polynomial, denoted = 0 , 1 , … , k−1 ⊤ in (1), so that the total

3
Polynomial whitening for high-dimensional data variation of the transformed data X A k = A k (X − ) , is minimized, subject to suitable constraints. For further intuition as to why we minimize the total variation, see Gillard et al. (2022). Let s j = trace( j ) , S (i,k) = s i , s i+1 , … , s i+k−1 , and define the matrix We seek to minimize the total variation of the transformed data X A k = A k (X − ) , which is given by: To ensure non-trivial solutions to the minimization of the total variation, we introduce a constraint. There are a number of options for this constraint; here we consider constraints of the form for a scalar value . This can be written in the above notation as A constraint of this form ensures that the minimal-variance polynomial matrix A k has similar qualities to −1∕2 , in the cases where this matrix exists. The constraint (2) will be revisited after the following theorem.

Constructing the minimal-variance polynomial
Theorem 1 Let X ∈ ℝ d×n be a d-dimensional dataset with n observations, having empirical mean and empirical covariance matrix .
is minimized, subject to the constraint ⊤ S ( ,k) = s −1∕2 , has coefficients given by ( Proof We will minimize 1 2 trace(D(X A k )) subject to the constraint (2), where the constant 1∕2 is introduced to simplify calculations. The Lagrange function L( , ) with Lagrange multiplier is given by We minimize the Lagrange function by differentiating with respect to and setting the result equal to 0, which gives: and we can therefore rearrange to find : Let = k . We can find the value of k by substituting (4) into the constraint Thus, the vector of coefficients of the polynomial (1) which minimizes trace D(X A k ) subject to the constraint (2) is given by (3). We call this polynomial the minimalvariance polynomial. ◻ When evaluating the polynomial, we recommend forming the matrix powers by iteratively multiplying by , or using Horner's method for polynomial evaluation. Both of these methods are outlined in Sect. 4.2 of Higham (2008).

Choice of the parameter ˛ in the constraint (2)
We studied the outcomes of polynomial whitening with different values of in the constraint (2). Theoretically, any value of will produce an alternative whitening matrix. Empirical investigations showed that the polynomial with = 1∕2 performed particularly well, in terms of data whitening success, stability and computational cost. Using this value of is equivalent to applying the constraint Our analysis found that, when using this constraint, data was approximately whitened using a relatively low value of k (when compared to the value of the dimension d of the dataset). Figure 1a considers a 50-dimensional dataset, with 5 eigenvalues greater than 1, 30 eigenvalues between 0 and 1, and 15 zero eigenvalues. Figure 1b considers a .
Polynomial whitening for high-dimensional data 150-dimensional dataset, with 5 eigenvalues greater than 1, 100 eigenvalues between 0 and 1, and 45 zero eigenvalues. The eigenvalues of these datasets are given in Appendix 1.1. These eigenvalues have been chosen to create a degenerate example which the Moore-Penrose pseudo-inverse would struggle to deal with well. We did so by setting roughly d∕3 eigenvalues equal to zero, and letting the nonzero eigenvalues taper towards zero, making the rank of the dataset unclear. We plot in blue dots the nonzero eigenvalues of the dataset on the horizontal axis, and the reciprocal square root of the nonzero eigenvalues on the vertical axis. We then plot the corresponding minimal-variance polynomial with = 1∕2 in degree k as follows. Find the coefficients = ( 0 , 1 , … , k−1 ) ⊤ of the minimal-variance polynomial using (3), and write the polynomial as in (1), replacing the matrix with a symbol t: We can then plot this polynomial p(t) for different values of t . In Fig. 1, we consider polynomials of degree 3, 4 and 5, (recalling that using the value k results in a (k − 1)-degree polynomial). The degree of these polynomials are much lower than the dimensionality of the datasets, yet still provide a good approximate fit to the inverse square root of the given eigenvalues. The constraint (5) with = 1∕2 works well in the case of non-degenerate data (when is essentially non-singular), but requires some simple tuning for degenerate or nearly-degenerate data, which has been applied here. This tuning will be discussed in Sect. 2.5.

Choice of the parameter k to determine the degree of the minimal-variance polynomial
The true inverse square root of a full-rank covariance matrix can be written as a (d − 1)-degree polynomial using the characteristic polynomial. The minimal-variance polynomial with parameter k forms a (k − 1)-degree polynomial, as defined in (1). As k increases, the polynomial can approximate the square root of the characteristic polynomial more accurately, but is more costly to compute. Therefore, the choice of the parameter k is often a trade-off between accuracy and cost. However, as k increases, so does the opportunity for instability in the polynomial, particularly when working in high dimensions (see Table 3 for an example of this, and Sect. 3.1 for more of a discussion on this topic). As such, keeping k relatively low is not only beneficial for cost, but for stability. Furthermore, choosing low values of k results in a good approximation for the inverse square root of the covariance matrix, as can be seen by the polynomial fit to the eigenvalues in Fig. 1. This will be further demonstrated in the numerical examples in Sect. 3.
In this paper, we often run the same experiments multiple times with different values of k , and use a problem-specific metric to identify the best value of k for that dataset. For example, in Sect. 3.1, we use the Wasserstein metric to compare the whitened data to the standard normal distribution, as well as a sum-of-squares-based metric. We then choose the value of k which produces the lowest values for these metrics. This is similar to methods used in many parameterized methods, such as using scree plots or silhouette scores to judge the best number of clusters to use in a clustering algorithm. It may often be best to apply minimal-variance whitening for multiple values of k to the dataset, and then inspect the empirical covariance matrix of the transformed data to see which value of k has performed best.

3
Polynomial whitening for high-dimensional data

Constraint adjustment for rank-deficient data
We provide an adaptation to the method when using = 1∕2 to make it suitable for use when is singular, as hinted at in discussions about the choice of . When using the constraint (5) with = 1∕2 , the polynomial aims to ensure that the trace of A k 1∕2 is equal to d . We propose that this trace should aim to equal r , the rank of the covariance matrix, in a similar way to the Moore-Penrose pseudo-inverse − having the property that trace ( − ) 1∕2 1∕2 = r . However, for matrices with many small eigenvalues, r is hard to calculate (Vidal and Favaro 2014), and often approximations of r are based on arbitrary eigenvalue thresholding or subjective elbow plots (Kishore Kumar and Schneider 2017).

Fig. 2
The minimal-variance polynomial with k = 5 fit to eigenvalues (blue, the same as in Fig. 1) before (red, dashed line) and after (green, solid line) adjustment for rank-deficiency, as described in Sect. 2.5 In cases where is not full-rank, we propose an adjustment to modify our constraint (2) without the need to calculate r directly. We will illustrate our adjustment using the two examples in Fig. 2, which plot the nonzero eigenvalues and nonzero reciprocal square root eigenvalues in the same way as in Fig. 1. The datasets used in Fig. 2 are the same as those in Fig. 1, details of which are given in Sect. 2.4 and Appendix 1.1.
We first calculate the minimal-variance polynomial using the constraint (5). This is shown in Figs. 2a and b as the red, dashed line. Although the polynomials take the correct shape, they are clearly placed too high and do not fit the plot of the inverse square root eigenvalues. We adjust the polynomial by multiplying by some constant c between 0 and 1, which ensures a better fit of the polynomial. Our method of choosing a value of c is as follows.
Let = { 1 , … , d } be the set of all eigenvalues of , and let ̃= { i ∈ ∶ i ≠ 0} be the set of all nonzero eigenvalues of . In the case of very large dimensions d , computation of eigenvalues i is certainly out of reach; in this case, as will be discussed in Sect. 2.6, we suggest to project the data to a low-dimensional space and use the set of eigenvalues for the low-dimensional version of the data. The constant c can be found in any number of ways which minimizes the distance between the polynomial p( ) and the target values 1∕ 1∕2 , for ∈̃ . We use the value c * from where w( ) is a suitable weight function. That is, we seek to minimize the weighted sum of squares between the polynomial and the reciprocal square root of the nonzero eigenvalues. The optimal value of the adjustment constant c is then found to be In Fig. 2 (and all other examples in this paper), we have used w( ) = , and in general we recommend this. However, the choice of w can be altered to give a different fit for the data given. If the user is more concerned about fitting the polynomial to the larger eigenvalues, they may decide to use w( ) = i with i > 1 , for example.
The adjusted polynomials (given by the green solid line) clearly fit the desired points much more successfully than the original polynomials. However, if this adjustment is not performed, the data transformed by the polynomial whitening matrix A k will still be approximately isotropic, so the adjustment is not necessary if equal variance is of variables is sufficient. This adjustment has been applied to all examples that follow in this paper.
This adjustment to the constraint can also be used to detect the singularity of a matrix. Let us first consider the case with d < N . If is full rank, and k is chosen appropriately, the value c * will be equal to (or very close to) 1, as the minimal-variance polynomial is aiming to make trace A k = d , which is correct in the case of full-rank . If the matrix is not full-rank, c * will be less than 1. To illustrate this, Table 1 Polynomial whitening for high-dimensional data dataset with N observations is generated using a covariance matrix generated with rank R . Further details on how these datasets were generated is given in Appendix 1.1. The empirical covariance matrix of the dataset has rank r = min(d, N, R) , and is used to find the minimal-variance polynomial matrix with k = 10 , and the constraint adjustment c * is given. In dataset 1, the empirical matrix has full rank r = d , so c * = 1 . In dataset 2, the 'true' covariance matrix has rank R = 50 , d = 100 and N = 1000 , therefore the empirical covariance matrix has rank r = min(d, N, R) = 50 . This produces a constraint adjustment value of c * = 0.50 < 1 , so we know the empirical covariance matrix is singular. We now consider cases with d ≥ N through the use of three examples. Dataset 3 in Table 1 has 100 dimensions and only 50 observations. The 'true' covariance matrix used to generate this dataset is full-rank R = 100 , but the empirical covariance matrix has rank r = min(d, N, R) = 50 . Therefore, using the empirical covariance matrix in the minimal-variance polynomial matrix gives adjustment value c * = 0.50 , informing us that this dataset is degenerate. Dataset 4 also has d = 100 , N = 50 , and the 'true' covariance matrix now has rank R = 50 . The adjustment value is therefore less than 1: c * = 0.50 . The final example we consider again has dimension d = 100 and number of observations N = 50 , but the 'true' covariance matrix has rank R = 30 . The empirical covariance matrix therefore has rank r = min(d, N, R) = 30 , and the adjustment value is c * = 0.30 . In all these examples, c * < 1 , as the empirical covariance matrix will never be full-rank in d < N examples.

Applications to extremely high-dimensional data
Given a dataset X with extremely high-dimension d , say d = 1, 000, 000 , finding the minimal-variance polynomial matrix can be too costly and time-intensive. We can instead sample some variables from X to produce a 'representative' dataset X in a much smaller dimension d . This representative dataset can be found through random samples of the variables in X , or projection to a lower dimensional space (see Bingham and Mannila (2001), Blum et al. (2014)). We can proceed with calculating the covariance matrix ̃ of X , and use ̃ to produce the minimal-variance polynomial alternative to ̃− 1∕2 : We can then replace the d -dimensional matrix ̃ in (6) with the d-dimensional covariance matrix to obtain the minimal-variance polynomial matrix A k . This can be used to whiten the original large dataset X , and is much cheaper than finding the minimal-variance polynomial matrix directly. For large datasets, it may be that we don't know the eigenvalues exactly, but can approximate the distribution of the eigenvalues. If this is the case, we can sample d eigenvalues from this distribution using the inverse cumulative distribution function. We will illustrate this using the Marchenko-Pastur distribution, as this distribution is known to model the eigenvalues of the sample covariance matrix of a random matrix as d, N → ∞ . Figure 3 considers an example with d = 10, 000 and N = 15, 000 , and the probability density function (PDF) of the Marchenko-Pastur distribution with these parameters is shown by the red line. The histogram represents a random sample of 300 eigenvalues, and shows such a sample models the distribution well.
Similar methods can be used in the case where d ≥ N . We can reduce d to a value smaller than N by sampling the real eigenvalues, if they are known or can be calculated. Alternatively, we can sample the eigenvalues from an assumed distribution as described above.

Data with d < N
We begin the numerical examples by whitening several synthetic and real datasets using the minimal-variance polynomial. The details of these datasets are given in Table 2. The four synthetic datasets (D1, D2, D3, D4) are sampled from a  (Vanschoren et al. 2013).
In some cases, it can be beneficial to rescale the data so that each variable has zero mean and unit variance, before finding the minimal-variance polynomial matrix. If rescaling the data provides less extreme eigenvalues in the covariance matrix, this scaling is likely to improve the performance of the polynomial whitening. The heatmaps in Fig. 4 show the covariance matrices of the datasets, and the distribution of the eigenvalues of these covariance matrices are given in  Figure 4 shows a lot of nonzero off-diagonal values in the heatmaps, indicating that these datasets are highly correlated. We can measure the proximity of the transformed data X A k ∼ N d (0, S) to the standard normal distribution N d (0, I) using the Wasserstein metric (Givens and Shortt 1984): where we divide by d here to account for the difference in the dimensions of each dataset.
The heatmaps in Fig. 5 show the covariance matrix of the transformed data X A k = A k X for each dataset, illustrating that the correlations between variables have been approximately whitened. The value of k used in these heatmaps is chosen as the value of k which gives the lowest Wasserstein score, as given in Table 3. The  Table 2 after minimal-variance polynomial whitening. The value of k used in constructing the minimal-variance polynomial is given in the caption for each dataset Table 3 The Wasserstein scores (7), denoted W A k X , which measure the distance between the polynomialwhitened dataset A k X and the standard normal distribution N(0, I) for each dataset Values in bold indicate the lowest Wasserstein score W A k X over all k for a given dataset

3
Polynomial whitening for high-dimensional data Wasserstein scores in Table 3 show that, in general, as the value of k increases, the transformed data is closer to the standard normal distribution, as desired. In some cases, such as the Musk dataset, higher values of k begin to show an increase in the Wasserstein score, indicating the whitening transformation is less successful than when using lower values of k . This is likely due to numerical instability, as the minimal-variance polynomial aims to fit itself to extremely small eigenvalues, causing erratic behaviour in the polynomial. As such, it is recommended to use lower values of k which provide a more reliable alternative to the inverse square root of the covariance matrix, or to compute several minimal-variance polynomial matrices for different k and use the one that best satisfies some metric, such as the Wasserstein score.
As indicated by the Wasserstein scores in Table 3 and the heatmaps in Fig. 5, we produce an effective alternative to the inverse square root of the covariance matrix using a polynomial of degree significantly lower than the dimension of the dataset.
The Wasserstein measure concerns itself with the diagonal values of the covariance matrix, as it is calculated using traces. We can consider it as a measure of standardization, rather than whitening. We therefore need to measure the extent to which the data has been decorrelated. The heatmaps in Fig. 5 show that the off diagonals of the covariance matrix of the transformed data are close to zero, indicating good decorrelation. Another way we can measure this is by considering the sum of squares of the off-diagonal entries of the covariance matrix of the transformed data. In Table 4, let SS A k X be the sum of squares of the off-diagonal entries of the covariance matrix of the whitened dataset A k X.
The sum of squares values in Table 4 decrease as k increases, until a certain value of k , much like the Wasserstein scores. Given we would like this value to be as small as possible, we see the value of k that gives the optimum sum of squares value for each dataset is close to value of k that gives the optimum Wasserstein score for each dataset. Therefore, when the data has been successfully standardized, it has also been decorrelated well. Table 13 in Appendix 3 shows the average time taken to produce the minimalvariance polynomial matrices for each dataset for each value of k considered, over Table 4 The sum of squares, denoted SS A k X , of the off diagonal values of the covariance matrix of the polynomial-whitened dataset A k X for each dataset Values in bold indicate the lowest value of SS A k X over all k for a given dataset 100 runs. The time taken increases as the dimensionality d of the dataset increases, and as the parameter k increases. However, the procedure to calculate the matrices only takes a matter of seconds, even for 1000-dimensional datasets. This time performance could be improved much further by implementing parallel computing methods.

Data with d ≥ N
It is increasingly common for data to have higher dimensionality than number of observations in many fields, such as genetic microarrays, medical imaging and chemometrics (Hall et al. 2005). Such data is clearly rank-deficient, with r ≤ N < d , and thus the sample covariance matrix of such data is always singular, rendering many multivariate data analysis methods unusable, including data whitening. Minimal-variance polynomial whitening is applicable in such cases, as illustrated by the following examples. We consider four synthetically generated datasets and four real datasets, detailed in Table 5. The first two synthetic datasets, E1 and E2, are sampled from a Gaussian distribution N d (0, ) , where the covariance matrices are produced as follows. Like the d < N case, we generate d eigenvalues, and produce a random d × d orthogonal matrix Q . Let L be the matrix with the eigenvalues on the diagonal and zeroes elsewhere, then let = Q ⊤ LQ . The third synthetic dataset, E3, is generated to copy the example in Wang and Fan (2017): a multivariate Gaussian with population covariance matrix with diagonal entries [50, 20, 10] + [1] * 47 . This creates a spiked eigenvalue model, which is of interest in HDLSS datasets . The fourth dataset uses a covariance matrix with eigenvalues generated from a random uniform distribution between 0 and 1, to produce a non-sparse set of eigenvalues. The madelon dataset was obtained from the UCI Machine Learning Repository (Dua and Graff 2017). The raw madelon dataset has 4400 observations, greater than the 500 features, so we sampled only the first 250 observations to create the madelon † dataset with d > N . The yeast dataset is a real genomic dataset with 2284 features and 17 observations (Tavazoie et al. 1999;Vanschoren et al. 2013). The third real dataset is a genomic dataset on colon cancer data (Alon et al. 1999), used by (Yata and Aoshima 2013) as an example of a spiked Successful whitening of these datasets would result in a covariance matrix with r eigenvalues equal to 1, and d − r eigenvalues equal to 0. We performed Moore-Penrose whitening on the four datasets in Table 5 by pre-multiplying the data by the Moore-Penrose inverse of the square root of the covariance matrix. We also performed minimal-variance polynomial whitening on the datasets as described in Sect. 2. Figure 6 compares the distribution of the eigenvalues of the covariance matrices after Moore-Penrose whitening and minimal-variance polynomial whitening. The eigenvalues are scaled such that the maximum eigenvalue is equal to 1. The first three synthetic datasets show that using minimal-variance whitening returns a dataset with eigenvalues only equal to 0 and 1, whereas using Moore-Penrose whitening gives a dataset with a spread of eigenvalues between 0 and 1. Figure 6d shows that minimal-variance whitening may not achieve perfect whitening, but that it is still more successful than the Moore-Penrose whitening method. Figure 6f considers the Yeast dataset, and shows that both Moore-Penrose whitening and minimal-variance whitening return only eigenvalues of value 0 or 1. However, the Moore-Penrose whitening gives one eigenvalue equal to 1, and the rest 0 (or very close to 0). When using the minimal-variance whitening, the dataset has rank equal to the original dataset ( r = 16 in this case). The madelon † , colon and DB-emails datasets are not whitened perfectly by either method, but the eigenvalues are much more dispersed when using Moore-Penrose whitening compared to minimal-variance polynomial whitening, whereas we seek eigenvalues only valued at 0 and 1, ideally.

Comparison to other whitening methods
Due to rotational freedom, there are infinitely many whitening matrices of the form W = Q −1∕2 , where Q is orthogonal and satisfies Q ⊤ Q = I d (Kessy et al. 2018).
Let us define some decompositions of the covariance matrix , beginning with = V 1∕2 PV 1∕2 , where V is the diagonal variance matrix and P is the correlation matrix. Let = U U ⊤ be the eigendecomposition of the covariance matrix, with U the matrix of eigenvectors and the diagonal matrix of eigenvalues. Analogously, define the eigendecomposition P = GOG ⊤ of the correlation matrix. We also define the Cholesky decomposition of the inverse covariance matrix LL ⊤ = −1 , when −1 exists.
Five whitening procedures are identified by Kessy et al. (2018) to be unique in fulfilling a given objective function. Most of these objective functions used in the paper are based on the cross-covariance matrix and the cross-correlation matrix between the original data X with covariance and the whitened data X W :

Fig. 6
Log-scale histograms, showing the eigenvalues of the covariance matrix after the data has been whitened by Moore-Penrose (MP) whitening (orange histogram), and minimal-variance polynomial (MV) whitening with k = 9 (green histogram), for each of the four datasets given in Table 5 Polynomial whitening for high-dimensional data In the following example, we will compare polynomial whitening to the three procedures from this paper that share our goal of whitening the data while changing as little else as possible. We do not consider the other methods in this paper, as these methods aim to maximize compression of variance into the first few variables of the whitened data. Although polynomial whitening performed relatively well in these scenarios, this is not the aim of our method. The three types of whitening we will consider alongside polynomial whitening are given below.
Mahalanobis whitening (MW) W = −1∕2 . Mahalanobis whitening is found to be the unique whitening procedure which maximizes trace( ) , the average crosscovariance between each variable of the original and the newly transformed data. This is equivalent to minimizing the total squared distance between the original data X and the whitened data X W , ensuring the whitened data is as similar as possible to the original data.
Mahalanobis-cor whitening (MCW) W = P −1∕2 V −1∕2 . Mahalanobis whitening can be affected by the differences in the scales of variables. To avoid this issue we may use a scale-invariant version, known as Mahalanobis-correlation whitening. The Mahalanobis-correlation whitening method maximizes the cross-correlation trace( ) between each variable of the standardized original data V −1∕2 X and the whitened data X W . Doing this is shown to be equivalent to minimizing the squared distance between V −1∕2 X and X W .
Cholesky whitening (CW) W = L ⊤ . Cholesky whitening is the only whitening procedure fulfilling the constraint of producing lower-triangular cross-covariance and cross-correlation matrices with positive diagonal entries. It does not result from fulfilling an objective function like the above methods, but rather from satisfying this constraint.
We evaluate the performance of these different whitening procedures by applying them to a dataset and considering the different objective functions in and . First, as in Kessy et al. (2018), we apply the whitening methods to the 4-dimensional Iris dataset (Fisher et al. 1936) in Table 6. Given the dataset's low dimension and well-conditioned covariance matrix, polynomial whitening (PW in the table) with k = d = 4 produces exactly the same results as Mahalanobis whitening. We also perform polynomial-cor whitening (PCW), where the data is standardized and polynomial whitening is performed using the correlation matrix P . This produces the same results as Mahalanobis-cor whitening. The polynomial whitening method is more effectively used when applied to higher dimensional datasets with singular or near-singular covariance matrices. As such, we repeat the above exercise with a different dataset. For the purposes of this example, we are unable to use a dataset which has a singular covariance matrix, as the Mahalanobis and Cholesky whitening methods are not usable in this case. We use the Wisconsin Breast Cancer dataset (Wolberg et al. 1992), which we have pre-standardized to give improved results from all methods. This dataset has dimension d = 32 and has a covariance matrix which could be considered ill-conditioned (see Appendix 1.4 for details on the eigenvalues). Table 7 shows that polynomial whitening outperforms Mahalanobis whitening, using both the covariance and correlation matrix.

The effect of different pre-processing methods on outlier detection algorithms
Outlier detection algorithms often require that data is pre-processed before the algorithm can be applied. It has been shown by Campos et al. (2016) that the normalization of datasets will often lead to a better performance of outlier detection algorithms.
Here we replicate a study described in Kandanaarachchi et al. (2020). The authors produced a collection of labelled benchmark datasets to be used for evaluating outlier detection algorithm performance. They evaluated the performance of various algorithms when used after applying different normalization methods to these datasets. Performance of an algorithm was measured using the area under the Receiver Operator Characteristic (ROC) curve, which compares the labels of an observation ('inlier' or 'outlier') produced by the algorithm to the 'true' labels. They found that two types of normalization method performed differently (dependent on data set and outlier detection method): 'Min-Max' normalization Each variable v of a dataset is normalized to only have values in the range [0, 1]: where min(v) and max(v) are the minimum and maximum values of the variable v , respectively.
' where median(v) and IQR(v) are the median and inter-quartile range of the variable v , respectively. We consider the following four different outlier detection methods from the Python package PyOD (Zhao et  Further details of each of these algorithms are provided in Campos et al. (2016). All of the above methods require a parameter choice K (different to the polynomial degree parameter k referred to throughout this paper) to set the so-called neighbourhood size, and a contamination value C to indicate how many observations the algorithm should label as outliers. We let K = 0.1 × N , where N is the number of observations in the dataset. The parameter C is equal to the percentage of outliers given by the 'true' labels.
For a dataset D , an outlier detection method o and a pre-processing method z , we denote the area under the ROC curve as AUC (D, o, z) . For each outlier detection method o listed above, we say a dataset D 'prefers' a pre-processing method z if AUC(D, o, z) ≥ AUC (D, o, y) for all other pre-processing methods y . We evaluate the AUC score for transformations A k D using (3) by taking the maximum AUC score over all k considered.
We tested the outlier detection methods with each pre-processing method on 7667 real datasets, as used in Kandanaarachchi et al. (2020). The datasets ranged from dimension 3 to dimension 359, and the number of observations in a dataset ranged from 44 to 5396. We impose no structural assumptions on the datasets for our method or the other normalization methods. Table 8 shows the percentage of datasets that prefer each pre-processing method for each of the given outlier detection algorithms. The results in this table indicate that the polynomial whitening method outperforms the two normalization methods.
The scatter graphs in Fig. 7 compare the minimal-variance polynomial whitening to the normalization methods considered individually. Each point represents a v − median(v) IQR(v) , dataset, and the diagonal line indicates those datasets where the two methods give equal AUC scores. Points below this line, in red, indicate that the minimal-variance whitening method outperformed the other method considered. A numerical breakdown of these scatter graphs is given in Table 9. Much like Table 8, Table 9 shows the percentage of datasets that prefer each pre-processing method, but shows a pairwise comparison. Table 10 shows the amount of datasets out of the total 7667 (and the percentage) for which the pre-processing methods produce strictly better results, for each outlier detection method. It is clear that the minimal-variance method performs as well as (and often better than) the techniques often used to preprocess datasets before applying common outlier detection methods.

Fig. 7
Scatter graphs plotting the AUC scores of outlier detection algorithms when performed using the minimal-variance polynomial whitening 'Min-Var' on the horizontal axis, and the AUC scores when using a-d 'Min-Max' or e-h 'Med-IQR' normalizations on the vertical axis. Points in red indicate a dataset where using Min-Var produced a better score than the alternative method, and points in blue indicate a dataset where using the alternative method produced a better score Table 9 The percentage of datasets for which the given pre-processing method (given in the column) produces AUC scores better than the alternative method in the adjacent column, for different outlier detection methods (given in the row). I.e. 34.4% of datasets produced higher AUC scores when using Min-Var than when using Min-Max, for the KNN outlier detection method. This differs from Table 8

3
Polynomial whitening for high-dimensional data

Principal component analysis
Principal component analysis (PCA) is a popular dimension-reduction technique, as it reduces a dataset to a chosen dimension p while retaining the greatest amount of variance from the original dataset as possible. PCA finds p linear combinations of the variables of the dataset, giving p new compressed variables with maximal variance. As such, it is highly sensitive to the variances of the variables in the dataset. If one variable is measured on a much larger scale than the others, this variable will be likely have much greater variance, and therefore be given much more weight in a linear combination than the other variables (Jolliffe and Cadima 2016). To prevent this, it is good practice to standardize the variables to ensure they are all measured on the same scale. We compare two methods of standardization prior to performing PCA: (Moore-Penrose) Mahalanobis standardization, which is most commonly used before PCA, and minimal-variance standardization. In Mahalanobis standardization, let the variable v i ∈ X have mean i and standard deviation i . We then consider the dataset made up of the variables for i ∈ {1, … , d} . If i = 0 , we use Moore-Penrose Mahalanobis (MPM) standardization, in which we find the square root of the Moore-Penrose inverse of the covariance matrix − , and then use ( − ) i,i (i.e. the i th diagonal entry of − ) in place of i .
In minimal-variance (MV) standardization, we find the minimal-variance polynomial matrix A k , and use the values on the diagonal of A k , denoted (A k ) i,i , in place of i : Note that this is different to minimal-variance whitening, in that we only use the diagonal of the minimal-variance polynomial matrix to perform the transformation. We do this to align our method with the Mahalanobis standardization method. Our method of comparing the different standardization methods for PCA is as follows. In the following sections, we consider 1000 generated datasets, each with K clusters. For each dataset, we consider three versions: let X be the original dataset, X MPM be the MPM standardized dataset and X MV be the MV standardized dataset. For each version, we find the data given by the PCA transformation for a given number of principal components, and then perform the K-means clustering algorithm (Lloyd 1982). This is repeated for 1000 different datasets, and the results are given in the following sections for different types of data.

Data with d < N
We first consider the impact of the different standardization methods on PCA for data with d < N . For each of the 1000 datasets, we generate 3 clusters from multivariate Gaussian distributions X (i) , i = {1, 2, 3} with dimension d = 100 , where the parameters (i) , (i) and N (i) denote the mean, covariance matrix and number of observations in cluster X (i) . The details of these parameters are given in Table 11. The eigenvalues of each (i) taper off towards zero gradually. This creates a degenerate dataset with a rank that is hard to identify, a situation which the Moore-Penrose inverse struggles to deal with well.
In this example, we make parameter choices based on the relative size of the eigenvalues of a dataset compared to the maximum eigenvalue. Let = { 1 , … , d } be the set of eigenvalues of a dataset, let max be the largest eigenvalue in , and let ̄ be the mean of the eigenvalues in .
Let p = p( ) be the number of principal components that we wish to reduce a dataset to using PCA. For each dataset, the parameter p( ) is chosen to be the number of eigenvalues in greater than the mean eigenvalue ̄: as commonly used in practice (Abdi and Williams 2010).
The parameter k = k( ) for the minimal-variance polynomial will chosen based on the number of scaled eigenvalues i = i ∕ max that are bigger than a given threshold t: In the examples that follow we use t = 0.1. The K-means clustering algorithm aims to assign each point within a dataset to a cluster, by estimating the distances from each point to the estimated centre-point of a cluster of points. For more information on the algorithm, see Jain (2010). We consider how well the K-means clustering algorithm performs after applying PCA, given the different standardization methods. We use the adjusted rand (AR) score (Hubert and Arabie 1985;Steinley 2004) of the cluster labels provided by K-means to judge how well the algorithm has found the correct clusterings. An AR score of 0 indicates random labellings, and an AR score of 1 means the clusters were perfectly labelled by the algorithm.
We also give the silhouette scores (Rousseeuw 1987) of the methods depending on the standardization methods. The silhouette score of a clustering indicates how well separated the clusters are. A score of 1 indicates well-distinguished clusters, whereas a score of -1 tells us that clusters have been incorrectly assigned. A higher silhouette score tells us that the standardization method and PCA have retained cluster structure well. Figure 8 shows that using the MPM standardization gives a slight improvement on using no standardization. However, using minimal-variance standardization before applying PCA and K-means clustering results in vastly better AR scores, as well as better silhouette scores.

Data with d ≥ N
We modify the above example slightly to help us consider the case where d ≥ N . In such circumstances, PCA can sometimes perform poorly due to difficulties in finding the eigenvectors of the covariance matrix correctly . As The latter two methods were formulated for performing dimension reduction on HDLSS data, and can avoid the difficulties sometimes faced by PCA in such settings. Examples given in these papers show promising results for eigenvalue estimation and dimension reduction in high dimensions. For more details on the implementation of these methods, see the papers referenced above. As in Sect. 3.4.1, we consider 1000 different datasets, each with d = 1000 and N = 430 . Each dataset is generated as a mixture of four multivariate Gaussian distributions X i ∼ N d i , i , i = {1, 2, 3, 4} . The population parameters of each cluster are given in Table 12. Figure 9 shows boxplots of the AR scores and silhouette scores of the labels given by K-means clustering, after applying one of the standardization methods and one of the three dimension reduction methods. Across all types of PCA, we see that MPM standardization gives very similar results to the datasets with no standardization. On the other hand, the MV standardization method provides a large improvement for all three methods of dimension reduction. The clustering results are better when MV standardization has been used, as indicated by the boxplots of AR scores in Fig. 9a, c and e. We see that cluster separation is also better when using MV standardization with dimension reduction, as the silhouette scores are much higher for all three dimension reduction methods.
Across the three dimension reduction techniques considered, the minimal-variance standardization method is clearly very useful in those cases where standardization would improve dimension reduction algorithms (or other multivariate data analysis methods), as it behaves similarly to the Moore-Penrose Mahalanobis standardization method, but does not struggle in cases where the rank is unclear and there are many small eigenvalues. 1 3 Polynomial whitening for high-dimensional data

Conclusion
We have developed a method of constructing polynomials in the empirical covariance matrix to provide an alternative to the inverse square root of a covariance matrix, particularly suitable to degenerate (or close to degenerate) data in high dimensions. The minimal-variance polynomial whitening method aims at minimizing the total variation in a transformed dataset, and in doing so it provides a dataset that has been decorrelated and standardized. We have demonstrated the potential applications of these polynomial matrices by considering whitening, outlier detection and principal component analysis for both d < N and d ≥ N cases. We have discussed and given recommendations for the choice of the parameter k which dictates the degree of the polynomial, as well as an alternative constraint and adjustments. We also suggested a method to reduce computational time in extremely high-dimensional cases, ensuring that this method can be applied in such scenarios.

Eigenvalues of the datasets in Section 3.1
Datasets in Section 3.1.1, with d < N The histograms in Fig. 10 give the distribution of the eigenvalues of the datasets used in Sect. 3.1.1, detailed in Table 2.

Datasets in Section 3.1.2, with d ≥ N
The histograms in Fig. 11 give the distribution of the eigenvalues of the datasets used in Sect. 3.1.2, detailed in Table 5.
Time to compute the minimal-variance polynomial in Section 3.1 Table 13 gives the time it took to calculate the minimal-variance polynomial (in seconds) for each dataset used in Sect. 3.1, for each different value of k.

Eigenvalues of the datasets in Section 3.2
In Sect. 3.2, we compare different whitening methods with the minimal-variance polynomial whitening method by applying them to the Iris dataset and the Wisconsin breast cancer dataset (the latter of which we have scaled to improve performance). The eigenvalues of these datasets are given below:  Fig. 11 Eigenvalues of the datasets used in Section 3.1.2. A dataset with a * in its caption has been rescaled such that each variable has zero mean and unit variance, and hence the eigenvalues of the correlation matrix are given