As J_H says, one of the ways to frame this problem is as a MIP/MILP. It might beat your (quintuple!) nested loop, but your mileage may vary based on input size and chosen solver, so you need to benchmark. There's a long list of gotchas:
- you'll either need to write a (maybe MPS) problem file generator and send off to a stand-alone solver, or pull in a Java FFI library that exposes one of the solver APIs for you.
- depending on which API you pick, you might have to reduce the list of expressions I demonstrate to matrix form.
- the problem cannot be expressed directly as ar_i + br_j = z because that's non-linear. You need to linearise via big-M constraints.
- the value of M should generally not be hard-coded, and leverage knowledge of the maximum value of arr.
Python equivalent
In PuLP/Python this looks like
import pulp
arr = (1, 2, 3, 6, 67)
Zarr = (1, 4, 30, 7, 88)
Z_MAX = 5
'''
a.arr[i] + b.arr[j] = z[k] for all k
what is a + b? a + b <= ZMAX
a >= 0, a int
b >= 0, b int
aarri + barrj = z
aarri >= a*1 - M(1 - seli)
'''
def solve(z: int) -> str | None:
isel = pulp.LpVariable.matrix(name='i', indices=range(len(arr)), cat=pulp.LpBinary)
jsel = pulp.LpVariable.matrix(name='j', indices=range(len(arr)), cat=pulp.LpBinary)
a = pulp.LpVariable(name='a', cat=pulp.LpInteger, lowBound=0, upBound=Z_MAX)
b = pulp.LpVariable(name='b', cat=pulp.LpInteger, lowBound=0, upBound=Z_MAX)
ari = pulp.LpVariable(name='ari', cat=pulp.LpContinuous)
brj = pulp.LpVariable(name='brj', cat=pulp.LpContinuous)
prob = pulp.LpProblem(name='diophantine', sense=pulp.LpMinimize)
prob.setObjective(a + b)
prob.addConstraint(name='z', constraint=ari + brj == z)
prob.addConstraint(name='iexcl', constraint=1 == pulp.lpSum(isel))
prob.addConstraint(name='jexcl', constraint=1 == pulp.lpSum(jsel))
M = 2*Z_MAX*max(arr)
for i, (arr_value, sel) in enumerate(zip(arr, isel)):
av = a*arr_value
relax = M*(1 - sel)
prob.addConstraint(name=f'lower_ari_{i}', constraint=ari >= av - relax)
prob.addConstraint(name=f'upper_ari_{i}', constraint=ari <= av + relax)
for j, (arr_value, sel) in enumerate(zip(arr, jsel)):
bv = b*arr_value
relax = M*(1 - sel)
prob.addConstraint(name=f'lower_brj_{j}', constraint=brj >= bv - relax)
prob.addConstraint(name=f'upper_brj_{j}', constraint=brj <= bv + relax)
# print(prob)
prob.solve()
if prob.status != pulp.LpStatusOptimal:
return None
ri = next(v for sel, v in zip(isel, arr) if sel.value() > 0.5)
rj = next(v for sel, v in zip(jsel, arr) if sel.value() > 0.5)
return f'{a.value()}*{ri} + {b.value()}*{rj} = {z}'
for z in Zarr:
print(solve(z))
LP construction
The above has the following problem description output, which shows the affine expressions necessary to constrain it:
diophantine:
MINIMIZE
1*a + 1*b + 0
SUBJECT TO
z: ari + brj = 1
iexcl: i_0 + i_1 + i_2 + i_3 + i_4 = 1
jexcl: j_0 + j_1 + j_2 + j_3 + j_4 = 1
lower_ari_0: - a + ari - 670 i_0 >= -670
upper_ari_0: - a + ari + 670 i_0 <= 670
lower_ari_1: - 2 a + ari - 670 i_1 >= -670
upper_ari_1: - 2 a + ari + 670 i_1 <= 670
lower_ari_2: - 3 a + ari - 670 i_2 >= -670
upper_ari_2: - 3 a + ari + 670 i_2 <= 670
lower_ari_3: - 6 a + ari - 670 i_3 >= -670
upper_ari_3: - 6 a + ari + 670 i_3 <= 670
lower_ari_4: - 67 a + ari - 670 i_4 >= -670
upper_ari_4: - 67 a + ari + 670 i_4 <= 670
lower_brj_0: - b + brj - 670 j_0 >= -670
upper_brj_0: - b + brj + 670 j_0 <= 670
lower_brj_1: - 2 b + brj - 670 j_1 >= -670
upper_brj_1: - 2 b + brj + 670 j_1 <= 670
lower_brj_2: - 3 b + brj - 670 j_2 >= -670
upper_brj_2: - 3 b + brj + 670 j_2 <= 670
lower_brj_3: - 6 b + brj - 670 j_3 >= -670
upper_brj_3: - 6 b + brj + 670 j_3 <= 670
lower_brj_4: - 67 b + brj - 670 j_4 >= -670
upper_brj_4: - 67 b + brj + 670 j_4 <= 670
VARIABLES
__dummy = 0 Continuous
0 <= a <= 5 Integer
ari free Continuous
0 <= b <= 5 Integer
brj free Continuous
0 <= i_0 <= 1 Integer
0 <= i_1 <= 1 Integer
0 <= i_2 <= 1 Integer
0 <= i_3 <= 1 Integer
0 <= i_4 <= 1 Integer
0 <= j_0 <= 1 Integer
0 <= j_1 <= 1 Integer
0 <= j_2 <= 1 Integer
0 <= j_3 <= 1 Integer
0 <= j_4 <= 1 Integer
This produces solutions
0.0*67 + 1.0*1 = 1
2.0*2 + 0.0*2 = 4
5.0*6 + 0.0*67 = 30
1.0*1 + 1.0*6 = 7
None
MPS problem format
If you can generate MPS (from Java or otherwise) looking like the following, then a few open-source solvers including CBC will accept it. This output is generated by PuLP.
*SENSE:Minimize
NAME MODEL
ROWS
N OBJ
E C0000000
E C0000001
E C0000002
G C0000003
L C0000004
G C0000005
L C0000006
G C0000007
L C0000008
G C0000009
L C0000010
G C0000011
L C0000012
G C0000013
L C0000014
G C0000015
L C0000016
G C0000017
L C0000018
G C0000019
L C0000020
G C0000021
L C0000022
COLUMNS
X0000000 OBJ 1.000000000000e+00
MARK 'MARKER' 'INTORG'
X0000001 C0000003 -1.000000000000e+00
X0000001 C0000004 -1.000000000000e+00
X0000001 C0000005 -2.000000000000e+00
X0000001 C0000006 -2.000000000000e+00
X0000001 C0000007 -3.000000000000e+00
X0000001 C0000008 -3.000000000000e+00
X0000001 C0000009 -6.000000000000e+00
X0000001 C0000010 -6.000000000000e+00
X0000001 C0000011 -6.700000000000e+01
X0000001 C0000012 -6.700000000000e+01
MARK 'MARKER' 'INTEND'
X0000002 C0000000 1.000000000000e+00
X0000002 C0000003 1.000000000000e+00
X0000002 C0000004 1.000000000000e+00
X0000002 C0000005 1.000000000000e+00
X0000002 C0000006 1.000000000000e+00
X0000002 C0000007 1.000000000000e+00
X0000002 C0000008 1.000000000000e+00
X0000002 C0000009 1.000000000000e+00
X0000002 C0000010 1.000000000000e+00
X0000002 C0000011 1.000000000000e+00
X0000002 C0000012 1.000000000000e+00
MARK 'MARKER' 'INTORG'
X0000003 C0000013 -1.000000000000e+00
X0000003 C0000014 -1.000000000000e+00
X0000003 C0000015 -2.000000000000e+00
X0000003 C0000016 -2.000000000000e+00
X0000003 C0000017 -3.000000000000e+00
X0000003 C0000018 -3.000000000000e+00
X0000003 C0000019 -6.000000000000e+00
X0000003 C0000020 -6.000000000000e+00
X0000003 C0000021 -6.700000000000e+01
X0000003 C0000022 -6.700000000000e+01
MARK 'MARKER' 'INTEND'
X0000004 C0000000 1.000000000000e+00
X0000004 C0000013 1.000000000000e+00
X0000004 C0000014 1.000000000000e+00
X0000004 C0000015 1.000000000000e+00
X0000004 C0000016 1.000000000000e+00
X0000004 C0000017 1.000000000000e+00
X0000004 C0000018 1.000000000000e+00
X0000004 C0000019 1.000000000000e+00
X0000004 C0000020 1.000000000000e+00
X0000004 C0000021 1.000000000000e+00
X0000004 C0000022 1.000000000000e+00
MARK 'MARKER' 'INTORG'
X0000005 C0000001 1.000000000000e+00
X0000005 C0000003 -6.700000000000e+02
X0000005 C0000004 6.700000000000e+02
MARK 'MARKER' 'INTEND'
MARK 'MARKER' 'INTORG'
X0000006 C0000001 1.000000000000e+00
X0000006 C0000005 -6.700000000000e+02
X0000006 C0000006 6.700000000000e+02
MARK 'MARKER' 'INTEND'
MARK 'MARKER' 'INTORG'
X0000007 C0000001 1.000000000000e+00
X0000007 C0000007 -6.700000000000e+02
X0000007 C0000008 6.700000000000e+02
MARK 'MARKER' 'INTEND'
MARK 'MARKER' 'INTORG'
X0000008 C0000001 1.000000000000e+00
X0000008 C0000009 -6.700000000000e+02
X0000008 C0000010 6.700000000000e+02
MARK 'MARKER' 'INTEND'
MARK 'MARKER' 'INTORG'
X0000009 C0000001 1.000000000000e+00
X0000009 C0000011 -6.700000000000e+02
X0000009 C0000012 6.700000000000e+02
MARK 'MARKER' 'INTEND'
MARK 'MARKER' 'INTORG'
X0000010 C0000002 1.000000000000e+00
X0000010 C0000013 -6.700000000000e+02
X0000010 C0000014 6.700000000000e+02
MARK 'MARKER' 'INTEND'
MARK 'MARKER' 'INTORG'
X0000011 C0000002 1.000000000000e+00
X0000011 C0000015 -6.700000000000e+02
X0000011 C0000016 6.700000000000e+02
MARK 'MARKER' 'INTEND'
MARK 'MARKER' 'INTORG'
X0000012 C0000002 1.000000000000e+00
X0000012 C0000017 -6.700000000000e+02
X0000012 C0000018 6.700000000000e+02
MARK 'MARKER' 'INTEND'
MARK 'MARKER' 'INTORG'
X0000013 C0000002 1.000000000000e+00
X0000013 C0000019 -6.700000000000e+02
X0000013 C0000020 6.700000000000e+02
MARK 'MARKER' 'INTEND'
MARK 'MARKER' 'INTORG'
X0000014 C0000002 1.000000000000e+00
X0000014 C0000021 -6.700000000000e+02
X0000014 C0000022 6.700000000000e+02
MARK 'MARKER' 'INTEND'
RHS
RHS C0000000 1.000000000000e+00
RHS C0000001 1.000000000000e+00
RHS C0000002 1.000000000000e+00
RHS C0000003 -6.700000000000e+02
RHS C0000004 6.700000000000e+02
RHS C0000005 -6.700000000000e+02
RHS C0000006 6.700000000000e+02
RHS C0000007 -6.700000000000e+02
RHS C0000008 6.700000000000e+02
RHS C0000009 -6.700000000000e+02
RHS C0000010 6.700000000000e+02
RHS C0000011 -6.700000000000e+02
RHS C0000012 6.700000000000e+02
RHS C0000013 -6.700000000000e+02
RHS C0000014 6.700000000000e+02
RHS C0000015 -6.700000000000e+02
RHS C0000016 6.700000000000e+02
RHS C0000017 -6.700000000000e+02
RHS C0000018 6.700000000000e+02
RHS C0000019 -6.700000000000e+02
RHS C0000020 6.700000000000e+02
RHS C0000021 -6.700000000000e+02
RHS C0000022 6.700000000000e+02
BOUNDS
FX BND X0000000 0.000000000000e+00
UP BND X0000001 5.000000000000e+00
FR BND X0000002
UP BND X0000003 5.000000000000e+00
FR BND X0000004
BV BND X0000005
BV BND X0000006
BV BND X0000007
BV BND X0000008
BV BND X0000009
BV BND X0000010
BV BND X0000011
BV BND X0000012
BV BND X0000013
BV BND X0000014
ENDATA
Matrix form
From Scipy or a Java equivalent, this is the sparse matrix form that you can use to define the problem constraints. It is trickier but more efficient to set up, particularly if you share common portions between your outer z loop iteration.
from functools import partial
from typing import Callable
import numpy as np
from scipy.optimize import milp, Bounds, LinearConstraint, OptimizeResult
import scipy.sparse as sp
def setup(arr: np.ndarray, ab_max: int) -> tuple[
Callable[[LinearConstraint], OptimizeResult], # binding to milp
sp.csc_matrix, # A matrix to constraint
np.ndarray, # lower and upper constraint bounds, z unpopulated
]:
n = len(arr)
M = 2*ab_max*arr.max()
'''
Columns:
a first term variable, int
b second term variable, int
ari = a*r_i, continuous
brj = b*r_j, continuous
isel (n) a binary assignment vector
jsel (n) b binary assignment vector
'''
# Common terms to big-M linearisation constraint rows
ab_rij = (
sp.kron(sp.eye_array(2), -arr[:, np.newaxis]),
sp.kron(sp.eye_array(2), np.ones((n, 1))),
)
msel = sp.diags_array(np.full(shape=2*n, fill_value=M))
A = sp.block_array((
( # ari + brj == z
[(0, 0)], [(1, 1)], None,
),
( # Kronecker: each term has one assignment
# sum(isel) == 1, sum(jsel) == 1
None, None, sp.kron(sp.eye_array(2), np.ones((1, n))),
),
( # Kronecker: lower bound for a*r_i, b*r_j
# ari >= a*arr_value - M*(1 - isel) -a*arr_value + ari - M*isel >= -M
# brj >= b*arr_value - M*(1 - jsel) -b*arr_value + brj - M*jsel >= -M
*ab_rij, -msel,
),
( # Kronecker: upper bound for a*r_i, b*r_j
# ari <= a*arr_value + M*(1 - isel) -a*arr_value + ari + M*isel <= M
# brj <= b*arr_value + M*(1 - jsel) -b*arr_value + brj + M*jsel <= M
*ab_rij, msel,
),
), format='csc')
# nan to be replaced with z
constraint_bounds = np.block([
[
np.array((np.nan, 1, 1)),
np.full(shape=2*n, fill_value=-M),
np.full(shape=2*n, fill_value=-np.inf),
],
[
np.array((np.nan, 1, 1)),
np.full(shape=2*n, fill_value=+np.inf),
np.full(shape=2*n, fill_value=+M)
],
])
c = np.zeros(2 + 2 + n + n)
c[:2] = 1 # minimize a + b
integrality = np.ones(shape=(2 + 2 + n + n), dtype=np.uint8)
integrality[2:4] = 0 # only ari/brj are continuous
bounds = Bounds(
# a b ari brj sel_ij
lb=( 0, 0, -np.inf, -np.inf) + (0,)*(2*n),
ub=(ab_max, ab_max, +np.inf, +np.inf) + (1,)*(2*n),
)
milp_bind = partial(
milp, c=c, integrality=integrality, bounds=bounds,
)
return milp_bind, A, constraint_bounds
def solve(
arr: np.ndarray,
milp_bind: Callable[[LinearConstraint], OptimizeResult],
A: sp.csc_matrix,
constraint_bounds: np.ndarray,
z: int,
) -> None | tuple[int, int, int, int]:
cbounds = constraint_bounds.copy()
cbounds[:, 0] = z
lba, uba = cbounds
result = milp_bind(
constraints=LinearConstraint(A=A, lb=lba, ub=uba),
)
if not result.success:
return None
(a, b, ari, brj), asel, bsel = np.split(
np.round(result.x).astype(int),
(4, 4 + len(arr)),
)
aval, = arr[asel.nonzero()]
bval, = arr[bsel.nonzero()]
return a, aval, b, bval
def demo() -> None:
arr = np.array((1, 2, 3, 6, 67))
milp_bind, A, constraint_bounds = setup(arr=arr, ab_max=5)
np.set_printoptions(linewidth=200)
print(A.astype(int).toarray())
for z in (1, 4, 30, 7, 88):
result = solve(
arr=arr, milp_bind=milp_bind, A=A, constraint_bounds=constraint_bounds, z=z,
)
if result is None:
print(f'z={z}: no solution')
else:
a, aval, b, bval = result
print(f'{a}*{aval} + {b}*{bval} = {z}')
if __name__ == '__main__':
demo()
[[ 0 0 1 1 0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 1 1 1 1 1 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 1 1 1 1 1]
[ -1 0 1 0 -670 0 0 0 0 0 0 0 0 0]
[ -2 0 1 0 0 -670 0 0 0 0 0 0 0 0]
[ -3 0 1 0 0 0 -670 0 0 0 0 0 0 0]
[ -6 0 1 0 0 0 0 -670 0 0 0 0 0 0]
[ -67 0 1 0 0 0 0 0 -670 0 0 0 0 0]
[ 0 -1 0 1 0 0 0 0 0 -670 0 0 0 0]
[ 0 -2 0 1 0 0 0 0 0 0 -670 0 0 0]
[ 0 -3 0 1 0 0 0 0 0 0 0 -670 0 0]
[ 0 -6 0 1 0 0 0 0 0 0 0 0 -670 0]
[ 0 -67 0 1 0 0 0 0 0 0 0 0 0 -670]
[ -1 0 1 0 670 0 0 0 0 0 0 0 0 0]
[ -2 0 1 0 0 670 0 0 0 0 0 0 0 0]
[ -3 0 1 0 0 0 670 0 0 0 0 0 0 0]
[ -6 0 1 0 0 0 0 670 0 0 0 0 0 0]
[ -67 0 1 0 0 0 0 0 670 0 0 0 0 0]
[ 0 -1 0 1 0 0 0 0 0 670 0 0 0 0]
[ 0 -2 0 1 0 0 0 0 0 0 670 0 0 0]
[ 0 -3 0 1 0 0 0 0 0 0 0 670 0 0]
[ 0 -6 0 1 0 0 0 0 0 0 0 0 670 0]
[ 0 -67 0 1 0 0 0 0 0 0 0 0 0 670]]
1*1 + 0*1 = 1
0*67 + 2*2 = 4
0*1 + 5*6 = 30
1*6 + 1*1 = 7
z=88: no solution
a,b
satisfying the constrains. Which one is the result? \$\endgroup\$