Find an algorithm for quickly computing results to a given equation

Question

There is an equation:

280a + 80b + 75c + 50d + 25e - 30f - 42g = R

• a, b, c, d, e, f, g - those are modifiers. They're strictly positive integer values.

• R - this is a solution range, the equation is valid as long as the result of its left side falls into that range.

---

The goal is:

• For each possible subset of the variable set find the smallest sum of modifiers values such that the result of the left side of the equation is within given solution range.

Subsets are defined by inclusion of modifiers. I.e., modifier can either == 0 or > 0.
For example:

• b, c, d, e, f, g == 0 and a > 0 is 1st subset

• a, c, d, e, f, g == 0 and b > 0 is 2nd subset

• c, d, e, f, g == 0 and a, b > 0 is 3rd subset

• etc...

Basically it's like counting in binary. Ergo there's a total of 128 subsets (technically 127 since a subset where all modifiers == 0 isn't needed as a solution).

Note: I know that e can replace c, d modifiers as both of these are divisible by it, but I still need those as separate modifiers for the sake of unique permutations.

Basically, the solution range is the input and the result should be a list of all valid subsets with the smallest sum of modifiers

Here's an example of a valid result:

• Let the range we input be from -2 to -2

Here's a list of smallest solutions for a couple different subsets (not all possible) with just b, c, d, e, f, g modifiers: https://imgur.com/a/ZlWeAlq (permutation == subset, I didn't use the correct terminology at the time).

Plugging these values in the equation will result in -2 on every subset (assume a modifier is 0, I didn't include it in my first solver).

---

I made a quick bruteforce solution in C# with a parallel loop just going over all possible values (within boundaries) until it finds all solutions, but it's painfully slow.

A couple algorithms I have considered:

1. Breadth-first search:

   • Add 1 of each modifier and check if any individual "node" is a solution. Then add 1 of each modifier to each previous node and see if any of the new nodes are solutions.

     1. I'll have to somehow optimize it to not check the same path multiple times since adding 2 modifiers in different order will lead to the same node

     2. This grows really fast in terms of memory, so I'll have to limit it to some range of values within which I allow this "branching", and then I'll only allow addition of identical modifiers until that "branching range" is close enough to the solution range. Something like:

        • if a solution range is smth like 500-550, I'll only allow modifiers with positive coefficients to be added to itself until the value on each node is within 100 of the solution range

2. Simplify an arbitrary working answer:

   • Calculate the minimum difference you can create given the selected modifier set

   • Calculate how many minimum differences added together it takes to get into a solution range which can give a large number of modifiers being used but its a working answer

   • Use a table to quickly simplify/neutralize the modifiers with each other until you get the minimum possible

3. Bruteforce, but with even more optimization:

   • I could add the same range check to a bruteforce algorithm so that it only tries every combination within a given range, as well as remove some unnecessary repetitions (my current solution sucks quite hard)

Another hopefully decent optimization I can think of right off the bat is to remove impossible subsets for a given range at the very start.

For example, if a range is positive, then subsets where only f, only g or only f, g modifiers are > 0 will never be a valid answer since both of these have negative coefficients (can't ever reach a positive range by adding negative values).

---

So my question is, which algorithm, either from the ones that I thought of or from other existing algorithms that I didn't think of, would be able to quickly and efficiently solve this problem?

I have yet to try to implement any of these since I'm not sure whether they will be good at all and if there are any existing algorithms that would be much better than any of the things I could think of. My bruteforce solution still exists and, if needed, I can link it. It is terribly slow tho and definitely not worth comparison in its current form.

Answer 1

Formulate this as one (highly-separable, but) large linear program. This reduces the amount of setup overhead at the expense of modest memory consumption, so long as you use sparse matrices correctly for your constraints.

The biggest drawback of this approach is that if any one permutation is infeasible, the entire problem is infeasible.

Set your optimization coefficients (c) to all 1. There will be 448 decision variables, one for each non-zero position in your permutation matrix (do not include zero-positions in the problem). All decision variables are integral and have a lower bound of 1. Each permutation has two constraint rows: one for the lower bound and one for the upper bound. Depending on which solver you use, you either need to (as below) represent everything in one inequality matrix, or separate into lower and upper constraint matrices. When implemented correctly this should execute in less than a second.

import numpy as np
import scipy.sparse
from scipy.optimize import linprog

coef = np.array((280, 80, 75, 50, 25, 30, 42))
n = len(coef)
m = 2**n - 1

perm_index = np.arange(m, dtype=np.uint8)
perm_bits = np.unpackbits(
    1 + perm_index[:, np.newaxis],
    axis=1, count=n, bitorder='little',
)  # 127x7
perm_mask = perm_bits.astype(bool).ravel()    # 889
perm_index = perm_index.repeat(n)[perm_mask]  # 448

low = 500
high = 600
c = np.ones_like(perm_index)  # 448

'''
280  .    .   . . . . (448)
.    80   .   . . . .
.    .    280 80 
.    .    .   . 75 . .
(127)
'''
sparse_coef = scipy.sparse.coo_array((
    # overload 4: data, (i, j)
    np.tile(coef, m)[perm_mask],  # data
    (
        perm_index,         # row indices
        np.arange(len(c)),  # column indices
    ),
)).tocsr()  # 127x448
A_ub = scipy.sparse.vstack((-sparse_coef, sparse_coef))  # 254x448
b_ub = np.concatenate((
    np.full(shape=m, fill_value=-low),
    np.full(shape=m, fill_value=high),
))  # 254

res = linprog(c=c, integrality=c, bounds=(1, None), A_ub=A_ub, b_ub=b_ub)
assert res.success, res.message

x_full = np.zeros(m*n, dtype=int)  # 889
x_full[perm_mask] = res.x
x_full = x_full.reshape((m, n))  # 127x7
print(x_full)

[[ 2  0  0  0  0  0  0]
 [ 0  7  0  0  0  0  0]
 [ 1  3  0  0  0  0  0]
 [ 0  0  7  0  0  0  0]
 [ 2  0  1  0  0  0  0]
 [ 0  6  1  0  0  0  0]
 [ 1  1  2  0  0  0  0]
 [ 0  0  0 10  0  0  0]
...
 [ 0  4  1  1  1  1  1]
 [ 1  1  1  1  1  1  1]]