2

Given arrays a, b, and c:

import numpy as np

a = np.array([100, 200, 300])
b = np.array([[1, 0, 0],
              [1, 0, 1],
              [0, 1, 1],
              [1, 1, 1]])
c = np.array([150, 300, 500, 650])

I'd like to optimize a such that each value minimizes the sum of the absolute difference defined in c_prime.

c_prime = c - np.sum(a*b, axis=1)
print(c_prime)
print(np.abs(c_prime).sum())

[  50 -100    0   50]
200

Manually... by changing the first element in a, c_prime starts to achieve the desired result.

a = np.array([150, 200, 300])

c_prime = c - np.sum(a*b, axis=1)
print(c_prime)
print(np.abs(c_prime).sum())

[   0 -150    0    0]
150

Now, my question, embarrassingly, is how can I achieve the desired result? I've tried scipy.optimize.minimize, but it's obvious this code misses the mark and the function may be conceptually incorrect altogether.

def f(x, b, c):
    return np.abs(c - np.sum(x*b, axis=1)).sum()

x0 = a
minimize(f, x0, args=(b,c))

      fun: 200.0
 hess_inv: array([[1, 0, 0],
       [0, 1, 0],
       [0, 0, 1]])
      jac: array([-1.,  0.,  1.])
  message: 'Desired error not necessarily achieved due to precision loss.'
     nfev: 327
      nit: 0
     njev: 63
   status: 2
  success: False
        x: array([100., 200., 300.])

Given the improved results from manually setting a[0] to 150 above, why do these results return a non-optimal x?

Improve this question
4
  • 1
    You dropped a sum call. Commented Sep 26, 2018 at 17:40
  • Thank you, the question has been updated. Commented Sep 26, 2018 at 17:48
  • Just an aside: I'm pretty sure there is a strict mathematical solution to this. Basically you are trying to find a such that sum(abs(B*a - c)) is minimized (here both a and c are 1-column vectors). That just seems strictly solvable to me. Maybe also ask the math stack exchange. Commented Sep 26, 2018 at 17:48
  • a=np.linalg.lstsq(b,c)[0][None,:] solves the matrix equation b*a-c=0 (by least squares), but it may be close enough to abs Commented Sep 26, 2018 at 18:02

1 Answer 1

Reset to default
3

The problem here is that your objective function is non-differentiable. SciPy is defaulting to BFGS optimization, which requires first derivatives of the objective function to exist.

I can think of 3 primary approaches to getting around this: use derivative-free optimization, use a differentiable approximation to your objective function, or transform the absolute values into constraints.

Almost every optimization method in scipy.optimize.minimize requires a differentiable objective function. A few don't, but even then, there's no guarantee that they'll find the minimum.

For example, specifying method='Nelder-Mead' results in successful optimization and a result of x: array([ 149.99998103, 349.99999851, 150.00000599]) in my test run, but as Paul Panzer points out in the comments, starting at x0=[1, 1, 1] results in convergence to a non-minimum. Nelder-Mead just does that sometimes; even with a differentiable objective function, it can converge to a non-stationary point.

A differentiable approximation to the objective function is easy and provides much better convergence properties, at the cost of a bit of error. For example, replacing np.abs with

def pseudoabs(x):
    return (x**2+0.1)**0.5

results in convergence to x: array([ 150.00000001, 350.00000011, 150.00000039]) with the default BFGS solver.

As for transforming the absolute values away, your problem is almost a standard linear programming problem, but with absolute values in the objective function. By introducing extra variables, it's possible to convert an absolute value into two new linear constraints. The idea is to replace an |x| term in the objective by an x' term subject to constraints x' >= x and x' >= -x.

Doing so would let you solve your problem with a standard linear programming solver like scipy.optimize.linprog, or with scipy.optimize.minimize if you prefer it.

Aside from that, there may be other algorithms for optimization problems of this form. I tried googling "l1 optimization", in the vein of "least squares optimization", but that turned up results for minimizing the l1 norm of the solution vector, which is the wrong vector for this case.

Improve this answer
Sign up to request clarification or add additional context in comments.

3 Comments

Doesn't seem to work for x0=[1, 1, 1]. Why not provide some guidance on how to get a good enough initial guess?
@PaulPanzer: After some more research, it looks like that's just a thing that sometimes happens with Nelder-Mead. I've expanded the answer to give more options.
Thanks for the descriptions of the possible solutions, they helped my understanding beyond the question as composed.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.