最优化第五次作业

最新推荐文章于 2024-07-18 15:53:23 发布

大写-凌祁

最新推荐文章于 2024-07-18 15:53:23 发布

阅读量785

点赞数 21

文章标签： python numpy 开发语言

本文链接：https://blog.csdn.net/qq_73553710/article/details/138979629

版权

$\min\:\sum_{i=1}^{11}\big(y_i-\tilde{f}(x,t_i)\big)^2,$
其中
$\tilde{f}(x,t)=\frac{x_1(t^2+x_2t)}{t^2+x_3t+x_4},$
数据 $t_i,y_i)$ 由表 5.5 给出。

表 5.5 第 4 题的数据

2	$t_i$	$y_i$	$i$	$t_i$	$y_i$
1	4.0000	0.1957	7	0.1250	0.0456
2	2.0000	0.1947	8	0.1000	0.0342
3	1.0000	0.1735	9	0.0833	0.0323
4	0.5000	0.1600	10	0.0714	0.0235
5	0.2500	0.0844	11	0.0625	0.0246
6	0.1670	0.0627

最速下降法求解

from scipy.optimize import approx_fprime
import numpy as np

# 定义数据点
t_data = np.array([4.0000, 2.0000, 1.0000, 0.5000, 0.2500, 0.1670, 0.1250, 0.1000, 0.0833, 0.0714, 0.0625])
y_data = np.array([0.1957, 0.1947, 0.1735, 0.1600, 0.0844, 0.0627, 0.0456, 0.0342, 0.0323, 0.0235, 0.0246])

# 定义模型函数
def model(x, t):
    x1, x2, x3, x4 = x
    return (x1 * (t**2 + x2*t)) / (t**2 + x3*t + x4)

# 定义要最小化的残差函数
def residuals(x, t, y):
    return model(x, t) - y

# 定义目标函数（总的误差平方和）
def objective_function(x, t, y):
    return np.sum((model(x, t) - y)**2)

# 使用数值梯度进行梯度计算
def gradient(x, t, y):
    eps = np.sqrt(np.finfo(float).eps)
    return approx_fprime(x, objective_function, eps, t, y)

# 最速下降法
def gradient_descent(x_init, t, y, alpha, max_iterations=10000, tolerance=1e-6):
    x = x_init
    for i in range(max_iterations):
        grad = gradient(x, t, y)
        x_new = x - alpha * grad
        
        # 检查收敛性
        if np.linalg.norm(x_new - x) < tolerance:
            break
        x = x_new
    
    # 计算最优参数的损失值
    loss = objective_function(x, t, y)
    
    return x, i, loss

# 初始猜测和学习率
x_init = np.array([0.2, 0.2, 0.1, 0.1])
alpha = 0.1

# 运行最速下降法
x_opt, num_iterations, loss = gradient_descent(x_init, t_data, y_data, alpha)

# 打印结果
print("Number of Iterations:", num_iterations)
print("Optimized Parameters:", x_opt)
print("Final Loss:", loss)

Number of Iterations: 3640
Optimized Parameters: [0.19292142 0.18828542 0.12222641 0.13473105]
Final Loss: 0.00030752241531585274

设置了较小的容忍度和较小的学习率得到如上结果，迭代次数和学习率相关，但是收敛结果基本一致。

牛顿法

import numpy as np
import matplotlib.pyplot as plt

# 拟合点
t_values = np.array([4.0000, 2.0000, 1.0000, 0.5000, 0.2500, 0.1670, 0.1250, 0.1000, 0.0833, 0.0714, 0.0625])
y_values = np.array([0.1957, 0.1947, 0.1735, 0.1600, 0.0844, 0.0627, 0.0456, 0.0342, 0.0323, 0.0235, 0.0246])

# 定义模型函数
def func(x, t):
    return (x[0] * (t**2 + x[1]*t)) / (t**2 + x[2]*t + x[3])

# 目标函数
def loss(x):
    return np.sum((y_values - func(x, t_values))**2)

# 目标函数的梯度
def grad_loss(x):
    epsilon = 1e-9
    grad = np.zeros_like(x)
    for i in range(len(x)):
        x1 = np.copy(x)
        x2 = np.copy(x)
        x1[i] += epsilon
        x2[i] -= epsilon
        grad[i] = (loss(x1) - loss(x2)) / (2 * epsilon)
    return grad

# 数值方法计算Hessian矩阵
def numerical_hessian(f, x, h=1e-5):
    n = len(x)
    H = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            x_ijp = np.array(x, dtype=float)
            x_ijp[i] += h
            x_ijp[j] += h

            x_ip = np.array(x, dtype=float)
            x_ip[i] += h

            x_jp = np.array(x, dtype=float)
            x_jp[j] += h

            x_orig = np.array(x, dtype=float)

            H[i, j] = (f(x_ijp) - f(x_ip) - f(x_jp) + f(x_orig)) / (h**2)
    return H

# 牛顿法
def newton_method(f, grad_f, hessian_f, x0, epsilon=1e-6, max_iter=1000):
    xk = np.array(x0, dtype=float)
    step_count = 0
    for k in range(max_iter):
        grad = grad_f(xk)
        if np.linalg.norm(grad) <= epsilon:
            break
        H = hessian_f(f, xk)
        dk = -np.linalg.solve(H, grad)  # 求解 dk
        xk += dk  # 更新位置
        step_count += 1
    return xk, step_count  # 返回位置和迭代次数

# 初始猜测
x0 = np.array([0.2, 0.1, 0.1, 0.1])

# 运行牛顿法
result, num_iterations = newton_method(loss, grad_loss, numerical_hessian, x0)

print("优化结果:", result)
print("迭代次数:", num_iterations)
print("最终损失:", loss(result))

优化结果: [0.19280693 0.19128247 0.12305651 0.1360624 ]
迭代次数: 5
最终损失: 0.00030750560384928163

牛顿法的求解结果一般，损失值很大，存在奇异矩阵的问题