diff options
Diffstat (limited to 'schofield.wl')
-rw-r--r-- | schofield.wl | 41 |
1 files changed, 26 insertions, 15 deletions
diff --git a/schofield.wl b/schofield.wl index eb0525c..908ca94 100644 --- a/schofield.wl +++ b/schofield.wl @@ -208,40 +208,48 @@ newSol[eqs_, oldSol_, newVars_, δ_:0, γ_:0, opts___] := FindRoot[ opts ] -levenbergMarquardtStep[M_, λ_, g_] := LinearSolve[M + λ DiagonalMatrix[Diagonal[M]], g] +levenbergMarquardtStep[M_, λ_, g_] := Quiet[LinearSolve[M + λ DiagonalMatrix[Diagonal[M]], g], LinearSolve::luc] -levenburgMarquardt[r_, β0_, λ0_ : 1, ν_ : 2] := +levenburgMarquardt[r_, β0_, λ0_ : 1, ν_ : 2, ε_ : 10^-15] := Module[ { n = Length[β0], β = β0, λ = λ0, newβ = β0, - x, rf, Jf, oldJ, oldr, newr, M, g, δ, oldC, newC + x, rf, Jf, oldJ, oldr, newr, M, g, δ, oldC, newC, + h = 0.1, fvv, α, a, v }, PrintTemporary["Compiling the Jacobian..."]; - rf = Compile[{{x, _Real, 1}}, + rf = Quiet[Compile[{{x, _Real, 1}}, Evaluate[r /. Thread[Rule[First /@ β, Part[x, #] & /@ Range[n]]]] - ]; - Jf = Compile[{{x, _Real, 1}}, - Evaluate[D[r, {First /@ β}] /. 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]]]] + ], Part::partd]; PrintTemporary["Beginning the algorithm."]; oldr = rf[β[[All, 2]]]; oldJ = Jf[β[[All, 2]]]; oldC = Re[Total[oldr^2]]; + newC = oldC; g = Re[Transpose[oldJ] . oldr]; M = Re[Transpose[oldJ] . oldJ]; - δ = levenbergMarquardtStep[M, λ / ν, g]; - PrintTemporary[Dynamic[oldC]] - While[Norm[δ] > 10^-15, - newβ[[All, 2]] -= δ; + 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]] + While[Norm[δ] > ε, + newβ[[All, 2]] += δ; newr = rf[newβ[[All, 2]]]; newC = Re[Total[newr^2]]; While[newC > oldC, - δ = levenbergMarquardtStep[M, λ, g]; + 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; newβ = β; - newβ[[All, 2]] -= δ; + newβ[[All, 2]] += δ; newr = rf[newβ[[All, 2]]]; newC = Re[Total[newr^2]]; λ *= ν; @@ -253,7 +261,10 @@ levenburgMarquardt[r_, β0_, λ0_ : 1, ν_ : 2] := oldJ = Jf[β[[All, 2]]]; g = Re[Transpose[oldJ] . oldr]; M = Re[Transpose[oldJ] . oldJ]; - δ = levenbergMarquardtStep[M, λ / ν, g]; + 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; ]; {newC, β} ] |