On 04/17/2014 10:07 AM, Matthias Felleisen wrote:
Neil, a two-layer design, using option contracts, could be a solution here.
Interesting! So we could have both math/matrix and (say)
math/unchecked-matrix, where the former enables the option contracts
(e.g. that check matrix condition numbers), and the latter doesn't. Each
could re-export the same implementation.
Enforcing low floating-point error using contracts requires solving a
few research-level problems, though.
1. Not only is floating-point error not compositional, it's not even
monotone. For example, for large arguments, `fllog` reduces error:
> (require math)
> (define a (flstep 1000.0 5)) ; 1000 approx., 5 ulps error
> (define b (flstep 1002.0 5)) ; 1002 approx., 5 ulps error
> (flulp-error (* a b) (* 1000 1002))
10.0 ; Error for * is additive
> (flulp-error (fllog (* a b))
(bigfloat->real (bflog (bf (* 1000 1002)))))
0.2789786104670918 ; Whoa, just rounding error
(Admittedly, decreasing error is less likely than increasing error
because most operations at least add rounding error.)
Pretending error is monotone would rule out useful high-error
approximations. Dealing with floating-point error's nonmonontonicity
would make it hard to determine blame, or when to raise exceptions.
2. A robust floating-point implementation usually approximates its
function in more than one way. One common strategy is to first try a
reasonably fast approximation that has low error on a large domain, and
switch to a slower but more accurate approximation if the fast one's
estimated error is high. Can contracts support this strategy? Can they
support it efficiently?
(Efficiently switching to another approximation is probably less of a
concern in linear algebra, where most operations are O(n^2) or O(n^3).)
3. To a user who expects flonums to be just like real numbers but maybe
with a bit of error, an approximation error contract violation - even if
it blames exactly the right expression - will seem capricious and
inexplicable. "It won't let me compute the dot product of this vector
and the result of solving this system of equations? What's wrong with
this library?"
4. How do we express the invariants of a function's floating-point
implementation as contracts?
5. A call site doesn't always need the absolute minimum error a
floating-point function can provide. Further, how much error it can
tolerate may be contextual. How do we express a call site's
approximation error requirements?
Again, I think these are interesting research problems. They may also be
some reasons why nobody has tried this yet (as far as I know). The first
three are the reasons I hesitate to use contracts to enforce low error
right now.
------------------------------------------------------------------------
On Apr 16, 2014, at 3:18 PM, Konrad Hinsen <konrad.hin...@fastmail.net> wrote:
Any linear-algebra library meant for "serious" work does not do any
non-essential work "just" for error checking. Performance is king.
Users are supposed to know the pitfalls and take care to avoid them.
Unfortunately, many don't.
It is incredible how often I find this scenario. Just yesterday I had
a colleague tell me (on a related subject) that "this is how you test
people who want to work in this field. If they can't figure out these
assumptions, they are not well suited. Flush them."
My proposal had been to create a type-like system to check some invariant
about geometric relations. His response in the end was "yeah, it may help
more novices to overcome this obstacle, but the field wont' care in the end.
Experts don't make these mistakes."
At this age, I am 1000% convinced that experts make these mistakes too, but
what can we say.
The most convincing argument I've made to myself is that *I* make these
mistakes even when I'm trying my best to be careful. Unfortunately, that
argument is convincing to others only if we can all agree I'm an expert,
and I still haven't convinced myself of that.
Neil ⊥
____________________
Racket Users list:
http://lists.racket-lang.org/users
