On the repeated inversion of a covariance matrix
M. de Jong





On the repeated inversion of a covariance matrix
    M. de Jong

  NWO-I, Nikhef, PO Box 41882, Amsterdam, 1098 DB Netherlands 

  Leiden University, Leiden Institute of Physics, PO Box 9504, Leiden, 2300 RA Netherlands

    
===========================================================================================================================================================================




  In many cases, the values of some model parameters are determined 
  by maximising the likelihood of a set of data points given the parameter values. 
  The presence of outliers in the data and correlations between data points
  complicate this procedure.
  An efficient procedure for the elimination of outliers is presented
  which takes the correlations between data points into account. 





§ INTRODUCTION


In general, the values of some model parameters are optimally determined 
by maximising the likelihood of a set of data points given the parameter values
(see e.g. <cit.>).
In this, a data point that has to low a probability to match the model can be considered an outlier.
The presence of such outliers in the data can readily be accommodated in
the probability density function of the data points.
It is, however, not straight forward to also take the correlations between the data points into account. 
If the underlying probability density functions are normal distributions, 
both the uncorrelated and correlated uncertainties of the data points 
can be incorporated in a single matrix.
The values of the model parameters can then be determined by minimising the χ^2:


    χ^2     =    ϵ^T  V^-1 ϵ


where
ϵ and ϵ^T are the (N × 1) and (1 × N) vectors
containing the distances between the model and the data points and
V is a (N × N) matrix.
The elements of V are set as follows.


    V_ii    =     (σ_i)^2                                
    
      V_ij    =    ∑_k∂ϵ_i/∂ u_k∂ϵ_j/∂ u_k(δ u_k)^2


where
σ_i refers to the uncorrelated uncertainty of data point i and 
u_k to some correlation parameter which itself has an an uncertainty δ u_k
The terms in the summation of equation <ref> are commonly referred
to as the covariances of the data points.
The matrix V is therefore often called the covariance matrix.
By construction, the matrix V is symmetric and can thus be inverted using an LDU decomposition.
The computation of V^-1 then requires 𝒪(N^3) operations. 
 


The presence of outliers cannot easily be incorporated in the covariance matrix.
So, it is desirable to remove the outliers and repeat the fit.
For this, a criterion is required to identify an outlier.
A common criterion is based on the value of the so-called standard deviation, D:


    D_k     ≡   |ϵ_k|/σ_k


A typical maximal allowed value of D is D_max  =   3
which –in the absence of correlations– corresponds to a probability to keep a good data point of about 0.997.
The removal of outliers could simply proceed by
removing the data point with the largest D and repeating the fit 
until there are no more data points with D > D_max. 
 


In case the correlations between the data points are strong,
this procedure may no longer be adequate.
In this scenario, the absolute values of some off-diagonal elements of V
are comparable to (or even larger than) the values of the diagonal elements.
The standard deviation is then no longer a good criterion because 
the distances and covariances of the other data points should also be taken into account.
A brute force procedure to identify outliers is to
1) remove data point k,
2) determine V',
3) invert V' and
4) minimise χ^2, 
for each data point k.
In this, V' corresponds to the (N-1)×(N-1) covariance matrix.
The change in χ^2 is then a good criterion to identify an outlier, i.e: 


    D_k     ≡   √(χ^2 - χ^2_k)


where χ^2_k corresponds to the χ^2 of the fit after removal of data point k.
As each inversion of V' requires 𝒪((N-1)^3) operations, 
this way of eliminating outliers requires 𝒪(N^4) operations.
For a large number of data points, the number of operations
needed to eliminate outliers may become too excessive.
 


An alternative exist, based on the fact that the removal of data point k
is equivalent to setting the corresponding uncertainty σ_k to infinity.
By doing so, the (N × N) covariance matrix V' can be decomposed as follows:


    V'     =      V   +   g  δ_k,k


In this, 
g is some arbitrary large value; in any case much larger than σ_k.
The matrix δ_i,j has 1 at row i and column j and 0 everywhere else.
To repeat the fit without data point k, it is required to invert matrix V'.
As a first step, the known inverse of the original matrix V is considered.
This is possible because V^-1 and V' have the same dimensions. 


    V' × V^-1      =         (V + g δ_k,k) × V^-1
    
          =         VV^-1 + g  δ_k,k V^-1
    
          =         I  +   g 
       ( [                      0                    3c⋯                      0;                      ⋮                                                                                           ⋮;                      0                    3c⋯                      0;               V^-1_k,1                      ⋯               V^-1_k,k                      ⋯ V^-1_k,N[0cm][l]←row k;                      0                    3c⋯                      0;                      ⋮                                                                                           ⋮;                      0                    3c⋯                      0;                        ])                                                            
    
       ≡     A


where I is the identity matrix.
It is obvious that the inverse of matrix V' is equal to the product V^-1 A^-1.
Consequently, the problem of inverting V' is reduced to inverting matrix A.
The inverse of matrix A is trivial, namely:


    A^-1      =         I  - g/1 + g V^-1_k,k( [                      0                    3c⋯                      0;                      ⋮                                                                                           ⋮;                      0                    3c⋯                      0;               V^-1_k,1                      ⋯               V^-1_k,k                      ⋯ V^-1_k,N[0cm][l]←row k;                      0                    3c⋯                      0;                      ⋮                                                                                           ⋮;                      0                    3c⋯                      0;                        ])


Now, one can let go g →∞ and multiply V^-1 and A^-1
to obtain the inverse of matrix V'.
As can be seem from equation <ref>,
the product of V^-1 and A^-1 requires 𝒪(N^2) operations.
Hence, this way of eliminating outliers only requires 𝒪(N^3) operations
which actually is equal to the number of operations needed to invert the covariance matrix for the first fit.
It is interesting to note that no additional memory is required for the computation of V'^-1.
 


A further reduction in the number of operations is possible when the outcome of the first fit is retained.
In that case, the distance ϵ_i between the model and data point i stays the same.
The χ^2 without data point k can then be expressed as:


    χ^2_k    =    ∑_i=1^N∑_j=1^N ϵ_i   V'^-1_i,j ϵ_j                   
    
        =    ∑_i=1^N∑_j=1^N ϵ_i   (V^-1A^-1)_i,j ϵ_j            
    
        =    χ^2 - 1/V^-1_k,k(∑_j=1^n V^-1_k,j ϵ_j)^2


As can be seen from equation <ref>,
the distances and covariances of the other data points are taken into account
but the inverse of V' is no longer required.
As a result, the summation takes only 𝒪(N) operations and
the elimination of outliers 𝒪(N^2) operations.




§ CONCLUSIONS


An efficient method is presented to eliminate outliers from a set of data points 
which takes the correlations between the data points into account.
The number of operations needed for this procedure is 
the same as that needed for the inversion of the covariance matrix.


plain


