Can/should one generate Ginibre ensembles of random matrices using low-discrepancy normal variates--and if so. how?

Question

I've been generating (via Mathematica) series of $4 \times 4$ "random density matrices with respect to Bures (minimal monotone) measure" https://arxiv.org/abs/0909.5094 [eq. (24)] and testing certain (``entanglement''-related) statistics based on them [https://arxiv.org/abs/quant-ph/9911058, https://arxiv.org/abs/quant-ph/0308037]. (Such matrices are Hermitian, positive-definite with unit trace--that is, unit-trace Wishart.) Each matrix ($\rho_{\frac{1}{2}}$ in the notation of the first ref.) requires 64 random unit normal variates for its generation.

Should I be able to speed convergence by using low-discrepancy series of normal variates, and if so how might that be effectively accomplished? (any Mathematica code/suggestions would be appreciated).

Answer 1

It now seems rather clear (but subject to more detailed assessment) that the answer to the first part ("can/should one,...") is yes.

As to the "how" part of the question, I've made use of the very elegant algorithm presented by Martin Roberts in his Mathematica.stackexchange answer https://mathematica.stackexchange.com/questions/143457/how-can-one-generate-an-open-ended-sequence-of-low-discrepancy-points-in-3d

now using $d = 36$ and $d = 64$ (setting, as suggested, the supplementary parameter $\alpha_0=\frac{1}{2}$) in the two cases ("two-re[al]bits" and "two-qu[antume"]bits) of particular interest. I seem to be getting much stronger convergence than when I simply use (unstructured) random normal variates.

To convert from the points ($p$) the Roberts algorithm yields, that are uniformly distributed over $[0,1]^d$, to the unit normal variates required in the random matrix generation, I employ the Mathematica command InverseCDF[NormalDistribution[0, 1],p].

The Mathematica code I am employing in the $d=36$ (two-rebit) case is

sp2 = x /. Solve[x^(37) == x + 1, x][[1]]; 
G = Array[1, 36]; 
Do[G[[i]] = 1/2 + i95/sp2^i, {i, 1, 36}]; 
rB = 0; 
Do[hw = TimeUsed[];
P = InverseCDF[NormalDistribution[0, 1], N[FractionalPart[G]]]; 
Y1 = (Orthogonalize[ArrayReshape[Take[P, {1, 16}], {4, 4}]] + 
     IdentityMatrix[4]).ArrayReshape[Take[P, {17, 36}], {4, 5}];
z = Partition[Y1.Transpose[Y1], {2, 2}]; 
If[PositiveDefiniteMatrixQ[
    ArrayFlatten@{{z[[1, 1]], z[[2, 1]]}, {z[[1, 2]], z[[2, 2]]}}] == 
   True, rB = rB + 1];
If[Mod[i95, 500000] == 0, 
  Print[{TimeUsed[] - hw, i95, rB, rB/i95, N[rB/i95, 20]}]; 
{i95,rB} >> LowDiscrepancyBuresTwoRebitsSave], {i95, 1, 1000000000}]

The ratio of the variable rb to the variable i95 is the quantity of interest (that is, the separability probability with respect to Bures measure of the two-rebit states).