infinite-newton-cg

from scipy import optimize
import numpy as np

def gauss(x, amp, mu, std):
    y = amp*np.exp(-(x - mu)**2/(2*std**2))
    return y

def gauss_ndim(xdata, params):
    amps, mus, sigmas = params[::3], params[1::3], params[2::3]
    data = np.zeros(len(xdata))
    for amp, mu, sigma in zip(amps, mus, sigmas):
        data += gauss(xdata, amp, mu, sigma)
    return data

def gauss_jac(params, x, y):
    jac = np.zeros_like(params)
    common = -gauss_ndim(x, params) + y
    for i in xrange(params.shape[0]):
        if i%3 == 0:
            amp = params[i]
        elif i%3 == 1:
            mu = params[i]
        elif i%3 == 2:
            sigma = params[i]
            amp_term = np.exp(-(-mu + x)**2/(2*sigma**2))
            jac[i-2] += sum(-2*common*amp_term)
            jac[i-1] += sum(amp*(2*mu - 2*x)*common*amp_term/sigma**2)
            jac[i] += sum(-2*amp*(-mu + x)**2*common*amp_term/sigma**3)
    return jac

def gauss_hess(params, x, y):
    H = np.zeros((params.shape[0], params.shape[0]))
    common = -gauss_ndim(x, params) + y
    for i in xrange(params.shape[0]):
        if i%3 == 0:
            amp = params[i]
        elif i%3 == 1:
            mu = params[i]
        elif i%3 == 2:
            sigma = params[i]
            for j in xrange(0, params.shape[0], 3):
                hess_amp, hess_mu, hess_sigma = params[j:j+3]
                # the diagonal in general refers to the hessian with partial
                # derivatives against itself (so a 3x3 or 4x4 area)
                diagonal = i == j+2
                H[i-2, j] = sum(2 * np.exp((-(-mu+x)**2)/(sigma**2)) if diagonal else 2 * np.exp((-(-mu+x)**2)/(2*sigma**2)) * np.exp((-(-hess_mu+x)**2)/(2*hess_sigma**2)))
                H[i-2, j+1] = sum(-2*amp*(mu-x)/(sigma**2) * np.exp((-(-mu+x)**2)/(sigma**2)) + 2*(mu-x)/(sigma**2)*common*np.exp(-((-mu+x)**2)/(2*sigma**2)) if diagonal else -2*hess_amp/(hess_sigma**2)*(hess_mu-x)*np.exp(-((-mu+x)**2)/(2*sigma**2))*np.exp(-((-hess_mu+x)**2)/(2*hess_sigma**2)))
                H[i-2, j+2] = sum(2*amp*(-mu+x)**2/(sigma**3)*np.exp(-(-mu+x)**2/(sigma**2)) - 2*(-mu+x)**2/(sigma**3)*common*np.exp(-(-mu+x)**2/(2*sigma**2)) if diagonal else 2*hess_amp*(-hess_mu+x)**2/(hess_sigma**3) * np.exp(-((-mu+x)**2)/(2*sigma**2)) * np.exp(-((-hess_mu+x)**2)/(2*hess_sigma**2)))
                H[i-1, j] = sum(-2*amp*(mu-x)/(sigma**2)*np.exp(-((-mu+x)**2)/(sigma**2))+2*(mu-x)/(sigma**2)*common*np.exp(-((-mu+x)**2)/(2*sigma**2)) if diagonal else 2*amp*(mu-x)/(2*sigma**2)*np.exp(-((-mu+x)**2)/(2*sigma**2)) * 2*hess_amp*(hess_mu-x)/(2*hess_sigma**2)*np.exp(-((-hess_mu+x)**2)/(2*hess_sigma**2)))
                H[i-1, j+1] = sum((amp/(2*sigma**4)*(2*mu-2*x)**2)*np.exp(-((-mu+x)**2)/(sigma**2)) + 2*amp/(sigma**2)*np.exp(-(-mu+x)**2/(2*sigma**2))*common - amp*((2*mu-2*x)**2)/(2*sigma**4)*np.exp(-((-mu+x)**2)/(2*sigma**2))*common if diagonal else 2*amp*(mu-x)/(2*sigma**2)*np.exp(-((-mu+x)**2)/(2*sigma**2)) * 2*hess_amp*(hess_mu-x)/(2*hess_sigma**2)*np.exp(-((-hess_mu+x)**2)/(2*hess_sigma**2)))
                H[i-1, j+2] = sum(-2*(amp**2)/(sigma**5)*((-mu+x)**2)*(mu-x)*np.exp(-((-mu+x)**2)/(sigma**2)) - 2*amp*2*(mu-x)/(sigma**3)*np.exp(-((-mu+x)**2)/(2*sigma**2))*common + 2*amp*(mu-x)*(sigma**5)*((-mu+x)**2)*np.exp(-((-mu+x)**2)/(2*sigma**2))*common if diagonal else -hess_amp/(hess_sigma**3)*((-hess_mu+x)**2)*np.exp(-((-hess_mu+x)**2)/(2*hess_sigma**2)) * 2*amp*(mu-x)/(sigma**2)*np.exp(-((-mu+x)**2)/(2*sigma**2)))
                H[i, j] = sum(2*amp*((-mu+x)**2)/(sigma**3)*np.exp(-((-mu+x)**2)/(sigma**2)) - 2*((-mu+x)**2)/(sigma**3)*np.exp(-((-mu+x)**2)/(2*sigma**2))*common if diagonal else 2*amp*((-mu+x)**2)/(sigma**3)*np.exp(-((-mu+x)**2)/(2*sigma**2)) * np.exp(-((-hess_mu+x)**2)/(2*hess_sigma**2)))
                H[i, j+1] = sum(-(amp**2)*((-mu+x)**2)*2*(mu-x)*np.exp(-((-mu+x)**2)/(sigma**2)) - 2*amp*2*(mu-x)/(sigma**3)*np.exp(-((-mu+x)**2)/(2*sigma**2))*common + amp*((-mu+x)**2)*2*(mu-x)/(sigma**5)*np.exp(-((-mu+x)**2)/(2*sigma**2))*common if diagonal else -hess_amp*2*(hess_mu-x)/(hess_sigma**2)*np.exp(-((-hess_mu+x)**2)/(2*hess_sigma**2)) * amp*(-mu+x)**2/(sigma**3)*np.exp(-((-mu+x)**2)/(sigma**2)))
                H[i, j+2] = sum(2*(amp**2)*((-mu+x)**4)/(sigma**6)*np.exp(-((-mu+x)**2)/(sigma**2)) + 6*amp*((-mu+x)**2)/(sigma**4)*np.exp(-((-mu+x)**2)/(2*sigma**2))*common - 2*amp*((-mu+x)**4)/(sigma**6)*np.exp(-((-mu+x)**2)/(2*sigma**2))*common if diagonal else 2*amp*((-mu+x)**2)/(sigma**3)*np.exp(-((-mu+x)**2)/(sigma**2)) * 2*hess_amp*((-hess_mu+x)**2)/(hess_sigma**3)*np.exp(-((-hess_mu+x)**2)/(hess_sigma**2)))
    return H

def gauss_func(guess, xdata, ydata):
    data = gauss_ndim(xdata, guess)
    # absolute deviation as our distance metric. Empirically found to give better results than
    # residual sum of squares for this data.
    return sum((ydata-data)**2)

amp, mu, std, mu2 = 1., 0., 1., 3.
one_gauss_params = np.array([amp, mu, std], dtype=np.float)
two_gauss_params = np.array([amp, mu, std, amp, mu2, std], dtype=np.float)
x = np.array(np.linspace(-5, 5, 101), dtype=np.float)
one_gauss = gauss_ndim(x, one_gauss_params)
two_gauss = gauss_ndim(x, two_gauss_params)
# passes
optimize.minimize(gauss_func, [0.8, 0.1, 0.8], (x, one_gauss), options={'maxiter': 1000, 'disp': True}, method='Newton-CG', jac=gauss_jac, hess=gauss_hess)
# infinite loop
optimize.minimize(gauss_func, [0.8, 0.1, 0.8, 0.8, 2.5, 1.], (x, two_gauss), options={'maxiter': 1000, 'disp': True}, method='Newton-CG', jac=gauss_jac, hess=gauss_hess)