Year: 2021

Abstract:
Implicit layers are computational modules that output the solution to some problem depending on the input and the layer parameters. The deep equilibrium model (DEQ) outputs a solution to a fixed point equation. On the other hand, deep declarative networks (DDNs) solve an optimisation problem in their forward pass, an arguably more intuitive, interpretable problem than finding a fixed point. We show that using a kernelised generalised linear model (kGLM) as an inner problem in a DDN yields a large class of commonly used DEQ architectures with a closed-form expression for the hidden layer parameters in terms of the kernel. The activation functions have interpretations in terms of the derivative of the log partition function and the Bayesian prior or kernel regulariser. Building on existing literature, we interpret DEQs as fine-tuned, unrolled classical algorithms, giving an intuitive justification for why DEQ models are sensible. We use our theoretical result to devise an initialisation scheme for DEQs that allows them to solve kGLMs in their forward pass at initialisation. We empirically show that this initialisation scheme improves training stability and performance over random initialisation.

heading: 1 INTRODUCTION
text: Implicit layers (Pineda, 1987; Almeida, 1990) have recently been subject to renewed attention (Kolter et al., 2020). In contrast with explicitly defined layers, implicit layers define a mapping in terms of a solution to some problem depending on the input and problem parameters. For example, deep equilibrium models (DEQs) consist of layers that output fixed points of parameterised functions (Bai et al., 2019). Deep declarative networks (DDNs) use declarative nodes which output solutions to optimisation problems (Gould et al., 2016; Amos & Kolter, 2017). Neural ordinary differential equations (NODEs) output solutions to ODEs (Chen et al., 2018; Dupont et al., 2019). Traditional explicit layers can always be represented as implicit layers (for example, see Proposition 4.10 of Gould et al. (2021)). Also, solutions to certain convex optimisation problems may be obtained via an iterative optimisation procedure such as Newton’s method or gradient descent, and as such, may be represented as fixed points of an iterative scheme. A correspondence between DDNs and DEQs is expected (but undiscovered), given the fundamental connection between fixed points of iterative maps and critical points of optimisation problems (Ryu & Boyd, 2016). This leads to two natural questions: (1) For a given optimisation problem, what is the corresponding DEQ architecture? (2) Can this correspondence be exploited for theoretical, conceptual, or practical benefit?
Contributions. (1) We prove an equivalence between a DEQ and a DDN with a classical statistical model — a kernelised generalised linear model (kGLM) — as the declarative node, as illustrated in Figure 1 and formally stated in Theorem 3 and Corollary 6. The weights of the DEQ layer have closed-form expressions in terms of the kernel. The surprise in our result is that the feature mapping involved in this correspondence is exactly the class of hidden layers that are most commonly used in practice. (2) We empirically demonstrate that initialising a DEQ as a DDN using the derived expression for the weights improves performance and stability over random parameter initialisation.
Notation. Let X ⊆ Rdx and Y ⊆ R be coordinate and target spaces. Define a Y−valued stochastic process {y(x, ω)}x∈X indexed by x ∈ X with outcome ω from a sample space Ω. LetX ∈ Rn×dx be a matrix such that the ith row ofX is some element xi ∈ Rdx . Similarly define X̃ ∈ Rñ×dx . Take n evaluations {yi}ni=1 = {y(xi, ω0)}ni=1 from a realisation of the stochastic process for every ith row in X , and form a corresponding matrix Y . Using the same realisation of the stochastic process, take
ñ evaluations {ỹi}ñi=1 = {y(x̃i, ω0)}ñj=1 for every ith row in X̃ . Call X and Y the inner coordinate and target data, and X̃ and Ỹ the outer coordinate and target data. X and X̃ encode the positions of corresponding elements in Y and Ỹ , and induce kernel matrices whose utility roughly corresponds with positional embeddings (Gehring et al., 2017) in transformers (Vaswani et al., 2017).
Example implication of result. Consider image inpainting a single 32 × 32 colour image. Take X = {1, . . . , 32}2 × {1, 2, 3}, and X ∈ Zn×3 to be a matrix consisting of (non-repeated) triples, where n < 322. Take X̃ ∈ Zñ×3 to be the matrix consisting of remaining triples, such that n+ ñ = 322 × 3. An incomplete image Y ∈ Rn together with its missing values Ỹ ∈ Rñ are jointly sampled from the same distribution from which CIFAR10 images are drawn. We may use kernel ridge regression or a Gaussian process with Gaussian likelihood to produce a per-pixel prediction of Ỹ given X,Y and X̃ . These models represent special cases of regularised or Bayesian kGLMs with closed-form predictors, but more generally an algorithm is required to compute the predictor— see Appendix A. We now move from a single image to a set of images. We wrap the kGLM with an outer minimisation loop over the loss L between a neural network Fγ output applied to the kGLM predictor and target Ỹ to obtain a DDN, as sketched in Figure 1. We view the DDN as producing a per-image prediction. We show that such a DDN is equivalent to a DEQ with a closed-form expression for the parameters depending on the kernel. When the kGLM has fixed hyperparameters, the DEQ has a fixed set of parameters (made precise in Theorem 3 and Corollary 6). Initialising the DEQ with these parameters improves training stability and performance, as shown in § 4. Another example involving kernel logistic regression for MNIST image segmentation is given in Appendix G.

heading: 2 DEEP EQUILIBRIUM MODELS
text: Given some Y ∈ Rn×dy containing row elements y ∈ Y, deep equilibrium models (DEQs) (Bai et al., 2019) compute an embedding z∗ ∈ Z ⊂ Rn×dz as a fixed point of some function and then output an affine transformation f∗(Y ) of the fixed point. More concretely, letting gθ(·;Y ) : Rn×dz → Rn×dz denote a parameterised function (e.g. a neural network) with parameters θ ∈ Θ,
f∗(Y ) = V z∗, z∗ = gθ(z ∗;Y ), z∗ ∈ Z ⊂ Rn×dz , (1)
where V ∈ Rdu×(n×dz) are referred to as readout parameters, and V may be thought of as a du×ndz matrix acting on a flattened vector with ndz entries. We will find it more natural to introduce a linear skip connection between Y and the output. We will also choose dy = dz = 1. That is, we consider
f∗(Y ) = V1z ∗ + V2Y, z ∗ = gθ(z ∗;Y ), z∗ ∈ Z ⊂ Rn (2)
Typically DEQs are trained in a supervised manner by running some iterative algorithm such as gradient descent on a loss that aims to match a neural network involving f∗(Y ) to some target Ỹ with respect to the DEQ and neural network parameters. Given a dataset {(Ys, Ỹs)}Ns=1, this procedure may be realised by running an iterative algorithm over the equality constrained optimisation problem
min γ,θ,V1,V2
N∑ s=1 L ( Ỹs,Fγ ( f∗(Ys) )) subject to z∗s = gθ(z ∗ s ;Ys), f ∗(Ys) = V1z ∗ s + V2Ys, s = 1, . . . , N
(3)
as if the algorithm were to solve the problem, even though it may be highly non-convex. Here L is some loss function and Fγ is a neural network with parameters γ. A large class of iterative algorithms requires the derivative with respect to the network parameters θ and V1, V2; these may be computed under mild conditions without backpropagating through the fixed point solution via implicit differentiation. We refer the reader to Bai et al. (2019) for details.
Dealing with existence and uniqueness of fixed points. In order for (1) and (2) to specify useful computational rules, they have to admit at least one fixed point. To avoid having to choose which fixed point should be returned with additional rules, it might be desirable that they admit exactly one fixed point. We give example configurations of fixed points in Figure 2. The classical Banach fixed point theorem (BFPT) gives sufficient conditions for the existence and uniqueness of fixed points. Theorem 1 (BFPT). Let (Z, d) be a complete metric space. A function H : Z → Z is said to be a contraction if there exists some 0 ≤ q < 1 such that for all z, z′ ∈ Z , d (H(z), H(z′)) ≤ qd(z, z′). Every contraction admits a unique fixed point.
The contraction property is suggestive of a simple algorithm (Hasselblatt & Katok, 2003) for finding the fixed point z∗ of a contraction H given an initial guess z0: While some termination condition is not met, Update the current estimate zr for the fixed point to beH(zr−1). This algorithm allows one to interpret the inner problem of (3) as an infinitely deep neural network with weights shared between each layer. More advanced algorithms are available, which we do not discuss here.
H is a contraction if the induced matrix norm ‖ · ‖ of the Jacobian of H is strictly less than 1 on Z , providing a useful test for checking whether a unique fixed point of H exists. More generally,
Proposition 2. Let Z be an open strictly convex set, Z its closure, H : Z → Z differentiable on Z and continuous on Z . If ‖DH‖ ≤ λ < 1 on Z , then H is a contraction on Z .
We provide a proof in Appendix B. While elegant, contractions are stronger than necessary and using them limits the set of admissible networks. It is therefore desirable to work with tighter conditions, or otherwise sidestep the issue. A consensus on how to deal with existence and uniqueness does not yet appear to have been reached. One may run the algorithm as if a unique fixed point exists (Bai et al., 2019), constrain or monitor the parameters during learning to ensure that a unique fixed point exists, add a penalisation term on an estimate of the spectral norm of the Jacobian (Roosta-Khorasani & Ascher, 2015) to encourage a contractive map (Bai et al., 2021), or design variants of DEQs (through appropriate restriction) that explicitly ensure that a unique fixed point exists (Winston & Kolter, 2020). Other conditions are given by El Ghaoui et al. (2021). We do not address this problem; we assume (in a precise sense, see Assumption 3) verifiable conditions for the BFPT.
Initialisation. For the vast majority of practical neural network architectures, learning involves applying an optimisation procedure (e.g. stochastic gradient descent or quasi-Newton’s method) to an empirical risk minimisation problem in the overparameterised network’s parameters outside of conditions that guarantee convergence to a minimum. Debate continues as to why such procedures lead to good generalisation performance (Belkin et al., 2019; Zhang et al., 2021). However, there are intuitive, empirical and theoretical reasons to suggest that appropriate parameter initialisation plays an important role (Glorot & Bengio, 2010; He et al., 2015; Poole et al., 2016; Hu et al., 2020). For DEQs, initialisations have not yet been studied in great detail. By interpreting DEQs as unrolled kGLMs, we find parameter initialisations as a corollary of our main theoretical result. See Figure 3.

heading: 3 EXACT FIXED PARAMETER EQUIVALENCE BETWEEN DDNS AND DEQS
text: We present our main result here and give a special case of a tractable Ising model in Appendix C.2. Our main result is a derivation of a fully connected or convolutional DEQ architecture as a model that is equivalent to a DDN that solves a kGLM as its inner problem.

heading: 3.1 SETTING: A KGLM INSIDE A DDN
text: Inner problem. We work with the exponential family in canonical form (Deisenroth et al., 2020, Section 6.6), with log partition function A : R → R and sufficient statistic T : R → R. We choose a strongly log-concave (Saumard & Wellner, 2014) prior pK/λ over evaluations f(X) of a kGLM,
pK/λ (f(X)) = q (f(X))φK/λ (f(X)) , φK/λ(v) ∝ e− λ 2 v >K−1v, (4)
where q is a log-concave function and φK/λ is a Gaussian density with covariance K/λ. The MAP f∗ of the kGLM with representer coefficients α∗ satisfies
α∗ = arg min α∈Rn −α>KT (Y ) + 1>A(Kα)− log pK/λ (Kα) , f∗(·;X,Y ) = k(·, X)α∗. (5)
For convenience, we define ρ :≡ (− log q)′◦(A′)−1, whenever it exists. When q is constant, pK/λ is Gaussian and ρ ≡ 0 if (A′)−1 exists. Note (5) is strongly convex in α and admits a unique solution.
Deep declarative network. We will use (5) as a declarative node (Gould et al., 2021) by wrapping it with an outer optimisation task in order to learn the parameters of a neural network by empirical risk minimisation. Let L denote a loss function (for example, mean squared error), and let Fγ define a generic neural network with parameters γ. Form a set of N quartets {(Xs, Ys, X̃s, Ỹs)}Ns=1, for example a set of N images. Define Ks = k(Xs, Xs). We consider the optimisation task
min γ
N∑ s=1 L ( Ỹs,Fγ ( f∗(X̃s;Xs, Ys) )) subject to α∗s = arg min
α∈Rn −α>KsT (Ys) + 1>A(Ksα)− log pK/λ (Ksα) , s = 1, . . . , N,
f∗s (·) = k(·, Xs)α∗s , s = 1, . . . , N, (6)
where we have used the strong convexity of the N inner problems to ensure a unique solution of the N inner problems. The outer problem is in general non-convex and highly nonlinear in γ.

heading: 3.2 MAIN RESULT: DERIVING A FIXED-PARAMETER DEQ
text: Restricting ourselves to strictly convex log partition functions or requiring that the prior (4) is Gaussian (i.e. q ≡ 1), z∗ may be represented as a solution to a fixed point iteration, whose form may be understood to be an embedding produced by a DEQ. Note that the log partition function A is always convex, and is strictly convex for so-called minimal exponential families.
Assumption 1. The log partition function A is strictly convex or q ≡ 1.
Under only Assumption 1, the equivalence between the embeddings z∗ of kGLMs and DEQs holds. Under Assumption 2, we may invert the mapping to obtain α∗ from the embeddings z∗, so that the DEQ and the kGLM are equivalent predictive models.
Assumption 2. The log partition function A is strictly convex and the kernel matrix K is strictly positive definite.
We now state our main result, the proof of which is given in Appendix D.
Theorem 3. Let α∗ be the minimiser of one of the inner optimisation tasks (5), inducing a function f∗(·;X,Y ) = k(·, X)α∗. Define z∗ := A′ (Kα∗). Under Assumption 1, z∗ is a fixed point
z∗ = σ ( W1z ∗ +W2T (Y ) +W3ρ(z ∗) ) , (7)
whereW1 = W3 = −λ−1K andW2 = λ−1K are parameter matrices, σ ≡ A′ and ρ ≡ (− log q)′◦ (A′)−1 are monotone non-decreasing. Furthermore, under Assumption 2, for any test index X̃ ∈ Rñ×dx and for at least one z∗ satisfying (7),
f∗(X̃;X,Y ) = V1 (z ∗ + ρ(z∗)) + V2T (Y ), (8)
where V1 = −λ−1k(X̃,X) and V2 = λ−1k(X̃,X).
Further conditions (Assumption 3) force the fixed point iteration to have a unique solution.
Assumption 3. Let Z be an open strictly convex set and Z its closure. Suppose A′ : R → Z is continuous on Z and A′ is differentiable on Z . Suppose
(A3.1) sup z∈Z |A′′(z)| < 1 and (A3.2) ‖K/λ‖ sup z∈Z ‖I + diagρ′(z)‖ ≤ 1.
Alternatively, swap the < in A3.1 with the ≤ in A3.2.
Note that when ρ is constant, i.e. p is Gaussian, A3.2 reduces to ‖K/λ‖ ≤ 1. Proposition 4. Under Assumptions 1 and 3, (7) admits a unique fixed point on Z .
The outer problem depends on the solution α∗s to the inner problem only through f ∗(X̃s;Xs, Ys), which under Assumption 2 in turn only depends on z∗s . We therefore have the following.
Corollary 5. Under Assumptions 1, 2 and 3, the optimisation problem (6) is equivalent to the constrained optimisation problem
min γ
N∑ s=1 L ( Ỹs,Fγ ( f∗(X̃s;Xs, Ys) )) subject to z∗s = σ ( W1sz ∗ s +W2sT (Ys) +W3sρ(z ∗ s ) ) , s = 1, . . . , N,
f∗(X̃s;Xs, Ys) = V1s ( z∗s + ρ(z ∗ s ) ) + V2sT (Ys), s = 1, . . . , N,
(9)
where V1s = −λ−1k(X̃s, Xs), V2s = λ−1k(X̃s, Xs), W1s = W3s = −λ−1Ks and W2s = λ−1Ks take the role of DEQ parameters and σ = A′ and ρ are monotone non-decreasing and take the role of DEQ activation functions. The fixed point condition for z∗s is met by exactly one element of Z .
The model in Corollary 5 allows the parameters W1s,W2s, V1s and V2s to vary with s whereas (3) does not. An additional assumption ensures these parameters are constant in s.
Corollary 6. In the same setting as Corollary 5, if Xs = X and X̃s = X̃ are constant in s, the parameters W1,W2,W3, V1, V2 are constant in s and we may write f∗(Ys) ≡ f∗s (X̃s;Xs, Ys).
Note that in (3), the outer problem is with respect to γ, θ, V1, V2 whereas in (9) the outer problem is only with respect to γ. The parameters V1, V2 and W1,W2 (which take the role of θ), are automatically determined according to X , X̃ and the kGLM kernel k.
Corollary 5 is suggestive of a DEQ model that we are yet to implement whose weights are a dynamic kernel function of the coordinate Xs. One parameterisation of such a model is to have an auxillary network that predicts a finite-dimensional feature mapping of the coordinates Xs.

heading: 3.3 REMARKS ON MATCHING DDNS AND DEQS
text: Special cases. Choosing pK/λ to be Gaussian, ρ ≡ 0. Additionally, if T is the identity and σ ≡ A′ ≡ tanh, we obtain a fully connected layer with tanh activations, as derived in Appendix C.2. Logistic sigmoidal activations are obtained when the kGLM is specialised as kernel logistic regression. When A′ and T are the identity, the kGLM becomes kernel ridge regression, with Gaussian process regression as the corresponding Bayesian model. Certain infinitely wide (nonlinear) neural networks are Gaussian processes (Neal, 1995). Our result shows that one may also construct Gaussian processes whose posterior predictive mean is an infinitely deep neural network with linear activations. The role of the coordinates X and X̃ differs between the two models.
Symbol matching. Theorem 3 says that an embedding (2) computed by a DEQ given Y is the same as an embedding (7) computed by a kGLM on coordinates X for a training set (X,Y ) if the hidden parameters W1, W2 and W3 are scalar multiples of the kernel matrix k(X,X).
Initialisation via optimisation warm-start. In a stronger setting, Corollary 6 shows a connection between the two optimisation tasks (3) and (9). Given a fixed DEQ architecture with a fixed point condition that may be expressed as (7) (and noting that this equation also includes convolutional architectures), we may initialise parameter values that ensure that the prediction is equivalent to a trained kGLM and identify the hyperparameters of the kGLM. This initialisation may be naive (i.e. assume a “reasonable” kernel function, coordinates and regularisation parameter for a wide range of tasks) or more informed (i.e. leverage some information about the task to suit the kernel and regularisation to the task being performed). Note that even without Assumptions 2 and 3, we may initialise the DEQ such that the feature embeddings z∗ are equivalent for at least one of the fixed points. We demonstrate both naive and informed approaches in § 4.
On augmentation and assumptions. Assumption 3 allows one to compute without concern using (7) since a unique fixed point is guaranteed to exist. DEQs are sometimes implemented without checking that a unique fixed point exists. This assumption may be replaced by any other assumption that is sufficient to ensure a unique fixed point. Using the derivative test on the open set Z instead of Z allows us to handle important edge cases, namely when A′ ≡ tanh (with derivative 1 at the origin) or A′ is Leaky ReLU with gradients in {0, 1} on both sides of the origin (with undefined
derivative at the origin). Assumption 2 allows one to map the kGLM embedding to the prediction. The skip connection present in (2) but not (1) is not a severe limitation from a practical perspective. The augmented model considers a coordinate X absent from the original DEQ. Note that any DEQ may be written as an augmented model by simply ignoring the coordinate X . Under only Assumption 1 and no others, the equivalence between the embeddings produced by kGLMs and DEQs is still valid for at least one of their respective fixed points.
Convolutional layers. The derived network architecture involves only matrix multiplication, however as is well-known, we may represent a convolution through a sparse matrix equation (Figure 4). Conveniently, the local structure expected in images may be encoded by a particular informed kernel choice in order to recover a convolutional architecture. We formally discuss how this may be done in 2D RGB images in Corollary 11, Appendix E through a particular choice of kernel k.
Other priors. One may replace the strongly log-concave prior pK/λ with, for example, the conjugate prior of the exponential family (which is only log-concave). Here we focus on the strongly log-concave prior to recover a form more similar to DEQs that are implemented in practice.
Spectral normalisation. Supposing ‖ · ‖ is the spectral norm, Assumption 3 is suggestive of a particular weight normalisation scheme. More precisely, if λ is the spectral norm of K, then A3.2 holds automatically. The technique of using layers with unit spectral norm is explored by Miyato et al. (2018) in the context of improving training stability of GANs. Our results indicate that such spectral normalisation also leads to improved stability in DEQs at initialisation, in the sense that it guarantees a unique fixed point by Proposition 4. Abusing terminology, we refer to the spectral norm of the linear transformation corresponding with a convolutional layer as the spectral norm of the convolutional layer. An easy to implement method for calculating the spectral norm of convolutional layers is available (Sedghi et al., 2018), which we use in our implementation and experiments.

heading: 4 EXPERIMENTS
text: DEQ parameters may be partitioned into two sets: the parameters of the layer that defines the fixed point iteration (θ = (W1,W2,W3, V1, V2)), and the remainder of the parameters. Leveraging our theoretical result, we give DEQs a warm-start by initialising θ such that the DEQ solves a kGLM.
We empirically demonstrate that such initialisation is superior to random intialisation. We stress that our comparison is between randomly initialised DEQs and kGLM initialised DEQs only. We are not interested in demonstrating that DEQs can achieve results competitive with state-of-the-art models, as this has already been demonstrated in previous studies (Bai et al., 2019; 2021). Our limited scope allows us to perform a more definitive empirical study with meaningful statistics. We give a detailed explanation of two of the tasks in the main text, and refer the reader to Appendix F for more experiments and Appendix G for an empirical demonstration of Theorem 3. For random initialisation, we initialise weights from a symmetric uniform distribution (the Pytorch default).
Smooth sequence-to-sequence task. Let X = [−2π − 2, 2π] and define a function h(x) = e−x 2/10 sin(x) + e−(x+9) 2 . Let Xs = X and X̃s = X̃ be 100 points on a uniform
grid in [−2π, 2π] and [−2π − 2, 2π − 2] respectively. Define sequences Ys and Ỹs to be joint evaluations of the sth realisation of a Gaussian process with covariance function ktrue(x, x ′) = 0.1 exp −‖x−x ′‖2 2 exp ( − sin2 ‖x− x′‖ ) , having means h(X) and h(X̃) respectively.
For informed kGLM initialisation, we pretend that we know the kernel and coordinate space. We choose k ≡ ktrue and form the required kernel matrices by taking Xs = X and X̃s = X̃ . For naive GLM initialisation, pretending we do not know the true underlying kernel or coordinate space, we default to knaive = exp −‖x−x′‖2 2 and chooseXs = X to be a uniform grid over the interval [−2π, 2π]. We initialise the hidden parameters W1 and W2 according to Theorem 3 and initialise V1 and V2 according to the same rule as random initialisation. When training using naive GLM, we update only the readout parameters for the first 10 epochs.
We generate a training and test set of size N = 20, 000 and 2, 000. Each sequence is evaluated on an n = 100 dimensional uniform grid. We choose gθ(z∗) = tanh (W1z∗ +W2Y + b) and train for 400 epochs using Adam with default hyperparameters. We repeat this for 100 trials using seeds 0 to 99 and collect the test error curves for each trial. See Figure 5(a) and (b).
Image denoising. We consider a convolutional architecture as described in Appendix E, with σ ≡ ReLU, T the identity, and ρ ≡ 0 applied to an image denoising task using the CIFAR10 dataset. The coordinate space X is the space of all 2D pixel and channel triples (i, j, c) ∈ {1, . . . , 32}2×{1, 2, 3}. Elements Ys are CIFAR10 images corrupted by i.i.d. additiveN (0, 0.2) noise, clipped to take values between 0 and 1. Elements Ỹs are the corresponding uncorrupted CIFAR10 images.
For random initialisation, we sample the filter weights i.i.d.∼ N (0,Var). For kGLM initialisation, we let λ = ‖K‖/C, where ‖ · ‖ denotes the spectral norm and C > 0 (see Assumption 3). Spectral norms are calculated using the fast and exact method described by Sedghi et al. (2018). The kernel is constructed randomly from the squared exponential kernel and described in Appendix F.
We are interested in separating the effects of the choice of scaling and initialisation scheme. For random initialisation, the scale is chosen through the variance, inducing a convolutional layer with random spectral norm. For kGLM initialisation, the scale C is equal to the spectral norm. We try all variances in logspace(10−3, 1, 25) and all C in logspace(10−2, 10, 25), which results in a grid of spectral norms covering roughly the same space for both initialisation schemes. We try 100 different random seeds for each configuration, resulting in 2× 100× 25 = 5000 different models. In Figure 6, we plot individually after training for 0, . . . , 5 epochs the test MSE against the spectral norm for both initialisation schemes and all random seeds. kGLM initialisation offers improved test error, training stability and reduced variance in test error. The curves for kGLM initialised models adhere to the shape predicted by classical generalisation theory, where the spectral norm
measures model complexity. For kGLM initialisation, spectral norms of 100 appear to represent a critical point; smaller values lead to stable training and larger values lead to unstable training. This contrasts with randomly initialised models, most of which have increasing test error with increasing spectral norm in their later epochs. This means it is easy to select an appropriate scale for kGLM initialisation, but not for random initialisation. In Figure 3, we plot sample targets, inputs and outputs for random and kGLM initialised DEQs without any training, showing that samples using kGLM initialisation are visually more similar to the expected output than random initialisation. Figure 3 and epoch 0 of Figure 6 provide a sanity check of Theorem 3. More plots are given in Appendix F.

heading: 5 DISCUSSION AND CONCLUSION
text: Related work. A review on connections between other optimisation-based iterative procedures and neural networks generally is given by Monga et al. (2021). Repeated applications of weight-tied layers can sometimes be seen as solving model-based problems, if the weights are chosen appropriately. Since the derived network implicitly minimises a conceptually tractable and principled model, it possesses a certain degree of interpretability. However, only the minimisation problems and not the forward pass computations of the network are interpretable, and as such the network does not satisfy the requirements of interpretability in many settings (Rudin, 2019). Monga et al. (2021) also discuss how disabling weight-tying and allowing parameters to be fine-tuned, thereby escaping the optimisation-based iterative procedure, can sometimes lead to improved performance, perhaps at the cost of interpretability, generalisation performance and theoretical guarantees.
Notable models analysed under this framework include ISTA (Beck & Teboulle, 2009) for solving LASSO, which when unrolled and fine-tuned yields LISTA (Gregor & LeCun, 2010), and ADMM for compressive sensing (Boyd et al., 2011) with an unrolled version called ADMM-CSNet (Yang et al., 2018). Our work considers a different class of iterative models to those surveyed by Monga et al. (2021), and recovers exactly the DEQ model. In these works, the architectures, tasks, model classes and activations are specific to a particular setting. More recently, Ramsauer et al. (2020) introduce a class of networks whose predictions minimise a continuous analogue of a Hopfield energy. This can be used to justify transformer architectures.
Conclusion. By inserting a kGLM as an inner problem inside a DDN, we derived a DEQ with a convolutional or fully connected implicit layer and a fixed closed-form expression for the weights. Such a result gives intuitive justification for computing with DEQ models: DEQs can solve classical statistical problems in their forward pass. Using this theoretical connection also allowed us to initialise DEQ parameters and obtain performance improvements over random initialisation.
ETHICS STATEMENT
Our work is mostly theoretical in nature and is unlikely to directly cause harm, discrimination or privacy violations. However, as our work may be used in future in applications, we encourage users of our work to adhere to ICLR’s general ethical principles. We welcome open enquiry and scrutiny of the theoretical results presented in this paper.
REPRODUCIBILITY STATEMENT
We have included code and the instructions required to reproduce our results in the supplementary material. The code will be made publically available upon acceptance.

heading: A KERNELISED GENERALISED LINEAR MODELS
text: A.1 KERNEL TRICK AND REPRESENTER THEOREM
The predictions of some classical algorithms depend on the training input data X ∈ Rn×dx only through the matrix product XX>. A widely exploited technique in machine learning, the kernel trick, replaces X with ψ(X) ∈ Rn×dψ , where dψ may be large or in fact infinite. Here ψ : Rdx → Rdψ and with an abuse of notation, when applied to a matrix it is understood to apply individually to each row so that ψ(X) ∈ Rn×dψ . It is easy to change the resulting predictor since the type of ψ(X)ψ(X)> ∈ Rn×n matches that of XX>, and scales with n and not dψ . One example is in kernel ridge regression (Saunders et al., 1998). Following our abuse of notation, we will use k(X,X) ∈ Rn×n to denote the matrix with ijth entry k(xi, xj) := ψ(xi)>ψ(xj). Any function of the form k together with the matrix X characterises the matrix ψ(X)ψ(X)>. Such a function is called a kernel function, and is a valid choice if and only if it corresponds with an inner product in feature space or equivalently, if it is positive semi-definite (PSD). The space of functions Hk that the kernel k characterises is called a reproducing kernel Hilbert space (RKHS).
Let k be a positive semi-definite (PSD) function inducing a reproducing kernel Hilbert space (RKHS) Hk. We will use K = k(X,X) ∈ Rn×n to denote the matrix with ijth entry k(xi, xj). We will make use of the representer theorem, which for appropriately regularised empirical risk minimisation problems, allows one to solve optimisation problems in nonparametric function space. Theorem 7 (Representer Theorem (Schölkopf et al., 2002)). Denote by R : [0,∞) → R a strictly monotonic increasing function, by X a set and by L : (X × R2)n → R⋃{∞} an arbitrary loss function. LetHk be an RKHS with kernel k. Then
f∗ ∈ arg min f∈Hk
L ({( xi, yi, f(xi) )}n i=1 ) +R (‖f‖Hk) =⇒ f∗(x) = n∑ i=1 αik(xi, x).
A.2 KERNELISED GENERALISED LINEAR MODELS
We work with the exponential family in canonical form (Deisenroth et al., 2020, Section 6.6), as shown in equation 10. As has also been noted by others (Canu & Smola, 2006; Cawley et al., 2007), the representer theorem allows one to extend appropriately regularised GLMs (McCullagh & Nelder, 1989) to kGLMs by replacing the linear predictor with an arbitrary element of an RKHS. Here we demonstrate this fact.
Generalised linear models use an exponential family for the likelihood of observing response yi. The exponential family in canonical form is
p(yi | ηi) = h(yi) exp (ηiT (yi)−A(ηi)) , (10) where h : R → R, A : R → R, and T : R → R are functions that must be chosen such that the resulting density (mass function) integrates (sums) to 1. The log-partition function A is always convex. Here ηi are called the natural parameters of the distribution. We work exclusively with univariate responses yi and canonical parameters ηi, but GLMs are constructed more generally. For GLMs with a canonical link function, one chooses ηi = x>i β = f(xi) to be a linear function of the predictor variables xi with parameters β ∈ Rdx . It is possible to extend GLMs to accommodate a non-linear predictor f in an RKHSHk with a small modification. This modification can be motivated by either a Bayesian or frequentist perspective by placing a Gaussian prior over the function values or regularising respectively, and then solving the resulting functional optimisation problem over an RKHS using the representer theorem. Following the former philosophy, the maximum a posteriori (MAP) estimate f∗ of p ( f(X) | Y ) with a particular prior with parameter λ satisfies
f∗ = arg min f∈Hk n∑ j=1 − log h(yi) +A (f(xi))− f(xi)T (yi) + λ 2 ‖f‖2Hk ,
so that by the representer theorem,
f∗(·) = k(·, X)α∗ where α∗ = arg min α∈Rn
1 >A (Kα)− α>KT (Y ) + λ
2 α>Kα, K = k(X,X).
As an example, in binary kernel logistic regression, using the representer theorem we seek to minimise
− log p (f(X) | Y ) = −Y >Kα∗ + 1> log ( 1 + eKα ∗ ) + λ
2 α∗>Kα∗,
where K = K(X,X) ∈ Rn×n is the kernel matrix with pqth entry Kpq = k(xp, xq). The solution can be found by applying Newton’s method, resulting in an iterative fixed-point numerical scheme for α∗. Details of this method can be found in a number of references (Zhu & Hastie, 2005). A detailed example of the Ising model is given in Appendix C.

heading: B DERIVATIVE TEST
text: We will require the mean value theorem. Theorem 8 (Mean value theorem). Let S : [a, b] → Rn be continuous on [a, b] and differentiable on (a, b). Then there exists some t ∈ (a, b) such that
‖S(b)− S(a)‖ ≤ ∥∥∥∥ ddtS(t) ∥∥∥∥ (b− a). We now restate the result to be proven. We follow the approach in theorem 2.2.16 of Hasselblatt & Katok (2003). Proposition 2. Let Z be an open strictly convex set, Z its closure, H : Z → Z differentiable on Z and continuous on Z . If ‖DH‖ ≤ λ < 1 on Z , then H is a contraction on Z . Proof. Let z1, z2 ∈ Z and define c(t) := z1 + t(z2 − z1) for t ∈ [0, 1] and S(t) := H ( S(t) ) . Then S : (0, 1) → Z by strict convexity. By the mean value theorem, there exists some t ∈ (0, 1) such that
‖H(z2)−H(z1)‖ = ‖S(1)− S(0)‖ ≤ ∥∥∥∥ ddtS(t) ∥∥∥∥ (1− 0) = ∥∥∥∥ ddtS(t) ∥∥∥∥
= ∥∥∥∥DH (S(t)) ddtc(t) ∥∥∥∥
= ‖DH (S(t)) (z2 − z1)‖ ≤ ‖DH (S(t))‖ ‖z2 − z1‖ ≤ λ ‖z2 − z1‖ < ‖z2 − z1‖
and therefore H is a contraction mapping.

heading: C WARM-UP: A TRACTABLE ISING MODEL
text: We demonstrate our main result for a special Ising model case. The general case is given in § 3.
C.1 AS A KGLM
Let x and x′ be two dimensional indices on the finite lattice, x ∈ {1, . . . , a}×{1, . . . , b} for integers a, b ∈ Z. The Ising model prescribes that
p (Y | f(X)) = Z−1 exp ( t−1 (∑ pq Jpqyipyiq + ∑ p fxipyip )) ,
where Z is the partition function, some normalising constant ensuring that p (Y | f(X)) is a valid probability mass function. Here the parameters Jpq determine interactions between the spin states yip and yiq , fxip controls an interaction between the spin state yip and an external magnetic field at lattice position xip and t is a temperature parameter. We let fx = f(x) belong to an RKHSHk with PSD kernel k. With a standard Gaussian prior over f , we have that
log p (f(X) | Y ) = t−1 (∑ p,q Jp,qyipyiq + ∑ p f(xip)yip ) − logZ − λ/2‖f‖2Hk .
Supposing no interactions, Jpq = 0 and Z = 2 cosh (t−1f(xip)) (Nguyen et al., 2017). By Theorem 7, letting K = t−1k(X,X) denote the matrix with pqth entry t−1k(xip, xiq),
α∗ = arg min α −Y >Kα+ 1> log (2 cosh (Kα)) + tλ 2 α>Kα.
The problem is strictly convex so we may find the solution by applying Newton’s method. Using that the derivative of log (2 cosh(·)) is tanh(·), stationarity of the solution implies that
0 = −KY +K tanh (Kα∗) + tλKα∗ := F (α∗). (11) The Jacobian ∂F∂α∗ is given by
J(α∗) := ∂F
∂α∗ = K>DK + tλK, where D : = diag
( sech2(Kα∗) ) ,
so that the cth Newton iteration requires solving a linear system, J(α∗(c−1)) ( α∗(c) − α∗(c−1) ) = K ( Y − tanh (Kα∗(c−1))− tλα∗(c−1) ) .
We employ the shorthand M(α∗(c)) = M(c) for any matrix function of α ∗ (c). The Newton update is typically rearranged as an Iteratively Re-Weighted Least Squares (IRWLS) form by substituting α∗(c−1) = J −1 (c−1)J(c−1)α ∗ (c−1), yielding
J(c−1)α ∗ (c) = KD(c−1)z(c−1),
where z := Kα + D−1(Y − µ) and µ = tanh (Kα∗). Once α∗ has been approximated to some tolerance using this numerical scheme, the predictor f∗ may be evaluated at points X̃ ∈ Rñ×dx using the representer theorem,
f(X̃) = tK̃α∗ ∈ Rn∗ , K̃ := t−1k(X̃,X). (12)
C.2 AS A DEQ
Rearranging (11), and applying tanh to both sides of the equation, we obtain
tanh (Kα∗) = tanh
( 1
tλ K (Y − tanh (Kα∗))
) .
Choosing z∗ = tanh (Kα∗), W1 = −(tλ)−1K, W2 = (tλ)−1K, we have z∗ = tanh (W1z∗ +W2Y ) . (13)
Note that α∗ being a unique fixed point of (11) implies that z∗ is a fixed point of (13), but not that the fixed point of (13) is necessarily unique. One condition that guarantees the uniqueness of the fixed point of (13) is that the induced matrix norm of K is less than or equal to 1. Without this assumption, (13) is still true but when used as a computational model, the output will in general depend on the fixed-point solver and initial conditions passed to that solver (and even still, might never return unstable fixed points). Finally by assuming that K is invertible and using (12), we obtain
f(Y ) = tK̃α∗ = tK̃K−1 tanh−1 z∗ = 1
λ K̃ (Y − z∗) = V1z∗ + V2Y,
where V1 = −(λ)−1K̃ and V2 = (λ)−1K̃.

heading: D PROOF OF MAIN RESULT
text: Theorem 3. Let α∗ be the minimiser of one of the inner optimisation tasks (5), inducing a function f∗(·;X,Y ) = k(·, X)α∗. Define z∗ := A′ (Kα∗). Under Assumption 1, z∗ is a fixed point
z∗ = σ ( W1z ∗ +W2T (Y ) +W3ρ(z ∗) ) , (7)
whereW1 = W3 = −λ−1K andW2 = λ−1K are parameter matrices, σ ≡ A′ and ρ ≡ (− log q)′◦ (A′)−1 are monotone non-decreasing. Furthermore, under Assumption 2, for any test index X̃ ∈ Rñ×dx and for at least one z∗ satisfying (7),
f∗(X̃;X,Y ) = V1 (z ∗ + ρ(z∗)) + V2T (Y ), (8)
where V1 = −λ−1k(X̃,X) and V2 = λ−1k(X̃,X).
Proof. Suppose A is strictly convex, implying A′ = σ strictly increasing and invertible. Q := − log q(·) is a convex function, and its derivative Q′ is monotone nondecreasing. Also note that Q′ ◦ (A′)−1 is monotone nondecreasing. Let K = k(X,X). From (5), stationarity at the minimum and the representer theorem implies
0 = −KT (Y ) +KA′(Kα∗) +KQ′(Kα∗) + λKα∗
Kα∗ = 1
λ (KT (Y )−KA′(Kα∗)−KQ′(Kα∗)) A′(Kα∗) = A′ ( 1
λ (KT (Y )−KA′(Kα∗)−KQ′(Kα∗)) ) A′(Kα∗) = A′ ( 1
λ
( KT (Y )−KA′(Kα∗)−KQ′ ◦ (A′)−1 ◦A′(Kα∗) )) (14)
z∗ = σ ( W1z∗ +W2T (Y ) +W3ρ(z∗) ) . (15)
Note a subtlety in moving from (14) to (15). Since A is strictly convex, α∗ is the unique fixed point of (14). However, z∗ is only a (not necessarily unique) fixed point of (15). If A is not strictly convex but q ≡ 1, the above again holds Q′ ≡ 0 and without the need to invert A′. Suppose Assumption 2. Since K is invertible, and σ is invertible (since A is strictly convex), then again using the representer theorem, new predictions may be formed as
K̃α∗ = K̃K−1σ−1z∗ = 1
λ K̃ (T (Y )− z∗ − ρ(z∗)) = V1 (z∗ + ρ(z∗)) + V2T (Y ),
where K̃ = k(X̃,X).
Proposition 4. Under Assumptions 1 and 3, (7) admits a unique fixed point on Z .
Proof. Let H be the mapping defined by iteration (15), H(z∗) = σ ( W1z∗ +W2T (Y ) +W3ρ(z∗) ) .
H is a contraction since the induced matrix norm of the Jacobian of H satisfies ‖DH‖ = ∥∥(W1 +W3diagρ′(z∗))diagσ′(W1z∗ +W2T (Y ))∥∥
≤ ‖K/λ‖ sup z∈Z ‖I + diagρ′(z)‖ sup z∈Z |A′′(z)|
< 1.
The result follows from Theorem 1.

heading: E CONVOLUTION AS MATRIX MULTIPLICATION
text: Let w, h ∈ Z>1 denote the pixel width and height of a space of images with c channels. Let Xwh = {1, . . . , h} × {0, . . . , w}, Xwhc = Xwh × {1, . . . , c} (16)
denote the corresponding pixel and pixel-channel spaces. LetX ∈ R(whc)×c denote a matrix, where every row is one of the pixel coordinates (i, j, c) ∈ Xwhc. Choose Xs = X for some set of N images {Ys}Ns=1. The images Ys are of size whc, and contain pixel intensities for the corresponding pixel coordinates in X .
Let x1, x2 and x3 denote the first, second and third coordinates of an element x ∈ Xwhc. Let h1(x, x
′) = 1(|x1 − x′1| ≤ r1), h2(x, x
′) = 1(|x2 − x′2| ≤ r2), where 1 is the indicator function for some truncating parameters r1, r2 > 0. Let κ0 : Xwhc × Xwhc → R denote any stationary PSD kernel, for example κ0(x, x) = e− 1 2‖x−x
′‖22 . Note that these functions each define a PSD kernel. Using the fact that kernels are closed under multiplication, define
κ(x, x′) := h1(x, x ′)h2(x, x ′)κ0(x, x ′). (17)
Here h1 and h2 ensure that only close pixels have non-zero kernel values.
Note that the matrix κ(X,X) is sparse when r1 or r2 are small and has an off-diagonal structure. More concretely choosing κ(x, x′) = e− 1 20‖[x1,x2]−[x ′ 1,x ′ 2]‖ 2 21(‖x−x′‖∞ ≤ 5) where ‖·‖∞ denotes the supremum norm, we obtain a whc×whcmatrix. Each block of size wh is repeated horizontally. This matrix represents c convolutional filters of filter length 5× 2 + 1 = 11 in both directions. We plot the induced 322 × 322 kernel matrix when w = h = 32 and c = 1 in Figure 4, as well as the corresponding operation as a convolution. When c > 1, this matrix is simply repeated vertically and horizontally c times. For rows in the middle of the matrix, there are 112 non-zero entries, corresponding to the indices included in the convolution operation when centered on the pixel of the corresponding row.
The proof is tedious and mostly centres around the notation of appropriate reshaping of arrays and zero-padding. In order to state our result, we first need to define the neural network convolution operation. Definition 9 (Single-channel convolution). Fix X = Xwh (in the sense of (16)) and let the rows of X be equal to the elements of X. Let ∗1 denote the single-channel image convolution operation (in the sense of convolutional networks).
More concretely, C ∗1 Y applies a filter C ∈ R(2r1+1)×(2r2+1) to a flattened image Y ∈ Rwh as follows. First flatten C to be an element C̃ of R(2r1+1)(2r2+1) and zero-pad Y to be an element Ỹ of R(2r1+w)(2r2+h). Then define for each j ∈ {1, wh} a vector C̃j with the same dimensionality as Ỹ i, such that C̃j,[j:j+(2r1+1)(2r2+1)] = C̃ and all other elements of C̃j are zero. Then for each j ∈ {1, wh}, define
(C ∗1 Y )j = C̃>j Ỹ = (CmatỸ )j so that C ∗1 Y ∈ Rwh, where Cmat ∈ Rwh×(2r1+w)(2r2+h) is a matrix with jth row equal to C̃j . Definition 10 (Multi-channel convolution). Let ∗ denote the muti-channel image convolution operation (in the sense of convolutional networks).
More concretely, C ∗ Y applies a set of filters C ∈ Rc1×c2×(2r1+1)×(2r2+1) to a flattened image Y ∈ Rwhc2 as
(C ∗ Y )j = c2∑ r=1 Cj,r,:,: ∗1 Yr ∈ Rwh,
so that C ∗ Y ∈ Rc1wh, where Yr ∈ Rwh represents the rth channel of Y .
Generally, we have the following corollary of Theorem 3.
Corollary 11. In the same setting as Theorem 3, under Definition 9, let k ≡ κ (in the sense of (17)) and let f∗(·;Y,X) = κ(·, X)α∗ be a global minimiser of (5). Under Assumption 1,
z∗ = σ ( C1 ∗ z∗ + C2 ∗ T (Y ) + C3 ∗ ρ(z∗) ) , (18)
where z∗ = A′(κ(X,X)α∗). C1 = C3 and C2 are convolutional filters with pqth entry −λ−1κ0(X>i,[r1+1,r2+1], Xi,[p,q]) and λ −1κ0(X > i,[r1+1,r2+1]
, Xi,[p,q]) respectively. σ = A′ and ρ = (− log q)′ ◦ (A′)−1 are monotone non-decreasing. Furthermore,
1. Under Assumption 3, (7) admits a unique fixed point on Z .
2. Under Assumption 2, for any test index X̃ ∈ Rñ×dx and for at least one z∗ satisfying (7),
f∗(X̃;Y,X) = V1 (z∗ + ρ(z∗)) + V2T (Y ), (19)
where V1 = −λ−1κ(X̃,X) and V2 = λ−1κ(X̃,X).
Proof. We begin with the case c = 1. By Definition 9, we have that
C1 ∗ z∗ = Cmatz̃∗ ∈ Rwh,
where z̃∗ ∈ R(2r1+w)(2r2+h) is a zero-padded z∗ and Cmat ∈ Rwh×(2r1+w)(2r2+h). The jth row of Cmat is equal to C̃j , where C̃j,[j:j+(2r1+1)(2r2+1)] = C̃ and all other elements of C̃j are zero. Here C̃ is a flattened version of C1.
By removing the zero entries from z̃∗ and the corresponding columns from Cmat, we may write
Cmatz̃∗ = Wmatz∗
for some W1mat ∈ Rwh×wh. The same holds for W2 and W3. We may therefore write σ ( C1 ∗ z∗ + C2 ∗ T (Y ) + C3 ∗ ρ(z∗) ) = σ ( W1matz∗ +W2matT (Y ) +W3matρ(z∗) ) .
On the other hand, Theorem 3 says that σ ( W1z∗ +W2T (Y ) +W3ρ(z∗) ) = z∗.
Choosing W1mat = W1 = −λ−1κ(X,X) and similarly for W2,W3, we obtain (18). This implies that the filter C1 has pqth entry −λ−1κ0(X>i,[r1+1,r2+1], Xi,[p,q]). The rest of the corollary follows from Theorem 3.

heading: F DETAILED EXPERIMENTS
text: Our experiments may be reproduced using the code provided in the supplementary material by following the README.md file.
Choice of kernel for initialisation. To intialise a convolutional network, we need to choose r1, r2 and κ0 according to (17). We choose r1 = r2 = 3 and κ0(x, x′) =∑c(c+1)/2 r=1 gr(x3, x ′ 3)e − 1 2`2r ‖x[1,2]−x′[1,2]‖ 2
, where gr(a, b) = 1 if a == ar and b == br or a == br and b == ar and 0 otherwise. We sample the squared lengthscales from a uniform distribution between 0 and 4.
The cost of kGLM initialisation is minuscule, and often even faster than random initialisation.
For a fully connected layer with hidden and readout weight matrices W and V of size n × n and ñ× n respectively, we require n2 + nñ evaluations of the kernel function and then a normalisation by the spectral norm (the default option for this calculation for a Python library such as Numpy will be the LAPACK divide and conquer algorithm, which will have a worst case of O(n3), but in practice the computation will be very fast). In contrast, random numbers (the usual method for initialising neural networks) requires n2 + nñ random number generations, which depending on software platform, is usually done through calls to the inverse transform sampling method, perhaps using a stochastic collocation Monte Carlo sampler. Either way, empirically we find that kGLM initialisation is actually much faster than random Gaussian initialisation for fully connected layers. See Figure 8.
For convolutional layers, which are basically just extremely sparse large fully connected layers, the kernel method induces an additional overhead associated with the various reshaping operations required to put the elements of the kernel matrix in the correct position in the convolutional filter layer. The calculation of the spectral norm is done through the method of Sedghi et al. (2018) at a cost of O ( w2c2(c + logw)), where c is the number of channels and w is the maximum of the pixel width/height of the image. Our implementation is far from optimised, but we find that even with this overhead, the initialisation costs less than 1 second on a laptop PC. This is a tiny fraction of the cost typically associated with training the network, which depends on the problem but will usually be much less than 1%. See Figure 8.
Scaling up parameter counts. The number of rows of the square kernel matrix is equal to the number of rows in X . This means, for example, in the setting of CIFAR10 where X ∈ R32∗32∗3×3, there are roughly 9.4 ∗ 106 elements in the kernel matrix. Using a convolutional representation, the sparse kernel matrix may be represented as a convolutional layer of size (3, 3, 32, 32). Crucially, the parameter count scales quadratically with the number of input channels, which in this case is 3.
One way to experiment on larger models is to scale up the size of the input data. To this effect, we take a dataset containing 3799 hyperspectral (25 channel) images (HSI) of roads (Lu et al., 2020) of total size (3799, 25, 192, 384). We consider a sequence of denoising tasks when the input data is a
subset of the 25 hyperspectral channels. The resulting plots are given in Figures 9, 10, 11 12. The results are summarised in Table 1

heading: G EXAMPLE EMPIRICAL INSTANTIATION OF EQUIVALENCE
text: We take binarised MNIST and for each image, independently corrupt 30% of bits by randomly flipping them. We are interested in obtaining a predictive model that can segment the original binarised image Ỹs from its noisy corrupt input Ys, together with the constant coordinate locations Xs = X̃s = X that may be represented as 784× 2 coordinate indices. The index 1 ≤ s ≤ 60, 000 runs over different images in the training dataset. We evaluate model performance on an independent testing dataset of size 10, 000. We empirically compare five related models.
1. Kernel logistic regression (KLR). Since each element of Ys ∈ R784 is binary, a natural kGLM model to use is kernel logistic regression (KLR) (Green & Yandell, 1985; Zhu & Hastie, 2005). With respect to the exponential family, such a model prescribes that the derivative of the log partition function A′(·) = σ(·) = (1 + e−·)−1 is the logistic sigmoid and the sufficient statistic T is the identity. As is common in the literature, we employ a zero-mean Gaussian prior so that ρ ≡ 0. Assumption 1 is satisfied under this setting (both A is strictly convex and q ≡ 1). The only remaining choices are the kernel function k and the regularisation parameter λ. In order to ensure that Assumption 2 is satisfied, we choose k to be the sum of a squared exponential kernel with lengthscale 1 and a white noise kernel with variance 10−8. This forces the resulting kernel matrix K to be strictly positive definite. In order to satisfy Assumption 3, we need to ensure that ‖K/λ‖ ≤ 1 (assumption A3.1 is automatically satisfied due to the setting of A in KLR). Accordingly, we set λ = 2‖K‖. Note that we do not actually need to ensure that Assumptions 2 and Assumptions 3 are satisfied for Theorem 3 to hold, but we do need these assumptions so that all of our theoretical results hold.
The KLR model is trained individually using every datapoint in {Ys}60,000s=1 as a dataset and a manually implemented iterated reweighted least squares algorithm consisting of updates according to (exact) Newton’s method. We allow Newton’s method to run for 5 steps, at which point the norm of the difference between subsequent values of α is typically on the order of 10−14. The output of the KLR model is a kernelised predictor f∗ that has been trained using dataXs and Ys, evaluated onXs and converted to the maximum a posteriori of the Bernoulli distribution through 1 { σ(f∗) > 0.5 } .
2. DEQ with kernel initial parameters (DEQ-kernel). Theorem 3 implies an equivalent DEQ model with fixed parameters. We use this model, namely
z∗ = σ ( W1z ∗ +W2Ys ) ,
f∗(Y ) = V1z ∗ + V2Y,
where W1 = V1 = −λ−1K and W2 = V2 = λ−1K are parameter matrices and σ is the logistic sigmoid function. To convert the kernelised predictor f∗(Y ) to the mean of the Bernoulli distribution, we must pass it again through the inverse link function,
σ(f∗(Y )) = σ(V1z ∗ + V2Y ) (Greyscale prediction)
and the resulting maximum a posteriori prediction is determined by examining whether the mean of the Bernoulli distribution is greater than 0.5,
1 { σ(f∗(Y )) > 0.5 } . (Binarised prediction)
3. DEQ with kernel initialisation and fine-tuned parameters (DEQ-kernel-trained). We interpret the DEQ with initial parameters model as a neural network, and update the parameters W1,W2, V1 and V2 using Adam to minimise the squared loss between the greyscale prediction and binarised ground truth. We also allow biases in both the hidden and readout layers. We train for 5 epochs using a batch size of 100, leading to 5× (60000/100) = 3000 gradient updates. 4. DEQ with random initialisation (DEQ-random). We take the same neural network architecture as with kernel initialisation, but randomly initialse the parameters from a Gaussian distribution using the Pytorch default.
5. DEQ with random initialisation and fine-tuned parameters (DEQ-random-trained). We take the same neural network architecture as with kernel initialisation, but randomly initialse the parameters from a Gaussian distribution using the Pytorch default. We also allow biases and train using the same optimiser and loss.
Observations. Remarkably (but not surprisingly, due to Theorem 3), the binarised prediction agrees exactly with KLR (an entirely independent piece of software, that uses an entirely different solver)
over all but 11 of 60, 000 MNIST datapoints each containing 784 pixels. For those 11 where disagreements were observed, only 1 of 784 pixels were different. The very small differences arise due to the tolerance of the solvers and the way tied predictions (i.e. those whose logistic sigmoid output is 0.5) are handled. In all 11 cases, the DEQ output a value of 0.5 and KLR output some value less than but very close to 0.5. See Table 2 for posterior probabilities, and Figure 13 for the corresponding binarised predictions. Figure 14 shows the test error as a function of training epoch for DEQ-kernel-trained and DEQ-random-trained.

Review summaries:
The authors draw a resemblance between two classes of models, namely deep declarative networks (DDNs) and deep equilibrium models (DEQs), under mild assumptions. In particular, they show that solving an empirical risk minimization problem of a kernelized generalized linear model (kGLM) over a RKHS in a DDN admits a fixed point that, when viewed with an outer minimization, resembles the optimization of a simple DEQ with a fully connected implicit layer. They use their result to motivate an initialization scheme for DEQs that matches the solution of a trained kGLM.

This paper describes a particular DEQ formulation that is motivated by a kGLM-based deep declarative network (DDN). The MAP of a kGLM admits a closed form solution and is guaranteed to exist due to the convexity of a kGLM's log probability. The paper further discusses making the solution unique through a Lipschitz constraint. While strongly convexity and Lipschitz conditions for making fixed points unique is known, the main advantage of this kGLM approach seems to be being able to have a closed form expression for the solution, which is useful as an initialization for training DEQs. Experimental results validate that careful initialization can benefit training compared to random initialization, and a kGLM approach is useful in cases where the kernel can include domain knowledge.

The paper discusses a new perspective that establishes the connection between optimization-based layers (i.e., DDNs) and the fixed-point forward computations of the deep equilibrium networks (DEQs). In particular, the paper shows that a kGLM optimization layer, under certain regularity and contractivity assumptions, can be written as an (arguably) simpler form of deep equilibrium model layer. The paper provides solid proof and a thorough discussion of the implication of this connection, such as 1) how we might want to parameterize the matrices; and 2) how this may induce a good initialization scheme for (this particular kind of) DEQ models which improves training stability and convergence.

The paper discusses the relationship between deep declarative networks(DDNs) and deep equilibrium models(DEMs). The paper starts with a specific set of DDNs with a kernelized generalize linear model and shows that some DEMs are equivalent to the DDNs. The paper discusses the insight from such observation and potential usage of DDN inspired initialization for weights for implicit models.