Solving LargeScale Optimization Problems Related to Bell’s Theorem
Abstract
Impossibility of finding local realistic models for quantum correlations due to entanglement is an important fact in foundations of quantum physics, gaining now new applications in quantum information theory. We present an indepth description of a method of testing the existence of such models, which involves two levels of optimization: a higherlevel nonlinear task and a lowerlevel linear programming (LP) task. The article compares the performances of the existing implementation of the method, where the LPs are solved with the simplex method, and our new implementation, where the LPs are solved with an innovative matrixfree interior point method. We describe in detail how the latter can be applied to our problem, discuss the basic scenario and possible improvements and how they impact on overall performance. Significant performance advantage of the matrixfree interior point method over the simplex method is confirmed by extensive computational results. The new method is able to solve substantially larger problems. Consequently, the noise resistance of the nonclassicality of correlations of several types of quantum states, which has never been computed before, can now be efficiently determined. An extensive set of data in the form of tables and graphics is presented and discussed. The article is intended for all audiences, no quantummechanical background is necessary.
keywords:
Quantum Information, LargeScale Optimization, Interior Point Methods, MatrixFree Methods.url]http://www.maths.ed.ac.uk/ gondzio/
url]http://www.gruca.org/ \cortext[cor1]Corresponding author. Tel. +48 602153844; fax +48 585232056. Postal address: Institute of Theoretical Physics and Astrophysics, University of Gdańsk, 80952 Gdańsk, ul. Wita Stwosza 57, Poland
1 Introduction
We discuss a new approach to tackling a certain class of large scale optimization problems arising in quantum information science and linked with Bell’s Theorem Bell (1964). These problems were described in detail in Kaszlikowski et al. (2000); Gruca et al. (2010) and are, in fact, twolevel optimization problems. The higher level problem is a nonconvex nonlinear optimization task. It requires the solution of a sequence of linear programming problems (LPs for short). Each of the LPs is a lower level task. We will focus on LP, as it has been identified as the major computational challenge. Reference Gruca et al. (2010) reports a range of problems which have been solved, with the small and medium scale LPs being solved using the GLPK myb (2011) simplex implementation. However, even for medium scale problems the CPU times were occasionally prohibitive. The solution time of a single problem with about 64,000 variables using GLPK sometimes exceeds 24 hours and the upper level optimization may require the solution of hundreds of such LP problems.
The solution of a single LP is a clear bottleneck in quantum information optimization. In this paper we report on the efforts to accelerate this part of optimization process. We make two contributions to achieve the goal:

we reformulate the problem by removing multiple redundant constraints from its original definition Gruca et al. (2010);

we replace the simplex method with a specialized variant of the matrixfree interior point method Gondzio (2012).
Although eliminating redundancy significantly reduces the number of LP constraints and produces more compact problem formulations, it is not sufficient to extend the applicability of the previous approach Gruca et al. (2010) which was based on the use of the simplex method. The reduced problems still defy both standard LP approaches: the simplex method and the interior point method. The use of the matrixfree variant of an interior point method provides a major advance because it allows the solution of these problems to be faster by orders of magnitude and guarantees a major reduction of the memory required to store them.
The aim of this article is twofold. Firstly, it is to describe the advantages of the interior point method over the simplex method in the context of the described problem and secondly, to include a complete formulation of the problem of computing the minimal critical visibility, accessible also to researchers outside the quantum information domain. In the following two sections we cover the second aim by discussing the formulation of the problems without going into excessive details on quantum physics in order to make the paper accessible to nonspecialists. We then discuss the original approach and the method along with its advantages and drawbacks, then we present our revised approach and a comparison between the two, thus covering the first aim of the article. Finally, we draw some conclusions and suggestions for potential future investigation.
2 Motivation
Quantum information science, also called quantum informatics, is an attempt to harness paradoxical aspects of quantum physics in the form of highly nonclassical effects in information transfer and processing, aiming at breaking classical limits (in the form of Shannon’s information theory or principles of ChurchTuring computation). It operates on quantum entangled states of particles and aims to provide protocols which are faster, more robust, secure or in any other way better than the ones provided by classical informatics. In order to do this it must harness the nonclassical properties of the states it operates on. This has led to a need for a measure of nonclassicality of quantum entangled states. Such a measure is also important in practical applications where obtaining a pure state is impossible since it always comes with a bit of noise which reduces its nonclassicality. For our purposes, classicality will be synonymous with local realism or local hidden variable models Fine (1982) (as it was shown that any localrealistic approach can be effectively simulated by local hidden variable models). In short, local realism requires that no interaction can occur between particles that exceeds the speed of light, and that particles must have a preexisting value for any measurement before the measurement is made. We will now present a formulation of the problem of computing a certain type of threshold for nonclassicality of quantum correlations.
3 Problem formulation
3.1 Critical visibility
Let us consider an experiment on a quantum entangled state of two particles emitted from a source in distinct directions. These particles travel towards Alice on one side and Bob on the other. The names Alice and Bob refer to the observers who both perform measurements of observables on their particles. The outcome of any single measurement is either 0 or 1. In the simplest case, Alice and Bob choose from two observables each: for Alice and for Bob
Quantum mechanics provides us with tools to calculate probabilities of measurement outcomes, which are denoted
(1) 
The above is a probability of a situation by which we denote that Alice obtains outcome while measuring observable and Bob obtains outcome while measuring observable . Under local realism there must also exist an underlying global probability distribution
(2) 
where and ( and ) are hypothetical values of measurements performed by Alice (Bob), if she (he) chooses to measure or ( or ). Both Alice and Bob have to choose exactly one observable for their measurements on a given pair of particles. The measurements of, say, and may be quantummechanically incompatible so are impossible to measure together (due to partial or full complementarity). Thus quantum formalism does not give us any method of deriving (2). As a matter of fact, quantum mechanics rules out the possibility of the existence of and hence is a nonlocalrealistic theory: it is a nonclassical theory. The practical question for quantum information applications is: where is the border between the classical and the nonclassical?
If we now assume local realism is correct then the following equations hold
(3) 
It can be shown that for some quantum entangled states the marginal sums (3.1) cannot be satisfied. This is contained in Bell’s Theorem. We will now discuss a method
If a strongly nonclassical state violates the above marginal sums (3.1) with a quantummechanical probability distribution:
(4) 
then we can admix to it a fraction of a fully classical state called white noise Laskowski (2007) to make it satisfy them. By white noise we understand completely random results (with a classical model of two coin tosses, one by Alice and one by Bob). The resulting state, with the noise admixture, will have the following probability distribution
(5) 
where is the visibility of the original state . If we admix of white noise to even the most nonclassical state, then we will always satisfy (3.1). On the other hand, if then, for some nonclassical states, we will see (3.1) violated. For such states there exists a threshold visibility, known as the critical visibility, above which we will see violation of the marginal sums. Let us now state the above as a set of linear equations (constraints) which have to be satisfied. From (3.1) and (5) we have:
where the top four lines form 16 equalities (each pair accounts for four of them as ). The orthant defined by the above 16 marginal sums, the probability summation constraint and the bounds on the visibility and probabilities forms the feasible region of our LP. The following cost function
(7) 
has a maximum on the feasible region equal to the critical visibility. If we adopt the traditional formulation of the LP problem:
(8) 
then, for the problem stated above, we have:
Let us observe that the above LP depends solely on the set of observables that the observers choose from during the measurements. In our experimental situation any observable can be (and in practice is) effectively parametrised by a pair of angles so our rightmost probabilities become
(10) 
The details of why such a parametrization takes place can be found in Nielsen and Chuang (2000) so are omitted here. However, what is important is that the critical visibility can be defined as the following function of the angles parameterising the observables.
(11) 
3.2 Minimal critical visibility
A single evaluation of this function is the lowerlevel task, requiring the solution of the LP problem described above with the maximum of the cost function as the critical visibility function’s value. Our ultimate goal is to minimize (11), thus arriving at the minimal critical visibility, the nonclassicality measure that we seek. The minimization of the critical visibility function is the higherlevel nonconvex nonlinear optimization task mentioned in Section 1.
The following sections contain much discussion, both about individual LPs and the higherlevel optimization. Hence, we would like to clarify certain vocabulary here. We will from now on use the expression “solve an LP” interchangeably with “calculate the value of the critical visibility function”, and the expression “minimize the critical visibility function” interchangeably with “execute the higherlevel optimization procedure”. We will also say “solve a sequence of LPs leading to the minimum of the critical visibility function”, as this is exactly what the higherlevel optimization procedure does. It requests a number of critical visibility function values (each of which requires the solution of a separate LP) in order to minimize it.
This work is mainly concerned with the method of obtaining numerical results rather than the results themselves. Therefore, for simplicity we will only analyse a single quantum entangled state called the GreenbergerHorneZeilinger (GHZ) state described in Greenberger et al. (1989); Scarani and Gisin (2001).
3.3 Redundancy in the problem formulation
One of our first findings was that the LPs constructed as described above contain a significant proportion of redundant constraints which increases with the problem size. For the 65536x65536 problem, where the first number is the number of rows excluding the summation of probabilities and the second is the number of variables excluding the visibility, only 6561 rows are needed for full formulation of the LP problem. We found that, when constructing the matrix , it is possible to skip the rows according to the following rule. If is the number of possible measurement outcomes (in our case ), then:
Rule 1
A row is redundant if its last column contains probabilities of measurements for which any of the observers obtained a result equal to while measuring an observable with index higher than .
In our example this means that one can skip rows 6, 8, 11, 12, 14, 15, 16. It is also possible to skip the last row representing the summation of probabilities to one since it is never needed. So, since 8 of 17 rows of matrix can be removed, the rank of the problem matrix () is 9. The percentage of redundant rows in matrix is, therefore, slightly less than 50%. This percentage, however, increases rapidly with the problem size and allows us to get rid of more than 90% of the rows for sufficiently large problems. This finding significantly sped up calculations both with the matrixfree interior point method and with the simplex method (GLPK).
Let us observe that, typically, both the number of rows and the number of variables in our LPs are multiples of 2 so, from now on, we will use ‘k’ to mean 1024, so the 65536x65536 problem will become problem 64kx64k. Thus the problem name refers to its origin (since it is informative for researchers in the field of quantum information science) but the size quoted for the problem is that of the LP after elimination of redundancy. Furthermore, since all the problems analysed in this paper are square, we will skip one of the dimensions and problem 64kx64k will also be referred to as problem 64k.
4 Original approach
The original approach adopted by the authors of Gruca et al. (2010) used the implementation of the simplex method Dantzig (1963) from the GLPK library myb (2011) for solving the lowerlevel LP task, and the downhill simplex method (DSM) Nelder and Mead (1965) from the SciPy package Jones et al. (2011) as the higherlevel nonlinear optimizer. As shown in Gruca et al. (2012); Gruca (2008) the critical visibility function is continuous but not differentiable so a robust higherlevel optimizer was needed. The DSM is a good fit since, in the vast majority of cases, it can be applied to any continuous function. It is worth mentioning that these two optimization techniques (the simplex method for LP and the DSM) are extremely powerful tools that have been used successfully in a plethora of applications. However, both of them have occasional drawbacks. The simplex method may struggle when applied to degenerate LP problems Hall and McKinnon (2004) and the DSM is not guaranteed to converge to even a local minimizer McKinnon (1998). The implementation has since been revised and now uses a newer version of the GLPK library (4.47 instead of 4.31) and a different implementation of DSM Press et al. (2007). This new code is called steamroller2. It exploits the property that the LP problems being solved in sequence are related to one another and differ only in the last column of , allowing the simplex solver to hotstart its solution procedure.
A typical execution of the method chooses random angles as a starting point for the higherlevel DSM optimization procedure. These angles are used to generate the first LP and its solution yields the value of the critical visibility function for the starting point. The DSM then modifies the angles as it sees fit and another LP is then generated and solved to calculate the critical visibility function value for the new angles. The procedure thus generates a sequence of LPs, each of which is solved by the simplex method. Finally, after a couple of hundred of such steps, the minimum of the continuous nondifferentiable critical visibility function is found.
Such an approach had a number of advantages: the two algorithms were a good fit since the DSM always converged and often returned the global minimum outright, and the simplex method could take advantage of the similarities between the subsequent LPs. They were also very precise: the results could be matched with known theoretical values (up to 15 significant digits). This is very much the limit of the machine’s precision.
It was found, however, that the calculations slowed down when the problem size increased and the GLPK simplex implementation struggled to solve the LP problems larger than 16k. The simplex method for linear programming is a nonpolynomial algorithm. Indeed, there does not exist a polynomial bound on the number of iterations that it may need to solve a given LP and occasionally it may need to perform a huge number of iterations. Although in practice such situations are rare Hall and McKinnon (2004) the LPs which model the critical visibility seem to challenge the simplex method. The solution of a single LP may require hours (or days) of computations so, since a sequence of hundreds of such LPs may have to be solved, the whole approach is questionable.
We have made a radical change in the solution approach and decided to use an interior point method (IPM) rather then the simplex method to solve the LPs. Such is the improvement of efficiency delivered by the use of the matrixfree interior point method on these LPs that we can afford to sacrifice the hotstart facility of the simplex method that would normally dictate its superiority. IPMs enjoy a polynomial worstcase iteration complexity. Indeed, an LP with variables is solved to optimality in iterations Wright (1997). These methods are generally believed to be wellsuited to the solution of very large scale optimization problems Gondzio (2012). IPMs converge to the optimal solution in very few iterations: they usually need only iterations to reach a solution of the problem with variables. However, a single iteration might be costly. A new variant, called the matrixfree IPM Gondzio (2012), removes this potential drawback. It replaces the exact Newton method with the inexact one and employs an iterative method based on conjugate gradients to solve the underlying systems of linear equations. In the next section we describe briefly the matrixfree interior point method and discuss in detail how it is applied to LPs which model the critical visibility.
5 New approach
We start this section with a discussion of a general purpose interior point method for linear programming and its special variant called the matrixfree interior point method. The reader familiar with IPMs may skip this part. The reader interested in more detail on the theory and implementation of IPMs should consult the excellent book of Wright Wright (1997) and the recent survey by Gondzio Gondzio (2012), respectively. The matrixfree IPM itself is introduced by Gondzio Gondzio (2012).
5.1 Matrixfree interior point method
Consider the following linear optimization problem in general form
(12) 
where and . An IPM approaches its solution by moving through the interior of the orthant defined by the nonnegativity constraints. This is achieved by removing the nonnegativity constraints and adding logarithmic barrier terms to the objective function to give the problem
(13) 
The parameter controls the force of the logarithmic barrier so influences the distance of iterates from the boundary of the positive orthant. This parameter is gradually reduced as the algorithm progresses, in order to allow it to get arbitrarily close to the optimum which, for LP problems, is always found on the boundary of the orthant Wright (1997). The first order optimality conditions for (13) are nonlinear equations of the form
(14)  
Here and are the dual variables (Lagrange multipliers) associated with the linear equality constraint and the linear inequalities , respectively, is the vector of ones and () is a diagonal matrix in with elements of the vector () along the diagonal. In interior point methods the system of equations (5.1) is solved using the Newton method Gondzio (2012). One IPM iteration calculates the Newton direction, making one step in this direction and reduces the parameter . The Newton direction is obtained by solving
(15) 
where is the identity matrix of size and represents the reduction of the parameter in every iteration of the algorithm. We can transform (15) by eliminating to get the following symmetric but indefinite augmented system
(16) 
where . Further elimination of gives the following system of normal equations whose matrix of coefficients is symmetric and positive definite:
(17) 
The componentwise products are driven to zero (recall the 3rd equation in (5.1)) and ultimately define an optimal splitting of indices into two subsets: for which and and for which and . Consequently, the diagonal scaling matrix is very illconditioned: when approaches zero, the elements go to infinity or zero for indices or , respectively. The matrices of linear systems (16) and (17) are therefore both very illconditioned. Different regularization techniques have been proposed to cure this illconditioning. Following Saunders Saunders (1996), Altman and Gondzio Altman and Gondzio (1999) propose regularizing both the primal and the dual formulation of the problem, replacing (16) by
(18) 
Here and are the dynamically chosen primal and dual positive definite diagonal regularization matrices and and are appropriately computed right hand side vectors. The regularized version of the normal equations is obtained by pivoting on the block of (18) to give
(19) 
The matrixfree interior point method Gondzio (2012) applies the inexact Newton method to (15). It replaces (17) by the regularized system (19) and applies the preconditioned conjugate gradients algorithm to solve (19) only approximately, and to rather loose accuracy. The resulting solution is used to compute the corresponding and , but clearly the complete Newton direction is only an approximate solution to (15). It is an inexact Newton direction for the first order optimality conditions (5.1).
To accelerate the convergence of the conjugate gradients method the system is preconditioned by a matrix to reduce its condition number so that
(20) 
For matrixfree IPM, a specially designed preconditioner Gondzio (2012) is identified as follows. Consider the rank partial Cholesky decomposition
(21) 
where is the Schur complement obtained by eliminating pivots, is a trapezoidal matrix containing the first columns of the Cholesky factor of and is a diagonal matrix containing the largest pivots of . The preconditioner
(22) 
is (21) with the Schur complement matrix replaced by its diagonal .
The rank of the partial Cholesky decomposition is carefully chosen to obtain an optimal tradeoff between precision and execution time. In general the approximation becomes more precise as the rank increases, but greater effort has to go into its computation and use. When the Cholesky decomposition is complete: there is no Schur complement matrix and so the conjugate gradients method terminates in one iteration. However, for reasons of efficiency, we always set . The choice of rank of the partial Cholesky used in the LP problems considered in this paper will be discussed in the following sections.
Another important feature of the matrixfree interior point method is its ability to allow for an implicit treatment of matrix Gondzio (2012). This is relevant when is very large and therefore requires a lot of memory to store. This is clearly the case for LPs arising in quantum information problems. The matrixfree IPM needs only the following results of operations performed with matrix .

Multiply the matrix and vector .

Multiply the matrix and vector .

Retrieve an arbitrary column of the matrix .

Retrieve the diagonal of the matrix .
Therefore, neither matrix nor matrix needs to be fully stored in order to solve (19). All that is needed is the implementation of the above operations. In the following sections we will discuss the implications of this feature and we will illustrate the practical behaviour of the method both in terms of the memory and CPU time requirements.
5.2 Implementation
As described above the matrixfree interior point method requires implementations of a number of operations which define the LP. Here is a full summary:

f_a_times_x_: the product of the matrix and vector , which is a parameter,

f_at_times_y_: the product of the matrix and vector , which is a parameter,

f_get_col_aat_: the th column of the matrix , where is a parameter,

f_get_diag_aat_: the diagonal of the matrix .
We implemented them all in a program integrating steamroller2 and the matrixfree interior point method of HOPDM software Gondzio (2012).
Performance optimization
The above functions are called very frequently during program execution and are responsible for a significant percentage of the CPU time used. Hence, it is essential that they are implemented efficiently. Our implementations are based on a rowwise storage of the matrix . For each row only the column indices of nonzeros are stored to save space (with the exception of the last column of the matrix, for which real values in the interval also need to be stored). This turned out to be a very efficient method of storage, with a 64k problem consuming less than 10 MB of memory. Indeed, after the redundancy of the problem has been removed, the memory required to store the problem grows sublinearly with its dimension. There are two reasons for this:

The vast majority of nonzero entries are equal to 1.

The percentage of redundant constraints increases with the problem size.
Our first implementations contained solely the definition of how each of the operations should be performed in the C programming language (without any optimization) and showed how much room for optimization there was. All results in this section are obtained from a development workstation and are, therefore, nonindicative in terms of absolute times and a production code would be much faster (see next section for production code results). Table 1 contains output from the UNIX operating system’s profiler gprof when solving a single 16k problem. After reduction of redundancy we have and , the latter being plus one variable representing the visibility. We report the number of calls of each routine, the CPU time (in seconds) per 1000 calls, and the percentage of overall solution time spent while executing a given routine.
operation  calls  time / 1000 calls  % of total time 

f_a_times_x_  12306  1.255  13.22 
f_at_times_y_  11307  3.766  36.46 
f_get_col_aat_  1000  0  0 
f_get_diag_aat_  10  0  0 
Only the costs of f_a_times_x_ and f_at_times_y_ are significant, accounting for about half of the execution time, so the following optimizations were performed to improve their efficiency.
Optimization of f_a_times_x_. Since all columns of apart from the
last contain either ones or zeros, the contributions of these entries
in a summation are either skipped (if zero) or added (if one)
when the intermediary sum builds a particular component of the
resulting vector. Thus the multiplication by one that would occur in
general is skipped. For sufficiently large problems (larger than
256) loop unrolling is performed (following an established
practice described, for example, by Sarkar and
Vivek Sarkar (2001) and commonly used in, for example,
LAPACK Anderson et al. (1992)). Each row contains the same number of
nonzero elements. These are divided into groups of 16 and added to
the intermediate sum all at once. This reduces the loop control
overhead and yields an almost twofold speedup of the f_a_times_x_ routine
(see Table 2).
operation  calls  time / 1000 calls  % of total time 

f_a_times_x_  11943  0.677  9.07 
f_at_times_y_  10944  2.404  29.50 
Optimization of f_at_times_y_. As for f_a_times_x_,
multiplication by one is skipped and loop unrolling is performed.
Here it is somewhat less efficient due to the rowwise storage
of matrix (see Table 2). The intermediate sums
have to be stored in the resulting dimensional vector itself and, despite
loop unrolling, they are frequently accessed at each step of the algorithm.
Further speedup was achieved when columnwise storage of matrix
was used. This enabled us to use intermediary sums and store them only once
in the resulting vector, just as in the case of f_a_times_x_.
As mentioned at the beginning of this section, storage is very efficient
so we can safely afford to store a duplicate copy of the problem matrix.
The implementation of loop unrolling for the columnwise
storage was more difficult due to the different number of nonzero
elements in columns, with even very large problems containing multiple
columns with one, three or five elements. After a considerable effort,
resulting in more than a hundred lines of loop unrolling code, we achieved
the significant speedup seen in Table 3.
The total resulting speedup of f_at_times_y_ after optimizations
is more than threefold compared to unoptimized version of the function.
operation  calls  time / 1000 calls  % of total time 

f_a_times_x_  11697  0.694  11.19 
f_at_times_y_  10698  0.955  14.08 
We do not provide total execution times as they are not informative. This is due to the order of additions and the associated numerical errors, causing the optimized and unoptimized implementations to differ in the last significant digits and, hence, gradually alter the path to optimality taken by the matrixfree interior point method. This is also the reason why the optimized and the unoptimized implementations have different numbers of function calls. Having run the program many times for different input we noticed a significant variation in the number of calls between the optimized and unoptimized functions and it is safe to assume that the difference does not favour either implementation. Neither of the two versions is numerically superior to the other in terms of precision and they would behave identically if no numerical errors were present.
Significant overall speedup of execution in the range 25–75% was observed after implementing the above improvements.
Precision
Obtaining sufficient precision was the major challenge of our investigations. It depends on a number of parameter settings and we have made a considerable effort to choose these parameters optimally. In the results reported below we analyse problems with so there are two possible outcomes for any measurement ( or ), 2 observables per each observer and increasing numbers of observers. In this setting the problem size grows fourfold each time an observer is added to the problem. So the minimal problem with 2 observers is 16, the problem with 3 observers is 64. However, after eliminating redundancy, the size of LP increases threefold with each additional observer.
The rank of partial Cholesky determines precision since, when it is large, it allows for more precise direction calculations. It cannot, however, be increased indefinitely lest the execution time exceed that of the simplexbased approach. We found a value of 100 to be optimal for our computations and we used it for all problem sizes. As a result, for LP problems of size 100 or less the preconditioner is exact so conjugate gradients is exact in one iteration and the interior point method uses the exact Newton direction. Since problems 256 and 1k (with 4 and 5 observers respectively) yield LPs of size 81 and 243 respectively, we only analyse problems with at least 5 observers since only for these problems is conjugate gradients nontrivial.
The matrixfree interior point method applies the preconditioned conjugate gradients algorithm to the system (19) and uses the partial Cholesky preconditioner (22). The solution of (19) is controlled by two parameters: the maximum number of conjugate gradient iterations allowed and the desired precision achieved by the algorithm. The use of high accuracy and/or setting a large conjugate gradients iterations limit might increase the execution time but these settings are needed to provide the accuracy required by (inexact) Newton directions. We vary these settings and start from a lower iteration limit at the first steps of IPM, gradually increasing it to 1000 towards the end of optimization. We also set the relative error tolerance in the CG algorithm to .
Another parameter which influences the precision of the method is the optimality tolerance used to terminate the IPM when
(23) 
where are defined in Section 5.1.
We set for all problems.
With the above settings we managed to obtain twodigit exact optimal solutions for all problems up to 16k that were attempted. We also achieved a precision of in the solutions for larger problems. Although these solutions are less accurate than those published in Gruca et al. (2010) they are sufficiently accurate to elicit the qualitatively important properties of most states. In practice the results are typically much more precise than this worstcase value, as may be seen in the tables in Section 6.
6 Performance comparison
We begin our comparison with individual LPs solved by the two approaches: the simplex method (GLPK) and the matrixfree interior point method (HOPDM). Results in this section are based on runs on a computer with an Intel^{®} Core™ i7 3.07GHz processor and 24 GB of RAM. Computational tests are performed using a Bell experiment with 2 observables per observer since the minimal critical visibility of the GHZ state is known and given in Mermin (1990) as
(24) 
where is the number of observers.
Table 4 reports the worstcase numerical error for results obtained with the matrixfree interior point method during an execution of the program which stopped after one hundred critical visibility function evaluations. This means that one hundred LPs were generated and solved successively. The angles used to generate LPs were provided by the higherlevel optimization method (DSM), as in a typical execution of the program. For a given number of observers, each line in Table 4 reports different execution results, one for each problem size. The reduced problem size corresponds to the number of rows in the constraint matrix after the redundancy is removed. For all problems solved with both solvers (GLPK and HOPDM) the solution results have been compared. The last column in the table contains the largest difference between the objective returned by GLPK (which have been confirmed repeatedly to be correct and precise) and by matrixfree HOPDM.
observers  problem name  reduced problem size  the largest difference 

5  1k  243  0.006535 
6  4k  729  0.006738 
7  16k  2187  0.003662 
In Table 5 we report the time required by the two solvers, the expected solution and the optimal objective function value (critical visibility) returned by HOPDM. The value returned by GLPK always matched the expected solution obtained analytically so it is skipped from the tables. Values calculated with (24) are given in the solution column of the table. The angles used to generate each LP were the same for both GLPK and HOPDM. They were also minimal for each problem for the GHZ state. Each line in Table 5 corresponds to one LP. Problems 64k and larger were scaled so that the returned value was in fact 8 times greater. This helped preserve the precision of the results, where the value of the solution is close to zero. As one can see from studying the table some of them are still off by a small fraction, which is the expected speedaccuracy tradeoff.
problem 




name  size  solution  time  time  solution  
4  256  81  0.354  0s  0s  0.354  
5  1k  243  0.250  0.02s  0.08s  0.250  
6  4k  729  0.177  0.93s  0.87s  0.177  
7  16k  2187  0.125  1m13s  11.88s  0.125  
8  64k  6561  0.088  6h51m  3m22s  0.088  
9  256k  19683  0.063  24h  28m38s  0.068  
10  1m  59049  0.044  unknown  1h34m  0.051  
11  4m  177147  0.031  unknown  9h15m  0.029 
The analysis of results collected in Table 5 reveals that the simplex method is faster for smaller problems but its solution time grows very quickly with problem size. For problem 16k (of size 2187) GLPK is already an order of magnitude slower than HOPDM.
Computations using GLPK for problem 256k were terminated after about 24 hours and computations for problem 1m were not even attempted as they cannot be expected to complete within any reasonable time. HOPDM times, on the other hand, increase roughly linearly with problem size, with successive factors being 10.9, 13.66, 17.00, 8.50, 3.29, 5.88. This confirms that there is no explosion of algorithm steps to take when solving the LP as the problem size increases. One can expect that the computation time will increase linearly with reduced problem size (and depend very little on the number of variables in the problem).
It should be noted, however, that the simplex method has a significant advantage over the interior point method when a sequence of similar LPs are solved. Successive LPs differ merely by the last column in the matrix which contains the angles and therefore one may expect that the solutions of two successive problems are similar. The simplex method in GLPK has an option which allows the optimal solution of one problem to be used to warmstart another and this technique proved very useful. Although, in general, it is possible to use the warmstarting facility of IPMs, this option is not available for the matrixfree variant of IPM. Therefore, another comparison which reflects this relative advantage of the simplex method is needed. The results collected in Table 6 demonstrate the performance of the complete approach which uses DSM to solve the upper level optimization problem and employs either the simplex method (GLPK) or the matrixfree interior point method (HOPDM) to solve the LPs at the lower level. For both variants of the program we report the number of LPs solved, the overall time of running the program and the optimal function value. We report only two significant digits of the result which corresponds to the precision requested from the DSM procedure.
problem  Simplex Method  Interior Point Method  

name  size  solution  # LPs  time  # LPs  time  solution  
2  16  9  0.71  44  0s  44  0.05s  0.71 
3  64  27  0.50  78  0s  78  0.13s  0.50 
4  256  81  0.35  156  0.07s  156  0.90s  0.35 
5  1k  243  0.25  231  1.67s  161  21s  0.25 
6  4k  729  0.18  251  56s  265  5m28s  0.18 
7  16k  2187  0.13  346  2h00m  330  56m52s  0.13 
8  64k  6561  0.088  unknown  24h  331  16h47m  0.090 
In both cases the same set of angles (drawn randomly beforehand) were provided to DSM as a starting point for the optimization process. Since GLPK and HOPDM return slightly different results for a single LP, the execution of DSM differs slightly (as one can deduce from the different number of LPs generated). Both, however, return the same correct result to within twodigits, except for problem 64k.
For the three smallest problems the number of LPs solved by both approaches are the same. This is because the number of constraints for first three problems in our test set (9, 27 and 81) are less than the rank of the partial Cholesky preconditioner so the latter is the complete Cholesky decomposition of . This perfect preconditioner allows the matrixfree IPM to produce very accurate optimal solutions. Analysis of results collected in Table 6 confirms an advantage of the simplex method approach over the matrixfree interior point method for smaller problems. However, it also indicates that for problem 16k the latter approach provides better efficiency.
The largest problem for which it was possible to obtain the solution with the simplex method was 16k. Larger problems are intractable by GLPK due to time constraints. It was possible to obtain a solution to a problem 4 times larger, that is problem of size 64k, with the IPM method. In order to obtain a reliable solution we had to make two adjustments to the computation method. Firstly, we had to relax the sensitivity of the DSM procedure to numerical errors of the optimized critical visibility function. This adjustment was necessary as the cumulative precision error introduced by inaccurate solutions to LPs misguides the DSM procedure and may eventually hamper its convergence. We set the parameter in our DSM implementation to compared to a previous setting of . is the fractional convergence tolerance, which we have found to return results correct on the second decimal place for the standard setting of used for all results in this paper so far. After relaxing this parameter by setting it to we can expect a less precise result, that is one potentially above the minimum, but also easier convergence. For details see Press et al. (2007). Having made this one adjustment we obtained a result , compared to the correct result of . We obtained a further improvement of the computational result by increasing the size of the partial Cholesky preconditioner rank from to . The result of obtained this way is reported in Table 6. This last problem raises a question of whether such a result, which is not precisely equal to the analytical predictions, is of use to the quantum information researcher. There are two reasons why this is so. The first one is that we can have confidence in this result. The inaccuracy comes from the lowerlevel optimization procedure feeding results to the higher level optimization task. Despite this inaccuracy the higher level optimization procedure did converge, which means that the result fulfilled all necessary (relaxed) convergence conditions. This indicates that the inaccuracy was kept within certain bounds throughout the optimization process. We can also have confidence in the result as due to the nature of the optimization process this final result is vitiated by the error of only the last LP computation. This means that the errors do not propagate in between the subsequent calls of the higher level optimization procedure. It is possible that, due to inaccuracy of the interior point method, the DSM procedure will fail to converge to a global minimum (but then again there is no guarantee it will converge to a global minimum in any case McKinnon (1998)). It is not possible, however, that due to the IPM inaccuracy DSM will converge to a point below the global minimum. In other words, the lowest result we can get will occur in a situation where we got a global minimum (point defined by angles), and for this point we got a critical visibility vitiated by the error of only the final LP calculation. For a researcher in the quantum information domain such a result can be harnessed in the following way. First, one would execute the full DSM optimization with the IPM method to arrive quickly at the global minimum (or somewhere nearby). Second, one would execute a single LP calculation using this point and a more precise LP method like the simplex method to confirm the result. If this is not possible for some reason, then at least the following conclusion can be made. The result obtained might be arbitrarily higher than the desired minimal critical visibility, depending on the convergence of the DSM procedure. But if it is lower, then it will not be by more than the inaccuracy of the single last LP calculation made. That is, the preceding LP calculation errors do not affect the final result in any way other than hampering the convergence of DSM. This allows the quantum information researcher to compare the nonclassicality of quantum entangled states to within the error bound introduced solely by the final LP computation. The above example suggests that in the case of problem size 64k the researcher can compare two states to within 0.002, which in many cases allows to elicit qualitative differences between states. It is worth adding that to protect ourselves from the risk of not getting the global minimum as a result of the DSM computation, we repeat it a number of times, be it with the simplex method or the IPM method as the underlying LP solver.
Let us highlight two other important points related to Table 6. Firstly, as would be clear to any researcher in the quantum information domain, there is a world of investigation opportunities available for problems of sizes between 16k and 64k (and potentially also slightly larger). It is possible to create problems of sizes in these bounds by adjusting the number of observers in an experiment (for simplicity this was the only variable in this section of this paper), the number of outcomes of any single measurement (which we assumed to be always 2), and the number of observables available for selection by any of the observers. The last of these was assumed to be always 2, whereas in practice not only can it vary for different experiments — it can also vary between different observers in the same experiment. With so many variables one can generate problems of almost any size in the range 16k – 64k. The majority of these problems were beyond computational reach before the introduction of the interior point method in the calculations. All of them can now be solved by HOPDM by making the above adjustments.
Secondly, the interior point method returns reasonably precise results for any single LP, but the inconsistencies involved are random leading to DSM convergence problems. Thus, global optimization does not demonstrate the biggest benefit of using HOPDM. The real power of the interior point method lies in single LP calculations, where the slight numerical error doesn’t impact the general contribution of the result. Deducing potential minimal angles is possible and frequently applied Mermin (1990). Therefore, even when full DSM minimization is not possible as in the case of problems larger than 64k, single LP calculations are still of paramount value for researchers in the field as they are a way to calculate the critical visibility for the deduced minimal angles. There is a substantial range of problems made available for numerical investigations thanks to the interior point method — see Table 5 and the next section for details.
7 Comparison of range of problems within reach
Finally let us present a broader range of singleLP problems that it is possible to solve with the simplex method and the interior point method. Table 7 and Figure 1 summarize this comparison. For each count of observers we solved the minimal problem and steadily increased the number of observables available to each of the observers by one (that is, for 2 observers we began with a 2x2 observable configuration, then tried 3x3, 4x4 and so on). We continued this process until we encountered a problem which either used too much RAM or exceed a 24 hour computation time threshold. We selected 24 hours as a threshold as we found that this is a maximal period which allows comfortable trial and error work with different states (used e.g. in Gruca et al. (2010)).
Simplex Method  Interior Point Method  






2  11x11  12x12  22m  13x13  
3  7x7x7  8x8x8  1h15m  9x9x9  
4  5x5x5x5  6x6x6x6  22h11m    
5  3x3x3x3x3  4x4x4x4x4  2h00m  5x5x5x5x5  
6  2x2x2x2x2x2  3x3x3x3x3x3  41m  4x4x4x4x4x4  
7  2x2x2x2x2x2x2  3x3x3x3x3x3x3  9h23m    
8  2x2x2x2x2x2x2x2  2x2x2x2x2x2x2x2  3m  3x3x3x3x3x3x3x3  
9    2x2x2x2x2x2x2x2x2  28m  3x3x3x3x3x3x3x3x3  
10    2x2x2x2x2x2x2x2x2x2  1h34m  3x3x3x3x3x3x3x3x3x3  
11    2x2x2x2x2x2x2x2x2x2x2  9h15m   
As one can see in both Table 7 and Figure 1, for each count of observers experiments with greater number of observables are within reach of the interior point method. By “within reach” we mean that it is possible to compute them on our production machine in 24 hours or less. Also, a greater number of observers is within reach of the interior point method than of the simplex method in the case of the minimal number 2 of observables. It should also be noted that, as with all other results presented in this section, we were limited to the 24 GB of RAM available on the machine we used for the computations. Thus, a number of problems can be expected to complete in a reasonable time, but it was not possible to verify it due to limited RAM. We indicated such problems in both the table and the figure. Our rule of thumb for distinguishing problems potentially possible to solve from the ones definitely beyond our reach, even with enough RAM, is very approximate and roughly based on the sequence of timeincrease factors (see Table 5). We assumed that if the preceding problem can be solved in about 2 hours or less, then the next problem could potentially be solved in 24 hours or less. Our estimation is obviously not precise and the results in the rightmost column of the table should be treated as directions of future investigation, when machines with more RAM are available.
8 Conclusions
Minimal critical visibility of a quantum entangled state is a very useful measure in quantum information. However its computation poses a substantial numerical challenge. At the heart of this challenge lie the very large linear programming problems, the solution of which is necessary each time the critical visibility itself is calculated. Previous work always used the simplex method to solve these LPs, an approach providing satisfactory precision in the computed results, but restricting researchers to a limited range of problems which can be solved. We proposed to replace the simplex method with an advanced matrixfree interior point method HOPDM and found that it enabled us to analyse problems much faster and extend our investigations to more complicated ones, which were not previously within our reach. We found that even though the results provided by the interior point method were less precise, they sufficed to analyse, order and compare quantum entangled states, even in complicated cases.
9 Acknowledgements
This work was carried out under the HPCEuropa2 project (project number: 228398) with the support of the European Commission Capacities AreaResearch Infrastructures Initiative. This work is also available as Technical Report ERGO2012004 (for other papers in this series see http://www.maths.ed.ac.uk/ERGO/). Jacek Gondzio was supported by EPSRC grant EP/I017127/1. Wiesław Laskowski and Marek Żukowski are supported by a MNiSW Grant IdP2011 00361 (Ideas Plus).
Footnotes
 journal: Journal of Computational and Applied Mathematics
 Cases where there are more possible measurement outcomes, more observers and more observables available to the observers can also be solved by our approach but are omitted here for simplicity.
 One could define other measures of nonclassicality of the correlations, e.g. see Gruca et al. (2012). However, as our aims are to present the basics behind the computational methods, we shall stick here to the simplest one.
References
 J. S. Bell, On the EinsteinPodolskyRosen paradox, Physics 1 (1964) 195–200.
 D. Kaszlikowski, P. Gnaciński, M. Żukowski, W. Miklaszewski, A. Zeilinger, Violations of local realism by two entangled dimensional systems are stronger than for two qubits, Phys. Rev. Lett. 85 (2000) 4418–4421.
 J. Gruca, W. Laskowski, M. Żukowski, N. Kiesel, W. Wieczorek, C. Schmid, H. Weinfurter, Nonclassicality thresholds for multiqubit states: Numerical analysis, Physical Review A 82 (2010) 012118.
 GNU Linear Programming Kit (GLPK), v4.47, GNU, http://www.gnu.org/software/glpk/, 2011.
 J. Gondzio, Matrixfree interior point method, Computational Optimization and Applications 51 (2012) 457–480.
 A. Fine, Hidden variables, joint probability, and the bell inequalities, Physical Review Letters 48 (1982) 291.
 J. Gruca, W. Laskowski, M. Żukowski, Nonclassicality of pure twoqutrit entangled states, Phys. Rev. A 85 (2012) 022118.
 W. Laskowski, Nonclassical correlations of systems of many qubits, Ph.D. thesis, Department of Mathematics, Physics and Informatics, University of Gdansk, Gdansk, Poland, 2007.
 M. Nielsen, I. Chuang, Quantum Computation and Quantum Information, Cambridge University Press, Cambridge, United Kingdom, 2000.
 D. M. Greenberger, M. A. Horne, A. Zeilinger, Bell’s Theorem, Quantum Theory, and Conceptions of the Universe, Kluwer Academic, Dordrecht, The Netherlands, p. 69, 1989.
 V. Scarani, N. Gisin, Spectral decomposition of Bell’s operators for qubits, Journal of Physics A: Mathematical and General 34 (2001) 6043.
 G. B. Dantzig, Linear Programming and Extensions, Princeton University Press, Princeton, 1963.
 J. A. Nelder, R. Mead, A simplex method for function minimization, Computer Journal 7 (1965) 308–313.
 E. Jones, T. Oliphant, P. Peterson, et al., SciPy: Open source scientific tools for Python, 2001–2011.
 J. Gruca, Nonclassical properties of multiqubit quantum states applicable in quantum information science, Master’s thesis, Department of Mathematics, Physics and Informatics, University of Gdansk, Gdansk, Poland, 2008.
 J. A. J. Hall, K. I. M. McKinnon, The simplest examples where the simplex method cycles and conditions where expand fails to prevent cycling, Mathematical Programming 100 (2004) 133–150.
 K. I. M. McKinnon, Convergence properties of the NelderMead simplex method in low dimensions, SIAM Journal on Optimization 9 (1998) 148–158.
 W. H. Press, S. Teukolsky, W. Vetterling, B. Flannery, Numerical Recipes, The Art of Scientific Computing, Cambridge University Press, Cambridge, United Kingdom, 3rd edition, 2007.
 S. J. Wright, PrimalDual InteriorPoint Methods, SIAM, Philadelphia, 1997.
 J. Gondzio, Interior point methods 25 years later, European Journal of Operational Research 218 (2012) 587–601.
 M. Saunders, Choleskybased methods for sparse least squares: the benefits of regularization, in: L. Adams, L. Nazareth (Eds.), Linear and Nonlinear Conjugate GradientRelated Methods, SIAM, Philadelphia, USA, 1996, pp. 92–100.
 A. Altman, J. Gondzio, Regularized symmetric indefinite systems in interior point methods for linear and quadratic optimization, Optimization Methods and Software 1112 (1999) 275–302.
 V. Sarkar, Optimized unrolling of nested loops, Int. J. Parallel Program. 29 (2001) 545–581.
 E. Anderson, Z. Bai, C. Bischof, J. Demmel, J. Dongarra, J. D. Croz, A. Greenbaum, S. Hammarling, A. McKenney, S. Ostrouchov, D. Sorensen, LAPACK Users’ Guide, Philadelphia, 1992.
 N. D. Mermin, Extreme quantum entanglement in a superposition of macroscopically distinct states, Phys. Rev. Lett. 65 (1990) 1838–1840.