Reenable and fix example 5.5 from Smiley et al. #153

gwater · 2020-04-26T21:12:52Z

No description provided.

gwater · 2020-04-26T21:19:29Z

fix for #86

gwater · 2020-04-26T22:18:07Z

unfortunately, the example takes to long to solve and some CI tests time out. So we might have to create a faster test for CI

gwater · 2020-04-26T22:24:59Z

implements the idea described in JuliaArrays/StaticArrays.jl#450

dpsanders · 2020-04-26T22:44:49Z

src/contractors.jl

    J = jacobian(X)
-
+    if contains_zero(det(J))


Isn't this very expensive? This should basically come out of the \ or LU factorization.

I cant imagine the determinant is as expensive as LU - its just a fixed number of multiplicaitons and additions without any memory allocations

As far as I understand, the determinant does exactly an LU decomposition.

(ah, maybe not if this is for small StaticArrays?)

https://github.com/JuliaLang/julia/blob/master/stdlib/LinearAlgebra/src/generic.jl#L1538

The fact that this check is even necessary is disturbing as we overload \ for intervals with custom methods in

https://github.com/JuliaIntervals/IntervalRootFinding.jl/blob/master/src/linear_eq.jl

Either something is broken there or somehow the wrong method is picked.

Maybe line 165 needs to have StaticVector instead of StaticArray?

We should definitely test if this is actually getting called.

maybe it makes more sense to call gaussian_elimination explicitely (instead of the \ operator) to make sure and obvious that no external algorithm should be called here

The fact that this check is even necessary is disturbing as we overload \ for intervals with custom methods in

https://github.com/JuliaIntervals/IntervalRootFinding.jl/blob/master/src/linear_eq.jl

Either something is broken there or somehow the wrong method is picked.

I can confirm that the internal method is not getting called – in fact it is broken for StaticArrays. And fixing the Array / Vector mixup is not sufficient either…

Are we just not actually importing \ from Base?

gwater · 2020-04-27T11:44:21Z

Ok, given that there apparently is no way to test whether LU decomposition will work, I dont see a solution. The fundamental problem remains that IntervalArithmetic tolerates silent failures in external algorithms. I spent a few hours generalizing (and cleaning up) IntervalRootFinding to support NumberIntervals and then the solution is easy:

function 𝒩(f::Function, jacobian::Function, X::AbstractVector{T}, α) where T # multidimensional Newton operator
    m = T.(mid.(X, α))
    J = jacobian(X)
    try
        return convert(typeof(X), m .- (J \ f(m)))
    catch
        return X .± Inf
    end
end

I've tested it: The exception is triggered (disturbingly often) and eventually all roots in Smiley 5.5 are discovered.

gwater · 2020-04-27T21:23:40Z

ok, this version uses the internal algorithm and finally finds all the roots for Smiley 5.5

gwater · 2020-04-28T09:25:26Z

not sure I trust the internal solver – for the matrix discussed above, it basically just returns this hard-coded value from

IntervalRootFinding.jl/src/linear_eq.jl

Line 90 in 80328b4

x .= -1e16..1e16

gwater · 2020-04-28T09:29:51Z

not sure I trust the internal solver – for the matrix discussed above, it basically just returns this hard-coded value from

IntervalRootFinding.jl/src/linear_eq.jl

Line 90 in 80328b4

x .= -1e16..1e16

It's also breaking the Julia naming convention – this function is marked as mutating but it doesn't

IntervalRootFinding.jl/src/linear_eq.jl

Line 102 in 80328b4

    
           function gauss_elimination_interval!(x::AbstractArray, A::AbstractMatrix, b::AbstractArray; precondition=true)

- fix type signature for ldiv operator - use similar instead of copy to get mutable arrays - remove arbitrary interval restriction

gwater · 2020-04-28T14:55:29Z

cleaned up the commit log a bit. this PR contains a few more changes than strictly necessary, lmk if I should separate them into other PRs

dpsanders

Many thanks for fixing this and sorry for the long delay.

I added a couple of comments but other than that LGTM!

dpsanders · 2020-06-24T19:39:13Z

test/roots.jl

@@ -62,7 +62,7 @@ end

    # Bisection
    rts = roots(f, X, Bisection, 1e-3)
-    @test length(rts) == 4


I do get 4 boxes here (even though there are only 2 roots).

Yes, it gives an incorrect answer. Marking it as @test_broken makes transparent that Bisection is unreliable and probably should never be used in practise. In addition, it will notify us if Bisection ever starts to give the correct result here. (I admit, this is not really related to the topic of this PR–I think I just noticed this inconsistency while I was working on it.)

I see what you mean, but I don't agree that the result is unreliable. It is reliably giving you a set of boxes that are guaranteed to contain all of the roots in the original box, and correctly telling you that it is unable to guarantee whether or not there are actually roots there.

The test is testing what the function returns, not what we would like it to ideally return ;)

Ok, that makes sense. But testing for the result 4 seems overly specific if any number above 2 would acceptable. Maybe it would be better to test that the exact solutions are contained in whatever boxes are returned? I've created a separate PR #156 for that.

dpsanders · 2020-06-24T19:45:54Z

src/linear_eq.jl

+    A = similar(A)
+    b = similar(b)
+    A .= _A
+    b .= _b


Is this different / better than _A = A; A = copy(_A)?

If I recall correctly, calling similar() on a StaticArray gives a MutableArray which is required for the setindex!() calls below.

dpsanders · 2020-06-26T00:28:12Z

Thanks a lot!

gwater force-pushed the dev_smiley55 branch from d034db2 to 49ded71 Compare April 26, 2020 21:16

dpsanders reviewed Apr 26, 2020

View reviewed changes

gwater force-pushed the dev_smiley55 branch from 49ded71 to 23d73a0 Compare April 27, 2020 21:14

gwater force-pushed the dev_smiley55 branch 2 times, most recently from 623fff3 to 39994d9 Compare April 28, 2020 08:36

gwater added 6 commits April 28, 2020 16:45

fix gaussian elimination for static arrays

d4008c4

- fix type signature for ldiv operator - use similar instead of copy to get mutable arrays - remove arbitrary interval restriction

use internal ldiv for Newton contractor

a2ce4a6

mark 2D bisection test as broken

f24d844

give CI tests more time

fdb0771

enable Krawczyk for Smiley examples / tests

3b4733c

bump patch version

779ba37

gwater force-pushed the dev_smiley55 branch from e6fab47 to 779ba37 Compare April 28, 2020 14:49

dpsanders approved these changes Jun 24, 2020

View reviewed changes

revert change to bisection test

7bc322e

gwater mentioned this pull request Jun 25, 2020

Make bisection test more expressive #156

Merged

dpsanders merged commit 9272e39 into JuliaIntervals:master Jun 26, 2020

Kolaru mentioned this pull request Jul 2, 2020

Use gauss_elimination_interval to solve A*b = x in contractors #111

Open

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Reenable and fix example 5.5 from Smiley et al. #153

Reenable and fix example 5.5 from Smiley et al. #153

gwater commented Apr 26, 2020

gwater commented Apr 26, 2020

gwater commented Apr 26, 2020

gwater commented Apr 26, 2020

dpsanders Apr 26, 2020

gwater Apr 26, 2020

dpsanders Apr 26, 2020

dpsanders Apr 26, 2020

dpsanders Apr 26, 2020

Kolaru Apr 27, 2020

dpsanders Apr 27, 2020

gwater Apr 27, 2020

gwater Apr 27, 2020

dpsanders Apr 27, 2020

gwater commented Apr 27, 2020

gwater commented Apr 27, 2020

gwater commented Apr 28, 2020 •

edited

Loading

gwater commented Apr 28, 2020

gwater commented Apr 28, 2020

dpsanders left a comment

dpsanders Jun 24, 2020

gwater Jun 24, 2020

dpsanders Jun 24, 2020 •

edited

Loading

gwater Jun 25, 2020

dpsanders Jun 24, 2020

gwater Jun 24, 2020

dpsanders commented Jun 26, 2020

Reenable and fix example 5.5 from Smiley et al. #153

Reenable and fix example 5.5 from Smiley et al. #153

Conversation

gwater commented Apr 26, 2020

gwater commented Apr 26, 2020

gwater commented Apr 26, 2020

gwater commented Apr 26, 2020

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

Choose a reason for hiding this comment

Choose a reason for hiding this comment

gwater commented Apr 27, 2020

gwater commented Apr 27, 2020

gwater commented Apr 28, 2020 • edited Loading

gwater commented Apr 28, 2020

gwater commented Apr 28, 2020

dpsanders left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dpsanders Jun 24, 2020 • edited Loading

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 Jun 26, 2020

gwater commented Apr 28, 2020 •

edited

Loading

dpsanders Jun 24, 2020 •

edited

Loading