# How to calculate tridiagonal approximate covariance matrix, for fast decorrelation?

Given a data matrix $X$ of say 1000000 observations $\times$ 100 features, is there a fast way to build a tridiagonal approximation $A \approx cov(X)$ ?
Then one could factor $A = L L^T$, $L$ all 0 except $L_{i\ i-1}$ and $L_{i i}$, and do fast decorrelation (whitening) by solving $L x = x_{white}$. (By "fast" I mean $O( size\ X )$.)

(Added, trying to clarify): I'm looking for a quick and dirty whitener which is faster than full $cov(X)$ but better than diagonal. Say that $X$ is $N$ data points $\times Nf$ features, e.g. 1000000$\times$ 100, with features 0-mean.

1) build $Fullcov = X^T X$, Cholesky factor it as $L L^T$, solve $L x = x_{white}$ to whiten new $x$ s. This is quadratic in the number of features.

2) diagonal: $x_{white} = x / \sigma(x)$ ignores cross-correlations completely.

One could get a tridiagonal matrix from $Fullcov$ just by zeroing all entries outside the tridiagonal, or not accumulating them in the first place. And here I start sinking: there must be a better approximation, perhaps hierarchical, block diagonal → tridiagonal ?

(Added 11 May): Let me split the question in two:

1) is there a fast approximate $cov(X)$ ?
No (whuber), one must look at all ${N \choose 2}$ pairs (or have structure, or sample).

2) given a $cov(X)$, how fast can one whiten new $x$ s ?
Well, factoring $cov = L L^T$, $L$ lower triangular, once, then solving $L x = x_{white}$ is pretty fast; scipy.linalg.solve_triangular, for example, uses Lapack.
I was looking for a yet faster whiten(), still looking.

• Do the columns have a natural ordering to them? Or do you want to find a tridiagonal approximation under some ("optimal") permutation of the columns? I'm assuming that when you say $A = \mathrm{Cov}(X)$ you're speaking of the covariance structure of the features. Can you confirm this? – cardinal May 10 '11 at 12:01
• No, there's no natural ordering, and yes, covariance of the 100 features. Methods that add up a full covariance matrix, then approximate it, would be >> O(size X); I'm looking for a fast simple approximation, which will necessarily be crude. – denis May 10 '11 at 12:07
• So, you want a tridiagonal approximation under some (to be determined by the data) permutation, yes? – cardinal May 10 '11 at 12:09
• added, tried to clarify. If a good (satisficing) permutation could be found in O(Nfeatures), yes, that would do. – denis May 10 '11 at 14:09
• There are approximations when the variables have additional structure, such as when they form a time series or realizations of a spatial stochastic process at various locations. These effectively rely on assumptions that let us relate the covariance between one pair of variables to that between other pairs of variables, such as between pairs separated by the same time lags. Calculations can be $O(Nf \log(Nf)$ in such cases. Absent such a model, I don't see how you can avoid computing all pairwise covariances. – whuber May 10 '11 at 16:10

## 2 Answers

Merely computing the covariance matrix--which you're going to need to get started in any event--is $O((Nf)^2)$ so, asymptotically in $N$, nothing is gained by choosing a $O(Nf)$ algorithm for the whitening.

There are approximations when the variables have additional structure, such as when they form a time series or realizations of a spatial stochastic process at various locations. These effectively rely on assumptions that let us relate the covariance between one pair of variables to that between other pairs of variables, such as between pairs separated by the same time lags. This is the conventional reason for assuming a process is stationary or intrinsically stationary, for instance. Calculations can be $O(Nflog(Nf)$ in such cases (e.g., using the Fast Fourier Transform as in Yao & Journel 1998). Absent such a model, I don't see how you can avoid computing all pairwise covariances.

On a whim, I decided to try computing (in R) the covariance matrix for a dataset of about the size mentioned in the OP:

z <- rnorm(1e8)
dim(z) <- c(1e6, 100)
vcv <- cov(z)


This took less than a minute in total, on a fairly generic laptop running Windows XP 32-bit. It probably took longer to generate z in the first place than to compute the matrix vcv. And R isn't particularly optimised for matrix operations out of the box.

Given this result, is speed that important? If N >> p, the time taken to compute your approximation is probably not going to be much less than to get the actual covariance matrix.