living in a discretized, computerized world

- “Computing: real numbers need not apply”
- “Floating-point-land: A
`R`

omance of many bits” - “Floating points: When you think you’re a mathematician”
- “I wrote the tests correctly though!”

- suppose you’ve coded up a really good example showing your functions work amazingly and you’ve checked your math (with a quick print statement)

`print(.1 + .2); print(.3)`

`## [1] 0.3`

`## [1] 0.3`

and you think “wow, I really like coding, and I’m like a good mathematician”, and then this happens…

`.1 + .2 == .3`

`## [1] FALSE`

And that’s when you go into pure mathematics cause no one likes computers anyway.

*Technically here’s a correction if you encounter this (but we’ll talk about why you did soon):*

`testthat::expect_equal(.1 + .2, .3, tolerance = 0) # says something useful`

```
## Error: 0.1 + 0.2 not equal to 0.3.
## 1/1 mismatches
## [1] 0.3 - 0.3 == 5.55e-17
```

```
testthat::expect_equal(.1 + .2, .3, tolerance = .Machine$double.eps)
# or
all.equal(.1 + .2, .3, tolerance = .Machine$double.eps)
```

`## [1] TRUE`

Suppose you came back to your computer and given that you’re into Stats, you’ve coded up some iterative algorithm to solve some linear algebra stuffs (you know matrix multiplication, inversion etc).

Say you have some code that can solve \(Ax = b\) where

\[ A = \begin{bmatrix} .78 &.563 \\ .913 & .659 \end{bmatrix} \quad \quad b = \begin{bmatrix} .217 \\ .254 \end{bmatrix} \]

Now, as mathematicians we know that \(x = A^{-1}b\) (aka we have a closed form solution), but for this example, we will just use the fact that we know the answer mathematically and we’d like to assess preformance of our algorithm through tests.

So, the set up:

```
A <- matrix(c(.78, .563,
.913, .659),
byrow = T, nrow = 2)
b <- c(.217, .254)
x_true <- solve(a = A,b = b); x_true
```

`## [1] 1 -1`

That is, we know the true solution is \(\begin{bmatrix} 1 \\ -1 \end{bmatrix}\).

So you decide that you want your answer (from your function) to be close to solving the equation, that is \(A\hat{x} - b\) is really close to zero (say \(l_2\) norm `< 1e-5`

). Note that, since we know the true solution (\(x\)), we will also be checking the **relative error** of the \(\hat{x}\) (your answer from your function), i.e. \[|x-\hat{x}|/|x|\]

The first function you write returns the following value for \(\hat{x}_1\):

```
x_hat1 <- c(.999,-.999)
res_1 <- b - A %*% x_hat1
res_1
```

```
## [,1]
## [1,] 0.000217
## [2,] 0.000254
```

`sprintf("l_2 norm of residual is %.6f", norm(res_1))`

`## [1] "l_2 norm of residual is 0.000471"`

`norm(x_true - x_hat1, type = "2")/norm(x_true, type = "2") `

`## [1] 0.0009999999`

A friend says they think they have a function that solved the problem, and sends over \(\hat{x}_2\).

```
x_hat2 <- c(.341, -.087)
res_2 <- b - A %*% x_hat2
res_2
```

```
## [,1]
## [1,] 1e-06
## [2,] 0e+00
```

`sprintf("l_2 norm of residual is %.6f", norm(res_2))`

`## [1] "l_2 norm of residual is 0.000001"`

`norm(x_true - x_hat2, type = "2")/norm(x_true, type = "2")`

`## [1] 0.7961941`

Dang it…

- R and floating points
- IEEE representation & underlying structure

- Understanding & Avoiding Errors (algebraic)
- Understanding & Avoiding Errors (linear algebra)
- Summary

*Introduction to floating points*

`print(.1 + .2); print(.3)`

`## [1] 0.3`

`## [1] 0.3`

`.1 + .2 == .3`

`## [1] FALSE`

What’s really happening?

\(0.1 + 0.2\):

```
## 0.1000000000000000055511151231257827021181583404541015625
## + 0.2000000000000000111022302462515654042363166809082031250
## ------------------------------------------------------------
## 0.3000000000000000444089209850062616169452667236328125000
```

vs. \(0.3\):

`## 0.2999999999999999888977697537484345957636833190917968750`

**why is “\(.3\)” is printed?**

`#?print.default`