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


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}]; 
    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).

| cite | improve this answer | |
  • $\begingroup$ To make your answer easier to read, I have suggested an edit to your answer that formats your mathematica code as a code block. $\endgroup$ – Martin Roberts Jul 23 '18 at 2:49

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.