|  Rapidly convergent coupled-cluster Monte Carlo using a Chebyshev projector
---|---
|  Zijun Zhao,∗a‡ Maria-Andreea Filip,a§ and Alex J W Thoma
|  The multi-reference coupled-cluster Monte Carlo (MR-CCMC) algorithm is a
determinant-based quantum Monte Carlo (QMC) algorithm that is conceptually
similar to Full Configuration Interaction QMC (FCIQMC). It has been shown to
offer a balanced treatment of both static and dynamic correlation while
retaining polynomial scaling, although scaling to large systems with
significant strong correlation remained impractical. In this paper, we
document recent algorithmic advances that enable rapid convergence and a more
black-box approach to the multi-reference problem. These include a
logarithmically scaling metric-tree based excitation acceptance algorithm to
search for determinants connected to the model space at the desired excitation
level and a symmetry-screening procedure for the reference space. We show
that, for moderately sized model spaces, the new search algorithm brings about
an approximately 8-fold acceleration of one MR-CCMC iteration, while the
symmetry screening procedure reduces the number of active model space
determinants at essentially no loss of accuracy. We also introduce a
stochastic implementation of an approximate wall projector, which is the
infinite imaginary time limit of the exponential projector, using a truncated
expansion of the wall function in Chebyshev polynomials. We show that this
wall-Chebyshev projector requires significantly fewer applications of the
Hamiltonian to achieve the same statistical convergence. We have applied these
methods to calculate the binding curve of the beryllium and carbon dimers with
large active spaces, obtaining binding curves that are close to FCIQMC
quality. Notably, the wall-Chebyshev projector presented here can be applied
to accelerate any projector-based QMC algorithm. This paper represents a
significant step in the algorithmic maturation of the MR-CCMC algorithm and
paves the way for large-scale calculations.

##

††footnotetext: a Yusuf Hamied Department of Chemistry, University of
Cambridge, Cambridge, UK. E-mail<EMAIL_ADDRESS>‡ Present
address: Department of Chemistry, Emory University, Atlanta, GA, USA.
††footnotetext: § Present address: Max Planck Institute for Solid State
Research, Stuttgart, Germany.

##  1 Introduction

Quantum Monte Carlo (QMC) methods have long provided a powerful alternative to
conventional electronic structure methods, by generating high accuracy results
at a fraction of the cost of standard approaches. The combination of
Variational Monte Carlo (VMC)1, 2 and Diffusion Monte Carlo (DMC)3, 4, 5 has
become a significant benchmarking approach in many areas of electronic
structure,6, 7, 8, 9, 10, but it is limited by the need to provide some
approximate nodal structure to avoid collapse onto bosonic solutions.
Fermionic Monte Carlo methods11, 12 have since been developed which act
directly in the anti-symmetrised Hilbert space of the electronic structure
problem, thereby removing the potential for bosonic solutions a priori.

First introduced in 2009 by Booth et al.11, full configuration interaction
quantum Monte Carlo (FCIQMC) can be variously described as a stochastic power
iteration algorithm or an iterative solution to the imaginary time
Schrödinger’s equation.

Here, we give a brief summary of theoretical underpinnings of FCIQMC by taking
the latter view. By applying a Wick rotation 13, $\tau\leftarrow it,\
\tau\in\mathbb{R}$, to the time-dependent Schrödinger’s equation
${i\dot{\Psi}=\hat{H}\Psi}$, one obtains the imaginary time Schrödinger’s
equation:

$\diffp{}{\tau}|\Psi(\tau)\rangle=-\hat{H}|\Psi(\tau)\rangle,\
\tau\in\mathbb{R}.$ (1)

It can be formally integrated to give

$|\Psi(\tau)\rangle=e^{-\tau(\hat{H}-S)}|\Psi(0)\rangle,$ (2)

with $S$ being the constant of integration, also known as the ‘shift’.

The reference wavefunction $|\Psi(0)\rangle$, commonly a Hartree–Fock (HF)
solution, can be expanded in the eigenbasis of the full Hamiltonian,
$\\{|\Psi_{i}^{\mathrm{FCI}}\rangle\\}$, leading to

$|\Psi(\tau)\rangle=\sum_{i}e^{-\tau(E_{i}-S)}c_{i}|\Psi_{i}^{\mathrm{FCI}}\rangle,$
(3)

with $\\{E_{i}\\}$ being the eigenspectrum of the full Hamiltonian and
$|\Psi(0)\rangle=c_{i}|\Psi_{i}^{\mathrm{FCI}}\rangle$. We can see that, if
$S=E_{0}$, in the limit of $\tau\rightarrow\infty$, we obtain the ground state
of the full Hamiltonian.

By discretising the projector in Equation 2 and further applying the first-
order Taylor series expansion, we obtain the ‘master equation’ of FCIQMC:

$|\Psi(\tau+\delta\tau)\rangle=[1-\delta\tau(\hat{H}-S)]|\Psi(\tau)\rangle.$
(4)

This equation can be projected onto the different determinants in the Hilbert
space to give

$\langle D_{\bm{\mathrm{i}}}|\Psi(\tau+\delta\tau)\rangle=\langle
D_{\bm{\mathrm{i}}}|\Psi(\tau)\rangle-\delta\tau\langle
D_{\bm{\mathrm{i}}}|\hat{H}-S|\Psi(\tau)\rangle,$ (5)

which gives an update equation for the corresponding FCI parameters
$c_{\bm{\mathrm{i}}}$, where
$|\Psi(\tau)\rangle=\sum_{\bm{\mathrm{i}}}c_{\bm{\mathrm{i}}}(\tau)|D_{\bm{\mathrm{i}}}\rangle$:

$c_{\bm{\mathrm{i}}}(\tau+\delta\tau)=c_{\bm{\mathrm{i}}}(\tau)-\delta\tau[(H_{\bm{\mathrm{ii}}}-S)c_{\bm{\mathrm{i}}}(\tau)+\sum_{\bm{\mathrm{j}}\neq\bm{\mathrm{i}}}H_{\bm{\mathrm{ij}}}c_{\bm{\mathrm{j}}}(\tau)],$
(6)

where $H_{\bm{\mathrm{ij}}}=\langle
D_{\bm{\mathrm{i}}}|\hat{H}|D_{\bm{\mathrm{j}}}\rangle$. This equation can be
viewed as describing the population dynamics of particles placed on the
different determinants and may be modelled by a stochastic process composed of
three steps:11

  * •

Spawning: Given determinant $D_{\bm{\mathrm{i}}}$, generate new particles on
determinant $D_{\bm{\mathrm{j}}}$ with probability ${p\propto\delta\tau
H_{\bm{\mathrm{ij}}}c_{\bm{\mathrm{j}}}(\tau)}$.

  * •

Death: Given determinant $D_{\bm{\mathrm{i}}}$, generate new particles on
determinant $D_{\bm{\mathrm{i}}}$ with probability
${p\propto\delta\tau(H_{\bm{\mathrm{ii}}}-S)c_{\bm{\mathrm{i}}}(\tau)}$

  * •

Annihilation: For a given determinant, cancel out particles carrying opposite
signs.

The formulation of CCMC closely matches that of FCIQMC, with the difference
that instead of residing on determinants, walkers reside on excitors,
$\hat{a}_{\bm{\mathrm{n}}}$, defined as
$\hat{a}_{\bm{\mathrm{n}}}|D_{\bm{\mathrm{0}}}\rangle=\pm|D_{\bm{\mathrm{n}}}\rangle$.
Replacing the FCI wavefunction by the coupled cluster ansatz in Equation 4 and
left-multiplying by $\langle D_{\bm{\mathrm{i}}}|$ gives

$\begin{split}\langle
D_{\bm{\mathrm{i}}}|\Psi^{\mathrm{CC}}(\tau+\delta\tau)\rangle=&\langle
D_{\bm{\mathrm{i}}}|\Psi^{\mathrm{CC}}(\tau)\rangle\\\ -&\delta\tau\langle
D_{\bm{\mathrm{i}}}|(\hat{H}-S)|\Psi^{\mathrm{CC}}(\tau)\rangle\end{split}$
(7)

The coupled cluster ansatz parametrises the wavefunction with cluster
amplitudes in a non-linear fashion. The mapping of CI coefficients to cluster
amplitudes can be done by a simple projection, which reveals contributions
from multiple clusters. For example in a CCSD wavefunction (i.e.,
${\hat{T}}={\hat{T}}_{1}+{\hat{T}}_{2}$)

$\langle
D_{ij}^{ab}|e^{{\hat{T}}}D_{\bm{\mathrm{0}}}\rangle=t_{ij}^{ab}+t_{i}^{a}t_{j}^{b}-t_{i}^{b}t_{j}^{a}$
(8)

with the negative sign arising from the fact that
${\hat{a}_{b}^{\dagger}\hat{a}_{i}\hat{a}_{a}^{\dagger}\hat{a}_{j}=-\hat{a}_{a}^{\dagger}\hat{a}_{i}\hat{a}_{b}^{\dagger}\hat{a}_{j}}$,
due to the anti-commutation relations of the second-quantised creation and
annihilation operators 14. Terms like $t_{ij}^{ab}$ are known as non-composite
cluster amplitudes, and the rest as composite cluster amplitudes. Here we make
the approximation that composite clusters have much smaller contributions than
non-composite ones, their changes will be negligible per time step, and hence
we can cancel out the $\mathcal{O}(\hat{T}^{2})$ contributions on both sides
to write

$t_{\bm{\mathrm{i}}}(\tau+\delta\tau)\rangle=t_{\bm{\mathrm{i}}}(\tau)-\delta\tau\langle
D_{\bm{\mathrm{i}}}|(\hat{H}-S)|\Psi^{\mathrm{CC}}(\tau)\rangle$ (9)

Compared to FCIQMC, an additional step needs to be performed for each Monte
Carlo iteration: the sampling of the exponential ansatz. For $N_{\mathrm{ex}}$
total walkers, also called excips in CCMC, $\mathcal{O}(N_{\mathrm{ex}})$
clusters are formed randomly by combining present excitors according to
specific biasing rules 15.

Finally, the intermediate normalisation 16 of the wavefunction is redefined to
give the CCMC ansatz:

$|\Psi_{\text{CCMC}}\rangle=N_{0}e^{{\hat{T}}/N_{0}}|D_{\bm{\mathrm{0}}}\rangle,$
(10)

which introduces the reference population as a new independent variable,
solving the problem that Equation 7 does not lead to a viable update equation
for $\langle D_{\bm{\mathrm{i}}}|=\langle D_{\bm{\mathrm{0}}}|$.

A multi-reference formulation of the CCMC algorithm (MR-CCMC) 17 has been
implemented, retaining a single-reference formalism, in common with the so-
called SRMRCC methods in Ref. 18. The flexibility of the CCMC algorithm allows
this multireference approach to be implemented with minimal code changes to
the single-reference algorithm, bypassing what would be a challenging process
in deterministic methods. Essentially, for a coupled cluster truncation level
$P$, the algorithm allows any number of determinants to become a ‘secondary
reference’, stores excitors that are within $P$ excitations from any
references (instead of just the HF determinant), and allows clusters to form
that are within $P+2$ excitations from any reference. The set of references is
commonly known as the model or reference space. To summarise, the algorithmic
modifications relative to single-reference CCMC are:

  * •

Store all the secondary references in some searchable data structure, and
additionally store the highest excitation level from the reference determinant
among the secondary references, $k_{\text{max}}$.

  * •

Cluster expansion: allow clusters with excitation level of up to
$k_{\text{max}}+P+2$ to form, instead of $P+2$ in the single-reference case.
Discard those that are not $P+2$ excitations away from some reference
determinant.

  * •

Spawning: for a randomly generated spawnee (i.e., $\langle D_{\bm{j}}|$),
check that it is within $P$ excitations of any secondary references.

  * •

Cloning/death: allow death on excitors that are within $P$ excitations from
any secondary references.

While this MR-CCMC method can treat systems which conventional single-
reference CC struggles with, this comes at an increased computational cost.
Comparing contributor excitation levels to all references becomes expensive as
the number of references grows, particularly when the contributor turns out to
lie outside of the desired space. Therefore, non-trivial computational effort
is expended on attempts that will not contribute to the overall estimators and
propagation, while also making successful steps more expensive than their
single-reference equivalents. In the rest of this paper, we will first
introduce the wall-Chebyshev projector, and show that it can be applied to
(mr-)CCMC and FCIQMC to reduce the number of times the Hamiltonian needs to be
applied to reach statistical convergence, thereby reducing the computational
cost. We will then introduce a suite of further modifications to the MR-CCMC
algorithm that accelerates the handling of the model space. We apply the
resulting algorithm to several benchmark systems to show that it offers a way
for mr-CCMC to become a competitive alternative to FCIQMC even in the strong
correlation regime.

##  2 The wall-Chebyshev projector

###  2.1 Motivation and theory

In all projector-based QMC methods, such as FCIQMC and CCMC, the linear
projector (Equation 4) is used. The first-order Taylor expansion turns out to
be a very reasonable approximation, since we demonstrate in Appendix A that
there is no benefit whatsoever in going to higher orders of the Taylor
expansion of the exponential projector. However, this does not mean that one
cannot devise more efficient projectors. An example is a projector based on a
Chebyshev expansion of the wall function, which was first proposed in Ref. 19
in the context of a deterministic projector-based selected CI algorithm.

The wall function is given by

${\text{wall}}(x)=\begin{cases}\infty,\ x<0\\\ 1,\ x=0\\\ 0,\ x>0\end{cases},$
(11)

and is physically motivated as the infinite time limit of the exponential
projector:

$\lim_{\tau\rightarrow\infty}e^{-\tau(x-S)}=\text{wall}(x-S),$ (12)

which can map any trial wave function $|\Phi_{0}\rangle$ to the exact ground
state $|\Psi_{0}\rangle$, if $\langle\Phi_{0}|\Psi_{0}\rangle\neq 0$ and
$E_{0}\leq S<E_{1}$.

While a Taylor series expansion does not exist for the discontinuous wall
function, an expansion in Chebyshev polynomials, like a Fourier expansion, is
trivial. The Chebyshev polynomials of the first kind, defined as
$T_{n}(\cos(\theta))=\cos(n\theta)$, form an orthogonal basis (with metric
$(1-x^{2})^{-1/2}$) for functions defined over $x\in[-1,1]$:

$\int^{1}_{-1}T_{m}(x)T_{n}(x)(1-x^{2})^{-1/2}\dl
x=\delta_{mn}\pi/(2-\delta_{m0}).$ (13)

To facilitate the following discussion, we define the spectral range, $R$, of
a Hamiltonian as $R=E_{N-1}-E_{0}$, where $E_{i}$ is the $i$-th eigenvalue of
the full Hamiltonian and $N$ is the size of the Hilbert space. Furthermore,
our energy range $\epsilon\in[E_{0},E_{N-1}]$ requires the application of an
affine transformation to the Chebyshev polynomials

$\epsilon=E_{0}+\frac{R}{2}(x+1),\ x\in[-1,1].$ (14)

We show in Appendix B that the $m$-th order Chebyshev expansion of the wall
function is

$g^{\text{wall-
Ch}}_{m}(\epsilon)=\frac{1}{1+2m}\sum_{k=0}^{m}(2-\delta_{k0})T_{k}(-x).$ (15)

For illustration purposes, we plot several orders of Chebyshev expansion in
Figure 1, where we can also observe the monotonic divergence to $+\infty$ for
$\epsilon<E_{0}$. The other tail also diverges to $\pm\infty$ depending on the
parity of the order $m$.

In this instance, the nodes of Equation 15 are analytically known (see
derivation in Appendix B) as

$a_{\nu}=E_{0}+\frac{R}{2}\left(1-\cos\frac{\nu\pi}{m+1/2}\right).$ (16)

This allows us to decompose the $m$-th order projector into a product of $m$
linear projectors, each with their own weight that ensures $g^{\text{wall-
Ch}}_{m}(E_{0})=1$:

$g^{\text{wall-Ch}}_{m}(\epsilon)=\prod_{\nu=1}^{m}\frac{\epsilon-
a_{\nu}}{E_{0}-a_{\nu}}.$ (17)

A decomposition for a fifth order Chebyshev expansion of the wall function can
be seen in Figure 2.

Fig. 1: The Chebyshev expansions of the wall function in an arbitrary range of
$[-75,5]$, compared to the linear projector with the maximal time step of
$\delta\tau=0.025$, and its corresponding exponential projector. Fig. 2: The
fifth order Chebyshev expansion of the wall function, shown here to decompose
into a product of $5$ linear projectors, each with their own effective time
steps.

###  2.2 Application to FCIQMC and CCMC

In FCIQMC and CCMC, the lowest eigenvalue estimate is the shift, $S$, and the
upper spectral bound can be a constant, estimated from the Gershgorin circle
theorem 20 as

$\widetilde{E}_{N-1}=H_{N-1,N-1}+{{{\sum}^{\prime}_{{\bm{\mathrm{j}}}\in\\{\mathbf{S},\mathbf{D}\\}}}}H_{N-1,{\bm{\mathrm{j}}}},$
(18)

where the sum is over all determinants connected to the highest determinant
(singles and doubles), and the ‘′’ restricts it to ${\bm{\mathrm{j}}\neq
N-1}$.

The action of the wall-Chebyshev projector on
${|\Psi^{(n,0)}\rangle=[g^{\text{wall-Ch}}(\hat{H})]^{n}|\Phi\rangle}$ is

$|\Psi^{(n+1,0)}\rangle=g^{\text{wall-Ch}}(\hat{H})|\Psi^{(n,0)}\rangle,$ (19)

which gives the wavefunction after $n+1$ applications of the projector. We can
additionally define the ‘intermediate’ wavefunctions as

$|\Psi^{(n,\mu)}\rangle=\left[\prod_{\nu=0}^{\mu}\frac{\hat{H}-a_{\nu}}{S-a_{\nu}}\right]|\Psi^{(n,0)}\rangle.$
(20)

We are now ready to derive the update equations for FCIQMC and CCMC. We start
with the slightly more involved derivation for CCMC . Projecting these
intermediate wavefunctions onto determinants, we have

$\langle
D_{\bm{\mathrm{i}}}|\Psi^{(n,\nu+1)}\rangle=\frac{1}{S-a_{\nu}}\langle
D_{\bm{\mathrm{i}}}|{\hat{H}-a_{\nu}}|\Psi^{(n,\nu)}\rangle$ (21)

It is important now to distinguish between $\tilde{t}_{\bm{\mathrm{i}}}$, the
projection of a wavefunction onto determinant $D_{\bm{\mathrm{i}}}$, and the
corresponding excitor amplitude, $t_{\bm{\mathrm{i}}}$, with the former
including unconnected (‘composite’) contributions. At convergence,

$\displaystyle\begin{split}\tilde{t}_{\bm{\mathrm{i}}}=-\frac{1}{a_{\nu}-S}\left[\sum_{{\bm{\mathrm{j}}}\neq{\bm{\mathrm{i}}}}H_{\bm{\mathrm{ij}}}\tilde{t}_{\bm{\mathrm{j}}}+(H_{\bm{\mathrm{ii}}}-a_{\nu})\tilde{t}_{\bm{\mathrm{i}}}\right]&\\\
\tilde{t}_{\bm{\mathrm{i}}}-t_{\bm{\mathrm{i}}}+t_{\bm{\mathrm{i}}}=-\frac{1}{a_{\nu}-S}\left[\sum_{{\bm{\mathrm{j}}}\neq{\bm{\mathrm{i}}}}H_{\bm{\mathrm{ij}}}\tilde{t}_{\bm{\mathrm{j}}}+(H_{\bm{\mathrm{ii}}}-a_{\nu})\tilde{t}_{\bm{\mathrm{i}}}\right]&\\\
t_{\bm{\mathrm{i}}}=t_{\bm{\mathrm{i}}}-\frac{1}{a_{\nu}-S}\left[\sum_{\bm{\mathrm{j}}\neq\bm{\mathrm{i}}}H_{\bm{\mathrm{ij}}}\tilde{t}_{\bm{\mathrm{j}}}+(H_{\bm{\mathrm{ii}}}-S)\tilde{t}_{\bm{\mathrm{i}}}\right]&\end{split}$
(22)

We may now convert the last equation into an update step,

$t_{\bm{\mathrm{i}}}(\tau+\delta\tau)=t_{\bm{\mathrm{i}}}(\tau)-\frac{1}{a_{\nu}-S}\left[\sum_{\bm{\mathrm{j}}\neq\bm{\mathrm{i}}}H_{\bm{\mathrm{ij}}}\tilde{t}_{\bm{\mathrm{j}}}(\tau)+(H_{\bm{\mathrm{ii}}}-S)\tilde{t}_{\bm{\mathrm{i}}}(\tau)\right]$
(23)

Comparing with the original update equations, which are given by

$t_{\bm{\mathrm{i}}}(\tau+\delta\tau)=t_{\bm{\mathrm{i}}}(\tau)-\delta\tau\left[\sum_{\bm{\mathrm{j}}\neq\bm{\mathrm{i}}}H_{\bm{\mathrm{ij}}}\tilde{t}_{\bm{\mathrm{j}}}(\tau)+(H_{\bm{\mathrm{jj}}}-S)\tilde{t}_{\bm{\mathrm{i}}}(\tau)\right],$
(24)

we reach the conclusion that the necessary modifications are

  1. 1.

Setting $\delta\tau=1$

  2. 2.

Applying the $m$ constituent linear projectors in the $m$-th order wall-
Chebyshev projector. For linear projector ${\nu\in\\{1,\dots,m\\}}$, scale
Hamiltonian elements in spawning and death by $1/(a_{\nu}-S)$ (‘Chebyshev
weights’).

The same analysis can be performed on FCIQMC, without the complication of
composite amplitudes, to obtain a similar set of update equations:

$c_{\bm{\mathrm{i}}}=c_{\bm{\mathrm{i}}}-\frac{1}{a_{\nu}-S}\left[\sum_{\bm{\mathrm{j}}\neq\bm{\mathrm{i}}}H_{\bm{\mathrm{ij}}}c_{\bm{\mathrm{j}}}+(H_{\bm{\mathrm{ii}}}-S)c_{\bm{\mathrm{i}}}\right].$
(25)

In terms of implementation, the two sets of update equations are nearly
identical, and can share the same code in large parts.

Analysis of the asymptotic rate of convergence (see Appendix C) shows that the
theoretical speedup of an order $m$ wall-Chebyshev projector relative to the
linear projector with largest allowed $\delta\tau$ is $(m+1)/3$. 19 Due to
blooms, the largest $\delta\tau$ is never reached in the conventional
propagator, so real speedups are expected to be larger.

###  2.3 The shift update procedure

The original shift update equation for CCMC and FCIQMC is given by

$S^{(n+A)}=S^{(n)}-\frac{\zeta}{A\delta\tau}\ln\left(\frac{N_{\mathrm{w}}^{(n+A)}}{N_{\mathrm{w}}^{(n)}}\right),$
(26)

where the update is performed every $A$ time steps, $\zeta$ is the shift
damping parameter, and $N_{\mathrm{w}}$ is the total number of walkers.

Due to the time step $\delta\tau$ being set to unity, the shift update
procedure is expected to become rather unresponsive to the changes in particle
population. As a consequence, there can be vastly uncontained spawning,
unchecked by the lower-than-expected deaths, resulting in unmanageable
population growths. To remedy this, initially, a scaled update procedure was
experimented with, setting $A=1$:

$S^{(n+1)}=S^{(n)}-{\zeta}\sum_{\nu=1}^{\mu}[a_{\nu}^{(n)}-S^{(n)}]\ln\frac{N_{\mathrm{w}}^{(n,\nu)}}{N_{\mathrm{w}}^{(n,\nu-1)}},$
(27)

which seemed attractive as it reduces to Equation 21 in the first-order case
where the sum only contains one term or if all Chebyshev weights are the same.
However, this was not successful in reigning in the population growth. We
believe this is because the intermediate wavefunctions in Equation 20 are ill-
behaved due to being generated by an effective time step potentially larger
than $\tau_{\mathrm{max}}$. A series of population changes that start and end
at $N^{(n)}$ and $N^{(n+1)}$ respectively can produce very different values of
shift update in Equation 27, and the shift update produced is very sensitive
to the unreliable intermediate values. Hence the population information from
these intermediate wavefunctions should not be used.

Another procedure that was more successful was to decrease the damping (by
increasing $\zeta$) of the shift updates, causing the shift to be more
responsive to the changes in populations, which in turn helps stabilise the
population. We also found it helpful to use the improved shift update
procedure outlined in Ref. 21, where an additional term is added to the shift-
update procedure:

$S^{(n+A)}\leftarrow-\frac{\xi}{A\delta\tau}\ln\left(\frac{N_{\mathrm{w}}^{(n+A)}}{N_{t}}\right),$
(28)

where $\xi$ is the ‘forcing strength’, and $N_{t}$ is the target population.
This has the effect of additionally stabilising the population by ‘pinning’ it
to the pre-set target population. A further proposal from the same paper,
arising from an argument from a scalar model of population dynamics, is for
critical damping to be achieved by setting $\xi=\zeta^{2}/4$. This is also
found to be helpful. Altogether, these modifications result in greatly
improved population control and we were able to obtain dynamics that can be
used in a reblocking analysis, as shown in Figure 8, for example.

In practice, we have also found that with increasing order of Chebyshev
projector, a larger target population is usually needed, otherwise the
calculation may exhibit a sign-problem-like divergence. This may be attributed
to the larger effective time steps that the higher order projectors use and is
documented elsewhere, for example, see Fig. 2 in Ref. 22.

##  3 Accelerating the MR-CCMC algorithm

###  3.1 Efficient cluster acceptance algorithm

In the spawning step of the MR-CCMC algorithm, we check that a spawnee is
wihtin $P$ excitations of any secondary reference. The same check needs to be
performed in the death step. The original MR-CCMC algorithm performed a linear
scan through the list of secondary references, which is clearly a
$\mathcal{O}(n_{\mathrm{ref}})$ operation, where $n_{\mathrm{ref}}$ is the
number of secondary references. The subroutine that decides whether a spawn is
accepted is the second most frequently called subroutine in the program, after
the excitation generator. Therefore, a linear search in this step can quickly
become prohibitively expensive in a moderate to large model space
($n_{\mathrm{ref}}>1000$, for example). We note that traditional data
structures and search algorithms, such as a binary search on a sorted list of
secondary references, would not work here, as the definition of ‘distance’ in
this case, i.e., excitation-rank, is non-Euclidean. A search algorithm in a
general metric space is therefore needed.

The excitation-rank distance between two Slater determinants is equivalent to
half the Hamming distance between the bit strings representing these two
determinants, and the Hamming distance is a well-known example of a discrete
metric 23. A data structure, known as the BK tree 24, is particularly well
suited for efficient searches in discrete metric spaces. The tree, an example
of which is given in Figure 3, is constructed only once at the beginning of
the calculation. Subsequently, the search can be performed in
$\mathcal{O}(\log n_{\mathrm{ref}})$ time using a recursive tree traversal
algorithm. The tree construction and search algorithms are pictured in Figure
4.

Fig. 3: The BK tree can conduct efficient nearest neighbour searches in a
discrete metric space, like the excitation rank. In this figure a BK tree
built from $20$ arbitrary determinants is shown. The topology of the tree is
not unique, and is dependent on the order the nodes were added to the tree.
Fig. 4: Flowcharts for the building and searching of a BK-tree.

###  3.2 Compression of the model space

Whereas many classical multireference coupled cluster (MRCC) methods work with
complete active spaces (CAS), the MR-CCMC algorithm as presented here is
highly flexible as to the shape of the model space, and as such can be
considered a general model space (GMS) method 18. This enables us to consider
arbitrary subsets of the CAS as the model space, and fine-tune the balance
between cost and accuracy. One of us has devised a compression method in this
vein 25. Here we briefly summarise its main thrust.

Two of us observed that 17 that for some $(ne,no)$ active space, the results
of a MR-CCSD calculation using all of the determinants in this active space as
references (i.e., a CAS MR-CCSD calculation) is qualitatively very similar to
the results of a MR-CCSD…$m$ calculation, where $m=n/2$, using only the
‘bottom’ and ‘top’ determinants (i.e., the aufbau and anti-aufbau determinant
respectively) of the CAS as the references. We term the latter calculation as
‘2r-CCSD…$m$’. Using this observation, we aim to algorithmically generate only
those determinants that are in the Hilbert spaces of both the CAS MR-CCSD and
2r-CCSD…$m$ calculations, which should provide us with a compressed set of
reference determinants that captures the most important determinants in the
CAS. It was shown that this set of determinants can be generated by
enumerating determinants of up to $(m-2)$–fold excitations from the bottom and
top determinants.

##  4 Computational details

###  4.1 Basis sets and point group symmetries

In this work we will study the carbon and beryllium dimers. For these systems,
Dunning’s cc-pVXZ basis sets are used 26. The required electron integrals are
generated by the Psi4 27 and PySCF 28 packages. The electron integrals are
generated in the $D_{2h}$ point group symmetry and transformed into the basis
of $\hat{L}_{z}$ eigenfunctions based on the TransLz.f90 script provided in
the NECI package 29, which we re-wrote to interface with PySCF. Unless
otherwise noted, the core $1s$ electrons of C2 are frozen. The ‘heat bath’
excitation generator 30 is used whenever possible, otherwise the renormalised
excitation generator 31 is used.

The use of the $\hat{L}_{z}$ eigenfunctions helps not only further reduce the
size of the relevant symmetry sector, but also helps distinguish low-lying
states that would descend to the same irreducible representation in $D_{2h}$.
For C2, this would be the ${}^{1}\Sigma_{g}^{+}$ state and the
${}^{1}\Delta_{g}$ states, which both descend into the ${}^{1}A_{g}$ state in
$D_{2h}$. The two states approach and cross each other 32, 33, which would
prove challenging, if not impossible, to distinguish in $D_{2h}$.

When employing the wall-Chebyshev projector, the upper spectral range estimate
obtained from the Gershgorin theorem (Equation 18) is scaled up by $10\%$ by
default to guarantee an upper bound on the spectral width of the Hamiltonian.

Quantum Monte Carlo calculations are carried out using the HANDE-QMC
package34.

###  4.2 Symmetry screening

The full $(8e,8o)$ CAS model space of C2 contains $(^{8}C_{4})^{2}=4900$ model
space determinants that have $M_{s}=0$. However, like most classical quantum
chemistry programs, HANDE performs calculations strictly in the symmetry
sector requested by the user, which in our case is the totally symmetric
irreducible representation, unless otherwise specified. This brings into
question the effect, if any, of including secondary references outside the
symmetry sector under consideration. We first note that the secondary
references are only used to define the shape of the spawning/death space, and
are not required to be occupied. It is conceivable that there are
determinants/excitors in the correct symmetry sector that are included in the
spawning space only by virtue of being connected (but of course not coupled by
the Hamiltonian) to a wrong-symmetry secondary reference. To quantify the
effect of these wrong-symmetry secondary references, we implemented a simple
symmetry-screening switch in HANDE, which when turned on, only processes
secondary references of the currently considered symmetry.

##  5 Results and discussion

###  5.1 BK-tree

For a C2 system with a full $(8e,8o)$ CAS as the model space without symmetry
screening ($4900$ references that preserve the $M_{s}=0$ symmetry), the BK-
tree search is benchmarked against a naïve linear search, which loops over all
secondary references and terminates when one of the references is within $P$
excitations of the target determinant. The validity of the BK-tree search is
separately established by performing a normal calculation with either search
algorithm using the same random number generator seed, and asserting that the
results are the same. Benchmarking results are given in Table 1.

Table 1: Timing comparison between the BK tree and naive search algorithms, for C2 using a $(8e,8o)$ CAS as the model space for a multireference CCMCSD calculation. | Overall timing / s | Time per spawning attempt / $\mu$s
---|---|---
BK tree | $809.28$ | $12.761$
Linear | $5995.16$ | $94.533$

An apparent $8\times$ speedup is observed. Without performing a full profiling
study, the actual reduction in time cost of the acceptance search is expected
to be greater than $8\times$ as a complex series of operations is performed
per spawning attempt on top of the acceptance search.

###  5.2 Model space compression

For the $(8e,8o)$ CAS used for C2, the compression method discussed in Section
3.2 yields a total of $722$ determinants in the compressed model space. Here
we show the results for the C2/cc-pVDZ system at separations of $0.9$ to $1.5$
Å. The performance of the compression scheme is shown in Figure 5. Here we
have employed the default quasi-Newton acceleration35 implemented in HANDE. We
can see that despite the almost $7$-fold reduction in the model space (and
consequently a similar reduction in computational cost), the errors are within
chemical accuracy ($1.6$ m$E_{h}$).

Fig. 5: The correlation energy for C2/cc-pVDZ at $r_{\mathrm{C-C}}=0.9-1.5$Å
with the compressed set of $721$ secondary references, relative to the full
CAS model space. We observe that, despite a $7$-fold reduction in the size of
the model space, the reductions in the correlation energy captured are much
smaller, making this an attractive trade-off. The stochastic error bars are
too small to be seen, due to the use of the semi-stochastic algorithm 25.

###  5.3 Symmetry screening

Results with symmetry screening turned on for the same C2 system are shown in
Figure 6. We can see that, although statistically significant, including the
symmetry-forbidden secondary references only captures a negligible amount of
additional correlation energy. The full set of $4900$ references are reduced
to $635$ symmetry-allowed references, while the compressed set of $722$
references, whose results are not shown here, are reduced to $109$ references.
Both of these bring about a roughly $8$-fold reduction of the model space
size.

Fig. 6: The differences in energy between the full set of secondary
references, $E_{\mathrm{full}}$, and the set of symmetry-screened secondary
references, $E_{\mathrm{sym}}$, in units of m$E_{h}$. The very small error
bars are due to the semi-stochastic algorithm 25.

###  5.4 Binding curve of the carbon dimer

We studied the ${}^{1}\Sigma_{g}^{+}$ state of the carbon dimer in the cc-pVDZ
basis. The carbon dimer is a challenging test case for electronic structure
methods, and the challenge is two-fold: firstly, as mentioned in Section 4.1,
the ${}^{1}\Sigma_{g}^{+}$ state becomes very nearly degenerate with the
exceptionally low-lying $\Delta_{g}$ state at stretched bond lengths, and both
states descend to the $A_{g}$ state in the commonly used $D_{2h}$ point group
symmetry, rendering it very challenging to distinguish both states without the
use of the $L_{z}$ symmetry, with one paper resorting to tracking individual
CI coefficients 32; secondly, there is an abundance of avoided crossings, and
specifically, the first excited ${}^{1}\Sigma_{g}^{+}$ state participates in
an avoided crossing with the ground state at a bond length of around $1.6$ Å
36, resulting in a change in the most highly weighted diabatic state (i.e.,
determinant). This makes MR-CCMC calculations based on RHF orbitals exhibit
long-imaginary-time instabilities for stretched geometries, which preclude
obtaining accurate estimators. The binding curve presented in Figure 7 used
the full $(8e,8o)$ CAS as the model space for a MR-CCMCSD calculation, with
orbitals obtained using PySCF from an $(8e,8o)$ state-average CASSCF
calculation over the lowest three ${}^{1}A_{g}$ states (in the $D_{2h}$ point
group). The orbital coefficients are still tagged with their corresponding
$D_{\infty h}$ irreducible representations, enabling us to perform the $L_{z}$
transformation. We ensured that the $\pi_{u}$ manifold was included in the
reference determinant that generates the Hilbert space and secondary
references.

Fig. 7: The binding curve of the ${}^{1}\Sigma_{g}^{+}$ state of C2/cc-pVDZ in
the range of $0.9$-$2.8$ Å separation. MR-CCMC calculations are based on
CASSCF orbitals and use an $(8e,8o)$ CAS as a model space, with clusters
truncated at the double excitation level.The FCIQMC result is from Ref. 37,
and the DMRG results are from 32.

The non-parallelity errors (NPE), defined here as the difference between the
maximal and minimal deviation from the DMRG energies (which are interpolated
with a cubic spline using the scipy.interpolate.CubeSpline class) or the
FCIQMC energies, and a summary in shown in Table 2.

| NPE (m$E_{\mathrm{h}}$) | Max AD (m$E_{\mathrm{h}}$) | Min AD (m$E_{\mathrm{h}}$)
---|---|---|---
FCIQMC | 10.0 | 5.8 (2.8) | 0.7 (2.4)
DMRG | 10.4 | 7.1 (1.9) | 0.4 (2.5)
Table 2: The non-parallelity error, maximal, and minimal deviations of the
carbon dimer binding curve calculated with MR-CCMCSD using a $(8e,8o)$ CAS
model space, compared to FCIQMC and DMRG. The numbers in parenthesis indicate
the bond length (in angstrom) at which the maximal/minimal absolute deviations
occur.

###  5.5 Chebyshev Results

Figure 8 shows an example of the power of the Chebyshev propagation, applied
to an MR-CCMCSD calculation for C2. The shoulder height is reached after
around $50$ iterations, and the calculation is equilibriated essentially
instantaneously, which means all that is left to do is collecting statistics.
On an Intel(R) Xeon(R) E5-2650 v2 CPU, the calculation in the figure was run
for only $2$ hours on $6$ physical cores. Without the Chebyshev projector, the
same calculation takes around $24$ hours with $12$ physical cores to give the
same statistical error bar.

Fig. 8: A fifth order Chebyshev projector applied to MR-CCMCSD for C2/cc-pVDZ
at $1.2$ Å separation, with a full $(8e,8o)$ CAS as the model space, with a
target population of $2\times 10^{6}$ and a shift damping parameter of $0.5$,
as compared to a linear projector with $\delta\tau=0.001$, a target population
of $1\times 10^{6}$, and a shift damping parameter of $0.05$. Both
calculations are carried out with a two-stage harmonic forcing shift update
scheme. The inset shows that the projected energy only barely stabilises
around the true value at the end of the calculation using the linear
projector. This means that the reblocking algorithm will disregard essentially
all of the iterations shown here, due to the shift and the projected energy
not agreeing.

Part of the benefit of our Chebyshev propagator algorithm is the automatic
determination of the effective time steps, or ‘Chebyshev weights’, via the
Gershgorin circle theorem, which in the first-order expansion limit reduces to
an automatic way of choosing a good time step, instead of using trial and
error.

The following example shows Be2, a smaller, modestly multi-reference system,
treated in the cc-pVQZ basis set. Referring to Figure 9, compared with the
linear projector with a guessed $\delta\tau=0.002$, the second-order Chebyshev
expansion shows a clear speed-up in convergence. In fact, the Chebyshev
calculation took $66$ applications of the Hamiltonian to reach the target
population of $3\times 10^{6}$, whereas the linear projector took $3034$
applications to reach the same target population.

Fig. 9: The second order projector and the default linear projector at
$\delta\tau=0.002$ applied to MR-CCMCSD for the Be2/cc-pVQZ system at $2.5$ Å
separation, with a symmetry-screened $(4e,8o)$ (full $2s,2p$ valence) CAS as
the model space. Both calculations have a target population of $3\times
10^{6}$.

However, the second-order Chebyshev projector is expected to be as efficient
as the linear projector with the maximum allowed time step (see Section 2.2).
To provide a fairer comparison, we present in Figure 10 the first, second and
fourth order Chebyshev projector applied to the Be2/cc-pVQZ at $2.5$ Å system.
The first order Chebyshev projector is equivalent to a linear projector with
$\delta\tau=3/(E_{N-1}-E_{0})$, which is $2/3$ of $\delta\tau_{\text{max}}$,
but this maximal time step is commonly found to give rise to destabilising
population blooms, so the first order Chebyshev projector gives an estimate
for a practical time step to use.

Fig. 10: MR-CCMCSD for Be2/cc-pVQZ at $2.5$ Å separation with the $(4e,8o)$
CAS model space run with the first, second and fourth order Chebyshev
projectors. Only the time evolution of the shift parameters is shown for
clarity.

Here we can most clearly see the reduction in time needed to reach the
shoulder height. It is worth bearing in mind that the three calculations
require different target populations to stabilise, with the first and second
order projectors having target populations of $3\times 10^{6}$, and the fourth
order projector having a target population of $5\times 10^{6}$, which slightly
increases the number of iterations required to reach the target population.

An even clearer example of the relative time reduction is shown in Figure 11,
where we simulate the effect of the first order Chebyshev projector by
compressing the imaginary time evolution of the linear projector with
$\delta\tau=0.001$ in Figure 8 accordingly.

Fig. 11: C2/cc-pVDZ at $1.2$ Å separation, with a full $(8e,8o)$ CAS model
space, run with the first (simulated) and fifth order Chebyshev projectors.
The latter requires significantly fewer Hamiltonian applications to reach the
target population.

As a demonstration of the applicability of the Chebyshev projector in
different correlation regimes, we have computed the binding curve of the
Be2/cc-pVTZ system using a $(4e,8o)$ CAS as a model space for MR-CCMCSD, using
the fifth-order Chebyshev projector. CCSD(T) and semistochastic heat-bath
configuration interaction with second-order perturbation correction (SHCI-PT2)
38 results are also shown. The CCSD(T) results are from Psi4,27, and the SHCI-
PT2 results are generated using the Dice plug-in to the PySCF package28. For
SHCI-PT2, the full (8e, 60o) space was correlated, and hence it can be
considered a surrogate for FCI results. The variational threshold was set to
$\epsilon_{1}=8\times 10^{-5}$ $E_{\mathrm{h}}$, and the PT2 threshold was set
to $\epsilon_{2}=1\times 10^{-8}$ $E_{\mathrm{h}}$, using $N_{d}=200$
deterministic determinants, with 5 PT2 iterations. Figure 12 shows the binding
curves. The cc-pVTZ basis is known to severely overbind the beryllium dimer,39
compared to the experimental value of $934.9\pm 2.5$ cm-1.40, 41, 42
Notwithstanding, the MR-CCMC method using the Chebyshev projector was able to
provide consistently better description of the binding curve than the ‘gold-
standard’ CCSD(T). The MR-CCMC results are close to but not qualitatively the
same as FCI results near the equilibrium, where static correlation
dominates,43, and are near-identical to FCI results in the stretched region,
where dynamic correlation dominates, and a compact coupled-cluster
representation of the wavefunction is beneficial.

Fig. 12: The binding curve of the beryllium dimer using the cc-pVTZ basis,
computed using the CCSD(T), MR-CCMCSD with a $(4e,8o)$ CAS, and SHCI-PT2
methods.

##  6 Conclusion

We present here a series of algorithmic changes that can be used to accelerate
the MR-CCMC algorithm in particular and QMC algorithms in general.

Specific to the MR-CCMC algorithm, we have introduced a BK-tree-based search
algorithm to verify whether proposed clusters and spawns are within the
accepted manifold for a given reference space. This reduces the scaling of
this step from $\mathcal{O}(n_{\mathrm{ref}})$ to $\mathcal{O}(\log
n_{\mathrm{ref}})$, which translates to an $8\times$ speed-up for the
molecular systems studied. We have also designed a compression method for the
reference space, which preserves only what we expect to be the most
significant reference determinants. This decreases the size of the space by
close to an order of magnitude. Finally, we have shown that only including
references that belong to the same symmetry sector as the desired solution is
also effective as a means to reduce the size of the reference space, while
introducing only negligible additional error to the results.

We have also developed a new projector based on the Chebyshev polynomial
expansion of the wall function, which significantly accelerates the
convergence of QMC calculations. In an example calculation on the Be2
molecule, this reduced the number of Hamiltonian applications needed to reach
the target population by two orders of magnitude. The wall-Chebyshev projector
is generally applicable to different QMC algorithms, including FCIQMC and
(mr-)CCMC approaches, so we believe that, together with many recent
developments in increasing the apparent scarcity of the Hamiltonian and
optimising the shift behaviour,44, 45, 21it is a promising step in expanding
the range of applications for these methods.

## Author contributions

ZZ, MAF and AJWT all contributed to the conceptualization of the project. ZZ
developed the necessary software and analysed the data, under supervision from
MAF and AJWT. The manuscript was written by ZZ and MAF, with input from all
authors.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

ZZ was in part funded by the U.S. Department of Energy under grant DE-
SC0024532. MAF thanks Corpus Christi College, Cambridge and the Cambridge
Trust for a studentship, as well as Peterhouse for funding through a Research
Fellowship. This work used the ARCHER2 UK National Supercomputing Service
(https://www.archer2.ac.uk).

## Notes and references

  * McMillan 1965 W. L. McMillan, _Phys. Rev._ , 1965, 138, A442.
  * Ceperley _et al._ 1977 D. Ceperley, G. V. Chester and M. H. Kalos, _Phys. Rev. B_ , 1977, 16, 3081.
  * Anderson 1975 J. B. Anderson, _J. Chem. Phys._ , 1975, 63, 1499–1503.
  * Anderson 1976 J. B. Anderson, _J. Chem. Phys._ , 1976, 65, 4121–4127.
  * Anderson 1980 J. B. Anderson, _J. Chem. Phys._ , 1980, 73, 3897.
  * Grossman 2002 J. C. Grossman, _J. Chem. Phys._ , 2002, 117, 1434–1440.
  * Nemec _et al._ 2010 N. Nemec, M. D. Towler and R. J. Needs, _J. Chem. Phys._ , 2010, 132, 034111.
  * Cox _et al._ 2014 S. J. Cox, M. D. Towler, D. Alfè and A. Michaelides, _J. Chem. Phys._ , 2014, 140, 174703.
  * Ganesh _et al._ 2014 P. Ganesh, J. Kim, C. Park, M. Yoon, F. A. Reboredo and P. R. C. Kent, _J. Chem. Theory Comput._ , 2014, 10, 5318–5323.
  * Della Pia _et al._ 2022 F. Della Pia, A. Zen, D. Alfè and A. Michaelides, _J. Chem. Phys._ , 2022, 157, 134701.
  * Booth _et al._ 2009 G. H. Booth, A. J. W. Thom and A. Alavi, _J. Chem. Phys._ , 2009, 131, 054106.
  * Thom 2010 A. J. W. Thom, _Phys. Rev. Lett._ , 2010, 105, 263004.
  * Wick 1954 G. C. Wick, _Phys. Rev._ , 1954, 96, 1124–1134.
  * Helgaker _et al._ 2014 T. Helgaker, P. Jorgensen and J. Olsen, _Molecular Electronic-Structure Theory_ , John Wiley & Sons, 2014.
  * Spencer and Thom 2016 J. S. Spencer and A. J. W. Thom, _J. Chem. Phys._ , 2016, 144, 084108\.
  * Shavitt and Bartlett 2009 I. Shavitt and R. J. Bartlett, _Many-Body Methods in Chemistry and Physics: MBPT and Coupled-Cluster Theory_ , Cambridge University Press, Cambridge, 2009.
  * Filip _et al._ 2019 M.-A. Filip, C. J. C. Scott and A. J. W. Thom, _J. Chem. Theory Comput._ , 2019, 15, 6625–6635.
  * Lyakh _et al._ 2012 D. I. Lyakh, M. Musiał, V. F. Lotrich and R. J. Bartlett, _Chem. Rev._ , 2012, 112, 182–243.
  * Zhang and Evangelista 2016 T. Zhang and F. A. Evangelista, _J. Chem. Theory Comput._ , 2016, 12, 4326–4337.
  * Geršgorin 1931 S. Geršgorin, _Bulletin de l’Académie des Sciences de l’URSS. Classe des sciences mathématiques et na_ , 1931, 749–754.
  * Yang _et al._ 2020 M. Yang, E. Pahl and J. Brand, _J. Chem. Phys._ , 2020, 153, 174103\.
  * Vigor _et al._ 2016 W. A. Vigor, J. S. Spencer, M. J. Bearpark and A. J. W. Thom, _J. Chem. Phys._ , 2016, 144, 094110.
  * Hamming 1950 R. W. Hamming, _The Bell System Technical Journal_ , 1950, 29, 147–160.
  * Burkhard and Keller 1973 W. A. Burkhard and R. M. Keller, _Commun. ACM_ , 1973, 16, 230–236.
  * Zhao 2022 Z. Zhao, _M.Sc. thesis_ , University of Cambridge, Cambridge, UK, 2022\.
  * Dunning 1989 T. H. Dunning, _J. Chem. Phys._ , 1989, 90, 1007–1023.
  * Smith _et al._ 2020 D. G. A. Smith, L. A. Burns, A. C. Simmonett, R. M. Parrish, M. C. Schieber, R. Galvelis, P. Kraus, H. Kruse, R. Di Remigio, A. Alenaizan, A. M. James, S. Lehtola, J. P. Misiewicz, M. Scheurer, R. A. Shaw, J. B. Schriber, Y. Xie, Z. L. Glick, D. A. Sirianni, J. S. O’Brien, J. M. Waldrop, A. Kumar, E. G. Hohenstein, B. P. Pritchard, B. R. Brooks, H. F. Schaefer, A. Y. Sokolov, K. Patkowski, A. E. DePrince, U. Bozkaya, R. A. King, F. A. Evangelista, J. M. Turney, T. D. Crawford and C. D. Sherrill, _J. Chem. Phys._ , 2020, 152, 184108.
  * Sun _et al._ 2018 Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova, S. Sharma, S. Wouters and G. K.-L. Chan, _WIREs Computational Molecular Science_ , 2018, 8, e1340.
  * Guther _et al._ 2020 K. Guther, R. J. Anderson, N. S. Blunt, N. A. Bogdanov, D. Cleland, N. Dattani, W. Dobrautz, K. Ghanem, P. Jeszenszki, N. Liebermann, G. L. Manni, A. Y. Lozovoi, H. Luo, D. Ma, F. Merz, C. Overy, M. Rampp, P. K. Samanta, L. R. Schwarz, J. J. Shepherd, S. D. Smart, E. Vitale, O. Weser, G. H. Booth and A. Alavi, _J. Chem. Phys._ , 2020, 153, 034107.
  * Holmes _et al._ 2016 A. A. Holmes, H. J. Changlani and C. J. Umrigar, _J. Chem. Theory Comput._ , 2016, 12, 1561–1571.
  * Booth _et al._ 2014 G. H. Booth, S. D. Smart and A. Alavi, _Mol. Phys._ , 2014, 112, 1855–1869.
  * Wouters _et al._ 2014 S. Wouters, W. Poelmans, P. W. Ayers and D. Van Neck, _Comput. Phys. Commun._ , 2014, 185, 1501–1514.
  * Sharma 2015 S. Sharma, _J. Chem. Phys._ , 2015, 142, 024107.
  * Spencer _et al._ 2019 J. S. Spencer, N. S. Blunt, S. Choi, J. Etrych, M.-A. Filip, W. M. C. Foulkes, R. S. T. Franklin, W. J. Handley, F. D. Malone, V. A. Neufeld, R. Di Remigio, T. W. Rogers, C. J. C. Scott, J. J. Shepherd, W. A. Vigor, J. Weston, R. Xu and A. J. W. Thom, _J. Chem. Theory Comput._ , 2019, 15, 1728–1742.
  * Neufeld and Thom 2020 V. A. Neufeld and A. J. W. Thom, _J. Chem. Theory Comput._ , 2020, 16, 1503–1510.
  * Varandas 2008 A. J. C. Varandas, _J. Chem. Phys._ , 2008, 129, 234103.
  * Booth _et al._ 2011 G. H. Booth, D. Cleland, A. J. W. Thom and A. Alavi, _J. Chem. Phys._ , 2011, 135, 084104.
  * Sharma _et al._ 2017 S. Sharma, A. A. Holmes, G. Jeanmairet, A. Alavi and C. J. Umrigar, _J. Chem. Theory Comput._ , 2017, 13, 1595–1604.
  * Guther _et al._ 2021 K. Guther, A. J. Cohen, H. Luo and A. Alavi, _J. Chem. Phys._ , 2021, 155, 011102.
  * Merritt _et al._ 2009 J. M. Merritt, V. E. Bondybey and M. C. Heaven, _Science_ , 2009, 324, 1548–1551.
  * Meshkov _et al._ 2014 V. V. Meshkov, A. V. Stolyarov, M. C. Heaven, C. Haugen and R. J. LeRoy, _J. Chem. Phys._ , 2014, 140, 064315.
  * Patkowski _et al._ 2009 K. Patkowski, V. Špirko and K. Szalewicz, _Science_ , 2009, 326, 1382–1384.
  * El Khatib _et al._ 2014 M. El Khatib, G. L. Bendazzoli, S. Evangelisti, W. Helal, T. Leininger, L. Tenti and C. Angeli, _J. Phys. Chem. A_ , 2014, 118, 6664–6673.
  * Ghanem _et al._ 2019 K. Ghanem, A. Y. Lozovoi and A. Alavi, _J. Chem. Phys._ , 2019, 151, 224108.
  * Ghanem _et al._ 2020 K. Ghanem, K. Guther and A. Alavi, _J. Chem. Phys._ , 2020, 153, 224115\.
  * Ayoub 1980 R. G. Ayoub, _Archive for History of Exact Sciences_ , 1980, 23, 253–277.
  * Riley _et al._ 2006 K. F. Riley, M. P. Hobson and S. J. Bence, _Mathematical Methods for Physics and Engineering: A Comprehensive Guide_ , Cambridge University Press, 3rd edn., 2006.
  * Dirichlet 1829 P. G. L. Dirichlet, _Journal für die reine und angewandte Mathematik_ , 1829, 4, 157–169.
  * Kosloff and Tal-Ezer 1986 R. Kosloff and H. Tal-Ezer, _Chem. Phys. Lett._ , 1986, 127, 223–230.

##  Appendix A Higher Taylor expansions of the exponential projector

In Ref. 19 it was proposed that there is no gain in going to higher order
Taylor expansions of the exponential projector, because all orders of
expansion have $\gamma=\tau$ (see Appendix C). The conclusion is correct, but
for a more subtle reason that we will now explain. The $m$-th order Taylor
series expansion of $g^{\text{exp}}$ is

$\sum_{k=0}^{m}\frac{1}{k!}(-\tau)^{k}(x-E_{0})^{k}.$ (29)

Eq. 6 in Ref. 19 requires that $g(E_{N-1})<1$, so, defining the spectral range
$R=E_{N-1}-E_{0}$, we have

$\left|\sum_{k=0}^{m}\frac{1}{k!}(-\tau)^{k}R^{k}\right|<1.$ (30)

The $m=1$ case leads to the familiar requirement in DMC, FCIQMC and CCMC that

$\tau_{\text{max}}<\frac{2}{E_{{N-1}}-E_{0}}$ (31)

and solving Equation 30 numerically shows that higher order expansions lead to
larger maximum allowed $\tau$. In fact, asymptotically, $\tau_{\text{max}}$
increases linearly with a gradient of $1/e\approx 0.368$ (see Figure 13). So,
naïvely we can expect the efficiency to increase linearly with $m$. To prove
the linearity, we note that $\tau_{\text{max}}R>1$, and we can approximate
$\sum_{k=0}^{m}|(-x)^{k}/{k!}|$ ($x\equiv\tau R$) with the leading order term
$|(-x)^{m}|/{m!}$, and so we are left with

$\displaystyle|(-x)^{m}|$ $\displaystyle<m!$ (32a) $\displaystyle m\ln x$
$\displaystyle\lessapprox m\ln m-m$ (32b) $\displaystyle x$
$\displaystyle<me^{-1},$ (32c)

where we used the Stirling’s formula in the second line.

Fig. 13: Calculated and fitted $\tau_{\text{max}}$ as a function of the order
of Taylor expansion of the exponential projector. There is indeed no gain
whatsoever in going to the second order Taylor expansion, but there is in
going to yet higher orders.

Therein lies the real reason for not using higher order expansions: a naive
implementation requires $m(m+1)/2$ applications of the Hamiltonian per
projection, and even a factorised implementation would require $m$
applications per projection, not to mention the lack of closed forms for the
roots of the $m$-th order expansion, due to the Abel-Ruffini theorem46, which
shows that no analytical solution can exist for $m\geq 5$. In any case, the
overall efficiency stays at best constant. Therefore, the conclusion that no
gains can be made is correct, although a more tortuous argument is needed.

##  Appendix B Properties of the wall-Chebyshev projector

Assuming the entire spectral range is re-scaled such that $x\in[-1,1]$, where
$x=2(E-E_{0}))/R-1$, and $R=E_{N-1}-E_{0}$ is the spectral range or the
Hamiltonian, the Chebyshev expansion coefficients of the wall function is

$\begin{split}c_{k}=&\int_{-1}^{1}\text{wall}(x+1)T_{k}(x)(1-x^{2})^{-1/2}\dl
x\\\ \equiv&\int^{1}_{-1}\text{wall}(x+1)T_{k}(x)\delta(x+1)\dl x\\\
=&(2-\delta_{k0})T_{k}(-1)=(2-\delta_{k0})(-1)^{k},\end{split}$ (33)

where the second line uses the fact that the wall function is zero everywhere
but at the lower bound, so the weight function $(1-x^{2})^{-1/2}$ has the same
action as the delta function centred at $-1$, $\delta(x+1)$. The last equality
exploits a well-known identity of the Chebyshev polynomials 47. We can then
write the $m$-th order expansion as

$\begin{split}g^{\text{wall-
Ch}}_{m}(x)&=\frac{1}{1+2m}\sum_{k=0}^{m}(2-\delta_{k0})(-1)^{k}T_{k}(x)\\\
&=\frac{1}{1+2m}\sum_{k=0}^{m}(2-\delta_{k0})T_{k}(-x),\end{split}$ (34)

where the last equality exploits the fact that $T_{k}$ has the same parity as
$k$ 47, and we scale the sum such that $g^{\text{wall-Ch}}_{m}(-1)=1$.

To obtain an analytical expression for the zeroes of the wall-Chebyshev
projector, we use the trigonometric definition of the Chebyshev polynomials.
Inverting the sign of the argument in Equation 34, we have

$\begin{split}g^{\text{wall-
Ch}}_{m}(-\cos\theta)&\propto\sum_{k=0}^{m}(2-\delta_{k0})T_{k}(\cos\theta)\\\
&=1+2\left[\sum_{k=0}^{m}\cos(k\theta)\right]\\\
&=\frac{\sin\left[{\left(m+1/2\right)\theta}\right]}{\sin\left(\theta/2\right)},\end{split}$
(35)

where the last equality is the Dirichlet kernel 48. The zeroes of
$g^{\text{wall-Ch}}_{m}(x)$ are then transparently

$a_{\nu}=-\cos\left(\frac{\nu\pi}{m+1/2}\right),\ \nu=1,2,\dots,m,$ (36)

where the negative sign is from account for the sign inversion in Equation 35.
In an arbitrary spectral range other than $[-1,1]$, these zeroes are

$a_{\nu}=E_{0}+\frac{R}{2}\left(1-\cos\frac{\nu\pi}{m+1/2}\right).$ (37)

Knowing its zeroes, we can decompose $g^{\text{wall-Ch}}_{m}(x)$ into a
product of linear projectors:

$g^{\text{wall-Ch}}_{m}(x)=\prod_{\nu=1}^{m}\frac{x-a_{\nu}}{E_{0}-a_{\nu}},$
(38)

where the numerators ensure $g^{\text{wall-Ch}}_{m}(E_{0})=1$.

##  Appendix C Convergence properties of generators

We summarise here some important properties of generators. The asymptotic rate
of convergence of a propagator is dominated by the slowest-decaying
eigencomponent:

$\mu=\lim_{n\rightarrow\infty}=\frac{||\Psi^{(n+1)}-\Psi_{0}||}{||\Psi^{(n)}-\Psi_{0}||}=\max_{i}\left|\frac{g(E_{i})}{g(E_{0})}\right|.$
(39)

Zhang and Evangelista 19 suggested that, in the common case that the first
excited state is the slowest-decaying component, and that the first excited
energy is small compared to the spectral range of $\hat{H}$, the above can be
approximated as

$\mu\approx\left|1+(E_{1}-E_{0})g^{\prime}(E_{0})\right|\equiv\left|1-\alpha\gamma\right|,$
(40)

where $\gamma=-g^{\prime}(E_{0})$ is the convergence factor for $g$. We now
derive the relation given in Ref. 49 that $\gamma$ is approximately the number
of times, $n$, that $g$ needs to be applied to achieve an error in the norm,
$\epsilon=||\Psi^{(n)}-\Psi_{0}||$ to the $N$-th decimal place:

$\displaystyle\epsilon=10^{-N}$ $\displaystyle\approx(1-\alpha\gamma)^{n}$
(41a) $\displaystyle 10^{-N}$ $\displaystyle\approx e^{-n\alpha\gamma}$ (41b)
$\displaystyle n$ $\displaystyle=\frac{N\ln 10}{\alpha\gamma}\equiv\kappa
N\cdot\frac{1}{\gamma},$ (41c)

where $\kappa=\ln 10/(E_{1}-E_{0})$ is the convergence prefactor, which is
inversely proportional to the first excited energy gap.


# EF-Calib: Spatiotemporal Calibration of
Event- and Frame-Based Cameras
Using Continuous-Time Trajectories

Shaoan Wang, Zhanhua Xin, Yaoqing Hu, Dongyue Li, Mingzhu Zhu, and Junzhi Yu
This work was supported in part by the National Natural Science Foundation of
China under Grant T2121002 and Grant 62233001, and in part by the Beijing
Natural Science Foundation under Grant 2022MQ05. _(Corresponding author:
Junzhi Yu.)_ Shaoan Wang, Zhanhua Xin, Yaoqing Hu, Dongyue Li, and Junzhi Yu
are with the State Key Laboratory for Turbulence and Complex Systems,
Department of Advanced Manufacturing and Robotics, College of Engineering,
Peking University, Beijing 100871, China (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>junzhi.yu@ia.ac.cn).Mingzhu Zhu is with the
Department of Mechanical Engineering, Fuzhou University, Fuzhou 350000, China
(e-mail: mzz@fzu.edu.cn).

###### Abstract

Event camera, a bio-inspired asynchronous triggered camera, offers promising
prospects for fusion with frame-based cameras owing to its low latency and
high dynamic range. However, calibrating stereo vision systems that
incorporate both event and frame-based cameras remains a significant
challenge. In this letter, we present EF-Calib, a spatiotemporal calibration
framework for event- and frame-based cameras using continuous-time
trajectories. A novel calibration pattern applicable to both camera types and
the corresponding event recognition algorithm is proposed. Leveraging the
asynchronous nature of events, a derivable piece-wise B-spline to represent
camera pose continuously is introduced, enabling calibration for intrinsic
parameters, extrinsic parameters, and time offset, with analytical Jacobians
provided. Various experiments are carried out to evaluate the calibration
performance of EF-Calib, including calibration experiments for intrinsic
parameters, extrinsic parameters, and time offset. Experimental results show
that EF-Calib achieves the most accurate intrinsic parameters compared to
current SOTA, the close accuracy of the extrinsic parameters compared to the
frame-based results, and accurate time offset estimation. EF-Calib provides a
convenient and accurate toolbox for calibrating the system that fuses events
and frames. The code of this paper will also be open-sourced at:
https://github.com/wsakobe/EF-Calib.

###### Index Terms:

Event camera, spatiotemporal calibration, continuous-time trajectory, time
offset estimation.

##  I Introduction

Recent years, there has been a growing interest among researchers in a novel
bio-inspired camera called the event camera [1]. Abandoning the frame-
triggered concept of conventional cameras, each pixel of the event camera can
be considered as responding independently and asynchronously to changes in
illumination, resulting in an extremely low-latency and high-dynamic-range
response pattern. These advantages offer competitive prospects for event
cameras in areas such as robotics [2], autonomous driving [3], VR/AR [4], and
camera imaging [5].

Figure 1: Overview diagram of EF-Calib. (a) The novel calibration pattern
consists of the concentric circle and crosspoint. (b) The stereo vision system
consists of an event camera and a frame-based camera. (c) The calibration
process of EF-Calib.

However, due to the imaging principle of event cameras, they can only react to
changes in illumination, making it difficult to capture absolute amounts of
illumination and RGB values as frame-based cameras do. This limitation weakens
the ability of event cameras to perceive and understand the environment.
Therefore, many recent studies have attempted to fuse events with images in
order to fully utilize the unique advantages of both modalities, as
illustrated in Fig. 1. Some novel SLAM systems achieve more robust
localization under fast motion by fusing events and frames [2, 8, 7, 6]. In
addition, several studies are exploring how to fuse events and frames for
object detection in challenging environments [3, 9]. In recent years, some
event-centric datasets with multiple sensors, including frame-based cameras,
have also been widely proposed [12, 11, 10].

It is important to note that calibrating the intrinsic and extrinsic
parameters of each camera is an indispensable step in the context of multi-
camera fusion. Classical camera-to-camera calibration schemes typically
require time synchronization, followed by the acquisition of each camera’s
parameters through the synchronous acquisition of images of the calibration
board with different poses [13]. However, due to the asynchronous nature of
the event camera, it is difficult to combine multiple events into single
“frames” and time-synchronize them with an image. In addition, events are only
generated if there is relative motion, so the event camera cannot capture a
stationary calibration board. In summary, a new calibration framework must be
designed for the system including event- and frame-based cameras.

To address the aforementioned issues, this letter proposes a novel
spatiotemporal calibration framework for event- and frame-based cameras. To
the best of our knowledge, EF-Calib is the first calibration framework to
achieve joint calibration of event- and frame-based cameras without requiring
any time synchronization. The main contributions of this paper are as follows:

  1. 1.

A novel spatiotemporal calibration framework for event- and frame-based
cameras is proposed. This framework can jointly obtain the intrinsic and
extrinsic parameters, as well as the time offset without requiring any
hardware synchronization.

  2. 2.

Leveraging the asynchronous and low-latency properties of the event camera,
the framework introduces a continuous-time trajectory to optimize its motion
trajectory, facilitating arbitrary timestamp alignment with the frame-based
camera.

  3. 3.

Extensive experiments are conducted in diverse scenarios to validate the
proposed calibration framework. The results indicate that the framework
achieves accuracy close to that of frame-based camera calibration methods and
consistently calibrates the time offset between the cameras.

The rest of the letter is organized as follows. Sec. II summarizes the related
works. Sec. III presents the preliminaries of event-based vision and
continuous-time trajectory. Sec. IV introduces the calibration framework. Sec.
V evaluates the calibration performance from different aspects. At last, Sec.
VI presents the conclusion of this letter.

##  II Related Works

For geometric vision, camera calibration is particularly crucial as it serves
as the initial step in processing the input image signal, with the quality of
calibration often dictating the performance of subsequent tasks. Traditional
camera calibration has undergone significant evolution. The most prevalent
calibration method today involves capturing images of a calibration pattern
with a known size, such as a checkerboard, from various viewpoints to identify
corresponding feature points [13]. Subsequently, the intrinsic and distortion
parameters of each camera, as well as the extrinsic parameters between
cameras, are automatically calculated.

Nevertheless, applying this static and discrete calibration method to event
cameras, which are triggered by changes in illumination or relative motion,
presents challenges. Initially, many open-source event camera calibration
toolkits utilized a synchronized blinking LED calibration board with a known
size or a blinking checkerboard pattern generated by an LED screen [14, 15,
16]. This allowed the event camera to identify features similar to a
conventional camera and utilize traditional calibration methods. However,
these toolkits require complex device preparation and are unsuitable for
calibrating extrinsic parameters between event- and frame-based cameras.

In recent years, there has been increased focus on designing new calibration
frameworks to facilitate event camera calibration using existing calibration
boards. Muglikar et al. [17] utilize deep learning-based image reconstruction
networks, such as E2VID [18], to record events generated by moving the
calibration board and then apply the reconstructed images to classical
calibration methods. However, these methods heavily rely on the quality of
image reconstruction and face challenges in achieving time synchronization
with conventional cameras. Another approach involves directly utilizing events
generated during camera motion for camera calibration. Huang et al. [19]
proposed a calibration framework based on a circular calibration board and
employed B-splines to optimize the movement trajectory, which is the most
similar method to the one proposed in this letter. However, they directly use
clustered asynchronous events as features for optimization, compromising sub-
pixel accuracy and being highly sensitive to noise. Additionally, their focus
is solely on the intrinsic calibration of event cameras, without addressing
extrinsic calibration between event cameras and frame-based cameras. Salah et
al. [20] also utilize circular calibration boards and introduce eRWLS to fit
circular features with sub-pixel accuracy. However, they compress events over
a period into a fixed timestamp to obtain a reference “frame” making this
method challenging to synchronize with a conventional camera. Furthermore, it
does not account for the deformation of circular features at different viewing
angles, leading to reduced sub-pixel localization accuracy. The calibration
framework proposed in this letter continues this concept and provides an
improvement to address the problems of these methods.

##  III Preliminaries

###  III-A Event-Based Vision

Unlike conventional cameras, each pixel of the event camera is independently
triggered and responds to changes in the logarithmic illumination signal
$L(\mathbf{u}_{k},t_{k})$. An event $(\mathbf{u}_{k},t_{k},p_{k})$ is
triggered when the change in logarithmic illumination received by a pixel
$\mathbf{u}_{k}=(x_{k},y_{k})$ exceeds a threshold value $C$, i.e.

$\Delta L(\mathbf{u}_{k},t_{k})\doteq
L(\mathbf{u}_{k},t_{k})-L(\mathbf{u}_{k},t_{k}-\Delta t)=p_{k}C$ (1)

where $\Delta t$ is the time since the last triggered event by the same pixel,
$p_{k}\in\\{-\text{1},+\text{1}\\}$ is the polarity of the event.

###  III-B Continuous-Time Trajectory Representation

Continuous-time trajectories are often represented utilizing a weighted
combination of the temporal basis functions [21], such as polynomial
functions, FFTs, and Bézier curves. In this letter, the uniform B-spline is
introduced as a representation of the continuous-time trajectory. B-splines
have the advantages of smoothness, local support, and analytic derivatives,
which are well-suited for representing the 6-DoF pose of the event camera [23,
22, 24]. Following the formulation of cumulative $k$th degree B-spline
$\mathcal{L}$, the event camera pose
$\mathbf{T}^{w}_{e}(\tau)\in\mathbb{SE}\text{3}$ at any time
$\tau\in[t_{i},t_{i+1})$, it can be represented by $N$ control points
$\mathbf{T}_{i}\in\mathbb{SE}\text{3},i\in[\text{0},\text{1},\ldots,N-\text{1}]$:

$\mathcal{L}:\mathbf{T}_{e}^{w}(\tau)=\mathbf{T}_{i}\cdot\prod_{j=\text{1}}^{k}\mathrm{Exp}\left(\tilde{\mathbf{B}}_{j}(\tau)\cdot\mathrm{Log}\left(\mathbf{T}_{i+j-1}^{-1}\mathbf{T}_{i+j}\right)\right)$
(2)

where $\tilde{\mathbf{B}}_{j}(\tau)$ is the cumulative basis function, which
is denoted by

$\tilde{\mathbf{B}}_{j}(\tau)=\tilde{\mathbf{M}}^{(k)}\mathbf{u}$ (3)

$\mathbf{u}=\begin{bmatrix}1&u&\cdots&u^{k}\end{bmatrix}^{T},u=\frac{\tau-
t_{i}}{t_{i+1}-t_{i}}$ (4)

where $\tilde{\mathbf{M}}^{(k)}\in\mathbb{R}^{(k+\text{1})\times(k+\text{1})}$
is the cumulative blending matrix of B-splines. Since the control points of
the B-splines are uniformly distributed on the time scale, the cumulative
blending matrix $\mathbf{M}^{(k)}$ is constant. In this letter, considering
the continuity and complexity, we use cubic B-splines to represent the camera
pose, i.e., $k=\text{3}$. The corresponding cumulative mixing matrix
$\tilde{\mathbf{M}}^{(3)}$ [25] is

$\tilde{\mathbf{M}}^{(\text{3})}=\dfrac{\text{1}}{\text{6}}\begin{bmatrix}\text{6}&\text{0}&\text{0}&\text{0}\\\
\text{5}&\text{3}&\text{-3}&\text{1}\\\
\text{1}&\text{3}&\text{3}&\text{-2}\\\
\text{0}&\text{0}&\text{0}&\text{1}\end{bmatrix}$ (5)

Commonly, to simplify the computation, some work decouples the rotation
$\mathbf{R}_{e}^{w}(\tau)\in\mathbb{SO}$3 and translation
$\mathbf{p}_{e}^{w}(\tau)\in\mathbb{R}^{\text{3}}$ into two independent cubic
B-splines, and the same process is carried out in this paper as well. Hence,
the continuous-time trajectory of the camera pose can be finally formulated as

$\mathbf{R}_{e}^{w}(\tau)=\mathbf{R}_{i}\cdot\prod_{j=\text{1}}^{3}\mathrm{Exp}\left(\tilde{\mathbf{B}}_{j}(\tau)\cdot\mathrm{Log}\left(\mathbf{R}_{i+j-1}^{-1}\mathbf{R}_{i+j}\right)\right)$
(6)

$\mathbf{p}_{e}^{w}(\tau)=\mathbf{p}_{i}+\sum_{j=\text{1}}^{3}\tilde{\mathbf{B}}_{j}(\tau)\cdot(\mathbf{p}_{i+j}-\mathbf{p}_{i+j-1})$
(7)

After decoupling the pose into two cubic B-splines, the corresponding analytic
derivatives [26] can also be derived

$\displaystyle\dot{\mathbf{R}_{e}^{w}}(\tau)$
$\displaystyle=\mathbf{R}_{e}^{w}(\tau)\cdot\left(\boldsymbol{\omega}^{(\text{3})}(\tau)\right)_{\wedge}$
(8)
$\displaystyle=\mathbf{R}_{i}\left(\dot{\mathbf{A}}_{\text{1}}\mathbf{A}_{\text{2}}\mathbf{A}_{\text{3}}+\mathbf{A}_{\text{1}}\dot{\mathbf{A}}_{\text{2}}\mathbf{A}_{\text{3}}+\mathbf{A}_{\text{1}}\mathbf{A}_{\text{2}}\dot{\mathbf{A}}_{\text{3}}\right)$

$\mathbf{v}_{e}^{w}(\tau)=\dot{\mathbf{p}}_{e}^{w}(\tau)=\mathbf{p}_{i}\cdot\sum_{j=\text{1}}^{\text{3}}\dot{\tilde{\mathbf{B}_{j}}}(\tau)\cdot(\mathbf{p}_{i+j}-\mathbf{p}_{i+j-1})$
(9)

where

$\mathbf{A}_{j}=\mathrm{Exp}\left(\tilde{\mathbf{B}}_{j}(\tau)\cdot\mathrm{Log}\left(\mathbf{R}_{i+j-1}^{-1}\mathbf{R}_{i+j}\right)\right)$
(10)

$\dot{\mathbf{A}}_{j}=\mathbf{A}_{j}\dot{\tilde{\mathbf{B}}}(\tau)_{j}\mathrm{Log}\left(\mathbf{R}_{i+j-1}^{-1}\mathbf{R}_{i+j}\right)$
(11)

$\dot{\tilde{\mathbf{B}_{j}}}(\tau)=\dfrac{1}{\Delta
t}\tilde{\mathbf{M}}^{(\text{3})}\begin{bmatrix}\text{0}\\\\[3.00003pt]
\text{1}\\\\[3.00003pt] \text{2}u\\\\[3.00003pt]
\text{3}u^{\text{2}}\end{bmatrix}$ (12)

Since this letter utilizes the uniform B-spline, $\Delta t$ is equal to the
time interval between any two consecutive knots, i.e. $\Delta
t=t_{i+1}-t_{i}$.

Figure 2: Flowchart of the proposed calibration framework

##  IV Methodology

###  IV-A Calibration Framework

For camera calibration, it is of utmost importance to accurately and robustly
identify features on the calibration pattern. However, the checkerboard
pattern [27], which is widely used, is difficult to apply to event cameras
because events disappear during parallel edge motion. Consequently,
calibration patterns for event cameras often employ circular features, but
these can produce blur in moving frames. To balance the characteristics of
event cameras and frame-based cameras, we designed a new calibration pattern
that combines isotropic circles and checkerboard crosspoints, as shown in Fig.
1(a). The center of each circle in this pattern coincides with the center of
the inner crosspoint. This hybrid pattern significantly enhances the
recognition efficiency and accuracy of the event camera while maintaining
compatibility with frame-based cameras.

Fig. 2 presents the flowchart of the proposed calibration framework. In this
letter, we divide the entire calibration process into two stages. The first
stage focuses on feature extraction and refinement of the calibration pattern.
The second stage is dedicated to optimizing the camera trajectory using
piecewise B-splines to achieve accurate calibration results. These two stages
will be elaborated in the following subsections.

###  IV-B Event-Based Feature Recognizer

Unlike frame-based camera, event camera only output asynchronous events during
relative motion, posing a challenge for robust calibration pattern
recognition. To address this, we propose an event-based calibration pattern
feature recognizer, illustrated in Fig. 3. First, we accumulate events over a
short period of time $\Delta t$ according to their polarity to obtain
“accumulation frames” that resemble traditional images. The introduction of
this “accumulation frames” can help us to recognize the feature plate using
some classical image processing algorithms. The following subsections describe
the recognition algorithm based on “accumulation frames” in detail.

Figure 3: The pipeline of event-based feature recognizer.

####  IV-B1 Noise Suppression

Event cameras often generate a large number of events from the static
background when they are in motion. These noisy events can adversely affect
the recognition of calibration boards, significantly reducing the operation
speed and accuracy of the recognizer. Therefore, after obtaining the
“accumulation frames”, a noise suppression module is designed to filter out
most of the events that are not related to the calibration board.

For circular features, the triggered events typically consist of two
semicircular arcs connecting regions of opposite polarity. However, many
structures in the background have straight edges, making them more likely to
have connected regions that resemble straight lines. To leverage this
property, we introduce a fast and accurate connected component labeling (CCL)
algorithm called BBDT, proposed by Grana et al. [28]. This algorithm merges
neighboring events with the same polarity to obtain all the connectivity
regions. Then, the magnitudes of the two principal components of each
connected region are calculated using PCA. For background-triggered connected
regions, the magnitude of the second principal component
$\|\textbf{PC}_{\text{2}}\|$ should be much smaller than the magnitude of the
first principal component $\|\textbf{PC}_{\text{1}}\|$, so that a large number
of noisy regions can be suppressed by the principal component magnitude ratio
$\beta_{PC}$, as given that
$\beta_{PC}=\|\textbf{PC}_{\text{1}}\|/\|\textbf{PC}_{\text{2}}\|<T_{PC}$,
where $T_{PC}$ is a threshold for the $\beta_{PC}$, and any region with
$\beta_{PC}$ higher than $T_{PC}$ is suppressed and not involved in subsequent
operations.

####  IV-B2 Feature Extraction

Following noise suppression, we proceed to extract potential circular features
from the remaining regions. Specifically, we identify two candidate regions of
opposite polarity based on their distance. Subsequently, we fit the elliptic
equation using all pixels contained within these regions, exploiting the fact
that circular features adhere to the elliptic model under a projective
transform. Then the fitting error $e_{fit}$ is calculated, excluding candidate
regions with a fitting error exceeding the fitting threshold $T_{fit}$.

Furthermore, additional geometric constraints are needed to eliminate the
remaining false positive regions. Specifically, the two candidate regions
constituting the same ellipse should demonstrate similar PCA magnitudes;
regions with a notable discrepancy in PCA magnitude fail to meet this
criterion. Additionally, the contributions of the two candidate regions to the
circumference of the ellipse should be close. In other words, the angular
range $\theta_{r}$ of the two candidate regions with respect to the center of
the ellipse should be close to 180∘.

Regions that successfully meet these geometric constraints are recognized as
accurately representing the elliptical features within the calibration
pattern. Consequently, depending on the distribution of these elliptical
features, their relative positions on the calibration pattern can be decoded
to correspond with the results identified in the frame-based camera.

Figure 4: Schematic of the moving ellipse model. The set of events belonging
to the same elliptical feature can be considered as a three-dimensional
oblique elliptical cylinder.

####  IV-B3 Feature Refinement

Neglecting the timestamps of events during elliptical fitting inevitably
introduces errors, thereby affecting the camera calibration precision.
However, by regarding the timestamps $t$ as a third dimension alongside the
pixel coordinates $[x,y]^{T}$, events can be conceptualized as points within a
three-dimensional space. Each “accumulation frames”, corresponding to a brief
time interval, allows for the assumption of solely translational movement with
speed $[v_{x},v_{y}]^{T}$, parallel to the pixel plane, for each elliptical
feature at any given moment within this interval, as demonstrated in Fig. 4.
Consequently, the moving ellipse model $\mathcal{F}$ is described by the
following representation:

$\mathcal{F}:ax(t)^{2}+\beta x(t)y(t)+\gamma y(t)^{2}+\delta x(t)+\epsilon
y(t)+\zeta=\text{0}$ (13)

In matrix form,

$\begin{cases}\mathcal{F}:\mathbf{P}(t)^{T}\mathbf{Q}\mathbf{P}(t)=\text{0}\\\
\mathbf{P}(t)=\mathbf{P}(t_{0})-\mathbf{V}\cdot(t-t_{0})\\\
\mathbf{Q}=\begin{pmatrix}\alpha&\beta/\text{2}&\delta/\text{2}\\\
\beta/\text{2}&\gamma&\epsilon/\text{2}\\\
\delta/\text{2}&\epsilon/\text{2}&\zeta\end{pmatrix}\end{cases}$ (14)

where $\mathbf{P}(t)=[x(t),y(t),\text{1}]^{T}$,
$\mathbf{V}=[v_{x},v_{y}]^{T}$, $t\in[t_{0},t_{0}+\text{2}\delta_{t}]$, and
$t_{0}$ represents the starting time of current “accumulation frame”. Here the
cost function for model optimization is defined as

$\arg\min_{\mathcal{F},\mathbf{V}}\sum_{i\in\\{e\\}}\left(\|\mathbf{P}_{i}^{T}\mathbf{Q}\mathbf{P}_{i}\|^{1}\right)$
(15)

By substituting the events into the aforementioned model, we utilize the
Levenberg-Marquardt algorithm to solve the model, thereby refining the
elliptical features.

###  IV-C Trajectory Optimization

Figure 5: Schematic of a piece-wise B-spline trajectory. The number of event
features inside the red box is insufficient; therefore, this segment of the
trajectory is omitted.

The refinement process converts the previously discrete elliptical features
into densely populated patterns within the corresponding time period of
associated events. This densely populated feature facilitates the optimization
of event camera poses within a continuous-time trajectory representation.
Nevertheless, maintaining continuous visibility of the entire calibration
plate during the calibration process poses a challenge. Incomplete calibration
patterns in a continuous event stream often result in unsuccessful
recognition. To mitigate this issue, the continuous-time trajectory is divided
into multiple segments based on the output from the recognizer. Consequently,
a piece-wise B-spline-based optimizer for event camera pose trajectories is
proposed.

First, based on the predefined knot interval $\Delta t$ and the results of the
recognizer, the features whose timestamps differ from the timestamps of other
features by more than $\Delta t$ are eliminated. In addition, segments
containing too few features are also excluded to ensure optimization accuracy,
and only the more desirable segments are preserved. Fig. 5 illustrates the
segmentation process of the trajectory. The entire calibration process
$\Lambda_{\mathcal{L}}$ is divided into a combination of $M$ segments of
trajectories $\mathcal{L}$:

$\Lambda_{\mathcal{L}}=\sum_{m}\\{\mathcal{L}_{m}\\}_{a_{m}}^{b_{m}},m=\text{1},\ldots,M$
(16)

where $a_{m}$ and $b_{m}$ are the starting and ending times corresponding to
the $m$th segment of B-splines, respectively. Each segment of B-splines is
optimized by only the features whose timestamps belong to its time period.

Figure 6: Schematic of sampling on a continuously moving ellipse model.

For $\mathcal{L}_{m}$, the corresponding state vector $\mathcal{X}_{e}$ is:

$\mathcal{X}_{e}=[\xi^{\text{1}}_{m},\ \xi^{\text{2}}_{m},\ \cdots,\
\xi^{N_{m}}_{m},\ K_{e},\ D_{e}]$ (17)

where $\xi_{m}^{i}$ is the control point of $\mathcal{L}_{m}$, $N_{m}$ is the
number of control points, and $K_{e}$ and $D_{e}$ are the intrinsic and
distortion parameters of the event camera. The corresponding visual residual
for the $i$th feature based on reprojection error is defined as:

$\displaystyle\mathbf{r}_{e}(\mathcal{X}_{e})=\sum_{k\in\mathcal{K}}\pi_{e}(\mathbf{R}_{w}^{e}(t_{i}+k\cdot\delta_{t_{i}})P_{i}^{w}+\mathbf{t}_{w}^{e}(t_{i}+k\cdot\delta_{t}))$
(18) $\displaystyle-(\begin{bmatrix}u_{i}^{e}\\\
v_{i}^{e}\end{bmatrix}+k\cdot\delta_{t}\mathbf{V}_{i})$

where $\mathbf{R}^{e}_{w}(\cdot)$ and $\mathbf{t}^{e}_{w}(\cdot)$ are derived
from the B-spline trajectory using equations (6) and (7), respectively. Since
the feature refinement yields a continuous moving ellipse model, the residuals
can be constructed by sampling the model at any time. Here, $\mathcal{K}$
denotes the partition of $\delta_{t_{i}}$, which defines the sampling interval
of the feature, as shown in Fig. 6. The function $\pi_{e}(\cdot)$ projects the
spatial point $P_{i}^{w}$ onto the “accumulation frame”.

The intrinsic and distortion parameters of the event camera, along with the
control points of the splines, are jointly optimized by minimizing the
following cost function:

$\arg\min_{\mathcal{X}_{e}}\left\\{\sum\rho(\|\mathbf{r}_{e}(\mathcal{X}_{e})\|^{2})\right\\}$
(19)

where $\rho(\cdot)$ is the Huber loss function.

###  IV-D Spatialtemporal Calibration

The final step involves jointly optimizing the event camera and the frame-
based camera, utilizing the previously optimized trajectories, to determine
the extrinsic parameters and time offset. The corresponding state vector
$\mathcal{X}_{f}$ is:

$\mathcal{X}_{f}=[K_{f},\ D_{f},\ \mathbf{T}^{f}_{e},\ t_{d}]$ (20)

where $\mathbf{T}^{f}_{e}$ is the transformation matrix between the two
cameras and $t_{d}$ is the difference between the real timestamps of the two
cameras, i.e., the time offset.

Similarly, define the visual residuals based on reprojection errors in
spatiotemporal calibration as:

$\mathbf{r}_{f}(\mathcal{X}_{f})=\pi_{f}\left(\mathbf{R}^{f}_{e}\left(\mathbf{R}^{e}_{w}(t_{i}+t_{d})P^{w}_{i}+\mathbf{t}^{e}_{w}(t_{i}+t_{d})\right)+\mathbf{t}^{f}_{e}\right)-\begin{bmatrix}u^{f}_{i}\\\\[6.0pt]
v^{f}_{i}\end{bmatrix}$ (21)

where $\pi_{e}(\cdot)$ projects the spatial points onto the image plane of the
frame-based camera.

From (21), the Jacobian $J_{t_{d}}$ of $\mathbf{r}_{f}$ w.r.t $t_{d}$ can be
obtained by the chain rule:

$J_{t_{d}}=\frac{\partial\mathbf{r}_{f}}{\partial P_{i}^{f}}\frac{\partial
P_{i}^{f}}{\partial t_{d}}=\frac{\partial\mathbf{r}_{f}}{\partial
P_{i}^{e}}(\frac{\partial
P_{i}^{f}}{\partial\mathbf{R}_{w}^{e}}\frac{\partial\mathbf{R}_{w}^{e}}{\partial
t_{d}}+\frac{\partial
P_{i}^{f}}{\partial\mathbf{t}_{w}^{e}}\frac{\partial\mathbf{t}_{w}^{e}}{\partial
t_{d}})$ (22)

where
$P_{i}^{f}=\mathbf{R}^{f}_{e}\left(\mathbf{R}^{e}_{w}(t_{i}+t_{d})P^{w}_{i}+\mathbf{t}^{e}_{w}(t_{i}+t_{d})\right)+\mathbf{t}^{f}_{e}$.

Based on (8), (9), and (22), the structure of $\partial P_{i}^{f}/\partial
t_{d}$ can be derived straightforwardly:

$\frac{\partial
P_{i}^{f}}{\partial\mathbf{R}_{w}^{e}}\frac{\partial\mathbf{R}_{w}^{e}}{\partial
t_{d}}=\mathbf{R}^{f}_{e}\dot{\mathbf{R}}^{e}_{w}(t_{i}+t_{d})P^{w}_{i}+\mathbf{R}^{f}_{e}\mathbf{v}^{e}_{w}(t_{i}+t_{d})$
(23)

To jointly optimize the intrinsic and extrinsic parameters of the frame-based
camera, as well as the time offset, the following cost function is minimized:

$\arg\min_{\mathcal{X}_{f}}\left\\{\sum\rho(\|\mathbf{r}_{f}(\mathcal{X}_{f})\|^{2})\right\\}$
(24)

Notably, the optimization problems in (15), (19), and (24) are solved using
Google Ceres111https://ceres-solver.org/.

##  V Experiments

In this section, several experiments are conducted to evaluate the performance
of EF-Calib, encompassing intrinsic calibration test, extrinsic calibration
test, and time offset calibration test. Additionally, an ablation study is
conducted to evaluate the contribution of several key modules within EF-Calib.

###  V-A System Setup

A real-world stereo vision system was designed as Fig. 1(b) shows. It contains
an event camera and a frame-based camera. The two cameras are integrated by a
slide, on which the baseline and viewing angle can be arbitrarily changed to
test the calibration performance of EF-Calib in different situations
comprehensively. The event camera utilized in this letter is the Inivation
DAVIS 346, featuring a resolution of 346$\times$260 and a maximum temporal
resolution of 1 µs. Additionally, this type of camera can also generate
regular frames at a frequency of 30 Hz under standard illumination conditions.
This configuration can be readily employed to compare EF-Calib with a high-
quality, frame-based calibration pipeline, such as the OpenCV calibration
toolbox [29]. The frame-based camera employed is the HikVision MV-CE013-80UM
industrial camera with a global shutter and a resolution of 1280$\times$1024
pixels. Note that no hardware synchronization was utilized in the stereo
vision system. This deliberate choice was made to provide a more rigorous
evaluation of the calibration capability of EF-Calib under real-world
conditions.

###  V-B Calibration Experiments

The calibration performance of a stereo vision system is usually affected by
the camera baseline and viewing angle. To fully evaluate the calibration
performance of EF-Calib, we conducted calibration experiments in three
settings and analyzed the corresponding calibration results separately. In the
first setting (Trial 1), the cameras are configured for a regular baseline. In
the second setting (Trial 2), the cameras are configured for a wide baseline.
In the third setup (Trial 3), the cameras are configured as a narrow baseline.
For each trial, the cameras are adjusted to obtain a reasonable viewing angle,
ensuring sufficient overlap of the camera field of view. For each trial,
images with event data are recorded simultaneously for sufficient time (about
40 s) to achieve converged calibration results.

####  V-B1 Intrinsic Calibration Test

We utilized only the event stream data from the above three trials and
completed the intrinsic calibration of the event camera respectively.
Previously, we completed the intrinsic calibration using OpenCV toolbox [29]
with the frame provided by DAVIS 346 and considered this calibration result as
the ground truth. In addition, EF-Calib is compared with two state-of-the-art
event camera intrinsic calibration methods [17, 20]. Note that the calibration
patterns used by the compared methods are the ones originally employed by
them: [17] utilizes a checkerboard pattern, while [20] employs an asymmetric
circular pattern.

Table I shows the intrinsic calibration results of each method. It can be seen
that the intrinsic parameter obtained by our method is closest to the ground
truth (GT) and the results corresponding to three trials are very stable. In
addition, Fig. 7 illustrates the plot of the intrinsic parameters over time.
It can be noticed that EF-Calib can get converged results in less than 20 s,
demonstrating the ease of use of our method.

Figure 7: Results of the intrinsic calibration test. TABLE I: Comparative Results of Intrinsic Calibration Test Methods | $f_{x}$ | $f_{y}$ | $c_{x}$ | $c_{y}$ | $k_{\text{1}}$ | $k_{\text{2}}$ | RPE
---|---|---|---|---|---|---|---
Frame-based (GT) | 413.84 | 413.80 | 157.42 | 132.25 | -0.38 | 0.31 | 0.13
E2Calib [17] | 417.56 | 417.24 | 159.86 | 132.35 | -0.36 | 0.09 | 0.41
E-Calib [20] | 404.66 | 403.99 | 159.81 | 132.69 | -0.37 | 0.32 | 0.33
EF-Calib (Trial 1) | 414.06 | 413.28 | 158.03 | 132.43 | -0.38 | 0.31 | 0.10
EF-Calib (Trial 2) | 414.79 | 413.85 | 158.74 | 131.90 | -0.38 | 0.34 | 0.12
EF-Calib (Trial 3) | 414.87 | 413.94 | 159.16 | 133.65 | -0.37 | 0.31 | 0.13

####  V-B2 Extrinsic Calibration Test

To evaluate the extrinsic calibration performance of the EF-Calib, we
calculated the errors corresponding to rotation and translation separately for
each trial, i.e.

$\displaystyle e_{t}$
$\displaystyle=\frac{1}{N}\sum_{i}^{N}\left\|\mathbf{t}^{e}_{w}(t_{i}+t_{d})-\mathbf{T}_{f}^{e}{\mathbf{t}^{f}_{w}}_{i}\right\|_{2}$
(25) $\displaystyle e_{r}$
$\displaystyle=\frac{1}{N}\sum_{i}^{N}\|\boldsymbol{\theta}(\mathbf{R}_{w}^{e}(t_{i}+t_{d}))-\boldsymbol{\theta}(\mathbf{R}_{f}^{e}{\mathbf{R}_{w}^{f}}_{i})\|_{2}$

where $N$ is the frame number and $\boldsymbol{\theta}(\cdot)$ represents the
Euler angle corresponding to the rotation matrix. Similarly to the intrinsic
calibration test, we also acquired 30 pairs of images containing a
checkerboard calibration board from the two cameras in different poses
simultaneously. These frames were calibrated using the OpenCV toolbox to
obtain the corresponding extrinsic parameters for both cameras. The computed
EF-Calib calibration errors were compared with the corresponding errors of the
frame-based extrinsic parameter calibration. From Table II, it can be seen
that EF-Calib can achieve the same level of error as the frame-based extrinsic
calibration, verifying its effectiveness in extrinsic calibration.

TABLE II: Comparative Results of Extrinsic Calibration Test Trial | Method |  $e_{t}$ (mm) |  $e_{r}$ (∘) | Frames
---|---|---|---|---
Trial 1 | Frame-based | 0.3499 | 0.0918 | 30
EF-Calib | 0.5336 | 0.1984 | 250
Trial 2 | Frame-based | 0.4199 | 0.1854 | 30
EF-Calib | 0.6572 | 0.3127 | 207
Trial 3 | Frame-based | 0.2541 | 0.0954 | 30
EF-Calib | 0.3638 | 0.2911 | 184

####  V-B3 Time Offset Calibration Test

Time offset estimation is crucial for multi-camera systems that lack hardware
triggering. In this test, we evaluate EF-Calib’s capability to calibrate time
offsets by manually adjusting the timestamp of each image frame. Specifically,
we compare the differences between the time offset calibration results
obtained with modified timestamps and the original results. This is achieved
by uniformly delaying or advancing the timestamps by a fixed time interval
$\Delta t_{d}$. Specifically, the timestamps were modified by $\pm$2.5 ms and
$\pm$5 ms for each trial, and the difference between the modified time delay
and the original time delay was calculated and compared to the delta. Fig. 8
illustrates the experimental results of the time offset calibration test, and
EF-Calib can accurately calibrate the time offset at different scenarios.

Figure 8: Results of the time offset calibration test. TABLE III: Ablation Study of EF-Calib Piece-wise trajectory |  Temporal calibration |  Feature refinement | Trial 1 | Trial 2 | Trial 3
---|---|---|---|---|---
$e_{t}$ | $e_{r}$ | $e_{t}$ | $e_{r}$ | $e_{t}$ | $e_{r}$
|  |  | 8.1633 | 2.8632 | 4.6642 | 2.6372 | 1.2161 | 1.4511
✔ |  |  | 8.0019 | 2.8187 | 3.8794 | 2.7859 | 0.4974 | 1.4322
✔ | ✔ |  | 0.9141 | 5.6910 | 0.6322 | 0.3239 | 0.3707 | 0.3065
✔ | ✔ | ✔ | 0.5336 | 0.1984 | 0.6572 | 0.3127 | 0.3638 | 0.2911

###  V-C Ablation Study

To thoroughly analyze and validate the performance and functionality of each
module within EF-Calib, we conducted an ablation study. Specifically, we
scrutinized three modules: piece-wise trajectories, temporal calibration, and
feature refinement. Table III demonstrates the impact of the introduction of
these three modules on the calibration error of the EF-Calib extrinsic
parameters. As can be seen from Table III, the introduction of all three
modules can greatly improve the extrinsic parameter calibration accuracy.

##  VI Conclusion

In this letter, we propose a spatiotemporal calibration framework called EF-
Calib, aiming to achieve joint calibration of intrinsic parameters, extrinsic
parameters, and time offset for event and frame-based cameras. First, we
design a novel calibration pattern that accommodates the heterogeneous nature
of event and frame-based representations, incorporating both circles and
crosspoints to facilitate simultaneous recognition by both camera types. A
corresponding event-based recognition algorithm is developed to ensure robust
and accurate feature recovery using this pattern. Additionally, to manage the
asynchronous characteristics of the event stream, we introduce a piece-wise
B-spline to continuously represent the pose trajectory of the event camera.
Finally, we provide the analytic Jacobian of the error term and implement the
joint calibration of intrinsic, extrinsic, and time offset for both camera
types. Experimental results demonstrate that EF-Calib outperforms current
state-of-the-art methods in intrinsic parameter estimation while also
achieving high accuracy in extrinsic parameter and time offset estimation.
These results demonstrate the spatiotemporal calibration capabilities of EF-
Calib and lay a robust foundation for the fusion of event and frame.

In the future, we aim to explore markerless online calibration based on EF-
Calib. Additionally, we plan to utilize EF-Calib to create novel visual
perception frameworks that fuse events and frames.

## References

  * [1] G. Gallego, T. Delbruck, G. Orchard, C. Bartolozzi, B. Taba, A. Censi, S. Leutenegger, A. Davison, J. Conradt, and K. Daniilidis, “Event-based vision: A survey,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 44, no. 1, pp. 154–180, Jan. 2020.
  * [2] S. Zhu, Z. Tang, M. Yang, E. Learned-Miller, and D. Kim, “Event camera-based visual odometry for dynamic motion tracking of a legged robot using adaptive time surface,” in _Proc. IEEE/RSJ Int. Conf. Intell. Rob. Syst._ , Detroit, MI, USA, 2023, pp. 3475–3482.
  * [3] Z. Zhou, Z. Wu, R. Boutteau, F. Yang, C. Demonceaux, and D. Ginhac, “RGB-Event fusion for moving object detection in autonomous driving,” in _Proc. IEEE Int. Conf. Robot. Autom._ , London, United Kingdom, Jul. 2023, pp. 7808–7815.
  * [4] J. Jiang, J. Li, B. Zhang, X. Deng, and B. Shi, “EvHandPose: Event-based 3D hand pose estimation with sparse supervision,” _IEEE Trans. Pattern Anal. Mach. Intell._ , early access, doi: 10.1109/TPAMI.2024.3380648, 2024.
  * [5] J. Han, Y. Yang, P. Duan, C. Zhou, L. Ma, C. Xu, T. Huang, I. Sato, B. Shi, “Hybrid high dynamic range imaging fusing neuromorphic and conventional images,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 45, no. 7, pp. 8553–8565, Jul. 2023.
  * [6] A. R. Vidal, H. Rebecq, T. Horstschaefer, and D. Scaramuzza, “Ultimate SLAM? Combining events, images, and IMU for robust visual SLAM in HDR and high-speed scenarios,” _IEEE Robot. Automat. Lett._ , vol. 3, no. 2, pp. 994–1001, Apr. 2018.
  * [7] M. S. Lee, J. H. Jung, Y. J. Kim, and C. G. Park, “Event-and frame-based visual-inertial odometry with adaptive filtering based on 8-DOF warping uncertainty,” _IEEE Robot. Automat. Lett._ , vol. 9, no. 2, pp. 1003–1010, Feb. 2024,
  * [8] W. Guan, P. Chen, Y. Xie, and P. Lu, “PL-EVIO: Robust monocular event-based visual inertial odometry with point and line features,” _IEEE Trans. Automat. Sci. Eng._ , early access, doi: 10.1109/TASE.2023.3324365, 2023.
  * [9] C. Luo, J. Wu, S. Sun, and P. Ren, “TransCODNet: Underwater transparently camouflaged object detection via RGB and event frames collaboration,” _IEEE Robot. Automat. Lett._ , vol. 9, no. 2, pp. 1444–1451, Feb. 2024.
  * [10] P. Chen, W. Guan, F. Huang, Y. Zhong, W. Wen, L. Hsu, P. Lu, “ECMD: An event-centric multisensory driving dataset for SLAM,” _IEEE Trans. Intell. Veh._ , vol. 9, no. 1, pp. 407–416, Jan. 2024.
  * [11] C. Creß, W. Zimmer, N. Purschke, B. N. Doan, S. Kirchner, V. Lakshminarasimhan, L. Strand, and A. Knoll, “TUMTraf event: Calibration and fusion resulting in a dataset for roadside event-based and RGB cameras,” _IEEE Trans. Intell. Veh._ , early access, doi: 10.1109/TIV.2024.3393749, 2024.
  * [12] L. Gao, Y. Liang, J. Yang, S. Wu, C. Wang, J. Chen, and L. Kneip, “VECtor: A versatile event-centric benchmark for multi-sensor SLAM,” _IEEE Robot. Automat. Lett._ , vol. 7, no. 3, pp. 8217–8224, Jul. 2022.
  * [13] Z. Zhang, “A flexible new technique for camera calibration,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 22, no. 11, pp. 1330–1334, Nov. 2000.
  * [14] “Calibration toolbox by RPG, University of Zurich,” https://github.com/uzh-rpg/rpg_dvs_ros/tree/master/dvs_calibration.
  * [15] G. Orchard, “Calibration toolbox by G. Orchard,” https://github.com/gorchard/DVScalibration.
  * [16] “Calibration toolbox by VLOGroup at TU Graz,” https://github.com/VLOGroup/dvs-calibration.
  * [17] M. Muglikar, M. Gehrig, D. Gehrig, and D. Scaramuzza, “How to calibrate your event camera,” in _Proc. IEEE/CVF Conf. Comput. Vis. Pattern Recognit. Workshops_ , Nashville, TN, USA, Jun. 2021, pp. 1403–1409.
  * [18] H. Rebecq, R. Ranftl, V. Koltun, and D. Scaramuzza, “High speed and high dynamic range video with an event camera,” _IEEE Trans. Pattern Anal. Mach. Intell._ vol. 43, no. 6, pp. 1964–1980, Jun. 2021.
  * [19] K. Huang, Y. Wang, and L. Kneip, “Dynamic event camera calibration,” in _Proc. IEEE/RSJ Int. Conf. Intell. Rob. Syst._ , Prague, Czech Republic, Sep. 2021, pp. 7021–7028.
  * [20] M. Salah, A. Abdulla, H. Muhammad, G. Daniel, A. Abdelqader, S. Lakmal, S. Davide, and Z. Yahya, “E-Calib: A fast, robust and accurate calibration toolbox for event cameras,” Jun. 2023, arXiv:2306.09078.
  * [21] J. Rehder, R. Siegwart, and P. Furgale, “A general approach to spatiotemporal calibration in multisensor systems,” _IEEE Trans. Robot._ , vol. 32, no. 2, pp. 383–398, Apr. 2016.
  * [22] J. Huai, Y. Zhuang, Y. Lin, G. Jozkow, Q. Yuan, and D. Chen, “Continuous-time spatiotemporal calibration of a rolling shutter camera-IMU system,” _IEEE Sensors J._ , vol. 22, no. 8, pp. 7920–7930, Apr. 2022.
  * [23] E. Mueggler, G. Gallego, H. Rebecq, and D. Scaramuzza, “Continuous-time visual-inertial odometry for event cameras,” _IEEE Trans. Robot._ , vol. 34, no. 6, pp. 1425–1440, Dec. 2018.
  * [24] A. Patron-Perez, S. Lovegrove, and G. Sibley, “A spline-based trajectory representation for sensor fusion and rolling shutter cameras,” _Int. J. Comput. Vis._ , vol. 113, no. 3, pp. 208–219, Feb. 2015.
  * [25] K. Qin, “General matrix representations for B-splines,” _Visual Comput._ , vol. 16, no. 3, pp. 177-186, 2000.
  * [26] C. Sommer, V. Usenko, D. Schubert, N. Demmel, and D. Cremers, “Efficient derivative computation for cumulative B-splines on lie groups,” in _Proc. IEEE/CVF Conf. Comput. Vis. Pattern Recognit._ , Seattle, WA, USA, Jun. 2020, pp. 11145–11153
  * [27] S. Wang, M. Zhu, Y. Hu, D. Li, F. Yuan, and J. Yu, “Accurate detection and localization of curved checkerboard-like marker based on quadratic form,” _IEEE Trans. Instrum. Meas._ , vol. 71, pp. 1–11, Jul. 2022.
  * [28] C. Grana, D. Borghesani, and R. Cucchiara, “Optimized block-based connected components labeling with decision trees,” _IEEE Trans. Imag. Process._ , vol. 19, no. 6, pp. 1596–1609, Jun. 2010.
  * [29] G. Bradski and A. Kaehler, _Learning OpenCV: Computer vision with the OpenCV library_. Sebastopol, CA, USA: O’Reilly, 2008, pp. 370–396.