2
$\begingroup$

I want to generate random positive hermitian matrices, I start with such a code

rho11R = Table[RandomReal[], {i, 1, 4}];
rho10R = Table[RandomComplex[], {i, 1, 4}];
Do[If[rho11R[[i]] - rho11R[[i]]^2 - 
     rho10R[[i]]*Conjugate[rho10R[[i]]] >= 0, 
   Nothing, {rho10R[[i]] = RandomComplex[], --i}], {i, 1, 
   4}];
rhoR = Table[{{rho11R[[i]], rho10R[[i]]}, {Conjugate[rho10R[[i]]], 
    1 - rho11R[[i]]}}, {i, 1, 4}]

This does not work, as the matrices I get are in general not positive. The problem is my If function, which does not seem to change the iteration variable i to i-1, when the condition for positivity is not fulfilled (so that I can get new random variable for rho10R[[I]] and again check the condition for positivity). How can I reach this?

Improve this question
$\endgroup$
2
  • 1
    $\begingroup$ Have you tried a While[] loop instead? $\endgroup$ – Somos Jul 19 '19 at 15:53
  • $\begingroup$ I assume that by positive you mean positive semi-definite, right? $\endgroup$ – Roman Jul 19 '19 at 17:04
4
$\begingroup$

This is one the cases where For can be helpful (generally it's just more complicated.) With For, you can manipulate the iteration variable, which you cannot do with Do.

Example:

For[
 i = 1,
 i <= 4,
 i++,
 If[
  i == 2,
  Print[i++],
  Print[i]
  ]
 ]

This prints 1, 2, 4.

Another option is to use While, as Somos wrote in a comment.

Improve this answer
$\endgroup$
4
$\begingroup$

A much easier way to construct random positive semi-definite Hermitian matrices is to start with Gaussian random matrices and Hermitian-square them:

randommatrix[n_Integer?Positive] := 
  RandomVariate[NormalDistribution[], {n, n, 2}].{1, I}

randomHermitian[n_Integer?Positive] :=
  (# + ConjugateTranspose[#])/2 &[ConjugateTranspose[#].# &[randommatrix[n]]]

In this way you don't need to reject anything.

Test:

randomHermitian[10] // Eigenvalues
(*    {71.4553, 53.6575, 46.3275, 31.8263, 21.4754,
       12.9687, 7.36107, 4.40568, 1.23665, 0.199904}    *)

There is of course the question of the measure (distribution) from which you pull the random matrices; you'd have to be more specific in your question to address this point.

Thanks to @mikado for pointing out that Hermitian symmetrization is needed to avoid generating matrices that are almost-but-not-quite Hermitian because of numerical precision limits.

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ While the matrices you generate are Hermitian by construction, they may not be exactly, because of rounding errors. It is good practice to force to Hermitian, by averaging with the conjugate transpose. $\endgroup$ – mikado Jul 19 '19 at 22:14
0
$\begingroup$

The following is a more "functional" way of doing acceptance/rejection 

mtx := RandomComplex[1+I, {2, 2}]

pmtx := Module[{m = mtx, m1}, m1 = Chop[m + ConjugateTranspose[m]]; 
  If[PositiveDefiniteMatrixQ[m1], m1, pmtx]]

You can then generate a list of positive definite matrices of any length e.g.

Table[pmtx, 20]

(Note: I'm not addressing the question of whether this gives an appropriate distribution of random matrices)

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.