套索回归和岭回归LASSO and Ridge Regression
这个脚本展示如何使用TensorFlow求解 =+y=Ax+b 的套索回归和岭回归。
我们使用iris数据集,特别地: y = Sepal Length, x = Petal Width
# import required libraries
import matplotlib.pyplot as plt
import sys
import numpy as np
import tensorflow as tf
from sklearn import datasets
from tensorflow.python.framework import ops
regression_type = 'LASSO'
# clear out old graph
ops.reset_default_graph()
#tf.set_random_seed(42)
np.random.seed(42)
# iris.data = [(Sepal Length, Sepal Width, Petal Length, Petal Width)]
iris = datasets.load_iris()
x_vals = np.array([x[3] for x in iris.data])
y_vals = np.array([y[0] for y in iris.data])
def model(x,w,b):
# Declare model operations
model_output = tf.add(tf.matmul(x, w), b)
return model_output
def loss1(x,y,w,b):
# Declare Deming loss function
if regression_type == 'LASSO':
# Declare Lasso loss function
# Lasso Loss = L2_Loss + heavyside_step,
# Where heavyside_step ~ 0 if A < constant, otherwise ~ 99
lasso_param = tf.constant(0.9)
heavyside_step = tf.truediv(1., tf.add(1., tf.exp(tf.multiply(-50., tf.subtract(w, lasso_param)))))
regularization_param = tf.multiply(heavyside_step, 99.)
loss = tf.add(tf.reduce_mean(tf.square(y - model(x,w,b))), regularization_param)
elif regression_type == 'Ridge':
# Declare the Ridge loss function
# Ridge loss = L2_loss + L2 norm of slope
ridge_param = tf.constant(1.)
ridge_loss = tf.reduce_mean(tf.square(w))
loss = tf.expand_dims(tf.add(tf.reduce_mean(tf.square(y - model(x,w,b))), tf.multiply(ridge_param, ridge_loss)), 0)
else:
print('Invalid regression_type parameter value',file=sys.stderr)
return loss
def grad1(x,y,w,b):
with tf.GradientTape() as tape:
loss_1 = loss1(x,y,w,b)
return tape.gradient(loss_1,[w,b])
# Declare batch size
batch_size = 50
# make results reproducible
seed = 13
np.random.seed(seed)
#tf.set_random_seed(seed)
# Declare batch size
learning_rate = 0.25 # Will not converge with learning rate at 0.4
iterations = 50
# Create variables for linear regression
w1 = tf.Variable(tf.random.normal(shape=[1,1]),tf.float32)
b1 = tf.Variable(tf.random.normal(shape=[1,1]),tf.float32)
optimizer = tf.optimizers.Adam(learning_rate)
# Training loop
loss_vec = []
for i in range(5000):
rand_index = np.random.choice(len(x_vals), size=batch_size)
rand_x = np.transpose([x_vals[rand_index]])
rand_y = np.transpose([y_vals[rand_index]])
x=tf.cast(rand_x,tf.float32)
y=tf.cast(rand_y,tf.float32)
grads1=grad1(x,y,w1,b1)
optimizer.apply_gradients(zip(grads1,[w1,b1]))
#sess.run(train_step, feed_dict={x_data: rand_x, y_target: rand_y})
temp_loss1 = loss1(x, y,w1,b1)[0,0].numpy()
#sess.run(loss, feed_dict={x_data: rand_x, y_target: rand_y})
loss_vec.append(temp_loss1)
if (i+1)%25==0:
print('Step #' + str(i+1) + ' A = ' + str(w1.numpy()) + ' b = ' + str(b1.numpy()))
print('Loss = ' + str(temp_loss1))
# Get the optimal coefficients
[slope] = w1.numpy()
[y_intercept] = b1.numpy()
# Get best fit line
best_fit1 = []
for i in x_vals:
best_fit1.append(slope*i+y_intercept)
# Plot the result
plt.plot(x_vals, y_vals, 'o', label='Data Points')
plt.plot(x_vals, best_fit1, 'r-', label='Best fit line', linewidth=3)
plt.legend(loc='upper left')
plt.title('Sepal Length vs Petal Width')
plt.xlabel('Petal Width')
plt.ylabel('Sepal Length')
plt.show()
# Plot loss over time
plt.plot(loss_vec, 'k-')
plt.title('L1 Loss per Generation')
plt.xlabel('Generation')
plt.ylabel('L1 Loss')
plt.show()
