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