How to define the derivative for Scipy.Optimize.Minimize

I am trying to use scipy.optimize.minimize to minimise a quadratic objective function: $f(x) =x^\top Q x$. As a start, I have successfully implemented this using the built-in Nelder-Mead Simplex algorithm, by defining a function:

def objective(x):
Q = np.asmatrix(DF.cov())       # Covariance matrix
x = np.asmatrix(x)
return x.transpose() * Q * x


which is passed to SciPy's minimise function:

minimize(objective, x0, method='Nelder-Mead',options={'xtol': 1e-6, 'disp': True})


Now, following the documentation at the above link, there are a number of other optimisation routines available, that I would like to try. Each of these require the calculation of the function derivative, $\nabla f(x)$, which must be written inside a python function similar to the above, and some require the Hessian $\nabla^2f(x)$. My initial attempt is below:

def der_objective(x):
Q = np.asmatrix(DF.cov())       # Covariance matrix
x = np.asmatrix(x).reshape(len(Q),1)
return 2 * Q * x


but this obviously returns the wrong dimensions. I'm not entirely sure how SciPy expects the result, and couldn't work it out from the Rosenbrock example in the tutorial here.

Could you please explain how the rosen_der example in the SciPy tutorial works (e.g. Why does it take the derivative with respect to a given $x_j$? Why does it return a vector of the same shape as $x$?) and how I should design my der_objective(x) function to achieve the same thing?

• As a side note, is there a factor of 2 missing from one of your functions? – origimbo Mar 19 '18 at 20:37
• Yes, that's right it should be 2 * Q * x. I will update the question to show this – Zac Mar 19 '18 at 23:58

If you ignore the mathematical formulae in the tutorial you link to, and just look at the call itself,

res = minimize(rosen, x0, method='BFGS', jac=rosen_der, ... options={'disp': True})

There are two python functions defined. One is the scalar functional, $f(\vec{x})$,

def rosen(x): ... """The Rosenbrock function""" ... return sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0)

and the other is the vector gradient function, $$\nabla f(\vec{x}) = \left(\left.\frac{\partial f}{\partial x_0}\right|_{\vec{x}},\ldots,\left.\frac{\partial f}{\partial x_{n-1}}\right|_{\vec{x}}\right)^T,$$

which is automatically the same size as $\vec{x}$.

def rosen_der(x): ... xm = x[1:-1] ... xm_m1 = x[:-2] ... xm_p1 = x[2:] ... der = np.zeros_like(x) ... der[1:-1] = 200*(xm-xm_m1**2) - 400*(xm_p1 - xm**2)*xm - 2*(1-xm) ... der[0] = -400*x[0]*(x[1]-x[0]**2) - 2*(1-x[0]) ... der[-1] = 200*(x[-1]-x[-2]**2) ... return der

In terms of the code, this is deliberately written in the same form as $\vec{x}$ so that the iterative minimisation routine can form objects like $$\vec{x}^{(n+1}=\vec{x}^n-\alpha \nabla f(\vec{x}^n)$$ so that the variable can move "downhill" towards a local minimum through successive iterations of the vector $\vec{x}$ (note just using the last answer is the method "steepest descent", which nobody would ever use choose to use outside of a classroom). Hopefully this isn't too much of a surprise to you. If it is, you might want to read up a bit more on the theory of gradient based optimisation before getting back to your coding.

In terms of your actual code, it should be enough to use basic arrays with numpy.dot to perform matrix multiplication,rather than writing everything as matrices (which numpy makes 2d), or if you really want to keep this style in your functions, use numpy.asarray and numpy.ravel when returning the result. At a worst case, you need to double check that your matrix $Q$ is $N\times N$ and that you make $x$ a column vector (which as I said before asmatrix doesn't do).

• Thanks. I actually have W = np.asmatrix(W).reshape(len(Mus),1) in my code ... which seems to get me a vector of the right shape, but didn't want to obscure the real question. I had been confused why, for each $x_i$ value in my vector, the vector gradient $\nabla f(x_i)$ appeared to be scalar... it seemed counter-intuitive for a scalar $x$ and a vector gradient to have the same dimension array. But I think now I realise that the $x$ vector is in fact one coordinate in an n-dimensional space.. and the gradient at this point can also be represented by an n-dimensional vector. Is this right? – Zac Mar 19 '18 at 23:55
• Coordinate isn't the best word to use (point might be appropriate), but that's the basic idea. Note that the computer science answer (a vector is a vector is a vector, i.e an array with a single axis) is somewhat simpler but more ambiguous than the functional analysis description. – origimbo Mar 20 '18 at 0:40
• Indeed; I struggle with this a little, coming from a Physics background but now studying scientific computing. I suppose $x$ is a coordinate vector in an n-d space. $f(x)$ is a map from $\mathcal{R}^N \rightarrow \mathcal{R}$. And $\nabla f(x)$ is a direction vector also in n-d space. Both $x$ and $\nabla f(x)$ are stored computationally in arrays of shape (N,1). But how does that code for rosen_der() produce an array of $\frac{\partial{f(x)}}{\partial{x_i}}$? for $i = 0,1,2..N$ ?? At least I think that's what the code (should) do, but can't see how to do it for my der_objective() function. – Zac Mar 20 '18 at 17:09
• That is what the code is doing. The return value is precisely that 1d array of length N. It's also what your routine is nearly doing. I strongly suspect your issues are down to using 2d matrices of size (N,1) when routines want the older vanilla numpy arrays of length N. These things don't automatically transform in the context lots of basic routines like numpy.dot. – origimbo Mar 20 '18 at 19:47