GitHub

Using automatic differentiation

This tutorial shows how to compute the power spectrum $P(k; \theta)$ and its (logarithmic) derivatives

\[\frac{\partial \lg P}{\partial \lg \theta_i}\]

using automatic differentiation with ForwardDiff.jl. This technique can also differentiate any other quantity.

1. Wrap the evaluation

First, we must decide which parameters $\theta$ the power spectrum $P(k; \theta)$ should be considered a function of. To do so, let us write a small wrapper function that calculates the power spectrum as a function of the parameters $(T_{\gamma 0}, \Omega_{c0}, \Omega_{b0}, N_\textrm{eff}, h, Y_p)$, following the Getting started tutorial:

using SymBoltz
M = ΛCDM(K = nothing)
pars = [M.γ.T₀, M.c.Ω₀, M.b.Ω₀, M.ν.Neff, M.g.h, M.b.YHe, M.h.m_eV, M.I.ln_As1e10, M.I.ns]
prob0 = CosmologyProblem(M, Dict(pars .=> NaN))

probgen = parameter_updater(prob0, pars)
P(k, θ) = spectrum_matter(probgen(θ), k; verbose = true, ptopts = (reltol = 1e-3,))
P (generic function with 1 method)

It is now easy to evaluate the power spectrum:

using Unitful, UnitfulAstro
θ = [2.7, 0.27, 0.05, 3.0, 0.7, 0.25, 0.02, 3.0, 0.95]
ks = 10 .^ range(-3, 0, length=100) / u"Mpc"
Ps = P(ks, θ)
100-element Vector{Unitful.Quantity{Float64, 𝐋^3, Unitful.FreeUnits{(Mpc^3,), 𝐋^3, nothing}}}:
 15091.151632943904 Mpc^3
 16046.243096773032 Mpc^3
 17053.12847431816 Mpc^3
 18113.60727346069 Mpc^3
 19228.362545205317 Mpc^3
 20398.792341501085 Mpc^3
 21625.665217421014 Mpc^3
 22909.643159974992 Mpc^3
 24250.70444047198 Mpc^3
 25648.66517770079 Mpc^3
     ⋮
   235.06475418563764 Mpc^3
   200.79387973248714 Mpc^3
   171.33410232783206 Mpc^3
   146.03469347296254 Mpc^3
   124.34637341877104 Mpc^3
   105.77058472703024 Mpc^3
    89.88876857175515 Mpc^3
    76.3115448431274 Mpc^3
    64.72921782278705 Mpc^3

This can be plotted with

using Plots
plot(log10.(ks/u"1/Mpc"), log10.(Ps/u"Mpc^3"); xlabel = "lg(k/Mpc⁻¹)", ylabel = "lg(P/Mpc³)", label = nothing)
Example block output

2. Calculate the derivatives

To get $\partial \lg P / \partial \lg \theta$, we can simply pass the wrapper function P(k, θ) through ForwardDiff.jacobian:

using ForwardDiff
lgP(lgθ) = log10.(P(ks, 10 .^ lgθ) / u"Mpc^3") # in log-space
dlgP_dlgθs = ForwardDiff.jacobian(lgP, log10.(θ))
100×9 Matrix{Float64}:
 -0.14814   -1.23995  -0.229644  -0.0098543  …  -0.00456279  3.0  -3.71642
 -0.162172  -1.23403  -0.228547  -0.0108805     -0.00453246  3.0  -3.65014
 -0.177443  -1.22755  -0.227343  -0.0119989     -0.00452908  3.0  -3.58385
 -0.194046  -1.22053  -0.226036  -0.0132178     -0.00455693  3.0  -3.51756
 -0.212126  -1.21298  -0.224641  -0.0145449     -0.0046202   3.0  -3.45128
 -0.231731  -1.20492  -0.223151  -0.0159878  …  -0.00472139  3.0  -3.38499
 -0.252872  -1.19613  -0.221522  -0.0175504     -0.00486125  3.0  -3.3187
 -0.275673  -1.18664  -0.219768  -0.0192415     -0.00503967  3.0  -3.25242
 -0.300201  -1.17644  -0.21788   -0.0210683     -0.00525571  3.0  -3.18613
 -0.326555  -1.16559  -0.215867  -0.0230406     -0.00550771  3.0  -3.11984
  ⋮                                          ⋱                    
 -5.9886     1.88713  -0.366274  -0.383656      -0.0253159   3.0   2.31565
 -6.0209     1.90387  -0.366372  -0.384704      -0.0253174   3.0   2.38194
 -6.05065    1.91964  -0.36668   -0.385783      -0.0253189   3.0   2.44823
 -6.07976    1.93503  -0.367155  -0.386891      -0.0253196   3.0   2.51451
 -6.10945    1.94966  -0.36726   -0.388046   …  -0.0253206   3.0   2.5808
 -6.13667    1.96434  -0.367485  -0.388952      -0.0253211   3.0   2.64709
 -6.16466    1.97805  -0.36763   -0.389994      -0.0253219   3.0   2.71337
 -6.19172    1.99126  -0.36773   -0.391052      -0.0253226   3.0   2.77966
 -6.21711    2.00426  -0.368018  -0.391939      -0.0253234   3.0   2.84595

The matrix element dlgP_dlgθs[i, j] now contains $\partial \lg P(k_i) / \partial \lg \theta_j$. We can plot them all at once:

plot(
    log10.(ks/u"1/Mpc"), dlgP_dlgθs;
    xlabel = "lg(k/Mpc⁻¹)", ylabel = "∂ lg(P) / ∂ lg(θᵢ)",
    labels = "θᵢ=" .* ["Tγ0" "Ωc0" "Ωb0" "Neff" "h" "YHe" "mh" "ln(10¹⁰As)" "ns"]
)
Example block output

Get values and derivatives together

The above example showed how to calculate the power spectrum values and their derivatives through two separate calls. If you need both, it is faster to calculate them simultaneously with the package DiffResults.jl:

using DiffResults

# Following DiffResults documentation:
Pres = DiffResults.JacobianResult(ks/u"1/Mpc", θ) # allocate buffer for values+derivatives for a function with θ-sized input and ks-sized output
Pres = ForwardDiff.jacobian!(Pres, lgP, log10.(θ)) # evaluate values+derivatives of lgP(log10.(θ)) and store the results in Pres
lgPs = DiffResults.value(Pres) # extract values
dlgP_dlgθs = DiffResults.jacobian(Pres) # extract derivatives

p1 = plot(
   log10.(ks/u"1/Mpc"), lgPs;
   ylabel = "lg(P/Mpc³)", label = nothing
)
p2 = plot(
   log10.(ks/u"1/Mpc"), dlgP_dlgθs;
   xlabel = "lg(k/Mpc⁻¹)", ylabel = "∂ lg(P) / ∂ lg(θᵢ)",
   labels = "θᵢ=" .* ["Tγ0" "Ωc0" "Ωb0" "Neff" "h" "YHe" "mh" "ln(10¹⁰As)" "ns"]
)
plot(p1, p2, layout=(2, 1), size = (600, 600))
Example block output

General approach

The technique shown here can be used to calculate the derivative of any SymBoltz.jl output quantity:

  1. Write a wrapper function output(input) that calculates the desired output quantities from the desired input quantities.
  2. Use ForwardDiff.derivative(output, input) (scalar-to-scalar), ForwardDiff.gradient(output, input) (vector-to-scalar) or ForwardDiff.jacobian(output, input) (vector-to-vector) to evaluate the derivative of output at the values input. Or use the similar functions in DiffResults to calculate the value and derivatives simultaneously.