Improments to Levenberg-Marquardt.

1 files changed, 32 insertions, 20 deletions

diff --git a/schofield.wl b/schofield.wl
index 908ca94..77dce96 100644
--- a/schofield.wl
+++ b/schofield.wl
@@ -63,10 +63,10 @@ ddη[h_][f_] := D[f, θ] / D[t[θ] / h[θ]^(1 / Δ[2]), θ]
 invDerivativeList[n_][f_][x_] := Module[
     {xp, dfs, fp, Pns},
   dfs = Rest[NestList[D[#, xp] &, f[xp], n]] /. xp -> x;
-  Pns = FoldList[{Pm, m} |-> 
+  Pns = FoldList[{Pm, m} \[Function]
         fp'[xp] D[Pm, xp] - (2 m - 1) fp''[xp] Pm, 1, Range[n - 1]] /. 
     Derivative[m_][fp][xp] :> dfs[[m]];
-  MapIndexed[{Pn, i} |-> Pn/dfs[[1]]^(2 i[[1]] - 1), Pns]
+  MapIndexed[{Pn, i} \[Function] Pn/dfs[[1]]^(2 i[[1]] - 1), Pns]
   ]
 
 dFdξLowList[n_, h_][m_] := Module[
@@ -99,7 +99,7 @@ dFdη[n_, h_][m_][tt_] := Module[{ff, hh}, Nest[ddη[hh], h[θ]^(-2 / Δ[]) (ff[
 dFdηList[n_, h_][m_][tt_] := Module[{ff, hh}, NestList[ddη[hh], h[θ]^(-2 / Δ[2]) (ff[θ] - t[θ]^2 Log[hh[θ]^2] / (8 π Δ[2])), m] /. θ -> tt /. Map[Derivative[#][ff][tt] -> Derivative[#][RF[n]][tt] &, Range[0, m]] /. hh -> h]
 
 ruleB[g_] := B -> (2 * OverBar[s] / π) * (- g'[θc] / t[θc]^Δ[2])
-ruleθ0[g_] := Around[0.18930, 0.00005] - Simplify[g[I θ0]/(-t[I θ0])^Δ[2]/I]
+ruleθ0[g_] := Simplify[g[I θ0]/(-t[I θ0])^Δ[2]/I] - Around[0.18930, 0.00005]
 ruleAL[g_] := AL -> Exp[Δ[2] t[θc]^(Δ[2] - 1) t'[θc] / (2 OverBar[s] / π g'[θc]) - t[θc]^Δ[2] g''[θc] / (4 OverBar[s] / π g'[θc]^2)] t[θc]^(1/8) OverBar[s] / (2 π) * g'[θc]
 ruleAH[g_] := AH / ((g[I θ0]/ I)^(2 / Δ[2]) * (-η[g]'[I θ0] / (2 θ0 I))^(5/6)) + Around[1.37, 0.02]
 
@@ -188,8 +188,8 @@ eq[n_, g_][m_, p_, q_] := Flatten[Join[{ruleθ0[g], ruleAH[g], g'[0] θc - 1}, e
 eqAround[n_, g_][m_, p_, q_] := Flatten[Join[{ruleθ0[g], ruleAH[g]}, eqLow[n, g][#] & /@ Range[0, m],eqMid[RF[n], g][#] & /@ Range[0, p], eqHigh[n, g] /@ Range[2, q, 2]]] //. rules[g]
 
 formResiduals[data_, functions_, δ_:10^(-15)] := If[Head[#1]===Around,
-  (#1["Value"] - #2) / Max[#1["Uncertainty"], δ],
-  (#1 - #2) / δ] & @@@ Thread[{data, functions}]
+  (#2-#1["Value"]) / Max[#1["Uncertainty"], δ],
+  (#2-#1) / δ] & @@@ Thread[{data, functions}]
 
 resLow[n_, g_, δ_][m_] := formResiduals[Gls[[;;m+1]], dFdξLowList[n, g][m], δ]
 
@@ -208,7 +208,16 @@ newSol[eqs_, oldSol_, newVars_, δ_:0, γ_:0, opts___] := FindRoot[
   opts
 ]
 
-levenbergMarquardtStep[M_, λ_, g_] := Quiet[LinearSolve[M + λ DiagonalMatrix[Diagonal[M]], g], LinearSolve::luc]
+levenbergMarquardtStep[M_, λ_, g_] := Quiet[
+  - LinearSolve[M + λ DiagonalMatrix[Diagonal[M]], g]
+  , LinearSolve::luc
+]
+
+levenbergMarquardtAcceleration[M_, λ_, v_, J_, rold_, rf_, h_, β_] := Quiet[
+  - Re[LinearSolve[M + λ DiagonalMatrix[Diagonal[M]], Transpose[J].(2/h((rf[β[[All,2]] + h v] - rold) / h - J . v))] / 2]
+  , LinearSolve::luc
+]
+
 
 levenburgMarquardt[r_, β0_, λ0_ : 1, ν_ : 2, ε_ : 10^-15] := 
  Module[
@@ -217,37 +226,37 @@ levenburgMarquardt[r_, β0_, λ0_ : 1, ν_ : 2, ε_ : 10^-15] :=
     β = β0,
     λ = λ0,
     newβ = β0,
-    x, rf, Jf, oldJ, oldr, newr, M, g, δ, oldC, newC,
-    h = 0.1, fvv, α, a, v
+    x, rf, Jf, oldJ, oldr, newr, M, g, δ, oldC, newC, δlast,
+    h = 0.1, fvv, α = 0.75, a, v, U, nSteps = 0
    },
   PrintTemporary["Compiling the Jacobian..."];
   rf = Quiet[Compile[{{x, _Real, 1}}, 
     Evaluate[r /. Thread[Rule[First /@ β, Part[x, #] & /@ Range[n]]]]
   ], Part::partd];
   Jf = Quiet[Compile[{{x, _Real, 1}}, 
-    Evaluate[-D[r, {First /@ β}] /. Thread[Rule[First /@ β, Part[x, #] & /@ Range[n]]]]
+    Evaluate[D[r, {First /@ β}] /. Thread[Rule[First /@ β, Part[x, #] & /@ Range[n]]]]
   ], Part::partd];
   PrintTemporary["Beginning the algorithm."];
   oldr = rf[β[[All, 2]]];
   oldJ = Jf[β[[All, 2]]];
+  U = First[SingularValueDecomposition[oldJ]];
   oldC = Re[Total[oldr^2]];
   newC = oldC;
   g = Re[Transpose[oldJ] . oldr];
   M = Re[Transpose[oldJ] . oldJ];
   v = levenbergMarquardtStep[M, λ / ν, g];
-  a = LinearSolve[oldJ, -2/h((rf[β[[All,2]] + h v] - oldr) / h - oldJ . v)];
-  If[2 Norm[a] / Norm[v] > α, a /= (2 / Norm[v])];
-  δ = v + a / 2;
-  PrintTemporary["Current cost value: ", Dynamic[oldC]]
+  a = levenbergMarquardtAcceleration[M, λ / ν, v, oldJ, oldr, rf, h, β];
+  δ = v + a;
+  δlast = δ;
+  PrintTemporary["Current cost value: ", Dynamic[oldC]];
   While[Norm[δ] > ε,
     newβ[[All, 2]] += δ;
     newr = rf[newβ[[All, 2]]];
     newC = Re[Total[newr^2]];
-    While[newC > oldC,
+    While[((1 - δ . δlast / Norm[δ] / Norm[δlast])newC > oldC) || (2 Norm[a] / Norm[v] > α),
       v = levenbergMarquardtStep[M, λ, g];
-      a = LinearSolve[oldJ, -2/h((rf[β[[All,2]] + h v] - oldr) / h - oldJ . v)];
-      If[2 Norm[a] / Norm[v] > α, a /= (2 / Norm[v])];
-      δ = v + a / 2;
+      a = levenbergMarquardtAcceleration[M, λ, v, oldJ, oldr, rf, h, β];
+      δ = v + a;
       newβ = β;
       newβ[[All, 2]] += δ;
       newr = rf[newβ[[All, 2]]];
@@ -259,13 +268,16 @@ levenburgMarquardt[r_, β0_, λ0_ : 1, ν_ : 2, ε_ : 10^-15] :=
     oldr = newr;
     β = newβ;
     oldJ = Jf[β[[All, 2]]];
+    U = First[SingularValueDecomposition[oldJ]];
     g = Re[Transpose[oldJ] . oldr];
     M = Re[Transpose[oldJ] . oldJ];
     v = levenbergMarquardtStep[M, λ / ν, g];
-    a = LinearSolve[oldJ, -2/h((rf[β[[All,2]] + h v] - oldr) / h - oldJ . v)];
-    If[2 Norm[a] / Norm[v] > α, a /= (2 / Norm[v])];
-    δ = v + a / 2;
+    a = levenbergMarquardtAcceleration[M, λ / ν, v, oldJ, oldr, rf, h, β];
+    δlast = δ;
+    δ = v + a;
+    nSteps += 1;
   ];
+  Print[nSteps];
   {newC, β}
  ]