Mailman 3 July 2008 - SciPy-Dev

SLSQP Constrained Optimizer Status
by Rob Falck 25 Nov '13

25 Nov '13

I'm currently implementing the Sequential Least Squares Quadratic Programming (SLSQP) optimizer, by Dieter Kraft, for use in Scipy. The Fortran code being wrapped with F2PY is here: http://www.netlib.org/toms/733 (its use within Scipy has been cleared) SLSQP provides for bounds on the independent variables, as well as equality and inequality constraint functions, which is a capability that doesn't exist in scipy.optimize. Currently the code works, although the constraint normals are being generated by approximation. I'm working on a way to pass these in. I think the most elegant way will be a single function that returns the matrix of constraint normals. For a demonstration of what the code can do, here is an optimization of f(x,y) = 2xy + 2x - x**2 - 2y**2 Example 14.2 in Chapra & Canale gives the maximum as x=2.0, y=1.0. The unbounded optimization tests find this solution. As expected, its faster when derivatives are provided rather than approximated. Unbounded optimization. Derivatives approximated. Elapsed time: 1.45792961121 ms Results [[1.9999999515712266, 0.99999996181577444], -1.9999999999999984, 4, 0, 'Optimization terminated successfully.'] Unbounded optimization. Derivatives provided. Elapsed time: 1.03211402893 ms Results [[1.9999999515712266, 0.99999996181577444], -1.9999999999999984, 4, 0, 'Optimization terminated successfully.'] The following example uses an equality constraint to find the optimal when x=y. Bound optimization. Derivatives approximated. Elapsed time: 1.384973526 ms Results [[0.99999996845920858, 0.99999996845920858], -0.99999999999999889, 4, 0, 'Optimization terminated successfully.'] I've tried to conform to the syntax used by the other optimizers in scipy.optimize. The function definition and doc string are below. If anyone is interested in testing it out, let me know. def fmin_slsqp( func, x0 , eqcons=[], f_eqcons=None, ieqcons=[], f_ieqcons=None, bounds = [], fprime = None, fprime_cons=None,args = (), iter = 100, acc = 1.0E-6, iprint = 1, full_output = 0, epsilon = _epsilon ): """ Minimize a function using Sequential Least SQuares Programming Description: Python interface function for the SLSQP Optimization subroutine originally implemented by Dieter Kraft. Inputs: func - Objective function (in the form func(x, *args)) x0 - Initial guess for the independent variable(s). eqcons - A list of functions of length n such that eqcons[j](x0,*args) == 0.0 in a successfully optimized problem f_eqcons - A function of the form f_eqcons(x, *args) that returns an array in which each element must equal 0.0 in a successfully optimized problem. If f_eqcons is specified, eqcons is ignored. ieqcons - A list of functions of length n such that ieqcons[j](x0,*args) >= 0.0 in a successfully optimized problem f_ieqcons - A function of the form f_ieqcons(x0, *args) that returns an array in which each element must be greater or equal to 0.0 in a successfully optimized problem. If f_ieqcons is specified, ieqcons is ignored. bounds - A list of tuples specifying the lower and upper bound for each independent variable [ (xl0, xu0), (xl1, xu1), ...] fprime - A function that evaluates the partial derivatives of func fprime_cons - A function of the form f(x, *args) that returns the m by n array of constraint normals. If not provided, the normals will be approximated. Equality constraint normals precede inequality constraint normals. args - A sequence of additional arguments passed to func and fprime iter - The maximum number of iterations (int) acc - Requested accuracy (float) iprint - The verbosity of fmin_slsqp. iprint <= 0 : Silent operation iprint == 1 : Print summary upon completion (default) iprint >= 2 : Print status of each iterate and summary full_output - If 0, return only the minimizer of func (default). Otherwise, output final objective function and summary information. epsilon - The step size for finite-difference derivative estimates. Outputs: ( x, { fx, gnorm, its, imode, smode }) x - The final minimizer of func. fx - The final value of the objective function. its - The number of iterations. imode - The exit mode from the optimizer, as an integer. smode - A string describing the exit mode from the optimizer. Exit modes are defined as follows: -1 : Gradient evaluation required (g & a) 0 : Optimization terminated successfully. 1 : Function evaluation required (f & c) 2 : Number of equality constraints larger than number of independent variables 3 : More than 3*n iterations in LSQ subproblem 4 : Inequality constraints incompatible 5 : Singular matrix E in LSQ subproblem 6 : Singular matrix C in LSQ subproblem 7 : Rank-deficient equality constraint subproblem HFTI 8 : Positive directional derivative for linesearch 9 : Iteration limit exceeded -- - Rob Falck

4 16

Question about Kaiser Implementation in firwin
by Buehler, Eric (AGRE) 04 Oct '08

04 Oct '08

Hello, I have been using the signal.firwin FIR design tool and I have a question about the implementation. Scipy-0.4.8 Numpy-0.9.6 I am generating a kaiser window with the following parameters: N = 128 Cutoff = ~4e7 width = .3 All of the special kaiser generation parameters in the function follow Oppenheim and Schafer well. Even the sinc function for linear phase is follows exactly. However, the last line of the firwin function is causing me some heartburn. filter_design.py 1538 win = get_window(window,N,fftbins=1) 1539 alpha = N//2 1540 m = numpy.arange(0,N) 1541 h = win*special.sinc(cutoff*(m-alpha)) 1542 return h / sum(h) Line 1542 of filter_design.py, "return h / sum(h)", normalizes the function where it doesn't seem necessary, at least in the kaiser window case. Without the normalization, the kaiser window already returns a value of 1 at the zero frequency point. This normalization scales all of the data, making the window difficult to use in the frequency domain. Can someone point me to the rationale for this line? Looking at the code, this seems to be a pretty recent change (within the last year/year and a half). Thanks, Eric Buehler eric <dot> buehler <at> smiths-aerospace <dot> com ****************************************** The information contained in, or attached to, this e-mail, may contain confidential information and is intended solely for the use of the individual or entity to whom they are addressed and may be subject to legal privilege. If you have received this e-mail in error you should notify the sender immediately by reply e-mail, delete the message from your system and notify your system manager. Please do not copy it for any purpose, or disclose its contents to any other person. The views or opinions presented in this e-mail are solely those of the author and do not necessarily represent those of the company. The recipient should check this e-mail and any attachments for the presence of viruses. The company accepts no liability for any damage caused, directly or indirectly, by any virus transmitted in this email. ******************************************

3 2

NNLS
by Uwe Schmitt 25 Aug '08

25 Aug '08

Hi, I'm new to this mailing list. I just wraped NNLS (non negative least sqaures) from NETLIB using ctypes and want to ask if anybody is interested in my code ? As far as I know scipy does not contain this algorithm. Greetings, Uwe -- Dr. rer. nat. Uwe Schmitt F&E Mathematik mineway GmbH Science Park 2 D-66123 Saarbrücken Telefon: +49 (0)681 8390 5334 Telefax: +49 (0)681 830 4376 uschmitt(a)mineway.de www.mineway.de Geschäftsführung: Dr.-Ing. Mathias Bauer Amtsgericht Saarbrücken HRB 12339

4 14

Removing globals from minpack wrapper, v3
by Yosef Meller 04 Aug '08

04 Aug '08

Before I enter exams period, one more patch in the series: Changes: * Now optimize.leastsq() and optimize.fsolve() should be completely thread-safe. * Rearranged the macros to reduce code duplication. Still to do: * chkder, but I don't know if it's significant there. * smjac_* in __minpack.h - what is that function used for? -- http://yosefm.imagekind.com/Eclectic

3 5

physical quantities: udunits?
by Darren Dale 01 Aug '08

01 Aug '08

I was wondering if anyone is familiar with the udunits package, would it be suitable as the foundation for physical quantities support in python? http://www.unidata.ucar.edu/software/udunits/ According to freshmeat, udunits is released under a BSD-style license.

5 14

[ANNOUNCE] ETS 3.0.0b1 released!
by Dave Peterson 01 Aug '08

01 Aug '08

Hello, I am very pleased to announce that ETS 3.0.0b1 has just been tagged and released! All ETS sub-projects, including Traits, Chaco, Mayavi, and Envisage, have been tagged and released in one form or another (alpha, beta, final) as part of this effort. All ETS projects have now been registered with PyPi and source tarballs have been uploaded! As a result, for the first time ever, you can now easy_install all of ETS via a simple command like: easy_install ETS Or, alternatively, you can get only specific projects with commands like: easy_install Mayavi Or easy_install Chaco We anticipate the final release of ETS 3.0, and all of its sub-projects at their current version number, in a matter of weeks. Please join us in helping make this final release the most robust release it can be by participating in this beta shakeout period. Discussions about the status of ETS and its sub-projects happens on the enthought-dev mailing list, who's home page is at: https://mail.enthought.com/mailman/listinfo/enthought-dev Thank you to all who have worked so hard on making this release possible! -- Dave Enthought Tool Suite --------------------------- The Enthought Tool Suite (ETS) is a collection of components developed by Enthought and open source participants, which we use every day to construct custom scientific applications. It includes a wide variety of components, including: * an extensible application framework * application building blocks * 2-D and 3-D graphics libraries * scientific and math libraries * developer tools The cornerstone on which these tools rest is the Traits package, which provides explicit type declarations in Python; its features include initialization, validation, delegation, notification, and visualization of typed attributes. More information is available for all the packages within ETS from the Enthought Tool Suite development home page at http://code.enthought.com/projects/tool-suite.php. Testimonials ---------------- "I set out to rebuild an application in one week that had been developed over the last seven years (in C by generations of post-docs). Pyface and Traits were my cornerstones and I knew nothing about Pyface or Wx. It has been a hectic week. But here ... sits in front of me a nice application that does most of what it should. I think this has been a huge success. ... Thanks to the tools Enthought built, and thanks to the friendly support from people on the [enthought-dev] list, I have been able to build what I think is the best application so far. I have built similar applications (controlling cameras for imaging Bose-Einstein condensate) in C+MFC, Matlab, and C+labWindows, each time it has taken me at least four times longer to get to a result I regard as inferior. So I just wanted to say a big "thank you". Thank you to Enthought for providing this great software open-source. Thank you for everybody on the list for your replies." — Gaël Varoquaux, Laboratoire Charles Fabry, Institut d’Optique, Palaiseau, France "I'm currently writing a realtime data acquisition/display application … I'm using Enthought Tool Suite and Traits, and Chaco for display. IMHO, I think that in five years ETS/Traits will be the most comonly used framework for scientific applications." — Gary Pajer, Department of Chemistry, Biochemistry and Physics, Rider University, Lawrenceville NJ

1 0

latest svn broken?
by Robert Cimrman 31 Jul '08

31 Jul '08

>>> import scipy Traceback (most recent call last): File "<stdin>", line 1, in ? File "/usr/lib/python2.4/site-packages/scipy/__init__.py", line 86, in ? ImportError: cannot import name Tester Looks like this update to scipy/__init__.py caused it: r4545 | alan.mcintyre | 2008-07-16 02:26:23 +0200 (Wed, 16 Jul 2008) | 2 lines Remove the testing package, and switched all test modules to use numpy.test instead. r.

2 2

[ANNOUNCE] Traits 3.0 has been released
by Dave Peterson 31 Jul '08

31 Jul '08

Hello, I am very pleased to announce that Traits 3.0 has just been released! All Traits projects have been registered with PyPi (aka The Cheeseshop) and each project's listing on PyPi currently includes a source tarball. In the near future, we will also upload binary eggs for Windows and Mac OS X platforms. Installation of Traits 3.0 is now as simple as: easy_install Traits The Traits projects include: http://pypi.python.org/pypi?:action=display&name=Traits&version=3.0.0 http://pypi.python.org/pypi?:action=display&name=TraitsGUI&version=3.0.0 http://pypi.python.org/pypi?:action=display&name=TraitsBackendQt&version=3.… http://pypi.python.org/pypi?:action=display&name=TraitsBackendWX&version=3.… The Traits project is at the center of all Enthought Tool Suite development and has changed the mental model used at Enthought for programming in the already extremely efficient Python programming language. We encourage everyone to join us in enjoying the productivity gains from using such a powerful approach. The Traits project allows Python programmers to use a special kind of type definition called a trait, which gives object attributes some additional characteristics: * Initialization: A trait has a default value, which is automatically set as the initial value of an attribute before its first use in a program. * Validation: A trait attribute's type is explicitly declared. The type is evident in the code, and only values that meet a programmer-specified set of criteria (i.e., the trait definition) can be assigned to that attribute. * Delegation: The value of a trait attribute can be contained either in the defining object or in another object delegated to by the trait. * Notification: Setting the value of a trait attribute can notify other parts of the program that the value has changed. * Visualization: User interfaces that allow a user to interactively modify the value of a trait attribute can be automatically constructed using the trait's definition. (This feature requires that a supported GUI toolkit be installed. If this feature is not used, the Traits project does not otherwise require GUI support.) A class can freely mix trait-based attributes with normal Python attributes, or can opt to allow the use of only a fixed or open set of trait attributes within the class. Trait attributes defined by a classs are automatically inherited by any subclass derived from the class. -- Dave

1 0

fmin_cobyla hangs on solaris x86 (sometimes)
by John Hunter 30 Jul '08

30 Jul '08

I have been using fmin_cobya so solve a constrained optimization problem, and for some initial parameters, the code hangs (unresponsive to CTRL-C) indefinitely. I have a free-standing script that replicates the problem, and it always hangs at the same point in the execution. Unfortunately, the problem is not repeatable on other platforms I have tried (linux xeon and os x x86). I'm testing all platforms with svn HEAD. To make matters worse, I spent some time trying to isolate where in the FORTRAN the code was hanging, using poor-man's insert print statements into cobyla2.f and trstlp.f, divide and conquer. What I find is that by the time I have inserted enough print statements to get to narrow down where it is occuring, the hang disappears (heisenbug). The same thing happens if I insert no print statements but set iprint=3 (the hang disappears). But I had isolated the freeze to trstlp.f before the hang mysteriously disappeared when I tried to narrow it down further. So I am not too hopeful that anyone will be able to help me, but I wanted to put this out there in case anyone has any ideas what might be going wrong. Fernando suggested offlist that is might be a gcc compiler problem, since we are using an old version here 3.4.1 and gcc isn't the most robust compiler on solaris x86. In any case, here is output of the code, which hangs at function call 18: johnh@flag:tmp> python port_opt_test.py call: 1 call: 2 call: 3 call: 4 call: 5 call: 6 call: 7 call: 8 call: 9 call: 10 call: 11 call: 12 call: 13 call: 14 call: 15 call: 16 call: 17 call: 18 and here is the example code: import numpy as np import scipy.optimize as optimize class OptimalPortfolio: """ Compute the optimal weighting for a risky portfolio given returns matrix R and investor risk tolerance factor q; see http://en.wikipedia.org/wiki/Modern_portfolio_theory """ def __init__(self, R, Sigma): """ R is a Nasset return vector Sigma is Nasset x Nasset covariance matrix """ self.numassets = len(R) self.R = np.matrix(R) self.R.shape = self.numassets, 1 self.Sigma = np.matrix(Sigma) self._cnt = 0 assert self.Sigma.shape==(self.numassets, self.numassets) def __call__(self, weights): """ evaluate the cost function with portfolio weights """ w = np.matrix(weights) w.shape = self.numassets, 1 self._cnt += 1 print 'call: %d'%self._cnt ret = float(0.5 * w.transpose() * self.Sigma * w - self.q * self.R.transpose() * w) return ret def optimal(self, w0, bounds=None, q=0.5, **kwargs): """ return the optimal portfolio weighting starting with initial guess w0 bounds, if not None, is a sequence of (min, max) bounds for each weight. q is the investor risk preference if the kwargs are passed on to scipy.optimize.fmin_cobyla, eg to control the maximum number of function calls or the verbosity of the print output. See help for details """ w0 = np.asarray(w0) self.q = float(q) def weights_leq1(w): 'the sum of the weights must be <= 1' ws = -(w.sum()-1) return ws def weights_geq1(w): 'the sum of the weights must be >=1' ws = w.sum()-1 return ws constraints = [weights_leq1, weights_geq1] if bounds is None: bounds = [] for w in w0: bounds.append((0,1)) if len(bounds)!=len(w0): raise ValueError('length of bounds must be equal to length of weights') class MinBound: def __init__(self, i, bound): self.i = i self.bound = bound def __call__(self, weights): v = weights[self.i]-self.bound return v class MaxBound: def __init__(self, i, bound): self.i = i self.bound = bound def __call__(self, weights): v = self.bound-weights[self.i] return v for i, (bmin, bmax) in enumerate(bounds): if bmin<0 or bmin>1: raise ValueError('bounds must be in 0..1') if bmax<0 or bmax>1: raise ValueError('bounds must be in 0..1') constraints.append(MinBound(i, bmin)) constraints.append(MaxBound(i, bmax)) xopt = optimize.fmin_cobyla(self, w0, constraints, **kwargs) return xopt q = 0.3 R = np.fromstring('\xf3\x8ci\xd6\xee\xdb[?#\x85phk\x0b[?\x89"\xa9\r6\xc7K?5.-\xa0\xff\x1a6?\x96\xd7X<\x87tU?\xc18h\xda\r\x84P?\xeb\x0f$\xf7&ga?u&Ho:\xbc[?\x1c\xfd|\xad.\xb7M?\xd8\x1a\x992\xc2\r??k)]\x94\x8c"U?\xed\x00$\xaf\x06\x8eg?', np.float64) C = np.fromstring('\xe6EU"\x0b 3?\x8f\x90b\xbeJ[\xdc>.w\xb7\xc8"\x9f\xda\xbe\x99\x87\x02o\xff\x8e\xd6\xbet\x89t\xf4\x7f\xa2\xa3>\xe8\x99\xef\xceK\x0b\xcc>\xf2\xc7\xc3\xedx\xe8\xfa>#um#\x0b\xf4\x03?\xd9\x10\xeac.\x87\x08?\xbe$\x15\xb6\xf9N#?%\xc2\xc2\xba\x9d\x06\x1e?\xbfJ\x01\x04p}\x17?\x8f\x90b\xbeJ[\xdc>r3\x9f\xebW\xd92?U\xa1/\xa9rQ\xd6>\x81\x94 ,\xd9X\x01\xbfa1\x9cv5-\xd4\xbe!\x84\x87\xd3\xa5\x96\xe1\xbe\xf3\xc1m\x9b\\/\xda\xbe8F\xdaF\xac\x9a\xd0\xbe\xd7\x15\xf7\x16\x0f\xe6\xc1\xbe\x84~\xc4\xfa\xda\xde\xb8>\x8f\xccyO\xf7\x99\xe9>Yn\x1e6j^\xd3\xbe9w\xb7\xc8"\x9f\xda\xbeU\xa1/\xa9rQ\xd6>\xeb\x91%,\x99\x1e3?h\xce\x12\x14E\x8c\xe3\xbeN%\xca\xc51\x1e\xf8\xbeb)\xf7,\x14@\xe8\xbe\xa1\xc8`\x8d>\xc0\xb8\xbe(\xd67{c"\xe3\xbe\x94/\x9b\x13\x87\xf0\xe5>K\x9a\xb1j\xde\x13\xe6\xbeV\xf2\xa7n3\xe5\xf0\xbe\\\xb1\xb6`,g\xe7\xbe\x99\x87\x02o\xff\x8e\xd6\xbe\x81\x94 ,\xd9X\x01\xbfh\xce\x12\x14E\x8c\xe3\xbe\x83_\x1atI\x892?\x88\xb2\xdb\xcbl\x9c\xf6>\x17P\xfa\xa6\xf8C\xf2>g\x97\xd3\xb8\xa6\x93\xef\xbe\xbd\x90\xe6\xb7\xe0I\xbf\xbe\xceA\x8a\x14=\xa1\xd9\xbel`\xd8C\x05b\xd3\xbe[\xc2\xa3z:\x1c\xea\xbe1\x7f V(\xe1\xaf>t\x89t\xf4\x7f\xa2\xa3>a1\x9cv5-\xd4\xbe?%\xca\xc51\x1e\xf8\xbe\x88\xb2\xdb\xcbl\x9c\xf6>g\x87+O\xc7D1?\xe2o\x95\xd0\xa8?-?\xd8\xf0]_\x88l\xf5\xbe\x96FD\xbd\x99h\xdc\xbe\x9e\x19\x14\xf0\x89-\xd4\xbe\xcf\xdfc\xceP\x1d\xde\xbe\xb5\xd0\x84o"\x10\xee\xbe\x9cl\xa3\x0c\x15"\xc3\xbe\xe8\x99\xef\xceK\x0b\xcc>!\x84\x87\xd3\xa5\x96\xe1\xbeb)\xf7,\x14@\xe8\xbe\x17P\xfa\xa6\xf8C\xf2>\xe2o\x95\xd0\xa8?-?H\xe8\x91 \x08\x9b1?\x97(\xe6\xed<\x89\xe5\xbe\xfd\x00\xd4\xb6\x92p\xc6\xbe.Y\xad\xa8\xd6\x9f\xe0\xbe\x8f\xb32U%_\xd7>\x01F\xa8\x18\n\xd6\xd2\xbe\xed\xda4u\xeb\xd8\xce>\xef\xc7\xc3\xedx\xe8\xfa>\xf3\xc1m\x9b\\/\xda\xbe\xa1\xc8`\x8d>\xc0\xb8\xbeg\x97\xd3\xb8\xa6\x93\xef\xbe\xef\xf0]_\x88l\xf5\xbe\x97(\xe6\xed<\x89\xe5\xbe)\x04z(\xe5x7?:u~~\xb9S"?\x9f&\xfdV\x07\x84\xd4>\x84\xb0c{\x94c\xf0\xbeSFj\xf4m\xc3\x93>\xedW\x84\x98SA\x07?#um#\x0b\xf4\x03?7F\xdaF\xac\x9a\xd0\xbe(\xd67{c"\xe3\xbe\xbd\x90\xe6\xb7\xe0I\xbf\xbe\x96FD\xbd\x99h\xdc\xbe\xfc\x00\xd4\xb6\x92p\xc6\xbe:u~~\xb9S"?fcJ\x9a\xa8\xc15?\x95\xd6\xb1\xdc\xc0w\xf1>>"\x86\xce}\xd6\xea>\xed\x8dhG\xd4\x90\xf5>\xba\xf2\xb7\xb6~\xfb\x05?\xd9\x10\xeac.\x87\x08?\xd7\x15\xf7\x16\x0f\xe6\xc1\xbe\x97/\x9b\x13\x87\xf0\xe5>\xceA\x8a\x14=\xa1\xd9\xbe\x9e\x19\x14\xf0\x89-\xd4\xbe.Y\xad\xa8\xd6\x9f\xe0\xbe\xb5&\xfdV\x07\x84\xd4>\x95\xd6\xb1\xdc\xc0w\xf1>h\x0bv\xf9\xc4)3?\xe5\x14\xc7-N\x87\x18?\xd8F\xd6\xe8\xe3\x82\x02?N\x84\xa7\xb6\x93f\xe1\xbe\xbe$\x15\xb6\xf9N#?\x84~\xc4\xfa\xda\xde\xb8>K\x9a\xb1j\xde\x13\xe6\xbek`\xd8C\x05b\xd3\xbe\xcf\xdfc\xceP\x1d\xde\xbe\x8f\xb32U%_\xd7>\x84\xb0c{\x94c\xf0\xbe>"\x86\xce}\xd6\xea>\xe5\x14\xc7-N\x87\x18?X*\xa7\n\xa8\xb36?\xd4\xf3\x15xRr\x19?\xe6C\x17\x83\xcd\x86\xea>/\xc2\xc2\xba\x9d\x06\x1e?\x8f\xccyO\xf7\x99\xe9>V\xf2\xa7n3\xe5\xf0\xbe[\xc2\xa3z:\x1c\xea\xbe\x89\xd0\x84o"\x10\xee\xbe\x01F\xa8\x18\n\xd6\xd2\xbeSFj\xf4m\xc3\x93>\xed\x8dhG\xd4\x90\xf5>\xceF\xd6\xe8\xe3\x82\x02?\xd4\xf3\x15xRr\x19?\x89\x1d\x9f\xd6\xaa<9?\xcaj\xabG\'u\x0b?\xbfJ\x01\x04p}\x17?Zn\x1e6j^\xd3\xbe\\\xb1\xb6`,g\xe7\xbe1\x7f V(\xe1\xaf>\x9cl\xa3\x0c\x15"\xc3\xbe\xee\xda4u\xeb\xd8\xce>\xedW\x84\x98SA\x07?\xba\xf2\xb7\xb6~\xfb\x05?N\x84\xa7\xb6\x93f\xe1\xbe\xe6C\x17\x83\xcd\x86\xea>\xcaj\xabG\'u\x0b?\x8a\xbcAnx\xa32?', np.float64) C.shape = len(R), len(R) w0 = np.ones(len(R)) / len(R) opt = OptimalPortfolio(R, C) w = opt.optimal(w0, q=q, iprint=0) print w

1 0

scipy.signal.convolve2d very slow, why is it there exactly?
by Joel Schaerer 28 Jul '08

28 Jul '08

Hi all, I had previously reported that signal.convolve2d was significantly slower than matlab for the following code: ------------- #!/usr/bin/python import scipy.io as io import scipy.signal aa=io.loadmat('input.mat') seq=aa["seq"] #seq_f=scipy.empty(seq.shape) seq_f=[] kernel=scipy.randn(21,21) for k in xrange(seq.shape[2]): seq_f.append(scipy.signal.convolve2d(seq[:,:,k],kernel,'same')) io.savemat("output_python2.mat",{"seq_f":seq_f}) ----------------- Turns out, I can get even better results than matlab by using fftconvolve instead of convolve2d. The generic signal.convolve gives intermediary results, a lot better than convolve2d but still slower than matlab. Here are the results of my measurements for anyone interested: convolve2d : 9 minutes convolve : 3 minutes matlab's conv2 : 1m2s fftconvolve : 35s So the question seems to be, why is convolve2d even there? It's very slow and limited to 2D arrays. It's misleading to newcommers who might very well conclude from its use that scipy is too slow to be useful.

1 0