Solving variant of the knapsack/money-changing problem

Question

I'm trying to solve a variant of the knapsack/changing money problem where I have a set of a few numbers and I'm trying to find the linear (integer) combinations of them which are close to a given value.

I have implemented it with a brute force method that generates all the possible combinations up to some max value of each and then sorts them and finds the desired value using a binary search (from Leonid's answer to another question). With some effort with Compile I was able to speed this up such that it is workable, but I was hoping there was a more elegant/efficient solution that would scale better. Unlike the common money changing problem, I am interested in the closest values themselves (not just exact values) and I don't care how many possible combinations there are.

Here is what I have, calling this function with a list (values) and max value (nmax) returns a second function which is used to search for the closest set (within a specified range):

makearray[values_, nmax_] := Module[
  {countlist, sums, sorted},
  countlist = Rest@Tuples[Range[0, nmax], Length@values];
  sums = dot[values, countlist];
  sorted = sort[sums];
  (* for my specific application duplicates are not desired so I remove  
     them with a custom compiled function that relies on the sorted list,  
     although they are rare with random data *)
  sorted = Pick[sorted, deleteduplicates[sorted], 1];
  Function[{desired, delta}, window[sorted, {desired - delta, desired + delta}]]
  ]

Here are the helper functions, if you don't have a C compiler installed simply remove CompilationTarget -> "C":

bsearchMin[list_List, elem_] := Module[
  {n0 = 1, n1 = Length[list], m},
  While[n0 <= n1, m = Floor[(n0 + n1)/2];
   If[list[[m]] == elem, While[list[[m]] == elem, m++];
    Return[m - 1]];
   If[list[[m]] < elem, n0 = m + 1, n1 = m - 1]];
  If[list[[m]] < elem, m, m - 1]]

bsearchMax[list_List, elem_] := Module[
   {n0 = 1, n1 = Length[list], m},
   While[n0 <= n1, m = Floor[(n0 + n1)/2];
    If[list[[m]] == elem, While[list[[m]] == elem, m--];
     Return[m + 1]];
    If[list[[m]] < elem, n0 = m + 1, n1 = m - 1]];
   If[list[[m]] > elem, m, m + 1]];

window[list_, {xmin_, xmax_}] := 
 With[{minpos = bsearchMax[list, xmin], 
   maxpos = bsearchMin[list, xmax]}, 
  Take[list, {minpos, maxpos}] /; ! MemberQ[{minpos, maxpos}, -1]]

window[__] := {}

dot = Compile[{{v1, _Real, 1}, {v2, _Real, 1}}, v1.v2, 
   RuntimeAttributes -> {Listable}, Parallelization -> True, 
   RuntimeOptions -> "Speed"];

sort = Compile[{{m, _Real, 1}}, Sort[m], RuntimeOptions -> "Speed", 
   CompilationTarget -> "C"];

deleteduplicates = Compile[{{v, _Real, 1}},
   Block[
    {i, len = Length[v], output},
    output = Table[1, {i, len}];
    For[i = 2, i < len, i++,
     If[Compile`GetElement[v, i] == Compile`GetElement[v, i - 1], 
      output[[i]] = 0]
     ];
    output
    ], RuntimeOptions -> "Speed", CompilationTarget -> "C"];

Now using the code lets start with a random list:

list = RandomReal[1,6];

solvef = makearray[list, 10]; // AbsoluteTiming

(*
0.4 seconds so this is workable
*)

solvef[10, 0.05]; // AbsoluteTiming

(*
0. seconds, basically instant
*)

Could invoke FrobeniusSolve in a table with targets being those numbers you regard as sufficiently close. — Daniel Lichtblau, Mar 11 '13 at 23:41

s0rce · Accepted Answer · 2013-03-12T01:40:17.533

Based on @DanielLichtblau's extremely useful comment I implemented a solution based on FrobeniusSolve. Required a bit of multipliying and rounding to get everything to work as integers but it works perfectly for my application.

frobeniussolver[values_, desired_, delta_] := 
  Flatten[
     Table[
        FrobeniusSolve[Round[values*100], x], {x, Round[(desired - delta)*100], Round[(desired + delta)*100], 1}],
         1]

frobeniussolver[RandomReal[1, 6], 10, 0.01]; // AbsoluteTiming

(*

0.137008 seconds!

*)

score 6 · Answer 2 · answered Mar 12 '13 at 01:09

Not really what you want, but you can find the closest hit using FindMinimum, as below.

n = 6;
target = 10;

SeedRandom[1111];
vals = RandomReal[1, n];

coeffs = Array[c, 6];
c1 = Map[# >= 0 &, coeffs];
diff = coeffs.vals - target;
c2 = {diff <= e, -diff <= e};
c3 = Element[coeffs, Integers];
constraints = Join[c1, c2, {c3}]; vars = Join[coeffs, {e}];

FindMinimum[{e, constraints}, vars]

This takes a few seconds or so. On your example I got around a factor of 10 smaller residual.

Solving variant of the knapsack/money-changing problem

2 Answers2

Linked