Calcul de la covariance en cours d'exécution (une passe)

Question

J'ai reçu un ensemble de vecteurs 3d (x, y, z), et je veux calculer la matrice de covariance sans stocker les vecteurs. Je vais le faire en C#, mais finalement je l'implémenterai en C sur un microcontrôleur, j'ai donc besoin de l'algorithme en lui-même, et non d'une bibliothèque.

Pseudocode serait aussi génial.

Answer 1

Je pense avoir trouvé la solution. Il est basé sur cet article sur how to calculate covariance manually et celui-ci sur calculating running variance. Et puis j'ai adapté l'algorithme dans ce dernier pour calculer la covariance au lieu de la variance, étant donné ma compréhension du premier article.

public class CovarianceMatrix 
{ 
    private int _n; 
    private Vector _oldMean, _newMean, 
        _oldVarianceSum, _newVarianceSum, 
        _oldCovarianceSum, _newCovarianceSum; 

    public void Push(Vector x) 
    { 
     _n++; 
     if (_n == 1) 
     { 
      _oldMean = _newMean = x; 
      _oldVarianceSum = new Vector(0, 0, 0); 
      _oldCovarianceSum = new Vector(0, 0, 0); 
     } 
     else 
     { 
      //_newM = _oldM + (x - _oldM)/_n; 
      _newMean = new Vector(
       _oldMean.X + (x.X - _oldMean.X)/_n, 
       _oldMean.Y + (x.Y - _oldMean.Y)/_n, 
       _oldMean.Z + (x.Z - _oldMean.Z)/_n); 

      //_newS = _oldS + (x - _oldM) * (x - _newM); 
      _newVarianceSum = new Vector(
       _oldVarianceSum.X + (x.X - _oldMean.X) * (x.X - _newMean.X), 
       _oldVarianceSum.Y + (x.Y - _oldMean.Y) * (x.Y - _newMean.Y), 
       _oldVarianceSum.Z + (x.Z - _oldMean.Z) * (x.Z - _newMean.Z)); 

      /* .X is X vs Y 
      * .Y is Y vs Z 
      * .Z is Z vs X 
      */ 
      _newCovarianceSum = new Vector(
       _oldCovarianceSum.X + (x.X - _oldMean.X) * (x.Y - _newMean.Y), 
       _oldCovarianceSum.Y + (x.Y - _oldMean.Y) * (x.Z - _newMean.Z), 
       _oldCovarianceSum.Z + (x.Z - _oldMean.Z) * (x.X - _newMean.X)); 

      // set up for next iteration 
      _oldMean = _newMean; 
      _oldVarianceSum = _newVarianceSum; 
     } 
    } 
    public int NumDataValues() 
    { 
     return _n; 
    } 

    public Vector Mean() 
    { 
     return (_n > 0) ? _newMean : new Vector(0, 0, 0); 
    } 

    public Vector Variance() 
    { 
     return _n <= 1 ? new Vector(0, 0, 0) : _newVarianceSum.DivideBy(_n - 1); 
    } 
} 

Answer 2

La formule est simple si vous avez Matrix et Vector classes à la main:

Vector mean; 
Matrix covariance; 
for (int i = 0; i < points.size(); ++i) { 
    Vector diff = points[i] - mean; 
    mean += diff/(i + 1); 
    covariance += diff * diff.transpose() * i/(i + 1); 
} 

Personnellement, je préfère toujours ce style plutôt que le calcul en deux passes. Le code est court et les résultats sont parfaits.

Matrix et Vector peuvent avoir une dimension fixe et peuvent être facilement codés à cette fin. Vous pouvez même réécrire le code en calculs discrets à virgule flottante et éviter de calculer la partie symétrique de la matrice de covariance.

Notez qu'il existe un vecteur produit externe sur la dernière ligne de code. Toutes les bibliothèques de vecteurs ne l'interprètent pas correctement.

Answer 3

Le code de l'émeu est élégant, mais nécessite une étape supplémentaire pour être correct:

Vector mean; 
Matrix covariance; 
for (int i = 0; i < points.size(); ++i) { 
    Vector diff = points[i] - mean; 
    mean += diff/(i + 1); 
    covariance += diff * diff.transpose() * i/(i + 1); 
} 

covariance = covariance/(points.size()-1); 

Notez la dernière étape de la normalisation de la covariance.