# Is Abs[z]^2 a bad way to calculate the square modulus of z?

For a numerical quantity z, Abs[z] returns the square root of the sum of the squares of the real and imaginary parts of z. If we then square it, do we lose numerical accuracy because we are squaring after taking a square root?
Is it better to use Re[z]^2+Im[z]^2?

-
There are two spelling mistakes in your post, but I don't dare to correct them ;-) – Sjoerd C. de Vries Oct 25 '12 at 7:31
Please do feel free to edit my posts, this was never an issue for me! (I corrected the mistakes, thank you :). – Joe Oct 25 '12 at 7:37
I don't know if this is a good figure of merit, but Accuracy[Abs[z]^2] and Accuracy[Re[z]^2+Im[z]^2] yield the same. Looking forward to answers from people who know more about numerics than I do. – sebhofer Oct 25 '12 at 9:29
Food for thought: compare Abs[(3 + 4 I) (1 + N[2^-52])]^2 // InputForm, (Re[#]^2 + Im[#]^2) &[(3 + 4 I) (1 + N[2^-52])] // InputForm, and N[(Re[#]^2 + Im[#]^2) &[(3 + 4 I) (1 + 2^-52)], 20]... do something similar using the factor 1 - 2^-52. – J. M. Oct 25 '12 at 9:32

Let's find out.

Here is a function to compare the results of evaluating f on a machine approximation of an exact argument z to its exact value:

test[f_] := Function[{z}, f[N[z]] -  f[z]];


We expect to find problems around the 52nd significant binary digit or so. Let's magnify these differences, then, by about $2^{60}$ and use logarithmic axes to study them. (In order to look at the logs, we neglect the signs of the errors and also add a small starting value to them to avoid taking the logarithm of $0$. Even so, the y-axes more or less indicate the number of binary digits lost out of $60$. Thus, we should expect log errors to be typically around $8$ when z has a modulus near $1$.) The errors will likely vary with the argument of the complex number z and will tend to be directly proportional to its modulus.

The following code--which is readily extensible to alternative procedures by modifying its first line--explores three ways of evaluating a squared modulus, using Norm, Abs, and the sum of squares of real and imaginary parts. To do this, it plots the errors at equally spaced points along semi circles of radii $3/4$, $1$, and $4/3$ (colored blue, red, and gold, respectively).

functions = {Norm[#]^2 &, Abs[#]^2 &, (Re[#]^2 + Im[#]^2) &};
TableForm[
{Table[ListPlot[
Table[Log2[1 + 2^60 Abs[test[g][\[Rho] Exp[I Pi \[Theta]]]]],
{\[Rho], {3/4, 1, 4/3}}, {\[Theta], 0, 1, 1/Prime[80]}],
PlotStyle -> PointSize[0.01], DataRange -> {0, Pi}, PlotRange -> {Full, {-1, 11}},
ImageSize -> 300],
{g, functions}]},
TableHeadings -> {{}, Style[TraditionalForm[#], 16] & /@ functions},
TableAlignments -> Center]


As expected, the error tends to be directly proportional to the modulus of the argument, as indicated by the approximately equal spacings between the three colors.

When it errs, Abs (shown in the middle column) tends to make relatively large errors, but it appears to err less often than the other two methods, which have similar (but not identical) error patterns. Moreover, the other two procedures have strong angular patterns in their errors: they tend to err more frequently for directions oriented along the axes, but may have slightly smaller errors there.

Abs probably uses a square root function that is optimized for accuracy: it works most of the time, but not all the time; and when it fails, it's just as bad--and sometimes slightly worse--than the alternatives.

Abs[#]^2 & behaves much like # Conjugate[#] & (not shown). This is a puzzle: the latter clearly does not optimize a square root, because it never uses one. This suggests that in the multiplication, Mathematica is exploiting the polar representation used to create points along the circles. Tracing the evaluation confirms this. Thus we learn that the accuracy of a calculation may depend not only on the value of the argument, but also on its form, even when the argument is only going to be evaluated to machine precision.

-
Awesome analysis! By slightly adjusting your code, I plotted the difference in the error between the functions Abs[z]^2 and Re[z]^2+Im[z]^2 when they are applied to the same argument, and indeed for most of the points you used Abs[z]^2 has a smaller error! This is really weird and counterintuitive to me, but at least now I know. – Joe Oct 25 '12 at 16:43
I would appreciate a bit of expounding on the last sentence of your answer. – Mr.Wizard Oct 26 '12 at 0:44
@Mr. W: First, let's establish an equality for two ways of expressing a number I will call $\theta$: Exp[I 1/19] - (Cos[1/19] + I Sin[1/19]) // FullSimplify returns 0. Now consider comparing these expressions to their machine-precision versions: With[{q = 1/19}, a = #~Join~N[#] & @ {Exp[I q], Cos[q] + I Sin[q]}]; Outer[(Abs[#1] - Abs[#2]) 2^60 &, a, a] // FullSimplify // MatrixForm. The result differs depending on how $\theta$ is expressed. – whuber Oct 26 '12 at 12:14