def f(w, X): return np.sum((X.dot(w)**2)) / len(X) def df_math(w, X): return X.T.dot(X.dot(w)) * 2. / len(X) def df_debug(w, X, epsilon=0.0001): res = np.empty(len(w)) for i in range(len(w)): w_1 = w.copy() w_1[i] += epsilon w_2 = w.copy() w_2[i] -= epsilon res[i] = (f(w_1, X) - f(w_2, X)) / (2 * epsilon) return res def direction(w): return w / np.linalg.norm(w) def gradient_ascent(df, X, initial_w, eta, n_iters = 1e4, epsilon=1e-8): w = direction(initial_w) cur_iter = 0 while cur_iter < n_iters: gradient = df(w, X) last_w = w w = w + eta * gradient w = direction(w) # 注意1:每次求一个单位方向 if(abs(f(w, X) - f(last_w, X)) < epsilon): break cur_iter += 1 return w
时间: 2023-05-26 10:07:42 浏览: 45
这是一个用于实现梯度上升算法的Python函数。其中,f(w, X)用来计算代价函数的值,df_math(w, X)用来计算代价函数的梯度,df_debug(w, X, epsilon)是用来调试df_math(w, X)函数的,direction(w)用来保证梯度向量的方向为最大增加方向,gradient_ascent(df, X, initial_w, eta, n_iters = 1e4, epsilon=1e-8)是梯度上升算法的实现函数。其中,eta为学习率,n_iters为迭代次数,epsilon为收敛判断的临界值。
相关问题
import numpy as np from sklearn import datasets from sklearn.linear_model import LinearRegression np.random.seed(10) class Newton(object): def init(self,epochs=50): self.W = None self.epochs = epochs def get_loss(self, X, y, W,b): """ 计算损失 0.5sum(y_pred-y)^2 input: X(2 dim np.array):特征 y(1 dim np.array):标签 W(2 dim np.array):线性回归模型权重矩阵 output:损失函数值 """ #print(np.dot(X,W)) loss = 0.5np.sum((y - np.dot(X,W)-b)2) return loss def first_derivative(self,X,y): """ 计算一阶导数g = (y_pred - y)*x input: X(2 dim np.array):特征 y(1 dim np.array):标签 W(2 dim np.array):线性回归模型权重矩阵 output:损失函数值 """ y_pred = np.dot(X,self.W) + self.b g = np.dot(X.T, np.array(y_pred - y)) g_b = np.mean(y_pred-y) return g,g_b def second_derivative(self,X,y): """ 计算二阶导数 Hij = sum(X.T[i]X.T[j]) input: X(2 dim np.array):特征 y(1 dim np.array):标签 output:损失函数值 """ H = np.zeros(shape=(X.shape[1],X.shape[1])) H = np.dot(X.T, X) H_b = 1 return H, H_b def fit(self, X, y): """ 线性回归 y = WX + b拟合,牛顿法求解 input: X(2 dim np.array):特征 y(1 dim np.array):标签 output:拟合的线性回归 """ self.W = np.random.normal(size=(X.shape[1])) self.b = 0 for epoch in range(self.epochs): g,g_b = self.first_derivative(X,y) # 一阶导数 H,H_b = self.second_derivative(X,y) # 二阶导数 self.W = self.W - np.dot(np.linalg.pinv(H),g) self.b = self.b - 1/H_bg_b print("itration:{} ".format(epoch), "loss:{:.4f}".format( self.get_loss(X, y , self.W,self.b))) def predict(): """ 需要自己实现的代码 """ pass def normalize(x): return (x - np.min(x))/(np.max(x) - np.min(x)) if name == "main": np.random.seed(2) X = np.random.rand(100,5) y = np.sum(X3 + X**2,axis=1) print(X.shape, y.shape) # 归一化 X_norm = normalize(X) X_train = X_norm[:int(len(X_norm)*0.8)] X_test = X_norm[int(len(X_norm)*0.8):] y_train = y[:int(len(X_norm)0.8)] y_test = y[int(len(X_norm)0.8):] # 牛顿法求解回归问题 newton=Newton() newton.fit(X_train, y_train) y_pred = newton.predict(X_test,y_test) print(0.5np.sum((y_test - y_pred)**2)) reg = LinearRegression().fit(X_train, y_train) y_pred = reg.predict(X_test) print(0.5np.sum((y_test - y_pred)**2)) ——修改代码中的问题,并补全缺失的代码,实现牛顿最优化算法
好的,我根据代码中的问题进行了修改和补全缺失的部分,如下所示:
```
import numpy as np
from sklearn import datasets
from sklearn.linear_model import LinearRegression
np.random.seed(10)
class Newton(object):
def __init__(self, epochs=50):
self.W = None
self.b = None
self.epochs = epochs
def get_loss(self, X, y, W, b):
"""
计算损失 0.5sum(y_pred-y)^2
input:
X(2 dim np.array):特征
y(1 dim np.array):标签
W(2 dim np.array):线性回归模型权重矩阵
b(1 dim np.array):线性回归模型偏置
output:损失函数值
"""
y_pred = np.dot(X, W) + b
loss = 0.5 * np.sum((y_pred - y) ** 2)
return loss
def first_derivative(self, X, y):
"""
计算一阶导数g = (y_pred - y)*x
input:
X(2 dim np.array):特征
y(1 dim np.array):标签
output:损失函数值
"""
y_pred = np.dot(X, self.W) + self.b
g = np.dot(X.T, np.array(y_pred - y))
g_b = np.mean(y_pred - y)
return g, g_b
def second_derivative(self, X, y):
"""
计算二阶导数 H_ij = sum(X.T[i]X.T[j])
input:
X(2 dim np.array):特征
y(1 dim np.array):标签
output:损失函数值
"""
H = np.dot(X.T, X)
H_b = 1
return H, H_b
def fit(self, X, y):
"""
线性回归 y = WX + b拟合,牛顿法求解
input:
X(2 dim np.array):特征
y(1 dim np.array):标签
output:拟合的线性回归
"""
self.W = np.random.normal(size=(X.shape[1]))
self.b = 0
for epoch in range(self.epochs):
g, g_b = self.first_derivative(X, y)
H, H_b = self.second_derivative(X, y)
H_inv = np.linalg.inv(H)
self.W = self.W - np.dot(H_inv, g)
self.b = self.b - H_inv @ np.array([g_b])
print("itration:{} ".format(epoch), "loss:{:.4f}".format(self.get_loss(X, y, self.W, self.b)))
def predict(self, X):
"""
线性回归预测
input:
X(2 dim np.array):特征
output:预测结果
"""
y_pred = np.dot(X, self.W) + self.b
return y_pred
def normalize(x):
return (x - np.min(x)) / (np.max(x) - np.min(x))
if __name__ == "__main__":
np.random.seed(2)
X = np.random.rand(100, 5)
y = np.sum(X ** 3 + X ** 2, axis=1)
print(X.shape, y.shape)
# 归一化
X_norm = normalize(X)
X_train = X_norm[:int(len(X_norm) * 0.8)]
X_test = X_norm[int(len(X_norm) * 0.8):]
y_train = y[:int(len(X_norm) * 0.8)]
y_test = y[int(len(X_norm) * 0.8):]
# 牛顿法求解回归问题
newton = Newton()
newton.fit(X_train, y_train)
y_pred = newton.predict(X_test)
print("Newton Method MSE: ", 0.5 * np.sum((y_test - y_pred) ** 2))
# 对比sklearn中的线性回归
reg = LinearRegression().fit(X_train, y_train)
y_pred = reg.predict(X_test)
print("Linear Regression MSE: ", 0.5 * np.sum((y_test - y_pred) ** 2))
```
在这个修改后的代码中,我主要做了以下几个修改和补充:
1. 在 `__init__` 函数中,添加了 `self.b = None`,用于存储线性回归模型的偏置;
2. 在 `get_loss` 函数中,修正了计算损失的公式,将 `b` 加入到了计算中;
3. 在 `fit` 函数中,对牛顿法的求解进行了修改,使用了矩阵求逆的方式来计算 Hessian 矩阵的逆;
4. 在 `predict` 函数中,添加了 `X` 作为输入参数,用于进行预测;
5. 在主函数中,修改了求解 MSE 的计算公式,以及添加了输出方法,用于输出牛顿法求解和 sklearn 中线性回归的 MSE 值。
希望这些修改能够帮助你。
请在不影响结果的条件下改变代码的样子:import numpy as np import matplotlib.pyplot as plt x1len = 21 x2len = 18 LEN = x1len + x2len POPULATION_SIZE = 100 GENERATIONS = 251 CROSSOVER_RATE = 0.7 MUTATION_RATE = 0.3 pop = np.random.randint(0,2,size=(POPULATION_SIZE,LEN)) def BinToX(pop): x1 = pop[:,0:x1len] x2 = pop[:,x1len:] x1 = x1.dot(2**np.arange(x1len)[::-1]) x2 = x2.dot(2**np.arange(x2len)[::-1]) x1 = -2.9 + x1*(12 + 2.9)/(np.power(2,x1len)-1) x2 = 4.2 + x2*(5.7 - 4.2)/(np.power(2,x2len)-1) return x1,x2 def func(pop): x1,x2 = BinToX(pop) return 21.5 + x1*np.sin(4*np.pi*x1) + x2*np.sin(20*np.pi*x2) def fn(pop): return func(pop); def selection(pop, fitness): idx = np.random.choice(np.arange(pop.shape[0]), size=POPULATION_SIZE, replace=True, p=fitness/fitness.sum()) return pop[idx] def crossover(IdxP1,pop): if np.random.rand() < CROSSOVER_RATE: C = np.zeros((1,LEN)) IdxP2 = np.random.randint(0, POPULATION_SIZE) pt = np.random.randint(0, LEN) C[0,:pt] = pop[IdxP1,:pt] C[0,pt:] = pop[IdxP2, pt:] np.append(pop, C, axis=0) return def mutation(idx,pop): if np.random.rand() < MUTATION_RATE: mut_index = np.random.randint(0, LEN) pop[idx,mut_index] = 1- pop[idx,mut_index] return best_chrom = np.zeros(LEN) best_score = 0 fig = plt.figure() for generation in range(GENERATIONS): fitness = fn(pop) pop = selection(pop, fitness) if generation%50 == 0: ax = fig.add_subplot(2,3,generation//50 +1, projection='3d', title = "generation:"+str(generation)+" best="+str(np.max(fitness))) x1,x2 = BinToX(pop) z = func(pop) ax.scatter(x1,x2,z) for idx in range(POPULATION_SIZE): crossover(idx,pop) mutation(idx,pop) idx = np.argmax(fitness) if best_score < fitness[idx]: best_score = fitness[idx] best_chrom = pop[idx, :] plt.show() print('最优解:', best_chrom, '| best score: %.2f' % best_score)
import numpy as np
import matplotlib.pyplot as plt
x1len = 21
x2len = 18
LEN = x1len + x2len
POPULATION_SIZE = 100
GENERATIONS = 251
CROSSOVER_RATE = 0.7
MUTATION_RATE = 0.3
pop = np.random.randint(0,2,size=(POPULATION_SIZE,LEN))
def BinToX(pop):
x1 = pop[:,0:x1len]
x2 = pop[:,x1len:]
x1 = -2.9 + x1*(12 + 2.9)/(np.power(2,x1len)-1)
x2 = 4.2 + x2*(5.7 - 4.2)/(np.power(2,x2len)-1)
return x1,x2
def func(x1,x2):
return 21.5 + x1*np.sin(4*np.pi*x1) + x2*np.sin(20*np.pi*x2)
def fn(pop):
x1,x2 = BinToX(pop)
return func(x1,x2);
def selection(pop, fitness):
idx = np.random.choice(np.arange(pop.shape[0]),size=POPULATION_SIZE, replace=True, p=fitness/fitness.sum())
return pop[idx]
def crossover(IdxP1,pop):
if np.random.rand() < CROSSOVER_RATE:
C = np.zeros((1,LEN))
IdxP2 = np.random.randint(0, POPULATION_SIZE)
pt = np.random.randint(0, LEN)
C[0,:pt] = pop[IdxP1,:pt]
C[0,pt:] = pop[IdxP2, pt:]
np.append(pop, C, axis=0)
return
def mutation(idx,pop):
if np.random.rand() < MUTATION_RATE:
mut_index = np.random.randint(0, LEN)
pop[idx,mut_index] = 1- pop[idx,mut_index]
return
best_chrom = np.zeros(LEN)
best_score = 0
fig = plt.figure()
for generation in range(GENERATIONS):
fitness = fn(pop)
pop = selection(pop, fitness)
if generation%50 == 0:
ax = fig.add_subplot(2,3,generation//50 +1, projection='3d',title = "generation:"+str(generation)+" best="+str(np.max(fitness)))
x1,x2 = BinToX(pop)
z = func(x1,x2)
ax.scatter(x1,x2,z)
for idx in range(POPULATION_SIZE):
crossover(idx,pop)
mutation(idx,pop)
idx = np.argmax(fitness)
if best_score < fitness[idx]:
best_score = fitness[idx]
best_chrom = pop[idx, :]
plt.show()
print('最优解:', best_chrom, '| best score: %.2f' % best_score)