JuliaIntervals
Replace signbit by interval_signbit, and let signbit(::Interval)=false
The idea is to define interval_signbit (which returns an interval) so copysign behaves properly (cf #273), and fix displaying complex intervals as proposed here.
cc @dpsanders @dlfivefifty @fp4code
Coverage increased (+0.007%) to 89.803% when pulling ce7cfd7c28bd2f4b1042ef2fc68773745961d999 on lb/show_complex into 478d673b61e1c7883bb6c05d08b72c53ada18dd7 on master.
Coverage increased (+0.007%) to 90.129% when pulling 202ea88ec89f66ddbd0af459c6f085ae93254b75 on lb/show_complex into d3a318c6356515ee2d2a0c47ab7af419dea106ea on master.
The question is if this affects anything somewhere else in some function that checks signbit to see if a number is negative.
We are coming up against limits of Julia's type system once again.
Isn't the new interval_signbit the same as the sign function for an Interval?
The question is if this affects anything somewhere else in some function that checks signbit to see if a number is negative.
Within IntervalArithmetics everything is consistent. I have only checked that IntervalRootFinding.jl works fine, and it does.
Incidentally, the standard does mentions (requires?) a sign function; the standard states "sign(x) is −1 if x<0; 0 if x=0; and 1 if x>0. It is discontinuous at 0 in its domain".
Isn't the new interval_signbit the same as the sign function for an Interval?
interval_signbit is identical to signbit in current master. As you noticed, it is almost the same as sign. The subtle difference is that sign returns the signs of the edges of the interval, while interval_signbit returns the signbit, in reversed order. Then, we get sign(-1 .. 1) == -1 .. 1 whereas interval_signbit(-1 .. 1) == 0 .. 1.
I introduced interval_signbit as a simple workaround, to allow having signbit(::Interval) = false which permits the correct display of complex intervals, and also have a function which does what signbit is currently doing in master; that's the reason that one method of copysign and flipsign were modified.
I have checked all packages in JuliaRegistries/General that have IntervalArithmetics.jl as a dependency. I went one by one, looking for signbit using GitHub search utility. I hope this is somewhat thorough.
I have found only two packages which somehow use signbit: https://github.com/JuliaIntervals/MPFI.jl and https://github.com/gwater/NumberIntervals.jl. In MPFI.jl, there is a test that involves an end of an interval, so it tests signbit(::Float64). In NumberIntervals.jl there are a couple of tests that do involve signbit of a NumberInterval; the function is overloaded here.
So I guess we should ask @gwater about this PR.
I have no objection to removing the signbit method from IntervalArithmetic. Since NumberIntervals imports it from Base and implements its own method, I doubt I even need to adapt anything.
Is this ready to merge? I confirm it fixes the show bug for complex intervals.
If I recall the reason I added the signbit overload was QR with intervals. Can we add a test for this?
@dlfivefifty I am in favor to add a test, can you please suggest one?
I just came up with the following, but do not know if this is what you have in mind:
julia> using IntervalArithmetic, LinearAlgebra
julia> A = [Interval(1) Interval(-1); Interval(1) Interval(1)]
2×2 Array{Interval{Float64},2}:
[1, 1] [-1, -1]
[1, 1] [1, 1]
julia> Q = qr(A)
QR{Interval{Float64},Array{Interval{Float64},2}}
Q factor:
2×2 LinearAlgebra.QRPackedQ{Interval{Float64},Array{Interval{Float64},2}}:
[-0.707107, -0.707106] [-0.707107, -0.707106]
[-0.707107, -0.707106] [0.707106, 0.707107]
R factor:
2×2 Array{Interval{Float64},2}:
[-1.41422, -1.41421] [-3.33067e-16, 4.4409e-16]
[0, 0] [1.41421, 1.41422]
julia> (1/sqrt(2)) ∈ Q.Q[2,2]
true
julia> sqrt(2) ∈ Q.R[2,2]
true
Would that be ok?
Does a qr decompositon actually make sense for interval matrices? What does it mean for two interval vectors to be orthogonal?
I do not know if qr makes sense or not for intervals, and you certainly have a good point about questioning the meaning of two interval vectors being orthogonal. My guess is that one could say that the interval vectors contain vectors that are indeed orthogonal.
@dlfivefifty Can you suggest a better test?
QR makes sense in that its the QR decomposition calculated via interval arithmetic. That is, Q and R are intervals containing the true QR decomposition.
I think the test is fine.
QR makes sense in that its the QR decomposition calculated via interval arithmetic. That is, Q and R are intervals containing the true QR decomposition.
I think the test is fine.
I'd like to remind you of this issue which is caused almost exactly by this test case (inverting a small matrix of Intervals through an external package which is unaware of Intervals). I fear, adding a test case will only further encourage users to expect results IntervalArithmetic.jl does not currently deliver (namely, safe matrix inversion through external libraries).
Don’t let the perfect be the enemy of the good: there are cases where interval QR works exceptionally well for inverting matrices, in particular, when A is diagonally dominant.
Also note that \ presumably uses Gaussian elimination which is well-known to be ill-conditioned in small neighbourhoods, which interval arithmetic picks up. QR is a very different beast: it’s always stable, interval arithmetic may be a bit pessimistic in its bounds but not in a catastrophically bad way.
I’d actually propose overriding \ to use QR for interval matrices.
I'm certainly not claiming that it's impossible to invert/decompose a matrix - the problem is that IntervalArithmetic.jl has no way of telling users when it's impossible (again, see the issue I referenced above). NumberIntervals.jl demonstrates how to implement interactions with interval-unaware algorithms safely.
Like I said: don’t let the perfect be enemy of the good. I encountered this issue because I needed it.
Like I said: don’t let the perfect be enemy of the good. I encountered this issue because I needed it.
I guess that is kind of my point: If you need to invert interval matrices, you might benefit from using the NumberIntervals wrapper because it protects you from cases which cannot be solved using Interval arithmetic. And there is no need to test functions from the standard library to guarantee that.
I actually quite like the solution in NumberIntervals. It's a clever way of introducing a 3-valued logic by returning missing when both true and false are correct answers for different parts of an interval (if I understand correctly).
I'm not sure what you mean by "cannot be solved" here.
As we saw in https://github.com/JuliaArrays/StaticArrays.jl/issues/450 sometimes external algorithms give somewhat unpredictable results when confronted with IntervalArithmetic.jl's intervals. Therefore I would say some matrices in conjunction with some external algorithms represent problems which IntervalAirthmetic.jl is unable to solve.
I actually quite like the solution in
NumberIntervals. It's a clever way of introducing a 3-valued logic by returningmissingwhen bothtrueandfalseare correct answers for different parts of an interval (if I understand correctly).
Thanks, it really copies all of IntervalArithmetics logic for set operations and reimplements number operations with the three-value logic you described. When I finally implemented it is was a lot simpler than I had imagined. (It lacks some of IntervalArithmetic's syntactic sugar though.)
I think @gwater has a point: Let's not include the proposed test, simply because it may give a false impression that things work "out of the box", when they may not.
What about checking if the changes proposed in this PR do not break the code which motivated the introduction of signbit? Could you share the code that motivated signbit, @dlfivefifty ?
The code that motivated it was calling QR on an interval valued matrix! So this test is the MWE.
I still don’t see nor believe there’s a problem in QR, which is very different from \
You also don’t want to falsely imply things don’t work out of the box, when they do...
I agree that we should include the qr test. And we should think about whether we should incorporate the use of three valued logic in the future.
We should also probably overload \ in the package as Sheehan suggests
@gwater How does the 3-valued logic affect performance though? Since that introduces a weak type instability, right?
@gwater How does the 3-valued logic affect performance though? Since that introduces a weak type instability, right?
I haven't benchmarked it yet but from what I hear these small Union's are quite performant.