3

I need to solve a system of linear diophantine equations with either numpy or sympy.

Is there any way to constrain numpy's linalg.solve/linalg.lstsq method to return only integer solutions? (probably not but thought I should ask)

I looked into Sympy's diophantine solver and it does not seem to be applicable to solving whole systems

The problem I am working on is something along the lines of

P1(X) + P2(Y) = TargetPro
F1(X) + F2(Y) = TargetFat
C1(X) + C2(Y) = TargetCarb

In this case, X,Y,Z would represent the approximate serving sizes, and P1/F1/C1 would be the pro/fat/carb profile respectively.

based on this paper https://www.math.uwaterloo.ca/~wgilbert/Research/GilbertPathria.pdf

it seems as though I could conduct a row-reduction to find the ref of this system (row-echelon form) and then plug it into sympy's solver.

Is there an easier way to go about it?

Here's a trivial example:

pro = [4,5]
fat = [1,2]
carb = [3,6]

A = np.array((pro, fat, carb))
b = np.array([22,12,21])

print(np.linalg.lstsq(A, b))

I expected to get an integer solution [3,2] and instead got [ 2.16666667, 2.66666667]

Both solutions are correct, but I want to bound my solutions to only integer solutions

smichr
  • 16,948
  • 2
  • 27
  • 34
Anthony Chung
  • 1,467
  • 2
  • 22
  • 44
  • 2
    Unless one of equations is redundant, a system of three linear equations in three variables has just one solution in real numbers. It this solution is in integers, you have what you want; if not, there is no solution in integers. –  May 08 '16 at 16:08
  • Could you make the title of your question more specific? A small example with the desired output would also help. The reader has to guess what kind of objects P1 and TargetFat are. – Han-Kwang Nienhuys May 08 '16 at 16:11
  • @Han-KwangNienhuys. I will provide an example shortly. working on it now – Anthony Chung May 08 '16 at 16:14
  • 1
    You should mention in the title/description that you are looking for a *least-squares* solution rather than an exact solution. Your example has no exact integer-valued solution. – ali_m May 08 '16 at 19:04
  • You now give an example where you have more equations than unknowns (3 versus 2). Is that generally the case? Also, the paper that you're referring to is about equations with not only integer solutions but also integer equations. Is that the case for you? Because then, the example of food types is a bit misleading. – Han-Kwang Nienhuys May 08 '16 at 19:20
  • In the example I provided, [3,2] is a valid solution. – Anthony Chung May 08 '16 at 19:20
  • Yes I simplified the example further, but I can have anywhere from 2x3 to a 7x3 matrix. The system of equations are composed of integers. – Anthony Chung May 08 '16 at 19:22
  • 2
    `[3, 2]` and `[ 2.16666667, 2.66666667]` might be "valid" solutions in the least-squares sense, but neither are exact solutions. `A·[3, 2] = [22, 7, 21]`, not `[22, 12, 21]` as your `b` vector is defined. – ali_m May 08 '16 at 20:09
  • welp. I multiplied my test case carelessly. you're right @ali_m haha for my purposes I'll use the lstsq method then – Anthony Chung May 08 '16 at 20:58
  • So, do your real systems have exact integer solutions or not? – ali_m May 08 '16 at 21:59
  • they might or might not. I'd like to grab integer solutions if possible but in all likelihood there will be some approximation necessary – Anthony Chung May 08 '16 at 22:10
  • Is it necessary to use scipy and numpy only? I think it'd be easier to formulate an integer program and use IP solver like PyGLPK directly. For instance: `max 0 subject to a11*x+a12*y = b1, a21*x + a22*y = b2 etc...`, if your system is solvable then the above IP is feasible and has optimal solution. – serge_k May 09 '16 at 06:56
  • I looked at the documentation here http://tfinley.net/software/pyglpk/ex_ref.html but I'm not sure how linear programming (seeking to max/min a target) would help me if I'm trying to find an exact (if possible) or approximate solution. Can you elaborate on how I could do that with PyGLPK? – Anthony Chung May 09 '16 at 23:51
  • Actually I showed an example of how to do it using IP. If you're trying to solve the system of linear equations you don't need to optimize any function (that's why the objective function is zero or any other constant) and the problem reduces to just cheaching wether the program is feasible or not and finding the solution. Plus, if you're not seeking to max/min a target then why are you using least-squares in your example? – serge_k May 10 '16 at 06:48
  • The coefficient matrix can be up to 3 x N, linalg.solve will only work with square matrices, so lst squares provides the closest approximation. Would IP give a closer approximation or should I just stick to trying to interpret least squares? – Anthony Chung May 10 '16 at 12:22
  • Furthermore, we can reframe it as a "minimize carb intake" while meeting the other nutritional requirements of fat/protein. Would that enable me to access the spectrum of feasible solutions in a LP program? – Anthony Chung May 10 '16 at 14:52

0 Answers0