<jupyter_start><jupyter_text>Covariance Non Trivial Example $ n $ - number of points$ m $ - number of dimensions$ h, i $ - indices goes over the rows of $X$ $ j, k$ - indices used to go over the columns of $X$$ X = [ x_1, x_2, ..., x_j, ... , x_m ] $What we want is the covariance matrix of $ X^T $ $$ \mu_j = \frac{1}{n} \sum^n_{i=0} x_{ij} $$ $$ C_{jk} = \frac{1}{n-1} \sum^n_{i=0} (x_{ij} - \mu_j)^T(x_{ik} - \mu_k)$$ If we were to expand $C_{jk} $ we could separate out different components$$ C_{jk} = \frac{1}{n-1} \left[ \sum^n_{i=0}{x_{ij} x_{ik}} - \mu_j \left(\sum^n_{i=0} {x_{ik}} \right) - \left( \sum^n_{i=0} {x_{ij}} \right) \mu_k  + n \mu_j \mu_k \right] $$ $$ C_{jk} = \frac{1}{n-1} \left[\sum_{i=0}^n x_{ij}x_{ji} - \mu_j ( n \mu_k ) - ( n \mu_j ) \mu_{k} + n \mu_{j} \mu_{k} \right] $$ $$ C_{ij} = \frac{1}{n-1} \left[ \sum^n_{i=0} x_{ij} x_{ji} - n \mu_j \mu_i \right] $$ If we substitute for the summations$$ A_{jk} = \sum^n_{i=0} {x_{ij} x_{ik}}$$We can write $C_{ij}$ as:$$C_{ij} = \frac{1}{n-1}\left[ A_{jk} - n\mu_j \mu_i \right]$$<jupyter_code># Here we are now going to show how to calculate C_ij for a data set
import numpy as np

X = np.array([
      [0.6787,    0.6948,    0.7094],
      [0.7577,    0.3171,    0.7547],
      [0.7431,    0.9502,    0.2760],
      [0.3922,    0.0344,    0.6797],
      [0.6555,    0.4387,    0.6551],
      [0.1712,    0.3816,    0.1626]])

# Where m is the number of dimensions or cols
m = 3

# Where n is the number of points or rows
n = 6
# We will start by calculating mu or the mean

Mu = np.mean(X,axis=0);
print(Mu)
# Now we can calculate C_ij

# Here we explicitly calculate the covariance matrix
Cov = np.zeros((m,m))
for j in range(0,m):
    for k in range(j,m):
        A_jk = 0.0
        for i in range(0,n):
            A_jk = A_jk + X[i][j]*X[i][k]
        
        Cov[j,k] = 1.0/(n-1) * ( A_jk - n*Mu[j]*Mu[k])
        Cov[k,j] = Cov[j,k]

print("\nUsing our algorithm\n")
print(Cov)

# To be safe we will compare with numpy builtin
print("\nUsing numpy\n")
print(np.cov(X.transpose()))<jupyter_output>Using our algorithm

[[ 0.05497947  0.037775    0.02970302]
 [ 0.037775    0.10060908 -0.03052289]
 [ 0.02970302 -0.03052289  0.06393645]]

Using numpy

[[ 0.05497947  0.037775    0.02970302]
 [ 0.037775    0.10060908 -0.03052289]
 [ 0.02970302 -0.03052289  0.06393645]]<jupyter_text>The Question remains how do we correctly append data to our covariance matrix.<jupyter_code># We will add the following new set of data

X_new = np.array([[0.0318, 0.7952, 0.4984],
   [0.2769, 0.1869, 0.9597],
   [0.0462, 0.4898, 0.3404],
   [0.0971, 0.4456, 0.5853],
   [0.8235, 0.6463, 0.2238]])

# This is 5 new points
n_new = 5

# For a total data set of 11 points
n_full  = n_new + n<jupyter_output><empty_output><jupyter_text>We will represent the new data with a $'$ and the full data set with a $\dagger$ symbol.We can start by calcuating a new mean for the total data set:$$ \mu^{\dagger}_j = \frac{1}{n^{\dagger}} \sum^{n^{\dagger}}_{i=0} x^{\dagger}_{ij} $$But remember we are assuming we don't have the full data set anymore just the covariance matrix and mean of the old data set and the also the new data set, so we can rewrite this as:$$ \mu^{\dagger}_{j} = \frac{1}{n^{\dagger}} \left( \mu_j n + \sum_{i=0}^{n'} x_{ij} \right) $$<jupyter_code>Mu_full = 1/(n_full) * ( Mu*n + np.sum(X_new,axis=0))

print(Mu_full)<jupyter_output>[0.4249     0.48914545 0.53137273]<jupyter_text>The next step will be to adjust $ \mathbf{A_{jk}} $$$ A^\dagger_{jk} = A_{jk} +  \sum^{n'}_{i=0} {x'_{ij} x'_{ik}} $$We can get $A_{jk}$ from the old data set by reversing $C_{jk}$ equation:$$ A_{jk} = C_{jk}\cdot(n-1) + n(\mu_j\mu_k) $$<jupyter_code># A more simplified version we can get rid of B entirely

C_full = np.zeros((m,m))
for j in range(0,m):
    for k in range(0,m):
        sum = 0.0
        for i in range(0,n_new):
            sum = sum + X_new[i][j]*X_new[i][k]
        
        A_jk = Cov[j][k]*(n-1) + n*Mu[j]*Mu[k]
        A_jk_full = A_jk + sum
        C_full[j,k] = 1.0/(n_full-1) * ( A_jk_full - Mu_full[j]*Mu_full[k]*n_full)
        
print("Our algorithm\n")
print(C_full)

# And comparing with what numpy would do if it had the full data set we have

X_full = np.concatenate((X, X_new))
print("\nFull dataset\n")
print(X_full)
print("\nNumpy solution\n")
print(np.cov(X_full.transpose()))<jupyter_output>Our algorithm

[[ 0.0981211   0.01732581  0.00371005]
 [ 0.01732581  0.07169848 -0.0343945 ]
 [ 0.00371005 -0.0343945   0.06386179]]

Full dataset

[[0.6787 0.6948 0.7094]
 [0.7577 0.3171 0.7547]
 [0.7431 0.9502 0.276 ]
 [0.3922 0.0344 0.6797]
 [0.6555 0.4387 0.6551]
 [0.1712 0.3816 0.1626]
 [0.0318 0.7952 0.4984]
 [0.2769 0.1869 0.9597]
 [0.0462 0.4898 0.3404]
 [0.0971 0.4456 0.5853]
 [0.8235 0.6463 0.2238]]

Numpy solution

[[ 0.0981211   0.01732581  0.00371006]
 [ 0.01732581  0.07169848 -0.0343945 ]
 [ 0.00371006 -0.0343945   0.06386179]]
{'fraction_non_alphanumeric': 0.1841411478401947, 'fraction_numerical': 0.14479821537213547, 'mean_word_length': 3.0760330578512396, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '559848', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2546, 'n_tokens_neox': 2161, 'n_words': 569, 'path': 'jupyter_notebooks/Covariance/UpdatingCovar.ipynb', 'repo_name': 'lanl/PANACEA', 'size': '4931'}