A semidefinite program solver for the conformal bootstrap

We introduce SDPB: an open-source, parallelized, arbitrary-precision semidefinite program solver, designed for the conformal bootstrap. SDPB significantly outperforms less specialized solvers and should enable many new computations. As an example application, we compute a new rigorous high-precision bound on operator dimensions in the 3d Ising CFT, Δσ = 0.518151(6), Δϵ = 1.41264(6).

With recent analytical work [34][35][36][37][38][39][40], we now understand in principle how to formulate many interesting bootstrap calculations, like studies of four-point functions of scalars, fermions, conserved currents, stress-tensors, and mixed correlators combining all of these ingredients.Such studies may shed light on the classification of critical points in condensed matter systems, the conformal windows of 4d gauge theories, the landscape of AdS string vacua, and more.Each study inevitably culminates in an optimization problem that only can be solved numerically at present.SDPB is a custom optimizer with three purposes: 1.To enable the next generation of bootstrap studies involving conserved currents, stress tensors, and other complex ingredients.
2. To help improve predictions for CFTs already isolated with the bootstrap, like the 3d Ising CFT, 3d O(N) vector models, and others.
3. To demonstrate the potential for dramatic improvement in current numerical bootstrap techniques and to encourage researchers in numerical optimization and computer science to contribute ideas and expertise.
Custom optimizers have improved bootstrap calculations in the past.For example in [5], a custom-written linear program solver enabled a new high-precision calculation of critical exponents in the 3d Ising CFT, surpassing what was possible with out-of-the-box solvers like Mathematica [41] (used in the original work [1]), GLPK [42], and CPLEX [43]. 1nfortunately, linear programming is not applicable to systems of multiple correlators or operators with spin.These more complicated cases can be attacked with semidefinite programming [17,9], for which previous studies have relied on the solvers SDPA [44,45] and SDPA-GMP [46].The study [9], in particular, pushed SDPA-GMP to its limits, with each optimization taking up to 2 weeks.By contrast, SDPB can perform the same optimization in 1-3 CPU-hours, or 4-12 minutes on a 16-core machine.
A general bootstrap problem can be approximated as a particular type of semidefinite program called a "polynomial matrix program" (PMP).SDPB implements a variant of the well-known primal-dual interior point method for semidefinite programming [47][48][49][50], specialized for PMPs.Specialization and parallelization are SDPB's advantages.We are optimistic that better designs and algorithms can be brought to bear in the future.
In section 2, we describe PMPs and the design of SDPB.This discussion is mostly independent of the conformal bootstrap and should be comprehensible to readers without a physics background.
In section 3, as an application of SDPB, we set a new world-record for precision of critical exponents in the 3d Ising CFT, using multiple correlators as in [9].Readers interested solely in physics can skip to this section.We conclude in section 4.

Polynomial Matrix Programs
SDPB solves the following type of problem, which we call a polynomial matrix program (PMP).Consider a collection of symmetric polynomial matrices P n j,11 (x) . . .P n j,1m j (x) . . . . . . . . .P n j,m j 1 (x) . . .P n j,m j m j (x) labeled by 0 ≤ n ≤ N and 1 ≤ j ≤ J, where each element P n j,rs (x) is a polynomial.Given b ∈ R N , we would like to maximize b • y over y ∈ R N , such that M 0 j (x) + N n=1 y n M n j (x) 0 for all x ≥ 0 and 1 ≤ j ≤ J. (2. 2) The notation M 0 means "M is positive semidefinite." As we review in section 3, a wide class of optimization problems from the conformal bootstrap can be written in this form.Conveniently, PMPs can be translated into semidefinite programs (SDPs) and solved using interior point methods.In the next few subsections, we perform this translation and describe an interior point algorithm for solving general SDPs.Subsequently, we make this algorithm more efficient by exploiting special structure in PMPs.

Translating PMPs into SDPs
where (2.4) Here, S K is the space of K × K symmetric real matrices, and Tr(A * Y ) denotes the vector (Tr(A 1 Y ), . . ., Tr(A P Y )) ∈ R P .The SDP (2.3) is similar to those treated by solvers like SDPA-GMP, except that it includes the variables y ∈ R N , called "free variables" because they are not constrained to be positive.The matrix C will eventually be set to zero, but we include it for generality.
The first step is to relate positive semidefiniteness of polynomial matrices to positive semidefiniteness of a single matrix Y .Let q 0 (x), q 1 (x), . . ., be a collection of polynomials with degrees 0, 1, . . .(for example, monomials q m (x) = x m ), and define the vector q δ (x) = (q 0 (x), . . ., q δ (x)).We call q m (x) a "bilinear basis" because products q m (x)q n (x) span the space of polynomials.In particular, any polynomial P (x) of degree d can be written where and Y 1 , Y 2 are (underdetermined) symmetric matrices For a symmetric m×m polynomial matrix M(x) of degree d, we can apply this construction to each element, where now Y 1 acts on R δ 1 +1 ⊗ R m and Y 2 similarly acts on R δ 2 +1 ⊗ R m , and each trace is over the first tensor factor.
The translation of PMPs into SDPs relies on the following theorem: Theorem 2.1.M(x) is positive semidefinite for x ∈ R + if and only if it can be written in the form (2.8) for some positive semidefinite Y 1 and Y 2 .
Proof.One direction is simple: choose a vector v ∈ R m and consider the pairing Suppose Y 1 , Y 2 are positive semidefinite and x ≥ 0. Then since Q δ 1 (x) ⊗ vv T and xQ δ 2 (x) ⊗ vv T are both positive semidefinite, it follows that v T M(x)v ≥ 0. The other direction is less trivial.It has been proven directly in [51] and also follows as a consequence of the "Biform Theorem" of [52], using the substitution x = y2 and results of [53]. 2   Theorem 2.1 lets us rewrite our PMP constraints (2.2) in terms of a collection of positive semidefinite matrices Y 1 , . . ., Y 2J .We equate each polynomial matrix M n j (x) to an expression of the form (2.8).Individual matrix elements P n j,rs can be isolated by taking the trace over R m j with symmetrized unit matrices This gives a set of polynomial equalities (2.10) Equality between polynomials of degree d j is equivalent to equality at d j + 1 points.Thus, evaluating (2.10) at points x 0 , . . ., x d j , we obtain a set of affine relations between the y n and Y j ,3 Let us group the Y j 's into a single block-diagonal positive semidefinite matrix (2.12) The constraints (2.11) now take the form where p runs over all tuples (j, r, s, k) satisfying 0 ≤ r ≤ s < m j , and 0 ≤ k ≤ d j .The matrices A p , B, C and vector c are given by c (j,r,s,k) = P 0 j,rs (x k ), (2.16) C = 0. (2.17) This completes the translation of our PMP into an SDP (2.3).Note that the matrices A p are far from generic.Exploiting this fact will help us solve PMPs much more efficiently than a generic SDP.SDPB is specifically designed to solve SDPs with constraint matrices of the form where B b (M) denotes the block-diagonal matrix with M in the b-th block and zeros everywhere else, and the sets blocks j are disjoint for different j.In the example above, we have (2.20)

Semidefinite Program Duality
Duality plays an important role in our interior point algorithm, so let us briefly review it.The problem (2.3) is related by duality to the following "primal" optimization problem: We refer to the problem (2.21) as P for "primal" and (2.3) as D for "dual."We say that P is "feasible" if there exist x, X satisfying the constraints (2.21).Similarly, D is feasible if there exist y, Y satisfying the constraints (2.3).The "duality gap" is defined as the difference in primal and dual objective functions, c • x − Tr(CY ) − b • y.
For our purposes, the statement of duality is as follows: Theorem 2.2 (Semidefinite Program Duality).Given a feasible point (x, X) of P and a feasible point (y, Y ) of D, the duality gap is nonnegative.If the duality gap vanishes, then (x, X) and (y, Y ) are each optimal solutions of P and D, and furthermore XY = 0.
Proof.Suppose we have feasible solutions (x, X) and (y, Y ).The duality gap is given by where nonnegativity follows because X and Y are positive semidefinite.Now suppose Tr(XY ) vanishes.Clearly this implies XY = 0 identically (this condition is called "complementarity").Because c • x is bounded from below by the dual objective Tr(CY ) + b • y and also equal to the dual objective, the point (x, X) must be optimal.Similarly, since Tr(CY ) + b • y is bounded from above by the primal objective c • x and also equal to the primal objective, the point (y, Y ) must be optimal as well.
Unlike in linear programming, there is no guarantee that either P or D will attain their respective optima, or that the duality gap will vanish.For this, we need additional regularity assumptions.One of them is Slater's condition, which says that the duality gap vanishes if there exist strictly feasible solutions to the primal and dual constraints -i.e.solutions where X, Y ≻ 0 are positive-definite.Slater's condition is generic in the sense that a small perturbation of a feasible but not strictly-feasible problem will typically produce a strictly-feasible problem.

An Interior Point Method
The idea behind primal-dual interior point methods is to solve both the primal and dual equations simultaneously to find an optimal point q = (x, X, y, Y ).As we saw in the previous subsection, the optimum is (generically) achieved by a pair of positive semidefinite matrices X, Y satisfying the "complementarity" condition XY = 0. Most algorithms work by deforming this condition to for some nonzero µ, where I is the identity matrix.The constraints together with (2.23) then have a unique family of solutions called the "central path:" q(µ) = (x(µ), X(µ), y(µ), Y (µ)) indexed by µ ∈ R + .By following the central path from positive µ towards µ = 0, we can find the optimum of the original problem.
In practice, instead of moving precisely along the central path, we use the following strategy.Consider an initial point q = (x, X, y, Y ) such that X, Y are positive semidefinite.Our goal is to decrease Tr(XY ) and move towards the constraint locus while maintaining positive semidefiniteness.
• Set µ = βTr(XY )/K for some β < 1, where K is the number of rows of X.
• Use Newton's method to compute a direction dq = (dx, dX, dy, dY ) towards the central path with the given µ.
• Take a step along dq, taking care not to violate the positive semidefiniteness of X, Y .This should result in a reduction of Tr(XY ) by roughly a factor of β.
This is essentially Newton's method with a moving target.
An important advantage of this method is that the initial starting point (x, X, y, Y ) need not satisfy any of the equality constraints in (2.3) and (2.21).As long as we start with positive semidefinite X, Y , and the problem is feasible, the above method will converge to a point where the equality constraints are satisfied.

Newton Search Direction
Let us describe a single Newton step in more detail.The direction dq is defined by replacing q → q + dq and solving the constraint equations (2.3, 2.21) and complementarity equation (2.23) at linear order in dq, The residues measure the failure of the current point to satisfy the constraints.These residues will decrease with each Newton step.The linearized equations (2.24) can then be written ) where Z = X −1 (P Y − R) and the "Schur complement" matrix S is defined by To find the search direction, we first solve (2.26) for dx, dy, and then plug into (2.27) and (2.28) to determine dX, dY .Naïvely applying (2.28) leads to a dY that is not necessarily symmetric, which would take us outside of the domain of definition of Y .Several solutions to this problem have been proposed [54].Our approach, following [49,55,56], will be to symmetrize dY by hand, replacing (2.28) with (2.30)

Mehrotra Predictor-Corrector Trick
The most expensive operations in the search direction calculation are forming the Schur complement matrix S and solving the Schur complement equation (2.26).We'd like to perform them as rarely as possible.A simple modification to the naïve Newton's method, due to Mehrotra [57], allows us to get closer to the central path while reusing S and most of the work done in solving (2.26).
The rough idea is to solve the constraint and complementarity equations at higher order.This proceeds in two steps, called the "predictor" and "corrector" steps, respectively.For the predictor step, we compute a direction as described above, which we call dq p = (dx p , dX p , dy p , dY p ).We then replace the linearized complementarity equation (2.24) with and re-solve to obtain a corrector direction dq c .In the corrector step, we may (optionally) use a smaller deformation parameter µ → µ c to get closer to Tr(XY ) = 0. Note that the replacement of (2.24) with (2.31) does not change the Schur complement matrix S ij , so it can be reused, together with any matrix decompositions performed in solving (2.26).

Termination Conditions
We say a point q is "primal feasible" if the residues p, P are sufficiently small.Similarly, the solution is "dual feasible" if the residue d is sufficiently small.The precise conditions are primal feasible: where primalErrorThreshold ≪ 1 and dualErrorThreshold ≪ 1 are parameters chosen by the user.
An optimal point should be both primal and dual feasible, and have (nearly) equal primal and dual objective values.Specifically, let us define dualityGap as the normalized difference between the primal and dual objective functions where dualityGapThreshold ≪ 1 is chosen by the user.

Complete Algorithm
Our complete interior point algorithm is as follows.
1. Choose an initial point q = (x, X, y, Y ) = (0, Ω P I, 0, Ω D I) where Ω P and Ω D are real and positive.This point probably does not satisfy the constraints.
2. Compute the residues (2.25) and terminate if q is simultaneously primal feasible, dual feasible, and optimal.(Sometimes we may wish to use a different termination criterion, see below.) 3. Let µ = Tr(XY )/K and compute the predictor deformation parameter µ p = β p µ where β p = 0 if q is primal and dual feasible; β infeasible otherwise. (2.36) Here, β infeasible ∈ (0, 1) is chosen by the user.
where α(M, dM) is the largest positive real number such that M + α(M, dM)dM is positive semidefinite, and γ ∈ (0, 1) is a parameter chosen by the user.
8. Replace and go to step 2. Note that the replacement (2.40) is guaranteed to preserve positive semidefiniteness of X and Y .
If the current point is close enough to a primal (or dual) feasible region, the step-length α P (α D ) in (2.39) can be exactly 1.When this occurs, the replacement (2.40) will exactly solve the primal (dual) equality constraints, up to numerical errors.This follows from linearity of the equality constraints, together with the fact that symmetrizing Y in (2.30) has no effect on the constraints.In cases where we only care about primal or dual feasibility, the iteration can be stopped here, see section 3.4.

Specialization to Polynomial Matrix Programs
As mentioned in section 2.4.2, the most expensive part of the search direction calculation is computing the Schur complement matrix S pq = Tr(A p X −1 A q Y ) and inverting (2.26) to obtain dx, dy.In this section, we will specialize to PMPs, and study how these calculations can be made more efficient.For similar optimizations, see [58].

Block Structure of the Schur Complement Matrix
Recall that for PMPs, the matrices A p are given by (2.18) with the index p running over tuples (j, r, s, k) satisfying 0 ≤ r ≤ s < m j , and 0 ≤ k ≤ d j .Since X and Y have the block structure (2.12), the Schur complement matrix S p 1 p 2 is block-diagonal: it has nonzero entries only if The size of the j-th block dim S (j) is equal to the number of choices for (r, s, k), Now consider equation (2.26), We could solve it using an LU (lower triangular × upper triangular) decomposition of T , but this is extremely expensive, taking cubic time in dim T = j dim S (j) + N.
We should use the block structure of T to our advantage.Let S = LL T be a Cholesky decomposition of S (which can be computed efficiently for each block S (j) = L (j) L (j)T ), and consider the decomposition The outer matrices on the right hand side are triangular, and can be solved efficiently by forward/backward-substitution.Meanwhile, the matrix Q ≡ B T L −T L −1 B typically has a much smaller dimension than S, dim Q = N ≪ dim T , so the middle block-matrix can be easily solved using a Cholesky decomposition. 4nfortunately, the decomposition (2.44) is numerically unstable when S is ill-conditioned.Indeed, suppose S has a very small eigenvalue (so L does too), while the full matrix T does not.Then quantities like L −1 B that appear in (2.44) will have large entries which must nearly cancel when recombined into a solution of (2.43).Near-cancellation reduces numerical precision.
These problems stem from the off-diagonal blocks of T , which ultimately come from the free variables y in our semidefinite program.Several authors have considered the problem of efficiently and stably solving semidefinite programs with free variables, with no obvious consensus [59].For example, [60] suggests eliminating free variables by explicitly solving the primal constraint B T x = b and taking appropriate linear combinations of the matrices A p .However this procedure destroys the sparsity structure of S, making it no longer block diagonal and forcing us to use an expensive full Cholesky decomposition.
Preserving the structure of S and T is paramount.The simplest way to stabilize (2.44) is to increase the precision of the underlying arithmetic.In practice, this appears to be necessary anyway for larger-scale bootstrap problems, see appendix A. For additional stabilization, we use an old trick of adding low-rank pieces to S to make it better-conditioned.Suppose u i are vectors, each with a single nonzero entry, such that S ′ ≡ S + i u i u T i = S + UU T has no small eigenvalues. 5Note that S ′ differs from S in only a few diagonal entries -in particular it has the same block structure.Now let us replace (2.43) with the system By solving for z and substituting back in, it is easy to see that (2.45) is precisely equivalent to (2.43).However, the advantage is that because S ′ is well-conditioned, a decomposition like (2.44) is numerically stable.Indeed, defining B ′ = (B U), we have where S ′ = L ′ L ′T .Because Q ′ is no longer necessarily positive semidefinite, we are forced to use an LU decomposition to invert the middle matrix in (2.46), which is slightly more expensive than Cholesky decomposition.Fortunately, Q ′ is usually small enough that this cost is inconsequential.If T itself is ill-conditioned, then this will manifest as instabilities when we try to LU decompose Q ′ .In this situation, there is little we can do to avoid imprecision.

Computing the Schur Complement Matrix
Now that we know what to do with the Schur complement matrix S, let us compute it efficiently.The fact that A p has low rank is helpful.This explains why we chose to evaluate the polynomial equalities (2.10) at discrete points x k in (2.11).Matching polynomial coefficients on each side, as done in [17,7,9], leads to higher-rank matrices A p .The helpfulness of low-rank constraints in solving SDPs, and their appearance in polynomial optimization, is well known [58].
Recall that X and Y are block diagonal (2.12), with each block X b acting on a tensor product of the form R δ+1 ⊗ R m .Let X (r,s) b ∈ R (δ+1)×(δ+1) denote the (r, s)-th block of X b in the second tensor factor, which acts on R δ+1 .
Since S is block diagonal, we need only compute elements with j 1 = j 2 = j.We have Thus, instead of performing repeated matrix multiplications to calculate Tr(A p 1 X −1 A p 2 Y ), we can precompute the bilinear pairings and plug them into (2.48).Whereas forming S is often the most expensive operation in lessspecialized solvers, the method outlined here makes it subdominant to other computations, see appendix B.1.

Computing Step Lengths
To find step lengths α P , α D , we must be able to compute α(M, dM) where M is a positive semidefinite matrix and α(M, dM) is the largest positive real number such that M + α(M, dM)dM is positive semidefinite.Let M = LL T be a Cholesky decomposition of In SDPB, we compute all the eigenvalues of L −1 dML −T using a QR decomposition and then simply take the minimum.Some solvers, like SDPA, implement the more efficient Lanczos method [61] for computing the minimum eigenvalue.In practice, the step-length calculation is a small part of the total running time, so we have neglected this optimization.

Implementation
SDPB is approximately 3500 lines of C++.It uses the GNU Multiprecision Library (GMP) [62] for arbitrary precision arithmetic, and MPACK [63] for arbitrary precision linear algebra.The relevant MPACK source files are included with SDPB, with some slight modifications: • The Cholesky decomposition routine Rpotrf has been modified to implement the stabilization procedure described in footnote 5.
• Some loops in the LU decomposition routine Rgetrf have been parallelized.
Previous experience shows that high-precision arithmetic is important for accurately solving bootstrap optimization problems.It is not fully understood why.The naïve reason is that derivatives ∂ m z ∂ n z g ∆,ℓ (z, z) of conformal blocks vary by many orders of magnitude relative to each other as ∆ varies.It is not possible to scale away this large variation, and answers may depend on near cancellation of large numbers.In practice, the matrix manipulations in our interior point algorithm "leak" precision, so that the search direction (dx, dX, dy, dY ) is less precise than the initial point (x, X, y, Y ).By increasing the precision of the underlying arithmetic, the search direction can be made reliable again.This strategy (which we adopt) comes at a cost of increased runtime and memory usage.Better strategies for dealing with numerical instabilities in bootstrap problems could bring enormous gains.
SDPB is parallelized with OpenMP [66].Because most matrices appearing in the interior point algorithm are block-diagonal, most computations are "embarrassingly parallel:" different blocks can be distributed to different threads.(The most prominent exception is the LU decomposition of Q ′ , which is why Rgetrf was modified.)Consequently, performance scales nearly linearly with the number of cores, as long as the number of matrix blocks is sufficiently large.This is usually the case for interesting bootstrap problems, where J (which sets the number of blocks) is typically much larger than the number of available cores.It should be possible to achieve favorable scaling up to dozens or even hundreds of cores using MPI and more careful memory management.Further scaling should be possible with more fine-grained parallelization.
SDPB is available online at https://github.com/davidsd/sdpbunder the MIT license.The source code is carefully commented and written for readability (to the extent that C++ code is ever readable).We hope this will encourage customization and improvement.

Example Application: 3d Ising Critical Exponents 3.1 A 3d Ising Optimization Problem
Bootstrap optimization problems can be naturally approximated as PMPs [17,9].In this section, we review the PMP for the system of correlators { σσσσ , σσǫǫ , ǫǫǫǫ } in the 3d Ising CFT.We will be brief.Much more detail is given in [9].Associativity of the Operator Product Expansion (OPE) for { σσσσ , σσǫǫ , ǫǫǫǫ } implies the consistency condition Here, O + runs over Z 2 -even operators of even spin and O − runs over Z 2 -odd operators of any spin.We have separated out the unit operator.∆ and ℓ are the dimension and spin of O, respectively.The object V −,∆,ℓ is a 5-vector and V +,∆,ℓ is a 5-vector of 2 × 2 matrices with entries that are functions of conformal cross-ratios u and v, The g (v, u) are conformal blocks, which are known special functions.
The OPE coefficients λ σσO , λ σǫO , λ ǫǫO and dimensions ∆ are not known a priori.Nonetheless, we can constrain them by understanding when it is possible for the terms in (3.1) to sum to zero.To do this, consider a 5-vector of functionals α = (α 1 , . . ., α 5 ), where each α i acts on the space of functions of u and v. Acting on (3.1) with α gives The bootstrap logic, in the spirit of [1], is as follows.First we make an assumption about which ∆, ℓ appear in (3.4).We then search for a functional α such that , for all Z 2 -even operators with even spin, α • V −,∆,ℓ ≥ 0, for all Z 2 -odd operators of any spin. (3.5) The OPE coefficients λ ijk are real in a unitary CFT.Thus, if α exists, then it is impossible to satisfy the consistency condition (3.1), and our assumption is ruled out.By making different assumptions and searching for functionals α, we can map out the space of allowed ∆.

Approximation as a PMP
The conditions (3.5) define a feasibility problem with an infinite number of semidefiniteness constraints (one for each ∆, ℓ).To obtain a PMP, we choose a particular type of functional where Although the bootstrap logic does not depend on the types of functionals considered, only functionals of the form (3.6) lead to a PMP.Other types of functionals require different optimization methods.
There are two important differences between (3.9) and our original PMP (2.2): 1.In (2.2) we have a finite number of positive semidefiniteness conditions j = 1, . . ., J, whereas here we have an infinite number since ℓ can be any nonnegative integer.In practice, we include spins ℓ up to some large but finite ℓ max .As long as ℓ max is large enough, a functional obtained by solving the problem with ℓ ≤ ℓ max should also satisfy positive semidefiniteness for spins ℓ > ℓ max .The proper choice of ℓ max depends on the problem at hand, see appendix A. See [27] for a more careful analysis.
2. In (2.2) we are trying to optimize an objective function b • y, whereas here we are only interested in feasibility.To determine feasibility, we can pick the trivial objective function b = 0 and run our interior point algorithm until it becomes dual feasible.

Setting Up SDPB
The natural objects entering our calculation are (approximate) derivatives of conformal blocks χ ℓ (∆)p ∆ 12 ,∆ 34 ;mn ℓ (∆), as opposed to just the polynomials p ∆ 12 ,∆ 34 ;mn ℓ (∆).Removing positive factors does not affect positive semidefiniteness, but it does affect the scaling of the resulting SDP.Restoring quantities to their "natural" size can improve numerical stability and performance. 6SDPB provides a few different ways to implement this rescaling.
As we saw in section 2.2, translating a PMP into an SDP requires a bilinear basis q m (x).SDPB allows a choice of bilinear basis q (j) m (x) for each j = 1, . . ., J.For bootstrap problems, we take q (j) m (x) to be orthogonal polynomials with respect to the norm where ℓ is the spin corresponding to j.The change of basis between orthogonal polynomials with respect to •, • (j) and monomials x m (used in previous bootstrap applications of SDP) is extremely ill-conditioned at high degree.So although the choice of q (j) m is unimportant in principle, it can have a dramatic effect on numerical stability. 7DPB also requires a set of sample points x (j) k at which to evaluate polynomials, as well as scaling factors s (j) k that modify the constraints (2.15, 2.16) as follows: Additionally, the A (j,r,s,k) are given by (2.18) with This s k -dependent rescaling gives an isomorphic, but potentially more numerically stable SDP.For bootstrap problems, it is natural to pick s k ) where ℓ corresponds to j.The x (j) k can be any sequence of distinct points.A natural choice are zeros of one of the q (j) m of sufficiently high degree.
To summarize, SDPB depends on the following input: • for each j = 1, . . ., J: polynomial matrices M 0 j (x), . . ., M N j (x) of maximum degree d j , -bilinear bases q SDPB reads this data in an XML format described in the manual.A Mathematica package that translates the above data from Mathematica expressions into XML is included with the source distribution.An example 2d bootstrap computation is also included.More details about SDPB's input and output formats and its various settings can be found in the manual.

Results
As an application of SDPB, let us improve upon the determinations of 3d Ising critical exponents in [5,9].We make an exclusion plot for the operator dimensions (∆ σ , ∆ ǫ ) as follows.Fix (∆ σ , ∆ ǫ ) and use SDPB to determine if the PMP (3.9) is dual-feasible, i.e. whether there exist (y, Y ) satisfying their associated constraints.If the PMP is dual-feasible, then the given (∆ σ , ∆ ǫ ) are excluded.If the PMP is not dual-feasible, then we cannot conclude anything about (∆ σ , ∆ ǫ ).By scanning over different points, we map out the excluded region in (∆ σ , ∆ ǫ ) space.
To determine dual feasibility, we use a vanishing objective function b = 0 and run SDPB with the option --findDualFeasible.This terminates the solver if (y, Y ) are found satisfying their constraints to sufficient precision (i.e. if dualError < dualErrorThreshold). 8In practice, if SDPB finds a primal feasible solution (x, X) after some number of iterations, then it will never eventually find a dual feasible one.Thus, we additionally include the option --findPrimalFeasible, which terminates the solver whenever primalError < primalErrorThreshold.The termination status of SDPB then determines whether a point is allowed or not: The precise SDPB options used for the computations in this work are described in appendix A.  1: Allowed region for a Z 2 -symmetric 3d CFT with two relevant scalars, computed using SDPB with the system of correlators σσσσ , σσǫǫ , and ǫǫǫǫ .The blue regions correspond to Λ = 19,27,35,43, in decreasing order of size.The larger black rectangle shows the current most precise Monte Carlo determinations of critical exponents in the 3d Ising CFT [67].The smaller black rectangle shows the estimate for (∆ σ , ∆ ǫ ) using c-minimization at Λ = 41 for the single correlator σσσσ [5].
In figure 1, we plot the allowed regions for different numbers of derivatives labeled by Λ = 19, 27, 35, 43, 9 corresponding to functionals α of dimension 275, 525, 855, and 1265, respectively. 10We focus on (∆ σ , ∆ ǫ ) near the 3d Ising CFT, leaving wider exploration to the future.The allowed region is an island that shrinks rapidly with increasing Λ. 11 The largest island, corresponding to Λ = 19 is the same as the allowed region in figure 5 of [9].We can estimate the point towards which the islands shrink as follows.Let (a Λ , b Λ ) be the bottom-left point of the Λ-allowed island, and similarly let (c Λ , d Λ ) be the top-right point.Define 9 n max = 10, 14, 18, 22 in the notation of [9]. 10 A performance analysis for different values of Λ is given in appendix B.2.
11 Each allowed region plotted in this work was computed by testing a grid of points and fitting curves to the boundary between allowed and disallowed gridpoints.The raw gridpoint data is available on request.and let r x , r y be the minima of E x , E y respectively.Our estimate is12 ≈ (0.5181478(5), 1.412617(4)), (3.16)where the errors are given by E x (r x ), E y (r y ).
The dimensions (∆ σ , ∆ ǫ ) were estimated in [5] using the conjecture that the 3d Ising CFT minimizes c ≡ ∆ 2 σ /λ σσTµν subject to the constraints of unitarity and crossing symmetry of σσσσ .This conjecture, called "c-minimization," is expected to be equivalent to the assumption that the 3d Ising CFT lives precisely at the kink in the dimension bound coming from the single correlator σσσσ .Although unproven, c-minimization's advantage is that it allows one to estimate (∆ σ , ∆ ǫ ) using a single scalar correlator σσσσ .Bootstrap computations for a single scalar correlator can be made relatively efficient with a modified primal simplex algorithm [5,70].
An advantage of multiple correlators is that it is possible to impose the condition that σ is the only relevant Z 2 -odd operator, causing the allowed region in (∆ σ , ∆ ǫ )-space to become a closed island, independent of auxiliary assumptions [9].Because our Λ = 43 island lies within the error bars of [5], our results verify c-minimization to the precision achieved in [5]. 13Our results also give further evidence for the conjecture that the 3d Ising CFT is the unique Z 2 -symmetric 3d CFT with two relevant scalars and ∆ σ 0.6.(The precise condition on ∆ σ , ∆ ǫ depends on the shape of the allowed region further away from the 3d Ising point.)It would be very interesting to prove these conjectures analytically, perhaps by showing that the island in figure 2 shrinks to a point as Λ → ∞.An alternative is that the conjectures are still true, but one needs information from other four-point functions to prove them.The analysis of [9] did not use permutation symmetry of the OPE coefficient λ σσǫ = λ σǫσ . 14Including this constraint leads to an additional modest reduction in the allowed region, 15 which we plot in figure 2 at Λ = 43.The resulting island gives a rigorous determination ∆ σ = 0.518151 (6), ∆ ǫ = 1.41264 (6), which is 5-10 times more precise than the Monte Carlo results of [67].We summarize the comparison to Monte Carlo in figure 3.

Discussion
With SDPB, we have significantly improved the precision of (∆ σ , ∆ ǫ ) in the 3d Ising CFT.Our numerics indicate that the window of allowed dimensions may shrink to a point in the limit of infinite computer time.In other words, they suggest that consistency of the correlators allowed region for Λ = 43, using three-point symmetry  { σσσσ , σσǫǫ , ǫǫǫǫ }, together with the assumption that σ and ǫ are the only relevant scalars in the theory, may uniquely fix the dimensions (∆ σ , ∆ ǫ ).This conjecture could be more tractable analytically than trying to solve the full CFT consistency conditions.
There are many more 3d Ising observables to explore.For example, the coefficient f ǫǫǫ should be computable, and it would be interesting to compare with the recent Monte Carlo prediction [72].It will also be important to consider larger systems of 3d Ising correlators.
However, SDPB should also enable wider exploration of new correlators and diverse theories.SDPB is already being used in several bootstrap studies that would have previously been difficult or impossible [71,73].An exciting direction that may now be accessible is studying a four-point function of stress-tensors in 3d CFTs.
In addition to the four-point function bootstrap, semidefinite programming has also recently been applied to the "modular bootstrap" in 2d CFTs [74,75].SDPB is equally applicable to modular bootstrap computations, since they too can be phrased in terms of polynomial matrix programs.
From the computing point of view, there are many opportunities for improvement.For example, it should be possible to parallelize SDPB up to hundreds of cores, which could lead to even more precise calculations, and (just as importantly) easier exploration.Very different algorithms, like Second Order Conic Programming (SOCP), cutting plane methods, or constrained nonlinear optimization may also be applicable.
The revival initiated in [1] is still young, and the technology (both analytical and numerical) is evolving rapidly.Current techniques are likely not maximally efficient, and it will be important to consider other methods, from new algorithms and optimization tools to conceptually different approaches.We are optimistic that much more will be possible.

A Choices and Parameters
The PMP for { σσσσ , σσǫǫ , ǫǫǫǫ } in the 3d Ising CFT depends on the following choices: • An integer Λ specifying which derivatives to include in the functional α.These are given by ∂ m z ∂ n z with m ≥ n and m + n ≤ Λ. 16 Because of the symmetry/antisymmetry of F ij,kl ±,∆,ℓ (u, v) under u ↔ v, some of these derivatives vanish identically.Including only non-vanishing derivatives, the dimension of α is • An integer κ controlling the accuracy of the approximation for conformal blocks (3.7).The positive prefactor is where r * = 3 − 2 √ 2 is the radius of the point z = z = 1 2 in the radial coordinates of [76,77].The dimensions ∆ * i are special values below the unitarity bound, so that the product i (∆ − ∆ * i ) is positive for unitary theories.The approximation (3.7) can be systematically improved by including more poles (∆ − ∆ * i ) −1 and increasing the degree of p ∆ ij ,∆ kl ;mn ℓ (∆).Our choice of poles is Smaller κ means smaller-degree polynomials and shorter runtimes.Larger κ is needed to get an accurate approximation for conformal blocks.We choose κ by computing bounds with successively larger values of κ until the results stabilize.Our final values are conservative: smaller κ may still give sufficient accuracy.Derivatives of conformal blocks were computed using the recursion relation in [9] to order r 90 (far greater accuracy than needed).
• A set of spins S = {ℓ 1 , . . ., ℓ L } to include.If not enough spins are included, the solver may find a functional α that violates a positive semidefiniteness constraint for some spin.Because derivatives of conformal blocks converge as a function of spin, in practice it is sufficient to include a finite number of spins to ensure α satisfies the constraints for all spins (as can be verified post-hoc by testing α on constraints that were not explicitly included).This sufficient number of spins grows with Λ.Our choices are given in (A.4).• The Schur complement matrix is dense, so each of its elements must be computed individually.Computing X −1 A q Y requires N + 2J(md) 3 multiplications, since it involves a block-diagonal dense matrix multiplication X −1 ×(A q Y ) (the linear term in N comes from multiplying 1 × 1 blocks).This must be repeated Jm 2 d 2 times: once for each q.Now, the matrices A p are typically sparse, with O(md) entries.Thus, evaluating all the traces Tr(A p X −1 A q Y ) requires approximately (N + 2J(md) 3 ) Jm 2 d 2 + 1 2 ( Jm 2 d 2 ) 2 (md) multiplications.These steps dominate the running time for the computations in [9].
• Because the Schur complement matrix is dense, it must be inverted using a full Cholesky decomposition, which requires 1 3 ( Jm 2 d 2 ) 3 multiplications.
For SDPB, computing S takes negligible time.The most important steps are in solving the Schur complement equation: • Computing the Cholesky decomposition S = LL T takes 1 3 J( m 2 d 2 ) 3 multiplications, since it can be done block-wise.
• Forming Q = (L −1 B) T (L −1 B) requires 1 2 N 2 (J m 2 d 2 ) multiplications.This step dominates the running time for the computations in this work.
• Computing the LU decomposition of Q requires 2  3 N 3 multiplications.A full analysis of the complexity in Λ would require determining the correct asymptotics of each quantity (including SDPB parameters like precision), which may depend on the computation at hand.In table Table 2: Runtimes for a single feasibility computation in the 3d Ising CFT using the correlators { σσσσ , σσǫǫ , ǫǫǫǫ }, as described in [9] and appendix A. Average runtimes are different depending on whether a spectrum is disallowed (dual feasible) or allowed (dual infeasible).All times are in CPU-hours (for SDPB, this means the actual runtime is multiplied by maxThreads, which was 16 for most of the computations in this work).Approximate SDPA-GMP times are from [9].All computations were performed on 3.3GHz 64-bit Intel Xeon Sandy Bridge processors.The column "mul./iter."gives the approximate number of multiplications per iteration, calculated according to the discussion in subsection B.1.(To estimate running time from the number of multiplications per iteration, one needs to take into account precision, which also increases with Λ.)The SDPA-GMP computations with Λ > 19 have not been attempted.

B.2 Running Time for 3d Ising Computations
Let us begin by translating the PMP (2.2) into a more standard semidefinite program of the following form: maximize Tr(CY ) + b • y over y ∈ R N , Y ∈ S K , such that Tr(A * Y ) + By = c, and Y 0,

Figure 2 :
Figure2: Allowed region for a Z 2 -symmetric 3d CFT with two relevant operators, computed with SDPB at Λ = 43.The light-blue region is a zoom of the smallest region in figure1.The darker-blue region additionally uses symmetry of the OPE coefficients λ σσǫ = λ σǫσ .The black rectangle shows the estimate for (∆ σ , ∆ ǫ ) using c-minimization at Λ = 41[5].