programming - Can I use Compile to speed up InverseCDF?

I'm repeatedly issuing the command

P = InverseCDF[NormalDistribution[0, 1], T]

where T is either a 36- or 64-vector of real numbers (points in the unit hypercubes $[0,1]^{36}$ or $[0,1]^{64}$, quasirandomly distributed according to the "golden-ratio" generalization formula of Martin Roberts given in his answer to How can one generate an open-ended sequence of low-discrepancy points in 3D?). (The points--after conversion by the indicated command--are employed for the generation of "two-rebit" and "two-qubit" $4 \times 4$ density matrices, randomly generated with respect to the Bures [minimal monotone] measure, in accordance with formulas (24) and (28), setting $x =\frac{1}{2}$, of "Random Bures mixed states and the distribution of their purity" https://arxiv.org/abs/0909.5094 .)

An analysis of mine shows that this is by far the most time-consuming step in my program.

So, can I use Compile so that the command "knows" that the vector is composed of reals? I presume/hope that, if so, this would lead to a considerable speed-up.

Answer

Out of curiosity, I tried to write my own version of inverse CDF for the normal distribution. I employ a qualitative approximation of the inverse CDF as initial guess and apply Newton iterations with line search until convergence.

This is the code:

f[x_] := CDF[NormalDistribution[0, 1], x];
finv[y_] := InverseCDF[NormalDistribution[0, 1], y];

p = 1/200;
q = 2/5;
g[x_] = N@Piecewise[{
     {finv[$MachineEpsilon], 0 <= x <= $MachineEpsilon},
     {Simplify[Normal@Series[finv[x], {x, 0, 1}], 0 < x < 1], $MachineEpsilon < x < p},
     {Simplify[PadeApproximant[finv[x], {x, 1/2, {7, 8}}]], q < x < 1 - q},
     {Simplify[Normal@Series[finv[x], {x, 1, 1}], 0 < x < 1], 1 - p < x < 1},
     {finv[1 - $MachineEpsilon], 1 - $MachineEpsilon <= x <= 1}
     },
    Simplify[PadeApproximant[finv[x], {x, 1/2, {17, 18}}]]

    ];
(*g[y_]:=Log[Abs[(1-Sqrt[1-y]+Sqrt[y])/(1+Sqrt[1-y]-Sqrt[y])]];*)

cfinv = Block[{T, S, Sτ}, With[{
     S0 = N[g[T]],
     ϕ00 = N[(T - f[S])^2],
     direction = N[Simplify[(T - f[S])/f'[S]]],
     residual = N[(T - f[Sτ])^2],
     σ = 0.0001,
     γ = 0.5,

     TOL = 1. 10^-15
     },
    Compile[{{T, _Real}},
     Block[{S, Sτ, ϕ0, τ, u, ϕτ},
      S = S0;
      ϕ0 = ϕ00;
      While[Sqrt[ϕ0] > TOL,
       u = direction;
       τ = 1.;
       Sτ = S + τ u;

       ϕτ = residual;
       While[ϕτ > (1. - σ τ) ϕ0, 
        τ *= γ;
        Sτ = S + τ u;
        ϕτ = residual;
        ];
       ϕ0 = ϕτ;
       S = Sτ;
       ];
      S

      ],
     CompilationTarget -> "C",
     RuntimeAttributes -> {Listable},
     Parallelization -> True,
     RuntimeOptions -> "Speed"
     ]
    ]
   ];

And here an obligatory test for speed and accuracy (tested on a Haswell Quad Core CPU):

T = RandomReal[{0., 1}, {1000000}];
a = finv[T]; // RepeatedTiming // First
b = cfinv[T]; // RepeatedTiming // First
Max[Abs[a - b]]

0.416

0.0533

3.77653*10^-12

So this one is almost three eight times faster than the built-in method.

I also expect a speedup should one replace the function g by a better but also reasonably quick approximation for the initial guess. First I tried an InterpolatingFunction for the (suitably) transformed "true" CDF, but that turned out to be way too slow.

Of course Newton's method has its problems on the extreme tails of the distribution (close to 0 and close 1) where the CDF has derivative close to 0. Maybe a secant method would have been more appropriate?

Using expansions of the inverse CDF at $0$, $1/2$ and $1$, I was able to come up with a way better initial guess function g.