[quanto]

> when I can get away with it I'd prefer to prove things with real numbers and assume magically they transfer to floating point.

True for some approaches, but numerical analysis does account for machine epsilon and truncation errors.

I am aware that Inria works with Coq as your link shows. However, the link itself does not answer my question. As a concrete example, how would you prove an implementation of a Kalman filter is correct?

[thesmtsolver]

There is nothing inherently difficult about practical implementations of continuous numbers for automated reasoning compared to more discrete mathematical structures. They are handleable by standard FOL itself.

See ACL2's support for floating point arithmetic.

https://www.cs.utexas.edu/~moore/publications/double-float.p...

SMT solvers also support real number theories:

https://shemesh.larc.nasa.gov/fm/papers/nfm2019-draft.pdf

Z3 also supports real theories:

https://smt-lib.org/theories-Reals.shtml

[kragen]

I'm curious about how you'd do that, too. I haven't tried doing anything like that, but I'd think you'd start by trying to formalize the usual optimality proof for the Kalman filter, transfer it to actual program logic on the assumption that the program manipulates real numbers, try to get the proof to work on floating-point numbers, and finally extract the program from the proof to run it.

https://youtu.be/_LjN3UclYzU has a different attempt to formalize Kalman filters which I think we can all agree was not a successful formalization.

[thesmtsolver]

It is really not that difficult. Here is a paper that formalizes a version of feed forward networks to prove properties about them.

https://arxiv.org/pdf/2304.10558

[kragen]

I only skimmed this paper, but it doesn't mention floating point (it's only modeling the FNN as a function on reals), and I don't think you can extract a working FNN implementation from an SMT-LIB or Z3 problem statement, so I think you have to take on faith the correspondence between the implementation you're running and the thing you proved correct.

But that's the actual problem we're trying to solve; nobody really doubts Kalman's proof of his filter's optimality.

So this paper is not a relevant answer to the question, "how would you prove an implementation of a Kalman filter is correct?"

[ajb]

In most cases the data a kalman filter is working on has some precision which is much lower than the available precision in the floating format you are using. The problem is inherently a statistical one, since the expected precision depends on the statistics of your data source.

So you would probably adopt some conservative approach in which you showed that the worst case floating point rounding error is << some quantile of error due to the data.

But, I think specialised tools are more commonly used than general process. Eg, see https://github.com/arpra-project/arpra

[kragen]

That's a start, but you might be able to do better than that; for example, you ought to be able to show that the floating-point rounding error is unbiased.

[dwohnitmok]

The link was to show Coq's floating point library.

> As a concrete example, how would you prove an implementation of a Kalman filter is correct?

This would be fun enough to implement if I had more time. But I need to spend time with my family. I'm usually loathe to just post LLM output, but in this case, since you're looking just for the gist, and I don't have the time to write it out in detail, here's some LLM output: https://pastebin.com/0ZU51vN0

Based on a quick skim, I can vouch for the overall approach. There's some minor errors/things I would do differently, and I'm extremely doubtful it typechecks (I didn't use a coding agent, just a chat window), but it captures the overall idea and gives some of the concrete syntax. I would also probably have a third implementation of the 1-d Kalman filter and would prefer some more theorems around the correctness of the Kalman filter itself (in addition to numerical stability). And of course a lot of the theorems are just left as proof stubs at the moment (but those could be filled in given enough time).

But it's enough to demonstrate the overall outline of what it would look like.

The overall approach I would have with floating point algorithms is exemplified by the following pseudo-code for a simple single argument function. We will first need some floating point machinery that establishes the relationship between floating point numbers and real numbers.

  # Calculates the corresponding real number from the sign, mantissa, and exponent
  floatToR : Float -> R
  
  # Logical predicate that holds if and only if x can be represented exactly in floating point
  IsFloat(x: R)

Then define a set of three functions, `f0` which is on the idealized reals, `f1` which is on floating point numbers as mechanically defined via binary mantissa, exponent, and sign, and `f2` which is on the finite subset of the reals which have an exact floating point representation.

  f0 : R -> R
  f1 : Float -> Float
  f2 : (x : R) -> IsFloat(x) -> R

The idea here is that `f0` is our "ideal" algorithm. `f1` is the algorithm implemented with exact adherence to low-level floating point operations. And `f2` is the ultimate algorithm we'll extract out to runnable code (because it might be extremely tedious, messy, and horrendous for runtime performance, to directly tie `f1`'s explicit representation of mantissa, exponent, and sign, which is likely some record structure with various pointers, to the machine 32-bit float).

Then I prove the exact correspondence between `f1` and `f2`.

  f1_equal_to_f2 : forall (x: R), (isFloatProof: IsFloat(x)) -> f1(floatToR(x)) = floatToR(f2(x, isFloatProof))

This gives the following corollary FWIW.

  f2_always_returns_floats : forall (x: R), IsFloat(x) -> IsFloat(f2(x))

This lets me throw away `f1` and work directly with `f2`, which lets me more easily relate `f2` and `f0` since they both work on reals and I don't need to constantly convert between `Float` and `R`.

So then I prove numerical stability on `f2` relative to `f0`. I start with a Wilkinson style backwards error analysis.

  f2_is_backwards_stable : forall (x: R), exists (y: R), IsFloat(x) -> f2(x) = f0(y) and IsCloseTo(y, x)

Then, if applicable, I would prove a forward error analysis (which requires bounding the input and then showing that overflow doesn't happen and optionally underflow doesn't happen)

  f2_is_forwards_stable : forall (x: R), IsFloat(x) -> IsReasonablyBounded(x) -> IsCloseTo(f2(x), f0(x))

And then given both of those I might then go on to prove some conditioning properties.

Then finally I extract out `f2` to runnable code (and make sure that compilation process after that has all the right floating point flags, e.g. taking into consideration FMA or just outright disabling FMA if we didn't include it in our verification, making sure various fast math modes are turned off, etc.).

[dwohnitmok]

Error in `f1_equal_to_f2`. Should be ` f1_equal_to_f2 : forall (x: Float), floatToR(f1(x)) = f2(floatToR(x))` with an implicit proof delivered to `f2` that `forall (x: R), IsFloat(floatToR(x))`.

Also there's some subtleties in how to set up the code extraction to make sure you don't inadvertently extract other real-valued functions, but this is already too long.