Use 10D problem from#108 (#242)

Author: Kolaru

benchmark/Project.toml

@@ -2,5 +2,6 @@
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
+IntervalRootFinding = "d2bf35a9-74e0-55ec-b149-d360ff49b807"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"

benchmark/benchmarks.jl

@@ -1,33 +1,32 @@
-include("../examples/smiley_examples.jl")
-using .SmileyExample22, .SmileyExample52, .SmileyExample54, .SmileyExample55
-
 using BenchmarkTools
 using ForwardDiff
 using IntervalArithmetic
 using IntervalRootFinding
-using StaticArrays
-
 import Random
+using StaticArrays
 
 const SUITE = BenchmarkGroup()
 
 Random.seed!(0)  # Seed the RNG to get consistent results
 tol = 1e-10
 
-include("dietmar_ratz_functions.jl")
+## Smiley and Chun examples
+include("../examples/smiley_examples.jl")
+using .SmileyExample22, .SmileyExample52, .SmileyExample54, .SmileyExample55
 
 S = SUITE["Smiley"] = BenchmarkGroup()
-for example in (SmileyExample22, SmileyExample52, SmileyExample54) #, SmileyExample55)
+
+for example in (SmileyExample22, SmileyExample52, SmileyExample54, SmileyExample55)
     s = S[example.title] = BenchmarkGroup()
     for contractor in (Newton, Krawczyk)
         s[string(contractor)] = @benchmarkable roots($(example.f), $(example.region) ; contractor = $contractor, abstol = $tol, infer_root_type = false)
     end
 end
 
 
+## Rastrigin function
 S = SUITE["Rastigrin stationary points"] = BenchmarkGroup()
 
-# Rastrigin function:
 const A = 10
 
 f(x, y) = 2A + x^2 - A*cos(2π*x) + y^2 - A*cos(2π*y)
@@ -43,22 +42,11 @@ for contractor in (Newton, Krawczyk)
 end
 
 
-S = SUITE["Linear equations"] = BenchmarkGroup()
-
-sizes = (2, 5, 10)
-
-for n in sizes
-    s = S["n = $n"] = BenchmarkGroup()
-    M = interval.(randn(n, n))
-    b = interval.(randn(n))
-
-    # s["Gauss seidel"] = @benchmarkable gauss_seidel_interval($M, $b)
-    # s["Gauss seidel contractor"] = @benchmarkable gauss_seidel_contractor($M, $b)
-    # s["Gauss elimination"] = @benchmarkable gauss_elimination_interval($M, $b)
-end
-
+## Dietmar-Ratz functions
+include("dietmar_ratz_functions.jl")
 
 S = SUITE["Dietmar-Ratz"] = BenchmarkGroup()
+
 X = interval(0.75, 1.75)
 
 for (k, dr) in enumerate(dr_functions)
@@ -70,3 +58,12 @@ for (k, dr) in enumerate(dr_functions)
         end
     end
 end
+
+# 10 dimensional problem
+include("../examples/10_dimensional.jl")
+
+S = SUITE["10 dimensional"] = BenchmarkGroup()
+
+for contractor in (Newton, Krawczyk)
+    S[string(contractor)] = @benchmarkable roots($f10d, $X10d_close ; contractor = $contractor, abstol = $tol, max_iteration = 1_000_000)
+end

examples/10_dimensional.jl

@@ -0,0 +1,41 @@
+using IntervalArithmetic
+using IntervalRootFinding
+using StaticArrays
+
+# 10 dimensional problem from https://discourse.julialang.org/t/solvers-fail-on-nonlinear-problem-which-has-solutions/20051
+function f10d(X)
+    p = [3600, 18, 18, 3600, 3600, 3600, 1800, 3600, 18, 18, 18, 1800, 18, 36, 11, 180, 0.7, 0.4, 30, 0.2, 4, 4.5, 0.4]
+    x1, x2, x3, x4, x5, x6, x7, x8, x9, x10 = X
+
+    return SVector(
+        -p[17] * x1 + -2 * p[1] * (x1 ^ 2 / 2) + 2 * p[2] * x2 + p[21] * p[18] * ((1 + p[19] * x7) / (p[20] + x7)),
+        -p[17] * x2 + p[1] * (x1 ^ 2 / 2) + -1 * p[2] * x2 + -1 * p[4] * x2 * x4 + p[9] * x3 + p[14] * x3 + -1 * p[6] * x2 * x7 + p[11] * x8,
+        -p[17] * x3 + p[4] * x2 * x4 + -1 * p[9] * x3 + -1 * p[5] * x3 * x4 + p[10] * x5 + -1 * p[14] * x3 + p[15] * x5 + p[7] * x8 * x4 + -1 * p[12] * x3 * x7,
+        -p[17] * x4 + -1 * p[4] * x2 * x4 + p[9] * x3 + -1 * p[5] * x3 * x4 + p[10] * x5 + -1 * p[7] * x8 * x4 + p[12] * x3 * x7 + p[16] * x9 + p[22] * p[18] * ((1 + p[19] * x7) / (p[20] + x7)),
+        -p[17] * x5 + p[5] * x3 * x4 + -1 * p[10] * x5 + -1 * p[15] * x5,
+        -p[17] * x6 + p[14] * x3 + p[15] * x5 + -1 * p[8] * x6 * x10 + p[13] * x9,
+        -p[17] * x7 + -1 * p[6] * x2 * x7 + p[11] * x8 + p[7] * x8 * x4 + -1 * p[12] * x3 * x7 + p[18] * ((1 + p[19] * x7) / (p[20] + x7)),
+        -p[17] * x8 + p[6] * x2 * x7 + -1 * p[11] * x8 + -1 * p[7] * x8 * x4 + p[12] * x3 * x7,
+        -p[17] * x9 + p[8] * x6 * x10 + -1 * p[13] * x9 + -1 * p[16] * x9,
+        x10 + x9 - p[23]
+    )
+end
+
+# Approximate solution
+sol10d = SVector(0.11584484298713685, 
+  0.9455230206055336,   
+  1.4935460680512724,   
+  0.014733756971650408, 
+  2.6673387628301497,   
+ 15.827943289607886,    
+  0.03777227940579977,  
+  5.088785613357363,    
+  0.3986099863617728,   
+  0.0013900136199035832
+)
+
+# This input is close enough to the solution to converge in less than 1e6 iterations
+X10d_close = SVector{10}(interval(x - 0.0002, x + 0.0002) for x in sol10d)
+
+# This input is too large, and requires more than 1e6 iterations to converge
+X10d_large = SVector{10}(fill(interval(0, 100), 10))

test/roots.jl

@@ -133,6 +133,19 @@ end
     end
 end
 
+@testset "10D roots" begin
+    include("../examples/10_dimensional.jl")
+
+    for contractor in newtonlike_methods
+        rts = roots(f10d, X10d_large ; contractor, max_iteration = 50_000)
+        @test any(X -> all(in_interval.(sol10d, X.region)), rts)
+
+        rts = roots(f10d, X10d_close ; contractor, max_iteration = 1_000_000)
+        @test length(rts) == 1
+        @test isunique(rts[1])
+    end
+end
+
 @testset "Dimension mismatch" begin
     f21(xy) = [xy[1]^2 - 2]
     f23(xy) = [xy[1]^2 - 2, xy[2]^2 - 3, xy[1] + xy[2]]