Fix rootfinding by enforcing bounds of initial x-values

[adrienmaillard]

• Tickets addressed: Closes #1336 #1468
• Review: By commit
• Merge strategy: Merge (no squash)

Description

There was a bug in the rootfinding algorithm. A heuristic value was provided without checking its bounds. Then the algorithm was using it and going off bound, causing it to find a solution that did not exist, violating the mutex conditions.

Verification

Validated against the clipper case, the bug disappeared.

Documentation

None.

Future work

None.

[adrienmaillard]

Clarification about why this happens as I have seen this is not clear. Let's say you want to find the root of f. Your x domain is [1,3], your y domain is [5,6]. You start at x0, provided as an argument of the rootfinding procedure. In our case, it's gonna be 1 (we take the lower bound of the x domain). The algorithm actually needs 2 starting values to really begin. we compute f(x0) = y0, let's say that f(1) = 8. Then we use a heuristic that assumes that the function is constant so we say x1 (the next x to be tested) is x1 = x0 - 8 (as f(x) is assumed to be x+8 and we want f(x) = 0). x1 = -7 which is outside of the x domain [1,3]. This x1 is passed to the nextValueAt as an init value. The domain of values found by nextValueAt was checked... except for the init value. So in the faulty case, it computes f(x1) and returns. f(x1) satisfies constraints on the y domain and rootfinding is done even though x1 is outside of the x domain. The x domain is determined by the application of constraints beforehand, of which mutex is one.

[Mythicaeda]

Clarification about why this happens as I have seen this is not clear. Let's say you want to find the root of f. Your x domain is [1,3], your y domain is [5,6]. You start at x0, provided as an argument of the rootfinding procedure. In our case, it's gonna be 1 (we take the lower bound of the x domain). The algorithm actually needs 2 starting values to really begin. we compute f(x0) = y0, let's say that f(1) = 8. Then we use a heuristic that assumes that the function is constant so we say x1 (the next x to be tested) is x1 = x0 - 8 (as f(x) is assumed to be x+8 and we want f(x) = 0). x1 = -7 which is outside of the x domain [1,3]. This x1 is passed to the nextValueAt as an init value. The domain of values found by nextValueAt was checked... except for the init value. So in the faulty case, it computes f(x1) and returns. f(x1) satisfies constraints on the y domain and rootfinding is done even though x1 is outside of the x domain. The x domain is determined by the application of constraints beforehand, of which mutex is one.

Why does the algorithm want the second value to be f(x)=0 instead of f(x_max)=y? If it took f(x_min) and f(x_max) as its first two values, then it shouldn't ever escape the fn's x-domain (and would remove the need for x_init, assuming of course that the caller isn't optimizing this). Even if x_init was kept, having the second point be the x-domain bound it's further from would still give the alg two points while ensuring it's in the bounds (while avoiding issues should x_init be x_min or x_max.

[dandelany]

@Mythicaeda and I talked this morning & confirmed that the change in this PR is a) safe and b) mitigates the issue reported by users, so I'd like to get this merged for our v2.14.0 release today. I've opened #1490 to continue the conversation from above, answer any remaining questions & determine if there's a deeper improvement to the algo to be made.

[adrienmaillard]

Why does the algorithm want the second value to be f(x)=0 instead of f(x_max)=y? If it took f(x_min) and f(x_max) as its first two values, then it shouldn't ever escape the fn's x-domain (and would remove the need for x_init, assuming of course that the caller isn't optimizing this). Even if x_init was kept, having the second point be the x-domain bound it's further from would still give the alg two points while ensuring it's in the bounds (while avoiding issues should x_init be x_min or x_max.

Our overall use of rootfinding is to find f(x) = y but we transpose this problem to a traditional rootfinding problem (there is one method with the y argument that transposes to another method without the y... because it's gonna be 0). The goal of traditional rootfinding is to find x such that f(x) = 0. Here we add bounds for x and f(x). The bounds for f(x) are just the wiggle room around 0 that we are willing to go with.

And yes it can escape because of the heuristic I talked about. This heuristic is good in most cases but it can sometimes get us outside of the bounds... not anymore though.

[dandelany]

Thanks for the fix @adrienmaillard and for the detailed review @Mythicaeda - discussion continued in #1490. I'm gonna go ahead and merge this for our 2.14 release today.