Quite often in scientific computing we're interested in finding the fixed-points of a function or transformation: those points where the value of the function f at a point x is equal to x itself, ie. f(x) = x. You can see how these might be useful – if the function f models any kind of real-world phenomenon, then a fixed-point suggests a sort of equilibrium or steady-state.
In fact, fixed-point iteration can be used in many numerical methods. Newton's method for finding the zeros of a differentiable function (which we'll look at in a later post) can be written in terms of a fixed-point computation, and many methods for solving ordinary differential equations apply these same fixed-point ideas.
In this post we'll code up a Scala function to (hopefully!) converge upon the fixed-points of a given transformation, and look at some examples; including a method for finding square roots of a number. We'll also look at the general technique of average-damping that we can apply to ensure/quicken convergence in particular cases.
[Inspired by the Structure and Interpretation of Computer Programs. (noticing a pattern yet?) The book is here, and the relevant Scheme solutions are for Exercises 1.35 and 1.36.]
The fixedPoint() Function
Our method to actually find the fixed
points relies on the following idea. Suppose we have a function f and
a starting point x. We want to find the point where f(x) = x, right?
So one way to do this would be to calculate f(x) and compare it to x:
if they're equal, then we're done. If not, we simply repeat that
process with f(x) as our next starting point (so we'd compare f(x) to
f(f(x)) the next time around, and so on.)
So instead we'll define x and f(x) to be “equal” if f(x) is within a small neighbourhood of x: that is, x-T < f(x) < x+T for some small tolerance T. For our case, we'll take T=0.0001, which gives us a decent notion of “equality”.
Next, how do we know that we're ever
going to converge here?
For the purposes of this post, most of our examples will magically work out okay! For the sake of pragmatism though, we'll include an upper bound on the number of iterations we run so we don't spiral off into an infinite loop. Most of our examples will converge rapidly, so we'll set the limit at 30 iterations. If we exhaust the number of iterations without converging, we'll return NaN
(Not-a-Number) to signal this.
Finally, it'd be cool to see the intermediate results of the computation for each step. So we'll print out the current value of our 'guess' at each point.
So with those choices made, our fixedPoint() function can be written. We'll define the function to take in the given function f to work on and a starting point firstGuess. We'll then update the value of our current guess with f(guess) at each iteration, and return the guess value if we're within the specified tolerance. If we pass 30 iterations, we'll return Java's Double.NaN value.
Alright! Let's check out some examples of fixed-points!
The Golden Ratio, Re-Revisited
I told you this thing pops up everywhere!
So let's use our fixedPoint() function to get yet another approximation for the golden ratio. Here we'll use an initial guess of 1.0:
fixedPoint(x => 1 + 1/x, 1.0)
//> 2.0
//> 1.5
//> 1.6666666666666665
//> 1.6
//> 1.625
//> 1.6153846153846154
//> 1.619047619047619
//> 1.6176470588235294
//> 1.6181818181818182
//> 1.6179775280898876
//> res0: Double = 1.6179775280898876
(1 + Math.sqrt(5))/2
//> res1: Double = 1.618033988749895
So we do end up converging to a decent approximation for the golden ratio.
A solution to xx=1000
We can use our fixedPoint() function to find the solution to xx=1000 by rewriting it as follows:
Here ln refers to the natural logarithm. We have access to this with Java's Math.log method.
Now note that ln(1.0) = 0, so if have 1.0 as our initial guess we'll get a division by zero:
fixedPoint(x => Math.log(1000) / Math.log(x), 1.0)
//> Infinity
//> 0.0
//> res2: Double = 0.0
We immediately get infinity followed by an answer of 0.0, which suggests that something went a little wrong in the calculation, as we'd expect. So we need to use a slightly different initial guess, say 1.5 instead:
fixedPoint(x => Math.log(1000) / Math.log(x), 1.5)
//> 17.036620761802716
//> 2.436284152826871
//> 7.7573914048784065
//> 3.3718636013068974
//> 5.683217478018266
//> 3.97564638093712
//> 5.004940305230897
//> 4.2893976408423535
//> 4.743860707684508
//> 4.437003894526853
//> 4.6361416205906485
//> 4.503444951269147
//> 4.590350549476868
//> 4.532777517802648
//> 4.570631779772813
//> 4.545618222336422
//> 4.562092653795064
//> 4.551218723744055
//> 4.558385805707352
//> 4.553657479516671
//> 4.55677495241968
//> 4.554718702465183
//> 4.556074615314888
//> 4.555180352768613
//> 4.555770074687025
//> 4.555381152108018
//> 4.555637634081652
//> 4.555468486740348
//> 4.555580035270157
//> res3: Double = 4.555580035270157
As you can see the convergence isn't quite as quick as for the golden ratio example, but we still get there in the end.
The Fixed Point of the Cosine Function
As one more easy example, we can apply our fixedPoint() method to find the fixed-point of the cos() function (since it's the last basic trigonometric function we haven't mentioned on the blog yet):
fixedPoint(Math.cos, 1.0)
//> 0.5403023058681398
//> 0.8575532158463934
//> 0.6542897904977791
//> 0.7934803587425656
//> 0.7013687736227565
//> 0.7639596829006542
//> 0.7221024250267077
//> 0.7504177617637605
//> 0.7314040424225098
//> 0.7442373549005569
//> 0.7356047404363474
//> 0.7414250866101092
//> 0.7375068905132428
//> 0.7401473355678757
//> 0.7383692041223232
//> 0.7395672022122561
//> 0.7387603198742113
//> 0.7393038923969059
//> 0.7389377567153445
//> 0.7391843997714936
//> 0.7390182624274122
//> 0.7391301765296711
//> res4: Double = 0.7391301765296711
Still quite a few iterations. We'll see a way to speed up the convergence shortly though, I promise!
Finding Square Roots
Right now we'll look at a more involved example though. We can use our fixedPoint() function to find the square root of a given number, using the following mathematical rearrangement:
So using this, we should be able to find the square root of 2:
fixedPoint(y => 2/y, 1.0)
//> 2.0
//> 1.0
//> 2.0
//> 1.0
//> 2.0
//> ... lots of bounces ...
//> ...
//> 1.0
//> 2.0
//> 1.0
//> res5: Double = NaN
Wow! What happened? That should have worked, right?
Actually, we have a bit of a problem here. Observe the following:
So we keep bouncing between the initial guess and the value of x divided by the initial guess, which exactly explains the output we get! We have an infinite oscillation between two values, so we never converge.
Fortunately, we can fix this using the technique of average-damping. The oscillation hints at our problem: we're bouncing too far on each iteration. So what we want to do is reduce the distance we jump (ie. We'll be damping the oscillation.)
Now we can see that the actual square root of x has to occur somewhere between the point y and the value x/y. (Why? Suppose √x is not between them; for instance, suppose √x < y and
√x < x/y. Then we have that y√x < x, ie. y < √x. So √x < y and, at the same time,
y < √x; a mathematical impossibility. So we have a contradiction. (The same argument applies to the case where √x > y and √x > x/y.) )
So using this, we can take our next guess at any point to be the average of the current guess and the value of the function at that guess. This will put us somewhere between them, which means that we're closer to the actual square root. If all goes according to plan, this should stop the oscillations since we won't be bouncing as far this time around.
Mathematically, this means we'll be solving the following equivalent fixed-point problem:
So you can see that we should still get the same answer for the square root if this works. So let's give it a try:
fixedPoint(y => (2/y + y)/2, 1.0)
//> 1.5
//> 1.4166666666666665
//> 1.4142156862745097
//> res6: Double = 1.4142156862745097
Math.sqrt(2)
//> res7: Double = 1.4142135623730951
Nice! So we actually arrived at an answer, and pretty fast too!
In fact, given that the solution to xx=1000 and the cosine fixed-point approximation also exhibited the same sort of back-and-forth bouncing on the way to the answer, we should be able to speed up the convergence for them by applying this same average-damping idea.
To help with this we'll define the higher-order function averageDamp() that takes in a function f, and returns a function that takes a value x and returns the average of x and the function f evaluated at x.
We'll then simply call this function as shown:
def averageDamp(f: Double => Double): Double => Double =
(x => (x + f(x))/2)
//> averageDamp: (f: Double => Double)Double => Double
fixedPoint(averageDamp(x => Math.log(1000) / Math.log(x)), 1.5)
//> 9.268310380901358
//> 6.185343522487719
//> 4.988133688461795
//> 4.643254620420954
//> 4.571101497091747
//> 4.5582061760763715
//> 4.555990975858476
//> 4.555613236666653
//> res8: Double = 4.555613236666653
fixedPoint(averageDamp(Math.cos), 1.0)
//> 0.7701511529340699
//> 0.7439782957302291
//> 0.7398792505064895
//> 0.7392146118880453
//> 0.7391062602582958
//> res9: Double = 0.7391062602582958
Awesome! The convergence was quite rapid indeed.
So now that this works, we've actually got a cool implementation for a square root method!
And this square root method will work quite well in practice, too! We can compare it to the Math.sqrt() method that already exists in the Java.lang.Math class:
sqrt(10)
//> res10: Double = 3.162277665175675
Math.sqrt(10)
//> res11: Double = 3.1622776601683795
sqrt(100)
//> res12: Double = 10.000052895642693
Math.sqrt(100)
//> res13: Double = 10.0
sqrt(12345)
//> res14: Double = 111.10805770848404
Math.sqrt(12345)
//> res15: Double = 111.1080555135405
This is a pretty good square-root finder, considering the tolerance value we used! In fact, many square root implementations (particularly in calculators) will use Newton's method. We'll see later that Newton's method can be written directly in terms of our fixed-point iteration, so we almost have the makings of an industrial-strength implementation of a sqrt() function right here!
Java Mathematician soooo cooool
ReplyDelete