What do We Have to Know About a Function in Order to Compute its Definite Integral?

Suppose that f is a continuous function from \lbrack 0,1\rbrack to \lbrack 0,1\rbrack and that we have a program which computes it. (Ignore for now exactly what it means to “compute” a real-valued function of the reals. Suffice it to say that almost every natural continuous function you come across is computable). If we want to compute \int_{0}^1 f(x)\,dx, say to within an accuracy of 0.001, how would we do it, and would we need any further information about f in order to do it?

The obvious thing to do is to compute a Riemann sum R(f,n)=\sum_{i=0}^{n-1} (1/n)\cdot f(i/n) for some large n. However, this could be arbitrarily far from the true value of \int_0^1 f(x)\,dx. For example, f(i/n) might be 0 for all 0 \leq  i \leq n-1, but f might curve sharply up in between each i/n and (i+1)/n so that its definite integral is arbitrarily close to 1.

However, since f is continuous on \lbrack 0,1\rbrack, it is uniformly continuous. This means that for all \epsilon > 0 there is a \delta > 0 such that whenever |x - y| < \delta, |f(x) - f(y)| < \epsilon. If we could compute a function \delta(\epsilon) such for all \epsilon > 0, whenever |x - y| < \delta(\epsilon), |f(x) - f(y)| < \epsilon, then we could compute the definite integral of f with arbitrary accuracy: If we want to know \int_0^1 f(x)\, dx to within 0.001, then choose n so that 1/n < \delta(0.001) and we can take R(f,n), since |R(f,n) - \int_0^1 f(x)\,dx| < 0.001.

So, one answer to the question “What extra information do we need in order to compute \int_0^1 f(x)\,dx?” is: a function \delta(\epsilon) witnessing f‘s uniform continuity.

In 1998, Alex Simpson showed, building on ideas of Ulrich Berger and others, that another answer is: Nothing!

For technical convenience, we switch from considering functions with domain and range \lbrack 0,1\rbrack to functions with domain and range \lbrack -1,1\rbrack, and try to compute the definite integral of such a function over \lbrack 0,1\rbrack, since that value will again be in \lbrack -1,1\rbrack.

For our representation of the reals in \lbrack -1,1\rbrack we will use functions \alpha from \mathbb{N}^+ = \{1,2,3,\ldots\} to \{-1,0,1\}, where we are thinking of \alpha as representing \sum_{i=1}^\infty \alpha(i) 2^{-i}. All reals in \lbrack -1,1\rbrack have representations. Some have more than one, but that won’t concern us.

We will use \lbrack -1,1\rbrack to denote both the interval of reals from -1 to 1 and the set of its representations. Hopefully this will not cause confusion.

So a computable function from \lbrack -1,1\rbrack to \lbrack -1,1\rbrack is a computable second-order functional which takes a function from \mathbb{N}^+ to \{-1,0,1\} and returns a function from \mathbb{N}^+ to \{-1,0,1\}. (There are some technicalities about the fact that a real can have more than one representation, but we will ignore them.)

Similarly, a predicate on \lbrack -1,1\rbrack is a function from \lbrack -1,1\rbrack to \{\text{true}, \text{false}\}. An amazing fact due to Ulrich Berger is that for any computable predicate P(x) we can decide whether or not there is an x\in \lbrack -1,1\rbrack such that P(x) holds. (To see why this is amazing, observe that the analogous statement for computable predicates on the integers is false, due to the Halting Problem.)

What’s more, we can define computable functionals \text{exists} and \text{witness} so that, if P(x) is a computable predicate, \mathop{\text{exists}}(P) returns \text{true} or \text{false} depending on whether or not \exists x\,P(x) and \mathop{\text{witness}}(P) returns a number x which satisfies P(x) if anything does.

Before we define those functions, I’ll give some notation. If \alpha is a function from \mathbb{N}^+ to \{-1,0,1\}, then we let 0 :: \alpha be defined by (0 :: \alpha)(0) = 0 and (0 :: \alpha)(i+1) = \alpha(i) for i\in\mathbb{N}^+. We define 1::\alpha and {-1} ::\alpha similarly.

Now define \text{exists} and \text{witness} by mutual recursion as follows: Let \mathop\text{exists}(P) = P(\mathop\text{witness}(P)) and let \mathop\text{witness}(P) do the following: it uses \text{exists} checks to see if there is an \alpha such that \alpha(0) = -1 and P(\alpha) holds and, if so, it returns {-1}::(\mathop\text{witness}(P')), where P'(\beta) = P(-1::\beta). If not, it uses \text{exists} to check if there is an \alpha such that \alpha(0) = 0 and P(\alpha) holds and if so, it returns 0::(\mathop\text{witness}(P')), where P'(\beta) = P(0::\beta). If both of those fail, it returns 1::(\mathop\text{witness}(P')), where P'(\beta) = P(1::\beta).

If you think about the above definitions for a while, it should probably become clear to you that they should work, if they terminate. But they seem way too crazy to terminate. In fact they do, given one condition which I’ll remark on below, and the crucial fact is that since P is computable and total, when given an \alpha, it only looks at finitely many values \alpha(n) before terminating.

The condition that allows them to terminate is that we specify that the evaluation is lazy. That is, when I say for example that \text{witness} returns {-1}::(\mathop\text{witness}(P')), I mean that it returns the function which on input {0} returns {-1} and only on non-zero input calculates \mathop\text{witness}(P') to determine its return value for that input. So, it only calls \text{witness} when it needs to. Since every call to \text{witness} involves a recursive call to \text{witness}, if you followed all of them immediately, the program could not possibly terminate.

Note also that we can use \text{exists} to define a \text{forall} as well.

Now let’s use this to compute the definite integral. Observe that if \alpha(0) = -1, then \alpha \in \lbrack -1,0\rbrack, if \alpha(0)=0, then \alpha\in \lbrack -1/2,1/2\rbrack, and if \alpha(0) = 1, then \alpha\in \lbrack 0,1\rbrack. Conversely, every real in \lbrack -1,0\rbrack has a representation \alpha such that \alpha(0) = -1, and similarly for the other two cases.

Therefore, if f is such that f(\alpha)(0) = -1 for all \alpha\in \lbrack -1,1\rbrack, it must be the case that \int_0^1 f(x)\, dx \in\lbrack -1,0\rbrack, and if we are trying to compute a representative \alpha for \int_0^1 f(x)\,dx, we may safely say that \alpha(0) = -1, and similarly for the other two cases. Furthermore, if f(\alpha)(0) = -1 for all \alpha \in \lbrack -1,1\rbrack. Then \int_0^1 f(x)\,dx = -1/2 + (1/2)\int_0^1 (2f(x) + 1)\,dx.

Next, observe that, for any f, by stretching it out and splitting it up, we may find f_0 and f_1 such that \int_0^{1/2} f(x)\,dx = (1/2)\int_0^1 f_0(x)\,dx, \int_{1/2}^1 f(x)\,dx = (1/2)\int_0^1 f_1(x)\,dx and the range of f on \lbrack 0, 1/2\rbrack is equal to the range of f_0 on \lbrack 0,1\rbrack and the range of f on \lbrack 1/2,1\rbrack is equal to the range of f_1 on \lbrack 0,1\rbrack

Therefore, the following definition for \text{integrate} (almost) works: Given f, first use \text{forall} to check if f(\alpha)(0) = -1 for all \alpha\in \lbrack 0,1\rbrack. If so, return -1::(\mathop\text{integrate}(2f + 1)). If not, check if f(\alpha)(0) = 0 for all \alpha \in\lbrack 0,1\rbrack or if f(\alpha)(0) = 1 for all \alpha\in\lbrack 0,1\rbrack and do the analogous thing. Otherwise, return (1/2)(\mathop\text{integrate}(f_0) + \mathop\text{integrate}(f_1)).

Since f is continuous, the process of splitting it up will eventually terminate, as will the algorithm.

Actually this algorithm doesn’t quite work, as there are some a couple of details that I’ve omitted, but it is essentially correct. For more information, see Alex Simpson’s paper, and for more on Berger’s work, see Martin Escardó’s blog post, both of which include actual runnable code.

5 thoughts on “What do We Have to Know About a Function in Order to Compute its Definite Integral?

  1. Yeah, it is pretty neat.

    I think that the idea is that a computable real number is (among other equivalent characterizations) a computable function which takes an n and returns a rational f(n) so that (f(1) - 1, f(1) + 1) \supset (f(2) - 1/2, f(2) + 1/2) \supset (f(3) - 1/3, f(3) + 1/3) \supset \cdots.

    Using a representation with digits \{-1,0,1\} is equivalent to this, but just using \{0,1\}, for example, isn’t. This is because if you’re given a method for computing the binary expansion of a real number x in \lbrack 0,1\rbrack, you could, for example, be able to confidently assert after looking at one bit either that x \geq 1/2 or that x \leq 1/2, which is not something that you’re guaranteed to be able to do with a computable real number in \lbrack 0,1\rbrack.

    Or have I completely misunderstood your objection?

  2. That is indeed the heart of my question.

    Integrating using with the usual “strict” binary representation actually requires a modulus of continuity. So this result is very surprising from that perspective. I had never seen such a practical example of the non-computable equivalence of these two representations of computable reals.

    I always thought that interval arithmetic was just a way to achieve better certainty, at a reasonable price. Now I understand that it is actually more powerful than floating-point arithmetic in a very concrete sense.

    I wonder if there is a more practical version of this. More precisely, is there a nice example of an integration problem for which there is a P-time algorithm to approximate it within an interval of length 1/2^n, but computing the n-th bit is NP-hard? (I suspect it’s not too difficult code a lot of things as integrals of trigonometric polynomials, but I don’t recall ever seeing that.)

  3. Hmm, that’s interesting. Do you have a reference for the necessity of a modulus of continuity for computing a definite integral with the regular binary representation?

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s