System of Linear Equations

[codecov-io]

Codecov Report

Merging #67 into master will increase coverage by 6.37%. The diff coverage is 85.13%.

@@            Coverage Diff             @@
##           master      #67      +/-   ##
==========================================
+ Coverage   56.14%   62.52%   +6.37%     
==========================================
  Files          10       10              
  Lines         431      467      +36     
==========================================
+ Hits          242      292      +50     
+ Misses        189      175      -14

Impacted Files	Coverage Δ
src/IntervalRootFinding.jl	4.65% <ø> (+2.61%)	:arrow_up:
src/linear_eq.jl	85.13% <85.13%> (ø)
src/contractors.jl	90% <0%> (-0.48%)	:arrow_down:
src/newton.jl	0% <0%> (ø)	:arrow_up:
src/krawczyk.jl	0% <0%> (ø)	:arrow_up:
src/bisect.jl
src/roots.jl	88.13% <0%> (+2.19%)	:arrow_up:
src/complex.jl	83.33% <0%> (+8.33%)	:arrow_up:

---

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update 25ebe4c...7746380. Read the comment docs.

[dpsanders]

Thanks! I don't see the benchmarks though.

[eeshan9815]

I manually benchmarked, and shared the results on #intervals in slack

[dpsanders]

Please include at least the benchmarking script (using BenchmarkTools) in a new directory perf.

[eeshan9815]

I have included the script and results for n = 1:25 cc @dpsanders

[dpsanders]

That's the kind of thing I had in mind, thanks. You could instead collect the data using @belapsed (I think) and output a table using Dataframes or something like that.

[eeshan9815]

Okay, I'll try it out! So which version of the function should I retain? I don't think MArray improved the performance a lot.

[dpsanders]

Try à version using IntervalBox and setindex.

[dpsanders]

IntervalBox is immutable like an SVector. So modifying it actually requires making a new copy and changing the relevant entry using setindex (without !)

[dpsanders]

This needs tests.

[eeshan9815]

@dpsanders Travis and AppVeyor both use the latest tagged version and don't have extended_div, and hence are failing.

[lbenet]

I tried to tag a new version but it was causing failures in ValidatedNumerics. I did some upgrades, but I think we should plan carefully how to make patches-versions, so they don't break everything. See this comment.

[dpsanders]

Actually the matrix version could be pretty efficient when using SMatrix and SVector everywhere. Have you tried that?

[eeshan9815]

Done. I think this version gives the best benchmarks.

[dpsanders]

So it seems just using normal arrays is fine?

[eeshan9815]

So it seems just using normal arrays is fine?

More or less, yes. Should I remove the other implementations and update the PR?

[dpsanders]

You should squash down the implementations as I suggested; I think you should leave all of them for now.

Could you please add the example from the Jaulin book in the examples folder?

Actually you should keep a separate repo, e.g. in your GitHub account or in JuliaIntervals, with a Jupyter notebook with usage examples and explanation.

[eeshan9815]

@dpsanders Is this accurate?

[eeshan9815]

@dpsanders I've added the gaussian elimination with preconditioning and benchmarking against Base.\

[eeshan9815]

ping @dpsanders

[dpsanders]

This and your other PRs need rebasing on current master, sorry!

[eeshan9815]

Done rebasing and tests pass on the current IntervalArithmetic master

[dpsanders]

There is still (or again) a conflict in test/runtests.jl.

[dpsanders]

(sorry, wrong button)

[eeshan9815]

Oh that occurred due to the merging of #70. I have rebased again.

[dpsanders]

@eeshan9815 Tests are still failing. Do you know what is happening?

[eeshan9815]

Tests fail as it uses extended_div which is not in the latest tagged version of IntervalArithmetic. Tests pass locally if the current master is used.

[dpsanders]

Ah great, thanks

[eeshan9815]

I have redefined \ for Intervals. Currently, it uses the gauss_elimination_interval method with precondition=true. This is due to the advantages of preconditioning such as reducing dependence mentioned in various books such as Global Optimisation using Interval Analysis by Eldon Hansen.

[eeshan9815]

Since \ is same as gauss_elimination_interval, for which tests exist already, is there a need to define tests for \ separately? Also, comparisions exist in the perf/linear_eq_tabular.txt file.

[dpsanders]

Right, but the tests are also for future-proofing, eg if we change the implementation of \ again. So you can just add a couple of tests to make sure the output of \ and gaussian elim are the same.

[dpsanders]

@eeshan9815 Did you actually check that the definition of \ in the package is being called? It doesn't seem to be. You need to import Base: \.

[dpsanders]

On the example in the testset "Linear Equations", the new definition seems to be worse than using Base's \.

[eeshan9815]

I have imported Base.\ in IntervalRootFinding.jland it seems to be using gauss_elimination_interval instead of Base.\ Also, preconditioning widens the solution set by a little bit. It's a trade-off mentioned in the Hansen book as well.

[dpsanders]

Hmm, I see, sorry. I don't know why it didn't work for me, then.

What's the advantage of preconditioning, then? I thought it was supposed to narrow, not widen, the result.

[dpsanders]

After checking out your PR, I get

julia> 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]
3×3 Array{IntervalArithmetic.Interval{Float64},2}:
 [4, 5]            [-1, 1]                    [1.5, 2.5]
    [-0.5, 0.5]   [-7, -5]                [1, 2]
    [-1.5, -0.5]       [-0.700001, -0.5]  [2, 3]

julia> b = [3..4, 0..2, 3..4]
3-element Array{IntervalArithmetic.Interval{Float64},1}:
 [3, 4]
 [0, 2]
 [3, 4]

julia> A \ b
3-element Array{IntervalArithmetic.Interval{Float64},1}:
 [-1.81928, 1.16873]
 [-0.41407, 1.72523]
  [0.700232, 3.42076]

julia> @which A \ b
\(A::AbstractArray{T,2} where T, B::Union{AbstractArray{T,1}, AbstractArray{T,2}} where T) in Base.LinAlg at linalg/generic.jl:820

which shows that the \ from Base is being used. Not sure why.

[dpsanders]

Ignore me, I seem not to have managed to checkout the latest version, apologies.

[dpsanders]

No, don't ignore me -- you defined \ in test/linear_eq.jl instead of in the source code!

[dpsanders]

No problem, thanks!

[dpsanders]

Running gauss_seidel seems to modify something it shouldn't? The result of A \ b shouldn't change in the following session:

julia> IntervalRootFinding.gauss_elimination_interval(A, b)
2-element Array{IntervalArithmetic.Interval{Float64},1}:
 [-130.228, 167.728]
  [-60.0001, 267.273]

julia> IntervalRootFinding.gauss_elimination_interval(A, b, precondition=false)
2-element Array{IntervalArithmetic.Interval{Float64},1}:
 [-120, 90]
 [-60, 240]

julia> IntervalRootFinding.gauss_seidel_interval(A, b, precondition=false)
2-element Array{IntervalArithmetic.Interval{Float64},1}:
 [-5e+15, 5.00001e+15]
 [-1e+16, 1e+16]

julia> A
2×2 Array{IntervalArithmetic.Interval{Float64},2}:
 [2, 3]  [0, 1]
 [1, 2]  [1, 3]

julia> b
2-element Array{IntervalArithmetic.Interval{Float64},1}:
   [0, 120]
 [-60, 240]

julia> A \ b
2-element Array{IntervalArithmetic.Interval{Float64},1}:
 [-1e+16, 1e+16]
 [-1e+16, 1e+16]

[dpsanders]

I think this needs more tests. E.g. the gauss_seidel methods give huge pointless answers for the single test case.

[dpsanders]

@eeshan9815 There seems to be a real test failure here.

[eeshan9815]

The test failure appears to be in a test for 2D roots, giving error messages similar to the ones faced when the SVector was initially made a data member of IntervalBox. It should work after merging #83

[eeshan9815]

The incorrect modification issue has also been resolved.

julia> IntervalRootFinding.gauss_elimination_interval(A, b)
2-element Array{IntervalArithmetic.Interval{Float64},1}:
 [-130.228, 167.728] 
  [-60.0001, 267.273]

julia> IntervalRootFinding.gauss_elimination_interval(A, b, precondition=false)
2-element Array{IntervalArithmetic.Interval{Float64},1}:
 [-120, 90]
 [-60, 240]

julia> IntervalRootFinding.gauss_seidel_interval(A, b, precondition=false)
2-element Array{IntervalArithmetic.Interval{Float64},1}:
 [-120.001, 90.0001] 
  [-60.0001, 240.001]

julia> A
2×2 Array{IntervalArithmetic.Interval{Float64},2}:
 [2, 3]  [0, 1]
 [1, 2]  [2, 3]

julia> b
2-element Array{IntervalArithmetic.Interval{Float64},1}:
  [0, 120]
 [60, 240]

julia> A \ b
2-element Array{IntervalArithmetic.Interval{Float64},1}:
 [-130.228, 167.728] 
  [-60.0001, 267.273]

[dpsanders]

@eeshan9815 Tests are now passing. Is this ready?

[eeshan9815]

Yes, I think so. I added 88 tests (mostly randomly generated) to make sure the methods are working.

[dpsanders]

Thanks very much!