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.
Continue reading: Iterative Square Root
Insights from the exploration of Square Root Finding Algorithms
The original text examined the workings of a square root finding algorithm using iterative convergence solutions. The core idea was applying Newton’s method of evaluation in different programming languages: Go, R, and APL. The algorithm loops through a calculation, adjusting the value of x, until the calculated root stops changing, thus finding the square root of any given number.
Long-term implications
Examining the process of iterative calculations and their convergence to an accurate value serves as a useful perspective on the practical implementation of mathematical concepts in programming. The iterative square root finding algorithm has wide-reaching implications for fields where precision is crucial, such as physical calculations in engineering, data approximation in machine learning, econometrics, or GPS satellite navigation.
It also underscores how different programming languages enable problem-solving in different ways. The iterations used in Go programming were replicated in R and APL – showing that these principles can be applied across programming languages, with varying degrees of complexity or conciseness.
Possible future developments
These processes and methods could become more relevant as computational requirements for data analysis continue to grow. Refining and optimizing such iterative calculation functions will contribute to the development of these fields. Furthermore, comparative studies of programming languages – such as illustrated here with Go, R, and APL – can highlight the efficiency, speed, or readability of different languages, guiding programmers towards the most suitable tool for their task.
Actionable Advice
- Improvement of existing functions: Existing iterative solutions can use these insights to relax numeric precision for performance gains or increase it for more robust results.
- Choosing the right language: comparing how different languages handle iterative calculations can guide programmers to make more informed decisions on choosing the best language for a given task.
- Constant learning and adaptation: As demonstrated, even common algorithms like square root calculation can have room for reinterpretation and improvement. Keeping an open mind for new approaches can lead to surprisingly elegant or efficient methods.
- Multiply Applicability: Once learned and understood, the iterative calculation method can be applied to ways outside finding square roots. The same logic can be used in other areas of computation where iterative improvement towards a converging value is required.
In conclusion, the article prompts us to keep digging deeper into programming constructs. Learnings from one area could potentially be transferred or adapted into others. It encourages adopting a cross-language perspective and constant learning as hallmarks of good programming practice.