@@ -258,14 +258,14 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::TrustRegion,
258
258
initial_trust_radius = convert (trustType, alg. initial_trust_radius)
259
259
if iszero (initial_trust_radius)
260
260
initial_trust_radius = convert (trustType, max_trust_radius / 11 )
261
- end
261
+ end
262
262
end
263
263
step_threshold = convert (trustType, alg. step_threshold)
264
264
shrink_threshold = convert (trustType, alg. shrink_threshold)
265
265
expand_threshold = convert (trustType, alg. expand_threshold)
266
266
shrink_factor = convert (trustType, alg. shrink_factor)
267
267
expand_factor = convert (trustType, alg. expand_factor)
268
-
268
+
269
269
# Parameters for the Schemes
270
270
p1 = convert (floatType, 0.0 )
271
271
p2 = convert (floatType, 0.0 )
@@ -322,7 +322,8 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::TrustRegion,
322
322
jac_cache, false , maxiters, internalnorm, ReturnCode. Default, abstol, prob,
323
323
radius_update_scheme, initial_trust_radius, max_trust_radius, step_threshold,
324
324
shrink_threshold, expand_threshold, shrink_factor, expand_factor, loss, loss_new,
325
- H, g, shrink_counter, du, u_tmp, u_gauss_newton, u_cauchy, fu_new, make_new_J, r, p1, p2, p3, p4, ϵ,
325
+ H, g, shrink_counter, du, u_tmp, u_gauss_newton, u_cauchy, fu_new, make_new_J, r,
326
+ p1, p2, p3, p4, ϵ,
326
327
NLStats (1 , 0 , 0 , 0 , 0 ))
327
328
end
328
329
@@ -338,9 +339,9 @@ function perform_step!(cache::TrustRegionCache{true})
338
339
339
340
# do not use A = cache.H, b = _vec(cache.g) since it is equivalent
340
341
# to A = cache.J, b = _vec(fu) as long as the Jacobian is non-singular
341
- linres = dolinsolve (alg. precs, linsolve, A = J, b = _vec (fu),
342
- linu = _vec (u_gauss_newton),
343
- p = p, reltol = cache. abstol)
342
+ linres = dolinsolve (alg. precs, linsolve, A = J, b = _vec (fu),
343
+ linu = _vec (u_gauss_newton),
344
+ p = p, reltol = cache. abstol)
344
345
cache. linsolve = linres. cache
345
346
@. cache. u_gauss_newton = - 1 * u_gauss_newton
346
347
end
@@ -374,7 +375,7 @@ function perform_step!(cache::TrustRegionCache{false})
374
375
375
376
# Compute the potentially new u
376
377
cache. u_tmp = u + cache. du
377
-
378
+
378
379
cache. fu_new = f (cache. u_tmp, p)
379
380
trust_region_step! (cache)
380
381
cache. stats. nf += 1
@@ -406,7 +407,7 @@ function trust_region_step!(cache::TrustRegionCache)
406
407
407
408
# Compute the ratio of the actual reduction to the predicted reduction.
408
409
cache. r = - (loss - cache. loss_new) / (dot (du, g) + dot (du, H, du) / 2 )
409
- @unpack r = cache
410
+ @unpack r = cache
410
411
411
412
if radius_update_scheme === RadiusUpdateSchemes. Simple
412
413
# Update the trust region radius.
@@ -446,19 +447,19 @@ function trust_region_step!(cache::TrustRegionCache)
446
447
end
447
448
448
449
# trust region update
449
- if r < 1 // 10 # cache.shrink_threshold
450
- cache. trust_r *= 1 // 2 # cache.shrink_factor
451
- elseif r >= 9 // 10 # cache.expand_threshold
450
+ if r < 1 // 10 # cache.shrink_threshold
451
+ cache. trust_r *= 1 // 2 # cache.shrink_factor
452
+ elseif r >= 9 // 10 # cache.expand_threshold
452
453
cache. trust_r = 2 * norm (cache. du) # cache.expand_factor * norm(cache.du)
453
- elseif r >= 1 // 2 # cache.p1
454
- cache. trust_r = max (cache. trust_r, 2 * norm (cache. du)) # cache.expand_factor * norm(cache.du))
454
+ elseif r >= 1 // 2 # cache.p1
455
+ cache. trust_r = max (cache. trust_r, 2 * norm (cache. du)) # cache.expand_factor * norm(cache.du))
455
456
end
456
457
457
458
# convergence test
458
459
if iszero (cache. fu) || cache. internalnorm (cache. fu) < cache. abstol
459
460
cache. force_stop = true
460
461
end
461
-
462
+
462
463
elseif radius_update_scheme === RadiusUpdateSchemes. NW
463
464
# accept/reject decision
464
465
if r > cache. step_threshold # accept
@@ -471,9 +472,9 @@ function trust_region_step!(cache::TrustRegionCache)
471
472
472
473
if r < 1 // 4
473
474
cache. trust_r = (1 // 4 ) * norm (cache. du)
474
- elseif (r > (3 // 4 )) && abs (norm (cache. du) - cache. trust_r)/ cache. trust_r < 1e-6
475
- cache. trust_r = min (2 * cache. trust_r, cache. max_trust_r)
476
- end
475
+ elseif (r > (3 // 4 )) && abs (norm (cache. du) - cache. trust_r) / cache. trust_r < 1e-6
476
+ cache. trust_r = min (2 * cache. trust_r, cache. max_trust_r)
477
+ end
477
478
478
479
elseif radius_update_scheme === RadiusUpdateSchemes. Hei
479
480
if r > cache. step_threshold
@@ -579,25 +580,24 @@ function dogleg!(cache::TrustRegionCache{true})
579
580
# Take intersection of steepest descent direction and trust region if Cauchy point lies outside of trust region
580
581
l_grad = norm (cache. g) # length of the gradient
581
582
d_cauchy = l_grad^ 3 / dot (cache. g, cache. H, cache. g) # distance of the cauchy point from the current iterate
582
- if d_cauchy >= trust_r
583
- @. cache. du = - (trust_r/ l_grad) * cache. g # step to the end of the trust region
583
+ if d_cauchy >= trust_r
584
+ @. cache. du = - (trust_r / l_grad) * cache. g # step to the end of the trust region
584
585
return
585
586
end
586
587
587
588
# Take the intersection of dogled with trust region if Cauchy point lies inside the trust region
588
- @. u_cauchy = - (d_cauchy/ l_grad) * cache. g # compute Cauchy point
589
+ @. u_cauchy = - (d_cauchy / l_grad) * cache. g # compute Cauchy point
589
590
@. u_tmp = u_gauss_newton - u_cauchy # calf of the dogleg -- use u_tmp to avoid allocation
590
-
591
+
591
592
a = dot (u_tmp, u_tmp)
592
- b = 2 * dot (u_cauchy, u_tmp)
593
+ b = 2 * dot (u_cauchy, u_tmp)
593
594
c = d_cauchy^ 2 - trust_r^ 2
594
- aux = max (b^ 2 - 4 * a * c, 0.0 ) # technically guaranteed to be non-negative but hedging against floating point issues
595
- τ = (- b + sqrt (aux)) / (2 * a) # stepsize along dogleg to trust region boundary
595
+ aux = max (b^ 2 - 4 * a * c, 0.0 ) # technically guaranteed to be non-negative but hedging against floating point issues
596
+ τ = (- b + sqrt (aux)) / (2 * a) # stepsize along dogleg to trust region boundary
596
597
597
598
@. cache. du = u_cauchy + τ * u_tmp
598
599
end
599
600
600
-
601
601
function dogleg! (cache:: TrustRegionCache{false} )
602
602
@unpack u_tmp, u_gauss_newton, u_cauchy, trust_r = cache
603
603
@@ -611,18 +611,18 @@ function dogleg!(cache::TrustRegionCache{false})
611
611
l_grad = norm (cache. g)
612
612
d_cauchy = l_grad^ 3 / dot (cache. g, cache. H, cache. g) # distance of the cauchy point from the current iterate
613
613
if d_cauchy > trust_r # cauchy point lies outside of trust region
614
- cache. du = - (trust_r/ l_grad) * cache. g # step to the end of the trust region
614
+ cache. du = - (trust_r / l_grad) * cache. g # step to the end of the trust region
615
615
return
616
616
end
617
-
617
+
618
618
# Take the intersection of dogled with trust region if Cauchy point lies inside the trust region
619
- u_cauchy = - (d_cauchy/ l_grad) * cache. g # compute Cauchy point
619
+ u_cauchy = - (d_cauchy / l_grad) * cache. g # compute Cauchy point
620
620
u_tmp = u_gauss_newton - u_cauchy # calf of the dogleg
621
621
a = dot (u_tmp, u_tmp)
622
- b = 2 * dot (u_cauchy, u_tmp)
622
+ b = 2 * dot (u_cauchy, u_tmp)
623
623
c = d_cauchy^ 2 - trust_r^ 2
624
- aux = max (b^ 2 - 4 * a * c, 0.0 ) # technically guaranteed to be non-negative but hedging against floating point issues
625
- τ = (- b + sqrt (aux)) / (2 * a) # stepsize along dogleg to trust region boundary
624
+ aux = max (b^ 2 - 4 * a * c, 0.0 ) # technically guaranteed to be non-negative but hedging against floating point issues
625
+ τ = (- b + sqrt (aux)) / (2 * a) # stepsize along dogleg to trust region boundary
626
626
627
627
cache. du = u_cauchy + τ * u_tmp
628
628
end
0 commit comments