Hello, I have written a very simple code that computes the gradient by finite differences of any general function. Keeping the same idea, I would like modify the code using numpy to make it faster. Any ideas? Thanks. def grad_finite_dif(self,x,user_data = None): assert len(x) == self.number_variables points=[] for j in range(self.number_variables): points.append(x.copy()) points[len(points)-1][j]=points[len(points)-1][j]+0.0000001 delta_f = [] counter=0 for j in range(self.number_variables): delta_f.append((self.eval(points[counter])-self.eval(x))/0.0000001) counter = counter + 1 return array(delta_f) -- View this message in context: http://old.nabble.com/Improvement-of-performance-tp28452458p28452458.html Sent from the Numpy-discussion mailing list archive at Nabble.com.
On Tue, May 4, 2010 at 4:06 PM, gerardob <gberbeglia@gmail.com> wrote:
Hello, I have written a very simple code that computes the gradient by finite differences of any general function. Keeping the same idea, I would like modify the code using numpy to make it faster. Any ideas? Thanks.
def grad_finite_dif(self,x,user_data = None): assert len(x) == self.number_variables points=[] for j in range(self.number_variables): points.append(x.copy()) points[len(points)-1][j]=points[len(points)-1][j]+0.0000001 delta_f = [] counter=0 for j in range(self.number_variables): delta_f.append((self.eval(points[counter])-self.eval(x))/0.0000001)
it looks like your are evaluating the same point several times self.eval(x)
counter = counter + 1 return array(delta_f)
That's what I used as a pattern for a gradient function #from scipy.optimize def approx_fprime(xk,f,epsilon,*args): f0 = f(*((xk,)+args)) grad = np.zeros((len(xk),), float) ei = np.zeros((len(xk),), float) for k in range(len(xk)): ei[k] = epsilon grad[k] = (f(*((xk+ei,)+args)) - f0)/epsilon ei[k] = 0.0 return grad Josef
-- View this message in context: http://old.nabble.com/Improvement-of-performance-tp28452458p28452458.html Sent from the Numpy-discussion mailing list archive at Nabble.com.
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
If your x data are equispaced I would do something like this def derive( func, x): """ Approximate the first derivative of function func at points x. """ # compute the values of y = func(x) y = func(x) # compute the step dx = x[1] - x[0] # kernel array for second order accuracy centered first derivative kc = np.array([-1.0, +0.0, +1.0]) / 2 / dx # kernel array for second order accuracy left and right first derivative kl = np.array([-3.0, +4.0, -1.0]) / 2 / dx kr = np.array([+1.0, -4.0, +3.0]) / 2 / dx # correlate it with the original array, # note that only the valid computation are performed derivs_c = np.correlate( y, kc, mode='valid' ) derivs_r = np.correlate( y[-3:], kr, mode='valid' ) derivs_l = np.correlate( y[:+3], kl, mode='valid' ) return np.r_[derivs_l, derivs_c, derivs_r] This is actually quite fast: on my machine (1.7GHz) i have this.
:x = np.linspace(0,2*np.pi, 1e6) :func = lambda x: np.sin(x) :timeit derive(func, x) 10 loops, best of 3: 177 ms per loop
I'm curious if someone comes up with something faster. Regards, Davide On 4 May 2010 22:17, <josef.pktd@gmail.com> wrote:
On Tue, May 4, 2010 at 4:06 PM, gerardob <gberbeglia@gmail.com> wrote:
Hello, I have written a very simple code that computes the gradient by
finite
differences of any general function. Keeping the same idea, I would like modify the code using numpy to make it faster. Any ideas? Thanks.
def grad_finite_dif(self,x,user_data = None): assert len(x) == self.number_variables points=[] for j in range(self.number_variables): points.append(x.copy())
points[len(points)-1][j]=points[len(points)-1][j]+0.0000001
delta_f = [] counter=0 for j in range(self.number_variables):
delta_f.append((self.eval(points[counter])-self.eval(x))/0.0000001)
it looks like your are evaluating the same point several times self.eval(x)
counter = counter + 1 return array(delta_f)
That's what I used as a pattern for a gradient function
#from scipy.optimize def approx_fprime(xk,f,epsilon,*args): f0 = f(*((xk,)+args)) grad = np.zeros((len(xk),), float) ei = np.zeros((len(xk),), float) for k in range(len(xk)): ei[k] = epsilon grad[k] = (f(*((xk+ei,)+args)) - f0)/epsilon ei[k] = 0.0 return grad
Josef
-- View this message in context: http://old.nabble.com/Improvement-of-performance-tp28452458p28452458.html Sent from the Numpy-discussion mailing list archive at Nabble.com.
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
playing devil's advocate I'd say use Algorithmic Differentiation instead of finite differences ;) that would probably speed things up quite a lot. On Tue, May 4, 2010 at 11:36 PM, Davide Lasagna <lasagnadavide@gmail.com> wrote:
If your x data are equispaced I would do something like this def derive( func, x): """ Approximate the first derivative of function func at points x. """ # compute the values of y = func(x) y = func(x) # compute the step dx = x[1] - x[0] # kernel array for second order accuracy centered first derivative kc = np.array([-1.0, +0.0, +1.0]) / 2 / dx # kernel array for second order accuracy left and right first derivative kl = np.array([-3.0, +4.0, -1.0]) / 2 / dx kr = np.array([+1.0, -4.0, +3.0]) / 2 / dx # correlate it with the original array, # note that only the valid computation are performed derivs_c = np.correlate( y, kc, mode='valid' ) derivs_r = np.correlate( y[-3:], kr, mode='valid' ) derivs_l = np.correlate( y[:+3], kl, mode='valid' ) return np.r_[derivs_l, derivs_c, derivs_r] This is actually quite fast: on my machine (1.7GHz) i have this.
:x = np.linspace(0,2*np.pi, 1e6) :func = lambda x: np.sin(x) :timeit derive(func, x) 10 loops, best of 3: 177 ms per loop I'm curious if someone comes up with something faster.
Regards, Davide
On 4 May 2010 22:17, <josef.pktd@gmail.com> wrote:
On Tue, May 4, 2010 at 4:06 PM, gerardob <gberbeglia@gmail.com> wrote:
Hello, I have written a very simple code that computes the gradient by finite differences of any general function. Keeping the same idea, I would like modify the code using numpy to make it faster. Any ideas? Thanks.
def grad_finite_dif(self,x,user_data = None): assert len(x) == self.number_variables points=[] for j in range(self.number_variables): points.append(x.copy())
points[len(points)-1][j]=points[len(points)-1][j]+0.0000001 delta_f = [] counter=0 for j in range(self.number_variables):
delta_f.append((self.eval(points[counter])-self.eval(x))/0.0000001)
it looks like your are evaluating the same point several times self.eval(x)
counter = counter + 1 return array(delta_f)
That's what I used as a pattern for a gradient function
#from scipy.optimize def approx_fprime(xk,f,epsilon,*args): f0 = f(*((xk,)+args)) grad = np.zeros((len(xk),), float) ei = np.zeros((len(xk),), float) for k in range(len(xk)): ei[k] = epsilon grad[k] = (f(*((xk+ei,)+args)) - f0)/epsilon ei[k] = 0.0 return grad
Josef
-- View this message in context: http://old.nabble.com/Improvement-of-performance-tp28452458p28452458.html Sent from the Numpy-discussion mailing list archive at Nabble.com.
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
On Tue, May 4, 2010 at 2:57 PM, Sebastian Walter <sebastian.walter@gmail.com> wrote:
playing devil's advocate I'd say use Algorithmic Differentiation instead of finite differences ;) that would probably speed things up quite a lot.
I would suggest that too, but aside from FuncDesigner[0] (reference in the end), I couldn't find any Automatic Differentiation tool that was easy to install for Python. To stay with simple solutions, I think that the "complex step" approximation gives you a very good compromise between ease of use, performance and accuracy. Here is an implementation (and you can take ideas from it to get rid of your "for" loops in your original code) import numpy as np def complex_step_grad(f, x, h=1.0e-20): dim = np.size(x) increments = np.identity(dim) * 1j * h partials = [f(x+ih).imag / h for ih in increments] return np.array(partials) **Warning**: you must convert your original real-valued function f: R^n -> R to the corresponding complex function f: C^n -> C. Use functions from the 'cmath' module. I strongly suggest that you take a look at the AutoDiff website[1] and at some references about the complex step [2][3] (or just Google "complex step" and "differentiation", both on normal Google and Google Scholar). [0] http://openopt.org/FuncDesigner [1] http://www.autodiff.org/ [2] http://doi.acm.org/10.1145/838250.838251 [3] http://mdolab.utias.utoronto.ca/resources/complex-step -- Guilherme P. de Freitas http://www.gpfreitas.com
I forgot to mention one thing: if you are doing optimization, a good solution is a modeling package like AMPL (or GAMS or AIMMS, but I only know AMPL, so I will restrict my attention to it). AMPL has a natural modeling language and provides you with automatic differentiation. It's not free, but there are trial licenses (60 days) and student a student edition (unlimited time, maximum of 300 variables). It is hooked to many great solvers, including KNITRO (commercial) and IPOPT (free), both for nonlinear programs. http://www.ampl.com And as this is a Python list, there is NLPy. I never tried it, but it seems to allow you to read model files written in AMPL, use AMPL's automatic differentiation capabilities, and still roll your own optimization algorithm, all in Python. It looks like it requires some other packages aside from Python, NumPy and AMPL, like a sparse linear solver and some other things. Worth taking a look. http://nlpy.sourceforge.net/how.html Best, Guilherme On Tue, May 4, 2010 at 5:23 PM, Guilherme P. de Freitas <guilherme@gpfreitas.com> wrote:
On Tue, May 4, 2010 at 2:57 PM, Sebastian Walter <sebastian.walter@gmail.com> wrote:
playing devil's advocate I'd say use Algorithmic Differentiation instead of finite differences ;) that would probably speed things up quite a lot.
I would suggest that too, but aside from FuncDesigner[0] (reference in the end), I couldn't find any Automatic Differentiation tool that was easy to install for Python.
To stay with simple solutions, I think that the "complex step" approximation gives you a very good compromise between ease of use, performance and accuracy. Here is an implementation (and you can take ideas from it to get rid of your "for" loops in your original code)
import numpy as np
def complex_step_grad(f, x, h=1.0e-20): dim = np.size(x) increments = np.identity(dim) * 1j * h partials = [f(x+ih).imag / h for ih in increments] return np.array(partials)
**Warning**: you must convert your original real-valued function f: R^n -> R to the corresponding complex function f: C^n -> C. Use functions from the 'cmath' module.
I strongly suggest that you take a look at the AutoDiff website[1] and at some references about the complex step [2][3] (or just Google "complex step" and "differentiation", both on normal Google and Google Scholar).
[0] http://openopt.org/FuncDesigner [1] http://www.autodiff.org/ [2] http://doi.acm.org/10.1145/838250.838251 [3] http://mdolab.utias.utoronto.ca/resources/complex-step
-- Guilherme P. de Freitas http://www.gpfreitas.com
-- Guilherme P. de Freitas http://www.gpfreitas.com
On Tue, May 4, 2010 at 8:23 PM, Guilherme P. de Freitas <guilherme@gpfreitas.com> wrote:
On Tue, May 4, 2010 at 2:57 PM, Sebastian Walter <sebastian.walter@gmail.com> wrote:
playing devil's advocate I'd say use Algorithmic Differentiation instead of finite differences ;) that would probably speed things up quite a lot.
I would suggest that too, but aside from FuncDesigner[0] (reference in the end), I couldn't find any Automatic Differentiation tool that was easy to install for Python.
To stay with simple solutions, I think that the "complex step" approximation gives you a very good compromise between ease of use, performance and accuracy. Here is an implementation (and you can take ideas from it to get rid of your "for" loops in your original code)
import numpy as np
def complex_step_grad(f, x, h=1.0e-20): dim = np.size(x) increments = np.identity(dim) * 1j * h partials = [f(x+ih).imag / h for ih in increments] return np.array(partials)
**Warning**: you must convert your original real-valued function f: R^n -> R to the corresponding complex function f: C^n -> C. Use functions from the 'cmath' module.
Interesting idea I tried it with some silly function that has special.gammaln and np.dot in it, and it seems to work without adjustments to the code. The precision is much better than simple forward differentiation, which I have only with 1e-5 accuracy. simple timing 40% slower than simple forward differentiation 50%-80% faster that forward and backward differentiation In [2] I didn't see anything about higher derivatives, so to get the Hessian I still had to do a finite difference (Jacobian) on the complex_step_grad. Even then the results look pretty good. Another recommendation especially to check whether the results are correct: http://pypi.python.org/pypi/Numdifftools/ is pure python, (optionally adaptive) finite difference method for Gradient, Jacobian and Hessian. And related: I was surprised that sympy knows what the analytical derivative of the gamma function is, and that the function is available in scipy.special, so one less reason not to use analytical derivatives. I will include complex_step_grad in the toolbox for Maximum Likelihood Estimation. Thanks for the information, Josef
I strongly suggest that you take a look at the AutoDiff website[1] and at some references about the complex step [2][3] (or just Google "complex step" and "differentiation", both on normal Google and Google Scholar).
[0] http://openopt.org/FuncDesigner [1] http://www.autodiff.org/ [2] http://doi.acm.org/10.1145/838250.838251 [3] http://mdolab.utias.utoronto.ca/resources/complex-step
-- Guilherme P. de Freitas http://www.gpfreitas.com _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
On Tue, May 4, 2010 at 9:23 PM, <josef.pktd@gmail.com> wrote:
In [2] I didn't see anything about higher derivatives, so to get the Hessian I still had to do a finite difference (Jacobian) on the complex_step_grad. Even then the results look pretty good.
Yes, the traditional complex step does not solve the second derivatives problem. I think there are papers that try to address that (I don't know the literature though). But I have seen a paper [0] that extends the notion of complex number to multi-complex numbers, and then using the algebra of those multi-complex numbers the authors obtain higher order derivatives. It lacks the convenience of the complex step (in the sense that complex numbers are already implemented), but it's something to keep in mind. [0] http://soliton.ae.gatech.edu/people/rrussell/FinalPublications/ConferencePap... But if we were to go that way (defining new numbers), maybe (I'm no expert in this area, some of the proponents of AD via dual numbers say it would be a better idea) it would be better to define dual numbers. They are like complex numbers, but their "imaginary component" (I don't think it's called that) d has the property that d^2 = 0. Using these numbers, one can obtain first and higher order derivatives "automatically" [1][2] (forward mode only, see refs.) [1] http://en.wikipedia.org/wiki/Automatic_differentiation#Automatic_differentia... [2] http://conal.net/papers/beautiful-differentiation/ (recent paper with references, related blog posts and video of a talk) There is a group of people in the Haskell community [3][4] that are working on Automatic Differentiation via dual numbers (that gives you the "forward mode" only, see refs.). It looks really interesting. I saw somewhere that one could write external modules in Haskell for Python... [3] http://www.haskell.org/haskellwiki/Automatic_Differentiation [4] http://hackage.haskell.org/packages/archive/fad/1.0/doc/html/Numeric-FAD.htm... Just throwing ideas out there. I found it really neat, and potentially very useful. Now, as for the reverse mode of AD, there seems to be no shortcut, and it would be really nice, as it seems that the computational complexity of the reverse mode is lower than the forward mode. There are libraries, like TAPENADE [5] and ADIFOR [6] that do it, though. They do it by source code transformation. The TAPENADE tool even accepts C or Fortran code via the web and returns a differentiated code (worth playing with). Really neat. [5] http://www-sop.inria.fr/tropics/ [6] http://www.mcs.anl.gov/research/projects/adifor Ok, now all we need is easy access to these tools from Python, and we are set! :) -- Guilherme P. de Freitas http://www.gpfreitas.com
Just to make this thread more useful for someone interested in these topics, this seems to be "the book" on automatic differentiation (it's one of the references in the autodiff website) "Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation" (2nd ed.), by Andreas Griewank and Andrea Walther. SIAM, 2008. Link: http://www.ec-securehost.com/SIAM/OT105.html
participants (5)
-
Davide Lasagna -
gerardob -
Guilherme P. de Freitas -
josef.pktd＠gmail.com -
Sebastian Walter