0

I am trying to solve an overdetermined linear system where the solution vector should sum to 1 and 0<=x<=1. I have tried using CVXPY to solve this, but sometimes the solution blatantly ignores the constraints. I also am having issues finding a good way to constrain the summation of x = 1. Any help would be great!

A = np.array([[7.958e+01 1.440e+00 6.971e+01 1.040e+00 3.733e+01 1.570e+01],
    [5.800e+00 3.200e-01 7.200e-01 1.020e+00 1.109e+01 4.980e+00],
    [1.460e+00 3.400e-01 6.600e-01 6.466e+01 4.050e+00 2.560e+00],
    [4.340e+00 5.256e+01 1.252e+01 1.427e+01 1.865e+01 4.231e+01],
    [9.200e-01 2.070e+00 3.610e+00 3.540e+00 5.100e+00 2.380e+00],
    [0.000e+00 6.000e-01 9.000e-02 2.100e-01 5.700e-01 1.280e+00],
    [1.450e+00 9.600e-02 2.800e-01 1.500e-01 1.478e+00 4.620e-01],
    [4.860e-01 2.700e-02 0.000e+00 3.000e-01 7.880e-01 2.330e-01]])

b = np.array([24.6875,  6.577,   2.2169, 59.9135,  3.3156,  1.0645,  0.7319,  
    0.3567])

# Define and solve the CVXPY problem.
x = cp.Variable(6)
objective = cp.Minimize(cp.sum_squares(A@x - b))
constraints = [0 <= x,  x <= 1, sum(x) == 1]
prob = cp.Problem(objective, constraints)
prob.solve()

#Print result.
print("\nThe optimal value is", prob.value)
print("The optimal x is")
print(x.value)
print("The norm of the residual is ", cp.norm(A @ x - b, p=2).value)

summ = sum(x.value)
print(summ)

>>>>
The optimal value is 370.4259696888074
The optimal x is
[-3.55262420e-18  7.02949975e-01  1.81476423e-01 -6.08113127e-16
-3.36614249e-16  1.15573602e-01]
The norm of the residual is  19.24645343144574
1.000000000000002

(This scenario represents the percent composition of a sample by input material. The input material (labeled material 1 through 4) has chemical composition measured in elements A through J. The sample composition is also measured in elements A through J. elemental compositions of materials and sample)

I have tried many different ways of summation, such as cp.sum_entries(x) and np.sum(x). I also do not know why my returned values for x are not between 0 and 1...

what's_python
  • 1
  • 1
  • If you're comfortable with Math behind least squares, you could set it up as a convex optimization problem (specifically quadratic programming, QP) with your constraints built into the problem. If I understand your question correctly, you need to augment a vanilla least squares cost function with non-negative and "sum to 1" constraints. Both of these are linear constraints and are easily incorporated into a standard QP. Please have a look at the package CVXOPT. – RDK Mar 09 '21 at 20:20
  • @RDK I looked into CVXOPT and began using CVXPY...but I am stuck again...see edits above – what's_python Mar 12 '21 at 04:41
  • Whats the problem? Looks good to me. Keep in mind, that the underlying solver (interior point) willl only approximate the solution (with a-priori defined tolerances; probably around 10^-10 or similar). Your simplex-constraint (sum=1) is violated by `0.000000000000002` and the biggest violation of the [0, 1] bounds is < 10^-16. Looks good to me. Remark: `cp.sum(x)` is now the canonical way (sum_entries has been renamed). `sum(x)` will work too. `np.sum(x)` should not. A simplex-algorithm (supporting quadratic objectives) would allow more nice (basic feasible) solutions, but this is pure luxury. – sascha Mar 12 '21 at 14:08
  • To demonstrate: change your prints to `print("The optimal x is: ", x.value.round(5))` and `print('sum: ', x.value.sum().round(5))` – sascha Mar 12 '21 at 14:12
  • Related: https://stackoverflow.com/q/33385898/6151828 – Roger Vadim Sep 29 '22 at 07:47

0 Answers0