Geometry of Program Synthesis

The challenge to efficiently learn arbitrary computable functions, that is, synthesise universal computation, seems a prerequisite to build truly intelligent systems. There is strong evidence that neural networks are necessary but not sufficient for this task (Zador, 2019; Akhlaghpour, 2020). Indeed, those functions efficiently learnt by a neural network are those that are continuous on a compact domain and that process memory in a rather primitive manner. This is a rather small subset of the complete set of computable functions, for example, those computable by a Universal Turing Machine. This suggests an opportunity to revisit the basic conceptual framework of modern machine learning by synthesising Turing machines themselves and re-evaluating our concept of programs. We will show that this leads to new unexplored directions for program synthesis.

In the classical model of computation a program is a code for a Universal Turing Machine (UTM). As these codes are discrete, there is no nontrivial “space” of programs. However this may be a limitation in the classical model, and not a fundamental feature of computation. It is worth bearing in mind that, in reality, there is no such thing as a perfect logical state: in any implementation of a UTM in matter a code $w\in\{0,1\}^{*}$ is realised by an equilibrium state of some physical system (say a magnetic memory). When we say the code on the description tape of the physical UTM “is” $w$ what we actually mean is, adopting the thermodynamic language, that the system is in a phase (a local minima of the free energy) including the microstate $c$ we associate to $w$. However, when the system is in this phase its microstate is not equal to $c$ but rather undergoes rapid spontaneous transitions between many microstates “near” $c$.

So in any possible physical realisation of a UTM, a program is realised by a phase of the physical system. Does this have any computational significance?

One might object that this is physics and not computer science. But this boundary is not so easily drawn: it is unclear, given a particular physical realisation of a UTM, which of the physical characteristics are incidental and which are essentially related to computation, as such. In his original paper Turing highlighted locality in space and time, clearly inspired by the locality of humans operating with eyes, pens and paper, as essential features of his model of computation. Other physical characteristics are left out and hence implicitly regarded as incidental.

In recent years there has been a renewed effort to synthesise programs by gradient descent (Neelakantan et al., 2016; Kaiser & Sutskever, 2016; Bunel et al., 2016; Gaunt et al., 2016; Evans & Grefenstette, 2018; Chen et al., 2018). However, there appear to be serious obstacles to this effort stemming from the singular nature of the loss landscape. When taken seriously at a mathematical level these efforts necessarily come up against this question of programs versus phases. When placed within the correct mathematical setting of singular learning theory (Watanabe, 2009), a clear picture arises of programs as singularities which complements the physical interpretation of phases. The loss function is given as the KL-divergence, an analytic function on the manifold of learnable weights, which corresponds to the microscopic Hamiltonian of the physical system. The Bayesian posterior corresponds to the Boltzmann distribution and we see emerging a dictionary between singular learning theory and statistical mechanics.

The purpose of this article is to (a) explain a formal mathematical setting in which it can be explained why this is the case (b) show how to explore this setting experimentally in code and (c) sketch the new avenues of research in program synthesis that this point of view opens up. The conclusion is that computation should be understood within the framework of singular learning theory with its associated statistical mechanical interpretation which provides a new range of tools to study the learning process, including generalisation and complexity.

We use Turing machines, but mutatis mutandis everything applies to other programming languages. Let $T$ be a Turing machine with tape alphabet $\Sigma$ and set of states $Q$ and assume that on any input $x\in\Sigma^{*}$ the machine eventually halts with output $T(x)\in\Sigma^{*}$. Then to the machine $T$ we may associate the set $\{(x,T(x))\}_{x\in\Sigma^{*}}\subseteq\Sigma^{*}\times\Sigma^{*}$. Program synthesis is the study of the inverse problem: given a subset of $\Sigma^{*}\times\Sigma^{*}$ we would like to determine (if possible) a Turing machine which computes the given outputs on the given inputs.

If we presume given a probability distribution $q(x)$ on $\Sigma^{*}$ then we can formulate this as a problem of statistical inference: given a probability distribution $q(x,y)$ on $\Sigma^{*}\times\Sigma^{*}$ determine the most likely machine producing the observed distribution $q(x,y)=q(y|x)q(x)$. If we fix a universal Turing machine $\mathcal{U}$ then Turing machines can be parametrised by codes $w\in W^{code}$ with $\mathcal{U}(x,w)=T(x)$ for all $x\in\Sigma^{*}$. We let $p(y|x,w)$ denote the probability of $\mathcal{U}(x,w)=y$ (which is either zero or one) so that solutions to the synthesis problem are in bijection with the zeros of the Kullback-Leibler divergence $K(w)$ between the true distribution and the model. So far this is just a trivial rephrasing of the combinatorial optimisation problem of finding a Turing machine $T$ with $T(x)=y$ for all $(x,y)$ with $q(x,y)>0$.

Smooth relaxation. One approach is to seek a smooth relaxation of the synthesis problem consisting of an analytic manifold $W\supseteq W^{code}$ and an extension of $K$ to an analytic function $K:W\rightarrow\mathbb{R}$ so that we can search for the zeros of $K$ using gradient descent. Perhaps the most natural way to construct such a smooth relaxation is to take $W$ to be a space of probability distributions over $W^{code}$ and prescribe a model $p(y|x,w)$ for propagating uncertainty about codes to uncertainty about outputs (Gaunt et al., 2016; Evans & Grefenstette, 2018). Supposing that such a smooth relaxation has been chosen together with a prior $\varphi(w)$ over $W$, smooth program synthesis becomes the study of the statistical learning theory of the triple $(p,q,\varphi)$.

There are at least two reasons to consider the smooth relaxation. Firstly, one might hope that stochastic gradient descent or techniques like Markov chain Monte Carlo will be effective means of solving the original combinatorial optimisation problem. This is not a new idea (Gulwani et al., 2017, §6) but so far its effectiveness for large programs has not been proven. Independently, one might hope to find powerful new mathematical ideas that apply to the relaxed problem and shed light on the nature of program synthesis.

Singular learning theory. All known approaches to program synthesis can be formulated in terms of a singular learning problem. Singular learning theory is the extension of statistical learning theory to account for the fact that the set of true parameters has in general the structure of an analytic space as opposed to an analytic manifold (Watanabe, 2007, 2009). It is organised around triples $(p,q,\varphi)$ consisting of a class of models $\{p(y|x,w):w\in W\}$, a true distribution $q(y|x)$ and a prior $\varphi$ on $W$.

Given a triple arising as above from a smooth relaxation (see Section 3 for the key example) the Kullback-Leibler divergence between the true distribution and the model is

We call points of $W_{0}=\{w\in W\,|\,K(w)=0\}$ solutions of the synthesis problem and note that

where $W_{0}\cap W^{code}$ is the discrete set of solutions to the original synthesis problem. We refer to these as the classical solutions. As the vanishing locus of an analytic function, $W_{0}$ is an analytic space over $\mathbb{R}$ (Hironaka, 1964, §0.1), (Griffith & Harris, 1978) and except in trivial cases it is not an analytic manifold.

The set of solutions fails to be a manifold due to singularities. We say $w\in W$ is a critical point of $K$ if $\nabla K(w)=0$ and a singularity of the function $K$ if it is a critical point where $K(w)=0$. Since $K$ is a Kullback-Leibler divergence it is non-negative and so it not only vanishes on $W_{0}$ but $\nabla K$ also vanishes, hence every point of $W_{0}$ is a singular point. Thus programs to be synthesised are singularities of analytic functions.

The geometry of $W_{0}$ depends on the particular model $p(y|x,w)$ that has been chosen, but some aspects are universal: the nature of program synthesis means that typically $W_{0}$ is an extended object (i.e. it contains points other than the classical solutions) and the Hessian matrix of second order partial derivatives of $K$ at a classical solution is not invertible - that is, the classical solutions are degenerate critical points of $K$. This means that singularity theory is the appropriate branch of mathematics for studying the geometry of $W_{0}$ near a classical solution. It also means that the Fisher information matrix is degenerate at a classical solution, so that the appropriate branch of statistical learning theory is singular learning theory (Watanabe, 2007, 2009). For an introduction to singular learning theory in the context of deep learning see (Murfet et al., 2020).

This geometry has important consequences for all synthesis tasks. In particular, given a dataset $D_{n}=\{(x_{i},y_{i})\}_{i=1}^{n}$ of input-output pairs from $q(x,y)$ and a predictive distribution $p^{*}(y|x,D_{n})=\int p(y|x,w)p(w|D_{n})dw$, the asymptotic form of the Bayes generalization error $B_{g}(n):=D_{KL}(q\|p^{*})$ is given in expectation by $\mathbb{E}[B^{*}_{g}]=\lambda$ where $\lambda$ is a rational number called the real log canonical threshold (RLCT) (Watanabe, 2009). Therefore, the RLCT is the key geometric invariant which should be estimated to determine the generalization capacity of program synthesis. We also show how this number measures the model complexity of program synthesis by relating it to the Kolmogorov complexity. One may interpret the RLCT as a refined measure of complexity.

The synthesis process. Synthesis is a problem because we do not assume that the true distribution is known: for example, if $q(y|x)$ is deterministic and the associated function is $f:\Sigma^{*}\rightarrow Q$, we assume that some example pairs $(x,f(x))$ are known but no general algorithm for computing $f$ is known (if it were, synthesis would have already been performed). In practice synthesis starts with a sample $D_{n}=\{(x_{i},y_{i})\}_{i=1}^{n}$ from $q(x,y)$ with associated empirical Kullback-Leibler distance

If the synthesis problem is deterministic and $u\in W^{code}$ then $K_{n}(u)=0$ if and only if $u$ explains the data in the sense that $\operatorname{step}^{t}(x_{i},u)=y_{i}$ for $1\leq i\leq n$. We now review two natural ways of finding such solutions in the context of machine learning.

Synthesis by stochastic gradient descent (SGD). The first approach is to view the process of program synthesis as stochastic gradient descent for the function $K:W\rightarrow\mathbb{R}$. We view $D_{n}$ as a large training set and further sample subsets $D_{m}$ with $m\ll n$ and compute $\nabla K_{m}$ to take gradient descent steps $w_{i+1}=w_{i}-\eta\nabla K_{m}(w_{i})$ for some learning rate $\eta$. Stochastic gradient descent has the advantage (in principle) of scaling to high-dimensional parameter spaces $W$, but in practice it is challenging to use gradient descent to find points of $W_{0}$ (Gaunt et al., 2016).

Synthesis by sampling. The second approach is to consider the Bayesian posterior associated to the synthesis problem, which can be viewed as an update on the prior distribution $\varphi$ after seeing $D_{n}$

where $Z^{0}_{n}=\int\varphi(w)\exp(-nK_{n}(w))dw$. If $n$ is large the posterior distribution concentrates around solutions $w\in W_{0}$ and so sampling from the posterior will tend to produce machines that are (nearly) solutions. The gold standard sampling is Markov Chain Monte Carlo (MCMC). Scaling MCMC to where $W$ is high-dimensional is a challenging task with many attempts to bridge the gap with SGD (Welling & Teh, 2011; Chen et al., 2014; Ding et al., 2014; Zhang et al., 2020). Nonetheless in simple cases we demonstrate experimentally in Section 7 that machines may be synthesised by using MCMC to sample from the posterior.

To demonstrate how program synthesis may be formalised in the context of singular learning theory we construct a very general synthesis model that can in principle synthesise arbitrary computable functions from input-output examples. This is done by propagating uncertaintly through the weights of a smooth relaxation of a universal Turing machine (UTM) (Clift & Murfet, 2018). Related approaches to program synthesis paying particular attention to model complexity have been given in (Schmidhuber, 1997; Hutter, 2004; Gaunt et al., 2016; Freer et al., 2014).

We fix a Universal Turing Machine (UTM), denoted $\mathcal{U}$, with a description tape (which specifies the code of the Turing machine to be executed), a work tape (simulating the tape of that Turing machine during its operation) and a state tape (simulating the state of that Turing machine). The general statistical learning problem that can be formulated using $\mathcal{U}$ is the following: given some initial string $x$ on the work tape, predict the state of the simulated machine and the contents of the work tape after some specified number of steps.

For simplicity, in this paper we consider models that only predict the final state; the necessary modifications in the general case are routine. We also assume that $W$ parametrises Turing machines whose tape alphabet $\Sigma$ and set of states $Q$ have been encoded by individual symbols in the tape alphabet of $\mathcal{U}$. Hence $\mathcal{U}$ is actually what we call a pseudo-UTM (see Appendix G). Again, treating the general case is routine and for the present purposes only introduces uninteresting complexity.

Let $\Sigma$ denote the tape alphabet of the simulated machine, $Q$ the set of states and let $L,S,R$ stand for left, stay and right, the possible motions of the Turing machine head. We assume that $|Q|>1$ since otherwise the synthesis problem is trivial. The set of ordinary codes $W^{code}$ for a Turing machine sits inside a compact space of probability distributions $W$ over codes

where $\Delta X$ denotes the set of probability distributions over a set $X$, see (H), and the product is over pairs $(\sigma,q)\in\Sigma\times Q$.¹¹1The space $W$ of parameters is clearly semi-analytic, that is, it is cut out of $\mathbb{R}^{d}$ for some $d$ by the vanishing $f_{1}(x)=\cdots=f_{r}(x)=0$ of finitely many analytic functions on open subsets of $\mathbb{R}^{d}$ together with finitely many inequalities $g_{1}(x)\geq 0,\ldots,g_{s}(x)\geq 0$ where the $g_{j}(x)$ are analytic. In fact $W$ is semi-algebraic, since the $f_{i}$ and $g_{j}$ may all be chosen to be polynomial functions. For example the point $\{(\sigma^{\prime},q^{\prime},d)\}_{\sigma,q}\in W^{code}$ encodes the machine which when it reads $\sigma$ under the head in state $q$ writes $\sigma^{\prime}$, transitions into state $q^{\prime}$ and moves in direction $d$. Given $w\in W^{code}$ let $\operatorname{step}^{t}(x,w)\in Q$ denote the contents of the state tape of $\mathcal{U}$ after $t$ timesteps (of the simulated machine) when the work tape is initialised with $x$ and the description tape with $w$. There is a principled extension of this operation of $\mathcal{U}$ to a smooth function

which propagates uncertainty about the symbols on the description tape to uncertainty about the final state and we refer to this extension as the smooth relaxation of $\mathcal{U}$. The details are given in Appendix H but at an informal level the idea behind the relaxation is easy to understand: to sample from $\Delta\operatorname{step}^{t}(x,w)$ we run $\mathcal{U}$ to simulate $t$ timesteps in such a way that whenever the UTM needs to “look at” an entry on the description tape we sample from the corresponding distribution specified by $w$.²²2Noting that this sampling procedure is repeated every time the UTM looks at a given entry.

The class of models that we consider is

where $t$ is fixed for simplicity in this paper. More generally we could also view $x$ as consisting of a sequence and a timeout, as is done in (Clift & Murfet, 2018, §7.1). The construction of this model is summarised in Figure 1.

As $\Delta\operatorname{step}^{t}$ is a polynomial function, $K$ is analytic and so $W_{0}$ is a semi-analytic space (it is cut out of the semi-analytic space $W$ by the vanishing of $K$). If the synthesis problem is deterministic and $q(x)$ is uniform on some finite subset of $\Sigma^{*}$ then $W_{0}$ is semi-algebraic (it is cut out of $W$ by polynomial equations) and all solutions lie at the boundary of the parameter space $W$ (Appendix B). However in general $W_{0}$ is only semi-analytic and intersects the interior of $W$ (Example E.2). We assume that $q(y|x)$ is realisable that is, there exists $w_{0}\in W$ with $q(y|x)=p(y|x,w_{0})$.

A triple $(p,q,\varphi)$ is regular if the model is identifiable, ie. for all inputs $x\in\mathbb{R}^{n}$, the map sending $w$ to the conditional probability distribution $p(y|x,w)$ is one-to-one, and the Fisher information matrix is non-degenerate. Otherwise, the learning machine is strictly singular (Watanabe, 2009, §1.2.1). Triples arising from synthesis problems are typically singular: in Example 3.5 below we show an explicit example where multiple parameters $w$ determine the same model, and in Example E.2 we give an example where the Hessian of $K$ is degenerate everywhere on $W_{0}$ (Watanabe, 2009, §1.1.3).

This important quantity will be estimated in Section 5. Intuitively, the more singular the analytic space $W_{0}$ of solutions is, the smaller the RLCT. One way to think of the RLCT is as a count of the effective number of parameters near $W_{0}$ (Murfet et al., 2020, §4). A thermodynamic interpretation is provided in Section 6. In Section 4 we relate the RLCT to Kolmogorov complexity and in Section 7 we estimate the RLCT of the synthesis problem detectA  using our estimation method.

Every Turing machine is the solution of a deterministic synthesis problem so Section 3 associates to any Turing machine a singularity of a semi-analytic space $W_{0}$. To indicate that this connection is not vacuous, we sketch how the length of a program is related to the real log canonical threshold of a singularity.

Let $q(x,y)$ be a deterministic synthesis problem for $\mathcal{U}$ which only involves input sequences in some restricted alphabet $\Sigma_{input}$, that is, $q(x)=0$ if $x\notin(\Sigma_{input})^{*}$. Let $D_{n}$ be sampled from $q(x,y)$ and let $u,v\in W^{code}\cap W_{0}$ be two explanations for the sample in the sense that $K_{n}(u)=K_{n}(v)=0$. Which explanation for the data should we prefer? The classical answer based on Occam’s razor (Solomonoff, 1964) is that we should prefer the shorter program, that is, the one using the fewest states and symbols.

Set $N=|\Sigma|$ and $M=|Q|$. Any Turing machine $T$ using $N^{\prime}\leq N$ symbols and $M^{\prime}\leq M$ states has a code for $\mathcal{U}$ of length $cM^{\prime}N^{\prime}$ where $c$ is a constant. We assume that $\Sigma_{input}$ is included in the tape alphabet of $T$ so that $N^{\prime}\geq|\Sigma_{input}|$ and define the Kolmogorov complexity of $q$ with respect to $\mathcal{U}$ to be the infimum $\mathfrak{c}(q)$ of $M^{\prime}N^{\prime}$ over Turing machines $T$ that give classical solutions for $q$.

Let $\lambda$ be the RLCT of the triple $(p,q,\varphi)$ associated to the synthesis problem (Definition 3.4).

The proof is deferred to Appendix C.

We have stated that the RLCT $\lambda$ is the most important geometric invariant in singular learning theory (Section 2) and shown how it is related to computation by relating it to the Kolmogorov complexity (Section 4). Now we explain how to estimate it in practice.

Given a sample $D_{n}=\{(x_{i},y_{i})\}_{i=1}^{n}$ from $q(x,y)$ let

be the negative log likelihood, where $S_{n}$ is the empirical entropy of the true distribution. We would like to estimate

where $Z^{\beta}_{n}=\int\varphi(w)\prod_{i=1}^{n}p(y_{i}|x_{i},w)^{\beta}dw$ for some inverse temperature $\beta$. If $\beta=\frac{\beta_{0}}{\log n}$ for some constant $\beta_{0}$, then by Theorem 4 of (Watanabe, 2013),

where $\{U_{n}\}$ is a sequence of random variables satisfying $\mathbb{E}[U_{n}]=0$ and $\lambda$ is the RLCT. In practice, the last two terms often vary negligibly with $1/\beta$ and so $\mathbb{E}^{\beta}_{w}[nL_{n}(w)]$ approximates a linear function of $1/\beta$ with slope $\lambda$ (Watanabe, 2013, Corollary 3). This is the foundation of the RLCT estimation procedure found in Algorithm 1 which is used in our experiments.

Each RLCT estimate $\hat{\lambda}(\mathcal{D}_{n})$ in Algorithm 1 was performed by linear regression on the pairs $\{(1/\beta_{i},\mathbb{E}^{\beta_{i}}_{w}[nL_{n}(w)])\}_{i=1}^{5}$ where the five inverse temperatures $\beta_{i}$ are centered on the inverse temperature $1/T$ where $T$ is the temperature reported for each experiment in Table 1 and Table 3.

Now that we have formulated program synthesis in the setting of singular learning theory, we introduce the thermodynamic dictionary that this makes available. Given a triple $(p,q,\varphi)$ associated to a synthesis problem, the physical system consists of the space $W$ of microstates and the random Hamiltonian

which depends on a sample $D_{n}$ from the true distribution of size $n$ and the inverse temperature $\beta$. The corresponding Boltzmann distribution is the tempered Bayesian posterior (2) for the synthesis problem

where $Z=\int_{W}\exp(-\beta H_{n}(w))dw$ is called the partition function. In thermodynamics we think of functions $f(w)$ of microstates as fundamentally unobservable; what we can observe are averages

possibly restricted to subsets $W^{\prime}\subseteq W$ (microstates compatible with some specified macroscopic condition). Taking program synthesis seriously leads us to systematically re-evaluate aspects of the theory of computation with this restriction (that only integrals are observable) in mind.

To give one important example: the energy of the Boltzmann distribution (tempered Bayesian posterior) at a fixed inverse temperature $\beta$ and value of $n$ is by definition

which is a random variable. The quantity $\frac{1}{\beta}\mathbb{E}^{\beta}_{w}\big{[}-\log\varphi(w)\big{]}$ is the contribution to the energy from interaction with the prior, and by (7) the most important contribution to the first summand is $\lambda T$ where $T=\frac{1}{\beta}$. We see the fundamental role that the RLCT plays in the thermodynamic picture.

We estimate the RLCT for the triples $(p,q,\varphi)$ associated to the synthesis problems detectA  (Example 3.5) and parityCheck  (Appendix F). Hyperparameters of the various machines are contained in Table 2 of Appendix D. The true distribution $q(x)$ is defined as follows: we fix a minimum and maximum sequence length $a\leq b$ and to sample $x\sim q(x)$ we first sample a length $l$ uniformly from $[a,b]$ and then uniformly sample $x$ from $\{A,B\}^{l}$.

We perform MCMC on the weight vector for the model class $\{p(y|x,w):w\in W\}$ where $w$ is represented in our PyTorch implementation by three tensors of shape $\{[L,n_{i}]\}_{1\leq i\leq 3}$ where $L$ is the number of tuples in the description tape of the TM being simulated and $\{n_{i}\}$ are the number of symbols, states and directions respectively. A direct simulation of the UTM is used for all experiments to improve computational efficiency (Appendix I). We generate, for each inverse temperature $\beta$ and dataset $D_{n}$, a Markov chain via the No-U-Turn sampler from (Hoffman & Gelman, 2014). We use the standard uniform distribution as our prior $\varphi$.

For the problem detectA  given in Example 3.5 the dimension of parameter space is $\dim W=30$. We use generalised least squares to fit the RLCT $\lambda$ (with goodness-of-fit measured by $R^{2}$), the algorithm of which is given in Section 5. Our results are displayed in Table 1 and Figure 3. Our purpose in these experiments is not to provide high accuracy estimates of the RLCT, as these would require much longer Markov chains. Instead we demonstrate how rough estimates consistent with the theory can be obtained at low computational cost. If this model were regular the RLCT would be $\dim W/2=15$.

We have developed a theoretical framework in which all programs can in principle be synthesised from input-output examples. This is done by propagating uncertainty through a smooth relaxation of a universal Turing machine. In approaches to program synthesis based on gradient descent there is a tendency to think of solutions to the synthesis problem as isolated critical points of the loss function $K$, but this is a false intuition based on regular models.

Since neural networks, Bayesian networks, smooth relaxations of UTMs and all other extant approaches to smooth program synthesis are strictly singular models (the map from parameters to functions is not injective) the set $W_{0}$ of parameters $w$ with $K(w)=0$ is a complex extended object, whose geometry is shown by Watanabe’s singular learning theory to be deeply related to the learning process. We have examined this geometry in several specific examples and shown how to think about programs from a geometric and thermodynamic perspective. It is our hope that this approach can assist in developing the next generation of synthesis machines.