Rotate covariance matrix

Question

I am generating 3D gaussian point clouds. I'm using the scipy.stats.multivariate.normal() function, which takes a mean value and a covariance matrix as arguments. It can then provide random samples using the rvs() method.

Next I want to perform a rotation of the cloud in 3D, but rather than rotate each point I would like to rotate the random variable parameters, and then regenerate the point cloud.

I'm really struggling to figure this out. After an rotation, the axes of variance will no longer align with the coordinate system. So I guess what what I want is to express variance along three arbitrary orthogonal axes.

Thank you for any help.

Final edit: Thank you, I got what I needed. Below is an example

cov = np.array([
   [ 3.89801357,  0.38668784,  1.47657614],
   [ 0.38668784,  0.87396495,  1.43575688],
   [ 1.47657614,  1.43575688, 15.09192414]])

rotation_matrix = np.array([
   [ 2.22044605e-16,  0.00000000e+00,  1.00000000e+00],
   [ 0.00000000e+00,  1.00000000e+00,  0.00000000e+00],
   [-1.00000000e+00,  0.00000000e+00,  2.22044605e-16]]) # 90 degrees around y axis

new_cov = rotation_matrix @ cov @ rotation_matrix.T # based on Warren and Paul's comments




rv = scipy.stats.multivariate_normal(mean=mean,cov=new_cov)

If you get an error

ValueError: the input matrix must be positive semidefinite

This page I found useful

*"So I guess what what I want is to express variance along three arbitrary orthogonal axes."* So you already know the axes, and you know the variance that you want along each axis? Is that your input? — Warren Weckesser, Oct 22 '18 at 05:04
Yes exactly. I have a rotation matrix, so I could calculate new x,y,z unit vectors, I just don't know how to modify the covariance matrix properly to plot it in the original coordinate system — skrei, Oct 22 '18 at 05:07
I'm not totally sure if this is the best site for this question. It's a good question, it just might be more mathematical than computational. If other people think similarly, [stats.se] might be a good place to ask it instead, but wait to see what others say first. — David Z, Oct 22 '18 at 05:16
If Σ is the unrotated covariance matrix, and R is the rotation matrix, the rotated covariance should be RΣR⁻¹. (David Z is right, this question would be better on the stats or math stackexchange sites.) — Warren Weckesser, Oct 22 '18 at 05:27
Btw, you probably know but just in case: the inverse R^-1 can be cheaply computed as R^T. — Paul Panzer, Oct 22 '18 at 05:58
Thanks everyone for the responses. The reason I posted here is because I don't have much of a formal math background and am more familiar with numpy syntax rather than symbolic syntax. Does RΣR⁻¹ refer to matrix multiplication of R and Σ, follwed by a multiplication of the inverse of R? — skrei, Oct 22 '18 at 17:52
Yes, those products are matrix multiplication. You can use `@` in Python 3. And as Paul Panzer pointed out, the inverse of the rotation matrix is just the transpose. So you can write: `new_cov = rotation_matrix @ cov @ rotation_matrix.T` — Warren Weckesser, Oct 22 '18 at 19:41

score 4 · Accepted Answer · answered Oct 22 '18 at 23:23

4

I edited the question with the answer, but again it is

new_cov = rotation_matrix @ cov @ rotation_matrix.T

answered Oct 22 '18 at 23:23

skrei

81
1
6

Rotate covariance matrix

1 Answers1