Numeric から Symbolic へ

　処理が重いなら WorkingPrecision オプションを Infinity から有限に変更する手もありますが，それは Mathematica 固有の解決策なので，ここでは一般的な方法を採ります．すなわち

Numeric に推定して，Symbolic に証明する

という道です．
　具体的には

g[r_] := Reduce[
　Exists[{a, b},ForAll[{x, y},Implies[(x - a)^2 + (y - b)^2 <= r^2, f[x, y] ] ] ],
　Reals];

のように r が前回の論理式を満たすか否かを判定する関数を作り

A = 0; B = 1; M = (A + B)/2;

から始めて

While[True,
　If[g[M],
　　A = M; Print["T", N[M, 50] ];,
　　B = M; Print["F", N[M, 50] ]; ];
　M = (A + B)/2;]

のように区間の 2 等分を繰り返すと

f[x_, y_] := x^2 <= y <= x;

の場合

F 0.50000000000000000000000000000000000000000000000000

F 0.25000000000000000000000000000000000000000000000000

F 0.12500000000000000000000000000000000000000000000000

T 0.062500000000000000000000000000000000000000000000000

F 0.093750000000000000000000000000000000000000000000000

T 0.078125000000000000000000000000000000000000000000000

T 0.085937500000000000000000000000000000000000000000000

F 0.089843750000000000000000000000000000000000000000000

T 0.087890625000000000000000000000000000000000000000000

…

といった列（ T は g[r] が True，F は g[r] が False となる r ）が出力され，程なく

0.088388347648318440550105545263106129910604492211059

に至りますので，その正体（？）を Wolfram|Alpha に尋ねてみると
http://www.wolframalpha.com/input/?i=0.088388347648318440550105545263106129910604492211059
と答えてくれますので，改めて，Mathematica に

Reduce[Implies[0 <= r <= Sqrt[2]/16,
r >= 0 &&
Exists[{a, b},
ForAll[{x, y},
Implies[(x - a)^2 + (y - b)^2 <= r^2, f[x, y]]]]], Reals]
Reduce[Implies[
r >= 0 &&
Exists[{a, b},
ForAll[{x, y}, Implies[(x - a)^2 + (y - b)^2 <= r^2, f[x, y]]]],
0 <= r <= Sqrt[2]/16], Reals]

と入力すれば，どちらも

True

と出力され，中心については

Reduce[ForAll[{x, y},
　Implies[(x - a)^2 + (y - b)^2 <= 1/128, f[x, y] ] ], Reals]

と入力すれば

b == 5/16 && a == 7/16

と出力されます．