redacted<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Mixture-of-Depths: Dynamically allocating compute in transformer-based
language models

David Raposo  Sam Ritter  Google DeepMind  Blake Richards  Google DeepMind
McGill University & Mila  Timothy Lillicrap  Google DeepMind  Peter Conway
Humphreys  Google DeepMind  Adam Santoro

###### Abstract

Transformer-based language models spread FLOPs uniformly across input
sequences. In this work we demonstrate that transformers can instead learn to
dynamically allocate FLOPs (or _compute_) to specific positions in a sequence,
optimising the allocation along the sequence for different layers across the
model depth. Our method enforces a total compute budget by capping the number
of tokens ($k$) that can participate in the self-attention and MLP
computations at a given layer. The tokens to be processed are determined by
the network using a top-$k$ routing mechanism. Since $k$ is defined _a priori_
, this simple procedure uses a static computation graph with known tensor
sizes, unlike other conditional computation techniques. Nevertheless, since
the identities of the $k$ tokens are fluid, this method can expend FLOPs non-
uniformly across the time and model depth dimensions. Thus, compute
expenditure is entirely predictable in sum total, but dynamic and context-
sensitive at the token-level. Not only do models trained in this way learn to
dynamically allocate compute, they do so efficiently. These models match
baseline performance for equivalent FLOPS and wall-clock times to train, but
require a fraction of the FLOPs per forward pass, and can be upwards of 50%
faster to step during post-training sampling.

##  1 Introduction

Not all problems require the same amount of time or effort to solve.
Analogously, in language modeling not all tokens and sequences require the
same time or effort to accurately make a prediction. And yet, transformer
models expend the same amount of compute per token in a forward pass. Ideally,
transformers would use smaller total compute budgets by not spending compute
unnecessarily.

Conditional computation is a technique that tries to reduce total compute by
expending it only when needed (Bengio, 2013; Bengio et al., 2013, 2016).
Various algorithms offer solutions to when and how much compute should be used
(Ainslie et al., 2023; Fedus et al., 2022; Bapna et al., 2020). However,
general formulations of this challenging problem may not work well with
existing hardware constraints since they tend to introduce dynamic computation
graphs (Graves, 2016; Dehghani et al., 2018). The most promising conditional
computation methods may instead be those that are harmonious with our current
hardware stack, which prioritizes static computation graphs, and known tensor
sizes that are selected to maximize hardware utilization.

Here we consider the problem of language modeling using a static compute
budget that can be made less than that used by a vanilla transformer. The
network must learn how to dynamically allocate the available compute by making
decisions per-token, in each layer, about where to spend compute from the
available budget. In our implementation total compute is user defined and
unchanging prior to training, rather than being a function of the network’s
on-the-fly decisions. Thus, hardware efficiency gains—such as reduced memory
footprint, or reduced FLOPs per forward pass—can be anticipated and exploited
ahead of time. As we will show, these gains can be had without sacrificing
overall performance.

We leverage an approach akin to Mixture of Experts (MoE) transformers, in
which dynamic token-level routing decisions are made across the network depth.
Departing from MoE, we choose to either apply a computation to a token (as
would be the case for a standard transformer), or pass it through a residual
connection (remaining unchanged and saving compute). Also in contrast to MoE,
we apply this routing to both forward MLPs and multi-head attention. Since
this therefore also impacts the keys and queries we process, the routing makes
decisions not only about which tokens to update, but also which tokens are
made available to attend to. We refer to this strategy as Mixture-of-Depths
(MoD) to emphasize how individual tokens pass through different numbers of
layers, or blocks, through the depth of the transformer (see figure 1).

The MoD technique also allows one to trade-off performance with speed. On the
one hand, one can train an MoD transformer that improves upon vanilla
transformers by as much as $1.5\%$ on the final log probability training
objective for equivalent training FLOPs (isoFLOP), and while taking an
equivalent amount of wall-clock time to train. On the other hand, one can
train an MoD transformer that achieves training loss parity with an isoFLOP
optimal vanilla transformer, but which uses a fraction of the FLOPs (upwards
of 50%) per forward pass, and hence is faster to step. Together, these results
imply that MoD transformers learn to route intelligently (i.e., skipping
computations that are unnecessary) since they can achieve equal or better log
probabilities per sequence despite a smaller FLOP footprint per forward pass.

Figure 1: Mixture-of-Depths Transformer. As in mixture-of-experts (MoE)
transformers we use a router to choose among potential computational paths.
But unlike in MoE transformers the possible choices are a standard block’s
computation (i.e., self-attention and MLP) or a residual connection. Since
some tokens take this second route, Mixture-of-Depths (MoD) transformers have
a smaller total FLOP footprint compared to vanilla or MoE transformers. On the
top right is depicted a trained model’s routing decisions for a short sequence
truncated to 64 tokens for visualization purposes. When examining the choices
one can find tokens processed by later blocks’ layers, despite passing through
relatively few total blocks throughout the model’s depth. This is a unique
feature of MoD compared to conventional halting-based, or "early-exit"
conditional computation, which instead engage blocks serially, or vanilla
transformers, which engage every block.

##  2 Background

The transformer architecture has become the workhorse of a revolution in
practical artificial intelligence, bringing unprecedented capabilities at the
cost of expensive training runs and serving procedures. This has spurred
tremendous interest in making transformer architectures more efficient (Tay et
al., 2020; Gupta and Agrawal, 2021). One of the promising approaches is
conditional computation, whereby learned mechanisms determine when and how to
expend computation. This terminology was introduced by Bengio (2013), and the
concept was explored further over the next several years (Bengio et al., 2013;
Cho and Bengio, 2014; Graves, 2016; Jernite et al., 2017; Bengio et al., 2016;
Wang et al., 2017).

A wide variety of recent work has developed conditional computation methods
for transformers. Some of this work focuses on "early exiting", that is,
learning to decide when to end computation on a given token, allowing the
token to skip any remaining transformer layers after the exit decision is made
(Elbayad et al., 2019; Liu et al., 2021; Schuster et al., 2022). In MoD,
unlike in early-exit methods, a token can skip middle layers, then be updated
via self-attention with tokens that that have gone through all the middle
layers. We speculate that this might be a useful property.

Other work has developed methods for iterating transformer layers with shared
weights for an adaptive number of steps (Simoulin and Crabbé, 2021; Dehghani
et al., 2018). Bolya et al. (2023) developed a method for choosing tokens to
merge when running inference on a trained vision transformer which notably
requires no learning. Lei et al. (2023) make use of conditional computation in
a fine tuning setting by building on adapter approaches (He et al., 2021) to
learn to skip blocks of frozen pre-trained weights in favor of running only a
small fine-tuned adapter.

CoLT5 (Ainslie et al., 2023) uses conditional routing to select whether a
given token will pass through a heavy or light pathway for each feedforward
layer. Further, they use the same routing mechanism to select whether a token
will attend to all other tokens or to a select few, as in Guo et al. (2022).
Like MoD, CoLT5 uses soft top-k for making routing decisions. However, CoLT5
focuses on a encoder-decoder setting, and thus does need to contend with the
problem of efficient sequential decoding given the non-causal nature of the
top-k operation. In contrast, our current work with MoD focuses on the
decoder-only setting, and so we propose a predictive router to enable
efficient inference for conditional computation in transformers.

One successful formulation of conditional computation is the the "mixture-of-
experts" layer (MoE) as introduced by Shazeer et al. (2017). Developed
initially in the context of LSTMs, later work showed compelling empirical
results for MoE with transformers (Lepikhin et al., 2020; Fedus et al., 2022;
Zoph et al., 2022). Unlike other conditional computation approaches that try
to conserve or expend additional compute, MoE transformers use conditional
logic to route tokens to one of many expert MLPs while keeping total compute
expenditure constant. Our mixture-of-depths method can be thought of as using
the routing logic from MoE transformers, but rather than having multiple
experts, MoD deploys a single expert which can be dynamically skipped.

##  3 Implementing Mixture-of-Depths Transformers

Our high-level strategy is as follows:

  * •

Set a static compute budget that is less than that of an equivalent vanilla
transformer by limiting the number of tokens in a sequence that can
participate in a block’s computations (i.e., self-attention and subsequent
MLP). For example, while a vanilla transformer might permit all the tokens in
a sequence to participate in self-attention, we might limit the number to 50%
of the tokens in a sequence. See section 3.1.

  * •

Use a per-block router to emit a scalar weight for each token, which expresses
the router’s preference for that token to participate in a block’s
computations or to route around it. See section 3.2.

  * •

Identify the top-$k$ scalar weights (per sequence, per block) to select those
tokens that will participate in a block’s computations. Since precisely $k$
tokens will participate in the block’s computations, the computation graph and
tensor sizes remain static throughout training; it is merely the tokens’
participation that is dynamic and context-sensitive, as determined by the
router. See section 3.3.

We then discuss some complications when sampling post-training in section 3.5.

###  3.1 Defining a compute budget

To enforce a total compute budget per forward pass we leverage the notion of
capacity, which defines the total number of tokens that comprise the input to
a given computation (e.g., the tokens participating in self-attention, a given
expert in MoE transformers, etc). For example, the self-attention and MLP in
each vanilla transformer block have a capacity of $T$—the total number of
tokens across the sequence and batch. MoE transformers, on the other hand, use
a capacity less than $T$ per expert MLP so as to more evenly divide the total
compute across each expert. But, since they use multiple experts per block,
their total capacity is approximately equal to that of a vanilla transformer.

Generally, it is the token capacity that determines the total FLOPs for
transformers that use conditional computation, rather than the outcomes of any
routing decisions. This is because static-graph implementations account for
the worst-case scenarios decisions; e.g., a computation’s inputs will be
padded to its capacity amount even if relatively few tokens actually end up
routing to it, and/or tokens will be dropped from the computation if the
capacity is exceeded.

We can achieve our goal of using a smaller compute budget per forward pass
compared to a vanilla transformer by lowering the capacity of the
computations. However, using a smaller compute budget haphazardly will result
in a performance degradation. We hypothesize that certain tokens might not
require as much processing as others, and these tokens can be identified
through learning. Therefore, if the network learns to choose the right tokens
to fill up its capacities, then it may preserve its performance. In the
following we describe routing schemes that can be used for this purpose.

###  3.2 Routing around transformer blocks

We consider the setting whereby we route tokens to one of two computational
paths: (1) self-attention and MLP blocks, and (2) a residual connection. The
latter is computationally cheap, and results in a block output that is
entirely determined by the value of its input. The former path is
computationally expensive.

The total number of FLOPs per forward pass will be fewer than that in a
vanilla transformer if we set the capacity for path (1) to be anything less
than $T$ (the total number of tokens across the sequence and batch). For
example, if we were to set a block’s capacity to $\frac{T}{2}$ (i.e., half the
number of tokens as would be the case in a vanilla transformer) then query-
times-key matrix multiplication during self-attention becomes $25\%$ as FLOP-
intensive as in a vanilla transformer ($(\frac{T}{2})^{2}$ vs. $T^{2}$).
Similar calculations can determine the FLOP-savings for the MLP.

Intuitively, the total FLOPs per forward pass decreases (and the time to
complete a forward pass decreases) in proportion to how aggressively we shrink
the blocks’ capacities. However, downstream performance will also be affected
by how aggressively we shrink the blocks capacities, and by the routing
algorithm we implement.

At one extreme, if we leave each block’s capacity at $T$ and route every token
to (rather than _around_) each block, then we recover a vanilla transformer.
At the other extreme, if we set each block’s capacity to $0$ and route all
tokens _around_ each block, then we’re left with a very fast model that
doesn’t engage with the vast majority of the transformer’s parameters, and
undoubtedly has poor downstream performance. We hypothesize that somewhere
between these two extremes is an optimal model that is faster than a vanilla
Transformer and performs as well, if not better, all while being faster to
step.

###  3.3 Routing schemes

Naively, one can leverage stochasticity to route tokens, akin to layer or
block “dropout”. We present this routing scheme as a control, and will show
that it significantly under-performs relative to vanilla transformers.

We hypothesize that _learned_ routing is preferable. Intuitively, the network
should be able to learn which tokens require more or less processing than
others. If we are correct that Transformers often expend more compute than
they need to make their predictions, then it is an empirical question as to
how aggressively we can shrink each block’s capacity, and hence, how many
tokens we can afford to route around each block.

Figure 2: Routing schemes. Tokens are funnelled to the computational path of
their choice when using token-choice routing (left). If a given path exceeds
its capacity (e.g., more than two tokens in this example) then surplus tokens
must be dropped (purple token). The exact token that is ultimately dropped
depends on the precise implementation in the underlying code. For example,
priority is often given to those tokens that come earlier in the sequence or
batch order. With expert-choice routing (middle), precisely $k$ (in this case,
two) tokens are chosen per path using a top-$k$ mechanism across the tokens’
router weights. Here, tokens are dropped if they are not among the top-$k$
with respect to any given path (orange token), and some tokens may even be
funnelled to multiple paths (yellow token). In this work we deploy expert-
choice routing (right). However, because we use just a single path, we
_leverage_ the implicit knowledge that tokens will be dropped if $k$ is less
than the sequence length so that we can route tokens away from the self-
attention and MLP computations, thus expending fewer FLOPs in a given forward
pass of the model.

There are two learned routing schemes we consider (see figure 2): token-choice
and expert-choice. In token-choice routing, a router produces per-token
probability distributions across computational paths (e.g., across expert
identities in MoE Transformers). Tokens are then shuttled to the path they
prefer—i.e., that with the highest probability—and auxiliary losses ensure
that all tokens don’t converge to the same path. Token-choice routing can have
load balancing problems since there isn’t a guarantee that tokens divide
themselves appropriately between the possible paths. “Expert choice routing”
flips this recipe on its head: rather than having tokens choose the path they
prefer, each path instead chooses the top-$k$ tokens based on the tokens’
preferences. This ensures a perfect load balance since $k$ tokens are
guaranteed to be shuttled to each path. However, it could result in over- or
under-processing of some tokens, since some tokens may be among the top-$k$
for multiple paths, or for none of them.

We decided to leverage expert-choice routing for a few reasons. First, it
obviates the need for an auxiliary balancing loss. Second, since the top-$k$
operation depends on the magnitude of the router weights, this routing scheme
allows for relative routing weights to help determine which tokens most need
the block’s computations; routers can try to _ensure_ that the most critical
tokens are among the top-$k$ by setting their weight appropriately, which is
not possible with token-choice routing schemes. For our specific use-case,
wherein one computational path is essentially a null operation, it might be
critical that important tokens are routed away from the null operation. Third,
because we only route through two paths, a single top-$k$ operation can
efficiently split the tokens into two mutually exclusive sets, one for each
computational path, preventing the over- or under-processing problem mentioned
above.

###  3.4 Routing implementation

As a reminder of the high-level intuition, each token is processed by a router
to produce a scalar weight, and the top-$k$ weights are then used to choose
the token identities that will route through a transformer’s block, which
comprises self-attention and the subsequent MLP.

Suppose we have the set of token embeddings in a sequence of length $S$ for a
given layer $l$; that is $X^{l}=\\{x_{i}^{l}|i\text{ is an integer, }1\leq
i\leq S\\}$. The router weight for a given token embedding is a scalar
produced as a result of a linear projection,
$r_{i}^{l}=w_{\theta}^{T}x_{i}^{l}$.

Our goal is to use these router weights to determine the output of a block’s
computation of each token. Suppose $P_{\beta}(R^{l})$ is the $\beta$-th
percentile of the set of router weights $R^{l}$, where $\beta=1-C/S$ and $C$
is the user-defined capacity per batch element (an integer $<S$ that defines
the number of tokens from a sequence that will be processed by a given
function). A block’s output for a given token is:

$\displaystyle
x_{i}^{l+1}=\begin{cases}r_{i}^{l}f_{i}(\tilde{X}^{l})+x_{i}^{l},&\text{if
}r_{i}^{l}>P_{\beta}(R^{l})\\\ x_{i}^{l},&\text{if
}r_{i}^{l}<P_{\beta}(R^{l})\\\ \end{cases}$ (1)

Here, $\tilde{X}^{l}$ is the set of tokens whose router values
$r_{i}^{l}>P_{\beta}(R^{l})$ (that is, the “top-k” tokens), and $f$ comprises
self-attention and the subsequent MLP. Note that the output for a given token
$x_{i}^{l+1}$ might depend on other tokens $x_{i\neq j}^{l}$ because of the
self-attention operation. The cardinality of $\tilde{X}^{l}$ is $C$ (or $k$):
the user-defined capacity. Therefore, the mixture-of-depths transformer
accrues compute savings relative to the baseline because the input to the
block’s computations $f$ comprise fewer tokens than usual ($C<S$), rendering
the self-attention and MLP less expensive.

Notably, we multiply the output of the function $f$ by the router weights.
This puts the router weights along the “gradient path”, thus subjecting them
to the forces of gradient descent through the course of the language modeling
task (We experimented with versions where the router weights are also included
along the computational path for those tokens that bypass the block’s
computations, but it seems to be sufficient—and implementationally simpler—to
only include the router weights along the computational path for those tokens
that do not bypass the block’s computations).

###  3.5 Sampling

While expert-choice routing has a number of advantages, it has one distinct
problem: the top-$k$ operation is non-causal. This means that whether a given
token’s routing weight is among the top-$k$ for the sequence depends on the
values of the routing weights for tokens that come after it, which we don’t
have access to when autoregressively sampling.

We tested two methods to work around this problem. The first introduces a
simple auxiliary loss that empirically affects the primary language modeling
objective by approximately $0.2-0.3\%$, but allows us to sample from the model
autoregressively. We use a binary cross-entropy loss wherein the router’s
outputs provide the logits, and the top-$k$ selections of these logits provide
the targets (i.e. 1 if a token was among the top-$k$, and 0 if not).
Intuitively, this loss centers the sigmoid of the router’s outputs around
$0.5$; those tokens that are selected among the top-k are pressured to produce
router outputs above $0.5$, and those not among the top-k will be pressured to
produce router outputs below $0.5$. The second method introduces a small
auxiliary MLP predictor (akin to a second router) that receives the same
inputs as the router (with a stop gradient), but whose output is a prediction
whether that token will be among the top-$k$ or not in the sequence. This
method does not affect the language modeling objective, and empirically does
not significantly impact the step speed.

Equipped with these new methods, we can sample autoregressively by choosing to
route tokens to or around a block based on the router’s output, which does not
depend on any information from future tokens. We provide empirical evidence
that this is a relatively easy auxiliary task that quickly achieves $99\%$
accuracy.

###  3.6 Training methods

All models use the same basic hyperparameter configurations (e.g. cosine
schedules equal to $1\times$ the training steps, 128 batch size, 2048 sequence
length) except for changes to the number of layers, heads, and embedding size
to produce differently sized models during isoFLOP analyses.

##  4 Results

###  4.1 Training, isoFLOP comparisons

Figure 3: MoD hyperparameter tuning. Variants of the MoD transformer were
trained for 6e18 FLOPs to determine the optimal hyperparameters for further
isoFLOP analyses. On the left plot, the grey box indicates models that perform
better than the isoFLOP optimal baseline. We found the best MoD variant to be
that which has the option to route every other block, and which uses a top-k
of $256$ (so, $256$, or 12.5% of the sequence’s tokens are processed by self-
attention and the subsequent MLP, while $1792$ tokens, or $87.5\%$ of the
sequence’s tokens route _around_ the block). Shown on the right are the
learning curves for a selection of models. Notably, model #3 achieves equal
performance to the isoFLOP optimal baseline but steps 66% faster, due to the
relatively fewer FLOPs needed per forward pass.

We first trained models with a relatively small FLOP budget (6e18) to
determine optimal hyperparameters (see figure 3). In general, we found that
MoD transformers drag the baseline isoFLOP curve "down and to the right". That
is, the optimal MoD transformer achieves a lower loss than the optimal
baseline, and also has more parameters. A fortunate consequence of this effect
is that there exist smaller MoD models that, while they are not themselves
isoFLOP optimal for their hyperparameter setting, are nevertheless as- or
better-performing than the optimal baseline model while being faster to step.
For example, a 220M parameter MoD (figure 3 model #3) variant slightly
outperforms the isoFLOP optimal baseline (also 220M, figure 3 model #1), but
is upwards of 60% faster to step during training. Crucially, when run on
equivalent hardware these two model variants take take approximately the same
amount of wall-clock time to train (figure 3).

We tested routing every block or every other block, using capacities from
12.5% to 95% of the total sequence. While routing every other block was
crucial for strong performance, we found that aggressive capacity reduction
was best (gradual improvements were observed when reducing the capacity down
to 12.5% of the total sequence, corresponding to 87.5% of tokens routing
around blocks, with performance degrading beyond this point). So, it seems the
network is robust to significant capacity reductions as long as there is
frequent opportunity for full capacity self-attention and MLP computations.

Learned routing is crucial, as MoD transformers that use stochastic routing
(implemented using a top-$k$ operation on router weights sampled from a
Gaussian distribution) perform drastically worse than both the baseline and
normal MoD transformer (figure 3).

Figure 4: isoFLOP analysis. We used the 12.5% capacity MoD variant to perform
an isoFLOP analysis for 6e18, 2e19, and 1e20 FLOPs, training models varying in
size from 60M to 3B parameters. Depicted on the right are the relative FLOPs
per forward pass (normalized to the isoFLOP optimal baseline). There exist MoD
variants that are both faster to step (by virtue of requiring fewer FLOPs per
forward pass) and better performing than the isoFLOP optimal baseline.

Depicted in figure 4 is an isoFLOP analysis for 6e18, 2e19, and 1e20 total
FLOPs. The trend that FLOP-optimal MoD transformers have more parameters than
the baseline continues for these larger FLOP budgets. Notably, there exist MoD
variants that are appreciably faster to step than the isoFLOP-optimal baseline
(measured as steps per second when training on equivalent hardware) while also
achieving a lower loss (in figure 4 we depict normalized FLOPs per forward
pass rather than wall-clock step time per se, but from our experiments the two
are tightly correlated. A similar plot can be produced showing relative wall-
clock step times and the same basic trend is present).

Step-wise speed gains come from two sources. First, the FLOP-per-parameter
ratio in MoD transformers is less than in the baselines because some
proportion of tokens are routed around blocks. So, for a given model size, a
transformer requires fewer FLOPs per forward pass. Second, since isoFLOP-
optimal MoD transformers are both bigger and achieve a lower loss than the
isoFLOP-optimal baseline, there exist smaller MoD variants that perform as
well or better than the isoFLOP-optimal baseline, and these variants are
faster to step because they are smaller. Altogether, then, there exist MoD
transformers that perform as well as isoFLOP-optimal baselines and are faster
to step, both because they use fewer FLOPs per parameter and because they use
fewer parameters.

Figure 4 also reveals another important finding: the optimal MoD transformer
is that which uses as many FLOPs per forward pass as the isoFLOP optimal
baseline. This finding allows one to directly predict which sized MoD
transformer will perform optimally for a given isoFLOP training budget: one
just needs to tune the model size for a given MoD configuration (i.e.,
capacity and routing frequency) to produce a model that uses as many FLOPs per
forward pass as the isoFLOP-optimal baseline, and they will have the optimally
performing MoD variant for that configuration. Empirically, we find that it is
better to add depth than to add width when adding FLOPs to the model.

Nevertheless, while the FLOPs per forward pass determines which model will be
the isoFLOP optimal, it does not predict whether the optimal loss will improve
upon the baseline (see figure 3. Namely, the optimal capacity appears to be
empirically determinable. We found that it is best to use 12.5% capacity
blocks, every other block.

We noticed that MoD transformers had memory savings relative to equivalently
sized baseline models at larger sizes, with some variants requiring fewer
total devices (i.e., a smaller TPU topology). We did not study this
extensively, but we anticipate that as one scales to larger models, these
savings could be an important consideration when choosing model variants to
train, and could have significant positive effects in regards to the KV cache
size during autoregressive sampling.

Figure 5: Routing analysis. We trained an MoD transformer that interleaved
$12.5\%$ capacity routing blocks with full-attention blocks. As expected, the
number of tokens that route to (rather than around) a block is sparse in
routing blocks, though the network does sometimes preferentially route certain
tokens to each block along its depth. This can be seen in the left figure that
depicts routing decisions, where we observe a vertical band of dark blue
towards the end of the sequence. As expected, the distribution of router
weights are as the auxiliary loss dictates: approximately $12.5\%$ of weights
are above 0.5 and $87.5\%$ are below (histogram, right).

Figure 5 shows the routing decisions for an MoD transformer trained with
interleaved routing blocks. Despite aggressive routing around the blocks,
transformers are able to achieve performance improvements relative to
baselines. We observe patterns that might warrant further study; namely, some
tokens appear to engage each block along the transformer’s depth, while others
decide to route around blocks whenever possible. Preliminary analyses suggest
that the tokens that engage with blocks more frequently are correlated with
output predictions that have higher entropy, which possibly corresponds to
predictions that are more difficult to make.

###  4.2 Auto-regressive Evaluation

Figure 6: Auto-regressive evaluation. Switching from the non-causal top-$k$
routing scheme in training to a causal predictor-based approach during auto-
regressive sampling leads to minimal performance degradation. This is perhaps
due to the ease of learning this prediction problem, which is upwards of 97%
accurate soon into training.

We evaluated MoD variants during auto-regressive sampling (see figure 6). Each
model was tested on exactly the same held-out data comprising $256000$
sequences ($500$M tokens). When switching from the top-$k$ routing method to
the predictor-based routing method we observed little performance degradation.
As in the training setting, there exist MoD variants that are better
performing than the isoFLOP-optimal baseline, while requiring fewer FLOPs per
forward pass. These results suggest that the compute savings offered by MoD
transformers should translate beyond the training setting.

###  4.3 Mixture-of-Depths-and-Experts (MoDE)

Figure 7: Mixture-of-Depths-and-Experts (MoDE). The MoD technique can be
implemented alongside MoE (together comprising MoDE models) in two
straightforward manners: staged, which first implements MoD machinery prior to
MoE machinery, and integrated, which uses one routing operation to funnel
tokens to either experts or no-op operations.

The MoD technique can be naturally integrated with MoE models (together
comprising MoDE models) in addition to vanilla transformers. In figure 7 we
present results showing that the performance improvments offered by MoD
compound with those of MoE. We tried two variants: in staged MoDE, which
routes tokens around or towards blocks prior to the self-attention step, and
integrated MoDE, which implements MoD routing by integrating “no-op” experts
among the conventional MLP experts. The former is advantageous because it
allows for tokens to skip the self-attention step, while the latter is
advantageous because it simplifies the routing machinery. We noticed that
implementing MoDE in the integrated manner was distinctly better than simply
reducing the capacity of experts in conventional MoE models, and relying on
token dropping to implement residual routing. We believe this is because with
the integrated MoDE machinery, tokens explicitly learn to choose the residual
path around the experts, as opposed to preferring an expert but being dropped
when implemented as a capacity reduction.

##  5 Discussion

Mixture-of-Depths transformers empirically demonstrate that one can improve on
isoFLOP-optimal baseline performance with models that use fewer FLOPs per
forward pass. This means that—for a given training FLOP budget—we can train
models that are both faster and better performing than their baseline
counterparts. Previously, to train models that are both faster and as- or
better-performing than isoFLOP-optimal models, one would have to use surplus
compute to _overtrain_ smaller models (notably, this overtraining technique is
still possible with MoD transformers, and speed gains should compound).

While MoD transformers require fewer FLOPs per forward pass, one cannot forego
FLOPs indiscriminately. Rather, it is crucial to use learned routing
decisions—much like in Mixture-of-Experts transformers—to determine whether a
token should participate in self-attention and the subsequent MLP (requiring
FLOPs), or not (saving FLOPs).We can then use any saved FLOPs by, for example,
making the model bigger or training it for longer. Our results show that
indeed FLOPs may be inefficiently used in vanilla transformer models, and that
there may be more efficient ways for them to be expended.

Learned routing mechanisms are sometimes non-causal; that is, information
about the future is used to determine a given token’s routing decision. This
is generally true for top-k routing mechanisms, which are useful because they
forego the need for auxiliary balancing losses. However, top-k routing
mechanisms present difficulties in post-training autoregressive sampling,
where it is impossible to use information about future token identities to
determine routing decisions. In this work we show that one can successfully
use a top-k routing scheme during training, but not require it during later
autoregressive sampling. Eiher a simple auxiliary classifier, or auxiliary
loss on the router, is sufficient to learn the top-$k$ routing decisions such
that it can mimic the top-$k$ decisions during autoregressive sampling, with
minimal to no performance degradation.

Intuitively, a token might learn to route around blocks because the prediction
being made at that step is easier, and hence, does not require as much
compute. However, this strategy is undoubtedly not all that the network
learns. If a token does not participate in self-attention at a certain block,
then later tokens will also not be able to attend to it. Thus, whether tokens
decide to route or not impacts both the current step’s prediction and future
predictions via causal self-attention, and how the network balances these
effects is guided by their influence on the overall language modeling
objective.

This insight opens the door to MoD variants that decouple the routing for
queries, keys and values. For example, perhaps a token would prefer to be
among the queries, but not the keys, for a given self-attention computation.
One can imagine extending this idea even further into the domain of "long-term
memory": perhaps there are tokens that would be extremely valuable as keys,
regardless of whether it is useful for them to also be among the queries at
the step of their occurrence. Learned routing could be a powerful mechanism
for deciding which tokens these might be, perhaps funnelling them into a long-
term memory buffer that is available during future self-attention. One
advantage of such an approach to long-term memory is that tokens decide once,
at the moment of "memory encoding", whether they should be retrieved in the
future. This is more computationally efficient than performing a full content-
based lookup across an entire memory buffer for each step in the future, and
could be one step towards drastically increasing the context-length available
for making a prediction.

Unlike MoE transformers that route between effectively the same computation
(usually MLPs), MoD transformers demonstrate the value of routing among
different types of computations. In this work the types were either the
conventional transformer block, or a null computation (functionally equivalent
to multiplying by zero). However, one can imagine extending this idea further
by routing between even more types of computation. For example, perhaps some
tokens are routed to "memory lookup" functions, and others are routed to "tool
use" functions. In general, the routing machinery we deployed provides a knob
for adjusting the types of computations available to the network and their
relative cost (in total FLOPs); if one wants to introduce an expensive
computation, then this can be offset by setting its capacity to some small
amount, and hence, by routing only a small number of tokens to it.

Altogether, MoD transformers are another tool one can use to tune a model’s
compute per forward pass (and hence inference time). The machinery used to
implement MoD is also generic, and opens the doors to many extensions and
integration with other techniques, such as MoE.

## References

  * Ainslie et al. (2023) J. Ainslie, T. Lei, M. de Jong, S. Ontañón, S. Brahma, Y. Zemlyanskiy, D. Uthus, M. Guo, J. Lee-Thorp, Y. Tay, Y.-H. Sung, and S. Sanghai.  Colt5: Faster long-range transformers with conditional computation, 2023\.
  * Bapna et al. (2020) A. Bapna, N. Arivazhagan, and O. Firat.  Controlling computation versus quality for neural sequence models.  _CoRR_ , abs/2002.07106, 2020.  URL https://arxiv.org/abs/2002.07106.
  * Bengio et al. (2016) E. Bengio, P.-L. Bacon, J. Pineau, and D. Precup.  Conditional computation in neural networks for faster models, 2016.
  * Bengio (2013) Y. Bengio.  Deep learning of representations: Looking forward, 2013.
  * Bengio et al. (2013) Y. Bengio, N. Léonard, and A. Courville.  Estimating or propagating gradients through stochastic neurons for conditional computation, 2013.
  * Bolya et al. (2023) D. Bolya, C.-Y. Fu, X. Dai, P. Zhang, C. Feichtenhofer, and J. Hoffman.  Token merging: Your vit but faster, 2023.
  * Cho and Bengio (2014) K. Cho and Y. Bengio.  Exponentially increasing the capacity-to-computation ratio for conditional computation in deep learning, 2014.
  * Dehghani et al. (2018) M. Dehghani, S. Gouws, O. Vinyals, J. Uszkoreit, and Ł. Kaiser.  Universal transformers.  _arXiv preprint arXiv:1807.03819_ , 2018.
  * Elbayad et al. (2019) M. Elbayad, J. Gu, E. Grave, and M. Auli.  Depth-adaptive transformer.  _CoRR_ , abs/1910.10073, 2019.  URL http://arxiv.org/abs/1910.10073.
  * Fedus et al. (2022) W. Fedus, B. Zoph, and N. Shazeer.  Switch transformers: Scaling to trillion parameter models with simple and efficient sparsity, 2022.
  * Graves (2016) A. Graves.  Adaptive computation time for recurrent neural networks.  _CoRR_ , abs/1603.08983, 2016.  URL http://arxiv.org/abs/1603.08983.
  * Guo et al. (2022) M. Guo, J. Ainslie, D. Uthus, S. Ontanon, J. Ni, Y.-H. Sung, and Y. Yang.  Longt5: Efficient text-to-text transformer for long sequences, 2022.
  * Gupta and Agrawal (2021) M. Gupta and P. Agrawal.  Compression of deep learning models for text: A survey, 2021.
  * He et al. (2021) J. He, C. Zhou, X. Ma, T. Berg-Kirkpatrick, and G. Neubig.  Towards a unified view of parameter-efficient transfer learning.  _arXiv preprint arXiv:2110.04366_ , 2021.
  * Jernite et al. (2017) Y. Jernite, E. Grave, A. Joulin, and T. Mikolov.  Variable computation in recurrent neural networks, 2017.
  * Lei et al. (2023) T. Lei, J. Bai, S. Brahma, J. Ainslie, K. Lee, Y. Zhou, N. Du, V. Y. Zhao, Y. Wu, B. Li, Y. Zhang, and M.-W. Chang.  Conditional adapters: Parameter-efficient transfer learning with fast inference, 2023.
  * Lepikhin et al. (2020) D. Lepikhin, H. Lee, Y. Xu, D. Chen, O. Firat, Y. Huang, M. Krikun, N. Shazeer, and Z. Chen.  Gshard: Scaling giant models with conditional computation and automatic sharding.  _arXiv preprint arXiv:2006.16668_ , 2020.
  * Liu et al. (2021) Z. Liu, Z. Xu, H.-J. Wang, T. Darrell, and E. Shelhamer.  Anytime dense prediction with confidence adaptivity.  _arXiv preprint arXiv:2104.00749_ , 2021.
  * Schuster et al. (2022) T. Schuster, A. Fisch, J. Gupta, M. Dehghani, D. Bahri, V. Q. Tran, Y. Tay, and D. Metzler.  Confident adaptive language modeling, 2022.
  * Shazeer et al. (2017) N. Shazeer, A. Mirhoseini, K. Maziarz, A. Davis, Q. Le, G. Hinton, and J. Dean.  Outrageously large neural networks: The sparsely-gated mixture-of-experts layer.  _arXiv preprint arXiv:1701.06538_ , 2017.
  * Simoulin and Crabbé (2021) A. Simoulin and B. Crabbé.  How many layers and why? An analysis of the model depth in transformers.  In _Proceedings of the 59th Annual Meeting of the Association for Computational Linguistics and the 11th International Joint Conference on Natural Language Processing: Student Research Workshop_ , pages 221–228, Online, Aug. 2021. Association for Computational Linguistics.  10.18653/v1/2021.acl-srw.23.  URL https://aclanthology.org/2021.acl-srw.23.
  * Tay et al. (2020) Y. Tay, M. Dehghani, D. Bahri, and D. Metzler.  Efficient transformers: A survey.  _CoRR_ , abs/2009.06732, 2020.  URL https://arxiv.org/abs/2009.06732.
  * Wang et al. (2017) X. Wang, F. Yu, Z. Dou, and J. E. Gonzalez.  Skipnet: Learning dynamic routing in convolutional networks.  _CoRR_ , abs/1711.09485, 2017.  URL http://arxiv.org/abs/1711.09485.
  * Zoph et al. (2022) B. Zoph, I. Bello, S. Kumar, N. Du, Y. Huang, J. Dean, N. Shazeer, and W. Fedus.  St-moe: Designing stable and transferable sparse expert models, 2022.


# A Deep Learning Approach for Active Anomaly Detection of Extragalactic
Transients

V. Ashley Villar Simons Junior Fellow, Department of Astronomy, Columbia University, New York, NY 10027, USA  Miles Cranmer Princeton University, Princeton, NJ 08544, USA  Edo Berger Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138-1516, USA  The NSF AI Institute for Artificial Intelligence and Fundamental Interactions  Gabriella Contardo Flatiron Institute, New York City, NY 10010, USA  Shirley Ho  Flatiron Institute, New York City, NY 10010, USA  Griffin Hosseinzadeh Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138-1516, USA  The NSF AI Institute for Artificial Intelligence and Fundamental Interactions  Joshua Yao-Yu Lin University of Illinois at Urbana-Champaign, Urbana, IL 61820, USA

###### Abstract

There is a shortage of multi-wavelength and spectroscopic followup
capabilities given the number of transient and variable astrophysical events
discovered through wide-field, optical surveys such as the upcoming Vera C.
Rubin Observatory. From the haystack of potential science targets, astronomers
must allocate scarce resources to study a selection of needles in real time.
Here we present a variational recurrent autoencoder neural network to encode
simulated Rubin Observatory extragalactic transient events using 1% of the
PLAsTiCC dataset to train the autoencoder. Our unsupervised method uniquely
works with unlabeled, real time, multivariate and aperiodic data. We rank
1,129,184 events based on an anomaly score estimated using an isolation
forest. We find that our pipeline successfully ranks rarer classes of
transients as more anomalous. Using simple cuts in anomaly score and
uncertainty, we identify a pure ($\approx 95$% pure) sample of rare transients
(i.e., transients other than Type Ia, Type II and Type Ibc supernovae)
including superluminous and pair-instability supernovae. Finally, our
algorithm is able to identify these transients as anomalous well before peak,
enabling real-time follow up studies in the era of the Rubin Observatory.

editorials, notices — miscellaneous — catalogs — surveys

††journal: ApJ

##  1 Introduction

Wide-field, optical surveys such as the Asteroid Terrestrial-impact Last Alert
System (ATLAS; Jedicke et al. 2012), the All-Sky Automated Survey for
SuperNovae (ASAS-SN; Shappee et al. 2014), the Panoramic Survey Telescope and
Rapid Response System 1 (Pan-STARRS1; Chambers et al. 2016) and the Zwicky
Transient Facility (ZTF; Bellm et al. 2018) have exponentially increased the
discovery rate of new transient events that vary on day to year timescales.
The upcoming Vera C. Rubin Observatory (Ivezić et al., 2019) and its decade-
long Legacy Survey of Space and Time (LSST) will greatly accelerate this
discovery rate to millions of new transients annually. However, a limited
fraction (likely $\lesssim 0.1$%) of all events can be followed up with
dedicated spectroscopic and multi-wavelength campaigns. Identifying transients
worthwhile of followup will be akin to finding needles in a haystack. Adding
to the challenge, we will need to identify such events quickly to capture
events pre- or near-peak to fully optimize the efficiency of follow-up
campaigns.

Over the past few years there have been several initial efforts aimed at
photometrically classifying transients to build pure samples of known
transient classes (Boone, 2019; Muthukrishna et al., 2019; Pasquet et al.,
2019; Hosseinzadeh et al., 2020; Sánchez-Sáez et al., 2020; Villar et al.,
2020). However, even the rarest transients known today, like superluminous
supernovae (SLSNe) and tidal disruption events (TDEs), will be discovered by
the thousands in the era of LSST (e.g., Villar et al. 2018; Bricman & Gomboc
2020). Detection and classification algorithms sensitive to anomalous
transients are essential in order to discover unexpected and even rarer
phenomena.

There is a growing literature on anomaly detection for astronomy applications.
For supernova (SN) light curves, Ishida et al. (2019); Pruzhinskaya et al.
(2019); Aleo et al. (2020); Martínez-Galarza et al. (2020) used isolation
forests and active anomaly discovery on archival datasets. Convolutional
autoencoders have recently been used to search for anomalies and glitches in
gravitational wave time series (Morawski et al., 2021). In the broader machine
learning literature, there has been increasing interest in active anomaly
detection (Chalapathy & Chawla, 2019). Typically these works focus on long and
well-sampled, single channel time series with anomalous periods of activities
(e.g., Zhang et al. 2019; Li et al. 2019), although some recent work has
focused on multivariate series (Zhao et al., 2020). There has been limited
focus on anomaly detection in irregularly sampled, aperiodic, multivariate
time series.

In this paper, we focus on out-of-distribution anomalies, which appear
distinct from all other known transients in some feature space. By taking a
completely data-driven approach to anomaly searches, our algorithm is agnostic
to physics and therefore sensitive to entirely unexpected phenomena. Our
anomaly detection pipeline is based on a variationl recurrent autoencoder
neural network (VRAENN) with no physical priors, and we search the learned
encoded space for out-of-distribution events. The paper is organized as
follows. In §2 we review the PLAsTiCC dataset used for training and anomaly
detection, and the breakdown of SN-like transients used in this study. In §3
we present our anomaly detection pipeline and the VRAENN architecture. We
discuss our results in §4 and conclude in §5.

Figure 1: Sample $grizY$ light curves of majority SN classes (Type Ia, Ibc,
II). Bold lines represent the 2D Gaussian Process mean function, with shaded
regions representing 68% confidence intervals. Even when entirely missing one
or more bands, our method is able to produce reasonable interpolated light
curves. Note that the y-axis, the magnitudes used to train the Gaussian
Process, is designed such that the light curve tends to zero as the light
curve approaches the survey magnitude limit.

##  2 Dataset and Preprocessing

In this study we use PLAsTiCC (Photometric LSST Astronomical Time-Series
Classification Challenge), a simulation of three years of Rubin Observatory
data that includes over 3.5 million transient events from eighteen unique
physical classes, extending to a redshift of $z\approx 1.5$ (Allam Jr et al.,
2018; Kessler et al., 2019). Each event is observed across six broadband
filters ($ugrizY$) following the LSST observing strategy at the time the
simulation was produced. Sample light curves are shown in Figure 1). Along
with the light curves, PLAsTiCC provides metadata including the redshift,
Milky Way reddening, physical parameters used to generate the model, and a
realistic photometric redshift estimate (see Kessler et al. 2019 for details).
For this observing strategy, each event is observed every few days (in any
filter) and roughly once a week in the same filter.

The PLAsTiCC dataset was originally created for a public data science (Kaggle)
competition111https://www.kaggle.com/c/PLAsTiCC-2018 to classify transients.
We re-purpose this dataset as a training set for anomaly detection in an LSST-
like data stream. Here, anomalous events will be determined by the metadata
(i.e., if the event comes from a rare astrophysical origin). We remove classes
observable only within the Milky Way (e.g., variable stars) and SNe with fewer
than three detections in any filter within 300 days of peak brightness. This
cut is not necessary, as our algorithm can take light curves of any length;
however, light curves with fewer points are very unlikely to be selected for
detailed follow-up in reality. Additionally, we only utilize transients from
the Wide-Fast-Deep (WFD) survey and remove events within the Deep Drilling
Fields. In total, our data set contains 1,129,184 extragalactic light curves
from thirteen classes:

  * •

Normal Type Ia SNe arise from the thermonuclear explosions of carbon-oxygen
white dwarfs. The models were generated using the standard SALT-II light curve
models (Guy et al., 2007), conditioned on $\approx 500$ light curves from the
Joint Lightcurve Analysis (Betoule et al., 2014). Type Ia SNe represent 54.2
percent of our dataset. We consider Type Ia SNe part of the majority classes.

  * •

Type II SNe222 In the original version, PLAsTiCC grouped normal Type II SNe
and Type IIn SNe into a single class. We separate these classes due to their
distinct physical origins and unique light curves. are the explosions of
massive stars that have retained their hydrogen envelopes. They are often
characterized by long plateaus in their optical light curves. The models were
generated from SED templates (Kessler et al., 2010; Anderson et al., 2014;
Galbany et al., 2016; Sako et al., 2018). Type II SNe make up 25.8 percent of
our dataset. We consider Type II SNe part of the majority classes.

  * •

Type Ibc SNe are core-collapse SNe of stars with stripped hydrogen (Ib) and
helium (Ic) envelopes. The models were generated using a combination of MOSFiT
(Villar et al., 2017; Guillochon et al., 2018) and SED templates (Kessler et
al., 2010). Type Ibc SNe make up 5.8 percent of our sample. We consider Type
Ibc SNe part of the majority classes.

  * •

Type I Superluminous SNe (SLSNe) are luminous, hydrogen-free events thought to
be powered by rapidly spinning, highly magnetized neutron stars. SLSN models
were produced using MOSFiT (Guillochon et al., 2018; Nicholl et al., 2017;
Villar et al., 2018). They make up 2.2 percent of our sample. We consider
SLSNe to be members of the minority classes.

  * •

Type Iax are irregular Type Ia SNe with typically lower luminosities and lower
velocities compared to normal Type Ia SNe (Li et al., 2003). The models were
generated using available data in the Open Supernova Catalog (Guillochon et
al., 2017). Type Iax SNe make up 1.8 percent of our dataset. We consider Type
Iax SNe to be members of the minority classes.

  * •

Type IIn SNe2 are core-collapse SNe mainly powered by the interaction of the
SN ejecta with circumstellar material (CSM). The models were generated using
MOSFiT (Villar et al., 2017; Guillochon et al., 2018; Jiang et al., 2020).
Type IIn SNe make up 1.8 percent of our sample. We consider Type IIn SNe to be
members of the minority classes.

  * •

Type Ia-91bg are fainter, faster and redder Type Ia SNe that make up $\approx
20$% of the volumetric Type Ia sample (Filippenko et al., 1992; Graur et al.,
2017), and $\approx 3$% of the observational sample (Li et al., 2011). The
model light curves are based on the spectral energy distribution (SED)
templates from Nugent et al. (2002). Type Ia-91bg SNe make up 1.2 percent of
our sample. We consider Type Ia-91bg SNe to be members of the minority
classes.

  * •

Tidal Disruption Events (TDEs) result from the tidal disruption of stars by
supermassive black holes (SMBH; Rees 1988). TDE models were generated using
MOSFiT (Guillochon et al., 2018; Mockler et al., 2019) and make up 0.6 percent
of our sample. We consider TDEs to be members of the minority classes.

  * •

Ca-rich Transients (CARTs) are intermediate luminosity transients whose
spectra appear rich in calcium (Kasliwal et al., 2012). CARTs are modeled
using MOSFiT, assuming they are powered by the radioactive decay of 56Ni. We
note that this is the same model used to generate Type Ibc SNe, but with a
distinct parameter space. CARTs make up 0.31 percent of our sample. We
consider CARTs to be members of the minority classes.

  * •

Intermediate Luminosity Optical Transients (ILOTs) are transients that are
brighter than novae but less luminous than SNe (Kasliwal, 2012). In this case,
we assume that ILOTs arise from CSM interaction with low-energy eruptions (or
explosions) of massive stars. ILOTs have been modelled using MOSFiT
(Guillochon et al., 2017; Villar et al., 2017; Jiang et al., 2020) and
represent 0.08 percent of our sample. We consider ILOTs to be members of the
minority classes.

  * •

Pair Instability Supernovae (PISNe) are the explosions of low-metallicity
massive stars ($M_{\mathrm{ZAMS}}\sim 130-260$ M⊙) that reach core
temperatures high enough to form electron-positron pairs (Kasen et al., 2011).
Compared to normal core-collapse SNe, PISNe have high kinetic energies and
larger ejecta masses. PISNe are modelled using MOSFiT assuming that they are
powered by the radioactive decay of 56Ni. PISNe make up 0.07 percent of our
sample. We consider PISNe to be members of the minority classes.

  * •

Kilonovae (KNe) arise from the formation of radioactive rapid neutron capture
elements in binary neutron star (and potentially neutron star - black hole)
mergers. The models are based on theoretical calculations (Kasen et al.,
2017). There are only two KNe in our sample. KNe are dim and short-lived,
making them nearly impossible for LSST to discover in the WFD strategy
explored here. We consider KNe to be members of the minority classes.

  * •

Active Galactic Nuclei (AGN) refers generally to galaxies with active SMBHs
from accreting gas. AGN have a wide range of observed behavior, although
typically they vary on timescales of weeks to years at the $\lesssim 10$%
level. AGN variability is modeled using a damped random walk as described in
MacLeod et al. (2011). AGN make up 6.3 percent of our dataset. We note that
AGN are the only class that is likely not representative of the true LSST
datastream; AGN will be much more numerous, with likely millions in the
complete sample. Many bright AGN (which are those represented in PLAsTiCC)
will be identified within the first year of LSST; however, transient bright
flares will be of interest to the community (Graham et al., 2017). Because AGN
are distinct from SN-like transients, we consider their effects separately in
an Appendix.

Figure 2: Breakdown of the various transient classes used in this study. Note
that AGN are likely highly underrepresented compared to the true breakdown in
the LSST datastream; however, AGN are exceptional in their light curve
properties and may be identified early in the survey. Parenthetical numbers
associated with each class represent the percentage breakdown. There are a
total of two kilonovae (KNe) in our sample.

The breakdown of the various classes is shown in Figure 2. The observational
rates of each class are a combination of the volumetric rates (intrinsic
rarity) and observational effects (e.g., luminosity function, duration). The
breakdown used here is designed to match what is empirically expected from a
wide-field survey such as LSST (see e.g., Perley et al. 2020; Villar et al.
2020). We note that this LSST cadence simulation detects just two KNe that
pass our cuts in the simulated observations, highlighting the need for target-
of-opportunity observations to better capture these rare events.

We define Type Ia SNe, Type II SNe and Type Ibc SNe together as the majority
classes because together they make up the bulk of the dataset ($\approx
86\%$). We consider all other classes (excluding AGN) to be minority classes,
with fraction of $\lesssim 2$% of the sample.

We pre-process the data as follows. Following Villar et al. (2020), we scale
the light curve magnitudes such that zero corresponds to the magnitude limit.
This is to aid the Gaussian process interpolation (which will tend towards
zero before and after the supernova). We correct the light curve of each event
for time dilation and convert to absolute magnitudes based on the provided
photo-$z$ values (see Villar et al. 2020 for details). The photo-$z$ values
are based on host galaxy association, and use a color-matched nearest neighbor
method presented in Graham et al. (2018). This method is trained on a
realistic sample of galaxies that would have spectroscopic redshifts
available. About 17% of the redshift estimates are outliers, defined as
$|z_{\mathrm{true}}-z_{\mathrm{phot}}|/(1+z_{\mathrm{phot}})>3\sigma_{\mathrm{IQR}}$
by Kessler et al. 2019, where $\sigma_{\mathrm{IQR}}$ is the typical error for
galaxies near $z_{\mathrm{true}}$. We make no cuts on redshift uncertainty,
and instead we account for redshift uncertainties via a simple Monte Carlo
method discussed in §3.

We additionally correct the light curves for Galactic reddening, as outlined
in Villar et al. (2020). Finally, we temporally shift each event such that the
observed time of peak brightness (in any filter) is considered $t=0$. This
means that as new data are taken during an event’s rise, the time of peak
luminosity will continue to shift until the true peak luminosity has been
observed or the event dims. For example, if a transient first peaks in
$g$-band 10 days post explosion and then reaches a brighter luminosity in
$r$-band at 15 days post explosion, the peak used by the neural network will
be at 15 days post explosion; this is independent of the bolometric peak
luminosity. We find that this prescription of phase (versus, for example, time
since first detection) leads to better performance in the autoencoder.

###  2.1 Interpolation using Gaussian Processes

The PLAsTiCC light curves are irregularly sampled across time and filters,
with no more than one filter observed at any given time. For the architecture
discussed below, we require a flux and error estimate for each filter at every
observation time. To produce this information, we use a 2D Gaussian Process
(GP) to interpolate the light curve over time and filter, with a multivariate
Gaussian kernel described by:

$\kappa(t_{i},t_{j},f_{i},f_{j};\sigma,l_{t}l_{f})=\\\
\sigma^{2}\times\exp\Big{[}-\frac{(t_{i}-t_{j})^{2}}{2l_{t}^{2}}\Big{]}\times\exp\Big{[}-\frac{d(f_{i},f_{j})^{2}}{2l_{f}^{2}}\Big{]}$
(1)

where $f$ represents the six ($ugrizY$) filters; $l_{t}$ and $l_{f}$ are
characteristic correlation length scales in time and filter, respectively; and
$d(f_{i},f_{j})$ is the Wasserstein-1 distance between each filter’s
normalized throughput which we optimize to each specific light curve. This
choice in the distance metric loosely measures the similarity between two
distributions. Mathematically, we treat each filter as a density function in
wavelength. This distance metric is minimized when filters overlap and
simplifies to difference between central wavelengths in the limit of
infinitely narrow passbands.

The GP interpolation both accounts for and produces error estimates for each
flux measurement. Flux uncertainties are utilized in both the encoding method
and anomaly scoring steps of our algorithm (see §3) in order to make our
algorithm robust to low-confidence outliers. Finally, we note that during
testing we implemented a similar (though less physically motivated)
interpolation scheme in Villar et al. (2020) and found that the interpolation
methods led to visually similar light curves.

We implement the GP pre-processing using sklearn, optimizing each light curve
independently via the minimize function from scipy, which uses the
Broyden–Fletcher–Goldfarb–Shanno optimization algorithm (Fletcher, 1986).
Sample light curves are shown in Figure 1. Even in cases of poorly sampled
light curves, or light curves in which a band is completely missing, the GP
produces reasonable flux and error estimates across all filters.

We note that our GP utilizes the complete light curve for interpolation. In
reality, only the light curve up to the most recent observation will be
available in real time. One may be concerned that because our GP has been
conditioned on the entire light curve, it has more information than what will
be available for real-time usage. This likely has no effect on our results, as
each observation heavily anchors the GP prediction (see Figure 1) and the
learned GP kernel sizes are similar to our priors (see e.g., Villar et al.
2020 for typical SN values).

##  3 VRAENN: Architecture, Training and Anomaly Detection

Figure 3: Summary of our anomaly detection pipeline. (1) We interpolate each
$ugrizY$ light curve using a 2D Gaussian process. (2) We train the variational
recurrent autoencoder in an unsupervised manner. Light curves are represented
as a time series, in which each epoch is represented by 11 features: 6 flux
values, 6 estimated error values, and one time. This time series (consisting
of $N$ points) is encoded. The encoding layer is repeated $N$ times, each time
appended with the associated time value. We use this network to encode each
light curve. (3) We use an isolation forest to assign an anomaly score for
each light curve using the encoded vectors. The isolation forest is
represented here in one subspace of the encoded space, with anomalous events
highlighted as purple stars.

Following pre-processing of the training set, our anomaly detection pipeline
includes two steps. First, we learn an encoded form of each light curve by
training a variational recurrent autoencoder neural network (VRAENN) on the
full dataset. Then, we use an isolation forest to rank each light curve’s
encoded form by an anomaly score. We utilize a simple Monte Carlo to estimate
the uncertainty on this score. In this section, we describe the VRAENN
architecture, the training process, the isolation forest and our error
estimation method. Our code, along with all chosen hyperparameters of our
model, is available via Github333https://github.com/villrv/vraenn

Our novel VRAENN architecture is well-suited to the problem of searching for
unknown, anomalous transients. While more traditional feature extraction via
model-fitting (e.g., Villar et al. 2019; Hosseinzadeh et al. 2020) or using
pre-defined quantities (Boone, 2019) has been successful for classification,
dozens of features are required. We want to avoid searching for anomalous
events in high-dimensional data, in which distance metrics are more
challenging to meaningfully define (Liu et al., 2012). Additionally, the
VRAENN architecture is built without any physical models, meaning that it is
sensitive to new and unexpected physical processes, which are observationally
distinct from known transients.

The architecture specifically used here has three additional benefits: (i) The
ability to handle unevenly sampled light curves across filters and time using
recurrent neurons; (ii) the ability to produce extrapolated and interpolated
light curves; and (iii) an insensitivity to noisy data and spurious outliers
due to the variational architecture and error estimation.

The VRAENN architecture, based on the model presented in Villar et al. (2020),
uses recurrent neurons to read in the GP light curve and estimated errors, and
encodes this light curve as a vector. This is achieved by encoding the light
curve into a series of smaller matrices until the information reaches a small
bottleneck layer of size $1\times 10$. This is an unoptimized choice of size,
although Villar et al. (2020) found a similar result after a hyperparameter
search of a similar architecture. This layer is known as the encoded layer,
and the layers preceding it are known as the “encoder”. Each of these layers
uses gated recurrent unit (GRU) neurons with a combination of hyperbolic
tangent, sigmoid, relu and linear activation functions. The GRU is memory-
efficient version of the Long short-term memory, the standard for RNNs (Cho et
al., 2014).

Rather than using a static vector as a bottlebeck, our VRAENN learns a
distribution of encoding vectors. In particular, the encoding layer consists
of two vectors: one that represents the mean of a multivariate Gaussian (MVG)
and one that represents its diagonal covariance matrix. For each SN passed
into the VRAENN, we randomly select an encoding from the MVG defined by this
learned mean and variance. The “variational” aspect of our architecture refers
to this process of learning a distribution of encodings rather than a singular
encoding. This is helpful in generating a smooth encoding space, in which most
events will cluster, allowing us to more easily pick out anomalous events.

Figure 4: Scatter plot of representative encoding space for various classes.
Events with high anomaly scores are shown as stars. Most events cluster near
the origin, while more anomalous classes are seen as clouds outside the main
distribution. The anomalous Type Ia, Ibc and II SNe are those with incorrect
photo-$z$ estimates.

Before being passed into the second half of the autoencoder (the “decoder”),
the encoded layer is repeated $N$ times, with each time appended as a phase
relative to maximum light. This can be thought of as evaluating a SN model at
specific times, where the model is specified by 10 free parameters, with an
11th parameter specifying the time. The unique repeat layer of our
architecture allows us to evaluate the light curve at times not included in
the real data; i.e., interpolating and extrapolating the light curve if
desired (although this feature is not used in this work).

The decoder then produces the light curves at the $N$ times specified in the
repeat layers. As with the encoder, the decoder uses GRU neurons with
hyperbolic tangent activation functions. Figure 3 illustrates our full
pipeline, including a schematic of the neural network architecture used. The
VRAENN is optimized using a loss function combining the log of the weighed
mean squared error and the standard Kullback–Leibler (KL) divergence (which
measured how well our MVG represents our latent variable space):

$\mathcal{L}=\log\sum_{i=1}^{N}\frac{\Big{[}F_{i,\mathrm{True}}(t,f)-F_{i,\mathrm{Predicted}}(t,f)\Big{]}^{2}}{N}\\\
+\sum_{i=1}^{N}-0.5\Big{(}1+\log(\sigma_{i}^{2})-\mu_{i}^{2}-\sigma_{i}^{2}\Big{)}$
(2)

We minimize the loss function using the Adam optimizer (Kingma & Ba, 2014)
with standard learning parameters, $\alpha=10^{-4}$, $\beta_{1}=0.9$,
$\beta_{2}=0.999$ using Keras (Chollet et al., 2015) with a TensorFlow backend
(Abadi et al., 2016). We train our VRAENN on 1% of the sample (12,159 events),
reflecting (for example) the small dataset that will be available within
months of LSST coming online.

A sampled subspace of our encoding vectors is shown in Figure 4. The majority
of events cluster near zero, with anomalous events (like SLSNe) forming a
cloud outside of the main distribution.

###  3.1 Scoring Anomalies with an Isolation Forest

Once our VRAENN is trained, we can encode any PLAsTiCC light curve, partial or
complete, as a $1\times 10$ vector. We then pass these encodings into an
isolation forest (Liu et al., 2012). The isolation forest works by generating
a series of decision trees over a random subset of attributes. Each tree
recursively splits the set. Out-of-distribution anomalous events will be
isolated with very few splits, while an average event will cluster with
similar events, even after many splits. The number of splits is inversely
related to an anomaly score. For the sake of interpretability, we then convert
this raw score to a percentile. We use sklearn to implement the isolation
forest, using 1000 base estimators.

We identify several sources of possible error in the anomaly score:

  * •

Flux uncertainty due to Poisson noise. This is provided by PLAsTiCC as a
standard deviation for each flux measurement.

  * •

Uncertainty in the flux estimates from the GP, which is also estimated by the
GP.

  * •

Photometric redshift error, reported as a standard deviation. This affects the
entirety of the light curve as both an overall multiplicative flux term and a
time dilation adjustment.

  * •

Model uncertainty from the neural network converting the light curve into an
encoded vector.

We account for the first three using a simple Monte Carlo method. For each
transient, we generate 10 light curves that have photo-$z$ and flux values
drawn from Gaussian distributions described by the reported mean and errors.
We do not account for noise from the neural network itself, which could be
accomplished via an ensemble of networks. We find that, in general, the error
estimated from this method is sufficient to eliminate anomalous events arising
from poor data quality and incorrect photo-$z$ estimates by making a cut on
the anomaly score uncertainty.

We show how the anomaly score uncertainty changes as a function of both time
and number of data points in Figure 5. We first examine how the error changes
as a function of number of observations. The error grows until the light curve
reaches $\approx 3$ observations and then decreases. The error is initially
small because the autoencoder, with limited information, returns the mean
encoding vector from the training set (i.e., only one number is needed to
encode a single observation). These light curves are not flagged as anomalous.
However, once the light curve reaches a sufficient number of data points, the
fractional uncertainty follows the expected trend of decreasing with
increasing number of observations. When translated into fractional error
versus time, the fractional error rises near peak (where most light curves
have few data points) and declines at late times.

##  4 Results and Discussion

Figure 5: Time evolution of the relative anomaly score error for a
representative subsample of our test set. In the upper panel, relative errors
are low before peak luminosity and subsequently rise; however, the lower panel
reveals that the relative error drops with increased number of data points.

Our anomaly classification algorithm generates a ranked list of anomalous
events given the full dataset. We are interested in which events the anomaly
detection pipeline ranks highly, but we are also interested in understanding
when and why an event becomes anomalous.

We first explore the anomaly scores of the full test set. In Figure 6, we show
the cumulative distribution function (CDF) for each class in our set. Even for
events that would not be labelled as anomalous by our pipeline, the anomaly
score distribution matches expectations. Type Ia-like SNe (including 91bg,
Iax) are the least anomalous on average. Type Ibc SNe and TDEs are largely
distributed evenly across scores, with a slight tail at the upper end. While
Type II SNe are, on average, more anomalous, only a small fraction have high
anomaly scores. CARTs and Type IIn SNe are also more anomalous on average.
PISNe, SLSNe, ILOTs and AGN most drastically cluster at the high end of
anomaly scores, indicating that they are the most likely to be classified as
anomalies with our algorithm. Only two KNe were in our sample; both had
moderately high anomaly scores around the 80th percentile. Reassuringly, the
events which prefer higher anomaly scores do not cluster in any obvious
section of observational phase space (e.g., luminosity or duration), implying
that our VRAENN has picked up on more fundamental features.

Figure 6:  CDF of anomaly score for various classes of transients. PISNe,
SLSNe, ILOTs, KNe and AGN are clustered towards the high end of anomaly score.
Figure 7: Top row: Breakdown of anomalous events, given the cuts listed in
each panels. The top number represents the anomaly percentile threshold, while
$\sigma_{\mathrm{A}}/\mathrm{A}$ refers to the anomaly score fractional error
threshold. The right-most panel shows optimized cuts on the anomaly score and
uncertainty to minimize contamination from majority classes. SLSNe
overwhelmingly make up the majority of our anomaly sample. Bottom row:
Relative abundance of anomaly samples shown in top rows. Each panel shows the
“abundance” (ratio between the number of “anomalous” events in each class and
the expected number based on input fraction. The minority classes are heavily
over-represented, while majority classes are underrepresented.

In practice, we will be limited to a small fraction ($\sim 0.1$%) of the LSST
transients for follow-up. We therefore investigate several thresholds on the
anomaly scores and uncertainties to search for anomalous events. The results
presented here are summarized as histograms of anomaly sample breakdown and
abundances in Figure 7. In short, our anomaly sample is pure, with $\lesssim
10$% contamination from majority (Type Ia, Ibc, and II SNe) classes. The
anomaly sample significantly over-represents minority classes.

Many of the non-anomalous transients (e.g., Type Ia SNe) with high anomaly
scores have fractional uncertainties on the anomaly scores of
$\sigma_{\mathrm{A}}/\mathrm{A}\gtrsim 0.1$ (due to large uncertainties on
photo-$z$ estimates). We first make a strict cut on the full sample, looking
only at events in the top 90, 95, and 99 percentile of anomaly scores, and
with $\sigma_{\mathrm{A}}/\mathrm{A}<0.01$. Within these cuts, we find that
SLSNe make up a majority of the remaining sample. For the 95 percentile
cutoff, the anomalous SLSNe that pass our cut make up $\approx 5$% of the
original sample. In contrast, just 0.03% of the original input sample for Type
Ia and Type II SNe remain in our 95 percentile anomaly score cutoff. In other
words, our 95 percentile cutoff removed 95% of the SLSN sample, but removed
99.97% of the Type Ia and Type II SNe samples (leading to SLSNe being over-
represented in the anomalous sample by a factor of about 4.5). Similarly,
PISNe make up a small fraction ($\approx 1$%) of our anomalous sample, but we
retain $\approx 5$% of the original sample in the 95 percentile cutoff. The
results are similar for higher percentiles, with an even larger bias towards
the more anomalous classes. This is highlighted as an abundance measurement in
Figure 7.

We then explore a higher cutoff in anomaly score fractional error, keeping
events with $\sigma_{\mathrm{A}}/\mathrm{A}<4\times 10^{-4}$ (a threshold
chosen by hand to maximize the ratio of rare events, shown in the right-most
column of Figure 7). Even without cuts on anomaly score itself, SLSNe dominate
the sample, followed by PISNe and Type Ia-like SNe.

Noting that stringent cuts in anomaly score uncertainty lead to pure samples
of minority classes, we investigate if incorrect photo-$z$ estimates are
responsible for false positive detection of anomalous events in majority
classes. We find that removing events with
$|z_{\mathrm{phot}}-z_{\mathrm{true}}|/z_{\mathrm{true}}>2$ does improve the
purity of anomalous samples for cuts with $A>99$% and $A>99.9$%, although has
a lesser effect for $A>95$% cut. For all cuts, we find that removing events
with incorrect photo-$z$ estimates drastically decrease the number of Type Ia
SNe which pass our anomaly thresholds. For example, in the $A>95$% sample, the
sample fraction of Type Ia SNe drops from $\approx 10$% to $\approx 3$%.
Reducing the number of outlier photo-$z$ estimates in LSST will likely
substantially improve the purity of our anomaly sample.

We investigate why some Type Ia SNe are within the most anomalous events even
with stringent cuts on anomaly score and uncertainty. All of the events which
pass our 99% cutoff for anomaly score have catastrophically incorrect
photo-$z$ values. Specifically, we find that these events are typically
injected at relatively low redshifts ($z_{\mathrm{true}}\lesssim 0.05$) yet
have reported photometric redshifts $z_{\mathrm{phot}}\approx 2$ with small
reported uncertainties—again highlighting the need for reliable photo-$z$
estimates.

Our anomaly detection algorithm is biased towards bright events, which begs
the question: Is our algorithm making a trivial cutoff on luminosity to search
for anomalies? We test this hypothesis by making such a cut. We rank each
event by peak luminosity (in any filter). If we keep the top 1% of brightest
events, 50% of the sample is made up of Type Ia SNe and only $\approx 4$% are
SLSNe. We can further make a cut on photometric redshift uncertainties. Even
with extremely aggressive cuts ($\sigma_{z}/z<0.05$) which remove 99.99% of
the sample, only $\approx 10$% of the sample are SLSNe. The plurality,
$\approx 40$%, are Type Ia SNe. Our detection algorithm, although biased
towards brighter events, is sensitive to all outliers regardless of
luminosity.

###  4.1 Anomaly Detection in Live Streaming Data

Our analysis thus far has focused on the full light curves, rather than real
time followup. We next turn to how the anomaly scores evolve over time for a
representative subsample of $\approx 10^{3}$ SNe. As shown in Figure 5, most
events are identified near peak luminosity ($t\approx 0$). We focus on how
anomaly scores vary for the majority classes versus the minority classes.
Events from the minority classes are much more likely to be triggered as
anomalous before peak. In our representative sample, $\approx 60$% of Type
Ia/Ibc/II SNe identified as anomalous are first marked as such after peak
magnitude. Type Ia and Type Ibc SNe are typically flagged just around peak,
while Type II SNe are flagged, on average, about a week post-peak. In
contrast, $\approx 65$% of anomalies from minority classes are identified as
anomalous before peak. They are flagged, on average, about one week before
peak. However, short-lived anomalous transients (such as ILOTs and CaRTs) are
flagged around or post-peak. We visualize our findings for a representative
sample of transients in Figure 8. The background histogram of the vector plot
in Figure 8 shows the overall density of the anomaly scores over time. Most
curves have just a few points around $t=0$ and have low anomaly scores. The
arrows show that the anomaly scores typically rise before peak but plateau
after peak. This implies that, on average, the scores are steady post-peak.

Finally, we investigate if light curves often begin as anomalous and then drop
to less anomalous over time. For this test, we use a 99th-percentile cutoff in
anomaly score. We show a selection of representative anomaly curves over time
in Figure 8; in this figure, grey curves are SNe which never reach the anomaly
threshold; blue curves are members of the minority classes which reach the
anomaly threshold; and orange curves are members of the majority classes which
reach the anomaly threshold. Light blue/orange curves represent members of the
minority/majority classes which reach the anomaly threshold but drop below
threshold before the end of the event. In contrast, dark curves are those
which remain anomalous until the end of the event. We find that about 7% of
events reach this threshold at least once (with, by definition, 1% remaining
anomalous by the final observation). As previously mentioned, we find that
members of the minority classes are typically flagged as anomalous before
peak, while members of the majority classes are flagged after peak; this is
regardless of whether or not those events drop below the anomaly threshold by
the final observation. The fact that the minority classes are flagged before
peak is useful in practice, as we are much more likely to follow events caught
before peak luminosity.

Figure 8: Evolution of anomaly scores as a function of time for a
representative set of SNe. Left: Vector plot showing the flow of anomaly
scores over time. Arrows represent the average gradients of anomaly scores.
Right: Anomaly score curves for a representative set of SNe. The green dotted
line is the 99th percentile threshold for final values (having an anomaly
score of $\mathrm{A}\approx 0.95$). Grey curves are SNe which never cross the
anomaly threshold. Blue and orange curves represent events from the minority
and majority classes, respectively. The light orange/blue curves drop below
the anomaly threshold before the final observation, while darker curves remain
above the anomaly threshold until the end. The colored vertical dashes above
the anomaly curves represents the times in which the anomalous events first
cross the anomaly threshold. Members of the minority class are much more
likely to cross this trigger threshold before maximum light, while erroneous
anomalies from the majority classes are triggered after peak luminosity.

##  5 Conclusions

We presented an anomaly detection pipeline for SN-like transients in an LSST-
like data stream. We re-purpose the PLAsTiCC dataset to train and test our
algorithm, allowing us to analyze the how, why, and when events are tagged as
anomalous. Our key results are as follows:

  * •

We present a novel variational recurrent autoencoder neural network
architecture that encodes SN-like light curves in real time into a low-
dimension vector. We train this neural network on 1% of 1,129,184 events from
the PLAsTiCC dataset.

  * •

We pair this neural network with an isolation forest to assign every transient
an anomaly score. We use a Monte Carlo method to estimate our uncertainty on
this score.

  * •

We examine the efficacy of our algorithm through a series of percentile and
uncertainty cuts. We find that our algorithm is successful in identifying
anomalous classes, especially luminous events.

  * •

We find that our algorithm is often limited by the photometric redshift
estimate. Catastrophically incorrect redshift estimates of Type Ia SNe are
especially challenging to remove from our anomaly samples.

  * •

We find that members of minority classes (i.e., SNe which are not Type Ia,
Type Ibc or Type II) are likely to be identified before peak luminosity. In
contrast, erroneously flagged members of the majority classes are more likely
to be flagged post-peak.

Much is left to be done to sufficiently prepare for the deluge of data which
will come in the new era of the Rubin Observatory. This algorithm must be
tested on real data with an active follow up campaign. Furthermore, it is
possible to use the anomaly score designed here in classification methods to
increase the purity of rare transient samples at the cost of completeness.

We thank Lehman Garrison for insightful conversations. VAV is supported by the
Simons Foundation through a Simons Junior Fellowship (#718240). The Berger
Time Domain group at Harvard is supported in part by NSF and NASA grants, as
well as by the NSF under Cooperative Agreement PHY-2019786 (The NSF AI
Institute for Artificial Intelligence and Fundamental Interactions
http://iaifi.org/). This research made use of the following software packages:
numpy (Harris et al., 2020), scipy (Virtanen et al., 2020), jupyter (Kluyver
et al., 2016), sklearn (Pedregosa et al., 2011), matplotlib (Hunter, 2007),
tensorflow (Abadi et al., 2016), and astropy (Astropy Collaboration et al.,
2018).

## References

  * Abadi et al. (2016) Abadi, M., Barham, P., Chen, J., et al. 2016, in 12th $\\{$USENIX$\\}$ Symposium on Operating Systems Design and Implementation ($\\{$OSDI$\\}$ 16), 265–283
  * Aleo et al. (2020) Aleo, P. D., Ishida, E. E., Kornilov, M., et al. 2020, Research Notes of the AAS, 4, 112
  * Allam Jr et al. (2018) Allam Jr, T., Bahmanyar, A., Biswas, R., et al. 2018, arXiv preprint arXiv:1810.00001
  * Anderson et al. (2014) Anderson, J. P., González-Gaitán, S., Hamuy, M., et al. 2014, The Astrophysical Journal, 786, 67
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Bellm et al. (2018) Bellm, E. C., Kulkarni, S. R., Graham, M. J., et al. 2018, Publications of the Astronomical Society of the Pacific, 131, 018002
  * Betoule et al. (2014) Betoule, M. e. a., Kessler, R., Guy, J., et al. 2014, Astronomy & Astrophysics, 568, A22
  * Boone (2019) Boone, K. 2019, The Astronomical Journal, 158, 257
  * Bricman & Gomboc (2020) Bricman, K., & Gomboc, A. 2020, The Astrophysical Journal, 890, 73
  * Chalapathy & Chawla (2019) Chalapathy, R., & Chawla, S. 2019, arXiv preprint arXiv:1901.03407
  * Chambers et al. (2016) Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv e-prints, arXiv:1612.05560.  https://arxiv.org/abs/1612.05560
  * Cho et al. (2014) Cho, K., van Merrienboer, B., Gulcehre, C., et al. 2014, arXiv e-prints, arXiv:1406.1078.  https://arxiv.org/abs/1406.1078
  * Chollet et al. (2015) Chollet, F., et al. 2015, Keras, GitHub.  https://github.com/fchollet/keras
  * Filippenko et al. (1992) Filippenko, A. V., Richmond, M. W., Branch, D., et al. 1992, The Astronomical Journal, 104, 1543
  * Fletcher (1986) Fletcher, R. 1986, Practical methods of optimization (Wiley)
  * Galbany et al. (2016) Galbany, L., Hamuy, M., Phillips, M. M., et al. 2016, The Astronomical Journal, 151, 33
  * Graham et al. (2017) Graham, M. J., Djorgovski, S., Drake, A. J., et al. 2017, Monthly Notices of the Royal Astronomical Society, 470, 4112
  * Graham et al. (2018) Graham, M. L., Connolly, A. J., Ivezić, Ž., et al. 2018, AJ, 155, 1, doi: 10.3847/1538-3881/aa99d4
  * Graur et al. (2017) Graur, O., Bianco, F. B., Modjaz, M., et al. 2017, The Astrophysical Journal, 837, 121
  * Guillochon et al. (2018) Guillochon, J., Nicholl, M., Villar, V. A., et al. 2018, The Astrophysical Journal Supplement Series, 236, 6
  * Guillochon et al. (2017) Guillochon, J., Parrent, J., Kelley, L. Z., & Margutti, R. 2017, The Astrophysical Journal, 835, 64
  * Guy et al. (2007) Guy, J., Astier, P., Baumont, S., et al. 2007, Astronomy & Astrophysics, 466, 11
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hosseinzadeh et al. (2020) Hosseinzadeh, G., Dauphin, F., Villar, V. A., et al. 2020, arXiv preprint arXiv:2008.04912
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * Ishida et al. (2019) Ishida, E. E., Kornilov, M. V., Malanchev, K. L., et al. 2019, arXiv preprint arXiv:1909.13260
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111, doi: 10.3847/1538-4357/ab042c
  * Jedicke et al. (2012) Jedicke, R., Tonry, J., Veres, P., et al. 2012, in AAS/Division for Planetary Sciences Meeting Abstracts# 44, 210–12
  * Jiang et al. (2020) Jiang, B., Jiang, S., & Villar, V. A. 2020, arXiv preprint arXiv:2008.10397
  * Kasen et al. (2017) Kasen, D., Metzger, B., Barnes, J., Quataert, E., & Ramirez-Ruiz, E. 2017, Nature, 551, 80
  * Kasen et al. (2011) Kasen, D., Woosley, S., & Heger, A. 2011, The Astrophysical Journal, 734, 102
  * Kasliwal (2012) Kasliwal, M. M. 2012, Publications of the Astronomical Society of Australia, 29, 482
  * Kasliwal et al. (2012) Kasliwal, M. M., Kulkarni, S., Gal-Yam, A., et al. 2012, The Astrophysical Journal, 755, 161
  * Kessler et al. (2010) Kessler, R., Conley, A., Jha, S., & Kuhlmann, S. 2010, arXiv preprint arXiv:1001.5210
  * Kessler et al. (2019) Kessler, R., Narayan, G., Avelino, A., et al. 2019, Publications of the Astronomical Society of the Pacific, 131, 094501
  * Kingma & Ba (2014) Kingma, D. P., & Ba, J. 2014, arXiv e-prints, arXiv:1412.6980.  https://arxiv.org/abs/1412.6980
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, ed. F. Loizides & B. Schmidt, IOS Press, 87 – 90
  * Li et al. (2019) Li, D., Chen, D., Shi, L., et al. 2019, arXiv e-prints, arXiv:1901.04997.  https://arxiv.org/abs/1901.04997
  * Li et al. (2003) Li, W., Filippenko, A. V., Chornock, R., et al. 2003, Publications of the Astronomical Society of the Pacific, 115, 453
  * Li et al. (2011) Li, W., Leaman, J., Chornock, R., et al. 2011, Monthly Notices of the Royal Astronomical Society, 412, 1441
  * Liu et al. (2012) Liu, F. T., Ting, K. M., & Zhou, Z.-H. 2012, ACM Transactions on Knowledge Discovery from Data (TKDD), 6, 1
  * MacLeod et al. (2011) MacLeod, C., Brooks, K., Ivezić, Ž., et al. 2011, The Astrophysical Journal, 728, 26
  * Martínez-Galarza et al. (2020) Martínez-Galarza, J. R., Bianco, F., Crake, D., et al. 2020, arXiv e-prints, arXiv:2009.06760.  https://arxiv.org/abs/2009.06760
  * Mockler et al. (2019) Mockler, B., Guillochon, J., & Ramirez-Ruiz, E. 2019, The Astrophysical Journal, 872, 151
  * Morawski et al. (2021) Morawski, F., Bejger, M., Cuoco, E., & Petre, L. 2021, arXiv e-prints, arXiv:2103.07688.  https://arxiv.org/abs/2103.07688
  * Muthukrishna et al. (2019) Muthukrishna, D., Narayan, G., Mandel, K. S., Biswas, R., & Hložek, R. 2019, Publications of the Astronomical Society of the Pacific, 131, 118002
  * Nicholl et al. (2017) Nicholl, M., Guillochon, J., & Berger, E. 2017, The Astrophysical Journal, 850, 55
  * Nugent et al. (2002) Nugent, P., Kim, A., & Perlmutter, S. 2002, Publications of the Astronomical Society of the Pacific, 114, 803
  * Pasquet et al. (2019) Pasquet, J., Pasquet, J., Chaumont, M., & Fouchez, D. 2019, Astronomy & Astrophysics, 627, A21
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825.  http://jmlr.org/papers/v12/pedregosa11a.html
  * Perley et al. (2020) Perley, D. A., Fremling, C., Sollerman, J., et al. 2020, The Astrophysical Journal, 904, 35
  * Pruzhinskaya et al. (2019) Pruzhinskaya, M. V., Malanchev, K. L., Kornilov, M. V., et al. 2019, Monthly Notices of the Royal Astronomical Society, 489, 3591
  * Rees (1988) Rees, M. J. 1988, Nature, 333, 523
  * Sako et al. (2018) Sako, M., Bassett, B., Becker, A. C., et al. 2018, Publications of the Astronomical Society of the Pacific, 130, 064002
  * Sánchez-Sáez et al. (2020) Sánchez-Sáez, P., Reyes, I., Valenzuela, C., et al. 2020, arXiv preprint arXiv:2008.03311
  * Shappee et al. (2014) Shappee, B., Prieto, J., Stanek, K., et al. 2014, in American Astronomical Society Meeting Abstracts# 223, Vol. 223, 236–03
  * Shemmer et al. (2018) Shemmer, O., Richards, G., Brandt, N., & et al. 2018, LSST AGN Science Collaboration Roadmap.  https://agn.science.lsst.org/sites/default/files/LSST_AGN_SC_Roadmap_v1p0.pdf
  * Villar et al. (2017) Villar, V. A., Berger, E., Metzger, B. D., & Guillochon, J. 2017, The Astrophysical Journal, 849, 70
  * Villar et al. (2018) Villar, V. A., Nicholl, M., & Berger, E. 2018, The Astrophysical Journal, 869, 166
  * Villar et al. (2019) Villar, V. A., Berger, E., Miller, G., et al. 2019, ApJ, 884, 83, doi: 10.3847/1538-4357/ab418c
  * Villar et al. (2020) Villar, V. A., Hosseinzadeh, G., Berger, E., et al. 2020, arXiv preprint arXiv:2008.04921
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: https://doi.org/10.1038/s41592-019-0686-2
  * Wu & Keogh (2020) Wu, R., & Keogh, E. J. 2020, arXiv e-prints, arXiv:2009.13807.  https://arxiv.org/abs/2009.13807
  * Zhang et al. (2019) Zhang, C., Song, D., Chen, Y., et al. 2019, in Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 33, 1409–1416
  * Zhao et al. (2020) Zhao, H., Wang, Y., Duan, J., et al. 2020, arXiv preprint arXiv:2009.02040

Here we investigate the AGN included in the PLAsTiCC simulation. In Figure 9,
we show the same anomaly samples shown in Figure 7 but including AGN. AGN
dominate the anomaly sample in all cuts. Our algorithm over-represents AGN in
our anomaly sample at the same rate as SLSNe and PISNe.

However, this is not likely not representative of reality. The PLAsTiCC AGN
are generated via a damped random walk with a structure function to define the
correlation between each filter, as described in MacLeod et al. (2011). We
find that, amongst the AGN identified as “anomalous” in this study, all AGN
have at least 20 data points spanning at least 60 days in duration. These
events would very likely be identified as AGN via other classification
methods, and be removed from our anomaly datastream. These events are an
example of a simple anomaly detection, in which simple filters would likely
also pick them out as distinct from SNe (see a summary of similar pitfalls in
Wu & Keogh 2020).

Finally, we note that being able to identify AGN in the LSST datastream is an
open problem (Shemmer et al., 2018). At first glance, our algorithm is
seemingly very successful at identifying AGN; however, these are likely AGN
which will be quickly identified by other means. Future studies will need to
identify extreme outbursts of AGN (not explicitly simulated by PLAsTiCC) and
AGN with substantial dust extinction (likely not identified in the LSST alert
stream unless they undergo a significant outburst). Searching for these
elusive AGN continue to be an open problem which require further development
of specialized classification methods.

Figure 9: Top row: Breakdown of anomalous events, given the cuts listed in
each panels. Cuts are identical to those of Figure 7, but now include AGN. AGN
dominate the samples in all cuts. Bottom row: Relative abundance of anomaly
samples shown in top rows. AGN are, like members of the minority class, highly
over-represented in the anomaly sample.