Ask Your Question
0

Random positive definite matrix

asked 2017-11-16 15:58:26 +0200

Deepak Sarma gravatar image

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 flag offensive close merge delete

2 Answers

Sort by ยป oldest newest most voted
0

answered 2017-11-16 19:58:21 +0200

dan_fulea gravatar image

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]
edit flag offensive delete link more

Comments

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.

Deepak Sarma gravatar imageDeepak Sarma ( 2017-11-17 06:52:45 +0200 )edit

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

as in https://en.wikipedia.org/wiki/Positive-definite_matrix#Simultaneous_diagonalization

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.

dan_fulea gravatar imagedan_fulea ( 2017-11-17 11:01:17 +0200 )edit

Got it. Thank you

Deepak Sarma gravatar imageDeepak Sarma ( 2017-11-17 11:55:45 +0200 )edit
1

answered 2017-11-16 17:00:47 +0200

vdelecroix gravatar image

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

edit flag offensive delete link more

Comments

Thank you.

Deepak Sarma gravatar imageDeepak Sarma ( 2017-11-17 06:54:29 +0200 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

Stats

Asked: 2017-11-16 15:58:26 +0200

Seen: 1,513 times

Last updated: Nov 16 '17