2
votes

Is there a way to generate a random positive semi-definite matrix with given eigenvalues and eigenvectors in Python?

I looked at this, but they do not allow to specify eigenvalues for matrix construction.

Context: I want to generate random multivariate Gaussians with controlled ellipticity and because the major/minor axes of the distribution have the length proportional to eigenvalues I want my covariance matrix to have them. Definiton could be found here (page 81).

1

1 Answers

4
votes

When you don't have the eigenvectors but only want some eigenvalues, you can list your desired eigenvalues and use a orthonormal matrix to jumble them up. Since congruence transformations don't change the inertia of a matrix (well up to numerical precision) you can use the Q matrix of the QR decomposition of a random matrix (or any other way to generate an orthonormal matrix).

import numpy as np
import scipy.linalg as la
des = [1, 0, 3, 4, -2, 0, 0]
n = len(des)
s = np.diag(des)
q, _ = la.qr(np.random.rand(n, n))
semidef = q.T @ s @ q
np.linalg.eigvalsh(semidef)

gives

array([-2.00000000e+00, -2.99629568e-16, -5.50063275e-18,  2.16993906e-16,
        1.00000000e+00,  3.00000000e+00,  4.00000000e+00])

When you actually have also the eigenvectors then you can simply construct the original matrix anyways which is the definition of eigenvalue decomposition.