How to create internally optimized expression for computing with high WorkingPrecision?

On the base of findings in this thread the procedure is as follows.

Assuming that residualVect is the residual vector and initGuess is the initial guess for FindMinimum, let us create an internally optimized residual vector via Experimental`CreateNumericalFunction:

resV = Experimental`CreateNumericalFunction[initGuess[[;; , 1]], 
   Sqrt[2]*residualVect, {Length[residualVect]}, WorkingPrecision -> 20];

Note that we have multiplied the true residual vector by Sqrt[2] in order to get the result consistent with the ordinary FindMinimum syntax where the first argument is equal to Total[residualVect^2].

Now FindMinimum can be used as follows:

FindMinimum[Null, initGuess, Method -> {"LevenbergMarquardt", 
   "Residual" :> resV[initGuess[[;; , 1]]], 
   "Jacobian" :> resV["Jacobian"@initGuess[[;; , 1]]]},
                             WorkingPrecision -> 20]

This approach gives more than one order of magnitude speedup in my application.

Thanks to the recent post by ilian who uncovered and explained the undocumented Optimization`FindFit`ObjectiveFunction, I have found Optimization`FindFit`ResidualFunction which seems to be ideally suited to the original goals of the question. With the latter we don't even need to construct the residual vector explicitly, we just specify the model function and the data and get the optimized NumericalFunction automatically! A little example from the linked thread:

data = BlockRandom[SeedRandom[12345];
   Table[{x, Exp[-2.3 x/(11 + .4 x + x^2)] + RandomReal[{-.5, .5}]}, {x, 
     RandomReal[{1, 15}, 20]}]];

residual = Optimization`FindFit`ResidualFunction[data, Exp[a x/(b + c x)], {a, b, c}, x]

Experimental`NumericalFunction[{a,b,c},Block[{x={2.69745,5.61891,11.9585,7.02236,4.1302,7.48274,11.3322,10.9067,12.0728,2.47999,6.32617,5.88977,8.63809,11.2304,11.0019,8.65308,4.26835,9.26918,3.04309,9.8615},Optimization`FindFit`y$588={0.972294,0.344088,0.767092,0.460704,0.947346,1.21116,0.439919,0.675085,1.17798,0.961116,0.720896,1.18664,0.605936,0.811138,0.861414,0.670646,0.770535,0.365288,0.649543,0.551157}},E^((a x)/(b+c x))-Optimization`FindFit`y$588],-NumericalFunctionData-]