summaries: The paper presents 3 main theoretical results: - An RNN with 40 unbounded-precision neurons is Turing-complete. - An RNN with 54 bounded-precision neurons and two stacks is Turing-complete. - An RNN with a finite number of bounded-precision neurons and no stacks can simulate a Turing machine with a bounded tape, where the maximum tape length is related to the number of RNN neurons.

summaries: This paper analyses the Turing completeness of recurrent neural networks, and proposes a memory-augmented RNN architecture reminiscent of existing stack RNNs.   Specifically:  First, a construction for simulating any Turing machine with an unbounded-precision RNN is presented, in which effectively the left and right sides of the Turing machine's tape (with respect to its head) are encoded as two 'stacks' in the RNN. This construction is analysed and requires less simulation time than that of Siegelmann and Sonntag [1] (3T instead of 4T+O(used tape length), where T is computation time of the simulated Turing machine). Additionally the authors analyse the exact precision of the neurons necessary (as function of the simulated Turing machine's computation length) for the simulation (as opposed to previous works, which simply note a requirement for infinite precision).  Second, noting that it is unavoidable that the precision must grow when simulating longer Turing machine computations (to maintain the contents of the tape), the authors propose augmenting the RNN with a dynamic memory module, such that the neurons themselves may have constant, bounded, precision. They theoretically analyse their proposed architecture and show how a Turing machine may be simulated in it.

summaries: This work proposes a dynamic-growing memory module for RNNs, which serves to simulate Turing machines of bounded precision.  First, the authors prove how to encode a Turing machine to an unbounded precision RNN. The encoding uses fractal encoding for symbols, like previous work, and it replaces each step of the Turing machine as 3 steps in the RNN. Leaning in the work of [14], they prove that there is a 40-neuron unbounded precision RNN to simulate any Turing machine (i.e., showing Turing-completeness). Next, they consider the previous unbounded-precision definitions with the memory module to obtain a bounded-precision RNN.The memory module is a stack of neurons, that requires two neurons in the RNN to control it: push and pop. Then it is proven that the proposed bounded-precision RNN utilizes 2 of these stack to simulate a Turing machine. Further, each stack is divided in groups of $p$ neurons, such that the most active neurons are in the RNN (near to each side of the memory head in the Turing machine). Leaning again in [14], they can show that there is a 54-neuron with $p$-precision RNN with 2 growing memory module, that can simulate any Turing machine. This proves Turing-completeness of the precision-bounded RNN. Finally, they remove the growing memory module and simulate the memory within the neurons of an RNN, showing that an infinite number of bounded-precision RNN can simulate any Turing machine. the work finishes with an interesting discussion of the results and limitations.

summaries: This theoretical work shows the existence of un-bounded precision RNNs for simulating turing machines, with growing/shrinking memory module. Although the paper seems quite interesting and well-written, I am afraid I don't have the right background to make an accurate judgment on the novelty and potential impact of the work. I share few thoughts below, similarly due to lack-of-knowledge and limited reviewing time I was not able to check the proofs.

abstractText: Previous works have proved that recurrent neural networks (RNNs) are Turingcomplete. However, in the proofs, the RNNs allow for neurons with unbounded precision, which is neither practical in implementation nor biologically plausible. To remove this assumption, we propose a dynamically growing memory module made of neurons of fixed precision. The memory module dynamically recruits new neurons when more memories are needed, and releases them when memories become irrelevant. We prove that a 54-neuron bounded-precision RNN with growing memory modules can simulate a Universal Turing Machine, with time complexity linear in the simulated machine’s time and independent of the memory size. The result is extendable to various other stack-augmented RNNs. Furthermore, we analyze the Turing completeness of both unbounded-precision and boundedprecision RNNs, revisiting and extending the theoretical foundations of RNNs.

authors:
Hava T Siegelmann
Eduardo D Sontag

title: On the computational power of neural nets

venue: Journal of computer and system sciences,

title: Computation beyond the turing limit

title: Neural networks and analog computation: beyond the Turing limit

venue: Springer Science & Business Media,

authors:
H Georg Kuhn
Heather Dickinson-Anson
Fred H Gage

title: Neurogenesis in the dentate gyrus of the adult rat: age-related decrease of neuronal progenitor proliferation

venue: Journal of Neuroscience,

authors:
Arturo Alvarez-Buylla
Jose Manuel Garcıa-Verdugo

title: Neurogenesis in adult subventricular zone

authors:
David Marr
David Willshaw
Bruce McNaughton

title: Simple memory: a theory for archicortex

venue: In From the Retina to the Neocortex,

authors:
MW Jung
BL McNaughton

title: Spatial selectivity of unit activity in the hippocampal granular layer. Hippocampus

authors:
Gerd Kempermann
H Georg Kuhn
Fred H Gage

title: Experience-induced neurogenesis in the senescent dentate gyrus

authors:
Georg Elias Müller
Alfons Pilzecker

title: Experimentelle beiträge zur lehre vom gedächtniss, volume 1

venue: JA Barth,

authors:
Edward Grefenstette
Karl Moritz Hermann
Mustafa Suleyman
Phil Blunsom

title: Learning to transduce with unbounded memory

venue: Advances in neural information processing systems,

authors:
Armand Joulin
Tomas Mikolov

title: Inferring algorithmic patterns with stack-augmented recurrent nets

authors:
Guo-Zheng Sun
C Lee Giles
Hsing-Hen Chen
Yee-Chun Lee

title: The neural network pushdown automaton: Model, stack and learning simulations

venue: arXiv preprint arXiv:1711.05738,

authors:
Ankur Arjun Mali
Alexander G Ororbia II
C Lee Giles

title: A neural state pushdown automata

venue: IEEE Transactions on Artificial Intelligence,

authors:
Ankur Mali
Alexander Ororbia
Daniel Kifer
Lee Giles

title: Recognizing long grammatical sequences using recurrent networks augmented with an external differentiable stack

venue: In International Conference on Grammatical Inference,

authors:
Joe Kilian
Hava T Siegelmann

title: The dynamic universality of sigmoidal neural networks

venue: Information and computation,

authors:
Alex Graves
Greg Wayne
Ivo Danihelka

title: Neural turing machines

venue: arXiv preprint arXiv:1410.5401,

authors:
Alex Graves
Greg Wayne
Malcolm Reynolds
Tim Harley
Ivo Danihelka
Agnieszka Grabska- Barwińska
Sergio Gómez Colmenarejo
Edward Grefenstette
Tiago Ramalho
John Agapiou

title: Hybrid computing using a neural network with dynamic external memory

authors:
Turlough Neary
Damien Woods

title: Four small universal turing machines

venue: Fundamenta Informaticae,

authors:
J Nicholas Hobbs
Hava Siegelmann

title: Implementation of universal computation via small recurrent finite precision neural networks

venue: In 2015 International Joint Conference on Neural Networks (IJCNN),

authors:
Fred C Hennie
Richard Edwin Stearns

title: Two-tape simulation of multitape turing machines

venue: Journal of the ACM (JACM),

title: Simple statistical gradient-following algorithms for connectionist reinforcement learning

venue: Machine learning,

heading: 1 Introduction

text:
Symbolic (such as Turing Machines) and sub-symbolic processing (such as adaptive neural networks) are two competing methods of representing and processing information, each with its own advantages. An ultimate way to combine symbolic and sub-symbolic capabilities is by enabling the running of algorithms on a neural substrate, which means a neural network that can simulate a Universal Turing Machine (UTM). Previous works [1, 2, 3] have shown that this is possible – there exists a recurrent neural network (RNN) that can simulate a UTM. These proofs assumed a couple of neurons with unbounded precision that equals the number of symbols used in the Turing tape. Here we provide an alternative simulation of a UTM by RNNs with bounded-precision neurons only.
The general idea works as follows. The Turing Machine’s tape is stored in a growing memory module, which is a stack of neurons with pushing and popping operations controlled by neurons in the RNN. The size of the growing memory module is determined by the usage of the Turing tape - it dynamically recruits new neurons when more memories are needed and releases them when memories become irrelevant. The neurons in the stack, except for the top neuron, are not regularly updated (and hence can be referred to as passive), saving computational cost for memories that are not in the focus of the computing and do not require change. Using growing memory modules, a 54-neuron bounded-precision RNN is constructed that can simulate any Turing Machine.
Our proposed growing memory modules are inspired by biological memory systems. The process of dynamically recruiting new neurons when more memories are necessary is also observed in biological memory systems. Neurogenesis is the process by which new neurons are produced in the central nervous system; it is most active during early development, but continues through life.
*Both authors contributed equally.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
In adult vertebrates, neurogenesis is known to occur in the dentate gyrus (DG) of the hippocampal formation [4] and the subventricular zone (SVZ) of the lateral ventricles [5]. Since DG is well-known in neuroscience for its role in pattern separation for memory encoding [6, 7], this suggests that biological memory systems also dynamically recruit new neurons. The rate of neurogenesis in adult mice has been shown to be higher if they are exposed to a wider variety of experiences [8]. This further suggests a role for self-regulated neurogenesis in scaling up the number of new memory that can be encoded and stored during one’s lifetime without catastrophic forgetting of previously consolidated memory. Besides the mechanism of recruiting new neurons, the process of storing neurons in growing memory modules also shares some similarities with biological memory consolidation, a process by which short-term memory is transformed into long-term memory [9, 10]. Compared to short-term memory, long-term memory is more long-lasting and robust to interference. This is similar to the neurons stored in growing memory modules - the values of these neurons (except the top neuron in the stack) remain unchanged and cannot be interfered by the RNN, providing a mechanism to store information stably.
Growing memory modules share similarities with other stack-augmented RNNs [11, 12, 13, 14, 15]. In neural stacks [11], the RNN outputs two continuous scalars that control the strength of pushing and popping operations, and the stack is made differentiable by adding a strength vector. In stack RNNs [12] and DiffStk-RNN [15], the RNN outputs the probability vector corresponding to pushing, popping, and no operation. In NNPDA [13], the RNN outputs a continuous scalar that controls the pushing and popping operations, with the minimum value corresponding to popping and the maximum value corresponding to pushing. NSPDA [14] uses a discrete-valued action neurons to control pushing and popping operations. In contrast to these models, growing memory modules have a simple design and do not need to be differentiable. However, it can be easily shown that growing memory modules can be simulated by these stack-augmented RNNs in linear time, and thus growing memory modules can be considered a generic type of stack-augmented RNNs. Therefore, our proof on the Turing completeness of an RNN with growing memory modules can be extended to stack-augmented RNNs in general, establishing their theoretical motivation.
A Turing-complete RNN that is fully differentiable was introduced in 1996 [16]; this feature is a prerequisite to have the network trainable by gradient descent. It was followed by the Neural Turing Machine (NTM) [17] and its improved version, the differentiable neural computer [18], which are both differentiable and trainable RNNs equipped with memory banks. Though inspired by Turing Machines, bounded-precision NTMs and differentiable neural computers are not Turing-complete due to the fixed-sized memory bank, but they can simulate space-bounded Turing Machines (see Section 5).
Simulation of a Turing Machine by an RNN with growing memory modules represents a practical and biologically inspired way to combine symbolic and sub-symbolic capabilities. All neurons in the RNN and growing memory modules have fixed precision. While the size of growing memory modules is linear in the number of symbols used in the Turing tape, the number of neurons in the RNN is still constant. Moreover, the neurons in growing memory modules (except the top neuron in the stack) are passive at most times. As a result, the time complexity of the simulation is linear in the simulated machine’s time and independent of the memory size. By showing how to simulate a Turing Machine with a bounded-precision RNN and thereby constructing a bounded-precision RNN that can run any algorithms, our paper proposes a practical method that combines symbolic and sub-symbolic capabilities.
The remainder of the paper is structured as follows. Section 2 describes the preliminary of the paper, including the definition of Turing Machines and RNNs. Section 3 revisits and extends theories relating to simulating a Turing Machine with unbounded-precision RNNs and shows the existence of a 40-neuron unbounded-precision RNN that is Turing-complete. Section 4 presents the growing memory modules and proves the existence of a 54-neuron bounded-precision RNN with two growing memory modules that is Turing-complete. Section 5 relates the number of neurons and the precision of RNNs when simulating Turing Machines. Section 6 concludes the paper.

heading: 2 Background and Notation

text:
A Turing Machine is a 7-tuple M = (Q,Σ,Γ, δ, q0, ♯, F ), where Q is a finite set of states, Σ is a finite set of input symbols, Γ is a finite set of tape symbols (note that Σ ⊂ Γ), δ : Q×Γ → Q×Γ×{L,R}
is the machine’s transition rule, q0 ∈ Q is the initial starting state, ♯ is the blank symbol (note that ♯ ∈ Γ but ♯ /∈ Σ), and F ⊂ Q is the set of final state. We only consider deterministic Turing Machines in this paper.
The instantaneous description (or the configuration) of a Turing Machine is typically defined as a tuple of state, tape, and the location of the read/write head. However, we use a slightly different definition in the paper. We define the left-tape symbols (or the left tape in short), denoted as sL, to be the string of symbols starting at the symbol under the read/write head and extending to the left. We define the right-tape symbols (or the right tape in short), denoted as sR, to be the string of symbols starting at the symbol to the right of the read/write head and extending to the right. The first symbol in both strings is the closest to the read/write head in both sL and sR (that is, the left tape is reversed in the representation) and the blank symbols at the two ends are omitted. Therefore, the length of sL and sR combined equals to the number of symbols used in the Turing tape in each step (that is, unbounded but not infinite). Since sL and sR encode both the tape and the location of the read/write head, we define the instantaneous description of a Turing Machine as a 3-tuple (q, sL, sR) ∈ (Q,Γ∗,Γ∗), where q denotes the state, sL denotes the left-tape symbols, and sR denotes the right-tape symbols. Though the two definitions are equivalent, this definition allows easy encoding of tape symbols into neurons’ values. The set of all possible instantaneous description is denoted as X := (Q,Γ∗,Γ∗). In each step of a Turing Machine, the symbol under the read/write head is read and, together with the state, determines the symbol to be written under the read/write head, the direction of moving the tape, and the next state. To be precise, the complete dynamic map of M, denoted as PM : X → X , is defined as follows: 1. Let x = (q, sL, sR) be the input configuration, sL,(1) denote the first symbol in sL and sR,(1) denote the first symbol in sR. The transition is defined by the 3-tuple (q′, y, d) = δ(q, sL,(1)); 2. Replace the state of the machine q with q′ and the first symbol in sL by y; 3. Move the symbol in sL,(1) to become the new sR,(1) if d = L, and move sR,(1) to become the new sL,(1) if d = R (if there are no symbols left in sL or sR for moving, append a blank symbol ♯ to it before moving). Denote the left-tape symbols and the right-tape symbols after 2. and 3. by s′L and s ′ R respectively. Then PM(x) = (q′, s′L, s′R) represents one transition of the Turing Machine M from one configuration to the next. The partial input-output function of M, denoted as P∗M : X → X , is defined by applying PM repeatedly until q ∈ F , and is undefined if it is not possible to have q ∈ F by applying PM repeatedly. A recurrent neural network (RNN) is a neural network consisting of n neurons. The value of neuron i at time t ∈ {1, 2, ...}, denoted as xi(t) ∈ Q (Q is the set of rational numbers), is computed by an affine transformation of the values of neurons in the previous state followed by an activation function σ, i.e. xi(t) = σ( ∑n j=1 wijxj(t− 1) + bi), where wij are the weights and bi is the bias; or in vector form:
x(t) = σ(Wx(t− 1) + b), (1)
where x(t) ∈ Qn, W ∈ Rn×n and b ∈ Rn. This defines a mapping TW,b : Qn → Qn which characterizes an RNN. For simplicity, we consider the saturated-linear function in this paper; that is:
σ(x) :=  0 if x < 0, x if 0 ≤ x ≤ 1, 1 if x > 1.
(2)
Thus, x(t) ∈ (Q ∩ [0, 1])n for all t > 0. We say that a neuron xi(t) has precision p in base b if for all t > 0, xi(t) can be expressed as∑p
l=1 a(l)(t)∏l
j=1 c(j)(t) for some strings a(t) ∈ {0, 1, ..., b}p and c(t) ∈ {1, ..., b}p.
For a string a, we use a(i) to denote the ith symbol in a and a(i:j) to denote the string a(i)a(i+1)...a(j). For a function f that maps from a set Y to a subset of Y, we use fn to denote the nth iterate of f where 1 ≤ n < ∞. For any two vectors x ∈ Rm,y ∈ Rn, we use x ⊕ y ∈ Rm+n to denote the concatenation of the two vectors. We use A∗ to denote all possible strings formed by elements from set A.
The notation is summarized in the table found in Appendix D.

heading: 3 Turing Completeness of Unbounded-Precision RNNs

text:
To simulate a Turing Machine M by an RNN, we first consider how to encode the instantaneous description (q, sL, sR) ∈ X by a vector of rational numbers. For the state q ∈ Q, we encode it with ⌈log2 |Q|⌉ binary values, denoted as ρ(q) : Q → {0, 1}⌈log2 |Q|⌉, with each possible combination of binary values representing a specific state.
Example. Let Q = {1, 2, 3, 4, 5, 6}. We can encode it by ⌈log2 |Q|⌉ = 3 binary values, with ρ(q)(1) = [0, 0, 0], ρ(q)(2) = [0, 0, 1], ρ(q)(3) = [0, 1, 0], ρ(q)(4) = [0, 1, 1], ρ(q)(5) = [1, 0, 0], and ρ(q)(6) = [1, 0, 1].
For the left tape sL ∈ Γ∗ and the right tape sR ∈ Γ∗, we use fractal encoding in two rational numbers, as recommended in [1, 2]. The fractal encoding, bearing similarity to Cantor sets, enables fast manipulation of the top symbols.
Without loss of generality, assume that the tape symbols Γ are encoded into numbers {1, 3, 5, ..., 2|Γ|− 1} and that the blank symbol ♯ is encoded by 1. Then, define the fractal encoding ρ(s) : Γ∗ → Q by:
ρ(s)(y) :=  |y|∑ i=1 y(i) (2|Γ|)i + 1 (2|Γ|)|y| · (2|Γ| − 1) . (3)
Example. Let Γ = {1, 3, 5, 7}, sL = (3, 5, 7, 3, 5), and ♯ = 1. Then ρ(s)(sL) = 38 + 5 82 + 7 83 + 3 84 +
5 85 + 1 86 + 1 87 + 1 88 + ... = 3 8 + 5 82 + 7 83 + 3 84 + 5 85 + 1 85·7 . The last term in (3) represents the infinite blank symbols of a tape.
This encoding requires the tape symbol neurons to have the same precision as the size of the active (non-blank) part of the tape in base 2|Γ|. As the tape in a Turing Machine has an unbounded size, it means that we require neurons with unbounded precision. This is different from infinite precision but still not applicable. Hence, we will discuss how to remove this unbounded-precision assumption in Section 4.
Finally, we encode the value of the top symbol in each tape neuron, sL,(1) and sR,(1), by a binary tuple ρ(r) : Γ → {0, 1}|Γ|−1:
ρ (r) i (y) := 1{y > 2i} (1 ≤ i ≤ |Γ| − 1). (4)
That is, the i coordinate of ρ(r)(y) is 1 if and only if the symbol y has a value larger than 2i, and is 0 otherwise.
Example. Let Γ = {1, 3, 5, 7}. Then ρ(r)(1) = [0, 0, 0], ρ(r)(3) = [1, 0, 0], ρ(r)(5) = [1, 1, 0], and ρ(r)(7) = [1, 1, 1].
Combining the above discussion, we define the encoding function of configurations ρ : X → Q2|Γ|+⌈log2 |Q|⌉+|Q||Γ|+5 by:
ρ(q, sL, sR) = ρ (q)(q)⊕ ρ(s)(sL)⊕ ρ(s)(sR)⊕ ρ(r)(sL,(1))⊕ ρ(r)(sR,(1))⊕ 0, (5)
where 0 is a zero vector of size |Q||Γ|+ 5. Let ρ−1 : ρ(X ) → X be the inverse function such that ρ−1(ρ(x)) = x for all x ∈ X (note that ρ is injective, i.e. ρ(x) = ρ(x′) implies x = x′); we call ρ−1 the decoder function.
It should be noted that ρ(q)(q) ⊕ ρ(s)(sL) ⊕ ρ(s)(sR) is sufficient to decode the instantaneous description. We include ρ(r)(sL,(1))⊕ ρ(r)(sR,(1))⊕ 0 only because it facilitates the full simulation on the RNN. This completes the construction of the encoding function.
Given an instantaneous description x ∈ X , we initialize the neurons in an RNN with values ρ(x). Then, it is possible to construct the parameters of RNN such that the update given by the RNN on these neurons is the same as the update given by the Turing Machine: Theorem 1. Given a Turing Machine M, there exists an injective function ρ : X → QN and an n-neuron unbounded-precision RNN TW,b : Qn → Qn, where n = 2|Γ|+ ⌈log2 |Q|⌉+ |Q||Γ|+ 5, such that for all instantaneous descriptions x ∈ X ,
ρ−1(T 3W,b(ρ(x))) = PM(x). (6)
Proof. A sketch of the proof is as follows. Neurons in the RNNs are grouped by their function. Tape neurons, initialized with ρ(s)(sL) and ρ(s)(sR), encode the tape in fractal encoding. Readout neurons, initialized with ρ(r)(sL,(1)) and ρ(r)(sR,(1)), encode the first symbol in the left and the right tape. State neurons, initialized with ρ(q)(q), encode the state. We need to update the values of these neurons to simulate one step of the Turing Machine. Three steps (or stages) of RNN are required to simulate one step of a Turing Machine. In the first stage, entry neurons, initialized with 0, compute the combination of the state and the symbol under the head. Since this combination fully determines the next transition, we use it to update the state neurons and the temporary tape neurons during stage two. Temporary tape neurons serve as a buffer for tape neurons when shifting the tape to the left or right. In stage three, we move the values from the temporary tape neurons to tape neurons to complete the update. The detailed proof can be found in Appendix A.
In other words, to simulate one step of a Turing Machine M, we can first encode the instantaneous description x by the encoder function ρ, apply the RNN three times T 3W,b, and decode the values back by ρ−1 to obtain PM(x), the instantaneous description after one step of the Turing Machine. Or equivalently, for any Turing Machine M, there exists an RNN such that every three steps of the RNN yield the same result as one step of the Turing Machine.
By applying Theorem 1 repeatedly, we simulate a Turing Machine with an RNN in a linear time. To be specific, the partial input-output function of an RNN, denoted as T ∗W,b : QN → QN , is defined by applying T 3W,b repeatedly until q ∈ F (where q is the state that the RNN simulates), and is undefined if it is not possible to have q ∈ F by applying T 3W,b repeatedly. Based on this definition and Theorem 1, it follows that:
Corollary 1.1. Given a Turing Machine M, there exists an injective function ρ : X → Qn and an n-neuron unbounded-precision RNN TW,b : Qn → Qn, where n = 2|Γ|+ ⌈log2 |Q|⌉+ |Q||Γ|+ 5, such that for all instantaneous descriptions x ∈ X , the following holds: If P∗M(x) is defined, then
ρ−1(T ∗W,b(ρ(x))) = P∗M(x), (7)
and if P∗M(x) is not defined, then T ∗W,b(ρ(x)) is also not defined. If P∗M(x) is defined and computed in T steps by M, then T ∗W,b(ρ(x)) is computed in 3T steps by the RNN.
Corollary 1.1 shares similarities with Theorem 1 in [1]. However, our theorem states that 3T , instead of 4T , is sufficient to simulate a Turing Machine. We also give the relationship between the number of neurons required by the RNN and the size of Q and Γ in the Turing Machine.
A small UTM with 6 states and 4 symbols, denoted by U6,4, was proposed [19] and can simulate any Turing Machine in time O(T 6), where T is the number of steps required by the Turing Machine (the one to be simulated) to compute the result. As U6,4 is also a Turing Machine, we apply Corollary 1.1 to simulate U6,4, leading to a Turing-complete RNN. Plugging in |Q| = 6 and |Γ| = 4, we obtain the following result:
Corollary 1.2. There exists a 40-neuron unbounded-precision RNN that can simulate any Turing Machine in O(T 6), where T is the number of steps required for the Turing Machine to compute the result.
It should be noted that [1] focused on capabilities and proved a Turing-complete RNN with 1058 neurons; [20] proposed a Turing-complete RNN with 52 neurons. Here we provide a plug-and-play formula to simulate any Turing Machine.

heading: 4 Turing Completeness of Bounded-Precision RNNs with Growing Memory

text:
In the following, we consider how to remove the assumption of unbounded neural precision (Section 3) without reducing computational capacity. If we assume all neurons to have precision bounded by p in base 2|Γ|, then the tape can be encoded by ⌈|sj |/p⌉ neurons using fractal encoding, where j ∈ {L,R}. We do this by encoding every p symbols of the tape into a single neuron. Since most of these neurons do not require updates, just like symbols far from the read/write head in a Turing tape, we propose to store them in a separate growing memory module organized in a stack-like manner:
Definition 1. A growing memory module is a stack of non-zero neurons with push and pop operations controlled by two neurons in an RNN, denoted as push neuron u(t) and pop neuron o(t), in the following way: for every step t after the RNN finished updating the values of neurons, (i) if u(t) > 0, then a new neuron with the value u(t) is pushed to the stack and u(t) is set to 0; (ii) if o(t) = 0 and the stack is not empty, then the top neuron is popped from the stack and o(t) is set to the value of the top neuron in the updated stack; (iii) if o(t) = 0 and the stack is empty, then o(t) is set to a default value c.
An RNN with a growing memory module is a mapping TW,b : (QN ,Q∗) → (QN ,Q∗) where the first element of the tuple corresponds to values of neurons in the RNN and the second element of the tuple corresponds to the stack of the growing memory module. We can equip an RNN with multiple growing memory modules, with each module having its own push and pop neurons controlled by the RNN. An RNN with two growing memory modules will be defined by a mapping TW,b : (QN ,Q∗,Q∗) → (QN ,Q∗,Q∗). We can view the growing memory module as a way of dynamically pointing to a sequence of non-zero neurons appended by one zero neuron. o(t) can be viewed as the pointer for the last non-zero neuron in the sequence, and u(t) can be viewed as the pointer for the zero neuron at the beginning of the sequence; see Figure 1.
With two growing memory modules (one for the left tape and one for the right tape), we can construct an RNN with bounded-precision neurons that can simulate any Turing Machine. We first describe how to encode the instantaneous description (q, sL, sR) ∈ X by a vector of rational numbers and two stacks, with which an RNN and its growing memory modules can be initialized. In the following discussion, we assume that all neurons have precision bounded by p ≥ 2 in base 2|Γ|; that is, each neuron can encode p symbols at most.
Both the state q and the top symbols sL,(1), sR,(1) are encoded with binary values as in Section 3. For the tape sj (j ∈ {L,R}), we define the fractal encoding ρ(s) : Γ∗ → Q by:
ρ(s)(y) := |y|∑ i=1 y(i) (2|Γ|)i , (8)
which is similar to (3) except for the encoding of infinite blank symbols. Then, we encode the tape sj into a stack of neurons as follows: First, encode the rightmost p symbols (the ones farthest from the read/write head) with ρt and push it to an empty stack, denoted as Mj . Then, encode the next
rightmost p symbols and push it to Mj again, and repeat until at least one and at most p symbols remain in the tape. Denote this encoding function for the tape as ρ(M) : Γ∗ → Q∗. The remaining symbols in the tape, denoted as sj,(1:h(|sj |)), where h(y) := ((y − 1) mod p) + 1, will be encoded with fractal encoding ρ(s) as well, but would appear in neurons inside the RNN.
The general idea is to let only the symbols closest to the read/write head (sj,(1:h(|sj |))) reside in the RNN. If the number of symbols residing in the RNN reaches 0 or p, then we pop from or push to the stack respectively to ensure that at least 1 and at most p symbols reside in the tape neurons. It is interesting to note that the kth neuron in the stack (from the top) requires at least kp steps of the Turing Machine before it may be updated, so the values of neurons near the bottom of the stack will not be changed for many steps. That is, the neurons in the stack, except for the top neurons, are passive.
Example. Let Γ = {1, 3, 5, 7}, sL = (3, 5, 7, 3, 5, 5, 3, 7) and p = 3. Then the number of symbols to remain in the RNN is 2 and they are encoded by ρ(s)(sL,(1:h(|sL|))) = 3 8 + 5 82 . The remaining six symbols are stored in the stack: ρ(M)(sL) = [ 78 + 3 82 + 5 83 , 5 8 + 3 82 + 7 83 ].
We encode in neuron h(|sj |), the number of symbols in sj (j ∈ {L,R}). In the above example, h(|sL|) = 2. These neurons help the RNN to know when pushing or popping operations are required. We use the encoding ρ(h) : {1, 2, ..., p} → Q, defined by:
ρ(h)(y) = y
p+ 1 . (9)
Together, define the encoding function ρ : X → (Q2|Γ|+⌈log2 |Q|⌉+|Q||Γ|+19,Q∗,Q∗) by:
ρ(q, sL, sR) =(ρ (q)(q)⊕ ρ(s)(sL,(1:h(|sL|)))⊕ ρ (s)(sR,(1:h(|sR|)))⊕ ρ (s)(sL,(h(|sL|)+1:h(|sL|)+p))⊕
ρ(s)(sR,(h(|sR|)+1:h(|sR|)+p))⊕ ρ (r)(sL,(1))⊕ ρ(r)(sR,(1))⊕ ρ(h)(h(|sL|))⊕
ρ(h)(h(|sR|))⊕ 0, ρ(M)(sL), ρ(M)(sR)), (10) where 0 is a zero vector of size |Q||Γ| + 15. The first element of the tuple is for initializing the neurons in the RNN, while the second and third element of the tuple is for initializing the two growing memory stack modules. All encoded values have precision of p. Similar to the previous section, ρ is injective and so we can define the decoder function ρ−1 : ρ(X ) → X . With the new encoding function ρ and growing memory modules, we can prove an alternative version to Theorem 1 that only requires bounded precision neurons: Theorem 2. Given a Turing Machine M, there exists an injective function ρ : X → (Qn,Q∗,Q∗) and an n-neuron p-precision (in base 2|Γ|) RNN with two growing memory modules TW,b : (Qn,Q∗,Q∗) → (Qn,Q∗,Q∗), where n = 2|Γ| + ⌈log2 |Q|⌉ + |Q||Γ| + 19 and p ≥ 2, such that for all instantaneous descriptions x ∈ X ,
ρ−1(T 3W,b(ρ(x))) = PM(x). (11)
Proof. The detailed proof is in Appendix B. To illustrate the construction of the RNN, the parameters for the neuron initialized with ρ(h)(h(|sL|)), called the left guard neuron and denoted as gL(t), will be described here.
We assume all neurons in the RNN are initialized with values from the encoder function ρ(x) at time t = 1. The guard neuron gL(t) encodes the number of left-tape symbols residing in the RNN. In three stages of an RNN, we need to update its value from gL(1) = h(|sL|)/(p + 1) to gL(4) = h(|s′L|)/(p+1), where s′L is the left tape after one step of the Turing Machine. First, notice that h(|s′L|) can be expressed as:
h(|s′L|) =  h(|sL|)− 1 if d = L and h(|sL|) ≥ 2, p if d = L and h(|sL|) = 1, h(|sL|) + 1 if d = R and h(|sL|) ≤ p− 1, 1 if d = R and h(|sL|) = p,
(12)
where d is the direction that the Turing Machine’s head is moved.
Example. Assume d = L and the size of the active symbols is h(|sL|) = 1, which means that the Turing Machine is moving left and there is only one symbol residing in the RNN. Since after the move, sL will have one symbol less and so the corresponding neuron sL(t) will encode no symbols. As a result, the top neuron of the left growing memory module would be popped out, and sL(t) would assume its value. This way the neuron sL(t) again encodes the active symbols (1 to p) of the left tape series, and h(|sL|) is set to p. (Alternatively, one may prove it directly by the definition of h and the fact that |s′L| = |sL| − 1.) An analogous process holds when the Turing Machine is moving right. We implement (12) with an RNN as follows. Define stage neurons as:
c1(t+ 1) = σ(1− c1(t)− c2(t)), (13) c2(t+ 1) = σ(c1(t)), (14)
with both neurons initialized to be zero. Define c(t) := [c1(t), c2(t), c3(t)] where c3(t) := 1 − c1(t) − c2(t), then c(1) = [0, 0, 1]; c(2) = [1, 0, 0]; c(3) = [0, 1, 0]. These stage neurons signal which one of the three stages that the RNN is in.
In the construction of the RNN, there exists a linear sum of neurons, denoted as d(t) = [dL(t), dR(t)], such that if the Turing Machine is moving left, d(1) = [0, 0]; d(2) = [1, 0]; d(3) = [0, 0]; and if the Turing Machine is moving right, d(1) = [0, 0]; d(2) = [0, 1]; d(3) = [0, 0]; this signals which direction the Turing Machine is moving to (formulas of d(t) appear in Appendix B).
Then consider the following update rule for the left guard neurons:
gL(t+ 1) = σ(gL(t) + (dR(t)− dL(t)− pg′L(t) + pg′′L(t))/(p+ 1)), (15) g′L(t+ 1) = σ((p+ 1)gL(t) + dR(t)− p− 2c2(t)− 2c3(t)), (16) g′′L(t+ 1) = σ(2− (p+ 1)gL(t)− dR(t)− 2c2(t)− 2c3(t)), (17)
where gL(1) = h(|sL|)/(p+1), g′L(1) = g′′L(1) = 0. It can be verified that gL(4) = h(|t′L|)/(p+1) as defined by (12), completing the proof for gL(t).
Example. Assume d = L and h(|sL|) = 1. Then gL(1) = gL(2) = 1/(p+ 1) and g′L(1) = g′′L(1) = g′L(2) = g ′′ L(2) = 0. On the second stage, gL(3) = σ(1/(p + 1) − 1/(p + 1)) = 0, g′L(3) = σ(1−p) = 0, and g′′L(3) = σ(2−1) = 1. On the third stage, gL(4) = σ(0+p/(p+1)) = p/(p+1) as required.
The full proof appears in Appendix B.
Similar to Corollary 1.1, it follows that: Corollary 2.1. Given a Turing Machine M, there exists an injective function ρ : X → (Qn,Q∗,Q∗) and an n-neuron p-precision (in base |2Γ|) RNN with two growing memory modules TW,b : (Qn,Q∗,Q∗) → (Qn,Q∗,Q∗), where n = 2|Γ| + ⌈log2 |Q|⌉ + |Q||Γ| + 19 and p ≥ 2, such that for all instantaneous descriptions x ∈ X , the following holds: If P∗M(x) is defined, then
ρ−1(T ∗W,b(ρ(x))) = P∗M(x), (18)
and if P∗M(x) is not defined, then T ∗W,b(ρ(x)) is also not defined. If P∗M(x) is defined and computed in T steps by M, then T ∗W,b(ρ(x)) is computed in 3T steps by the RNN.
Finally, applying Corollary 2.1 to U6,4, we obtain: Corollary 2.2. There exists a 54-neuron p-precision (in base 8) RNN with two growing memory modules that can simulate any Turing Machine in O(T 6), where T is the number of steps required for the Turing Machine to compute the result and p ≥ 2.
The architecture of the Turing-complete 54-neuron RNN (fully described in the proof of Theorem 2) is depicted in Figure 2.

heading: 4.1 Relationship of the Growing Memory Modules with Stack-augmented RNNs

text:
The proposed growing memory module belongs to the generic class of stack-augmented RNNs, which refers to any RNNs augmented with a stack-like mechanism. Many different forms of stackaugmented RNNs have been proposed [11, 12, 13, 14, 15]. Given the simplicity of the design,
the growing memory module represents the foundation form of these stack-augmented RNNs. It is easy to show that these stack-augmented RNNs can simulate the growing memory module in linear time. For example, the growing memory modules can be simulated by the neural stack [11] as follows. For pushing operation, set vt in the neural stack to u(t) in the growing memory module, and dt in the neural stack to 1{u(t) > 0} in the growing memory module. For popping operation, set ut in the neural stack to 1{o(t) = 1} in the growing memory module. Therefore, the proof for the Turing completeness of bounded-precision RNNs with growing memory modules can extend to other stack-augmented RNNs. That is, bounded-precision stack-augmented RNNs are also Turing-complete.
Different from other stack-augmented RNNs, the proposed growing memory modules use a simple mechanism to control pushing and popping, as only the top neurons in the stack are included in the RNNs. This allows theories relating to growing memory modules to be easily extended to other forms of RNNs.

heading: 5 Bounded-Precision RNNs and Space-Bounded Turing Machines

text:
As discussed above, it is inefficient to update all neurons that encode tape information, but if one still wants to remove the growing memory modules and focuses purely on a bounded-precision RNN only, then the resulting network can simulate space-bounded Turing Machines (that is, Turing Machines with a bounded-size tape) only.
Theorem 3. Given a Turing Machine M with a bounded tape of size F , there exists an injective function ρ : X → Qn and an n-neuron p-precision (in base |2Γ|) RNN TW,b : Qn → Qn, where n = O(⌈F/p⌉) and p ≥ 2, such that for all instantaneous descriptions x ∈ X ,
ρ−1(T 3W,b(ρ(x))) = PM(x). (19)
The proof, which is constructive, can be found in Appendix C. The general idea of the proof is to implement the growing memory module in Section 4 by an RNN as well and place all neurons inside the RNN. The theorem shows that the number of neurons required to simulate a space-bounded Turing Machine correlates with the tape size.
To simulate a Turing Machine with an unbounded tape, we would need to add neurons to the RNN once the read/write head reaches the end of the memory. To be specific, we say an RNN has an
unbounded number of neurons if the RNN either has an infinite number of neurons to start with or can increase the number of neurons after each update step depending on the neurons’ values. The UTM that simulates any Turing Machine in a fast way was described in [21], which does so with only O(T log T ) slowdown. While this UTM has multiple tapes, Theorem 3 can be generalized to multiple-tape Turing Machines easily. We now obtain:
Corollary 3.1. There exists an unbounded-neuron bounded-precision RNN that can simulate any Turing Machine in O(T log T ), where T is the number of steps required for the Turing Machine to compute the result.

heading: 6 Discussion and Conclusion

text:
To construct a Turing-complete RNN, we have to incorporate some encoding for the unbounded number of symbols on the Turing tape. This encoding can be done by: (a) unbounded precision of some neurons (Theorem 1), (b) an unbounded number of neurons (Theorem 3), or (c) a separate growing memory module (Theorem 2). The main contribution of this paper is spelling out the details of (c), which provides a practical way to construct an RNN that runs any given algorithms. We prove the Turing completeness of a 40-neuron unbounded-precision RNN, which is the smallest Turing-complete RNN to date. We analyze the relationship between the number of neurons and the precision of an RNN when simulating a Turing Machine. Most importantly, we propose a 54-neuron bounded-precision RNN with growing memory modules that is Turing-complete, and this proof of Turing completeness can be extended to stack-augmented RNNs in general.
This paper focuses on the computational capabilities and representability of symbolic and subsymbolic processing via stack-augmented RNNs; it does not engage yet with methods to train them. Since growing memory modules are not differentiable, we cannot train them directly by the frequently used error backpropagation algorithm. One may want to construct a differentiable version of the modules, or alternatively use a different learning rule (e.g., REINFORCE [22]) to deal with the discrete pushing and popping operations. Understanding methods to train growing memory modules, incorporate symbolic information into the sub-symbolic representation efficiently, and retrieve both symbolic and non-symbolic information are the next steps towards the goal of combining symbolic and sub-symbolic capabilities in an adaptive and applicable manner.

title: Turing Completeness of Bounded-Precision Recurrent Neural Networks


summaries: In this work, the authors consider linear quantile regression, paying particular attention to the "coverage" properties of the resulting predictor. The traditional goal of quantile regression is to learn a quantile of the conditional distribution of output $Y$, conditioned on the inputs $X$. Denoting the predictor returned by any quantile regression algorithm by $\hat{f}(\cdot)$, one can use this to construct a predictive interval for a freshly-drawn $Y$. The authors consider the one-sided interval $C(X) = (-\infty, \hat{f}(X)]$, and the event $\\{ Y \in C(X) \\}$. The authors call the probability of this event (over the draw of a pair $(X,Y)$) the coverage of $\hat{f}$. Their main results elucidate conditions under which the traditional pinball-loss based ERM quantile regression solution is biased below the desired quantile level, that is, the predictive interval tends to be "too small," due to $\hat{f}$ under-estimating the desired quantile.  In particular, they show that in the limit where the sample size and dimensionality are proportional, the coverage (the aforementioned probability, which depends on $\hat{f}$ and thus the training data) converges below the desired quantile, with this error depending directly on the ratio of dimension to sample size $d/n$. Furthermore, they show (under fixed $n$ and $d$) how for a certain family of data distributions and a linear quantile model, given any estimator of the weights, the downward bias term scales with the weight estimation error. Their empirical tests also effectively highlight the phenomena described formally.

summaries: This paper presents a theoretical study on the coverage of quantile regression algorithms. More specifically, the authors prove theoretically that linear quantile regression exhibits an inherent under-coverage bias, with an order d/n difference in the marginal coverage regardless of the noise distribution (but Gaussian for the covariates). They also show that the primary source of this under-coverage bias is the estimation error in the linear coefficients, while the estimation error in the bias term can have either an under- or over-coverage effect. Numerical experiments further support their findings.

summaries: Quantifying predictive uncertainty is a crucial task for high-stakes prediction problems. The paper focus on quantile regression, a common approach for quantifying predictive uncertainty in regression problems. Authors consider realizable linear setting when the true model belongs to the considered functional class. With a focus on the under-parameterized setting, they theoretically justify the under-coverage of the learnt quantile regressors. Established theoretical results are backed up with supporting empirical evidence.

summaries: This paper gives theoretical insights into the under-coverage phenomenon often observed in quantile regression, i.e., for a continuous distribution, the probability of observing  a true label below a learned  $\alpha$-quantile is often below $\alpha$, which is consider a failure in estimating the uncertainty. Two important results are shown. The first one provides a characterization of the under-coverage in the setting of a Gaussian linear model in the under-parameterized high-dimensional proportional limit (i.e. $d,n \rightarrow \infty$ , $d/n\rightarrow\kappa$ with small $\kappa$) independently of the noise distribution. The experimental results show that the formulas match pretty well the observed coverage in realistic situations of dimension around 100, on synthetic and real data. The second important result shows that the source of the under-coverage stems essentially from the error in estimating the high-dimensional linear coefficient, under some assumptions on the underlying source.

abstractText: Estimating the data uncertainty in regression tasks is often done by learning a quantile function or a prediction interval of the true label conditioned on the input. It is frequently observed that quantile regression—a vanilla algorithm for learning quantiles with asymptotic guarantees—tends to under-cover than the desired coverage level in reality. While various fixes have been proposed, a more fundamental understanding of why this under-coverage bias happens in the first place remains elusive. In this paper, we present a rigorous theoretical study on the coverage of uncertainty estimation algorithms in learning quantiles. We prove that quantile regression suffers from an inherent under-coverage bias, in a vanilla setting where we learn a realizable linear quantile function and there is more data than parameters. More quantitatively, for α > 0.5 and small d/n, the α-quantile learned by quantile regression roughly achieves coverage α− (α− 1/2) · d/n regardless of the noise distribution, where d is the input dimension and n is the number of training data. Our theory reveals that this under-coverage bias stems from a certain highdimensional parameter estimation error that is not implied by existing theories on quantile regression. Experiments on simulated and real data verify our theory and further illustrate the effect of various factors such as sample size and model capacity on the under-coverage bias in more practical setups.

authors:
C. Achilles
H.P. Bain
F. Bellott
J. Boyd-Zaharias
J. Finn
J. Folger
J. Johnston
E. Word

title: Tennessee’s student teacher achievement ratio (star) project

venue: Harvard Dataverse,

authors:
A. Angelopoulos
S. Bates
J. Malik
M.I. Jordan

title: Uncertainty sets for image classifiers using conformal prediction

venue: arXiv preprint arXiv:2009.14193,

authors:
Y. Bai
S. Mei
H. Wang
C. Xiong

title: Don’t just blame over-parametrization for overconfidence: Theoretical analysis of calibration in binary classification

venue: arXiv preprint arXiv:2102.07856,

authors:
R.F. Barber
E.J. Candes
A. Ramdas
R.J. Tibshirani

title: Conformal prediction under covariate shift

venue: arXiv preprint arXiv:1904.06019,

title: The limits of distribution-free conditional predictive inference

title: Predictive inference with the jackknife+

venue: The Annals of Statistics,

authors:
S. Bates
A. Angelopoulos
L. Lei
J. Malik
M.I. Jordan

title: Distribution-free, risk-controlling prediction sets

venue: arXiv preprint arXiv:2101.02703,

authors:
M. Bayati
A. Montanari

title: The dynamics of message passing on dense graphs, with applications to compressed sensing

venue: IEEE Transactions on Information Theory,

authors:
M. Bayati
M. Lelarge
A. Montanari

title: Universality in polytope phase transitions and message passing algorithms

venue: Annals of Applied Probability,

authors:
E. Begoli
T. Bhattacharya
D. Kusnezov

title: The need for uncertainty quantification in machine-assisted medical decision making

venue: Nature Machine Intelligence,

authors:
E.J. Candès
P. Sur

title: The phase transition for the existence of the maximum likelihood estimate in high-dimensional logistic regression

authors:
M. Cauchois
S. Gupta
A. Ali
J.C. Duchi

title: Robust validation: Confident predictions even when distributions shift

venue: arXiv preprint arXiv:2008.04267,

authors:
M. Cauchois
S. Gupta
J.C. Duchi

title: Knowing what you know: valid and validated confidence sets in multiclass and multilabel prediction

venue: Journal of Machine Learning Research,

authors:
A. Christmann
I. Steinwart

title: How svms can estimate quantiles and the median

venue: In Advances in neural information processing systems,

authors:
D. Donoho
A. Montanari

title: High dimensional robust m-estimation: Asymptotic variance via approximate message passing

venue: Probability Theory and Related Fields,

authors:
D.L. Donoho
A. Maleki
A. Montanari

title: Message-passing algorithms for compressed sensing

venue: Proceedings of the National Academy of Sciences,

authors:
D. Dua
C. Graff

title: UCI machine learning repository, 2017

venue: URL http://archive.ics. uci.edu/ml

authors:
N. El Karoui
D. Bean
P.J. Bickel
C. Lim
B. Yu

title: On robust regression with highdimensional predictors

authors:
Y. Gal
Z. Ghahramani

title: Dropout as a bayesian approximation: Representing model uncertainty in deep learning

venue: In international conference on machine learning,

title: The predictive sample reuse method with applications

venue: Journal of the American statistical Association,

authors:
C. Guo
G. Pleiss
Y. Sun
K.Q. Weinberger

title: On calibration of modern neural networks

venue: arXiv preprint arXiv:1706.04599,

authors:
C. Gupta
A. Podkopaev
A. Ramdas

title: Distribution-free binary classification: prediction sets, confidence intervals and calibration

venue: arXiv preprint arXiv:2006.10564,

authors:
G. Huang
Y. Li
G. Pleiss
Z. Liu
J.E. Hopcroft
K.Q. Weinberger

title: Snapshot ensembles: Train 1, get m for free

venue: arXiv preprint arXiv:1704.00109,

authors:
X. Jiang
M. Osl
J. Kim
L. Ohno-Machado

title: Calibrating predictive model estimates to support personalized medicine

venue: Journal of the American Medical Informatics Association,

authors:
C. Jung
C. Lee
M.M. Pai
A. Roth
R. Vohra

title: Moment multicalibration for uncertainty estimation

venue: arXiv preprint arXiv:2008.08037,

title: Asymptotic behavior of unregularized and ridge-regularized high-dimensional robust regression estimators: rigorous results

venue: arXiv preprint arXiv:1311.2445,

authors:
A. Kendall
Y. Gal

title: What uncertainties do we need in bayesian deep learning for computer vision

venue: arXiv preprint arXiv:1703.04977,

authors:
D. Kivaranovic
K.D. Johnson
H. Leeb

title: Adaptive, distribution-free prediction intervals for deep networks

venue: In International Conference on Artificial Intelligence and Statistics,

authors:
R. Koenker
G. Bassett Jr.

title: Regression quantiles

venue: Econometrica: journal of the Econometric Society,

authors:
R. Koenker
K.F. Hallock

title: Quantile regression

venue: Journal of economic perspectives,

authors:
A. Kumar
P. Liang
T. Ma

title: Verified uncertainty calibration

venue: arXiv preprint arXiv:1909.10155,

authors:
B. Lakshminarayanan
A. Pritzel
C. Blundell

title: Simple and scalable predictive uncertainty estimation using deep ensembles

venue: arXiv preprint arXiv:1612.01474,

title: Classification with confidence

authors:
J. Lei
M. G’Sell
A. Rinaldo
R.J. Tibshirani
L. Wasserman

title: Distribution-free predictive inference for regression

venue: Journal of the American Statistical Association,

authors:
L.T. Liu
M. Simchowitz
M. Hardt

title: The implicit fairness criterion of unconstrained learning

venue: In International Conference on Machine Learning,

title: Bayesian methods for adaptive models

venue: PhD thesis, California Institute of Technology,

authors:
W.J. Maddox
P. Izmailov
T. Garipov
D.P. Vetrov
A.G. Wilson

title: A simple baseline for bayesian uncertainty in deep learning

venue: Advances in Neural Information Processing Systems,

authors:
X. Mai
Z. Liao
R. Couillet

title: A large scale analysis of logistic regression: Asymptotic performance and new insights

venue: IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP),

authors:
A. Malinin
M. Gales

title: Predictive uncertainty estimation via prior networks

venue: arXiv preprint arXiv:1802.10501,

authors:
A. Malinin
B. Mlodozeniec
M. Gales

title: Ensemble distribution distillation

venue: arXiv preprint arXiv:1905.00076,

title: Quantile regression forests

authors:
R. Michelmore
M. Kwiatkowska
Y. Gal

title: Evaluating uncertainty quantification in end-toend autonomous driving control

venue: arXiv preprint arXiv:1811.06817,

title: Last iterate of sgd converges (even in unbounded domains), 2020; Accessed: May, 2021

authors:
Y. Ovadia
E. Fertig
J. Ren
Z. Nado
D. Sculley
S. Nowozin
J.V. Dillon
B. Lakshminarayanan
J. Snoek

title: Can you trust your model’s uncertainty? evaluating predictive uncertainty under dataset shift

venue: arXiv preprint arXiv:1906.02530,

title: Inductive conformal prediction: Theory and application to neural networks

venue: In Tools in artificial intelligence. Citeseer,

title: Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods

title: Approximate tests of correlation in time-series

venue: Journal of the Royal Statistical Society: Series B (Methodological),

authors:
Y. Romano
E. Patterson
E.J. Candès

title: Conformalized quantile regression

venue: arXiv preprint arXiv:1905.03222,

authors:
E. Shabat
L. Cohen
Y. Mansour

title: Sample complexity of uniform convergence for multicalibration

venue: arXiv preprint arXiv:2005.01757,

authors:
G. Shafer
V. Vovk

title: A tutorial on conformal prediction

authors:
I. Steinwart
A. Christmann

title: Estimating conditional quantiles with the help of the pinball

venue: loss. Bernoulli,

title: A framework to characterize performance of lasso algorithms

venue: arXiv preprint arXiv:1303.7291,

title: Cross-validatory choice and assessment of statistical predictions

authors:
P. Sur
E.J. Candès

title: A modern maximum-likelihood theory for high-dimensional logistic regression

authors:
I. Takeuchi
Q.V. Le
T.D. Sears
A.J. Smola

title: Nonparametric quantile estimation

authors:
C. Thrampoulidis
S. Oymak
B. Hassibi

title: Regularized linear regression: A precise analysis of the estimation error

venue: In Conference on Learning Theory,

authors:
C. Thrampoulidis
E. Abbasi
B. Hassibi

title: Precise error analysis of regularized m-estimators in high dimensions

title: Bias and confidence in not quite large samples

venue: Ann. Math. Statist.,

title: Asymptotic statistics, volume 3

venue: Cambridge university press,

title: High-dimensional probability: An introduction with applications in data science, volume 47

title: Conditional validity of inductive conformal predictors

venue: In Asian conference on machine learning,

title: Cross-conformal predictors

venue: Annals of Mathematics and Artificial Intelligence,

authors:
V. Vovk
A. Gammerman
G. Shafer

title: Algorithmic learning in a random world

authors:
V. Vovk
I. Nouretdinov
V. Manokhin
A. Gammerman

title: Cross-conformal predictive distributions

venue: In Conformal and Probabilistic Prediction and Applications,

title: Determination of sample sizes for setting tolerance limits

venue: The Annals of Mathematical Statistics,

title: Statistical prediction with special reference to the problem of tolerance limits

venue: The annals of mathematical statistics,

authors:
B. Zadrozny
C. Elkan

title: Transforming classifier scores into accurate multiclass probability estimates

venue: In Proceedings of the eighth ACM SIGKDD international conference on Knowledge discovery and data mining,

heading: 1 Introduction

text:
This paper is concerned with the problem of uncertainty estimation in regression problems. Uncertainty estimation is an increasingly important task in modern machine learning applications—Models should not only make high-accuracy predictions, but also have a sense of how much the true label may deviate from the prediction. This capability is crucial for deploying machine learning in the real world, in particular in risk-sensitive domains such as medical AI [15, 29], self-driving cars [47], and so on. A common approach for uncertainty estimation in regression is to learn a quantile function or a prediction interval of the true label conditioned on the input, which provides useful distributional information about the label. Such learned quantiles are typically evaluated by their coverage, i.e., probability that it covers the true label on a new test example. For example, a learned 90% upper quantile function should be an actual upper bound of the true label at least 90% of the time.
Algorithms for learning quantiles date back to the classical quantile regression [35], which estimates the quantile function by solving an empirical risk minimization problem with a suitable loss function that depends on the desired quantile level α. Quantile regression is conceptually simple, and is
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
theoretically shown to achieve asymptotically correct coverage as the sample size goes to infinity [34] or approximately correct coverage in finite samples under specific modeling assumptions [46, 60, 56]. However, it is observed that quantile regression often under-covers than the desired coverage level in practice [53]. Various alternative approaches for constructing quantiles and confidence intervals are proposed in more recent work, for example by aggregating multiple predictions using Bayesian neural networks or ensembles [24, 37], or by building on the conformal prediction technique to construct prediction intervals with finite-sample coverage guarantees [68, 66, 39, 53]. However, despite these advances, a more fundamental understanding on why vanilla quantile regression exhibits this under-coverage bias is still lacking.
This paper revisits quantile regression and presents a first precise theoretical study on its coverage, in a new regime where the number of samples n is proportional to the dimension d, and the ratio d/n is small (so that the problem is under-parametrized). Our main result shows that quantile regression exhibits an inherent under-cover bias under this regime, even in the well-specified setting of learning a linear quantile function when the true data distribution follows a Gaussian linear model. To the best of our knowledge, this is the first rigorous theoretical justification of the under-coverage bias. Our main contributions are summarized as follows.
• We prove that linear quantile regression exhibits an inherent under-coverage bias in the wellspecified setting where the data is generated from a Gaussian linear model, and the number of samples n is proportional to the feature dimension d with a small d/n (Section 3). More quantitatively, quantile regression at nominal level α ∈ (0.5, 1) roughly achieves coverage α− (α− 1/2)d/n regardless of the noise distribution. To the best of our knowledge, this is the first rigorous characterization of the under-coverage bias in quantile regression.
• Towards understanding the source of this under-coverage bias, we disentangle the effect of estimating the bias and estimating the linear coefficient on the coverage of the learned linear quantile (Section 4). We show that the estimation error in the bias can have either an undercoverage or over-coverage effect, depending on the noise distribution. In contrast, the estimation error in the linear coefficient always drives the learned quantile to under-cover, and we show this effect is present even on broader classes of data distributions beyond the Gaussian linear model.
• We perform experiments on simulated and real data to test our theory (Section 5). Our simulations show that the coverage of quantile regression in Gaussian linear models agrees well with our precise theoretical formula as well as the α− (α− 1/2)d/n approximation. On real data, we find quantile regression using high-capacity models (such as neural networks) exhibits severe under-coverage biases, while linear quantile regression can also have a mild but non-negligible amount of under-coverage, even after we remove the potential effect of model misspecification.
• On the technical end, our analysis builds on recent understandings of empirical risk minimization problems in the high-dimensional proportional limit with a small d/n, and develops new techniques such as a novel concentration argument to deal with an additional learnable variable in learning linear models with biases, which we believe could be of further interest (Section 6).

heading: 1.1 Related work

text:
Algorithms for uncertainty estimation in regression The earliest methods for uncertainty estimation in regression adopted subsampling methods (bootstrap) or leave-one-out methods (Jackknife) for assessing or calibrating prediction uncertainty [52, 63, 58, 25]. More recently, a growing line of work builds on the idea of conformal prediction [55] to design uncertainty estimation algorithms for regression. These algorithms provide confidence bounds or prediction intervals by post-processing any predictor, and can achieve distribution-free finite-sample marginal coverage guarantees utilizing exchangeability of the data [50, 66, 39, 53, 33, 67, 69, 68, 11]. Further modifications of the conformal prediction technique can yield stronger guarantees such as group coverage [10] or coverage under distribution shift [9] under additional assumptions. Our under-coverage results advocate the necessity of such post-processing techniques, and are complementary in the sense that we provide understandings on the more vanilla quantile regression algorithm. Quantiles and prediction intervals can also be obtained by aggregating multiple predictors, such as using Bayesian neural networks [41, 24, 32, 44, 42] or ensembles [37, 49, 28, 45]. These methods offer an alternative approach for uncertainty estimation, but do not typically come with coverage guarantees.
Theoretical analysis of quantile regression Linear quantile regression with the pinball loss dates back to the late 1970s [34]. The same work proved the asymptotic normality of the regression coefficients in the n → ∞, fixed d limit. Takeuchi et al. [60] studied non-parametric quantile regression using kernel methods, and provided generalization bounds (with the pinball loss) based on the Rademacher complexity. Meinshausen [46] studied non-parametric quantile regression using random forest and showed its consistency under proper assumptions. Christmann and Steinwart [19], Steinwart et al. [56] established a “self-calibration” inequality for the quantile loss, which, when combined with standard generalization bounds, can be translated to an estimation error bound for quantile regression. These works all focus on bounding the parameter or function estimation error, which can be translated to bounds on the coverage bias, but does not tell the sign of this coverage bias as we do in this paper. We also remark that conformalization can be used in conjunction with quantile regression to correct its coverage bias [53].
Uncertainty quantification for classification For classification problems, two main types of uncertainty quantification methods have been considered: outputting discrete prediction sets with guarantees of covering the true (discrete) label [70, 71, 38, 7, 12, 18, 17], or calibrating the predicted probabilities [51, 72, 73, 37, 26]. The connection between prediction sets and calibration was discussed in [27]. The sample complexity of calibration has been studied in a number of theoretical works [36, 27, 54, 30, 40, 8]. Our work is inspired by the recent work of Bai et al. [8], which showed that logistic regression is over-confident even if the model is correctly specified and the sample size is larger than the dimension.
High-dimensional behaviors of empirical risk minimization There is a rapidly growing literature on limiting characterizations of convex optimization-based estimators in the n ∝ d regime [21, 13, 23, 31, 57, 61, 20, 62, 43, 59, 16]. Our analysis builds on results for unregularized Mestimator derived in [62] and generalizes theirs in certain aspects (see also [23, 20, 31]).

heading: 2 Preliminaries

text:
In this paper we focus on the problem of learning quantiles. Suppose we observe a training dataset {(xi, yi)}ni=1 drawn i.i.d. from some joint distribution P on Rd × R, where xi ∈ Rd is the input features and y ∈ R is the real-valued response (label). Let F (t|x) := P(Y ≤ t|X = x) denote the conditional CDF of Y |X. Our goal is to learn the α-(conditional) quantile of Y |X:
q?α(x) := inf {t ∈ R : F (t|x) ≥ α}. For example, q?0.95(x) is the ground truth 95% quantile of the true conditional distribution Y |X, and can be seen as the “ideal” 95% upper confidence bound for the label y given the features x. Throughout this paper we work with upper quantiles, that is, α ∈ (0.5, 1) (some typical choices are α ∈ {0.8, 0.9, 0.95}); by symmetry our results hold for learning lower quantiles as well.
Coverage For any learned quantile function f̂ : Rd → R, the marginal coverage (henceforth “coverage”) of f̂ is the probability of y ≤ f̂(x) on a new test example (x, y):
Coverage(f̂) := P(x,y) ( y ≤ f̂(x) ) = Ex [ P ( y ≤ f̂(x)|x )] . (1)
For learning the α-quantile (α > 0.5), we usually expect Coverage(f̂) ≈ α, i.e. f̂(x) covers the label y on approximately α proportion of the data, under the ground truth data distribution.
We say that f̂ has under-coverage if Coverage(f̂) < α and over-coverage if Coverage(f̂) > α. Note that these two notions are not symmetric: Over-coverage means that the learned upper quantile f̂(x) is overly conservative (higher than enough), and is typically tolerable; In contrast, under-coverage means that f̂(x) fails to cover y with α probability, and is typically considered as a failure. We remark that while there exist more fine-grained notions of coverage such as conditional coverage [10], the (marginal) coverage is still a basic requirement for any quantile learning algorithm.
Quantile regression We consider quantile regression, a standard method for learning quantiles from data [35]. Quantile regression estimates the true quantile function qα(·) via the pinball loss [34, 56]
`α(t) = −(1− α)t1 {t ≤ 0}+ αt1 {t > 0} . (2)
Note that in the special case of α = 0.5, we have `0.5(t) = |t|/2, and thus the pinball loss strictly generalizes the absolute loss (for learning medians) to learning any quantile. Given the training dataset and any function class {fθ : θ ∈ Θ} (e.g. linear models or neural networks), quantile regression solves the (unregularized) empirical risk minimization (ERM) problem
θ̂ = arg min θ∈Θ
R̂n(θ) := 1
n n∑ i=1 `α(yi − fθ(xi)). (3)
(We take θ̂ as any minimizer of R̂n when the minimizer is non-unique.) Let R(θ) := E[R̂n(θ)] denote the corresponding population risk. It is known that the population risk over all (measurable) functions is minimized at the true quantile q?α = arg minf R(f) under minimal regularity conditions (for completeness, we provide a proof in Appendix B.1).

heading: 3 Quantile regression exhibits under-coverage

text:
We analyze quantile regression in the vanilla setting where the input distribution is a standard Gaussian and y follows a linear model of x:
y = w>? x + z, where x ∼ N(0, Id), z ∼ Pz. (4)
Above, w? ∈ Rd is the ground truth coefficient vector, and the noise z ∼ Pz is independent of x. The Gaussian input assumption is required only for technical convenience in the high-dimensional limiting analysis, and we believe it is not strictly required for the same result to hold1 (an extension to more general input distributions can also be found in Section 4). The noise distribution Pz is required to satisfy the following smoothness assumption, but can otherwise be arbitrary: Assumption A (Smooth density). The noise distribution Pz has a smooth density φz ∈ C∞(R) (with corresponding CDF Φz), with bounded derivatives: supt∈R |φ (k) z (t)| < ∞ for any k ≥ 0. We further assume that φz(zα) > 0, where α ∈ (0.5, 1) is our pre-specified quantile level, and zα := inf {t ∈ R : Φz(t) ≥ α} is the α-quantile of Pz .
Under the above model, it is straightforward to see that the true α-conditional quantile of y|x is also a linear model (with bias):
q?α(x) = w > ? x + zα. (5)
Given the training data {(xi, yi)}ni=1, we learn a linear quantile function f̂(x) = ŵ>x + b̂ via quantile regression:
(ŵ, b̂) = arg min w,b
R̂n(w, b) := 1
n n∑ i=1 `α(yi − (w>xi + b)), (6)
where `α is the pinball loss in (2). As our linear function class realizes the true quantile function (5), the population risk is minimized at the true quantile: arg minw,bR(w, b) = (w?, zα).
We are now ready to state our main result, which shows that quantile regression exhibits an inherent under-coverage bias even in this vanilla realizable setting. Theorem 1 (Quantile regression exhibits under-coverage bias). Suppose the data is generated from the linear model (4) and the noise satisfies Assumption A. Let f̂(x) = ŵ>x + b̂ be the output of quantile regression (6) at level α ∈ (0.5, 1). Then, in the limit of n, d → ∞ and d/n → κ where κ ∈ (0, κ0] for some small κ0 > 0, for the coverage (1), we have (
p→ denotes convergence in probability)
Coverage(f̂) p→ α− Cα,κ for some Cα,κ > 0.
That is, the limiting coverage of the learned quantile function is less than α. Further, for small enough κ we have the local linear expansion
Cα,κ = (α− 1/2)κ+ o(κ). (7) 1Our results can be extended directly to any correlated Gaussian input x ∼ N(0,Σ) by the transform x̃ = Σ−1/2x and w̃? = Σ1/2w?. We believe our results also hold for i.i.d. sub-Gaussian inputs by the universality principle (e.g. [14]).
Theorem 1 builds on the precise characterization of ERM problems in the high-dimensional proportional limit [62], along with new techniques over existing work for dealing with the unique challenges in quantile regression (such as analyzing the additional learnable bias b). An overview of the main technical steps is provided in Section 6, and the full proof is deferred to Appendix C.
Implications Theorem 1 can be illustrated by the following numeric example. Suppose we perform quantile regression at α = 0.9, where the data follows the linear model (4), and our κ = d/n = 0.1 (so that the sample size is 10x number of parameters). Then Theorem 1 shows that, even in this realizable, under-parametrized setting, the coverage of the learned quantile f̂ is going to be roughly 0.9− Cα,κ when n, d are large, and further Cα,κ ≈ (α− 1/2)κ = 0.04. Thus the actual coverage is around 0.9− 0.04 = 0.86, and such a 4% under-coverage bias can be rather non-negligible in reality. We remark that a symmetric conclusion of Theorem 1 also holds for lower quantiles, and we further expect similar results also hold for learning prediction intervals, where the coverage is defined as the two-sided coverage of the prediction interval formed by the learned {lower quantile, upper quantile}. To the best of our knowledge, this offers a first precise theoretical understanding of why practically trained quantiles or prediction intervals often under-cover than the desired coverage level [53].
Comparison against existing theories An important feature of the under-coverage bias shown in Theorem 1 is that it only shows up in the n, d proportional regime, and is not implied by existing theories on quantile regression. Classical asymptotic theory only shows asymptotic normality √ n([ŵ, b̂] − [w?, zα]) → N(0,V) in the n → ∞, fixed d limit [34, 64]. Under this limit, Coverage(f̂) is consistent at α with O(1/ √ n) deviation. Christmann and Steinwart [19], Steinwart et al. [56] consider the finite n, d setting and establish self-calibration inequalities (similar to strong
convexity) that bounds the quantile estimation error by the square root excess loss √ R(f̂)−R(q?α). Combined with standard generalization theories (e.g. via Rademacher complexities) and Lipschitzness, this can be turned into a bound on |Coverage(f̂) − α|, but does not tell the sign (positive or negative) of the coverage bias.
Large κ; extension to over-parametrized learning While Theorem 1 requires a small κ = d/n, the approximation formula (7) suggests that the under-coverage should get more severe as κ—the measure of over-parametrization in this problem—gets larger. We confirm this trend experimentally in our simulations in Section 5.1.
As an extension to Theorem 1, we also show theoretically that the under-coverage bias indeed becomes even more severe in over-parametrized learning, under the same linear model (4): When d/n > Õ(1), and the noise Pz is sub-Gaussian and symmetrically distributed about 0, the convergence point of the gradient descent path on the quantile regression risk R̂n is the minimum-norm interpolator of the data, which has coverage 0.5± Õ(1/ √ d) with high probability (see Appendix D for the formal statement and the proof). Notably, this 0.5 coverage does not depend on α and exhibits a severe under-coverage.

heading: 4 Understanding the source of the under-coverage bias

text:
In this section, we take steps towards a deeper understanding of how the under-coverage bias shown in Theorem 1 happens. Recall that the quantile regression returns f̂(x) = ŵ>x+ b̂ where (ŵ, b̂) is a solution to the ERM problem (6) and estimates the true parameters (w?, zα). Our main approach in this section is to disentangle the effect of the two sources—the estimation error in b̂ and the estimation error in w—on the coverage of f̂ .
We show that the estimation error in b̂ can have either an under-coverage or an over-coverage effect, depending on the noise distribution (Section 4.1). In contrast, the estimation error in ŵ always has an under-coverage effect; this holds not only for the linear model assumed in Theorem 1, but also on more general data distributions (Section 4.2). In the setting of Theorem 1, this under-coverage effect of ŵ is always strong enough to dominate the effect of b̂, leading to the overall under-coverage.

heading: 4.1 Effect of estimation error in b̂

text:
To study the effect of b̂, we use the quantity b̂ − zα as a measure for its effect on the coverage— Recall that the true quantile is q?α(x) = w > ? x + zα, thus having b̂ < zα means that b̂ contributes to under-coverage, whereas b̂ > zα means b̂ contributes to over-coverage. (This can be seen more straightforwardly in the easier case where we know w? and only output b̂ to estimate zα.)
The following corollary shows that, under the same settings of Theorem 1, the error b̂− zα can be understood precisely. The proof can be found in Appendix E.1.
Corollary 2 (Effect of b̂ on coverage depends on noise distribution). Under the same settings as Theorem 1, for any α ∈ (0.5, 1), as n, d→∞ with d/n→ κ ∈ (0, κ0], we have
(a) The learned bias b̂ from quantile regression (6) converges to the following limit:
b̂− zα p→ Cbα,κ = b̄0κ+ o(κ),
where b̄0 has a closed-form expression:
b̄0 := −α(1− α)φ′z(zα)− (2α− 1)φ2z(zα)
2φ3z(zα) . (8)
(b) For any α ∈ (0.5, 1), when Pz is the Gaussian distribution (with arbitrary scale), we have b̄0 < 0 in which case Cbα,κ < 0 for small enough κ. Conversely, for any α ∈ (0.5, 1), there exists some noise distribution Pz for which b̄0 > 0, in which case Cbα,κ > 0 for small enough κ.
Corollary 2 shows that the sign ofCbα,κ in the limiting regime (and thus the effect of b̂ on the coverage) depends on b̄0, which in turn depends on the noise distribution Pz . For common noise distributions such as Gaussian we have Cbα,κ < 0 at small κ, but there also exists Pz such that C b α,κ > 0. Note that the second claim in part (b) follows directly from (8): we can always design the density φz by varying φz(zα) and φ′z(zα) so that b̄0 > 0. Overall, this result shows that the under-coverage bias in Theorem 1 cannot be simply explained by the under-estimation error in b̂.

heading: 4.2 Effect of estimation error in ŵ; relaxed data distributions

text:
We now show that the primary source of the under-coverage is the estimation error in ŵ, which happens not only on the linear data distribution assumed in Theorem 1, but also on a broader class of data distributions. We consider the following relaxed data distribution assumption
y = µ?(x) + σ?(x)z, (9)
where the noise z ∼ Pz . We do not put structural assumptions on (µ?, σ?), except that we assume the true α-quantile is still a linear function of x, that is, there exists (w?, b?) for which
q?α(x) = µ?(x) + σ?(x)zα = w > ? x + b?. (10)
Since here we are interested in the effect of estimating w?, for simplicity, we assume that we know b? and only estimate w? via some estimator ŵ. We now collect our assumptions and state the result.
Assumption B (Relaxed data distribution). The data is distributed as model (9) with a linear αquantile function (10). Further, the data distribution satisfies the following regularity conditions:
(a) The distribution of x ∈ Rd is symmetric about 0, has a lower bounded covariance E[xx>] γId, and is K-sub-Gaussian, for constants γ,K > 0.
(b) The variance function σ?(·) is bounded and symmetric: For all x ∈ Rd we have σ ≤ σ?(x) ≤ σ for some constants σ, σ > 0, and σ?(x) = σ?(−x).
(c) The noise density φz is continuously differentiable and symmetric about 0, i.e. φz(t) = φz(−t) for all t ∈ R. Further, φz is uni-modal, i.e. φ′z(t)|t<0 > 0 and φ′z(t)|t>0 < 0.
Theorem 3 (Estimation error in ŵ leads to under-coverage on a family of data distributions). Under the relaxed data distribution assumption (Assumption B), for any α > 3/4, there exists constants c, r0 > 0 such that for any learned quantile estimate f̂(x) = ŵ>x + b? with small estimation error ‖ŵ −w?‖2 ≤ r0, we have
Coverage(f̂) ≤ α− cγ/σ2 · ‖ŵ −w?‖22 ,
that is, the learned quantile under-covers by at least Ω(‖ŵ −w?‖22). Above, c > 0 is an absolute constant, and r0 > 0 depends on (γ, σ,K,Φz, α) but not (n, d).
Implications; proof intuition Theorem 3 shows that, for a broad class of data distributions, any estimator ŵ will under-cover by at least Ω(‖ŵ −w?‖22). If particular, any estimator satisfying ‖ŵ −w?‖2 Õ( √ d/n) (e.g. from standard generalization theory) will under-cover by Õ(d/n). This confirms that the estimation error in the (bulk) regression coefficient ŵ is the primary source of the under-coverage bias, under assumptions that are more general than Theorem 1 in certain aspects (such as the distribution of x and y|x). We remark that as opposed to Theorem 1, Therorem 3 does not give an end-to-end characterization of any specific algorithm, but assumes we have an estimator ŵ with a small error.
At a high-level, Theorem 3 follows from the fact that any estimator ŵ>x+ b? must be lower than the true quantile w>? x + b? for some x and higher for some other x. Averaging the coverage indicator over x, such under-coverage and over-coverage cancel out on the first-order if x has a symmetric distribution, but aggregate to yield a under-coverage effect on the second-order as long as Φz(t)|t>0 is concave (which holds if Pz is unimodal). The proof of Theorem 3 can be found in Appendix E.2.

heading: 5 Experiments

heading: 5.1 Simulations

text:
Setup We first test our Theorem 1 via simulations. We generate data from the linear model (4) in d = 100 dimensions with ‖w?‖2 = 1 and noise distribution Pz = N(0, 0.25). We vary κ = d/n ∈ {0.02, 0.04, . . . , 0.5} where κ determines a sample size n, and vary α ∈ {0.5, 0.52, ..., 0.98}. For each combination of (α, κ), we generate 8 random problem instances, and solve the quantile regression ERM problem (6) on each instance via (sub)-gradient descent. We evaluate the coverage of the learned quantile f̂ (thanks to the linear model (4), the coverage can be computed exactly without needing to introduce a test set). Additional details about the setup can be found in Appendix F.1.
Results Figure 1 plots the coverage of the learned quantiles. Observe that quantile regression exhibits under-coverage consistently across different values of (α, κ). Figure 1a shows that at fixed κ, the amount of under-coverage gets more severe at a higher α, which is qualitatively consistent with our approximation formula (α− 1/2)κ. Figure 1b further compares simulations with our analytical formula α − Cα,κ (found numerically through solving the system of equations 12), as well as the
local linear approximation α− (α− 1/2)κ claimed in Theorem 1. Note that the simulations agree extremely well with the analytical formula. The approximation α− (α− 1/2)κ is also very accurate for almost all κ at α = 0.8, and accurate for small κ at α = 0.9, 0.95. These verify our Theorem 1 and suggests it holds at rather realistic values of the dimension (d = 100).

heading: 5.2 Real data experiments

text:
Datasets and models We take six real-world regression datasets: community and crimes (Community) [2], bike sharing (Bike) [1], Tennessee’s student teacher achievement ratio (STAR) [6], as well as the medical expenditure survey number 19 (MEPS_19) [3], number 20 (MEPS_20) [4], and number 21 (MEPS_21) [5]. All datasets are pre-processed to have standarized features and randomly split into a 80% train set and 20% test set.
To go beyond linear quantile functions, we perform quantile regression with one of the following four models as our fθ: linear model (Linear), a 3-layer MLP (two non-linear layers) with width 64 (MLP-3-64), 512 (MLP-3-512), and a variant of the width-512 MLP where all representation layers are frozen and only the last linear layer is trained (MLP-freeze-3-512). All linear layers include a trainable bias. We minimize the α-quantile loss (3) via momentum SGD with batch size 64. For each setting, we average over 8 random seeds where each seed determines the train-validation split, model initialization, and SGD batching. In our real experiments we fix α = 0.9. (Results at α ∈ {0.8, 0.95} as well as additional experimental setups can be found in Appendix F.2).
Results Table 1 reports the coverage of the learned quantile functions (evaluated on the test sets). Observe that all MLPs exhibit under-coverage compared with the nominal level 90%. Additionally, the amount of under-coverage correlates well with model capacity—the two vanilla MLPs undercovers more severely than the MLP-freeze and the linear model. Notice that the linear model does not have a notable under-coverage on most datasets—we believe this is a consequence of d/n being small on these datasets. The only exception is the Community dataset with the highest d/n ≈ 1/16, on which the linear model does under-cover mildly by roughly 1%.

heading: 5.3 Linear quantile regression on pseudo-labels

text:
To further test the coverage of linear quantile regression on real data distributions, we make two modifications: (1) We subset the training data by fixing d and reducing n, so as to test the coverage across differerent values of κ = d/n; (2) We compare linear quantile regression on both true labels yi, and pseudo-labels y pseudo i generated from estimated linear models. These pseudo-labels are generated by first fitting a linear model ŵ ∈ Rd (with square loss) on the training data, and then generating a new label using the fitted linear model ŵ:
ypseudoi = ŵ >xi + σ̂zi,
where σ̂ is estimated as √
Ê(x,y)[(y − ŵ>x)2] on a separate hold-out split, and zi ∼ N(0, 1). The motivation for the pseudo-labels is to make sure that the data comes from a true linear model, removing the potential effect of model misspecification.
Table 2 shows that on the MEPS_20 dataset, linear quantile regression exhibits under-coverage at relatively large values of κ (0.1, 0.2, 0.5) for both kinds of labels. Also, there is no notable difference
between pseudo-labels and true labels. This provides evidence that our theory on linear quantile regression may hold broadly on real-world data distributions.

heading: 6 Proof overview of Theorem 1

text:
Closed-form expression for coverage Our first step is to obtain a closed-form expression for the coverage. Recall that
Coverage(f̂) := P(x,y)(y ≤ f̂(x)) = P(x,z)(〈w?,x〉+ z ≤ 〈ŵ,x〉+ b̂).
As x is standard Gaussian, and the random variable z has cumulative distribution function Φz , standard calculation then yields the closed form expression (Lemma B.1)
Coverage(f̂) = EG∼N(0,1)[Φz(‖ŵ −w?‖2G+ b̂)].
Concentration of ‖ŵ − w?‖2 and b̂ We generalize results from recent advances in highdimensional M-estimator in linear models [23, 20, 31, 62] to show that ‖ŵ −w?‖2 and b̂ obtained by quantile regression 6 concentrates around fixed values in the high-dimensional limit. We show that, in the limit of d, n→∞ and d/n→ κ, the following concentration happens:
‖ŵ −w?‖2 p→ τ?(κ), and b̂ p→ b?(κ). (11)
Above, τ? and b? are determined by the solutions of a system of nonlinear equations with three variables (τ, λ, b):  τ2κ = λ2 · E(G,Z)∼N(0,1)×Pz [e ′ `αb (τG+ Z;λ)2], τκ = λ · E(G,Z)∼N(0,1)×Pz [e ′ `αb (τG+ Z;λ)G],
0 = E(G,Z)∼N(0,1)×Pz [e ′ `αb (τG+ Z;λ)],
(12)
where e`(x; τ) = minv 12τ (x − v) 2 + `(v) and `αb = ` α(t − b) is the shifted pinball loss (2). (See Theorem C.1 for the formal statement.) This is established via two main steps: We first build on the results of Thrampoulidis et al. [62] to show that a variant of the risk minimization problem with a fixed bias b concentrates around the solution to the first two equations in (12). We then develop a novel concentration argument to deal with the additional learnable bias b in the minimization problem (6), which introduces the third equation in (12) that will be used in characterizing the limiting value of the minimizer b̂.
The concentration (11) implies that Coverage(f̂) also converges to the following limiting coverage value (Lemma C.3):
Coverage(f̂) p→ EG∼N(0,1)[Φz(τ?(κ)G+ b?(κ))] =: α− Cα,κ. (13)
Calculating the limiting coverage via local linear analysis In this final step, as another technical crux of the proof, we further evaluate the small κ approximation of coverage value (13), and determine the sign of Cα,κ. This is achieved by a local linear analysis on the solutions of the aforementioned system of equations at small κ (Lemma C.2) in a similar fashion as in [8], and a precise analysis on the interplay between the concentration values τ?, b?, and the noise density φz . Combining these calculations yields that Cα,κ/κ = −(α−1/2) +o(1) for small enough κ (Lemma C.4). As α > 1/2, this establishes Theorem 1. All details on these analyses can be found in our proofs in Appendix C.

heading: 7 Conclusion

text:
This paper presents a first theoretical justification of the under-coverage bias in quantile regression. We prove that quantile regression suffers from an inherent under-coverage bias even in well-specified linear settings, and provide a precise quantitative characterization of the amount of the under-coverage bias on Gaussian linear models. Our theory further identifies the high-dimensional estimation error in the regression coefficient as the main source of this under-coverage bias, which holds more generally on a broad class of data distributions. We believe our work opens up several interesting directions for future work, such as analyzing non-linear quantile regression, as well as analyzing other notions of uncertainty in regression problems.
Funding transparency statement
YB, HW, CX are funded through employment with Salesforce.

title: Understanding the Under-Coverage Bias in Uncertainty Estimation