# Small integer programming problem

**URL:** <https://discourse.julialang.org/t/small-integer-programming-problem/121647>\
**Category:** Optimization (Mathematical)\
**Tags:** question\
**Created:** [October 23, 2024, 2:15pm UTC](https://discourse.julialang.org/t/small-integer-programming-problem/121647 "2024-10-23T14:15:14Z")\
**Posts on this page:** 8\
**Page:** 2

<div class="post-metadata">

**Author:** ![Tortar](https://avatars.discourse-cdn.com/v4/letter/t/6bbea6/32.png) [@Tortar](https://discourse.julialang.org/u/Tortar)\
**Post date:** [October 24, 2024, 6:21pm UTC](https://discourse.julialang.org/t/small-integer-programming-problem/121647/21 "2024-10-24T18:21:23Z")

</div>

Just some thoughts, both of this should hold

v-u \> 0

n\_1(v-u) - 1 \leq f(n\_1) = n\_1 v - \lceil n\_1 u \rceil \leq n\_1(v-u)

so that for all n\_1 \geq \frac{1}{v-u} then f(n\_1) \geq 0, so `next_lowest` shouldn’t be necessary to check from that point on. Same for n\_2. Which means that n \leq \max{(\lceil\frac{1}{v-u}\rceil, \lceil\frac{1}{t-s}\rceil)}. What about starting from this value rather than n=1? Didn’t think deeply about this idea actually, so I’m probably wrong, but I think these results should lead to some improvements of the algorithm given by @barucden

---

<div class="post-metadata">

**Author:** ![barucden](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/barucden/32/26154_2.png) [@barucden](https://discourse.julialang.org/u/barucden)\
**Post date:** [October 24, 2024, 6:34pm UTC](https://discourse.julialang.org/t/small-integer-programming-problem/121647/22 "2024-10-24T18:34:19Z")

</div>

Your upper bound is correct. However, we are looking for the smallest n, which can be much smaller than the upper bound you identified.

If we started from that upper bound n\_\text{max} and then checked n\_\text{max} - 1, n\_\text{max} - 2, \ldots, we would eventually find some n for which the condition holds. But there could another n' \< n for which the condition still holds and we wouldn’t know for sure until we checked all the remaining n-1, n-2, \ldots

Say that P(n) = 1 if \lceil nu \rceil \leq nv for given u\<v and otherwise P(n)=0. P(n) is not monotonic. You can find n such that P(n)=1 but P(n-1) = P(n+2) = 0. So you cannot even use binary search or something like that.

---

<div class="post-metadata">

**Author:** ![Benny](https://avatars.discourse-cdn.com/v4/letter/b/49beb7/32.png) [@Benny](https://discourse.julialang.org/u/Benny)\
**Post date:** [October 24, 2024, 10:15pm UTC](https://discourse.julialang.org/t/small-integer-programming-problem/121647/23 "2024-10-24T22:15:58Z")

</div>

> [@barucden](#):
>
> ```julia
> # Improvement: do larger steps
> n += 1
> 
> ```

I was going to suggest that once an `i,j` pair was identified for `n==i*a+j*b`, you could increment by `a` or `b` to the next candidate `n` and only have to check for `n == l*c+k*d`. But after working out an example, there appears to be an unbounded memory usage and `min` check for backtracking from increments of the larger number when iterating candidate n:

```julia
a=4 < b=9
4i+9j, i, j
4, 1, 0
8, 2, 0
9, 0, 1 new increment path, checked min(4*3+9*0, 4*0+9*1)
12, 3, 0 backtrack to previous path, n+=3
13, 1, 1 n+=1 because 1 == 9-2*4
16, 4, 0 backtrack, n+=3
17, 2, 1
18, 0, 2 new increment path, checked min(4*5+9*0, 4*3+9*1, 4*0+9*2)
20, 5, 0 backtrack, n+=2
21, 3, 1
22, 1, 2

```

which as far as I know is only avoidable when either `a` or `b` is a multiple of the other, making the diverging increment paths redundant because the sums are equal. When not backtracking, the increments of n are `1 == 9-2*4` when there are 4s to spare, and backtracking becomes less frequent and decreases the increments with more paths, so it’s not saving increments in the long run, either. I imagine other coprimes pairs would have a much larger set of increments, so I can’t simplify this.

---

<div class="post-metadata">

**Author:** ![rocco\_sprmnt21](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rocco_sprmnt21/32/20127_2.png) [@rocco\_sprmnt21](https://discourse.julialang.org/u/rocco_sprmnt21)\
**Post date:** [October 24, 2024, 11:20pm UTC](https://discourse.julialang.org/t/small-integer-programming-problem/121647/24 "2024-10-24T23:20:29Z")

</div>

```julia
a, b = 9, 35
c, d = 6, 25

gg = gcd(b,d)

k = b ÷ gg
h = d ÷ gg

h*b + 0*a
k*d + 0*c

[[[h1,j]'*[b,a] for h1 in h-j:-1:0] for j in 0:h]
[[[k1,j]'*[d,c] for k1 in k-j:-1:0] for j in 0:k]

g = gcd(a,c)

k = a ÷ g
h = c ÷ g

h*a + 0*b
k*c + 0*d

[[[h1,j]'*[a,b] for h1 in h-j:-1:0] for j in 0:h]
[[[k1,j]'*[c,d] for k1 in k-j:-1:0] for j in 0:k]

```

I propose this example to illustrate the idea.  
We start from an upper limit given by  
n0 = h_a+0_b = k_c+0_d

the other potentially valid values ​​are to be chosen within the vectors provided by the following formulas

```julia
[[[h1,j]'*[b,a] for h1 in h-j:-1:0] for j in 0:h]

[[[k1,j]'*[d,c] for k1 in k-j:-1:0] for j in 0:k]

```

So it would be a matter of comparing, at worst, two ordered vectors of (h+2)(h+1)/2 and (k+2)(k+1)/2 elements

---

<div class="post-metadata">

**Author:** ![tchr](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tchr/32/15138_2.png) [@tchr](https://discourse.julialang.org/u/tchr)\
**Post date:** [October 25, 2024, 7:34am UTC](https://discourse.julialang.org/t/small-integer-programming-problem/121647/25 "2024-10-25T07:34:21Z")

</div>

Not sure if relevant or a fruitful direction, but just figured I’d note that a relaxed version of the problem, with coefficients i, j, k, l in \mathbb{Z} (instead of \mathbb{Z}\_{\geq0}) is solvable by linear algebra, using the Smith normal form.

In particular, the original problem can also be phrased as:

\mathbf{A}\mathbf{v} = \begin{bmatrix} a & b & 0 & 0 \\ 0 & 0 & c & d\end{bmatrix}\begin{bmatrix} i \\ j \\ k \\ l \end{bmatrix} = n\begin{bmatrix} 1 \\ 1 \end{bmatrix}

The integer solutions to this equation, if they exist, can be expressed using the generalized inverse \mathbf{A}^g = \mathbf{T}^{-1}\boldsymbol{\Lambda}^g\mathbf{S}^{-1} constructed from the Smith normal form \mathbf{S}\mathbf{\Lambda}\mathbf{T} = \mathbf{A} with unimodular matrices \mathbf{S} and \mathbf{T} and a diagonal matrix of invariant factors \mathbf{\Lambda}. In particular, the solutions are:

\mathbf{v} = \begin{bmatrix} i \\ j \\ k \\ l \end{bmatrix} = n \mathbf{A}^g \begin{bmatrix} 1 \\ 1\end{bmatrix} + (\mathbf{1} - \mathbf{A}^g\mathbf{A}) \mathbf{t},

with \mathbf{t}\in \mathbb{Z}^4.  
One maybe useful thing in the original context, is that integer solutions exist if and only if n \mathbf{A}^g \begin{bmatrix} 1 \\ 1\end{bmatrix} is integral.

For one of the variants noted above, a = 9, b = 35, c = 6, d = 25, this procedure then gives the full set of solutions for i, j, k, l \in \mathbb{Z}:

\begin{bmatrix} i \\ j \\ k \\ l \end{bmatrix} = \begin{bmatrix} 4 \\ -1 \\ -4 \\ 1 \end{bmatrix} n + \begin{bmatrix} -35t\_1 - 140t\_2 \\ 9t\_1 + 36t\_2 \\ 25t\_3 + 100t\_4 \\ -6t\_3 - 24t\_4 \end{bmatrix}.

From the perspective of this relaxed version of the problem, the original problem then comes down to finding t\_i\in \mathbb{Z} such that the above is nonnegative.  
(Here, the minimal-n solution seems to be n=18 with t\_1, t\_2, t\_3, t\_4 = -6, 2, -5, 2 and i, j, k, l = 2, 0, 3, 0).  
This is the intersection of a cone-constraint and a lattice, i.e., a monoid - and a basis for the solutions should then be in principle obtainable from an associated Hilbert basis.

( **But** , not very clear to me that this is practically very useful for the original un-relaxed problem; computing the Smith normal form is a bit much, unless the `gcd` relations can somehow help getting it analytically, and computing the Hilbert basis is even worse. Still, sort of fun to see that part of what makes the original problem interesting is the non-negative constraint)

---

<div class="post-metadata">

**Author:** ![sgaure](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sgaure/32/14779_2.png) [@sgaure](https://discourse.julialang.org/u/sgaure)\
**Post date:** [October 25, 2024, 7:41am UTC](https://discourse.julialang.org/t/small-integer-programming-problem/121647/26 "2024-10-25T07:41:36Z")

</div>

> [@ambiso](#):
>
> There’s probably some connection to lattice reduction.
> 
> Here’s a first attempt using sagemath:

That’s my first thought as well. The package Nemo.jl has the LLL algorithm, so it should be possible to solve the shortest vector problem, and other related lattice problems there.

---

<div class="post-metadata">

**Author:** ![GunnarFarneback](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gunnarfarneback/32/1827_2.png) [@GunnarFarneback](https://discourse.julialang.org/u/GunnarFarneback)\
**Post date:** [October 25, 2024, 7:51am UTC](https://discourse.julialang.org/t/small-integer-programming-problem/121647/27 "2024-10-25T07:51:04Z")

</div>

> [@tchr](#):
>
> computing the Smith normal form is a bit much, unless the `gcd` relations can somehow help getting it analytically,

I worked out the Smith normal form for `A = [a b -c -d]` and it turned out, not all that surprisingly, to basically coincide with Bézout’s identity. I’m fairly certain that’s the case here as well. (As noted in other posts, Julia’s `gcdx` gives you what’s needed to solve for integer solutions for a fixed `n`.)

---

<div class="post-metadata">

**Author:** ![tchr](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tchr/32/15138_2.png) [@tchr](https://discourse.julialang.org/u/tchr)\
**Post date:** [October 25, 2024, 5:49pm UTC](https://discourse.julialang.org/t/small-integer-programming-problem/121647/28 "2024-10-25T17:49:58Z")

</div>

> [@tchr](#):
>
> This is the intersection of a cone-constraint and a lattice, i.e., a monoid - and a basis for the solutions should then be in principle obtainable from an associated Hilbert basis.

Following on this thread a bit, I thought it might be fun to compute this Hilbert basis for the specific problem of `a, b, c, d = 9, 35, 6, 25`. From this, we can get _all_ possible non-negative solutions - not just the minimum one (but also the minimum one).

We can compute the Hilbert basis using Normaliz, accessed via PyCall and PyNormaliz (the Normaliz.jl package isn’t working afaik):

```julia
using PyCall
const PyNormaliz = pyimport("PyNormaliz")

# set up problem from where it was left by generalized inverse approach
v = [4, -1, -4, 1]
Q = [-35 -140 0 0; 9 36 0 0; 0 0 25 100; 0 0 -6 -24]
S = hcat(v, Q) # solution space in ijkl is spanned by S*[n, t1, t2, t3, t4]

# compute Hilbert basis for monoid S*[n, t1, t2, t3, t4] >= 0
C = PyNormaliz.Cone(inequalities = S)
C.Compute("HilbertBasis", "DualMode")
nts = transpose(C.HilbertBasis())
ijkls = collect(eachcol(S*nts))
ns = nts[1,:]

```

The possible solutions in n and i,j,k,l are then given by `ns` and `ijkls`:

```julia
julia> println.(Ref("n = "), ns, Ref(": ijkl = "), ijkls)
n = 18: ijkl = [2, 0, 3, 0]
n = 62: ijkl = [3, 1, 2, 2]
n = 79: ijkl = [1, 2, 9, 1]
n = 81: ijkl = [9, 0, 1, 3]
n = 105: ijkl = [0, 3, 5, 3]
n = 106: ijkl = [4, 2, 1, 4]
n = 114: ijkl = [1, 3, 19, 0]
n = 125: ijkl = [10, 1, 0, 5]
n = 140: ijkl = [0, 4, 15, 2]
n = 149: ijkl = [1, 4, 4, 5]
n = 150: ijkl = [5, 3, 0, 6]
n = 175: ijkl = [0, 5, 0, 7]
n = 175: ijkl = [0, 5, 25, 1]
n = 210: ijkl = [0, 6, 35, 0]
n = 225: ijkl = [25, 0, 0, 9]

```

Any solution to ai + bj = ck + dl = n (here for `a, b, c, d = 9, 35, 6, 25`) can now be obtained from the non-negative integer span of the above `ijkls`.

[Previous page](https://discourse.julialang.org/t/small-integer-programming-problem/121647.md?page=1)