Units of computer storage¶

• Humans use decimal digits (why?)
• Computers use binary digits (why?)
  • Bit = binary digit (coined by statistician John Tukey).
  • byte = 8 bits.
  • KB = kilobyte = $10^3$ bytes; KiB = kibibyte = $2^{10}$ bytes.
  • MB = megabyte = $10^6$ bytes; MiB = mebibyte = $2^{20}$ bytes.
  • GB = gigabyte = $10^9$ bytes. Typical RAM size.
  • TB = terabyte = $10^{12}$ bytes. Typical hard drive size.
  • PB = petabyte = $10^{15}$ bytes.

• (cont'd)
  • EB = exabyte = $10^{18}$ bytes.
  • ZB = zetabyte = $10^{21}$ bytes. Size of all healthcare data in 2018 is ~33 ZB.

Julia function Base.summarysize shows the amount of memory (in bytes) used by an object.

varinfo() function prints all variables in workspace and their sizes.

Storage of Characters¶

• Plain text files are stored in the form of characters: .jl, .r, .c, .cpp, .ipynb, .html, .tex, ...
• ASCII (American Code for Information Interchange): 7 bits, only $2^7=128$ characters.

• Extended ASCII: 8 bits, $2^8=256$ characters.

• Unicode: UTF-8, UTF-16 and UTF-32 support many more characters including foreign characters; last 7 digits conform to ASCII.

  • UTF-8 is the current dominant character encoding on internet.

• Julia supports the full range of UTF-8 characters. You can type many Unicode math symbols by typing the backslashed LaTeX symbol name followed by tab.

See https://docs.julialang.org/en/v1/manual/unicode-input/

Integers: fixed-point number system¶

• Fixed-point number system is a computer model for integers $\mathbb{Z}$.
  • Remember that computer memory is finite whereas the cardinality of $\mathbb{Z}$ is (countably) infinite.
  • Any representation of numbers in computer has to be an approximation.

• Number of bits M to represent integers and method of representing negative numbers vary from system to system.
  • The integer type in R has $M=32$ (packages such as 'bit64' support 64 bit integers).
    • https://www.r-bloggers.com/r-in-a-64-bit-world/
  • C has (unsigned) char, int, short, long (and long long), whose sizes depend on the machine.
  • Matlab has (u)int8, (u)int16, (u)int32, (u)int64.

• Julia has even more integer types:

Unsigned integers¶

• Model for $\mathbb{N} \cup \{0\}$.
• For unsigned integers, the range is $[0,2^M-1]$.
• Julia functions typemin(T) and typemax(T) give the lowest and highest representable number of a type T respectively

Signed integers¶

• Model of $\mathbb{Z}$. Can do subtraction.

• First bit ("most significant bit" or MSB) is the sign bit.

  • 0 for nonnegative numbers
  • 1 for negative numbers

• Two's complement representation for negative numbers

  • Set the sign bit to 1
  • Negate (0->1, 1->0) the remaining bits
  • Add to 1 to the result
  • Two's complement representation of a negative integer $x$ is the same as the unsigned integer $2^M - x$.

• Two's complement representation respects modular arithmetic nicely. Addition of any two signed integers are just bitwise addition, possibly modulo $2^M$
  • $M=4$ case:

• Range of representable integers by $M$-bit signed integer is $[-2^{M-1},2^{M-1}-1]$:

BigInt¶

Julia BigInt type is arbitrary precision.

Overflow in integer arithmetic¶

R reports NA for integer overflow. Julia outputs the result according to modular arithmetic.

Reals: floating-point number system¶

Floating-point number system is a computer model for $\mathbb{R}$.

• Most computer systems adopt the IEEE 754 standard, established in 1985, for floating-point arithmetics.

For the history, see an interview with William Kahan.

• In the scientific notation, a real number is represented as

Humans use the base $b=10$ and digits $d_i=0, 1, \dotsc, 9$.
In computer, the base is $b=2$ and the digits $d_i$ are 0 or 1. The exponent $e$ is between $e_{\min}$ and $e_{\max}$.

• Normalized vs denormalized numbers. For example, decimal number 18 is

or, equivalently, $$ +0.1001 \times 2^5 \quad (\text{denormalized}).$$

• In the floating-number system, computer stores
  • sign bit
  • the fraction (or mantissa, or significand) of the normalized representation (except some special cases)
  • the actual exponent plus a bias
    • Because $e_{\min}$ is likely negative and $e_{\max}$ is positive, the exponent requires a sign. IEEE 754 handles the sign of the exponent by subtracting the bias from the value of the exponent evaluated as an unsigned integer.

• Julia provides Float16 (half precision, implemented in software), Float32 (single precision), Float64 (double precision), and BigFloat (arbitrary precision).

Half precision (Float16)¶

{width="500" fig-align="center"}

::: {style="font-size: 50%;"}

https://en.wikipedia.org/wiki/Half-precision_floating-point_format

:::

• In Julia, Float16 is the type for half precision numbers.

• MSB is the sign bit.

• 10 significand bits (fraction=mantissa), hence $p=11$ (why?)

• 5 exponent bits: $e_{\max}=15$, $e_{\min}=-14$, bias=15 = $01111_2$ for encoding:

  • $e_{\min} = \mathbf{00001_2} - 01111_2 = -14_{10}$
  • $e_{\max} = \mathbf{11110_2} - 01111_2 = 15_{10}$

• $e_{\text{min}}-1$ and $e_{\text{max}}+1$ are reserved for special numbers.

• range of magnitude: $10^{\pm 4}$ in decimal because $\log_{10} (2^{15}) \approx 4$.

• Precision: $\log_{10}2^{11} \approx 3.311$ decimal digits. $$ (value) = (-1)^{b_{15}}\times 2^{(\sum_{j=1}^5 b_{15-j}2^{5-j}) - 15} \times \left( 1 + \sum_{i=1}^{10}\frac{b_{10-i}}{2^i}\right) $$

Single precision (Float32)¶

{width="800" fig-align="center"}

::: {style="font-size: 50%;"}

https://en.wikipedia.org/wiki/Single-precision_floating-point_format

:::

• In Julia, Float32 is the type for single precision numbers.

• MSB is the sign bit.

• 23 significand bits ($p=24$).

• 8 exponent bits: $e_{\max}=127$, $e_{\min}=-126$, bias=127.

• $e_{\text{min}}-1$ and $e_{\text{max}}+1$ are reserved for special numbers.

• range of magnitude: $10^{\pm 38}$ in decimal because $\log_{10} (2^{127}) \approx 38$.

• precision: $\log_{10}(2^{24}) \approx 7.225$ decimal digits.

Double precision (Float64)¶

{width="800" fig-align="center"}

::: {style="font-size: 50%;"}

https://en.wikipedia.org/wiki/Double-precision_floating-point_format

:::

• Double precision (64 bits = 8 bytes) numbers are the dominant data type in scientific computing.

• In Julia, Float64 is the type for double precision numbers.

• MSB is the sign bit.

• 52 significand bits ($p=53$).

• 11 exponent bits: $e_{\max}=1023$, $e_{\min}=-1022$, bias=1023.

• $e_{\text{min}}-1$ and $e_{\text{max}}+1$ are reserved for special numbers.

• range of magnitude: $10^{\pm 308}$ in decimal because $\log_{10} (2^{1023}) \approx 308$.

• precision to the $\log_{10}(2^{53}) \approx 15.95$ decimal point.

Special floating-point numbers¶

• Exponent $e_{\max}+1$ plus a zero mantissa means $\pm \infty$.

• Exponent $e_{\max}+1$ plus a nonzero mantissa means NaN. NaN could be produced from 0 / 0, 0 * Inf, ...

• In general NaN ≠ NaN bitwise. Test whether a number is NaN by isnan function.

• Exponent $e_{\min}-1$ with a zero mantissa represents the real number 0 ("exact zero").

• Exponent $e_{\min}-1$ with a nonzero mantissa are for numbers less than $b^{e_{\min}}$ (subnormal numbers).
• Numbers are denormalized in the range $(0,b^{e_{\min}})$ -- gradual underflow.
• For example, in half-precision, $e_{\min}=-14$ but $2^{-24}$ is represented by $0.0000000001_2 \times 2^{-14}$.
• If floating point numbers are all normalized, the spacing between 0 and $b^{e_{\min}}$ is just $b^{e_{\min}}$, whereas the spacing between $b^{e_{\min}}$ and $b^{e_{\min}+1}$ is $b^{e_{\min}-p+1}$. With subnormal numbers, the spacing between 0 and $b^{e_{\min}}$ can be $b^{e_{\min}-p+1}$, which is consistent with the spacing just above $b^{e_{\min}}$.

Rounding¶

• Rounding is necessary whenever a number has more than $p$ significand bits. Most computer systems use the default IEEE 754 round to nearest mode (also called ties to even mode). Julia offers several rounding modes, the default being RoundNearest.

• "Round to nearest, ties to even" rule: rounds to the nearest value; if the number falls midway, it is rounded to the nearest value with an even least significant digit (i.e., zero; default for IEEE 754 binary floating point numbers)

• For example, the number 1/10 cannot be represented accurately as a (binary) floating point number:

Errors¶

Rounding (more fundamentally, finite precision) incurs errors in floating porint computation. If a real number $x$ is represented by a floating point number $[x]$, then

• Absolute error: $| [x] - x |$

• Relative error: $\frac{| [x] - x |}{|x|}$

(if $x \neq 0$).

Of course, we want to ensure that the error after a computation is small.

Machine epsilons¶

• Floating-point numbers do not occur uniformly over the real number line

{width="800" fig-align="center"}

::: {style="font-size: 50%;"}

Source: What you never wanted to know about floating point but will be forced to find out

:::

• Same number of representible numbers in $(2^i, 2^{i+1}]$ as in $(2^{i+1}, 2^{i+2}]$. Within an interval, they are uniformly distributed.

• Machine epsilons are the spacings of numbers around 1:
  • $\epsilon_{\max}$ = (smallest positive floating point number that added to 1 will give a result different from 1) = $\frac{1}{2^p} + \frac{1}{2^{2p-1}}$
  • $\epsilon_{\min}$ = (smallest positive floating point number that subtracted from 1 will give a result different from 1) = $\frac{1}{2^{p+1}} + \frac{1}{2^{2p}}$.
  • That is, $1 + \epsilon_{\max}$ is the number in the "midway" between 1 and the floating point number right next to it, etc.

{width="1000" fig-align="center"}

::: {style="font-size: 50%;"}

Source: Computational Statistics, James Gentle, Springer, New York, 2009.

Caution: the definition of $\epsilon_{\max}$ and $\epsilon_{\min}$ in this book is different from the lecture note.

:::

• Any real number in the interval $\left[1 - \frac{1}{2^{p+1}}, 1 + \frac{1}{2^p}\right]=[1.111\dotsb1|1_2 \times 2^{-1}, 1.000\dotsb 0|1 \times 2^0]$ is represented by a floating point number $1 = 1.00\dotsc 0_2 \times 2^0$ (assuming the "ties to even" rule: consider $p=2$ case).

• Adding $\frac{1}{2^p}$ to 1 won't reach the next representable floating point number $1.00\dotsc 01_2 \times 2^0 = 1 + \frac{1}{2^{p-1}}$. Hence $\epsilon_{\max} > \frac{1}{2^p} = 1.00\dotsc 0_2 \times 2^{-p}$.

• Adding the floating point number next to $\frac{1}{2^p} = 1.00\dotsc 0_2 \times 2^{-p}$ to 1 WILL result in $1.00\dotsc 01_2 \times 2^0 = 1 + \frac{1}{2^{p-1}}$, hence $\epsilon_{\max} = 1.00\dotsb 01_2 \times 2^{-p} = \frac{1}{2^p} + \frac{1}{2^{p+p-1}}$.

• Subtracting $\frac{1}{2^{p+1}}$ from 1 results in $1-\frac{1}{2^{p+1}} = \frac{1}{2} + \frac{1}{2^2} + \dotsb + \frac{1}{2^p} + \frac{1}{2^{p+1}}$, which is represented by the floating point number $1.00\dotsb 0_2 \times 2^{0} = 1$ by the "ties to even" rule. Hence $\epsilon_{\min} > \frac{1}{2^{p+1}}$.

• The smallest positive floating point number larger than $\frac{1}{2^{p+1}}$ is $\frac{1}{2^{p+1}} + \frac{1}{2^{2p}}=1.00\dotsc 1_2 \times 2^{-p-1}$. Thus $\epsilon_{\min} = \frac{1}{2^{p+1}} + \frac{1}{2^{2p}}$.

Machine precision¶

• Machine epsilon is often called the machine precision.

• If a positive real number $x \in \mathbb{R}$ is represented by $[x]$ in the floating point arithmetic, then $$ [x] = \left( 1 + \sum_{i=1}^{p-1}\frac{d_{i+1}}{2^i}\right) \times 2^e. $$

Thus $x - \frac{2^e}{2^p} \le [x] \le x + \frac{2^e}{2^p}$, and $$ \begin{split} \frac{| x - [x] |}{|x|} &\le \frac{2^e}{2^p|x|} \le \frac{2^e}{2^p}\frac{1}{[x]-2^e/2^p} \\ &\le \frac{2^e}{2^p}\frac{1}{2^e(1-1/2^p)} \quad (\because [x] \ge 2^e) \\ &\le \frac{2^e}{2^p}\frac{1}{2^e}(1 + \frac{1}{2^{p-1}}) \\ &= \frac{1}{2^p} + \frac{1}{2^{2p-1}} = \epsilon_{\max}. \end{split} $$ That is, the relative error of the floating point representation $[x]$ of real number $x$ is bounded by $\epsilon_{\max}$.

In Julia, eps(x) gives the distance between consecutive representable floating point values at x, i.e.,

eps(x) == max(x-prevfloat(x), nextfloat(x)-x)

which is roughly twice the $\epsilon_{\max}$.

• In R, the variable .Machine contains numerical characteristics of the machine. double.eps and double.neg.eps are roughly twice our $\epsilon_{\max}$ and $\epsilon_{\min}$, respectively.

Overflow and underflow¶

• For double precision, the range is $\pm 10^{\pm 308}$. In most situations, underflow (magnitude of result is less than $10^{-308}$) is preferred over overflow (magnitude of result is larger than $10^{+308}$). Overflow produces $\pm \inf$. Underflow yields zeros or subnormal numbers.

• Example: the logit link function is

The former expression can easily lead to Inf / Inf = NaN, while the latter expression leads to gradual underflow.

• floatmin and floatmax functions gives largest and smallest non-subnormal number represented by the given floating point type.

• BigFloat in Julia offers arbitrary precision.

Catastrophic cancellation¶

The result of computation is just the digits that represented the rounding.

• Scenario 1 (benign cancellation): Addition or subtraction of two numbers of widely different magnitudes: $a+b$ or $a-b$ where $a \gg b$ or $a \ll b$. We lose the precision in the number of smaller magnitude. Consider $$ a = x.xxx ... \times 2^{30}, \quad b = y.yyy... \times 2^{-30} $$

What happens when computer calculates $a+b$? We get $a+b=a$!

Analysis: suppose we want to compute $x + y$, $x, y > 0$. Let the relative error of $x$ and $y$ be $$ \delta_x = \frac{[x] - x}{x}, \quad \delta_y = \frac{[y] - y}{y} . $$ What the computer actually calculates is $[x] + [y]$, which in turn is represented by $[ [x] + [y] ]$. The relative error of this representation is $$ \delta_{\text{sum}} = \frac{[[x]+[y]] - ([x]+[y])}{[x]+[y]} . $$ Recall that $|\delta_x|, |\delta_y|, |\delta_{\text{sum}}| \le \epsilon_{\max}$.

We want to find a bound of the relative error of $[[x]+[y]]$ with respect to $x+y$. Since $|[x]+[y]| = |x(1+\delta_x) + y(1+\delta_y)| \le |x+y|(1+\epsilon_{\max})$, we have $$ \small \begin{split} | [[x]+[y]]-(x+y) | &= | [[x]+[y]] - [x] - [y] + [x] - x + [y] - y | \\ &\le | [[x]+[y]] - [x] - [y] | + | [x] - x | + | [y] - y | \\ &= |\delta_{\text{sum}}([x]+[y])| + |\delta_x x| + |\delta_y y| \\ &\le \epsilon_{\max}(x+y)(1+\epsilon_{\max}) + \epsilon_{\max}x + \epsilon_{\max}y \\ &\approx 2\epsilon_{\max}(x+y) \end{split} $$ because $\epsilon_{\max}^2 \approx 0$.

Thus $$ \frac{| [[x]+[y]]-(x+y) |}{|x+y|} \le 2\epsilon_{\max} $$ approximately.

• Scenario 2 (catastrophic cancellation): Subtraction of two nearly equal numbers eliminates significant digits. $a-b$ where $a \approx b$. Consider $$ a = x.xxxxxxxxxx1ssss \quad b = x.xxxxxxxxxx0tttt $$

The result is $1.vvvvu...u$ where $u$ are unassigned digits.

Analysis: Let $$ [x] = 1 + \sum_{i=1}^{p-2}\frac{d_{i+1}}{2^i} + \frac{1}{2^{p-1}}, \quad [y] = 1 + \sum_{i=1}^{p-2}\frac{d_{i+1}}{2^i} + \frac{0}{2^{p-1}} . $$

• $[x]-[y] = \frac{1}{2^{p-1}} = [[x]-[y]]$.

• The true difference $x-y$ may lie anywhere in $(0, \frac{1}{2^{p-2}}]$.

• Relative error

is unbounded -- no guarantee of any significant digit!

• Implications for numerical computation
  • Rule 1: add small numbers together before adding larger ones
  • Rule 2: add numbers of like magnitude together (paring). When all numbers are of same sign and similar magnitude, add in pairs so each stage the summands are of similar magnitude

• Rule 3: avoid substraction of two numbers that are nearly equal
  • Example: in solving quadratic equation $ax^2 + bx + c = 0$ with $b > 0$ and $b^2 - 4ac > 0$, the roots are
  $$ x = \frac{-b \pm \sqrt{b^2 - 4ac}}{2a} . $$ If $b \approx \sqrt{b^2 - 4ac}$, then a scenario 2 occurs. A resolution is to compute only one root $x_1 = \frac{-b - \sqrt{b^2 - 4ac}}{2a}$ by the formula, and find the other root $x_2$ using the relation $x_1x_2 = c/a$.

Algebraic laws¶

Floating-point numbers may violate many algebraic laws we are familiar with, such associative and distributive laws. See the example in the Machine Epsilon section and HW1.

Coditioning¶

Consider solving a linear system $Ax=b$:

whose solution is $(x_1, x_2) = (1.000, 1.000)$.

If we perturb $b$ by 0.001, i.e., solve

then the solution changes to $(x_1, x_2) = (0.000, 3.000)$.

If we instead perturb $A$ by 0.001, i.e., solve

then this time the solution changes to $(x_1, x_2) = (2.000, -1.000)$.

In other words, an input perturbation of order $10^{-3}$ yield an output perturbation of order $10^0$. Thats 3 orders of magnutide of relative change!

Floating point representation $[x]$ of a real number $x$ may introduce such input perturbation easily. The perturbation of output of a problem with respect to the input is called conditioning.

• Abstractly, a problem can be viewed as function $f: \mathcal{X} \to \mathcal{Y}$ where $\mathcal{X}$ is a normed vector space of data and $\mathcal{Y}$ is a normed vector space of solutions.
  • The problem of solving $Ax=b$ for fixed $b$ is $f: A \mapsto A^{-1}b$ with $\mathcal{X}=\{M\in\mathbb{R}^{n\times n}: M \text{ is invertible} \}$ and $\mathcal{Y} = \mathbb{R}^n$.
  • The combination of a problem $f$ with a given data $x$ is called a problem instance, or simply problem unless no confusion occurs.

• A well-conditioned problem (instance) is one such that all small perturbations of $x$ lead to only small changes in $f(x)$.

• An ill-conditioned problem is one such that some small perturbation of $x$ leads to a large change in $f(x)$.

• The (relative) condition number $\kappa=\kappa(\theta)$ of a problem is defined by $$ \kappa = \lim_{\delta\to 0}\sup_{\|\delta \theta\|\le \delta}\frac{\|\delta f\|/\|f(\theta)\|}{\|\delta \theta\|/\|\theta\|} , $$

where $\delta f = f(\theta + \delta \theta) - f(\theta)$.

• For the problem of solving $Ax=b$ for fixed $b$, $f: A \mapsto A^{-1}b$, it can be shown that the condition number of $f$ is $$ \kappa = \|A\|\|A^{-1}\| =: \kappa(A) , $$

where $\|A\|$ is the matrix norm induced by vector norm $\|\cdot\|$, i.e., $$ \|A\| = \sup_{x\neq 0} \frac{\|Ax\|}{\|x\|}. $$

• To see this, recall $f(A) = A^{-1}b = x$. Since $b = Ax$ and the LHS is fixed, if a perturbation of the input $A \to A + \delta A$ yields a perturbation of the output $x \to x + \delta x$ or $f \to f + \delta f$, then we have $$ 0 = (\delta A) x + A \delta x + (\delta A) \delta x $$

or $$ \delta x = -A^{-1}(\delta A)x + o(\Vert \delta A \Vert) . $$ That is, $\delta f = f(A + \delta A) - f(A) = - A^{-1}(\delta A)A^{-1}b + o(\Vert \delta A \Vert)$.

• Then $$ \begin{aligned} \kappa &= \lim_{\delta\to 0}\sup_{\|\delta A\|\le \delta}\frac{(\Vert A^{-1}(\delta A)A^{-1}b\Vert + o(\Vert \delta A \Vert)\Vert) / \Vert A^{-1}b \Vert}{\Vert \delta A \Vert / \Vert A \Vert} \\ &= \lim_{\delta\to 0}\sup_{\|\delta A\|\le \delta} \frac{\Vert A^{-1}(\delta A)A^{-1}b\Vert / \Vert A^{-1}b \Vert}{\Vert \delta A \Vert / \Vert A \Vert} + o(1) \\ &= \lim_{\delta\to 0}\sup_{\|\delta A\|\le \delta} \frac{\Vert A^{-1}(\delta A)A^{-1}b\Vert / \Vert A^{-1}b \Vert}{\Vert \delta A \Vert / \Vert A \Vert} . \end{aligned} $$

• From the property of a matrix norm $\|AB\| \leq \|A\|\|B\|$, $$ \frac{\Vert A^{-1}(\delta A)A^{-1}b\Vert / \Vert A^{-1}b \Vert}{\Vert \delta A \Vert / \Vert A \Vert} \le \frac{\Vert A^{-1} \Vert \Vert \delta A \Vert \Vert A^{-1}b\Vert / \Vert A^{-1}b \Vert}{\Vert \delta A \Vert / \Vert A \Vert}

% = \Vert A^{-1} \Vert \Vert A \Vert , $$ hence $\kappa \le \Vert A^{-1} \Vert \Vert A \Vert$.

• It can be shown that for any invertible $A$ and $b$, there exits a perturbation $\delta A$ such that

• It follows that $\kappa \ge \Vert A^{-1} \Vert \Vert A \Vert$.

• If Euclidean norm is used, then

the ratio of the maximum and minimum singular values of $A$.

In the above problem, the condition number is matrix $A$ (w.r.t. Euclidean norm) is

Further readings¶

• Section II.2, Computational Statistics by James Gentle (2009).

• Sections 1.5 and 2.2, Applied Numerical Linear Algebra by James W. Demmel (1997).

• What every computer scientist should know about floating-point arithmetic by David Goldberg (1991).