7
$\begingroup$

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
$\endgroup$
3
  • 2
    $\begingroup$ Around was reimplemented in V12.3. This information may not be valuable for you, but it may also have fixed the issues you are having. $\endgroup$ Commented Sep 2, 2021 at 1:30
  • $\begingroup$ @CATrevillian, could you test WM code in V12.3, do the results match? $\endgroup$
    – I.M.
    Commented Sep 2, 2021 at 1:36
  • 1
    $\begingroup$ I'm running 12.3.1 and my results match what you show for WM. $\endgroup$
    – Mark R
    Commented Sep 2, 2021 at 1:54

1 Answer 1

7
$\begingroup$

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.

$\endgroup$
3
  • 1
    $\begingroup$ Thanks, this explains the difference in WM code. Looks like one should be careful using Around in expressions. Still need to wrap my head around which case is 'correct'. For torch approximate formula is used sigma_f^2 = grad(f) var(x) grad(f) , this might explain the difference, since in WM it's more like using variance-covariance $\endgroup$
    – I.M.
    Commented Sep 2, 2021 at 4:56
  • 1
    $\begingroup$ Regarding the issues with Abs: The issue is that the complex derivative of Abs is not well-defined, hence the issues. You can simply use /.Abs->RealAbs to replace it with RealAbs, which does not have these problems $\endgroup$
    – Lukas Lang
    Commented Sep 2, 2021 at 11:59
  • $\begingroup$ Gotcha. Interesting that it doesn't have a downvalue defined on Reals (without 0), though. I think that would be useful! $\endgroup$
    – thorimur
    Commented Sep 2, 2021 at 19:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.