# Fast Whitening Transformation

#### Developed by Gabriel Hoffman

#### Run on 2025-07-18 08:21:55.634263

Data whitening is a widely used preprocessing step to remove correlation structure since statistical models often assume independence (Kessy, et al. 2018). The typical procedures transforms the observed data by an inverse square root of the sample correlation matrix (Figure 1). For low dimension data (i.e. \(n > p\)), this transformation produces transformed data with an identity sample covariance matrix. This procedure assumes either that the true covariance matrix is know, or is well estimated by the sample covariance matrix. Yet the use of the sample covariance matrix for this transformation can be problematic since 1) the complexity is \(\mathcal{O}(p^3)\) and 2) it is not applicable to the high dimensional (i.e. \(n \ll p\)) case since the sample covariance matrix is no longer full rank.

Here we use a probabilistic model of the observed data to apply a whitening transformation. Our Gaussian Inverse Wishart Empirical Bayes (GIW-EB) 1) model substantially reduces computational complexity, and 2) regularizes the eigen-values of the sample covariance matrix to improve out-of-sample performance.

Figure 1: Intuition for data whitening transformation. Given the singular value decomposition of \(Y/\sqrt{n}\) after centering columns is \(UDV^T\), the steps in the transformation are: A) Original data, B) Data rotated along principal components, C) Data rotated and scaled, D) Data rotated, scaled and rotated back to original axes. Green arrows indicate principal axes and lengths indicate eigen-values.

# Basic usage

    library(decorrelate)
    library(Rfast)

    n <- 500  # number of samples
    p <- 200  # number of features

    # create correlation matrix
    Sigma <- autocorr.mat(p, 0.9)

    # draw data from correlation matrix Sigma
    Y <- rmvnorm(n, rep(0, p), sigma = Sigma * 5.1, seed = 1)
    rownames(Y) <- paste0("sample_", 1:n)
    colnames(Y) <- paste0("gene_", 1:p)

    # eclairs decomposition implements GIW-EB method: *E*stimate
    # *c*ovariance/correlation with *l*ow *r*ank and *s*hrinkage
    ecl <- eclairs(Y)

    # decorrelate data using eclairs decomposition
    Y_whitened <- decorrelate(Y, ecl)

    # the same whitening can be performed with one command where the eigen-value
    # shrinkage is performed internally
    Y_whitened2 <- whiten(Y)

Here plot A) the eigen-values of the covariance matrix before and after shrinkage. B) The correlation of the observed and C) whitened data.

    oldpar <- par(mfrow = c(1, 3))

    # plot shrinkage of eigen-values
    plot(ecl)

    # correlation between variables in observed data
    image(cor(Y), axes = FALSE, main = "Correlation of observed data")

    # decorrelate data using eclairs decomposition
    image(cor(Y_whitened), axes = FALSE, main = "Correlation of whitened data")

    par(oldpar)

# Advanced usage

The `decorrelate` package has advanced features to examine details of the whitening transformation.

## Directly compute whitening matrix

While `eclairs()`, `decorrelate()`, and `whiten()` perform the probabilistic whitening transformation efficiently without directly computing the whitening matrix, `getWhiteningMatrix()` can directly compute the matrix.

    # compute whitening matrix from eclairs decomposition
    W <- getWhiteningMatrix(ecl)

    # transform observed data using whitening matrix
    Z <- tcrossprod(Y, W)

    # evalute difference between whitened computed 2 ways
    max(abs(Z - Y_whitened))

    ## [1] 4.374279e-14

The difference between function is due only to machine precision.

## Explicit covariance or correlation

The full covariance for correlation matrix implied by the `eclairs()` decomposition can be computed explicitly. Note that computing and storing these matries is \(O(p^2)\), it may not feasable for large datasets.

    # compute correlation matrix from eclairs
    getCor(ecl)

    # compute covariance matrix from eclairs
    getCov(ecl)

## Sample from multivariate normal

The form of the `eclairs()` decomposition can be used to efficiently sample from a multivariage normal distribution with specified covariance.

    # draw from multivariate normal
    n <- 1000
    mu <- rep(0, ncol(Y))

    # using eclairs decomposition
    X.draw1 <- rmvnorm_eclairs(n, mu, ecl)

## Low-rank models

A low rank `eclairs()` decomposition can be computed more efficiently when \(k\) is small relative to \(min(n,p)\). Importantly, the emprical Bayes estimate of the shrinkage parameter \(\lambda\) can still be computed accurately for sufficiently large \(k\). Note that the low rank method trades computational efficientcy for accuracy in the whitening transform.

    # use low rank decomposition with 50 components
    ecl <- eclairs(Y, k = 60)

    # decorrelate data using eclairs decomposition
    Y_whitened <- decorrelate(Y, ecl)

In this case, the low rank whitening produces transformed features that are approximately independent. The approximation improves as the rank increases.

## Computing condition number

Compute the condition number (i.e. the ratio between the largest and smallest eigen-value) of the correlation/covariance matrix from the `eclairs()` decomposition.

    kappa(ecl)

    ## [1] 88.64936

## Removing correlation vs covariance

By default `eclairs()` computes the covariance between columns by using the default `compute = "covariance"`. Running `decorrelate()` using this removes the covariance between columns. Setting `compute = "correlation"`, evaluates and then removes correlation between columns while retaining the variance.

    library(clusterGeneration)

    # generate covariance matrix, where the diagonals (i.e. variances) vary
    Sigma <- genPositiveDefMat(p, rangeVar = c(1, 1e+06))$Sigma

    Y <- rmvnorm(n, rep(0, p), sigma = Sigma)

    # examine variances of the first 5 variables
    apply(Y, 2, var)[1:5]

    ## [1]  6.279349  6.171918  2.320736  1.306302 10.851267

    # transform removes covariance between columns so variance of transformed
    # features are *approximately* equal
    ecl_cov <- eclairs(Y, compute = "covariance")
    Z1 <- decorrelate(Y, ecl_cov)

    # variance are *approximately* equal
    apply(Z1, 2, var)[1:5]

    ## [1] 0.9757741 0.9929288 0.9928508 0.9924159 0.9902491

    # transform removes **correlation** between columns but variables are not
    # scaled
    ecl_cor <- eclairs(Y, compute = "correlation")
    Z2 <- decorrelate(Y, ecl_cor)

    # variances are not standardized
    apply(Z2, 2, var)[1:5]

    ## [1]  6.173119  6.131217  2.296492  1.286559 10.779611

# Session Info

    ## R version 4.4.2 (2024-10-31)
    ## Platform: aarch64-apple-darwin23.6.0
    ## Running under: macOS Sonoma 14.7.1
    ##
    ## Matrix products: default
    ## BLAS:   /Users/gabrielhoffman/prog/R-4.4.2/lib/libRblas.dylib
    ## LAPACK: /opt/homebrew/Cellar/r/4.5.1/lib/R/lib/libRlapack.dylib;  LAPACK version 3.12.1
    ##
    ## locale:
    ## [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
    ##
    ## time zone: America/New_York
    ## tzcode source: internal
    ##
    ## attached base packages:
    ## [1] stats     graphics  grDevices utils     datasets  methods   base
    ##
    ## other attached packages:
    ##  [1] Rfast_2.1.5.1           RcppParallel_5.1.10     zigg_0.0.2
    ##  [4] RcppZiggurat_0.1.8      Rcpp_1.1.0              clusterGeneration_1.3.8
    ##  [7] MASS_7.3-65             latex2exp_0.9.6         colorRamps_2.3.4
    ## [10] cowplot_1.2.0           ggplot2_3.5.2           mvtnorm_1.3-3
    ## [13] decorrelate_0.1.6.4
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] sass_0.4.10        generics_0.1.4     stringi_1.8.7      lattice_0.22-7
    ##  [5] digest_0.6.37      magrittr_2.0.3     evaluate_1.0.4     grid_4.4.2
    ##  [9] RColorBrewer_1.1-3 fastmap_1.2.0      jsonlite_2.0.0     Matrix_1.7-3
    ## [13] formatR_1.14       scales_1.4.0       codetools_0.2-20   jquerylib_0.1.4
    ## [17] cli_3.6.5          rlang_1.1.6        CholWishart_1.1.4  withr_3.0.2
    ## [21] cachem_1.1.0       yaml_2.3.10        tools_4.4.2        parallel_4.4.2
    ## [25] dplyr_1.1.4        corpcor_1.6.10     vctrs_0.6.5        R6_2.6.1
    ## [29] lifecycle_1.0.4    stringr_1.5.1      irlba_2.3.5.1      pkgconfig_2.0.3
    ## [33] pillar_1.11.0      bslib_0.9.0        gtable_0.3.6       glue_1.8.0
    ## [37] whitening_1.4.0    xfun_0.52          tibble_3.3.0       tidyselect_1.2.1
    ## [41] knitr_1.50         farver_2.1.2       htmltools_0.5.8.1  rmarkdown_2.29
    ## [45] labeling_0.4.3     compiler_4.4.2