Sunday, January 28, 2024

Exponentiating Functions

If we compose a function F(x) with itself, we get F(F(x)) or F∘F. If we compose it again, we get F(F(F(x))) or F∘F∘F. Rather than writing ‘F’ over and over again, we can abuse exponential notation and write (F∘F∘F) as F3, where the superscript indicates how many times we compose the function. F1 is, naturally, just F. F0 would be zero applications of F, which would be the function that leaves its argument unchanged, i.e. the identity function.

The analogy with exponentiation goes deeper. We can quickly exponentiate a number with a divide and conquer algorithm:

(defun my-expt (base exponent)
  (cond ((zerop exponent) 1)
        ((evenp exponent) (my-expt (* base base) (/ exponent 2)))
        (t (* base (my-expt base (- exponent 1))))))
The analagous algorithm will exponentiate a function:
(defun fexpt (f exponent)
  (cond ((zerop exponent) #'identity)
        ((evenp exponent) (fexpt (compose f f) (/ exponent 2)))
        (t (compose f (fexpt f (- exponent 1))))))

The function (lambda (x) (+ 1 (/ 1 x))) takes the reciprocal of its input and adds one to it. What happens if you compose it with itself? We can rewrite it as a linear fractional transformation (LFT) (make-lft 1 1 1 0) and try it out:

> (make-lft 1 1 1 0)
#<LFT 1 + 1/x>

> (compose * (make-lft 1 1 1 0))
#<LFT (2x + 1)/(x + 1)>

> (compose * (make-lft 1 1 1 0))
#<LFT (3x + 2)/(2x + 1)>

> (compose * (make-lft 1 1 1 0))
#<LFT (5x + 3)/(3x + 2)>

> (compose * (make-lft 1 1 1 0))
#<LFT (8x + 5)/(5x + 3)>

> (compose * (make-lft 1 1 1 0))
#<LFT (13x + 8)/(8x + 5)>

> (compose * (make-lft 1 1 1 0))
#<LFT (21x + 13)/(13x + 8)>
Notice how the coefficients are Fibonacci numbers.

We can compute Fibonacci numbers efficiently by exponentiating the LFT 1 + 1/x.

> (fexpt (make-lft 1 1 1 0) 10)
#<LFT (89x + 55)/(55x + 34)>

> (fexpt (make-lft 1 1 1 0) 32)
#<LFT (3524578x + 2178309)/(2178309x + 1346269)>
Since we’re using a divide and conquer algorithm, raising to the 32nd power involves only five matrix multiplies.

Fixed points

If an input x maps to itself under function f, we say that x is a fixed point of f. So suppose we have a function f with a fixed point x. We consider the function f’ which ignores its argument and outputs x. If we compose f with f’, it won’t make difference that we run the result through f again, so f’ = f∘f’ = f. You can find a fixed point of a function by composing the function with its fixed point. Unfortunately, that only works in a lazy language, so you have two options: either choose a finite number of compositions up front, or compose on demand.

You can approximate the fixed point of a function by exponentiating the function to a large number.

(defun approx-fixed-point (f) 
  (funcall (fexpt f 100) 1))

> (float (approx-fixed-point (make-lft 1 1 1 0)))
1.618034

> (float (funcall (make-lft 1 1 1 0) *))
1.618034

Alternatively, we could incrementally compose f with itself as needed. To tell if we are done, we need to determine if we have reached a function that ignores its input and outputs the fixed point. If the function is a LFT, we need only check that the limits of the LFT are equal (up to the desired precision).

(defun fixed-point (lft)
  (let ((f0   (funcall lft 0))
        (finf (funcall lft ’infinity)))  
    (if (and (numberp f0)
             (numberp finf)
             (= (coerce f0 ’single-float) 
                (coerce finf ’single-float)))
        (coerce f0 ’single-float)
        (fixed-point (compose lft lft)))))

> (fixed-point (make-lft 1 1 1 0))
1.618034    

No comments: