Navigator Function for the Conformal Bootstrap

Current numerical conformal bootstrap techniques carve out islands in theory space by repeatedly checking whether points are allowed or excluded. We propose a new method for searching theory space that replaces the binary information"allowed"/"excluded"with a continuous"navigator"function that is negative in the allowed region and positive in the excluded region. Such a navigator function allows one to efficiently explore high-dimensional parameter spaces and smoothly sail towards any islands they may contain. The specific functions we introduce have several attractive features: they are everywhere well-defined, can be computed with standard methods, and evaluation of their gradient is immediate due to an SDP gradient formula that we provide. The latter property allows for the use of efficient quasi-Newton optimization methods, which we illustrate by navigating towards the 3d Ising island.


Introduction and summary
Over the last decade, the numerical conformal bootstrap program 1 has relied on the idea [4] that for any point in CFT parameter space it is possible to check if the point is allowed or excluded by constructing positive linear functionals. In this work we will dramatically upgrade this idea, replacing the binary information "allowed/excluded" by a continuous measure of success, called a "navigator function." For excluded points, the navigator function will tell us how far we are from the allowed region. Minimizing the navigator, we will be able to quickly find the allowed region, starting from an excluded point. For allowed points, the navigator will tell us how far inside the allowed region we are, and navigator minima will be excellent predictors for the position of an actual CFT.
To describe what we have in mind in some detail, let X be an infinite-dimensional vector containing all parameters characterizing a CFT (i.e. all operator dimensions and OPE coefficients, bundled together). We split it as X = (x, y) where x ∈ R k are parameters we are especially interested in, and y contains all the rest. We also select a finite subset of the infinitely many bootstrap equations.
Most bootstrap computations performed so far proceeded in what one may call "oracle mode." 2 One picks a sequence of trial vectors x 1 , x 2 , . . . and asks for each of them if there is any y such that X = (x i , y) satisfies the selected subset of bootstrap equations. A bootstrap solver such as SDPB [5,6] provides an answer: "allowed" or "excluded". By trying many x i 's, one maps out the allowed region. 3 Thus, we compute the characteristic function χ R of the allowed region R (i.e. χ R (x) = 1 for x ∈ R and χ R (x) = 0 otherwise). Experience shows that the boundary of the allowed region ∂R is typically smooth, apart from isolated points (kinks). This can guide the choice of future trial points and speed up the computation. 4 By trying many points, one zooms in on the boundary ∂R of the allowed region. Importantly, a single oracle query does not provide any information about whether one is close to or far from ∂R. Rather, one knows that one is close to ∂R if one can find two nearby trial points x i and x i such that they are on two different sides of the boundary.
We will modify this setup so that a single SDPB run computes a continuous function N (x), called a navigator, which will give a more nuanced measure of success than simply "allowed/excluded." To be maximally useful, the navigator should have the following properties: • N (x) is continuous and differentiable; 1 Ssee [1] for a thorough review, and [2,3] for pedagogical introductions. 2 In technical jargon referred to as "feasibility mode." 3 Other typical bootstrap computations are OPE coefficient optimizations. Sometimes these computations allow to zoom in on actual CFTs, as e.g. c-minimization is conjectured to lead to the 3d Ising CFT [7]. 4 Other speed-up tricks include the cutting surface algorithm [8], which allows in some cases to use a single oracle computation to rule out not just one point but a large swath of the parameter space.
• N (x) > 0 outside the allowed region R, and N (x) < 0 inside R. In particular, N (x) = 0 on the boundary ∂R; 5 • N (x) should be defined not just in a tiny neighborhood of the allowed region but globally; • The allowed region R should be a basin of attraction of the navigator function from a sizable neighborhood of R.
Assuming these nice properties, the navigator value will allow us to guess how far we are from the allowed region. We will also be able to reach the allowed region by starting from some initial trial point x 0 and by minimizing the navigator until we reach a point with negative N (x). We'd like to be optimistic and hope that the navigator has no local minima away from the allowed region where such a search may get stuck. The idea of replacing the binary information of "oracle mode" with continuous information from solving an optimization problem is not completely new [9,10]. Notably Ref. [10] emphasized the power of this idea to quickly determine the boundary of the allowed region once its approximate position is known, replacing bisection with the secant method. 6 A crucial difference here is our requirement that the navigator should be defined in a wide region and not only near the boundary, which will greatly increase the list of potential applications. This requirement is non-trivial and the early navigator avatars [9,10] don't satisfy it (see Section 2.1.3).
In this paper we will lay down the systematic theory of navigator functions by showing three important results: 1. First, we will show that navigators satisfying all of the above properties can indeed be found for a generic bootstrap problem. We will present both the general principle of their existence, and several explicit constructions (see Section 2). Please scroll down to Fig. 1 for a concrete navigator example in the mixed σbootstrap setup used to isolate the 3d Ising model. It has all the nice properties, and in particular a single minimum (within the range we show), located within the 3d Ising island. See Section 3 for more beautiful navigator plots.
2. Our navigators can be evaluated using standard conformal bootstrap software such as SDPB. In practical applications that we have in mind, it's important to know not just the navigator but also its gradient. Our second important result is a general "SDP gradient formula," Eq. (4.16). This formula shows that navigator gradient can be evaluated essentially for free once the navigator value has been computed using SDPB. 3. We foresee that one of the most important navigator applications will be to quickly look for allowed points, i.e. to "sail towards the Ising island," by minimizing the navigator. Naive minimization strategies, such as the gradient descent, are inefficient, getting stuck in narrow "valleys" of the navigator surface. Our third important result is to demonstrate how a quasi-Newton method-the BFGS algorithm [11]-successfully overcomes these difficulties (Section 5). This algorithm finds first the allowed region, and then the navigator minimum, in a relatively small number of steps.
The paper is structured as follows. Section 2 will explain our two main navigator constructions: the GFF-navigator and the Σ-navigator. (A third construction is in App. B). In Section 3 we will show various plots of these navigators, to gain intuition about their shape. In Section 4 we will derive the SDP gradient formula. In Section 5 we will describe the BFGS algorithm and its bounding-box modification, to look for an allowed point and the navigator minimum, and show that it performs well in realistic multiple-correlator setups. In Section 6 we describe another possible navigator application: extremizing operator dimension within the allowed region. This represents an attractive alternative to the tiptop algorithm recently introduced for this purpose in the feasibility setup [12]. In Section 7 we conclude. Appendix C shows how one can also evaluate the navigator Hessian, in addition to the gradient, provides numerical tests of these procedures.

Navigator function
Our motivation to look for the navigator function, and its desired properties, have already been described in the introduction. The crucial requirement is that the navigator should be finite. Indeed, a navigator which is negative inside the allowed region and equals +∞ outside would be rather useless for the purposes we have in mind, such as looking for an allowed point starting from an excluded one. Furthermore, once a finite navigator is constructed, other nice properties turn out to also be satisfied.
How to get a robustly finite navigator is one of the main ideas of our paper (see Section 2.1.3 for an account of naive attempts which fail). Although the idea is general, we will start in Section 2.1 by presenting it in the simplest single-correlator setup. We will then move on to more realistic multiple-correlator problems.

Single-correlator problems
Consider the simplest bootstrap setup: scalar gap maximization in a single 4pt function of four identical scalars [4]. Thus we are solving the bootstrap equation  (v, u). Here ∆ φ is the external scalar dimension which for simplicity is considered fixed (although see footnote 7). The set S(∆ * ) is given by: S(∆ * ) = {(∆, ) : = 0 and ∆ ∆ * , or = 2, 4, . . . and ∆ + d − 2} . (2. 2) The variables to be solved for in (2.1) are the set of appearing pairs (∆, ) and the corresponding coefficients p ∆, . We are interested to know what is the maximal ∆ * such that (2.1) has a solution.
We would like to define a navigator function N (∆ * ) such that it is negative if a solution exists and is positive if it does not exist. To this end we will consider a modified problem of the form We just added an extra term in the l.h.s. with a fixed function M (u, v) and a new parameter λ. The function M (u, v) will be chosen so that the following crucial property holds: For any ∆ * , problem (2.3) has a solution with some λ = λ 0 (∆ * ) > 0. (2.4) Given this property, the navigator function will be defined as the minimal value of λ such that (2.3) has a solution: 7 N (∆ * ) = min λ such that (2.3) has a solution. (2.5) Property (2.4) then guarantees that the navigator is bounded from above, as we have N (∆ * ) λ 0 (∆ * ). We also see that the navigator is monotonically non-decreasing in the ∆ * direction, negative in the allowed region and positive outside. 8 This described construction does not formally guarantee other nice properties of the navigator that we wish to have (that N (∆ * ) is differentiable, strictly negative in the allowed region, has no local minima outside the allowed region where minimization can get stuck etc.) It also does not guarantee that the navigator is finite inside the allowed region (it may be −∞ there). Nevertheless, explicit navigator functions constructed below using this idea will have all these additional nice properties, by inspection.
We will now give two examples of functions M (u, v) that have the required property (2.4). 7 Although in this section we consider ∆ φ fixed, it is trivial to relax this and consider the navigator as a function of both ∆ φ and ∆ * , defined by the same Eq. (2.5). The zero set of N (∆ φ , ∆ * ) is then a curve which is the upper bound on ∆ * as a function of ∆ φ . We will not develop this idea further here but we will encounter analogous situations below in the multiple-correlator context. 8 Note that for any ∆ * the set of λ's for which (2.3) has a solution is a connected subset of the real axis. This follows from the fact that a convex linear combination of solutions is again a solution.

GFF-navigator
We know that for any ∆ φ , Eq. (2.1) has a Generalized Free Field (GFF) solution with the spectrum ∆ = 2∆ φ + 2n + , n 0, = 0, 2, 4, . . ., corresponding to operators of schematic form φ∂ n φ. The GFF-navigator is obtained by taking M (u, v) to be the first term in this solution: The GFF solution to crossing provides a solution to (2.3) with λ = 1 as long as all GFF operators besides φ 2 belong to S(∆ * ), which will be the case for ∆ * 2∆ φ + 2. Hence N (∆ * ) 1 for any ∆ * in this range.
Note that having a finite navigator in the range ∆ * 2∆ φ + 2 is sufficient for the problem at hand, since the boundary of the allowed region for (2.1) is known to satisfy this condition. Alternatively, higher GFF operators which do not satisfy gap assumptions may be added to the r.h.s. of Eq. (2.6). See App. A for this tweak of the GFF-navigator, important for bootstrap problems with additional gaps in the spectrum.

Σ-navigator
Another possibility, called the Σ-navigator, results from choosing: where (∆ i , i ) are any n spectrum points in S(∆ * ), c i > 0 some fixed positive coefficients, and n is a sufficiently large number. Since the coefficients c i are, apart from being positive, essentially arbitrary, there is a lot of freedom in choosing the Σ-navigator. Consider Eq. (2.3) with this M (u, v). In practice, in the numerical conformal bootstrap we analyze this equation in Taylor expansion around some point, i.e. we replace functions of u, v by vectors of Taylor coefficients of some finite length n 0 . Denoting vectors by boldface symbols, we have We claim that this equation will generically have a solution with some positive λ as long as the number of terms n in (2.7) is n n 0 . Indeed, generically the vectors F ∆ i , i are not expected to be linearly independent. Thus the equation will have a solution as longs as x i are allowed to have either sign. We rewrite this solution as For sufficiently large positive λ = λ 0 all the coefficients x i +λ 0 c i 0 so this is a solution to (2.8), proving the above claim. Hence, by the general arguments, the navigator is bounded from above by λ 0 .
In the described construction the number of terms n in (2.7) may have to be increased with the number of conformal block derivatives used in the numerical analysis. Alternatively, we may replace the sum in (2.7) by an integral with a positive continuous measure in some interval of ∆'s. Then the same navigator may be used independently of the number of derivatives.

Dual picture
In the dual approach to the numerical conformal bootstrap, the problem of computing the navigator (2.5) is formulated as follows: Our construction guarantees that the choices (2.6) or (2.7) lead to this problem having a solution bounded from above for any ∆ * . From this dual formulation we can see that the Σ-navigator is guaranteed to be finite also in the allowed region (i.e. it cannot be −∞ there). That's because for any ∆ * there is always some functional which satisfies the positivity condition in (2.11). Rescaling this functional we may make it also satisfy the normalization condition. This provides a finite lower bound for the Σ-navigator. For the GFF-navigator this argument clearly fails if ∆ * 2∆ φ . In this case there is no functional α satisfying both the normalization and the positivity conditions. Thus the GFF-navigator equals −∞ for ∆ * 2∆ φ . 9 This is not so problematic in practice, since this range is anyway deep inside the allowed region for the single-correlator problem. In principle the GFF-navigator could become −∞ even for ∆ * somewhat above 2∆ φ , but we have not seen this happen.
It is instructive to compare the above dual formulation with how one computes the maximal allowed value p max ∆ 0 , 0 of the squared OPE coefficient for an operator (∆ 0 , 0 ) present in the spectrum [13,9]: Comparing (2.12) with (2.11), one may wonder if one could perhaps define a navigator simpler than in our proposals, namely as for some appropriate choice of (∆ 0 , 0 ) in S(∆ * ). E.g. what if one tries 0 = 0 and ∆ 0 a little above the boundary of the allowed region? It turns out however that such simpleminded choices of functional normalization are inadequate. Namely, they give a finite navigator only in a rather small neighborhood of the boundary of the allowed region, which moreover gets smaller and smaller as one increases the number of derivatives used in the conformal bootstrap computation. 10 If one already knows quite well where the boundary is (e.g. via bisection), then using this navigator one can quickly determine it even more precisely. But if one starts far away from the boundary, this navigator would not help. Our Σ-navigator proposal shows that to get a robustly bounded navigator one needs to modify this idea by normalizing not on a single conformal block in the allowed region as in (2.12) but on a positive linear combination of many blocks as in (2.7). Analogously, one could have hoped to get a bounded navigator by normalizing the functional to −1 on a single conformal block in the region outside S(∆ * ). But again, one finds that choosing 0 = 0 and ∆ 0 a little below the boundary of the allowed region gives a navigator which is finite only in a small neighborhood of the boundary of the allowed region. Instead, our GFF-navigator proposal shows that if ∆ 0 is lowered all the way to 2∆ φ , which is quite a bit lower than the boundary of the allowed region, then the navigator becomes robustly bounded from above.

Multiple-correlator problems
We will now discuss how the navigator function construction generalizes to bootstrap problems involving several correlation functions. The main idea will be the same: we just need to add a new term so that crossing can always be obeyed, and minimize its coefficient.
We will consider the example of three 4pt functions σσσσ , σσ and where σ and are an odd and even scalars in a Z 2 -invariant CFT (such as the critical 3d Ising model). This system of correlators leads to 5 independent crossing relations [14]: 14) 10 Ref. [10] considered an early version of navigator function corresponding to normalizing one particular component of the functional to 1. This navigator prototype suffered from the same problem of being finite only in a small region. We are grateful to Tom Hartman and Amir Tajdini for enlightening communications concerning their findings, which sparked our search for a robust navigator function.
where V −,∆, is a 5-vector of functions while V +,∆, is a 5-vector of 2 × 2 symmetric matrices of functions of u, v: (2.16) See [14] for the expressions of the functions F ij,kl ±,∆, (u, v). The first sum in (2.14) runs over the Z 2 -even operators O + in the OPEs σ × σ and × (whose spin is necessarily even), while the second sum in (2.14) is over all Z 2 -odd operators O − in the OPE σ × (which can have any spin).
As usual, we will treat separately the unit operator contribution Furthermore, we will group the contributions of and σ using the relation λ σσ = λ σ σ . We will work in d = 3 and assume that all other scalars apart from and σ are irrelevant, so all remaining O ± will satisfy the spectrum restrictions: Then we can write (2.14) as If the point (∆ σ , ∆ ) is allowed, this equation must have a solution with P ∆ ,0 , P ∆, 0, p ∆, 0. As discovered in [14], 11 this condition gives rise to an allowed region in the (∆ σ , ∆ ) plane consisting of a small island containing the 3d Ising CFT and a larger detached "continent." We will first discuss how this can be reproduced using a twoparameter navigator N (∆ σ , ∆ ). See Section 2.2.1 below for how to include the third parameter θ parametrizing the ratio of the OPE coefficients λ σσ /λ . Analogously to (2.3), we consider the modification of (2.20) adding to the l.h.s. an extra term λ M where λ ∈ R and M is a particular 5-vector of functions of u, v: In general M will also have some dependence on ∆ σ and ∆ (just like all the other vectors in the equation). We will be looking for solutions of (2.21) with P ∆ ,0 , P ∆, 0 and p ∆, 0. Analogously to (2.4) and (2.5), the navigator N (∆ σ , ∆ ) is defined as the minimal λ such that a solution exists: N (∆ σ , ∆ ) = min λ such that (2.21) has a solution, (2.22) while M has to be chosen such that there is always some solution for a sufficiently large λ. This then provides an upper bound for the navigator and in particular guarantees that N < +∞. The GFF-navigator idea from Section 2.1.1 generalizes to the present multiplecorrelator setup. Indeed, we always have a GFF solution to crossing in which σ and are independent GFFs. The vector M is constructed from the contributions of (unitnormalized) operators 1 √ 2 : σ 2 : ∈ σ × σ, 1 With this M , Eq. (2.21) has a solution with λ = 1, P ∆ ,0 = 0 and P ∆, and p ∆, coming from the rest of the GFF spectrum in the σ × σ, × , σ × OPE. This guarantees that N GFF (∆ σ , ∆ ) 1. 12 To describe Σ-navigators we choose two finite sets R ± ⊂ S ± of (∆, ) pairs, and the linear equation where the variables X ∆, and x ∆, don't have to satisfy any positivity requirement. As in Section 2.1.2, the boldface symbols mean that we have switched to working at some finite order in Taylor expansion. Taking into account the structure of V 0,0 , V ±,∆, , V +,∆, , and the fact that the functions F ij,kl ±,∆, (u, v) are generically linearly independent (as follows from their expressions in [14]), Eq. (2.24) has a solution as long as R ± include sufficiently many points. 13 We won't need to know anything about the solution apart from the fact that it exists.
So let us pick any two such sets R ± with sufficiently many points, and define with some strictly positive fixed coefficients C ∆, 0, c ∆, > 0. For any such M Σ , Eq. (2.21) has a solution with some positive λ, by the same argument as in Section 2.1.2. Hence the corresponding Σ-navigator defined via (2.22) will be bounded from above.
As a final comment, we would like to recall another problem with the feasibilitymode searches which is resolved by our navigators. Feasibility-mode SDPB runs may not converge due to precision issues for points that can already be excluded using the bootstrap of crossing equations involving only a subset of the correlators [16]. E.g. this sometimes happens for points outside the 3d Ising island which are excluded by a singlecorrelator constraint. The navigators presented in this section converge in all the cases we tested, including the exact Ising setup that does exhibit this problem when run in feasibility-mode. Thus, navigators also provide a more robust method of checking the feasibility of any point.

Including the angles
As shown in [15], the allowed region in the 3-correlator bootstrap can be further reduced by treating the P ∆ ,0 term in (2.20) differently from the other P ∆, . This is possible since we are assuming is non-degenerate. Writing λ σσ = λ cos θ, λ = λ sin θ, p = λ 2 0, we can then specialize Eq. (2.20) as The original numerical implementation of this setup [15] involved scanning over the angle θ in addition to ∆ σ and ∆ , which was computationally laborious. Significant progress in reducing the computational cost has been recently achieved via the cutting surface algorithm [8].
In this paper we will show how this setup can be analyzed even more efficiently using the navigator function. The construction is almost the same as above. We simply add to the l.h.s. of (2.20) the term λ M and define the navigator N (∆ σ , ∆ , θ) as the minimal value of λ for which the so modified equation has a solution with p 0, P ∆, 0, p ∆, 0. We can choose M GFF as in (2.23), or M Σ as in (2.25), with R ± ⊂ S ± . The numerical results will be shown below.

Dual picture
The primal definition of the navigator function given above was convenient for clarifying the condition under which the navigator is bounded from above. For the actual numerical computation, we translate the primal definition to an equivalent dual formulation. As an example, for the 2-parameter navigator N (∆ σ , ∆ ), Eq. (2.22), the dual definition takes the form: For the 3-parameter navigator N (∆ σ , ∆ , θ) from Section 2.2.1 we have to simply replace condition (2.29) with (see (2.27)) We recall that the above dual problems can be then transformed into a polynomial matrix problem using rational approximations of conformal blocks expanded up to some finite derivative order around the z =z = 1/2 point. This polynomial matrix problem is then transformed into a semidefinite programming problem, which can be solved by SDPB [5,6].
In App. B we describe an alternative construction of the navigator function, which turns the feasibility problem into an optimization problem not at the level of crossing equations, but after the problem has already been dualized and translated into an SDP. We have not used that construction in this work, but it may turn out useful in future applications.

Visualizing the GFF-navigator
In the previous section we provided a formal definition of navigator functions. Their actual numerical evaluation can be performed using SDPB. Since navigator evaluation involves maximization, it will be comparable in cost to an OPE coefficient maximization, and more expensive than say testing feasibility of a point. Of course, we hope that this extra cost will be offset due to additional information provided by the navigator. And indeed, in subsequent sections we will see that complicated bootstrap tasks can be achieved with relatively few navigator evaluations.
Before we go to those applications, in this section we will explicitly visualize the various navigator functions of Section 2. We will do this to get some intuition about their "shape," and to check that they are sufficiently well behaved to allow application of minimization algorithms. Visualization will be done by performing fine scans in all variables. We emphasize again that in realistic applications we will not need to perform such expensive visualization scans.
We will focus on the 2-and 3-parameter GFF-navigators N (∆ σ , ∆ ) and N (∆ σ , ∆ , θ) from Sections 2.2 and 2.2.1. Numerical evaluation is done using the dual formulations given in Section 2.2.2, where we need to put M = M GFF from Eq. (2.23). We will not show plots for the Σ-navigators, although we have checked that they behave similarly to the GFF-navigators. 2-parameter case. We start with Fig. 1 showing N (∆ σ , ∆ ) in an extended region around the 3d Ising island at the derivative order Λ = 11. We can see from it that the region of negative navigator value matches in size the Λ = 11 allowed region of [14], Figs. 3 and 4. 14 On this scale the navigator is observed to be smooth (see however below) and approaching its predicted asymptotic value N max = 1 far away from allowed regions. There is clearly a valley coming from the top right of Fig. 1(left), narrowing to a tight gorge as it approaches its minimum inside the island. The surface has only one local minimum located in the plotted region and, as expected, it is inside the island. This feature will be essential when we discuss navigator minimization strategies in Section 5. Indeed, local minima in the disallowed region would have required more computationally expensive optimization methods than the BFGS algorithm discussed there.
In addition to the island, the allowed region found in [14] also included a detached "continent" at larger values of ∆ σ , beyond the range of Fig. 1. This continent is of course also found to be a region of negative navigator. Our navigator minimization strategies will use a bounding box, see Section 5.2, to make sure that we sail to the island and not to the continent. 3-parameter case. To get an idea of the shape of N (∆ σ , ∆ , θ), we will show two-dimensional slices for fixed values of one of the 3 parameters. Thus, in Fig. 2 we fix θ = 0.96926 (the central value from [15]), and let (∆ σ , ∆ ) vary in a region close to the navigator minimum. The surface shape is similar to the two-parameter navigator surface in Fig. 1. 15 Furthermore, in Fig. 3 we show 2d slices of the 3-parameter navigator arising for a fixed ∆ σ and ∆ . Although the precise shapes here are somewhat different, all three 2d slice surfaces are found to be smooth at this scale and free of local minima in the disallowed region (i.e. where the navigator is positive). This is a good sign that optimization algorithms should be able to quickly converge towards the Ising island given a reasonably precise initial guess.
Variation with Λ. Here we will explore how navigator shape changes with the derivative order Λ. By design, the navigator function monotonically increases pointwise with Λ, i.e. N Λ 2 (x) N Λ 1 (x) for Λ 2 > Λ 1 . This generalizes the fact that the allowed region shrinks with Λ. It is interesting to know how this increase happens. E.g. does the navigator surface move up with Λ uniformly or not? To answer this question, we show in Fig. 4 the 2d slice of the 3-parameter navigator at fixed θ = 0.96926 with Λ = 19, comparing it to Λ = 11 from Fig. 2. We see that the navigator surface has indeed moved up, but in non-uniform fashion. Most notably, the surface along one of the nearly flat "valley" directions gets lifted up much more than near the minimum. As a result, the minimum became more pronounced, which is a good sign.

Derivative of the navigator
The visualizations show navigator functions that are seemingly smooth and free of local minima. Both these properties would be very helpful for the numerical minimization algorithms, but they did not automatically follow from the definition of the navigator functions and we cannot guarantee that they hold in other setups. In fact, in the course of our investigations we found that even the navigator function under consideration is not entirely smooth: more precisely, we believe that it is not everywhere C 2 .
Our evidence is provided in figure 5. In this figure we consider a GFF navigator function with Λ = 11 for ∆ σ = 0.51831848513294, as a function of ∆ . (The chosen values of ∆ σ and ∆ are in the vicinity of the minimum that we found using the techniques described below. Notice that the navigator is negative along the entirety of the cross-section in figure 5 and so we are inside the Ising island. We also imposed the OPE relation λ σσ = λ σ σ but left the ratio λ /λ σ σ unspecified.) We plot both the navigator function itself as well as its first derivative in the ∆ direction. The kink in the latter plot strongly suggests that there is a discontinuity in the second derivative of the navigator. Indeed, the straight lines on either side of the kink allow us to reliably estimate the second derivative with finite differences: we find the value to be 767.762901557722(1) on the left and 219229.421457(1) on the right. Furthermore, using the two points closest to the kink we can estimate that the third derivative would have to be at least 10 23 if the navigator function were smooth, which seems highly unlikely.
Although we have only shown a single cross section plot, it is likely that the nonsmoothness persists along a line (segment) in the (∆ σ , ∆ ) plane. It would be interesting to understand its origin and whether there is a connection with the physics of the problem. Some preliminary investigations indicate that the discontinuity might be due to rearrangements of the extremal spectrum, but a detailed investigation is beyond the scope of this work.
Fortunately we will see below that the jump in the second derivative does not appear to inhibit the functioning of our minimization algorithm. We will comment more on this in the section 5.3.

Gradient at primal-dual optimality
In order to find points x where N (x) < 0 we will use a numerical minimization algorithm. The convergence rate of such algorithms is significantly improved if we also provide it with derivative information. In this section we therefore outline a procedure to compute the gradient ∇N (x).
Naively, one might think that gradient evaluation would involve computational overhead. For example, evaluating it via finite differences would require k additional SDPB runs where k is the number of variables on which the navigator depends. However this naive expectation is wrong: the main result of this section will be that ∇N (x) can be evaluated at negligible computational cost if we have already evaluated the function N (x) itself. The underlying reason is that the evaluation of N (x) is an extremization problem, and at extremality the first-order variation can be computed using only the original, unperturbed solution. This remains true even for constrained minimization problems, as is the case for us, when solved via primal-dual algorithms such as in SDPB, because primal and dual variables play the role of each other's Lagrange multipliers. To explain this in more detail we first have to introduce the semidefinite programming problem that underlies the computation of N (x).

Semidefinite programming reminder
Now we will explain how to compute the gradient of the objective in the above setup. As mentioned above, the evaluation of N (x) is computationally analogous to an OPE extremization problem that is often encountered in numerical bootstrap studies. Let us recall that, using a rational approximation for conformal blocks [17], these extremization problems become semidefinite programs with a particular structure of the constraint matrices. We will use the notation of [5], using which the problem can be written as: with S K the space of symmetric matrices of size K. Note that c ∈ R P is a vector, B ∈ (R n ) P a rectangular matrix, and the A * = (A 1 , . . . , A P ) ∈ (S K ) P is a vector of matrices. 16 In the language of convex optimization the program (4.1) is called a dual program (D), and the corresponding primal program P is given by: 17 We need a few more definitions. A vector x is said to be primal feasible if all the conditions in (4.2) are obeyed, even if optimality is not necessarily achieved. In the same vein a pair (y, Y ) can be dual feasible if it obeys all the conditions in (4.1). The duality gap is defined as the difference between the objectives: If x is primal feasible and (y, Y ) is dual feasible, then the duality gap is nonnegative: by the positive semidefiniteness of X and Y . So for any primal feasible point x the value of c T x provides an upper bound for the dual optimum, and similarly for any dual feasible point (y, Y ) the value of b T y provides a lower bound for the primal optimum. Now suppose one finds primal and dual feasible points with D(x, y) = 0. Then clearly both the primal and dual problem have been solved and brought to extremality, because neither objective has any room left to improve. It is a non-trivial fact of life that this condition is not only sufficient but also necessary for optimality in a generic semidefinite program (see [5] and references therein for details). In other words, rather than solving the primal or dual extremization problem, we can equivalently solve and then the optimal value of (4.1) and (4.2) is given by b T y = c T x. Notice that the fourth equation in (4.5) states that XY = 0 as a matrix equation. We call this the complementarity condition, and it follows from the vanishing duality gap, i.e. Tr(XY ) = 0, together with X, Y 0.

SDP gradient formula
Suppose we have found a primal-dual optimal point (x, y, X, Y ) such that the equations (4.5) are solved. To compute the gradient of the objective we change the parameters in the problem a little bit, and ask how the objective will change. So we need to investigate the corresponding linearized problem. The change in the solution must obey the linearized version of the optimality equations (4.5): (4.8) Our goal will be to compute the change in the dual objective, which is given by: In fact, since the duality gap remains zero we find d(c T x) = d(b T y) and one could equally well have computed the change in the primal objective. We start by showing a useful auxiliary result. The dX Y + X dY = 0 in (4.8) implies of course that Tr(dX Y ) + Tr(X dY ) = 0. We claim that a stronger result is true, namely that the two terms vanish independently: The proof is as follows. If XY = 0 and X, Y 0 then X and Y must have some zero eigenvalues. We can choose a basis where X is an upper block matrix, with X 11 0. Then any symmetric Y obeying XY = 0 must look like Now it becomes clear that the condition dX Y + X dY = 0 implies that X 11 dY 11 = 0 and dX 22 Y 22 = 0, which in turn implies (4.10). Let us return to the change in the dual objective as given in equation (4.9). Using the linearized optimality equations it can be written as: At this point we recall that x T A * = X. Moreover we have just shown Tr(XdY ) = 0. So the term proportional to dY in (4.15) vanishes, and we obtain: This "SDP gradient formula" constitutes one of the main points of our paper. It shows that the variation of the objective function of semidefinite programs (4.1) and (4.2) can be computed just from the variation of the data (db, dc, dB, dA * ) provided that we know the primal-dual solution (x, y, X, Y ). A remarkable fact is that we have eliminated all the dependence on (dx, dy, dX, dY ) from this formula.
In this work we will apply Eq. (4.16) to the navigator function. Once the navigator has been evaluated for some parameter values, Eq. (4.16) computes the gradient at the same point with negligible extra computational cost (see Section 4.2.1 below for how we organized the computation in practice). It's worth pointing out that this observation holds also for more familiar conformal bootstrap problems such as the OPE coefficient maximization. Such problems have been analyzed for years using primal-dual methods, but the existence of the "SDP gradient formula" has never been suspected by the people in the bootstrap community.
There is one important caveat to the preceding derivation. Although the solution (dx, dy, dX, dY ) to the linearized optimality equations does not appear in equation (4.16), we did need to assume that it existed in the intermediate steps. On the other hand, it is not guaranteed that the equations in (4.8) always have a solution. Fortunately this question has been analyzed in the semidefinite programming literature: for example, the paper [18] proves that the linearized equations for the semidefinite programs considered here will have a solution if X + Y 0, which in our notation is equivalent to the requirement that Y 22 0 rather than just Y 22 0. The paper [19] shows that this is in fact a generic property of the optimal matrices in semidefinite programs. We therefore expect the navigator functions to be generically C 1 , which is also confirmed experimentally by the smooth plots shown in the previous section.

Practical details for navigator gradient evaluation
As was shown in the previous section the change of the objective under a small perturbation of a bootstrap problem can be computed using only the solution to the initial problem (x, y, X, Y ) and the differences between the data (db, dc, dB, dA * ) defining the SDP. In this short technical section we describe the precise workflow using available codes. In order to be able to compute the gradient, one should first run SDPB on the original problem specified by (b, c, B, A * ) using the option --writeSolution="x,y,X,Y" to save the full solution to a file. The values (db, dc, dB, dA * ) can either be obtained by taking the difference between the perturbed and unperturbed bootstrap problem on the level of the polynomial matrix problem (PMP) and converting that to an SDP using pvm2sdp or by first converting both PMP's to SDP's and taking the difference between the resulting (b, c, B, A) and (b , c , B , A ). For the computations in this paper we did the latter. A dedicated tool called approx objective that takes one optimal checkpoint containing (x, y, X, Y ), one unperturbed SDP-file and one perturbed SDP-file 18 as input and outputs the corresponding change in the objective has been packaged with SDPB as of version 2.5.
When converting the PMP to an SDP, we can choose to keep the bilinear basis, sample points, and sample scalings (see [5]) the same for both the perturbed and the unperturbed SDP. With such choices, we automatically have dA * = 0 and the gradient formula simplifies to d(b T y) = db T y + dc T x − x T dB y.

Lagrangian perspective
In this section we give an alternative derivation of the SDP gradient formula. This derivation may look like a trick, but it provides an interesting perspective on why we were able to eliminate the variation (dx, dy, dX, dY ) from the change d(b T y) in the dual objective.
Consider the following Lagrange function: with µ > 0 a parameter, and it is understood that X 0. As is readily verified, the stationarity equations of this Lagrangian with respect to x, y and Y yield exactly the primal and dual feasibility conditions, i.e. the first three conditions in (4.5). Demanding stationarity with respect to X yields: with I the identity matrix. We can think of the last term −µ log det X in (4.17) as a barrier function that guarantees that X, Y 0. In the limit µ ↓ 0 the barrier disappears and we recover the original complementarity condition. 19 As is well known (e.g. [20]), the barrier function − log det X is convex.
We denote by (x(µ), y(µ), X(µ), Y (µ)) the stationary point of the Lagrange function, i.e. the solution of all the feasibility conditions and of the deformed complementarity condition XY = µI. Apart from degenerate situations, this solution exists; it is also unique. 20 The value of the Lagrange function at this solution is given by: since the constraints multiplying y and Y are obeyed by assumption. Furthermore, since the Lagrange function is stationary with respect to (x, y, X, Y ), its variation with respect to the parameters (b, c, B, A * ) is immediate: The modified complementarity condition (4.18) is also at the heart of primal-dual interior point algorithms as used in SDPB. Keeping µ finite and therefore X, Y strictly positive is useful to avoid getting stuck at the boundary where X and Y are singular. In the course of the algorithm the value of µ is then gradually reduced to zero in order to obtain a solution that obeys the original complementarity condition. See [5] for details. 20 Let M be the convex set of all On this set we consider the convex and smooth function t(x) = c T x − µ log det X. Generically the sublevel sets of this function are bounded. Indeed, any unbounded direction inside M can be parametrized as x 0 + λx with x 0 ∈ M , andx T B = 0 and with λ → ∞. If the original primal minimization problem is bounded, we have c Tx 0 for all directions. Generically we have a stronger condition c Tx > 0, in which case we eventually exit all sublevel sets for any such direction. Therefore t(x) must have a minimum inside M , and by convexity it is unique. At this point ∇t(x) = c − µ Tr A * X −1 is orthogonal to M . But the directions orthogonal to M are spanned by the gradient of the constraints, so by the columns of the matrix B. There must then exists some coefficients y such that ∇t(x) = By, which solve the last remaining equation. Now we can ask what happens if we take µ very small. Since the original semidefinite program is assumed to have a solution, we expect (x(µ), y(µ), X(µ), Y (µ)) to smoothly approach that solution as we send µ ↓ 0. Clearly, the variation of the Lagrangian (4.20) at the stationary point will then approach the right-hand side of our previous result (4.16). What of the value of the Lagrangian (4.19) itself? We know that X becomes singular and so det(X) will diverge. However, for Y = µX −1 to remain finite the eigenvalues of X cannot vanish faster than linearly with µ. We conclude that −µ log det X = O(µ log µ), the additional term in equation (4.19) vanishes in the limit, and so lim µ↓0 L(µ) = c T x. Together with equation (4.20) this reproduces (4.16).
This derivation elucidates the absence of (dx, dy, dX, dY ) from the variation of the objective. To summarize, the point is to replace the original constrained problem with an unconstrained one, involving a barrier function times a regulator µ. The unconstrained variation involves only variation of the data, and not of the solution itself. The constrained variation is recovered in the µ ↓ 0 limit and also has this property.
In Appendix C we push this logic one step further and explain how it can be used to compute the second variation of the objective, which one may call the "SDP Hessian formula." Also there we provide numerical tests of the SDP gradient and Hessian formulas. Having the Hessian as opposed to just the gradient could further speed up the minimization algorithms to be described in the next section, allowing to use Newton rather than quasi-Newton methods, but exploring this is postponed to future work.

Navigator minimization
A central task in the numerical bootstrap is the search for a feasible point. This corresponds to finding a point where the navigator function is negative. In addition, we may also be interested in finding the minimum of the navigator function, since its location might be close to the true CFT (we will show shortly that this indeed seems to be the case).
Given an n-dimensional search space, the search for a local minimum of the navigator function N (x) is a standard optimization problem. As explained in Section 4, the gradient of the navigator function is cheap to compute. Quasi-Newton methods can make good use of this cheap gradient. Recall that Newton's method requires computing a gradient and a Hessian at each point in the search. By contrast, quasi-Newton methods approximate the Hessian using gradient information. 21 In this work, we use the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm ( [21], Sec. 6.1) which is a well documented and widely used quasi-Newton method.
The BFGS algorithm maintains an approximation to the Hessian, which it updates using gradient information at each step. This update enforces positive-definiteness of the Hessian. Thus, it can only provide a truthful representation of the Hessian if the objective function is convex. In the examples studied in this paper, we have found that the navigator function is convex close to its minimum. However, this is not true further away from the minimum (for example, the GFF navigator tends to its asymptotic value 1 in a concave manner far away from allowed regions). This can lead to failure of the BFGS algorithm or less than optimal convergence. Therefore it is helpful to compose the navigator function with a monotonic function so that it becomes convex in a larger region but maintains the same minima. For example, if the maximal value of the navigator N (x) is N max (e.g. N max = 1 by construction for the GFF navigator), we can instead minimize Note that f (x) < 0 if and only if N (x) < 0, so that the allowed region is unchanged after this transformation. It's also easy to show that f (x) is convex in a larger region than N (x). 22 Intuitively, the main idea is that where N (x) approaches its asymptotic limit and hence is not convex, f (x) instead grows and has a chance to be convex.
x, then f (x) grows as |x| a , which is convex for a > 1. 23 To see how this works in practice, consider the GFF-navigator plotted in Fig. 1, which is clearly not convex. Fig. 6 shows the result of applying to it transformation (5.1) with N max = 1. We can see that the fractional linear transformation indeed improves the convexity of the objective function fed to BFGS. The function in Fig. 6  still not globally convex, but it is locally convex, or close to it, in a much larger region than the original function in Fig. 1. We will see below that this transformation indeed results in more appropriate step lengths in the initial line searches and that BFGS has a higher rate of success of finding the Ising model minimum, even when starting in regions of relative flatness of the untransformed navigator.
In our studies we will use the standard implementation of the BFGS algorithm which can be found in the SciPy library [22], with some minor modifications. In Section 5.1, we review the BFGS method. We describe our modifications and their motivation in Section 5.2. We will see in Section 5.3 that the resulting algorithm gives good results when applied to the 3d Ising model case. Finally, we will comment in Section 5.4 on further possible improvements on our modified BFGS algorithm.

BFGS algorithm
Let f (x) be the objective function to be minimized. The BFGS algorithm attempts to minimize f (x) by taking successive steps . . , where step k is taken using the information from an approximated quadratic model of the function at x k . This approximate quadratic model is where B k is an approximation to the Hessian at x k . After BFGS takes the k th step x k − → x k+1 , it determines the approximate Hessian at x k+1 by updating the one at x k using only gradient information at x k and x k+1 . For the full updating formula, see (6.19) of [21], Sec. 6. The minimum of the quadratic model (5.2) is the so-called "Newton step" In Newton's method, at each iteration the Newton step would be taken, so that x k+1 = x k + p k . In BFGS, the Newton step is replaced by a line search in the direction of p k . An exact line search would correspond to In practice, one uses an inexact line search, which means that one looks for an "approximate" minimum of φ(α) at α > 0. It turns out that a rather rough approximation is sufficient for good performance of the algorithm. A typical termination criterion is the "strong Wolfe conditions:" The first condition enforces that the function decreases sufficiently. The parameter µ controlling this is usually chosen to be very small. We used the default value µ = 10 −4 implemented in SciPy. The second condition demands that the gradient decreases sufficiently. This is usually called the curvature condition, and it guarantees that the BFGS update to the Hessian maintains positive-definiteness, 24 which in turn implies that p k+1 will be a decrease direction, allowing the algorithm to proceed. The parameter η controlling the demanded decrease is usually chosen somewhat below 1. We used the default value η = 0.9, with satisfactory results. The SciPy BFGS algorithm used in this paper relies on the Moré-Thuente line search algorithm [24], a standard and robust algorithm for finding points obeying the strong Wolfe conditions. Once an "accepted" point, i.e. a point obeying these conditions, is found, the Hessian is updated, 25 and the BFGS algorithm proceeds with its next step. The algorithm terminates once the norm of the gradient gets smaller than some value g tol supplied by the user (we used g tol = 10 −5 ).
It's worth pointing out that in the line searches, the Newton step α = 1 is used as the initial guess. Once the Hessian has been well-approximated (as may happen towards the end of the minimization run), the first step α = 1 will usually be accepted, and a convergence similar to that of Newton's method is expected. On the contrary, α = 1 may not be a good guess at the beginning of the run unless we have an idea of the typical size of the region in which the minimum is expected to lie. This is provided via a bounding box in our modified BFGS algorithm described below.

Modified BFGS algorithm
The BFGS algorithm requires an initial guess for the Hessian at the first step B 0 . This guess is usually taken to be the identity, which does not take into account the different scalings of the different variables. This is often okay because the BFGS algorithm recovers scale information after a sufficient number of steps have been taken, i.e. once the Hessian approximation becomes accurate in all directions. However, we still found that if some idea of the scale of the problem is known, e.g. if we have a vague idea about the location of the allowed region, it is best to incorporate this information into the initial Hessian. By setting a well-scaled initial Hessian, an appropriate step length in the initial line searches can be achieved. (Recall that the initial line search step is always α = 1 in the direction of p k , and the length as well as the direction of this search clearly depends on B k via (5.3).) This will ensure that the BFGS algorithm explores the vicinity of the starting point rather than a much larger space-a crucial feature in cases where we are interested in one specific nearby local minimum. For example, when we want to study the 3d Ising model, we are not interested in studying the navigator in the big allowed "continent" found at large ∆ σ [14]. We found that the BFGS algorithm may end up in this much larger feasible region unless an appropriately scaled initial Hessian is supplied.
One trick to set an appropriately scaled initial Hessian (based on [21], p.142) is the following: Compute the gradient at the initial point, and set B 0 to Then, from (5.3), the initial Newton step α = 1 will result in probing the function at . Hence the parameters α i 0 have the meaning of the characteristic desired |∆x i | during the initial step of the first line search. Alternatively, one could use the procedure described in Appendix C to explicitly compute the initial Hessian. However we do not advise this, since the Hessian for a point far away from the minimum could very well not provide an accurate scale for the problem, nor is the Hessian far away from the minimum likely to provide a more accurate starting point for the approximation of the Hessian at the minimum than an appropriately scaled diagonal matrix. 26 Apart from specifying the initial Hessian, some minor modifications have to be made in order to apply the BFGS algorithm to conformal bootstrap problems. Firstly, the navigator function is naturally defined only in certain regions and not globally. Consider for example the case of the 3d Ising model. The navigator function is naturally defined only for ∆ σ and ∆ above the unitarity bound. Similarly, when demanding the existence of exactly one relevant parity odd and even singlet, we restrict the domain of N (∆ σ , ∆ ) to values ∆ σ or ∆ below 3. Additionally, as discussed above, one might only be interested in minima or negative values that are located in a certain region around the starting point.
Hence, we ask the user to provide a bounding box for the search space, past which we do not allow the search to move. This constraint is implemented by altering the line search such that if a step outside of the bounding box would be taken, the maximal step in the same direction within the boundaries is taken instead. If this point on the edge is accepted, i.e. obeys the strong Wolfe conditions, we check whether the new search direction points inside or outside of the bounding box. If the new search direction points outside of the bounding box, but the gradient descent direction lies inside, the search direction is taken to be the gradient descent direction for the next step. If neither the initial search direction nor the negative gradient lie inside the bounding box, the search is terminated. The user should then either try a different initial point or change the bounding box.
The search direction p k is updated to p k+1 if x k+1 is at the boundary then if p k+1 points back inside the bounding box then continue else if −∇f (x k+1 ) points back inside the bounding box then p k+1 = −∇f (x k+1 ) else return x k+1 and the termination message "Out of the bounding box" end end end Optional: if f (x k+1 ) < 0 then return x k+1 and the termination message "Found a negative point" end end return x k and the termination message "Minimum found: gradient is smaller than the tolerance" end The provided bounding box is also used to specify the desired step lengths in the initial Hessian of (5.7). We found satisfactory results by setting the desired step lengths in each direction to be 20% of the supplied bounding box.
It is fair to ask how the user will know which bounding box to specify. We assume that the user has some idea of the range of parameters they want to explore. Results obtained at lower derivative order Λ can also be used for guidance, as well as estimates of CFT data coming from other methods such as RG or Monte Carlo simulations.
The BFGS algorithm including these modifications is summarized as Algorithm 1.

Minimization results
To illustrate the effectiveness of our minimization algorithm, we apply it to the classic conformal bootstrap problem of finding an allowed point corresponding to the 3d Ising model using the system of correlators { σσσσ , σσ , }, which contains the lowest dimensional Z 2 -odd scalar σ and the lowest dimensional Z 2 -even scalar , under the assumption that those operators are the only relevant ones, as described in Section 2.2. We will see that the navigator function enables us to locate an allowed point with a relatively small number of SDPB calls. Finding an allowed point naively by checking feasibility for a dense grid of points covering the search space would take orders of magnitude more SDPB calls.
Of course, in a decade of feasibility searches many useful tricks have been found to speed them up. 28 Still, we foresee that navigator-function methods will offer even better performance. They should eventually allow computations in more complicated setups involving an even higher number of parameters to scan over, such as e.g. bootstrapping the full system of σ, , 4pt functions, which were not possible to treat so far via feasibility-based methods.

2-parameter searches
We start with the 2-parameter case which is easier to visualize. So we minimize N (∆ σ , ∆ ) of Section 2.2. We use Λ = 11 and the bounding box [0.510, 0.530] × [1.30, 1.50], i.e. the same range as in Fig. 1. Running our algorithm for 10 different starting points chosen at random within this bounding box, the number F C of function calls to reach a point of negative navigator value was 9 F C 31, while F C = 19.3 on average. All runs terminated at essentially the same point (with an error controlled by g tol ) x f = (∆ σ , ∆ ) = (0.5182861212(4), 1.41521640889(6)), N (x f ) = −0.00267253307546000(2), (5.8) 28 E.g. for Ising and O(N ) we can use the fact that they live close to the kink in a single correlator bound, for Ising we can use c-minimization [7], OPE scans can be replaced with the cutting surface algorithm [8], etc.
where the tiny error bars show the largest difference observed between different runs. We conclude that the minimum is unique and all the runs terminate near it. . This run took 29 function evaluations to reach the island, and 66 function evaluations to reach the minimum within the specified g tol (see Fig. 8).
Only the first 38 points are marked, the rest being too closely spaced to be distinguishable.
A representative run is shown in Fig. 7, where the numbered points correspond to the path taken by our modified BFGS algorithm. Convergence rate in this run is illustrated in Fig. 8(left) where we plot the navigator values N i returned by subsequent function calls, until the negative navigator region is reached. This plot can be correlated with the navigator shape in Fig. 1, which features an arrow-shaped valley around the Ising island (see Section 3). Thus, Fig. 8(left) shows a period of modest progress in the minimization of N (∆ σ , ∆ ), in some sense looking for the the valley. This is followed by a period of fast decrease once the valley is found (at around i = 25).
Another way to evaluate the convergence rate is shown in Fig. 8(right), where we plot for the same run the distance x i − x f between the point x i and the eventually found minimum x f . This measure of convergence is appropriate also for the region where N (x) < 0. This plots show a period of greatly accelerated convergence towards the end of the run. Indeed, we expect Newton-like, i.e. superlinear, 29  Finally, we observe that minimum (5.8) of the Λ = 11 navigator function gives a good prediction for the actual location of the Ising model, as compared to a generic point in the Ising island. Indeed, the distance between this minimum and the best prediction from [15] (3-parameter scan at Λ = 43) is only ∼ 10% of the size of the Λ = 11 island, see Fig. 9.
These runs are shown in Fig. 10. Eighteen of them successfully converged to the 29 Recall that superlinear means i+1 = o( i ) where i is the error after step i. The Newton method has quadratic convergence, i+1 = O( 2 i ), while for the BFGS only weaker theorems showing superlinear convergence are available [21]. One-dimensional bisection in this notation has linear convergence, i+1 α i with α < 1. 30 Once the navigator function is negative, we often observe a bit of a plateau, for example between iterations 30-60 in Fig. 8. At that point we are relatively close to the minimum, but the long sequences of blue dots indicate that the BFGS quadratic model of the navigator function is not yet accurate. It is quite possible that this is caused by the non-C 2 locus that we identified above, and it would certainly be interesting to investigate this further. Either way, the algorithm eventually recovers and then continues to converge rapidly to the minimum.  same minimum inside the Λ = 19 Ising island: x f = (∆ σ , ∆ , θ) = (0.5181536110(7), 1.412692879(8), 0.969334757(6)), N (x f ) = −0.0000208827730 (5). (5.9) A typical successful run is shown separately in Fig. 11. Figure 11: A typical BFGS run from Fig. 10 (only a part of the bounding box is shown). First, the search is seen to be looking for the "valley", and once it has found it, converges rapidly to the Ising island.
Two runs terminated at the boundary of the bounding box with both the subsequent BFGS search direction and the gradient pointing outside, according to the safe-guarding procedure (see Algorithm 1). This suggests that these points were being attracted by a minimum outside the bounding box. By inspection, these runs started close to the edge of the bounding box in regions where the navigator surface is non-convex even after applying transformation (5.1).
Limiting to the successful runs, it took on average 50.3 function calls to reach the negative navigator region. Of course, the Λ = 19 island is orders of magnitude smaller in all directions than the bounding box. This demonstrates our point that the navigator minimization method is capable of finding a small isolated island given even a rough estimate of its location. We will comment in Section 5.4 on an iterative way to speed up high-Λ calculations.
Using the run in Fig. 11 as an example, we show its rate of convergence in Fig. 12, following the same conventions as in Fig. 8. Comparing Figs. 11 and 12, it's easy to reconstruct what is going on. The initial line searches are spent finding the bottom of the valley. Once this is found, the algorithm quickly manages to follow the valley towards the negative navigator region. Similar plots for six more runs from Fig. 10 are collected in App. E.2.
It's worth pointing out that in both Fig. 8(right) and Fig. 12 we see two periods of accelerated convergence: one when the negative region is approached and another towards the end of the run. The slower rate of convergence in between might be due to the function exhibiting some local concavity, or due to a large change in the local Hessian. We have not investigated this in detail.

Other algorithms and possible improvements
We have shown in the previous section that navigator minimization using our modified BFGS algorithm offers a robust and efficient method for finding an allowed point. However, there are bound to be avenues for improvement. We will remark on some potential improvements in this section. We hope that the algorithm presented here sets a good benchmark to which future algorithms will be compared. In order to efficiently find an allowed point at high values of Λ, one could imagine an iterative procedure where the navigator minimum point at a lower derivative order is used as an initial point for a minimization run at a higher derivative order (perhaps reducing the bounding box, or inheriting the Hessian estimate from the lower Λ BFGS run). This is expected to perform well for two reasons. First, because the navigator minimum provides an excellent estimate of the position of the Ising model, see Section 5.3.1, and hopefully also of other CFTs. Second, because of the accelerated convergence properties of the BFGS algorithm after reaching the convex region around the minimum (see Figs. 8 and 12). Using the exact Hessian computed as explained in Appendix C may be also especially beneficial in the convex region around the minimum.
In the above, we did not make use of the fact that the minimum of the navigator occurs close to N (x) = 0, and that in some cases, one will only be interested in reaching any negative point rather than the minimum. This information could be e.g. incorporated in the initial guess for the Hessian, by scaling the identity matrix such that the initial step aims towards zero of N (x) in a first-order expansion around the initial point (instead of scaling it so that the initial step explores some percentage of the bounding box, as done here) As discussed before, we have also found that the navigator function is not globally convex. We have found that in our case, this problem can be mitigated by minimizing another related, more convex function instead, Eq. (5.1). Even in regions of nonconvexity of this transformed function, we have found that line searches provided robustness to the algorithm. Still, other bootstrap problems may require more care when dealing with non-convexity. In such cases, algorithms where the updating formula for the approximate Hessian does not enforce positive-definiteness could be advantageous, see [25].
We have opted in our algorithm to constrain the search space in a rudimentary way via the bounding box, and found this to be adequate for our needs. With that being said, there exist a myriad of other algorithms for constrained optimization that could offer more robustness with the way they deal with constraints. Here we mention two included in SciPy: L-BFGS-B, a bounded limited memory version of the BFGS method optimized for dealing with problems with search spaces with a large number of dimensions, and SLSQP, allowing general, as opposed to box, constraints. See [21] for more information on constrained optimization.
Similarly, there are many unexplored avenues for parallelization. One could imagine parallelizing the line search, or using an inherently parallel optimization algorithm, in the spirit of particle swarms [26]. Particle swarm algorithms that we have seen do not make use of gradient information. Since we have gradients for free (Section 4), it would be desirable to develop a similar algorithm taking advantage of the gradients.

An application: exploring the tip of an island
In order to connect the Ising island to physical observables it is important to know its extreme points. For example, the left-and rightmost point of the island provide a rigorous lower and upper bound on the critical exponent η = 2∆ σ + 2 − d. In previous applications such bounds were often found by simply mapping out the entire island, using a higher-dimensional analogue of a binary search based on a Delaunay triangulation, and then locating its extremal points. A more systematic triangulation algorithm, suitable for parallelization, was introduced in [12] and used to determine the instability of the O(3) fixed point.
In future bootstrap applications one might want to study more complicated systems of correlators and this inevitably means the introduction of new parameters. If we wish to locate the extremal point of an island in such a higher-dimensional space then any triangulation algorithm based on a sequence of feasible and infeasible points will scale extremely poorly. A constrained optimization algorithm based on a navigator function is much less sensitive to the dimensionality of the parameter space and will perform much better. We therefore expect that the use of a navigator function is essential for the high-precision determination of critical exponents in the future.
In the next section we present a simple algorithm inspired by these general ideas. We will then maximize ∆ σ in the Ising island as an illustration.

A constrained optimization algorithm
Suppose we want to locate an extremal point of the allowed region in the direction specified by a vector n. The problem is then: over all x such that N (x) 0 . (6.1) We will use optimality conditions where the latter equation sets to zero all components of the gradient orthogonal to n. We propose to work towards a solution of these equations in a manner inspired by the quasi-Newton method from Section 5. We will now explain the full algorithm (see Algorithm 2 below for a summary). As in Section 5 we will use a quadratic model around a point x k : 3) The function and the gradient at x k are assumed known, while B k can be either the exact Hessian at x k (computed as explained in Appendix C), or an approximation like the one obtained from the BFGS method. In the following we will assume that B k 0. Substituting the quadratic model in (6.2) we find a simple system involving one quadratic and many linear equations, which can be solved exactly, yielding two solutions. 31 These are real if x k in the allowed region, so that N (x k ) < 0, and by continuity also in some domain outside the feasible region. In this case the surface N (2) (x) = 0 is an ellipsoid, and the second condition in (6.2) picks out the extremal points of this ellipsoid along the n direction. Some distance away from the allowed region the ellipsoid shrinks to zero size and the solutions become complex-conjugate. We denote by x # the real solution which has the largest value of n T x, when the solutions are real. When the solutions are complex conjugate, we let x # denote their real part (and then x # turns out to simply correspond to the minimum of the model function).
Denote p k = x # − x k ; this is our search direction. The next point x k+1 is then found using a line search along p k starting from x k . We use the initial step length α = 1, however the rest of the line search algorithm is not the Moré-Thuente algorithm used in BFGS. This should not be surprising since we are now solving a different problem. Instead of minimizing N (x) we would now like to maximize n T x while moving along a trajectory remaining close to the boundary of the allowed region (but not exactly along the boundary). One could think that a safe choice would be to remain always inside the allowed region (a sort of interior point algorithm). We have found however that a much faster algorithm results if we allow the algorithm to choose points on both sides of the boundary. To make sure that the algorithm does not veer off too much away from the boundary, we impose with a parameter λ rel > 0. Clearly λ rel < 1 would be safer but might slow down the algorithm in the later stages. We found it advantageous to use λ rel somewhat above 1, e.g. λ rel = 2 works well. So (6.4) is our line search termination condition. In practice, this condition with λ rel = 2 is not very constraining and the initial step α = 1 is almost always accepted if we start with a good initial Hessian. (E.g. in the run shown in Section 6.2 this happened for 100% of the steps.) In the cases that the initial step α = 1 does not obey Eq. (6.4), we proceed as follows. We construct cubic polynomial approximation P (α), fitted to match the value and gradient at the initial point and the previous line search point. If x k is in the feasible region we choose the next α by solving P (α) = 0, and if not then by minimizing P (α). Iterating this, eventually we find an α such that x k+1 = x k + αp k satisfies (6.4).
Once we have accepted x k+1 , we construct a new quadratic model around this point. In particular, if the approximate Hessian is used, then B k is updated as in BFGS. However, the update is carried out only if the curvature condition is obeyed at x k+1 ; as explained in footnote 24 this is sufficient to ensure that B k+1 0. If the curvature condition is not satisfied, then the Hessian is not updated.
We then repeat the process. The algorithm terminates if the conditions (6.2) are obeyed within a certain tolerance.

The tip of the Ising island
As an example, let's apply the above algorithm to find the maximal value of ∆ σ within the Ising allowed island N (∆ σ , ∆ ) 0 where N is the 2-parameter navigator for the Ising 3-correlator setup at derivative order Λ = 11. 32 The search was started from the navigator minimum reached via a BFGS run, and with the initial Hessian approximation B 0 inherited from BFGS, which is expected to be close to the true Hessian. The algorithm path is shown in Fig. 13. The algorithm took 17 steps to reach the tip of the island, i.e. the point with maximal ∆ σ . Termination condition max(|N (∆ σ , ∆ )|, |∂ ∆ N (∆ σ , ∆ )|) g tol was satisfied with g tol = 10 −27 .
For comparison, Fig. 13 also show the blue allowed region obtained from the Delaunay triangulation method. We finely sampled the zoomed-in region around the very tip Input: A navigator function N (x), a vector n indicating the maximizing direction, a precision goal g tol and a line search parameter λ rel . Output: The final point x f . begin Use Algorithm 1 to construct a feasible point x 0 and Hessian estimate B 0 and their gradients if N (x k ) < 0 then find α such that P (α) = 0 else find α such that P (α) is minimized end end of the island, with a total number of sampled points being around 480 33 . In contrast, our algorithm takes only 10 steps to locate the maximal ∆ σ point more accurately than the triangulation resolution. The line search never had to be activated, the initial try α = 1 having been accepted in 100% of the steps.
In Fig. 14 we show the convergence rate towards the minimum. These plots demonstrate superlinear convergence towards the end of the run, as should be expected from this type of algorithm.
We would like to warn the reader about a difference in spirit between our Algorithms 1 and 2. Algorithm 1 for navigator minimization is backed up by decades of experience in numerical optimization, and should be widely applicable without major modifications. On the other hand, Algorithm 2 is our own custom-made procedure. It served well the purpose to demonstrate the point that the navigator can be used to find extremal values of allowed parameters, but it has a somewhat tentative character and is expected to evolve more in the future.
For example, the Ising island is admittedly a simple model with a convex island and a single local maximum of ∆ σ . If the island does not have such a nice shape, Algorithm 2 can get stuck in a local optimum instead of the global optimum. In more realistic cases it is therefore important to have a rough idea of the shape of the island, and then perhaps an admixture of triangulation-based methods and navigator methods might be the best approach.

Conclusions and future directions
We have presented in this work a powerful alternative to the scanning-based approach employed so far in the numerical conformal bootstrap program. This came from the realization that there exist functions, for which we have coined the term "navigator functions," which measure how far a given point is from the boundary between allowed and disallowed regions and can thus be used to efficiently find an allowed point as well as the boundary of an allowed region. We have explicitly constructed two such navigator functions. It was shown that the computation of these navigator functions can be written as a semi-definite programming problem of the same form as an OPE maximization. Adding the generalized free field solution to the crossing equation has led us in Section 2.1.1 to the definition of the GFF navigator. The Σ-navigator was introduced in Section 2.1.2 as an another equally valid option. 34 With the help of such functions, we have shown it is possible to quickly locate allowed regions in parameter space by ways of minimization. We have presented in Algorithm 1 a modified BFGS algorithm which does so quite efficiently. To prove this, we set out to study the canonical bootstrap problem of the 3d Ising model. First we showed that the navigator is (C 1 ) smooth and has no local minima in the disallowed region. With both a two-dimensional search space at Λ = 11, and a three-dimensional search space at Λ = 19, we have shown that it took on average a few dozen SDPB calls to find the Ising island (19.3 for the former, 50.3 for the latter), starting only from very 34 While the GFF-navigator is naturally normalized, the Σ-navigator has its own set of advantages. It is actually easier to set up, since one does not have to work out the GFF OPE coefficients. In addition, there is not one but infinitely many Σ-navigators, corresponding to different choices of terms in the r.h.s. of (2.7) or (2.25), and this flexibility may prove useful in the future. At present, we see no definite reason to prefer one or the other navigator. For comparison, we performed some of the reported computations using both navigators (e.g. section 6.2), and they performed equally well.
conservative estimates of the parameters. This is competitive with previous methods for isolating islands and bounding CFT data. Moreover, these previous methods suffered from exponential scaling with the dimensionality of the search space. This constituted a major bottleneck for the kinds of problems that could be tackled: realistically only setups with a handful of free parameters could be considered. We expect that the scaling of the minimization-based navigator method with the number of parameters will greatly outperform scanning methods.
Crucially, efficient minimization of a navigator function, for example with the BFGS algorithm presented in this paper, requires the knowledge of derivatives of the navigator function. We have derived the "SDP gradient formula," Eq. (4.16), which gives the variation of the objective function of an SDP as only a function of the variation of the SDP input parameters around the point where the derivative is requested. This means that computing derivatives does not require additional SDPB runs, making one function and gradient evaluation in a BFGS run just about equivalent in cost to one OPE maximization.
We also tested the efficiency of the navigator method to search for extremal parameter values allowed by the bootstrap constraints. So, we presented in Section 6 a way to find optimal bounds on CFT data using a custom-made constrained-optimization routine. The algorithm was able to walk in and around the allowed region and converge in 17 steps to the maximal allowed ∆ σ , determining it to an accuracy of ∼ 10 −35 . 35 A similar triangulation-based search only achieves an accuracy of 10 −6 even after testing over 400 points, see Fig. 14. Again we expect that the increase in performance can only become greater as the dimensionality of the search space increases.
We feel that the applications shown in this paper demonstrate only a small part of the power the navigator method, and we are hopeful that the future will show it to be a great addition to the toolbox of all bootstrap enthusiasts.
We would like to conclude by mentioning here some of the ideas that we are going to start exploring immediately ourselves using this new tool. Indeed, these applications, out of reach of traditional bootstrap techniques, were among our chief motivations to start thinking hard about the navigator function.
One class of situations where navigator is going to be useful is when we know a solution to bootstrap constraints for some value of a parameter (such as space dimension d or the symmetry group rank N ) and we would like to perform a deformation in this parameter. We imagine doing this by considering a navigator function depending on the dimensions of several exchanged operators, and imposing sparsity of the exchanged spectrum. Among other things, this should allow a more robust determination of critical values of parameters when bootstrap solutions disappear, than the more traditional approach of looking for kinks and trying to see when those kinks get rounded off. One long-standing problem which could benefit from this approach is determining the upper critical dimension of the 3-state Potts model. Including exchanged operator dimensions among the arguments of the navigator function could also provide a useful (and more rigorous) alternative to estimating the spectrum using the extremal functional method [9,27,7].
The use of the navigator function to quickly find extremal allowed values (Section 6) will benefit all cutting-edge bootstrap computations. One problem on our to-do list is to bootstrap the system of correlators in O(3) symmetric CFTs involving lowest scalar primaries in vector (φ), scalar (s), rank-2 tensor (t), and rank-4 tensor (t 4 ) O (3) representations. This setup extends that of [12] by including t 4 as an external operator. The physics interest in doing so is that it will allow access to the OPE coefficient λ t 4 ,t 4 ,t 4 , and other data needed to study the RG flow leading from the O(3) fixed point to the cubic fixed point in conformal perturbation theory (see [12], Section 5). The parameter space for this problem is 13-dimensional (4 ∆'s and 9 OPE coefficients), out of reach of traditional approaches, but we expect that the navigator function will put it within reach.
As a final example, we expect that the navigator functions will allow an exploration of hybrid methods where the numerical bootstrap data is complemented with analytical data at high spins obtained from the light-cone bootstrap, as suggested in Section 9.1 of [28]. We imagine a navigator function depending on many parameters accurately parametrizing one or more Regge trajectories. In this context a navigator function will be very useful not only to localize an allowed point, but also because the minimum of the navigator offers a natural "most feasible point" that can be used to compare different parametrizations. Although this method is not entirely rigorous, it might lead to more precise estimates of the numerical bounds.
formance Cluster, partially supported by a grant from the Gordon and Betty Moore Foundation. This work also used the Extreme Science and Engineering Discovery Environment (XSEDE) Comet Cluster at the San Diego Supercomputing Center (SDSC) through allocation PHY190023, which is supported by National Science Foundation grant number ACI-1548562. The computations in this paper were partially run on the Symmetry cluster of Perimeter institute and on the Hopper cluster of theÉcole Polytechnique.

A Tweaks of the GFF-navigator
As mentioned in Section 2.1.1 and footnote 12, the GFF-navigator definition has to be tweaked in presence of additional GFF operators violating gap assumptions. These modifications will be discussed here. In addition we will explain how to deal with the case where the navigator function depends on the magnitude of a squared OPE coefficient.
A relevant example in the single-correlator setup of Section 2.1.1 is to assume a gap in the scalar spectrum above ∆ * . E.g. suppose that all further scalars above the one at ∆ * are required to be above ∆ gap . This corresponds to changing the constraint ∆ ∆ * for = 0 in (2.2) to "∆ = ∆ * or ∆ ∆ gap ." We can still define the navigator by the same Eq. (2.5). In this case we don't in general expect the navigator to be monotonic in the ∆ * direction. For large ∆ gap , definition (2.6) of M GFF will have to be modified, including all scalar GFF conformal block contributions below ∆ gap : where c n are explicitly known coefficients (c 0 = 2). These are contributions of GFF operators of schematic form φ n φ.
For the 3-correlator setup, let us discuss how the GFF-navigator definition (2.23) should be modified in the case of gap assumptions in the spectrum of 1 operators. As a concrete example, let us define the navigator N (∆ σ , ∆ , c T ) where c T is the 2pt function coefficient of the canonically normalized stress-tensor. The c T parametrizes the OPE coefficients of the corresponding unit-normalized ∆ = 3, = 2 primary O + as: where K d is a known d-dependent constant. To isolate the stress tensor, we need to impose a gap assumption on the higher-dimension = 2 O + operators. We will assume that all of them have ∆ ∆ gap where ∆ gap > 3 is some fixed parameter. E.g. let us choose ∆ gap = 5, which allows the 3d Ising CFT. 36 36 Recall that the 3d Ising CFT has ∆ T = 5.50915(44) [28].
To define the GFF navigator, we will proceed analogously to (A.1) and include in M GFF additional terms corresponding to all GFF primaries violating the gap assumptions. In the case at hand, we have to check the spin-2 GFF operators of the schematic form σ∂∂ n σ and ∂∂ n . For ∆ gap = 5 and ∆ σ , ∆ around the 3d Ising island, only the n = 0 operators of this form are below the gap. So we take where c(∆ φ ) is the (explicitly known) coefficient of the ∆ = 2∆ φ + 2, = 2 conformal block in the decomposition of the GFF 4pt function φφφφ . The first two lines in (A.4) are the analogue of (A.1). The last line is an additional small modification needed due to the presence of OPE coefficient parameter c T among navigator function variables. It is the negative of the stress tensor contribution in (A.3). Including this piece into M GFF is needed to guarantee that problem (A.3) has a solution with λ = 1 for any fixed value of c T . This in turn guarantees that the navigator N (∆ σ , ∆ , c T ) is bounded from above by 1 for any value of its arguments.

B Feasibility as optimization
In this appendix we discuss the problem of finding a navigator function from the more abstract semidefinite programming perspective. We will assume the reader is familiar with the semidefinite programming terminology of Section 4.1 in the main text. As we review there, a general numerical bootstrap problem of the opimization type can be formulated as the dual problem given in Eq. (4.1) on p.20. For a feasibility problem, on the other hand, the question is merely whether there exist any y and Y that obey the constraints. In that case the standard approach is to set b = 0 in (4.1) and run SDPB until one of two termination conditions are met: • If a dual feasible point (y, Y ) is found, terminate with 'success'; • If a primal feasible point x is found and c T x < 0, then terminate with 'failure'.
The last termination condition is explained by the duality gap: if b = 0 then D(x, y) = c T x, which can only be negative (for a primal feasible x) if no dual feasible point exists. 37 The above two termination conditions correspond to the binary oracle output discussed in Section 1: "success" means that the point is excluded (CFT does not exist), while "failure" means that the point is allowed (CFT may exist).
To pass from this to a navigator function, we need to reformulate the feasibility search as an optimization problem. The commonly adopted approach to do so is to use slack variables that relax the constraints. As discussed in the main text, in the context of the conformal bootstrap one can add an additional term to the crossing equations, in such a way that these equations can always be obeyed if the coefficient of this extra term is positive. The minimization of the coefficient of this term is then a potential navigator function: if it is positive we are in the 'success' region and if negative we are in the 'failure' region.
We will now describe an alternative navigator function construction, which does not rely on the physical intuition of the crossing symmetry equations. Instead, we will start with a general feasibility semidefinite program of the type (4.1) with b = 0, and transform it into an optimization SDP.
As a first attempt, consider replacing the condition Y 0 in (4.1) with a maximization problem: with I the identity matrix. With this transformation the 'success' and 'failure' cases mentioned above respectively correspond to ν > 0 and ν < 0 at optimality, and (in the conventions of the main text) we can take ν at optimality as a candidate navigator function.
Unfortunately the modification (B.1) is not guaranteed to give a finite navigator in the "success" region. E.g. suppose there exists a Y 0 such that Tr(A * Y ) = 0. In the 'success' region one can add this Y to any feasible solution Y with arbitrarily large coefficient. This would then imply that ν → +∞ at optimality. We therefore cannot exclude a divergence in this candidate navigator function unless we know that the program does not allow such Y .
To guarantee boundedness in the 'success' region, we apply the same idea, but on the primal side, that its, by modifying the primal problem (4.2). For simplicity, let us 37 With b = 0 the primal problem is completely homogeneous in the sense that the constraints are invariant under rescalings x → λx with non-negative λ. In particular, there is an obviously primal feasible point x = 0. Since this point teaches us nothing about dual feasibility, the inequality in the second termination condition has to be strict. Furthermore, if we were to ignore the above termination conditions and run the program to optimality then we would either find x → 0 (in the 'success' case) or x diverges such that c T x → −∞ (in the 'failure' case). We thank Petr Kravchuk for a discussion of these issues. first assume that there always exists an x such that meaning we only need to introduce a slack variable for the positive semidefiniteness condition. In that case the right problem to solve is: This is a standard primal semidefinite programming problem, and we can run it to optimality without special termination conditions. The value of ν at optimality is the navigator function. In the 'success' region it is guaranteed to be positive (and finite) and then it is likely to be as good a navigator function as the ones used in the main text.
The dual version of the program in (B.3) is: As usual, the introduction of free variables on one side yields additional constraints on the other side. In this case the trace condition on Y guarantees the boundedness of the problem, and the parameter ξ allows for the re-scaling of a feasible (y, Y ) such that this constraint can be met.
Let us also discuss boundedness (from below) in the 'failure' region of (B.3). We do not have a first-principles argument for boundedness everywhere: 38 for the same reasons as above, the navigator function of (B.3) diverges in the 'failure' region if there exists a x which obeys Fortunately, in conformal bootstrap applications this is unlikely. To see this, recall that the formulation (4.1) with c and b arises only after eliminating one component of y from a normalization condition n T y = 1 for some normalization vector n, which is typically the identity operator. Reinstating this normalization condition as a separate constraint to (4.1) one finds that unboundedness of the modification (B.4) can really only occur if there is a solution to the crossing symmetry equations (with positive coefficients) without an identity operator. Although this is known to be the case for problems in d = 2 and d = 1, it is an unlikely possibility in most numerical bootstrap problems and then (B.4) is also bounded in the 'failure' region. The corresponding navigator therefore obeys the same manifest properties as those used in the main text. Finally let us consider the case where the equality constraints in the primal problem cannot obviously be met. In that case not all is lost: one can simply replace them with with 1 = (1, 1, . . . 1) a constant vector, and proceed by minimizing ν + i λ i . As before, a positive value at optimality means that no feasible point exists and so we still have a good candidate for a navigator function in the 'success' region. The navigator functions introduced in this appendix are more general since they work for any feasibility problem of the type described in Eq. (4.1) with b = 0. On the other hand, for numerical conformal bootstrap applications they offer little upside compared to the GFF and Σ-navigators discussed in the main text. Furthermore they also suffer from a practical disadvantage. To see this, note that the GFF and Σ-navigators are readily implemented with the usual conformal bootstrap software: programs like sdp2input or pvm2sdp can be used to translate the problems into a format acceptable by SDPB, which e.g. involves setting up matrices B and A * , and SDPB then does the rest of the computation. Unfortunately this workflow does not quite work for the navigator function described in Eq. (B.4). The main problem is that SDPB is meant to solve problems where the matrices A p have rank one and the constraint Tr(Y ) = 1 is not of this form. 39

C Comments on variations of the objective
In section 4, we found a simple formula (4.16) for the linear-order variation in the objective function under changing the SDP. In this appendix, we give a formula for the quadratic-order variation as well, and explain how it can be computed easily using machinery already present in SDPB. We also present numerical checks of both the linear and quadratic variations, determining how their errors scale with the duality gap. 40 39 One can probably impose the trace constraint in an SDPB compatible way, by extending y with spurious variablesŷ. One then needs to set these equal to the diagonal components of Y in the sense thatŷ 1 = Y 11 ,ŷ 2 = Y 22 , etc. This can be done by including one additional equation for each diagonal value of Y by extending b, c, B and A. Finally, by extending these quantities by one more entry we can impose the trace constraint by demanding iŷ i = 1. Alternatively one can use this equation to eliminate one of these extra components instead. It is unclear whether such an altered semi-definite problem still corresponds to any polynomial matrix problem. 40 The quadratic variation of the objective could be used to compute the Hessian of the navigator function, enabling the use of Newton's method for finding allowed points and extremizing CFT data. We leave possible applications of the quadratic variation to future work.

C.1 A formula for the quadratic variation
Consider changing an SDP by (b, c, B, A) → (b, c, B, A) + (db, dc, dB, dA). For simplicity, we assume dA = 0. (In practice, we can ensure this by keeping constant the "bilinear basis" and "sample scalings" discussed in [5].) The linear-order change in the objective at optimality is where L is the Lagrange function (4.17).
As explained in section 4.3, dL is independent of (dx, dy, dX, dY ) because the variation of the Lagrange function with respect to (x, y, X, Y ) vanishes at optimality. The same reasoning implies that the quadratic variation in the objective should be linear in (dx, dy, dX, dY ). To compute it, we will work at finite µ. Afterwards, we consider the µ → 0 limit of the resulting expression and assess the size of finite-µ corrections.
For brevity, let us write s = (b, c, B) and z = (x, y, X, Y ). Given a change s → s+ds, the solution changes as z → z +dz +d 2 z +. . . , where dz and d 2 z are linear and quadratic in ds, respectively, and ". . . " represent higher order terms in ds. The quadratic change in the Lagrange function is Here, s and z are multidimensional and we suppress indices for brevity. The first term on the first line vanishes by the optimality equations ∂L ∂z = 0, and the last term vanishes because L is linear in s. The remaining two terms are proportional to each other. To see this, note that under changing s → s + ds, the shifted optimality equations become 0 = ∂L(s, z) ∂z Contracting (C.3) with dz and plugging this result into (C.2), we find The variations dx, dy can be computed from the linearized optimality equations (C.3), which are written in more detail in (4.8). After some rearrangement, we find where S pq = Tr(A p X −1 A q Y ) is the so-called Schur complement matrix. This is precisely the equation solved by SDPB in its main optimization algorithm, with a modified right-hand side. Consequently, it is straightforward to adapt SDPB to determine dx, dy and compute dL and d 2 L. We have implemented this computation in a program approx objective packaged with SDPB as of version 2.5. 41

C.2 Possible sources of error
We note two possible sources of error in the results for dL and d 2 L -one conceptual and one practical: The formulas for dL and d 2 L were derived assuming finite µ (so that the optimization problem is well-posed). Is the µ → 0 limit of these expressions well-behaved? How big are the finite-µ corrections?
As with the objective function itself, we expect errors in dL and d 2 L to be of order O(µ log µ), provided the SDP is generic. This expectation comes from thinking about L as a function to be optimized b T y + c T x − x T By + Tr (X − x T A * )Y , plus a barrier function −µ log det X that imposes that X is positive semidefinite. Near a smooth point on the boundary of the positive-semidefinite cone, the barrier function effectively moves the boundary of the cone by a smoothly-varying amount proportional to µ.
As we vary the parameters (b, c, B, A), the optimal solution with µ = 0 moves along the boundary of the positive semidefinite cone. Similarly, the optimal solution with finite µ moves along the "effective" boundary a distance µ away. As long as the boundary is smooth, derivatives of the finite-µ objective will differ from derivatives of the µ = 0 objective by O(µ log µ) (the size of the barrier function).
(E 2 ) Errors from XY = µI. One of the optimality equations (4.18) is XY = µI. Under normal operation, SDPB does not attempt to solve this equation with high precision. Instead, it performs repeated Newton steps toward solutions of XY = µ (i) I with values µ (i) that change with each iteration. This turns the equation XY = µI into a kind of moving target. Solutions computed by SDPB will generally have nonzero (but small) XY − µI.
It is not a-priori obvious how large errors resulting from nonzero XY − µI will be. (We show a numerical example in figure 16.) However, they can be mitigated with a simple strategy: After SDPB terminates with a primal-dual optimal solution, we can perform a few extra iterations toward a solution of XY = µI. In practice, this can be done by running SDPB from the most recent checkpoint with the options listed in table 1 (in addition to whatever other options were used in the optimization). Because the locus XY = µI is called the "central path," we call these extra iterations "centering iterations." option explanation --maxIterations=n Control the number of iterations. We take n = 10 below.

--stepLengthReduction=1
Take full Newton steps instead of decreasing the step size. --infeasibleCenteringParameter=1 Ensure that µ stays (nearly) constant instead of changing µ → βµ with each iteration. This option is only effective if SDPB has both a primal-and dualinfeasible internal state.

C.3 Numerical checks
To describe our numerical checks of the expressions for dL and d 2 L, we need some quick definitions. Given an SDP s, let f (s) be the optimal value of its objective. We also define g(s, ds) ≡ f (s) + dL + d 2 L = f (s) + ∂L(s, z) ∂s ds + 1 2 ∂L(s, z) ∂s∂z ds dz, (C. 6) where z is the optimum of s, and dz (which is linear in ds) is the solution to equation (C.3). Note that g is arbitrarily nonlinear in its first argument, but quadratic in its second argument -in fact, g(s 0 , s − s 0 ) provides a quadratic approximation to f (s) around a given s 0 : 42 Consider now a family of SDP's s(∆) depending smoothly on a parameter ∆. Consider a sequence of values ∆ 0 + δ∆ converging to ∆ 0 , and let us write s 0 = s(∆ 0 ). Equation (C.7) with s = s(∆ 0 + δ∆) implies that h(δ∆) ≡ f (s(∆ 0 + δ∆)) − g(s 0 , s(∆ 0 + δ∆) − s 0 ) ∼ O(δ∆ 3 ), (C. 8) where we used that s(∆) depends locally smoothly on ∆. We can use this to check our expressions for dL and d 2 L: we compute h(δ∆) for several values of δ∆ and check whether it decreases cubically in δ∆. . . , 25}, and derivative order Λ = 11. We see that the difference between the true objective and its quadratic approximation is cubic in δ∆. The optimizations for this plot were computed with a duality gap threshold of 10 −30 , and 10 centering iterations.
In figure (15), we plot the ratio h(δ∆)/δ∆ 3 for a one-parameter family of SDP's describing the GFF navigator function for correlators of σ and in the 3d Ising model. For small δ∆, the ratio h(δ∆)/δ∆ 3 approaches a constant. This is a strong check of our results for dL and d 2 L and our ability to compute them accurately: cubic dependence of h(δ∆) on δ∆ requires delicate cancellations between the true objectives of s(∆ 0 + δ∆) and s 0 , the linear correction dL, and the quadratic correction d 2 L. The SDPB computations in figure 15 were performed with duality gap threshold D = 10 −30 , with 10 centering iterations. Evidently these choices effectively remove both sources of error (E 1 ) and (E 2 ) in this example. 43 In figures 16 and 17, we show the effects of (E 1 ) and (E 2 ) on dL and d 2 L. Figure 16 was produced with no centering iterations, so it shows the effects of both (E 1 ) and (E 2 ). In that case, the relative error in dL scales approximately as µ 0.8 , and the relative error in d 2 L scales as µ 0.235 . These numbers presumably are not universal: they depend on the whole history of the optimization procedure in SDPB, and are not uniquely determined by the final solution. Figure 17 was produced with 10 centering iterations. In that case, the errors in dL and d 2 L both scale linearly with µ, and are much smaller overall. This is strong evidence that centering iterations effectively mitigate (E 2 ), and it also supports our estimate of the size of finite-µ effects.   Figure 17: Errors for dL and d 2 L as a function of the duality gap D with the same setup as figure 17, but where for each optimization we perform 10 centering iterations of SDPB. The errors now decrease linearly with D (which is proportional to µ). This is consistent with our naive estimate µ log µ in section C.2. (To detect the logarithm log µ, we would need more data and a more careful fit.)  Table 2: Parameters used to setup the SDPs, along with the SDPB parameters. The definition of these can be found in [5] (where order was 90 and keptPoleOrder was κ).

D Parameters for numerics
The computation of the navigator function can be translated to the form of a semidefinite program (SDP), to solve which we use the arbitrary precision solver SDPB [5,6]. We used simpleboot [29], PyCFTBoot [30], and sdpb-haskell 44 to setup the SDPs. The parameters used for the computations are presented in Table 2. We used the same conformal block normalization as [1]. For the Λ = 19 results in Section 5.3, we used the Python package PyCFTBoot [30] to setup the SDP, with parameters (k max , l max , n max , m max ) = (28,28,1,9). The parameters (n max , m max ) control the number of derivatives used in the (a, b) coordinates (see [30] for more details). This choice results in the same navigator value as taking (z,z) derivatives up to Λ = 19.
To numerically implement the BFGS Algorithm 1, we have used the BFGS algorithm minimize(method='BFGS') of Python's SciPy library, with the additional modifications of the rescaling of the initial Hessian and the implementation of the bounding box. All parameters used were the default ones, both for the Moré and Thuente line search SciPy implements and the actual BFGS algorithm.

E Further plots
Here we collect plots like Figs. 7 and 8 for six additional runs of our modified BFGS algorithm, for both the two parameter Λ = 11 case discussed in Sec. 5.3.1, and the 44 https://gitlab.com/davidsd/sdpb-haskell 45 The computations presented in Sections 5.3 and 6 were set up using a version of simpleboot where the definition of keptPoleOrder was slightly different. Here the poles were kept without modifying the residue to better approximate the contribution of discarded poles and thus the blocks were less accurate than those used in [5]. three parameter Λ = 19 case discussed in Sec. 5.3.2   Figure 19: This plot is analogous to Fig. 8(left). It shows navigator values N i at the i-th function call for the 6 runs from Fig. 18, and with the same color code for the dots.