# DATA 442/642 - Advanced Machine Learning Midterm
## [Yunting Chiu](https://www.linkedin.com/in/yuntingchiu/) 

# Exercise 1 (10 points)
Show that the $\ell 1$ norm is a convex function (as all norms), yet it is not strictly convex. In contrast, show that the squared Euclidean norm is a strictly convex function.

## Answer
A function $f$ is **convex** if
$$
f(\lambda x + (1 − \lambda )y) \leq \lambda f(x) + (1 − \lambda )f(y)
$$
for all x, y $\in$ dom f and all $\lambda \in$ [0, 1].

If the inequality holds strictly (i.e. $<$ rather than $\leq$) for all $\lambda \in$ (0, 1) and x not equal to y, then we say that $f$ is **strictly convex**.\
The norm of a vector is a generalization of the concept of length. The $l_1$ norm of a vector $x \in R^n$ is defined as:
$$
||x||_1 := \sum_{i=1}^{n}|x_i|
$$
According to Minkowski's Inequality, we have the following inequality for $l_1$ norm:
$$
\sum_{i}^{l}(|\lambda x_i + (1- \lambda)y_i|) \leq \lambda \sum_{i}^{l}|X_i| + (1 - \lambda) \sum_{i}^{l}|y_i|
$$
which is not strictly convex. All $x_i$ and all $y_i$ are non-negative then there is no reason for strict inequality.

By constract, the squared Euclidean norm is a strictly convex function because of the follwing inequality 
$$
\sum_{i}(|\lambda x_i + (1-\lambda)y_i|)^2 < \lambda \sum_{i}|X_i|^2+ (1 - \lambda) \sum_{i}|y_i|^2
$$
which gives
$$
\lambda ^2 \sum_{i}x^2_i + (1-\lambda)^2\sum_{i}y^2_i + 2\lambda (1-\lambda)x_iy_i  < \lambda \sum_{i}x^2_i + (1-\lambda)^2\sum_{i}y^2_i
$$
then
$$
2 \lambda(1- \lambda) x_iy_i < \lambda(1-\lambda)(\sum_{i}x^2_i + \sum_{i}y^2_i)
$$
So we have 
$$
\sum_{i}(x_i-y_i) > 0, for x \neq y
$$
Which is strictly convex.

# Exercise 2 (10 points)
Let the observations resulting from an experiment be $x_n, n = 1, 2,..., N$. Assume that they are independent and that they originate from a Gaussian PDF with mean $\mu$ and standard deviation $\sigma$. Both, the mean and the variance, are unknown. Prove that the maximum likelihood (ML) estimates of these quantities are given 

$$
\hat{\mu}ML = \frac{1}{N}\sum_{n=1}^{N}X_n,\space\hat{\sigma}^2ML =  \frac{1}{N}\sum_{n=1}^{N}X_n(X_n - \hat{\mu}ML)^2
$$
## Answer
In statistics, maximum likelihood estimation (MLE) is a method of estimating the parameters of an assumed probability distribution. The log-likelihood function is given by
$$
L(\mu, \sigma^2) = - \frac{N}{2}ln(2 \pi) - \frac{N}{2}ln(\sigma^2) - \frac{1}{2\sigma^2}\sum_{n=1}^{N}(x_n - \mu)^2
$$
Why we start with log-likehood? Because we can simply the notation. Taking the derivative with respect to $\mu, \sigma^2$, and equating it to zero we can get the following equations. We have two unknowns, so we need to derivate two equations.
$$
 \frac{\partial L}{\partial \mu} = \frac{1}{\sigma^2}\sum_{n=1}^{N}(x_n - \mu) = \frac{1}{\sigma^2}(\sum_{n=1}^{N}x_n \sum_{n=1}^{N} \mu)=0
$$
$$
 \frac{\partial L}{\partial \sigma^2} = -\frac{N}{2}\frac{1}{\sigma^2} + \frac{1}{2\sigma^4}\cdot \sum_{n=1}^{N}(x_n - \mu)^2 = 0
$$
We know the $\sigma^4$ is $\sigma^2 \cdot \sigma^2$, so we can slove it in the above equation.

# Exercise 3 (15 points)
For the regression model where the noise vector $\eta= [ \, \eta_1,..., \eta_N	] \,^T$ comprises samples from zero mean Gaussian random variable, with covariance matrix $\Sigma_n$, show that the Fisher information matrix is given by

$$
I(\theta) = X^T\Sigma^{-1}_{n}X,
$$
where $X$ is the input matrix.

## Answer
A Gaussian distribution is defined as
$$
P(y, \theta) = \frac{1}{2 \pi^\frac{N}{2}\cdot \sqrt{det(\Sigma_{n}})} \cdot{exp(-\frac{1}{2}(y - x \cdot \theta)^T \cdot \Sigma_{n}(y - x \cdot \theta)}
$$
We take natural log of the equation and get
$$
ln (p(y, \theta)) = -\frac{1}{2}(y - x \cdot \theta)^T \cdot \Sigma_{n}^{-1}(y - x \cdot \theta) + C
$$
where $C$ is a constant.

We take the first derivative from the above equation:
$$
\frac{\partial ln (p(y, \theta))}{\partial \theta} = -\frac{1}{\theta^2}(x - \mu)
$$
Then, we take the second derivative
$$
\frac{\partial^2 ln (p(y, \theta))}{\partial \theta^2} = - x^T \cdot \Sigma_{n}^{-1} \cdot x
$$

Therefore, the Fisher information matrix is defined as follows:
$$
I(\theta) = - E(\frac{\sigma^2\ell_n P (y, \theta)}{\sigma \theta^2}) = x^T \cdot \Sigma_{n}^{-1} \cdot x
$$


# Exercise 4 (20 points)
Consider the regression problem described in one of our labs. Read the same audio file, then add white Gaussian noise at a 15 dB level and randomly “hit” 10% of the data samples with outliers (set the outlier values to 80% of the maximum value of the data samples).


Recall from Lab 6 of the kernal ridge regression


```python
from sklearn.svm import SVR
from datetime import datetime
import numpy as np
import math
import soundfile as sf
import matplotlib.pyplot as plt
import sys
import os
path = "/content/drive/MyDrive/American_University/2021_Fall/DATA-642-001_Advanced Machine Learning/GitHub/Labs/06/BladeRunner.wav"
```


```python
# This is a helping funvtion that contains different types of kernels
def kappa(x, y, kernel_type, kernel_params):
  value = None
  if kernel_type == 'gaus':
    sigma = kernel_params[0]
    N = len(x)
    norm = sum((x-y)**2)
    value = np.exp(-norm/(sigma**2))
  elif kernel_type == 'gaus_c':
    sigma = kernel_params[0]
    N = len(x)
    exponent = sum( (x-y.conj())**2 )
    # value = 2*(math.exp( -np.real(exponent)/(sigma**2) ))
    value = 2*np.real(np.exp(-exponent/(sigma**2)))
  elif kernel_type == 'linear':
    value = 0
    N = len(x)
    for i in range(0, N):
      value = value + x[i]*y[i].conj()
  elif kernel_type == 'poly':
    d = kernel_params[0]
    value = (1 + x*y.transpose())**d
    # value = ( (1 + x*y.transpose())/( math.sqrt(np.real(x*x.transpose())* np.real(y*y.transpose()) ) ) )**d;
  elif kernel_type == 'poly_c':
    d = kernel_params[0]
    value = 2*np.real((1 + x*y.transpose())**d)
  # value = 2*np.real( ( (1 + x*y.transpose())/( math.sqrt(np.real(x*x.transpose()) * np.real(y*y.transpose()) ) ) )**d )
  return value
```


```python
def kernel_regression_l2_l2_unbiased(X, y, params, kernel_type, kernel_params):
  lmbda = params[0]
  # build kernel matrix
  [d, N] = X.shape
  # for n=1:N
  #   for m=1:N
  #     K(n,m) = kappa(X(:,n), X(:,m), kernel_type, kernel_params);
  #   end;
  # end;
  if kernel_type == 'gaus':
    par = kernel_params
    norms = np.zeros(shape=(N, N))
    for i in range(0, N):
      T = X - X[:, i] # bsxfun(@minus,X,X(:,i))
      norms[i, :] = np.sum(T**2, axis=0)
    K = np.exp(-norms/(par**2))
  elif kernel_type == 'gaus_c':
    par = kernel_params
    norms = np.zeros(shape=(N, N))
    for i in range(0, N):
      T = X - X[:, i].conj() # bsxfun(@minus,X,conj(X(:,i)))
      norms[i, :] = np.sum(T**2, axis=0)
    K = 2*np.real(np.exp(-norms/(par**2)))
  else:
    K = np.zeros(shape=(N, N))
    for i in range(0, N):
      for j in range(0, N):
        K[i, j] = kappa(X[:, i], X[:, j], kernel_type, kernel_params)

  I = np.eye(N)
  A = lmbda*I+K
  c = y
  
  # Solve A*x=c
  sol = np.linalg.solve(A, c)
  return sol
```


```python
# -----------------------------------------------------------------
# Kernel Ridge Regression
# on an audio sequence corrupted by
# Gaussian noise and outliers
# You need a .wav file to run the experiment...
# Python3 required packages: numpy, soundfile, matplotlib
# -----------------------------------------------------------------
def kernelridge_11_19():
  np.random.seed(0)
  # --------------------------------------------------------------------
  # Reading wav file. x corresponds to time instances (is., x_i in [0,1])
  # fs is the sampling frequency
  # Replace the name "BladeRunner.wav" with the name of the file
  # you intend to use.
  # --------------------------------------------------------------------
  N = 100
  samples = 1000
  indices = range(0, samples,int(samples/N))
  start = 100000
  [data, fs] = sf.read(path)
  sound = np.array(data[start:(start+samples+1), :], dtype=np.float32)
  y = np.reshape(sound[indices, 0], newshape=(len(indices), 1))
  Ts = 1/fs # sampling period
  x = np.array(range(0, samples)).conj().transpose()*Ts # time instances of sampling
  x = x[indices]
  x = np.reshape(x, newshape=(x.shape[0], 1))
  #print(x)
  # -------------------------------------------------------
  # Add white Gaussian noise
  snr = 15 # dB
  y = py_awgn(y, snr)
  # add outliers
  O = 0.8*np.max(np.abs(y))
  percent = 0.1
  M = int(math.floor(percent*N))
  out_ind = np.random.choice(N, M, replace=False)
  outs = np.sign(np.random.randn(M, 1))*O
  y[out_ind] = y[out_ind] + outs
  M = len(y)
  # ----------Code for unbiased L2 Kernel Ridge Regression (KRR-L2)-----------
  C = 0.0001
  kernel_type = 'gaus'
  kernel_params = 0.004
  sol = kernel_regression_l2_l2_unbiased(x.conj().transpose(), y, [C], kernel_type, kernel_params)
  a0 = sol[0:N]
  # Generate regressor
  t = np.array([range(0, samples)]).conj().transpose()*Ts # Here we generate all 1000 points that will be used for prediction
  #t = np.array([range(0, samples+1)]).conj().transpose()*Ts
  M2 = len(t)
  print(M2)
  z0 = np.zeros(shape=(M2, 1))
  for k in range(0, M2):
    z0[k] = 0
    for l in range(0, N):
      z0[k] = z0[k] + a0[l]*kappa(x[l], t[k], kernel_type,[kernel_params])
  # ------------------end unbiased KRR-L2-------------------------------------
  # For unbiased KRR-L2
  plt.figure(1)
  # plot(x,y);
  plt.xlabel('time in sec')
  plt.ylabel('amplitude')
  plt.plot(t, z0, 'r', lw=1)
  plt.plot(x, y, '.', markeredgecolor=[0.3, 0.3, 0.3], markersize=5)
  title = 'unbiased KRR-L2 C= %s' % str(C)
  plt.title(title)
  plt.show()

```


```python
def py_awgn(input_signal, snr_dB, rate=1.0):
  """ Addditive White Gaussian Noise (AWGN) Channel.
  Parameters
  __________
  input_signal : 1D ndarray of floats
  Input signal to the channel.
  snr_dB : float
  Output SNR required in dB.
  rate : float
  Rate of the a FEC code used if any, otherwise 1.
  Returns
  _______
  output_signal : 1D ndarray of floats
  Output signal from the channel with the specified SNR.
  """
  avg_energy = np.sum(np.dot(input_signal.conj().T, input_signal)) / input_signal.shape[0]
  snr_linear = 10 ** (snr_dB / 10.0)
  noise_variance = avg_energy / (2 * rate * snr_linear)
  if input_signal.dtype is np.complex:
    noise = np.array([np.sqrt(noise_variance) * np.random.randn(input_signal.shape[0]) * (1 + 1j)], ndmin=2)
  else:
    noise = np.array([np.sqrt(2 * noise_variance) * np.random.randn(input_signal.shape[0])], ndmin=2)
  output_signal = input_signal + noise.conj().T
  return output_signal

if __name__ == '__main__':
  kernelridge_11_19()
```

Construct the support vector regression model


```python
N = 100
samples = 1000
indices = range(0, samples,int(samples/N))
start = 100000
[data, fs] = sf.read(path)
sound = np.array(data[start:(start+samples+1), :], dtype=np.float32)
y = np.reshape(sound[indices, 0], newshape=(len(indices), 1))
Ts = 1/fs # sampling period
x = np.array(range(0, samples)).conj().transpose()*Ts # time instances of sampling
x = x[indices]
x = np.reshape(x, newshape=(x.shape[0], 1))
#print(x.shape)

#-------------------------------------------------------
# Add white Gaussian noise
snr = 15 # dB
y = py_awgn(y, snr)

# add outliers
O = 0.8*np.max(np.abs(y)) # set the outlier values to 80%
percent = 0.1
M = int(math.floor(percent*N))
out_ind = np.random.choice(N, M, replace=False)
outs = np.sign(np.random.randn(M, 1))*O
y[out_ind] = y[out_ind] + outs
y = y.ravel()
M = len(y)
# print(M)
#print(x.shape)
#print(y.shape)
#-------------------------------------------------------
# support vector regression
#C=1
#epsilon=0.003
#kernel_type='gaus'
#kernel_params=0.004

#-------------------------------------------------------

def svm_regression(sigma, epsilon, C):
  # fit regression model
  # C - float, default=1.0: Regularization parameter. The strength of the regularization is inversely proportional to C. Must be strictly positive. The penalty is a squared l2 penalty.
  # epsilon - float, default=0.1 Epsilon in the epsilon-SVR model. It specifies the epsilon-tube within which no penalty is associated in the training loss function with points predicted within a distance epsilon from the actual value.
  svr_rbf = SVR(kernel='rbf', C=C, gamma = 1/(2*(sigma**2)), epsilon = epsilon)
  svr_lin = SVR(kernel='linear', C=C, gamma = 1/(2*(sigma**2)), epsilon = epsilon)
  svr_poly = SVR(kernel='poly', C=C, gamma = 1/(2*(sigma**2)), epsilon = epsilon, degree=2)
  y_rbf = svr_rbf.fit(x, y).predict(x)
  y_lin = svr_lin.fit(x, y).predict(x)
  y_poly = svr_poly.fit(x, y).predict(x)

  ###############################################################################
  # look at the results
  plt.scatter(x, y, c='k', label='data')
  #plt.hold('on')
  plt.plot(x, y_rbf, c='g', label='RBF kernel')
  plt.plot(x, y_lin, c='r', label='Linear kernel')
  plt.plot(x, y_poly, c='b', label='Polynomial kernel')
  plt.xlabel('time in sec')
  plt.ylabel('amplitude')
  plt.title('Support Vector Regression')
  plt.legend()
  plt.show()
  plt.tight_layout()

```

(a) Find the reconstructed data samples obtained by the support vector regression. Employ the Gaussian kernel with $\sigma$ = 0.004 and set $\epsilon$ = 0.003 and $C$ = 1. Plot the fitted curve of the reconstructed samples together with the data used for training.


```python
print(svm_regression(sigma = 0.004, epsilon = 0.003, C = 1))
```

(b) Repeat step (a) using $C$ = 0.05, 0.1, 0.5, 5, 10, 100.



```python
print(svm_regression(sigma = 0.004, epsilon = 0.003, C = 0.05))
```


```python
print(svm_regression(sigma = 0.004, epsilon = 0.003, C = 0.1))
```


```python
print(svm_regression(sigma = 0.004, epsilon = 0.003, C = 0.5))
```


```python
print(svm_regression(sigma = 0.004, epsilon = 0.003, C = 5))
```


```python
print(svm_regression(sigma = 0.004, epsilon = 0.003, C = 10))
```


```python
print(svm_regression(sigma = 0.004, epsilon = 0.003, C = 100))
```

(c) Repeat step (a) using $\epsilon$ = 0.0005, 0.001, 0.01, 0.05, 0.1.



```python
print(svm_regression(sigma = 0.004, epsilon = 0.0005, C = 1))
```


```python
print(svm_regression(sigma = 0.004, epsilon = 0.001, C = 1))
```


```python
print(svm_regression(sigma = 0.004, epsilon = 0.01, C = 1))
```


```python
print(svm_regression(sigma = 0.004, epsilon = 0.05, C = 1))
```


```python
print(svm_regression(sigma = 0.004, epsilon = 0.1, C = 1))
```

(d) Repeat step (a) using $\sigma$ = 0.001, 0.002, 0.01, 0.05, 0.1.



```python
print(svm_regression(sigma = 0.001, epsilon = 0.003, C = 1))
```


```python
print(svm_regression(sigma = 0.002, epsilon = 0.003, C = 1))
```


```python
print(svm_regression(sigma = 0.01, epsilon = 0.003, C = 1))
```


```python
print(svm_regression(sigma = 0.05, epsilon = 0.003, C = 1))
```


```python
print(svm_regression(sigma = 0.1, epsilon = 0.003, C = 1))
```

(e) Comment on the results.

We experiment with various Gaussian kernal, epsilon, and regularization parameters. After that, we conclude that the polynomial RBF kernel fits the dataset better. When we change the value of C, the model does not change vastly, and the model doesn't perform well in the linear and polynomial kernel. When we increase the value of espilon above 0.01, the predicted lines of SVR models almost touch the 0 in the y-axis, meaning that using a greater epsilon in SVR models with this dataset is not a good idea. Also, if the value of Gaussian kernal above 0.05, the fitted curve would not performs well.

# Exercise 5 (15 points)
Show, using Lagrange multipliers, that the $\ell_2$ minimizer in equation (9.18) from the textbook accepts the closed form solution


$$
\hat{\theta} = X^T(XX^T)^{−1}y
$$

Now, show that for the system $y = X\theta$ with $X \in R^{n×l}$ and $n > l$ the least squares solution is given by:

$$
\hat{\theta} = (X^TX)^{−1}X^Ty
$$


## Answer - $\ell>n$
In regression, we want to select those components of $\theta$ that have the most say in the formation of the output variable. The solution here is given by $\theta$. The $\ell_2$ minimizer tries to minimize

$$
min ||\theta||^{2}_{2}
$$

Such that


$$
x^{T}_{n}\theta = y_n, n = 1, 2, ..., N \to y = x \cdot \theta
$$


where $N<\ell$ (features < observations) is the assumption. Because it penalizes very large coeffients. The corresponding Lagrangian can be calculated as follows:


$$
L(\theta) = \theta^{T}\theta +  \sum_{n=1}^{N}\lambda_n(y_n-X^{T}_n\theta)
$$

Where $\lambda$ is the parameter that controls the weight of the constant.
We take the first order of the above equation and set that to zero. This will give us the minimum:


$$
L'(\theta) = 0 \to \theta = X^T\lambda
$$


where$\lambda = [ \, \lambda_1, \lambda_2, ..., \lambda_N	] \,^T$, and $X$ being the input matrix


$$
X = \begin{bmatrix}
    x_1^T \\
    \vdots  \\
    x_N^T 
\end{bmatrix}
$$


We plug the solution into the set of constraints, which can be written
as $ y = X\theta$, now we get:


$$
X\theta = y = XX^T\lambda \to \lambda = (XX^T)^{-1} y
$$


So in ridge regression we have a closed form solution. Therefore, the final result is


$$
\theta = X^T(XX^T)^{-1}y
$$

## Answer - $n>\ell$

We know that
$$
y = X \cdot \theta
$$
based on the above, we can get, 
$$
||y - x \cdot \theta ||^2_2 = (y - x \cdot \theta)^T \cdot (y - x \cdot \theta)^T
$$
then,
$$
= y^T\cdot y-y^T\cdot x\cdot\theta-\theta^T\cdot x^T\cdot y+\theta^T\cdot x^T\cdot x\cdot\theta
$$
then,
$$
(-y^T\cdot x)^T - (x^T\cdot y) + 2 \cdot x^T \cdot x \cdot \theta = 0
$$
So we can get, 
$$
\theta = X^T(XX^T)^{-1}y
$$

# Exercise 6 (10 points)
Show that the null space of a full rank $N × l$ matrix $X$ is a subspace of imensionality $l − N$ , for $N < l$.

## Answer
Because $N \cdot \ell$ is a full rank matrix where $N < \ell$. We also know `rank(x) = dim(rank(x)) = N`. So we have `dim(null(x)) + dim(rank(x)) = l`. Therefore, `dim(null(x)) = l-n`, where `dim` is a dimension of a matrix `x`.

# Exercise 7
Generate in Python a sparse vector $\theta \in R^l, l = 100$, with its first five components taking random values drawn from a normal distribution with mean zero and variance one and the rest being equal to zero. Build, also, a sensing matrix $X$ with $N = 30$ rows having samples normally distributed, with mean zero and variance $1/N$, in order to get 30 observations based on the linear regression model $y = X\theta$. Then perform the following tasks.


```python
import numpy as np
import matplotlib.pyplot as plt
import math
import random

np.random.seed(1234)
N = 30 # 30 rows
k = 5 # first five components taking random values drawn from a normal distribution with mean zero and variance one and the rest being equal to zero
l = 100 # 100 columns

theta = np.zeros((l, 1))
theta[0:k] = np.random.standard_normal(size = (k, 1)) # draw samples from a standard Normal distribution (mean=0, stdev=1)
# print(theta)
X = np.dot(np.random.standard_normal(size = (N, l)),(1/math.sqrt(N)))
y = np.dot(X, theta)
print("X size is {}".format(X.shape))
print("y size is {}".format(y.shape))
#print(theta)
```

    X size is (30, 100)
    y size is (30, 1)


(a) Use a LASSO implementation to reconstruct $\theta$ from $y$ and $X$.

Lasso regression performs L1 regularization, which adds a penalty equal to the absolute value of the magnitude of coefficients. This type of regularization can result in sparse models with few coefficients; Some coefficients can become zero and eliminated from the model. Larger penalties result in coefficient values closer to zero, which is the ideal for producing simpler models.

L1 tends to shink coefficients to zero. As we can see below, the cofficients are almost zero.


```python
# https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Lasso.html
from sklearn.linear_model import Lasso
from sklearn.linear_model import LinearRegression

lasso_reconst = Lasso(fit_intercept=True, normalize=True, max_iter=1e6)
lasso_reconst.fit(X, y)
error_lasso = np.linalg.norm(np.reshape(np.array(lasso_reconst.coef_), newshape=theta.shape)-theta)
print(lasso_reconst.coef_, lasso_reconst.intercept_)
print(error_lasso)
```

    [ 0. -0.  0. -0. -0. -0. -0.  0.  0. -0. -0. -0.  0.  0.  0.  0. -0. -0.
     -0. -0. -0.  0.  0.  0. -0. -0.  0.  0. -0.  0.  0.  0. -0. -0. -0.  0.
      0.  0.  0.  0. -0.  0.  0.  0.  0. -0.  0.  0. -0.  0.  0. -0.  0. -0.
      0.  0. -0. -0.  0. -0.  0.  0. -0. -0. -0.  0.  0.  0. -0. -0. -0.  0.
     -0.  0.  0.  0. -0.  0. -0.  0.  0.  0.  0. -0. -0.  0.  0.  0. -0. -0.
      0.  0.  0. -0. -0. -0. -0.  0. -0. -0.] [0.03136091]
    2.0761316931119906


When we compare the coefficients in the linear model, they are non-zero.


```python
linear_reconst = LinearRegression()
linear_reconst.fit(X, y)
#error_lin = np.linalg.norm(theta - linear_reconst.coef_)
print(linear_reconst.coef_)
#print(error_lin)
```

    [[ 0.1504216  -0.24852887  0.2851724  -0.08915576 -0.10327576 -0.12618162
      -0.04602372  0.05473597  0.07021556 -0.05015795 -0.0767972  -0.00251141
       0.08548682  0.19020808  0.06857647  0.03895134 -0.09698147 -0.10210427
      -0.22613006 -0.04752304 -0.0878581   0.07235256  0.04932862 -0.00278751
      -0.11543169 -0.00846009  0.12466218 -0.01432902 -0.04694153  0.09194229
       0.02842385  0.05432163 -0.02859287 -0.10724279 -0.06568882  0.12864594
       0.01044759  0.16077251  0.10993205  0.12337903 -0.04621973  0.08477134
       0.04870411  0.11392832  0.11209879 -0.02476117 -0.02252337 -0.00092738
      -0.11909358  0.07232582 -0.01891333 -0.14706331 -0.03034476 -0.08839947
       0.0574367   0.08681011 -0.02308903 -0.02174072  0.13930187 -0.0745217
       0.02068142  0.0311157  -0.03109188 -0.16113841 -0.01404451  0.07568557
       0.04640705  0.06853807 -0.09811168 -0.19712532 -0.00386257  0.07060244
      -0.11120965 -0.01141535 -0.02502231  0.12010114 -0.05180121  0.2208544
      -0.10120215  0.02234681  0.04724423 -0.07606814  0.07811075 -0.00621995
      -0.01792998 -0.00034544  0.06737222  0.01680416 -0.09043721 -0.03464077
       0.02525259 -0.00389582  0.01578966 -0.03366412 -0.01596172 -0.105041
      -0.19233346  0.01089446  0.05992144  0.02080798]]


(b) Repeat the experiment 500 times, with different realizations of $X$, in order to compute the probability of correct reconstruction (assume the reconstruction is exact when $||y =X\theta|| < 10^{−8}$).


The probability is 0 if we computer $||y = X\theta|| < 10^{-8}$ with the experiment 500 times.



```python
rep = 500 # repeat the experiment 500 times
error = np.zeros((rep, 1))
print(error.shape)

for i in range(0, rep):
  X = np.random.randn(N,l)*(1/math.sqrt(N))
  #X = np.dot(np.random.standard_normal(size = (N, l)),(1/math.sqrt(N))) # random values
  y = np.dot(X, theta)
  model = Lasso(alpha=0.0001, fit_intercept=False, normalize=False, max_iter=1e6).fit(X, y)
  values = model.coef_
  errorX = np.linalg.norm(np.reshape(np.array(values), newshape = theta.shape)-theta)
  #print(errorX)
  if errorX < 10**(-8):
    error[i] = 1
  else:
    error[i] = 0

probRandn = np.sum(error)/rep
print('Random Sensing Mtx: '+ str(probRandn))
```

    (500, 1)
    Random Sensing Mtx: 0.0


(c) Repeat the same experiment (500 times) with matrices of the form


$$
\begin{equation}
            \nonumber X(i, j) = \left\{
            \begin{array}{l l}
              +  \sqrt{\frac{\sqrt{p}}{N}}, & \quad\text{with probability}  \frac{1}{2\sqrt{p}}\\
              0, & \quad\text{with probability}  1- \frac{1}{\sqrt{p}}\\
              - \sqrt{\frac{\sqrt{p}}{N}}, &\quad\text{with probability}  \frac{1}{2\sqrt{p}}
            \end{array} \right.
            \end{equation}
$$


for p equal to 1, 9, 25, 36, 64 (make sure that each row and each column of $X$ has at least a nonzero component). Give an explanation why the probability of reconstruction falls as $p$ increases (observe that both the sensing matrix and the unknown vector are sparse).

## Explanation
The range of error terms are going to wider if p increased. With the larger error value, the estimated results are more likely to be bias, so the probability of reconstruction falls as p increases.



```python
p = [1, 9, 25, 36, 64]
Error = np.zeros((rep, 1))
eachMinErr = 9999
eachMaxErr = 0

# construct sparse sensing matrices
for pval in p:
  for i in range(rep):
    OK = False
    while not OK:
      kk = np.zeros((N*l, 1))
      P = np.random.permutation(N*l)
      #numofzeros = round(N*l*(1-1/math.sqrt(pval)))
      numofzeros = int(np.around(N*l*(1-1/math.sqrt(pval))))
      P[(N*l-numofzeros+1):] = 0
      #kk[P] = np.multiply(np.matmul(math.sqrt((math.sqrt(pval)/N)), np.ones((len(P),1))), np.sign(np.random.randn(len(P),1)))
      kk[P] = np.sqrt((math.sqrt(pval)/N))*np.multiply(np.ones(shape=(len(P), 1)), np.sign(np.random.randn(len(P), 1)))
      X = np.zeros((N, l))
      X[:(l*N)] = np.reshape(kk, newshape = (N, l))
      #OK = np.linalg.matrix_rank(X) == N # Return matrix rank of array using SVD method, also check if there is a full rank
      # Check if it is full rank
      if np.linalg.matrix_rank(X) == N:
        OK = True
      else:
        OK = False
    y = np.dot(X, theta)
    L1_model = Lasso(alpha=0.0001, fit_intercept=False, normalize=False, max_iter=1e6).fit(X, y)
    values = L1_model.coef_
    #errorX = np.linalg.norm(L1_model.coef_-theta)
    errorX = np.linalg.norm(np.reshape(np.array(values),newshape=theta.shape) - theta)
    #print(errorX)
    if errorX < 10 ** (-8):
      Error[i] = 1
    else:
      Error[i] = 0
    eachMinErr = min(eachMinErr, errorX)
    eachMaxErr = max(eachMaxErr, errorX)
  print("min err: {}, max err:{}".format(eachMinErr, eachMaxErr))
  probSparse = np.sum(Error)/rep
  #print(pval)
  print('Sparse Sensing Mtx, p = '+str(pval)+', the probability is '+str(probSparse))
```

    min err: 0.0039041880894823524, max err:1.0101686575393911
    Sparse Sensing Mtx, p = 1, the probability is 0.0
    min err: 0.0039041880894823524, max err:1.5537180405710938
    Sparse Sensing Mtx, p = 9, the probability is 0.0
    min err: 0.003875605385666494, max err:1.5537180405710938
    Sparse Sensing Mtx, p = 25, the probability is 0.0
    min err: 0.0038363421335191514, max err:1.6001763493778818
    Sparse Sensing Mtx, p = 36, the probability is 0.0
    min err: 0.0034991983600240503, max err:1.7485271633545891
    Sparse Sensing Mtx, p = 64, the probability is 0.0


# Testing Zone


```python
"""
kk = np.zeros((N*l, 1))
X = np.zeros((N, l))
print(kk.shape)
print(X.shape)
print(X[:l*N].shape)
a = np.multiply(np.dot(math.sqrt((math.sqrt(9)/N)), np.ones((len(P),1))) , np.sign(np.random.rand(len(P),1)))
print(a.shape)
"""
```




    '\nkk = np.zeros((N*l, 1))\nX = np.zeros((N, l))\nprint(kk.shape)\nprint(X.shape)\nprint(X[:l*N].shape)\na = np.multiply(np.dot(math.sqrt((math.sqrt(9)/N)), np.ones((len(P),1))) , np.sign(np.random.rand(len(P),1)))\nprint(a.shape)\n'




```python
"""
P = np.random.permutation(N*l)
print(P.shape)
"""
```




    '\nP = np.random.permutation(N*l)\nprint(P.shape)\n'



# Output


```python
# should access the Google Drive files before running the chunk
%%capture
!sudo apt-get install texlive-xetex texlive-fonts-recommended texlive-plain-generic 
!jupyter nbconvert --to html "/content/drive/MyDrive/American_University/2021_Fall/DATA-642-001_Advanced Machine Learning/GitHub/Midterm/submit/midterm_Yunting.ipynb"
```

# References:
- https://www.statisticshowto.com/lasso-regression/ 
- https://scikit-learn.org/stable/modules/generated/sklearn.svm.SVR.html
- https://ai.stanford.edu/~gwthomas/notes/convexity.pdf
- https://en.wikipedia.org/wiki/Convex_function
- https://math.stackexchange.com/questions/1356842/is-l-2-norm-a-strictly-convex-function
- https://cims.nyu.edu/~cfgranda/pages/OBDA_spring16/material/convex_optimization.pdf
- https://math.stackexchange.com/questions/80139/why-is-the-l-p-norm-strictly-convex-for-1p-infty
- https://ai.stanford.edu/~gwthomas/notes/convexity.pdf
- https://scikit-learn.org/stable/modules/generated/sklearn.svm.SVR.html


a6d89cf3afafec9f046e8a3a38974f140f4857dd
ipynb
Jupyter Notebook
Midterm/submit/midterm_Yunting.ipynb
twyunting/DATA-642_Advanced_Machine_Learning
fffc2621f0dfe86c315a75b6331da9b54a11dd1a
MIT
Midterm/submit/midterm_Yunting.ipynb
twyunting/DATA-642_Advanced_Machine_Learning
fffc2621f0dfe86c315a75b6331da9b54a11dd1a
MIT
Midterm/submit/midterm_Yunting.ipynb
twyunting/DATA-642_Advanced_Machine_Learning
fffc2621f0dfe86c315a75b6331da9b54a11dd1a
MIT
Qwen/Qwen-72B
1. YES
2. YES