test function list¶

In this part, we will show gradient descent optimization on various 2-dimensioanl test function from this wiki page wiki/test-function-for-optimization.

We will plot the function landscape around the optimal point and also the optimization path.

We choose the initial point for optimization uniformly from the range $[-1, 1]$, so it could be biased to the local optimal point around the origin.

In [1]:

import numpy as np
import matplotlib.pyplot as plt
import torch

try:
    import numqi
except ImportError:
    %pip install numqi
    import numqi

np_rng = np.random.default_rng()
hf_uniform_para = lambda *x: torch.nn.Parameter(torch.tensor(np_rng.uniform(-1, 1, size=x), dtype=torch.float64))

import numpy as np import matplotlib.pyplot as plt import torch try: import numqi except ImportError: %pip install numqi import numqi np_rng = np.random.default_rng() hf_uniform_para = lambda *x: torch.nn.Parameter(torch.tensor(np_rng.uniform(-1, 1, size=x), dtype=torch.float64))

In [2]:

def hf_demo(model, num_repeat=3, xlim=None, ylim=None, x_optim=None, tag_logscale=False):
    callback = numqi.optimize.MinimizeCallback(extra_key='path', print_freq=0)
    theta_optim = numqi.optimize.minimize(model, tol=1e-10, num_repeat=num_repeat, callback=callback)
    path = np.stack(callback.state['path'])
    print(f'optimal theta: {theta_optim.x}')
    print(f'optimal loss: {theta_optim.fun}')

    hf_model = numqi.optimize.hf_model_wrapper(model)
    xdata = np.linspace(*xlim, 100)
    ydata = np.linspace(*ylim, 100)
    if len(x_optim)==2:
        zdata = np.array([[hf_model(np.array([x, y]), tag_grad=False) for x in xdata] for y in ydata])
    else:
        zdata = np.array([[hf_model(np.concatenate([np.array([x, y]),x_optim[2:]]), tag_grad=False) for x in xdata] for y in ydata])
    if tag_logscale:
        zdata = np.log(np.maximum(1e-4,zdata))
    fig,ax = plt.subplots()
    contour_set = ax.contourf(xdata, ydata, zdata, levels=15, cmap='winter')
    fig.colorbar(contour_set)
    ax.plot([x_optim[0]], [x_optim[1]], 'rx', label='optimal')
    ax.plot(path[:,0], path[:,1], '-', color='orange', label='path')
    ax.set_xlim(*xlim)
    ax.set_ylim(*ylim)
    ax.legend()
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    if tag_logscale:
        ax.set_title('log scale plot of loss function')
    fig.tight_layout()

def hf_demo(model, num_repeat=3, xlim=None, ylim=None, x_optim=None, tag_logscale=False): callback = numqi.optimize.MinimizeCallback(extra_key='path', print_freq=0) theta_optim = numqi.optimize.minimize(model, tol=1e-10, num_repeat=num_repeat, callback=callback) path = np.stack(callback.state['path']) print(f'optimal theta: {theta_optim.x}') print(f'optimal loss: {theta_optim.fun}') hf_model = numqi.optimize.hf_model_wrapper(model) xdata = np.linspace(*xlim, 100) ydata = np.linspace(*ylim, 100) if len(x_optim)==2: zdata = np.array([[hf_model(np.array([x, y]), tag_grad=False) for x in xdata] for y in ydata]) else: zdata = np.array([[hf_model(np.concatenate([np.array([x, y]),x_optim[2:]]), tag_grad=False) for x in xdata] for y in ydata]) if tag_logscale: zdata = np.log(np.maximum(1e-4,zdata)) fig,ax = plt.subplots() contour_set = ax.contourf(xdata, ydata, zdata, levels=15, cmap='winter') fig.colorbar(contour_set) ax.plot([x_optim[0]], [x_optim[1]], 'rx', label='optimal') ax.plot(path[:,0], path[:,1], '-', color='orange', label='path') ax.set_xlim(*xlim) ax.set_ylim(*ylim) ax.legend() ax.set_xlabel('x') ax.set_ylabel('y') if tag_logscale: ax.set_title('log scale plot of loss function') fig.tight_layout()

Rastrigin function¶

$$ f(x)=An+\sum_{i=1}^n [x_i^2-A\cos(2\pi x_i)], A=10$$

$$ f(0,\cdots,0)=0 $$

In [3]:

class Rastrigin(torch.nn.Module):
    def __init__(self, n):
        super().__init__()
        self.x = hf_uniform_para(n)
        self.A = 10
        # solution [0,0,...,0], 0

    def forward(self):
        x = self.x
        loss = self.A*x.shape[0] + (x*x - self.A*torch.cos(2*np.pi*x)).sum()
        return loss

n = 2
model = Rastrigin(n)
x_optim = np.zeros(n)
hf_demo(model, num_repeat=10, xlim=(-5.12, 5.12), ylim=(-5.12, 5.12), x_optim=x_optim)
# when n is large, it's almost impossible to find the global minimum

class Rastrigin(torch.nn.Module): def __init__(self, n): super().__init__() self.x = hf_uniform_para(n) self.A = 10 # solution [0,0,...,0], 0 def forward(self): x = self.x loss = self.A*x.shape[0] + (x*x - self.A*torch.cos(2*np.pi*x)).sum() return loss n = 2 model = Rastrigin(n) x_optim = np.zeros(n) hf_demo(model, num_repeat=10, xlim=(-5.12, 5.12), ylim=(-5.12, 5.12), x_optim=x_optim) # when n is large, it's almost impossible to find the global minimum

[round=0] min(f)=0.9949590570932898, current(f)=0.9949590570932898
[round=1] min(f)=0.0, current(f)=0.0
[round=2] min(f)=0.0, current(f)=0.9949590570932898
[round=3] min(f)=0.0, current(f)=1.9899181141865796
[round=4] min(f)=0.0, current(f)=16.91420284955502
[round=5] min(f)=0.0, current(f)=7.105427357601002e-15
[round=6] min(f)=0.0, current(f)=0.9949590570932898
[round=7] min(f)=0.0, current(f)=0.9949590570932898
[round=8] min(f)=0.0, current(f)=1.9899181141865867
[round=9] min(f)=0.0, current(f)=0.9949590570932898
optimal theta: [ 4.51037240e-12 -6.44942176e-10]
optimal loss: 0.0

Ackley function¶

$$ f(x,y)=-20\mathrm{exp}[-0.2\sqrt{0.5(x^2+y^2)}] - \mathrm{exp}[0.5(\cos(2\pi x)+\cos(2\pi y))] + e + 20 $$

$$ f(0,0)=0 $$

In [4]:

class Ackley(torch.nn.Module):
    def __init__(self):
        super().__init__()
        self.theta = hf_uniform_para(2)
        # solution [0,0] 0

    def forward(self):
        x,y = self.theta
        tmp0 = -20*torch.exp(-0.2*torch.sqrt(0.5*(x*x + y*y)))
        tmp1 = -torch.exp(0.5*(torch.cos(2*np.pi*x)+torch.cos(2*np.pi*y))) + np.e + 20
        ret = tmp0 + tmp1
        return ret

model = Ackley()
x_optim = np.zeros(2)
hf_demo(model, num_repeat=10, xlim=(-5, 5), ylim=(-5, 5), x_optim=x_optim)

class Ackley(torch.nn.Module): def __init__(self): super().__init__() self.theta = hf_uniform_para(2) # solution [0,0] 0 def forward(self): x,y = self.theta tmp0 = -20*torch.exp(-0.2*torch.sqrt(0.5*(x*x + y*y))) tmp1 = -torch.exp(0.5*(torch.cos(2*np.pi*x)+torch.cos(2*np.pi*y))) + np.e + 20 ret = tmp0 + tmp1 return ret model = Ackley() x_optim = np.zeros(2) hf_demo(model, num_repeat=10, xlim=(-5, 5), ylim=(-5, 5), x_optim=x_optim)

[round=0] min(f)=2.5799275570298725, current(f)=2.5799275570298725
[round=1] min(f)=2.5799275570298725, current(f)=3.574451877257676
[round=2] min(f)=2.5799275570298725, current(f)=2.5799275570298725
[round=3] min(f)=2.579927557029869, current(f)=2.579927557029869
[round=4] min(f)=2.579927557029869, current(f)=3.5744518772576797
[round=5] min(f)=2.579927557029869, current(f)=2.5799275570299756
[round=6] min(f)=9.048761739904876e-11, current(f)=9.048761739904876e-11
[round=7] min(f)=9.048761739904876e-11, current(f)=2.5799275570298903
[round=8] min(f)=9.048761739904876e-11, current(f)=3.5744518772576797
[round=9] min(f)=9.048761739904876e-11, current(f)=2.5799275570298725
optimal theta: [-1.42936591e-11  2.86218824e-11]
optimal loss: 9.048761739904876e-11

Rosenbrock function¶

$$ f(x)=\sum_{i=1}^{n-1}[100(x_{i+1}-x_i^2)^2 + (1-x_i)^2] $$

$$ f(1,\cdots,1)=0 $$

In [5]:

class Rosenbrock(torch.nn.Module):
    def __init__(self, n):
        super().__init__()
        self.theta = hf_uniform_para(n)
        # solution [1,1,...,1] 0

    def forward(self):
        tmp0 = self.theta[1:] - self.theta[:-1]**2
        tmp1 = 1-self.theta[:-1]
        ret = 100*torch.dot(tmp0, tmp0) + torch.dot(tmp1,tmp1)
        return ret

n = 2
model = Rosenbrock(n)
x_optim = np.ones(n)
hf_demo(model, num_repeat=10, xlim=(-2, 2), ylim=(-1, 3), x_optim=x_optim, tag_logscale=True)

class Rosenbrock(torch.nn.Module): def __init__(self, n): super().__init__() self.theta = hf_uniform_para(n) # solution [1,1,...,1] 0 def forward(self): tmp0 = self.theta[1:] - self.theta[:-1]**2 tmp1 = 1-self.theta[:-1] ret = 100*torch.dot(tmp0, tmp0) + torch.dot(tmp1,tmp1) return ret n = 2 model = Rosenbrock(n) x_optim = np.ones(n) hf_demo(model, num_repeat=10, xlim=(-2, 2), ylim=(-1, 3), x_optim=x_optim, tag_logscale=True)

[round=0] min(f)=1.2639860158827894e-19, current(f)=1.2639860158827894e-19
[round=1] min(f)=1.2639860158827894e-19, current(f)=3.4143011335172966e-19
[round=2] min(f)=1.2639860158827894e-19, current(f)=1.623849586387913e-19
[round=3] min(f)=1.2639860158827894e-19, current(f)=1.309454909852618e-13
[round=4] min(f)=1.2639860158827894e-19, current(f)=2.7376182986765596e-19
[round=5] min(f)=1.2639860158827894e-19, current(f)=3.01206762634152e-13
[round=6] min(f)=1.2639860158827894e-19, current(f)=2.0083017631151303e-17
[round=7] min(f)=1.2639860158827894e-19, current(f)=2.1534881151585027e-14
[round=8] min(f)=8.463749318267302e-20, current(f)=8.463749318267302e-20
[round=9] min(f)=8.463749318267302e-20, current(f)=2.9760011914445767e-15
optimal theta: [1. 1.]
optimal loss: 8.463749318267302e-20

Beale function¶

$$ f(x,y)=(1.5-x+xy)^2 + (2.25-x+xy^2)^2 + (2.625-x+xy^3)^2 $$

$$ f(3,0.5)=0 $$

In [6]:

class Beale(torch.nn.Module):
    def __init__(self) -> None:
        super().__init__()
        self.theta = hf_uniform_para(2)
        # [3,0.5] 0

    def forward(self):
        x,y = self.theta
        ret = (1.5 - x + x*y)**2 + (2.25 - x + x*y**2)**2 + (2.625 - x + x*y**3)**2
        return ret

model = Beale()
x_optim = np.array([3, 0.5])
hf_demo(model, num_repeat=10, xlim=(-4.5, 4.5), ylim=(-4.5, 4.5), x_optim=x_optim, tag_logscale=True)

class Beale(torch.nn.Module): def __init__(self) -> None: super().__init__() self.theta = hf_uniform_para(2) # [3,0.5] 0 def forward(self): x,y = self.theta ret = (1.5 - x + x*y)**2 + (2.25 - x + x*y**2)**2 + (2.625 - x + x*y**3)**2 return ret model = Beale() x_optim = np.array([3, 0.5]) hf_demo(model, num_repeat=10, xlim=(-4.5, 4.5), ylim=(-4.5, 4.5), x_optim=x_optim, tag_logscale=True)

[round=0] min(f)=9.550644975334525e-15, current(f)=9.550644975334525e-15
[round=1] min(f)=4.40736957365878e-18, current(f)=4.40736957365878e-18
[round=2] min(f)=1.740336010493108e-20, current(f)=1.740336010493108e-20
[round=3] min(f)=1.740336010493108e-20, current(f)=2.792038894125406e-19
[round=4] min(f)=1.740336010493108e-20, current(f)=1.0806746258653203e-15
[round=5] min(f)=1.740336010493108e-20, current(f)=2.0943463350959024e-16
[round=6] min(f)=1.740336010493108e-20, current(f)=3.449445256081479e-14
[round=7] min(f)=2.2398965846651986e-26, current(f)=2.2398965846651986e-26
[round=8] min(f)=2.2398965846651986e-26, current(f)=0.4532781718454318
[round=9] min(f)=2.2398965846651986e-26, current(f)=1.0153810449367475e-14
optimal theta: [3.  0.5]
optimal loss: 2.2398965846651986e-26

Goldstein-Price function¶

$$ f(x,y)=[1+(x+y+1)^2(19-14x+3x^2-14y+6xy+3y^2)][30+(2x-3y)^2(18-32x+12x^2+48y-36xy+27y^2)]-3 $$

$$ f(0,-1)=0 $$

comment: the function value is shifted by a constant $3$ for good plotting.

In [7]:

class GoldsteinPrice(torch.nn.Module):
    def __init__(self) -> None:
        super().__init__()
        self.theta = hf_uniform_para(2)
        # [0,-1] 0

    def forward(self):
        x,y = self.theta
        tmp0 = 1 + (x+y+1)**2 * (19-14*x+3*x*x-14*y+6*x*y+3*y*y)
        tmp1 = 30 + (2*x-3*y)**2 * (18-32*x+12*x*x+48*y-36*x*y+27*y*y)
        ret = tmp0*tmp1 - 3 #shift 3 for good plotting
        return ret

model = GoldsteinPrice()
x_optim = np.array([0, -1])
hf_demo(model, num_repeat=10, xlim=(-2, 2), ylim=(-3, 1), x_optim=x_optim, tag_logscale=True)

class GoldsteinPrice(torch.nn.Module): def __init__(self) -> None: super().__init__() self.theta = hf_uniform_para(2) # [0,-1] 0 def forward(self): x,y = self.theta tmp0 = 1 + (x+y+1)**2 * (19-14*x+3*x*x-14*y+6*x*y+3*y*y) tmp1 = 30 + (2*x-3*y)**2 * (18-32*x+12*x*x+48*y-36*x*y+27*y*y) ret = tmp0*tmp1 - 3 #shift 3 for good plotting return ret model = GoldsteinPrice() x_optim = np.array([0, -1]) hf_demo(model, num_repeat=10, xlim=(-2, 2), ylim=(-3, 1), x_optim=x_optim, tag_logscale=True)

[round=0] min(f)=27.0, current(f)=27.0
[round=1] min(f)=0.0, current(f)=0.0
[round=2] min(f)=0.0, current(f)=7.105427357601002e-15
[round=3] min(f)=0.0, current(f)=8.704148513061227e-14
[round=4] min(f)=-1.0658141036401503e-14, current(f)=-1.0658141036401503e-14
[round=5] min(f)=-1.0658141036401503e-14, current(f)=27.0
[round=6] min(f)=-1.0658141036401503e-14, current(f)=3.197442310920451e-14
[round=7] min(f)=-6.394884621840902e-14, current(f)=-6.394884621840902e-14
[round=8] min(f)=-6.394884621840902e-14, current(f)=-2.4868995751603507e-14
[round=9] min(f)=-6.394884621840902e-14, current(f)=27.00000000000142
optimal theta: [-8.08005565e-10 -9.99999999e-01]
optimal loss: -6.394884621840902e-14

Booth function¶

$$ f(x,y)=(x+2y-7)^2 + (2x+y-5)^2 $$

$$ f(1,3)=0 $$

In [8]:

class Booth(torch.nn.Module):
    def __init__(self) -> None:
        super().__init__()
        self.theta = hf_uniform_para(2)
        # [1,3] 0

    def forward(self):
        x,y = self.theta
        tmp0 = x+2*y-7
        tmp1 = 2*x+y-5
        ret = tmp0*tmp0 + tmp1*tmp1
        return ret

model = Booth()
x_optim = np.array([1, 3])
hf_demo(model, num_repeat=10, xlim=(-10, 10), ylim=(-10, 10), x_optim=x_optim, tag_logscale=True)

class Booth(torch.nn.Module): def __init__(self) -> None: super().__init__() self.theta = hf_uniform_para(2) # [1,3] 0 def forward(self): x,y = self.theta tmp0 = x+2*y-7 tmp1 = 2*x+y-5 ret = tmp0*tmp0 + tmp1*tmp1 return ret model = Booth() x_optim = np.array([1, 3]) hf_demo(model, num_repeat=10, xlim=(-10, 10), ylim=(-10, 10), x_optim=x_optim, tag_logscale=True)

[round=0] min(f)=7.328366864799852e-17, current(f)=7.328366864799852e-17
[round=1] min(f)=8.757725601558686e-19, current(f)=8.757725601558686e-19
[round=2] min(f)=5.3945992048351445e-22, current(f)=5.3945992048351445e-22
[round=3] min(f)=5.3945992048351445e-22, current(f)=6.687337534864696e-22
[round=4] min(f)=5.3945992048351445e-22, current(f)=1.5742554806826965e-21
[round=5] min(f)=5.3945992048351445e-22, current(f)=4.841600505028755e-18
[round=6] min(f)=2.5935379980951205e-26, current(f)=2.5935379980951205e-26
[round=7] min(f)=1.5777218104420236e-30, current(f)=1.5777218104420236e-30
[round=8] min(f)=1.5777218104420236e-30, current(f)=3.362628554741069e-17
[round=9] min(f)=1.5777218104420236e-30, current(f)=7.057151235828982e-24
optimal theta: [1. 3.]
optimal loss: 1.5777218104420236e-30

Bukin Function N.6¶

$$ f(x,y)=100\sqrt{|y-0.01x^2|}+0.01|x+10| $$

$$ f(-10,1)=0 $$

In [9]:

class BukinFunctionN6(torch.nn.Module):
    def __init__(self) -> None:
        super().__init__()
        self.theta = hf_uniform_para(2)
        # [-10,1] 0

    def forward(self):
        x,y = self.theta
        # possible numerical unstable here, also hard for Adam
        ret = 100*torch.sqrt(torch.abs(y-0.01*x*x)) + 0.01*torch.abs(x+10)
        return ret

model = BukinFunctionN6()
x_optim = np.array([-10, 1])
hf_demo(model, num_repeat=10, xlim=(-15, -5), ylim=(-4, 6), x_optim=x_optim, tag_logscale=True)

class BukinFunctionN6(torch.nn.Module): def __init__(self) -> None: super().__init__() self.theta = hf_uniform_para(2) # [-10,1] 0 def forward(self): x,y = self.theta # possible numerical unstable here, also hard for Adam ret = 100*torch.sqrt(torch.abs(y-0.01*x*x)) + 0.01*torch.abs(x+10) return ret model = BukinFunctionN6() x_optim = np.array([-10, 1]) hf_demo(model, num_repeat=10, xlim=(-15, -5), ylim=(-4, 6), x_optim=x_optim, tag_logscale=True)

[round=0] min(f)=7.995237311598105, current(f)=7.995237311598105
[round=1] min(f)=7.995237311598105, current(f)=20.744564182075024
[round=2] min(f)=7.995237311598105, current(f)=47.44056601472836
[round=3] min(f)=7.995237311598105, current(f)=63.387719797292455
[round=4] min(f)=4.140171204777714, current(f)=4.140171204777714
[round=5] min(f)=4.140171204777714, current(f)=nan
[round=6] min(f)=0.25744936598052937, current(f)=0.25744936598052937

[round=7] min(f)=0.25744936598052937, current(f)=480.7121964002795
[round=8] min(f)=0.25744936598052937, current(f)=1081.1035234675178
[round=9] min(f)=0.25744936598052937, current(f)=90.13261498640846
optimal theta: [9.32673574e-02 5.28605102e-05]
optimal loss: 0.25744936598052937

Matyas function¶

$$ f(x,y)=0.26(x^2+y^2)-0.48xy $$

$$ f(0,0)=0 $$

In [10]:

class Matyas(torch.nn.Module):
    def __init__(self) -> None:
        super().__init__()
        self.theta = hf_uniform_para(2)
        # [-0,0] 0

    def forward(self):
        x,y = self.theta
        ret = 0.26*(x*x + y*y) - 0.48*x*y
        return ret

model = Matyas()
x_optim = np.array([0, 0])
hf_demo(model, num_repeat=10, xlim=(-10, 10), ylim=(-10, 10), x_optim=x_optim, tag_logscale=True)

class Matyas(torch.nn.Module): def __init__(self) -> None: super().__init__() self.theta = hf_uniform_para(2) # [-0,0] 0 def forward(self): x,y = self.theta ret = 0.26*(x*x + y*y) - 0.48*x*y return ret model = Matyas() x_optim = np.array([0, 0]) hf_demo(model, num_repeat=10, xlim=(-10, 10), ylim=(-10, 10), x_optim=x_optim, tag_logscale=True)

[round=0] min(f)=1.397269878372716e-30, current(f)=1.397269878372716e-30
[round=1] min(f)=1.397269878372716e-30, current(f)=2.3899057235360275e-19
[round=2] min(f)=1.3816327809443071e-31, current(f)=1.3816327809443071e-31
[round=3] min(f)=8.581943832189534e-33, current(f)=8.581943832189534e-33
[round=4] min(f)=8.581943832189534e-33, current(f)=2.9902758688533968e-31
[round=5] min(f)=8.581943832189534e-33, current(f)=1.080369661603464e-31
[round=6] min(f)=8.581943832189534e-33, current(f)=3.7028444079564495e-15
[round=7] min(f)=8.581943832189534e-33, current(f)=1.067427412377183e-31
[round=8] min(f)=8.581943832189534e-33, current(f)=9.444760447275019e-31
[round=9] min(f)=8.581943832189534e-33, current(f)=1.31425459404985e-32
optimal theta: [4.44089210e-16 4.71844785e-16]
optimal loss: 8.581943832189534e-33

Levi function N.13¶

$$ f(x,y)=\sin^2(3\pi x)+ (x-1)^2(1+\sin^2(3\pi y)) + (y-1)^2(1+\sin^2(2\pi y)) $$

$$ f(1,1)=0 $$

In [11]:

class LeviFunction(torch.nn.Module):
    def __init__(self) -> None:
        super().__init__()
        self.theta = hf_uniform_para(2)
        # [0,0] 0

    def forward(self):
        x,y = self.theta
        ret = torch.sin(3*np.pi*x)**2 + (x-1)**2*(1+torch.sin(3*np.pi*y)**2) + (y-1)**2*(1+torch.sin(2*np.pi*y)**2)
        return ret

model = LeviFunction()
x_optim = np.array([1, 1])
hf_demo(model, num_repeat=10, xlim=(-10, 10), ylim=(-10, 10), x_optim=x_optim)

class LeviFunction(torch.nn.Module): def __init__(self) -> None: super().__init__() self.theta = hf_uniform_para(2) # [0,0] 0 def forward(self): x,y = self.theta ret = torch.sin(3*np.pi*x)**2 + (x-1)**2*(1+torch.sin(3*np.pi*y)**2) + (y-1)**2*(1+torch.sin(2*np.pi*y)**2) return ret model = LeviFunction() x_optim = np.array([1, 1]) hf_demo(model, num_repeat=10, xlim=(-10, 10), ylim=(-10, 10), x_optim=x_optim)

[round=0] min(f)=0.10987365835883793, current(f)=0.10987365835883793
[round=1] min(f)=7.281733255089467e-16, current(f)=7.281733255089467e-16
[round=2] min(f)=7.281733255089467e-16, current(f)=0.10987365835883804
[round=3] min(f)=7.281733255089467e-16, current(f)=0.6316425103184409
[round=4] min(f)=7.281733255089467e-16, current(f)=0.10987365835883871
[round=5] min(f)=7.281733255089467e-16, current(f)=1.182316538532437
[round=6] min(f)=3.3861083352848305e-18, current(f)=3.3861083352848305e-18
[round=7] min(f)=3.3861083352848305e-18, current(f)=4.952205341584722
[round=8] min(f)=3.3861083352848305e-18, current(f)=2.75270922028409
[round=9] min(f)=3.3861083352848305e-18, current(f)=2.7465302928244633
optimal theta: [1. 1.]
optimal loss: 3.3861083352848305e-18

Himmelblau's function¶

$$ f(x,y)=(x^2+y-11)^2 + (x+y^2-7)^2 $$

$$ f(3,2)=f(-2.805118,3.131312)=f(-3.779310,-3.283186)=f(3.584428,-1.848126)=0 $$

In [12]:

class Himmelblau(torch.nn.Module):
    def __init__(self) -> None:
        super().__init__()
        self.theta = hf_uniform_para(2)
        # [3,2],[-2.805118,3.131312],[-3.779310,-3.283186],[3.584428,-1.848126] 0

    def forward(self):
        x,y = self.theta
        tmp0 = x*x + y - 11
        tmp1 = x + y*y - 7
        ret = tmp0*tmp0 + tmp1*tmp1
        return ret

model = Himmelblau()
x_optim = np.array([[3, 2],[-2.805118,3.131312],[-3.779310,-3.283186],[3.584428,-1.848126]])
hf_demo(model, num_repeat=10, xlim=(-5, 5), ylim=(-5, 5), x_optim=x_optim[0], tag_logscale=True)

class Himmelblau(torch.nn.Module): def __init__(self) -> None: super().__init__() self.theta = hf_uniform_para(2) # [3,2],[-2.805118,3.131312],[-3.779310,-3.283186],[3.584428,-1.848126] 0 def forward(self): x,y = self.theta tmp0 = x*x + y - 11 tmp1 = x + y*y - 7 ret = tmp0*tmp0 + tmp1*tmp1 return ret model = Himmelblau() x_optim = np.array([[3, 2],[-2.805118,3.131312],[-3.779310,-3.283186],[3.584428,-1.848126]]) hf_demo(model, num_repeat=10, xlim=(-5, 5), ylim=(-5, 5), x_optim=x_optim[0], tag_logscale=True)

[round=0] min(f)=1.6298243412751323e-18, current(f)=1.6298243412751323e-18
[round=1] min(f)=3.733537809572998e-19, current(f)=3.733537809572998e-19
[round=2] min(f)=1.4181712529275272e-22, current(f)=1.4181712529275272e-22
[round=3] min(f)=1.4181712529275272e-22, current(f)=5.762989524783656e-18
[round=4] min(f)=1.4181712529275272e-22, current(f)=3.603815021392727e-18
[round=5] min(f)=1.4181712529275272e-22, current(f)=6.136043407070982e-20
[round=6] min(f)=1.4181712529275272e-22, current(f)=1.0212519845220034e-15
[round=7] min(f)=1.4181712529275272e-22, current(f)=3.741436603324366e-18
[round=8] min(f)=1.4181712529275272e-22, current(f)=2.306131702168999e-16
[round=9] min(f)=1.4181712529275272e-22, current(f)=1.0010855576225206e-21
optimal theta: [3. 2.]
optimal loss: 1.4181712529275272e-22

Three-hump camel function¶

$$ f(x,y)=2x^2-1.05x^4+\frac{x^6}{6}+xy+y^2 $$

$$ f(0,0)=0 $$

In [13]:

class ThreeHumpCamel(torch.nn.Module):
    def __init__(self) -> None:
        super().__init__()
        self.theta = hf_uniform_para(2)
        # [0,0] 0

    def forward(self):
        x,y = self.theta
        ret = 2*x*x - 1.05*x*x*x*x + x*x*x*x*x*x/6 + x*y + y*y
        return ret

model = ThreeHumpCamel()
x_optim = np.array([0, 0])
hf_demo(model, num_repeat=10, xlim=(-5, 5), ylim=(-5, 5), x_optim=x_optim, tag_logscale=True)

class ThreeHumpCamel(torch.nn.Module): def __init__(self) -> None: super().__init__() self.theta = hf_uniform_para(2) # [0,0] 0 def forward(self): x,y = self.theta ret = 2*x*x - 1.05*x*x*x*x + x*x*x*x*x*x/6 + x*y + y*y return ret model = ThreeHumpCamel() x_optim = np.array([0, 0]) hf_demo(model, num_repeat=10, xlim=(-5, 5), ylim=(-5, 5), x_optim=x_optim, tag_logscale=True)

[round=0] min(f)=3.702127835467326e-16, current(f)=3.702127835467326e-16
[round=1] min(f)=3.702127835467326e-16, current(f)=2.5250534692635775e-15
[round=2] min(f)=3.4623859820660514e-19, current(f)=3.4623859820660514e-19
[round=3] min(f)=3.4623859820660514e-19, current(f)=7.067218075686402e-17
[round=4] min(f)=3.4623859820660514e-19, current(f)=1.5776743117754285e-16
[round=5] min(f)=3.4623859820660514e-19, current(f)=1.6450376582766694e-18
[round=6] min(f)=3.4623859820660514e-19, current(f)=1.1538084294666314e-15
[round=7] min(f)=1.0973534534997989e-20, current(f)=1.0973534534997989e-20
[round=8] min(f)=1.0973534534997989e-20, current(f)=3.278758156209364e-14
[round=9] min(f)=1.0973534534997989e-20, current(f)=2.6319834144319464e-17
optimal theta: [1.13766593e-11 9.79795770e-11]
optimal loss: 1.0973534534997989e-20

Easom function¶

$$ f(x,y)=-\cos(x)\cos(y)\exp[-(x-\pi)^2-(y-\pi)^2] $$

$$ f(\pi,\pi)=-1 $$

In [14]:

class Easom(torch.nn.Module):
    def __init__(self) -> None:
        super().__init__()
        self.theta = hf_uniform_para(2)
        # [pi,pi] -1

    def forward(self):
        x,y = self.theta
        ret = -torch.cos(x)*torch.cos(y)*torch.exp(-(x-np.pi)**2-(y-np.pi)**2)
        return ret

model = Easom()
x_optim = np.array([np.pi, np.pi])
hf_demo(model, num_repeat=10, xlim=(-100, 100), ylim=(-100, 100), x_optim=x_optim)

class Easom(torch.nn.Module): def __init__(self) -> None: super().__init__() self.theta = hf_uniform_para(2) # [pi,pi] -1 def forward(self): x,y = self.theta ret = -torch.cos(x)*torch.cos(y)*torch.exp(-(x-np.pi)**2-(y-np.pi)**2) return ret model = Easom() x_optim = np.array([np.pi, np.pi]) hf_demo(model, num_repeat=10, xlim=(-100, 100), ylim=(-100, 100), x_optim=x_optim)

[round=0] min(f)=-1.0, current(f)=-1.0
[round=1] min(f)=-1.0, current(f)=-6.8661887399191835e-15
[round=2] min(f)=-1.0, current(f)=-1.0
[round=3] min(f)=-1.0, current(f)=-1.0
[round=4] min(f)=-1.0, current(f)=-0.6074873577867838
[round=5] min(f)=-1.0, current(f)=-0.8198623458484934
[round=6] min(f)=-1.0, current(f)=-8.110223891561197e-05
[round=7] min(f)=-1.0, current(f)=-1.0
[round=8] min(f)=-1.0, current(f)=-0.9993947224509269
[round=9] min(f)=-1.0, current(f)=-1.0
optimal theta: [3.14159265 3.14159265]
optimal loss: -1.0

Cross-in-tray function¶

$$ f(x,y)=-0.0001\left[\left|\sin(x)\sin(y)\exp\left(\left|100-\frac{\sqrt{x^2+y^2}}{\pi}\right|\right)\right|+1\right]^{0.1} $$

$$ f(1.34941,-1.34941)=f(1.34941,1.34941)=f(-1.34941,1.34941)=f(-1.34941,-1.34941)=-2.06261 $$

In [15]:

class CrossInTray(torch.nn.Module):
    def __init__(self) -> None:
        super().__init__()
        self.theta = hf_uniform_para(2)
        # [1.34941,-1.34941],[1.34941,1.34941],[-1.34941,1.34941],[-1.34941,-1.34941] -2.06261

    def forward(self):
        x,y = self.theta
        tmp0 = torch.abs(100 - torch.sqrt(x*x + y*y)/np.pi)
        tmp1 = torch.abs(torch.sin(x)*torch.sin(y)*torch.exp(tmp0))
        ret = -0.0001*tmp1**0.1
        return ret

model = CrossInTray()
x_optim = np.array([[1.34941,-1.34941],[1.34941,1.34941],[-1.34941,1.34941],[-1.34941,-1.34941]])
hf_demo(model, num_repeat=10, xlim=(-10, 10), ylim=(-10, 10), x_optim=x_optim[0])

class CrossInTray(torch.nn.Module): def __init__(self) -> None: super().__init__() self.theta = hf_uniform_para(2) # [1.34941,-1.34941],[1.34941,1.34941],[-1.34941,1.34941],[-1.34941,-1.34941] -2.06261 def forward(self): x,y = self.theta tmp0 = torch.abs(100 - torch.sqrt(x*x + y*y)/np.pi) tmp1 = torch.abs(torch.sin(x)*torch.sin(y)*torch.exp(tmp0)) ret = -0.0001*tmp1**0.1 return ret model = CrossInTray() x_optim = np.array([[1.34941,-1.34941],[1.34941,1.34941],[-1.34941,1.34941],[-1.34941,-1.34941]]) hf_demo(model, num_repeat=10, xlim=(-10, 10), ylim=(-10, 10), x_optim=x_optim[0])

[round=0] min(f)=-2.062611870822727, current(f)=-2.062611870822727
[round=1] min(f)=-2.062611870822736, current(f)=-2.062611870822736
[round=2] min(f)=-2.062611870822736, current(f)=-2.0626118708226504
[round=3] min(f)=-2.062611870822736, current(f)=-2.0626118708226393
[round=4] min(f)=-2.062611870822736, current(f)=-2.0626118708227343
[round=5] min(f)=-2.062611870822736, current(f)=-2.062611870822733
[round=6] min(f)=-2.062611870822736, current(f)=-2.0626118708223564
[round=7] min(f)=-2.062611870822736, current(f)=-2.062611870822689
[round=8] min(f)=-2.0626118708227374, current(f)=-2.0626118708227374
[round=9] min(f)=-2.0626118708227374, current(f)=-2.062611870822492
optimal theta: [ 1.34940669 -1.34940652]
optimal loss: -2.0626118708227374

Eggholder function¶

$$ f(x,y)=-(y+47)\sin\sqrt{|x/2+(y+47)|}-x\sin\sqrt{|x-(y+47)|} $$

$$ f(512,404.2319)=-959.6407 $$

In [16]:

class Eggholder(torch.nn.Module):
    def __init__(self) -> None:
        super().__init__()
        self.theta = hf_uniform_para(2)
        # [512,404.2319] -959.6407

    def forward(self):
        x,y = self.theta
        tmp0 = -(y+47)*torch.sin(torch.sqrt(torch.abs(x/2 + (y+47))))
        tmp1 = -x*torch.sin(torch.sqrt(torch.abs(x-(y+47))))
        ret = tmp0 + tmp1
        return ret

model = Eggholder()
x_optim = np.array([512, 404.2319])
hf_demo(model, num_repeat=10, xlim=(-1000, 1000), ylim=(-1000, 1000), x_optim=x_optim)

class Eggholder(torch.nn.Module): def __init__(self) -> None: super().__init__() self.theta = hf_uniform_para(2) # [512,404.2319] -959.6407 def forward(self): x,y = self.theta tmp0 = -(y+47)*torch.sin(torch.sqrt(torch.abs(x/2 + (y+47)))) tmp1 = -x*torch.sin(torch.sqrt(torch.abs(x-(y+47)))) ret = tmp0 + tmp1 return ret model = Eggholder() x_optim = np.array([512, 404.2319]) hf_demo(model, num_repeat=10, xlim=(-1000, 1000), ylim=(-1000, 1000), x_optim=x_optim)

[round=0] min(f)=-66.84371732946161, current(f)=-66.84371732946161
[round=1] min(f)=-66.84371732946161, current(f)=-66.84371732943485
[round=2] min(f)=-66.84371732946161, current(f)=-66.84371732946006
[round=3] min(f)=-66.84371732946161, current(f)=-66.84371732944952
[round=4] min(f)=-66.84371732946398, current(f)=-66.84371732946398
[round=5] min(f)=-66.84371732946398, current(f)=-66.84371732946384
[round=6] min(f)=-66.843717329464, current(f)=-66.843717329464
[round=7] min(f)=-66.843717329464, current(f)=-66.84371732946312
[round=8] min(f)=-66.843717329464, current(f)=-66.84371732946397
[round=9] min(f)=-66.843717329464, current(f)=-66.84371732946396
optimal theta: [ 8.45693446 15.65091813]
optimal loss: -66.843717329464

Holder table function¶

$$ f(x,y)=-\left|\sin(x)\cos(y)\exp\left(\left|1-\frac{\sqrt{x^2+y^2}}{\pi}\right|\right)\right| $$

$$ f(8.05502,9.66459)=f(-8.05502, 9.66459)=f(8.05502, -9.66459)=f(-8.05502, -9.66459)=-19.2085 $$

In [17]:

class HolderTable(torch.nn.Module):
    def __init__(self) -> None:
        super().__init__()
        self.theta = hf_uniform_para(2)
        # [8.05502,9.66459] -19.2085

    def forward(self):
        x,y = self.theta
        tmp0 = torch.sin(x)*torch.cos(y)
        tmp1 = torch.exp(torch.abs(1 - torch.sqrt(x*x + y*y)/np.pi))
        ret = -torch.abs(tmp0*tmp1)
        return ret

model = HolderTable()
x_optim = np.array([[8.05502, 9.66459], [-8.05502, 9.66459], [8.05502, -9.66459], [-8.05502, -9.66459]])
hf_demo(model, num_repeat=10, xlim=(-10, 10), ylim=(-10, 10), x_optim=x_optim[1])

class HolderTable(torch.nn.Module): def __init__(self) -> None: super().__init__() self.theta = hf_uniform_para(2) # [8.05502,9.66459] -19.2085 def forward(self): x,y = self.theta tmp0 = torch.sin(x)*torch.cos(y) tmp1 = torch.exp(torch.abs(1 - torch.sqrt(x*x + y*y)/np.pi)) ret = -torch.abs(tmp0*tmp1) return ret model = HolderTable() x_optim = np.array([[8.05502, 9.66459], [-8.05502, 9.66459], [8.05502, -9.66459], [-8.05502, -9.66459]]) hf_demo(model, num_repeat=10, xlim=(-10, 10), ylim=(-10, 10), x_optim=x_optim[1])

[round=0] min(f)=-1.7329722438960702, current(f)=-1.7329722438960702
[round=1] min(f)=-1.7329722438960702, current(f)=-1.7329722438960684
[round=2] min(f)=-1.7329722438960702, current(f)=-1.73297224389607
[round=3] min(f)=-1.7329722438960702, current(f)=-1.73297224389607
[round=4] min(f)=-1.7329722438960702, current(f)=-1.7329722438960702
[round=5] min(f)=-1.7329722438960702, current(f)=-1.732972243896069
[round=6] min(f)=-1.7329722438960702, current(f)=-1.7329722438960684
[round=7] min(f)=-1.7329722438960702, current(f)=-1.7329722438960578
[round=8] min(f)=-1.7329722438960702, current(f)=-1.7329722438960702
[round=9] min(f)=-1.7329722438960704, current(f)=-1.7329722438960704
optimal theta: [-1.26262726e+00 -4.30338145e-11]
optimal loss: -1.7329722438960704

McCormick function¶

$$ f(x,y)=\sin(x+y)+(x-y)^2-1.5x+2.5y+1 $$

$$ f(-0.54719,-1.54719)=-1.9133 $$

In [18]:

class McCormick(torch.nn.Module):
    def __init__(self) -> None:
        super().__init__()
        self.theta = hf_uniform_para(2)
        # [-0.54719,-1.54719] -1.9133

    def forward(self):
        x,y = self.theta
        tmp0 = torch.sin(x+y)
        tmp1 = (x-y)**2
        tmp2 = -1.5*x + 2.5*y + 1
        ret = tmp0 + tmp1 + tmp2
        return ret

# out of domain
model = McCormick()
x_optim = np.array([-0.54719,-1.54719])
hf_demo(model, num_repeat=10, xlim=(-2.5, 4), ylim=(-3, 4), x_optim=x_optim)

class McCormick(torch.nn.Module): def __init__(self) -> None: super().__init__() self.theta = hf_uniform_para(2) # [-0.54719,-1.54719] -1.9133 def forward(self): x,y = self.theta tmp0 = torch.sin(x+y) tmp1 = (x-y)**2 tmp2 = -1.5*x + 2.5*y + 1 ret = tmp0 + tmp1 + tmp2 return ret # out of domain model = McCormick() x_optim = np.array([-0.54719,-1.54719]) hf_demo(model, num_repeat=10, xlim=(-2.5, 4), ylim=(-3, 4), x_optim=x_optim)

[round=0] min(f)=-1.9132229549810353, current(f)=-1.9132229549810353
[round=1] min(f)=-1.9132229549810362, current(f)=-1.9132229549810362
[round=2] min(f)=-1.9132229549810367, current(f)=-1.9132229549810367
[round=3] min(f)=-1.9132229549810367, current(f)=-1.9132229549810362
[round=4] min(f)=-1.9132229549810367, current(f)=-1.9132229549810362
[round=5] min(f)=-1.9132229549810367, current(f)=-1.9132229549810362
[round=6] min(f)=-1.9132229549810367, current(f)=-1.913222954981028
[round=7] min(f)=-1.9132229549810367, current(f)=-1.913222954981027
[round=8] min(f)=-11.338000915750419, current(f)=-11.338000915750419
[round=9] min(f)=-11.338000915750419, current(f)=-1.9132229549810236
optimal theta: [ -9.97197551 -10.97197551]
optimal loss: -11.338000915750419

Schaffer function N.2¶

$$ f(x,y)=0.5+\frac{\sin^2(x^2-y^2)-0.5}{[1+0.001(x^2+y^2)]^2} $$

$$ f(0,0)=0 $$

In [19]:

class SchafferFunctionN2(torch.nn.Module):
    def __init__(self) -> None:
        super().__init__()
        self.theta = hf_uniform_para(2)
        # [0,0] 0

    def forward(self):
        x,y = self.theta
        tmp0 = torch.sin(x*x-y*y)
        tmp1 = 1 + 0.001*(x*x + y*y)
        ret = 0.5 + (tmp0*tmp0 - 0.5)/(tmp1*tmp1)
        return ret

model = SchafferFunctionN2()
x_optim = np.array([0, 0])
hf_demo(model, num_repeat=10, xlim=(-50, 50), ylim=(-50, 50), x_optim=x_optim)

class SchafferFunctionN2(torch.nn.Module): def __init__(self) -> None: super().__init__() self.theta = hf_uniform_para(2) # [0,0] 0 def forward(self): x,y = self.theta tmp0 = torch.sin(x*x-y*y) tmp1 = 1 + 0.001*(x*x + y*y) ret = 0.5 + (tmp0*tmp0 - 0.5)/(tmp1*tmp1) return ret model = SchafferFunctionN2() x_optim = np.array([0, 0]) hf_demo(model, num_repeat=10, xlim=(-50, 50), ylim=(-50, 50), x_optim=x_optim)

[round=0] min(f)=0.0, current(f)=0.0
[round=1] min(f)=0.0, current(f)=1.34692257347524e-12
[round=2] min(f)=0.0, current(f)=4.440892098500626e-16
[round=3] min(f)=0.0, current(f)=4.5075054799781356e-14
[round=4] min(f)=0.0, current(f)=2.3314683517128287e-13
[round=5] min(f)=0.0, current(f)=4.993783164763954e-13
[round=6] min(f)=0.0, current(f)=1.6084911180769268e-12
[round=7] min(f)=0.0, current(f)=1.7674750552032492e-13
[round=8] min(f)=0.0, current(f)=1.021405182655144e-14
[round=9] min(f)=0.0, current(f)=1.9710257870286796e-09
optimal theta: [-3.86757880e-08 -2.40706185e-08]
optimal loss: 0.0

Schaffer function N.4¶

$$ f(x,y)=0.5+\frac{\cos^2(\sin|x^2-y^2|)-0.5}{[1+0.001(x^2+y^2)]^2} $$

$$ f(0,1.25313)=f(0,-1.25313)=f(1.25313,0)=f(-1.25313,0)=0.292579 $$

In [20]:

class SchafferFunctionN4(torch.nn.Module):
    def __init__(self) -> None:
        super().__init__()
        self.theta = hf_uniform_para(2)
        # [0,1.25313] 0.292579

    def forward(self):
        x,y = self.theta
        tmp0 = torch.cos(torch.sin(torch.abs(x*x-y*y)))
        tmp1 = 1 + 0.001*(x*x + y*y)
        ret = 0.5 + (tmp0*tmp0 - 0.5)/(tmp1*tmp1)
        return ret

model = SchafferFunctionN4()
x_optim = np.array([[0, 1.25313], [0, -1.25313], [1.25313,0], [-1.25313,0]])
hf_demo(model, num_repeat=10, xlim=(-50, 50), ylim=(-50, 50), x_optim=x_optim[0])

class SchafferFunctionN4(torch.nn.Module): def __init__(self) -> None: super().__init__() self.theta = hf_uniform_para(2) # [0,1.25313] 0.292579 def forward(self): x,y = self.theta tmp0 = torch.cos(torch.sin(torch.abs(x*x-y*y))) tmp1 = 1 + 0.001*(x*x + y*y) ret = 0.5 + (tmp0*tmp0 - 0.5)/(tmp1*tmp1) return ret model = SchafferFunctionN4() x_optim = np.array([[0, 1.25313], [0, -1.25313], [1.25313,0], [-1.25313,0]]) hf_demo(model, num_repeat=10, xlim=(-50, 50), ylim=(-50, 50), x_optim=x_optim[0])

[round=0] min(f)=0.34584416439334076, current(f)=0.34584416439334076
[round=1] min(f)=0.2925786320360779, current(f)=0.2925786320360779
[round=2] min(f)=0.29257863203598067, current(f)=0.29257863203598067
[round=3] min(f)=0.29257863203598067, current(f)=0.2925786320359848
[round=4] min(f)=0.29257863203598067, current(f)=0.2925786320359811
[round=5] min(f)=0.29257863203598067, current(f)=0.29257863203600587
[round=6] min(f)=0.29257863203598067, current(f)=0.29387375881132394
[round=7] min(f)=0.29257863203598067, current(f)=0.29257863203598233
[round=8] min(f)=0.2925786320359806, current(f)=0.2925786320359806
[round=9] min(f)=0.2925786320359806, current(f)=0.2925786320362329
optimal theta: [-2.03812623e-07 -1.25313183e+00]
optimal loss: 0.2925786320359806

Styblinski-Tang function¶

$$ f(x)=\frac{1}{2}\sum_{i=1}^n[x_i^4-16x_i^2+5x_i] $$

$$ f(-2.903534,-2.903534,\cdots,-2.903534)=-39.16599n $$

In [21]:

class StyblinskiTangFunction(torch.nn.Module):
    def __init__(self, n) -> None:
        super().__init__()
        self.theta = hf_uniform_para(n)
        # [-2.903534]*n -39.16599*n

    def forward(self):
        x = self.theta
        ret = 0.5*(x*x*x*x - 16*x*x + 5*x).sum()
        return ret

n = 2
model = StyblinskiTangFunction(n)
x_optim = np.array([-2.903534]*n)
hf_demo(model, num_repeat=10, xlim=(-5, 5), ylim=(-5, 5), x_optim=x_optim)

class StyblinskiTangFunction(torch.nn.Module): def __init__(self, n) -> None: super().__init__() self.theta = hf_uniform_para(n) # [-2.903534]*n -39.16599*n def forward(self): x = self.theta ret = 0.5*(x*x*x*x - 16*x*x + 5*x).sum() return ret n = 2 model = StyblinskiTangFunction(n) x_optim = np.array([-2.903534]*n) hf_demo(model, num_repeat=10, xlim=(-5, 5), ylim=(-5, 5), x_optim=x_optim)

[round=0] min(f)=-78.33233140754214, current(f)=-78.33233140754214
[round=1] min(f)=-78.33233140754214, current(f)=-64.19561235905516
[round=2] min(f)=-78.33233140754214, current(f)=-50.058893310567896
[round=3] min(f)=-78.33233140754277, current(f)=-78.33233140754277
[round=4] min(f)=-78.33233140754277, current(f)=-78.33233140754173
[round=5] min(f)=-78.33233140754277, current(f)=-78.33233140754263
[round=6] min(f)=-78.33233140754277, current(f)=-64.19561235905536
[round=7] min(f)=-78.33233140754277, current(f)=-50.058893310567896
[round=8] min(f)=-78.33233140754277, current(f)=-64.19561235905533
[round=9] min(f)=-78.33233140754277, current(f)=-50.05889331056785
optimal theta: [-2.90353407 -2.90353399]
optimal loss: -78.33233140754277