Iteratively updating the decomposition of a covariance matrix

Question

Say I have a set of data points in $\mathbb{R}^d$. At each time step, I add a point to the set, or remove one. I then fit a multivariate normal distribution to the points in the set, and compute the entropy of this distribution. I also run some other computations, which involve evaluating the density at various points.

The naive way to do this would be to re-compute the mean and covariance matrix from scratch at each time step. Of course, this would be computationally wasteful. An improvement would be to maintain running sums, and use these to efficiently compute the mean and covariance matrix at each time step. But, computing the entropy and the density requires a decomposition of the covariance matrix (e.g. eigendecomposition, cholesky decompositon, etc.). The 'improved' approach would still require re-computing this decomposition at each time step.

What's needed is a way to iteratively update the decomposition of the covariance matrix, given the data point added to or removed from the set. Can you recommend a good way to do this?

If the data points had constant mean, this problem could be formulated in terms of rank one updates. There are some papers describing this kind of approach. Unfortunately, this isn't the case; adding/removing a point changes both the mean and covariance matrix.

Answer 1

I think you can achieve part of what you want by using an incremental SVD and/or an online PCA algorithm. Given a known decomposition we update it to take into account a new data-point.

In terms of theoretical background, I would suggest you look at computer vision literature. They have nice application papers that explain the relevant mechanics without getting to heavy into the numerics side of things. The papers "Incremental singular value decomposition of uncertain data with missing values" (Brand, 2002) and "Incremental learning for robust visual tracking" (Ross et al. 2008) are well-cited papers than can offer a good starting point. The whole idea is the "forgetting vector" f encapsulates how much we value our recent points compared to the ones we have already observed. Setting f to something very small will lead us to obtain essentially a new eigen-decomposition for $n+1$ points given we already have seen $n$ points.

In terms of implementations in R I would suggest you look at the package onlinePCA. It offers an online PCA implementation using incremental SVD (function incRpca); it also offers the updateCovariance and updateMean functions if you want to take matters in your own hands more (for PCA-purposes incremental SVD will be much faster though).

The second part of the task as hand is removing the influence of existing data-points from the decomposition already computed. In that case the procedure we need is decremental SVD or SVD downdating.

In terms of theoretical background, decremental SVD is a much less popular problem than incremental SVD and I have no experience on it. Starting references appear to be the papers: "Merging and splitting eigenspace models" (Hall et al. 2000) and "Efficiently downdating, composing and splitting singular value decompositions preserving the mean information" (Melenchon and Martinez, 2007); I came across a paper that implements among other things decremental LDA - Incremental and decremental LDA learning with applications (Pang et al. 2010) which might prove useful too.

I have not come across any decremental SVD routines. For a particular problem if there are at most one or two consecutive down-dates regarding the most recent points (eg. we just read in corrupted points), we might even get away by simply storing the past couple decompositions and retrieving them. The Melenchon and Martinez paper outlines a relatively simple implementation if you are inclined to try your hand (in which case please publish it!).

Update on the matter: I looked around the Eigen library offers a rank update for their $LDL^T$ decomposition (ie. a variant of the standard Cholesky decomposition) and it allows for down-dates. That means we can remove from a decomposition a previously-added vector. Eigen can be very easily linked with R through the excellent RcppEigen package. If it possible to work with the Cholesky decomposition of the covariance matrix this functionality might be fully resolve the issue presented.