[SciPy-user] optimization using fmin_bfgs with gradient information

[Sebastian Walter]

Hi,

thanks for the function :).
 It is now part of pyadolc's unit test.
http://github.com/b45ch1/pyadolc/blob/3cfff80c062d43e812379f4606eda0ebaaf4e82c/tests/complicated_tests.py
 , Line 246

I added two functions: one is providing the gradient by finite
differences and the other by using automatic differentiation.
The finite differences gradient has very poor accuracy, only the first
three digits are correct.


Sebastian



2009/7/19 Ernest Adrogué <eadrogue at gmx.net>:
> 19/07/09 @ 17:41 (+0200), thus spake Sebastian Walter:
>> Hello,
>> I would like to add your objective function as part of the unit test
>> for on of the automatic differentiation tools I've been developing:
>> PYADOLC (download at http://github.com/b45ch1/pyadolc).
>>
>> Unfortunately, the description of the objective function is
>> incomplete: there are some references to self.n and so on in the code.
>> Would it be possible for you to post the full description? That would
>> be great :)
>
> Yes, of course.
> You can use this version that works independently. The original
> one feeds on data from a dataset, where as this one uses randomly-generated
> data. Example usage:
>
> In [6]: m = Model(random_obs())
>
> In [7]: x0 = m.guess()
>
> In [8]: x0
> Out[8]:
> [1.3414634146341464,
>  1.4634146341463414,
>  1.6341463414634145,
>  1.5609756097560976,
>  1.4878048780487805,
>  1.8048780487804879,
>  ........
>  1.3902439024390243,
>  1.6097560975609757,
>  1.5,
>  0]
>
> In [9]: m.fobj(x0)
> Out[9]: 1281.027571347247
>
>
>
> ------------------------------
> import random
> import math
>
> class Observation(object):
>    def __init__(self, data):
>        self.ht = data[0]
>        self.at = data[1]
>        self.hg = data[2]
>        self.ag = data[3]
>
> def random_obs():
>    names = 'abcdefghijklmnopqrstuvwxyz'
>    random_data = [(i,j,random.randint(0,5),random.randint(0,5))
>                   for i in names for j in names if i != j]
>    return [Observation(i) for i in random_data]
>
> class Model(object):
>    def __init__(self, observations):
>        self.observations = tuple(observations)
>        self.names = [i.ht for i in observations]
>        self.names.extend([i.at for i in observations])
>        self.names = tuple(set(self.names))
>        self.n = len(self.names)
>    def guess(self):
>        a, b = [], []
>        total_hg = sum([i.hg for i in self.observations])
>        for j in self.names:
>            sh_j = sum([i.hg for i in self.observations if i.ht == j])
>            ca_j = sum([i.hg for i in self.observations if i.at == j])
>            a.append(sh_j/math.sqrt(total_hg))
>            b.append(ca_j/math.sqrt(total_hg))
>        return a + b + [1.5, 0]
>    def tau(self, mu1, mu2, rho, x, y):
>        if x == 0 and y == 0:
>            t = 1 - mu1 * mu2 * rho
>        elif x == 0 and y == 1:
>            t = 1 + mu1 * rho
>        elif x == 1 and y == 0:
>            t = 1 + mu2 * rho
>        elif x == 1 and y == 1:
>            t = 1 - rho
>        else:
>            t = 1
>        return t
>    def fobj(self, x):
>        n = self.n
>        y = [abs(i) for i in x[:-2]]
>        y.insert(0, abs(n - sum(y[:n-1])))
>        a = dict(zip(self.names, y[:n]))
>        b = dict(zip(self.names, y[n:]))
>        g = abs(x[-2])
>        r = x[-1]
>        pseudo_loglikelihood = 0
>        for m in self.observations:
>            x = m.hg
>            y = m.ag
>            mu1 = a[m.ht] * b[m.at] * g
>            mu2 = a[m.at] * b[m.ht]
>            tau = self.tau(mu1, mu2, r, m.hg, m.ag)
>            mu1 = mu1 > 0 and mu1 or 1e-10
>            mu2 = mu2 > 0 and mu2 or 1e-10
>            tau = tau > 0 and tau or 1e-10
>            pseudo_loglikelihood += math.log(tau)
>            pseudo_loglikelihood += m.hg * math.log(mu1) - mu1
>            pseudo_loglikelihood += m.ag * math.log(mu2) - mu2
>        return -pseudo_loglikelihood
>
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>