Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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}},
    {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]
    ], 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
share|improve this question
Could invoke FrobeniusSolve in a table with targets being those numbers you regard as sufficiently close. – Daniel Lichtblau Mar 11 '13 at 23:41
up vote 7 down vote accepted

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_] := 
        FrobeniusSolve[Round[values*100], x], {x, Round[(desired - delta)*100], Round[(desired + delta)*100], 1}],

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


0.137008 seconds!

share|improve this answer

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

n = 6;
target = 10;

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.

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.