How can I enforce both x ≥ 0 and Ax ≥ b when solving a least squares problem in Python?

You won’t find a specialized function that directly “forces” both non-negativity on x and the extra condition Ax ≥ b. In other words, neither np.linalg.lstsq nor scipy.optimize.nnls supports additional linear inequality constraints on Ax.

What you can do instead is reformulate your problem as a linear (or quadratic) programming problem. For example, if you have an objective (say, minimizing the norm of x or some cost function) and you need to enforce x ≥ 0 and Ax ≥ b, you can use a general-purpose optimizer like SciPy’s linprog (for linear objectives) or a quadratic programming solver if your objective is quadratic.

Below is an example using SciPy’s linprog. Note that linprog uses constraints of the form A_ub x ≤ b_ub, so the constraint Ax ≥ b must be reformulated as –Ax ≤ –b.

PYTHONCopied
import numpy as np
from scipy.optimize import linprog

# Define your arrays (example shapes used, replace with your actual A and b)
# A is 37 x 5, b is a 37 element vector
# (Your v1, v2, v3, v4, v5, and b here)
v1 = np.array([0.00000961, 0.0000011, 0.0000011, 0.000015, 0.00000884, 0.00000286, 0,
               0.00000006, 0, 0.000196, 0, 0.0000071, 0.000000023, 0.000131, 0.00038,
               0, 0, 0.00000161, 0, 0, 0.0000069, 0.00027, 0.000005, 0, 0,
               0.00054, 0.00475, 0.000000002, 0.00036, 0.0000032, 0, 0, 0.033, 0.0015,
               0.02, 0.00052, 0.207]).T
v2 = np.array([0, 0.0000064, 0.00000135, 0.000121, 0.0000177, 0.00000348, 0,
               0.00000006, 0, 0, 0, 0.0000833, 0, 0.000525, 0.00062,
               0, 0, 0.0000114, 0, 0, 0.0000458, 0.00168, 0.0000193, 0, 0,
               0.00376, 0.00705, 0.000000072, 0.00018, 0.0000327, 0, 0, 0.085,
               0.492, 0.258, 0.0628, 0.161]).T
v3 = np.array([0, 0.000009, 0.00000193, 0.0000196, 0.00000899, 0,
               0, 0.00000444, 0, 0, 0, 0.0000021, 0.000000056, 0.000664,
               0.00123, 0, 0, 0.00000841, 0, 0, 0.0000502, 0.00171,
               0.0000106, 0, 0, 0.00352, 0.0148, 0.000000032, 0.00005,
               0.0000365, 0, 0, 0.155, 0.0142, 0.216, 0.00366, 0.624]).T
v4 = np.array([0, 0.00000763, 0.00000139, 0.00000961, 0.0000135, 0.00000119,
               0, 0.00000056, 0, 0, 0, 0, 0, 0, 0.00054, 0, 0, 0.00000626,
               0, 0, 0.0000472, 0.00177, 0.0000492, 0, 0, 0.00523, 0.00429,
               0, 0.00002, 0.0000397, 0, 0, 0.106, 0.069, 0.169, 0.0122,
               0.663]).T
v5 = np.array([0.00000241, 0.00000113, 0.00000347, 0.0000118, 0.0000037,
               0.00000147, 0, 0.00000062, 0, 0.000934, 0, 0, 0, 0.000005,
               0.00254, 0, 0, 0.00000053, 0, 0, 0.000016, 0.00033,
               0.0000092, 0, 0, 0.00055, 0.00348, 0, 0.00053, 0.0000039,
               0, 0, 0.041, 0.0149, 0.0292, 0.00178, 0.0442]).T

b = np.array([0.00657, 0.00876, 0.00949, 0.1168, 0.0365, 0.01241, 0.000219,
              0.00292, 0, 0.657, 0.000146, 0.1095, 0.000876, 4.015, 8.76,
              16.79, 0.0002555, 0.00657, 0.0292, 0.001095, 0.0584, 3.066,
              0.01679, 0.0003285, 0, 5.11, 24.82, 0.0004015, 10.95, 0.0803,
              0, 0.00219, 204.4, 569.4, 365, 146, 2007.5]).T

# Build the matrix A (37 x 5)
A = np.column_stack((v1, v2, v3, v4, v5))

# Determine the number of variables (foods)
n = A.shape[1]

# Choose an objective. For example, you might want to minimize the sum of x.
# If you have a different notion of cost, adjust accordingly.
c = np.ones(n)  # Minimizing sum(x), subject to x >= 0 and Ax >= b

# The constraint Ax >= b can be transformed to -Ax <= -b.
A_ub = -A
b_ub = -b

# x must be nonnegative
bounds = [(0, None)] * n

# Run the linear programming solver
res = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method='highs')

if res.success:
    x_solution = res.x
    print("x solution:", x_solution)
    print("Ax:", A @ x_solution)
else:
    print("No feasible solution found:", res.message)

Points to Consider

• Objective Function: In the example above, the objective is to minimize the sum of x, but you can change it to any linear function depending on your goals.
• Feasibility: With 5 foods (variables) and 37 nutrient constraints, it’s possible that there might not be any x ≥ 0 that satisfies Ax ≥ b exactly. In such cases, you may need to relax or reconsider the constraints, or allow for some slack in the model.
• Alternative Approaches:
  • If you want to minimize a quadratic objective (like minimizing the squared error) under Ax ≥ b constraints, you could use a quadratic programming solver from packages like cvxopt or OSQP.
  • General constrained optimization can also be performed with scipy.optimize.minimize, though you’ll have to define the inequality constraints in the proper format.

In summary, while there isn’t a single, one-line function in SciPy that handles both the nonnegativity on x and the Ax ≥ b constraint, you can handle these constraints by formulating your problem as a linear (or quadratic) programming problem as shown above.