Return to Answer

The following implementation is basically Newton's mwthod. It assumes that towards the end of the optimization, the relation between data points and closest reference points does not change anymore. In that regime, the objective function is quadratic and Newton converges in a single step. The Newton search direction is also a good choice outside of this regime, and so we can employ Armijo line search to globalize Newton's method. (Newton's method is notorious for wild behavior outsize a basin of attraction. One of the reasons (when the Hessian is positive-definite) is that the size of the step is typically too large.)

The following implementation is basically Newton's method. It assumes that towards the end of the optimization, the relation between data points and closest reference points does not change anymore. In that regime, the objective function is quadratic and Newton converges in a single step. The Newton search direction is also a good choice outside of this regime, and so we can employ Armijo line search to globalize Newton's method. (Newton's method is notorious for wild behavior outsize a basin of attraction. One of the reasons (when the Hessian is positive-definite) is that the size of the step is typically too large.)

ClearAll[findShift]; findShift[data_, reference_, {a_, b_, n_Integer}] := Module[{shift, d, h, Y, nf}, (* n+1 = the number of shifting operations per dimension.*) (* a = leftmost shift in each direction*) (* b rightmost shift in each direction*) (* The step size.*) h = (b - a)/n; (*Build only one nearest function for all shifting operations.*) nf = Nearest[data -> "Distance"]; d = Dimensions[data][[-1]]; (*a vector to store the shifts*) shift = ConstantArray[0., d]; Y = reference; Do[ Y[[All, direction]] -= (a - h); shift[[direction]] = MinimalBy[ (*You can parallize this loop if you want.*) Table[ (*The trick is to shift the reference,  not the data! *) (*This allows us to reuse the Nearest \ function nf.*) Y[[All, direction]] -= h; {s, Total[nf[Y]^2, 2]} , {s, a, b, h} ], Last ][[1, 1]]; (*Apply optimal shift in current direction.*) Y[[All, direction]] += b - shift[[direction]]; , {direction, 1, d}]; (*Finally, apply the shifts to the data once.*) data + ConstantArray[shift, Length[data]] ] 

ClearAll[findShift]; findShift[data_, reference_, {a_, b_, n_Integer}] := Module[{shift, d, h, Y, nf}, (* n+1 = the number of shifting operations per dimension.*) (* a = leftmost shift in each direction*) (* b rightmost shift in each direction*) (* The step size.*) h = (b - a)/n; (*Build only one nearest function for all shifting operations.*) nf = Nearest[data -> "Distance"]; d = Dimensions[data][[-1]]; (*a vector to store the shifts*) shift = ConstantArray[0., d]; Y = reference; Do[ Y[[All, direction]] -= (a - h); shift[[direction]] = MinimalBy[ (*You can parallize this loop if you want.*) Table[ (*The trick is to shift the reference, not the data! *) (*This allows us to reuse the Nearest function nf.*) Y[[All, direction]] -= h; {s, Total[nf[Y]^2, 2]} , {s, a, b, h} ], Last ][[1, 1]]; (*Apply optimal shift in current direction.*) Y[[All, direction]] += b - shift[[direction]]; , {direction, 1, d}]; (*Finally, apply the shifts to the data once.*) data + ConstantArray[shift, Length[data]] ] 

Okay, after I have seen that Newton's method does not work well here, here my reformulation of OP's algorithm. This implementation finds the same shifts, but it is far more efficient for large datasets.

ClearAll[findShift]; findShift[data_, reference_, {a_, b_, n_Integer}] := Module[{shift, d, h, Y, nf}, (* n+1 = the number of shifting operations per dimension.*) (* a = leftmost shift in each direction*) (* b rightmost shift in each direction*) (* The step size.*) h = (b - a)/n; (*Build only one nearest function for all shifting operations.*) nf = Nearest[data -> "Distance"]; d = Dimensions[data][[-1]]; (*a vector to store the shifts*) shift = ConstantArray[0., d]; Y = reference; Do[ Y[[All, direction]] -= (a - h); shift[[direction]] = MinimalBy[ (*You can parallize this loop if you want.*) Table[ (*The trick is to shift the reference, not the data! *) (*This allows us to reuse the Nearest \ function nf.*) Y[[All, direction]] -= h; {s, Total[nf[Y]^2, 2]} , {s, a, b, h} ], Last ][[1, 1]]; (*Apply optimal shift in current direction.*) Y[[All, direction]] += b - shift[[direction]]; , {direction, 1, d}]; (*Finally, apply the shifts to the data once.*) data + ConstantArray[shift, Length[data]] ] 

And here is a a usage example.

reference = RandomReal[{-1, 1}, {100, 3}]; data1 = RandomReal[{-1, 1}, {100000, 3}]; data1shifted2 = findShift[data1, reference, {-1., 1., 20}]; // AbsoluteTiming // First 

0.027269

OP's code requires 5.87231 seconds for the same task on my machine. So this is an improvement of factor 200.

Newton's method is not very effective for this type of problem since there are many local minima. I keep the following only as reference.