<jupyter_start><jupyter_text>HW13.1 DFPBy ZincCat用DFP优化$f(\mathbf{x})=\frac{x_{1}^{4}}{4}+\frac{x_{2}^{2}}{2}-x_{1} x_{2}+x_{1}-x_{2}$<jupyter_code>%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt
from matplotlib.colors import LogNorm

# 设置随机种子
np.random.seed(19890817)<jupyter_output><empty_output><jupyter_text>DFP<jupyter_code># 定义一些辅助函数
def f(x):
    # 计算函数值
    return x[0]**4/4 + x[1]*x[1]/2 - x[0]*x[1] + x[0] - x[1]

def gradient_f(x):
    # 计算函数梯度
    return np.array([x[0]*x[0]*x[0] - x[1] + 1, x[1] - x[0] - 1])

def linesearch(f, x, g, d, alpha=0.4, beta=0.8):
    # backtrack linesearch
    t = 1.0
    value = f(x)
    while f(x + t*d) > value + alpha*t*np.dot(g, d):
        t *= beta
    return t
alpha = 0.4
beta = 0.8
maxIter = 100

# 设置变量初值
# x = np.array([0, 0], dtype=np.float64) #x0
x = np.array([1.5, 1], dtype=np.float64) #x0
H = np.eye(2) #H_0
g = gradient_f(x)
d = -g
eta = 1e-5
eps = 1e-18
timestep = 0

Plot = True
xlog1 = []
xlog2 = []

while True:
    print(timestep, "th iteration, x=", x, ", f(x)=", f(x))
    print("H=", H)
    print("g=", g)
    if Plot:
        xlog1.append(x[0])
        xlog2.append(x[1])
    if np.linalg.norm(g) < eta: # 算法终止条件
        break
    alpha = linesearch(f, x, g, d, alpha, beta)
    dx = alpha*d
    x += dx
    dg = gradient_f(x) - g
    g += dg
    t = H@dg
    if np.abs(np.dot(dg, t)) < eps: #避免除0
        break
    H = H + np.outer(dx, dx)/np.dot(dx, dg) - np.outer(t, t)/np.dot(dg, t) # DFP
    d = -H@g
    timestep += 1
    if timestep >= maxIter:
        break
print("result: x:", x, "\nf(x):",f(x))<jupyter_output><empty_output><jupyter_text>后面的代码仅用于绘图<jupyter_code>xmin, xmax, xstep = -4.5, 4.5, .2
ymin, ymax, ystep = -4.5, 4.5, .2
x = np.meshgrid(np.arange(xmin, xmax + xstep, xstep), np.arange(ymin, ymax + ystep, ystep))
z = f(x)+1
zmin, zmax = np.min(z), np.max(z)
minima1 = np.array([1, 2]).reshape(-1, 1)
minima2 = np.array([-1, 0]).reshape(-1, 1)
fig = plt.figure(figsize=(8, 5))
ax = plt.axes(projection='3d', elev=50, azim=-50)
ax.plot_surface(*x, z, edgecolor='none', norm = LogNorm(), alpha=.8, cmap=plt.get_cmap('rainbow'))
ax.plot(*minima1, f(minima1)+1, 'r*', markersize=10)
ax.plot(*minima2, f(minima2)+1, 'r*', markersize=10)

ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_zlabel('$f$')

ax.set_xlim((xmin, xmax))
ax.set_ylim((ymin, ymax))
plt.savefig('131func.pdf')
plt.show()
a = plt.contourf(*x, z, 10, cmap=plt.get_cmap('rainbow'))
b = plt.contour(*x, z, levels = np.logspace(-2,3,100))
plt.xlabel('$x_1$'); plt.ylabel("$x_2$")
plt.title("Contour plot of loss function for $x_1 = 1.5, x_2=1$")
plt.plot(xlog1, xlog2, 'r--')
plt.savefig('131cont1.51.pdf')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.22017673048600883, 'fraction_numerical': 0.053387334315169364, 'mean_word_length': 3.68448275862069, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '458050', 'license': 'MIT', 'n_tokens_mistral': 1399, 'n_tokens_neox': 1355, 'n_words': 255, 'path': '13/131.ipynb', 'repo_name': 'zinccat/Convex-Analysis-homework', 'size': '2716'}