Variance-Covariance Matrix¶

Miftahur Rahman, Ph.D.
March 26, 2022

The motivation of writing a quick review article on "$\bf Variance-Covariance \ Matrix$" is to study the protection of active satellites and International Space Station (ISS) from impacts by space debris objects. Collision avoidance for the orbiting satellites and ISS serves both adaptation and mitigation. A major collision such as the Iridium-Cosmos collision in early 2009 created much additional debris and thus, as the debris population increases collision avoidance becomes urgently necessary. The Joint Space Operations Center of the U.S. conducting conjunction assessments for all operational spacecraft in Earth's orbits, regardless of ownership. Any prediction of a close approach, typically within 1 km, will be notified with the spacecraft owner or the operator immediately without any cost.

Predictions of the collision in between of the satellites, satellite and debris are probabilistic due to inherent uncertainties in space surveillance measurements, the dynamic state of the atmosphere and, the instability of at least one of the conjuncting objects. Typical maneuver probability thresholds for collision avoidance are 1 in 10,000 for human space flight and 1 in 1000 (or more) for robotic satellites.

Collision avoidance maneuvers normally integrate well with standard satellite operations. Collision avoidance maneuvers are typically very small, that is, involve changes in velocity of less than $1 \frac {m}{s}$ and in most cases can be conducted in a manner that does not waste propellant resources. As an example, collision avoidance maneuvers performed by the ISS and robotic satellites are almost always result in a small increase in orbital altitude and thus simply constitute an unscheduled antidrag maneuver.

More than 99% of the risk to operational spacecraft comes from collisions with objects too small to track routinely, that is, objects smaller than 5–10 cm. Therefore, there is no alternative to an improvement in the orbital debris environment that could drastically reduce collision risks for operational spacecraft.

Covariance matrices, are compact summaries of the variability and covariation present in data. The covariance matrix and vector of means constitute sufficient statistics for models that assume a multivariate normal distribution.

$\bf Variance$ is a measure of the variability or spread in a set of data. Mathematically, it is the average squared deviation from the mean data. We use the following formula to compute population variance.

$$\sigma_x^2=Var(X) = \sum_{i=0}^{N} \frac{(X_i-\bar X)^2}{N}=\sum \frac{x_i^2}{N}$$

where

N is the total number of data points in a data set,
$\bar X$ is the mean of the N number of data set,
$X_i$ is the ith raw data in the set of data,
$x_i$ is the ith deviation of data in the set of data,
Var(X) is the variance of all the data points in the data set.

$\bf Covariance$ is a measure of the extent to which corresponding elements from two sets of ordered data move in the same direction. We use the following formula to compute population covariance.

$$\sigma(x,y)=Cov(X, Y) = \sum \frac { (X_i - \bar X ) ( Y_i - \bar Y )}{N} = \sum \frac{x_iy_i}{N}$$

where

N is the number of data points in each set of data,
$\bar X$ is the mean of the N data points in the first data set,
$X_i$ is the ith raw data in the first set of data,
$x_i$ is the ith deviation of data in the first set of data,
$\bar Y$ is the mean of the N data in the second data set,
$Y_i$ is the ith raw data in the second set of data,
$y_i$ is the ith deviation data in the second set of data,
Cov(X,Y) is the covariance of corresponding data in the two sets of data.

Variance-Covariance Matrix¶

$\bf Variance \ and \ covariance$ are often displayed together in a variance-covariance matrix. The variances appear along the diagonal and covariances appear in the off-diagonal elements, as shown below.

$$\begin{equation*} \bf V = \begin{bmatrix} \sum \frac{x_1^2}{N} & \sum \frac{x_1x_2}{N} & ..... \sum \frac{x_1x_r}{N}\\ \sum \frac{x_2x_1}{N} & \sum \frac{x_2^2}{N} & ..... \sum \frac{x_2x_r}{N}\\ .....\\ \sum \frac{x_rx_1}{N} & \sum \frac{x_rx_2}{N} & ..... \sum \frac{x_r^2}{N}\\ \end{bmatrix} \end{equation*}$$

where

$\bf V$ is a r x r variance-covariance matrix,
N is the number of data points in each of the r data sets,
$x_i$ is a deviation data from the ith data set,
$\sum \frac{x_i^2}{N}$ is the variance of elements from the ith data set,
$\sum \frac{x_ix_j}{N}$ is the covariance for elements from the ith and jth data sets.

Example of Variance-Covariance Matrix¶

Suppose $\bf A$ is an n x k matrix holding ordered sets of raw data of the scores on k tests for n students, as shown in Example-1.

Starting with the raw data of matrix $\bf A$, you can create a variance-covariance matrix to show the variance within each column and the covariance between columns. Here's how.

Transform the raw scores from matrix X into deviation scores for matrix x.

$\bf a$ = $\bf A$ - $\bf 1$ $\bf 1'$ $\bf A$ ( 1 / n )

where

1 is an n x 1 column vector of ones 1' is the transpose of vector 1
$\bf a$ is an n x k matrix of deviation scores: $a_{11}$, $a_{12}$, . . . , $a_{nk}$
$\bf A$ is an n x k matrix of raw scores: $A_{11}$, $A_{12}$, . . . , $A_{nk}$

Compute $\bf a'a$, the k x k deviation sums of squares and cross products matrix for $\bf a$.

Then, divide each term in the deviation sums of squares and cross product matrix by n to create the variance-covariance matrix. That is,

$$\bf V = \bf a'a ( 1 / n )$$

where

$\bf V$ is a k x k variance-covariance matrix, $\bf a'$ is the transpose of matrix $\bf a$, $\bf a'a$ is the deviation sums of squares and cross product matrix, n is the number of scores in each column of the original matrix $\bf A$,

Example-1¶

The table below displays scores on PHY, MAT, and ENG tests for 5 students.

Student	PHY	MAT	ENG
1	90	60	90
2	90	90	30
3	60	60	60
4	60	60	90
5	30	30	30

Note that data from the table can be represented in matrix $\bf A$, where each column in the matrix shows scores on a test and each row shows scores for a student.

$$\begin{equation*} \bf A = \begin{bmatrix} 90 & 60 & 90\\ 90 & 90 & 30\\ 60 & 60 & 60\\ 60 & 60 & 90\\ 30 & 30 & 30\\ \end{bmatrix} \end{equation*}$$

Given the data represented in matrix $\bf A$, compute the variance of each test and the covariance between the tests.

The solution involves a three-step process.

First, we transform the raw scores in matrix $\bf A$ to deviation scores in matrix $\bf a$, using the transformation formula described at how to transform raw scores to deviation scores.

$$\bf a = \bf A - 11'\bf A \left (\frac{1}{n}\right)$$

where

1 is an 5 x 1 column vector of ones,
$\bf a$ is an 5 x 3 matrix of deviation scores: $a_{11}$, $a_{12}$, . . . , $a_{53}$
$\bf A$ is an 5 x 3 matrix of raw scores: $A_{11}$, $A_{12}$, . . . , $A_{53}$
n is the number of rows in matrix $\bf A$

$$\begin{equation*} \bf a = \begin{bmatrix} 90 & 60 & 90\\ 90 & 90 & 30\\ 60 & 60 & 60\\ 60 & 60 & 90\\ 30 & 30 & 30\\ \end{bmatrix}- \begin{bmatrix} 1 & 1 & 1 & 1 & 1\\ 1 & 1 & 1 & 1 & 1\\ 1 & 1 & 1 & 1 & 1\\ 1 & 1 & 1 & 1 & 1\\ 1 & 1 & 1 & 1 & 1\\ \end{bmatrix} \begin{bmatrix} 90 & 60 & 90\\ 90 & 90 & 30\\ 60 & 60 & 60\\ 60 & 60 & 90\\ 30 & 30 & 30\\ \end{bmatrix} \left (\frac{1}{5}\right) \end{equation*}$$ $$\begin{equation*} \bf a = \begin{bmatrix} 90 & 60 & 90\\ 90 & 90 & 30\\ 60 & 60 & 60\\ 60 & 60 & 90\\ 30 & 30 & 30\\ \end{bmatrix}- \begin{bmatrix} 66 & 60 & 60\\ 66 & 60 & 60\\ 66 & 60 & 60\\ 66 & 60 & 60\\ 66 & 60 & 60\\ \end{bmatrix} \end{equation*}$$ $$\begin{equation*} \bf a = \begin{bmatrix} 24 & 0 & 30\\ 24 & 30 & -30\\ -6 & 0 & 0\\ -6 & 0 & 30\\ -36 & -30 & -30\\ \end{bmatrix} \end{equation*}$$

Then, to find the deviation score sums of squares matrix, we compute a'a, as shown below.

$$\begin{equation*} \bf a'a = \begin{bmatrix} 24 & 24 & -6 & -6 & -36\\ 0 & 30 & 0 & 0 & -30\\ 30 & -30 & 0 & 30 & -30\\ \end{bmatrix}- \begin{bmatrix} 66 & 60 & 60\\ 66 & 60 & 60\\ 66 & 60 & 60\\ 66 & 60 & 60\\ 66 & 60 & 60\\ \end{bmatrix} \end{equation*}$$ $$\begin{equation*} \bf a'a = \begin{bmatrix} 2520 & 1800 & 900\\ 1800 & 1800 & 0\\ 900 & 0 & 3600\\ \end{bmatrix} \end{equation*}$$

And finally, to create the variance-covariance matrix, we divide each element in the deviation sum of squares matrix by n, as shown below.

$$\begin{equation*} \bf V = \frac{a'a}{5} = \begin{bmatrix} 2520 & 1800 & 900\\ 1800 & 1800 & 0\\ 900 & 0 & 3600\\ \end{bmatrix} \left (\frac{1}{5}\right) \end{equation*}$$ $$\begin{equation*} \bf V = \begin{bmatrix} 504 & 360 & 180\\ 360 & 360 & 0\\ 180 & 0 & 720\\ \end{bmatrix} \end{equation*}$$

We can interpret the variance and covariance statistics in matrix $\bf V$ to understand how the various test scores vary and covary.

Shown along the diagonal, we see the variance of scores for each test. The English test has the biggest variance (720); and the MAT test, the smallest (360). So we can say that English test scores are more variable than MAT test scores. The covariance is displayed in the off-diagonal elements of matrix $\bf V$.

The covariance between PHY and MAT is positive (360), and the covariance between PHY and ENG is positive (180). This means the scores tend to covary in a positive way. As scores on PHY go up, scores on ENG and MAT also tend to go up; and vice versa.

The covariance between MAT and ENG, however, is zero. This means there tends to be no predictable relationship between the movement of MAT and ENG scores. If the covariance between any tests had been negative, it would have meant that the test scores on those tests tend to move in opposite directions. That is, students with relatively high scores on the first test would tend to have relatively low scores on the second test.

Steps to Create a Covariance Matrix using Python¶

Now we are going to test the Example-1 using Python libraries.

Step 1: Using the Same Dataset¶

Student	PHY	MAT	ENG
1	90	60	90
2	90	90	30
3	60	60	60
4	60	60	90
5	30	30	30

Step 2: Get the Population Covariance Matrix using Python¶

To get the population covariance matrix (based on N), you’ll need to set the bias to True in the code below.

This is the complete Python code to derive the population covariance matrix using the numpy package:

Step 3: Get a Visual Representation of the Matrix¶

You can use the seaborn and matplotlib packages in order to visually represent the covariance matrix.

Here is the complete code that you can apply in Python:

Derive the Sample Covariance Matrix¶

To get the sample covariance (based on N-1), you’ll need to set the bias to False in the code below.

Here is the code based on the numpy package:

You can also use the pandas package in order to get the sample covariance matrix.

You may then apply the following code using pandas:

Finally, you can visually represent the covariance matrix using the seaborn and matplotlib packages:

Geometric and Intuitive Explanation¶

Following is a geometric and intuitive explanation of the covariance matrix and the way it describes the shape of a data set. We will describe the geometric relationship of the covariance matrix with the use of linear transformations and eigendecomposition.

Before we get started, we shall take a quick look at the difference between covariance and variance. Variance measures the variation of a single random variable (like the height of a person in a population), whereas covariance is a measure of how much two random variables vary together (like the height of a person and the weight of a person in a population). The formula for variance is given by

$$\sigma_x^2=\frac{1}{n-1}\sum_{i=1}^n(x_i-\bar x)^2$$

where n is the number of samples (e.g. the number of people) and $\bar x$ is the mean of the random variable $\bf x$ (represented as a vector). The covariance $\sigma (x,y)$ of two random variables x and y is given by

$$\sigma(x,y)=\frac{1}{n-1}\sum_{i=1}^n(x_i-\bar x)(y_i-\bar y)$$

with n samples. The variance $\sigma_x^2$ of a random variable x can be also expressed as the covariance with itself by $\sigma(x,x)$.

Covariance Matrix¶

With the covariance we can calculate entries of the covariance matrix, which is a square matrix given by

$C_{i,j}=\sigma(x_i,x_j)$ where $C \in \mathbb{R}^{d \times d}$ and d describes the dimension or number of random variables of the data (e.g. the number of features like height, width, weight, …). Also the covariance matrix is symmetric since $\sigma(x_i,x_j) = \sigma(x_j,x_i)$. The diagonal entries of the covariance matrix are the variances and the other entries are the covariances. For this reason, the covariance matrix is sometimes called the variance-covariance matrix. The calculation for the covariance matrix can be also expressed as

$$C= \frac{1}{n-1}\sum_{i=1}^n(X_i-\bar X)(X_i-\bar X)^T$$

where our data set is expressed by the matrix $X \in \mathbb{R}^{n \times d}$. Following from this equation, the covariance matrix can be computed for a data set with zero mean with $C = \frac{XX^T}{n-1}$ by using the semi-definite matrix $XX^T$.

In this article, we will focus on the two-dimensional case, but it can be easily generalized to more dimensional data. Following from the previous equations the covariance matrix for two dimensions is given by

$$\begin{equation*} C = \begin{bmatrix} \sigma(x,x) & \sigma(x,y) \\ \sigma(y,x) & \sigma(y,y) \\ \end{bmatrix} \end{equation*}$$

We want to show how linear transformations affect the data set and in result the covariance matrix. First we will generate random points with mean values $\bar x,\bar y$ at the origin and unit variance $\sigma_x^2=\sigma_y^2=1$ which is also called white noise and has the identity matrix as the covariance matrix.

This case would mean that x and y are independent (or uncorrelated) and the covariance matrix C is

$$\begin{equation*} C = \begin{bmatrix} \sigma_x^2 & 0 \\ 0 & \sigma_y^2 \\ \end{bmatrix} \end{equation*}$$

We can check this by calculating the covariance matrix as follows

The above results approximatelly gives us our expected covariance matrix with variances $\sigma_x^2=\sigma_y^2=1$.

Linear Transformations of the Data Set¶

Now, we will see how linear transformations affect the data and the covariance matrix C. We will transform our data with the following scaling matrix.

$$\begin{equation*} S = \begin{bmatrix} s_x & 0 \\ 0 & s_y \\ \end{bmatrix} \end{equation*}$$

where the transformation simply scales the x and y components by multiplying them by $s_x$ and $s_y$ respectively. The covariance matrix C of the transformed data set will simply be

$$\begin{equation*} C = \begin{bmatrix} (s_x\sigma_x)^2 & 0 \\ 0 & (s_y\sigma_y)^2 \\ \end{bmatrix} \end{equation*}$$

which means that we can extract the scaling matrix from our covariance matrix by calculating $S=\sqrt C$ and the data is transformed by $Y=SX$. .

Now we will apply a linear transformation in the form of a transformation matrix T to the data set which will be composed of a two dimensional rotation matrix R and the scaling matrix S as follows

$$T = RS$$

where the rotation matrix R is given by

$$\begin{equation*} C = \begin{bmatrix} cos\theta & -sin\theta \\ sin\theta & cos\theta \\ \end{bmatrix} \end{equation*}$$

where $\theta$ is the angle of rotation. The transformed data is then calculated by Y = TX or Y = RSX.

Eigen Decomposition of the Covariance Matrix¶

Eigen Decomposition is one connection between a linear transformation and the covariance matrix. An eigenvector is a vector whose direction remains unchanged when a linear transformation is applied to it. It can be expressed as

$$Av = \lambda v$$

where v is an eigenvector of A and $\lambda$ is the corresponding eigenvalue. If we put all eigenvectors into the columns of a Matrix $\bf V$ and all eigenvalues as the entries of a diagonal matrix L we can write for our covariance matrix C with the following equation

$$CV = VL$$

where the covariance matrix can be represented as

$$C=VLV^{-1}$$

which can be also obtained by $\bf Singular \ Value \ Decomposition$. The eigenvectors are unit vectors representing the direction of the largest variance of the data, while the eigenvalues represent the magnitude of this variance in the corresponding directions. This means V represents a rotation matrix and $\sqrt L$ represents a scaling matrix. From this equation, we can represent the covariance matrix C as

$$C = RSSR^{-1}$$

where the rotation matrix R = V and the scaling matrix $S=\sqrt L$. Hence,

$$C = RSSR^{-1}=TT^{-1}$$

because $T^T=(RS)^T=S^TR^T=SR^{-1}$ due to the properties $R^{-1}=R^T$ since R is orthogonal and $S=S^T$ since S is a diagonal matrix. This enables us to calculate the covariance matrix from a linear transformation. In order to calculate the linear transformation of the covariance matrix, one must calculate the eigenvectors and eigenvectors from the covariance matrix C. This can be done by calculating

$$T = V\sqrt L$$

where V is the previous matrix where the columns are the eigenvectors of C and L is the previous diagonal matrix consisting of the corresponding eigenvalues. The transformation matrix can be also computed by the Cholesky decomposition with $Z = L^{-1}(X-\bar X)$ where L is the Cholesky factor of $C = LL^T$.

We can see the basis vectors of the transformation matrix by showing each eigenvector v multiplied by $\sigma = \sqrt \lambda$. By multiplying $\sigma$ with 3 we cover approximately 99.7% of the points according to the three sigma rule if we would draw an ellipse with the two basis vectors and count the points inside the ellipse.

We can now get from the covariance the transformation matrix T and we can use the inverse of T to remove correlation (whiten) the data.