One possibility for the situation where one gets a local min is as follows. Intersect the region with a sphere of radius slightly larger than the given distance, centered at the query point. Now one can repeat with a better chance of getting the global min. Moreover one can get an approximation by meshing the subregion, extracting vertices, and using ordinary Nearest to find the closest vertex to the query point. I show this below.

Here is the original setup and the local min found.

pt = {-0.5, 4};
region = ImplicitRegion[y - x^2 == 0, {x, y}];
pt1 = RegionNearest[region, pt]
dist = EuclideanDistance[pt1, pt]

(* Out[64]= {1.83399714916, 3.36354556229}

Out[65]= 2.41921825051 *)

If we use a subregion comprised of points at most slightly further away, then we get a better minimal distance/closer point. We could of course go in the direction of being more exclusive, by taking a radius that excludes the point we found, and for some purposes that might be better. But I'll leave that for others to consider.

region2 = 
  ImplicitRegion[
   y - x^2 == 0 && ({x, y} - pt).({x, y} - pt) <= 1.01*dist^2, {x, y}];
pt2 = RegionNearest[region2, pt]
dist2 = EuclideanDistance[pt2, pt]

(* Out[68]= {-1.90556930822, 3.63119438843}

Out[69]= 1.45314922129 *)

Here is the meshing approach mentioned above.

ptlist = MeshCoordinates[DiscretizeRegion[region2]];
pt3 = First[Nearest[ptlist, pt]]
dist3 = EuclideanDistance[pt3, pt]

(* Out[71]= {-1.92080363722, 3.68948661277}

Out[72]= 1.4543388667 *)

One will note that it does give a good approximation to the globally closest point. This method might be slow, due to meshing algorithm internals, or too coarse, due to quality of mesh, so caveat emptor.

After writing the above I realized some of this is built in. One can just discretize the original region and do RegionNearest on that.

pt4 = RegionNearest[DiscretizeRegion[region], pt]
dist4 = EuclideanDistance[pt4, pt]

(* Out[74]= {-1.90493346777, 3.63315045973}

Out[75]= 1.45203878531 *)

Again, there might be speed issues and or artifacts from quality of the discretization such as a result that is too close (if discretization goes outside the boundary of the actual region) or too far (if discretization is too coarse to capture sufficient detail).

I do not know if this is regarded as a bug or feature of using approximate values/methods. I have passed the question on to people who might know.

Global optimization is quite difficult, so I'd say that the documentation statement " NMinimize may sometimes only find a local minimum" is overly optimistic.

What I've done is to find a subregion containing the global min, which might or might not be the local one. I also keep the local min and in fact do it this way to make sure I do not have an empty region. I suppose I could allow it to be empty, and let region code tell me that, and deduce from there that the local min was also either a global one, or at least quite close.

It may have made it linear but that hardly means there is no (local) min there. We are minimizing distance from a point, and that doesn't care so much whether the "edges" are locally linear or actual curves.

What matters is that with a discretization we now have, internally, a simplicial complex in the sense of computational geometry. This has considerable structure, and internal code for RegonNearest is better equipped to handle that than it is to deal with a general implicit region.

I think it would make more sense for Minimize and the related functions to rationalize inputs that are not infinite precision, rather than switching to the numerical version.

Hi Philip,

Indeed, this seems a little strange but expected in some way. Either giving the starting point or the curve with some limited precision will give the same kinda error. Probably it will transfer to different algorithm once either the start or the curve has some precision.

For this curve it is quite straight-forward but you can imagine much harder cases where the curve makes loops and crossing et cetera. I think they have a general algorithm that looks for a minimum.

This does not change of course the fact that it is indeed wrong, please submit an bug-report with just:

{xcorrect, ycorrect} = RegionNearest[ImplicitRegion[y - x^2 == 0, {x, y}], {-1/2, 4}] // N
{xwrong, ywrong} = RegionNearest[ImplicitRegion[y - x^2 == 0, {x, y}], {-0.5, 4}] 

to https://www.wolfram.com/support/contact/email/ ("product feedback").

NB:

{xcorrect, ycorrect} = RegionNearest[ImplicitRegion[y - x^2 == 0, {x, y}], {1/2, 4}] // N
{xwrong, ywrong} = RegionNearest[ImplicitRegion[y - x^2 == 0, {x, y}], {0.5, 4}] 

does work correctly because the numerical minimalization is starting on the 'right' side.