System of Linear Equations #67

eeshan9815 · 2018-05-14T14:03:22Z

No description provided.

codecov-io · 2018-05-14T14:21:16Z

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%)`	⬆️
src/linear_eq.jl	`85.13% <85.13%> (ø)`
src/contractors.jl	`90% <0%> (-0.48%)`	⬇️
src/newton.jl	`0% <0%> (ø)`	⬆️
src/krawczyk.jl	`0% <0%> (ø)`	⬆️
src/bisect.jl
src/roots.jl	`88.13% <0%> (+2.19%)`	⬆️
src/complex.jl	`83.33% <0%> (+8.33%)`	⬆️

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 · 2018-05-14T23:24:10Z

Thanks! I don't see the benchmarks though.

eeshan9815 · 2018-05-15T06:09:58Z

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

dpsanders · 2018-05-15T06:59:08Z

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

eeshan9815 · 2018-05-15T09:36:18Z

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

dpsanders · 2018-05-15T10:05:53Z

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 · 2018-05-15T12:05:11Z

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 · 2018-05-15T12:17:54Z

Try à version using IntervalBox and setindex.

dpsanders · 2018-05-15T12:19:01Z

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

dpsanders · 2018-05-16T06:28:27Z

src/linear_eq.jl

+function gauss_seidel_interval_static1{T, N}(A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}; precondition=true, maxiter=100)
+
+    n = size(A, 1)
+    x = @MVector fill(-1e16..1e16, n)


@SVector?

I realized only x has to be modified, so I tried a version where A and b are SVectors, and that showed improvement.

Please try a version where x is an SVector and you use setindex to "modify" it (i.e. replace x with a new SVector).

Done.
I think this version gives the best benchmarks.

dpsanders · 2018-05-16T06:29:39Z

This needs tests.

eeshan9815 · 2018-05-17T14:15:35Z

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

lbenet · 2018-05-17T16:00:24Z

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 · 2018-05-17T21:32:07Z

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

eeshan9815 · 2018-05-17T23:17:48Z

Done.
I think this version gives the best benchmarks.

dpsanders · 2018-05-17T23:23:27Z

src/linear_eq.jl

+"""
+function gauss_seidel_interval_static2!{T, N}(x::SVector{N, Interval{T}}, A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}; precondition=true, maxiter=100)
+
+    precondition && ((M, r) = preconditioner_static2(A, b))


What happens is precondition is false? Where are M and r defined then?

dpsanders · 2018-05-18T13:57:51Z

So it seems just using normal arrays is fine?

eeshan9815 · 2018-05-18T14:01:10Z

So it seems just using normal arrays is fine?

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

dpsanders · 2018-05-18T19:10:55Z

perf/linear_eq.jl

+        t2 = @belapsed gauss_seidel_interval_static($mA, $mb)
+        t3 = @belapsed gauss_seidel_interval_static1($sA, $sb)
+        t4 = @belapsed gauss_seidel_interval_static2($sA, $sb)
+        t5 = @belapsed gauss_seidel_contractor($A, $b)


Can you try this one with the static vectors too?

gauss_seidel_contractor with static vectors? Okay

dpsanders · 2018-05-18T19:11:40Z

perf/linear_eq.jl

+        t3 = @belapsed gauss_seidel_interval_static1($sA, $sb)
+        t4 = @belapsed gauss_seidel_interval_static2($sA, $sb)
+        t5 = @belapsed gauss_seidel_contractor($A, $b)
+        df[Symbol("n = $n")] = [t1, t2, t3, t4, t5]


Why do you need Symbol here? Can't you just use a string?

julia> df["test"] = [1, 2, 3, 4] ERROR: MethodError: no method matching setindex!(::DataFrames.DataFrame, ::Array{Int64,1}, ::String)

The DataFrames package takes a symbol as the index.

dpsanders · 2018-05-18T19:12:08Z

perf/linear_eq.jl

+        t3 = @belapsed gauss_seidel_interval_static1($sA, $sb)
+        t4 = @belapsed gauss_seidel_interval_static2($sA, $sb)
+        t5 = @belapsed gauss_seidel_contractor($A, $b)
+        df[Symbol("n = $n")] = [t1, t2, t3, t4, t5]


Why don't you in fact just make n another column?

dpsanders · 2018-05-18T19:12:55Z

perf/linear_eq_tabular.txt

+│ 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  │


Actually n seems to be redundant since it's just the same as the row.

The row numbers need not be sequential due to the stack operations, That's why I kept the column with the value of n

dpsanders · 2018-05-18T19:13:40Z

src/linear_eq.jl

+
+    n = size(A, 1)
+
+    for iter in 1:maxiter


You should add @inbounds to the outer for loop in each of these functions.

dpsanders · 2018-05-18T19:17:18Z

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.

dpsanders · 2018-05-18T19:14:48Z

src/linear_eq.jl

+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_static1!{T, N}(x::MVector{N, Interval{T}}, A::SMatrix{N, N, Interval{T}}, b::SVector{N, Interval{T}}; precondition=true, maxiter=100)


There are several versions of the function whose code is identical.
Just make a single version and do much looser type annotations as AbstractArray etc.

(Or no type annotations at all)

eeshan9815 · 2018-05-19T13:34:07Z

@dpsanders Is this accurate?

eeshan9815 · 2018-06-29T16:50:19Z

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 · 2018-06-29T16:54:12Z

Ah great, thanks

eeshan9815 · 2018-06-30T01:10:21Z

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.

dpsanders

Thanks! Could you please add Tests for \ and comparisons to Base.\

eeshan9815 · 2018-06-30T19:41:40Z

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 · 2018-06-30T20:59:37Z

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 · 2018-07-01T03:17:28Z

@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 · 2018-07-01T03:18:48Z

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

eeshan9815 · 2018-07-01T06:41:00Z

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 · 2018-07-01T15:08:49Z

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 · 2018-07-01T15:11:19Z

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 · 2018-07-01T15:13:39Z

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

dpsanders · 2018-07-01T15:15:33Z

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

dpsanders · 2018-07-01T15:16:46Z

test/linear_eq.jl

+    end
+end
+
+\(A::StaticMatrix{Interval{T}}, b::StaticArray{Interval{T}}; kwargs...) where T = gauss_elimination_interval(A, b, kwargs...)


These should be moved to the source code.

I think you can define a single method if you use AbstractMatrix and AbstractVector.

These should be moved to the source code.

I'm terribly sorry, I didn't realize.

I think you can define a single method if you use AbstractMatrix and AbstractVector

I tried AbstractArrays earlier, but it doesn't work since it matches its type with the \ for StaticArrays or Arrays instead.

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

Preconditioning doesn't always give a better result. Specifically, for M-matrices, the non-preconditioned version computes the exact hull whereas preconditioning causes some widening.

dpsanders · 2018-07-01T22:47:15Z

No problem, thanks!

dpsanders · 2018-07-09T21:41:20Z

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 · 2018-07-09T22:10:21Z

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

dpsanders · 2018-07-09T22:14:43Z

src/linear_eq.jl

+            Z = extended_div(Y, A[i, i])
+            x[i] = hull((x[i] ∩ Z[1]), x[i] ∩ Z[2])
+        end
+        if all(x .== x¹)


This should probably use a tolerance instead.
Can you write this as all(==, x, x1)?

Doing that gives this error -

julia> A = [2..3 0..1;1..2 2..3] 2×2 Array{IntervalArithmetic.Interval{Float64},2}: [2, 3] [0, 1] [1, 2] [2, 3] julia> all(==, A, A) ERROR: ArgumentError: reduced dimension(s) must be integers Stacktrace: [1] reduced_indices(::Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}, ::Array{IntervalArithmetic.Interval{Float64},2}) at ./reducedim.jl:35 [2] reducedim_initarray(::Array{IntervalArithmetic.Interval{Float64},2}, ::Array{IntervalArithmetic.Interval{Float64},2}, ::Bool, ::Type{Bool}) at ./reducedim.jl:73 [3] mapreducedim at ./reducedim.jl:242 [inlined] [4] all(::Function, ::Array{IntervalArithmetic.Interval{Float64},2}, ::Array{IntervalArithmetic.Interval{Float64},2}) at ./reducedim.jl:583 [5] macro expansion at /home/eeshan/.julia/v0.6/Atom/src/repl.jl:118 [inlined] [6] anonymous at ./<missing>:?

dpsanders · 2018-07-17T02:16:29Z

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

eeshan9815 · 2018-07-17T05:58:54Z

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 · 2018-07-17T11:12:52Z

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 · 2018-07-17T20:20:27Z

@eeshan9815 Tests are now passing. Is this ready?

eeshan9815 · 2018-07-17T20:57:44Z

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

dpsanders · 2018-07-17T21:09:01Z

Thanks very much!

eeshan9815 force-pushed the linear-eq branch from b10e219 to 8ae13da Compare May 14, 2018 14:12

eeshan9815 force-pushed the linear-eq branch from e812548 to 50d532c Compare May 15, 2018 09:56

dpsanders reviewed May 16, 2018

View reviewed changes

eeshan9815 force-pushed the linear-eq branch 2 times, most recently from db737c7 to e3f58c1 Compare May 17, 2018 23:17

dpsanders reviewed May 17, 2018

View reviewed changes

dpsanders reviewed May 18, 2018

View reviewed changes

dpsanders reviewed May 19, 2018

View reviewed changes

Import and redefine \ for intervals

3ed3dfd

dpsanders reviewed Jun 30, 2018

View reviewed changes

eeshan9815 changed the title ~~[WIP] System of Linear Equations~~ System of Linear Equations Jun 30, 2018

Add tests for \

45fbb24

dpsanders reviewed Jul 1, 2018

View reviewed changes

Fix declaration in wrong file

a9065f1

dpsanders reviewed Jul 9, 2018

View reviewed changes

Add more tests and bug fixes

7746380

dpsanders merged commit 5b7a8ab into JuliaIntervals:master Jul 17, 2018

System of Linear Equations #67

System of Linear Equations #67

Conversation

eeshan9815 commented May 14, 2018

codecov-io commented May 14, 2018 • edited Loading

Codecov Report

dpsanders commented May 14, 2018

eeshan9815 commented May 15, 2018

dpsanders commented May 15, 2018

eeshan9815 commented May 15, 2018 • edited Loading

dpsanders commented May 15, 2018

eeshan9815 commented May 15, 2018

dpsanders commented May 15, 2018

dpsanders commented May 15, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dpsanders commented May 16, 2018

eeshan9815 commented May 17, 2018

lbenet commented May 17, 2018

dpsanders commented May 17, 2018

eeshan9815 commented May 17, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dpsanders commented May 18, 2018

eeshan9815 commented May 18, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dpsanders commented May 18, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

eeshan9815 commented May 19, 2018

eeshan9815 commented Jun 29, 2018

dpsanders commented Jun 29, 2018

eeshan9815 commented Jun 30, 2018

dpsanders left a comment

Choose a reason for hiding this comment

eeshan9815 commented Jun 30, 2018

dpsanders commented Jun 30, 2018

dpsanders commented Jul 1, 2018

dpsanders commented Jul 1, 2018

eeshan9815 commented Jul 1, 2018

dpsanders commented Jul 1, 2018

dpsanders commented Jul 1, 2018

dpsanders commented Jul 1, 2018

dpsanders commented Jul 1, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dpsanders commented Jul 1, 2018

dpsanders commented Jul 9, 2018

dpsanders commented Jul 9, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dpsanders commented Jul 17, 2018

eeshan9815 commented Jul 17, 2018 • edited Loading

eeshan9815 commented Jul 17, 2018

dpsanders commented Jul 17, 2018

eeshan9815 commented Jul 17, 2018

dpsanders commented Jul 17, 2018

codecov-io commented May 14, 2018 •

edited

Loading

eeshan9815 commented May 15, 2018 •

edited

Loading

eeshan9815 commented Jul 17, 2018 •

edited

Loading