I was wondering how I can generate about 50 random 2x2 symmetrical positive definite matrices in PYTHON. Is there a block of code using scikit or numpy anyone can guide me to?
3 Answers
Strictly diagonally dominant matrices are positive definite. So you could generate a random A, compute AA= A'A and then increase the elements on the diagonal to make sure that
AA[i,i] > sum( abs(AA[i,j]), j != i ),
e.g., compute the sum on the right-hand side and then add 1 to it and assign the result to AA[i,i].
More simply, you could compute A'A + alpha*I for some alpha > 0 of your choice and where I is the
identity matrix. All eigenvalues of this matrix are >= alpha, which make it "safely" positive definite.
You can also use sklearn to do it easier and then make matrix values random. See: https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_spd_matrix.html
A real symmetric 2-by-2 matrix
/ a b \
\ b c /
is positive definite iff its trace and determinant are positive.
As the determinant is ac-b^2 this happens iff a and c are positive and -m < b < m where m is the geometric mean \sqrt ac.
You can therefore draw positive a and c and a factor f with -1<f<1 and then use
/ a f * sqrt(ac) \
\ f * sqrt(ac) c /
Note that this method can easily be vectorized.
"Random" means different things to different people but I'll take a stab.
As a previous answer noted, a real matrix X is positive semi-definite if (and only if) it can be factored as X=M'M, i.e. it is a Gram matrix.
One simple way to make X random, while also guaranteeing positive definiteness, is to use random normal entries for M. See also, the Wishart distribution.
Related