[SOLVED] Random positive semi-definite matrix with given eigenvalues and eigenvectors

Random positive semi-definite matrix with given eigenvalues and eigenvectors

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).

Solution

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.