How to simulate a low precision calculation?

Questions
,
1

I recently encountered some unexpected behavior in my numerical simulation, and I suspect that may come from the round off error in the calculation. I plan to test it using low precision calculations, but I can’t figure out how to do that.

For example, a Fourier transform

f[x_] = Cos[2 \[Pi] x];
N0 = 10;
ls = Most@Array[f, N0, {0, 1}];

using the default machine precision I get

Fourier[ls]

{-6.4763*10^-17 + 0. I, 1.5 + 5.55112*10^-17 I, 1.60247*10^-17 + 2.77556*10^-17 I, 3.23815*10^-17 - 8.01234*10^-18 I, -1.60247*10^-17 - 2.77556*10^-17 I, -1.60247*10^-17 + 2.77556*10^-17 I, 3.23815*10^-17 + 8.01234*10^-18 I, 1.60247*10^-17 - 2.77556*10^-17 I,1.5 - 5.55112*10^-17 I}

but using fixed lower precision, I get

Block[{$MaxPrecision = 8, $MinPrecision = 8}, 
 Fourier[SetPrecision[ls, 8]]]

{-4.3981121*10^-35, 1.5000000, 1.8230938*10^-34, 2.1989581*10^-35, 1.8230938*10^-34, 1.8230938*10^-34, 2.1989581*10^-35, 1.8230938*10^-34, 1.5000000}

It appears that fix 8 precision has smaller errors than the default machine precision. I would expect using fixed precision at 8 would give larger errors than machine precision. So why does it behave like this, and what are the correct ways to simulate the low precision calculations?

2

Let’s write our own DFT function for a comparison
I’ve taken the code for this from the documentation for Fourier. Let’s fist compute the values symbolically

(* Calculate the DFT symbolically *)
x = ls; n = Length[x]; u /: Subscript[u, r_] := x[[r]]; 
FullSimplify[Table[(1/Sqrt[n])*Sum[Subscript[u, r]*E^(2*Pi*I*(r - 1)*((s - 1)/n)), {r, 1, n}], {s, 1, n}]]

{0, 3/2, 0, 0, 0, 0, 0, 0, 3/2}

We can run this same code with different precisions and compare to the “correct” results above.

(*With precision = 8*)
Block[{$MaxPrecision = 8, $MinPrecision = 8},
x = SetPrecision[ls, 8]; n = Length[x]; u /: Subscript[u, r_] := x[[r]]; 
Table[(1/Sqrt[n])*Sum[Subscript[u, r]*E^(2*Pi*I*(r - 1)*((s - 1)/n)), {r, 1, n}], {s, 1, n}]
]

{-4.4041855*10^-35, 1.5000000, 1.8255231*10^-34, 2.2020928*10^-35, 
 1.8255231*10^-34, 1.8255231*10^-34, 2.2020928*10^-35, 1.8255231*10^-34, 1.5000000}

(*With boring old machine precision*)
x = N[ls]; n = Length[x]; u /: Subscript[u, r_] := x[[r]]; 
Table[(1/Sqrt[n])*Sum[Subscript[u, r]*E^(2*Pi*I*(r - 1)*((s - 1)/n)), {r, 1, n}], {s, 1, n}]

{-1.11022*10^-16, 1.5 + 0. I, -9.25186*10^-17 + 0. I, 
 0. + 0. I, -1.11022*10^-16 + 0. I, -1.11022*10^-16 + 0. I, 
 0. + 0. I, -9.25186*10^-17 + 0. I, 1.5 + 0. I}

So the built-in compiled version and a version we built in Mathematica do the same thing. If we run with a bunch of different precisions, we can see how the results differ based on the accuracy of the numbers used (X axis is precision used. Y axis is how off the actual results the calculations end up being. Closer to 0 is better ):

ex1.png
ex1.png360×203 2.87 KB
ex2.png
ex2.png461×255 3.81 KB

Higher precision is generally good, but answering “Why is not always the case that that is true” can be difficult. After all, if I asked “why does 0.1+0.2 give 0.30000000000000004 in floating point computation?” what answer could be given except to walk step-by-step through how the addition plays out in floating point in all the painful detail?

Higher precision usually means better answers, yes, but there’s nothing that guarantees this that I know of.

3

It appears that fix 8 precision has smaller errors than the default machine precision. I would expect using fixed precision at 8 would give larger errors than machine precision. So why does it behave like this

As mentioned on ref/Precision, machine precision is considered lower than any explicit precision like 8. The machine precision computation is carried out directly using the computer’s native floating point arithmetic, hence without any error tracking. This has the advantage of being very fast. On the other hand, when using precision 8, the arbitrary-precision routines are used instead and the error is tracked to keep the result is as accurate as possible given the input precision and operations performed.

4

Thanks for the answers. I tried to make the example in Sean’s answer into a fortran (at the end), and by compile with
ifort test.f90or ifort -r8 test.f90
I can set it to use single (fixed precision at 8) or double precision (fixed precision at 16).And the results are expected, that double precision gives smaller errors.
In Sean’s example, it shows that even fixed precision at 1 would give a smaller error compared to machine precision. This seems to suggested that Mathematica uses higher precision in the calculation even we set set $MaxPrecision or $MinPrecision to a small number (for example 1). This is very confusing for me, because I though setting $MinPrecision and $MaxPrecision t0 8 would minic the single precision calculation, and setting $MinPrecision and $MaxPrecision to 16 would minic the double precision calculation. So could you explain the differences between a calculation by setting setting $MinPrecision and $MaxPrecision to 16 and a machine precision calculation? And is there ways to minic the round off error in low precision calculations?
----------------------
test program

program main

 implicit none



 real :: Pi=4.0*atan(1.0)

 integer,parameter :: n=9

 integer r,s

 real :: x(n)

 complex :: ft(n)

 complex :: ci=(0.,1.)





 do r=1,n 

  x(r)=cos(2.*(r-1.)*Pi/n)

 end do



 do s=1,n

  ft(s)=(0.,0.)

  do r=1,n

   ft(s) = ft(s)+ x(r)*exp(2.*Pi*ci*(r - 1.)*(s - 1.)/n)

  end do

  ft(s)=ft(s)/sqrt(real(n))

 end do



 do r=1,n

  write (*,*) ft(r)

 end do



end program

single precision, compiled with

ifort test.f90

results:

 (1.9868216E-08,0.0000000E+00)

 (1.500000,3.2782555E-07)

 (2.5828680E-07,-4.1723251E-07)

 (7.9472862E-08,-1.7881393E-07)

 (1.1920929E-07,-1.7881393E-07)

 (1.3907750E-07,-4.5696893E-07)

 (-4.5696893E-07,-2.5828680E-07)

 (-6.7551929E-07,-1.3907750E-07)

 (1.500000,1.5695890E-06)

double precision, compiled with

ifort -r8 test.f90

results:

 (-2.960594732333751E-016,0.000000000000000E+000)

 (1.50000000000000,2.220446049250313E-016)

 (-1.850371707708594E-017,-1.110223024625157E-016)

 (-1.665334536937735E-016,-3.700743415417188E-017)

 (1.850371707708594E-016,1.295260195396016E-016)

 (-6.661338147750939E-016,-3.515706244646329E-016)

 (4.625929269271486E-016,-7.401486830834377E-017)

 (4.903485025427775E-016,7.771561172376096E-016)

 (1.50000000000000,-2.090920029710711E-015)
5

The result is decidedly not of high precision. If you check the InputForm you will see it is 8 as expected. What it is, however, is high accuracy on the zeros.

{-4.404185501080837196522866226126013806214324`8.*^-35,
1.49999999999999999999999999999999982511`8.,
1.8253664026743891567492525318255642007953505`8.*^-34,
2.202092750540418598261433113063006903107162`8.*^-35,
1.8251704869492521088312577424026938265157147`8.*^-34,
1.8251704869492521088312577424026938265157147`8.*^-34,
2.202092750540418598261433113063006903107162`8.*^-35,
1.8253664026743891567492525318255642007953505`8.*^-34,
1.49999999999999999999999999999999982511`8.}

So we know these are zero to 35 or so places (to the extent that working in fixed precision does not mess things up), but we do not know the actual digit values at all well.