This example shows a hazard of floating-point, but it does not show that the FP error is worse than modeling error.

Here, the problem arises because the equational system is ill-conditioned. Yet you present the system as if the coefficients are known exactly.

In reality, this ill-conditioned model can arise for two reasons:

+ The equations represent a physical system that is ill-conditioned -- a slight change in the system can produce a huge change in its output. So where did those coefficients, with their 8 digits, come from? A minor error measuring the coefficients will throw off your "correct" solution. In this case, modeling error is every bit as important as floating-point error.

+ The physical system is not poorly conditioned, but the model is. In this case, by definition, modeling error overwhelms floating-point error.

So the essay still advances the correct perspective -- if you are modeling an ill-conditioned system, you had better spend a lot of time worrying about the effects of modeling errors, before you get hung up worrying about FP errors. And if you know the system is ill-conditioned, then poor conditioning in the model is one of the things that an experienced analyst will be prepared for.

> In reality, this ill-conditioned model can arise for two reasons

> + The equations represent a physical system that is ill-conditioned

> + The physical system is not poorly conditioned, but the model is

In my experience, there are other ways this kind of problem can come up where it would have nothing to do with the model or the system being ill-conditioned. I'll give one example below, but I expect there are more I'm not thinking of.

One that I recall encountering is, say you have a model, and you compute the time 't' at which some event occurs (say, the time at which a condition x < h is violated) in some part of your model. Once you find this 't', you plug it back into your full model and evaluate the model at that time. Now you try to compute the time of the next event and repeat. But there's no guarantee that x = h (or even x >= h) at this time 't' here, even though that's the condition you solved for... and if x < h again, you can get stuck here re-evaluating the same event forever. And notice this had nothing to do with the system being ill-conditioned; it was purely due to FP round-off making x just a tiny bit smaller than h.

Now, there are two ways to fix this example. The usual hack (once you notice the implicit x = h assumption, which only obvious because I just told it to you, not when you're actually coding) is to just set a tolerance parameter for |x - h| being small... but introducing tolerances doesn't always solve the problem and in addition can introduce more messes of its own, since tolerances become quite difficult to reason about and deal with in the presence of feedback (they can cancel or compound themselves, for example) and you don't want to have to keep figuring out reasonable tolerances for every condition. On the other hand, there's another (nicer) solution: simply force x = h here, because that's literally the condition you solved for. The trouble is, this is a very unnatural approach, and hence very, very far from obvious when you're coding it for the first time. You have to introduce a mechanism for forcing the state of one part of your model externally, while the rest is computed internally. And you have to explicitly invoke this mechanism. It's something that makes sense after the fact, but it's very difficult to foresee because it makes no sense when you're thinking about the problem purely mathematically (i.e. in \mathbb{R}), and yet it can bite you just when you think your code that switches to finite-precision is fine (as it has bitten me).

My main point is, I wouldn't assume you can predict when these sorts of problems come up. They're not always due to poor conditioning or due to poor modeling, or poor coding for that matter. They can be products of reasonable, competent work and still crop up when you really don't anticipate them, no matter how "nicely-behaved" the inputs to your solver are.

> And notice this had nothing to do with the system being ill-conditioned; it was purely due to FP round-off making x just a tiny bit smaller than h.

This doesn't sound right: the comparison function x(t) >= h is extremely ill-conditioned at the point t where x(t)=h, having condition number Inf.

That is not the system being ill-conditioned. The system would be well-conditioned. That's just a comparison that is assumed to hold true inside the simulator because it's trying to advance its time to that of the next event. It's internal to the simulation algorithm.

I think it is: your code ends up taking different paths based on the result of a comparison, which is ill-conditioned. The underlying system might be well-conditioned, but the question you are asking about it is not.

I'm not disputing whether there is an ill-conditioning present here; I am saying that it is not the system that is ill-conditioned or otherwise somehow particularly pathological, and that this is obviously not due to a modeling error. Please consider re-reading my comments since my point is something other than what you are disputing -- you seem to be disputing a point that I was never making.

It is not the system that is ill-conditioned or pathological, nor is it modelling error, I agree. It is the comparison x(t) >= h that is ill-conditioned, and you'd be wrong to rely on its result (much like catastrophic cancellation, it is a bad idea to compare approximately equal numbers with errors in them, no matter where the errors come from—the errors could come from measurements, and the same logic would apply). We might have different ideas about what in your example is specified as given (with full precision) and what is specified approximately? But I do not think I misunderstood anything you wrote, please consider re-reading my comments too.

If you try to solve a linear system with a condition number of 3.73e16, it's not really correct to say that the problem lies with floating-point errors, because the real problem is that you tried to solve an extremely ill-conditioned problem. The wrong answer you get is backward stable, in fact—it is the exact solution of a nearby perturbed system. One way to think of it is that the phrase "the answer you should get" is somewhat ambiguous. I don't think your example goes against the post's main points.

> the real problem is that you tried to solve an extremely ill-conditioned problem

This makes it sound like you're saying just not to solve problems when they're too annoying, but I somehow don't think that's what you're trying to say.

In a way that's technically true, because solving ill-conditioned problems with noisy data is both difficult and the results are often meaningless. But you almost always have some choice about which of a multitude of subproblems you try and solve to reach your actual goal. It's like the difference between executing an algorithm and designing an algorithm. So instead of solving ill-conditioned subproblems, you spend your effort getting to well-conditioned subproblems instead.

I think a common example of an ill-conditioned problem that's too annoying to solve is weather prediction past 2 weeks in advance, because it's a chaotic system, so I don't think what you're saying is wrong.

That's only a place you would "never expect" because you've carefully hidden the gotcha. Stated less obtusely, you're finding the intersection between two lines which are parallel to JUST above the precision of an IEEE double. This is a cute example, not a real world situation. In the real world, precision loss in this circumstance is acceptable and any real system that accepted a "parallel" condition like this would have a check for it.

For a similar situation: atan2() exists because of a coverage issue in the problem area and not a bug in the idea of a tangent or a float.

I would think that atan2() exists because atan() doesn't invert tan() correctly for some purposes -- if you're trying to find the angle associated with a point then atan(y/x) is just the wrong answer.

The fact that it avoids the division by 0 is a convenient side effect of also being more correct.

(I'm just guessing; I don't know the actual history of it.)

The point was more that problems that need to use atan() in a context where the line may point near the Y axis don't flip out over "floating point precision", they just use atan2().

Problems like the one in the grandparent that need to intersect nearly-parallel lines (or subspaces) use special techniques better suited to the problem and don't just feed it to a linear solver.

How does that contradict what the parent was saying? Parent said problems like this can come up where you least expect it, and you rebutted that "this is not a real world situation" because "any real system that accepted a 'parallel' condition like this would have to check for it"? If anything you're confirming the point, not refuting it.

This is a carefully constructed (malformed) problem presented in an obfuscated way. It didn't "just happen" that those ~20 bit numbers happen to produce slopes that differ only after ~57 bits, come on. If you're going to do that kind of trickery, you can come up with an equivalent that blows up any fixed precision limit and/or explodes the representation size for any arbitrary precision computation system.

It's a parlor trick, not a serious problem.

in MATLAB:

>> A = [0.25510582 0.52746197; 0.80143857 1.65707065];

>> B = [0.79981812; 2.51270273];

>> X = A\B

X =

   0.272534072955269

   1.384541698505404

>> A * X

ans =

   0.799818120000000

   2.512702730000000

>> A * [-1;2]

ans =

   0.799818120000000

   2.512702730000000

Yeah, floating point error is the least of your worries... until it's not. But usually, it is.

Wolfram has no problems with it.

Other commentators have already shown that this is easily solvable with systems like MATLAB, Octave, or Wolfram. I will add, though, that I believe that the potential pitfalls of the above are still more due to modeling error than floating point.

Even if we don't relate the above system to a physical problem, the issue is that asking for a solution to `Ax=b` is generally an ill-defined problem. We always have to ask three questions:

1. Does a solution exist?

2. Does a unique solution exist?

3. If a solution exists, but it's not unique, then which solution do we want?

For a system `Ax=b`, I contend the first one is generally the hardest to answer. One way to view (1) is to ask whether there exists an `x` such that `|| Ax-b || < eps`. Here, we have a floating point question, but I don't think it's a floating point error issue. Simply, since we're working with floating point numbers, we have limited accuracy and have to define what we believe zero to be. Depending on the scaling, `eps` could be really big or really small and knowing the correct scaling is a modeling issue. Now, if we don't believe that there exists an `x` such that `|| Ax -b ||`, then we can solve something like a minimum norm problem and let `x=argmin || Ax - b ||` or we could do something else. Again, this is a modeling issue and it depends what we want.

Now, if we believe a solution exists and is unique, then we could just factorize the system with something like an LU. However, even here, we have a modeling issue because we have to decide what we mean by unique in floating point arithmetic. Again, this is not a floating point error issue, it's really an issue of what we want to define zero to be in floating point arithmetic. One way to define uniqueness is to look at the Hessian of the least-squares objective, `J(x) = 0.5 || Ax -b ||^2`, which is just `A'A`. If the smallest eigenvalue is smaller than `eps`, then we're rank deficient and the optimization problem `argmin J(x)` does not have a unique solution. Now, this relates to the condition number. Since `A'A` is symmetric, `cond(A)^2=cond(A'A)=|lambda_max(A'A)|/|lambda_min(A'A)|`. Basically, we have the scaling issue again. If the largest eigenvalue of `A'A` is much larger than the smallest, then we should probably just call the smallest eigenvalue effectively zero because we're limited by floating point arithmetic in our ability to correctly model differences.

These leads us to (3). If we have a solution, but it's not unique, then we need to pick what solution we want. Generally, this becomes the minimum norm solution or `x=argmin{ ||x|| : Ax = b }`. There's a lot of ways to calculate this, but using a QR factorization is my favorite. For the parent problem, this solution is `[0.594352357167031;1.228894994185025]`, which is precisely what MATLAB/Octave gives for `A\b`. The solution `[-1;2]` satisfies `Ax=b`, but it's not minimal norm. I guess we could call it the aesthetically pleasing solution.

Anyway, again, all of this relates to floating point, but I contend it's not a floating point error issue, but more of an issue of how do we model zero in floating point arithmetic. This relates to the scaling of the problem, but once we settle on that, then we can run through questions 1-3 to find the solution we want.