Precision loss in FourierDCT

I have high precision data that I want to do Fourier transform with, but it looks like the precision is lost dramatically.

data = CloudGet["https://www.wolframcloud.com/objects/2135516f-ec65-44ad-8fe1-e87f6ab0c15e"];

Precision[data]
(* 50. *)

Block[{$MaxExtraPrecision = ∞}, FourierDCT[data]] // Precision (* 38.4857 *) Block[{$MaxExtraPrecision = ∞}, FourierDCT[data, 1]] // Precision
(* 14.9341 *)


So why do I get a large amount of precision loss, especially in the type-I DCT? Is this a bug?

• It might be worthwhile to export the data to a file, upload in another program. i.e. NumPy and SciPy have FFT tools. – wgwz Feb 19 '15 at 16:57
• Check the Accuracy of the coefficients. I believe the absolute error is stable, but if you're approximating an analytic function, then you can expect the Fourier coefficents to vanish quickly. Therefore the precision gets worse. – Michael E2 Dec 10 '17 at 22:53

I've noticed the fall-off of precision of Fourier coefficients with all the Chebyshev approximation work I've done over the last few years. I left a comment a few months ago hinting at what I thought was going on, namely that each coefficient has about the same absolute error.

I'm not an expert on the FFT, but I've glanced at Kaneko 1970, Ramos 1971, Schatzman 1996, which discuss floating-point round-off error in the FFT. It seems traditional to discuss the error in terms of the RMS error, $$E_\text{RMS}(v) = \left(\sum_{j=1}^N |v_j|^2 \Big/ N \right)^{1/2}$$ and the relative error in these terms, $${E_\text{RMS}\left(v_{\text{approx}} - v_{\text{exact}} \right) \over E_\text{RMS}\left(v_{\text{exact}} \right)}\,.$$ This relative error is equivalent to the relative error in the 2-norm, since the factors of $N^{-1/2}$ cancel. This is different from Precision[v], which is essentially equivalent to the infinity norm, since Precision returns the precision of the element of v with the worst precision.

I have two things to explain I think: why the precision loss is to be expected and how to understand the errors in the results of Fourier and its ilk. The second informs the first.

On the errors of FourierDCT[v, 2]

I'll focus mainly on one function, the DCT-II. I suggested in my comment that Accuracy could be used to understand the error. The literature seems to prefer the RMS error, in part probably because it is easier to prove bounds on it. I can see from my own examples that they are correlated (some are given below), and from the examples it is not hard to understand why. Accuracy[] is the negative logarithm of the absolute error. In the examples at least, the Accuracy[] of each term stays roughly constant, which means the absolute error estimates of the terms are roughly the same. Hence the root-mean-square will be approximately equal to the ordinary mean of the absolute errors and also not far from largest or smallest errors, except for an occasional outlier.

On the other hand, Precision[] is the negative logarithm of the relative error. If the absolute errors of the coefficients $c_j$ stay roughly constant, say $|\Delta c_j| = e$, then Precision[cj = FourierDCT[data, 2]] will be the negative logarithm of the maximum of $${e \over |c_0|}, {e \over |c_0|}, \dots, {e \over |c_n|}\,.$$ If $c_j \rightarrow 0$, then the precision gets worse. The range of the coefficients determines most of the precision loss. If $${\mathop{\text{Max}} |c_j| \over \mathop{\text{Min}} |c_j|} = 10^p\,,$$ then I would expect that the Precision[cj] to be $p$ digits less than Precision[data], minus any loss due to round-off error. This is just what we see:

Block[{$MaxExtraPrecision = ∞}, dct2 = FourierDCT[data, 2]]; First@Differences@MinMax@RealExponent@dct2 + Precision[dct2] Precision[data] - % (* 49.3543 0.645744 <-- loss due to round-off error *)  One can see that the Accuracy[] is roughly constant and the magnitude of the coefficients is decreasing, which lead to decreasing in Precision[]: ListPlot[Accuracy /@ dct2, PlotRange -> {45.8, 50.3}] ListPlot[RealExponent@dct2]  The error of DCT-I Here's the same analysis applied to the DCT-I: Block[{$MaxExtraPrecision = ∞}, dct1 = FourierDCT[data, 1]];
First@Differences@MinMax@RealExponent@dct1 + Precision[dct1]
Precision[data] - %
(*
49.089
0.910985  <-- loss due to round-off error
*)

ListPlot[Accuracy /@ dct1, PlotRange -> {45.8, 50.3}]
ListPlot[RealExponent@dct1]


Other examples

Here are some examples in which the data is generated by the same process as for Chebyshev interpolation. The series coefficients are computed at a working precisions of 100, to be considered accurate; 50, to compare with the OP's example; and MachinePrecision, since different code is used for it than for arbitrary precision. The "error" is the relative RMS error mentioned above. Things to note:

• The machine-precision computation (column 1) is highly accurate. So is the 50-digit computation (column 2). (Compared to the 100-digit computation.)
• The Max, Median, and Min precision loss of the coefficients cj are in column 3. The max is equal to 100 - Precision[cj]; it should be compare to column 4, which is the precision loss due to the range of the coefficients cj. The difference in all cases is less than 1 and is due to numerical error of the FFT, which I've been calling round-off error. (The range of the exponents of cj may be observed in the plots in column 6.)
• The Accuracy range (column 5) shows that the accuracy of the coefficients is nearly constant in all examples.
• The last column shows what the Precision[] loss is when computing the Chebyshev interpolant. It is not particularly important; I was using it as a consistency check.

(* SetPrecision[x, 3] is for output formatting only *)
ClearAll[errornorm];
errornorm[approx_, exact_] := SetPrecision[Norm[approx - exact]/Norm[exact], 3];

(* abscissae for Chebyshev interpolation *)
With[{n = 2^10},
abscissae = Sin[Pi/2 Range[n, -n, -2]/n]];

ds = Dataset@Association@Table[
Block[{data, cj, eval, nonzero, prec0 = 100},
data = f@abscissae;
cj = FourierDCT[N[data, prec0], 1];
eval = Sqrt[2/(Length@abscissae - 1)] cj;
eval[[{1, -1}]] /= 2;
nonzero = DeleteCases[cj, c_ /; c == 0];
f[\[FormalX]] -> <|
"MP error" -> StandardForm@
errornorm[SetPrecision[FourierDCT[N[data], 1], prec0], cj],
"WP-50 error" -> StandardForm@
errornorm[SetPrecision[FourierDCT[N[data, 50], 1], prec0], cj],
"c(j) prec loss" -> Column[prec0 -
Through[{Min, Median, Max}[Precision /@ nonzero]]],
"Δlog10|c(j)|" -> First@Differences@MinMax@RealExponent@nonzero,
"acc range" -> Column@MinMax[Accuracy /@ cj],
"coeff plot" ->
ListPlot[RealExponent@cj, PlotRange -> {-101., 1.8}],
"eval prec loss" -> Column[
prec0 - Through[
{Min, Median, Max}[
Precision /@ (Cos[
Outer[Times,
RandomReal[{0, Pi}, 100, WorkingPrecision -> prec0],
Range[0, Length[eval] - 1]]].eval)]]
]
|>
],
{f, {Exp[(# - 1)] &, Sqrt[(101/100 + #)/2] &, Sqrt[(1 + #)/2] &,
1/(64 #^2 + 1) &, Exp@Sin[25 #] &}}]

• Thanks for linking to those papers, I haven't encountered them yet! – J. M.'s torpor Apr 1 '18 at 0:07
• @J.M. You're welcome. I'm not sure they're the best refs, though. – Michael E2 Apr 1 '18 at 2:17

(This got too long for a comment.)

I got very curious about the lowered precision in the results for the DCT, so as a check, I used the relation between FFT and DCT-I/II to generate the transformed sequences:

Block[{$MaxExtraPrecision = ∞}, FourierDCT[data, 1]] // Precision 14.9341 Block[{$MaxExtraPrecision = ∞},
Re[Take[Fourier[Join[data, data[[-2 ;; 2 ;; -1]]]], Length[data]]]] // Precision
15.0717


and

Block[{$MaxExtraPrecision = ∞}, FourierDCT[data, 2]] // Precision 38.4857 Block[{$MaxExtraPrecision = ∞},
Re[Take[Exp[2 π I Range[0, 2 Length[data] - 1]/(4 Length[data])]
Fourier[Join[data, Reverse[data]]]/Sqrt[2], Length[data]]]] // Precision
38.2994


I would not have found this of note if the precision results of using direct DCT and the FFT-based method were different (in fact, I would have expected an even lower precision result), but they are quite close. I am now curious as to how FFTW (which IIRC is used under the hood) would deal with this if directly fed this sequence.