The authors present empirical distributions for the halting time (measured by the number of iterations to reach a given accuracy) of optimization algorithms applied to two random systems: spin glasses and deep learning. Given an algorithm, which we take to be both the optimization routine and the form of the random landscape, the fluctuations of the halting time follow a distribution that, after centering and scaling, remains unchanged even when the distribution on the landscape is changed. We observe two main classes, a Gumbel-like distribution that appears in Google searches, human decision times, QR factorization and spin glasses, and a Gaussian-like distribution that appears in conjugate gradient method, deep network with MNIST input data and deep network with random input data. This empirical evidence suggests presence of a class of distributions for which the halting time is independent of the underlying distribution under some conditions.
1 INTRODUCTION
In this paper we discuss both the presence and application of universality in optimization algorithms. More precisely, in order to optimize an energy functional when the functional itself and the initial guess are random, we consider the following iterative algorithms: conjugate gradient for solving a linear system, gradient descent for spin glasses, and stochastic gradient descent for deep learning.
A bounded, piecewise differentiable random field (See Adler & Taylor (2009) for an account on the connection of random fields and geometry), where the randomness is non-degenerate, yields a landscape with many saddle points and local minima. Given such a landscape and a moving particle that takes steps to reach a low-energy level, an essential quantity is the time the particle takes until it stops which we call the halting time. Many useful bounds on the halting time are known for convex cases, where the stopping condition produces a halting time that is, essentially, the time to find the minimum. In non-convex cases, however, the particle knows only the information that can be calculated locally. And a locally measurable stopping condition, such as the norm of the gradient at the present point, or the difference in altitude with respect to the previous step, can lead the algorithm to locate a local minimum. This feature allows the halting time to be calculated in a broad range of non-convex, high-dimensional problems. A prototypical example of such a random field is the class of polynomials with random coefficients. Spin glasses and deep learning cost functions are then special cases of such fields that yield different landscapes. Polynomials with random coefficients are not only a broad class of functions, but also they are hard to study mathematically in any generality. Therefore, in order to capture essential features of such problems, we focus on their subclasses that are well studied (spin glasses) and practically relevant (deep learning cost functions).
The halting time in such landscapes, when normalized to mean zero and variance one (subtracting the mean and dividing by the standard deviation), appears to follow a distribution that is independent of the input data, in other words it follows a universal distribution: the fluctuations are universal. In statistical mechanics, the term “universality” is used to refer to a class of systems which, on a certain macroscopic scale, behave statistically the same while having different statistics on a microscopic scale. An example of such a law is the central limit theorem, which states that the sums of observations tend to follow the same distribution independent of the distribution of the individual observations, as long as contribution from individual observations is reasonably small. It may fail to hold, if the microscopic behavior is not independent, does not have a finite second-moment, or if we consider something different than the sum. This work’s focus is an attempt to put forward the cases
where we see universality. But in this spirit, we show a degenerate case in which halting time fails to follow a universal law.
A rather surprising example of halting time universality is in the cases of observed human decision times and GoogleTM query times. In Bakhtin & Correll (2012) the time it takes a person make a decision in the presence of visual stimulus is shown to have universal fluctuations. The theoretically predicted curve in this experiment follows a Gumbel-like distribution. In addition, we randomly sampled words from two different dictionaries and submitted search queries. The time it takes Google to present the results are recorded. The normalized search times closely follow the same Gumbel-like curve.
In the cases we observe, we find two main universality classes: (1) A Gumbel-like distribution that appears in Google searches, human decision times, QR factorization and spin glasses, and (2) a Gaussian-like distribution that appears in conjugate gradient algorithm and deep learning. To the best of our knowledge, our work along with the accompanying references in this introduction are the first ones to address the question of observing and classifying the distribution of the halting time.
1.1 DEFINITION OF UNIVERSALITY
Definition 1.1. An algorithm A consists of both a random cost function F (x, w) where x is a given random input and an optimization routine that seeks to minimize F with respect to w.
To each algorithm we attach a precise -dependent halting criteria for the algorithm. The halting time, which is a random variable, is the time it takes to meet this criteria. Within each algorithm there must be an intrinsic notion of dimension which we denote by N . The halting time T ,N,A,E depends on , N , the choice of algorithm A, and the ensemble E (or probability distribution). We use the empirical distribution of T ,N,A,E to provide heuristics for understanding the qualitative performance of the algorithms.
The presence of universality in an algorithm is the observation that for sufficiently large N and = (N), the halting time random variable satisfies
τ ,N,A,E := T ,N,A,E − E[T ,N,A,E ]√
Var(T ,N,A,E) ≈ τ∗A, (1)
where τ∗A is a continuous random variable that depends only on the algorithm. The random variable τ ,N,A,E is referred to as the fluctuations and when such an approximation appears to be valid we
say that N and (and any other external parameters) are in the scaling region. For example, in Section 1.2, A is the QR eigenvalue algorithm, N is the size of the matrix, is a small tolerance and E is given by a distribution on complex Hermitian (or real symmetric) matrices.
Some remarks must be made:
• A statement like (1) is known to hold rigorously for a few algorithms (see Deift & Trogdon (2016; 2017)) but in practice, it is verified experimentally. This was first done in Pfrang et al. (2014) and expanded in Deift et al. (2014) for a total of 8 different algorithms.
• The random variable τ∗A depends fundamentally on the functional form of F . And we only expect (1) to hold for a restricted class of ensembles E.
• T ,N,A,E is an integer-valued random variable. For it to become a continuous distribution limit must be taken. This is the only reason N must be large — in practice, the approximation in (1) is seen even for small to moderate N .
Universality in this sense is a measure of stability in an algorithm. For example, it is known from the work of Kostlan (1988) that halting time for the power method to compute the largest eigenvalue (in modulus) of symmetric Gaussian matrices has infinite expectation and hence this type of universality is not believed to be present. One could use this to conclude that the power method is a naı̈ve method for these matrices. Yet, it is known that the power method is much more efficient on positive-definite matrices where universality can be shown Deift & Trogdon (2017). Therefore, we have evidence that the presence of universality is a desirable feature of a numerical method.
1.2 EXAMPLE: DEMONSTRATION OF UNIVERSALITY IN THE QR ALGORITHM
To give some context, we discuss the universality in the solution of the eigenvalue problem with the classical QR algorithm. Historically, this was first noticed in Pfrang et al. (2014). In this example the fundamental object is the QR factorization (Q,R) = QR(A) where A = QR, Q is orthogonal (or unitary) and R is upper-triangular with positive diagonal entries. The QR algorithm applied to a Hermitian N ×N matrix A is given by the iteration
A0 := A,
(Qj , Rj) := QR(Aj),
Aj+1 := RjQj .
Generically, Aj → D as j → ∞ where D is a diagonal matrix whose diagonal entries are the eigenvalues of A. The halting time in Pfrang et al. (2014) was set to be the time of first deflation, T ,N,A,E(A), as:
min{j : √ N(N − k)‖Aj(k + 1 : N, 1 : k)‖∞ <
for some 1 ≤ k ≤ N − 1}.
Here ‖A‖∞ refers to the maximum entry of a matrix A in absolute value and the notation A(i : j, k : l) refers to the submatrix of A consisting of entries only in rows i, i+1, . . . , j and in columns k, k + 1, . . . , l. Thus the halting time for the QR algorithm is the time at which at least one offdiagonal block is appropriately small. Next, we have to discuss choices for the randomness, or ensembleE, by choosing different distributions on the entries ofA. Four such choices for ensembles are, Bernoulli ensemble (BE), Gaussian orthogonal ensemble (GOE), Gaussian unitary ensemble (GUE), Quartic unitary ensemble (QUE):
BE A is real-symmetric with iid Bernoulli ±1 entries on and below the diagonal. GOE A is real-symmetric with iid standard normal entries below the diagonal. The entries on the diagonal are iid normal with mean zero and variance two. GUE A is complex-Hermitian with iid standard complex normal entries below the diagonal. The
entries on the diagonal are iid complex normal mean zero and with variance two.
QUE A is complex-Hermitian with probability density ∝ e−trA4dA. See Deift (2000) for details on such an ensemble and Olver et al. (2015) for a method to sample such a matrix. Importantly, the entries of the matrix below the diagonal are correlated.
Here we have continuous and discrete, real and complex, and independent and dependent ensembles but nevertheless we see universality in Figure 2 where we take N = 150 and = 10−10.
Remark 1.1. The ensembles discussed above (GOE, GUE, BE and QUE) exhibit eigenvalue repulsion. That is, the probability that two eigenvalues are close1 is much smaller than if the locations of the eigenvalues were just given by iid points on the line. It turns out that choosing a random matrix with iid eigenvalues breaks the universality that is observed in Figure 2. See Pfrang et al. (2014) for a more in-depth discussion of this. Remark 1.2. To put the QR algorithm in the framework, let B = UAU∗ define F (A,U) by
min{j : √ N(N − k)‖B(k + 1 : N, 1 : k)‖∞ <
for some 1 ≤ k ≤ N − 1} We then use the QR algorithm to minimize F with respect to unitary matrices U using the initial condition U = I . If A is random then F (A,U) represents a random field on the unitary group.
1.3 CORE EXAMPLES: SPIN GLASS HAMILTONIANS AND DEEP LEARNING COST FUNCTIONS
A natural class of random fields is the class of Gaussian random functions on a high-dimensional sphere, known as p-spin spherical spin glass models in the physics literature (in the Gaussian process literature they are known as isotropic models). From the point of view of optimization, minimizing the spin glass model’s Hamiltonian is fruitful because a lot is known about its critical points. This allows us to experiment with questions regarding whether the local minima and saddle points, due to the non-convex nature of landscapes, present an obstacle in the training of a system. Such observations on the Hamiltonian doesn’t imply that it is a cost function or a simplified version of a cost function. Rather, the features that both systems have in common hint at a deeper underlying structure that needs to be discovered.
In recent years Dauphin et al. (2014) attacked the saddle point problem of non-convex optimization within deep learning. In contrast, Sagun et al. (2014) and the experimental second section of Choromanska et al. (2014) jointly argue that if the system is large enough, presence of saddle points is not an obstacle, and add that the local minimum practically gives a good enough solution within the limits of the model. However, Sagun et al. (2014) and Choromanska et al. (2014) hold different perspectives on what the qualitative similarities between optimization in spin glasses and deep learning might imply. The latter asserts a direct connection between the two systems based on these similarities. On the contrary, the former argues that these similarities hint at universal behaviors that are generically observed in vastly different systems rather than emphasizing a direct connection.
1By close, we mean that their distance is much less than O(1/N) where N is the size of the matrix.
In line with the asymptotic proof in Auffinger et al. (2013), the local minima are observed to lie roughly at the same energy level in spherical spin glasses. Auffinger et al. (2013) also gives asymptotic bounds on the value of the ground state and the exponential behavior of the average of the number of critical points below a given energy level. It turns out, when the dimension is large, the bulk of the local minima tend to have the same energy which is slightly above the global minimum. This level is called the floor level of the function. Simulations of the floor in spin glass can be found in Sagun et al. (2014). Sagun et al. (2014) also exhibits floor in a specially designed MNIST experiment: A student network is trained by the outputs of a pre-trained teacher network. Zero cost is achievable by the student, but the stochastic gradient descent cannot find zeros. It also does not have to because the floor level already gives a decent performance.
• Given data (i.e., from MNIST) and a measure L(x`, w) for determining the cost that is parametrized by w ∈ RN , the training procedure aims to find a point w∗ that minimizes the empirical training cost while keeping the test cost low. Here x` ∈ Z for ` ∈ {1, ..., S}, where Z is a random (ordered) sample of size S from the training examples. Total training cost is given by
F (Z,w) = LTrain(w) = 1
S S∑ `=1 L(x`, w). (2)
• Given couplings x(·) ∼ Gaussian(0, 1) that represent the strength of forces between triplets of spins. The state of the system is represented by w ∈ SN−1( √ N) ⊂ RN .
The Hamiltonian (or energy) of the simplest complex2 spherical spin glass model is given by:
F (x(·), w) = HN (w) = 1
N N∑ i,j,k xijkwiwjwk. (3)
The two functions are indeed different in two major ways. First, the domain of the Hamiltonian is a compact space and the couplings are independent Gaussian random variables whereas the inputs for (2) are not independent and the cost function has a non-compact domain. Second, at a fixed point w, variance of the function LTrain(w) is inversely proportional to the number of samples, but the variance of HN (w) is N . As a result a randomly initialized Hamiltonian can take vastly different values, but a randomly initialized cost tend to have very similar values. The Hamiltonian has macroscopic extensive quantities: its minimum scales with a negative constant multiple ofN . In contrast, the minimum of the cost function is bounded from below by zero. All of this indicates that landscapes with different geometries (glass-like, funnel-like, or another geometry) might still lead to similar phenomena such as existence of the floor level, and the universal behavior of the halting time.
1.4 SUMMARY OF RESULTS
We discuss the presence of universality in algorithms that are of a very different character. The conjugate gradient algorithm, discussed in Section 2.1, effectively solves a convex optimization problem. Gradient descent applied in the spin glass setting (discussed in Section 2.2) and stochastic gradient descent in the context of deep learning (MNIST, discussed in Section 2.3) are much more complicated non-convex optimization processes. Despite the fact that these algorithms share very little geometry in common, we demonstrate three things they share:
• A scaling region in which universality appears and performance is good. • Regions where the computation is either ineffective or inefficient. • A moment-based indicator for finding the universality class.
22-spin spherical spin glass, sum of xijwiwj terms, has exactly 2N critical points. When p ≥ 3, p−spin model has exponentially many critical points with respect to N . For the latter case, complexity is a measure on the number of critical points in an exponential scale. Deep learning problems are suspected to be complex in this sense.
2 EMPIRICAL OBSERVATION OF UNIVERSALITY
2.1 THE CONJUGATE GRADIENT ALGORITHM
The conjugate gradient algorithm (Hestenes & Stiefel, 1952) for solving the N × N linear system Ax = b, when A = A∗ is positive definite, is an iterative procedure to find the minimum of the convex quadratic form:
F (A, y) = 1
2 y∗Ay − y∗b,
where ∗ denotes the conjugate-transpose operation. Given an initial guess x0 (we use x0 = b), compute r0 = b−Ax0 and set p0 = r0. For k = 1, . . . , N ,
1. Compute rk = rk−1 − ak−1Apk−1 where3 ak−1 = 〈rk−1, rk−1〉/〈pk−1, Apk−1〉. 2. Compute pk = rk + bk−1pk−1 where bk−1 = 〈rk, rk〉/〈rk−1, rk−1〉. 3. Compute xk = xk−1 + ak−1pk−1.
If A is strictly positive definite xk → x = A−1b as k → ∞. Geometrically, the iterates xk are the best approximations of x over larger and larger affine Krylov subspaces Kk,
‖Axk − b‖A = minx∈Kk‖Ax− b‖A, Kk = x0 + span{r0, Ar0, . . . , Ak−1r0}, ‖x‖2A = 〈x,A−1x〉,
as k ↑ N . The quantity one monitors over the course of the conjugate gradient algorithm is the norm ‖rk‖:
T ,N,CG,E(A, b) := min{k : ‖rk‖ < }.
In exact arithmetic, the method takes at mostN steps: In calculations with finite-precision arithmetic the number of steps can be much larger than N and the behavior of the algorithm in finite-precision arithmetic has been the focus of much research (Greenbaum, 1989; Greenbaum & Strakos, 1992). What is important for us here is that it may happen that ‖rk‖ < but the true residual r̂k := b−Axk (which typically differs from rk in finite-precision computations) satisfies ‖r̂k‖ > .
Now, we discuss our choices for ensembles E of random data. In all computations, we take b = (bj)1≤j≤N where each bj is iid uniform on (−1, 1). We construct positive definite matrices A by A = XX∗ where X = (Xij)1≤i≤N, 1≤j≤M and each Xij ∼ D is iid for some distribution D. We make the following three choices for D, Positive definite Bernoulli ensemble (PDE), Laguerre orthogonal ensemble (LOE), Laguerre unitary ensemble (LUE):
3We use the notation ‖y‖2 = 〈y, y〉 = ∑
j |yj | 2 for y = (y1, y2, . . . , yN ) ∈ CN .
PBE D a Bernoulli ±1 random variable (equal probability). LOE D is a standard normal random variable. LUE D is a standard complex normal random variable.
The choice of the integer M , which is the inner dimension of the matrices in the product XX∗, is critical for the existence of universality. In Deift et al. (2014) and Deift et al. (2015) it is demonstrated that universality is present when M = N + bc √ Nc and the -accuracy is small, but fixed. Universality is not present when M = N and this can be explained by examining the distribution of the condition number of the matrix A in the LUE setting (Deift et al., 2015). We demonstrate this again in Figure 3(a). We also demonstrate that universality does indeed fail for M = N in Figure 3(b).
2.2 SPIN GLASSES AND GRADIENT DESCENT
The gradient descent algorithm for the Hamiltonian of the p-spin spherical glass will find a local minimum of the non-convex function (3). Since variance of HN (w) is typically of order N , a local minimum has size N . More precisely, by Auffinger et al. (2013), the energy of the floor level where most of local minima are located is asymptotically at −2 √ 2/3N ≈ −1.633N and the ground state is around−1.657N . The algorithm starts by picking a random element w of the sphere with radius √ N , SN−1( √ N), as a starting point for each trial. We vary the environment for each trial and introduce ensembles by setting x(·) ∼ D for a number of choices of distributions. For a fixed dimension N , accuracy that bounds the norm of the gradient, and an ensemble E: (1) Calculate the gradient steps: wt+1 = wt − ηt∇wH(wt), (2) Normalize the resulting vector to the sphere: √ N w t+1
||wt+1|| ← wt+1, and (3) Stop when the norm of the gradient size is below and record T ,N,GD,E . This procedure is repeated 10,000 times for different ensembles (i.e. different choices forD). Figure 4 exhibit the universal halting time presenting evidence that τ ,N,GD,E is independent of the ensemble.
2.3 DIGIT INPUTS VS. RANDOM INPUTS IN DEEP LEARNING
A deep learning cost function is trained on two drastically different ensembles. The first is the MNIST dataset, which consists of 60,000 samples of training examples and 10,000 samples of test examples. The model is a fully connected network with two hidden layers, that have 500 and 300 units respectively. Each hidden unit has rectified linear activation, and a cross entropy cost is attached at the end. To randomize the input data we sample 30K samples from the training set each time we set up the model and initialize the weights randomly. Then we train the model by the stochastic gradient descent method with a minibatch size of 100. This model gets us about 97%
accuracy without any further tuning. The second ensemble uses the same model and outputs, but the input data is changed from characters to independent Gaussian noise. This model, as expected, gets us only about 10% accuracy: it randomly picks a number! The stopping condition is reached when the average of successive differences in cost values goes below a prescribed value. As a comparison we have also added a deep convolutional network (convnet), and we used the fully connected model with a different stopping condition: one that is tied to the norm of the gradient. Figure 5 demonstrates universal fluctuations in the halting time in all of the four cases.
3 CONCLUSIONS
What are the conditions on the ensembles and the model that lead to such universality? What constitutes a good set of hyperparameters for a given algorithm? How can we go beyond inspection when tuning a system? How can we infer if an algorithm is a good match to the system at hand? What is the connection between the universal regime and the structure of the landscape? This research attempts to exhibit cases where one can extract answers to these questions in a robust and quantitative way. The examples we have presented clearly exhibit universality. The normalized moment analysis, presented in the Appendix, gives a quantitative way to test for universality. And we further believe that an algorithm that exhibits universality is running in a scaling region of “high performance”: universality is a measure of insensitivity to initial data which is a beneficial property of a numerical method. Establishing this claim is a difficult task, beyond the scope of this primarily empirical work.
More specifically, the current work gives empirical evidence that within an appropriate scaling region, the halting time can often be approximated as
T ,N,A,E ≈ µ+ στ∗A, where τ∗A is a mean-zero, variance one universal distribution. If this holds, a simple estimate of the mean µ = µ ,N,A,E and the standard deviation σ = σ ,N,A,E using a few samples will give a good a priori estimate of algorithm run time
P (|T ,N,A,E − µ| ≥ σ`) ≈ P (|τ∗A| ≥ `) . If τ∗A has (or is just conjectured to have) exponential tails, for example, this can be quite useful.
This work also validates the broad claims made in Deift et al. (2015) that universality is present in all or nearly all (sensible) computation. Future work will be along the lines of using these heuristics to identify when we have universality, to identify the different kinds of landscapes, and to guide both algorithm development and algorithm tuning. Furthermore, one would like theoretical estimates for the mean µ ,N,A,E and the standard deviation σ ,N,A,E .
ACKNOWLEDGMENTS
We thank Percy Deift for valuable discussions and Gérard Ben Arous for his mentorship throughout the process of this research. The first author thanks very much to Uğur Güney for his availability for support and valuable contributions in countless implementation issues. This work was partially supported by the National Science Foundation under grant number DMS-1303018 (TT).
Summary
For several algorithms, previous research has shown that the halting time follows a two-parameter distribution (the so-called universal property investigated by the authors). In this work, the authors extend the investigation to new algorithms (spin-glass, gradient descent in deep learning).
An algorithm is considered to satisfy the universality property when the centered/scaled halting time fluctuations (empirical distribution of halting times) depend on the algorithm but do not depend on the target accuracy epsilon, an intrinsic measure of dimension N, the probability distribution/random ensemble. (This is clear from Eq 1 where on the left the empirical halting time distribution depends on epsilon, N, A, E and on the right, the approximation only depends on the algorithm)
The authors argue that empirically, the universal property is observed when both algorithms (spin glass and deep learning) perform well and that it is not observed when they do not perform well.
A moment-based indicator is introduced to assess whether universality is observed.
Review
This paper presents several problems.
page 2: “[…] for sufficiently large N and eps = eps(N)”
The dependence of epsilon on N is troubling.
page 3: “Universality is a measure of stability in an algorithm […] For example […] halting time for the power method […] has infinite expectation and hence this type of universality is *not* present. One could use this to conclude that the power method is naive. Therefore the presence of universality is a desirable feature of a numerical method”
No. An algorithm is naive if there are better ways to answer the problem. One could not conclude from a halting time with infinite expectation (e.g. solving a problem extremely quickly 99% of the time, and looping forever in 1% of cases) or infinite variance, that the algorithm is naive. 
Moreover the universal property is more restrictive than having a finite halting time expectation. Even if in many specific cases, having a finite halting time expectation is a desirable property, showing that the presence of universality is desirable would require a demonstration that the other more restrictive aspects are also desirable. 
Also, the paragraph only concerns one algorithm. why would the conclusions generalise to all numerical methods ? 
Even if the universality property is arguably desirable (i.e. event if the conclusion of this paragraph is assumed correct), the paragraph does not support the given conclusion.
Comparing Eq 1 and figures 2,3,4,5 
From Eq 1, universality means that the centered/scaled halting time fluctuations (which depend on A, epsilon, N, E) can be approximated by a distribution that only depends on A (not on epsilon, N, E) but in the experiments only E varies (figures 2,3,4,5). The validity of the approximation with varying epsilon or N is never tested
The ensembles/distributions parameter E (on which halting fluctuations should not depend) and the algorithm A (on which halting fluctuations are allowed to depend) are not well defined, especially w.r.t. the common use of the words. In the optimisation setting we are told that the functional form of the landscape function is part of A (in answer to the question of a reviewer) but what is part of the functional form ? what about computations where the landscape has no known functional form (black box) ?
The conclusion claims that the paper “attempts to exhibit cases” where one can answer 5 questions in a robust and quantitative way.
Question 1: “What are the conditions on the ensembles and the model that lead to such universality ?” The only quantitative way would be to use the moments based indicator however there is only one example of universality not being observed which concerns only one algorithm (conjugate gradient) and one type of failure (when M = N). This does not demonstrate robustness of the method.
Question 2: “What constitutes a good set of hyper parameters for a given algorithm ?”
The proposed way to choose would be to test whether universality is observed. If it is then the hyper parameters are good, if not the hyper parameters are bad. The correspondance between bad hyper-parameters and observing no universality concerns only one algorithm and one type of failure. Other algorithms may fail in the universal regime or perform well in the non universal regime. The paper does not show how to answer this question in a robust way.
Question 3: "How can we go beyond inspection when tuning a system ? "
The question is too vague and general and there is probably no robust and quantitative way to answer it at all.
Question 4: "How can we infer if an algorithm is a good match to the system at hand ? "
The paper fails to demonstrate convincingly that universality is either a good or robust way to approach the very few studied algorithms. The suggested generalisation to all systems and algorithms is extremely far fetched.
Question 5: "What is the connection between the universal regime and the structure of the landscape ?"
 Same as before, the question is extremely vague and cannot be answered in a robust or quantitative way at all. The fact that what corresponds to A and what corresponds to E is not clear does not help.
In the conclusion it is written that the paper validates the claim that universality is present in all or nearly all sensible computation. It does not. The paper does not properly test whether universality is present (only 1 parameter in 3 that should not vary is tested). The paper does not properly test whether universality is lost when the computation is no longer sensible (only one failure case tested). Finally the experiments do not apply to all or nearly all computations but only to very few  specific algorithms.