2
$\begingroup$

I want an EFFICIENT way of solving the Frobenius Equation

a1 n1  + a2 n2 + a3 n3 + ... + aK nK = U

where n1, n2, ..., nK are restricted to be either 0 or 1. I also want to count the number of such solutions as a function of U for a fixed set of a1,a2, ....., aK.

I also want an efficient way of solving simultaneous Frobenius equations.

a1 n1  + a2 n2 + a3 n3 + ... + aK nK = U

b1 n1  + b2 n2 + b3 n3 + ... + bK nK = W
Improve this question
$\endgroup$
  • 2
    $\begingroup$ It looks like your problem is known as subset sum problem. A quick look at what the Wikipedia page has to say about it suggest that no efficient algorithms to solve it are known. $\endgroup$ – Kiro Feb 27 '19 at 12:56
  • 7
    $\begingroup$ Yours is an interesting problem. However, I wanted to leave you a comment regarding how your question reads. I may be overly sensitive to tone, but your wording sounds rude and demanding to me. You should consider sharing any efforts you have already made towards a solution, if any. Perhaps you could also propose a better algorithm with whose implementation we could help you. Remember though: this is a community of volunteers that donate their time. $\endgroup$ – MarcoB Feb 27 '19 at 13:07
  • $\begingroup$ The C code I used at mathematica.stackexchange.com/questions/192030 could be modified and generalized to solve this problem. It's not efficient in a computer-science sense of the word, but it gets the job done. I'm wondering though if there isn't a built-in function in Mathematica that does this job. $\endgroup$ – Roman Feb 27 '19 at 14:55
  • 2
    $\begingroup$ @Kiro Regarding "subset sum", yes, that's exactly what this is. General NP-ness notwithstanding, there are effective methods that work well for certain classes, based on a particular parameter that I cannot recall off the top of my head. Which is another reason a representative example would be useful here. $\endgroup$ – Daniel Lichtblau Feb 27 '19 at 16:41
  • $\begingroup$ Thanks for all the answers. Sorry for the demanding tone. Also the simultaneous Frobenius equation solutions are also very much appreciated.. $\endgroup$ – Quasar Supernova Feb 28 '19 at 5:11
13
$\begingroup$

A concrete example would be helpful. Here is a simple way to count solutions, illustrated by example.

Generate random set of smallish positive integers.

nlist = Union[RandomInteger[{8, 80}, 12]]

(* Out[1713]= {9, 16, 19, 35, 41, 42, 58, 59, 68, 74, 78} *)

Form a product of binomials.

prod = Times @@ (1 + x^nlist)

(* Out[1714]= (1 + x^9) (1 + x^16) (1 + x^19) (1 + x^35) (1 + x^41) (1 + 
   x^42) (1 + x^58) (1 + x^59) (1 + x^68) (1 + x^74) (1 + x^78) *)

Find out how many ways we can get 93 as a subset sum.

SeriesCoefficient[prod, {x, 0, 93}]

(* Out[1715]= 5 *)

--- edit ---

We can adapt this to get all the solutions, though in a way that for large problems will not be so efficient as using some of the other methods presented.

Set up as before, but now we throw in a new monomial for each power.

nlist = {9, 16, 19, 35, 41, 42, 58, 59, 68, 74, 78};
prod = Times @@ (1 + (Thread[x[Range[Length[nlist]]]]*t)^nlist);

We take the series coefficient in the common variable t, expand, group terms and factors into lists and sublists, and the powers in the new variables are the values selected for the subset sums. (That may have been terse, but it strikes me as better than some other explanations I have attempted in past.)

Apply[List, Expand[SeriesCoefficient[prod, {t, 0, 93}]]] /. 
  Times -> List /. x[_]^n_. :> n

(* Out[82]= {{16, 35, 42}, {16, 19, 58}, {35, 58}, {9, 16, 68}, {19, 74}} *)

--- end edit ---

Improve this answer
$\endgroup$
  • $\begingroup$ CoefficientList[prod, x] if you want all counts in one go. This is very fast. $\endgroup$ – Roman Mar 6 '19 at 6:26
  • $\begingroup$ @Roman Correct, assuming all coefficients are nonnegative. $\endgroup$ – Chip Hurst Mar 6 '19 at 13:59
  • $\begingroup$ Also when nlist grows larger, SeriesCoefficient becomes substantially faster. $\endgroup$ – Chip Hurst Mar 6 '19 at 14:24
8
$\begingroup$

For the first problem, if you have enough memory available you could just generate all subsets and count how many times they form a certain sum:

aa = {7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20};
Rest@BinCounts[Total /@ Subsets[aa]]

{1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 3, 3, 4, 4, 4, 4, 5, 6, 7, 8, 10, 10, 11, 12, 13, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 36, 38, 42, 45, 49, 51, 55, 58, 63, 66, 71, 75, 81, 84, 89, 94, 99, 103, 109, 114, 120, 125, 130, 136, 141, 146, 151, 158, 163, 168, 173, 179, 183, 189, 192, 199, 202, 207, 210, 216, 218, 222, 225, 229, 230, 234, 235, 238, 239, 240, 241, 242, 242, 242, 242, 241, 240, 239, 238, 235, 234, 230, 229, 225, 222, 218, 216, 210, 207, 202, 199, 192, 189, 183, 179, 173, 168, 163, 158, 151, 146, 141, 136, 130, 125, 120, 114, 109, 103, 99, 94, 89, 84, 81, 75, 71, 66, 63, 58, 55, 51, 49, 45, 42, 38, 36, 32, 30, 28, 26, 24, 22, 20, 18, 16, 14, 13, 12, 11, 10, 10, 8, 7, 6, 5, 4, 4, 4, 4, 3, 3, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1}

@DanielLichtblau's method is much faster than this though, and gives the same result when used with CoefficientList instead of SeriesCoefficient. My slower method, on the other hand, can be used to find out which subsets sum to a certain number.

For the problem with several simultaneous equations, Solve works but may be too slow:

aa = {1, 1, 1, 1}; A = 2;
bb = {1, 2, 3, 4}; B = 5;
nn = Array[n, 4];
nmax = 1;
nn /. Solve[Join[{aa.nn == A, bb.nn == B}, Thread[0 <= nn <= nmax]],
        nn, Integers]

{{0, 1, 1, 0}, {1, 0, 0, 1}}

A brief test suggests that the execution time of this solution scales polynomially with the number of unknowns.

As an efficiency example, this 100-variable search runs in 68 seconds on my laptop, and generates 33'401 solutions:

M = 100;
aa = ConstantArray[1, M]; A = 10;
bb = Range[M]; B = 100;
nn = Array[n, M];
nmax = 1;
Timing[F = nn /. Solve[
  Join[{aa.nn == A, bb.nn == B}, Thread[0 <= nn <= nmax]],
  nn, Integers];]
Improve this answer
$\endgroup$
5
$\begingroup$

Note that it's hard to know how 'efficient' you need without supplying specific examples.

If you're only interested in a single instance, you can reduce this problem to the 0-1 knapsack problem (or equivalently 0-1 ILP) by calling KnapsackSolve. Unfortunately I don't see an immediate (built in) way to find more than one solution with this approach.

(* Values from Danny's answer *)
A = {9, 16, 19, 35, 41, 42, 58, 59, 68, 74, 78};
b = 93;

If no solution exists, KnapsackSolve will still return a result -- a nearby RHS that's satisfiable. We can first use Danny's answer to verify that there is a solution. If there's not, we stop.

SeriesCoefficient[Times @@ (1+x^A), {x, 0, b}]

5

A single solution is then:

ksol = KnapsackSolve[Thread[{A, A, 1}], b]

{0.004771, {0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0}}

If you don't want to call SeriesCoefficient beforehand you can verify here:

A.ksol == b

True

If you need all solutions and your inputs are small enough, you can call FrobeniusSolve and pick the results you want.

sol = FrobeniusSolve[A, b]; // AbsoluteTiming
Pick[sol, UnitStep[Max /@ sol - 2], 0]

{0.040838, Null}

{{0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0}, {0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0}, 
  {0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0}, {0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0}, 
  {1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0}}

We can also recurse brute force, pruning the search space along the way. Here's one such solution. I didn't spend time making a faster implementation as I don't know the how large the intended inputs are.

binaryFrobenius[{}, _] = {};

binaryFrobenius[list_, 0] := {ConstantArray[0, Length[list]]}

binaryFrobenius[{k_, rest___}, n_] /; n < k := binaryFrobenius[{rest}, n]

binaryFrobenius[{k_, rest___}, n_] := Join[
  Prepend[0] /@ binaryFrobenius[{rest}, n], 
  Prepend[1] /@ binaryFrobenius[{rest}, n - k]
]

binaryFrobenius[A, b] // AbsoluteTiming

{0.000681, {{0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0}, {0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0}, 
   {0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0}, {0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0}, 
   {1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0}}}

Edit

Here's my attempt at leveraging a SAT solver.

I've managed to reduce subset sum to SAT, but not necessarily in polynomial space though. This method works, but is slow. I wonder how one can perform this reduction more efficiently.

For each ni we'll have ai boolean variables. Therefore this scales with the coefficients themselves, which is not ideal.

groupedvars = MapIndexed[Thread[x[First[#2], Range[#1]]]&, A];
vars = Flatten[groupedvars];

Assert that all variables corresponding to the same ni are the same:

coeffs = Equivalent @@@ groupedvars;

Assert that exactly b variables are true:

constraint = BooleanCountingFunction[{b}, Total[A]] @@ vars;

Verify that we only have 5 solutions:

SatisfiabilityCount[And @@ Prepend[coeffs, constraint]]

5

The explicit solutions:

insts = SatisfiabilityInstances[And @@ Prepend[coeffs, constraint], vars, All];

res = Boole[insts][[All, Prepend[Most[Accumulate[A] + 1], 1]]]

{{1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0}, {0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0}, 
  {0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0}, {0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0}, 
  {0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0}}

res.A

{93, 93, 93, 93, 93}
Improve this answer
$\endgroup$
2
$\begingroup$

For the simultaneous equations case, you can again use SeriesCoefficient as in Daniel's answer to get the number of solutions. For instance:

SeedRandom[1]
a = Range[10];
b = RandomSample[1;;30, 10];

Using SeriesCoefficient:

SeriesCoefficient[Times @@ (1 + x^a y^b), {x, 0, 12}, {y, 0, 37}]

2

Using Solve to confirm:

Solve[
    a.x == 12 && b.x == 37 && x ∈ RegionProduct @@ Table[Point[{{0},{1}}], 10],
    x
]

{{x -> {0, 0, 1, 1, 1, 0, 0, 0, 0, 0}}, {x -> {1, 1, 1, 0, 0, 1, 0, 0, 0, 0}}}

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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