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.
 A: 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(Nf\,\log(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.
A: 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.
A: Extra two cents:
Algorithmically speaking, I don' think there are any faster algorithms to do this for generic $X$. If there were, they must have been already implemented in the programs so far. However, from a software-engineering perspective, speeds can differ dramatically between implementations (e.g., legacy Blas, Goto Blas, Intel's MKL, and OpenBlas).
Depending on the application scenario, you can code your implementations specific to your use cases, but this requires do code-juggling in a more native language such as Fortran, C, and C++.  Now to do fast implementations, one thing definitely needed is to take advantage of new CPU features such as AVX512, which requires bits of ASM knowledge.
Also, depending on how crude you mean, one possible approximation is to use Monte Carlo integration.  For example, for a  1000000 x 1 column $x$ (say, 1000000 observations of something), the literal computation for mean will be $\sum_{i=1}^{1000000}x_i/1000000$, but for an approximation, you can take a subsmaple of only 1000 to get a crude estimation of mean, say $\sum_{j=1}^{1000}x_{i(j)}/1000$, where $x_{i(j)}$ is a sub sample of the original x. The idea applies to computing your covariance. Say, you have two columns x and y; rather than doing $\sum_{i=1}^{1000000}x_i*y_i$, you can approximate it with only a subsample (e.g, 1000) $\sum_{j=1}^{1000}x_{i(j)}*y_{i(j)}/1000 *1000000 $. If the ordering of your rows are randomly enough, you don't need to explicitly sample these 1000 indices, you may take it by a regular step or simply take the first 1000 rows.
For your back-substitution to solve $Lx=x_w$, if L is re-used many times, one bit of minor improvement is to explicitly store the diagonal elements of $L$ as their inversion (e.g., $1/Lii$); this will avoid doing the division but rather doing multiplication in the back-substitution stage. Nowadays, this saving is of no real sigfifance, given the speed for float division and multiplication should be roughly the same. But that is not the case for old CPUS because division is used to be much slower than multiplication. Again, this is just an minor thing and wouldn't expedite  the computation dramatically, but depending on user cases, it may help a little bit.
