14

Here I'm not interested in accuracy (see 13614) but rather in raw speed. You'd think that for a complex machine-precision number z, calculating Abs[z]^2 should be faster than calculating Abs[z] because the latter requires a square root whereas the former does not. Yet it's not so:

s = RandomVariate[NormalDistribution[], {10^7, 2}].{1, I};
Developer`PackedArrayQ[s]
(* True *)
Abs[s]^2; // AbsoluteTiming // First
(* 0.083337 *)
Abs[s]; // AbsoluteTiming // First
(* 0.033179 *)

This indicates that Abs[z]^2 is really calculated by summing the squares of real and imaginary parts, taking a square root (for Abs[z]), and then re-squaring (for Abs[z]^2).

Is there a faster way to compute Abs[z]^2? Is there a hidden equivalent to the GSL's gsl_complex_abs2 function? The source code of this GSL function is simply to return Re[z]^2+Im[z]^2; no fancy tricks.

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
Roman
  • 47,322
  • 2
  • 55
  • 121
  • 1
    Here's an even slower way: (Re[#]^2 + Im[#]^2) & /@ s. And even slower still: Total[ReIm[#]^2] & /@ s – bill s May 10 '19 at 14:24

2 Answers2

22

There's Internal`AbsSquare:

s = RandomVariate[NormalDistribution[], {10^7, 2}].{1, I};
foo = Internal`AbsSquare[s]; // AbsoluteTiming // First
murf = Abs[s]^2; // AbsoluteTiming // First
(*
  0.022909
  0.063441
*)

foo == murf
(*  True  *)
Michael E2
  • 235,386
  • 17
  • 334
  • 747
  • 1
    Ah yes precisely what I was looking for, many thanks Michael! Is there a repository of such tricks? – Roman May 10 '19 at 14:25
  • 1
    @Roman I was just looking. I thought there was a post about useful Internal` functions, but I couldn't find it just now. The context contains some useful numerical functions like Log1p and Expm1. Statistics`Library` also contains some nice, well-programmed functions. – Michael E2 May 10 '19 at 14:31
  • 3
    @MichaelE2 https://mathematica.stackexchange.com/questions/805/what-are-some-useful-undocumented-mathematica-functions – Chris K May 10 '19 at 14:31
  • 1
    @ChrisK That must be it! Thanks. – Michael E2 May 10 '19 at 14:32
  • I’m going digging for Internal functions tonight! Is there some large difference in precisely how this function operated faster than the others? Just optimized for it? – CA Trevillian May 11 '19 at 01:53
  • 1
    @CATrevillian I would have thought it was in the MKL (Intel Math Kernel Library), but I didn't find it there. I guess I don't know. – Michael E2 May 11 '19 at 03:10
  • 1
    @CATrevillian it's such a simple operation (at least in the GSL source code) that I'd look for it in the icc/gcc/clang optimization capabilities directly, not in the MKL. – Roman May 11 '19 at 13:10
  • 2
    @CATrevillian the Abs function likely operates by calling hypot which is safer but slower than just taking the square root of AbsSquare. AbsSquare, on the other hand, is likely expressed directly as re*re+im*im in assembly, or a vectorized version thereof: with the Compiler Explorer I find that all x86-64 compilers express this function as mulsd xmm0,xmm0; mulsd xmm1,xmm1; addsd xmm0, xmm1 even at -O3 level, which is a literal translation without frills or optimizations. – Roman May 14 '19 at 11:01
0

for v5.2, s Conjugate[s] is fast too, ref the pic:

enter image description here

infoage
  • 9
  • 3
  • 3
    On my computer, Re[s*Conjugate[s]] is about five to ten times slower than Internal`AbsSquare[s]. What is your $Version and what CPU do you have? – Roman Jul 08 '20 at 12:27
  • 1
    Hi, people here generally like users to post code as Mathematica code instead of just images or TeX, so they can copy-paste it. It makes it convenient for them and more likely they will engage with your posts. You may find this meta Q&A helpful. -- BTW, have you seen RandomVariate[NormalDistribution[], {10^7, 2}]? It's much faster on my machine. Ditto for RandomComplex[]. – Michael E2 Jul 08 '20 at 12:48
  • @Roman Re[] is unnecessary, though it's very fast. My version is very old, it's v5.2. So there's no Internal`AbsSquare[]. – infoage Jul 08 '20 at 18:14
  • @MichaelE2 Thanks, man. My version is v5.2. This code is so simple that I had no motivation to paste text version at that moment. Sorry. – infoage Jul 08 '20 at 18:18
  • 1
    Maybe it's worth adding the version info to your answer. It turns out that I don't have the Statistics`NormalDistribution package (in V12.1.1), I suppose because it's been replaced by top-level statistics functions some versions ago. – Michael E2 Jul 08 '20 at 19:08