Potential Titles

Motivation 1: Testing (raw floating point errors)

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

Motivation 2: Matrix Problems and Stability

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…

Outline

  1. R and floating points
    • IEEE representation & underlying structure
  2. Understanding & Avoiding Errors (algebraic)
  3. Understanding & Avoiding Errors (linear algebra)
  4. Summary

Part I

Introduction to floating points

Introduction

print(.1 + .2); print(.3)
## [1] 0.3
## [1] 0.3
.1 + .2 == .3
## [1] FALSE

What’s really happening?

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