Student_exam_oriented_ex_3_2

!wget --no-cache -O init.py -q https://raw.githubusercontent.com/jdariasl/OTBD/main/content/init.py
import init; init.init(force_download=False)

#!sudo apt install cm-super dvipng texlive-latex-extra texlive-latex-recommended

import numpy as np
from local.lib.data import load_data
import scipy as sc
import matplotlib.pyplot as plt

#!pip install cvxpy

import cvxpy as cp

Exercise

• Algorithm: FISTA

• Problem: LASSO

\(\underset{{\bf{w}}}{\min}f({\bf{w}})=\underset{{\bf{w}}}{\min}\left(\frac{1}{n}\left\Vert {\bf{X}}{\bf{w}}-{\bf{y}}\right\Vert _{2}^{2}+\lambda\left\Vert {\bf{w}}\right\Vert _{1}\right)\)

• Weather forecast dataset

  • 23 features: max / min air temperature, forecast of next day temperatures, precipitations, wind speed…

    • We use only 20 features: \(\bf{X}\) is a \(1000\times 20\) matrix containing 1000 dataset entries.

  • Target: to predict the right next day temperature (bias correction of forecast model)

    • Thus, \({\bf{y}}\) is a \(1000\times1\) vector containing the regression target (i.e., actual next day maximum temperature)

#load data
X,y = load_data("regression", 4)
n,d = X.shape
#data normalisation
X = sc.stats.zscore(X)
y = sc.stats.zscore(y)

# Constant parameters
lamb = 0.1  #regularisation parameter
Niter=600   # Number of iterations for each algorithm

#cvx_solver
def solver_cvx(n,X,Y,lamb,objective_fn):
    n_columns = X.shape[1]
    w = cp.Variable(n_columns)
    lambd = cp.Parameter(nonneg=True)
    lambd.value = lamb
    problem = cp.Problem(
        cp.Minimize(objective_fn(n, X, Y, w, lambd))
    )
    problem.solve()
    return w.value

# Get the optimum value for comparison purposes 
#===================================
loss_fn = lambda n, X, Y, w: (1/n)*cp.pnorm(X @ w - Y, p=2)**2
reg_L1 = lambda w: cp.pnorm(w, p=1)
loss_lasso = lambda n, X, Y, w, lambd: loss_fn(n, X, Y, w) + lambd * reg_L1(w)


# Solution of the empirical risk using CVX
#=========================================
w_lasso_cvx=solver_cvx(n,X,y,lamb,loss_lasso)

w = cp.Variable(w_lasso_cvx.shape[0])
w.value = w_lasso_cvx
f_cvx=loss_lasso(n,X,y,w_lasso_cvx,lamb).value

print(f'The loss function f at the optimum takes the value {f_cvx}')

The loss function f at the optimum takes the value 0.3388340670437767

eta = 0.1  # Constant step size

#Function that estimates the loss for several w at once.
f = lambda n, X, Y, w, lambd: (1/n)*np.sum((X@w - np.kron(Y.reshape(-1,1),np.ones((1,Niter+1))))**2, axis=0) + lambd*np.sum(np.abs(w),axis=0)


# FISTA
w_fista = np.zeros((d,Niter+1))
v_fista = np.zeros((d,Niter+1))

for k in range(Niter):
    #Complete the code including the updating formulas. Keep the weight values for all the iterations
    
    v_fista = ...

f_fista=f(n,X,y,v_fista,lamb)

plt.rcParams.update({
    "text.usetex": True,
    "font.family": "Helvetica"
})
t = range(Niter+1)
plt.plot(t, 10*np.log10((f_fista-f_cvx)**2+np.finfo(float).eps), color = 'b',label = 'FISTA')
plt.grid()
plt.legend()
plt.xlabel('Iteration')
plt.ylabel(r'$10\log(f({\bf{x}})-f*)^2)$ (MSE)')
plt.show()