# Experiments with Variable Format Half Precision

For the past month I have been working with the variable format 16-bit floating point arithmetic that I described in this post. It has been frustrating work. I have found that the limited precision and limited range of half precision make it barely usable for the kind of experiments with matrix computation that I like to do. In this post I will describe a few of these experiments.

I am not saying that half precision floating point is not useful. Organizations including Intel and Google are adopting both the IEEE 754 standard *fp16* and the new *bfloat16* formats for work in machine learning and image processing, link. I am just not recommending either format for general engineering and scientific computing.

### Contents

#### vfp16

I have revised the `help` entry for my variable format half precision object. Here is the updated version.

```
doc vfp16
```

Here is the list of functions currently available. As you can see, this is far from a complete set. The workhorse functions `mldivide` (backslash), `lu`, `qr` and `svd` are available, but I have not yet done anything about complex arithmetic, so `eig` is not on the list.

methods(vfp16)

Methods for class vfp16: abs fix lu plus single tiny binary ge max power size tril ctranspose gt minus precision sqrt triu diag hex mldivide qr subsasgn uint disp horzcat mrdivide realmax subsindex uminus display isfinite mtimes realmin subsref vertcat double le ne rem sum vfp16 eps log2 norm round svd eq lt packed sign times

#### vfp16 Anatomy

This is a repeat of the chart I had in my previous post, except now I am using the `vfp16` object itself to compute the parameters.

`p`is the precision, the number of bits in the fraction.`eps`is the distance from 1 to the next larger`vfp16`number.`realmax`is the largest`vfp16`number.`realmin`is the smallest normalized`vfp16`number.`tiny`is the smallest subnormal`vfp16`number.

disp(' p eps realmax realmin tiny') for p = 1:15 vfp16('precision',p) x = vfp16(1); y = [eps(x) realmax(x) realmin(x) tiny(x)]; fprintf('%4d %12.4g %12.4g %12.4g %12.4g\n',p,double(y)) end

p eps realmax realmin tiny 1 0.5 Inf 0 0 2 0.25 Inf 0 0 3 0.125 Inf 0 0 4 0.0625 1.742e+308 2.225e-308 1.391e-309 5 0.03125 1.32e+154 2.983e-154 9.323e-156 6 0.01563 1.149e+77 3.454e-77 5.398e-79 7 0.007813 3.39e+38 1.175e-38 9.184e-41 8 0.003906 1.841e+19 2.168e-19 8.47e-22 9 0.001953 4.291e+09 9.313e-10 1.819e-12 10 0.0009766 6.55e+04 6.104e-05 5.96e-08 11 0.0004883 255.9 0.01563 7.629e-06 12 0.0002441 16 0.25 6.104e-05 13 0.0001221 4 1 0.0001221 14 0 2 Inf 0.0001221 15 0 4 -0 0.0001221

For precision equal to 1, 2 or 3 there are 12 or more bits devoted to the exponent and so `realmax`, `realmin` and `tiny` cannot be expressed in MATLAB's conventional double precision arithmetic This effectively means that `vfp16` cannot take advantage of the extended range.

On the other hand, for precision equal to 11 or greater, `eps` is less than `realmin`, so underflow rather than roundoff dominates accuracy. Especially if gradual underflow is `off`, computation with these precisions behaves badly. Consequently I keep `p` in the range `4 <= p <= 10`.

#### Gradual Underflow

Before we had the IEEE standard, code that was portable across different architectures had to discover machine accuracy. One way this was done involved the following computation. Let's use the standard `fp16` which has `p = 10` and `'subnormals'` turned `'on'`.

format short vfp16('fp16') a = vfp16(4)/vfp16(3) b = a-1 c = b+b+b e = 1-c

a = 1.3330 b = 0.3330 c = 0.9990 e = 9.7656e-04

With exact arithmetic, `b` would be 1/3, `c` would be 1, and `e` would be zero. But, unless we have a ternary computer, there must be a roundoff error in the division by 3. That's the only roundoff error. So `b` is not exactly 1/3, `c` is not exactly 1, and `e` is equal to the spacing of floating numbers near, but larger than, 1.0. In other words, `e` is equal to `eps(1)`.

Let's do that again with precision equal to 11 and gradual underflow turned off.

vfp16('precision',11,'subnormals','off') a = vfp16(4)/vfp16(3) b = a-1 c = b+b+b e = 1-c

a = 1.3335 b = 0.3335 c = 1.0005 e = 0

Now the division by 3 rounds up instead of down and `c` is slightly larger than one, but the computation of `e` underflows and is flushed to zero. This is a simple example, but it demonstrates the importance of gradual underflow when we have a limited exponent range.

#### Fused Multiply Add

Fused multiply add makes only a modest difference in the roundoff error in the residual `R` from LU decomposition of a random 10-by-10 matrix, as measured by number of nonzeros in `R`.

vfp16('subnormals','on') form = {'fp16','bfloat16'}; fma = {'off','on'}; rng(17) for j = 1:2 vfp16(form{j}) A = vfp16(randn(10,10)); for k = 1:2 vfp16('fma',fma{k}) [L,U,p] = lu(A); R = L*U - A(p,:); subplot(2,2,2*(j-1)+k) spy(double(R)) title([form{j} ', fma ' fma{k}]) end end vfp16('subnormals','on')

#### QR with Column Pivoting, Frank Matrix

The Frank matrix $F_n$ is an interesting test matrix. Here is $F_6$.

```
F6 = gallery('frank',6)
```

F6 = 6 5 4 3 2 1 5 5 4 3 2 1 0 4 4 3 2 1 0 0 3 3 2 1 0 0 0 2 2 1 0 0 0 0 1 1

We are going to compute the QR decomposition with column pivoting and look at $r_n$, the last diagonal element of $R$. It turns out that $F_n$ is close to a matrix of rank $n-1$ and this is forecast by $r_n$ decaying to zero exponentially as $n$ increases. Here is the $R$ from the QR decomposition (with column pivoting) of $F_6$; you can see the beginnings of exponential smallness in the last diagonal element.

[Q,R,E] = qr(F6); R

R = -8.1240 -6.7700 -5.1698 -6.8931 -3.4466 -1.7233 0 -3.8944 0.5136 0.6847 0.3424 0.1712 0 0 3.6068 2.4986 2.7742 1.3871 0 0 0 -1.6652 0.5547 0.2773 0 0 0 0 1.0001 1.0000 0 0 0 0 0 0.0053

When the QR decomposition is computed by our half precision arithmetic this exponential decay cannot be tracked for long. In the following logarithmic plot the black line is $r_n$. The other lines are the values obtained with `vfp16` arithmetic of varying precision. Although $r_{15}$ is less than $10^{-10}$, none of the other values is less than $10^{-4}$.

The broken lines are caused by exact zeros which cannot be plotted on a logarithmic plot.

Rn = zeros(1,15); rn = zeros(10,15); for n = 1:15 A = gallery('frank',n); [~,R,E] = qr(A); Rn(n) = R(n,n); for p = 4:10 vfp16('precision',p) X = vfp16(A); [~,r,e] = qrp(X); rn(p,n) = r(n,n); end end xqrp_print xqrp_plot

n 2 3 4 ... 13 14 15 Rn 4.47e-01 2.43e-01 -9.41e-02 ... -1.41e-09 1.09e-10 -7.82e-12 rn p 4 -0.4375 -0.2813 0.0547 ... -0.0547 -0.1875 -0.1406 5 -0.4844 -0.2969 0.0977 ... -0.0332 0.0000 0.0332 6 -0.4297 -0.2266 0.1074 ... -0.0308 0.0154 -0.0308 7 -0.4453 -0.2402 0.1016 ... 0.0000 0.0078 -0.0078 8 -0.4434 -0.2402 0.0879 ... -0.0118 0.0000 0.0000 9 -0.4473 -0.2412 0.0879 ... -0.0020 0.0000 -0.0039 10 -0.4468 -0.2429 0.0952 ... -0.0039 -0.0039 0.0000

#### SVD, Hilbert Matrix

Now for two players familiar to readers of this blog, the Hilbert matrix and the Singular Value Decomposition. Here is $H_5$.

```
format rat
H5 = hilb(5)
```

H5 = 1 1/2 1/3 1/4 1/5 1/2 1/3 1/4 1/5 1/6 1/3 1/4 1/5 1/6 1/7 1/4 1/5 1/6 1/7 1/8 1/5 1/6 1/7 1/8 1/9

The matrix is symmetric, positive definite, so its singular values are equal to its eigenvalues. See how the singular values $H_6$ are decaying rapidly.

format short g H6 = hilb(6); sigma6 = svd(H6)

sigma6 = 1.6189 0.24236 0.016322 0.00061575 1.2571e-05 1.0828e-07

This experiment is like the QR of the Frank matrix, except the logarithmic plot of the smallest singular value is almost perfectly linear and has a steeper slope. You can see that the red line from the IEEE standard with p = 10 does the best job of tracking $\sigma_n$.

sn = zeros(10,13); for n = 2:13 A = hilb(n); S = svd(A); Sn(n) = S(n); for p = 4:10 vfp16('precision',p) X = vfp16(A); s = svd(X); sn(p,n) = s(n); end end xsvd_hilbert_print xsvd_hilbert_plot

n 2 3 4 ... 11 12 13 Sn 6.57e-02 2.69e-03 9.67e-05 ... 3.40e-15 1.11e-16 3.79e-18 sn p 4 6.99e-02 4.81e-03 1.89e-03 ... 4.44e-04 2.11e-04 3.13e-04 5 6.91e-02 1.90e-03 2.24e-03 ... 2.20e-04 2.05e-04 5.45e-05 6 6.56e-02 2.36e-03 0.00e+00 ... 1.18e-04 2.44e-05 1.33e-04 7 6.72e-02 3.09e-03 9.89e-05 ... 1.22e-04 1.70e-05 6.97e-05 8 6.46e-02 2.79e-03 1.74e-04 ... 4.46e-05 2.36e-05 2.43e-05 9 6.59e-02 2.91e-03 3.12e-05 ... 8.70e-06 6.78e-06 1.27e-05 10 6.55e-02 2.61e-03 3.67e-05 ... 9.33e-07 8.26e-06 3.45e-06

#### SVD, Parter Matrix

Now for something different. I discussed what the gallery calls the Parter matrix five years ago. Here is $P_5$.

format rat P5 = gallery('parter',5)

P5 = 2 -2 -2/3 -2/5 -2/7 2/3 2 -2 -2/3 -2/5 2/5 2/3 2 -2 -2/3 2/7 2/5 2/3 2 -2 2/9 2/7 2/5 2/3 2

Most of the singular value are converging to $\pi$. I'm not sure where the smallest singular value is going, but it is certainly not going to zero exponentially. Now the computations with `vfp16` with any precision in the range `4 <= p <= 10` does a good job.

Sn = zeros(1,13); sn = zeros(10,13); for n = 2:13 A = gallery('parter',n); S = svd (A); Sn(n) = S(n); for p = 4:10 vfp16('precision',p) X = vfp16(A); s = svd(X); sn(p,n) = s(n); end end xsvd_parter_print xsvd_parter_plot

n 2 3 4 ... 11 12 13 Sn 1.7370 1.6053 1.5217 ... 1.2814 1.2640 1.2485 sn p 4 1.7430 1.5907 1.4846 ... 1.3287 1.3204 1.3339 5 1.7430 1.5502 1.4985 ... 1.2493 1.2373 1.2223 6 1.7227 1.5997 1.5093 ... 1.2799 1.2582 1.2521 7 1.7327 1.6050 1.5283 ... 1.2874 1.2639 1.2488 8 1.7407 1.6038 1.5229 ... 1.2843 1.2641 1.2478 9 1.7380 1.6040 1.5216 ... 1.2799 1.2628 1.2474 10 1.7361 1.6051 1.5206 ... 1.2822 1.2646 1.2479

#### Cleve's Laboratory

I will update Cleve_s Laboratory in the MATLAB Central File Exchange to Version 4.30 to include my changes to `@vfp16`.

**Category:**- Eigenvalues,
- Matrices,
- Numerical Analysis,
- Precision

## Comments

To leave a comment, please click here to sign in to your MathWorks Account or create a new one.