Eigenvalue spread of `random_psd_operator` is very skewed

[Newtech66]

The issue

The eigenvalue spread of a PSD matrix generated by random_psd_operator has a strange distribution that I didn't expect:

from toqito.rand import random_psd_operator
import matplotlib.pyplot as plt
import seaborn as sns

eigs_rpsd = []
for _ in range(10):
    eigs_rpsd.extend(eigvalsh(random_psd_operator(2 ** 6, is_real=True)).tolist())
plt.title('random_psd_operator(N)')
ax = sns.histplot(data=eigs_rpsd, kde=True)
ax.lines[0].set_color('crimson')
plt.xlabel('Eigenvalues')
plt.show()

Expected behavior

I'm not sure what I should expect really: when we say "random PSD operator", we have to define what random really means. Perhaps drawing from a probability measure over PSD operators, like scipy.stats.wishart? In my case, I thought the eigenvalues would be uniformly distributed.

Wishart has the following distribution for reference:

from scipy.stats import wishart
import matplotlib.pyplot as plt
import seaborn as sns

eigs_wishart = []
for _ in range(10):
    eigs_wishart.extend(eigvalsh(wishart.rvs(2 ** 8, np.eye(2 ** 8) / 2 ** 8)).tolist())
plt.title('Wishart(N, I / N)')
ax = sns.histplot(data=eigs_wishart, kde=True)
ax.lines[0].set_color('crimson')
plt.xlabel('Eigenvalues')
plt.show()

[vprusso]

Thanks for flagging, @Newtech66. Indeed, we likely want to parameterize the probability distribution we are using and replace the rand_mat = gen.random((dim, dim)) distribution with something more customizable for a better and more appropriate distribution. If you'd be inclined to put in a fix for that, that would be a great improvement! If not, that's okay too, and thanks in any case for flagging, @Newtech66 !