[SciPy-dev] fmin_cobyla hangs on solaris x86 (sometimes)

John Hunter jdh2358@gmail....
Wed Jul 30 09:16:46 CDT 2008


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


More information about the Scipy-dev mailing list