Error in tests when changing bissection lines

[matbesancon]

In the bisection function in utils,

mid = (left + right) / 2
if left == mid || right == mid
    return mid
end

We should be able to use isapprox(mid, left) instead of equality since it should be enough as convergence criterion, but the change yields errors in the tests, for example:

Exponential Cone Projections: Test Failed at /home/mbesancon/Documents/projects_julia/MathOptSetDistances.jl/test/projections.jl:126
  Expression: _test_proj_exp_cone_help(x, atol; dual = false)

and

Exponential Cone Projections: Test Failed at /home/mbesancon/Documents/projects_julia/MathOptSetDistances.jl/test/projections.jl:129
  Expression: _test_proj_exp_cone_help(x, atol; dual = true)
Stacktrace:
 [1] top-level scope at /home/mbesancon/Documents/projects_julia/MathOptSetDistances.jl/test/projections.jl:129
 [2] top-level scope at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/Test/src/Test.jl:1115
 [3] top-level scope at /home/mbesancon/Documents/projects_julia/MathOptSetDistances.jl/test/projections.jl:83

I left it as-is but might be missing something here.

[matbesancon]

@tjdiamandis would you have time to take a look? I'm not sure why the tests fail with isapprox, it could be a sign

[tjdiamandis]

Sure, I can take a look later today.

[tjdiamandis]

This error appears due to the exponential of the solution to the root finding problem, i.e. x comes from the bisection search and we return ((x - 1)*r + s)/(x^2 - x + 1) * [x; 1.0; exp(x)]. Even with the scaling factor, the exp can turn a small numerical error in the solution to a larger error in the projection. Seems to happen when x is greater than ~10.

For example, the bisection search involved in the projection of

v = [0.06487712363946631, -0.6752808310718523, 0.5195365780773017]

produces solutions that differ by 1.27e-8, the resulting scale factor ((x - 1)*r + s)/(x^2 - x + 1) differs by only 6.88e-12, yet the third term introduces error since the exponential of the two solutions (both around 11.42) differ by 1.16e-3. This ultimately introduces a norm error above the 1e-7 tolerance in the tests.

I'd advocate keeping accuracy to floating point precision. The search stops anyway when the abs(f_mid) < tol, so at most the floating point precision stopping criterion adds ~25 extra function evals / iterations. This case should only occur if the function changes very rapidly near the solution, in which case it makes sense to do anyway. Let me know what you think.

Note: the "ground truth" is from SCS, which I believe uses a Newton–Raphson method to solve the root finding problem with an abstol on the solution of 1e-8. So we could alternatively consider loosening the tolerance on abs(f_mid). I also may implement this method anyway in an attempt to deal with the numerical errors I'm getting in the power cone #35.

[matbesancon]

Thanks a lot for looking into it. I'll try with tighter tolerances for the first check than tol

[matbesancon]

Increasing the number of tests run yields errors in the exp projection, some are errors between the projected vector and the JuMP SCS solution, another one is more problematic:

v = [-0.9181930672817855, 0.0006532070373337695, -1.2227350895604592]

hits the point

Failure to find bracketing interval for exp cone projection.

of case 4

[matbesancon]

I was also trying to change the test to:

isapprox(norm(x - px) ^ 2, objective_value(model), atol = tol)

but get very high differences, maybe I got something wrong

[odow]

Couldn't you just change the function to:

function _bisection(
    f::Function,
    left,
    right;
    max_iters::Int = 500,
    atol::Float64 = 1e-14,
)
    f_left, f_right = f(left), f(right)
    for i = 1:max_iters
        if signbit(f_left) == signbit(f_right)
            error(
                "Interval `[$(left), $(right)]` became non-bracketing at " *
                "iteration $(i)."
            )
        end
        mid = 0.5 * (left + right)
        f_mid = f(mid)
        if abs(f_mid) < atol
            return mid
        elseif signbit(f_mid) == signbit(f_left)
            left, f_left = mid, f_mid
        else
            right, f_right = mid, f_mid
        end
    end
    return nothing
end

I don't see why you need the extra check?

[tjdiamandis]

The upper bound evaluates to ub = -1405.6692821767872 for

v = [-0.9181930672817855, 0.0006532070373337695, -1.2227350895604592]

but h(ub) = -Inf. This is problematic, since the theoretical results suggest that h should be strictly increasing on this interval (lb, ub), where lb = -Inf here.

Edit: the upper bound uses r/s, and s is small here. If I increase s by a factor of 3, I no longer get an error. Maybe it's numerical?

Edit 2: It looks like it's numerical: (r - x*s)*exp(-x) evaluates to Inf. Will need to handle this case for negative, large magnitude x.

[tjdiamandis]

I don't see why you need the extra check?

Agreed. I think I added the early stopping (abs(f(mid)) < tol) in later on. I'd advocate removing the if left == mid || right == mid condition

[matbesancon]

Thanks both, I'll remove the first check and re-arrange the iterations as @odow suggested. For the interval finding issue I'll experiment a bit

[matbesancon]

Mmmh, that's gonna be tricky:

julia> h(ub+0.0000000000005)
Inf

julia> h(ub+0.0000000000001)
-Inf

[matbesancon]

This seems indeed to be a precision error:

julia> ub = s > 0 ? big(r)/big(s) : Inf
-1405.669282176787092512159672269710140073349064255701228886964385615708308712058

julia> h(ub)
1.222735089560459176283302440424449741840362548828125

Which sounds more sensible as a result

[matbesancon]

Ok the numerics are extremely tricky with this one, even with big floats, the lower bound quickly takes a huge dive:

julia> h(ub-0.0000000000000000000000000000000000000000000000000000000000000000000000001)
-1.041713637010071423028633113342996800418033643807315172561975501815968710900029e+528

julia> h(ub)
1.222735089560459176283302440424449741840362548828125

[matbesancon]

This also applies when converting v to big before defining h

[tjdiamandis]

@matbesancon do you think we can go ahead and close this how that #43 is merged?

[matbesancon]

yup agreed, thanks!