<jupyter_start><jupyter_text>Unsupervised Learning Methods   Problem Set 1   Optimization and Clustering<jupyter_code>import numpy             as np
import matplotlib.pyplot as plt
import matplotlib

matplotlib.rc('font', **{'size' : 16})
A = set()
A.add('3')
A.update(['3'])
A.update(['4'])
A<jupyter_output><empty_output><jupyter_text>1. Optimization  1.3 Descent Methods  Gradient descent* Let $\boldsymbol{Y}\in\mathbb{R}^{d\times d}$ be a blurred versionof the image $\boldsymbol{X}\in\mathbb{R}^{d\times d}$:$$\boldsymbol{Y}=\boldsymbol{H}\boldsymbol{X}\boldsymbol{H}^{T}$$where $\boldsymbol{H}\in\mathbb{R}^{d\times d}$ is a known (separable)blurring matrix.  Load the data ($\boldsymbol{Y}$ and $\boldsymbol{H}$):<jupyter_code>from scipy.io import loadmat

data = loadmat('data.mat')
Y    = data['Y']
H    = data['H']

Y.shape, H.shape<jupyter_output><empty_output><jupyter_text>Plot $Y$ and $H$:<jupyter_code>_, ax = plt.subplots(1, 2, figsize=(12, 6))
ax[0].imshow(Y, cmap='gray')
ax[1].matshow(H)
ax[0].set_title('$Y = H X H^T$')
ax[1].set_title('$H$')
plt.show()<jupyter_output><empty_output><jupyter_text>* Our goal is to find $\boldsymbol{X}$ that minimizes:$$\boldsymbol{X}^{\star}=\arg\min_{\boldsymbol{X}}f\left(\boldsymbol{X}\right)=\arg\min_{\boldsymbol{X}}\left\Vert \boldsymbol{Y}-\boldsymbol{H}\boldsymbol{X}\boldsymbol{H}^{T}\right\Vert _{F}^{2}$$  1.3.11. Find a closed form expression for $\boldsymbol{X}^{\star}$ as a function of $\boldsymbol{Y}$ and $\boldsymbol{H}$.2. Compute $\boldsymbol{X}^{\star}$ and plot it.<jupyter_code>'''
Your code here...
'''

_, ax = plt.subplots(1, 2, figsize=(12, 6))
ax[0].imshow(Y,     cmap='gray')
ax[1].imshow(Xstar, cmap='gray')
ax[0].set_title('$Y = H X H^T$')
ax[1].set_title('$X^{\star}$')
plt.show()<jupyter_output><empty_output><jupyter_text>------------------------------------------------------------------------------------------------------* To avoid matrix inversion, we will find $\boldsymbol{X}^{\star}$ using gradient descent. 1.3.21. Implement the following functions where$$f\left(\boldsymbol{X}\right)=\left\Vert \boldsymbol{Y}-\boldsymbol{H}\boldsymbol{X}\boldsymbol{H}^{T}\right\Vert _{F}^{2}$$```python==================================================def f(X):    '''    f(X) = ||Y - HXH^T||_F^2 is the objective function    '''    pass====================================================================================================def Df(X):    '''    Df(X) = ∇f(X) is the objective gradient    '''    pass====================================================================================================def ApplyGradientDescent(f, Df, X0, lr, ε):    '''    Apply gradient descent.    Args:        f  - Objective function        Df - Objective gradient        X0 - Initial point        lr - Learning rate (step size)        ε  - Stopping criterion value    Output:        X  - The converged X    Stopping criterion:        Stop iterate when f(X) < ε    '''    pass==================================================```2. Set a reasonable starting point $\boldsymbol{X}^{\left(0\right)}$.3. Find $\boldsymbol{X}$ such that $f\left(\boldsymbol{X}\right)<\epsilon$ where $\epsilon=10^{-4}$.    * Time your function.    * Plot the final obtained $\boldsymbol{X}$.<jupyter_code>'''
Your code here...
'''<jupyter_output><empty_output><jupyter_text>------------------------------------------------------------------------------------------------------ Momentum* Momentum is an effective method to improve convergence rate.* Gradient descent update rule:$$\boldsymbol{x}^{\left(t\right)}=\boldsymbol{x}^{\left(t-1\right)}-\mu\nabla f\left(\boldsymbol{x}^{\left(t-1\right)}\right)$$* Momentum update rule:$$\boldsymbol{x}^{\left(t\right)}=\boldsymbol{x}^{\left(t-1\right)}-\mu\nabla f\left(\boldsymbol{x}^{\left(t-1\right)}\right)+\beta\left(\boldsymbol{x}^{\left(t-1\right)}-\boldsymbol{x}^{\left(t-2\right)}\right)$$where $\beta\in[0,1)$ and typically $0.9\leq\beta\leq0.99$. 1.3.3Show that $\boldsymbol{x}^{\left(t\right)}$ can be expressed by:$$\boldsymbol{x}^{\left(t+1\right)}=\boldsymbol{x}^{\left(t\right)}-\mu\sum_{s=0}^{t}\beta^{t-s}\nabla f\left(\boldsymbol{x}^{\left(s\right)}\right)$$ ------------------------------------------------------------------------------------------------------ 1.3.41. Implement the following function:```pythondef ApplyGradientDescentMomentum(f, Df, X0, lr, β, ε):    '''    Apply gradient descent with momentum.    Args:        f  - Objective function        Df - Objective gradient        X0 - Initial point        lr - Learning rate (step size)        β  - Momentum decaying factor        ε  - Stopping criterion value    Output:        X  - The converged X    Stopping criterion:        Stop iterate when f(X) < ε    '''    pass```2. Set (the same) reasonable starting point $\boldsymbol{X}^{\left(0\right)}$.3. Find $\boldsymbol{X}$ such that $f\left(\boldsymbol{X}\right)<\epsilon$ where $\epsilon=10^{-4}$.    * Time your function.    * Plot the final obtained $\boldsymbol{X}$.     * Compare the run-time performance with and without momentum.        (play with the hyper-parameters).        * **2% Bonus:**              Plot the objective value $f(\boldsymbol{X})$ as a function of the iterations.              You may need to slightly modify your code.<jupyter_code>'''
Your code here...
'''<jupyter_output><empty_output><jupyter_text>------------------------------------------------------------------------------------------------------ 2. Clustering 2.1 K-Means 2.1.3Implement the K-Means algorithm.```pythondef KMeans(mX, K):    '''    Apply K-Means    Args:        mX   - Input data,         mX.shape = d, N        K    - Number of clusters, positive integer    Output:        vIdx - Index of the cluster each point belongs to, vIdx.shape = N        mMu  - Clusters' centers,                          mMu.shape  = d, K            '''    pass```**Notes**:* The algorithm halts when no update occur.  * **4% Bonus**: Implement the K-Menas++ initialization: https://en.wikipedia.org/wiki/K-means%2B%2B------------------------------------------------------------------------------------------------------<jupyter_code>'''
Your code here...
'''<jupyter_output><empty_output><jupyter_text>Super-pixels1. Load the `Fruits.jpeg` image and covert it to NumPy ndarray `mI`:$$\boldsymbol{I}\in\mathbb{R}^{\text{height}\times\text{width}\times3}$$2. Use the LAB color space (instead of RGB) (`mLAB`).3. Create a new "image" `mX` such that:$$\boldsymbol{X}\in\mathbb{R}^{\text{height}\times\text{width}\times5}$$where:    * The first 3 channels are the LAB image.    * The 4th channel is the $x$ position.    * The 5th channel is the $y$ position.4. Apply K-Means to the pixels of `mX` (set a reasonable $K$).5. Create a mask image `mMask` such that:$$\boldsymbol{M}\in\mathbb{R}^{\text{height}\times\text{width}}$$where each pixel in `mMask` is the cluster index of the corresponding pixel in `mI`.6. Plot the segmentation (Superpixels) map.**Tip:** Try different weights for the LAB image and the XY position.<jupyter_code>from PIL                  import Image
from skimage              import color
from skimage.segmentation import mark_boundaries

oImage = Image.open('Fruits.jpeg')
oImage
mI   = np.array(oImage)
mLAB = color.rgb2lab(mI)

'''
Your code here...
'''

mO = mark_boundaries(mI, mMask, color=(0,1,1))

plt.figure(figsize=(12, 12))
plt.imshow(mO)
plt.show()<jupyter_output><empty_output><jupyter_text>------------------------------------------------------------------------------------------------------ 2.2 GMMThe GMM algorithm aims to maximize the (log) likelihood objective:$$\arg\max_{\left\{ \left(w_{k},\boldsymbol{\mu}_{k},\boldsymbol{\Sigma}_{k}\right)\right\} _{k=1}^{K}}f=\arg\max_{\left\{ \left(w_{k},\boldsymbol{\mu}_{k},\boldsymbol{\Sigma}_{k}\right)\right\} _{k=1}^{K}}\sum_{i=1}^{N}\log\left(\sum_{k=1}^{K}w_{k}\mathcal{N}\left(\boldsymbol{x}_{i};\boldsymbol{\mu}_{k},\boldsymbol{\Sigma}_{k}\right)\right)$$ 2.2.3Implement the GMM algorithm.```pythondef GMM(mX, K, ε=1e-6, maxIter=100):    '''    Apply GMM    Args:        mX      - Input data,                   mX.shape = d, N        K       - Number of clusters,           positive integer        ε       - Stopping criterion threshold, positive real        maxIter - Maximum number of iterations, positive integer    Output:        vIdx - Index of the cluster each point belongs to, vIdx.shape = N        vW   - The weight of each Gaussian,                vW.shape   = K        mMu  - Centers,                                    mMu.shape  = d, K        mSig - Covariance,                                 mSig.shape = d, d, K    '''    pass```**Notes**:* Stopping criteria:    * The objective $f$ has changed less than $\epsilon$.    * Maximum number of iteration.* The index $s$ of the point $\boldsymbol{x}_i$ is defined by:$$\boldsymbol{x}_{i}\in\mathcal{D}_{s} \iff s=\arg\max_{k}w_{k}\mathcal{N}\left(\boldsymbol{x}_{i},\boldsymbol{\mu}_{k},\boldsymbol{\Sigma}_{k}\right)$$* **2% Bonus:** Use [broadcasting](https://numpy.org/doc/stable/user/basics.broadcasting.html) and avoid `for` loops in step II.<jupyter_code>'''
Your code here...
'''<jupyter_output><empty_output><jupyter_text>------------------------------------------------------------------------------------------------------Consider the following data:<jupyter_code>import scipy.stats
import scipy as sp


np.random.seed(0)

N1    = 250
N2    = 150
N3    = 200

vMu1  = np.array([0, 0])
vMu2  = np.array([2, 0.5])
vMu3  = np.array([4, 1])

mSig1 = .5 * np.array([[1.00, 1.25],
                       [1.25, 2.00]])

mSig2 = .5 * np.array([[ 1.00, -1.25],
                       [-1.25,  2.00]])

mSig3 = .5 * np.array([[1.00, 1.25],
                       [1.25, 2.00]])

mX1 = sp.stats.multivariate_normal.rvs(mean=vMu1, cov=mSig1, size=N1)
mX2 = sp.stats.multivariate_normal.rvs(mean=vMu2, cov=mSig2, size=N2)
mX3 = sp.stats.multivariate_normal.rvs(mean=vMu3, cov=mSig3, size=N3)
mX  = np.r_[mX1, mX2, mX3].T

plt.figure(figsize=(6,5))
plt.scatter(mX[0,:], mX[1,:], s=50, edgecolors='k', color='b')
plt.axis('equal')
plt.axis('tight')
plt.show()<jupyter_output><empty_output><jupyter_text>2.2.4 Compare K-Means to GMM * Set $K=3$ and apply K-Means and GMM to the above dataset.* Repeat several times.<jupyter_code>'''
Your code here...
'''<jupyter_output><empty_output><jupyter_text>------------------------------------------------------------------------------------------------------ 2.4 DBSCAN 2.4.1Implement the DBSCAN algorithm.```pythondef DBSCAN(mX, Z, r):    '''    Apply DBSCAN    Args:        mX   - Input data,                                   mX.shape = d, N        Z    - Number of points required to be a core point, positive integer        r    - Neighborhood radius,                          positive real.    Output:        vIdx - Index of the cluster each point belongs to,   vIdx.shape = N    '''    pass```**Notes**:* Noise points should have index `-1`.* Implement an auxiliary function to compute connected components (using BFS or DFS).<jupyter_code>'''
Your code here...
'''<jupyter_output><empty_output><jupyter_text>2.4.2 Comparison1. Consider the following datasets:<jupyter_code>from sklearn import datasets

np.random.seed(0)

N   = 1000
A   = np.array([[0.6, -0.6],
                [-0.4, 0.8]]);
mX1 = datasets.make_circles(n_samples=N, noise=0.02)      [0]
mX2 = datasets.make_moons  (n_samples=N, noise=0.05)      [0]
mX3 = datasets.make_blobs  (n_samples=N, random_state=170)[0] @ A
mX4 = datasets.make_blobs  (n_samples=N, random_state=170, cluster_std=[.8, 2, .4])[0] 
mX5 = np.load('clusterable_data.npy')

lDatasets = [mX1, mX2, mX3, mX4, mX5]

fig, ax = plt.subplots(1, 5, figsize=(18, 4))
for ii, ax in enumerate(fig.axes):
    mXi = lDatasets[ii]
    ax.scatter(*mXi.T, c='lime', s=15, edgecolor='k')
    ax.axis('equal')
    
plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>2. Apply and compare your clustering methods to these datasets.    * **2% Bonus:**      Add to your comparison the hierarchical and HDBSCAN methods.      You may use existing libraries for these methods.    **Note:** The datasets are transposed.<jupyter_code>'''
Your code here...
'''<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.2493052149746608, 'fraction_numerical': 0.022233120810854993, 'mean_word_length': 2.7222391238211134, 'pattern_counts': {'":': 0, '<': 61, '<?xml version=': 0, '>': 57, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '589819', 'license': 'MIT', 'n_tokens_mistral': 4027, 'n_tokens_neox': 3788, 'n_words': 950, 'path': 'ProblemSet1/ProblemSet1_Original.ipynb', 'repo_name': 'alonbar/unsupervised-learning', 'size': '12234'}