Linear Hull for solving system of linear equations

examples/linear_eq.jl

@@ -0,0 +1,32 @@
+# Examples from Luc Jaulin, Michel Kieffer, Olivier Didrit and Eric Walter - Applied Interval Analysis
+
+using IntervalArithmetic, IntervalRootFinding, StaticArrays
+
+A = [4..5 -1..1 1.5..2.5; -0.5..0.5 -7.. -5 1..2; -1.5.. -0.5 -0.7.. -0.5 2..3]
+sA = SMatrix{3}{3}(A)
+mA = MMatrix{3}{3}(A)
+
+b = [3..4, 0..2, 3..4]
+sb = SVector{3}(b)
+mb = MVector{3}(b)
+
+p = fill(-1e16..1e16, 3)
+
+rts = gauss_seidel_interval!(p, A, b, precondition=true) # Gauss-Seidel Method; precondition=true by default
+rts = gauss_seidel_interval!(p, sA, sb, precondition=true) # Gauss-Seidel Method; precondition=true by default
+rts = gauss_seidel_interval!(p, mA, mb, precondition=true) # Gauss-Seidel Method; precondition=true by default
+
+rts = gauss_seidel_interval(A, b, precondition=true) # Gauss-Seidel Method; precondition=true by default
+rts = gauss_seidel_interval(sA, sb, precondition=true) # Gauss-Seidel Method; precondition=true by default
+rts = gauss_seidel_interval(mA, mb, precondition=true) # Gauss-Seidel Method; precondition=true by default
+
+rts = gauss_seidel_contractor!(p, A, b, precondition=true) # Gauss-Seidel Method (Vectorized); precondition=true by default
+rts = gauss_seidel_contractor!(p, sA, sb, precondition=true) # Gauss-Seidel Method (Vectorized); precondition=true by default
+rts = gauss_seidel_contractor!(p, mA, mb, precondition=true) # Gauss-Seidel Method (Vectorized); precondition=true by default
+
+rts = gauss_seidel_contractor(A, b, precondition=true) # Gauss-Seidel Method (Vectorized); precondition=true by default
+rts = gauss_seidel_contractor(sA, sb, precondition=true) # Gauss-Seidel Method (Vectorized); precondition=true by default
+rts = gauss_seidel_contractor(mA, mb, precondition=true) # Gauss-Seidel Method (Vectorized); precondition=true by default
+
+rts = gauss_elimination_interval!(p, A, b, precondition=true) # Gaussian Elimination; precondition=true by default
+rts = gauss_elimination_interval(A, b, precondition=true) # Gaussian Elimination; precondition=true by default

perf/linear_eq.jl

@@ -0,0 +1,64 @@
+using BenchmarkTools, Compat, DataFrames, IntervalRootFinding, IntervalArithmetic, StaticArrays
+
+function randVec(n::Int)
+    a = randn(n)
+    A = Interval.(a)
+    mA = MVector{n}(A)
+    sA = SVector{n}(A)
+    return A, mA, sA
+end
+
+function randMat(n::Int)
+    a = randn(n, n)
+    A = Interval.(a)
+    mA = MMatrix{n, n}(A)
+    sA = SMatrix{n, n}(A)
+    return A, mA, sA
+end
+
+function benchmark(max=10)
+    df = DataFrame()
+    df[:Method] = ["Array", "MArray", "SArray", "Contractor", "ContractorMArray", "ContractorSArray"]
+    for n in 1:max
+        A, mA, sA = randMat(n)
+        b, mb, sb = randVec(n)
+        t1 = @belapsed gauss_seidel_interval($A, $b)
+        t2 = @belapsed gauss_seidel_interval($mA, $mb)
+        t3 = @belapsed gauss_seidel_interval($sA, $sb)
+        t4 = @belapsed gauss_seidel_contractor($A, $b)
+        t5 = @belapsed gauss_seidel_contractor($mA, $mb)
+        t6 = @belapsed gauss_seidel_contractor($sA, $sb)
+        df[Symbol("$n")] = [t1, t2, t3, t4, t5, t6]
+    end
+    a = []
+    for i in 1:max
+        push!(a, Symbol("$i"))
+    end
+    df1 = stack(df, a)
+    dfnew = unstack(df1, :variable, :Method, :value)
+    dfnew = rename(dfnew, :variable => :n)
+    println(dfnew)
+    dfnew
+end
+
+function benchmark_elimination(max=10)
+    df = DataFrame()
+    df[:Method] = ["Gauss Elimination", "Base.\\", "Linear Hull"]
+    for n in 1:max
+        A, mA, sA = randMat(n)
+        b, mb, sb = randVec(n)
+        t1 = @belapsed gauss_elimination_interval($A, $b)
+        t2 = @belapsed gauss_elimination_interval1($A, $b)
+        t3 = @belapsed linear_hull($A, $b)
+        df[Symbol("$n")] = [t1, t2, t3]
+    end
+    a = []
+    for i in 1:max
+        push!(a, Symbol("$i"))
+    end
+    df1 = stack(df, a)
+    dfnew = unstack(df1, :variable, :Method, :value)
+    dfnew = rename(dfnew, :variable => :n)
+    println(dfnew)
+    dfnew
+end

perf/linear_eq_tabular.txt

@@ -0,0 +1,53 @@
+Gauss Seidel
+
+│ Row │ variable │ Array       │ Contractor  │ MArray      │ SArray1     │ SArray2     │
+├─────┼──────────┼─────────────┼─────────────┼─────────────┼─────────────┼─────────────┤
+│ 1   │ n = 1    │ 1.1426e-5   │ 2.4173e-5   │ 1.5803e-5   │ 1.014e-5    │ 9.512e-6    │
+│ 2   │ n = 2    │ 2.4413e-5   │ 4.0504e-5   │ 3.9829e-5   │ 2.4083e-5   │ 2.4083e-5   │
+│ 3   │ n = 3    │ 4.5194e-5   │ 6.4343e-5   │ 9.1236e-5   │ 4.5533e-5   │ 5.1221e-5   │
+│ 4   │ n = 4    │ 7.2639e-5   │ 9.3671e-5   │ 0.00016058  │ 7.0972e-5   │ 6.9937e-5   │
+│ 5   │ n = 5    │ 0.000103281 │ 0.000127489 │ 0.000262814 │ 0.000104823 │ 0.000102729 │
+│ 6   │ n = 6    │ 0.000141315 │ 0.000169992 │ 0.000416266 │ 0.000144021 │ 0.000139086 │
+│ 7   │ n = 7    │ 0.000195001 │ 0.000217059 │ 0.000680848 │ 0.000198164 │ 0.000191145 │
+│ 8   │ n = 8    │ 0.000247193 │ 0.00027535  │ 0.00122559  │ 0.000251089 │ 0.000241972 │
+│ 9   │ n = 9    │ 0.000306262 │ 0.000338028 │ 0.0020223   │ 0.000310961 │ 0.000299153 │
+│ 10  │ n = 10   │ 0.00037781  │ 0.000414073 │ 0.00294134  │ 0.000386235 │ 0.000367335 │
+│ 11  │ n = 11   │ 0.000465017 │ 0.000489114 │ 0.00585246  │ 0.000473425 │ 0.000447176 │
+│ 12  │ n = 12   │ 0.00055403  │ 0.000587152 │ 0.0110894   │ 0.000573395 │ 0.00054274  │
+│ 13  │ n = 13   │ 0.000659994 │ 0.000692402 │ 0.0160662   │ 0.000686098 │ 0.000643391 │
+│ 14  │ n = 14   │ 0.000777718 │ 0.000793626 │ 0.0184595   │ 0.000800443 │ 0.000767576 │
+│ 15  │ n = 15   │ 0.000910174 │ 0.000952207 │ 0.0239422   │ 0.000972855 │ 0.000899402 │
+│ 16  │ n = 16   │ 0.00107785  │ 0.00111443  │ 0.0303295   │ 0.00113467  │ 0.0010867   │
+│ 17  │ n = 17   │ 0.00122326  │ 0.00125759  │ 0.0398127   │ 0.00134961  │ 0.00125675  │
+│ 18  │ n = 18   │ 0.00167176  │ 0.00159463  │ 0.043795    │ 0.00181409  │ 0.00162584  │
+│ 19  │ n = 19   │ 0.0018632   │ 0.00185617  │ 0.0549107   │ 0.00208529  │ 0.00196587  │
+│ 20  │ n = 20   │ 0.0020604   │ 0.00210595  │ 0.0635598   │ 0.00232704  │ 0.00212374  │
+
+10×7 DataFrames.DataFrame
+│ Row │ n  │ Array       │ Contractor  │ ContractorMArray │ ContractorSArray │ MArray      │ SArray      │
+├─────┼────┼─────────────┼─────────────┼──────────────────┼──────────────────┼─────────────┼─────────────┤
+│ 1   │ 1  │ 1.0188e-5   │ 2.5276e-5   │ 0.000781355      │ 0.000764627      │ 1.6959e-5   │ 1.1216e-5   │
+│ 2   │ 10 │ 0.000377853 │ 0.000407817 │ 0.00122723       │ 0.00123259       │ 0.0029247   │ 0.000393896 │
+│ 3   │ 2  │ 2.609e-5    │ 4.3597e-5   │ 0.000802467      │ 0.000796962      │ 4.1039e-5   │ 2.6654e-5   │
+│ 4   │ 3  │ 4.7362e-5   │ 6.317e-5    │ 0.000832029      │ 0.000821289      │ 9.1695e-5   │ 4.7007e-5   │
+│ 5   │ 4  │ 7.0378e-5   │ 9.1502e-5   │ 0.00085429       │ 0.000844048      │ 0.000163604 │ 7.1166e-5   │
+│ 6   │ 5  │ 0.000104813 │ 0.000129112 │ 0.000898626      │ 0.000889969      │ 0.000268756 │ 0.000112476 │
+│ 7   │ 6  │ 0.000142316 │ 0.000168891 │ 0.000941422      │ 0.000922727      │ 0.000430442 │ 0.000146906 │
+│ 8   │ 7  │ 0.000192073 │ 0.000218861 │ 0.00108794       │ 0.000997821      │ 0.00073304  │ 0.000198962 │
+│ 9   │ 8  │ 0.000242832 │ 0.000273341 │ 0.00107159       │ 0.00106769       │ 0.00126407  │ 0.000250956 │
+│ 10  │ 9  │ 0.000308256 │ 0.000333722 │ 0.00116203       │ 0.00116428       │ 0.00196451  │ 0.000314499 │
+
+Gauss Elimination
+
+│ Row │ n  │ Base.\\    │ Gauss Elimination │ Linear Hull │
+├─────┼────┼────────────┼───────────────────┼─────────────┤
+│ 1   │ 1  │ 1.9415e-6  │ 8.40867e-6        │ 3.94457e-6  │
+│ 2   │ 10 │ 7.1184e-5  │ 7.8812e-5         │ 4.6999e-5   │
+│ 3   │ 2  │ 3.34387e-6 │ 1.012e-5          │ 5.30567e-6  │
+│ 4   │ 3  │ 4.548e-6   │ 1.0955e-5         │ 6.5832e-6   │
+│ 5   │ 4  │ 7.256e-6   │ 1.3845e-5         │ 1.1115e-5   │
+│ 6   │ 5  │ 1.0878e-5  │ 1.7847e-5         │ 1.4564e-5   │
+│ 7   │ 6  │ 1.5747e-5  │ 2.1526e-5         │ 1.9468e-5   │
+│ 8   │ 7  │ 2.296e-5   │ 3.109e-5          │ 2.7597e-5   │
+│ 9   │ 8  │ 3.4512e-5  │ 4.0778e-5         │ 3.7466e-5   │
+│ 10  │ 9  │ 4.8627e-5  │ 5.7028e-5         │ 5.3766e-5   │

perf/slopes.jl

@@ -0,0 +1,116 @@
+using BenchmarkTools, Compat, DataFrames, IntervalRootFinding, IntervalArithmetic, StaticArrays
+
+function benchmark_results()
+    f = [] # Example functions from Dietmar Ratz - An Optimized Interval Slope Arithmetic and its Application
+    push!(f, x->((x + sin(x)) * exp(-x^2)))
+    push!(f, x->(x^4 - 10x^3 + 35x^2 - 50x + 24))
+    push!(f, x->((log(x + 1.25) - 0.84x) ^ 2))
+    push!(f, x->(0.02x^2 - 0.03exp(-(20(x - 0.875))^2)))
+    push!(f, x->(exp(x^2)))
+    push!(f, x->(x^4 - 12x^3 + 47x^2 - 60x - 20exp(-x)))
+    push!(f, x->(x^6 - 15x^4 + 27x^2 + 250))
+    push!(f, x->(atan(cos(tan(x)))))
+    push!(f, x->(asin(cos(acos(sin(x))))))
+
+    s = interval(0.75, 1.75)
+    df = DataFrame()
+
+    df[:Method] = ["Automatic Differentiation", "Slope Expansion"]
+    for n in 1:length(f)
+
+        t1 = ForwardDiff.derivative(f[n], s)
+        t2 = slope(f[n], s, mid(s))
+        df[Symbol("f" * "$n")] = [t1, t2]
+    end
+    a = []
+    for i in 1:length(f)
+        push!(a, Symbol("f" * "$i"))
+    end
+    df1 = stack(df, a)
+    dfnew = unstack(df1, :variable, :Method, :value)
+    dfnew = rename(dfnew, :variable => :Function)
+    println(dfnew)
+    dfnew
+end
+
+function benchmark_time()
+    f = []
+    push!(f, x->((x + sin(x)) * exp(-x^2)))
+    push!(f, x->(x^4 - 10x^3 + 35x^2 - 50x + 24))
+    push!(f, x->((log(x + 1.25) - 0.84x) ^ 2))
+    push!(f, x->(0.02x^2 - 0.03exp(-(20(x - 0.875))^2)))
+    push!(f, x->(exp(x^2)))
+    push!(f, x->(x^4 - 12x^3 + 47x^2 - 60x - 20exp(-x)))
+    push!(f, x->(x^6 - 15x^4 + 27x^2 + 250))
+    push!(f, x->(atan(cos(tan(x)))))
+    push!(f, x->(asin(cos(acos(sin(x))))))
+
+    s = interval(0.75, 1.75)
+    df = DataFrame()
+    df[:Method] = ["Automatic Differentiation", "Slope Expansion"]
+    for n in 1:length(f)
+
+        t1 = @belapsed ForwardDiff.derivative($f[$n], $s)
+        t2 = @belapsed slope($f[$n], $s, mid($s))
+        df[Symbol("f" * "$n")] = [t1, t2]
+    end
+    a = []
+    for i in 1:length(f)
+        push!(a, Symbol("f" * "$i"))
+    end
+    df1 = stack(df, a)
+    dfnew = unstack(df1, :variable, :Method, :value)
+    dfnew = rename(dfnew, :variable => :Function)
+    println(dfnew)
+    dfnew
+end
+
+struct SlopesMulti
+    f::Function
+    x::IntervalBox
+    c::Vector
+    sol::Matrix{Interval}
+end
+
+function benchmark_multi()
+
+    rts = SlopesMulti[]
+    f(x, y) = SVector(x^2 + y^2 - 1, y - 2x)
+    f(X) = f(X...)
+    X = (-6..6) × (-6..6)
+    c = [0.0, 0.0]
+    push!(rts, SlopesMulti(f, X, c, [-6..6 -6..6; -2.. -2 1..1]))
+
+    function g(x)
+        (x1, x2, x3) = x
+        SVector(    x1^2 + x2^2 + x3^2 - 1,
+                    x1^2 + x3^2 - 0.25,
+                    x1^2 + x2^2 - 4x3
+                )
+    end
+
+    X = (-5..5)
+    XX = IntervalBox(X, 3)
+    cc = [0.0, 0.0, 0.0]
+    push!(rts, SlopesMulti(g, XX, cc, [-5..5 -5..5 -5..5; -5..5 0..0 -5..5; -5..5 -5..5 -4.. -4]))
+
+    function h(x)
+        (x1, x2, x3) = x
+        SVector(    x1 + x2^2 + x3^2 - 1,
+                    x1^2 + x3 - 0.25,
+                    x1^2 + x2 - 4x3
+                )
+    end
+
+    XXX = IntervalBox(1..5, 2..6, -3..7)
+    ccc = [3.0, 4.0, 2.0]
+    push!(rts, SlopesMulti(h, XXX, ccc, [1..1 6..10 -1..9; 4..8 0..0 1..1; 4..8 1..1 -4.. -4]))
+
+    for i in 1:length(rts)
+        println("\nFunction $i")
+        println("Slope Expansion: ")
+        println(DataFrame(slope(rts[i].f, rts[i].x, rts[i].c)))
+        println("\nJacobian: ")
+        println(DataFrame(ForwardDiff.jacobian(rts[i].f, rts[i].x)))
+    end
+end

perf/slopes_tabular.txt

@@ -0,0 +1,82 @@
+Results
+
+│ Row │ Function │ Automatic Differentiation │ Slope Expansion       │
+├─────┼──────────┼───────────────────────────┼───────────────────────┤
+│ 1   │ f1       │ [-5.44573, 0.886249]      │ [-2.79917, 0.0521429] │
+│ 2   │ f2       │ [-87.6875, 77.0625]       │ [-43.875, 38.25]      │
+│ 3   │ f3       │ [-0.474861, 0.787211]     │ [-0.159199, 0.432842] │
+│ 4   │ f4       │ [-2.97001, 21.0701]       │ [0.04, 0.326667]      │
+│ 5   │ f5       │ [2.63258, 74.8333]        │ [6.03135, 33.2205]    │
+│ 6   │ f6       │ [-94.5871, 115.135]       │ [-38.9908, 65.5595]   │
+│ 7   │ f7       │ [-279.639, 167.667]       │ [-146.852, 67.0665]   │
+│ 8   │ f8       │ [-∞, ∞]                   │ [1, 2]                │
+│ 9   │ f9       │ [-∞, ∞]                   │ [1.3667, ∞]           │
+
+Time
+<Using abstract types>
+│ Row │ Function │ Automatic Differentiation │ Slope Expansion │
+├─────┼──────────┼───────────────────────────┼─────────────────┤
+│ 1   │ f1       │ 8.847e-6                  │ 2.4192e-5       │
+│ 2   │ f2       │ 3.0109e-5                 │ 7.6821e-5       │
+│ 3   │ f3       │ 6.1046e-6                 │ 9.447e-6        │
+│ 4   │ f4       │ 1.6145e-5                 │ 3.8827e-5       │
+│ 5   │ f5       │ 8.29233e-6                │ 2.2186e-5       │
+│ 6   │ f6       │ 3.0895e-5                 │ 7.9524e-5       │
+│ 7   │ f7       │ 3.0847e-5                 │ 7.6188e-5       │
+<After using concrete types>
+│ Row │ Function │ Automatic Differentiation │ Slope Expansion │
+├─────┼──────────┼───────────────────────────┼─────────────────┤
+│ 1   │ f1       │ 8.34067e-6                │ 2.1679e-5       │
+│ 2   │ f2       │ 2.9766e-5                 │ 7.194e-5        │
+│ 3   │ f3       │ 6.0278e-6                 │ 7.827e-6        │
+│ 4   │ f4       │ 1.7114e-5                 │ 3.5235e-5       │
+│ 5   │ f5       │ 8.36567e-6                │ 2.0747e-5       │
+│ 6   │ f6       │ 3.0061e-5                 │ 7.273e-5        │
+│ 7   │ f7       │ 3.0198e-5                 │ 7.2014e-5       │
+│ 8   │ f8       │ 7.43125e-6                │ 8.046e-6        │
+│ 9   │ f9       │ 7.118e-6                  │ 7.8705e-6       │
+
+
+Multi-Dimensional Case:
+Function 1
+Slope Expansion:
+│ Row │ x1       │ x2      │
+├─────┼──────────┼─────────┤
+│ 1   │ [-6, 6]  │ [-6, 6] │
+│ 2   │ [-2, -2] │ [1, 1]  │
+
+Jacobian:
+│ Row │ x1        │ x2        │
+├─────┼───────────┼───────────┤
+│ 1   │ [-12, 12] │ [-12, 12] │
+│ 2   │ [-2, -2]  │ [1, 1]    │
+
+Function 2
+Slope Expansion:
+│ Row │ x1      │ x2      │ x3       │
+├─────┼─────────┼─────────┼──────────┤
+│ 1   │ [-5, 5] │ [-5, 5] │ [-5, 5]  │
+│ 2   │ [-5, 5] │ [0, 0]  │ [-5, 5]  │
+│ 3   │ [-5, 5] │ [-5, 5] │ [-4, -4] │
+
+Jacobian:
+│ Row │ x1        │ x2        │ x3        │
+├─────┼───────────┼───────────┼───────────┤
+│ 1   │ [-10, 10] │ [-10, 10] │ [-10, 10] │
+│ 2   │ [-10, 10] │ [0, 0]    │ [-10, 10] │
+│ 3   │ [-10, 10] │ [-10, 10] │ [-4, -4]  │
+
+Function 3
+Slope Expansion:
+│ Row │ x1     │ x2      │ x3       │
+├─────┼────────┼─────────┼──────────┤
+│ 1   │ [1, 1] │ [6, 10] │ [-1, 9]  │
+│ 2   │ [4, 8] │ [0, 0]  │ [1, 1]   │
+│ 3   │ [4, 8] │ [1, 1]  │ [-4, -4] │
+
+Jacobian:
+│ Row │ x1      │ x2      │ x3       │
+├─────┼─────────┼─────────┼──────────┤
+│ 1   │ [1, 1]  │ [4, 12] │ [-6, 14] │
+│ 2   │ [2, 10] │ [0, 0]  │ [1, 1]   │
+│ 3   │ [2, 10] │ [1, 1]  │ [-4, -4] │

src/IntervalRootFinding.jl

@@ -18,12 +18,16 @@ export
     derivative, jacobian,  # reexport derivative from ForwardDiff
     Root, is_unique,
     roots, find_roots,
-    bisect, newton1d
+    bisect, newton1d, slope,
+    slope_newton1d, linear_hull,
+    gauss_seidel_interval, gauss_seidel_interval!,
+    gauss_seidel_contractor, gauss_seidel_contractor!,
+    gauss_elimination_interval, gauss_elimination_interval!
 
 export isunique, root_status
 
 
-import IntervalArithmetic.interval
+import IntervalArithmetic: interval, wideinterval
 
 
 
@@ -58,7 +62,8 @@ include("complex.jl")
 include("contractors.jl")
 include("roots.jl")
 include("newton1d.jl")
-
+include("linear_eq.jl")
+include("slopes.jl")
 
 
 gradient(f) = X -> ForwardDiff.gradient(f, SVector(X))