How use convex optimization for combinatorial ksum problem

Not a lot time ago I have been asked is it possible to solve k-sum problem via convex optimization.

Convex optimization can be used as a good heuristic for it.

Also it is rather interesting that 3-SUM problem is open problem in Computer Science. And nobody knows is it possible to solve it faster then ~N^2.  Here I would like to give extra details and the source code for whom it can be interesting for:

   - Brute force algorithm to solve k-sum problem. But it's exponential in variable "size of subset"

   - Convex-optimization based way with L1 heuristic and Iterated reweighting heuristic 

Approach was based on convex relaxation for solve cardinality problem. For details please look on S.Boyd book on convex optimization.

Complexity of one time evaluate "solve" on convex optimization problem with n variables and m constraints via various interior-point methods is rougly equal to 20-100 calls of procedure which takes time equal to: 

max {n^3, n^2*m, time to evaluate for objective and constraint functions it's values, derivatives and Hessians } 

And for record - for convex optimization problem there is no analytical solution in general (and most often).

Problem Setup

I’d like to write how it is possible to solve the following combinatorial problem: Which  k numbers among n give a number L?

Theoretical approach

Link to PDF and Link to MS Word docx if you have problems with reading embeded document.

Python Source Code

#!/home/bruzzo/anaconda2/bin/python2.7# Brute force algorithm and convex-optimization algorithm for solve k-summ combinatorial problem#   How to launch this scripts:#   1. Install cvxpy http://www.cvxpy.org/en/latest/install/index.html (this tool have been created by: Steven Diamond, Eric Chu, Stephen Boyd)#   2. If your OS support reading of sha-bang that tune sha-bang in the first line of this scriptimport cvxpy as cvx import numpy as np import matplotlib.pyplot as plt import random, time #======================================================================================#                 HELP FUNCTION FOR BRUTE FORCE IMPLEMENTATION#======================================================================================def maybeCarry(x, maxValue):     for i in xrange(len(x)):         if x[len(x) - i - 1] == maxValue:             x[len(x) - i - 1] = 0             x[len(x) - i - 1 - 1] = x[len(x) - i - 1 - 1] + 1         else:             returndef timeToFinish(x, maxValue):     for i in xrange(len(x)):         if x[i] != maxValue - 1:             return False     return Truedef isIncreasingSequence(x):     for i in xrange(len(x)-1):         if x[i+1] <= x[i]:             return False     return True#======================================================================================#                   PROBLEM CONFIGURATION#====================================================================================== n = 15   # Set size k = 7    # Subset size lo = -20 # Lower bound on items in array (not used in algorithms) hi = 40  # High  bound on items in array (not used in algorithms)# Target sums and various accumulators SumToGet = [24, 51, 55, 73, 73, 76, 89, 100] ReconstructedSum = [] TimeBruteForce = [] TimeCvxWay = []print "Size of array: ", n print "Subset size: ", k print "Items are between: [%i, %i]" % (lo, hi)   # Fill test array a = [] random.seed(12)for i in xrange(n):     a.append(random.randint(lo, hi))#======================================================================================#                        BRUTE FORCE ALGORITHM#======================================================================================for L in SumToGet:     usedIndicies = [0 for i in xrange(k)]     t0 = time.time()     while True:            if (timeToFinish(usedIndicies, n)):             break          if isIncreasingSequence(usedIndicies):             items = map(lambda i: a[i], usedIndicies)             sumItems = sum(items)             if sumItems == L:                 print "Found items with summ equal to ", L, " indicies: ", usedIndicies, " items: ", items                 break # Early break          usedIndicies[-1] += 1         maybeCarry(usedIndicies, n)     spendTime = time.time() - t0     print "", "BRUTE FORCE ALGORITHM COMPLETED IN: %s SECONDS" %(str(spendTime))     TimeBruteForce += [spendTime]#======================================================================================# Conex-optimization algorithm based on L1 Heuristic for cardinality problem#======================================================================================        x = cvx.Variable(n)   # Vector optimization variable     A = np.matrix(a)     w = np.ones(n)      # Create constraints     constraints = [x >= 0,                    x <= 1,                    cvx.sum_entries(x) == k,                    A*x - L == 0]      # Coefficient for iterated L1 weight heuristic     eps = 0.1      # Formulate problem and solve it and make reweighting step. 5 is heurisric number mentioned by S.Boyd in EE364A in context of iterated reweighting     t1 = time.time()     for i_reweight in xrange(5):         prob = cvx.Problem(cvx.Minimize(cvx.norm(cvx.diag(w)*x, 1)), constraints)         prob.solve()         for i in xrange(len(w)):             w[i] = 1.0/(eps+x.value[i])      printReport = True     if printReport:         print "status:", prob.status, "('optimal' - the problem was solved. Other statuses: 'unbounded', 'infeasible')"         print "optimal value (inf for objection): ", prob.value         print "optimal variables (x): ", x.value      spendTime = time.time() - t1     print "", "SOLVING VIA SERIES OF CONVEX OPTIMIZATION PROBLEMS: %s SECONDS" %(str(spendTime))      indicies = np.argsort(x.value.T)     indicies = indicies[0,::-1]     indicies = indicies[0,0:k]     print "indicies:", indicies     print "sum:", np.sum(A[0,indicies])     ReconstructedSum += [np.sum(A[0,indicies])]     TimeCvxWay += [spendTime]#==================================================================================================================================================== fig = plt.figure() ax = fig.add_subplot(111) plt.plot(range(1,len(ReconstructedSum)+1), SumToGet, 'r--', label='sum to get(target)') plt.plot(range(1,len(ReconstructedSum)+1), ReconstructedSum, 'r', label='reconstructed sum(via convex relaxation)') plt.xlabel('experiment number') plt.ylabel('Sum_target and Sum_reconstructed') plt.legend(loc='upper left') plt.grid()  fig = plt.figure() ax = fig.add_subplot(111) plt.plot(range(1,len(ReconstructedSum)+1), TimeBruteForce, 'b--', label='Brute force time to calc') plt.plot(range(1,len(ReconstructedSum)+1), TimeCvxWay, 'b', label='Cvx time to calc') plt.xlabel('experiment number') plt.ylabel('seconds') plt.legend(loc='upper left') plt.grid()  plt.show()#====================================================================================================================================================# Copyright (c) 2017, Konstantin Burlachenko (burlachenkok@gmail.com). All rights reserved.#====================================================================================================================================================

Result plots

Last update: Konstnatin Burlachenko, 10SEP2017