<jupyter_start><jupyter_text>AMATH 515 Homework 1**Due Date: 01/22/2020***Homework Instruction*: Please follow order of this notebook and fill in the codes where commented as `TODO`.<jupyter_code># TODO: Put your name and UW ID here (for homework grading purposes)
UW_ID = "1772371"
FIRST_NAME = "Philip"
LAST_NAME = "Pham"
# require numpy module
import numpy as np
from numpy.linalg import norm
from numpy.linalg import solve
import matplotlib.pyplot as plt
# import the supplementary function for homework 1
import sys
sys.path.insert(0, './')
from hw1_supp import *<jupyter_output><empty_output><jupyter_text>Gradient Descent SolverRecall the gradient descent algorithm we learned from the class and complete the gradient descent solver.<jupyter_code>def optimizeWithGD(x0, func, grad, step_size, tol=1e-6, max_iter=1000, use_line_search=False):
    """
    Optimize with Gradient Descent
    
    input
    -----
    x0 : array_like
        Starting point for the solver.
    func : function
        Input x and return the function value.
    grad : function
        Input x and return the gradient.
    step_size : float or None
        If it is a float number and `use_line_search=False`, it will be used as the step size.
        Otherwise, line search will be used
    tol : float, optional
        Gradient tolerance for terminating the solver.
    max_iter : int, optional
        Maximum number of iteration for terminating the solver.
    use_line_search : bool, optional
        When it is true a line search will be used, other wise `step_size` has to be provided.
        
    output
    ------
    x : array_like
        Final solution
    obj_his : array_like
        Objective function value convergence history
    err_his : array_like
        Norm of gradient convergence history
    exit_flag : int
        0, norm of gradient below `tol`
        1, exceed maximum number of iteration
        2, line search fail
        3, others
    """
    # safeguard
    if not use_line_search and step_size is None:
        print('Please specify the step_size or use the line search.')
        return x0, np.array([]), np.array([]), 3
    
    # initial step
    x = np.copy(x0)
    g = grad(x)
    #
    obj = func(x)
    err = norm(g)
    #
    obj_his = np.zeros(max_iter + 1)
    err_his = np.zeros(max_iter + 1)
    #
    obj_his[0] = obj
    err_his[0] = err
    
    # start iterations
    iter_count = 0
    while err >= tol:
        if use_line_search:
            step_size = lineSearch(x, g, g, func)
        #
        # if line search fail step_size will be None
        if step_size is None:
            print('Gradient descent line search fail.')
            return x, obj_his[:iter_count+1], err_his[:iter_count+1], 2
        #
        # gradient descent step
        x -= step_size * g

        # update function and gradient
        g = grad(x)
        #
        obj = func(x)
        err = norm(g)
        #
        iter_count += 1
        obj_his[iter_count] = obj
        err_his[iter_count] = err
        #
        # check if exceed maximum number of iteration
        if iter_count >= max_iter:
            print('Gradient descent reach maximum number of iteration.')
            return x, obj_his[:iter_count+1], err_his[:iter_count+1], 1
    #
    return x, obj_his[:iter_count+1], err_his[:iter_count+1], 0<jupyter_output><empty_output><jupyter_text>Newton's SolverRecall the Newton's method we learned from the class and complete the Newton's solver.*Remark*: This is a simplified version of Newton's method which do not evolve line search with constant step size 1.<jupyter_code>def optimizeWithNT(x0, func, grad, hess, tol=1e-6, max_iter=100):
    """
    Optimize with Newton's Method
    
    input
    -----
    x0 : array_like
        Starting point for the solver.
    func : function
        Input x and return the function value.
    grad : function
        Input x and return the gradient.
    hess : function
        Input x and return the Hessian matrix.
    tol : float, optional
        Gradient tolerance for terminating the solver.
    max_iter : int, optional
        Maximum number of iteration for terminating the solver.
        
    output
    ------
    x : array_like
        Final solution
    obj_his : array_like
        Objective function value convergence history
    err_his : array_like
        Norm of gradient convergence history
    exit_flag : int
        0, norm of gradient below `tol`
        1, exceed maximum number of iteration
        2, others
    """
    # initial step
    x = np.copy(x0)
    g = grad(x)
    H = hess(x)
    #
    obj = func(x)
    err = norm(g)
    #
    obj_his = np.zeros(max_iter + 1)
    err_his = np.zeros(max_iter + 1)
    #
    obj_his[0] = obj
    err_his[0] = err
    
    # start iteration
    iter_count = 0
    while err >= tol:
        # Newton's step
        #####
        x -= solve(H, g)

        # update function, gradient and Hessian
        g = grad(x)
        H = hess(x)
        #
        obj = func(x)
        err = norm(g)
        #
        iter_count += 1
        obj_his[iter_count] = obj
        err_his[iter_count] = err
        #
        # check if exceed maximum number of iteration
        if iter_count >= max_iter:
            print('Gradient descent reach maximum number of iteration.')
            return x, obj_his[:iter_count+1], err_his[:iter_count+1], 1
    #
    return x, obj_his[:iter_count+1], err_his[:iter_count+1], 0<jupyter_output><empty_output><jupyter_text>Simple TestConsider a simple test problem,$$\min_x~~f(x) \triangleq \frac{1}{2} \|x - b\|^2$$We could easily calculate the gradient and Hessian of $f$,$$\nabla f(x) = x - b, \quad \nabla^2 f(x) = I$$where $I$ is the identity matrix, and we know that $x^* = b$ is the solution.Please complete the solvers above, use the following function as a simple test.*Remark*: Having simple test is always a good idea. It can easily detect the bugs in the code and save a lot of time and emotion energy when use the code in the more complex context.<jupyter_code># create b
b = np.array([1.0, 2.0, 3.0])
# define test function
def test_func(x):
    return 0.5*sum((x - b)**2)
# define test gradient
def test_grad(x):
    return x - b
# define test Hessian
def test_hess(x):
    return np.eye(b.size)
# test gradient descent
x0_gd = np.zeros(b.size)
#
x_gd, obj_his_gd, err_his_gd, exit_flag_gd = optimizeWithGD(x0_gd, test_func, test_grad, 1.0)
# check if solution is correct
err_gd = norm(x_gd - b)
#
if err_gd < 1e-6:
    print('Gradient Descent: OK')
else:
    print('Gradient Descent: Err')
# test Newton's method
x0_nt = np.zeros(b.size)
#
x_nt, obj_his_nt, err_his_nt, exit_flag_nt = optimizeWithNT(x0_nt, test_func, test_grad, test_hess)
# check if solution is correct
err_nt = norm(x_nt - b)
#
if err_nt < 1e-6:
    print('Newton: OK')
else:
    print('Newont: Err')<jupyter_output>Newton: OK<jupyter_text>Logistic RegressionRecall the logistic regression model from the class,$$f(x) = \sum_{i=1}^m\left\{\log\left[1+\exp(\langle a_i,x\rangle)\right]-b_i\langle a_i,x\rangle\right\} + \frac{\lambda}{2}\|x\|^2.$$Calculate the gradient and Hessian of the function and define them below using the given data.<jupyter_code># fix the random seed to generate same data set
np.random.seed(123)
# set dimension and create the data
m_lgt = 500
n_lgt = 50
A_lgt = 0.3*np.random.randn(m_lgt, n_lgt)
x_lgt = np.random.randn(n_lgt)
b_lgt = sampleLGT(x_lgt, A_lgt)
lam_lgt = 0.1<jupyter_output><empty_output><jupyter_text>Define functions<jupyter_code># define function, gradient and Hessian
def lgt_func(x):
    logits = np.matmul(A_lgt, x)
    return np.sum(np.log(1 + np.exp(logits)) - b_lgt * logits) + np.sum(np.square(x)) * lam_lgt / 2.
#
def lgt_grad(x):
    logits = np.matmul(A_lgt, x)
    exp_logits = np.exp(logits)
    return np.sum(A_lgt.T * exp_logits / (1 + exp_logits), axis=1) - np.matmul(A_lgt.T, b_lgt) + lam_lgt * x
#
def lgt_hess(x):
    logits = np.matmul(A_lgt, x)
    exp_logits = np.exp(logits)
    rescaled_A_T = A_lgt.T * np.sqrt(exp_logits) / (1 + exp_logits)
    return np.matmul(rescaled_A_T, rescaled_A_T.T) + np.eye(n_lgt) * lam_lgt
# uniform step size for gradient descent or should we use line search?
step_size_lgt = 1./np.sort(np.linalg.eig(np.matmul(A_lgt.T, A_lgt) / 4. + np.eye(n_lgt) * lam_lgt)[0])[-1]
use_line_search = False
step_size_lgt<jupyter_output><empty_output><jupyter_text>Apply gradient descent solver<jupyter_code>x0_lgt_gd = np.zeros(n_lgt)
x_lgt_gd, obj_his_lgt_gd, err_his_lgt_gd, exit_flag_lgt_gd = \
    optimizeWithGD(x0_lgt_gd, lgt_func, lgt_grad, step_size_lgt, use_line_search=use_line_search)
# plot result
fig, ax = plt.subplots(1, 2, figsize=(12,5))
ax[0].plot(obj_his_lgt_gd)
ax[0].set_title('function value')
ax[1].semilogy(err_his_lgt_gd)
ax[1].set_title('norm of gradient')
fig.suptitle('Gradient Descent on Logistic Regression')
# fig.savefig('logistic_gradient_descent.pdf')<jupyter_output><empty_output><jupyter_text>Apply Newton's solver<jupyter_code>x0_lgt_nt = np.zeros(n_lgt)
x_lgt_nt, obj_his_lgt_nt, err_his_lgt_nt, exit_flag_lgt_nt = \
    optimizeWithNT(x0_lgt_nt, lgt_func, lgt_grad, lgt_hess)
# plot result
fig, ax = plt.subplots(1, 2, figsize=(12,5))
ax[0].plot(obj_his_lgt_nt)
ax[0].set_title('function value')
ax[1].semilogy(err_his_lgt_nt)
ax[1].set_title('norm of gradient')
fig.suptitle('Newton\'s Method on Logistic Regression')
# fig.savefig('logistic_newton.pdf')<jupyter_output><empty_output><jupyter_text>Poisson RegressionRecall the logistic regression model from the class,$$f(x)=\sum_{i=1}^m\left[\exp(\langle a_i,x\rangle)-b_i\langle a_i,x\rangle\right]+\frac{\lambda}{2}\|x\|^2$$Calculate the gradient and Hessian of the function and define them below using the given data.<jupyter_code># fix the random seed to generate same data set
np.random.seed(123)
# set dimension and create the data
m_psn = 500
n_psn = 50
A_psn = 0.3*np.random.randn(m_psn, n_psn)
x_psn = np.random.randn(n_psn)
b_psn = samplePSN(x_psn, A_psn)
lam_psn = 0.1<jupyter_output><empty_output><jupyter_text>Define functions<jupyter_code># define function, gradient and Hessian
def psn_func(x):
    logits = np.matmul(A_psn, x)
    return np.sum(np.exp(logits)) - np.matmul(b_psn, np.matmul(A_psn, x)) + lam_psn * np.sum(np.square(x)) / 2.
#
def psn_grad(x):
    logits = np.matmul(A_psn, x)
    return np.sum(A_psn.T * np.exp(logits), axis=1) - np.matmul(A_psn.T, b_psn) + lam_psn * x
#
def psn_hess(x):
    logits = np.matmul(A_psn, x)
    A_psn_scaled_T = A_psn.T * np.exp(logits / 2)
    return np.matmul(A_psn_scaled_T, A_psn_scaled_T.T) + np.eye(n_psn) * lam_psn
# uniform step size for gradient descent or should we use line search?
step_size_psn = None
use_line_search = True<jupyter_output><empty_output><jupyter_text>Apply gradient descent solver<jupyter_code>x0_psn_gd = np.zeros(n_lgt)
x_psn_gd, obj_his_psn_gd, err_his_psn_gd, exit_flag_psn_gd = \
    optimizeWithGD(x0_psn_gd, psn_func, psn_grad, None, use_line_search=use_line_search)
# plot result
fig, ax = plt.subplots(1, 2, figsize=(12,5))
ax[0].plot(obj_his_psn_gd)
ax[0].set_title('function value')
ax[1].semilogy(err_his_psn_gd)
ax[1].set_title('norm of gradient')
fig.suptitle('Gradient Descent on Poisson Regression')
# fig.savefig('poisson_gradient_descent.pdf')<jupyter_output><empty_output><jupyter_text>Apply Newton's solver<jupyter_code>x0_psn_nt = np.zeros(n_lgt)
x_psn_nt, obj_his_psn_nt, err_his_psn_nt, exit_flag_psn_nt = \
    optimizeWithNT(x0_psn_nt, psn_func, psn_grad, psn_hess)
# plot result
fig, ax = plt.subplots(1, 2, figsize=(12,5))
ax[0].plot(obj_his_psn_nt)
ax[0].set_title('function value')
ax[1].semilogy(err_his_psn_nt)
ax[1].set_title('norm of gradient')
fig.suptitle('Newton\'s Method on Poisson Regression')
# fig.savefig('poisson_newton.pdf')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11042529804684197, 'fraction_numerical': 0.015726726980637526, 'mean_word_length': 3.4871016691957513, 'pattern_counts': {'":': 0, '<': 50, '<?xml version=': 0, '>': 52, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '429464', 'license': 'MIT', 'n_tokens_mistral': 4483, 'n_tokens_neox': 4290, 'n_words': 1122, 'path': 'hw1/515Hw1_Coding.ipynb', 'repo_name': 'ppham27/amath515', 'size': '11827'}