Why does the following Solve take so long?

Question

Why does the following take Mathematica so long to reply to?

Solve[BetaRegularized[0.5, k, n - k - 1] == 0.035, k]

where $n = 1000$

By the way, In case there is no exact solution, how can I request Mathematica to find the nearest solution?

Answer 1

Have a look for FindRoot. Per default, it tries to use Newton's method. However, the derivatives behave rather badly, so a secant method might be a more robust choice. Brent's method (mixture of bisection and secant method) works very well, here:

n = 1000;
b = 0.035;
sol = FindRoot[
  BetaRegularized[0.5, x, n - x - 1] == b, {x, 0.4 n, 0.6 n}, 
  Method -> "Brent"]
F[x] - b /. sol

{x -> 528.115}

-3.84415*10^-15

Note that 0.4 n and 0.6 n are supposed to bracket the root of the equation.

Answer 2

Solve originally was principally a solver of polynomial equations $p(x)=0$ and some equations like $p(f(x))=0$ for suitable transcendental functions $f$. But beginning around V7, I think, efforts have been made to extend the reach of Solve so that it might live up to its name. Many transcendental equations can be solved if the roots can be isolated, and if in no other way, their roots can be represented by exact, symbolic Root objects.

It often helps to have exact coefficients.

Here Solve works, if the root is isolated and a sufficiently high, but finite, working precision is used:

Block[{n = 1000},
 Solve[BetaRegularized[1/2, k, 1000 - k - 1] == 7/200 && 
    2 n/5 <= k <= 3 n/5, k, Reals,
  WorkingPrecision -> 400]
 ]
(*  {{k -> 528.1147041429262349889...}}  *)

Alternatively, one can construct the corresponding Root object:

(sol = Block[{n = 1000},
    With[{f = BetaRegularized[1/2, k, n - k - 1] - 7/200 & /. k -> #},
     Root[{f,
       k /. FindRoot[f[k], {k, n/2, 2 n/5, 3 n/5}, 
         WorkingPrecision -> 400]}]
     ]]) // AbsoluteTiming
(*
  {20.7221, 
   Root[{-(7/200) + BetaRegularized[1/2, #1, 999 - #1] &, 
     528.11470414292623498897865816253955484889... + 0.*10^-403 I}]}
*)

Remarks on the numerics of the equation

We can see that at machine precision, the function is constant for n = 1000 until the middle of its domain 0 < k < 999:

Block[{n = 1000},
 With[{f = 
    BetaRegularized[0.5, k, n - k - 1] - 0.035 & /. k -> # // 
     Rationalize},
  InputForm /@ Table[f[N@k], {k, 0, 1000, 100}]
  ]]
(*
  {0.965`, 0.965`, 0.965`, 0.965`,
   0.9649999998691152`, 0.4523748657765868`, -0.03499999991364527`,
   -0.035`, -0.035`, -0.035`, -0.035`} 
*)

A locally constant function generally stumps solvers that use local methods. One should consider RootSearch or findAllRoots; or perhaps System`TRootsDump`GuessRoots (or System`TRootsDump`GuessRealRoots) from this related Q&A: Faster variant of Reduce for finding zeros of holomorphic functions in a region

Perhaps it is not surprising that the derivative is difficult to evaluate at machine precision:

Block[{n = 1000},
 With[{f = BetaRegularized[0.5, k, n - k - 1] - 0.035 & /. k -> # // Rationalize},
  Table[f'[N@k], {k, Join[{1}, 100 Range[9], {998}]}]
  ]]
(*
  {-4.46485*10^155, -7.43223*10^246, -8.32217*10^275, 1.77529*10^275, 
   4.69065*10^254, -2.14015*10^216, 3.1861*10^253, -1.13516*10^276, 
   2.6163*10^275, 2.96348*10^246, -4.46485*10^155}
*)

You need ~500 digits to get accurate values throughout the whole domain; only ~400 are needed for k between 400 and 600, which perhaps explains why WorkingPrecision -> 400 works above.

Block[{n = 1000},
 With[{f = BetaRegularized[0.5, k, n - k - 1] - 0.035 & /. k -> # // Rationalize},
  N@Table[
    f'[SetPrecision[k, 500]], {k, Join[{1}, 100 Range[9], {998}]}]
  ]]
(*
  {-2.79399*10^-300, -5.31457*10^-162, -7.31643*10^-86, 
   -6.31716*10^-38, -5.41614*10^-11, -0.0252418, -3.60875*10^-11, 
   -2.70387*10^-38, -1.82439*10^-86, -5.86766*10^-163, -2.79399*10^-300} 
*)

Answer 3

I think there is something strange about the derivative of BetaRegularized.

First let's plot your function.

Show[
 Plot[BetaRegularized[0.5, k, 1000 - k - 1], {k, 450, 550}, 
  PlotStyle -> Black],
 ListLinePlot[{{450, 0.035}, {550, 0.035}}, PlotStyle -> Red]
 ]

I defined a function

f[k_] := BetaRegularized[0.5, k, 1000 - k - 1] - 0.035

and attempted to locate k where f[k] was zero.

I tried some numerical methods (e.g, FindMinimum, FindRoot) and had no luck.

I computed the derivative symbolically and numerically

D[f[k], k] /. k -> 520
(* -2.34977*10^206 *)

the result is inconsistent with the plot.

I defined an approximate numerical derivative

df[k_] := (f[k + 10^-6] - f[k])/10^-6

df[520]
(* -0.0108823 *)

This made much more sense.

I then used FixedPoint and a simple Newton's method to get an approximation of k.

FixedPoint[# - f[#]/df[#] &, 530, 20]
(* 528.115 *)

Here is a plot on a magnified scale near the solution

Show[
 Plot[BetaRegularized[0.5, k, 1000 - k - 1], {k, 527, 529}, 
  PlotStyle -> Black],
 ListLinePlot[{{527, 0.035}, {529, 0.035}}, PlotStyle -> Red],
 Graphics[{
   Blue,
   Opacity[0.6],
   PointSize -> 0.025,
   Point[{{528.115, 0.035}}]
   }]
 ]

So my guess is that there is something strange going on with the derivative.

Note that I am completely unfamiliar with BetaRegularized so you may very well get a better answer.