# Solving variant of the knapsack/money-changing problem

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[CompileGetElement[v, i] == CompileGetElement[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

## 2 Answers

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!

*)


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

n = 6;
target = 10;

SeedRandom;
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.