Covariant-matrix of a matrix¶

Types of variance and covariance¶

  • Population variance: Used when you have data for the entire population. $$\sigma^2 = \frac{1}{N} \sum_{i=1}^{N} (x_i - \mu)^2$$
  • Sample variance: Used when you only have a sample of the population. To avoid bias, divide by $n-1$ instead of $n$. $$s^2 = \frac{1}{n-1} \sum_{i=1}^{n} (x_i - \bar{x})^2$$ The same thing goes for covariance.
  • Population covariance $$\text{Cov}(X, Y) = \frac{1}{N} \sum_{i=1}^{N} (x_i - \mu_X)(y_i - \mu_Y)$$
  • Sample covariance $$\text{Cov}(X, Y) = \frac{1}{n-1} \sum_{i=1}^{n} (x_i - \bar{x})(y_i - \bar{y})$$ When estimating variance or covariance from a sample, using $n$ underestimates the true variability (bias). So instead of dividing by $n$, we divide by $n-1$. This is called Bessel's correction.


There are 2 popular ways of finding the covariant-matrix of a matrix:

1. The direct method¶

If we want to find the covariant matrix w.r.t the column: $$ \begin{bmatrix} Var(r_1) & Cov(r_1, r_2) & \cdots & Cov(r_1, r_n) \\ Cov(r_2, r_1) & Var(r_2) & \cdots & Cov(r_2, r_n) \\ \vdots & \vdots & \ddots & \vdots \\ Cov(r_n, r_1) & Cov(r_n, r_2) & \cdots & Var(r_n) \end{bmatrix} $$ If we want to find the covariant matrix w.r.t the row, we just transpose the input matrix and do the same operation.

2. The matrix formulation method¶

  • If we want to find the covariant matrix w.r.t the column, calculate the mean of each column and subtract that value from each of the columns. Let the resultant matrix be $X$. $$\text{Covariance matrix with Bessel's correction}=\frac{1}{N-1}X \ X^T$$
  • If we want to find the covariant matrix w.r.t the row, calculate the mean of each row and subtract that value from each of the rows. Let the resultant matrix be $X$. $$\text{Covariance matrix with Bessel's correction}=\frac{1}{n-1}X^T \ X$$
In [1]:
matrix = [
 [2, 5, -7, 8],
 [-5, 9, -1, -3],
 [-1, 2, -2, 3],
 [-8, -6, 7, 10]
]

Function to find covariance between 2 arrays¶

In [2]:
from typing import List

def covariance(arr1: List[float], arr2: List[float], ddof: int = 1) -> float:
    """
    Calculate the covariance between two arrays using only standard Python.
    
    Parameters:
    arr1, arr2 : list
        Input arrays of the same length
    
    Returns:
    float
        Covariance between arr1 and arr2
    """
    if len(arr1) != len(arr2):
        raise ValueError("Arrays must have the same length")
    if len(arr1) == 0:
        raise ValueError("Arrays cannot be empty")
    
    n = len(arr1)
    mean1 = sum(arr1) / n
    mean2 = sum(arr2) / n
    
    # Calculate covariance: sum((x_i - mean1) * (y_i - mean2)) / (n - 1)
    return sum((x - mean1) * (y - mean2) for x, y in zip(arr1, arr2)) / (n - ddof)

Function to multiply 2 matrices¶

In [3]:
def matrix_multiply(A: List[List[float]], B: List[List[float]]) -> List[List[float]]:
    """
    Multiply two matrices A and B using standard Python.
    
    Parameters:
    A : list of lists
        First matrix (m x n)
    B : list of lists
        Second matrix (n x p)
    
    Returns:
    list of lists
        Resulting matrix (m x p)
    """
    # Check if matrices can be multiplied
    if not A or not B or len(A[0]) != len(B):
        raise ValueError("Invalid matrix dimensions for multiplication")
    
    rows_A = len(A)
    cols_A = len(A[0])
    cols_B = len(B[0])
    
    # Initialize result matrix with zeros
    result = [[0 for _ in range(cols_B)] for _ in range(rows_A)]
    
    # Perform matrix multiplication
    for i in range(rows_A):
        for j in range(cols_B):
            for k in range(cols_A):
                result[i][j] += A[i][k] * B[k][j]
    
    return result

Function to multiply a scalar to a matrix¶

In [4]:
def scalar_matrix_multiply(scalar, matrix):
    """
    Multiply a scalar by each element of a matrix.
    
    Parameters:
    scalar : number (int or float)
        The scalar value to multiply
    matrix : list of lists
        Input matrix (m x n)
    
    Returns:
    list of lists
        Resulting matrix (m x n) with each element multiplied by scalar
    """
    if not matrix or not matrix[0]:
        raise ValueError("Matrix cannot be empty")
    
    rows = len(matrix)
    cols = len(matrix[0])
    
    # Initialize result matrix
    result = [[0 for _ in range(cols)] for _ in range(rows)]
    
    # Perform scalar multiplication
    for i in range(rows):
        for j in range(cols):
            result[i][j] = round(scalar * matrix[i][j], 2)
    
    return result

Function to transpose a matrix¶

In [5]:
def matrix_transpose(matrix: List[List[float]]) -> List[List[float]]:
    """
    Calculate the transpose of a matrix.
    
    Parameters:
    matrix : list of lists
        Input matrix (m x n)
    
    Returns:
    list of lists
        Transposed matrix (n x m)
    """
    if not matrix or not matrix[0]:
        raise ValueError("Matrix cannot be empty")
    
    rows = len(matrix)
    cols = len(matrix[0])
    
    # Initialize result matrix with dimensions (cols x rows)
    result = [[0 for _ in range(rows)] for _ in range(cols)]
    
    # Compute transpose
    for i in range(rows):
        for j in range(cols):
            result[j][i] = matrix[i][j]
    
    return result

Direct method to find covariance¶

In [6]:
def mean(arr: List[int]):
    return sum(arr) / len(arr)

def variance(arr: List[int], ddof: int = 1):
    n = len(arr)
    mean_ = sum(arr) / n
    return sum([(mean_ - i) ** 2 for i in arr]) / (n - ddof)

def covariance_matrix_direct(matrix_: List[List[float]], ddof_ = 1, axis = 0) -> List[List[float]]:
    """ Parameters
    ddof_: 0 if bessel's correction is to be applied, otherwise 1
    axis: 0 if we want the correlation between rows, otherwise (for column) 1
    """
    n = len(matrix)
    m = len(matrix_[0])

    if not n == m:
        raise ValueError("Number of rwos and column of matrix must be equal")

    if not (axis == 0 or axis == 1):
        raise ValueError("Axis must be either 0 or 1")
    
    matrix_t = (matrix_transpose(matrix_) if axis == 1 else matrix_)
    covariance_matrix = []
    
    for i in range(n):
        arr = []
        for j in range(n):
            if i == j:
               arr.append(round(variance(matrix_t[i], ddof = ddof_), 2))
            else:
               arr.append(round(covariance(matrix_t[i], matrix_t[j], ddof = ddof_), 2))
        covariance_matrix.append(arr[:])
    
    return covariance_matrix

Matrix formulation method to find covariance¶

In [7]:
def covariance_matrix_formulation(matrix_: List[List[float]], ddof_: int = 1, axis = 0) -> List[List[float]]:
    """ Parameters
    ddof_: 0 if bessel's correction is to be applied, otherwise 1
    axis: 0 if we want the correlation between rows, otherwise (for column) 1
    """
    n = len(matrix_)
    m = len(matrix_[0])

    if not n == m:
        raise ValueError("Number of rwos and column of matrix must be equal")

    if not (axis == 0 or axis == 1):
        raise ValueError("Axis must be either 0 or 1")
    
    mean_vector = (
        [mean([matrix_[i][j] for i in range(n)]) for j in range(n)] if axis == 1
        else [mean([matrix_[i][j] for j in range(n)]) for i in range(n)]
    )
    
    matrix_t = [row[:] for row in matrix]  # deep-copy of the matrix
    
    for row in range(n):
        for column in range(n):
            matrix_t[row][column] -= mean_vector[(column if axis == 1 else row)] 
    
    if axis == 1:
        covariance_matrix = matrix_multiply(matrix_transpose(matrix_t), matrix_t)
    else:
        covariance_matrix = matrix_multiply(matrix_t, matrix_transpose(matrix_t))
    return scalar_matrix_multiply(1.0 / (n - ddof_), covariance_matrix)

Driver code for Direct Method with Bessel correction¶

W.r.t column¶

In [8]:
covariance_matrix_direct(matrix, axis = 1)
Out[8]:
[[19.33, 13.67, -24.0, 0.67],
 [13.67, 40.33, -27.5, -28.67],
 [-24.0, -27.5, 33.58, 8.17],
 [0.67, -28.67, 8.17, 33.67]]

W.r.t row¶

In [9]:
covariance_matrix_direct(matrix, axis = 0)
Out[9]:
[[42.0, 6.0, 14.0, -7.0],
 [6.0, 38.67, 5.33, -17.0],
 [14.0, 5.33, 5.67, 3.5],
 [-7.0, -17.0, 3.5, 82.25]]

Driver code for Direct Method without Bessel correction¶

W.r.t column¶

In [10]:
covariance_matrix_direct(matrix, ddof_ = 0, axis = 1)
Out[10]:
[[14.5, 10.25, -18.0, 0.5],
 [10.25, 30.25, -20.62, -21.5],
 [-18.0, -20.62, 25.19, 6.12],
 [0.5, -21.5, 6.12, 25.25]]

W.r.t row¶

In [11]:
covariance_matrix_direct(matrix, ddof_ = 0, axis = 0)
Out[11]:
[[31.5, 4.5, 10.5, -5.25],
 [4.5, 29.0, 4.0, -12.75],
 [10.5, 4.0, 4.25, 2.62],
 [-5.25, -12.75, 2.62, 61.69]]

Driver code for Matrix Formulatiom Method with Bessel correction¶

W.r.t column¶

In [12]:
covariance_matrix_formulation(matrix, axis = 1)
Out[12]:
[[19.33, 13.67, -24.0, 0.67],
 [13.67, 40.33, -27.5, -28.67],
 [-24.0, -27.5, 33.58, 8.17],
 [0.67, -28.67, 8.17, 33.67]]

W.r.t row¶

In [13]:
covariance_matrix_formulation(matrix, axis = 0)
Out[13]:
[[42.0, 6.0, 14.0, -7.0],
 [6.0, 38.67, 5.33, -17.0],
 [14.0, 5.33, 5.67, 3.5],
 [-7.0, -17.0, 3.5, 82.25]]

Driver code for Matrix Formulatiom Method without Bessel correction¶

W.r.t column¶

In [14]:
covariance_matrix_formulation(matrix, ddof_ = 0, axis = 1)
Out[14]:
[[14.5, 10.25, -18.0, 0.5],
 [10.25, 30.25, -20.62, -21.5],
 [-18.0, -20.62, 25.19, 6.12],
 [0.5, -21.5, 6.12, 25.25]]

W.r.t row¶

In [15]:
covariance_matrix_formulation(matrix, ddof_ = 0, axis = 0)
Out[15]:
[[31.5, 4.5, 10.5, -5.25],
 [4.5, 29.0, 4.0, -12.75],
 [10.5, 4.0, 4.25, 2.62],
 [-5.25, -12.75, 2.62, 61.69]]

Driver code for covariance using numpy with bessel correction¶

W.r.t column¶

In [16]:
import numpy as np

covariance_matrix = np.cov(matrix, rowvar = False)
covariance_matrix
Out[16]:
array([[ 19.33333333,  13.66666667, -24.        ,   0.66666667],
       [ 13.66666667,  40.33333333, -27.5       , -28.66666667],
       [-24.        , -27.5       ,  33.58333333,   8.16666667],
       [  0.66666667, -28.66666667,   8.16666667,  33.66666667]])

W.r.t row¶

In [17]:
covariance_matrix = np.cov(matrix, rowvar = True)
covariance_matrix
Out[17]:
array([[ 42.        ,   6.        ,  14.        ,  -7.        ],
       [  6.        ,  38.66666667,   5.33333333, -17.        ],
       [ 14.        ,   5.33333333,   5.66666667,   3.5       ],
       [ -7.        , -17.        ,   3.5       ,  82.25      ]])

Driver code for covariance using numpy without bessel correction¶

W.r.t column¶

In [18]:
covariance_matrix = np.cov(matrix, ddof = 0, rowvar = False)
covariance_matrix
Out[18]:
array([[ 14.5   ,  10.25  , -18.    ,   0.5   ],
       [ 10.25  ,  30.25  , -20.625 , -21.5   ],
       [-18.    , -20.625 ,  25.1875,   6.125 ],
       [  0.5   , -21.5   ,   6.125 ,  25.25  ]])

W.r.t row¶

In [19]:
covariance_matrix = np.cov(matrix, ddof = 0, rowvar = True)
covariance_matrix
Out[19]:
array([[ 31.5   ,   4.5   ,  10.5   ,  -5.25  ],
       [  4.5   ,  29.    ,   4.    , -12.75  ],
       [ 10.5   ,   4.    ,   4.25  ,   2.625 ],
       [ -5.25  , -12.75  ,   2.625 ,  61.6875]])