2025.7.28 Himmelblau関数を共役勾配法で解く

関数 search_min_alpha 、刻み幅の設定に注意が必要。

こちらも初期条件のチューニングが必要

パラメータチューニングさえ行えば、それなりに安定して解を求めることができるようになった。

code:p.py

import numpy as np

import torch as pt

from torch.linalg import inv, eigvals

from torch.autograd.functional import jacobian, hessian

import matplotlib.pyplot as plt

def himmelblau(x):

x, y = x

return (x**2 + y - 11)**2 + (x + y**2 - 7)**2

def plot_himmelblau(levels=30):

x = pt.linspace(-5, 5, 100)

y = pt.linspace(-5, 5, 100)

xx, yy = pt.meshgrid(x, y, indexing='xy')

zz = himmelblau((xx, yy)) # タプルで与えることに注意

plt.contour(xx, yy, zz, levels=levels)

# plt.pcolormesh(xx, yy, zz)

plt.colorbar()

def search_min_alpha(func, x, dx):

h = 0.001 # 十分に小さくないと、谷を飛び越す可能性がある

alpha = h

# eps = 0.0001

loop_count = 0

loop_max = 10000

while True:

f = func(x)

x_next = x + dx*h

f_next = func(x_next)

if loop_count >= loop_max:

# print('# searchを打ち切りました', loop_count, alpha, end='')

break

if pt.abs(f_next - f) < eps:

# print('# search is converged at', loop_count, x, pt.linalg.norm(dx).item(), alpha, end='')

break

if f_next < f: # 減少しているので、xを先に進める

print('+', end='')

x = x_next

alpha += h

else: # 増加したので、xは進めずにhを半分にする

print('/', end='')

h = h / 2

loop_count += 1

return alpha

def armijo_backtracking(func, x, p, alpha=0.1, rho=0.1, c=0.0001):

# x ... 現在の点

# p ... 探索の向き

# c ... 1だとxにおける接線（大きい減衰を要求）、0だと

loop_count = 0

grad = - jacobian(func, x)

eps = 0.001

fx = func(x)

while True:

x_next = x + alpha*p

fx_next = func(x_next)

if fx_next <= fx + c * alpha * grad @ p:

myprint('# linear search is converged at', loop_count, x, alpha)

break

alpha *= rho

myprint('/', end='')

if alpha <= eps:

myprint('# searchを打ち切りました', loop_count, alpha)

break

loop_count += 1

return alpha

######################

def solve_gd():

print('# Gradient Descent method')

x = pt.tensor(x0, dtype=pt.float)

x_hist = []

loop_count = 0

while True:

x_hist.append(x)

dx = - jacobian(himmelblau, x)

# alpha = search_min_alpha(himmelblau, x, dx)

alpha = armijo_backtracking(himmelblau, x, dx)

x = x + alpha * dx

fn_dx = pt.linalg.norm(dx) # Frobenius Norm of ds

if fn_dx < eps:

print('converged at', loop_count)

hess = hessian(himmelblau, x)

break

if pt.isnan(fn_dx):

print('発散した')

break

if loop_count >= loop_max:

print('計算を打ち切った', loop_count, dx)

print('収束していない？解と勾配のノルム', x, dx, pt.linalg.norm(dx))

break

loop_count += 1

x_hist = pt.stack(x_hist, dim=0)

return x_hist

def solve_nt():

print('# Newton method')

x = pt.tensor(x0, dtype=pt.float)

x_hist = []

loop_count = 0

while True:

x_hist.append(x)

grad = jacobian(himmelblau, x)

hess = hessian(himmelblau, x)

dx = - inv(hess)@grad

x = x + dx

fn_dx = pt.linalg.norm(dx) # Frobenius Norm of ds

if fn_dx < eps:

print('converged at', loop_count)

hess = hessian(himmelblau, x)

break

if pt.isnan(fn_dx):

print('発散した')

break

if loop_count >= loop_max:

print('計算を打ち切った', loop_count, dx)

print('収束していない？解と勾配のノルム', x, pt.linalg.norm(dx))

break

loop_count += 1

x_hist = pt.stack(x_hist, dim=0)

return x_hist

def solve_cggd(meth):

# meth

# 1 ... Fletcher–Reeves / フレッチャー・リーブス

# 2 ... Polak–Ribiere / ポラック・リビエール

# 3 ... Hestenes-Stiefel /

# 4 ... Dai–Yuan

# dx ... x_{n} (now)

# dx_p ... x_{n-1} (previous)

print('# Nonlinear Conjugate Gradient Method ... ', end=' ')

if meth == 1:

print(' Fletcher–Reeves')

def beta(dx, dx_p, s_p):

return (dx @ dx) / (dx_p @ dx_p)

# return (dx.T @ dx) / (dx_p.T @ dx_p)

elif meth == 2:

print(' Polak–Ribiere')

def beta(dx, dx_p, s_p):

return dx @ (dx - dx_p) / (dx_p @ dx_p)

# return dx.T @ (dx - dx_p) / (dx_p.T @ dx_p)

elif meth == 3:

print(' Hestenes-Stiefel')

def beta(dx, dx_p, s_p):

return - dx @ (dx - dx_p) / s_p @ (dx - dx_p)

# return - dx.T @ (dx - dx_p) / s_p.T @ (dx - dx_p)

elif meth == 4:

print(' Dai–Yuan')

def beta(dx, dx_p, s_p):

return - (dx @ dx) / (s_p @ (dx - dx_p))

# return - (dx.T @ dx) / (s_p.T @ (dx - dx_p))

x = pt.tensor(x0, dtype=pt.float)

x_hist = []

loop_count = 0

# 0-step

x_hist.append(x)

dx = - jacobian(himmelblau, x)

# alpha = search_min_alpha(himmelblau, x, dx)

alpha = armijo_backtracking(himmelblau, x, dx)

x = x + alpha * dx

s = dx

while True:

x_hist.append(x)

dx_p = dx

dx = - jacobian(himmelblau, x)

s_p = s

b = max(0, beta(dx, dx_p, s_p))

if loop_count % n != 0:

s = dx

else:

s = dx + b * s_p # CGの更新

# alpha = search_min_alpha(himmelblau, x, dx)

alpha = armijo_backtracking(himmelblau, x, s)

x = x + alpha * s

# print(loop_count, pt.linalg.norm(dx).item(), dx @ dx, dx, s, b.item())

fn_dx = pt.linalg.norm(dx) # frobenius norm of dx

if fn_dx < eps:

print('converged at', loop_count)

hess = hessian(himmelblau, x)

break

if pt.isnan(fn_dx):

print('発散した')

break

if loop_count >= loop_max:

print('計算を打ち切った', loop_count, dx)

print('収束していない？解と勾配のノルム', x, dx, pt.linalg.norm(dx))

break

loop_count += 1

x_hist = pt.stack(x_hist, dim=0)

return x_hist

alpha = 0.001

loop_max = 1000

eps = 0.001

x_hist = solve_gd()

x_hist = solve_cggd(1)

plot_himmelblau(levels=100)

plt.legend()

plt.show()