# 3. Hands-on experiment:  Computing a geometric embedding for the Sparsest Cut Problem via Semidefinite Programming (20 pts)


The Uniform Sparsest Cut problem (USC) aims to find a bipartition $(S, \bar{S})$ of the nodes of a graph $G = (V, E)$, $\vert V \vert = p$, which minimizes the quantity 
\begin{equation*}
\label{sceq}
       \frac{E(S, \bar{S})}{\lvert S \rvert \; \lvert \bar{S} \rvert}, 
\end{equation*}
where $E(S, \bar{S})$ is the number of edges connecting $S$ and $\bar{S}$, and $\lvert S \rvert $ is the number of nodes in $S$.

Computing such a bipartition is NP-hard and intense research has gone into designing efficient approximation algorithms for this problem. In the seminal work of Arora, Rao and Vazirani (2009), an $\mathcal{O}(\sqrt{\log p})$ approximation algorithm is proposed for solving USC, which relies on finding a well-spread $\ell_2^2$ geometric representation of $G$ where each node $i\in V$ is mapped to a vector $\mathbf{v}_i$ in $\mathbb{R}^p$. In this experimental section we focus on solving the SDP that computes this geometric embedding. 

### SDP formulation
The canonical formulation of the SDP is
\begin{align}
\mathbf{X}^\star \in \arg \min_{ \mathbf{X} } \bigg\{ \langle \boldsymbol{C}, \mathbf{X} \rangle : &~p \; \mathrm{Tr}(\mathbf{X}) -  \mathrm{Tr}(\mathbf{1}_{p\times p}\mathbf{X}) = \frac{p^2}{2}, \\\\ 
& \; \mathbf{X}_{i,j} + \mathbf{X}_{j,k} - \mathbf{X}_{i, k} - \mathbf{X}_{j,j} \leq 0, \; \forall\ i \neq j \neq k \neq i\in V, \\\\
&~\underbrace{\mathrm{Tr}(\mathbf{X}) \leq p, ~\mathbf{X} \in \mathbb{R}^{p\times p} ,~\mathbf{X}\succeq 0}_{\mathcal{X}} \bigg\},
\end{align}
where $\boldsymbol{C}$ represents the Laplacian of graph $G$ and $\mathbf{X}_{i,j} = \langle \mathbf{v}_{i}, \mathbf{v}_j\rangle$ gives the geometric embedding of the nodes. We can rewrite the above SDP as


\begin{equation}
    \label{eq:reform-sdp}
    \min_{\mathbf{X} \in \mathcal{X}} f(\mathbf{X}) + g(A(\mathbf{X})) ~~~~~~~~\text{ subject to } ~~~~B_{i,j,k}(\mathbf{X}) \in \mathcal{K}, \;\forall\ i \neq j \neq k \neq i\in V,
\end{equation}
where $f(\mathbf{X}) = \langle \boldsymbol{C}, \mathbf{X} \rangle$ and  $g(\cdot) = \delta_{\left\{\frac{p^2}{2}\right\}}(\cdot)$ is the indicator function of singleton $\left\{\frac{p^2}{2}\right\}$.



```python
import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse.linalg import eigsh, svds, eigs
from math import sqrt
import scipy.io as sio
import random
import numpy.matlib
from helpers import *
import time
from scipy.sparse import isspmatrix
import itertools
from datetime import datetime

# fix the seed
random.seed( 3 )
```


```python
# Load data
dataset = 'data/102n-insecta-ant-colony4-day10.mat'
# dataset = 'data/55n-insecta-ant-colony1-day37.mat'
# dataset = 'data/25mammalia-primate-association-13.mat'

Problem = sio.loadmat(dataset)
C = np.double(Problem['LAPL']) # euclidean distance matrix
p = np.int(Problem['d']) # number of data points
opt_val = Problem['opt_val'] # optimum value 
X_true = Problem['X_true']
nchoosek_inds = list(itertools.combinations(range(p), 3))
```


#### <span style="font-variant:small-caps;">3.(c) *(10 points)*</span>

We will now observe the behavior of HCGM on three graphs from the Network Repository dataset \cite{nr}: 
* **G1 (25 nodes)**: *mammalia-primate-association-13*, 
* **G2 (55 nodes)**: *55n-insecta-ant-colony1-day37*
* **G3 (102 nodes)**: *insecta-ant-colony4-day10*

$\triangleright$ Fill in the missing code and run the algorithm.

$\triangleright$ Add resulting plots to your report.

$\triangleright$ What do you notice about the running times of the algorithm for the three problem instances?

$\triangleright$ What are the potential bottlenecks to applying this method to large graphs?



```python
def grad_F(X, beta):
    # Add the equality constraint
    [grad_val, displacement] = get_equality_constr_grad(X, beta)
    feas_eq = abs(displacement)
   
    # Add the triangle constraints. This needs to be scaled by d to match
    # the order of the equality constr
    n = len(nchoosek_inds)
    feas_ineq = 0
    for l in range(n):
        index_triplet = nchoosek_inds[l]
        [update_vector, row_idxs, col_idxs, part_feas_normsq] = get_triangle_constr_grad(X, beta, index_triplet)
        
        grad_val[row_idxs, col_idxs] = grad_val[row_idxs, col_idxs] + update_vector
        grad_val[col_idxs, row_idxs] = grad_val[col_idxs, row_idxs] + update_vector
        
        feas_ineq = feas_ineq + part_feas_normsq
    
    # Total gradient
    grad_val = grad_val + C
    feas_ineq = sqrt(feas_ineq)
    
    return (grad_val, feas_eq, feas_ineq)



def get_equality_constr_grad(X, beta):
    # X is a p x p matrix
    displacement = p * np.trace(X) - np.sum(X) - p**2/2
    grad_val = (p/beta) * (p * np.eye(p) - np.ones((p, p))) * displacement
    
    return (grad_val, displacement)



def get_triangle_constr_grad(X, beta, index_triplet):
    # X needs to be in a matrix shape, not vectorized
    #print(index_triplet)
    i = index_triplet[0]
    j = index_triplet[1]
    k = index_triplet[2]
    t1 = max(X[i, j] + X[j, k] - X[i, k] - X[j, j], 0.)
    t2 = max(X[i, k] + X[i, j] - X[j, k] - X[i, i], 0.)
    t3 = max(X[i, k] + X[j, k] - X[i, j] - X[k, k], 0.)

    grad_update = (1. / beta) * np.array([t1 + t2 - t3, # for Xij 
                                -t1 + t2 + t3, # for Xik
                                t1 - t2 + t3, # for Xjk
                                -t2/2., # for Xii - ./2 because we update these twice when we symmetrize
                                -t1/2., # for Xjj - ./2 because we update these twice when we symmetrize
                                -t3/2.]) # for Xkk - ./2 because we update these twice when we symmetrize

    update_row_idxs = [i, i, j, i, j, k]
    update_col_idxs = [j, k, k, i, j, k]
    part_feas_normsq = 2. * (t1**2 + t2**2 + t3**2)
    
    return (grad_update, update_row_idxs, update_col_idxs, part_feas_normsq)
```


```python
# Plotting function
def plot_func(cur_iter, feasibility1,feasibility2, objective, X, X_true):
    date = datetime.now().strftime("%y%m%d_%H%M%S")

    plt.figure(figsize=(12, 8))
    plt.subplot(221)
    plt.loglog(cur_iter, feasibility1)#, 'go--', linewidth=2, markersize=12))
    plt.xlabel('iteration',fontsize=15)
    plt.ylabel('$\|A(X) - b\|$',fontsize=15)
    plt.grid(True)

    plt.subplot(222)
    plt.loglog(cur_iter, feasibility2)
    plt.xlabel('iteration',fontsize=15)
    plt.ylabel('dist$(B(X), (-\infty, 0])$',fontsize=15)
    plt.grid(True)
    plt.savefig(f'./results1_{date}.pdf', bbox_inches='tight')
    plt.show()

    #plt.subplot(223)
    obj_res = np.reshape(np.abs(objective - opt_val)/opt_val, (len(objective),))
    plt.figure(figsize=(12, 8))
    plt.loglog((cur_iter), (obj_res))
    plt.xlabel('iteration',fontsize=15)
    plt.ylabel('$(f(X) - f^*)/f^*$',fontsize=15)
    plt.title('Relative objective residual',fontsize=15)
    plt.grid(True)
    plt.savefig(f'./results2_{date}.pdf', bbox_inches='tight')
    plt.show()
    
    plt.figure(figsize=(12, 8))
    plt.subplot(221)
    plt.imshow(X)
    plt.title('HCGM solution',fontsize=15)
    plt.colorbar()
    
    plt.subplot(222)
    plt.imshow(X_true)
    plt.title('True solution',fontsize=15)
    plt.colorbar()
    plt.savefig(f'./results3_{date}.pdf', bbox_inches='tight')
    plt.show()
```


```python
def HCGM(kappa=p, maxit=np.int(1e3), beta0=1e2):
    # Initialize
    X = np.zeros((p,p))
    
    feasibility1 = [] # norm(A1(X)-b1)/norm(b1)
    feasibility2 = [] # dist(X, \mathcal{K})
    objective    = [] # f(x)
    cur_iter    = [] 
    t    = [] 
    
    #u = np.zeros((N,1))
    iter_track = np.unique(np.ceil(np.power(2, np.linspace(0,20,50))))
    
    start = time.time()
    
    for iteration in range(1, maxit+1):
        
        # Update Step Size
        gamma = 2 / (iteration + 1)
        
        # Update beta
        beta_ = beta0 / np.sqrt(iteration + 1)
        
        # Write down the vk to use in the lmo (eigenvalue routine)
        (vk, feas_eq, feas_ineq) = grad_F(X, beta_)
        # the above needs to be symmetric
        
        
        # Linear minimization oracle
        q, u = eigsh(vk, k=1, tol=1e-16, which='SA')
        if q >= 0:
            X_bar = 0
        else:
            u = sqrt(kappa)*u
            X_bar = np.outer(u,u)
        
        # Update X
        X = (1.0 - gamma) * X + gamma * (X_bar)
                
        if any(iteration == iter_track) or iteration==maxit:
            feasibility1.append(feas_eq)
            feasibility2.append(feas_ineq)
            objective.append(np.sum(C.flatten()*X.flatten()))
            cur_iter.append(iteration)
            t.append(time.time()-start)
            print('{:03d} | {:.4e}| {:.4e}| {:.4e}|'.format(iteration, feasibility1[-1], feasibility2[-1],objective[-1]))
            
    return X, feasibility1, feasibility2, objective, cur_iter, t
```


```python
print("\nDataset = {}. \n".format(dataset))
print("\n f_opt = {}. \n".format(opt_val))

X_HCGM, f1_HCGM, f2_HCGM, obj_HCGM, iter_HCGM, time_HCGM = HCGM(kappa=p, maxit=np.int(1e3), beta0=5e1)

print("\nHCGM finished. Running time: {} seconds.".format(time_HCGM[-1]))
```

    
    Dataset = data/102n-insecta-ant-colony4-day10.mat. 
    
    
     f_opt = [[669.56459553]]. 
    
    001 | 5.2020e+03| 0.0000e+00| 1.3078e+03|
    002 | 5.2020e+03| 1.4711e+02| 4.3593e+02|
    003 | 1.7340e+03| 4.9038e+01| 1.2256e+03|
    004 | 1.7340e+03| 3.1488e+02| 7.3536e+02|
    005 | 1.0404e+03| 1.8893e+02| 1.1127e+03|
    006 | 1.0404e+03| 7.3615e+01| 7.9480e+02|
    008 | 7.4314e+02| 9.7239e+01| 8.1120e+02|
    010 | 5.7800e+02| 7.3494e+01| 9.4545e+02|
    013 | 4.0015e+02| 2.6638e+01| 1.0002e+03|
    017 | 3.0600e+02| 1.5553e+01| 8.8891e+02|
    023 | 2.2617e+02| 1.2735e+01| 8.2143e+02|
    030 | 1.7938e+02| 6.9966e+00| 7.2529e+02|
    040 | 1.3338e+02| 4.0708e+00| 7.0069e+02|
    053 | 9.8151e+01| 3.6254e+00| 7.1116e+02|
    070 | 7.5391e+01| 3.2214e+00| 6.7731e+02|
    093 | 5.5935e+01| 2.8352e+00| 6.8735e+02|
    123 | 4.2293e+01| 2.0105e+00| 6.8107e+02|
    163 | 3.1914e+01| 1.7190e+00| 6.7766e+02|
    217 | 2.3972e+01| 1.6077e+00| 6.7498e+02|
    287 | 1.8125e+01| 1.3839e+00| 6.7301e+02|
    381 | 1.3654e+01| 1.2024e+00| 6.7176e+02|
    505 | 1.0301e+01| 1.0463e+00| 6.7097e+02|
    670 | 7.7758e+00| 9.1208e-01| 6.6846e+02|
    889 | 5.8515e+00| 7.9142e-01| 6.7013e+02|
    1000 | 5.2072e+00| 7.4784e-01| 6.6869e+02|
    
    HCGM finished. Running time: 5032.813604831696 seconds.



```python
plot_func(iter_HCGM, f1_HCGM, f2_HCGM, obj_HCGM, X_HCGM, X_true)
```

1a3ba78dd2bdb6ab557c16915d3c2fa2b7040be8
ipynb
Jupyter Notebook
HW3/code/part3/Test_Sparsest_Cut.ipynb
3Represents/EE-556_MathsOfData
91790e214f2cd2f27a08f343ed89d050bc377d1e
CC0-1.0
HW3/code/part3/Test_Sparsest_Cut.ipynb
3Represents/EE-556_MathsOfData
91790e214f2cd2f27a08f343ed89d050bc377d1e
CC0-1.0
HW3/code/part3/Test_Sparsest_Cut.ipynb
3Represents/EE-556_MathsOfData
91790e214f2cd2f27a08f343ed89d050bc377d1e
CC0-1.0
Qwen/Qwen-72B
1. YES
2. YES