Unexpected results from Brent minimizer.

[alexgian]

Unexpected results from Brent minimizer:

I came across this problem when trying to find the nearest point on a parametric curve (in this case an ellipse) to a random given point in its vicinity. I tried minimizing the distance formula, the difference between the two points.

Obviously you'd expect such distances around an ellipse to be symmetrical about the axes, but I found a large section of the ellipse, pizza-slice shaped, gave unlikely to very unlikely results. I tried the code in GJS' scmultils, and the answers were correct there. Checked the default values for eps tolerances, adjusted them, but the results were still way out.

Here's the code, both scmutils and Emmy.

(define (fparam t)
  (up (* 6 (cos t)) (* 3 (sin t))))

(define ((dist-from param-func point) t)
  (abs (- point (param-func t))))

(define (constrain-func point)
  (fparam (first (minimize (dist-from fparam point) 0 2pi))))

(constrain-func (up 1 2))
(constrain-func (up 1 -2))
(constrain-func (up -1 1))
(constrain-func (up -1 -1))
(constrain-func (up -4 1))
(constrain-func (up -4 -1))
(constrain-func (up 5 1))
(constrain-func (up 5 -1))
(constrain-func (up 6 2))
(constrain-func (up 6 -2))

gives

1 ]=> #| fparam |#

1 ]=> #| dist-from |#

1 ]=> #| constrain-func |#

1 ]=> #|
(up 1.0875785060182517 2.9503039247358527)
(up 1.0874951745568635 -2.9503116041070787)
(up -1.197540201106018 2.939638475507446)
(up -1.1974119961407907 -2.939651531708228)
(up -4.5555647684665725 1.9523594469452918)
(up -4.555564768466676 -1.9523594469452314)
(up 5.351607713907495 1.3564931695778815)
(up 5.351607054569401 -1.356493819879307)
(up 5.35694860003839 1.351212575480515)
(up 5.35699007907589 -1.351171463275709)
|#

All correct and symmetric as it should be. OTOH with the Emmy code:

(def tau (* 2 pi))

(defn fparam [t]
  [(* 6 (cos t))  (* 3 (sin t))])

;; distance from a point to that parametric function
(defn dist-from [parametric-function [x y]]
  (fn [t]  (emmy.env/abs (- [x y] (parametric-function t) ))))

;; get the point on curve at minimal distance
(defn constrain-func [[x y]]
    (fparam ((emmy.env/brent-min (dist-from fparam [x y]) 0 tau
                                 {:absolute-threshold (sqrt emmy.value/machine-epsilon) :relative-threshold 1.0e-5})
 :result)))

running in REPL...

(constrain-func (up 1 2))
(constrain-func (up 1 -2))
(constrain-func (up -1 1))
(constrain-func (up -1 -1))
(constrain-func (up -4 1))
(constrain-func (up -4 -1))
(constrain-func (up 5 1))
(constrain-func (up 5 -1))
(constrain-func (up 6 2))
(constrain-func (up 6 -2))

; second values of each pair are incorrect
=> [1.0875839178752857 2.9503034259877]
=> [1.7364359622122425 2.871619323184048]
=> [-1.1975658217144087 2.939635866168695]
=> [-1.506227489636631 2.9039317635519146]
; this is the only correct pair
=> [-4.5555412374816395 1.9523731734484224]
=> [-4.55554123748163 -1.952373173448428]
; totally goofy
=> [5.3516073498306564 1.3564935286648523]
=> [5.999999999999998 6.31257342499718E-8]
=> [5.356961624488731 1.3511996663849362]
=> [5.999999999999998 6.31257342499718E-8]

Symmetry is just not there. However, in every case the first value of each pair matches pretty well with the scmutils results.

Basically, I wanted to apply this function to an emmy-viewers.mafs/moveable-point to constrain it to move on an ellipse that is specified parametrically. This is a bit more elegant than using x-y coordinates and conditionals (which kinda works, but feels klutzy).

[sritchie]

I used this opportunity to create a nice visual debugging notebook in Maria, thanks to @mhuebert: https://2.maria.cloud/gist/e2e9fa24033f5b50d6cc7af44f22243f

The difference with scmutils was the initial condition for Brent. Apache Math-Commons used the midpoint between the left and right ends of the interval, while the original algorithm, and scmutils, used a golden section cut. It doesn't matter usually, but in this case the difference pushed you to a different local minimum.

I'm going to fix / change this to match scmutils in #175 . The real thing you want to use is #114 , and that would be an awesome contribution @alexgian if you're willing to work on it.

[alexgian]

Does this mean that in different circumstances choosing the midpoint between left and right would have been the better choice? Or is a Golden Section cut always preferred? If the the former is the case, then it needs be an option, I think. Yeah, I'll look at #114, Sam.

[sritchie]

@alexgian I think this is an edge case for the algorithm and there's no general answer here. I did provide a new option for :initial-guess, but I adjusted the default to the golden cut in #175 .

I'll close this now since this is fixed in #175 , and that's under review now by @littleredcomputer .

Thanks @alexgian !