Lambda the Ultimate

You're thinking of exact and approximate numbers as separate types for dispatch, etc. They are not. They are simply different values for ordinary numbers.

The operations are implemented on the representation of ordinary numbers, which is the same for exact and approximate numbers.

The 'exact' bit of the result for an operation like addition or multiplication is simply the 'exact' bit of argument 1 AND the 'exact' bit of argument 2. The 'exact' bit of the results for operations based on Taylor series, irrational numbers, etc, is just #FALSE regardless of whether the arguments are exact.

There's no advantage in doing type dispatch on the value of one varying bit among 32 bytes of other varying bits.

Floating point accuracy is in fact quite good. That's not the issue, in itself. There's a qualitative difference in how much you can trust subsequent calculations based on a number, between knowing that something is the exact answer and knowing that something is very close to the exact answer. The problem with roundoff accuracy is that it is amplified by repetition; the more rounding calculations you do with rounded numbers, the less you can trust your result's accuracy. I'm therefore regarding it as important to know exactly when the switch from exact to approximate answers occurs.

I've already implemented that normalizing step, BTW. It turns out to be surprisingly complex (and therefore slow). :-(

Bear

Complexity of certain operations was one of my first impressions. Does your normalizing step always round to the nearest representable value? If not, is addition commutative?

More food for thought: Have you settled whether -0 behaves differently from +0? It seems that infinities are allowed. Is +∞ equal to -∞?

When "approximately 1" is compared to "exactly 1", does it come out less than, greater than, equal, or something else? (When different values compare equal, it can wreck Haskell containers like Data.Map.)

Are not always amplified by repetition. Sometimes the round off errors cancel each other. The gold example of this is Newton's method, but it happens in other algorithms as well. As a side note, this is also why interval arithmetic fails, as it says every iteration of Newton's method is less accurate than before, when the reality is that it's more accurate (and generally, lots more accurate).

And it doesn't matter when round off error is introduced, once it's introduced, you have problems. And if your algorithm isn't numerically stable (i.e. if errors don't tend to cancel instead of accumulate), you're screwed. And you can't avoid round off error in any finite-bit representation.

In any successive approximation algorithm, as long as the bits lost due to error in each iteration are less than the bits of accuracy gained in that iteration, the inaccuracy-inducing errors can be ignored except that you may need more iterations to reach your desired accuracy.

This is because of the feedback in each iteration; some function that senses an approximate distance between the current approximation and the accurate answer is called, and the distance reported by this function becomes more accurate the closer the approximation is to the actual answer. The crux of this is that it measures the distance based on the current approximation, which automatically takes into account all rounding errors so far affecting it. The inaccuracy introduced by rounding errors is worked out the same way as the inaccuracy in the initial guess.

IOW, the interval arithmetic that fails is a misapplication of interval arithmetic in the first place. If you do it correctly, you call the feedback function separately for the upper bound and lower bound of the current approximation, perform the correction based on each (resulting in two intervals), combine them into a new interval which includes both of the two intervals you've computed, and observe that the new upperbound and lowerbound are closer to each other than those of the previous iteration.

Of course, this is getting a bit dangerously close to the "sufficiently smart compiler" argument, in that you can't simply perform your operations linearly (as with any other type of number); instead, you have to call functions separately on the upperbound and lowerbound. IOW, to do interval arithmetic right you have to support it using algorithmic transformations in the compiler rather than just by extending basic mathematical functions to accept intervals and then using the same algorithms you use with single-valued numbers.

Ray

With an iterative algorithm (where the output of one iteration is fed into the next iteration), the question becomes: do errors tend to cancel, or tend to accumulate? That's the key question: can the system recognize that errors are tending to cancel? The default assumption has to be that errors tend to accumulate- as uncorrelated errors will tend to accumulate. So you have to detect correlations not just within an operation, but between (arbitrarily distant) operations. If you don't detect these correlations, then your estimate of the error bounds of the computation will be much wider than the real ones.

An example of what I mean: summing a long list of numbers is unstable. Having a very large value early in the sequence will cause subsequent small numbers to be rounded off and lose precision. But sorting the list of numbers and adding them smallest to largest is stable. So any system to detect the true error bounds of the calculation has to be able to recognize that we're summing the numbers in increasing order, and shrink the respective bounds.

And if you errors are accumulating, you're screwed- it's just a question of when. If you have an N-bit representation of the number, and are losing 1 bit per iteration, after N iterations you no longer have any (non-error) information left. More bits just increases the N at which you're screwed. Interval arithmetic just says you have no meaningful information left.

But one has to ask- is this the fault of the number representation, or a fault of the algorithm being used?

But one has to ask- is this the fault of the number representation, or a fault of the algorithm being used?

Can I say, "sort of both and sort of neither?"

I'll start with algorithms.

You can, in theory, automagically transform any valid converging numerical algorithm into a valid converging interval algorithm. So in some technical sense, the algorithm is "not the problem" in that the interval algorithm can be derived from (and therefore expressed by) any expression of the numerical algorithm.

But in practice, doing that is hard to implement in a compiler. It's far easier to implement interval arithmetic by just implementing alternate versions of the lowest-level math functions (no fancy analysis of user-written algorithms needed) and then apply the relatively simpler numerical algorithms to intervals -- bottoming out in type dispatch only at the lowest level. And that is what most implementations do, resulting in very rapid loss of precision.

A few of the better implementations provide more complete interval-aware math libraries, in which converging functions (like Taylor series, etc) are provided in separate numerical and interval versions and type dispatch occurs at the highest possible level. To get valid converging interval functions here, you need to call library functions instead of expecting any converging numerical algorithms you write yourself to be valid on intervals. I'm not aware of *any* implementations that do the actual analysis to correctly detect converging numerical algorithms written by the programmer and transform them into converging interval algorithms -- though several symbolic-algebra packages such as Mathematica and Maxima certainly have all the infrastructure you'd need to do so.

So in practical terms, if you want to do correct interval mathematics I think you have to write interval algorithms separately from your numerical algorithms, and explicitly compute the new upperbound and new lowerbound separately. So even though they are technically only transformations of the numerical algorithms, and could theoretically be expressed by the same code, you still have to worry about it and write them.

Representation is also part of the difficulty. Fractional numbers appear in source code mostly in base ten or as ratios, and IEEE floats are represented inside the computer mostly in normalized base two. There is no accurate representation of any negative power of ten (or most other possible denominators) in base two, thus usually no completely accurate representation of most source code fractions inside the computer. Finally, rounding is almost always part of any operation on these numbers, and every time you round, you lose information.

Changing the representation can palliate the problem somewhat, but can't alleviate it. In making the numerical representation A/B * C^10 * D^2, I'm providing a way to directly represent base-ten fractions and rationals with no rounding, and a way to represent the results of most operations on numbers likely to appear in source code with no rounding. So the results of simple calculations in six or twelve terms are likely to be exact. But, as you observe, with any finite representation, it's just a question of when -- how many operations can you do until you're looking at an approximation? How many before you're losing information? Fifty? A hundred? Sooner or later, and likely sooner, the results are going to stop being exact.

Ideally, source code should say what kind of precision the output requires, and the compiler should either find a way to provably provide that amount of precision by marshalling the available representation sizes and operator precisions for intermediate results, or issue a warning and provide error bounds if it can't figure out how. But that's hard, so we're busily specifying one particular method of doing it rather than actually saying what we want done.

By specifying how many bits of the operands and intermediate terms to keep, as most source code is doing now, we are (over)specifying a method, but neither specifying the desired resulting precision nor checking whether the method we specify actually gives the desired resulting precision. And, that's kind of a problem.

Ray

When does your representation stop being precise and start being approximate? Consider the square root function- is it's return value exact, or approximate? The square root of 4 can be computed exactly, the square root of 2, not so much. This would be a good question for you to answer- how does the square root function work?

This highlights what I think of as a core problem with language design, which goes well beyond number representation (and is thus on topic for this forum in a way that arguing about the behavior of interval arithmetic isn't), although it's on display here- the introduction of "do what I mean" magic.

What do I mean by "do what I mean" magic? It's an attempt by the language designer to predict what the programmer meant by some construct, by a complicated (and tending to grow more complicated over time) set of rules. The difference is that the language is trying not only to understand what the programmer explicitly stated, but what the programmer actually meant. Fairly quickly, the rules become complicated enough that the programmer is unlikely to be able to remember them all, or correctly apply them to any sufficiently complex situation, at which point the language decides what to do "by magic". Examples of this (in non-numerical terms) include C++'s function overloading rules and most of perl.

Unfortunately, trying to suss out what the programmer meant is a mug's game. Not the least of which because the programmer themselves often do not have a clear idea of what they mean.

Consider the string "0.3333" entered into a program. What does the programmer mean by that string? Does the programmer mean 3333/10000, or does the programmer mean 1/3? Of course it's obvious, right? OK, so we add some rules. How about the number 0.142857? It's obvious that I mean 1/7, right? Then how about the number 3.14159, that's obviously pi. When you get to the point where 2.9987355 "obviously means" pi - 1/7, you start to see the problem.

At a certain point you have to give up trying to guess what the programmer meant, and just do something. In your terms, at a certain point you have to switch over to approximate results- with errors and round offs and the whole baggage train of problems. Adding "do what I mean" magic never works- it's much more worthwhile to have a simple to understand (or at least well understood) system, which the programmer can hold in their heads and apply.

IEEE floating point is at least well understood at this point. Those programmers who care to learn it already have, and that is knowledge that can be applied to just about any language. Those programmers who do not care to learn how numbers work, and expect the computer to just do what they mean, will be disappointed- the only question is when will they be disappointed?

Most of the motivation for a decimal datatype, as I understand it, is the desire to be able to do precise financial computations, with predictable rounding in the event that they cannot be precise (as with annuities, for example). This is not a small sector of the problem space, although it is also not very high status, but I think the world is generally catching up, so that decimal types are implemented or proposed for most programming languages, and not just COBOL.

The theory is that programmers (should be) clever enough to figure out whether they want a decimal or binary floating point number.

Your point about numeric literals in programs is extremely well-taken. It's unfortunate that 1/7 does not give you the right answer in most programming languages, but we've probably all gotten used to writing 1.0/7.0 and M_PI_4 or math.pi, or whatever.

Consider the square root function- is it's return value exact, or approximate? The square root of 4 can be computed exactly, the square root of 2, not so much. This would be a good question for you to answer- how does the square root function work?

You're right, it's instructive.

The square root function, as specified, returns exact results if given an exact argument whose square root is representable, returns an inexact result but introduces no additional error if given an inexact argument whose square root is representable, and gives an inexact result accurate to within practical limits (a factor of one over the maximum denominator) otherwise.

It starts by creating a "square-normalized" representation of the argument, by multiplying the numerator by 1, 2, 10, or 20, and adding 0 or 1 to the binary power and decimal power, as needed to achieve a representation in which the binary and decimal powers are even. If an overflow results from this multiplication or either of these additions, the attempt to produce an exact result is aborted.

Our possible result must then have binary and decimal powers equal to half those respective powers of the square-normalized form.

Next it uses Hensel's Lemma to find the integer square roots of the numerator and denominator of the square-normalized form, if they exist. If they exist, then they are the numerator and denominator of the square root and an exact square root is returned. However, it marks the result as an inexact number if the argument is inexact.
If they don't exist, then the attempt to produce an exact result is aborted.

If the attempt to produce an exact result has been aborted, it continues refining the result using Newton's method, likely for only one or two iterations, and returns an inexact result instead.

Consider the string "0.3333" entered into a program. What does the programmer mean by that string? Does the programmer mean 3333/10000, or does the programmer mean 1/3? Of course it's obvious, right? OK, so we add some rules. How about the number 0.142857? It's obvious that I mean 1/7, right? Then how about the number 3.14159, that's obviously pi. When you get to the point where 2.9987355 "obviously means" pi - 1/7, you start to see the problem.

Well, I'm presuming, where anything the programmer says is valid code with an unambiguous semantic meaning, that the programmer means exactly what he says. If he meant 1/3, he should have typed "1/3" instead of "0.3333" and if he meant 1/7, he should have typed "1/7" instead of "0.142857." All of these would be read as exact numbers, whether in thirds, sevenths, ten-thousandths, or millionths.

Pi, however, is another issue. It has no exact representation. The closest available value (marked inexact, of course) is provided as a constant named "π." π however is only used as a constant, so it is never provided as output, nor is output ever shown as a multiple of π.

A programmer who thinks he means to take the usual shortcuts in representation or use the most general algorithm for something when a more precise algorithm exists or is more efficient, may suppose that the runtime is refusing to do what he asks, but that programmer will be wrong.

I intend to include "constraint relaxation clauses" that permit inexact methods or representations with less-stringent specification to be used within their dynamic extent, and "requirement clauses" that specify constraints such as a maximum acceptable average time for operations within their dynamic extent, but you'd have to use both (ie, permit inexact operations/representations AND provide a reason why the exact operations/representations are unacceptable) to get, eg, inexact square roots on IEEE reals, where exact square roots on the extended numeric representation are possible -- and even then, the exact method and representations would probably be used if the specified timings could be achieved that way.

Ray

You can, in theory, automagically transform any valid converging numerical algorithm into a valid converging interval algorithm. So in some technical sense, the algorithm is "not the problem" in that the interval algorithm can be derived from (and therefore expressed by) any expression of the numerical algorithm.

Is that strictly true? Even in theory, in addition to the practical cases that you discuss, it seems reasonable that there would be valid numerical algorithms that converge for reasons that we do not understand. If we have no understanding of why they are convergent then the transformation into a convergent interval algorithm would require unknown steps. I can think of no guarantee that these unknown steps should exist - transforming a program from one set of operations onto another set of operations while preserving a property of the calculated result can be a hard problem.

William Kahan has a nice example in a paper somewhere; I can try and find it if you want although he has rather a lot of papers on IEEE floating point behaviour so it would take a while. It is only about three lines long, according to a reasonable understanding of stability and error bounds it *should* be divergent, yet there is an unusual correlation in the values that causes the calculation to remain convergent.

This is not unusual in numerical analysis: the behaviour of pure mathematical objects is difficult to model, but in many cases there is centuries of previous work to build upon. As soon as we introduce state / steps / iteration it becomes significantly harder and there is a much shorter literature to draw upon.

Even in theory, in addition to the practical cases that you discuss, it seems reasonable that there would be valid numerical algorithms that converge for reasons that we do not understand. If we have no understanding of why they are convergent then the transformation into a convergent interval algorithm would require unknown steps. I can think of no guarantee that these unknown steps should exist.

This ... is true. We can transform numerical algorithms into interval algorithms only when the numerical algorithms and the reasons why they converge are well-understood.

I presumed, without thinking thoroughly, that programming, by intent, would use only algorithms whose convergence properties are well-understood. But that just isn't true. Programmers have a long history of invoking Magic Which They Do Not Understand.

Trying to imagine an interval algorithm that converges on a region in something like, say, the Collatz Fractal, or how it could be directly constructed from the corresponding numerical algorithms, is a brain-breaking experience.

Ray

I'd be interested in that example of Kahan's if it turns out to be easy for you to find. Thanks!

But I did stumble across an interesting discussion in this section of the wiki page on rounding. The argument that Kahan is quoted in making is somewhat stronger than what I said. The final paragraph of the section is somewhat frustrating as it cannot provide an explicit example as the argument is not constructive - it looks like an application of Rice's theorem.

Thanks.

Thanks.

I have a hypothesis on duplicate posts we see lately. When I last did it, clicking "preview comment" button seemed to actually post, and then the "post comment" button posted again as a duplicate.

That's the problem ;)

Summing smaller numbers first is better because of the reason you gave (loss of precision) but it's not a matter of errors cancelling out. I also wouldn't describe sorting before summing in a single pass as the difference between stable and unstable. If you add up enough numbers, you can eventually have a total that's much larger than the numbers you're adding and run into the same issues, even if they're sorted. You can do better in theory in those situations by doing the addition hierarchically, but in practice, a linear running sum is going to be kept on the FPU at 80 bits even though your numbers are all only doubles, which make many of the loss of precision issues around summing numbers moot unless you're summing trillions of numbers. Edit: OK, my mental math was off for that last calculation

Granted. Most of the examples I can recall of algorithms with simple changes to the algorithm that greatly increase/decrease the precision of the result (due to error cancellation/accumulation) require a fair bit more explanation.

Point of order: Intel FPUs have 80-bit internal results; most others don't. Code that assumes 80-bit internal results is highly non-portable (even between different x64 compilers, if there are spills to variables declared as 64-bit double).

I mean, in theory, hardware could contain a 2048-bit fixed-point accumulator and sum numbers in any order precisely before rounding, but the existence of that hardware doesn't affect algorithms that can't assume it :)

So what does "stable" mean in this context? For me, a *usefully* stable algorithm would have the property of similar errors (relative to epsilon*result) regardless of internal precision. If you drop down to 32-bit floats, or even 16- or 8-bit floats, does it still work? Algorithms that are fine *given enough internal precision* don't fall into that class of stable.

breaks down in several applications:

* Cryptography and similar places, where hundreds if not thousands of bits are required.

* Financial applications, where exact decimal fractions are required. A common workaround if using IEEE floats is to scale values by some power of ten, as the number of decimal places is frequently bounded.

IEEE nowadays defines decimal floating point systems, though I'm aware of any hardware that implements decimal floats.

At least, POWER7 has decimal floating point hardware support, and I think POWER6 did too. (The IBM 360 had decimal fixed point, so it might be something lingering in the corporate culture.) It may not be coming to a desktop near you any time soon, but the software emulations seem to be acceptably performant and the language integration means that hardware support could be transparently added.

I'm not sure why you would expect floating point representations of any kind to be helpful in cryptography, where the multikilobit numbers are generally integers, and would not, as far as I know, benefit in any way from having an attached exponent. There is no shortage of languages which provide reasonably usable exact bignum datatypes, and at least a couple of reasonably used C libraries.

In short, there has not been a revolution in numeric representation, but there are ongoing incremental improvements. Despite these improvements, nothing has really changed the issue that:

1) Numerical analysis is hard, and the results are sometimes surprising and often counter-intuitive.

nor the resulting cultural feature:

2) Many programmers instinctively distrust floating point computations, at least in part because the preceding makes it seem like black magic.

in number-theoretic crypto algorithms, where you are using very large integers FTMP. My post was in response to a "doubles are good enough for most things" comment above--there are things where 64-bits of precision are not enough.

Most financial applications are not computationally-intensive, particularly in comparison to most forms of scientific computing, so software implementation of decimal numbers is probably good enough for most business applications. (I'm sure someone will posit a counter-example, of course.)

Agreed on the general observations about numerical analysis. One interesting thing I've long noted about PLT culture is a bias, in some quarters, against the various shortcuts taken by hardware numerics--in particular, with the fixed-width ints (with modulo arithmetic) found in many mainstream languages, often being scorned.

My preference is for a good numeric tower, with many different types of representations available for what application you like--hardware ints and floats/doubles (particularly IEEE); software bignums, rationals, fixed-point numbers, and decimals, complex numbers and quaternions, symbolic representations of numbers, suspended computations, etc. Getting this right seems to be a difficult problem...

One interesting thing I've long noted about PLT culture is a bias, in some quarters, against the various shortcuts taken by hardware numerics--in particular, with the fixed-width ints (with modulo arithmetic) found in many mainstream languages, often being scorned.

I guess that's directed at me, because I'm more-or-less implementing something else and that's what I started the thread about.

The hardware fixed-width integer (or float) has the same relationship to numbers that the fixed-width byte array has to strings. At some level a more useful implementation of strings has to be implemented in terms of them, or in terms of dynamically allocated structures containing them, but if fixed-width byte arrays are what a high-level language generally uses for strings, then the language is deficient, and its strings won't behave in the ways people want strings to behave.

In the same way, numbers must be implemented in terms of what the hardware provides, but what the hardware provides has limitations that I do not want to tie a higher-level language to. Also, I do not want the language's semantics to depend on a particular hardware, because new hardware shouldn't change the semantics or behavior of existing programs.

I'm implementing "numbers" that are a product of rational, decimal-exponent, and binary-exponent values because those are the kinds of numbers people use a lot of.

I'm of the opinion that the "normal" numerics that fifth-grade math students spend time learning, with fractions and decimal exponents and answers that are exactly right, require far fewer axioms to describe and reason about, and are therefore more appropriate in a high-level language and drastically easier to prove things about.

Bignums are even simpler than my fixed-width representation, of course, but I'm making them a separate type because they contain a trap for the unwary; bignums can easily balloon to consume all memory (and all compute power) in an iterative calculation, and I don't want it to be hard to avoid that particular bad behavior.

So my numerics are an attempt to provide something useful above the level of hardware-dependent semantics and below the level of burning all resources available for standard iterative refinement algorithms.

I want to provide some guarantees, like knowing when the answer is an approximation, like returning a value that doesn't require additional rules (two's complement integers, normalized and denormalized floating-point representation, rounding errors, and modulo arithmetic) to explain, and like keeping sources of error as few and as small as I can.

Fixed-width hardware FP as provided by IEEE754 is better than what has gone before in terms of hardware floating-point numerics -- but it is still quite bizarre as compared to the way math works in grade-school math classes.

For example, in grade-school math classes, when you add 3/4 to 3/21, you get 75/84, not some string of digits behind a decimal point that don't after all exactly add up to 75/84 and which, in memory, are represented by a set of binary charges and a complicated set of rules about them which in fact says they don't after all add up to exactly the decimal fraction represented by that string of digits behind a decimal point.

Imagine doing proofs about programs that use the hardware-provided version of fractional numbers, and think about the set of axioms involved, and it's bewilderingly complex. It becomes hard to tell whether the proofs you do about them are in fact correct. So I'm trying to create numbers that behave the way numbers are expected to behave by normal grade-school math rules, until an approximation *has* to be made, and then *warn* that the result is an approximation rather than the right answer.

Ray

In cryptography the type of primitive determines the number of bits. For public key crypto it is common to use modulo exponentiation as the underlying operation so a large enough group size to guarantee security does require thousands of bits. But this is not the only choice of operation to implement that primitive, or the only interesting primitive. Elliptic curves use a different operation and the group size is normally just a few hundred bits. Symmetric encryption and hashes also operate on much smaller values.

Even for large bit-sizes there are ways of implementing multi-word arithmetic. Residue-number systems are one method of encoding, although they prove to be somewhat inefficient. Dekker style splits use the IEEE754 error guarantees to use a sequence of floating point values as adjacent words in a multi-word value. The exponents of each value in the sequence decrease by the word size, each word uses a half-ULP error bound and the subsequent word is the error term to correct below that precision.

Although Dekker-style splits were invented for arithmetic over rationals they work just as well for embedding integers into sets of FP values. Recently some students used them to implement Blum-Blum-Shub on a GPU. They have been used in RSA implementation on a GPU previously.

Edit: Now that I look more closely at the indentation, I think that it may be the case that what Andrew Moss believes is "not so true" was Scott Johnson's comment, and not mine. Unfortunately, I can't re-indent this comment, so I'll just leave it where it is. To Andrew: if we're in agreement -- and quite possibly we are -- that's great.

You can use FP primitives to implement multiword integer arithmetic, and that is a good strategy on hardware where FP units are wider, faster, more parallel, or otherwise better than integer units. But that's an internal implementation detail, invisible to the program which uses multiword integer arithmetic.

So there's no reason for anyone to expect an FP implementation to provide multiword integer arithmetic, and I continue to find it odd to complain that IEEE754 "breaks down" in the cryptography application. In fact, it's the reverse, as Andrew Moss says: the IEEE754 error guarantees *allow* the implementer of a multiword integer arithmetic library to make use of FP if it is a reasonable alternative.

I believe it is an indentation thing that makes it look the other way around, we are in agreement :)

You mention that you want to see better support for error models and I agree with you. But the best notion of error is often a bound on the difference from the "exact number" intended. So I'd love to be able to specify the exact answer that I intend using limits and other not-naively computable expressions and then have the have the IDE help me find an efficient approximation. (This is different from the notion of an "exact number" flag at runtime, which may be all you were disagreeing with)

1. Approximation for the sake of efficiency, e.g. floating point numbers.
2. Uncertainty of inputs, e.g. inexact measurement.

Solutions for dealing with these two will probably be different, since in the second case you want to track the inexactness, but in the first case tracking that would defeat the efficiency. Probabilistic programming may be a good way to deal with #2 (much superior to interval arithmetic!)

I think we can track inexactness without defeating efficiency. One way is to model it as a parallel computation that represents how errors will accumulate, which could then be calculated (lazily or in parallel) if desired (or dropped as dead code if not desired).

Another technique is to track error and uncertainty models in the type-system, which might prove a very useful approach for controlling the amount of error. As I understand it, floats often cause trouble for losing information; having a static understanding of what information might be lost is quite useful.

As somebody who has dabbled in numerical algorithms I can tell you that it's not that easy. The actual accuracy of a floating point operations cannot easily be modeled statically or even dynamically independently of the real computation as it highly depends on the actual values involved. For example 100.0001 - 100.0000 may be highly inaccurate, while 1.0001 - 1.0000 may be very accurate (depending on the number of bits of precision used).

The only practical way I know is to run the computation with higher precision and seeing if you still got the same result. Languages can certainly help by easily supporting that, and tooling to show a side by side execution of the different precisions to determine at which point the lower precision begins to diverge from the higher precision. Efficiency suffers tremendously once you go beyond hardware precision so that's only practical during development.

A language that addressed this in its type system might force you to write alternative algorithms that are easier to reason about, or at least be explicit about accepting the risk where problems might occur.

If the type system was sufficiently expressive, you might be able to reason about relative sizes of values even without knowing the exact values.

I think managing error in numerical algorithms would be a lot less difficult with better tools.

A formal approach seems highly unlikely to work. Even proving numerical accuracy or reasoning about it by hand is impossibly hard in most cases, let alone trying to do it formally in a type system. In fact in most cases it's not just hard, but actually impossible since all but the most trivial numerical algorithms DO NOT work for all inputs, and it's not actually known on which inputs they are accurate to which degree, since floating point semantics makes this impossibly hard to determine. Saying lets solve this with a type system is like saying "proving Fermats last theorem is hard, lets do a formal proof in Coq". It's the "???" in "1. Lets formalize! 2. ??? 3. Profit!" that ruins the plan.

As I said above, you can make assumptions about inputs more explicit. It isn't clear to me why you think reasoning by hand would be easier. Formal methods are often quite useful in cases where by hand effort has inconsistent success.

Those ??? problems are often well worth tackling. It just means there is no obvious or well known solution. The only impossibility I'm sure of is improving a situation without changing it.

Because formal proof is virtually always harder than informal proof.

The reason you find numerical analysis "hard" isn't because you can't perform an informal proof; it's because your informal proofs are inadequate, and you're bad at doing more formal proofs by hand.

If we had to prove type safety or security of large programs while "reasoning about it by hand", I think we'd find that quite difficult. As far as I can tell, this is a similar situation.

That wouldn't be hard in the same sense that numerical analysis is hard. It would just be a lot of work due to the size of the programs. In contrast, numerical analysis on very small programs is already hard. Even numerical analysis on code that does not have any loops and does not have any conditionals is hard.

Meh. I disagree with your assumption that 'size of the program' is the only thing that makes hand reasoning about safety and security 'hard'. Even safety for code that does not have any loops and does not have any conditionals can be hard. If you remove access to the ambient program (and thus assumptions about values) it can be even harder. The same is true for numerical analysis, except you don't currently have even an effective mental model (ontology, types, framework... abstractions) to discipline your code and support reasoning about potential numerical errors.

Sure, reasoning about security may be hard for small programs in some cases, but in those cases type systems or other forms of automated reasoning would certainly not help much. The kind of analysis that type systems do can be quite trivially be done by a human with a piece of paper. It's just a lot of work.

It would really help if you put some weight behind your assertions by providing concrete ideas on how types or anything else make analysis of floating point behavior feasible.

The kind of analysis that type systems do can be quite trivially be done by a human with a piece of paper. It's just a lot of work.

Assuming the type system is already understood by the human, the human can structure code to support reasoning with those types. This makes it easy for a human to later perform "the kind of analysis that type systems do", though it can still be difficult with for the more expressive type systems. But if you take away the human's initial knowledge of type and structure, the resulting code will not be organized in a manner that the type system can handle. In these circumstances, the kind of analysis that type systems do will become very hard, and often the only option will be ad-hoc abstract interpretation.

You argue that numerical analysis is very hard. And I think you are right to believe so. You also claim that the cause for this difficulty is independent of the tooling, i.e. essential to the problem. But you have not justified your pessimism.

type systems or other forms of automated reasoning would certainly not help much

Even for small problems, dependent types, ghost types, hoare logic, first-class abstract types (or sealers/unsealers), and other frameworks for structure and automated reasoning can help a great deal for safety and security. I do not share your certainty.

put some weight behind your assertions by providing concrete ideas on how types or anything else make analysis of floating point behavior feasible

I don't have solid ideas on how to approach numerical analysis specifically, just a strong feeling that it is approachable. I do have intuitions on how to approach it, some of which I have mentioned:

• Types that refine absolute ranges for values
• Dependent types that track absolute error for values
• Types that can track relative values of numbers (e.g. x is 20% to 80% of y)
• Types that track whether a value was the result of an unregulated subtraction, and thus whether it may be dangerous to use in a later division or multiplication

You also claim that the cause for this difficulty is independent of the tooling, i.e. essential to the problem. But you have not justified your pessimism.

I think I have? Nobody knows for certain how accurate any but the most trivial numerical algorithms are. There is some knowledge about the stability of certain algorithms (e.g. the condition number of a matrix is related to how accurately you can invert it), but this is very much a simplification of actual floating point behavior. Inverting a matrix is just a couple of lines of code. What makes you think that a type system would suddenly provide tremendous insight into the behavior of those couple of lines of code that is impossible for humans to get? Or take it even simpler: summing a list of numbers. The resulting error due to floating point is very much dependent on the actual numbers involved, and in which order they appear in the list. What would the type of such a summation function be in your type system? I'd expect you can't do anything significantly more useful than "the error is what the error is", i.e. "the error is floatsum(xs) - infiniteprecisionsum(xs)".

Even for small problems, dependent types, ghost types, hoare logic, first-class abstract types (or sealers/unsealers), and other frameworks for structure and automated reasoning can help a great deal for safety and security. I do not share your certainty.

Well I certainly agree if you ignore the first half of that sentence and take the last half on its own. If you take the complete sentence I'd disagree ;) If a program is small but hard to analyze that is going to be because it has tricky behavior, like a Collatz sequence program. But in those cases a type system won't help you analyze it. Type systems and other automated checking facilities are great for doing huge quantities of trivial reasoning. This is even the case with dependent types. The theorem you're proving there might be non-trivial, but constructing that proof requires ingenuity on the part of the human. The automated checking of the proof is still a trivial process. Types don't (yet) give you magic proving powers.

A matrix could be validated as an invertible matrix (with constrained relative error) before you even try to invert it. This could be based on actual number values in the matrix. Similarly, for summing a list of numbers, we could have a type system constrain the actual numbers involved, and potentially even the relative order they appear in the list.

You seem to be forgetting that type systems will constrain how you write code, and provide answers based on those constraints. When you ask "how will numerical type systems help with FOO, which I wrote without using the numerical type system", you're doing it wrong. Of course a type system won't help you unless you use it while developing and designing your code.

If a program is small but hard to analyze that is going to be because it has tricky behavior

Yes. But behaviors aren't tricky all for the same reasons. There are many cases where type systems can force developers to recognize, reify, and address potential concerns.

The theorem you're proving there might be non-trivial, but constructing that proof requires ingenuity on the part of the human.

I have not suggested that type systems eliminate need for 'ingenuity on the part of the human'. Indeed, it often some ingenuity to get clever or tricky concepts approved by the type system.

The automated checking of the proof is still a trivial process.

It's trivial when automated. By hand, those automated proofs are often non-trivial. Even for programs whose runtime aspect is very small.

Types don't (yet) give you magic proving powers.

Type systems support proofs by two mechanisms. First, automation reduces errors and improves maintenance and confidence far beyond anything you could achieve by hand. Second, they force you to write code that is more easily provable (even if you subtracted the type system after writing it). It isn't magic. But it is advanced technology.

You seem to be forgetting that type systems will constrain how you write code, and provide answers based on those constraints. When you ask "how will numerical type systems help with FOO, which I wrote without using the numerical type system", you're doing it wrong. Of course a type system won't help you unless you use it while developing and designing your code.

Well yes, but that's moving the goalposts. Sometimes you are just given a list of numbers, and you want to know their sum. Having a function that only accepts certain kinds of lists doesn't help me. Similarly, sometimes you want to invert a matrix.

But even when you allow yourself to move the goalposts, it's still very hard. For example can you give a concrete constraint that you'd place on a matrix that would let you say something useful about the floating point error after inverting it? Of course you could say "it must be the identity matrix", but keep in mind that your matrix inversion function should still have *some* usefulness.

It's trivial when automated. By hand, those automated proofs are often non-trivial. Even for programs whose runtime aspect is very small.

Do you really not understand what I meant? Surely you must see that doing type checking by hand is "non trivial" in a completely different way than proving interesting theorems. Somehow our discussions always end up being about definitions, and I get the feeling that you're intentionally trying to find definitions or take things out of context so that you can interpret what I write as a false statement :( It would be much more interesting to discuss the things, instead of the names of the things.

Type systems support proofs by two mechanisms. First, automation reduces errors and improves maintenance and confidence far beyond anything you could achieve by hand. Second, they force you to write code that is more easily provable (even if you subtracted the type system after writing it). It isn't magic. But it is advanced technology.

I agree with this, but this doesn't address the point.

A typed language will not prevent you from summing arbitrary lists or inverting arbitrary matrices. But it also won't promise you good data from your results. Rejecting safety for convenience is an option that must be addressed based on context.

I do not believe I have moved any goalposts. But I might have different goalposts than you do from the start. I believe that "improvement without change" is a contradiction, and I quite expect that addressing numerical safety (and addressing how code behaves in the presence of error-models in general) might require changes in how we code.

can you give a concrete constraint that you'd place on a matrix that would let you say something useful about the floating point error after inverting it?

A practical constraint is that the matrix is normalized. This allows developing inversion code that corrects floating point errors. And we can extract information to recover from normalization. But there may be weaker constraints that can achieve the same goals (e.g. allowing developers to track epsilon error) without the efficiency overheads.

(The problem is then shifted to normalizing the matrix.)

doing type checking by hand is "non trivial" in a completely different way than proving interesting theorems

I think they are not so completely different. Sure, we can separate the frustrations of 'finding a proof' from those of 'validating a proof' after your code has been written, especially if it's still fresh in memory. But I think you can't easily separate them during development. Or, a few years later, during maintenance.

With more expressive type systems, and especially on code you don't already know, I think you will find the hand-proof is not merely rote where you'll make only mechanical errors.

and I quite expect that addressing numerical safety (and addressing how code behaves in the presence of error-models in general) might require changes in how we code.

So are you saying we should completely abandon all numerical algorithms, since pretty much none of them can be analyzed with respect to floating point behavior? You say that the type system will only work for certain kinds of code. What makes you think the intersection of that with the practical problems we want to solve is nonempty? It would help if you'd give at least 1 example of something that is both useful and a human armed with your type system can practically analyze to prove that it is accurate to some degree.

A practical constraint is that the matrix is normalized.

I don't know what you mean by normalized. From a Google search is is apparent that there is no accepted definition of normalized matrix. Can you elaborate?

With more expressive type systems, and especially on code you don't already know, I think you will find the hand-proof is not merely rote where you'll make only mechanical errors.

No, what I said was that CHECKING the types is rote, not coming up with the proof (for example coming up with a proof that a certain algorithm is accurate to degree X is HARD, but once you have formalized such a proof in some type system, then checking those types is rote). That's exactly my point: type systems automate the trivial stuff, and don't (yet) significantly help prove non trivial stuff (and they certainly don't help much with coming up with new theorems of the kind "matrix inversion algorithm X is accurate to degree Y" -- but it seems that we do agree on this point since you insist on coming up with new algorithms that are supposedly easier to analyze).

Numerical programmers aren't used to needing to justify their code for floating point performance. Best effort is often good enough, and they can attack the few cases where it empirically proves otherwise.

I do not suggest that we abandon existing programs. I'd favor an incremental migration strategy, where existing programs are typed in a manner that indicates they make no useful guarantees about how errors grow or propagate. The type system in this case is not be about limiting how code is composed (except for first-class functions). Rather, it's more about understanding error.

Unless the type system is stupidly inexpressive, I think there is not much risk of it failing to address problems. Numerical programmers can often provide assumptions/conditions/contracts/data bounds for which they can justify their algorithms are good, and we can reify these conditions and their justifications in the code.

But there are social issues with a type model. The system will often hinder expressing the most efficient algorithm with any useful guarantees. Programmers will sometimes grow frustrated being unable to achieve the guarantees on errors they want (due to their newfound awareness of error).

I don't know what you mean by normalized.

I understand them to be a standard concept. (cf. Normal matrix)

No, what I said was that CHECKING the types is rote

I was disagreeing with you. Hand-proof, when CHECKING the types of more expressive type systems, and especially on code you don't already know, is not merely rote.

To say that "type systems automate the trivial stuff" seems quite circular. I think you only call it trivial because type systems automate it.

Unless the type system is stupidly inexpressive, I think there is not much risk of it failing to address problems. Numerical programmers can often provide assumptions/conditions/contracts/data bounds for which they can justify their algorithms are good, and we can reify these conditions and their justifications in the code.

Again, this is a false assumption. Numerical programmers in almost all cases CANNOT provide assumptions/conditions/etc for which they can justify that their algorithms work correctly with floating point. This is the whole point!

I understand them to be a standard concept. (cf. Normal matrix)

Well normal != normalized. I don't see how you can guarantee anything about numerical matrix inversion by restricting to normal matrices, since a normal matrix isn't necessarily invertible in the first place, even with infinite precision! I can help you a little: check out the wikipedia page on condition number. Now that does allow you to say something about the accuracy of matrix inversion, but only in a very much idealized model. In particular it does NOT model floating point arithmetic.

I was disagreeing with you. Hand-proof, when CHECKING the types of more expressive type systems, and especially on code you don't already know, is not merely rote.

Meh, again uselessly arguing over definitions like "rote" by taking things out of context. The point is that type systems do not magically let you come up with new theorems regarding the degree of accuracy of numerical algorithms under floating point arithmetic (whether you are allowed to change those algorithms or not), something which large numbers of numerical analysis researches have not succeeded in.

Numerical programmers in almost all cases CANNOT provide assumptions/conditions/etc for which they can justify that their algorithms work correctly with floating point

Are you trying to argue that this is an acceptable situation? Shouldn't programmers be able to justify their code?

I think your claim doesn't hold in any case. At least in the fields I've worked in that use a lot of numerical methods (e.g. computer vision), developers can understand error introduced by their code with an understanding of ranges on input, limits on matrix size, avoiding division, etc..

I don't see how you can guarantee anything about numerical matrix inversion by restricting to normal matrices

Maybe so. I didn't think too hard on it. My intuition is (1) we can constrain error introduced during the process of inversion that might otherwise arise from having very high or very low numbers in the matrix, (2) Since the inversion of a normal matrix is also a normal matrix, we can use this to help recover errors introduced (perhaps through an iterative algorithm).

type systems do not magically let you come up with new theorems

Indeed. There is no magic involved. Just pervasive tools and an encouraging environment.

It is clear that you believe numerical analysis with floating point arithmetic is a case that will defy all attempts to tame it. I think there is much we can do to improve our ability to track, control, and reason about errors - not just from floating point representations, but error models for all measured values. I doubt I can convince you, and I know none of your arguments have convinced me. So I think it is time to stop arguing the point.

something which large numbers of numerical analysis researches have not succeeded in

Historically, having large numbers of people working on a problem and not succeeding is a common precedent to success.

Shouldn't programmers be able to justify their code?

Is it an acceptable situation that mathematicians haven't been able to prove the Collatz conjecture? Ideally they already would have. Meanwhile the world works just fine, it's not high priority. Is it a problem that you can't prove handwriting recognition algorithms correct (since it wouldn't even be clear what the theorem would be, much less the proof) if they work in practice?

I think your claim doesn't hold in any case. [...] developers can understand error introduced by their code with an understanding of ranges on input, limits on matrix size, avoiding division, etc..

If you think my claim does not hold please provide counterexamples. In practice programmers do just fine by heuristic methods and testing, but that's something different than formal proof.

rror introduced during the process of inversion that might otherwise arise from having very high or very low numbers in the matrix

A matrix being normal has nothing to do with the sizes of the elements in it. You are probably being (understandably) confused between one of the many meanings of normal & normalized.

Since the inversion of a normal matrix is also a normal matrix, we can use this to help recover errors introduced

I don't follow this implication. Everything is easy in the abstract, but giving a concrete method is hard.

I doubt I can convince you, and I know none of your arguments have convinced me. So I think it is time to stop arguing the point.

You could easily convince me that floating point analysis is practically tractable with a single example, but I agree that in the absence of that further discussion is pointless.

Historically, having large numbers of people working on a problem and not succeeding is a common precedent to success.

Survivorship bias. It's an even more common precedent to failure, especially in this case since pretty much nobody is even trying to do such a thing since it is intractable, and instead use heuristic methods like analyzing stability and doing testing that mostly work in practice but do not give formal proof that the result is accurate.

If a numerical programmer wants to make a guarantee they can't formally prove, they can just put it in the documentation. If they need to make a formal assumption that others have not managed to formally prove, they can take it as a parameter, and users of their code can work on figuring out how to pass it in. In the latter approach, their code may not be practically useful for a long time, so they should probably find a way not to rely on that assumption after all. Perhaps they could weaken their guarantees or reduce their assumption to the minimum they actually need.

If the domain the programmer cares about is extremely difficult to formalize, as you claim for floating point numerical algorithms, then the limit outcome will be that it isn't really formalized at all. The algorithms will still be possible to express using basic types like Float, but most of these algorithms' benefits will have to be taken for granted based on the documentation, unit tests, load tests, or other assurance mechanisms. Even in this pessimistic outcome, we're no worse off than usual!

Personally, I think if floating point analysis is so hard to study formally, that's a pretty good reason to work with simpler representations. How about a format xxxxyyyy which is treated semantically as e^xxxx.yyyy? This way the distribution of representable values is about the same as floating point. We lose the ability to represent exact integers, but we gain a cheap multiplication operator that's isomorphic to modulo addition. Looks like this already exists under the name logarithmic fixed-point.

Also known as logarithmic number system. I agree with Jules Jacobs that formal guarantees for numerical calculation are not very plausible. Nevertheless, it may still be possible to provide average-case guarantees, and rely on the law of large numbers to ensure that total error stays reasonable. In this case, formal methods may help ensure that our numerical algorithms are reasonably well-behaved, and avoid known sources of error accumulation.

If you're replying to me, I've considered them but they're not what I want, I think.

I have considered constructive reals. In some trivial way, I have them already in the way that any language that supports lazy evaluation has them. But to be useful, there has to be support for a final step whose semantics are something like, "now calculate an actual value for this promise, exactly or to within some known margin of error." And it's not easy to build that step.

If you could tell from an expression what accuracy you needed in subexpressions in order to get an answer within some margin of error, you could backward-chain through the expression to find out what margins of error you needed at each calculation, then forward-chain knowing what margins of error you needed to reach in each calculation.

But it's really hard to tell that from an expression given a standard floating-point form (or even an alternate form like the one I've proposed), without first having an approximation of the value so you can know things like the order of magnitude of each operand. So before you can backward-chain to establish required margins of error, you have to solve it approximately once to (hopefully!) establish the binary magnitude of each operand. For efficiency, you want to memoize this just in case you can use the results of the preliminary calculation in the final calculation. For accuracy, you want to memoize it so you can check against the final computation and verify that your estimated orders of magnitude were correct.

And theoretical maximum error in rounding tends to accumulate in computation a lot faster than actual error does, so the final calculation is still likely to be carried out to far more precision than the particular instance of defense against accumulated error requires.

The alternate representation form I've described above eases this pain somewhat by allowing more calculations to be done exactly; taking the maximum rounding error to a smaller power puts constraints on accumulated error that will save a lot of effort in large calculations. But it still has the problem of maximum rounding error taken to an exponent by long computations being MUCH greater than average rounding error taken to the same exponent.

Anyway, I haven't implemented that yet. The above is just a discussion of the problem space.

Ray