Iterative Square Root
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 root finding algorithm and wanted to dig a bit deeper into how it works, with a diversion into some APL.
This was the toot from Jim Gardner
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 to see why it works for this application, so I thought it would be neat to walk through 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 iterative method for solving equations of this type (not all equations, mind you – I have an entire chapter of my PhD thesis discussing exactly why it can’t be used to solve the equations I was solving that my supervisor insisted it could). It works by using the slope (derivative) at a point to guide towards a solution. The formula for 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 than
relying on exact numerical equality. The while(TRUE)
construct is equivalent
to 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, which is 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 until
the output no longer changes (is equal). A classic example is the
golden ratio; take the reciprocal
then add 1 until it converges, i.e.
\[ x_{n+1} = 1+\frac{1}{x_n} \]
which you can try for yourself here
1+∘÷⍣=1 1.618033989
In this, +∘÷
is the (tacit) function created by
composing (∘
) ‘addition of 1’ (1+
, a partial application of a function) and
‘reciprocal’ (÷
), which is iterated until it no longer changes (within
machine precision).
Iterating until convergence is exactly what we want, since we’re looking for the value 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 the
function 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 left
argument, we get
16{⍵-(((⍵×⍵)-⍺)÷(2×⍵))}⍣=1 4
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 25 5 sqrt 81 9
Taking the simplification above, we can write this a bit shorter as
sqrt←{⍵{(⍵+(⍺÷⍵))÷2}⍣=1} sqrt 144 12
As clean as the Go code looks, I think there’s a certain beauty to being able to 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 use the comment section below, or hit me up on Mastodon.
devtools::session_info()
## ─ Session info ─────────────────────────────────────────────────────────────── ## setting value ## version R version 4.3.3 (2024-02-29) ## os Pop!_OS 22.04 LTS ## system x86_64, linux-gnu ## ui X11 ## language (EN) ## collate en_AU.UTF-8 ## ctype en_AU.UTF-8 ## tz Australia/Adelaide ## date 2024-05-29 ## pandoc 3.1.11 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/x86_64/ (via rmarkdown) ## ## ─ Packages ─────────────────────────────────────────────────────────────────── ## package * version date (UTC) lib source ## blogdown 1.18 2023-06-19 [1] CRAN (R 4.3.2) ## bookdown 0.36 2023-10-16 [1] CRAN (R 4.3.2) ## bslib 0.6.1 2023-11-28 [3] CRAN (R 4.3.2) ## cachem 1.0.8 2023-05-01 [3] CRAN (R 4.3.0) ## callr 3.7.3 2022-11-02 [3] CRAN (R 4.2.2) ## cli 3.6.2 2023-12-11 [3] CRAN (R 4.3.2) ## crayon 1.5.2 2022-09-29 [3] CRAN (R 4.2.1) ## devtools 2.4.5 2022-10-11 [1] CRAN (R 4.3.2) ## digest 0.6.34 2024-01-11 [3] CRAN (R 4.3.2) ## ellipsis 0.3.2 2021-04-29 [3] CRAN (R 4.1.1) ## evaluate 0.23 2023-11-01 [3] CRAN (R 4.3.2) ## fastmap 1.1.1 2023-02-24 [3] CRAN (R 4.2.2) ## fs 1.6.3 2023-07-20 [3] CRAN (R 4.3.1) ## glue 1.7.0 2024-01-09 [3] CRAN (R 4.3.2) ## htmltools 0.5.7 2023-11-03 [3] CRAN (R 4.3.2) ## htmlwidgets 1.6.2 2023-03-17 [1] CRAN (R 4.3.2) ## httpuv 1.6.12 2023-10-23 [1] CRAN (R 4.3.2) ## icecream 0.2.1 2023-09-27 [1] CRAN (R 4.3.2) ## jquerylib 0.1.4 2021-04-26 [3] CRAN (R 4.1.2) ## jsonlite 1.8.8 2023-12-04 [3] CRAN (R 4.3.2) ## knitr 1.45 2023-10-30 [3] CRAN (R 4.3.2) ## later 1.3.1 2023-05-02 [1] CRAN (R 4.3.2) ## lifecycle 1.0.4 2023-11-07 [3] CRAN (R 4.3.2) ## magrittr 2.0.3 2022-03-30 [3] CRAN (R 4.2.0) ## memoise 2.0.1 2021-11-26 [3] CRAN (R 4.2.0) ## mime 0.12 2021-09-28 [3] CRAN (R 4.2.0) ## miniUI 0.1.1.1 2018-05-18 [1] CRAN (R 4.3.2) ## pkgbuild 1.4.2 2023-06-26 [1] CRAN (R 4.3.2) ## pkgload 1.3.3 2023-09-22 [1] CRAN (R 4.3.2) ## prettyunits 1.2.0 2023-09-24 [3] CRAN (R 4.3.1) ## processx 3.8.3 2023-12-10 [3] CRAN (R 4.3.2) ## profvis 0.3.8 2023-05-02 [1] CRAN (R 4.3.2) ## promises 1.2.1 2023-08-10 [1] CRAN (R 4.3.2) ## ps 1.7.6 2024-01-18 [3] CRAN (R 4.3.2) ## purrr 1.0.2 2023-08-10 [3] CRAN (R 4.3.1) ## R6 2.5.1 2021-08-19 [3] CRAN (R 4.2.0) ## Rcpp 1.0.11 2023-07-06 [1] CRAN (R 4.3.2) ## remotes 2.4.2.1 2023-07-18 [1] CRAN (R 4.3.2) ## rlang 1.1.3 2024-01-10 [3] CRAN (R 4.3.2) ## rmarkdown 2.25 2023-09-18 [3] CRAN (R 4.3.1) ## rstudioapi 0.15.0 2023-07-07 [3] CRAN (R 4.3.1) ## sass 0.4.8 2023-12-06 [3] CRAN (R 4.3.2) ## sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.3.2) ## shiny 1.7.5.1 2023-10-14 [1] CRAN (R 4.3.2) ## stringi 1.8.3 2023-12-11 [3] CRAN (R 4.3.2) ## stringr 1.5.1 2023-11-14 [3] CRAN (R 4.3.2) ## urlchecker 1.0.1 2021-11-30 [1] CRAN (R 4.3.2) ## usethis 2.2.2 2023-07-06 [1] CRAN (R 4.3.2) ## vctrs 0.6.5 2023-12-01 [3] CRAN (R 4.3.2) ## xfun 0.41 2023-11-01 [3] CRAN (R 4.3.2) ## xtable 1.8-4 2019-04-21 [1] CRAN (R 4.3.2) ## yaml 2.3.8 2023-12-11 [3] CRAN (R 4.3.2) ## ## [1] /home/jono/R/x86_64-pc-linux-gnu-library/4.3 ## [2] /usr/local/lib/R/site-library ## [3] /usr/lib/R/site-library ## [4] /usr/lib/R/library ## ## ──────────────────────────────────────────────────────────────────────────────
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.