{'fraction_non_alphanumeric': 0.03898190427870793, 'fraction_numerical': 0.013867749027566378, 'mean_word_length': 4.313974387778027, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stable Invariant Models via Koopman Spectra

Weight-tied models have attracted attention in the modern development of neural networks. The deep equilibrium model (DEQ) represents infinitely deep neural networks with weight-tying, and recent studies have shown the potential of this type of approach. DEQs are needed to iteratively solve root-finding problems in training and are built on the assumption that the underlying dynamics determined by the models converge to a fixed point. In this paper, we present the stable invariant model (SIM), a new class of deep models that in principle approximates DEQs under stability and extends the dynamics to more general ones converging to an invariant set (not restricted in a fixed point). The key ingredient in deriving SIMs is a representation of the dynamics with the spectra of the Koopman and Perron--Frobenius operators. This perspective approximately reveals stable dynamics with DEQs and then derives two variants of SIMs. We also propose an implementation of SIMs that can be learned in the same way as feedforward models. We illustrate the empirical performance of SIMs with experiments and demonstrate that SIMs achieve comparative or superior performance against DEQs in several learning tasks.


Introduction
A feedforward neural network learns a representation by explicitly iterating a number of layer-by-layer computations. Each layer performs a transformation of outputs from the previous layer, which is typically characterized by different sets of parameters among the layers. However, several recent studies have shown that models with weight-tying, i.e. the ones employing the same transformation in each layer, achieve results competitive with stateof-the-art performances (Bai et al., 2019b;Dabre & Fujita, 2019;Dehghani et al., 2019). Motivated from this fact, Bai et al. (2019a) recently proposed the deep equilibrium model (DEQ), which is equal to running an infinitely deep feedforward model with weight-tying instead of using a finite number of layers. The models compute a representation by finding a fixed point (or an equilibrium point) with root-finding in practice, and are thus regarded as an instance of the so-called implicit-depth models such as neural ordinary differential equations (Chen et al., 2018). The following studies on DEQs have shown capability in this type of approach in several learning tasks (Bai et al., 2020;Winston & Kolter, 2020).
The forward pass of DEQs and their variants involves solving root-finding problems, which can lead to a high computational time and tends to be unstable (regarding, for example, the sensitivity in hyper-parameter tuning and initialization). Thus, DEQs sometimes require extensive and time-consuming tuning to achieve strong performance and convergence to solutions (Linsley et al., 2020). Additionally, DEQs assume that the underlying dynamics determined by the models converge to a fixed point. However, as is often reported in papers on sequential neural models such as recurrent neural networks and also known in the scientific studies of brain activity, a broader class of convergence such as nonlinear oscillations could convey preferred capabilities in learning (Selverston & Moulins, 1985;Townley et al., 2000;Chang et al., 2019;Kag et al., 2020).
In this study, we propose a novel class of deep models, referred to as the stable invariant model (SIM). The key insight behind SIMs is to interpret the underlying dynamics regarding DEQs through the Koopman operator. The Koopman operator is a linear operator over functions defined on latent states of dynamics (Koopman, 1931;Mezić, 2005). Because of its linearity, we can capture inherent temporal and spatial patterns of dynamics by the representation of the spectra, i.e. eigenvalues and eigenvectors, of the corresponding Koopman operator. We first show that the spectra of the Koopman operator clarify stable dynamics with DEQs and then develop two variants of SIMs. The first models, single-tier SIMs, principally approximate DEQs under the stability. The resulting models, somewhat surprisingly, consist of only threestep transformations, although they approximate DEQs that are infinitely deep (see Figure 1). Moreover, the second models, two-tier SIMs, extend the dynamics to broader ones converging to an invariant set (e.g. a set of points, curve, and more general manifold) by employing the connection between the Koopman and its conjugate (Perron-Frobenius) operators. We further provide a practical scheme to implement SIMs so that they can be learned in the same manner as feedforward models. Finally, we illustrate the behaviors of SIMs with numerical experiments in supervised learning tasks. We demonstrate that our models achieve competitive or superior performances compared to DEQs with less computational time. The remainder of this paper is organized as follows. First, in Section 2, we briefly review DEQs, and the Koopman and Perron-Frobenius operators. In Section 3, we propose SIMs along with the description of their resulting architectures and characteristics. We describe the related works in Section 4 and investigate the empirical performance of our models in three learning tasks in Section 5. We conclude this paper in Section 6. The details of some equation derivations and experiments are presented in Appendix A and Appendix B, respectively.

Deep Equilibrium Models
One of the core ideas in DEQs is weight-tying, i.e. the same set of parameters is shared across the layers of a deep network. Formally, DEQs consider an L-layer weight-tied transformation with shared parameters θ: where x ∈ R D is the input to the model, z l ∈ R d is the hidden state of the l-th layer, and f θ : R d+D → R d is a continuous function. DEQs suppose that stacking such layers infinitely tends to a fixed point: The forward pass of DEQs uses root-finding algorithms to directly compute the fixed point z * by solving the equation z * = f θ (z * , x). Then, z * is transformed to an output y by a function h as y = h(z * ). We can train DEQs with backpropagation by computing the gradient of the fixed point through implicit differentiation (Krantz & Parks, 2013). A series of transformations in Eq. (1) can be viewed as a discrete-time nonlinear dynamical system where the hidden state z l is a state vector at step l. The underlying dynamics are determined by the transformation f θ , which is affected by the input x at every step.

Koopman and Perron-Frobenius Operators
We briefly overview the Koopman and Perron-Frobenius operators. Please see other references, e.g. (Mauroy et al., 2020), for more information.
Definitions. Consider a discrete-time nonlinear dynamical system: z t+1 = f (z t ), defined on a state space S ⊂ R d , where z t ∈ S is the state vector at time t, and f is a (possibly, nonlinear) state-transition function. The Koopman operator (Koopman, 1931;Mezić, 2005), which we denote here by K, is an infinite-dimensional linear operator that acts on a scalar observable function g : S → C, mapping g to a new function Kg, given by where z ∈ S is a state vector, and • denotes the composition of g with f . K acts linearly on the observable g, although the dynamics described by f may be nonlinear. The Perron-Frobenius operator P describes the transition of the density over the state. That is, if we denote the density by p t over z, then P is defined by where δ is the Dirac delta function. It should be noted that, from the definition, the Koopman (or Perron-Frobenius) operator defined on a space always has its unique conjugate as the Perron-Frobenius (or Koopman) defined on the corresponding space.
Koopman Spectrum. Because K and P are linear operators, they have eigenvalues and the associated eigenfunctions although their numbers are not generally finite. Let λ j ∈ C and φ j : S → C for j = 1, 2, . . . be the j-th eigenvalues and eigenfunctions of K, respectively. A Koopman eigenvalue λ j and the eigenfunction φ j satisfy the relation This leads to an expression of the dynamics: where the coefficients v j ∈ C for j = 1, 2, . . . are referred to as the Koopman modes associated with g.
Finite-Dimensional Approximation. Let ϕ j : S → C be N basis functions (j = 1, 2, . . . , N ). If g exists on the subspace spanned by {ϕ j }, then g is represented as a linear combination of the basis functions, i.e. g(z) = w ϕ(z), where we denote the concatenation of the basis functions as a vector-valued one ϕ = (ϕ 1 , . . . , ϕ N ) : S → C N , and w ∈ C N is the coordinate of g on the subspace. By projecting the action of K onto the span of {ϕ j }, we approximate K with another linear operator K N . This approximates the Koopman operator K and satisfies where K ∈ C N ×N is referred to as the Koopman matrix. The above shows that a temporal evolution of g with K N is represented by applying K to the coordinate w. Therefore, the Koopman matrix owns the one-to-one correspondence to K N . Moreover, the Koopman eigenvalues, eigenfunctions, and modes of K N can be obtained from the eigenvalues, right-eigenvectors, and left-eigenvectors of K, respectively.
Additionally, the Koopman matrix provides an approximation of the dynamics through the basis functions: where A is the transpose of the Koopman matrix, i.e. A = K . Eq. (3) indicates that a temporal evolution with f can be approximated by a temporal evolution with the finite and linear dynamics described by A over a lifted space with ϕ. The equation holds if K = K N . It should be noted that, if g is a real-valued function, then it is sufficient to consider real-valued ones for the corresponding quantities (such as ϕ j , w and K).

Stable Invariant Models
In this section, we introduce our proposed SIMs. The models are motivated by two problems concerning DEQs. First, DEQs work under the assumption that the underlying dynamics are stable, i.e. Eq. (2) holds for any input and initial state. However, it is in general hard to assess whether nonlinear dynamics satisfy the assumption. We address this problem by approximately specifying the stable dynamics using the Koopman spectrum, which leads our first single-tier SIMs by identifying the convergent behavior for a fixed point.
The second problem is that DEQs only consider convergence to a fixed point. In the literature on dynamical systems, one often considers more general convergence to an invariant set: for a dynamical system on a state space S, a set S ⊂ S is said to be an invariant set if any trajectory starting in S remain in S. The convergence to an invariant set allows for the dynamics that oscillate on the set. Notable examples include limit cycles, tori, and other nonlinear oscillations (Strogatz, 2015). We also refer to dynamics as stable if the dynamics converge to a bounded invariant set for any input and initial state and construct our second two-tier SIMs by incorporating only the dynamics converging to an invariant set via the Koopman and Perron-Frobenius operators.

Approximating DEQs
We begin by considering the representations of DEQs with the Koopman operator. However, it should be noted that f θ Eq. (1) has another vector x as an input differently from f in Section 2.2. Although there are several ways to define the Koopman operator for such a case, we consider the following Koopman operator K acting on the observable g of both the hidden state z and input x (Proctor et al., 2018): The action of this operator is restricted so that x is maintained at the same point. Even if x is injected, the properties of the Koopman operator discussed in Section 2.2 still hold.
Hereafter, we consider a real-valued observable g. First, let ϕ j (z, x) : R d+D → R be N real-valued basis functions (j = 1, . . . , N ). The concatenation is given by a vector-valued function ϕ. Here, we assume that there exists the inverse function ϕ −1 , and the action of the Koopman operator always stays in the subspace spanned by {ϕ j }. Then, if Eq. (1) converges to a fixed point z * as in Eq.
(2), we can approximate z * with x by the Koopman operator of Eq. (4) as where the last approximation follows the finite-dimensional approximation of the Koopman operator through the lifted dynamics as described in Eq. (3). By dividing the inverse function ϕ −1 into two parts corresponding to z and x, i.e. ϕ −1 = (ϕ −1 z , ϕ −1 x ), the fixed point z * can be approximated more directly as

Convergent Behavior via Koopman Spectra
The approximation (5) assumes that DEQs converge to a fixed point. However, whether a DEQ converges to a fixed point depends on the behavior of the underlying dynamics with the DEQ. The approximation (5) allows us to characterize the convergent behavior of a DEQ via eigenvalues of the corresponding A.
First, we denote by λ j ∈ C for j = 1, 2, . . . , N the eigenvalues of A. It should be noted that the eigenvalues can be complex values even though A is a real matrix (because it is not necessarily symmetric). For any real matrix A, there exists a nonsingular matrix U ∈ R N ×N that consists of the generalized-eigenvectors including the ordinal eigenvectors of A. We represent U and U −1 with N vectors, respectively, as U = (u 1 , . . . , u N ) and U −1 = (v 1 , . . . , v N ) , where u j is associated with a generalized-eigenvector of λ j if λ j is real, and the real or imaginary part of a generalized-eigenvector of λ j if λ j is non-real. Additionally, let ρ(A) be the spectral radius of A, i.e. ρ(A) := max{|λ 1 |, . . . , |λ N |}. We can then classify the convergent behavior of the lifted dynamics into the following four cases: (i). If ρ(A) < 1, then the dynamics converge to the origin: (ii). If ρ(A) = 1, all the eigenvalues with |λ j | = 1 take the values 1, and their corresponding eigenvectors are linearly independent, then the lifted dynamics converge to a fixed point. That is, if we denote by J 1 = {j | λ j = 1} the index set of such eigenvalues, then we have (iii). If ρ(A) = 1 and the eigenvectors with eigenvalues of |λ j | = 1 are linearly independent, then the lifted dynamics do not converge to a point but oscillates in the state space. More concretely, if we denote Re where i is the imaginary unit, and ∆ j = arctan(β j /α j ). J 3 is the set of index pairs whose eigenvalues are conjugated. It should be noted that there always exists a conjugate eigenvalue for every non-real eigenvalue when A is real. (iv). Otherwise, at least one element among the states of the lifted dynamics diverges.
We describe the derivations of the four cases in Appendix A. Figure 2 shows the complex plane and coordinates of the eigenvalues described in the above cases. The red and green dots, and the blue line denote the locations of the eigenvalues that correspond to J 1 , J 2 , and J 3 , respectively. The gray area shows the area where the absolute values of the eigenvalues are less than 1. The first three cases approximately correspond to the stable dynamics with the DEQ. Case (i) converges to the origin regardless of x, which is useless for any learning problems. In contrast, case (ii) converges to a fixed point that reflects x. Case (iii) does not converge to any fixed point but oscillates on a manifold in the state space, where the terms corresponding to J 2 and J 3 respectively include coefficients such as (−1) l and cos(l∆ j ) that neither diverge nor converge. Lastly, case (iv) diverges and the corresponding dynamics are not stable.

Model Description
We derive SIMs based on the above analysis of DEQs from the perspective of the Koopman spectrum. We first describe the variant that approximates DEQs by leveraging case (ii). Moreover, we present another variant of SIMs that incorporates a broader class of dynamics by utilizing case (iii).
Fixed Point. Case (ii) represents the dynamics that converge to a fixed point as in Eq. (6). Because DEQs assume that the underlying dynamics converge towards a fixed point, which are covered by case (ii), we can approximate the fixed point of a DEQ by plugging Eq. (6) into Eq. (5): whereÛ ∈ R N ×K andV ∈ R K×N consist of the vectors that are the rows of U and V whose indices belong to J 1 and K is the size of J 1 . Interestingly, this approximation implies that a fixed point of DEQs can be represented as a finite-depth model that consists of three-step transformations, under the assumption that the corresponding Koopman operator can be approximated by a finite set of basis functions. We refer to the right-hand side of Eq. (8) single-tier SIMs. This comes from the fact that a state of the dynamics is lifted to another tier with ϕ.
Invariant Set. Case (iii) represents the dynamics that converge to an invariant set given by Eq. (7). Such areas perturbed by noises are known as almost invariant sets in physics, and can be characterized by spectra with eigenvalues 1 of the corresponding Perron-Frobenius operator (Billings & Schwartz, 2008;Froyland & Padberg, 2009). That is, the convergent trajectories that are characterized by the spectra of case (iii) of the Koopman operator can be given as the ones in case (ii) of the corresponding Perron-Frobenius operator. Now, let ψ be the (finite-dimensional) basis functions of an embedding of the state into some Hilbert space that encodes the density defining the corresponding Perron-Frobenius operator. The above insight allows us to represent the extended equilibrium of the states as follows: where J 1 is the index set of the spectra of the Perron-Frobenius operator analogous to J 1 , andÛ andV are the corresponding matrices. This new point z reflects the convergent behavior of the dynamics over an invariant set. We refer to this type of models as two-tier SIMs, which come from the structure involving two-tier maps by ϕ and ψ. Two-tier SIMs are also represented as finite-depth models. Figure 3 shows the overview of SIMs. The black arrows denote the first transformation with the basis functions ϕ. The blue and red arrows correspond to the transformations of the single-tier and two-tier SIMs, respectively. Each model considers the convergent behavior of the dynamics in the first-tier or second-tier space, and then returns to the original space with the inverse functions.
Although DEQs are defined by the transformation f θ , SIMs require the specification of the basis functions ϕ and ψ. One natural question may be how the two approaches and convergence properties relate to each other. Generally, if two dynamical systems converge to similar fixed points or invariant sets, the corresponding dynamics also have similar topological properties. However, the actual transformations are determined by many factors (e.g. the training algorithms and choice of basis functions), and the convergence regions of the two approaches will be different. Investigating how f θ can be associated with ϕ and ψ is an important future work.

Practical Implementation
The key to encode high expressiveness in the dynamics with SIMs is to utilize rich basis functions. Therefore, we leverage neural networks with structures to construct the components in our models (8) and (9) as follows.
First, we utilize a neural network whose inputs are only x, which we denote by µ NN (x), to approximate ϕ(z 0 , x). This is because the initial state z 0 can be basically assumed to be zero. Although we can prepare a non-zero z 0 (e.g. the original implementation of DEQs allows to use the final state of the previous step in the training loop), the effectiveness of such a biased initial state is unclear. Therefore, we assume that z 0 is always zero and encoded in the basis functions, and thus drop z 0 from the input of ϕ.
Next, becauseÛ andV in Eq. (8) are simply real matrices, we can deal with them as two consecutive linear layers althoughÛ andV respectively consist of linearly independent vectors. Although we could consider the constraint, we ignore it because the vectors in a learned matrix will rarely be linearly dependent.V has another implicit constraint:V originates from U −1 and thus is affected byÛ . However,V is also influenced by the vectors of U which correspond to eigenvalues with an absolute value less than 1. Although these vectors disappear by taking the limit and do not appear in the model,V will have degrees of freedom thanks to the vectors. From this observation, we modelV independently fromÛ .
Further, we prepare another neural network for the inverse functions, ϕ −1 z and ψ −1 . This network is, principally, required to ensure its invertibility corresponding to the respective part of ϕ and ψ. However, modeling a neural network to be invertible is costly and even difficult unless the dimension of the lifted space is identical to that of the original space (Papamakarios et al., 2021). The invertibility will also require SIMs to consider ϕ −1 x that is omitted from Eq. (5). This is because it is necessary to fully construct ϕ −1 including ϕ −1 x to guarantee the invertibility even though the output of ϕ −1 x is never used in subsequent transformations. In this paper, we focus on the practical aspect in our implementation and simply approximate ϕ −1 z in Eq. (8) or the composition ϕ −1 z •ψ −1 in Eq. (9) by a neural network ν NN without restriction of invertibility. This approach follows the same manner as common encoderdecoder models; a encoder model is usually designed independently of the corresponding decoder model despite their close connection of going back and forth between original and latent spaces.
where ω j ∈ R N for j = 1, . . . , M/2 are random vectors drawn from a probability distribution P (ω), to avoid the increase of the computational cost along the sample size when using reproducing kernels. Because ψ is the M -dimensional real-valued function,Û andV are defined as M × K and K × M real matrices, respectively, where K is the size of J 1 . We model those matrices in the same way asÛ andV .
Putting the above pieces together, we model Eqs. (8) and (9) by z single-tier and z two-tier , respectively, as follows: (10) In the end, those implementations of SIMs are realized as feedforward models. The implementations no longer require any root-finding algorithms for the forward pass and implicit differentiation for the backward pass. In the following experiments, we will instantiate more specific examples of the implementations for each task. However, the implementations do not limit themselves to such particular forms and have the flexibility to take various configurations by changing µ NN and ν NN . A practical merit of DEQs compared to feedforward models is the memory efficiency: the required memory does not depend on the number of transformations by f θ owing to root-finding and implicit differentiation. Because we realize SIMs as feedforward models, they do not fully inherit this efficiency. However, DEQs also model f θ with feedforward architectures such as Transformers (Vaswani et al., 2017;Bai et al., 2019a) and hence need the memory for f θ , which could be potentially large. Additionally, we notice that DEQs are infinitely deep. The naive approximation with a feedforward model may need to increase the depth of the feedforward model. Our derivation showed that it is critical to approximate the DEQs according to the form of Eq. (10). If the basis functions of SIMs could be represented by shallower models, it will save a lot of memory compared to naive approximations. We can regard SIMs as an intermediate approach that can approximate infinitely deep models while avoiding the simple dependence on depth.

Relation to Implicit Neural Representation
Implicit neural representation has recently gained interest in the machine learning community. The aim is to represent complex natural signals (e.g. images, 3D objects, and audio waves) with a function modeled by a neural network (Sitzmann et al., 2020;Tancik et al., 2020).
SIMs have an interesting connection to implicit neural representation by (Tancik et al., 2020). Tancik et al. (2020) proposed to use a Fourier feature mapping as a pre-processing. If the basis functions modeled by µ NN of singletier SIMs are modeled by an RFF ψ, then the implementation of the SIMs is the same as the Fourier feature mapping by (Tancik et al., 2020) up to the termÛV : As in Eq. (10), z 0 is omitted from ψ(x) because it is supposed to be zero and does not affect the output of ψ. Alternatively, this RFF only model can be interpreted as an instance of two-tier SIMs when µ NN is an identity mapping.

Related Works
The origin of DEQs dates back to the work on recurrent backpropagation (RBP) (Almeida, 1987;Pineda, 1987). Those studies proposed to utilize early implicit-depth models and have been applied to other studies, e.g. graph neural networks (Gori et al., 2005;Scarselli et al., 2008). Recently, Liao et al. (2018) revisited the RBP algorithm and improved it with the conjugate gradient method and Neumann series. Bai et al. (2019a) introduced a perspective of the use of a fixed point as a replacement for depth and called the approaches the DEQ. This study also proved the universality of a single DEQ layer and developed a practical quasi-Newton method that works for large-scale sequential tasks. Subsequent studies have reported the theoretical analysis (Kawaguchi, 2021;Pabbaraju et al., 2021) and proposed improved architectures (Bai et al., 2020;Xie et al., 2021). Several relevant studies focused on the stability issue of DEQs. Winston & Kolter (2020) proposed monotone DEQs (monDEQs) with guaranteed convergence to a fixed point based on monotone operators. Bai et al. (2021) also proposed to stabilize the training using Jacobian regularization.
The Koopman operator has been known for a long time as a tool for analyzing dynamics in physics (Koopman, 1931). It has recently received attention that the spectral properties of the Koopman operator play an important role in revealing global characteristics of the underlying dynamics (Mezić, 2005). Because the Koopman operator is linear but infinite-dimensional, a major challenge is how to compute the spectra of the Koopman operator in practice. Dynamic mode decomposition (DMD) (Schmid, 2010) has gained popularity as a data-driven approach to computing a reasonable finitedimensional approximation of the Koopman spectra (Rowley et al., 2009). While the original DMD algorithm supposes that the basis functions are linear, several subsequent studies have been proposed to utilize nonlinear basis functions (Williams et al., 2015), a kernel-based approach (Kawahara, 2016), a Bayesian formulation (Takeishi et al., 2017a), and learning from data (Takeishi et al., 2017b).
The Koopman operator is also applied to machine learning and the related fields. Dogra & Redman (2020) leveraged the representation with the Koopman operator to accelerate the training of neural networks. Manojlović et al. (2020) analyzed the dynamics in the training of neural networks with the Koopman operator and proposed a method to characterize the architectures of neural networks with the Koopman spectrum. Takeishi & Kawahara (2021) proposed a neural network to learn stable dynamical systems by utilizing a map obtained by, for example, an eigenfunction of the Koopman operator (although it is not described explicitly in the paper) by extending stable deep dynamics models by Manek & Kolter (2019).

Experiments
We evaluated SIMs through experiments on three supervised learning tasks. For all tasks, we first trained models with candidates of hyperparameters on training data and evaluated them on validation data. We then selected the best hyper-parameter, re-trained the model with the best one, and evaluated the trained model on test data. We implemented our proposed models and training algorithms using Pytorch (Paszke et al., 2019). For the basic building blocks and optimizers, we used the existing implementation in Pytorch. We tuned the hyper-parameters using Tune (Liaw et al., 2018). In all experiments, we set the search algorithm to the random search and trial scheduler to the asynchronous successive halving algorithm (ASHA) (Li et al., 2020). We conducted the experiments on an internal computing server with Intel Xeon Bronze 3104 CPUs and NVIDIA V100 GPUs.
We set µ NN as a task-specific neural network for each task. For ψ, we used a normal distribution to sample random vectors. For ν NN , we used a three-layer fully connected network (FCN) with ReLU activation (Nair & Hinton, 2010) for all tasks in common. We used a linear function as h to fit the dimension of the hidden state to the output size and then the softmax function if the addressed task is a classification problem. For SIMs, we set K = N/2 and K = M/2, respectively.
The details of the datasets, model architectures, and training algorithms are described in Appendix B.

Copy Memory Task
We first report the results of the copy memory task (Hochreiter & Schmidhuber, 1997) to evaluate the effectiveness of SIMs against DEQs. The goal of the copy memory task is to predict a sequence of digits of length T + 20 from another input sequence of the same length. The first ten elements of the input sequence consist of digits randomly drawn from 1 to 8, the subsequent T − 1 elements are filled with 0, and the last eleven elements are all 9. Given this input, the first T + 10 elements of the output sequence are 0 and the last ten elements are the same as the first ten ones of the input sequence. Hence, this task evaluates how well a model can remember the first elements of an input. In the experiment, we set T to 500 and followed the experimental procedure of (Bai et al., 2018).
For SIMs, we applied the temporal convolutional network (TCN) architecture to model µ NN by following the implementation of (Bai et al., 2018): the architecture includes 1D dilated causal convolution, ReLU activation, and residual connection (He et al., 2016). Because the method of applying the TCN is rather complicated compared to other tasks, we explain the case of the single-tier SIM as an example. Formally, if we denote x ∈ {0, 1, . . . , 9} 520 as an input sequence, the TCN outputs the sequence of states: where W = (w 1 , . . . , w 520 ) ∈ R N ×520 and w j ∈ R N denotes a state on the lifted space at the j-th position in the sequence. Each state feeds the common architecture: where z single,j ∈ R d is the hidden state of the j-th position. The prediction at the j-th position of the output sequence is obtained from z single,j . We can interpret this architecture as follows; for this task, the model has different µ NN for each position, but it transforms all the positions to the hidden states at once by the TCN, and then the same subsequent transformation is applied to all the positions. Although we could use different ν NN andÛV for each position, the above architecture can retain memory in practice while keeping the model size small. For the DEQ, we applied the Universal Transformer (Dehghani et al., 2019) as the base function f θ and mostly adopted the original implementation of (Bai et al., 2019a). However, the public source code was optimized for the tasks of language modeling; hence, we mainly modified the following two parts of the implementation. First, the implementation uses the adaptive softmax function (Baevski & Auli, 2019;Grave et al., 2017) to address the large word vocabulary. Because the copy memory task considers only ten digits, we did not use this architecture in the experiment. Second, the implementation uses the memory padding and nonzero initial hidden states that employ the final states of the previous step in a training loop. We tested the two techniques but observed that the test loss was considerably worse even though the training loss was fine. Hence, we omit the techniques and instead used the empty memory and zero initial hidden states in the experiment. Table 1 shows the test cross-entropy loss per sequence and the number of learnable parameters for each model. Although the DEQ achieved the lowest loss, the single-tier SIM obtained was comparable with almost the same number of learnable parameters. Figure 4 (a) shows the progress of the training losses along the run-time for 20 epochs when we re-trained the models. The training speed depends on the batch size. During the hyperparameter search, 1, 1, and 10 were selected as the batch size for the DEQ, single-tier SIM, and two-tier SIM, respectively. Hence, the two-tier SIM was the fastest to complete the training. For the DEQ and single-tier SIM, the single-tier model was ten times faster than the DEQ under the same batch size. Figure 4 (b) shows the progress of the training losses along the epochs. We can observe that the single-tier SIM converged faster than the DEQ.

Image Classification
We next report the results of the image classification to compare SIMs to monDEQs. We followed the experiment of image classification conducted in (Winston & Kolter, 2020). We prepared the CIFAR-10 (Krizhevsky, 2009), SVHN (Netzer et al., 2011), andMNIST (LeCun et al., 1998) datasets, which contain images in 10 different classes and evaluated the classification performance in the standard setting. Following (Winston & Kolter, 2020), we evaluated the models on test data three times with different initialization and reported the averaged performance.
For SIMs, we employed convolutional neural networks for µ NN . Following the VGG models (Simonyan & Zisserman, 2015), µ NN consists of two convolutional layers each of which has two convolution filters with ReLU activation and batch normalization and one max pooling. The output of the convolution layers is additionally transformed by a linear layer to fit the dimension of the lifted space. We constrained the number of learnable parameters of SIMs to be comparable to the one of monDEQs in (Winston & Kolter, 2020). Table 2 lists the means of the test accuracy over three runs of monDEQs and SIMs with different initialization. SIMs showed comparable or better performances against monDEQs for all datasets. Particularly, the performance was improved for the CIFAR-10 dataset even though the SIMs have a similar number of learnable parameters to monDEQ ones.

Image Regression
Finally, we report the results of the image regression task, which is an example of implicit neural representation tasks. In Section 3.4, we found that SIMs have a close connection to the work of (Tancik et al., 2020). The purpose of the experiment is to verify 1) that the RFF only model actually works well for this task, and 2) how well the other types of SIMs perform. The goal of the task is to obtain a neural network where the input is a 2-D pixel coordinate and the output is its 3-D RGB value. Following (Tancik et al., 2020), we evaluate SIMs with 32 datasets, where 16 are natural images 1 and the rest are text images 2 . This task considers each of the images to be one dataset. For an image, we picked 1/4 pixels as training data and other 1/4 as test data. Additionally, we prepared another 1/4 pixels as validation data to select the hyper-parameters.
We compare three instances of SIMs. The first two models employ an FCN for µ NN : it consists of two linear layers with ReLU activation and one linear layer to fit the dimension of the lifted space. The last one is the RFF only model (11). Table 3 shows the mean of the test PSNR over 16 images of two types of datasets. The results indicate that the RFF only model was better than the single-tier SIM. This result is consistent with the original work on (Tancik et al., 2020). The two-tier SIM was also effective although the test PSNRs were slightly worse than the RFF only model. Figure 5 illustrates examples of prediction for two images. Although the single-tier SIM produced blurred images, the RFF only model can generate sharper ones. Although the images of the two-tier SIM are a little blurry, the detail can be recognized compared to the single-tier SIM.

Discussion
The first two results indicate that SIMs can provide time-effective alternatives to DEQs. This is because SIMs are implemented as feedforward models, which do not require implicit differentiation even though they can approximate DEQs. The results also showed that two-tier SIMs did not necessarily outperform single-tier SIMs for the first two tasks. One of the reasons would be that even single-tier SIMs suffice for the tasks where the models can capture the essential feature of sequences and images. The third result demonstrated that SIMs are also effective in the task of implicit neural representation. Particularly, we observed that the two-tier SIM also obtained moderate performance. The work of (Tancik et al., 2020) does not cover the form of the two-tier models, and we believe that this result is an interesting finding.

Conclusions
In this study, we considered DEQs from the viewpoint of the Koopman operator, which enables us to identify stable dynamics described by DEQs via the representation of the spectra. This perspective yielded our proposed SIMs that approximated DEQs and exploited more general dynamics converging to an invariant set. Despite having such noteworthy properties, the resulting models can be represented as simple feedforward models, which will provide new insights into the studies on DEQs.
We clarify the convergent behavior of the lifted dynamics in Section 3. It should be noted that the following result is not novel; it employs a known consequence of linear algebra.
First, we define the l-th step of the lifted dynamics as where ϕ 0 := ϕ(z 0 , x). To reveal the property of ϕ l , we decompose A l with the eigenvalues and eigenvectors. We begin with the decomposition of A. If A is diagonalizable, it can be represented by the set of the eigenvalues and corresponding eigenvectors. However, if A has repeated eigenvalues, A is not necessarily diagonalizable. Although the Jordan canonical form can be used in such a case, this form is constructed by complex matrices if A contains complex eigenvalues. To represent A by real matrices for convenience, we consider the real Jordan canonical form: any real square matrix can be written as where U is defined in Section 3.2, and J ∈ R N ×N is the following block diagonal matrix: where R ≤ N corresponds to the number of linearly independent eigenvectors. Each block J r ∈ R sr×sr (r = 1, . . . , R) is associated with s r repeated eigenvalues with a value λ (r) if λ (r) is real or s r /2 repeated eigenvalues with a value λ (r) and their s r /2 complex conjugates with a value λ (r) if λ (r) is nonreal. If λ (r) is real, J r is also associated with one ordinal eigenvector of λ (r) and s r − 1 non-ordinal generalized-eigenvectors. If λ (r) is nonreal, J r is associated with a complex conjugate pair of ordinal eigenvectors of λ (r) and λ (r) and (s r − 1)/2 complex conjugate pairs of non-ordinal generalizedeigenvectors. It should also be noted that there always exists a conjugate eigenvalue for every nonreal eigenvalue when A is real. J r takes one of the following two forms: Here, C r and I are defined by As a special case, J r,R can be a scalar λ (r) when s r = 1, and J r,C can be C r when s r = 2. In this case, J r,R and J r,C are only associated with one ordinal eigenvector and a complex conjugate pair of ordinal eigenvectors, respectively. Moreover, the matrix C r can be rewritten as a polar form: where r (r) = |λ (r) | = α 2 (r) + β 2 (r) , ∆ (r) = arctan(β (r) /α (r) ).
The real Jordan canonical form can represent ϕ l as Each J l r can be written as cos ∆ (r) sin ∆ (r) − sin ∆ (r) cos ∆ (r) l 5 . If |λ (r) | = 1, λ (r) is nonreal, and s r = 2, then J l r neither converges nor diverges and takes the form of 6 . If |λ (r) | = 1, λ (r) is nonreal, and s r > 2, then the non-block-diagonal parts of J l r (= J l r,C ), e.g. l 1 C l−1 r , diverge. Hence, if at least one block falls into cases 4 or 6 , the lifted dynamics diverges: the complexity of the non-diagonal elements or non-block-diagonal parts is at least O(l). This case falls into case (iv).
Case (ii) means all the blocks are either case 1 or 2 , and each block in case 2 corresponds to one of the eigenvalues in J 1 . Hence, the eigenvectors with eigenvalues of λ j = 1 are linearly independent. In the end, J l converges to a diagonal matrix P J 1 ∈ R N ×N where the diagonal elements corresponding to J 1 are 1 and the rest are 0, and we obtain Eq. (6): Case (iii) means that all the blocks are either case 1 , 2 , 3 , or 5 . Each block in cases 3 and 5 corresponds to one of the eigenvalues of J 2 and the complex conjugate pairs of the eigenvalues of J 3 , respectively. Hence, the eigenvectors with eigenvalues of |λ j | = 1 are linearly independent. In addition to P J 1 , if we denote P J 2 ∈ R N ×N as a diagonal matrix where the diagonal elements corresponding to J 2 are -1 and the rest are 0 and P J 3 ∈ R N ×N as a block diagonal matrix where the diagonal blocks corresponding to the pairs of J 3 take the form of Eq. (A.3) and the rest are zero, we obtain Eq. (7): Lastly, if ρ(A) > 1, at least one block of J l diverges. This case falls into case (iv).

Appendix B. Experimental Details
We show the detailed settings and configurations of the experiments in Section 5. It should be noted that we used the default values of Pytorch for the arguments not mentioned below. Table B.4 shows that the settings of Tune and ASHA for each task. #trials denotes the number of trials in running a search algorithm. Scope denotes the method of selecting the best hyper-parameters. We used last that selects the best one by comparing the last performance at the end of training. #max t denotes the number of maximum epochs or iterations in a trial. ASHA periodically stops trials when the specified metric is poor and reduces them by a factor of the reduction factor. However, any trial is run until the grace period. Those settings are partially different for each task but the same as compared methods in a task. Table B.5 shows the detailed architecture of the reverse model ν NN in SIMs. The first column denotes the layers of ν NN , which transform the input in order from the top to the bottom layers. "× 2" means that the same transformation has been repeated twice. The second column denotes detailed information about the layers. The variable #hidden units is set to a different value for each task.
Appendix B.1. Copy Memory Task  Table B.6 lists the statistics of the dataset in the copy memory task. #train, #valid, and #test denote the number of training, validation, and test data, respectively. We first trained the models with the candidates of hyper-parameters on the training data and selected the best hyper-parameter on the validation data. We then used both training and validation data to re-train the model with the best hyper-parameter and evaluated it on the test data. A data point of this dataset consists of the pair of an input sequence of length 520 and an output sequence of length 520. We generated all the data points by following the procedure of (Bai et al., 2018). Table B.7 shows the model architecture of the TCN. The first building block is the temporal convolution layer where we adopted the architecture of (Bai et al., 2018): it mainly consists of the 1D dilated causal convolution, ReLU activation, and residual connection, and down-sampling with 1D convolution is applied to the residual connection (He et al., 2016) except for the first layer. We stacked this layer eight times and then connected a linear layer to fit the dimension of the lifted space. Table B.8 lists the configurations of the compared models. The upper and lower rows denote the configuration of the model architectures and training algorithms, respectively. The configuration with the set notation (e.g. [1e-3, 1e+3) and {1, 2, 5, 10, 20}) means that the corresponding variable was tuned as a hyper-parameter within the range of the set. For each trial, ASHA sampled a candidate from the set uniformly at random. For the configuration with a real interval, it was sampled in the logarithmic space with base 10. In this task, we tuned the bandwidth of the RFF, batch size, and learning rate. The configuration of the DEQ indicates that the variables of the implementation of (Bai et al., 2019a): we set it to ensure almost the same number of learnable parameters as SIMs. For the training algorithm of the DEQ, we adopted Adam (Kingma & Ba, 2015) with the step-wise cosine decaying schedule for the learning rate as in (Bai et al., 2019a). For SIMs, we used Adam with a constant learning rate schedule and slightly larger to improve the stability of the training algorithm. We also used gradient clipping with the default values of the DEQ and TCN. For all compared models, we minimized the cross-entropy loss, i.e. the negative log-likelihood with the softmax function. Although every model contains the Dropout  Table B.9 shows the statistics of three datasets in the image classification task. We pre-processed these datasets in almost the same way as (Winston & Kolter, 2020). The difference from (Winston & Kolter, 2020) is to additionally prepare the validation data. We initially used 90% of the prepared training data as training data and the remaining 10% as validation data. We then used both of them to train models with the best hyper-parameters for evaluating the performance on the test data. Table B.10 shows the architecture of µ NN in SIMs. We first used a combination of convolution and max-pooling layers, which follows the VGG architecture (Simonyan & Zisserman, 2015) and then applied a linear layer to fit the dimension of the lifted space. Table B.11 lists the configuration of SIMs. In this task, we tuned the bandwidth of the RFF, the batch size, the number of epochs, the learning rate, and the learning rate schedule. We set the model architectures of SIMs as almost the same size as that of the compared monDEQs. We also applied the training algorithm of the monDEQs to SIMs: the learning rate schedules  are also implemented in the same way as the monDEQs. As in the monDEQs, we minimized the cross-entropy loss for SIMs.
Appendix B.3. Image Regression Table B.12 shows the statistics of a dataset in the image regression task. All the images have 512 × 512 pixels, and we used equally spaced 1/4 (65,536) pixels from an image as training data, another 1/4 pixels as validation data, and other 1/4 pixels as test data. It should be noted that all the images are the same size and thus have the same statistics. Table B.13 shows the architecture of µ NN in SIMs. We used a neural network that consists of two fully connected layers and one additional linear layer to fit the dimension of the lifted space. Table B.14 shows the configuration of SIMs. In this task, we tuned the number of hidden units of µ NN , the first and second lift dimensions (i.e. N  and M ), the bandwidth of the RFF, and the learning rate. Following (Tancik et al., 2020), we applied Adam with the constant learning rate and performed the full batch gradient descent which uses all data points (pixels) for the update of parameters. #iterations denotes the number of steps in the gradient descent and was set to 2000 as in (Tancik et al., 2020). We minimized the mean squared loss for all models. The difference from (Tancik et al., 2020) is that we did not apply the sigmoid function before the output for each model: we observed such a setting improved the performance.