Around workaround

Question

I'm testing Around for error propagation (v12.2), but it seems to be inconsistent. Error results depend on the form. Am I missing something here?

Edit Same results in wolfram cloud, which has v12.3.

ClearAll[get$amplitude] ;
get$amplitude[frequency_, signal_] := Block[
  {c, s},
  c = Cos[2*Pi*frequency*Range[Length[signal]]] ;
  s = Sin[2*Pi*frequency*Range[Length[signal]]] ;
  c = 2.0*Dot[signal, c] ;
  s = 2.0*Dot[signal, s] ;
  Divide[{Norm[{c, s}], Sqrt[c^2 + s^2], Sqrt[c*c + s*s]}, Length[signal]]
] ;

frequency = 0.12 ;
signal = 0.5*Sin[2*Pi*frequency*Range[256]] ;
signal = Around[signal, 0.1] ;
TableForm[Map[FullForm, get$amplitude[frequency, signal]]]
(* Around[0.501507305064027`,0.008852097698515098`] *)
(* Around[0.501507305064027`,0.008852097698515098`] *)
(* Around[0.501507305064027`,0.006259378310345857`] *)

Here is a torch implementation, I was expecting results to match exactly. What can be the reason behind this mismatch? Looks like Around doesn't use automatic differentiation, I thought it uses Series under the hood.

import torch
from math import pi

pi = torch.tensor(pi, dtype=torch.float64)
length = 256
time = torch.linspace(start=1.0, end=length, steps=length, dtype=torch.float64)
signal = 0.5 * torch.sin(2.0 * pi * 0.12 * time)
signal.requires_grad_(True) ;

def get_amplitude(length, frequency, time, signal):
  c = torch.cos(2.0 * pi * frequency * time)
  s = torch.sin(2.0 * pi * frequency * time)
  c = 2.0 * torch.dot(c, signal)
  s = 2.0 * torch.dot(s, signal)
  return torch.sqrt(c*c + s*s) / length

amplitude = get_amplitude(256, 0.12, time, signal)
amplitude.backward()

sigma = 0.1*0.1 + torch.zeros(length, dtype=torch.float64)
sigma = torch.diag(sigma)
print(amplitude.item())
print(torch.sqrt(torch.dot(signal.grad, sigma @ signal.grad)).item())
# 0.501507305064027
# 0.008852294777182603

Answer 1

It seems to be because (since Around represents statistical uncertainty) it tries to figure out which expressions are correlated and which ones are uncorrelated. In general, Around[x, a]^2 is different than Around[x, a] * Around[x, a], since in the former, we're drawing from the same distribution, but in the latter, we're drawing from two uncorrelated distributions.

Around[x, a]^2

(* Out: Around[x^2, 2 Abs[a x]] *)

Around[x, a] * Around[x, a]

(* Out (upon simplifying): Around[x^2, Sqrt[2] Abs[a x] *)

Likewise for addition:

2 Around[x, a]

(* Out: Around[2 x, 2 a] *)

Around[x, a] + Around[x, a]

(* Out (upon simplifying): Around[2 x, Sqrt[2] Abs[a]] *)

It's a bit clunky, but to avoid this, it seems like you need to use AroundReplace:

AroundReplace[z + z, {z -> Around[x, a]}]

(* Out (upon simplifying): Around[2 x, 2 Abs[a]] *)

AroundReplace[2 z, {z -> Around[x, a]}]

(* Out (upon simplifying): Around[2 x, 2 Abs[a]] *)

Note that these actually seem to agree "anyway" because z + z gets simplified to 2 z first. But the documentation assures us that AroundReplace "does what we want" and considers two instances of the same replaced variable as coming from the same distribution, and tells us to contrast the following as an example:

expr = p q/(p + q);
rules = {p -> Around[p0, \[Delta]p], q -> Around[q0, \[Delta]q]};
$Assumptions = (p | q | p0 | q0 | \[Delta]p | \[Delta]q) \[Element] 
   Reals;
AroundReplace[expr, rules] // Simplify
ReplaceAll[expr, rules] // Simplify

(Note that AroundReplace supports VectorAround when we want to consider the uncertainties of the entries in a vector to be correlated with each other.)

So, a workaround for your particular function might consist of giving it this definition for a list of Around's:

ClearAll[get$amplitude]
get$amplitude[frequency_, signal : {___Around}] :=
 Block[
  {c, s, c0, s0},
  c = Cos[2*Pi*frequency*Range[Length[signal]]] ;
  s = Sin[2*Pi*frequency*Range[Length[signal]]] ;
  c = 2.0*Dot[signal, c] ;
  s = 2.0*Dot[signal, s] ;
  AroundReplace[
   Unevaluated[
    {Norm[{c0, s0}], Sqrt[c0^2 + s0^2], Sqrt[c0*c0 + s0*s0]}/Length[signal]], {c0 -> c, s0 -> s}]
]

One problem is that Mathematica apparently doesn't know that Derivative[1][Abs][x] is 1 for x > 0:

frequency = 0.12;
signal = 0.5*Sin[2*Pi*frequency*Range[256]];
signal = Around[signal, 0.1];
TableForm[Map[FullForm, get$amplitude[frequency, signal]]]

(* Out: *)

Around[0.501507305064027`,Power[Plus[Times[1.734130926267436`*^-9,Power[Derivative[1][Abs][0.6057819465394019`],2]],Times[0.00007835789953313003`,Power[Derivative[1][Abs][128.38444091337783`],2]]],Rational[1,2]]]
Around[0.501507305064027`,0.008852097698515098`]
Around[0.501507305064027`,0.008852097698515098`]

@Lukas Lang identified a better way to resolve the issue in a comment: append /. {Abs -> RealAbs} to the result (as Mathematica assumes `Abs` is a function on the compelxes, which is not holomorphic, and so doesn't have a simple derivative that behaves nicely), and it'll work itself out.