Return to Revisions

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.)

This improves two parts of OP's implementation:

1. It updates all coordinate directions at once.

2. It builds only one NearestFunction for the reference data and applies is multiple times. This is an advantage because everytime we call Nearest, a certain space partition tree hast to be build. So while the actual data queries are very fast, there is a certain overhead for building the seach index. By calling Nearest only once, we have to pay the overhead only once, of course.

Here is the actual code:

(*We have to generate only on NearestFunction for the reference data; then we can apply it repeatedly.*) nf = Nearest[reference]; X = data1; (*Some hard-coded parameters for the algorithm.*) \[Sigma] = 0.5; \[Gamma] = 0.5; TOL = 1. 10^-8; maxbacktrackings = 100; maxiterations = 100; Y = nf[X, 1][[All, 1]]; (*Objective: sum of squares*) F = 1/2 Total[(Y - X)^2, 2]; (*This is the derivative of the objective.*) DF = Total[X - Y]; (*This is the Newton search direction.*) u = -DF/Length[X]; residualsquared = u.u; iter = 0; While[residualsquared > TOL^2 && iter < maxiterations, iter++; (*Armijo backtracking line search "with quadratic fit".*) (*Eventually, step size t = 1 should work because we use Newton's method.*) t = 1.; Xt = X + t ConstantArray[DF, Length[X]]; Yt = nf[Xt, 1][[All, 1]]; Ft = 1/2 Total[(Yt - Xt)^2, 2]; slope = DF.u; backtrackings = 0; While[Ft > F - t \[Sigma] slope && backtrackings < maxbacktrackings, backtrackings++; tnew = -\[Sigma] t^2 slope/(Ft - F - t slope); t = Max[\[Gamma] t, tnew]; Xt = X + t ConstantArray[u, Length[X]]; Yt = nf[Xt, 1][[All, 1]]; Ft = 1/2 Total[(Yt - Xt)^2, 2] ]; (*Now we have found a feasible step size t. Updating data*) X = Xt; Y = Yt; F = Ft; (*This is the derivative of the objective.*) DF = Total[X - Y]; (*This is the Newton search direction.*) u = -DF/Length[X]; residualsquared = u.u; ]; 

The result (the optimally translated data) is now stored in X.

On the example data, this converges with two iterations.

Edit

The algorithm minimizes the sum of squared distances from each reference point to its clostests point in the dataset. If you want to minimize the sum of squared distances of each data point from its closest point in the references, the following should do:

Z = reference; X = data1; nf = Nearest[X]; (*Some hard-coded parameters for the algorithm.*) \[Sigma] = 0.5; \[Gamma] = 0.5; TOL = 1. 10^-8; maxbacktrackings = 20; maxiterations = 100; Y = nf[Z, 1][[All, 1]]; (*Objective:sum of squares*) F = 1/2 Total[(Z - Y)^2, 2]; Print["initial fit" -> F]; (*This is the derivative of the objective.*) DF = Total[Y - Z]; (*This is the Newton search direction.*) u = -DF/Length[Z]; residualsquared = u.u; iter = 0; While[ residualsquared > TOL^2 && iter < maxiterations, iter++; (*Armijo backtracking line search "with quadratic \ fit".*)(*Eventually,step size t= 1 should work because we use Newton's method.*)t = 1.; t = 1.; Xt = X + t ConstantArray[u, Length[X]]; nft = Nearest[Xt]; Yt = nft[Z, 1][[All, 1]]; Ft = 1/2 Total[(Z - Yt)^2, 2]; slope = DF.u; backtrackings = 0; While[Ft > F + t \[Sigma] slope && backtrackings < maxbacktrackings, backtrackings++; tnew = -\[Sigma] t^2 slope/(Ft - F - t slope); t = Max[\[Gamma] t, tnew]; Xt = X + t ConstantArray[u, Length[X]]; nft = Nearest[Xt]; Yt = nft[Z, 1][[All, 1]]; Ft = 1/2 Total[(Z - Yt)^2, 2]; ]; (*Now we have found a feasible step size t.Updating data*) X = Xt; Y = Yt; F = Ft; nf = nft; (*This is the derivative of the objective.*) DF = Total[Y - Z]; (*This is the Newton search direction.*) u = -DF/Length[Z]; residualsquared = u.u; ]; Print["final fit" -> F]; Print["residual" -> Sqrt[residualsquared]];