# Variance of product of a matrix and vector using Python

I am trying to calculate the variance of the product of matrix A and vector b. As it was mathematically discussed on another post https://math.stackexchange.com/questions/2365166/what-is-the-variance-of-a-constant-matrix-times-a-random-vector, I have tried to do the calculation in python.

To recap, we have $$$$Y = \bf{Ab}$$$$ where $$$$\bf{b} \sim N(\mu_b, \bf{I}\sigma_b^2)$$$$

Therefore $$$$Var(Y) = \bf{A}\sigma_b^2\bf{A}^T$$$$ I have tried this code with

import numpy as np
import statistics

A = np.arange(12).reshape(4, 3)   # some 4 by 3 array
b = np.array([2,4,5])              # some 3 by 1 vector

Y = np.dot(A,b)
statistics.variance(Y)
>> 1815

Vb = np.dot(np.identity(b.shape[0]), statistics.variance(b))
np.dot(A.dot(Vb), np.transpose(A))
>> array([[ 10.,  28.,  46.,  64.],
[ 28., 100., 172., 244.],
[ 46., 172., 298., 424.],
[ 64., 244., 424., 604.]])



However, the statistic.variance(Y) outputs a different result compared to np.dot(A.dot(Vb), np.transpose(A)).

I wonder how should I do it properly in python?

To test formulas like these numerically, you need to compute sample statistics, like this:

    import numpy as np
A = np.arange(12).reshape(4, 3)  # Note: A is a fixed matrix
b_avg = np.array([2, 4, 5])  # b is random, with a mean vector
b_var = 1e0 * np.eye(b_avg.size)  # and a (co)variance matrix

# generate a random sample of b vectors
n_sample = 10**4
b_sample = np.random.multivariate_normal(
mean=b_avg, cov=b_var, size=n_sample).T

# compute sample moments (each column is a sample)
b_avg_smp = b_sample.mean(axis=1)
b_var_smp = np.cov(b_sample)

# compare sample moments to population moments
print('b vector: statistics')
print('Averages (population vs. sample):',
b_avg.tolist(), b_avg_smp.round(1).tolist(), sep='\n')
print('Variances (population vs. sample):',
b_var.tolist(), b_var_smp.round(1).tolist(), sep='\n')


This will give something like this:

b vector: statistics
Averages (population vs. sample):
[2, 4, 5]
[2.0, 4.0, 5.0]
Variances (population vs. sample):
[[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]
[[1.0, -0.0, 0.0], [-0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]


i.e. the sample statistics are consistent with the specified $$b$$ distribution.

Then we can transform the $$b$$ sample into a $$y$$ sample, and do the analogous comparison, using your analytically computed $$y$$ moments, i.e.

    # transform b moments into y moments (@ is matrix multiply)
y_avg = A @ b_avg
y_var = A @ b_var @ A.T

# transform b sample into y sample, and compute sample moments
y_sample = A @ b_sample
y_avg_smp = y_sample.mean(axis=1)
y_var_smp = np.cov(y_sample)

# compare population vs. sample moments
print('y vector: statistics')
print('Averages (population vs. sample):',
y_avg.tolist(), y_avg_smp.round(1).tolist(), sep='\n')
print('Variances (population vs. sample):',
y_var.tolist(), y_var_smp.round(1).tolist(), sep='\n')


This will give something like this:

y vector: statistics
Averages (population vs. sample):
[14, 47, 80, 113]
[13.9, 46.8, 79.7, 112.7]
Variances (population vs. sample):
[[5.0, 14.0, 23.0, 32.0], [14.0, 50.0, 86.0, 122.0], [23.0, 86.0, 149.0, 212.0], [32.0, 122.0, 212.0, 302.0]]
[[5.0, 14.0, 23.1, 32.1], [14.0, 50.1, 86.1, 122.1], [23.1, 86.1, 149.0, 212.0], [32.1, 122.1, 212.0, 302.0]]


i.e. the sample $$y$$ statistics are consistent with the analytically computed $$y$$ distribution.

(Alternatively, the uncertainties package will automatically do linearized error propagation for Gaussian random variables.)

• Thank you very much, @GeoMatt22, I got the idea – user2842390 Jan 9 at 2:29

Your are mixing some concepts. It's hard to explain but i will try.

statistics.variance(b) gives you the variance of n samples in variable b in 1 dimension. The variance matrix for b is a 1x1 matrix.

Vb is the variance matrix for 3 dimensional sample with several samples and zero covarinace between the dimensions. But it is not the variance matrix for b.

The definition of A for b is wrong. It works because np.dot uses automatically b.T, but conceptually is wrong.

I hope this help to understand it.

• Thanks. Would you suggest any simulations as I showed in the code above to demonstrate how to calculate the variance of Y as in the scenario with Y=Ab where b ~ $N(\mu,\bf{I\sigma_b^2})$ ? Thank you! – user2842390 Jan 8 at 16:38