[This article was first published on rstats on Irregularly Scheduled Programming, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.I saw a toot celebrating a short, clean implementation of a square rootfinding algorithm and wanted to dig a bit deeper into how it works, with adiversion into some APL.
This was the toot from JimGardner
Doubly pleased with myself.
Been doing the Tour of Go. Got to the section where you make a square root function, which should return once the calculated value stops changing. Struggled for ages. Trimmed and trimmed. Until finally… this!
The calculation for z was given, and I don’t understand it at all. But I don’t care. It was a total mess when I started and has turned out quite neat. I’m very satisfied.
But why “doubly pleased”? Because I’ve been solely using Neovim so far for Go!!
package main import ( "fmt" ) func Sqrt(x float64) float64 {z := 1.0 for { y := z z -= (z*z - x) / (2 * z) if z == y { return z } }} func main() { fmt.Println(Sqrt(16)) }
It’s a nice, not-too-complicated algorithm to play with, and I agree it’s hard tosee why it works for this application, so I thought it would be neat to walkthrough that.
What we’re trying to solve here is the function (y = x^2) which we could write as(f(x) = x^2 – y) for which we want the value (x) where (f(x) = 0).
Newton’s Method is an iterativemethod for solving equations of this type (not all equations, mind you – I have anentire chapter of my PhD thesis discussing exactly why it can’t be used tosolve the equations I was solving that my supervisor insisted it could). It worksby using the slope (derivative) at a point to guide towards a solution. The formulafor the updated value (x_{n+1}) given some guess (x_n) is
[x_{n+1} = x_n – \frac{f(x_n)}{f'(x_n)}]where (f'(x)) is the derivative of the function (f) at the point (x). For (f(x) = x^2 – y)the derivative is (f'(x) = 2x) so we can substitute this and (f(x)) into the above formula
[x_{n+1}=x_n−\frac{x_n^2-y}{2x_n}]This is what the Go code calculates; given an initial guess of (x_n = 1) it calculates the next value as
x = x - (x*x - y) / (2 * x)
where, here, y is the value we’re finding the square root of.
In R this could be written as
SQRT <- function(x) { z <- 1 while (TRUE) { y <- z z <- z - (z*z - x)/(2*z) if (abs(y -z) < 0.0001) return(z) }}
(since base::sqrt is already defined) where I’ve used a tolerance rather thanrelying on exact numerical equality. The while(TRUE) construct is equivalentto Go’s for {} syntax; an infinite loop.
R actually has another way to write that which is even closer; repeat {}
SQRT <- function(x) { z <- 1 repeat { y <- z z <- z - (z*z - x)/(2*z) if (abs(y -z) < 0.0001) return(z) }}
One might notice that this approach requires essentially squaring a value, whichis hardly expensive, but we can simplify and cancel out (x_n), so
[x_{n+1}=\frac{x_n-\frac{y}{x_n}}{2}]in which case we have
SQRT <- function(x) { z <- 1 repeat { y <- z z <- (z + x/z)/2 if (abs(y -z) < 0.0001) return(z) }}
One of the reasons I wanted to dig into this was the fact that it’s a convergence…
In APL the power operator (⍣ aplwi)applies a function some specified number of times, so
f ⍣n x
applies f to x n times, i.e. (f⍣3)x produces f(f(f(x))).
It can also be used as ⍣= where it will continue to apply the function untilthe output no longer changes (is equal). A classic example is thegolden ratio; take the reciprocalthen add 1 until it converges, i.e.
[x_{n+1} = 1+\frac{1}{x_n}]
which you can try for yourselfhere
1+∘÷⍣=11.618033989
In this, +∘÷ is the (tacit) function created bycomposing (∘) ‘addition of 1’ (1+, a partial application of a function) and‘reciprocal’ (÷), which is iterated until it no longer changes (withinmachine precision).
Iterating until convergence is exactly what we want, since we’re looking for thevalue satisfying
[x_n = x_{n+1}=\frac{x_n-\frac{y}{x_n}}{2}]APL uses ⍵ as the right argument placeholder and ⍺ as the left, so thefunction we want to apply repeatedly to the right argument is
{⍵-(((⍵×⍵)-⍺)÷(2×⍵))}
If we provide 1 as the right argument (the start value) and 16 as the leftargument, we get
16{⍵-(((⍵×⍵)-⍺)÷(2×⍵))}⍣=14
You can try this out yourself at tryapl.org(link should load that expression).
We can turn that into a function, once again using the argument placeholder
sqrt←{⍵{⍵-(((⍵×⍵)-⍺)÷(2×⍵))}⍣=1} sqrt 255 sqrt 819
Taking the simplification above, we can write this a bit shorter as
sqrt←{⍵{(⍵+(⍺÷⍵))÷2}⍣=1} sqrt 14412
As clean as the Go code looks, I think there’s a certain beauty to being ableto write this in just 20 characters. It’s not for everyone, I get that.
I love these opportunities to learn a bit more about languages.
If you have comments, suggestions, or improvements, as always, feel free to usethe comment section below, or hit me up onMastodon.
devtools::session_info() ```
``` To leave a comment for the author, please follow the link and comment on their blog: rstats on Irregularly Scheduled Programming.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.Continue reading: Iterative Square Root