Skip to content

Instantly share code, notes, and snippets.

Show Gist options
  • Save hobson/ad451869dd11439cc8ff to your computer and use it in GitHub Desktop.
Save hobson/ad451869dd11439cc8ff to your computer and use it in GitHub Desktop.
r"""Faulty positive semidefinite matrix generator
Using a normally distributed random matrix it's pretty
rare that you'll get a semidefinite matrix from the product,
even for low dimensions. For 3 dimensions you have about a 0.001%% yield
>>> pos_semi_matrices = generate_positive_semidefinite_matrices(num_samples=1000)
6.7% were positive semidefinite and not positive definite
45.5% were positive definite
52.2% were positive semidefinite
47.8% were neither positive definite nor positive semidefinite
"""
import numpy as np
def generate_positive_semidefinite_matrices(num_dimensions=3, num_samples=1000):
np.random.seed(2)
positive_semidefinite_matrices = []
pos_semidef = np.zeros(num_samples).astype(bool)
pos_def = np.zeros(num_samples).astype(bool)
for i in range(num_samples):
B = np.random.randn(num_dimensions - 1, num_dimensions)
A = np.dot(B.T, B)
eigvals = np.linalg.eigvals(A)
# is it positive semidefinite?
pos_semidef[i] = np.all(eigvals >= 0)
pos_def[i] = np.all(eigvals > 0)
if pos_semidef[i]:
positive_semidefinite_matrices.append(A)
num_samples = float(num_samples)
print("{: 6.1%} were positive semidefinite and not positive definite".format(
np.sum(~pos_def & pos_semidef) / num_samples))
print("{: 6.1%} were positive definite".format(
np.sum(pos_def) / num_samples))
print("{: 6.1%} were positive semidefinite".format(
np.sum(pos_semidef) / num_samples))
print("{: 6.1%} were neither positive definite nor positive semidefinite".format(
np.sum(~pos_semidef & ~pos_def) / num_samples))
return positive_semidefinite_matrices, pos_semidef, pos_def
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment