Algorithm about finding a combination that the sum is closest to a given number

Given a matrix $$P\in \mathbb{R}^{n*k}$$ (just for ease of notation, no matrix or linear algebra is actually needed; bound to $$(0,1)$$ if necessary), select one number from each row and compute the sum. Find the choice that the produced sum is closest to a given number $$x$$.

$$\arg\min_{1\leq j_1,\cdots,j_n\leq k}{|x - \sum_{i=1}^n{P_{ij_i}}|}$$

An direct algorithm with time complexity $$k^n$$ is trivial. For $$n=2$$, I could come up with a $$k\log k$$ algorithm using binary search. But for $$k>2$$, all I could think of is to transform the matrix to a $$R^{2*k^{n/2}}$$ equivalent, then applying algorithm for $$k=2$$. This could be a cutdown, but the complexity is still exponential.

Could I do better than this? Is there any relevant topic, algorithm or data structure that I'm missing?

• Is the fact that the values are elements of a matrix even relevant? This seems equivalent to the elements just being in a list. Some relevant search terms are subset sum problem, which is equivalent to a knapsack problem. There are lots of algorithms for these problems. Commented Jun 1, 2023 at 3:54
• Sorry for the potential confusion. I just found it easier to express with a matrix. Otherwise, no matrix or linear algebra is involved. Commented Jun 1, 2023 at 4:36
• There might be multiple solutions. Commented Jun 1, 2023 at 5:51
• @whpowell96 the matrix is relevant insofar as there's a fixed stride in which selection groups are exclusive. This could technically be represented as a list with definitions for such strides, but that would be awkward. Commented Jun 2, 2023 at 23:40
• @Reinderien Thanks for pointing out. I fixed my notation. A bound is added in case of necessity. Commented Jun 3, 2023 at 0:36

A typical - though probably not maximally efficient - approach is mixed-integer linear programming. Use binary assignment variables for your element selection, and minimize the absolute error from your total.

Note that your "closest to a given number" is not well-defined. Here I show the Manhattan norm. If you need another norm (Frobenius, etc.) then you will have to use a different approach than LP.

from time import monotonic

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

n = 20
k = 30
rand = default_rng(seed=0)
P = rand.random((n, k))
x = rand.uniform(low=0.3*n, high=0.7*n)

# Decision variables: nk selection of each cell, binary
# Auxiliary variable: sum error absolute upper bound, continuous
c = np.zeros(n*k + 1)
c[-1] = 1  # minimize sum error

integrality = np.ones(n*k + 1, dtype=int)
integrality[-1] = 0

bounds = np.zeros((n*k + 1, 2))
bounds[:-1, 1] = 1
bounds[-1, 1] = np.inf

A_eq = scipy.sparse.bmat((
(   # exactly one element selected per row (kronecker)
scipy.sparse.kron(
scipy.sparse.eye(n),
np.ones(k),
),
scipy.sparse.csr_array((n, 1)),
),
))
A_ub = np.block([
[ P.ravel(), -1],  #  R.assign - x <= sumerr: R.assign - sumerr <= x
[-P.ravel(), -1],  # -R.assign + x <= sumerr: -R.assign - sumerr <= -x
])

start = monotonic()
result = linprog(
c=c, integrality=integrality, bounds=bounds, method='highs',
A_eq=A_eq, b_eq=np.ones(n),
A_ub=A_ub, b_ub=(x, -x),
)
stop = monotonic()
print(stop - start, 's')
print(result.message)

assignments = result.x[:-1].reshape((n, k))
xappr = (P * assignments).sum()
# print('P:')
# print(P)
# print('Assignments:')
# print(assignments)
print(xappr, '~', x, ' err:', xappr-x)

2.6089999999385327 s
Optimization terminated successfully. (HiGHS Status 7: Optimal)
6.04941432208226 ~ 6.049414630274225  err: -3.0819196528142356e-07

• Thanks for your answer and effort. But there may be a bit misunderstanding in your statement and code (or my description is not clear enough). The target $x$ is a decimal scalar, and the sum of chosen numbers is also a decimal scalar (see the equation). For two scalars, the closeness is clearly defined. Commented Jun 3, 2023 at 7:42
• @obfish Re. scalars - that's already what this code does. Re. closeness of two scalars: incorrect. There are many norms to choose from, for example least squares which is common in optimization but not expressible in an LP context. Commented Jun 4, 2023 at 2:08