# Automatic Differentiable Numerical Renormalization Group

Jonas B. Rigo<EMAIL_ADDRESS>School of Physics, University College
Dublin, Belfield, Dublin 4, Ireland  Centre for Quantum Engineering, Science,
and Technology, University College Dublin, Belfield, Dublin 4, Ireland  Andrew
K. Mitchell<EMAIL_ADDRESS>School of Physics, University College
Dublin, Belfield, Dublin 4, Ireland  Centre for Quantum Engineering, Science,
and Technology, University College Dublin, Belfield, Dublin 4, Ireland

###### Abstract

Machine learning techniques have recently gained prominence in physics,
yielding a host of new results and insights. One key concept is that of
backpropagation, which computes the exact gradient of any output of a program
with respect to any input. This is achieved efficiently within the
differentiable programming paradigm, which utilizes automatic differentiation
(AD) of each step of a computer program and the chain rule. A classic
application is in training neural networks. Here, we apply this methodology
instead to the numerical renormalization group (NRG), a powerful technique in
computational quantum many-body physics. We demonstrate how derivatives of NRG
outputs with respect to Hamiltonian parameters can be accurately and
efficiently obtained. Physical properties can be calculated using this
differentiable NRG scheme—for example, thermodynamic observables from
derivatives of the free energy. Susceptibilities can be computed by adding
source terms to the Hamiltonian, but still evaluated with AD at precisely zero
field. As an outlook, we briefly discuss the derivatives of dynamical
quantities and a possible route to the vertex.

##  I Introduction

The Numerical Renormalization Group (NRG) Wilson (1975); Bulla  _et al._
(2008) is a standard tool in computational physics for solving a certain class
of quantum many-body problem, known as quantum impurity models. These comprise
a few interacting quantum degrees of freedom, coupled to one or more non-
interacting baths. The most famous example is the Kondo model, which involves
a single spin-$\tfrac{1}{2}$ ‘impurity’ coupled to a single fermionic bath
Hewson (1997). The low-energy physics is non-perturbative and controlled by a
single emergent energy scale (the Kondo temperature, $T_{\rm K}$), below which
the impurity is dynamically screened by a spatially-extended entanglement
‘Kondo cloud’ of surrounding conduction electrons Mitchell  _et al._ (2011).
Generalized quantum impurity models describe the scattering from magnetic
impurities in metals Hewson (1997); Costi  _et al._ (2009), semiconductor
quantum dot devices Goldhaber-Gordon  _et al._ (1998); Keller  _et al._
(2014); Iftikhar  _et al._ (2015); *mitchell2016universality;
*iftikhar2018tunable; Pouse  _et al._ (2021), single molecule transistors
Liang  _et al._ (2002); Mitchell  _et al._ (2017), and are the local effective
models within dynamical mean field theory (DMFT) for correlated materials
Kotliar  _et al._ (2006); Georges  _et al._ (1996); Stadler  _et al._ (2015).

NRG is a tensor network method Weichselbaum  _et al._ (2009) that exploits the
renormalization group structure of such quantum impurity models Wilson (1975).
At its core, the method involves the iterative numerical diagonalization of a
set of renormalized effective Hamiltonians. There is no statistical element,
as with quantum Monte Carlo impurity solvers Gull  _et al._ (2011).

In this paper, we apply the ‘differentiable programming’ methodology to NRG,
and demonstrate some advantages of its use for practical physics applications.
We show that the core eigensolver routine in NRG is well-suited to Automatic
Differentiation (AD), an emerging programming paradigm that has its origin in
Deep Learning (DL) Goodfellow  _et al._ (2016). In DL, hugely complicated
neural networks (NN) feed their output into a loss function, which must be
optimized with respect to the NN parameters (weights). This can be
accomplished efficiently via gradient descent using AD programming: the
gradient of the loss-function with respect to the NN weights can be computed
at about the cost of the evaluation of the NN Bartholomew-Biggs  _et al._
(2000). The remarkable aspect of this technique is that the computed gradient
is not approximated, but is numerically exact. Without AD, derivatives can of
course be approximated using finite difference (FD) derivatives. However,
without care, FD derivatives can be inaccurate; when treated rigorously to
guarantee convergence, FD can become computationally expensive. The efficiency
and accuracy of AD owes to the fact that if each step of a computer program
can be differentiated, then outputs can be differentiated with respect to
inputs via the chain rule. Evaluating derivatives is trivial and cheap since
the program itself is differentiated.

AD has been made more accessible through the introduction of AD libraries such
as jax Bradbury  _et al._ (2018), tensorflow Abadi  _et al._ (2015) or PyTorch
Paszke  _et al._ (2019). The advantage for computational physics is that when
formulated using an AD library, exact derivatives can be obtained ‘for free’
Tamayo-Mendoza  _et al._ (2018); Zhang  _et al._ (2019); Xie  _et al._ (2020);
McGreivy  _et al._ (2021); Ponsioen  _et al._ (2021); Dick and Fernandez-Serra
(2021); Coopmans  _et al._ (2021). Recent applications in physics include
optimal control in quantum systems Vargas-Hernández  _et al._ (2021),
mitigating the sign problem in Monte Carlo Wan  _et al._ (2020), and
optimizing tensor networks Liao  _et al._ (2019). In addition to its use for
optimization problems, AD has particular appeal in physics because physical
quantities are related to derivatives of generating functionals Tu  _et al._
(2021).

Here we formulate a differentiable NRG ($\partial$NRG) scheme. The key step is
an explicit expression for the derivative of NRG eigenvalues and eigenvectors:
the exact gradient of the entire eigensolver routine at any given NRG step is
determined symbolically. This eliminates the need for full backpropagation,
and when combined with AD for the other program elements, provides highly
efficient access to derivatives of program outputs with respect to model
parameters. For example, thermodynamic observables may be obtained from
derivatives of the free energy; and the physical response of a system to
perturbations can be simply and cheaply determined.

##  II Forward- and backward-mode AD

The AD approach is similar to symbolic differentiation, with standard
derivative rules applied algorithmically Linnainmaa (1976); Goodfellow  _et
al._ (2016). Regarding a program as a function $f$ that maps a given input $x$
to an output $y$, we may write $f(x)=y$. This is typically achieved in
practice by concatenating several primitive functions $f^{(i)}$, viz:

$f=f^{(n)}\circ f^{(n-1)}\circ...\circ f^{(1)}\;,$ (1)

where “$\circ$” denotes the composition of two functions. We denote the number
primitive functions comprising $f$ as $\operatorname{\textrm{ops}}(f)$. The
flow of data between primitive functions can be expressed as a compute graph
Tinhofer  _et al._ (2012); for simplicity, we consider here only programs
whose graph has a chain topology as in Eq. (1). To compute the derivative of
$f$ with respect to $x$ one can apply the chain rule,

$\displaystyle\frac{\partial f}{\partial x}$ $\displaystyle=\frac{\partial
f^{(n)}}{\partial f^{(n-1)}}\frac{\partial f^{(n-1)}}{\partial
f^{(n-2)}}...\frac{\partial f^{(1)}}{\partial x}\;.$ (2)

Eq. (2) can be evaluated algorithmically in either the forward mode, or the
backward mode. In the forward mode, we apply the recursion

$\frac{\partial f^{(i)}}{\partial x}=\frac{\partial f^{(i)}}{\partial
f^{(i-1)}}\delta_{i-1}\equiv\delta_{i},~{}\textrm{with}~{}\delta_{0}=1,~{}f^{(0)}=x$
(3)

starting with $i=1$ and increasing up to $i=n$, with $f^{(n)}\equiv f$ such
that $\delta_{n}=\partial f/\partial x$ is the desired derivative. In backward
mode (also known as backpropagation) we use the recursion

$\frac{\partial f}{\partial f^{(i)}}=\frac{\partial f^{(i+1)}}{\partial
f^{(i)}}\bar{\delta}_{i+1}\equiv\bar{\delta}_{i},~{}\textrm{with}~{}\bar{\delta}_{n}=1,~{}f^{(0)}=x$
(4)

starting from $i=n-1$ and decreasing down to $i=0$, with
$\bar{\delta}_{0}=\partial f/\partial x$ the desired derivative.

These forward and backward mode AD recursions are symbolically equivalent
methods for calculating derivatives. In the practical implementation of AD the
intermediate derivatives $\delta_{i}$ or $\bar{\delta}_{i}$ are evaluated
numerically and stored in memory. This provides a way to compute numerically
exact derivatives, with the true analytical result recovered when the value of
$\delta_{i}$ is not truncated Griewank and Walther (2008).

By contrast, the FD approximation reads,

$\frac{\partial f}{\partial x}\approx\frac{f(x+h)-f(x)}{h}\equiv
D_{h}[f](x)\;,$ (5)

with $h$ being the FD value. The exact limit $h\rightarrow 0^{\pm}$ cannot be
taken numerically, and so $\partial_{x}f$ must be approximated for finite $h$
and the convergence $D_{h}[f](x)\rightarrow\partial_{x}f$ checked explicitly.
This may require many function evaluations (runs of the whole program) since
for very small $h$, $f(x)$ and $f(x+h)$ may be numerically indistinguishable
(truncation error). Furthermore, convergence should be checked for both the
forwards ($h>0$) and backwards ($h<0$) difference quotient. The situation gets
worse for higher-order derivatives, and numerical convergence of FD
approximations can be costly and/or fraught Grossmann  _et al._ (2007).

Above we considered scalar functions of the form
$f:\mathbb{R}\rightarrow\mathbb{R}$. However, both AD and FD approaches can be
straightforwardly generalized to compute the Jacobian of vector-valued
functions defined over a vector space
$f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$. The FD and forward/backward mode
AD scale differently with respect to the input dimension $n$ and output
dimension $m$. The FD exhibits the most expensive scaling of computational
cost with $\sim 2nm$ function evaluations required to compute the Jacobian of
$f$ in the single-shot case. By contrast, AD methods do not require the
function to be evaluated in order to compute the derivative, as with the
symbolic approach (although since many derivatives are related to their anti-
derivative, function evaluations can be reused to reduce computational cost at
the expense of increased memory cost Baydin  _et al._ (2018)). Instead, the
computational cost of AD methods is controlled by
$\operatorname{\textrm{ops}}(f)$. Specifically, the forward mode scales as
$\sim c_{F}n\operatorname{\textrm{ops}}(f)$, while the backward mode scales as
$\sim c_{B}m\operatorname{\textrm{ops}}(f)$, with $c_{F},c_{B}<6$ (although
the forward mode typically outperforms the backward mode unless $m\ll n$
Baydin  _et al._ (2018)). The computational cost of the FD approach is
therefore typically much higher than either of the AD modes; in this work we
use forward mode AD unless stated otherwise. A comparison of the performance
and precision of the different methods, as applied to NRG, is presented in
Appendix A.

Through AD one can obtain the numerically exact derivative of any program with
a single run, as long as all operations comprising the program are themselves
differentiable. Thus, one can differentiate a given physics solver program
simply by utilizing known derivatives of its constituent operations.

On the other hand, in certain cases it may be possible to find the symbolic
derivative of the entire solver with respect to certain input parameters. In
the following we show how the latter can be achieved for NRG through analysis
of the backpropagation chain.

##  III NRG

Wilson’s NRG is a numerical solver for quantum impurity models, of the type
$\hat{H}=\hat{H}_{\rm imp}+\hat{H}_{\rm bath}+\hat{H}_{\rm hyb}$, where
$\hat{H}_{\rm imp}$ is the ‘impurity’ Hamiltonian describing a few local,
interacting quantum degrees of freedom, while $\hat{H}_{\rm
bath}=\sum_{\sigma,\mathbf{k}}\epsilon_{\mathbf{k}}^{\phantom{\dagger}}\hat{c}^{\dagger}_{\sigma\mathbf{k}}\hat{c}_{\sigma\mathbf{k}}^{\phantom{\dagger}}$
describes a non-interacting fermionic bath. The coupling between them is given
by $\hat{H}_{\rm
hyb}=\sum_{\sigma,\mathbf{k}}V_{\mathbf{k}}^{\phantom{\dagger}}(\hat{d}^{\dagger}_{\sigma}\hat{c}_{\sigma\mathbf{k}}^{\phantom{\dagger}}+\hat{c}^{\dagger}_{\sigma\mathbf{k}}\hat{d}_{\sigma}^{\phantom{\dagger}})$,
where we have assumed for simplicity here that a single impurity orbital
$\hat{d}_{\sigma}$ couples to the bath.

Due to the impurity interactions, such a model is a genuine quantum many body
problem, and in general has no exact solution Hewson (1997). NRG treats
instead a discretized version of the model, which becomes computationally
tractable through a process of iterative diagonalization and truncation. This
constitutes an RG procedure, in which useful physical information can be
extracted at each step, as progressively lower energy scales are probed. Full
details can be found in Refs. Wilson (1975); Bulla  _et al._ (2008); here we
introduce only the elements necessary to formulate $\partial$NRG.

The first step is the logarithmic discretization of the bath and mapping to a
Wilson chain. The continuous local density of states of the uncoupled bath at
the impurity position is divided up into intervals of decreasing width,
defined by the points $x_{n}^{\pm}=\pm D\Lambda^{-n}$, where $D$ is the half-
bandwidth, $\Lambda>1$ the discretization parameter, and $n=0,1,2,...$. The
spectrum is discretized by replacing the continuous density in each interval
by a single pole of the same weight and at the average position. A semi-
infinite tight-binding chain (the Wilson chain) is then defined such that the
spectral function at the end is the same as the discretized bath spectrum
(this is achieved in practice by the Lanczos algorithm). The discretized bath
Hamiltonian reads,

$\hat{H}_{\rm bath}~{}\mapsto~{}\hat{H}_{\rm bath}^{\rm
disc}=\sum_{\sigma}\sum^{\infty}_{n=0}t_{n}(\hat{c}^{\dagger}_{\sigma
n+1}\hat{c}_{\sigma n}+\hat{c}^{\dagger}_{\sigma n}\hat{c}_{\sigma n+1})\;,$
(6)

where we have assumed particle-hole symmetry here for simplicity. The original
model is recovered as $\Lambda\rightarrow 1$. The specific form of the Wilson
chain hopping parameters $t_{n}$ depends on details of the dispersion relation
$\epsilon_{\mathbf{k}}$. However, for a metallic system, $t_{n}\sim
D\Lambda^{-n/2}$ for large $n$.

For a Wilson chain of total length $N$, adding an extra site $N+1$ is
therefore always a small perturbation. The NRG exploits this energy scale
separation down the chain through the iterative diagonalization scheme.
Starting from the impurity, the chain is built up by successively adding new
Wilson chain sites. At each step, the system is diagonalized, and high-energy
states discarded. One defines a sequence of rescaled Hamiltonians, comprising
the impurity and the first $N$ sites of the Wilson chain,

$\displaystyle\hat{H}_{N}=\Lambda^{(N-1)/2}\left[\hat{H}_{\text{imp}}+V\sum_{\sigma}(\hat{d}^{\dagger}_{\sigma}\hat{c}_{\sigma
0}+\hat{c}^{\dagger}_{\sigma 0}\hat{d}_{\sigma})\right.$
$\displaystyle\left.+\sum_{\sigma}\sum^{N-1}_{n=0}t_{n}(\hat{c}^{\dagger}_{\sigma
n+1}\hat{c}_{\sigma n}+\hat{c}^{\dagger}_{\sigma n}\hat{c}_{\sigma
n+1})\right]\;,$ (7)

with $V$ the effective impurity-bath coupling. The full (discretized)
Hamiltonian is recovered in the limit
$\hat{H}=\lim_{N\rightarrow\infty}\Lambda^{-(N-1)/2}\hat{H}_{N}$. The sequence
of Hamiltonians $\hat{H}_{N}$ satisfy the recursion relation,

$\hat{H}_{N+1}=\sqrt{\Lambda}\hat{H}_{N}+\Lambda^{N/2}\sum_{\sigma}t_{N}(\hat{c}^{\dagger}_{\sigma
N+1}\hat{c}_{\sigma N}+\hat{c}^{\dagger}_{\sigma N}\hat{c}_{\sigma N+1})\;.$
(8)

At each step $N$, the Hamiltonian $\hat{H}_{N}$ is diagonalized to find the
eigenvectors (states) $\\{|i\rangle_{N}\\}$ and eigenvalues (energies)
$\\{E_{N;i}\\}$ that satisfy the Schrödinger equation,

$\hat{H}_{N}\ket{i}_{N}=E_{N;i}\ket{i}_{N}\;.$ (9)

We denote excitation energies relative to the ground state as $\Delta
E_{N;i}=E_{N;i}-E_{N;0}$, and $\hat{H}_{\rm imp}\equiv\hat{H}_{-1}$. To
construct $\hat{H}_{N+1}$ we add another Wilson chain site. The Fock space of
$\hat{H}_{N+1}$ is spanned by basis states,

$\ket{i;s}_{N+1}=\ket{i}_{N}\otimes\ket{s}\;,$ (10)

comprising the tensor product of eigenstates $|i\rangle_{N}$ of $\hat{H}_{N}$,
and the added chain site $N+1$ denoted $|s\rangle$. Matrix elements of
$\hat{H}_{N+1}$ in this basis read,

$\hat{H}_{N+1}(is;i^{\prime}s^{\prime})=\;_{N+1}\\!\bra{i;s}\hat{H}_{N+1}\ket{i^{\prime};s^{\prime}}_{N+1}\;.$
(11)

Finally, $\hat{H}_{N+1}$ is diagonalized to find the new eigenbasis at step
$N+1$, viz:

$\ket{j}_{N+1}=\sum_{is}U_{N+1}(j,is)\ket{i;s}_{N+1}\;.$ (12)

Eq. (11) can be simplified using the energies $E_{N;i}$ of $\hat{H}_{N}$ and
the tensor $\eta_{ss^{\prime}\sigma}=\bra{s^{\prime}}\hat{c}_{\sigma
N+1}\ket{s}$ (which is independent of $N$), viz

$\displaystyle\hat{H}_{N+1}(is;i^{\prime}s^{\prime})=$
$\displaystyle\sqrt{\Lambda}\delta_{ss^{\prime}}\delta_{ii^{\prime}}\times
E_{N;i}$ (13)
$\displaystyle+(-1)^{s}\Lambda^{N/2}t_{N}\sum_{\sigma}\eta_{ss^{\prime}\sigma}\eta_{N;ii^{\prime}\sigma}+{\rm
H.c.}\;,$

where
$\eta_{N;jj^{\prime}\sigma}=\sum_{\bar{i}\bar{s},\bar{i}^{\prime}\bar{s}^{\prime}}U^{\dagger}_{N}(j,\bar{i}\bar{s})~{}U_{N}(j^{\prime},\bar{i}^{\prime}\bar{s}^{\prime})\times\delta_{\bar{i}\bar{i}^{\prime}}\eta_{\bar{s}^{\prime}\bar{s}\sigma}$
and $s=\\{0,-1,+1,2\\}$ for
$|s\rangle=\\{|0\rangle,|\downarrow\rangle,|\uparrow\rangle,|\uparrow\downarrow\rangle\\}$.

When following these steps, the dimension of the Fock space grows
exponentially with the length of the chain. This is avoided in NRG by
truncating the Fock space at each step, discarding high-energy states. In
practice, one retains $M_{K}$ of the lowest energy eigenstates of
$\hat{H}_{N}$ at each step. The NRG approximation Wilson (1975); Weichselbaum
and von Delft (2007) is that the states and energies of $\hat{H}_{N}$
approximate those of the full $\hat{H}$. This approximation is justified by
the exponential decay of the Wilson chain coefficients $t_{n}$, which means
that high-energy states discarded at a given iteration do not become low-lying
states at a later iteration. Convergence of the NRG calculation can be checked
_post hoc_ by increasing $M_{K}$. Hereafter it is to be understood that the
condition $i\leq M_{K}$ applies to Eq. (10) (with the eigenvalues sorted
lowest to highest). The dimensionality of the Fock space is therefore constant
at each step, and the NRG calculation scales _linearly_ with chain length $N$.
With increasing $N$ the physics on successively lower energy scales is
uncovered.

The entire process of going from one iteration to the next can be summarized
as an RG transformation,

$\hat{H}_{N+1}=R[\hat{H}_{N}]\;.$ (14)

The full (discretized, renormalized) Hamiltonian $\hat{H}_{N}$ can therefore
be constructed iteratively starting from the impurity Hamiltonian
$\hat{H}_{\rm imp}$. Regarding the latter as a function of $n$ physical model
parameters $\\{\theta\\}$, the entire NRG algorithm can be viewed as a
function,

$\displaystyle f(\\{\theta\\})\equiv R^{\circ N+1}[\hat{H}_{\rm
imp}(\\{\theta\\})]$ $\displaystyle
f:\mathbb{R}^{n}\rightarrow\mathbb{C}^{M_{K}\times M_{K}}\;,$ (15)

where we have used the shorthand notation for an $N\\!+\\!1$-fold
concatenation of functions $R^{\circ N+1}=R\circ R\circ...\circ R$. All
physical quantities of interest for the original quantum impurity model are
obtained from the eigenvalues and eigenvectors of the set of
$\\{\hat{H}_{N}\\}$ Wilson (1975); Bulla  _et al._ (2008); Weichselbaum and
von Delft (2007).

Note that an operator $\hat{O}$, acting only on impurity degrees of freedom,
can be expressed in the in the eigenbasis of $\hat{H}_{N}$. We denote this
matrix as
$\hat{O}_{N}(i;i^{\prime})=\;_{N}\\!\bra{i}\hat{O}\ket{i^{\prime}}_{N}$. From
Eq. (12) it follows that,

$\displaystyle\hat{O}_{N+1}(j;j^{\prime})$
$\displaystyle=\sum_{is,i^{\prime}s^{\prime}}\\!\\!{\vphantom{\sum}}^{\prime}U^{\dagger}_{N+1}(j,is)~{}U_{N+1}(j^{\prime},i^{\prime}s^{\prime})\times\hat{O}_{N+1}(is;i^{\prime}s^{\prime})$
$\displaystyle=\sum_{is,i^{\prime}s^{\prime}}\\!{\vphantom{\sum}}^{\prime}U^{\dagger}_{N+1}(j,is)~{}U_{N+1}(j^{\prime},i^{\prime}s^{\prime})\times\hat{O}_{N}(i;i^{\prime})\delta_{s,s^{\prime}}$
(16)

where the primed sum implies that $i,i^{\prime}\leq M_{K}$, and the second
line follows from the fact that $\hat{O}$ does not act on states $|s\rangle$
of site $N+1$. Eq. (III) can be used iteratively to transform
$\hat{O}_{-1}(i;i^{\prime})$, evaluated explicitly in the impurity basis, into
the eigenbasis of any later iteration.

Symbolically, we denote this process,

$\hat{O}_{N+1}=P[\hat{O}_{N}]\;,$ (17)

which we refer to as propagating forward the operator.

##  IV Derivative of eigenvalues and eigenvectors in NRG

As outlined in the previous section, NRG outputs the eigenvalues and
eigenvectors for the set of $\hat{H}_{N}$. Since all physical quantities for
the quantum impurity model are obtained from these, calculating their
derivatives using Eq. (2) requires the derivatives of the eigenvalues and
eigenvectors of $\hat{H}_{N}$. From Eq. (III), we regard NRG as a function
defined over an $n$ dimensional domain, with an $M_{K}\times M_{K}$
dimensional image. Since $n$ denotes the number of model parameters in
$\hat{H}_{\rm imp}$, and $M_{K}\times M_{K}$ is the Fock space dimension of
$\hat{H}_{N}$, we have $M_{K}\times M_{K}\gg n$ and therefore forward mode AD
is significantly faster Griewank and Walther (2008). In the following we
consider only forward mode AD.

For a Hermitian matrix $\hat{A}$ (such as the Hamiltonian $\hat{H}_{N}$) with
non-degenerate eigenvectors $|i\rangle$ and distinct eigenvalues $\lambda_{i}$
satisfying $\hat{A}|i\rangle=\lambda_{i}|i\rangle$, we may express the
differentials through Seeger  _et al._ (2017); Giles (2008),

$\displaystyle\operatorname{\text{d}\\!}\lambda_{i}$
$\displaystyle=\bra{i}\operatorname{\text{d}\\!}\hat{A}\ket{i}\;,$ (18a)
$\displaystyle\operatorname{\text{d}\\!}\ket{i}$ $\displaystyle=\sum_{i\neq
j}\frac{\bra{j}\operatorname{\text{d}\\!}\hat{A}\ket{i}}{\lambda_{i}-\lambda_{j}}\ket{j}\;.$
(18b)

These expressions are of course well-known in the context of non-degenerate
perturbation theory Feynman (1939).

However, the assumption implicit in Eq. (18) that $\hat{A}$ must be free of
eigenvalue degeneracies is rather restrictive. Methods treating the general
degenerate case are much more computationally involved and may suffer from
numerical instabilities Dailey (1989); Giles (2008). Of course in physics,
where $\hat{A}$ is the Hamiltonian of some system, energy eigenvalue
degeneracies are common. This may be because of underlying non-Abelian
symmetries (for example SU(2) spin symmetry), which endow the spectrum of the
Hamiltonian with a multiplet structure. In this case, the Hamiltonian becomes
block diagonal in the associated conserved quantum numbers, and the
diagonalization Eq. (12) can be done separately in each block. Alternatively,
NRG can be formulated directly in multiplet space by using the Wigner-Eckart
theorem Weichselbaum (2012). Either approach removes the problem of such
degeneracies in Eq. (18). However, this may not fully solve the problem, since
there could be accidental degeneracies or emergent symmetries Wilson (1975);
Hewson (1997); Keller  _et al._ (2014) that lead to additional energy
eigenvalue degeneracies in physical systems.

To overcome this, we make a simple approximation. To the diagonal entries of
$\hat{H}_{N}$ we add noise, with random variables drawn from a normal
distribution of width $\sigma$. The noise lifts the eigenvalue degeneracy,
meaning that derivatives can be obtained using Eq. (18). Note that the
smallest terms in the rescaled Hamiltonian $\hat{H}_{N}$ are $\mathcal{O}(1)$
Wilson (1975); Bulla  _et al._ (2008), so $\sigma\ll 1$ constitutes a small
perturbation. Care should be taken to add the noise in such a way to respect
bare symmetries.

With a straight application of AD via Eq. (3), the derivatives of eigenvalues
and eigenvectors for $\hat{H}_{N}$ require the evaluation of Eq. (18) at every
NRG step (and hence noise must be added at every NRG step). For small enough
noise width $\sigma$, physical properties at a given iteration should be
unaffected. Indeed, in Appendix A we show for the AIM that NRG remains highly
accurate when $\sigma\lesssim 10^{-6}$. At larger noise levels, errors may
propagate through the iterative process (and will snowball if they introduce
RG relevant perturbations), so care must be taken to avoid this with such an
AD approach. On the other hand, forward mode AD for the AIM requires
$\sigma\gtrsim 10^{-6}$ to stabilize the calculation of derivatives using Eq.
(18). Smaller noise levels introduce more severe numerical instabilities
because of the recursive nature of the derivative calculation in AD via Eq.
(3) – derivatives at one iteration depend on those of previous iterations.
This stability-accuracy tradeoff in AD is analyzed in Appendix A.

Below, we derive an alternative formulation used in $\partial$NRG, which
utilizes the symbolic derivative of the entire eigensolver for $\hat{H}_{N}$.
This has the advantage that Eq. (18) need only be applied once, at the NRG
step for which derivatives are required. Not only is this far less
computationally costly, but it also avoids accumulating errors introduced by
noise, which can then be added at a much lower level.

###  IV.1 Formulation in $\partial$NRG

In the context of NRG, one can formulate a differential recursion relation
based on the RG transformation Eq. (14),

$\operatorname{\text{d}\\!}\hat{H}_{N}=\operatorname{\text{d}\\!}\hat{H}_{N-1}\frac{\operatorname{\text{d}\\!}R(\hat{H}_{N-1})}{\operatorname{\text{d}\\!}\hat{H}_{N-1}}\;.$
(19)

This allows derivatives of $\hat{H}_{N}$ to be computed in AD forward or
backward mode. However, one can also reformulate the problem in a much more
simple and efficient way. Taking derivatives with respect to model parameters
$\theta$ of $\hat{H}_{\rm imp}$, we may write from Eq. (III) the operator
identity
$\partial_{\theta}\hat{H}_{N}\simeq\Lambda^{(N-1)/2}~{}\partial_{\theta}\hat{H}_{\rm
imp}$ (where $\simeq$ follows from the NRG approximation). Since
$\partial_{\theta}\hat{H}_{\rm imp}$ only involves impurity operators, it can
be trivially evaluated in the impurity basis, $[\partial_{\theta}\hat{H}_{\rm
imp}]_{-1}(k;k^{\prime})={}_{-1}\langle k|\partial_{\theta}\hat{H}_{\rm
imp}|k^{\prime}\rangle_{-1}$. Eq. (17) can then be used to propagate this
operator forward into the basis of $\hat{H}_{N}$. It therefore follows that,

$[\partial_{\theta}\hat{H}_{N}]_{N}\simeq\Lambda^{(N-1)/2}~{}P^{\circ
N+1}([\partial_{\theta}\hat{H}_{\rm imp}]_{-1})\;.$ (20)

From Eq. (18) we then obtain,

$\displaystyle\partial_{\theta}E_{N;i}$
$\displaystyle\simeq\Lambda^{(N-1)/2}~{}P^{\circ
N+1}([\partial_{\theta}\hat{H}_{\rm imp}]_{-1})(i;i)$ (21a)
$\displaystyle\partial_{\theta}\ket{i}_{N}$
$\displaystyle\simeq\Lambda^{(N-1)/2}\sum_{i\neq j}\frac{P^{\circ
N+1}([\partial_{\theta}\hat{H}_{\rm
imp}]_{-1})(j;i)}{E_{N;i}-E_{N;j}}\ket{j}_{N}$ (21b)

Eq. (21) is the main result of this work. It shows how the derivative of
energies and states of the NRG Hamiltonians $\hat{H}_{N}$ can be obtained by
simply forward-propagating the impurity operators
$\partial_{\theta}\hat{H}_{\rm imp}$. Derivatives of physical quantities can
then be related to those obtained via Eq. (21). The above can also be
generalised to deal with derivatives $\partial_{\phi}$ of other Hamiltonian
parameters $\phi$ (for example $V$ or $t_{n}$ in Eq. (III)) and higher order
derivatives (see Appendix B).

To calculate first-order derivatives of eigenvalues and eigenvectors in AD
forward mode (assuming the simplest NRG implementation) one must perform seven
additional operations for every NRG step. The $\partial$NRG approach involves
only one additional operation per NRG step and does not require any
propagation of derivatives related to the chain rule (see Appendix C). Hence
$\partial$NRG, utilizing Eq. (21), is far more efficient than standard forward
mode AD. Performance benchmarking is demonstrated in A.

Derivatives with respect to multiple parameters can be obtained
straightforwardly by propagating forward any $[\partial_{\theta}\hat{H}_{\rm
imp}]_{-1}$ of interest using Eq. (III) in $\partial$NRG. Therefore,
$\partial$NRG scales linearly with the number of derivatives $n$, but is
independent of $\operatorname{\textrm{ops}}(f)$, unlike forward mode AD.

In practice, the use of Eq. (21) again requires lifting eigenvalue
degeneracies by adding a small diagonal noise term. We found in numerical
tests for the AIM that $\sigma\simeq 10^{-11}$ is sufficient for stabilizing
the numerics in $\partial$NRG, without noticeably impacting any measurable
physical properties – see Appendix A.

If derivatives at only step $N$ are required, such a noise term need only be
introduced at that step (rather than adding noise at each step) and the RG
flow is unaffected. This is typically the situation for ground-state
properties or low-temperature thermodynamics. On the other hand, if
information is required from every NRG step (for example in the calculation of
dynamical quantities Weichselbaum and von Delft (2007)) noise must be
introduced at each iteration. Unlike with the straight AD implementation,
$\partial$NRG allows us to mitigate the possibility of snowballing errors
introduced by propagating the noise terms – if this level of accuracy should
be required (e.g. in the vicinity of a quantum critical point). This is
because the eigensystem derivatives at different $N$ via Eq. (21) are
independent in $\partial$NRG. Noise may be added to $\hat{H}_{N}$ for the
purpose of evaluating Eq. (21), but the pristine $\hat{H}_{N}$ without noise
can be used for the main NRG recursion. This gives the most accurate and
reliable results, at the cost of an additional matrix diagonalization. See
Appendices A and C for performance comparisons (and note that this more
rigorous approach is still faster than FD and straight AD in many
circumstances, and certainly more accurate).

##  V Application to Anderson Impurity Model

We illustrate the use of Eq. (21) by applying the $\partial$NRG scheme to the
paradigmatic Anderson impurity model (AIM), for which

$\hat{H}_{\text{imp}}=\epsilon\hat{n}+U\hat{n}_{\uparrow}\hat{n}_{\downarrow}\;,$
(22)

where $\hat{n}=\sum_{\sigma}\hat{n}_{\sigma}$ and
$\hat{n}_{\sigma}=\hat{d}^{\dagger}_{\sigma}\hat{d}^{\phantom{\dagger}}_{\sigma}$.
With $\hat{H}_{\rm imp}$ so defined, the exact (discretized) $\hat{H}_{N}$ is
given by Eq. (III). In NRG, $\hat{H}_{N}$ is approximated through the RG
procedure Eq. (14), $\hat{H}_{N}=R^{\circ N+1}(\hat{H}_{\text{imp}})$.

For a given $\hat{H}_{N}$, we obtain the partition function
$\mathcal{Z}_{N}(\bar{\beta})=\sum_{i}e^{-\bar{\beta}E_{N;i}}$ and free energy
$F_{N}(\bar{\beta})=-\frac{1}{\beta}\log[\mathcal{Z}_{N}(\bar{\beta})]$. In
the original Wilsonian formulation Wilson (1975); Bulla  _et al._ (2008), the
effective inverse temperature $\beta\equiv 1/k_{\rm B}T$ is related to the NRG
iteration number (Wilson chain length) $N$ via
$\beta=\Lambda^{(N-1)/2}\bar{\beta}$, with $\bar{\beta}=\mathcal{O}(1)$ in
practice. With this definition, the NRG free energy at inverse temperature
$\bar{\beta}$ is a good approximation to the true free energy of the original
(undiscretized) model at inverse temperature $\beta$ Wilson (1975); Bulla  _et
al._ (2008), $F(\beta)\simeq F_{N}(\bar{\beta})$. The corresponding
differential follows as,

$\operatorname{\text{d}\\!}F_{N}=\frac{\Lambda^{-(N-1)/2}}{\mathcal{Z}_{N}}\sum_{i}e^{-\bar{\beta}E_{N;i}}\operatorname{\text{d}\\!}E_{N;i}\;.$
(23)

For the AIM, we may use Eq. (21) to obtain derivatives of the free energy with
respect to the impurity parameters $\epsilon$ and $U$. For example,

$\partial_{\epsilon}F_{N}=\frac{1}{\mathcal{Z}_{N}}\sum_{i}e^{-\bar{\beta}E_{N;i}}P^{\circ
N+1}([\hat{n}]_{-1})(i;i)\;.$ (24)

Since
$\partial_{\epsilon}F_{N}(\bar{\beta})=\langle\hat{n}\rangle_{H_{N},\bar{\beta}}\simeq\langle\hat{n}\rangle_{H,\beta}$,
Eq. (24) is precisely equivalent to the standard Wilsonian approach to
calculating local thermodynamic expectation values in NRG Bulla  _et al._
(2008). The above illustration demonstrates that $\partial$NRG is analytically
equivalent to the well-known result for such local thermodynamic quantities.

Another commonly computed quantity for such models is the local impurity
magnetic susceptibility at zero field,

$\displaystyle\chi(T)=\int^{\beta}_{0}\operatorname{\text{d}\\!}\tau\langle\hat{S}_{z}(\tau)\hat{S}_{z}\rangle-\beta\langle\hat{S}_{z}\rangle^{2}$

where $\hat{S}_{z}=\tfrac{1}{2}(\hat{n}_{\uparrow}-\hat{n}_{\downarrow})$ is
the impurity spin operator and $\tau$ is imaginary time. This can be
alternatively obtained by adding a source term $B\hat{S}z$ to the Hamiltonian
Eq. 22, and then taking the second-order derivative of the free energy with
respect to $B$, evaluated at $B=0$. That is, we can write
$\chi(T)=\frac{\partial^{2}}{\partial B^{2}}F\big{|}_{B=0}$. Since
$\partial$NRG is able to deal with second (and higher) order derivatives (see
Appendix B) we calculate
$\chi_{H_{N},\bar{\beta}}=\frac{\partial^{2}}{\partial
B^{2}}F_{N}(\bar{\beta})\big{|}_{B=0}$. This is obtained automatically in
$\partial$NRG, but from Eq. (29) it can also be expressed as,

$\displaystyle\chi_{H_{N},\bar{\beta}}=-\frac{\bar{\beta}}{\mathcal{Z}^{2}_{N}}\left(\sum_{i}e^{-\bar{\beta}E_{i;N}}P^{\circ
N+1}([\hat{S}_{z}]_{-1})(i;i)\right)^{2}$
$\displaystyle+\frac{\bar{\beta}}{\mathcal{Z}_{N}}\sum_{i}e^{-\bar{\beta}E_{i;N}}|P^{\circ
N+1}([\hat{S}_{z}]_{-1})(i;i)|^{2}$
$\displaystyle-\frac{2}{\mathcal{Z}_{N}}\sum_{i}e^{-\bar{\beta}E_{i;N}}\sum_{i\neq
j}\frac{|P^{\circ N+1}([\hat{S}_{z}]_{-1})(j;i)|^{2}}{E_{N;i}-E_{N;j}}\;.$
(25)

This form of $\chi_{H_{N},\bar{\beta}}$ is equivalent to the Lehmann
representation of $\chi(T)$ evaluated in NRG Hamiltonian $\hat{H}_{N}$ at
effective inverse temperature $\bar{\beta}=\beta\Lambda^{-(N-1)/2}$ Bulla  _et
al._ (2008). Since all degeneracies are lifted in $\partial$NRG, convergence
factors typically used in the Lehmann representation of dynamical quantities
are not needed here. Note that $\chi_{H_{N},\bar{\beta}}$ can be obtained at
exactly zero magnetic field, so no symmetries are broken and no limit is taken
numerically. Since the derivatives of primitive program functions are
symbolic, source terms added to $\hat{H}_{\rm imp}$ may be evaluated at zero
coupling constant and still yield finite derivatives. However, we note that
such symmetry-breaking source terms cannot be added to the Hamiltonian if non-
Abelian quantum numbers are implemented from the outset. In the case of the
magnetic susceptibility, we can therefore utilize conserved total $S_{z}$, but
not conserved total $S$, when adding the source term $B\hat{S}_{z}$ (even
though in the end we set $B=0$). Numerical results are presented in Appendix
A, and reproduce precisely the results of standard NRG using Eq. (V).

Similarly, $\partial$NRG may be used to obtain the charge susceptibility,
which in the wide-band limit is given simply by
$\chi_{C}(T)\simeq\partial^{2}F_{N}(\bar{\beta})/\partial\epsilon^{2}=\partial_{\epsilon}\langle\hat{n}\rangle_{H_{N},\bar{\beta}}$
(and can be obtained by replacing $\hat{S_{z}}$ by $\hat{n}$ in Eq. (V)).
Indeed, Maxwell relations provide non-trivial connections between physical
quantities obtained as derivatives. This was exploited recently in Ref. Han
_et al._ (2021) to extract the fractional entropy of multi-channel Kondo
states in quantum dot experiments, comparing with NRG calculations of charge
derivatives.

The above examples for the AIM demonstrate that the $\partial$NRG algorithm is
equivalent to known expressions for certain derived quantities. However, the
power of $\partial$NRG is that it does this automatically within a generalized
framework, and works equally well for any derivative in any quantum impurity
model.

Figure 1: Ground state energy $E_{N;0}$ of the AIM (green line) and its
normalized derivative $\partial_{U}E_{N;0}/\nu_{0}$ obtained by $\partial$NRG
(blue line) and by FD (red circles) for iteration $N=5$ (top panel) and $N=40$
(bottom panel). Results plotted as a function of $U$, with $\epsilon=-U/2$ and
constant $J\equiv 8V^{2}/U=0.3$. We define
$\nu_{0}=\partial_{U}E_{N;0}|_{U\rightarrow 0}$ as the derivative in the limit
$U\to 0$. Figure 2: Free energy derivatives of the AIM obtained by
$\partial$NRG. Inset shows
$\partial_{\epsilon}F_{N}\simeq\langle\hat{n}\rangle_{T}$; main panel shows
$\partial_{U}(\partial_{\epsilon}F_{N})\simeq\partial_{U}\langle\hat{n}\rangle_{T}$,
comparing AD (blue lines) with FD (red points). Same parameters used as for
Fig. 1, but with $N=20$ and $\bar{\beta}=0.9$ yielding an effective
temperature $T/D\simeq 3\times 10^{-5}$.

##  VI Numerical Results

For the following numerical demonstrations, we implemented a basic NRG code in
jax Bradbury  _et al._ (2018), using the included AD routines to obtain
derivatives via Eq. (21). The code is available open-source at Ref. rig
(2021). For simplicity we did not exploit quantum numbers here, and so
eigenstate degeneracies were removed by adding Gaussian noise with a small
variance $\sigma\simeq 10^{-11}$ (FD derivatives were calculated for
comparison without noise, $\sigma=0$). In the following we set the conduction
electron half-bandwidth to $D=1$ and use an NRG discretization parameter
$\Lambda=3$. Further details on the numerical calculations and the finite
difference derivatives can be found in Appendix D and Appendix E respectively.

First we use $\partial$NRG to compute the derivative of the ground-state
energy $E_{N;0}$ of the AIM with respect to the interaction $U$. AD results in
Fig. 1 (blue lines) are compared with FD approximations (points), as a
function of $U$ at iteration $N=5$ (upper panel) and $N=40$ (lower panel),
normalized by their respective $U\to 0$ derivatives. The green lines show the
variation of $E_{N;0}$ itself. The results show the non-trivial effect of
renormalization going from $N=5$ to $N=40$ at intermediate $U$, as well as the
saturation of the ground state derivatives at both large and small $U$. In
this case we see excellent agreement between AD and FD results (although the
former are far less computationally expensive to obtain).

Fig. 2 demonstrates the use of $\partial$NRG to obtain thermodynamic
quantities from derivatives of the NRG free energy. The inset shows the
impurity occupation $\langle\hat{n}\rangle\simeq\partial_{\epsilon}F_{N}$ at a
temperature $T/D\sim 10^{-5}$ (corresponding to $N=20$ and $\bar{\beta}=0.9$)
for the same systems as in Fig. 1. Since $\epsilon=-U/2$, the AIM possesses an
exact particle-hole symmetry and hence is at half-filling,
$\langle\hat{n}\rangle=1$. The AD results (blue line) satisfy this exact
result precisely, while the FD results (points) show some numerical error.
More interestingly, the $\partial$NRG framework allows to obtain higher
derivatives with equal ease (Eq. (29) is used instead of Eq. (18a)), as shown
in the main panel of Fig. 2. Here we calculate the corresponding second
derivatives $\partial^{2}F_{N}/\partial
U\partial\epsilon\simeq\partial_{U}\langle\hat{n}\rangle$, which again show
non-trivial behaviour as a function of interaction strength $U$. The FD
approximations agree well, but are much more costly to obtain Grossmann  _et
al._ (2007), requiring for every point five executions of the entire NRG code
per second derivative, and an expensive convergence test. $\partial$NRG
requires only a single function evaluation (see Appendix C).

Figure 3: Top panel: NRG energy level flow diagram for the AIM. The lowest 32
rescaled excitation energies $\Delta E_{N;i}$ above the ground state at
iteration $N$ are plotted as a function of $N$ (for even $N$ only). Plotted
for model parameters $U=0.01$, $\epsilon=-U/2$, and $J\equiv 8V^{2}/U=0.2$,
with the free orbital (FO), local moment (LM) and strong coupling (SC) regimes
marked. Bottom panel: the corresponding derivatives with respect to $U$
obtained by $\partial$NRG.

As a final application of $\partial$NRG for the AIM, we turn to the RG energy
level flow diagram shown in the top panel of Fig. 3. The excitation spectrum
of the effective rescaled $\hat{H}_{N}$ plotted against iteration number $N$
shows the well-known flow between fixed points Wilson (1975); Bulla  _et al._
(2008); Hewson (1997). In the example shown, the crossover scales are well-
separated, such that we see distinct level structures associated with the free
orbital (FO), local moment (LM), and strong coupling (SC) fixed points as
marked on the diagram. At a fixed point, the level structure does not change
with $N$. Indeed, the RG structure of the problem and resulting universality
implies that the fixed point level structure is the same, independent of the
microscopic model parameters $U$ ($\equiv-2\epsilon$) and $V$ – only the
crossover scales between fixed points are affected.

This is demonstrated in the bottom panel of Fig. 3, where we plot the
derivatives with respect to $U$ of the NRG energy levels $E_{N;i}$, again as a
function of iteration number $N$. As expected, the derivatives vanish at the
fixed points (the level structure does not depend on $U$ at the fixed points);
but there is a strong dependence along the crossovers, since the crossover
scales depend on $U$.

##  VII Outlook: dynamics
and the vertex

The above numerical results for the AIM are provided as a demonstration proof-
of-principle. Future useful applications exploiting the full power of
$\partial$NRG may be found for more complex models, situations involving
higher-order derivatives, optimization techniques requiring exact gradients,
or in cases where partial derivatives of NRG outputs with respect to inputs
are difficult to obtain by standard FD means. An example of the latter is the
derivative of frequency-dependent dynamical quantities, such as impurity
Green’s functions or the conduction electron scattering T-matrix, with respect
to bare model parameters.

In the context of dynamical mean field theory (DMFT) for correlated materials
Georges  _et al._ (1996); Kotliar  _et al._ (2006); Sheridan  _et al._ (2021),
model machine learning techniques Rigo and Mitchell (2020) could be applied to
optimize simpler effective models with respect to more complicated microscopic
ones, by comparing their Green’s functions. Given the non-trivial dependence
of such dynamical quantities at different energies on model parameters, the
exact gradient within the AD framework becomes an essential ingredient for
gradient descent optimization.

Refs. Kugler  _et al._ (2021); Lee  _et al._ (2021) recently uncovered the
analytic structure of the full local vertex and proposed a scheme to compute
it within NRG. The vertex is an important object, entering for example in
extensions of DMFT beyond the local limit Rohringer  _et al._ (2018). It is
possible that some partial information on the vertex could also be obtained by
$\partial$NRG. This is inspired by the result in Ref. van Loon  _et al._
(2020) for the AIM which, within the Matsubara formalism, connects the vertex
at zero bosonic frequency $F^{\rm loc}_{\nu\nu^{\prime}}(\omega=0)$ to the
functional derivative of the interaction self-energy $\Sigma_{\nu}$ with
respect to the hybridization $\Delta_{\nu^{\prime}}$, viz:

$\frac{\delta\Sigma_{\nu}}{\delta\Delta_{\nu^{\prime}}}=T[G_{\nu^{\prime}}]^{2}F^{\textrm{loc}}_{\nu\nu^{\prime}}(\omega=0)\;,$
(26)

where $G_{\nu^{\prime}}$ is the single-particle impurity Matsubara Green’s
function. For a precise definition and discussion of $F^{\rm
loc}_{\nu\nu^{\prime}}(\omega=0)$, see Ref. van Loon  _et al._ (2020); Kugler
_et al._ (2021). Since $G_{\nu^{\prime}}$ and $\Sigma_{\nu}$ can be calculated
within standard NRG Osolin  _et al._ (2013), we argue that such an object is
obtainable within $\partial$NRG. We speculate that a Keldysh version of Eq. 26
may similarly provide access to certain information about the real-frequency
vertex at finite temperatures.

In order to calculate such derivatives of dynamical quantities, further code
development is required, since the full-density-matrix NRG method would need
to be implemented using AD libraries such as jax Bradbury  _et al._ (2018) and
integrated within the $\partial$NRG scheme described in this paper. We leave
this for future work.

##  VIII Conclusion

NRG is the gold-standard method of choice for solving generalized quantum
impurity models Wilson (1975); Bulla  _et al._ (2008). In this work, we
introduce a new variant of the standard algorithm – $\partial$NRG – which
makes use of the differential programming paradigm to automatically and
efficiently obtain derivatives of NRG outputs with respect to input model
parameters.

We make use of the AD jax library Bradbury  _et al._ (2018), together with a
bespoke routine based on Eq. (21) which allows the derivatives of Hamiltonian
eigenvalues and eigenvectors to be obtained at any iteration in an accurate
and inexpensive way. Our fully commented code is available open-source rig
(2021).

We demonstrated the use of $\partial$NRG by application to the Anderson
impurity model, for which we obtained the derivative of NRG energy levels with
respect to model parameters, calculated thermodynamic quantities from
derivatives of the NRG free energy, and susceptibilities from derivatives of
Hamiltonian source terms. $\partial$NRG may be useful for machine learning
applications involving NRG Rigo and Mitchell (2020) for which gradient descent
optimization requires derivatives of a loss function; for adaptive broadening
schemes Lee and Weichselbaum (2016); or for optimal control protocols Vargas-
Hernández  _et al._ (2021). Richer physical information may be obtained from
derivatives of dynamical quantities. This also opens the door to automatic
differentiable DMFT, with $\partial$NRG as the impurity solver. Finally, we
note that the $\partial$NRG methodology is compatible with interleaved-NRG
(iNRG) Mitchell  _et al._ (2014); Stadler  _et al._ (2016) and generalizations
utilising full symmetries Weichselbaum (2012). This is left for future work.

###### Acknowledgements.

_Acknowledgments.–_ The authors thank Lei Wang, Luuk Coopmans, Jan von Delft,
Seung-Sup Lee, and Emma Minarelli for the helpful discussions. We acknowledge
funding from the Irish Research Council Laureate Awards 2017/2018 through
grant IRCLA/2017/169.

Figure 4: (a) Dependence of $T_{K}$ on the noise variance $\sigma$ at
different $M_{K}$. The number of kept states ranges from $150$ to $500$ in
steps of $50$. The AIM model parameters used are
$\epsilon=-0.15,~{}U=0.3,~{}V=0.1~{}B=0$. (b) Zero-field impurity
magnetization $\langle\hat{S}_{z}\rangle$ computed via derivatives of the free
energy using $\partial$NRG and forward mode AD, as a function of Wilson chain
length $N$ (for $M_{K}=500$). The noise variance is $\sigma=1\times
10^{-8},3\times 10^{-8},1\times 10^{-7},3\times 10^{-7}...,1\times 10^{-2}$.
Same model parameters as in (a).

##  Appendix A Numerical benchmarking: accuracy and speedup

In this Appendix we demonstrate that: (i) NRG and $\partial$NRG are not
affected by the addition of Gaussian noise with a small variance; and that
(ii) the real-world performance of $\partial$NRG exceeds that of basic AD and
FD in terms of both accuracy and speed.

A degeneracy-free NRG Hamiltonian is required for application of Eq. (18).
However, physical systems often do have eigenvalue degeneracies, arising for
example from underlying symmetries. Although one can partially mitigate this
problem by exploiting Abelian and non-Abelian quantum numbers, some
degeneracies may remain. Here we consider the ‘worst case’ scenario in which
no symmetries are exploited.

First we examine the effect of adding Gaussian noise with mean $\mu=0$ and
variance $\sigma$ to the diagonal elements of $\hat{H}_{N}$. We use the Kondo
temperature $T_{K}$ as one figure-of-merit for assessing the effect of the
noise term. For simplicity we define the Kondo temperature through the
impurity entropy via $S_{\rm imp}(T=T_{K})=0.5$ (suitably between the local
moment and strong coupling fixed point values). We compute $T_{K}$ for
$10^{-8}\leq\sigma\leq 10^{-2}$ and different numbers of kept states $100\leq
M_{K}\leq 500$ (using fixed $\Lambda=3$). Fig. 4(a) shows clearly that for
$\sigma\leq 10^{-6}$ the Kondo temperature does not noticeably depend on the
noise level (we have confirmed explicitly down to $\sigma=10^{-15}$ that the
results are fully converged for each $M_{K}$). Other physical quantities
computed in standard NRG show similar behaviour.

However, the effect of adding noise is more pronounced in derivatives. In Fig.
4(b) we compare
$\partial_{B}F_{N}(\bar{\beta})|_{B=0}=\langle\hat{S}_{z}\rangle_{H_{N},\bar{\beta}}$
as a function of Wilson chain length $N$, as computed with $\partial$NRG and
straight AD for different $\sigma$. Since the impurity magnetization is
evaluated at zero field, $SU(2)$ spin symmetry implies the exact result
$\langle\hat{S}_{z}\rangle_{H_{N},\bar{\beta}}=0$. However, this spin symmetry
is not enforced, and so we see deviations due to the noise. As Fig. 4(b)
shows, $\partial$NRG accurately approximates the exact result even at
relatively high $\sigma$. By contrast, AD yields derivatives that strongly
depend on $\sigma$ and have much higher error than $\partial$NRG at a given
$\sigma$. Indeed, derivatives are numerically not computable for all $N$ in
standard forward mode AD for $\sigma\lesssim 10^{-6}$; at larger noise levels,
AD derivatives are computable but the accuracy can become poor, especially at
later iterations due to the propagation and accumulation of errors through the
derivative chain rule.

Therefore, although NRG for the AIM is insensitive to noise for $\sigma\leq
10^{-6}$, the AD approach is not stable at these noise levels. At higher noise
levels, AD is stabilized, but the accuracy suffers. On the other hand,
$\partial$NRG avoids such problems: a lower noise level can be used since each
derivative calculation is independent, and highly accurate results can be
obtained.

Figure 5: (a) Magnetic susceptibility $\chi$ calculated with $\partial$NRG for
different noise variances $\sigma$. (b) $\chi$ computed using forward mode AD
and $10^{-8}\leq\sigma\leq 10^{-2}$. (c) $\chi$ computed with $\partial$NRG
using the pole broadening approach, with $10^{-15}\leq\delta^{2}\leq 10^{0}$.
All calculations are performed with $M_{K}=500$ and with the same AIM model
parameters as in Fig. 4.

Within $\partial$NRG, the magnetization calculation via the analog of Eq. (24)
shown in Fig. 4(b) does not involve eigenvector derivatives, and is therefore
particularly stable. By contrast, the AD calculation of any derivative is
built up recursively via Eq. (3) and therefore does involve the computation of
eigenvector derivatives at each and every step. This contributes to the
relative performance gain in $\partial$NRG.

In Fig. 5 we examine the magnetic susceptibility $\chi(T)$. The calculation of
this quantity, which is obtained automatically in $\partial$NRG via the second
derivative of the NRG free energy, is formally equivalent to the analytic
expression, Eq. (V) (and does involve eigenvector derivatives). Gaussian noise
of width $\sigma$ is added to the diagonal of $\hat{H}_{N}$ at iteration $N$.
The calculated $\chi_{H_{N},\bar{\beta}}$ at $\bar{\beta}=0.9$ (used
hereafter) is then interpreted as the true $\chi(T)$ at inverse temperature
$\beta=\Lambda^{(N-1)/2}\bar{\beta}$. The numerical results from $\partial$NRG
show that $\chi_{H_{N},\bar{\beta}}$ is obtained reliably at all $N$ (and
hence all $T$) for $\sigma$ as low as $10^{-11}$. Only at smaller $\sigma$
does the method break down: derivatives are then not computable for earlier
iterations/higher temperatures (the correct low-$T$ behavior is however still
captured). As demonstrated in Fig. 4 and confirmed in Fig. 5(a), highly
accurate results are obtained for $\sigma<10^{-6}$. Therefore in
$\partial$NRG, we have a wide window $10^{-11}<\sigma<10^{-6}$ over which
numerical results are fully converged and stable.

By contrast, in Fig. 5(b) we show results for the same quantity obtained by
straight AD. The calculation is not numerically stable for small $\sigma$, but
very inaccurate for large $\sigma$ (typical of the breakdown for higher-order
derivatives). As such there is no reliable regime for which robust results can
be obtained by straight AD.

In Fig. 5(c) we examine the feasibility of an alternative approach to the
eigenvalue-degeneracy problem that avoids adding Gaussian noise. The method,
proposed by Liao et al in Ref. Liao  _et al._ (2019), consists of
reformulating Eq. (18b),

$\operatorname{\text{d}\\!}\ket{i}=\sum_{i\neq
j}\frac{\bra{j}\operatorname{\text{d}\\!}\hat{A}\ket{i}}{\lambda_{i}-\lambda_{j}}\ket{j}\equiv\sum_{i\neq
j}\bra{j}\operatorname{\text{d}\\!}\hat{A}\ket{i}f(\lambda_{i}-\lambda_{j})\ket{j}\;,$
(27)

with $f(x)=1/x$. Divergences induced by eigenvalue degeneracies can be avoided
by replacing the function $f(x)$ with the approximate broadened form
$f(x)\approx x/(x^{2}+\delta^{2})$ with $\delta\ll 1$. This broadening
approximation is known to distort somewhat the overall results, but has the
advantage that divergences are removed without the need to lift degeneracies.
Fig. 5(c) shows numerical results for $\chi(T)$ obtained in this way, as a
function of $T$ for different broadenings $\delta$. Although the method always
yields a computable result, the true result (dashed line) is only obtained in
the $\delta\to 0$ limit at the lowest temperature scales. Finite $\delta$ at
intermediate $T$ appears to yield rather poor results. We conclude that
$\partial$NRG is indeed the method of choice in this context, in terms of
accuracy.

Figure 6: Representative calculation runtime in seconds, vs number of kept
states $M_{K}$, comparing different automatic derivative techniques and
computed quantities. In all cases shown, blue lines correspond to
$\partial$NRG with single diagonalization, orange lines for $\partial$NRG with
double diagonalization, red lines for forward mode AD, and green lines for
backward mode AD. (a) Calculation at $N=40$ of
$\langle\hat{S}_{z}\rangle_{H_{N},\bar{\beta}}$ (solid lines) compared with
the calculation of both $\langle\hat{S}_{z}\rangle_{H_{N},\bar{\beta}}$ and
$\langle\hat{n}\rangle_{H_{N},\bar{\beta}}$ (dashed lines), obtained via the
first derivatives of the NRG free energy. (b) Calculation for all $N\leq 40$
of $\langle\hat{S}_{z}\rangle_{H_{N},\bar{\beta}}$ requiring first derivatives
(solid lines), and $\chi_{H_{N},\bar{\beta}}$ requiring second derivatives
(dashed lines). Model parameters as in Fig. 4.

We now turn to an analysis of the real-world performance of $\partial$NRG in
terms of calculation time, comparing with straight AD implemented in jax
Bradbury  _et al._ (2018). The figure of merit is the runtime measured in
seconds $[s]$ of both algorithms run on the same machine (in this case an AMD
Threadripper 2950X platform running Python 3.8.10 with the `OMP_NUM_THREAD =
16` flag; further package specifications can be found in sum (2019)). We do
not utilize any quantum numbers here.

In Fig. 6(a) we compare the runtime for the calculation of (i)
$\langle\hat{S}_{z}\rangle_{H_{N},\bar{\beta}}$ (solid lines), and (ii)
$\langle\hat{S}_{z}\rangle_{H_{N},\bar{\beta}}$ and
$\langle\hat{n}\rangle_{H_{N},\bar{\beta}}$ together (dashed lines); using
$\partial$NRG (blue), forward mode AD (red), and backward-mode AD (green), at
the last NRG iteration, $N=40$. We see that for all $M_{K}$, $\partial$NRG
performs best (with the relative speedup becoming more pronounced at larger
$M_{K}$). For straight AD, we see that forward mode beats backward mode for
the calculation of a single derivative; but since backward mode AD scales with
the number of inputs rather than outputs, it will overtake the forward mode
when many derivatives are calculated.

In Fig. 6(b) we compare runtimes using $\partial$NRG with a single
diagonalization of $\hat{H}_{N}$ (blue), $\partial$NRG with double
diagonalization of $\hat{H}_{N}$ (orange), and forward mode AD (red). We
calculate $\langle\hat{S}_{z}\rangle_{H_{N},\bar{\beta}}$ (solid lines), and
$\chi_{H_{N},\bar{\beta}}$ (dashed lines) at _all_ iterations $N\leq 40$. For
the orange lines, _two_ diagonalizations of $\hat{H}_{N}$ are performed at
each step (with and without noise), which thereby eliminates error propagation
from the added noise and provides the most accurate calculation (this may or
may not be needed in practice, depending on the situation). We see that
$\partial$NRG with a single diagonalization per step is the fastest in all
cases. For the simpler quantity
$\langle\hat{S}_{z}\rangle_{H_{N},\bar{\beta}}$ (which involves only a first-
derivative), the $\partial$NRG using two diagonalizations and forward mode AD
have similar runtimes, although AD is slightly faster. However, for
$\chi_{H_{N},\bar{\beta}}$, both versions of $\partial$NRG are considerably
faster. This is because the calculation of $\chi_{H_{N},\bar{\beta}}$ requires
a second-derivative, which is much more computationally expensive in AD, but
almost as cheap in $\partial$NRG. The relative performance gain for
$\partial$NRG also increases if several derivatives are computed.

In conclusion, in a typical setting, $\partial$NRG is considerably more
efficient than AD (often by orders of magnitude) – while at the same time
being more accurate.

##  Appendix B Higher order derivatives

Eq. (18) allows us to take the derivatives of eigenvalues and eigenvectors of
some Hermitian, non-degenerate Hamiltonian $\hat{H}$. Consider now a
Hamiltonian that can be linearly decomposed as,

$\hat{H}(\\{\theta\\})=\sum_{i}\theta_{i}\hat{h}_{i}\;,$ (28)

where $\hat{h}_{i}$ is an operator defined over the same Hilbert space as
$\hat{H}$, and with the corresponding coupling constant $\theta_{i}$ a real
scalar. In this case, all second order derivatives of the Hamiltonian must
vanish, $\frac{\partial^{2}\hat{H}}{\partial\theta_{i}\partial\theta_{j}}=0$.
It is then possible to derive a closed-form expression for the second-order
derivatives of eigenvalues and eigenvectors with respect to the couplings
$\theta_{\mathbf{1}},\theta_{\mathbf{2}}\in\\{\theta\\}$, viz:

$\displaystyle\partial_{\theta_{\mathbf{1}}}\partial_{\theta_{\mathbf{2}}}E_{i}$
$\displaystyle=2\sum_{i\neq j}\frac{\langle
i|\partial_{\theta_{\mathbf{2}}}\hat{H}|j\rangle\langle
j|\partial_{\theta_{\mathbf{1}}}\hat{H}|i\rangle}{E_{i}-E_{j}}\;,$ (29a)
$\displaystyle\partial_{\theta_{\mathbf{1}}}\partial_{\theta_{\mathbf{2}}}|i\rangle$
$\displaystyle=\sum_{i\neq j}\left[\Pi_{ij}|j\rangle+\frac{\langle
j|\partial_{\theta_{\mathbf{1}}}\hat{H}|i\rangle}{\lambda_{i}-\lambda_{j}}\sum_{j\neq
k}\frac{\langle
k|\partial_{\theta_{\mathbf{2}}}\hat{H}|j\rangle}{\lambda_{j}-\lambda_{k}}|k\rangle\right]\;,$
(29b)

where,

$\displaystyle\Pi_{ij}=-\frac{\langle
i|\partial_{\theta_{\mathbf{2}}}\hat{H}|i\rangle-\langle
j|\partial_{\theta_{\mathbf{2}}}\hat{H}|j\rangle}{(\lambda_{i}-\lambda_{j})^{2}}\langle
j|\partial_{\theta_{\mathbf{1}}}\hat{H}|i\rangle+\sum_{j\neq k}\frac{\langle
j|\partial_{\theta_{\mathbf{2}}}\hat{H}|k\rangle}{\lambda_{j}-\lambda_{k}}\langle
k|\partial_{\theta_{\mathbf{1}}}\hat{H}|i\rangle+\sum_{k\neq i}\frac{\langle
k|\partial_{\theta_{\mathbf{2}}}\hat{H}|i\rangle}{\lambda_{i}-\lambda_{k}}\langle
j|\partial_{\theta_{\mathbf{1}}}\hat{H}|k\rangle\;.$ (30)

With these formulae one can compute physical observables depending on second-
derivatives, such as susceptibilities, using $\partial$NRG.

##  Appendix C Comparing differentiation methods

Forward Primal Trace |  Forward Tangent Trace
---|---
$v_{-3}=\Lambda$ |  $\operatorname{\text{d}\\!}v_{-3}=0$
$v_{-2}=\\{t_{n}\\}$ |  $\operatorname{\text{d}\\!}v_{-2}=0$
$v_{-1}=\eta$ |  $\operatorname{\text{d}\\!}v_{-1}=0$
$v_{0}=H_{\rm imp}$ |  $\operatorname{\text{d}\\!}v_{0}=\operatorname{\text{d}\\!}H_{\rm imp}$
$v_{1}=\textrm{EIGENVALUES}(v_{0})$ |  $[\operatorname{\text{d}\\!}v_{1}]_{i}=[v_{2}]^{\dagger}_{i}\cdot\operatorname{\text{d}\\!}v_{0}\cdot[v_{2}]_{i}$
$v_{2}=\textrm{EIGENVECTORS}(v_{0})$ |  $[\operatorname{\text{d}\\!}v_{2}]_{i}=\sum_{j\neq i}\frac{[v_{2}]^{\dagger}_{j}\cdot\operatorname{\text{d}\\!}v_{0}\cdot[v_{2}]_{i}}{[v_{1}]_{i}-[v_{1}]_{j}}[v_{2}]_{j}$
$v_{3}=v_{-1}^{T}$ |  $\operatorname{\text{d}\\!}v_{3}=0$
FOR $n$ = 0 to $N$ |  FOR $n$ = 0 to $N$
$\qquad v_{7n+4}=v_{7n}$ |  $\qquad\operatorname{\text{d}\\!}v_{7n+4}=\operatorname{\text{d}\\!}v_{7n}$
$\qquad[v_{7n+5}]_{iiss}=v_{-3}^{1/2}[v_{7n+1}]_{l}$ |  $\qquad[\operatorname{\text{d}\\!}v_{7n+5}]_{iiss}=v_{-3}^{1/2}[\operatorname{\text{d}\\!}v_{7n+2}]_{l}$
$\qquad[v_{7n+6}]_{ii^{\prime}ss^{\prime}}=(-1)^{i^{\prime}}[v_{-2}]_{n}\sum_{\sigma}[v_{7n+3}]_{\sigma ii^{\prime}}[v_{-1}]_{\sigma ss^{\prime}}$ |  $\qquad[\operatorname{\text{d}\\!}v_{7n+6}]_{ii^{\prime}ss^{\prime}}=(-1)^{i^{\prime}}[v_{-2}]_{n}\sum_{\sigma}[\operatorname{\text{d}\\!}v_{7n+3}]_{\sigma ii^{\prime}}[v_{-1}]_{\sigma ss^{\prime}}$
$\qquad[v_{7n+7}]_{ii^{\prime}ss^{\prime}}=[v_{7n+5}]_{ii^{\prime}ss^{\prime}}+[v_{7n+6}]_{ii^{\prime}ss^{\prime}}+[v_{7n+6}]^{\dagger}_{ii^{\prime}ss^{\prime}}$ |  $\qquad[\operatorname{\text{d}\\!}v_{7n+8}]_{i}=[v_{7n+9}]^{\dagger}_{i}\cdot\operatorname{\text{d}\\!}v_{7n}\cdot[v_{7n+9}]_{i}$
$\qquad v_{7n+8}=\textrm{EIGENVALUES}(v_{7n+7})$ |  $\qquad[\operatorname{\text{d}\\!}v_{7n+8}]_{i}=[v_{7n+9}]^{\dagger}_{i}\cdot\operatorname{\text{d}\\!}v_{7n}\cdot[v_{7n+9}]_{i}$
$\qquad v_{7n+9}=\textrm{EIGENVECTORS}(v_{7n+7})$ |  $\qquad[\operatorname{\text{d}\\!}v_{7n+9}]_{i}=\sum_{j\neq i}\frac{[v_{7n+9}]^{\dagger}_{j}\cdot\operatorname{\text{d}\\!}v_{7n}\cdot[v_{7n+9}]_{i}}{[v_{7n+8}]_{i}-[v_{7n+8}]_{j}}[v_{7n+9}]_{j}$
$\qquad[v_{7n+10}]_{\sigma ii^{\prime}}=\sum_{ss^{\prime}}\sum_{jj^{\prime}}[v_{7n+9}]^{\dagger}_{i;sj}[v_{7n+9}]_{i;s^{\prime}j^{\prime}}[v_{-1}]_{\sigma ss^{\prime}}$ |  $\qquad[\operatorname{\text{d}\\!}v_{7n+10}]_{\sigma ii^{\prime}}=\sum_{ss^{\prime}}\sum_{jj^{\prime}}\left([\operatorname{\text{d}\\!}v_{7n+9}]^{\dagger}_{i;sj}[v_{7n+9}]_{i;s^{\prime}j^{\prime}}\right.$
$\qquad~{}~{}~{}~{}~{}~{}~{}+\left.[v_{7n+9}]^{\dagger}_{i;sj}[\operatorname{\text{d}\\!}v_{7n+9}]_{i;s^{\prime}j^{\prime}}\right)[v_{-1}]_{\sigma
ss^{\prime}}$
$y_{0}=v_{7N+7}$ |  $\operatorname{\text{d}\\!}y_{0}=\operatorname{\text{d}\\!}v_{7N+7}$
Table 1: Wengert list for NRG with forward mode AD. Left column: the primal trace prescribes all steps to be performed to execute an NRG calculation with $N$ iterations. Right column: the tangent trace prescribes all steps to compute the derivative the NRG Hamiltonian $\hat{H}_{N}$. Forward Primal Trace |  Forward Tangent Trace
---|---
$v_{-4}=\operatorname{\text{d}\\!}H_{\rm imp}$ |  $dv_{-4}=\operatorname{\text{d}\\!}H_{\rm imp}$
$v_{-3}=\Lambda$ |  $\operatorname{\text{d}\\!}v_{-3}=0$
$v_{-2}=\\{t_{n}\\}$ |  $\operatorname{\text{d}\\!}v_{-2}=0$
$v_{-1}=\eta$ |  $\operatorname{\text{d}\\!}v_{-1}=0$
$v_{0}=H_{\rm imp}$ |  $dv_{-4}=0$
$v_{1}=\textrm{EIGENVALUES}(v_{0})$ |
$v_{2}=\textrm{EIGENVECTORS}(v_{0})$ |
$v_{3}=v_{-1}^{T}$ |
$[v_{4}]_{ii^{\prime}}=\sum_{ss^{\prime}}\sum_{jj^{\prime}}[v_{2}]^{\dagger}_{i;sj}[v_{2}]_{i;s^{\prime}j^{\prime}}[v_{-4}]_{jj^{\prime}}\delta_{ss^{\prime}}$ |
FOR $n$ = 0 to $N$ |
$\qquad v_{8n+5}=v_{8n}$ |
$\qquad[v_{8n+6}]_{iiss}=v_{-3}^{1/2}[v_{8n+1}]_{l}$ |
$\qquad[v_{8n+7}]_{ii^{\prime}ss^{\prime}}=(-1)^{i^{\prime}}[v_{-2}]_{n}\sum_{\sigma}[v_{8n+3}]_{\sigma ii^{\prime}}[v_{-1}]_{\sigma ss^{\prime}}$ |
$\qquad[v_{8n+8}]_{ii^{\prime}ss^{\prime}}=[v_{8n+6}]_{ii^{\prime}ss^{\prime}}+[v_{8n+7}]_{ii^{\prime}ss^{\prime}}+[v_{8n+7}]^{\dagger}_{ii^{\prime}ss^{\prime}}$ |
$\qquad v_{8n+9}=\textrm{EIGENVALUES}(v_{8n+8})$ |
$\qquad v_{8n+10}=\textrm{EIGENVECTORS}(v_{8n+8})$ |
$\qquad[v_{8n+11}]_{\sigma ii^{\prime}}=\sum_{ss^{\prime}}\sum_{jj^{\prime}}[v_{8n+10}]^{\dagger}_{i;sj}[v_{8n+10}]_{i;s^{\prime}j^{\prime}}[v_{-1}]_{\sigma ss^{\prime}}$ |
$\qquad[v_{8n+12}]_{ii^{\prime}}=\sum_{ss^{\prime}}\sum_{jj^{\prime}}[v_{8n+10}]^{\dagger}_{i;sj}[v_{8n+10}]_{i;s^{\prime}j^{\prime}}[v_{8n+4}]_{jj^{\prime}}\delta_{ss^{\prime}}$ |
$y_{0}=v_{8N+8}$ |  $\operatorname{\text{d}\\!}y_{0}=v_{8N+12}$
Table 2: Wengert list for $\partial$NRG. Left column: the forward primal trace
for $\partial$NRG outputs the NRG Hamiltonain $\hat{H}_{N}$ and its derivative
$d\hat{H}_{N}$. Since derivatives are computed via Eq. (21) in the main
algorithm, there is no corresponding forward tangent trace in the Right
column.

The basic NRG algorithm is described in the main text, with the key steps
contained in Eq. (6)-(13). For more details, see Ref. Bulla  _et al._ (2008).
Assuming that degeneracies in the NRG Hamiltonian $\hat{H}_{N}$ are lifted, we
can use Eq. (18) to compute the derivatives of eigenvectors and eigenvalues of
$\hat{H}_{N}$. The Wengert list Wengert (1964); Baydin  _et al._ (2018) with
the forward primal trace (function evaluation) and the forward tangent trace
(function differentiation) can then be compiled, as shown in Table 1.
Similarly we compile the forward primal trace and forward tangent trace for
the $\partial$NRG algorithm, see Table 2. For $\partial$NRG, we have one
additional step in the primal trace for the main algorithm (corresponding to
Eq. (III)), but a trivial tangent trace.

We can now compare the efficiency of the two approaches. For forward mode AD
(fAD) we have $\operatorname{\textrm{ops}}(\textrm{fAD})=7N$, and to compute
forward primal and tangent traces,
$2\times\operatorname{\textrm{ops}}(\textrm{fAD})$ operations are required.
For $\partial$NRG by contrast,
$\operatorname{\textrm{ops}}(\textrm{$\partial$NRG})=8N$, but no other step is
required to compute the derivative of $\hat{H}_{N}$. However, this does not
mean that $\partial$NRG is twice as fast as fAD because the bottleneck
eigensolver routines Demmel  _et al._ (2008) appear only in the forward primal
trace. Nonetheless $\partial$NRG still provides a considerable performance
advantage, as established by explicit benchmarks in Appendix A

Figure 7: FD derivative errors $\Delta(h)$ vs FD value $h$ for the calculation
of $\partial_{U}\langle\hat{n}\rangle_{H_{N},\bar{\beta}}$ at $N=40$. Shown
for AIM model parameters $\epsilon=-U/2$, $V=0.1$ and $B=0$ with (a) $U=0.01$,
and (b) $U=100$. In all cases $M_{K}=500$.

##  Appendix D Details of numerics - NRG calculations

We used standard thermodynamic NRG Wilson (1975) to calculate the free energy,
impurity occupation and eigenvalues of $\hat{H}_{N}$. An NRG discretization
parameter $\Lambda=3$ was used, and $M_{K}=200$ states were retained at each
step for Figs. (1) and (3), while $M_{K}=600$ states were retained for Fig.
(2). No quantum numbers were used in this demonstration calculation.

##  Appendix E Details of finite difference calculations

We calculate the finite difference (FD) derivative of a function $f$ via,

$D_{h}[f](x)=\frac{f(x+h)-f(x)}{h}\;,$ (31)

where $h$ is the FD value. The FD derivative is connected to the definition of
the derivative by taking the limit, $\frac{\partial f}{\partial x}=\lim_{h\to
0}D_{h}[f](x)$.

In practice we compute FD derivatives for a set of finite $\\{h_{i}\\}$, and
check for convergence. However, care must be taken due to the trade-off
between round-off and truncation errors Baydin  _et al._ (2018). The
truncation error arises due to the finite $h>0$, which is required for the
numerical evaluation of $D_{h}[f](x)$, and which would diminish in the formal
limit $h\rightarrow 0$. However for very small $h$, the difference in the
function evaluations for $f(x)$ and $f(x+h)$ cannot be distinguished
numerically due to inevitable round-off errors in any finite precision
arithmetic. This leads to increasing errors as $h\rightarrow 0$.

We quantify the FD error as,

$\Delta(h)=\left|D_{h}[f](x)-\frac{\partial f}{\partial x}\right|\;.$ (32)

This quantity is plotted in Fig. 7 for the derivative of the AIM impurity
occupation with respect to the impurity interaction,
$\partial_{U}\langle\hat{n}\rangle_{H_{N},\bar{\beta}}$. We have used the
$\partial$NRG derivative as the true value of $\frac{\partial f}{\partial x}$
in Eq. (32). Fig. 7 shows rather typical behavior, with truncation errors
dominating at large $h$ and round-off errors dominating at small $h$, with a
stability plateau in between where the derivative should be evaluated.
However, comparison of the upper and lower panels (which correspond to model
parameters $U=10^{-2}$ and $U=10^{2}$ respectively), demonstrates an important
limitation of FD differentiation: the optimal $h$ is not fixed but depends on
input model parameters. Reliable results therefore require such a convergence
analysis for each new set of model parameters and for each new derivative.
This becomes computationally costly. The situation becomes significantly worse
when considering higher-order derivatives.

## References

  * Wilson (1975) Kenneth G Wilson, “The renormalization group: Critical phenomena and the kondo problem,” Reviews of modern physics 47, 773 (1975).
  * Bulla  _et al._ (2008) Ralf Bulla, Theo A Costi, and Thomas Pruschke, “Numerical renormalization group method for quantum impurity systems,” Reviews of Modern Physics 80, 395 (2008).
  * Hewson (1997) Alexander Cyril Hewson,  _The Kondo problem to heavy fermions_ , 2 (Cambridge university press, 1997).
  * Mitchell  _et al._ (2011) Andrew K Mitchell, Michael Becker, and Ralf Bulla, “Real-space renormalization group flow in quantum impurity systems: local moment formation and the kondo screening cloud,” Physical Review B 84, 115120 (2011).
  * Costi  _et al._ (2009) TA Costi, Lars Bergqvist, A Weichselbaum, J Von Delft, T Micklitz, A Rosch, P Mavropoulos, Peter H Dederichs, Francois Mallet, Laurent Saminadayar, _et al._ , “Kondo decoherence: Finding the right spin model for iron impurities in gold and silver,” Physical review letters 102, 056802 (2009).
  * Goldhaber-Gordon  _et al._ (1998) David Goldhaber-Gordon, Hadas Shtrikman, D Mahalu, David Abusch-Magder, U Meirav, and MA Kastner, “Kondo effect in a single-electron transistor,” Nature 391, 156–159 (1998).
  * Keller  _et al._ (2014) AJ Keller, S Amasha, I Weymann, CP Moca, IG Rau, JA Katine, Hadas Shtrikman, G Zaránd, and D Goldhaber-Gordon, “Emergent su (4) kondo physics in a spin–charge-entangled double quantum dot,” Nature Physics 10, 145–150 (2014).
  * Iftikhar  _et al._ (2015) Z Iftikhar, Sébastien Jezouin, A Anthore, U Gennser, FD Parmentier, A Cavanna, and F Pierre, “Two-channel kondo effect and renormalization flow with macroscopic quantum charge states,” Nature 526, 233–236 (2015).
  * Mitchell  _et al._ (2016) Andrew K Mitchell, LA Landau, L Fritz, and E Sela, “Universality and scaling in a charge two-channel kondo device,” Physical review letters 116, 157202 (2016).
  * Iftikhar  _et al._ (2018) Z Iftikhar, A Anthore, AK Mitchell, FD Parmentier, U Gennser, A Ouerghi, A Cavanna, C Mora, P Simon, and F Pierre, “Tunable quantum criticality and super-ballistic transport in a “charge” kondo circuit,” Science 360, 1315–1320 (2018).
  * Pouse  _et al._ (2021) Winston Pouse, Lucas Peeters, Connie L Hsueh, Ulf Gennser, Antonella Cavanna, Marc A Kastner, Andrew K Mitchell, and David Goldhaber-Gordon, “Exotic quantum critical point in a two-site charge kondo circuit,” arXiv preprint arXiv:2108.12691 (2021).
  * Liang  _et al._ (2002) Wenjie Liang, Matthew P Shores, Marc Bockrath, Jeffrey R Long, and Hongkun Park, “Kondo resonance in a single-molecule transistor,” Nature 417, 725–729 (2002).
  * Mitchell  _et al._ (2017) Andrew K Mitchell, Kim GL Pedersen, Per Hedegård, and Jens Paaske, “Kondo blockade due to quantum interference in single-molecule junctions,” Nature communications 8, 1–10 (2017).
  * Kotliar  _et al._ (2006) Gabriel Kotliar, Sergej Y Savrasov, Kristjan Haule, Viktor S Oudovenko, O Parcollet, and CA Marianetti, “Electronic structure calculations with dynamical mean-field theory,” Reviews of Modern Physics 78, 865 (2006).
  * Georges  _et al._ (1996) Antoine Georges, Gabriel Kotliar, Werner Krauth, and Marcelo J Rozenberg, “Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions,” Reviews of Modern Physics 68, 13 (1996).
  * Stadler  _et al._ (2015) KM Stadler, ZP Yin, J Von Delft, G Kotliar, and A Weichselbaum, “Dynamical mean-field theory plus numerical renormalization-group study of spin-orbital separation in a three-band hund metal,” Physical review letters 115, 136401 (2015).
  * Weichselbaum  _et al._ (2009) Andreas Weichselbaum, Frank Verstraete, Ulrich Schollwöck, J Ignacio Cirac, and Jan von Delft, “Variational matrix-product-state approach to quantum impurity models,” Physical Review B 80, 165117 (2009).
  * Gull  _et al._ (2011) Emanuel Gull, Andrew J Millis, Alexander I Lichtenstein, Alexey N Rubtsov, Matthias Troyer, and Philipp Werner, “Continuous-time monte carlo methods for quantum impurity models,” Reviews of Modern Physics 83, 349 (2011).
  * Goodfellow  _et al._ (2016) Ian Goodfellow, Yoshua Bengio, and Aaron Courville,  _Deep learning_ (MIT press, 2016).
  * Bartholomew-Biggs  _et al._ (2000) Michael Bartholomew-Biggs, Steven Brown, Bruce Christianson, and Laurence Dixon, “Automatic differentiation of algorithms,” Journal of Computational and Applied Mathematics 124, 171–190 (2000).
  * Bradbury  _et al._ (2018) James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, and Qiao Zhang, “JAX: composable transformations of Python+NumPy programs,”  (2018).
  * Abadi  _et al._ (2015) Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dandelion Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng, “TensorFlow: Large-scale machine learning on heterogeneous systems,”  (2015), software available from tensorflow.org.
  * Paszke  _et al._ (2019) Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala, “Pytorch: An imperative style, high-performance deep learning library,” in  _Advances in Neural Information Processing Systems 32_, edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett (Curran Associates, Inc., 2019) pp. 8024–8035.
  * Tamayo-Mendoza  _et al._ (2018) Teresa Tamayo-Mendoza, Christoph Kreisbeck, Roland Lindh, and Alán Aspuru-Guzik, “Automatic differentiation in quantum chemistry with applications to fully variational hartree–fock,” ACS central science 4, 559–566 (2018).
  * Zhang  _et al._ (2019) Shi-Xin Zhang, Zhou-Quan Wan, and Hong Yao, “Automatic differentiable monte carlo: Theory and application,” arXiv preprint arXiv:1911.09117 (2019).
  * Xie  _et al._ (2020) Hao Xie, Jin-Guo Liu, and Lei Wang, “Automatic differentiation of dominant eigensolver and its applications in quantum physics,” Physical Review B 101, 245139 (2020).
  * McGreivy  _et al._ (2021) Nick McGreivy, Stuart R Hudson, and Caoxiang Zhu, “Optimized finite-build stellarator coils using automatic differentiation,” Nuclear Fusion 61, 026020 (2021).
  * Ponsioen  _et al._ (2021) Boris Ponsioen, Fakher F Assaad, and Philippe Corboz, “Automatic differentiation applied to excitations with projected entangled pair states,” arXiv preprint arXiv:2107.03399 (2021).
  * Dick and Fernandez-Serra (2021) Sebastian Dick and Marivi Fernandez-Serra, “Using differentiable programming to obtain an energy and density-optimized exchange-correlation functional,” arXiv preprint arXiv:2106.04481 (2021).
  * Coopmans  _et al._ (2021) Luuk Coopmans, Di Luo, Graham Kells, Bryan K Clark, and Juan Carrasquilla, “Protocol discovery for the quantum control of majoranas by differentiable programming and natural evolution strategies,” PRX Quantum 2, 020332 (2021).
  * Vargas-Hernández  _et al._ (2021) Rodrigo A Vargas-Hernández, Ricky TQ Chen, Kenneth A Jung, and Paul Brumer, “Fully differentiable optimization protocols for non-equilibrium steady states,” arXiv preprint arXiv:2103.12604 (2021).
  * Wan  _et al._ (2020) Zhou-Quan Wan, Shi-Xin Zhang, and Hong Yao, “Mitigating sign problem by automatic differentiation,” arXiv preprint arXiv:2010.01141 (2020).
  * Liao  _et al._ (2019) Hai-Jun Liao, Jin-Guo Liu, Lei Wang, and Tao Xiang, “Differentiable programming tensor networks,” Physical Review X 9, 031041 (2019).
  * Tu  _et al._ (2021) Wei-Lin Tu, Huan-Kuang Wu, Norbert Schuch, Naoki Kawashima, and Ji-Yao Chen, “Generating function for tensor network diagrammatic summation,” Physical Review B 103, 205155 (2021).
  * Linnainmaa (1976) Seppo Linnainmaa, “Taylor expansion of the accumulated rounding error,” BIT Numerical Mathematics 16, 146–160 (1976).
  * Tinhofer  _et al._ (2012) Gottfried Tinhofer, Rudolf Albrecht, Ernst Mayr, Hartmut Noltemeier, and Maciej M Syslo,  _Computational graph theory_ , Vol. 7 (Springer Science & Business Media, 2012).
  * Griewank and Walther (2008) Andreas Griewank and Andrea Walther,  _Evaluating derivatives: principles and techniques of algorithmic differentiation_ (SIAM, 2008).
  * Grossmann  _et al._ (2007) Christian Grossmann, Hans-Görg Roos, and Martin Stynes,  _Numerical treatment of partial differential equations_ , Vol. 154 (Springer, 2007).
  * Baydin  _et al._ (2018) Atilim Gunes Baydin, Barak A Pearlmutter, Alexey Andreyevich Radul, and Jeffrey Mark Siskind, “Automatic differentiation in machine learning: a survey,” Journal of machine learning research 18 (2018).
  * Weichselbaum and von Delft (2007) Andreas Weichselbaum and Jan von Delft, “Sum-rule conserving spectral functions from the numerical renormalization group,” Physical review letters 99, 076402 (2007).
  * Seeger  _et al._ (2017) Matthias Seeger, Asmus Hetzel, Zhenwen Dai, Eric Meissner, and Neil D Lawrence, “Auto-differentiating linear algebra,” arXiv preprint arXiv:1710.08717 (2017).
  * Giles (2008) Mike Giles, “An extended collection of matrix derivative results for forward and reverse mode automatic differentiation,” (2008).
  * Feynman (1939) Richard Phillips Feynman, “Forces in molecules,” Physical review 56, 340 (1939).
  * Dailey (1989) R Lane Dailey, “Eigenvector derivatives with repeated eigenvalues,” AIAA journal 27, 486–491 (1989).
  * Weichselbaum (2012) Andreas Weichselbaum, “Non-abelian symmetries in tensor networks: A quantum symmetry space approach,” Annals of Physics 327, 2972–3047 (2012).
  * Han  _et al._ (2021) Cheolhee Han, Andrew K Mitchell, Zubair Iftikhar, Yaakov Kleeorin, Anne Anthore, Frédéric Pierre, Yigal Meir, and Eran Sela, “Extracting entropy of exotic quasiparticles from conductance measurements,” arXiv preprint arXiv:2108.12878 (2021).
  * rig (2021) https://github.com/JonasRigo/AD-Numerical-Renormalization-Group (2021).
  * Sheridan  _et al._ (2021) Evan Sheridan, Christopher Rhodes, Francois Jamet, Ivan Rungger, and Cedric Weber, “Data-driven dynamical mean-field theory: An error-correction approach to solve the quantum many-body problem using machine learning,” Physical Review B 104, 205120 (2021).
  * Rigo and Mitchell (2020) Jonas B Rigo and Andrew K Mitchell, “Machine learning effective models for quantum systems,” Physical Review B 101, 241105 (2020).
  * Kugler  _et al._ (2021) Fabian B Kugler, Seung-Sup B Lee, and Jan von Delft, “Multipoint correlation functions: spectral representation and numerical evaluation,” Physical Review X 11, 041006 (2021).
  * Lee  _et al._ (2021) Seung-Sup B Lee, Fabian B Kugler, and Jan von Delft, “Computing local multipoint correlators using the numerical renormalization group,” Physical Review X 11, 041007 (2021).
  * Rohringer  _et al._ (2018) G Rohringer, H Hafermann, A Toschi, AA Katanin, AE Antipov, MI Katsnelson, AI Lichtenstein, AN Rubtsov, and K Held, “Diagrammatic routes to nonlocal correlations beyond dynamical mean field theory,” Reviews of Modern Physics 90, 025003 (2018).
  * van Loon  _et al._ (2020) Erik GCP van Loon, Friedrich Krien, and Andrey A Katanin, “Bethe-salpeter equation at the critical end point of the mott transition,” Physical Review Letters 125, 136402 (2020).
  * Osolin  _et al._ (2013) Žiga Osolin _et al._ , “Padé approximant approach for obtaining finite-temperature spectral functions of quantum impurity models using the numerical renormalization group technique,” Physical Review B 87, 245135 (2013).
  * Lee and Weichselbaum (2016) Seung-Sup B Lee and Andreas Weichselbaum, “Adaptive broadening to improve spectral resolution in the numerical renormalization group,” Physical Review B 94, 235127 (2016).
  * Mitchell  _et al._ (2014) Andrew K Mitchell, Martin R Galpin, Samuel Wilson-Fletcher, David E Logan, and Ralf Bulla, “Generalized wilson chain for solving multichannel quantum impurity problems,” Physical Review B 89, 121105 (2014).
  * Stadler  _et al._ (2016) KM Stadler, AK Mitchell, J von Delft, and A Weichselbaum, “Interleaved numerical renormalization group as an efficient multiband impurity solver,” Physical Review B 93, 235101 (2016).
  * sum (2019) https://github.com/JonasRigo/Summer-Wine/wiki/Build-BLAS-on-a-TR2-NUMA-platform (2019).
  * Wengert (1964) Robert Edwin Wengert, “A simple automatic derivative evaluation program,” Communications of the ACM 7, 463–464 (1964).
  * Demmel  _et al._ (2008) James W Demmel, Osni A Marques, Beresford N Parlett, and Christof Vömel, “Performance and accuracy of lapack’s symmetric tridiagonal eigensolvers,” SIAM Journal on Scientific Computing 30, 1508–1526 (2008).