# Instant Continual Learning of Neural Radiance Fields

Ryan Po Zhengyang Dong Alexander W. Bergman Gordon Wetzstein
Stanford University
{rlpo, awb, leozdong<EMAIL_ADDRESS>
###### Abstract

Neural radiance fields (NeRFs) have emerged as an effective method for novel-
view synthesis and 3D scene reconstruction. However, conventional training
methods require access to all training views during scene optimization. This
assumption may be prohibitive in continual learning scenarios, where new data
is acquired in a sequential manner and a continuous update of the NeRF is
desired, as in automotive or remote sensing applications. When naively trained
in such a continual setting, traditional scene representation frameworks
suffer from catastrophic forgetting, where previously learned knowledge is
corrupted after training on new data. Prior works in alleviating forgetting
with NeRFs suffer from low reconstruction quality and high latency, making
them impractical for real-world application. We propose a continual learning
framework for training NeRFs that leverages replay-based methods combined with
a hybrid explicit–implicit scene representation. Our method outperforms
previous methods in reconstruction quality when trained in a continual
setting, while having the additional benefit of being an order of magnitude
faster.

##  1 Introduction

High-quality reconstruction and image-based rendering of 3D scenes is a long-
standing research problem spanning the fields of computer vision [23, 36],
computer graphics [7, 18], and robotics [3, 15, 41]. Recently, the
introduction of Neural Radiance Fields (NeRFs) [39] has led to substantial
improvements in this area through the use of differentiable 3D scene
representations supervised with posed 2D images. However, NeRFs require access
to all available views of the 3D scene during training, a condition that is
prohibitive for automotive and remote sensing applications, among others,
where data is sequentially acquired and an updated 3D scene representation
should be immediately available. In such conditions, the scene representation
must be trained in a continual setting, where the model is given access to a
limited number of views at each stage of training, while still tasked with
reconstructing the entire scene.

Figure 1: Continual learning of NeRFs. Conventionally, NeRFs are trained with
access to all training views. However, for continual learning scenarios we
must train on batches of input views without access to previously seen data
(top). When trained in these settings, conventional methods suffer from
catastrophic forgetting, leading to poor reconstructions (center). In
contrast, our method reconstructs the entire scene with high quality (bottom).
Figure 2: Problem overview. (a) Continual learning setting for training NeRFs.
Instead of training the scene representation over all input views at once, the
model is given 2D views of the scene in sequential batches. During a
particular stage in training, the model is only given access to the most
recently captured views. (b) Training NeRF is the continual setting leads to
catastrophic forgetting. Previously learned 3D scene content is corrupted
after training on newly captured views.

When trained in a continual setting, NeRFs suffer from catastrophic forgetting
[17], where previously learned knowledge is forgotten when trained on new
incoming data. Recent work [53, 13] has shown promise in tackling catastrophic
forgetting through replay-based methods. Such approaches aim to alleviate
forgetting by storing information from previous tasks either explicitly or in
a compressed representation, then revisiting this information during training
of subsequent tasks. Existing methods [69, 48] have seen success through the
application of replay-based techniques in conjunction with NeRF for addressing
the task of simultaneous mapping and localization (SLAM) [3], however such
methods either suffer from memory scalability or latency issues.

In this work, we tackle the task of continually learning NeRFs by leveraging
the benefits of replay-based techniques. Specifically, we acknowledge that a
trained NeRF itself is a compressed representation of all previously observed
2D views. By freezing a copy of the scene representation after the training of
each task, we essentially have access to pseudo ground truth RGB values for
all previously seen data by querying this oracle. We also modify the
underlying neural scene representation architecture motivated by one key
insight: catastrophic forgetting is a fundamental problem faced by neural
networks. Therefore, the fully implicit (MLP) representation used by NeRF is
fundamentally ill-suited for the task of continual learning. We minimize the
reliance of our underlying scene model on the decoder neural network by using
a hybrid implicit–explicit representation. By replacing the frequency encoding
in NeRF with a multi-resolution hash encoding [40], we greatly reduce the size
of the decoder multilayer perceptron (MLP), minimizing the effects of
catastrophic forgetting.

As an additional benefit, our method is also an order of magnitude faster than
previous replay-based methods [13]. This enables fast continual scene fitting,
as our method can learn additional 3D scene content from new input views in as
little as 5 seconds (see Section 5.4 for details).

##  2 Related Work

#### Neural radiance fields.

Scene representation networks [51] and neural rendering [58, 59] have emerged
as a family of techniques enabling effective 3D scene reconstruction. Given a
set of images and corresponding ground truth camera poses, neural radiance
fields (NeRFs) [39], for example, optimizes a underlying scene representation
by casting rays, sampling the scene volume and aggregating sampled color and
density values to synthesize an image. The success of NeRFs has spawned a line
of works on improving the quality and efficiency of the method [5, 4, 11, 20,
24, 6, 32, 33, 38, 40, 45, 55, 56, 61, 62, 65, 66, 68], while extending the
method to a range of applications [63, 12, 34, 43, 19, 22, 42, 53, 69]. NeRFs
leverage a neural implicit representation (NIR) [50] in the form of a simple,
yet effective multi-layer pereceptron (MLP) to represent the 3D scene. Many
follow-up works improve on the underlying NIR, enabling features such as real-
time rendering [45, 67, 8] and faster training [40, 33, 66, 10]. A key
limitation for the training of NeRFs is the assumption that all input images
of the target scene are available during training. In scenarios such as
autonomous vehicle or drone footage captures, this assumption no longer holds
as data is sequentially acquired and an updated 3D representation should be
immediately available. NeRFs trained on sequential data suffer from
catastrophic forgetting [46]. Our method overcomes this limitation, providing
a high quality reconstruction of the entire scene, while imparting minimal
computational and memory overhead.

Figure 3: Memory replay through NeRF distillation. Scene representation is
sequentially trained on sequentially acquired views. After each stage of
training, a frozen copy of the scene parameters is stored. While optimizing
for the next set of incoming images, the frozen network is queried to obtain
pseudo ground truth values. The current network $(\Phi_{i},\Theta_{i})$ is
trained on a mixed objective that minimizes photometric loss with respect to
ground truth images from the current task, and pseudo ground truth values for
previous tasks (Equation 4).

#### Continual learning.

Continual learning is a long-standing problem in the field of machine
learning, where partial training data is available at each stage of training.
As mentioned above, NeRFs trained in a continual learning setting suffers from
catastrophic forgetting [46]. Existing work in this field fall into three main
categories [29]: parameter regularization [31, 60, 1, 25], parameter isolation
[2, 64, 37, 16] and data replay [26, 44, 47, 49, 35, 9]. Parameter isolation
methods aim at combating catastrophic forgetting by attempting to learn a sub-
network for each task, while parameter regularization methods identify
parameters important for preserving old knowledge and penalizing changes to
them. Finally, data replay methods preserve previous knowledge by storing a
subset of previous training data. Subsequent tasks are then optimized over old
and new incoming data. Our proposed method leverages a self-distillation
method similar to previous data replay approaches, storing pseudo ground truth
values for all previous training with minimal memory usage.

#### SLAM & continual learning of NeRFs.

Works in the field of simultaneous mapping and localization (SLAM) [3] aim at
reconstructing a 3D scene from a continuous stream of images, similar to the
continual learning setting. Recent works [53, 69, 48] combine NIRs and
traditional SLAM-based methods with promising results. These methods fall
under the data replay category, as they approach the task of continual
learning by explicitly storing key-frames from previous image streams. Storing
data explicitly can be expensive, and designing an appropriate importance
heuristic for selecting key-frames is non-trivial. In contrast, our approach
stores previous data as an implicitly defined generator, greatly reducing
memory overhead.

##  3 Continual Learning of NeRFs

Before we explain the details of our proposed method, it is important to first
formally establish the task of continual learning of NeRFs. We consider the
scenario where $t$ sets of image data come in sequentially, represented by
$\\{\mathcal{I}_{1},\dots,\mathcal{I}_{t}\\}$. Each image data set is
represented by $\mathcal{I}_{i}=(\mathbf{I}_{i},\mathbf{R}_{i})$, where
$\mathbf{I}_{i}$ represents the per-pixel RGB values of the image data and
$\mathbf{R}_{i}$ represent the camera rays corresponding to each image pixel.
Note that $\mathbf{R}_{i}$ can either be explicitly stored as values in
$\mathbb{R}^{6}$ (ray origin and direction) or implicitly through camera
extrinsic and intrinsic matrices.

The objective of our optimization remains the same, we wish to minimize
reconstruction loss across all provided ground truth views in
$\\{\mathcal{I}_{1},\dots,\mathcal{I}_{t}\\}$. However, the training procedure
differs from conventional NeRF training. Training is performed sequentially as
illustrated in Figure 2a. At a given stage of training, our model is only
given access to a subset of all of the RGB images (visualized in Figure 2b),
but access to ray information from all previous tasks. Formally, at time step
$i$, the model is able to access $\mathbf{I}_{i}$ and
$\\{\mathbf{R}_{1},\dots,\mathbf{R}_{i}\\}$. Note that this formulation is
slightly different from prior works such as MEIL-NeRF [13] where access to ray
information is also constrained. However, we believe this constraint is
unwarranted since ray information can be stored implicitly for every input
view with only 6 scalar values111Camera extrinsic matrices can be implicitly
stored in the form $(t_{x},t_{y},t_{z},r_{x},r_{y},r_{z})$, where
$t_{x},t_{y},t_{z}$ represents the position of the camera optical center and
$r_{x},r_{y},r_{z}$ the orientation of the camera., assuming all input views
share the same camera intrinsics. Similar to prior work [13], our method is
based on self distillation [21], therefore we also assume that we have access
to a frozen copy of our trained representation from the previous task.

##  4 Method

In this section, we first provide a brief recap behind the formulation of
NeRFs [39], then introduce our solution to catastrophic forgetting in the
context of training NeRFs in a continual setting. There are two main
contributors to our solution: namely, the use of a hybrid feature
representation (Section 4.2) and task specific network distillation (Section
4.3).

###  4.1 NeRF preliminaries

Neural radiance fields (NeRFs) [39] represent a 3D scene through an implicit
function from a point in 3D space $\mathbf{x}=(x,y,z)$ along with a
corresponding viewing direction $\mathbf{d}=(\theta,\phi)$ to a density value
$\mathbf{\sigma}$ and RGB color $\mathbf{c}=(r,g,b)$. Conventionally, NeRFs
are represented with an MLP characterized by its parameters $\Theta$, giving
the mapping

$F_{\Theta}:(\mathbf{x},\mathbf{d})\mapsto(\mathbf{\sigma},\mathbf{c}).$ (1)

Novel views of the 3D scene are generated through volume rendering [28] of the
5D radiance field. Given an image pixel with the corresponding ray
$\mathbf{r}=(\mathbf{r}_{o},\mathbf{r_{d}})$, by sampling points
$\mathbf{x}_{i}$ along this ray and evaluating the radiance field values
$(\mathbf{\sigma}_{i},\mathbf{c}_{i})$ at these points, the color associated
with this ray can be recovered. With $\mathbf{N}$ sampled points, the RGB
color of a ray $\mathbf{r}$ is obtained by

$\hat{\mathbf{C}}(\mathbf{r};\Theta)=\sum_{i=1}^{N}T_{i}(1-\text{exp}(-\sigma_{i}\delta_{i}))\mathbf{c}_{i},$
(2)

where $\delta_{i}$ represents the distance between the $i^{th}$ and
$(i+1)^{th}$ sampled point and $T_{i}$ represents the accumulated
transmittance from $r_{o}$ to the current sample point, given by
$T_{i}=\text{exp}\left(-\sum_{j=1}^{i-1}\sigma_{j}\delta_{j}\right)$.

###  4.2 Multi-resolution hash encoding

Prior work [40] has found success in replacing the fully implicit $F_{\Theta}$
with a hybrid representation, leading to faster convergence rates along with
better memory and computational efficiency. Hybrid representations map 3D
coordinates to an explicitly defined feature space before passing these
features into a significantly smaller implicit MLP decoder to obtain density
and RGB values. We leverage these explicit feature mappings to alleviate the
effects of catastrophic forgetting.

Table 1: Quantitative results: unconstrained setting. PSNR of different continual learning methods. Every method is trained on each task until convergence, which differs by method. Approximate training time for all 10 tasks is listed next to each method. For each scene, we mark the best performing methods with gold , silver  and bronze  medals. Results marked with * are trained in a non-continual setting, where ground truth data from all tasks are available during scene optimization. These results serve as an upper bound for scenes trained in a continual setting. Our method consistently out-performs all baselines while taking significantly less time to converge. | ScanNet | Tanks & Temples | TUM RGB-D
---|---|---|---
Method |  0101 |  0146 |  0160 |  Truck |  Caterpillar |  Family |  Desk 0 |  Desk 1
NeRF-Incre (2 hours) |  13.70 |  13.20 |  17.31 |  16.88  |  15.36  |  22.96  |  13.05  |  14.03
iNGP-Incre (10 min) |  16.51  |  16.64 |  19.98 |  13.49 |  14.55 |  21.15 |  12.70 |  14.65
iNGP + EWC (10 min) |  16.11 |  17.32  |  20.16  |  12.50 |  13.61 |  19.28 |  12.50 |  10.85
MEIL-NeRF (2 hours) |  24.32  |  26.82  |  28.93  |  22.74  |  20.89  |  26.57  |  20.79  |  19.80
Ours (10 min) |  25.72  |  27.87  |  30.28  |  22.71  |  22.51  |  29.33  |  20.65  |  20.34
NeRF* (2 hours) |  26.15 |  28.48 |  30.88 |  24.80 |  23.14 |  29.33 |  22.35 |  20.88
iNGP* (10 min) |  26.00 |  28.43 |  31.16 |  24.22 |  24.02 |  31.14 |  20.95 |  20.73

#### Multi-resolution feature grids.

Following Instant-NGP [40], we map 3D coordinates to explicit features
arranged into $L$ levels, each level containing a maximum of $T$ features,
with each feature having a dimensionality of $F$. Each level stores features
corresponding to vertices of a 3D grid with fixed resolution. Consider a
single feature level $l$: the queried 3D coordinate $\mathbf{x}$ is first
scaled to match the native resolution of $l$, and the neighboring $2^{3}$
vertices from the fixed resolution 3D grid are identified. Each vertex of
interest is mapped to an entry in the $l^{th}$ level feature array and the
final feature value corresponding to $\mathbf{x}$ is obtained through tri-
linear interpolation. This feature value is then passed into an implicit
function represented by an MLP, mapping from feature space to density and RGB
values.

#### Forgetting in explicit features.

Consider the case where $T$ matches the total number of vertices at each grid
resolution, such that a 1:1 mapping exists between grid vertices and feature
embeddings. In the continual setting, features are only updated when the
corresponding voxel is visible in the training views of the current task,
whereas other features remain constant, unaffected by catastrophic forgetting.
This is in stark contrast to the global updates observed in fully implicit
representations such as in NeRF [39]. In NeRF, each network parameter
influences radiance and density values at every point in 3D space, and
training on new data points overwrites information learned in the entire
scene, even for regions not visible in the current training views.

#### Hashed feature tables.

In an effort to lower memory usage at higher grid resolutions, Instant-NGP
proposes a hashed encoding scheme. At finer levels, a hash function
$h:\mathcal{Z}^{d}\mapsto\mathcal{Z}_{T}$ is used to index into the feature
array, effectively acting as a hash table. Following prior work [40], we use a
spatial hash function of the form

$h(\mathbf{x})=\left(\bigoplus_{i=1}^{d}x_{i}\pi_{i}\right)\;\text{mod}\;T,$
(3)

where $\bigoplus$ represents the bit-wise XOR operator and $\pi_{i}$ are
unique, large primes.

In contrast to dense feature grids, hashed feature tables suffer from
catastrophic forgetting in the feature space due to hash collisions. Consider
a single task $\mathcal{I}_{i}$. A vertex $v_{1}$ visible in $\mathcal{I}_{i}$
may share the same hash table entry as another vertex $v_{2}$ that is not
visible in $\mathcal{I}_{i}$. During training, the training objective will
only optimize the shared hash table entry for the current task
$\mathcal{I}_{i}$, learning the correct feature value for $v_{1}$, while
forgetting any information learnt for $v_{2}$. The effects of forgetting are
dependent on the frequency of hash collisions between grid vertices, which
increases as the hash table size $T$ decreases.

###  4.3 Memory replay through NeRF distillation

Catastrophic forgetting results from a misalignment between the current and
cumulative training objectives. Replay-based approaches [26, 44, 47, 49, 35,
9] combat network forgetting by storing information from previous tasks either
explicitly or implicitly through a generative model. Consider a NeRF with
explicit feature embeddings trained on a set of tasks
$\\{\mathcal{I}_{1},\dots,\mathcal{I}_{i}\\}$, with feature and MLP parameters
characterized by $(\Phi_{i},\Theta_{i})$. We can then treat
$(\Phi_{i},\Theta_{i})$ as a generator for 2D image data found in tasks
$\\{\mathcal{I}_{1},\dots,\mathcal{I}_{i}\\}$. Let $\hat{\mathbf{R}}_{i}$ be
the union of all ground truth rays in the first $i$ tasks. The ground truth
RGB value corresponding to a ray $\mathbf{r}\in\hat{\mathbf{R}}_{i}$ can then
be approximated by $\hat{\mathbf{C}}(\mathbf{r};\Phi_{i},\Theta_{i})$
following Eq. 2.

We approach continual learning in a self-distillation manner [21]. When
training on the subsequent task $\mathcal{I}_{i+1}$, we no longer have access
to ground truth image data from previous tasks. However, as explored in prior
work [13], by saving network parameters $(\Phi_{i},\Theta_{i})$ we effectively
have access to pseudo ground truth values for all rays in
$\hat{\mathbf{R}}_{i}$. We can then modify our training objective to minimize
photometric loss for all rays in tasks
$\\{\mathcal{I}_{1},\dots,\mathcal{I}_{i+1}\\}$, rather than just
$\mathcal{I}_{i+1}$. The modified training objective is then given by

$\displaystyle\mathcal{L}(\Phi,\Theta)_{i+1}$
$\displaystyle=\sum_{\mathbf{r}\in\mathcal{I}_{i+1}}||\hat{\mathbf{C}}(\mathbf{r};\Phi,\Theta)-\mathbf{C}(\mathbf{r})||^{2}$
$\displaystyle+\sum_{\mathbf{r}\notin\mathcal{I}_{i+1}}||\hat{\mathbf{C}}(\mathbf{r};\Phi,\Theta)-\hat{\mathbf{C}}(\mathbf{r};\Phi_{i},\Theta_{i})||^{2}.$
(4)

During each task, we still sample rays uniformly over all previous and current
tasks. However, for previous tasks where ground truth RGB values are no longer
available, we instead query the frozen network to obtain a pseudo ground truth
value. Figure 3 shows a visualization of the replay-based distillation method.

Figure 4: Qualitative results: unconstrained setting. We show reconstructed
views from a previously supervised (forgotten) task across different methods.
Our method consistently outperforms all other baselines in visual quality.
NeRF trained in a continual setting suffers from catastrophic forgetting, as
illustrated by poor early-task reconstruction results. Parameter
regularization through EWC aids in alleviating forgetting effects, however,
reconstruction results still suffer from severe visual artefacts. MEIL-NeRF
adopts a similar replay approach as our method, using a frozen copy of the
scene representation as guidance when training future tasks. However, the
fully implicit representation in MEIL-NeRF forgets high-frequency detail from
earlier tasks. In contrast, our method is able to retain high-frequency
details for earlier tasks through the use of explicit features.

##  5 Experiments

To highlight the effectiveness of our method in overcoming catastrophic
forgetting, we compare our method against existing continual learning methods
[25, 13]. We describe baseline methods in Section 5.1, datasets used in
Section 5.2 and experimental settings in Section 5.3.

###  5.1 Baselines

#### NeRF and iNGP.

We train NeRFs under the continual setting using frequency and multi-
resolution hash encodings, referring to these baselines as NeRF-Incre and
iNGP-Incre respectively. For our hash encoding experiments, we used a feature
grid of $L=16$ levels, a hash table size of $T=2^{17}$, a feature dimension of
$F=2$ and grid resolutions ranging from 16 to 512. We also scale the original
NeRF representation [39] to have 8 fully connected layers with 512 channels
each, matching the total number of trainable parameters as the hash encoding
models.

#### Elastic Weight Consolidation.

Elastic Weight Consolidation (EWC) [25] is a form of feature regularization
method for alleviating catastrophic forgetting. Let $\Phi_{A}$ be the set of
hashed feature embeddings learned on task $\mathcal{I}_{A}$. Consider a
subsequent task $\mathcal{I}_{B}$. EWC modifies the training objective to the
following:

$\mathcal{L}(\Phi)=\mathcal{L}_{B}(\Phi)+\frac{\lambda}{2}F(\Phi-\Phi_{A})^{2}.$
(5)

$\mathcal{L}_{B}$ represents the training objective on task $\mathcal{I}_{B}$
and $F$ is an estimation of the diagonal of the Fischer information matrix
given by the squared gradients of parameters $\Phi_{A}$ with respect to the
training objective $\mathcal{L}_{A}$. Intuitively, $\Phi_{A}$ is recorded as a
set of reference parameters. Deviation from these reference parameters are
penalized, weighted on their importance relative to the training objective. We
implement EWC on top of an iNGP backbone as a baseline method by fixing the
trained network parameters after each training task as the reference
parameters.

#### MEIL-NeRF.

Recently, MEIL-NeRF [13] also proposed the use of memory replay through
network distillation for alleviating catastrophic forgetting effects in NeRFs.
However, MEIL-NeRF uses the original fully implicit NeRF representation as a
backbone, which limits reconstruction quality and convergence speed. We
include continual learning results following the general implementation of
MEIL-NeRF. While MEIL-NeRF uses an additional ray generator network for
sampling previous rays from previous tasks, this additional step leads to
significant degradation in reconstruction results while providing minimal
memory savings; we therefore omit this step and sample ground truth rays
instead. MEIL-NeRF also explores using Charbonnier penalty function, we
consider changes to the penalty function a tangential area of exploration, and
choose to train both our method and MEIL-NeRF using the loss function detailed
in Equation 4.

###  5.2 Datasets

We compare methods on the task of continual scene fitting using the Tanks &
Temples [27], ScanNet [14] and TUM RGB-D datasets [52]. Data for each scene is
represented by a trajectory of ground truth camera poses and corresponding RGB
images, with each trajectory containing 100–300 images depending on scene. We
emulate the setting of continual learning by partitioning each trajectory into
10 temporally sequential tasks.

| ScanNet | Tanks & Temples | TUM RGB-D
---|---|---|---
Method |  0101 |  0146 |  0160 |  Truck |  Caterpillar |  Family |  TUM 1 |  TUM 2
Ours (1 s) |  19.61 |  22.18 |  23.84 |  19.19 |  17.24 |  23.24 |  15.05 |  16.65
Ours (5 s) |  24.10  |  26.13  |  28.37  |  21.93  |  20.59  |  26.29  |  19.35  |  19.02
Ours (30 s) |  25.54  |  27.84  |  30.45  |  23.97  |  22.62  |  29.21  |  21.09  |  20.42
MEIL–NeRF (30 s) |  18.85 |  21.41 |  22.72 |  18.11 |  16.93 |  21.78 |  16.05 |  15.96
MEIL–NeRF (1 min) |  20.65 |  22.76 |  24.39 |  19.38 |  18.41 |  23.19 |  17.44 |  16.19
MEIL-NeRF (10 min) |  24.32  |  26.82  |  28.93  |  22.74  |  20.89  |  26.57  |  20.79  |  19.80
Table 2: Quantitative results: time constrained. We show reconstruction PSNR
for our method and MEIL-NeRF trained on a fixed time limit per task. Our
method converges to better results at a much faster rate. Our method trained
for only 5s per task outperforms MEIL-NeRF trained for 1 min per task and is
competitive with MEIL-NeRF trained for 10 min per task. Given its rapid
convergence, our method uniquely enables real-time continual scene
reconstruction. Figure 5: Qualitative results: time constrained. We show
reconstructed views from an earlier supervised (forgotten) task for our method
and MEIL-NeRF trained for fixed times per task. Our method consistently
outperforms MEIL-NeRF given equal time budget. With only 5s per task, our
method already reconstructs the scene with reasonable fidelity, illustrating
that our method is well-suited for real-time continual scene fitting.

###  5.3 Experimental settings

We evaluate our method in two separate settings: an unconstrained setting
where each method is trained on every task until convergence, and a
constrained setting where each task is trained on a fixed time budget. The
unconstrained setting aims at testing the upper-bound performance of each
method, while the constrained setting mimics a real-time continual scene
reconstruction setting. Each model is trained on a single RTX 6000 GPU, with a
ray batch size of 1024. For the unconstrained settings, we trained methods
using a hash encoding for 1 minute per task and methods built on fully
implicit NeRFs for 10 minutes per task.

###  5.4 Results

#### Unconstrained setting.

We show quantitative results of each method for the unconstrained setting in
Table 1. Methods are evaluated using peak signal-to-noise ratio (PSNR),
averaged over all images in the test trajectory. We also provide quantitative
results of the fully implicit NeRF and hash-encoding representations trained
in a non-continual setting. These results serve as an upper bound for their
continual learning counterparts. Quantitatively, our method consistently
outperforms all baselines in reconstruction quality. While performance of
MEIL-NeRF comes close to our method for certain scenes, our method takes
significantly less time to train due to the convergence properties of the hash
encoding representation. Results from our method also come very close to the
theoretical upper bound set by the results obtained from non-continual
training, further illustrating the effectiveness of our method.

Figure 4 shows qualitative results from the unconstrained setting. Naively
training NeRF under the continual setting leads to catastrophic forgetting, as
earlier views contain heavy artefacts. Parameter regularization through EWC
helps alleviate forgetting for certain scenes, however, reconstruction quality
is still limited. MEIL-NeRF produces visually pleasing results, but
reconstruction of earlier views lack high-frequency details. In contrast, our
method is able to retain these high frequency details, as the underlying
multi-resolution hash encoding stores high-frequency features explicitly,
allowing high frequency details to be retained during training.

#### Time-constrained setting.

We evaluate our method against MEIL-NeRF in the time-constrained setting. We
trained both methods on each task for a fixed period of time and show
reconstruction PSNR averaged over all views along the test trajectory in Table
2. Our method trained on 30 seconds per task out performs MEIL-NeRF, even when
trained for 10 minutes per task. More importantly, our method trained for just
5 seconds produces results comparable to MEIL-NeRF at convergence. Qualitative
results in Figure 5 show that our method provides reasonable scene
reconstruction quality at much shorter training times, illustrating that our
method is uniquely suited for the task of real-time continual scene fitting.

###  5.5 Ablation study

#### Degradation of early tasks.

We evaluate reconstruction PSNR of the second task (of ten total) over the
course of training using different methods in Figure 6. Conventional methods
naively trained in the continual setting (NeRF-Incre & iNGP-Incre) experience
severe degredation due to catasatrophic forgetting. MEIL-NeRF succesfully
alleviates forgetting effects through self-distillation, however, forgetting
effects are still observed after training for many tasks. In contrast, our
method is able to maintain high PSNR for previous tasks even after training
for many tasks.

Figure 6: Reconstruction quality of early-supervised tasks. Reconstruction
PSNR of task 2 over the course of training. Our method successfully alleviates
degradation effects from catastrophic forgetting, and consistently outperforms
all other baselines.

###  5.6 Applications: autonomous vehicle data

Our method is well-suited for scenarios such as autonomous vehicle captures
and drone footage, where data is sequentially acquired and an updated 3D scene
representation should be immediately available. To illustrate this, we train
our method on data obtained from the Waymo open dataset [54]. A single
trajectory in the Waymo dataset consists of a video stream from 5 calibrated
cameras mounted at the top of the vehicle. Similar to the experimental
settings described in Section 5, we split each trajectory into 10 temporally
sequential tasks. We show qualitative results using our method and iNGP-Incre
in Figure 7, training each task for 30 seconds for a total of 5 minutes. Our
method recovers meaningful geometry and reconstructs earlier views with much
higher quality.

##  6 Discussion

#### Limitations and future work.

Our method relies on ground truth camera poses to perform scene fitting.
Although prior works have explored simultaneous optimization of camera poses
and scene parameters for NeRFs, they either rely on good initializations [32]
or specific constraints on the distribution of camera poses [30]. Simultaneous
estimation of camera poses in the setting of continual learning for NeRFs has
also yet to be explored. It may be fruitful to explore this direction further
in relation to methods for SLAM [3].

We chose to use multi-resolution hash encodings [40] to leverage its fast
convergence properties and explicitly defined features to combat forgetting.
Alternate representations, such as triplanes [8] and TensoRF [10] can also be
explored as potential substitutes, potentially further increasing robustness
to catastrophic forgetting through more structured encodings.

Our method uses a frozen version of the scene representation network trained
on previous tasks as a pseudo ground truth oracle. Querying the network for
pseudo ground truth values requires volume rendering through the scene, adding
computational overhead to the training process. A potential direction of
exploration is to find other forms of compression, such as 2D coordinate
networks [50], to act as the pseudo ground truth oracle. Additionally, if any
single oracle network is not of sufficient quality, this will continue to
affect downstream training on subsequent tasks.

Figure 7: Qualitative results on Waymo open dataset. Our method recovers
earlier training views at higher quality than training NeRFs naively in a
continual learning setting.

#### Conclusion.

In this work, we aim to extend the practical viability of NeRFs, specifically
in the continual setting, where training data is sequentially captured and a
3D representation needs to be immediately available. By combining multi-
resolution hash encodings and replay methods through network distillation, our
approach alleviates the effects of catastrophic forgetting observed in the
continual learning of NeRFs. While previous approaches struggle with quality
and speed, our method is able to produce visually compelling reconstruction of
earlier tasks while being an order of magnitude faster than existing methods.

##  7 Acknowledgements

We thank Geoff Burns, Abhinav Modi, Jing Cui, and Hamid Izadi Nia for
invaluable discussions and feedback. This project was in part supported by
Rivian, the Samsung GRO program and a PECASE from the ARO. Ryan Po is
supported by the Stanford Graduate Fellowship.

## References

  * [1] Rahaf Aljundi, Francesca Babiloni, Mohamed Elhoseiny, Marcus Rohrbach, and Tinne Tuytelaars.  Memory aware synapses: Learning what (not) to forget.  ArXiv, abs/1711.09601, 2017.
  * [2] Rahaf Aljundi, Punarjay Chakravarty, and Tinne Tuytelaars.  Expert gate: Lifelong learning with a network of experts.  2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 7120–7129, 2016.
  * [3] Tim Bailey and Hugh F. Durrant-Whyte.  Simultaneous localization and mapping (slam): part ii.  IEEE Robotics & Automation Magazine, 13:108–117, 2006.
  * [4] Jonathan T. Barron, Ben Mildenhall, Matthew Tancik, Peter Hedman, Ricardo Martin-Brualla, and Pratul P. Srinivasan.  Mip-nerf: A multiscale representation for anti-aliasing neural radiance fields.  2021 IEEE/CVF International Conference on Computer Vision (ICCV), pages 5835–5844, 2021.
  * [5] Jonathan T. Barron, Ben Mildenhall, Dor Verbin, Pratul P. Srinivasan, and Peter Hedman.  Mip-nerf 360: Unbounded anti-aliased neural radiance fields.  2022 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 5460–5469, 2021.
  * [6] Alexander W. Bergman, Petr Kellnhofer, and Gordon Wetzstein.  Fast training of neural lumigraph representations using meta learning.  In NeurIPS, 2021.
  * [7] Fabio Bruno, Stefano Bruno, G. De Sensi, Maria Laura Luchi, Stefania Mancuso, and Maurizio Muzzupappa.  From 3d reconstruction to virtual reality: A complete methodology for digital archaeological exhibition.  Journal of Cultural Heritage, 11:42–49, 2010.
  * [8] Eric Chan, Connor Z. Lin, Matthew A. Chan, Koki Nagano, Boxiao Pan, Shalini De Mello, Orazio Gallo, Leonidas J. Guibas, Jonathan Tremblay, S. Khamis, Tero Karras, and Gordon Wetzstein.  Efficient geometry-aware 3d generative adversarial networks.  2022 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 16102–16112, 2021.
  * [9] Arslan Chaudhry, Marc’Aurelio Ranzato, Marcus Rohrbach, and Mohamed Elhoseiny.  Efficient lifelong learning with a-gem.  ArXiv, abs/1812.00420, 2018.
  * [10] Anpei Chen, Zexiang Xu, Andreas Geiger, Jingyi Yu, and Hao Su.  Tensorf: Tensorial radiance fields.  In European Conference on Computer Vision, 2022.
  * [11] Anpei Chen, Zexiang Xu, Fuqiang Zhao, Xiaoshuai Zhang, Fanbo Xiang, Jingyi Yu, and Hao Su.  Mvsnerf: Fast generalizable radiance field reconstruction from multi-view stereo.  2021 IEEE/CVF International Conference on Computer Vision (ICCV), pages 14104–14113, 2021.
  * [12] Jianchuan Chen, Ying Zhang, Di Kang, Xuefei Zhe, Linchao Bao, and Huchuan Lu.  Animatable neural radiance fields from monocular rgb video.  ArXiv, abs/2106.13629, 2021.
  * [13] Jaeyoung Chung, Kang-Ho Lee, Sungyong Baik, and Kyoung Mu Lee.  Meil-nerf: Memory-efficient incremental learning of neural radiance fields.  ArXiv, abs/2212.08328, 2022.
  * [14] Angela Dai, Angel X. Chang, Manolis Savva, Maciej Halber, Thomas A. Funkhouser, and Matthias Nießner.  Scannet: Richly-annotated 3d reconstructions of indoor scenes.  2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 2432–2443, 2017.
  * [15] Hugh F. Durrant-Whyte and Tim Bailey.  Simultaneous localization and mapping: part i.  IEEE Robotics & Automation Magazine, 13:99–110, 2006.
  * [16] Chrisantha Fernando, Dylan S. Banarse, Charles Blundell, Yori Zwols, David R Ha, Andrei A. Rusu, Alexander Pritzel, and Daan Wierstra.  Pathnet: Evolution channels gradient descent in super neural networks.  ArXiv, abs/1701.08734, 2017.
  * [17] Robert M. French.  Catastrophic forgetting in connectionist networks.  Trends in Cognitive Sciences, 3:128–135, 1999.
  * [18] Guy Gafni, Justus Thies, Michael Zollhofer, and Matthias Nießner.  Dynamic neural radiance fields for monocular 4d facial avatar reconstruction.  2021 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 8645–8654, 2020.
  * [19] Zekun Hao, Arun Mallya, Serge J. Belongie, and Ming-Yu Liu.  Gancraft: Unsupervised 3d neural rendering of minecraft worlds.  2021 IEEE/CVF International Conference on Computer Vision (ICCV), pages 14052–14062, 2021.
  * [20] Peter Hedman, Pratul P. Srinivasan, Ben Mildenhall, Jonathan T. Barron, and Paul E. Debevec.  Baking neural radiance fields for real-time view synthesis.  2021 IEEE/CVF International Conference on Computer Vision (ICCV), pages 5855–5864, 2021.
  * [21] Geoffrey E. Hinton, Oriol Vinyals, and Jeffrey Dean.  Distilling the knowledge in a neural network.  ArXiv, abs/1503.02531, 2015.
  * [22] Jeffrey Ichnowski, Yahav Avigal, Justin Kerr, and Ken Goldberg.  Dex-nerf: Using a neural radiance field to grasp transparent objects.  In Conference on Robot Learning, 2021.
  * [23] Joel Janai, Fatma Güney, Aseem Behl, and Andreas Geiger.  Computer vision for autonomous vehicles: Problems, datasets and state-of-the-art.  ArXiv, abs/1704.05519, 2017.
  * [24] Petr Kellnhofer, Lars Jebe, Andrew Jones, Ryan P. Spicer, Kari Pulli, and Gordon Wetzstein.  Neural lumigraph rendering.  2021 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 4285–4295, 2021.
  * [25] James Kirkpatrick, Razvan Pascanu, Neil C. Rabinowitz, Joel Veness, Guillaume Desjardins, Andrei A. Rusu, Kieran Milan, John Quan, Tiago Ramalho, Agnieszka Grabska-Barwinska, Demis Hassabis, Claudia Clopath, Dharshan Kumaran, and Raia Hadsell.  Overcoming catastrophic forgetting in neural networks.  Proceedings of the National Academy of Sciences, 114:3521 – 3526, 2016.
  * [26] Georg S. W. Klein and David William Murray.  Parallel tracking and mapping for small ar workspaces.  2007 6th IEEE and ACM International Symposium on Mixed and Augmented Reality, pages 225–234, 2007.
  * [27] Arno Knapitsch, Jaesik Park, Qian-Yi Zhou, and Vladlen Koltun.  Tanks and temples.  ACM Transactions on Graphics (TOG), 36:1 – 13, 2017.
  * [28] Kiriakos N. Kutulakos and Steven M. Seitz.  A theory of shape by space carving.  International Journal of Computer Vision, 38:199–218, 1999.
  * [29] Matthias De Lange, Rahaf Aljundi, Marc Masana, Sarah Parisot, Xu Jia, Ale _vs_. Leonardis, Gregory G. Slabaugh, and Tinne Tuytelaars.  A continual learning survey: Defying forgetting in classification tasks.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 44:3366–3385, 2019.
  * [30] Axel Levy, Mark J. Matthews, Matan Sela, Gordon Wetzstein, and Dmitry Lagun.  Melon: Nerf with unposed images using equivalence class estimation.  ArXiv, abs/2303.08096, 2023.
  * [31] Zhizhong Li and Derek Hoiem.  Learning without forgetting.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 40:2935–2947, 2016.
  * [32] Chen-Hsuan Lin, Wei-Chiu Ma, Antonio Torralba, and Simon Lucey.  Barf: Bundle-adjusting neural radiance fields.  2021 IEEE/CVF International Conference on Computer Vision (ICCV), pages 5721–5731, 2021.
  * [33] Lingjie Liu, Jiatao Gu, Kyaw Zaw Lin, Tat-Seng Chua, and Christian Theobalt.  Neural sparse voxel fields.  ArXiv, abs/2007.11571, 2020.
  * [34] Stephen Lombardi, Tomas Simon, Gabriel Schwartz, Michael Zollhoefer, Yaser Sheikh, and Jason M. Saragih.  Mixture of volumetric primitives for efficient neural rendering.  ACM Transactions on Graphics (TOG), 40:1 – 13, 2021.
  * [35] David Lopez-Paz and Marc’Aurelio Ranzato.  Gradient episodic memory for continual learning.  In NIPS, 2017.
  * [36] Xinzhu Ma, Zhihui Wang, Haojie Li, Wanli Ouyang, and Pengbo Zhang.  Accurate monocular 3d object detection via color-embedded 3d reconstruction for autonomous driving.  2019 IEEE/CVF International Conference on Computer Vision (ICCV), pages 6850–6859, 2019.
  * [37] Arun Mallya and Svetlana Lazebnik.  Packnet: Adding multiple tasks to a single network by iterative pruning.  2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 7765–7773, 2017.
  * [38] Ricardo Martin-Brualla, Noha Radwan, Mehdi S. M. Sajjadi, Jonathan T. Barron, Alexey Dosovitskiy, and Daniel Duckworth.  Nerf in the wild: Neural radiance fields for unconstrained photo collections.  2021 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 7206–7215, 2020.
  * [39] Ben Mildenhall, Pratul P. Srinivasan, Matthew Tancik, Jonathan T. Barron, Ravi Ramamoorthi, and Ren Ng.  Nerf.  Communications of the ACM, 65:99 – 106, 2021.
  * [40] Thomas Müller, Alex Evans, Christoph Schied, and Alexander Keller.  Instant neural graphics primitives with a multiresolution hash encoding.  ACM Transactions on Graphics (TOG), 41:1 – 15, 2022.
  * [41] Raul Mur-Artal, José M. M. Montiel, and Juan D. Tardós.  Orb-slam: A versatile and accurate monocular slam system.  IEEE Transactions on Robotics, 31:1147–1163, 2015.
  * [42] Michael Oechsle, Songyou Peng, and Andreas Geiger.  Unisurf: Unifying neural implicit surfaces and radiance fields for multi-view reconstruction.  2021 IEEE/CVF International Conference on Computer Vision (ICCV), pages 5569–5579, 2021.
  * [43] Sida Peng, Junting Dong, Qianqian Wang, Shang-Wei Zhang, Qing Shuai, Xiaowei Zhou, and Hujun Bao.  Animatable neural radiance fields for modeling dynamic human bodies.  2021 IEEE/CVF International Conference on Computer Vision (ICCV), pages 14294–14303, 2021.
  * [44] Sylvestre-Alvise Rebuffi, Alexander Kolesnikov, G. Sperl, and Christoph H. Lampert.  icarl: Incremental classifier and representation learning.  2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 5533–5542, 2016.
  * [45] Christian Reiser, Songyou Peng, Yiyi Liao, and Andreas Geiger.  Kilonerf: Speeding up neural radiance fields with thousands of tiny mlps.  2021 IEEE/CVF International Conference on Computer Vision (ICCV), pages 14315–14325, 2021.
  * [46] Anthony V. Robins.  Catastrophic forgetting, rehearsal and pseudorehearsal.  Connect. Sci., 7:123–146, 1995.
  * [47] David Rolnick, Arun Ahuja, Jonathan Schwarz, Timothy P. Lillicrap, and Greg Wayne.  Experience replay for continual learning.  In Neural Information Processing Systems, 2018.
  * [48] Antoni Rosinol, John J. Leonard, and Luca Carlone.  Nerf-slam: Real-time dense monocular slam with neural radiance fields.  ArXiv, abs/2210.13641, 2022.
  * [49] Hanul Shin, Jung Kwon Lee, Jaehong Kim, and Jiwon Kim.  Continual learning with deep generative replay.  In NIPS, 2017.
  * [50] Vincent Sitzmann, Julien N. P. Martel, Alexander W. Bergman, David B. Lindell, and Gordon Wetzstein.  Implicit neural representations with periodic activation functions.  ArXiv, abs/2006.09661, 2020.
  * [51] Vincent Sitzmann, Michael Zollhöfer, and Gordon Wetzstein.  Scene representation networks: Continuous 3d-structure-aware neural scene representations.  Advances in Neural Information Processing Systems, 32, 2019.
  * [52] Jürgen Sturm, Nikolas Engelhard, Felix Endres, Wolfram Burgard, and Daniel Cremers.  A benchmark for the evaluation of rgb-d slam systems.  2012 IEEE/RSJ International Conference on Intelligent Robots and Systems, pages 573–580, 2012.
  * [53] Edgar Sucar, Shikun Liu, Joseph Ortiz, and Andrew J. Davison.  imap: Implicit mapping and positioning in real-time.  2021 IEEE/CVF International Conference on Computer Vision (ICCV), pages 6209–6218, 2021.
  * [54] Pei Sun, Henrik Kretzschmar, Xerxes Dotiwalla, Aurelien Chouard, Vijaysai Patnaik, Paul Tsui, James Guo, Yin Zhou, Yuning Chai, Benjamin Caine, Vijay Vasudevan, Wei Han, Jiquan Ngiam, Hang Zhao, Aleksei Timofeev, Scott M. Ettinger, Maxim Krivokon, Amy Gao, Aditya Joshi, Yu Zhang, Jonathon Shlens, Zhifeng Chen, and Dragomir Anguelov.  Scalability in perception for autonomous driving: Waymo open dataset.  2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 2443–2451, 2019.
  * [55] Matthew Tancik, Ben Mildenhall, Terrance Wang, Divi Schmidt, Pratul P. Srinivasan, Jonathan T. Barron, and Ren Ng.  Learned initializations for optimizing coordinate-based neural representations.  2021 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 2845–2854, 2020.
  * [56] Matthew Tancik, Pratul P. Srinivasan, Ben Mildenhall, Sara Fridovich-Keil, Nithin Raghavan, Utkarsh Singhal, Ravi Ramamoorthi, Jonathan T. Barron, and Ren Ng.  Fourier features let networks learn high frequency functions in low dimensional domains.  ArXiv, abs/2006.10739, 2020.
  * [57] Matthew Tancik, Ethan Weber, Evonne Ng, Ruilong Li, Brent Yi, Justin Kerr, Terrance Wang, Alexander Kristoffersen, Jake Austin, Kamyar Salahi, Abhik Ahuja, David McAllister, and Angjoo Kanazawa.  Nerfstudio: A modular framework for neural radiance field development.  ArXiv, abs/2302.04264, 2023.
  * [58] Ayush Tewari, Ohad Fried, Justus Thies, Vincent Sitzmann, Stephen Lombardi, Kalyan Sunkavalli, Ricardo Martin-Brualla, Tomas Simon, Jason Saragih, Matthias Nießner, et al.  State of the art on neural rendering.  In Computer Graphics Forum, volume 39, pages 701–727. Wiley Online Library, 2020.
  * [59] Ayush Tewari, Justus Thies, Ben Mildenhall, Pratul Srinivasan, Edgar Tretschk, Wang Yifan, Christoph Lassner, Vincent Sitzmann, Ricardo Martin-Brualla, Stephen Lombardi, et al.  Advances in neural rendering.  In Computer Graphics Forum, volume 41, pages 703–735. Wiley Online Library, 2022.
  * [60] A. Triki, Rahaf Aljundi, Matthew B. Blaschko, and Tinne Tuytelaars.  Encoder based lifelong learning.  2017 IEEE International Conference on Computer Vision (ICCV), pages 1329–1337, 2017.
  * [61] Dor Verbin, Peter Hedman, Ben Mildenhall, Todd E. Zickler, Jonathan T. Barron, and Pratul P. Srinivasan.  Ref-nerf: Structured view-dependent appearance for neural radiance fields.  2022 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 5481–5490, 2021.
  * [62] Qianqian Wang, Zhicheng Wang, Kyle Genova, Pratul P. Srinivasan, Howard Zhou, Jonathan T. Barron, Ricardo Martin-Brualla, Noah Snavely, and Thomas A. Funkhouser.  Ibrnet: Learning multi-view image-based rendering.  2021 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 4688–4697, 2021.
  * [63] Yiheng Xie, Towaki Takikawa, Shunsuke Saito, Or Litany, Shiqin Yan, Numair Khan, Federico Tombari, James Tompkin, Vincent Sitzmann, and Srinath Sridhar.  Neural fields in visual computing and beyond.  Computer Graphics Forum, 41, 2021.
  * [64] Ju Xu and Zhanxing Zhu.  Reinforced continual learning.  ArXiv, abs/1805.12369, 2018.
  * [65] Gengshan Yang, Minh Vo, Natalia Neverova, Deva Ramanan, Andrea Vedaldi, and Hanbyul Joo.  Banmo: Building animatable 3d neural models from many casual videos.  2022 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 2853–2863, 2021.
  * [66] Alex Yu, Sara Fridovich-Keil, Matthew Tancik, Qinhong Chen, Benjamin Recht, and Angjoo Kanazawa.  Plenoxels: Radiance fields without neural networks.  2022 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 5491–5500, 2021.
  * [67] Alex Yu, Ruilong Li, Matthew Tancik, Hao Li, Ren Ng, and Angjoo Kanazawa.  Plenoctrees for real-time rendering of neural radiance fields.  2021 IEEE/CVF International Conference on Computer Vision (ICCV), pages 5732–5741, 2021.
  * [68] Alex Yu, Vickie Ye, Matthew Tancik, and Angjoo Kanazawa.  pixelnerf: Neural radiance fields from one or few images.  2021 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 4576–4585, 2020.
  * [69] Zihan Zhu, Songyou Peng, Viktor Larsson, Weiwei Xu, Hujun Bao, Zhaopeng Cui, Martin R. Oswald, and Marc Pollefeys.  Nice-slam: Neural implicit scalable encoding for slam.  2022 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), pages 12776–12786, 2021.

Appendix

##  Appendix A Additional Results

###  A.1 Degradation of early task

We show additional qualitative results on the degradation of earlier views
under the continual learning setting for each method in Supp. Fig. 9. Naive
methods such as NeRF-Incre and EWC suffer from catastrophic forgetting.
Although the early view is reconstructed well when the current view is part of
the training views of the current task, reconstruction quality quickly
degrades as the training task moves to other regions of the 3D scene. Our
method and MEIL-NeRF [13] maintains relatively high quality throughout all
stages of training.

Figure 8: Reconstructed view of earlier task during different stages of
training. Quantitative results shown in Figure 6 are obtained from evaluating
PSNR of the right-most column.

###  A.2 Ablation on # of tasks

We perform an ablation over the number of tasks used for the continual
learning setting, showing results in Table 3. As the number of tasks
increases, the continual learning setting becomes more challenging, causing
performance of our method to drop minimally as number of tasks increases.
Naive method of NeRF-Incre performs poorly in the continual learning setting
even for lower number of tasks.

| # of Tasks
---|---
Method | 2 | 5 | 10 | 20
NeRF-Incre | 12.30 | 16.40 | 13.70 | 13.43
Ours | 24.29 | 24.27 | 24.13 | 23.87
Table 3: Ablation on # of tasks in continual learning setting. Results
obtained from 0101 scene of the ScanNet dataset [14].

###  A.3 Reconstruction quality for individual tasks

Table 4 shows reconstruction PSNR for scene ScanNet-0101 for every single task
after training has concluded. As expected, non-replay baselines (NeRF-Incre,
iNGP-Incre, EWC) suffer from catastrophic forgetting, where earlier tasks are
reconstructed poorly, but having high reconstruction PSNR for the final task.
Replay based methods such as MEIL-NeRF and ours maintain high PSNR across all
tasks, with our method achieving better performance across all tasks.

Method | Task 1 | Task 2 | Task 3 | Task 4 | Task 5 | Task 6 | Task 7 | Task 8 | Task 9 | Task 10
---|---|---|---|---|---|---|---|---|---|---
NeRF-Incre | 13.92 | 13.17 | 11.01 | 10.94 | 12.04 | 13.84 | 15.96 | 18.73 | 17.43 | 25.01
iNGP-Incre | 14.86 | 15.02 | 12.89 | 15.15 | 16.50 | 18.21 | 18.40 | 21.07 | 18.68 | 25.37
EWC | 13.79 | 14.07 | 12.66 | 14.66 | 16.50 | 18.09 | 18.46 | 20.64 | 17.98 | 25.34
MEIL-NeRF | 23.84 | 24.29 | 22.95 | 22.85 | 24.10 | 24.72 | 25.47 | 26.31 | 25.57 | 24.68
Ours | 24.19 | 26.33 | 24.74 | 24.32 | 26.12 | 26.08 | 26.78 | 27.15 | 25.96 | 24.92
Table 4: Reconstruction PSNR of individual tasks after completion of training
over all tasks in the continual learning setting.

##  Appendix B Code and Implementation Details.

We implement our method as a modular component of the Nerfstudio pipeline
[57]. Our code includes custom data parsers for processing datasets to be
suited for the continual learning setting described in the paper. All code and
custom data formats will be made available upon acceptance of the paper.

Figure 9: Nerfstudio interface for continual learning of NeRFs. Current view
shows training for an initial task from the Waymo open dataset [54].


[2]Charles University, Sokolovská 83 186, 75 Praha 8, Czech Republic<EMAIL_ADDRESS><EMAIL_ADDRESS>)
This work was supported by Charles University Research programs No. PRIMUS/21/SCI/009 and UNCE/SCI/023, and by the Magica project ANR-20-CE29-0007 funded by the French National Research Agency.

A new method for solving non-autonomous ordinary differential equations is proposed, the method achieves spectral accuracy.
It is based on a new result which expresses the solution of such ODEs as an element in the so called $\star$-algebra.
This algebra is equipped with a product, the $\star$-product, which is the integral over the usual product of two bivariate distributions.
Expanding the bivariate distributions in bases of Legendre polynomials leads to a discretization of the $\star$-product and this allows for the solution to be approximated by a vector that is obtained by solving a linear system of equations.
The effectiveness of this approach is illustrated with numerical experiments.

§ INTRODUCTION

A new method for solving non-autonomous ordinary differential equations that achieves spectral accuracy <cit.> is proposed.
Consider a smooth function $\tilde{f}\in C^{\infty}$ that is analytic, then the ODE studied here is

\begin{equation}\label{eq:ODE}
\frac{d}{dt} \tilde{u}(t) = \tilde{f}(t) \tilde{u}(t), \quad \tilde{u}(-1) = 1, \quad t\in\mathcal{I}, \quad \mathcal{I}:=\left[-1,1\right].
\end{equation}

This method forms an essential building block in the development of a numerical method to solve the matrix ODE

\begin{equation*}
\frac{d}{dt} \tilde{U}(t) = \tilde{A}(t) \tilde{U}(t), \quad \tilde{U}(-1) = I, \quad t\in\mathcal{I},
\end{equation*}

for sparse large-to-huge matrix functions $\tilde{A}(t)$.
The matrix ODE appears in many applications, e.g., in nuclear magnetic resonance spectroscopy (NMR) [3].
In NMR $\tilde{A}(t) = - 2\pi \imath H(t)$, where $H(t)$ is the Hamiltonian of the system describing the dynamics of the nuclear spins of some sample in a magnetic field.
Hamiltonians appearing in NMR are of size $2^\ell \times 2^\ell$ for a system with $\ell$ spins and is usually sparse since spins only interact with close neighbors.

The new method for the scalar ODE (<ref>) is based on expressing this ODE and its solution $\tilde{u}(t)$ in the so called $\star$-algebra, which is equipped with the $\star$-product [5].
Section <ref> introduces the $\star$-product, which is an integral over two bivariate distributions, and the $\star$-algebra generated by this product.
In this algebra, the solution is given by a closed form expression [6].
In principle, this expression can be computed symbolically, however, in general, the computation is too complex for practical purposes.

Therefore, a numerical procedure is proposed that computes a discretization of this solution in the matrix algebra, equipped with the usual matrix-matrix product.
The key to going from the $\star$-algebra to the matrix algebra is finding a suitable discretization of the $\star$-product, which can be based on a quadrature rule [1] or on the expansion of the bivariate distributions in a basis of orthonormal polynomials.
The latter approach is followed in this paper, a basis of orthonormal Legendre polynomials is chosen and the resulting discretization is discussed in Section <ref>.
Section <ref> describes how to obtain an accurate approximation inside the matrix algebra and illustrates the method with some numerical experiments.

§ $\STAR$-PRODUCT

The solution $\tilde{u}(t)$ can be expressed as a closed form expression in the $\star$-algebra.
First, all functions and the ODE (<ref>) are represented in the space
$\mathcal{D}(\mathcal{I})$, which is spanned by all distributions $d$ of the form $d(t,s) = \tilde{d}(t,s)\Theta(t-s) + \sum_{i=0}^N \tilde{d}_i(t,s) \delta^{(i)}(t-s)$, with $N\geq 0$ and $\tilde{d},\tilde{d}_i$ are smooth bivariate functions over $\mathcal{I}\times \mathcal{I}$.
In $\mathcal{D}(\mathcal{I})$ the ODE becomes

\begin{equation}\label{eq:ODE_bivariate}
\frac{d}{dt} u(t,s) = f(t,s) u(t,s),\quad u(s,s) = 1,\quad t,s\in\left[-1,1\right].
\end{equation}

The function is $f(t,s)=\tilde{f}(t)\Theta(t-s)$, where the Heaviside function $\Theta(t-s) := \begin{cases}
1, \quad \text{if }t\geq s,\\
0, \quad \text{if } t<s\end{cases}$ enforces the starting time, which is given by the parameter $s$.
The space of smooth functions $\tilde{f}(t)$ multiplied with the Heaviside function $\Theta(t-s)$ is denoted by $C^\infty_{\Theta} := \{f: f(t,s) = \tilde{f}(t)\Theta(t-s), \quad \tilde{f}\in C^\infty\}$.
For $f,g\in C^{\infty}_{\Theta}$ the $\star$-product, defined as

\begin{equation}\label{eq:starProd}
f(t,s)\star g(t,s) := \int_{-1}^1 g(t,\tau) f(\tau,t) d\tau \in C^{\infty}_{\Theta},
\end{equation}

is closed.
In the larger space $\mathcal{D}(\left[-1,1\right])$ the $\star$-product is also closed and an inverse for this product, $f(t,s)^{-\star}$ such that $f(t,s)\star f(t,s)^{-\star} = \delta(t-s)$, exists under certain conditions on $f(t,s)$ [5].
The corresponding identity element is the Dirac impulse $\delta(t-s) = \begin{cases}
1, \quad \text{if }t= s,\\
0, \quad \text{else}
\end{cases}.$
The elements and operations which compose the $\star$-algebra are given in Table <ref>, for details we refer to [5].

Operation/element     Properties
$f(t,s) \star g(t,s) = \int_{-1}^1 g(t,\tau) f(\tau,t) d\tau$                  $\star$-product is closed in $\mathcal{D}(\mathcal{I})$

$f + g$               addition is closed in $\mathcal{D}(\mathcal{I})$

$1_{\star} = \delta(t-s)$ (Diract delta)      identity element for the $\star$-product

$f^{\star-1}(t,s)$                        $\star$-inverse [5]

$R_\star(f)(t,s):=(1_{\star}- f)^{\star-1}(t,s)$        $\star$-resolvent [5]

Summary of elements and operations which compose the $\star$-algebra. The functions $f,g,q\in \mathcal{D}(\left[-1,1\right])$.

In the $\star$-algebra the solution $\tilde{u}(t)$ to (<ref>) is given by evaluating the solution $u(t,s)$ in $s=-1$:

\begin{equation*}
\tilde{u}(t) = u(t,s)\vert_{s=-1}, \quad \text{with }u(t,s) = \Theta(t-s)\star R_{\star}(\tilde{f}(t)\Theta(t-s)).
\end{equation*}

Computing $u(t,s)$ symbolically is usually too complex, therefore we will discretize the problem and compute an approximation to the solution by numerical computation.

§ FROM $\STAR$-ALGEBRA TO MATRIX ALGEBRA

The key to going from the $\star$-algebra to the matrix algebra is discretizing the $\star$-product.
One way to discretize the $\star$-product is by the use of quadrature rules [1].
Another way is by expansion in a basis of orthonormal polynomials (ONPs), which is the topic of this paper.
A natural choice of ONPs is the sequence of Legendre polynomials, these polynomials are discussed in Section <ref>.
Using these polynomials as a basis, Section <ref> describes how distributions living in $C^{\infty}_\Theta$ can be expanded as a series and Section <ref> provides details on how to compute the coefficients in this series.
The discretization of the $\star$-product that follows from the Legendre basis expansion and the resulting matrix algebra are the topic of Section <ref>.

§.§ Legendre polynomials

The sequence of Legendre polynomials $\{p_k\}_k$ satisfies the orthogonality conditions

\begin{equation*}
\int_{-1}^1 p_k(x) p_{\ell}(x) dx \begin{cases}
=0,\quad \text{if } k\neq \ell\\
\neq 0,\quad \text{if } k=\ell
\end{cases}.
\end{equation*}

We choose to normalize this sequence such that $\int_{-1}^1 p^2_k(x)dx = 1$.
These orthonormal Legendre polynomials satisfy the property stated in Lemma <ref>, which is paramount to efficiently computing the Legendre series expansion of funtions in $C^\infty_\Theta$.

Consider orthonormal Legendre polynomials $\{p_\ell(t)\}_{l}$.
Then for $-1\leq\tau\leq 1$ the following equality holds, for $\ell>0$,

\begin{equation*}
\int_{-1}^{\tau} p_\ell(\rho)d\rho = \frac{1}{\sqrt{2\ell+1}} \left(\frac{1}{\sqrt{2\ell+3}}p_{\ell+1}(\tau) - \frac{1}{\sqrt{2\ell-1}} p_{\ell-1}(\tau)  \right)
\end{equation*}

and for $\ell=0 $

\begin{equation*}
\int_{-1}^{\tau} p_0(\rho)d\rho = \frac{1}{\sqrt{3}} p_1(\tau) + p_0(\tau).
\end{equation*}

The expansion of a given function $\tilde{f}(x)$ in a basis of Legendre polynomials is given by the series

\begin{equation*}
\tilde{f}(x) = \sum_{k=0}^{\infty} \tilde{f}_k p_k(x), \quad \text{with coefficients }\tilde{f}_k = \int_{-1}^1 \tilde{f}(\tau) p_k(\tau) d\tau.
\end{equation*}

For functions finite and continuous on $\mathcal{I}$ the Legendre series expansion is uniformly convergent [9].
All functions we encountered in our applications are entire functions, i.e., analytic on the whole complex plane.
The Legendre series expansions of entire functions converges faster than geometric [8].
For more details on the rate of convergence for analytic and differentiable functions see [9].

Consider the truncated Legendre series $\hat{f}_N(x) := \sum_{k=0}^N \tilde{f}_k p_k(x)$.
An upper bound for the error of $\hat{f}_N(x)$ to $\tilde{f}(x)$ can be obtained by noting that $\vert p_k(x)\vert \leq \sqrt{\frac{2k+1}{2}}$ on $\mathcal{I}$:

\begin{equation*}
\Vert \tilde{f}(x) - \hat{f}_N(x)  \Vert_\infty = \sum_{k=N+1}^\infty \hat{f}_k p_k(x) \leq \sum_{k=N+1}^\infty \vert \hat{f}_k\vert \sqrt{\frac{2k+1}{2}}.
\end{equation*}

The fast decay of the magnitude of the coefficients of smooth functions cancels out the square root growth as $k$ increases.
Thus, if the series $\hat{f}_N(x)$ contains all the coefficients above machine precision appearing in the Legendre series expansion of $\tilde{f}(x)$, it represents $\tilde{f}(x)$ up to high accuracy.
To be able to use the FFT, we will consider interpolating Legendre series instead of truncated Legendre series, the coefficients $\{f_k\}_{k=0}^{N-1}$ of the interpolating Legendre series can be computed using  [2] at a complexity of $\mathcal{O}(N\log^2 (N))$.
The accuracy of the interpolating Legendre series is expected to be close to the truncated Legendre series <cit.>, i.e., $\vert\tilde{f}_k-f_k\vert$ is small for $k=0,1,\dots N-1$.

§.§ Expansion of distributions living in $C^{\infty}_{\Theta}$

The expansion of $\tilde{f}(t) \Theta(t-s) = f(t,s)\in C^{\infty}_{\Theta}$ in Legendre bases is given by

\begin{equation*}
f(t,s) = \sum_{k=0}^{\infty} \sum_{\ell=0}^{\infty} f_{k,\ell} p_k(t) p_\ell(s), \quad \text{for } t\neq s, \quad \text{with coefficients } f_{k,\ell} = \int_{-1}^1 \int_{-1}^1 f(\tau,\rho) p_k(\tau) p_\ell(\rho) d\tau d\rho.
\end{equation*}

The coefficients form the coefficient matrix $F :=\left[f_{k,\ell}\right]_{k,\ell=1}^\infty$ which represents $f(t,s)$ in the bases of Legendre polynomials:

\begin{equation*}
f(t,s) \approx \begin{bmatrix}
p_0(t)& p_1(t) &p_2(t) &\dots
\end{bmatrix} F \begin{bmatrix}
\vdots
\end{bmatrix}.
\end{equation*}

In the sequel, we will work with a truncation of this double series, we consider the $M\times M$ leading principal submatrix $F_M$ of $F$, such that $f(t,s)$ is represented by

\begin{equation*}
f(t,s) \approx \sum_{k=0}^{M-1} \sum_{\ell=0}^{M-1} f_{k,\ell} p_k(t) p_\ell(s)= \begin{bmatrix}
p_0(t)& p_1(t)& p_2(t)& \dots & p_{M-1}(t)
\end{bmatrix} F_M \begin{bmatrix}
\vdots\\
\end{bmatrix}.
\end{equation*}

Symbolic computation of the coefficients $f_{k,\ell}$ can be slow.
In general, a numerical computation of the coefficients is needed, this can be achieved by, e.g., using a Gauss-Legendre quadrature rule to discretize the double integral.
The number of nodes required to achieve accuracy close to machine precision depends on the given function $\tilde{f}(t)$.
However, there is a more straightforward and more efficient approach, which is the topic of next section.

§.§ Computing basis coefficient matrices

The coefficient matrix $F_M$ can be computed up to high precision.
The procedure to compute $F_M$ requires as input only the function $\tilde{f}(t)$, size of the matrix $M$ and the chosen accuracy for the entries of $F_M$.
The coefficients are obtained in two steps.
First $\tilde{f}(t)\in C^{\infty}$ is represented, up to machine precision, by its interpolating Legendre series $\tilde{f}(t) \approx \sum_{k=0}^{N-1} {f}_k p_k(x)$.
This series is obtained by using  [2], which automatically chooses the number of terms $N$.
Second, since the series can be expressed as $f(t,s) \approx \sum_{d=0}^{N-1} \hat{f}_d \sum_{k=0}^{\infty}  \sum_{\ell=0}^{\infty} p_d(t) \Theta(t-s) p_k(t) p_\ell(s)$, it suffices to compute the coefficients for Legendre polynomials of degree $d$ in $C^\infty_{\Theta}$, i.e., $p_d(t)\Theta(t-s)$:

\begin{equation*}
b_{k,\ell}^{(d)} := \int_{-1}^1 \int_{-1}^1 p_d(\tau) \Theta(\tau-\rho) p_k(\tau) p_\ell(\rho) d\rho d\tau.
\end{equation*}

The corresponding basis coefficient matrices $\{B^{(d)}\}_{d=0}^{N-1}$ are infinite matrices $B^{(d)} := \left[b_{k,\ell}^{(d)} \right]_{k,l=0}^\infty$.
Note that these basis coefficient matrices do not change for different $f(t,s)$, we only require the expansion coefficients $\{{f}_d\}_{d=0}^{N-1}$ of $\tilde{f}(t)$ in a basis of orthonormal Legendre polynomials in order to compute $F\approx F^{(N)} = \sum_{d=0}^{N-1} \hat{f}_d B^{(d)}$, where superscript $(N)$ denotes the number of terms used in the Legendre expansion of $\tilde{f}(t)$.
This approximation is convenient since, as is shown in the following, it is possible to express $b_{k,\ell}^{(d)}$ analytically.

Consider integers $a,b,c\geq 0$ and set $s:=\frac{a+b+c}{2}$.
The integral of the product of three orthonormal Legendre polynomials is

\begin{equation*}
\mathcal{P}_{a,b,c} := \int_{-1}^{1} {p}_a(\rho) {p}_b(\rho) {p}_c(\rho) d\rho
\begin{cases}
0\quad \textnormal{if } a+b+c \textnormal{ odd or } s<\max(a,b,c), \textnormal{ or }a<\vert b-c \vert,\\ \frac{2(2a+1)(2b+1)(2c+1)}{\sqrt{8}(a+b+c)}\begin{pmatrix}
\end{pmatrix} \begin{pmatrix}
\end{pmatrix} \begin{pmatrix}
\end{pmatrix} \begin{pmatrix}
\end{pmatrix}^{-1} \enspace \textnormal{else}.
\end{cases}
\end{equation*}

Let $p_d(t)$ be the orthonormal Legendre polynomial of degree $d$ and ${\mathcal{P}}_{a,b,c}$ as in Property <ref>.
Then the coefficients $b_{k,l}^{(d)}$ of $p_d(t)\Theta(t-s)$ expanded in orthonormal Legendre bases $\{p_k(t)\}_k$ and $\{p_\ell(s)\}$, i.e., $p_d(t)\Theta(t-s) = \sum_{k=0}^\infty \sum_{\ell=0}^\infty b_{k,\ell}^{(d)} p_k(t) p_\ell(s)$ are given by

\begin{equation}\label{eq:LegBasisCoeffs}
b^{(d)}_{k,\ell}==  \frac{1}{\sqrt{2\ell+1}} \left[ \frac{1}{\sqrt{2\ell+3}}\mathcal{P}_{d,k,\ell+1} - \frac{1}{\sqrt{2\ell-1}}\mathcal{P}_{d,k,\ell-1}\right].
\end{equation}

By orthonormality of the Legendre polynomials, we can write the coefficients as

\begin{align*}
b^{(d)}_{k,\ell} &= \int_{-1}^1  p_d(\tau) p_k(\tau) \left(\int_{-1}^1\Theta(\tau-\rho)  p_\ell(\rho) d\rho\right) d\tau = \int_{-1}^1  p_d(\tau) p_k(\tau) \underbrace{\left(\int_{-1}^\tau p_\ell(\rho) d\rho\right)}_{\text{Apply Lemma \ref{prop:int_shift_Legendre}}} d\tau \\
&=  \frac{1}{\sqrt{2\ell+1}} \left[ \frac{1}{\sqrt{2\ell+3}}\int_{-1}^1  p_d(\tau) p_k(\tau) p_{\ell+1}(\tau) d\tau - \frac{1}{\sqrt{2\ell-1}}\int_{-1}^1  p_d(\tau) p_k(\tau)  p_{\ell-1}(\tau)  d\tau\right]\\
&=  \frac{1}{\sqrt{2\ell+1}} \left[ \frac{1}{\sqrt{2\ell+3}}\mathcal{P}_{d,k,\ell+1} - \frac{1}{\sqrt{2\ell-1}}\mathcal{P}_{d,k,\ell-1}\right].
\end{align*}

An immediate consequence of this formula is that $B^{(d)}$ is a banded matrix with bandwidth $(d+1)$ and therefore $F^{(N)}$ has bandwidth $N$.
In the sequel, the truncated coefficient matrix $F^{(N)}_M\in\mathbb{C}^{M\times M}$ is used, which is the $M\times M$ leading principal submatrix of $F^{(N)}$.
Choosing an adequate value of $M$ is outside the scope of this paper and is the subject of ongoing research.

§.§ A finite basis

Once $f(t,s)$ is replaced by its truncated coefficient matrix $F^{(N)}_M$, the symbolic operations in the $\star$-algebra are replaced by operations involving matrices.
Consider $f(t,s),g(t,s)$ and their truncated coefficient matrices $F^{(N)}_M, G^{(N)}_M$, respectively.
The coefficient matrix of the result of the $\star$-product $q(t,s) = f(t,s)\star g(t,s)$ is then approximately given by the coefficient matrix $Q^{(N)}_M = F^{(N)}_M G^{(N)}_M$, i.e., the usual matrix matrix product.
This is obtained by plugging in the truncated series in the $\star$-product:

\begin{align*}
&q(t,s)=f(t,s)\star g(t,s) = \int_{-1}^1 f(t,\tau) g(\tau,s) d\tau \\
&\approx \int_{-1}^1 \left(\sum_{k=0}^{M-1} \sum_{\ell=0}^{M-1} f_{k,\ell} p_k(t) p_\ell(\tau)\right) \left(\sum_{k=0}^{M-1} \sum_{\ell=0}^{M-1} g_{k,\ell} p_k(\tau) p_\ell(s)\right) d\tau\\
&= \begin{bmatrix}
p_0(t) & \dots & p_M(t)
\end{bmatrix}F^{(N)}_M
\begin{bmatrix}
\int_{-1}^1 p_0(\tau)p_0(\tau)d\tau & \dots & \int_{-1}^1 p_0(\tau)p_\infty(\tau)d\tau\\
\int_{-1}^1 p_1(\tau)p_0(\tau) d\tau& \dots & \int_{-1}^1 p_1(\tau)p_\infty(\tau)d\tau\\
\vdots & & \vdots \\
\int_{-1}^1 p_{M-1}(\tau) p_0(\tau)d\tau & \dots &\int_{-1}^1  p_M(\tau) p_{M-1}(\tau)d\tau
\end{bmatrix}
G^{(N)}_M \begin{bmatrix}
\vdots\\
\end{bmatrix} \\
&= \begin{bmatrix}
p_0(t) & \dots & p_{M-1}(t)
\end{bmatrix}
\underbrace{F^{(N)}_M G^{(N)}_M}_{=Q^{(N)}_M} \begin{bmatrix}
\vdots\\
\end{bmatrix}.
\end{align*}

Table <ref> describes the matrix algebra for the coefficient matrices of functions in $\mathcal{D}$($\mathcal{I}$).
This matrix algebra is, in some sense, the discretization of the $\star$-algebra.
The coefficient matrix for the expansion of $\Theta(t-s)$ in the Legendre bases is denoted by $H_M\in\mathcal{C}^{M\times M}$.
This matrix appears in the coefficient matrix $U^{(N)}_M = H_M (I_M-F_M^{(N)})^{-1}$ which represents the series approximating the solution $u(t,s) = \Theta(t-s)\star (1_{\star}- f)^{\star-1}(t,s)$.
This expression is used in the next section to approximate the solution $\tilde{u}(t)$.

$f(t,s)$     $F_M$
$\star$-operation/elements     matrix operation/elements
$q(t,s) = f(t,s) \star g(t,s)$                  $Q_M = F_M G_M$

$f + g$               $F_M+G_M$

$1_{\star} := \delta(t-s)$       $I_M$, identity matrix

$f^{\star-1}(t,s)$                        $F_M^{-1}$

$R_\star(f)(t,s):=(1_{\star}- f)^{\star-1}(t,s)$        $R(F_M) := (I_M-F_M)^{-1}$
Solution to ODE (<ref>)     Approximate solution
$u(t,s) = \Theta(t-s)\star (1_{\star}- f)^{\star-1}(t,s)$     $U_M = H_M (I_M-F_M)^{-1}$

Matrix algebra for the coefficient matrices of functions in $\mathcal{D}(\mathcal{I})$ and corresponding operations and elements in the $\star$-algebra.

§ APPROXIMATION IN MATRIX FRAMEWORK

For the solution $u(t,s)$ of (<ref>) the truncated coefficient matrix in the Legendre bases is given by $U_M^{(N)} = H_M (I_M-F_M^{(N)})^{-1}$.
Since the solution of interest $\tilde{u}(t) = u(t,s)\vert_{s=-1}$, its coefficients in Legendre basis $u^{(N)} = \begin{bmatrix}
u_0 & u_1 & \dots & u_{M-1}
\end{bmatrix}$ can be computed by solving the linear system $ (I_M-F_M^{(N)}) y =\begin{bmatrix}
p_0(-1) &  p_1(-1) & \dots & p_{M-1}(-1)
\end{bmatrix}^\top$ for $y$ and forming the product $u^{(N)}=H_M y$.
The first $L$ coefficients in $u\in\mathbb{C}^{M}$ can be computed up to high accuracy, which leads to the approximate series $\tilde{u}(t) \approx \sum_{k=0}^{L} u_k p_k(t)=:\hat{u}_L(t)$.

Only the first $L$ coefficients of $u^{(N)}$ are accurate because of the use of a truncated basis.
The linear system for $y$ involves an $N$-banded matrix $(I_M-F_M^{(N)})$, which is a truncation of infinite matrices.
Thanks to the bandedness of $(I_M-F_M^{(N)})$, the first $M-K$ coefficients in $y$ are computed accurately, i.e., as if these coefficients were computed using the infinite matrices.
Here, $K$ corresponds to the numerical bandedness of the matrix $(I_M-F_M^{(N)})^{-1}$.
The multiplication with the tridiagonal $H_M$ then results in the first $L = M-K-1$ coefficients in $u$ being accurately computed.

We illustrate this using an example, consider $\tilde{f}(t) = \cos(4t)$, which can be accurately represented in Legendre basis by the series $\sum_{d=0}^{21}f_d p_d(t)$.
Therefore, its coefficient matrix $F^{(22)}_{101}\in\mathbb{C}^{101\times 101}$ has band size $N=22$, Figure <ref> shows the entries of $F^{(22)}_{101}$ that are larger than machine precision $\epsilon_{\textrm{mach}} \approx 2.2204 e-16$.
The coefficient matrix $U_{101}^{(22)} = H_{101} (I_{101}-F_{101}^{(22)})^{-1}$ is also shown on this figure and we observe that it has a numerical band size of $K+1=29$ in the trailing part of the matrix.
This suggests that the first $M-K-1 = 71$ of $u^{(22)}$ are accurate.
Figure <ref> verifies that the first 71 computed coefficients $\{u_k\}_k$ by comparing them to the coefficients of the Legendre expansion of the known exact solution $\tilde{u}(t)$ and by comparing the Legendre series $\hat{u}_n(t):=\sum_{k=0}^n u_k p_k(t)$ to the exact solution $\tilde{u}(t)$ in 1000 equispaced nodes in $\mathcal{I}$.

Entries in $F^{(22)}_{101}$ (left) for $f(t)=\cos(4t)$ and $U^{(22)}_{101}$ (right) with an amplitude larger than machine precision, $\vert f_{i,j}\vert \leq \epsilon_{\textrm{mach}}$ and $\vert u_{i,j}\vert \leq \epsilon_{\textrm{mach}}$.

Left: Magnitude of coefficients representing the solution to the ODE with $\tilde{f}(t) = \cos(4t)$. The coefficients of the Legendre series of the exact solution $\tilde{u}(t)$ are shown as ${\color{blue}\circ}$ and the coefficients $u = H_M(I_M-F_M)^{-1}y$ obtained by the method described in this paper as ${\color{red}\ast}$.
Right: error of the series $\hat{u}_n(t) = \sum_{k=0}^{n} u_k p_k(t)$ compared to $\tilde{u(t)}$ measured in infinity norm for increasing $n$.

Consider now a more oscillatory function $\tilde{f}(t) = -2\pi \imath( 0.1 + \cos(6\pi(t+1))  + \cos(12\pi(t+1)))$ for which we repeat the experiment.
We choose $M = 601$ and the functions is represented by $\tilde{f}(t) \approx \sum_{d=0}^{74} f_d p_d(t)$.
Numerically we determine that the band of $(I_{601}-F^{75}_{601})^{-1}$ is $K = 196$, which suggests that $M-K-1 = 404$ coefficients are computed accurately.
This is verified in Figure <ref>, where the computed coefficients are compared with the exact Legendre coefficients and the accuracy of $\hat{u}_n(t)$ is compared to $\tilde{u}(t)$ in the infinity norm.

Left: Magnitude of coefficients representing the solution to the ODE with $\tilde{f}(t) = -2\pi \imath( 0.1 + \cos(6\pi(t+1))  + \cos(12\pi(t+1)))$. The coefficients of the Legendre series of the exact solution $\tilde{u}(t)$ are shown as ${\color{blue}\circ}$ and the coefficients $u^{(N)} = H_M(I_M-F_M^{(N)})^{-1}y$ obtained by the method described in this paper as ${\color{red}\ast}$.
Right: error of the series $\hat{u}_n(t) = \sum_{k=0}^{n} u_k p_k(t)$ compared to $\tilde{u(t)}$ measured in infinity norm for increasing $n$.

These numerical experiments illustrate that the proposed method is capable of solving scalar ODEs up to high accuracy.
This method can be generalized to the matrix ODE, where similar numerical behavior has been observed, however treating this case is outside the scope of this report.
In order to develop a numerical algorithm that can compete with the state-of-the-art methods, the coefficient computation must be performed very efficiently, an a priori estimate of the required size of basis $M$ and a procedure to automatically truncate the series $\sum_{k=0}^{n} u_k p_k(t)$ at an appropriate value of $n$.

[1]
S., Cipolla, S., Pozza, M., Redivo-Zaglia, and N., Van Buggenhout, A Lanczos-type procedure for tensors. Numer Algor (2022). Published online: <https://doi.org/10.1007/s11075-022-01351-6>.

[2]
T. A., Driscoll, N., Hale, and L. N., Trefethen, Chebfun guide, (Pafnuty Publications, Oxford, 2014). <www.chebfun.org/docs/guide/>.

[3]
S., Hafner, and H. W., Spiess, Advanced solid-state NMR spectroscopy of strongly dipolar coupled spins under fast magic angle spinning.
Concepts Magn. Reson., 10, 99-128 (1998). <https://doi.org/10.1002/(SICI)1099-0534(1998)10:2<99::AID-CMR3>3.0.CO;2-Q>.

[4]
J., Gillis, J., Jedwab, and D., Zeilberger, A Combinatorial Interpretation of the Integral of the Product of Legendre Polynomials.
SIAM J. Math. Anal. 19:6, 1455-1461 (1988). <https://doi.org/10.1137/0519109>.

[5]
P-L., Giscard, and S. Pozza, Lanczos-Like Algorithm for the Time-Ordered Exponential: The $\ast$-Inverse Problem.
Linear Algebra Appl 65, 807–827 (2020). <https://doi.org/10.21136/AM.2020.0342-19>.

[6]
P-L., Giscard, and S. Pozza, Tridiagonalization of systems of coupled linear differential equations with variable coefficients by a Lanczos-like method.
Appl Math 624, 153–173 (2021). <https://doi.org/10.1016/j.laa.2021.04.011>.

[7]
P-L., Giscard, and S. Pozza, A Lanczos-like method for non-autonomous linear ordinary differential equations.
Boll Unione Mat Ital (2022). Published online: <https://doi.org/10.1007/s40574-022-00328-6>.

[8]
L. N., Trefethen, Approximation Theory and Approximation Practice. (SIAM, Philadelphia, PA, 2013).

[9]
H., Wang, and S., Xiang, On the convergence rates of Legendre approximation.
Math. Comp. 81, 861-877 (2012). <https://doi.org/10.1090/S0025-5718-2011-02549-4>.