System of Linear Equations (#67)

* Add Gauss-Seidel method

* Use StaticArrays and add benchmarking

* Add benchmark script

* Add another version with StaticArrays

* contractor approach

* Add tests

* replace deepcopy with copy

* Add test dependency

* another method with a static vector x

* bug fix

* Display benchmark data in a table

* add tabular results

* make results more readable

* generalized the various implementation

* Add Gauss Elimination with preconditioning

* Add benchmarking data

* Add tests and fix typos

* Add examples and bug fix

* Import and redefine \ for intervals

* Add tests for \

* Fix declaration in wrong file

* Add more tests and bug fixes

Author: eeshan9815

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,63 @@
+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.\\"]
+    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)
+        df[Symbol("$n")] = [t1, t2]
+    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_results.txt

@@ -0,0 +1,49 @@
+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
+
+5×3 DataFrames.DataFrame
+│ Row │ n │ Base.\\    │ Gauss Elimination │
+├─────┼───┼────────────┼───────────────────┤
+│ 1   │ 1 │ 1.84e-6    │ 8.127e-6          │
+│ 2   │ 2 │ 3.1615e-6  │ 1.0029e-5         │
+│ 3   │ 3 │ 4.53986e-6 │ 1.1211e-5         │
+│ 4   │ 4 │ 6.8655e-6  │ 1.4342e-5         │
+│ 5   │ 5 │ 1.0773e-5  │ 1.7725e-5         │

perf/linear_eq_tabular.txt

@@ -9,7 +9,7 @@ using ForwardDiff
 using StaticArrays
 
 
-import Base: ⊆, show, big
+import Base: ⊆, show, big, \
 
 import Polynomials: roots
 
@@ -18,7 +18,10 @@ export
     derivative, jacobian,  # reexport derivative from ForwardDiff
     Root, is_unique,
     roots, find_roots,
-    bisect, newton1d, quadratic_roots
+    bisect, newton1d, quadratic_roots,
+    gauss_seidel_interval, gauss_seidel_interval!,
+    gauss_seidel_contractor, gauss_seidel_contractor!,
+    gauss_elimination_interval, gauss_elimination_interval!
 
 export isunique, root_status
 
@@ -43,6 +46,7 @@ include("contractors.jl")
 include("roots.jl")
 include("newton1d.jl")
 include("quadratic.jl")
+include("linear_eq.jl")
 
 
 gradient(f) = X -> ForwardDiff.gradient(f, X[:])

src/IntervalRootFinding.jl

@@ -0,0 +1,160 @@
+"""
+Preconditions the matrix A and b with the inverse of mid(A)
+"""
+function preconditioner(A::AbstractMatrix, b::AbstractArray)
+
+    Aᶜ = mid.(A)
+    B = inv(Aᶜ)
+
+    return B*A, B*b
+
+end
+
+function gauss_seidel_interval(A::AbstractMatrix, b::AbstractArray; precondition=true, maxiter=100)
+
+    n = size(A, 1)
+    x = similar(b)
+    x .= -1e16..1e16
+    gauss_seidel_interval!(x, A, b, precondition=precondition, maxiter=maxiter)
+    return x
+end
+"""
+Iteratively solves the system of interval linear
+equations and returns the solution set. Uses the
+Gauss-Seidel method (Hansen-Sengupta version) to solve the system.
+Keyword `precondition` to turn preconditioning off.
+Eldon Hansen and G. William Walster : Global Optimization Using Interval Analysis - Chapter 5 - Page 115
+"""
+function gauss_seidel_interval!(x::AbstractArray, A::AbstractMatrix, b::AbstractArray; precondition=true, maxiter=100)
+
+    precondition && ((A, b) = preconditioner(A, b))
+
+    n = size(A, 1)
+
+    @inbounds for iter in 1:maxiter
+        x¹ = copy(x)
+        for i in 1:n
+            Y = b[i]
+            for j in 1:n
+                (i == j) || (Y -= A[i, j] * x[j])
+            end
+            Z = extended_div(Y, A[i, i])
+            x[i] = hull((x[i] ∩ Z[1]), x[i] ∩ Z[2])
+        end
+        if all(x .== x¹)
+            break
+        end
+    end
+    x
+end
+
+function gauss_seidel_contractor(A::AbstractMatrix, b::AbstractArray; precondition=true, maxiter=100)
+
+    n = size(A, 1)
+    x = similar(b)
+    x .= -1e16..1e16
+    x = gauss_seidel_contractor!(x, A, b, precondition=precondition, maxiter=maxiter)
+    return x
+end
+
+function gauss_seidel_contractor!(x::AbstractArray, A::AbstractMatrix, b::AbstractArray; precondition=true, maxiter=100)
+
+    precondition && ((A, b) = preconditioner(A, b))
+
+    n = size(A, 1)
+
+    diagA = Diagonal(A)
+    extdiagA = copy(A)
+    for i in 1:n
+        if (typeof(b) <: SArray)
+            extdiagA = setindex(extdiagA, Interval(0), i, i)
+        else
+            extdiagA[i, i] = Interval(0)
+        end
+    end
+    inv_diagA = inv(diagA)
+
+    for iter in 1:maxiter
+        x¹ = copy(x)
+        x = x .∩ (inv_diagA * (b - extdiagA * x))
+        if all(x .== x¹)
+            break
+        end
+    end
+    x
+end
+
+function gauss_elimination_interval(A::AbstractMatrix, b::AbstractArray; precondition=true)
+
+    n = size(A, 1)
+    x = similar(b)
+    x .= -1e16..1e16
+    x = gauss_elimination_interval!(x, A, b, precondition=precondition)
+
+    return x
+end
+"""
+Solves the system of linear equations using Gaussian Elimination,
+with (or without) preconditioning. (kwarg - `precondition`)
+Luc Jaulin, Michel Kieffer, Olivier Didrit and Eric Walter - Applied Interval Analysis - Page 72
+"""
+function gauss_elimination_interval!(x::AbstractArray, a::AbstractMatrix, b::AbstractArray; precondition=true)
+
+    if precondition
+        ((a, b) = preconditioner(a, b))
+    else
+        a = copy(a)
+        b = copy(b)
+    end
+    n = size(a, 1)
+
+    p = zeros(x)
+
+    for i in 1:(n-1)
+        if 0 ∈ a[i, i] # diagonal matrix is not invertible
+            p .= entireinterval(b[1])
+            return p .∩ x
+        end
+
+        for j in (i+1):n
+            α = a[j, i] / a[i, i]
+            b[j] -= α * b[i]
+
+            for k in (i+1):n
+                a[j, k] -= α * a[i, k]
+            end
+        end
+    end
+
+    for i in n:-1:1
+        sum = 0
+        for j in (i+1):n
+            sum += a[i, j] * p[j]
+        end
+        p[i] = (b[i] - sum) / a[i, i]
+    end
+
+    p .∩ x
+end
+
+function gauss_elimination_interval1(A::AbstractMatrix, b::AbstractArray; precondition=true)
+
+    n = size(A, 1)
+    x = fill(-1e16..1e16, n)
+
+    x = gauss_elimination_interval1!(x, A, b, precondition=precondition)
+
+    return x
+end
+"""
+Using `Base.\``
+"""
+function gauss_elimination_interval1!(x::AbstractArray, a::AbstractMatrix, b::AbstractArray; precondition=true)
+
+    precondition && ((a, b) = preconditioner(a, b))
+
+    a \ b
+end
+
+\(A::StaticMatrix{Interval{T}}, b::StaticArray{Interval{T}}; kwargs...) where T = gauss_elimination_interval(A, b, kwargs...)
+\(A::Matrix{Interval{T}}, b::Array{Interval{T}}; kwargs...) where T = gauss_elimination_interval(A, b, kwargs...)

src/linear_eq.jl

@@ -0,0 +1,40 @@
+using IntervalArithmetic, StaticArrays, IntervalRootFinding
+
+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
+
+@testset "Linear Equations" begin
+
+    A = [[2..3 0..1;1..2 2..3], ]
+    b = [[0..120, 60..240], ]
+    x = [[-120..90, -60..240], ]
+
+    for i in 1:10
+        rand_a = randMat(i)[1]
+        rand_b = randVec(i)[1]
+        rand_c = rand_a * rand_b
+        push!(A, rand_a)
+        push!(b, rand_c)
+        push!(x, rand_b)
+    end
+    for i in 1:length(A)
+        for precondition in (false, true)
+            for f in (gauss_seidel_interval, gauss_seidel_contractor, gauss_elimination_interval, \)
+                @test all(x[i] .⊆ f(A[i], b[i]))
+            end
+        end
+    end
+end

test/linear_eq.jl

@@ -7,3 +7,4 @@ include("roots.jl")
 include("newton1d.jl")
 include("quadratic.jl")
 include("test_smiley.jl")
+include("linear_eq.jl")