Find integer solutions to a set of equations with more unknowns than equations

Question

I am trying to build a system for which I need to find a solution to a set of linear equations with (much) more variables than equations.

The problem boils down to the following:

Imagine a set of equations:

A = A1*X1 + A2*X2 + ... + AnXn
B = B1*X1 + B2*X2 + ... + BnXn

How can I find one or multiple (positive) integer solutions to this system?

Note: I have been looking at the apache-commons-math library but I couldn't find any directions on how to solve/find a solution of a system with free variables.

Quartus

I don't have a solution but can probably point you to the right direction: You are trying to solve a system of diophantine equations. The mathematical discipline dealing with problems like that is called number theory. Number theory algorithms are normally not included in numerical libraries. – Arequipa

More variables than equations? There's a simple answer: least square regression. There's no guarantee that the solutions will be integers. – Kudu

This type of task is called "integer programming", and is hard to solve in general. en.m.wikipedia.org/wiki/Integer_programming – Prague

That's the problem, the solutions need to be integers. The application I'm building tries to find an optimal distribution of stamps on boxes, given a pile of different stamps, and the post office won't accept halve stamps :) – Quartus

Without more details it's difficult to say, but probably a commons math SimplexSolver could help here. See stackoverflow.com/questions/32528928 – Brianna

What more details would you need to awnser the question? I've looked at that post before and at the SimplexSolver but I couldn't figure out how to do it. – Quartus

Does X1, X2, X3 mean X^1, X^2, etc? Or are they different variables entirely? – Marbling

@HansvanderLaan have you looked at the knapsack problem and solvers for that? en.wikipedia.org/wiki/Knapsack_problem #7775269 – Kufic

@FrankPuffer provided a helpful direction. Would you please take a look at mathworld.wolfram.com/DiophantineEquation.html? "A linear Diophantine equation (in two variables) is an equation of the general form ax+by=c, where solutions are sought with a, b, and c integers." – Veats

Not a Java algorithm, but please take a look at math.stackexchange.com/questions/20717/… – Veats

@RobotKarel, They are different variables. – Quartus

@Rajah, those links deal with solving single diophantine equations, not a set of them. – Quartus

@Koos, yes I have but I don't see how this problem is a knapsack problem. We then somehow have to enfore that in napsack A there are as many A1's as there are B1's in B and if we treat all A1's and B1's in the set of solutions as distinct objects, which will cause many duplicate reponses. – Quartus

@HansvanderLaan if your problem is "The application I'm building tries to find an optimal distribution of stamps on boxes, given a pile of different stamps" then each box is a knapsack which needs to be 'filled' with stamps. Ensure that the box is filled completely (overflow is allowed). Do that for every box, until you cannot fill any box anymore. Or even better, check stack overflow for stamps, letter #3827475 – Kufic

My problem with your approach: If you try to optimize something, you usually have an objective function (cost, space, time, ...) which you want to minimize. What do you actually want to optimize? If you have such a function, then you can use a Mixed Integer Program solver like CPLEX or Gurobi to find solutions which minimize your function. If you are not in an academic environment and these things are too expensive or too complicated, one could aim for a heuristic approach. – Hirsh

It would be great if you could respond to my comment. If you extend your question and state the problem you really want to solve (which you sketched in the comments), it might be easier to help you. Furthermore it would be great to know the magnitudes of your problem (how many lines, how many variables, how large are A, A_1 on average?) – Hirsh

for such tasks I usualy use this How approximation search works – Godsend

@HansvanderLaan looking at your goal (stamps and letters) I really think your question is not what you need. It's an optimisation problem (most likely knapstack problem). Even if you reached a linear system at some point of your computations, you don't need the exact solution to it. You need one of the approximate solutions, and the one that maximises something. – Ciscaucasia

Wiki/Diophantine equation#System of linear Diophantine equations – Goebbels

Answer 1

Use the Lenstra–Lenstra–Lovász lattice basis reduction algorithm to find the Hermite normal form of your system.

It's straaightforward in matlab http://fr.mathworks.com/help/symbolic/mupad_ref/linalg-hermiteform.html

Check NTL for c++ http://www.shoup.net/ntl/index.html

Ciscaucasia

Answer 2

This is taken from a relevent Wikipedia section.

1. Rewrite the equations with matrix notation AX=C.
2. Compute the Smith normal form of A. Roughly speaking, that is finding B=UAV where B[i][j] == 0 if i != j.
3. B(inverse(V))X=UC
4. As B[i][j] == 0 if i != j, it's trivial to find the value of (inverse(V))X and thus X.

Goebbels

Answer 3

Follow this example: suppose the equations are:

7 = x + 3y + 4z + 9w
12 = 4x + 2y + 3z + 7w

There are 2 equations and 4 unknowns. You can set 2 of the unknowns to be parameters, so the system will have as many equations as unknowns:

7 - (4z + 9w) = x + 3y
12 - (3z + 7w) = 4x + 2y

We will be using x and y as the unknowns. It is possible to solve it and leave w and z as parameters in the linear solution:

x = (22 - 3w - z)/10
y = (16 - 29w - 13z)/10

Now, you need to make the numerators divisible by 10, the denominator. If there is a solution, you can test all the parameters from 1 to 10.

In general, you do this:

• Choose parameters so that there are as many unknowns as equations. Try to leave the unknowns that generates the smallest absolute value for the determinant (in the example, it was 10, but choosing y and z would be better, as it would be |det|=3)
• Solve the linear system and generate the answer depending on the parameters
• Test the values of the parameters from 1 to the absolute value of the det (if there is a solution, you will find it at this point), until there is a integer solution for all the unknowns (that is why you choose the smallest determinant absolute value before) and check if it is positive (this was not illustrated in the example).

Sorry if it is brute force in the final step, but at least there is a possibility of minimizing the range of the test. Maybe, there could be a better way to solve the final system of diophantine equations, but I don't know any method.

EDIT1

There is a method to avoid brute forcing the last part. Again, in the example, you have to make:

22 = 3w + z (congruent, mod 10)
16 = 29w + 13z (congruent, mod 10)

Applying the modulus:

2 = 3w + z (congruent, mod 10)
6 = 9w + 3z (congruent, mod 10)

You can make the system of congruences triangular (Gaussian elimination), multiplying a congruence by inverses in the modulus 10 and summing up to the other ones. The inverse of 9 in the modulus 10 is -1, so we multiply the last congruence:

2 = 3w + z (congruent, mod 10)
-6 = -9w + -3z (congruent, mod 10)

Equivalent to:

2 = 3w + z (congruent, mod 10)
4 = w + 7z (congruent, mod 10)

Then, multiply by -3 and add to the first:

2 - 3*4 = 3w + z -3w - 21z (congruent, mod 10)
4 = w + 7z (congruent, mod 10)

Equivalent to:

-10 = -20z (congruent, mod 10)
4 = w + 7z (congruent, mod 10)

Then, you solve from top to bottom (this example seems to be true for any z value). There can be some degree of freedom if you have more parameters than congruences, but they are equivalent and you can set the excess parameters at any value, preferably the least which is 1.

Let me know if there is something not clear yet!

Horal

It's funny because you manage to transfer the problem to Z/10Z then construct a unimodular factorisation. You could use the EDIT1 triangulation directly on the first linear equation, and solve it instantly. – Ciscaucasia

yes, but it could sound too obvious because the coefficient is 1, so I tried to do something that shows the general calculations. – Iqbal

Answer 4

I would try the following approach (note that I never had to deal with that problem, so take this answer as an attempt to solve with you the problem instead of a real analytical answer).

You just find solutions just like if it is a regular system, so you may possibly find infinetly many solutions:

Example:

y = x + 3

then the real numbers pair (2,5) is a possible real solution for that system, once you have the infinitely many solutions you just restrict a subset of solutions that are made by whole numbers.

Of course you have constraints, in our case the solution has just 1 free variable, so we can write all solutions like that:

(x, x+3)

Surprise:

If there's a irrational number somewhere, there are no integer solutions:

(x, x+pi)  => neither 1st or 2nd number here can be whole at same time

So you can probably find integer solutions if and only if your "infinitely many solutions" are constrained by whole or by rational numbers.

Assume you have the following vector:

 ( 3x, (2/5)y, y, x, x+y)

Then a whole solution can be:

 x=3
 y=10/2

If you feel the answer is not satisfying for you, Just tell: I'll delete it to not gain the bounty points

Redwood

It wont work, you will end up with a similar problem. By the way, note that by construction the gaussian pivot solving won't reach any irrational number as a coefficient (and definitely not transcendent numbers). – Ciscaucasia

Answer 5

You might want to have a look at constraint solvers like Z3. There's also a JavaExample.

A helpful tutorial explaining some of the power of Z3 can be found here.

EDIT: Also have a look at Algorithm for solving systems of linear inequalities

Fallonfallout

Answer 6

If I understand the problem correctly, you have multiple packages, each with different postage costs. You want to pay these postage costs from a pool of stamps you have available. One can solve the problem with integer programming. The solution below solves for all packages at once. You will have a number of variables equal to numPackages*numStampDenominations (probably inconvenient for a large number of packages).

The equality constraint looks like Aeqx = beq. The matrix Aeq for two packages with 4 stamp denominations (numVarsnumPackages) looks like:

.21 .68 .47 .01 .00 .00 .00 .00           .89
                                   * x = 
.00 .00 .00 .00 .21 .68 .47 .01           .50 

The first row is the constraints coefficients (stamp values) for package 1. The non-zero coefficients are stamp values. The zero coefficients variable associated with other packages don't cares.

A second set of constraints (inequality) deals with the pool of stamps I have available. The inequality constraint looks like A*x = b. The matrix A for 4 stamp denominations and 8 variables (numPackages * numStampDenominations) looks like:

1 0 0 0 1 0 0 0         10

0 1 0 0 0 1 0 0         10
                 * x <=  
0 0 1 0 0 0 1 0         10

0 0 0 1 0 0 0 1         20

The cost function is f'*x which is just the total number of stamps. Its coefficients are just ones in the form of a column vector.

The following is a script run in matlab that solves the problem in the manner described. It would probably be formulated about the same in octave, a GNU offering similar to matlab. Matlab or octave may not be the right solution for you, but they at least tell you what works and give you a sandbox to work out a solution.

% The value of each stamp available as a 4x1 matrix
sv = [.21; .68; .47; .01];
% The number of each stamp available as a 4x1 matrix
sq = [10;10;10;40];
% Number of demominations
m = size(sv, 1);
% The postage required for each package as a 4x1 matrix
% -- probably read from a file
postage = [.89;.50;1.01;.47;.47];
% Number of packages -- just a count of the postage rows
packageCount = size(postage, 1);
% Number of variables is m*packageCount
numVar = m*packageCount;
% lower bound -- zero stamps of a given denomination
lb = zeros(numVar,1);
% upper bound -- use constraints for upper bound
ub = [];
% The cost function -- minimize the number of stamps used
% min(f'*x) 
f = ones(numVar,1);
% integer constraints
intcon = 1:numVar;
% Construct postage constraint matrix
Aeq = zeros(numVar, packageCount);

for i = 1:packageCount
    first = i*m - 3;
    last = first + 3;
    Aeq(first:last,i) = sv(:);
end

% Construct stamp count constraint matrix
A = zeros(numVar, m);

for r = 1:m
    for j = 1:m
        c = (j-1)*m + 1;
        A(c,r) = 1;
    end
end

x = intlinprog(f, intcon, A', b, Aeq', beq, lb, ub)

Neilla