As pointed out by @SevC_10 in his answer, you are missing dx parameter.

I like to show case the use of sympy for derivation operations, I find it much easier in many cases.

Copyimport sympy
import numpy as np

x = sympy.Symbol('x')

f = sympy.exp(x) # my function e^x
df = f.diff() # y' of the function = e^x

f_lambda = sympy.lambdify(x, f, 'numpy')
df_lambda = sympy.lambdify(x, yprime, 'numpy') # use lambdify

print(f_lambda(np.ones(5)))

# array([2.71828183, 2.71828183, 2.71828183, 2.71828183, 2.71828183])

print(df_lambda(np.ones(5)))

# array([2.71828183, 2.71828183, 2.71828183, 2.71828183, 2.71828183])

print(f_lambda(np.zeros(5)))

# array([1., 1., 1., 1., 1.])

print(df_lambda(np.zeros(5)))

# array([1., 1., 1., 1., 1.])


print(f_lambda(np.array([0, 1, 2, 3, 4])))
# array([ 1.        ,  2.71828183,  7.3890561 , 20.08553692, 54.59815003])

print(df_lambda(np.array([0, 1, 2, 3, 4])))
# array([ 1.        ,  2.71828183,  7.3890561 , 20.08553692, 54.59815003])

As the documentation for derivative says:

derivative(func, x0, dx=1.0, n=1, args=(), order=3)
    Find the n-th derivative of a function at point x0.

    Given a function, use a central difference formula with spacing `dx` to
    compute the n-th derivative at `x0`.

You didn't specify dx, so it used the default of 1, which is too large here. For example:

In [1]: from scipy.misc import derivative

In [2]: derivative(lambda x: 3*x**4 + 2*x**3 - 10*x**2 + 15*x - 2, -2, dx=1)
Out[2]: -39.0

In [3]: derivative(lambda x: 3*x**4 + 2*x**3 - 10*x**2 + 15*x - 2, -2, dx=0.5)
Out[3]: -22.5

In [4]: derivative(lambda x: 3*x**4 + 2*x**3 - 10*x**2 + 15*x - 2, -2, dx=0.1)
Out[4]: -17.220000000000084

In [5]: derivative(lambda x: 3*x**4 + 2*x**3 - 10*x**2 + 15*x - 2, -2, dx=0.01)
Out[5]: -17.0022000000003

In [6]: derivative(lambda x: 3*x**4 + 2*x**3 - 10*x**2 + 15*x - 2, -2, dx=1e-5)
Out[6]: -17.000000001843318

Alternatively, you could increase the order:

In [7]: derivative(lambda x: 3*x**4 + 2*x**3 - 10*x**2 + 15*x - 2, -2, dx=1, order=5)
Out[7]: -17.0

Taking numerical derivatives is always a little troublesome.

The Savitzky-Golay filter uses a constant delta (the spacing of the samples,) and the default value of the delta in the filter implementation is 1, according to https://docs.scipy.org/doc/scipy-0.16.1/reference/generated/scipy.signal.savgol_filter.html. Your data set has irregular deltas, not 1, so the result from the Savgol filter is incorrect.

I would write a simple wrapper, something along the lines of

def partial_derivative(func, var=0, point=[]):
    args = point[:]
    def wraps(x):
        args[var] = x
        return func(*args)
    return derivative(wraps, point[var], dx = 1e-6)

Demo:

>>> partial_derivative(foo, 0, [3,1])
6.0000000008386678
>>> partial_derivative(foo, 1, [3,1])
2.9999999995311555

You have four options

1. Finite Differences
2. Automatic Derivatives
3. Symbolic Differentiation
4. Compute derivatives by hand.

Finite differences require no external tools but are prone to numerical error and, if you're in a multivariate situation, can take a while.

Symbolic differentiation is ideal if your problem is simple enough. Symbolic methods are getting quite robust these days. SymPy is an excellent project for this that integrates well with NumPy. Look at the autowrap or lambdify functions or check out Jensen's blogpost about a similar question.

Automatic derivatives are very cool, aren't prone to numeric errors, but do require some additional libraries (google for this, there are a few good options). This is the most robust but also the most sophisticated/difficult to set up choice. If you're fine restricting yourself to numpy syntax then Theano might be a good choice.

Here is an example using SymPy

In [1]: from sympy import *
In [2]: import numpy as np
In [3]: x = Symbol('x')
In [4]: y = x**2 + 1
In [5]: yprime = y.diff(x)
In [6]: yprime
Out[6]: 2⋅x

In [7]: f = lambdify(x, yprime, 'numpy')
In [8]: f(np.ones(5))
Out[8]: [ 2.  2.  2.  2.  2.]

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).

FFT returns a complex array that has the same dimensions as the input array. The output array is ordered as follows:

1. Element 0 contains the zero frequency component, F0.

2. The array element F1 contains the smallest, nonzero positive frequency, which is equal to 1/(Ni Ti), where Ni is the number of elements and Ti is the sampling interval.

3. F2 corresponds to a frequency of 2/(Ni Ti).

4. Negative frequencies are stored in the reverse order of positive frequencies, ranging from the highest to lowest negative frequencies.

5. For an even number of points, the frequencies corresponding to the returned complex values are: 0, 1/(NiTi), 2/(NiTi), ..., (Ni/2–1)/(NiTi), 1/(2Ti), –(Ni/2–1)/(NiTi), ..., –1/(NiTi) where 1/(2Ti) is the Nyquist critical frequency.

6. For an odd number of points, the frequencies corresponding to the returned complex values are: 0, 1/(NiTi), 2/(NiTi), ..., (Ni–1)/2)/(NiTi), –(Ni–1)/2)/(NiTi), ..., –1/(NiTi)

Using this information we can construct the proper vector of frequencies that should be used for calculating the derivative. Below is a piece of self-explanatory Python code that does it all correctly. Note that the factor 2$\pi$N cancels out due to normalization of FFT.

from scipy.fftpack import fft, ifft, dct, idct, dst, idst, fftshift, fftfreq
from numpy import linspace, zeros, array, pi, sin, cos, exp, arange
import matplotlib.pyplot as plt


N = 100
x = 2*pi*arange(0,N,1)/N #-open-periodic domain                                                   

dx = x[1]-x[0]
y = sin(2*x)+cos(5*x)
dydx = 2*cos(2*x)-5*sin(5*x)


k2=zeros(N)

if ((N%2)==0):
    #-even number                                                                                   
    for i in range(1,N//2):
        k2[i]=i
        k2[N-i]=-i
else:
    #-odd number                                                                                    
    for i in range(1,(N-1)//2):
        k2[i]=i
        k2[N-i]=-i

dydx1 = ifft(1j*k2*fft(y))

plt.plot(x,dydx,'b',label='Exact value')
plt.plot(x,dydx1, color='r', linestyle='--', label='Derivative by FFT')
plt.legend()
plt.show()