求信赖域算法加局部二次近似求解 关于b级数i从1到nxi[yi-(e∧bxi/1+e∧bxi)])的最大值(b是P元向量)的Python代码
时间: 2023-12-03 07:46:51 浏览: 91
信赖域算法加局部二次近似求解的Python代码如下:
```python
import numpy as np
from scipy.optimize import minimize
def trust_region(func, grad, hess, x0, radius=1.0, tol=1e-6):
x = x0
fval = func(x)
gradval = grad(x)
hessval = hess(x)
while np.linalg.norm(gradval) > tol:
# solve local quadratic model
p = -np.linalg.solve(hessval, gradval)
# calculate predicted and actual reduction
pred = -0.5 * np.dot(gradval, p)
actual = func(x + p) - fval
# calculate rho
rho = actual / pred
if rho < 0.25:
radius *= 0.25
elif rho > 0.75 and np.linalg.norm(p) == radius:
radius = min(2.0 * radius, 100.0)
if rho > 0.1:
x += p
fval = func(x)
gradval = grad(x)
hessval = hess(x)
if np.linalg.norm(p) < tol or np.linalg.norm(gradval) < tol:
break
return x
def b_series(b, x, y):
n = len(x)
s = 0.0
for i in range(n):
s += np.power(np.e, b.dot(x[i])) * (y[i] - np.power(np.e, np.dot(b, x[i]))) / (1 + np.power(np.e, b.dot(x[i])))
return -s
def b_series_grad(b, x, y):
n = len(x)
s = np.zeros(len(b))
for i in range(n):
ebx = np.power(np.e, b.dot(x[i]))
s += ebx * x[i] * (y[i] - ebx) / np.power(1 + ebx, 2)
return -s
def b_series_hess(b, x, y):
n = len(x)
s = np.zeros((len(b), len(b)))
for i in range(n):
ebx = np.power(np.e, b.dot(x[i]))
s += np.outer(x[i], x[i]) * ebx * (ebx - y[i]) / np.power(1 + ebx, 2)
return -s
# example usage
x = np.array([[1, 2], [3, 4], [5, 6], [7, 8]])
y = np.array([2, 4, 6, 8])
b0 = np.zeros(2)
radius = 1.0
tol = 1e-6
b_opt = trust_region(lambda b: b_series(b, x, y), lambda b: b_series_grad(b, x, y), lambda b: b_series_hess(b, x, y), b0, radius, tol)
print("Optimal b:", b_opt)
print("Optimal value:", -b_series(b_opt, x, y))
```
其中,`trust_region` 函数是信赖域算法的实现,`b_series` 函数是目标函数,`b_series_grad` 函数是目标函数关于变量的梯度,`b_series_hess` 函数是目标函数关于变量的黑塞矩阵。
阅读全文