# Random positive definite matrix

How can I get a random positive definite (or positive semi definite) matrix of given order? I know that there is a inbuilt function random_matrix with an additional feature algorithm to set with it. But there we can get some special matrices like 'echelon_form', 'orthogonal' ', 'echelonizable', 'diagonalizable'.... but positive definite command is not in built there. Some other important matrix classes are not there. So how to obtain random matrix of those classes?

edit retag close merge delete

Sort by ยป oldest newest most voted

Alternatively, built in method random_matrix can be used, while explicitly declaring the (ring, size, and) random eigenvalues, and their dimensions. Here, the sample code uses a poisson distribution with parameter lambda=7 of the given size. (Feel free to fill in in an other way the diagonal.)

from scipy.stats import poisson
# my way to get some random eigenvalues involves the above...
# use an other way to generate diag, if this is too artificial

N = 10
# we will produce an NxN positive definite random matrix with integer eigenvalues

diag = list( poisson.rvs( 7, size=N ) )
diag = [ 1+d for d in diag ]    # all entries are now > 0
evs  = list( set( diag ) )
dims = [ diag.count( ev ) for ev in evs ]

print "diag =", diag
print "evs  =", evs
print "dims =", dims

X = random_matrix( ZZ
, N
, algorithm  ='diagonalizable'
, eigenvalues=evs
, dimensions =dims )
print X


This gives this time:

diag = [9, 6, 11, 12, 3, 3, 10, 8, 10, 5]
evs  = [3, 5, 6, 8, 9, 10, 11, 12]
dims = [2, 1, 1, 1, 1, 2, 1, 1]
[  15    8    4    2   -2   -4    0   14   66  -24]
[  24   54   80  -24 -128  108   28    0 1012  384]
[  44   64   96  -20 -132  100   28   28 1152  348]
[ -28  -44  -64   25   98  -77  -21  -14 -834 -272]
[ -17  -28  -42   11   76  -52  -14   -7 -549 -186]
[   2   16   36  -14  -62   70   14  -14  452  202]
[  -2  -16  -36   14   62  -58   -4   14 -448 -206]
[  -8    8   32  -16  -60   62   14  -25  368  248]
[  -6  -12  -20    6   32  -27   -7    0 -247  -96]
[  -5   -4   -2   -1    1    2    0   -7  -33   17]

more

Thank you. One more thing, what to do if I don't want my matrix to be 'diagonalizable'. I tried to remove the algorithm, but there were some errors.

( 2017-11-17 06:52:45 +0200 )edit

If $A$ is a positive definite (square, symmetric, real valued) $N\times N$ matrix,

then it induces a Hilbert space structure on $\mathbb R^N$, seen as space of column vectors, i.e. matrices $N\times 1$, by setting: $$\langle x,y\rangle _A = x' Ay\ .$$ A Hilbert basis exists, and can be used to diagonalize. In fact, any symmetric matrix can be diagonalized.

( 2017-11-17 11:01:17 +0200 )edit

Got it. Thank you

( 2017-11-17 11:55:45 +0200 )edit

You could start from a random positive diagonal and conjugate

sage: A = diagonal_matrix(RDF, [RDF.random_element(1,5) for _ in range(5)])
sage: B = random_matrix(RDF, 5)
sage: C = B.transpose() * A * B
sage: C.is_positive_definite()
True


(The above is not a serious algorithm for integer matrices since all eigenvalues would be integral...)

more

Thank you.

( 2017-11-17 06:54:29 +0200 )edit