Analysis of the Karmarkar-Karp Differencing Algorithm
aa r X i v : . [ c s . NA ] O c t EPJ manuscript No. (will be inserted by the editor)
Analysis of the Karmarkar-Karp Differencing Algorithm
Stefan Boettcher and Stephan Mertens , Department of Physics, Emory University, Atlanta GA 30322-2430, U.S.A. Institut f¨ur Theoretische Physik, Otto-von-Guericke Universit¨at, PF 4120, 39016 Magdeburg, Germany Santa Fe Institute, 1399 Hyde Park Road, Santa Fe, New Mexico 87501, U.S.A.October 27, 2018
Abstract
The Karmarkar-Karp differencing algorithm is the best known polynomial time heuristic for the numberpartitioning problem, fundamental in both theoretical computer science and statistical physics. We analyze the per-formance of the differencing algorithm on random instances by mapping it to a nonlinear rate equation. Our analysisreveals strong finite size effects that explain why the precise asymptotics of the differencing solution is hard to establishby simulations. The asymptotic series emerging from the rate equation satisfies all known bounds on the Karmarkar-Karp algorithm and projects a scaling n − c ln n , where c = / ( ) = . ... . Our calculations reveal subtle relationsbetween the algorithm and Fibonacci-like sequences, and we establish an explicit identity to that effect. PACS.
Consider a list of n positive numbers. Replacing the two largestnumbers by their difference yields a new list of n − n − differencing , and the procedure that itera-tively selects the two largest numbers for differencing is knownas largest differencing method or LDM. This method was intro-duced in 1982 by Karmarkar and Karp [1] as an algorithm forsolving the number partitioning problem (N PP ): Given a list a , a , . . . , a n of positive numbers, find a partition, i.e. a subset A ⊂ { , . . . , n } such that the discrepancy D ( A ) = (cid:12)(cid:12)(cid:12) (cid:229) i ∈ A a i − (cid:229) i A a i (cid:12)(cid:12)(cid:12) , (1)is minimized. Obviously, LDM amounts to deciding iterativelythat the two largest numbers will be put on different sides ofthe partition, but to defer the decision on what side to put eachnumber. The final number then represents the discrepancy.Despite its simple definition, the N PP is of considerableimportance both in theoretical computer science and statisti-cal physics. The N PP is NP-hard, which means (a) that no al-gorithm is known that is essentially faster than exhaustivelysearching through all 2 n partitions, and (b) that the N PP is com-putationally equivalent to many famous problems like the Trav-eling Salesman Problem or the Satisfiability Problem [2]. Infact, the N PP is one of Garey and Johnson’s six basic NP-hard problems that lie at the heart of the theory of NP-completeness[3], and it is the only one of these problems that actually dealswith numbers. Hence it is often chosen as a base for NP-hard-ness proofs of other problems involving numbers, like bin pack-ing, multiprocessor scheduling [4], quadratic programming orknapsack problems. The N PP was also the base of one of thefirst public key crypto systems [5].In statistical physics, the significance of the N PP resultsfrom the fact that it was the first system for which the localREM scenario was established [6, 7]. The notion local REMscenario refers to systems which locally (on the energy scale)behaves like Derrida’s random energy model [8, 9]. It is con-jectured to be a universal feature of random, discrete systems[10]. Recently, this conjecture has been proven for several spinglass models [11, 12] and for directed polymers in random me-dia [13].Considering the NP-hardness of the problem it is no sur-prise that LDM (which runs in polynomial time) will generallynot find the optimal solution but an approximation. Our initialquestion asks for the quality of the LDM solution to N PP , andto address this question we will focus on random instances ofthe N PP where the numbers a j are independent, identically dis-tributed (i.i.d.) random numbers, uniformly distributed in theunit interval. Let L n denote the output of LDM on such a list.Yakir [14] proved that the expectation E [ L n ] is asymptoticallybounded by n − b ln n ≤ E [ L n ] ≤ n − a ln n , (2)where a and b are (unknown) constants such that b ≥ a ≥
12 ln 2 = . . . . . (3)In this contribution we will argue that b = a = . Stefan Boettcher, Stephan Mertens: Analysis of the Karmarkar-Karp Differencing Algorithm
The paper is organized as follows. We start with a compre-hensive description of the differencing algorithm, a simple (butwrong) argument that yields the scaling (2) and a presentationof simulation data that seems to violate the asymptotic bound(3). In section 3 we reformulate LDM in terms of a stochasticrecursion on parameters of exponential variates. This recursionwill then be simplified to a deterministic, nonlinear rate equa-tion in section 4. A numerical investigation of this rate equationreveals a structure in the dynamics of LDM that can be used asan Ansatz to simplify both the exact recursions and the rateequation. This will lead to a simple, Fibonacci like recursion(section 5) and to an analytic solution of the rate equation (sec-tion 6). In both cases we can derive the asymptotics includingthe corrections to scaling, and we claim that a similar asymp-totic expansion holds for the original LDM. The latter claim iscorroborated by fitting the asymptotic expansion to the avail-able numerical data on LDM. (c) (d)(a) (b)15 6 783 12 14 4
Figure 1.
The differencing algorithm in action.
The differencing scheme as described in the introductiongives the value of the discrepancy, but not the actual partition.For that we need some additional bookkeeping, which is mosteasily implemented in terms of graphs (Fig. 1). The algorithmmaintains a list of rooted trees where each root is labeled witha number. The algorithm starts with n trees of size one and theroots labeled with the numbers a i . Then the following steps areiterated until a single rooted tree of size n remains:1. Among all roots, find those with the largest ( x ) and secondlargest ( y ) label.2. Join nodes x and y with an edge, declare node x as the rootof the new tree and relabel it with x − y .After n − ( , , , , ) .The final two coloring corresponds to the partition ( , , ) ver-sus ( , ) with discrepancy 2. Note that the optimum partition ( , , ) versus ( , ) achieves discrepancy 0.Technically, LDM boils down to deleting items from andinserting items into a sorted list of size n . This can be donein time O ( n ln n ) using an advanced data structure like a heap[15]. Hence LDM is very efficient, but how good is it? Aswe have already seen in the example, LDM can miss the op-timal partition. And for random instances, the corridor (2) is farabove the true optimum, which is known to scale like Q ( √ n − n ) [7]. Yet LDM yields the best results that can be achieved inpolynomial time. Many alternative algorithms have been inves-tigated in the past [16, 17], but they all produce results worsethan (2). The few algorithms that can actually compete with theKarmarkar-Karp procedure use the same elementary differenc-ing operation [18, 19]. It seems as if the differencing schememarks an inherent barrier for polynomial time algorithms.The following argument explains the scaling (2). The typ-ical distance between adjacent pairs of the n numbers in theinterval [ , ] is n − . Hence after n / n / [ , n − ] . The typi-cal distance between pairs is now 2 n − . After another round of n / n / [ , n − ] . In general, after 2 k differencing operations we areleft with n / k numbers in the range [ , ( k ) n − k ] . Reducing theoriginal list to a single number requires k = log n differencingoperations, and applying the above argument all the way downsuggests that E [ L n ] (cid:181) n − c ln n (4)with c =
12 ln 2 = . . . . . (5)As we will see, this is the right scaling, yet the argument can-not be correct. This follows from the fact that it predicts thesame scaling for the paired differencing method (PDM). Herein each round all pairs of adjacent numbers are replaced by theirdifference in parallel. This method, however, yields an averagediscrepancy of order Q ( n − ) [20]. Yet, our analysis below sug-gests that (4) and (5) indeed describe the asymptotic behaviorcorrectly, although a far more subtle treatment is required.An obvious approach to find the quality of LDM are sim-ulations. We ran LDM on random instances of varying size n ,and Figure 2 shows the results for E [ L n ] . Apparently ln E [ L n ] scales like ln n , in agreement with (2) and (4). A linear fitseems to yield c ≃ . c ≥ / n = is too small to see the trueasymptotic behavior. This may be the reason why Monte Carlostudies of LDM never have been published.A plot of the probability density function (pdf) of L n / E [ L n ] reveals a data collapse varying values of n (Fig. 3). Apparentlythe complete statistics of L n is asymptotically dominated by asingle scale n − c ln n .Some technical notes about simulating LDM are appro-priate. Differencing means subtracting numbers over and overagain. The numerical precision must be adjusted carefully to tefan Boettcher, Stephan Mertens: Analysis of the Karmarkar-Karp Differencing Algorithm 3
50 100 150ln n - l n E [ L n ] Figure 2.
Results of LDM applied to n random i.i.d. numbers, uni-formly drawn from the unit interval. Each data point represents be-tween 10 (large n ) and 10 samples (small n ). The solid line is thelinear fit − ln E [ L n ] = . + .
65 ln n . L n / E[L n ] -3 -2 -1 n = 10 n = 10 n = 10 n = 10 Figure 3.
Probability density function of L n / E [ L n ] . support this and to be able to represent the final discrepancyof order n − c ln n . We used the freely available GMP library [21]for the required multiple precision arithmetic and ran all simu-lations on ℓ -bit integers where the number of bits ranges from ℓ =
40 (for n =
20) to ℓ =
300 for n = . · . The integerdiscrepancies were then rescaled by 2 − ℓ . The pseudo randomnumber generator was taken from the TRNG library [22]. A common problem in the average-case analysis of algorithmslike LDM is that numbers become conditioned and cease to beindependent as the algorithm proceeds. Lueker [20] proposedto use exponential instead of uniform variates to cope with thisproblem. Let X , . . . , X n + be i.i.d. random exponentials withmean 1 and consider the partial sums S k = (cid:229) ki = X i . Then the joint distribution of the ratios S k / S n + , k = , . . . , n , is the sameas that of the order statistics of n i.i.d. uniform variates from [ , ] [23]. As a consequence, LDM will produce the same dis-tribution of data no matter whether it is run on uniform variatesor on S k / S n + . Let ˆ L n denote the result of LDM on the par-tial sums S , S , . . . , S n . Since the output of LDM is linear in itsinput, we have ˆ L n D = S n + L n , (6)where S n + is the sum of n + X D = Y indicates that the random variable X and Y have the same distribution. The probability density of S n + isthe gamma density g n + ( s ) = s n n ! e − s . (7)Taking expectations of both sides of (6) we getE [ L n ] = E (cid:2) ˆ L n (cid:3) n + . (8)This allows us to derive the asymptotics of E [ L n ] from theasymptotics of E (cid:2) ˆ L n (cid:3) .Exponential variates are well suited for the analysis of LDMbecause the sum and difference of two exponential variates areagain exponential variates. Once started on exponential vari-ates, LDM keeps working on exponentials all the time. Thisallows us to express the operation of LDM in terms of a recur-sive equation for the parameters of exponential densities [14].We start with the following Lemma: Lemma 1
Let X and X be independent exponential randomvariables with parameter l and l , resp.. The probability ofthe event X < X is given by P ( X < X ) = l l + l . (9) Furthermore, conditioned on the event X < X , the variablesX and X − X are independent exponentials with parameters l + l (for X ) and l for X − X . The proof of Lemma 1 consists of trivial integrations of theexponential densities and is omitted here.Next we consider generalized partial sums of exponentials,described by n -tuples ( l , l , . . . , l n ) . This n -tuple is shorthand for the sequence of partial sums ( X , X + X , . . . , n (cid:229) i = X i ) with X i = EXP ( l i ) .Now let us look at the result of one iteration of LDM on ( l , l , . . . , l n ) . The two largest numbers are removed and re-placed by their difference X n which is an EXP ( l n ) variate.Lemma 1 tells us, that the probability that this number is thesmallest in the list isP ( X n < X ) = l n l + l n , Stefan Boettcher, Stephan Mertens: Analysis of the Karmarkar-Karp Differencing Algorithm and conditioned on that event, the smallest number is an EXP ( l + l n ) variate and the increment to the 2nd smallest number X − X n is an independent EXP ( l ) variate. Conditioned on X n < X we get another l -tuple as the input for the next iteration: X n < X ⇒ ( l + l n , l , l , . . . , l n − ) The probability that X n ≥ X isP ( X n ≥ X ) = l l + l n , and in this case X is an EXP ( l n + l ) variate, whereas thedifference X n − X is an EXP ( l n ) variate. Now the probabilitythat the new number X n is second in the new list readsP ( X n ≥ X ∩ X n < X + X ) = P ( X n ≥ X ∩ X n − X < X )= l l + l n l n l + l n and conditioned on that event the input for the new iteration is ( l + l n , l + l n , l , . . . , l n − ) . This argument can be iterated to calculate the probability of X n becoming the k -th number in the new list. Denoting the partialsums by S k we getP ( X n ≥ S k − ∩ X n < S k ) = l n l k + l n k − (cid:213) i = l i l i + l n (10)for k = , . . . , n − ( l + l n , . . . , l k + l n , l k , l k + , . . . , l n − ) . (11)The final case is that X n becomes the largest number in the newlist. This happens with probabilityP ( X n ≥ S n − ) = n − (cid:213) i = l i l i + l n (12)and leads to the list ( l + l n , . . . , l n − + l n , l n ) . (13)In all cases we stay within the set of instances given by partialsums of independent exponentials, and we can apply Eqs. (10)to (13) recursively until we have reduced the original problemto a ( l , l ) -instance which tells us that the final difference isan EXP ( l ) variate.Fig. 4 shows the result of this analysis on the input ( , , , ) ,our original problem with n =
4. We have to explore the treethat branches according to the position that is taken by thenew number inserted in the shortened list. The numbers writ-ten on the edges of the tree are the probabilities for the cor-responding transition. Note that we have combined the twobranches emerging from the root that both lead to a ( , , ) -configuration into a single one by adding their probabilities. Inthe end we get p ( x ) =
23 e − x +
23 e − x (1 , , , , ,
1) (2 , , ,
2) (3 , ,
2) (3 ,
12 1213 23 13 23
Figure 4.
Statistics of LDM on n =
4. The final difference is dis-tributed according to p ( x ) = · e − x + · − x k \ n
23 1324 41120 49180 4312520
13 16 518 18 5273456
724 772 10734320 307938880
118 53360 1492100
772 4863595443200 Table 1.
Coefficients a ( n ) k in (14). for the probability density function (pdf) of ˆ L . In general, thepdf of ˆ L n is a sum of exponentials, p n ( x ) = (cid:229) k a ( n ) k k e − kx (14)where a ( n ) k is the probability of LDM returning an EXP ( k ) -variate. For small values of n , this probabilities can be calcu-lated by expanding the recursions explicitly (Table 1), but forlarger values of n this approach is prohibited by the exponen-tial growths of the number K ( n ) of branches that have to beexplored.Alternatively we can explore the tree of l -tuples by walk-ing it randomly. Given a tuple ( l . . . , l n ) , we generate a ran-dom integer 1 ≤ k ≤ n − ( k ≤ ℓ ) = ( − (cid:213) ℓ j = l j l j + l n ( ℓ < n − ) ( ℓ = n − ) (15) tefan Boettcher, Stephan Mertens: Analysis of the Karmarkar-Karp Differencing Algorithm 5 l / E[l ] Figure 5.
Probability density function of l / E [ l ] . and using this random k we generate a new tuple of size n − l is the parameter forthe statistics of ˆ L . The probability density of l / E [ l ] is shownin Fig. 5). Again the data collapse corroborates the claim thatthe statistics of LDM is dominated by a single scale.
We can turn the exact recursions from Sec. 3 into a set of rateequations for the time-evolution of the average l -tuple. Let l ti denote the value of l i after t iterations, such that (cid:0) l t , l t , . . . , l tn − t (cid:1) → (cid:0) l t + , l t + , . . . , l t + n − t − (cid:1) . (16)As explained in Sec. 3, at “time” t a number k , 1 ≤ k ≤ n − − t is chosen with probabilityP t ( k ≤ ℓ ) = − (cid:213) ℓ j = l tj l tj + l tn − t ( ℓ < n − − t ) ( ℓ = n − − t ) . (17)Depending on the choice of k , Eqs. (11) and (13) suggestthat l t + i only takes on one of two possible values. For 1 ≤ i < n − t −
1, these are l t + i = ( l ti + l tn − t ( i ≤ k ≤ n − t − ) l ti − ( ≤ k < i ) , (18)whereas for i = n − t −
1, the two values are l t + n − t − = ( l tn − t ( k = n − t − ) l tn − t − ( ≤ k < n − t − ) . (19)We introduce the shorthand P ti = i − (cid:213) j = l tj l tj + l tn − t , (20) for the probability of k ≥ i at iteration t . On average, the evolu-tion of l ti is given by the rate equation l t + i = l ti − (cid:0) − P ti (cid:1) + (cid:0) l ti + l tn − t (cid:1) P ti , (21)for all 1 ≤ i < n − − t , and at the upper boundary l t + n − ( t + ) = l tn − − t (cid:0) − P tn − − t (cid:1) + l tn − t P tn − − t . (22)These equations are defined on the triangular domain 0 ≤ t ≤ n −
1, 1 ≤ i ≤ n − t . The initial conditions are l t = i = ( ≤ i ≤ n ) . (23)As described in Sec. 3, the process terminates at t = n − l n − characterizing the exponential variate for the final differ-ence in LDM. Yet, Eq. (22) for t = n − l n − = l n − ,reflecting the final, trivial differencing step, and it will proveconceptually advantageous to focus on the asymptotic proper-ties of l n − instead.Since the rate equation is an approximation to the exact re-cursion, we need to check how accurate it is. We have solvedthe rate equations (21-23) numerically up to n = · . Fig. 8shows ln (cid:0) l n − ( n + ) (cid:1) ln n from the rate equation versus 1 / ln n . If l n − were calculatedas an average from the exact recursion, it should be equal to − ln E [ L n ] ln n from the direct simulation of LDM. Fig. 8 shows this quantity,too. Apparently the error introduced by approximating the ex-act recursion by the rate equation vanishes for n → ¥ , and ourconjecture is that the rate equation and the exact recursion areasymptotically equivalent. Judging from our numerical studiesbelow, see Tab. 2, both asymptotic series have a relative differ-ence of size ln ln n / ln ( n ) .The time to solve the rate equation numerically scales like O (cid:0) n (cid:1) , so it is actually more efficient to simulate LDM di-rectly, not least because the sampling for the latter can be doneefficiently on a parallel machine. For analytic approaches, how-ever, the rate equation is more convenient.The initial probabilities decay exponentially, P i = − i , (24)which implies that only the first values l , l , . . . increase. Ev-erywhere else, P i is essentially zero, and those entries will notincrease until the first term of (21) has copied the values fromthe low index boundary. Hence we expect a “wavefront” of in-creased l -values to travel with a velocity one index per timestep toward the upper boundary, which in turn travels with thesame velocity towards the lower boundary. As can be seen fromFig. 6, this traveling wavefronts of increasing heights are a hall-mark of the rate equation for all times t . We will use this intu-itive picture for an Ansatz to analyze both the exact recursionand the rate equation in the next two sections. Stefan Boettcher, Stephan Mertens: Analysis of the Karmarkar-Karp Differencing Algorithm i t log ( λ i,t ) Figure 6.
Contour plot on a logarithmic scale for the numerical solu-tion l ti of the rate equations (21-23) at n = l ti ≃ t above that. The solution rises by about a decade betweeneach repeat of a band color. Note the ever more-rapid alternation be-tween narrowing and widening bands, signifying regions of rapid gaininterrupted by extended plateaus. The regular banded structure alongdiagonals t − i = const justifies the similarity solution in Eq. (38). Theonly notable exceptions occur in asymptotically diminishing regionsnear i = t = n / , n / , n / ,... . Both the exact recursion and the rate equations yield l t + = l t + l tn − t (25)for the lower boundary that we are ultimately interested in. Thisrecursion connects the lower and the upper boundaries at i = i = n − t . Unfortunately, l tn − depends in a complicatedway on entries of the l -tuple at different times and differentplaces. However, Fig. 6 suggests a similarity Ansatz l ti = l t − xi − x , (26)which makes the upper boundary readily available: l t + = l t + l t − n + ( ≤ t < n − ) l t = ( t ≤ ) (27)Note that we have extended the initial conditions l ti = l n − ofthis recursion in terms of the corresponding values in smallersystems, which leads to a simple recursion in n . To derive thisrecursion it is convenient to visualize (27) in terms of pathsin a right-angled triangle D n (Fig. 7). The hypotenuse of D n represents t and ranges from − n + n −
1, the height is n − n =
8. Thefinal recursion reads l = l + l , and the two terms on the right hand side correspond to twopaths: one that connects 6 with 7 along the hypotenuse, theother connects 5 with 7 along the path that is “reflected” at the − − − − − − − − − − − − − − Figure 7.
Proof of the Fibonacci recursion: The number of differentpaths from the leftmost point to the rightmost point in the triangle for n is the sum of the number of paths in the corresponding triangle ofsize [ n / ] (top) plus the number of paths in the triangle of size n − right leg of D . In our case a reflected path moves diagonallyupward until it touches the right leg above point i . From thereit moves downward to point i +
1. This peculiar “law of refrac-tion” implies that only every second point of the left half of thehypotenuse is connected to the right half by a reflected path.We can apply the recursion again and write l = l + l = l + l + l + l Here we have connected 6 with 5 along the hypotenuse andwith 3 along a reflected path, and similarly for 5. We iteratethis path finding process until all paths end on the left half ofthe hypotenuse (negative t ). Here the paths collect the initialvalues l t =
1, hence l equals the number of different pathsthat connect the points − , − , . . ., − −
7. The rule for path finding then is: if youare on an even index, move one unit to the right. If you areon an odd index, there are two branches: one to the right, theother 45 degrees upward and reflected down to the hypotenuse.Obviously, l equals the number of different paths that con-nects the leftmost point of D to the rightmost point accordingto this rules. Let T n ( i ) denote the number of paths that connectthe point i with n − D n . Then we have T n ( − n + ) = l n − . tefan Boettcher, Stephan Mertens: Analysis of the Karmarkar-Karp Differencing Algorithm 7 Now, starting at − n +
1, we have two choices: move upwardfor a reflection that will take us to point 1 or move along thehypotenuse to point − n + T n ( − n + ) = T n ( ) + T n ( − n + ) . As we can see in Fig. 7 (top), the number of paths from 1 to n − D n / .Hence T n ( ) = T n / ( − n / + ) . Similarly, the number of paths from − n + n − T n ( − n + ) = T n − ( − n + ) , and all three equations yield T n ( − n + ) = T n / ( − n / + ) + T n − ( − n + ) . The derivation of a corresponding equation for odd values of n is straightforward. If we define F ( n ) : = T n ( − n + ) = l n − , (28)the recursion for T n translates into the Fibonacci like recursion F ( n ) = F ( n − ) + F ([ n / ]) F ( ) = [ x ] refers to the integer part of x . The resulting sequenceis known as A033485 in [24]. The generating function g ( z ) = (cid:229) n F ( n ) z n satisfies the functional equation g ( z ) ( − z ) = z + ( + z ) g ( z ) , (30)and is given by g ( z ) = ( − z ) − (cid:213) k ≥ ( − z k ) − ! . (31) F ( n ) can be evaluated numerically for values of n that are largerthan the values feasible for simulations of LDM or for solvingthe rate equation. The bottleneck for calculating F ( n ) is mem-ory, not CPU time, since n / F ( n ) .With 3 GByte of memory, we managed to calculate F ( n ) for n ≤ · . We will derive the asymptotics of F ( n ) in the nextsection.Fig 8 shows F ( n ) within the same scaling as the simula-tions of LDM and the numerical solution of the rate equation.Apparently the similarity Ansatz does not capture the full com-plexity of the LDM algorithm or the rate equation. Yet it yieldsa very similar qualitative behavior. And in the next section wewill show that lim n → ¥ ln F ( n ) ln n =
12 ln 2 . (32) l n [ Z ( n + )] / l n n rate eq.LDMF(n) (Fibonacci)f(n) (cont. limit) Figure 8.
Four models of LDM: Direct simulation ( Z = / E [ nL n ] ),rate equation ( Z = l n − ), the Fibonacci model Z = F ( n ) from (29)and the similarity solution Z = f ( n ) of the continuous rate equation,given by (49). The dashed line represents (50). All dotted lines arenumerical fits of the type (51). To analyze the rate equations (21-23), it is convenient to con-sider the continuum limit for n → ¥ . Asymptotically, a con-tinuum solution may differ from the discrete problem in cor-rections of order 1 / n . As we will see, such corrections are in-accessible, as the asymptotic expansion is a series in terms of1 / ln ( n ) .We rewrite Eq. (21) in terms of discrete differences, l t + i − l ti = − (cid:0) l ti − l ti − (cid:1) + (cid:0) l ti − l ti − + l tn − t (cid:1) P ti . (33)Setting t = sn ( ≤ s ≤ ) , i = xn ( ≤ x ≤ − s ) , (34) l ti = y ( x , s ) , we obtain for large n n (cid:20) ¶¶ x + ¶¶ s (cid:21) y ( x , s ) = P ( x , s ) (cid:20) ¶ n ¶ x y ( x , s ) + y ( − s , s ) (cid:21) , (35)where we have set P ti → P ( x , s ) = exp (cid:26) n Z x d x ln a ( x , s ) (cid:27) (36)with a ( x , s ) = y ( x , s ) y ( x , s ) + y ( − s , s ) ≤ . (37)The left-hand side of Eq. (35), as well as the numerical solu-tion of the full rate equations (21-23) displayed in Fig. 6, againsuggest a similarity Ansatzy ( x , s ) = g ( s − x ) . (38) Stefan Boettcher, Stephan Mertens: Analysis of the Karmarkar-Karp Differencing Algorithm
This Ansatz yields immediately for Eq. (35):0 = P ( x , s ) (cid:20) − n g ′ ( s − x ) + g ( s − ) (cid:21) . (39)For almost all x >
0, the right-hand side vanishes by virtue of P ( x , s ) →
0, as indicated by Eq. (36) for a < n → ¥ .Correspondingly, P ( x = − s , s ) = P ( x = , s ) = s , hence we are leftwith 1 n g ′ ( s ) = g ( s − ) , (40)which can be interpreted as the continuous version of (27).From the initial conditions of the discrete problem in (23) itis clear that y ( x , ) =
1. For the similarity solution, this impliesthat g ( s ) = , ( − ≤ s ≤ ) . (41)Integrating (40), we formally obtain g ( s ) = g ( ) + n Z s d xg ( x − ) . (42)Thus, we can evaluate the integral for 0 ≤ s ≤ to get g ( s ) = + ns , (cid:18) ≤ s ≤ (cid:19) . (43)We can continue this process for ≤ s ≤ , i. e., 0 ≤ s − ≤ , exactly the domain of validity of (43), to obtain g ( s ) = g ( ) + n Z d xg ( x − ) + n Z s d xg ( x − ) , = + ns + n ( s − ) (cid:18) ≤ s ≤ (cid:19) . (44)The emergent pattern is best represented by defining g k ( s ) = g ( s ) , (cid:16) − − k ≤ s ≤ − − k (cid:17) , (45)for k = , , , . . . , where Eqs. (41-44) represent k =
0, 1 and 2.In general, we find that g k + ( s ) = g k (cid:16) − − k (cid:17) + n Z s − − k d xg k ( x − ) , (46)which is solved by g k ( s ) = k (cid:229) j = n j j ! 2 ( j ) (cid:0) j − s − j − + (cid:1) j . (47)For any n , we are interested in g ( s → ) ∼ lim t → n − l t , hence g ( ) = lim k → ¥ g k (cid:16) − − k (cid:17) = ¥ (cid:229) j = n j j ! 2 ( j ) , (48)which concludes our solution of (40). The sum for g ( ) stilldepends on n , hence we define f ( n ) = ¥ (cid:229) j = n j j ! 2 ( j ) . (49) Z f F l n − E [ nL n ] − c -1.44 -1.45 -1.22 -1.24 c -1.00 -1.42 -3.06 -3.86 c Table 2.
Parameters for (51) used in Fig. 8.
Now f ( n ) can be evaluated numerically for very large valuesof n . Fig. 8 shows the result for n ≤ . Don’t try this athome unless you have a computer algebra system. Interest-ingly, ln f ( n ) / ln n asymptotically approaches a value that isextremely close to 1 / [ f ( n )( n + )] ln n ≃
12 ln 2 + n (cid:18) ln ln 2 + + (cid:19) + n (cid:18) ln 2 + − ln ln 22 ln 2 (cid:19) − ln ln n ln n − ln ln n ln n + ln ln n ln n
12 ln 2 , (50)which is the dashed line in Fig. 8. The dotted lines are numer-ical least square fits of the ln ln n terms of this scaling, i.e., fitsof the formln [ Z ( n )( n + )] ln n ≃
12 ln 2 + n (cid:18) ln ln 2 + + (cid:19) + n (cid:18) ln 2 + − ln ln 22 ln 2 (cid:19) + ln ln n ln n c + ln ln n ln n c + ln ln n ln n c . (51)with values for c i as shown shown in Table 2. Note that the se-ries (49) as a solution of (40) and the first terms of the asymp-totic expansion (50) have been derived independently in thecontext of dynamical systems [25]. The numerical data supports the claim that the complete statis-tics of LDM is dominated by a single scale ∼ n − c ln n , not justthe expectation as described in (2). The available data is notsufficient to pin down the precise asymptotic scaling, however.In fact a naive extrapolation of the available data even con-tradicts the known asymptotic bound (3). With its O ( n ln n ) complexity, LDM is a very efficient algorithm, but probing theasymptotics requires ln n to be large. This discrepancy of scaleseliminates simulations as a means to study the asymptotics ofLDM and calls for alternative approaches.We have taken a step in the direction of a rigorous asymp-totic analysis by mapping the differencing algorithm onto a rateequation. The structure seen in the evolution of this rate equa-tion (Fig. 6) suggests a similarity Ansatz (26). With the helpof this Ansatz we could reduce the exact recursion in l -spaceto the Fibonacci model (29). The asymptotics of this modelcan be calculated, and it agrees with (2) and (3). The same tefan Boettcher, Stephan Mertens: Analysis of the Karmarkar-Karp Differencing Algorithm 9 Ansatz plugged into the rate equation even allows us to calcu-late the first terms of an asymptotic expansion (50). Althoughour Ansatz does not yield a proof, the extracted asymptoticbehavior satisfies all previous constraints and provides a con-sistent interpretation of the numerical results. Hence, our rateequations pave the way for further systematic investigations.
We appreciate stimulating discussions George E. Hentschel and CrisMoore. S.M. enjoyed the hospitality of the Cherry L. Emerson Cen-ter for Scientific Computation at Emory University, where part of thiswork was done. Most simulations were run on the Linux-Cluster T
INA at Magdeburg University. S.M. was sponsored by the European Com-munity’s FP6 Information Society Technologies programme, contractIST-001935, EVERGROW.
A Asymptotic Analysis
To evaluate the series (49) we apply Laplace’s saddle-pointmethod for sums as described on p. 304 of Ref. [26]. For a j = n j j ! 2 ( j ) = e f j , the saddle point is determined by D f j = f j − f j − =
0, i. e.,0 = D ln ( a j ) = ln (cid:0) a j / a j − (cid:1) , or1 = a j a j − = nj j − . (52)Hence, we obtain a moving ( n -dependent) saddle point at j ∼ ln n ln 2 − ln (cid:0) ln n ln2 (cid:1) ln 2 + + ln (cid:0) ln n ln2 (cid:1) ln 2 ln n − n + . . . , (53)including terms to the order needed to determine f ( n ) up tothe correct prefactor. We keep the 1 / ln ( n ) -corrections, since f j contains terms like j ln ( n ) . In particular, it is f j = j ln n − ln j ! − j ( j − ) . (54)As the saddle point j is large for large n , we can replace j !by its Stirling-series [26]. Then, we expand around the saddlepoint by substituting j = j + h , keeping only terms to order in h and those that are non-vanishing for n → ¥ . We find f j + h ∼ ln n −
12 ln ( p ) − ln 22 h ( h + ) + C ( n ) , with log-polynomial corrections C ( n ) ∼ − ln n ln 2 (cid:20) ln (cid:18) ln n ln 2 (cid:19) − − ln 22 (cid:21) + (cid:20)
12 ln 2 ln (cid:18) ln n ln ( ) (cid:19) − ln (cid:18) ln n ln ( ) (cid:19)(cid:21) . (55)We finally obtain for the asymptotic expansion of (50): f ( n ) ∼ Z ¥ − ¥ d h exp (cid:0) f j + h (cid:1) ∼ √ ln 2 exp (cid:26) ln n + C ( n ) (cid:27) (56) References
1. Narendra Karmarkar and Richard M. Karp. The differencingmethod of set partitioning. Technical Report UCB/CSD 81/113,Computer Science Division, University of California, Berkeley,1982.2. Stephan Mertens and Cristopher Moore.
The Nature of Computa-tion . Oxford University Press, Oxford, 2008.3. Michael R. Garey and David S. Johnson.
Computers andIntractability. A Guide to the Theory of NP-Completeness .W.H. Freeman, New York, 1997.4. Heiko Bauke, Stephan Mertens, and Andreas Engel. Phasetransition in multiprocessor scheduling.
Phys. Rev. Lett ,90(15):158701, 2003.5. R. C. Merkle and M. E. Hellman. Hiding informations and signa-tures in trapdoor knapsacks.
IEEE Transactions on InformationTheory , 24:525–530, 1978.6. Stephan Mertens. Random costs in combinatorial optimization.
Phys. Rev. Lett. , 84(6):1347–1350, February 2000.7. Christian Borgs, Jennifer Chayes, and Boris Pittel. Phase tran-sition and finite-size scaling for the integer partitioning problem.
Rand. Struct. Alg. , 19(3–4):247–288, 2001.8. Bernard Derrida. Random-energy model: Limit of a family ofdisordered models.
Phys. Rev. Lett. , 45:79–82, 1980.9. Bernard Derrida. Random-energy model: An exactly solvablemodel of disordered systems.
Phys. Rev. B , 24(5):2613–2626,1981.10. Heiko Bauke and Stephan Mertens. Universality in the levelstatistics of disordered systems.
Phys. Rev. E , 70:025102(R),2004.11. Anton Bovier and Irina Kurkova. Local energy statistics in dis-ordered systems: a proof of the local REM conjecture.
Commun.Math. Phys. , 263:513–533, 2006.12. Anton Bovier and Irina Kurkova. Local energy statistics in spinglasse.
J. Stat. Phys. , 126:933–949, 2007.13. Irina Kurkova. Local energy statistics in directed polymers.
Elec-tronic Journal of Probability , 13:5–25, 2008.14. Benjamin Yakir. The differencing algorithm LDM for partition-ing: a proof of a conjecture of Karmarkar and Karp.
Math. Oper.Res. , 21(1):85–99, 1996.15. Jon Kleinberg and ´Eva Tardos.
Algorithm Design . Addison Wes-ley, Boston, 2006.16. David S. Johnson, Cecicilia R. Aragon, Lyle A. McGeoch, andCatherine Schevron. Optimization by simulated annealing: anexperimental evaluation; part II, graph coloring and number par-titioning.
Operations Research , 39(2):378–406, May-June 1991.17. W. Ruml, J.T. Ngo, J. Marks, and S.M. Shieber. Easily searchedencodings for number partitioning.
Journal of Optimization The-ory and Applications , 89(2):251–291, May 1996.18. Richard E. Korf. A complete anytime algorithm for number par-titioning.
Artificial Intelligence , 106:181–203, 1998.19. Robert H. Storer, Seth W. Flanders, and S. David Wu. Problemspace local search for number partitioning.
Annals of OperationsResearch , 63(4):465–487, 1996.20. George S. Lueker. A note on the average-case behavior of asimple differencing method for partitioning.
Oper. Res. Lett. ,6(6):285–287, 1987.21. The GNU Multiple Precision Arithmetic Library. http://gmplib.org .22. TRNG - portable random number generators for parallel comput-ing. http://trng.berlios.de/ .23. William Feller.
An Introduction to Probability and Its Applica-tions , volume 2. John Wiley & Sons, 2 edition, 1972.0 Stefan Boettcher, Stephan Mertens: Analysis of the Karmarkar-Karp Differencing Algorithm24. N.J.A. Sloane. The on-line encyclopedia of integer sequences. .25. Cristopher Moore and Porus Lakdawala. Queues, stacks, andtranscendentality at the transition to chaos.
Physica D , 135:24–40, 2000.26. Carl M. Bender and Steven A. Orszag.