Ch2.3.2 Loss of Significance

• Round-off error often occurs near last digit.

pi

[1] 3.141593

round(pi,3)

[1] 3.142

round(pi,4)

[1] 3.1416

• For R, round-off at a “5” will typically “go towards the even digit” which is one of the IEEE standards.
• But the rounding depends on underlying number being stored rather than number being displayed (before the rounding).

x <- c(2.3456,0.0045,0.0035,0.0345,0.0355,5.5555)
round(x,3)

[1] 2.346 0.004 0.004 0.034 0.035 5.556

• Consider this subtraction problem:

1/3 - 0.33333333333333

[1] 3.330669e-15

• Now 1/3 is stored as 53 binary digits, which is around 15 decimal digits of significance.
• Thus the floating point representation does not have significant digits after the 15th.
  
• Our subtraction result must be the round-off error.

• This is called loss of significance, or cancellation error.
• Cancellation will occur when subtracting two close numbers when at least one is not perfectly represented in binary.
• This error can intensify if the subtraction is followed by a multiplication, for which the error will increase in magnitude.

(1/3 - 0.33333333333333)*1000

[1] 3.330669e-12

• Round-off error can change intermediate results.
• Subtraction is not necessarily commutative:

20.55 - 19.2 - 1.35

[1] 1.332268e-15

20.55 - 1.35 - 19.2

[1] 0

• Loss of significance is characterized by the change in relative error is much larger than change in absolute error.
• Generally, the relative error should be insignificant compared to the absolute error.
• Cancellation flips this around by creating a situation where the change in absolute error is very small, comparatively.
• This hinders our ability to understand the error when performing error analysis (an important step for numerical analysts).

abs(x - xp),  abs((x-xp)/x)

• Formulas can often be rearranged to remedy this, as we saw with the various summation algorithms.
• The quadratic formula provides an excellent example of this.

\[ x = \frac{-b \pm \sqrt{b^2-4ac}}{2a} \]

• Conjugate of radical expression:

\[ \begin{align*} x &= \frac{-b \pm \sqrt{b^2-4ac}}{2a} \\ &= \frac{-b \pm \sqrt{b^2-4ac}}{2a} \times \frac{-b \mp \sqrt{b^2-4ac}}{-b \mp \sqrt{b^2-4ac}}\\ &= \frac{b^2 - (b^2-4ac)}{2a(-b \mp \sqrt{b^2-4ac})}\\ &= \frac{-2c}{b \pm \sqrt{b^2-4ac}} \end{align*} \]

• For the two roots \( x_1 \) and \( x_2 \), we now have two formulas each to choose from.
• Under certain conditions, cancellation error will arise.
• The correct choice will avoid cancellation error.

\[ \begin{align*} x_1 &= \frac{-b + \sqrt{b^2-4ac}}{2a} = \frac{-2c}{b + \sqrt{b^2-4ac}}\\ \\ x_2 &= \frac{-b - \sqrt{b^2-4ac}}{2a} = \frac{-2c}{b - \sqrt{b^2-4ac}} \end{align*} \]

• Consider \( x^2 + 3x + 1 = 0 \)
• Which formula is best for \( x_1 \)? Which is best for \( x_2 \)?

\[ \begin{align*} x_1 &= \frac{-b + \sqrt{b^2-4ac}}{2a} = \frac{-2c}{b + \sqrt{b^2-4ac}}\\ \\ x_2 &= \frac{-b - \sqrt{b^2-4ac}}{2a} = \frac{-2c}{b - \sqrt{b^2-4ac}} \end{align*} \]

• Consider \( x^2 + 3x + 1 = 0 \)
• Answer: Regular formulas will work, since \( b^2 \) not large compared to \( 4ac \), so cancellation error will not occur.
  

\[ \begin{align*} x_1 &= \frac{-b + \sqrt{b^2-4ac}}{2a} = \frac{-2c}{b + \sqrt{b^2-4ac}}\\ \\ x_2 &= \frac{-b - \sqrt{b^2-4ac}}{2a} = \frac{-2c}{b - \sqrt{b^2-4ac}} \end{align*} \]

• Consider \( x^2 + 100x + 1 = 0 \).
• Which formula is best for \( x_1 \)? Which is best for \( x_2 \)?

\[ \begin{align*} x_1 &= \frac{-b + \sqrt{b^2-4ac}}{2a} = \frac{-2c}{b + \sqrt{b^2-4ac}}\\ \\ x_2 &= \frac{-b - \sqrt{b^2-4ac}}{2a} = \frac{-2c}{b - \sqrt{b^2-4ac}} \end{align*} \]

• Consider \( x^2 + 100x + 1 = 0 \)
• Answer: Since \( b > 0 \) and \( b^2 \gg 4ac \), choose second formula for \( x_1 \) and first formula for \( x_2 \) to avoid cancellation error.

\[ \begin{align*} x_1 &= \frac{-b + \sqrt{b^2-4ac}}{2a} = \frac{-2c}{b + \sqrt{b^2-4ac}}\\ \\ x_2 &= \frac{-b - \sqrt{b^2-4ac}}{2a} = \frac{-2c}{b - \sqrt{b^2-4ac}} \end{align*} \]

• Consider \( x^2 - 100x + 1 = 0 \).
• Which formula is best for \( x_1 \)? Which is best for \( x_2 \)?

\[ \begin{align*} x_1 &= \frac{-b + \sqrt{b^2-4ac}}{2a} = \frac{-2c}{b + \sqrt{b^2-4ac}}\\ \\ x_2 &= \frac{-b - \sqrt{b^2-4ac}}{2a} = \frac{-2c}{b - \sqrt{b^2-4ac}} \end{align*} \]

• Consider \( x^2 - 100x + 1 = 0 \)
• Answer: Since \( b < 0 \) and \( b^2 \gg 4ac \), choose first formula for \( x_1 \) and second formula for \( x_2 \) to avoid cancellation error.

\[ \begin{align*} x_1 &= \frac{-b + \sqrt{b^2-4ac}}{2a} = \frac{-2c}{b + \sqrt{b^2-4ac}}\\ \\ x_2 &= \frac{-b - \sqrt{b^2-4ac}}{2a} = \frac{-2c}{b - \sqrt{b^2-4ac}} \end{align*} \]

quadratic <- function (b2 , b1 , b0) {
t1 <- sqrt (b1^2 - 4*b2*b0)
t2 <- 2*b2
x1 <- (-b1 + t1)/t2
x2 <- (-b1 - t1)/t2
return (c(x1,x2))   }

quadratic2 <- function (b2 , b1 , b0) {
t1 <- sqrt (b1^2 - 4*b2*b0)
t2 <- 2*b0
x1 <- t2/(-b1 - t1)
x2 <- t2/(-b1 + t1)
return (c(x1 , x2)) }

• Actual: \( x_1 = 1, ~~ x_2 = 1.000000028975958... \)

b2 <- 94906265.625
b1 <- 189812534.000
b0 <- 94906268.375
print(quadratic(b2, b1, b0), digits = 20)

[1] -1.0000000144879793 -1.0000000144879793

print(quadratic2(b2, b1, b0), digits = 20)

[1] -1.000000014487979 -1.000000014487979

quadraticP3 <- function(a,b,c) {
d <- sqrt(b^2 - 4*a*c)
x1 <- 0 
x2 <- 0
test <- (4*a*c)/b**2 #If b^2 >> 4ac, test = small.
eps = 0.1
if (test > eps) { #Reqular formulas 
x1 <- (-b+d)/(2*a) 
x2 <- (-b-d)/(2*a)
Loop <- 1  }   else {
if (b > 0) 
{ x1 = 2*c/(-b-d)     #New formula
x2 = (-b-d)/(2*a) #Regular formula
Loop <- 2 } else  
{ x1 = (-b+d)/(2*a) #Regular formula
x2 = 2*c/(-b+d)     #New formula
Loop <- 3  } }
return (c(x1,x2,Loop)) }

• Actual: \( x_1 = 1, ~~ x_2 = 1.000000028975958... \)

b2 <- 94906265.625
b1 <- 189812534.000
b0 <- 94906268.375
print(quadraticP3(b2, b1, b0), digits = 20)

[1] -1.0000000144879793 -1.0000000144879793  1.0000000000000000

print(quadratic(b2, b1, b0), digits = 20)

[1] -1.0000000144879793 -1.0000000144879793

\( x^2 + 3x + 1 = 0 \)

a <- 1
b <- 3
c <- 1
print(quadraticP3(a, b, c), digits = 20)

[1] -0.3819660112501051 -2.6180339887498949  1.0000000000000000

print(quadratic(a, b, c), digits = 20)

[1] -0.3819660112501051 -2.6180339887498949

\( x^2 + 100x + 1 = 0 \)

a <- 1
b <- 100
c <- 1
print(quadraticP3(a, b, c), digits = 20)

[1]  -0.010001000200050014 -99.989998999799951207   2.000000000000000000

print(quadratic(a, b, c), digits = 20)

[1]  -0.010001000200048793 -99.989998999799951207

\( x^2 - 100x + 1 = 0 \)

a <- 1
b <- -100
c <- 1
print(quadraticP3(a, b, c), digits = 20)

[1] 99.989998999799951207  0.010001000200050014  3.000000000000000000

print(quadratic2(a, b, c), digits = 20)

[1] 99.989998999812158331  0.010001000200050014

\( x^2 + 10^{15} x + 1 = 0 \)

a <- 1
b <- 10^15
c <- 1
print(quadraticP3(a, b, c), digits = 20)

[1] -1.0000000000000001e-15 -1.0000000000000000e+15  2.0000000000000000e+00

print(quadratic2(a, b, c), digits = 20)

[1] -1.0000000000000001e-15                     Inf

\( 6.972 x^2 + 43.06 x + 0.8431 = 0 \)

a <- 6.972
b <- 43.06
c <- 0.8431
quadraticP3(a, b, c)

[1] -0.01964212 -6.15649098  2.00000000