<jupyter_start><jupyter_text>1. Data processing- Download the Diabete dataset from https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/binary/diabetes- Load the data using sklearn.- Preprocess the data.  1.1. Load the data<jupyter_code>from sklearn import datasets
import numpy

x_sparse, y = datasets.load_svmlight_file('diabetes')
x = x_sparse.todense()

print('Shape of x: ' + str(x.shape))
print('Shape of y: ' + str(y.shape))<jupyter_output>Shape of x: (768, 8)
Shape of y: (768,)<jupyter_text>1.2. Partition to training and test sets<jupyter_code># partition the data to training and test sets
n = x.shape[0]
n_train = int(numpy.ceil(n * 0.8))
n_test = n - n_train

rand_indices = numpy.random.permutation(n)
train_indices = rand_indices[0:n_train]
test_indices = rand_indices[n_train:n]

x_train = x[train_indices, :]
x_test = x[test_indices, :]
y_train = y[train_indices].reshape(n_train, 1)
y_test = y[test_indices].reshape(n_test, 1)

print('Shape of x_train: ' + str(x_train.shape))
print('Shape of x_test: ' + str(x_test.shape))
print('Shape of y_train: ' + str(y_train.shape))
print('Shape of y_test: ' + str(y_test.shape))<jupyter_output>Shape of x_train: (615, 8)
Shape of x_test: (153, 8)
Shape of y_train: (615, 1)
Shape of y_test: (153, 1)<jupyter_text>1.3. Feature scaling Min-max normalization and standardization are two popular feature scaling methods.- Min-max normalization scales the features to the interval $[0, 1]$.- Standardization makes the features have zero mean and unit variance.<jupyter_code># Min-Max Normalization
import numpy

d = x.shape[1]
xmin = numpy.min(x, axis=0).reshape(1, d)
xmax = numpy.max(x, axis=0).reshape(1, d)
xnew = (x - xmin) / (xmax - xmin)

print(xnew)

print('max = ')
print(numpy.max(xnew, axis=0))

print('min = ')
print(numpy.min(xnew, axis=0))
# Standardization
import numpy

d = x.shape[1]
mu = numpy.mean(x, axis=0).reshape(1, d)
sig = numpy.std(x, axis=0).reshape(1, d)
xnew = (x - mu) / sig

print('xnew = ')
print(xnew)

print('mean = ')
print(numpy.mean(xnew, axis=0))

print('std = ')
print(numpy.std(xnew, axis=0))<jupyter_output>xnew = 
[[ 0.63994726  0.84832379  0.14964075 ...  0.20401252  0.46849198
   1.4259954 ]
 [-0.84488505 -1.12339636 -0.16054575 ... -0.68442195 -0.36506078
  -0.19067191]
 [ 1.23388019  1.94372388 -0.26394125 ... -1.10325559  0.60439732
  -0.10558415]
 ...
 [ 0.3429808   0.00330087  0.14964075 ... -0.73518952 -0.68519336
  -0.27575966]
 [-0.84488505  0.1597866  -0.47073225 ... -0.24020459 -0.37110101
   1.17073215]
 [-0.84488505 -0.8730192   0.04624525 ... -0.20212882 -0.47378505
  -0.87137393]]
mean = 
[[-7.74843153e-17  3.61400724e-18 -1.32724416e-17  7.76288755e-17
  -5.49329101e-17  5.12683067e-15  1.92438658e-15  2.19297959e-16]]
std = 
[[1. 1. 1. 1. 1. 1. 1. 1.]]<jupyter_text>In this tutorial, we use the standardization to trainsform both training and test features<jupyter_code># Standardization
import numpy

# calculate mu and sig using the training set
d = x_train.shape[1]
mu = numpy.mean(x_train, axis=0).reshape(1, d)
sig = numpy.std(x_train, axis=0).reshape(1, d)

# transform the training features
x_train = (x_train - mu) / sig

# transform the test features
x_test = (x_test - mu) / sig


print('test mean = ')
print(numpy.mean(x_test, axis=0))

print('test std = ')
print(numpy.std(x_test, axis=0))<jupyter_output>test mean = 
[[ 0.04282502  0.14644089  0.13515084  0.07852463  0.12107582  0.05411205
  -0.04668225  0.01295459]]
test std = 
[[0.98823714 1.05366321 1.01518686 0.97945751 1.12619899 1.01806036
  0.82103702 0.86354754]]<jupyter_text>Add a dimension of all ones<jupyter_code>n_train, d = x_train.shape
x_train = numpy.concatenate((x_train, numpy.ones((n_train, 1))), axis=1)

n_test, d = x_test.shape
x_test = numpy.concatenate((x_test, numpy.ones((n_test, 1))), axis=1)

print('Shape of x_train: ' + str(x_train.shape))
print('Shape of x_test: ' + str(x_test.shape))<jupyter_output>Shape of x_train: (615, 9)
Shape of x_test: (153, 9)<jupyter_text>2. Logistic regression modelThe objective function is $Q (w; X, y) = \frac{1}{n} \sum_{i=1}^n \log \Big( 1 + \exp \big( - y_i x_i^T w \big) \Big) + \frac{\lambda}{2} \| w \|_2^2 $.<jupyter_code># Calculate the objective function value
# Inputs:
#     w: d-by-1 matrix
#     x: n-by-d matrix
#     y: n-by-1 matrix
#     lam: scalar, the regularization parameter
# Return:
#     objective function value (scalar)
def objective(w, x, y, lam):
    n, d = x.shape
    yx = numpy.multiply(y, x) # n-by-d matrix
    yxw = numpy.dot(yx, w) # n-by-1 matrix
    vec1 = numpy.exp(-yxw) # n-by-1 matrix
    vec2 = numpy.log(1 + vec1) # n-by-1 matrix
    loss = numpy.mean(vec2) # scalar
    reg = lam / 2 * numpy.sum(w * w) # scalar
    return loss + reg
    
# initialize w
d = x_train.shape[1]
w = numpy.zeros((d, 1))

# evaluate the objective function value at w
lam = 1E-6
objval0 = objective(w, x_train, y_train, lam)
print('Initial objective function value = ' + str(objval0))<jupyter_output>Initial objective function value = 0.6931471805599453<jupyter_text>3. Numerical optimization  3.1. Calculate the full gradientThe gradient at $w$ is $g = - \frac{1}{n} \sum_{i=1}^n \frac{y_i x_i }{1 + \exp ( y_i x_i^T w)} + \lambda w$<jupyter_code># Calculate the gradient
# Inputs:
#     w: d-by-1 matrix
#     x: n-by-d matrix
#     y: n-by-1 matrix
#     lam: scalar, the regularization parameter
# Return:
#     g: g: d-by-1 matrix, full gradient
def gradient(w, x, y, lam):
    n, d = x.shape
    yx = numpy.multiply(y, x) # n-by-d matrix
    yxw = numpy.dot(yx, w) # n-by-1 matrix
    vec1 = numpy.exp(yxw) # n-by-1 matrix
    vec2 = numpy.divide(yx, 1+vec1) # n-by-d matrix
    vec3 = -numpy.mean(vec2, axis=0).reshape(d, 1) # d-by-1 matrix
    g = vec3 + lam * w
    return g<jupyter_output><empty_output><jupyter_text>3.2. Gradient descent<jupyter_code># Gradient descent for solving logistic regression
# Inputs:
#     x: n-by-d matrix
#     y: n-by-1 matrix
#     lam: scalar, the regularization parameter
#     stepsize: scalar
#     max_iter: integer, the maximal iterations
#     w: d-by-1 matrix, initialization of w
# Return:
#     w: d-by-1 matrix, the solution
#     objvals: a record of each iteration's objective value
def grad_descent(x, y, lam, stepsize, max_iter=100, w=None):
    n, d = x.shape
    objvals = numpy.zeros(max_iter) # store the objective values
    if w is None:
        w = numpy.zeros((d, 1)) # zero initialization
    
    for t in range(max_iter):
        objval = objective(w, x, y, lam)
        objvals[t] = objval
        print('Objective value at t=' + str(t) + ' is ' + str(objval))
        g = gradient(w, x, y, lam)
        w -= stepsize * g
    
    return w, objvals
# example
lam = 1E-6
stepsize = 1.0
w, objvals_gd = grad_descent(x_train, y_train, lam, stepsize)<jupyter_output>Objective value at t=0 is 0.6931471805599453
Objective value at t=1 is 0.5857422010224674
Objective value at t=2 is 0.5404055860105454
Objective value at t=3 is 0.5151957841771834
Objective value at t=4 is 0.4990843613352312
Objective value at t=5 is 0.4879904229793065
Objective value at t=6 is 0.4799953399338087
Objective value at t=7 is 0.4740522631313984
Objective value at t=8 is 0.46953259598819874
Objective value at t=9 is 0.46603377304313864
Objective value at t=10 is 0.4632859399675407
Objective value at t=11 is 0.4611019171637835
Objective value at t=12 is 0.45934833927905316
Objective value at t=13 is 0.4579280772961021
Objective value at t=14 is 0.45676907664354227
Objective value at t=15 is 0.4558170315694747
Objective value at t=16 is 0.4550304454901687
Objective value at t=17 is 0.45437722328184754
Objective value at t=18 is 0.453832273543992
Objective value at t=19 is 0.4533757919269716
Objective value at t=20 is 0.4529920128789861
Objective value at t=21 is 0.45266828925[...]<jupyter_text>3.3. Stochastic gradient descent (SGD)Define $Q_i (w) = \log \Big( 1 + \exp \big( - y_i x_i^T w \big) \Big) + \frac{\lambda}{2} \| w \|_2^2 $.The stochastic gradient at $w$ is $g_i = \frac{\partial Q_i }{ \partial w} = -\frac{y_i x_i }{1 + \exp ( y_i x_i^T w)} + \lambda w$.<jupyter_code># Calculate the objective Q_i and the gradient of Q_i
# Inputs:
#     w: d-by-1 matrix
#     xi: 1-by-d matrix
#     yi: scalar
#     lam: scalar, the regularization parameter
# Return:
#     obj: scalar, the objective Q_i
#     g: d-by-1 matrix, gradient of Q_i
def stochastic_objective_gradient(w, xi, yi, lam):
    d = xi.shape[0]
    yx = yi * xi # 1-by-d matrix
    yxw = float(numpy.dot(yx, w)) # scalar
    
    # calculate objective function Q_i
    loss = numpy.log(1 + numpy.exp(-yxw)) # scalar
    reg = lam / 2 * numpy.sum(w * w) # scalar
    obj = loss + reg
    
    # calculate stochastic gradient
    g_loss = -yx.T / (1 + numpy.exp(yxw)) # d-by-1 matrix
    g = g_loss + lam * w # d-by-1 matrix
    
    return obj, g
# SGD for solving logistic regression
# Inputs:
#     x: n-by-d matrix
#     y: n-by-1 matrix
#     lam: scalar, the regularization parameter
#     stepsize: scalar
#     max_epoch: integer, the maximal epochs
#     w: d-by-1 matrix, initialization of w
# Return:
#     w: the solution
#     objvals: record of each iteration's objective value
def sgd(x, y, lam, stepsize, max_epoch=100, w=None):
    n, d = x.shape
    objvals = numpy.zeros(max_epoch) # store the objective values
    if w is None:
        w = numpy.zeros((d, 1)) # zero initialization
    
    for t in range(max_epoch):
        # randomly shuffle the samples
        rand_indices = numpy.random.permutation(n)
        x_rand = x[rand_indices, :]
        y_rand = y[rand_indices, :]
        
        objval = 0 # accumulate the objective values
        for i in range(n):
            xi = x_rand[i, :] # 1-by-d matrix
            yi = float(y_rand[i, :]) # scalar
            obj, g = stochastic_objective_gradient(w, xi, yi, lam)
            objval += obj
            w -= stepsize * g
        
        stepsize *= 0.9 # decrease step size
        objval /= n
        objvals[t] = objval
        print('Objective value at epoch t=' + str(t) + ' is ' + str(objval))
    
    return w, objvals
# example
lam = 1E-6
stepsize = 0.1
w, objvals_sgd = sgd(x_train, y_train, lam, stepsize)<jupyter_output>Objective value at epoch t=0 is 0.5151549707339808
Objective value at epoch t=1 is 0.5206758924832361
Objective value at epoch t=2 is 0.5043274140498304
Objective value at epoch t=3 is 0.5017007852231665
Objective value at epoch t=4 is 0.5061405326180962
Objective value at epoch t=5 is 0.4999228760733248
Objective value at epoch t=6 is 0.49963857212327895
Objective value at epoch t=7 is 0.49163172119009857
Objective value at epoch t=8 is 0.493069056246553
Objective value at epoch t=9 is 0.4907411452841133
Objective value at epoch t=10 is 0.4863129094437151
Objective value at epoch t=11 is 0.48762427080878673
Objective value at epoch t=12 is 0.484085909137715
Objective value at epoch t=13 is 0.48001497432115775
Objective value at epoch t=14 is 0.4817474833510589
Objective value at epoch t=15 is 0.4793509087801017
Objective value at epoch t=16 is 0.4767657378561972
Objective value at epoch t=17 is 0.4771544041110307
Objective value at epoch t=18 is 0.4771127791611159
Objective value at e[...]<jupyter_text>Compare GD with SGDPlot objective function values against epochs.<jupyter_code>import matplotlib.pyplot as plt
%matplotlib inline

fig = plt.figure(figsize=(6, 4))

epochs_gd = range(len(objvals_gd))
epochs_sgd = range(len(objvals_sgd))

line0, = plt.plot(epochs_gd, objvals_gd, '--b', LineWidth=4)
line1, = plt.plot(epochs_sgd, objvals_sgd, '-r', LineWidth=2)
plt.xlabel('Epochs', FontSize=20)
plt.ylabel('Objective Value', FontSize=20)
plt.xticks(FontSize=16)
plt.yticks(FontSize=16)
plt.legend([line0, line1], ['GD', 'SGD'], fontsize=20)
plt.tight_layout()
plt.show()
fig.savefig('compare_gd_sgd.pdf', format='pdf', dpi=1200)<jupyter_output><empty_output><jupyter_text>3.3. Mini-batch SGDDefine $Q_I (w) = \frac{1}{b} \sum_{i \in I} \log \Big( 1 + \exp \big( - y_i x_i^T w \big) \Big) + \frac{\lambda}{2} \| w \|_2^2 $, where $I$ is a set containing $b$ indices randomly drawn from $\{ 1, \cdots , n \}$ without replacement.The stochastic gradient at $w$ is $g_I = \frac{\partial Q_I }{ \partial w} = \frac{1}{b} \sum_{i \in I} \frac{- y_i x_i }{1 + \exp ( y_i x_i^T w)} + \lambda w$.  4. Prediction<jupyter_code># Predict class label
# Inputs:
#     w: d-by-1 matrix
#     X: m-by-d matrix
# Return:
#     f: m-by-1 matrix, the predictions
def predict(w, X):
    xw = numpy.dot(X, w)
    f = numpy.sign(xw)
    return f
# evaluate training error
f_train = predict(w, x_train)
diff = numpy.abs(f_train - y_train) / 2
error_train = numpy.mean(diff)
print('Training classification error is ' + str(error_train))
# evaluate test error
f_test = predict(w, x_test)
diff = numpy.abs(f_test - y_test) / 2
error_test = numpy.mean(diff)
print('Test classification error is ' + str(error_test))<jupyter_output>Test classification error is 0.28104575163398693
{'fraction_non_alphanumeric': 0.13045486030939737, 'fraction_numerical': 0.11944893404140691, 'mean_word_length': 3.588276836158192, 'pattern_counts': {'":': 0, '<': 36, '<?xml version=': 0, '>': 36, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '404056', 'license': 'MIT', 'n_tokens_mistral': 6002, 'n_tokens_neox': 5093, 'n_words': 1395, 'path': 'LectureNotes/Logistic/code/Logistic.ipynb', 'repo_name': 'RachelZ82/DeepLearning', 'size': '12993'}