Changing iteration variable in Do loop

Question

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?

Answer 1

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.

Answer 2

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.

Answer 3

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)