# Contrasting covariance calculation using R, Matlab, Pandas, NumPy cov, NumPy linalg.svd

I wish there were a kind of Rosetta Stone for covariance (and other) stat methods supplied by popular statistics libraries and software, so that we could understand why what is ostensibly the same method gives different results in different libraries. I am aware that most of the discrepancies are due to scaling and normalizing factors, but these factors are usually not explained in the documentation sets of these libraries, and some of the discrepancies are pretty significant. If approaches to normalization and scaling are widely understood, why is it so difficult to write Python code that manually reproduces the output of these methods? Shouldn't the algorithms used in these popular libraries be more transparent?

Test example 1.

Input matrix A = $\begin{pmatrix} 1 & 2 \\ 3 & 4 \end{pmatrix}$

1. Output of cov(A) in R: $\begin{pmatrix} 2 & 2 \\ 2 & 2 \end{pmatrix}$

2. Output in R when input is $A^\top$: $\begin{pmatrix} 0.5 & 0.5 \\ 0.5 & 0.5 \end{pmatrix}$

3. Output of cov(A) in Matlab: $\begin{pmatrix} 2 & 2 \\ 2 & 2 \end{pmatrix}$

4. Output in Pandas: $\begin{pmatrix} 2 & 2 \\ 2 & 2 \end{pmatrix}$

5. Output in SciKit Learn Empirical Covariance: $\begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix}$
What sort of normalization or scaling is SciKit doing here?

6. Output in NumPy using np.cov(): $\begin{pmatrix} 0.5 & 0.5 \\ 0.5 & 0.5 \end{pmatrix}$

7. Output from Python code that supposedly reflects what np.cov() is doing (according to this forum post): $\begin{pmatrix} +2 & -2 \\ -2 & +2 \end{pmatrix}$
NOTE the sign difference in the off-diagonal elements!

8. Output in NumPy using np.linalg.svd() to calculate covariance: $\begin{pmatrix} 10 & -14 \\ -14 & 20 \end{pmatrix}$
The values here differ from Matlab by more than a constant factor or a square.

Code used to generate input matrices for Test 1:

1. R code for A = matrix( + c(1, 3, 2, 4), + nrow=2, + ncol=2)

2. R code for input of A.transpose = matrix(+ c(1, 2, 3, 4), + nrow=2, + ncol=2)

3. Matlab code for input A = [1, 2; 3, 4]

4. Panda code: DataFrame.cov(np.vstack(([1, 2], [3, 4])))

5. SciKit Learn empirical covariance code:

X = np.vstack(([1, 2], [3, 4]))
fit = EmpiricalCovariance().fit(X)
fit.covariance_

6. NumPy code using np.cov():

X = np.vstack(([1, 2], [3, 4]))
np.cov(X)

7. Python code that supposedly reflects what's happening in np.cov():

X0 = np.vstack(([1, 2], [3, 4]))
X = X0 - X0.mean(axis=0)
N = X.shape[1]
fact = float(N - 1)
C = np.dot(X, X.T) / fact


NOTE that putting the transpose X.T first instead of second in np.dot(X.T, X) gives the same result as Matlab and R without the negative signs in the off-diagonal elements, but this is supposed to be an outer product!! What is backwards about the NumPy conventions such that we need to put the transpose term first in order to get an outer product?

1. NumPy code to obtain covariance using np.linalg.svd() according to this post:

X0 = np.vstack(([1, 2], [3, 4]))
U, s, V = np.linalg.svd(X0, full_matrices = 0)
D = np.dot(np.dot(V,np.diag(s**2)),V.T)
Dadjust = D / (X0.shape[0] - 1)
print (Dadjust)


Again, the values here differ from Matlab by more than a constant factor or a square. Can anyone help me unravel this?

Test example 2.

Input matrix A = $\begin{pmatrix} 1 & 2 \\ 3 & 4 \\ 22 & 44 \end{pmatrix}$

1. Output of cov(A) in R: $\begin{pmatrix} 134.3333 & 274.3333 \\ 274.3333 & 561.3333 \end{pmatrix}$

2. Output in R when input is $A^\top$: $\begin{pmatrix} 0.5 & 0.5 & 11 \\ 0.5 & 0.5 & 11 \\ 11.0 & 11.0 & 242 \end{pmatrix}$

3. Output in Matlab: $\begin{pmatrix} 134.3333 & 274.3333 \\ 274.3333 & 561.3333 \end{pmatrix}$

4. Output in Pandas: $\begin{pmatrix} 134.3333 & 274.3333 \\ 274.3333 & 561.3333 \end{pmatrix}$

5. Output in SciKit Learn Empirical Covariance: $\begin{pmatrix} 89.55555556 & 182.88888889 \\ 182.88888889 & 374.22222222 \end{pmatrix}$

6. Output in NumPy using np.cov(): $\begin{pmatrix} 0.5 & 0.5 & 11 \\ 0.5 & 0.5 & 11 \\ 11 & 11 & 242 \end{pmatrix}$

7. Output from Python code that supposedly reflects what np.cov() is doing (according to this forum post): $\begin{pmatrix} 273.88888889 & 229.22222222 & -503.11111111 \\ 229.22222222 & 192.55555556 & -421.77777778 \\ -503.11111111 & -421.77777778 & 924.88888889 \end{pmatrix}$

8. Output in NumPy using np.linalg.svd() to calculate covariance: $\begin{pmatrix} 247 & 491 \\ 491 & 978 \end{pmatrix}$

Code used to generate input matrices for Test 2:

1. R code for A = matrix(+ c(1, 3, 22, 2, 4, 44), + nrow=3, + ncol=2)

2. R code for covariance with input of $A^\top$:

Atranspose = t(A)
cov(Atranspose)

3. Matlab code for input A = [1, 2; 3, 4; 22, 44]

4. Panda code: DataFrame.cov(np.vstack(([1, 2], [3, 4], [22, 44])))

5. SciKit Learn empirical covariance code:

X = np.vstack(([1, 2], [3, 4], [22, 44])
fit = EmpiricalCovariance().fit(X)
fit.covariance_

6. NumPy code using np.cov():

X = np.vstack(([1, 2], [3, 4], [22, 44])
np.cov(X)

7. Python code that supposedly reflects what's happening in np.cov():

X0 = np.vstack(([1, 2], [3, 4], [22, 44])
X = X0 - X0.mean(axis=0)
N = X.shape[1]
fact = float(N - 1)
C = np.dot(X, X.T) / fact

8. NumPy code to obtain covariance using np.linalg.svd() according to this post:

X0 = np.vstack(([1, 2], [3, 4], [22, 44])
U, s, V = np.linalg.svd(X0, full_matrices = 0)
D = np.dot(np.dot(V,np.diag(s**2)),V.T)
Dadjust = D / (X0.shape[0] - 1)
print (Dadjust)

• People speak different language, and different software have different default setting and syntax. Why that is a problem? And could you elaborate your question? Are you asking why there are "inconsistencies" among software? Feb 22, 2017 at 21:47
• Welcome to Cross Validated! @hxd1011, I'm not sure this is a problem, per se, but it would be useful to know why two superficially-similar programs do different things. Some of these might be technical (row vs. columns) and others appear to be related to the underlying math, but if you don't know, it can definitely be confusing. Feb 22, 2017 at 22:11
• Octave also gives the same results as Matlab with the function cov(A). To answer your question about the difficulty in implementing the function into Python, it may be helpful to look at the code the other packages use to make the useful function cov(). As Octave and R are open source, it is relatively easy to look at the code behind these functions. Feb 22, 2017 at 22:24
• This seems potentially on topic here, potentially not. I'm not sure. I think this is ambiguous. Feb 23, 2017 at 16:27
• (7) and (8) clearly are programming errors: if you would like an explanation, please show details of the calculations. The rest of the differences are readily understood by referring to the documentation. The R documentation alone has enough information to understand why such differences would crop up and even to predict them quantitatively.
– whuber
Feb 24, 2017 at 0:46

## 1 Answer

Note that numpy.cov() considers its input data matrix to have observations in each column, and variables in each row, so to get numpy.cov() to return what other packages do, you have to pass the transpose of the data matrix to numpy.cov().

The Python code that you linked can be used to simulate what other packages do, but it contains some errors: N should be the number of rows, not columns, and you have to perform the matrix multiplication in the other order:

import numpy as np

def cov(X0):
print "\n==\nMatrix:"
print X0
X = X0 - X0.mean(axis=0)
N = X.shape[0]                # !!!
fact = float(N - 1)
print "Covariance:"
print np.dot(X.T, X) / fact   # !!!

X0 = np.vstack(([1, 2], [3, 4]))
cov(X0)
cov(X0.T)

X0 = np.vstack(([1, 2], [3, 4], [22, 44]))
cov(X0)
cov(X0.T)


With these fixes, the covariance behaves as expected:

==
Matrix:
[[1 2]
[3 4]]
Covariance:
[[ 2.  2.]
[ 2.  2.]]

==
Matrix:
[[1 3]
[2 4]]
Covariance:
[[ 0.5  0.5]
[ 0.5  0.5]]

==
Matrix:
[[ 1  2]
[ 3  4]
[22 44]]
Covariance:
[[ 134.33333333  274.33333333]
[ 274.33333333  561.33333333]]

==
Matrix:
[[ 1  3 22]
[ 2  4 44]]
Covariance:
[[   0.5    0.5   11. ]
[   0.5    0.5   11. ]
[  11.    11.   242. ]]


As for the numpy.linalg.svd() code, you need to center the data matrix by subtracting off the variable means, and the multiplication involving the V matrix must be performed in the other order. With these changes you will replicate everybody else's behavior:

import numpy as np

def cov(X0):
print "\n==\nMatrix:"
print X0
X = X0 - X0.mean(axis=0)
U, s, V = np.linalg.svd(X, full_matrices = 0)
D = np.dot(np.dot(V.T,np.diag(s**2)),V)
Dadjust = D / (X0.shape[0] - 1)
print "Covariance:"
print (Dadjust)

• Key difference between np.cov() and the corresponding cov() methods in Matlab, R and Pandas is, as you say, NumPy cov() considers rows to be observations instead of columns. You are right, thank you. I hope we'll get an answer shedding light on why SciKit Learn's Empirical Covariance method returns different values or why the SVD likewise fails to produce the results promised by the post I linked to. Feb 22, 2017 at 22:40
• @JohnStrong I haven't looked at the SVD approach, but the reason why SciKit Learn returns a different result is that the Empirical Covariance is dividing by N instead of N-1. You'll get everybody else's covariance if you multiply the output from sklearn by $N/(N-1)$. Feb 22, 2017 at 22:50
• Thank you for clearing up what the empirical variance routine is doing, grand_chat!! Feb 23, 2017 at 13:32
• Just noticed the follow-up explaining what I was missing with regard to the little SVD algorithm. Very much appreciated. Thanks for taking time. Feb 23, 2017 at 13:35