<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 <center>Structural Analysis and Visualization of Networks</center>
<center>Home Assignment #4
Step1:  Problems
Task 1* (For those who have not done that during the seminar)
On this seminar your are asked to implement simple community detection algorightm. It is called Markov Cluster Algorithm (MCL).
Markor Clustering Algorithm
Input
Step2:  The agorithm is desinged to transform the graph connectivity in such a way as to disconnect different communities and concentrate connectivity within one.
At each iteration the algorithm cacluates a new transition kernel, based on the $p$-hop importance score $u_{ij}$ of connectivity of nodes $i$ and $j$ in the original kernel.
The score is a nonlinear transformation of the $p$-hop connectivity probability $\pi^p_{sd} = \mathbb{P}( s\rightarrow_p d )$, which reflects the chances of reaching any node $d$ from a vertex $s$ in $p$-hops. 
The procedure consists of "pruning" and "inflation" steps. At the pruning step, the algorithm zeroes the $p$-hop transition prbabilities lower that a $\theta$. The rationale is, that if the random walk from $i$ is unlikely to end up in vertex $j$ in $p$ hops, then the $i\sim j$ path is not likely to connect vertices within one community. Thus the pairs of nodes $(i,j)$ with $p$-hop transition probability $\pi_{ij}^p$ lower than the threshold are forcefully disconnected. However, if the nodes $i$ and $j$ belong to the same community, the higher is the chance that a random walk reaches $j$ from $i$. Therefore zeroing negligible in effect severs weak between community connections.
At next step the remaining non-zero $\pi_{ij}^p$ are damped by a power transformation with parameter $\alpha>1$ to get the importance score $u_{ij} = \big(\pi^p_{ij}  \big)^\alpha 1_{\pi_{ij}^p\geq \theta}$. Finally the renormalisation step, makes the scores into one-hop transition probabilities for the next iteration. The a new transition kernel is thus $\pi_{ij}' \propto u_{ij}$ with the constraint $\sum_j \pi_{ij}' = 1$. The matrix $M$ in the algorithm is a transposed stochastic transition matrix
Step3:  Below are procedure for running the test, collecting the results and visualizing and printing them.
Step4:  Create $10\times 7$ grid of $(\alpha, p)$ pairs, with $\alpha$ spaced evenly apart from $1.1$ to $10$ on the $\log$-scale, and $p$ running from $1$ to $20$ with step $3$.
Step5:  This grid will be used throughout this task.
Let's load the first netork and run the MCL community detection procedure with default parameters.
Step6:  Below are the communities detected by the MCL with default parameters $\alpha=2$, $p=2$ and $\theta=10^{-8}$ for various leves of noise in the source data.
Step7:  Case
Step8:  Run the mcl for the first network.
Step9:  Run the mcl for the second network.
Step10:  And finally for the last one.
Step11:  Analysis
Display the number of communities detected (empty cells represent failure to detect anything, empty rows are omitted).
Step12:  As expected the higher the $\alpha$ the less perceptive the alogrithm is to the macrostructure. In contrast, high scouting range $p$ permits longer paths within a community, thereby making the iterative procedure absorb smaller communities into larger ones.  
Case
Step13:  Run the experiment
Step14:  Show the resulting tables
Step15:  The conclusions about the effects of $\alpha$ and $p$ match the $\theta=10^{-8}$ case. Higher $\theta$ narrowed the area, where the algorithm was able to detect any community.
Case
Step16:  Show the results
Step17:  The significance threshold $\theta$ severely reduced the range of $\alpha$-$p$ pairs, for which MCT detects anything at all. Within that region, however, parameters affect the detection similarly to the previous cases. 
Intersting cases
Let's have a look at some "intersting" parameter pairs for the case of $\theta = 10^{-8}$.
Step18:  Summary
In conclusion
Step19:  Residual code from the failed attempt.
Step20:  Maybe the all-pars shortest path matrix could yield some insight into the effects of the different parameter settings.
Step21:  <hr />
Task 2
Criterion for partitioning
Lets consider a problem of partitioning a  similarity matrix $A=\big(a_{ij}\big)_{i,j\in I}$ in two clusters non overlapping clsuters $A, B\subseteq I$ with $A\uplus B = I$. The matrix $A$ is a symmetric matrix of nonnegative elements.
The cut, or in other words, the between-cluster similarity, is defined as the sum of all edges across the boundary of the partition in a weighted graph with adjacency matrix $A$
Step22:  Load Yahoo Music network. Edges in this network appear if enough number of users have given ratings to both music bands. Note, that edges are weighted with similarity of the ratings.
Implement multilevel spectral recursive partitioning algorithm that was described during the lecture
Visualize community structure of the network and output some of the dense clusters (with interpretation, if you can)
You can load .mat files with the following commands
Step23:  Consider a graph $G=(V,E)$ given by a similarity matrix $A$. The basic idea of recursive spectral clsutering is to continue splitting $V$ in two clusters until either a cluster is depleted or a clique is detected. To cut a long story short, using the equivalent reformulation of the spectral clustering problem, the goal is to find an eigenvector corresponding to the second largest eigenvalue.
Step24:  Results
Run the recursive spectral clustering on the Yahoo music similarity data to uncover the hidden cluster strucutre.
Step25:  Recursive spectral clustering does indeed reveal certain connectivity structure in the network even for a restriction on depth of the recusion tree to at most seven levels ($\leq$128 clusters).
Now compute the density of each cluster in the full recursive tree of clusters.
Step26:  Oviously denser cluster have less nodes
Step27:  Lets visualize the densest clusters
Step28:  It is actaully quite difficult to analyse how well the clustering preformed, since we were not given the information on the artist's genre of music. However moderate googling revealed the genres of some of the densest clusters (#rank)
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import networkx as nx
from matplotlib import pyplot as plt
%matplotlib inline
from scipy.io import loadmat
import warnings
warnings.filterwarnings( 'ignore' )
Explanation: <center>Structural Analysis and Visualization of Networks</center>
<center>Home Assignment #4: Community Detection Algorithms
<center>Student: Nazarov Ivan</center>
<hr /> General Information
Due Date: 19.03.2015 23:59 <br >
Late submission policy: -0.2 points per day <br >
Please send your reports to &#108;&#101;&#111;&#110;&#105;&#100;&#46;&#101;&#46;&#122;&#104;&#117;&#107;&#111;&#118;&#64;&#103;&#109;&#97;&#105;&#108;&#46;&#99;&#111;&#109; and &#115;&#104;&#101;&#115;&#116;&#97;&#107;&#111;&#102;&#102;&#97;&#110;&#100;&#114;&#101;&#121;&#64;&#103;&#109;&#97;&#105;&#108;&#46;&#99;&#111;&#109; with message subject of the following structure:<br > [HSE Networks 2015] Nazarov Ivan HA4
Support your computations with figures and comments. <br >
If you are using IPython Notebook you may use this file as a starting point of your report.<br >
<br >
<hr \>
Preamble
The toolbox includes Numpy for computation, NetworkX for graph manipulation and algorithms, matplotlib for visualization and Scipy.IO for processing matlab objects.
End of explanation
def mcl_iter( A, p = 2, alpha = 2, theta = 1e-8, rel_eps = 1e-4, niter = 10000 ) :
## Convert A into a transition kernel: M_{ij} is the probability of making a transition from i to j.
    M = np.multiply( 1.0 / A.sum( axis = 1, dtype = np.float64 ).reshape(-1,1), A )
    i = 0 ; status = -1
    while i < niter :
        M_prime = M.copy( )
## Expansion step: M_{ij} is the probability of reaching a vertex j from i in p hops.
        M = np.linalg.matrix_power( M, p )
## Pruning: make paths with low transition probability into almost surely unused.
        M[ np.abs( M ) < theta ] = 0
## Inflation step: dampen the probabilites
        M = np.power( M, alpha )
## Renormalisation step: make the matrix into a stochastic transition kernel
        N = M.sum( axis = 1, dtype = np.float64 )
## If a nan is encountered, then abort
        if np.any( np.isnan( N ) ) :
            status = -2
            break
        M = np.multiply( 1.0 / N.reshape(-1,1), M )
## Convergence criterion is the L1 norm of relative divergence of transition probabilities
        if np.sum( np.abs( M - M_prime ) / ( np.abs( M_prime ) + rel_eps ) ) < rel_eps :
            status = 0
            break
## Advance to the next iteration
        i += 1
    return ( M, (status, i) )
Explanation: Problems
Task 1* (For those who have not done that during the seminar)
On this seminar your are asked to implement simple community detection algorightm. It is called Markov Cluster Algorithm (MCL).
Markor Clustering Algorithm
Input: Transition matrix $T = D^{-1}A$
Output: Adjacency matrix $M^$
1. Set $M = T$
2. repeat:
    3. Expansion Step: $M = M^p$ (usually $p=2$)
    4. Inflation Step: Raise every entry of $M$ to the power $\alpha$ (usualy $\alpha=2$)
    5. Renormalize: Normalize each row by its sum
    6. Prunning: Replace entries that are close to $0$ by pure $0$
7. until $M$ converges
8. $M^ = M$
As a result you should get a cluster matrix s.t. elements of the cluster correspont to nonzero elements of the columns of the matrix. 
Run this method for network 1, 2 and 3.
Play with the parameters ($p$, $\alpha$, zero tolerance), analyse the results
Solution
Though the specification above clearly states the order of operations within each iteration, the algorithm implemented below mpved "pruning" to before "inflation". First of all this uncovered more straightforward connection of the pruning parameter with transition probability. Secondly, as a side effect the final step within each iteration became "renormalization", which produces a row-stocastic matrix for the next iteration.
End of explanation
def extract_communities( M, lengths = True ) :
## It is extected that the MCL matrix detects communities in columns
    C = list( ) ; i0 = 0
    if np.any( np.isnan( M ) ) :
        return C
## Find all indices of nonzero elements
    r, c = np.where( M )
## Sort them by the column index and find the community sizes
    r = r[ np.argsort( c ) ]
    u = np.unique( c, return_counts = True )
    if np.sum( u[ 1 ] ) > M.shape[ 1 ] :
        return C
    if lengths :
        return u[ 1 ]
## Columns indices of nonzero entries are ordered, so we just need to
##  sweep across the sizes
    for s in u[ 1 ] :
## Row indices for a column with a nonzero element are the indices of
##  nodes in the community.
        list.append( C, r[ i0:i0+s ] )
        i0 += s
    return C
Explanation: The agorithm is desinged to transform the graph connectivity in such a way as to disconnect different communities and concentrate connectivity within one.
At each iteration the algorithm cacluates a new transition kernel, based on the $p$-hop importance score $u_{ij}$ of connectivity of nodes $i$ and $j$ in the original kernel.
The score is a nonlinear transformation of the $p$-hop connectivity probability $\pi^p_{sd} = \mathbb{P}( s\rightarrow_p d )$, which reflects the chances of reaching any node $d$ from a vertex $s$ in $p$-hops. 
The procedure consists of "pruning" and "inflation" steps. At the pruning step, the algorithm zeroes the $p$-hop transition prbabilities lower that a $\theta$. The rationale is, that if the random walk from $i$ is unlikely to end up in vertex $j$ in $p$ hops, then the $i\sim j$ path is not likely to connect vertices within one community. Thus the pairs of nodes $(i,j)$ with $p$-hop transition probability $\pi_{ij}^p$ lower than the threshold are forcefully disconnected. However, if the nodes $i$ and $j$ belong to the same community, the higher is the chance that a random walk reaches $j$ from $i$. Therefore zeroing negligible in effect severs weak between community connections.
At next step the remaining non-zero $\pi_{ij}^p$ are damped by a power transformation with parameter $\alpha>1$ to get the importance score $u_{ij} = \big(\pi^p_{ij}  \big)^\alpha 1_{\pi_{ij}^p\geq \theta}$. Finally the renormalisation step, makes the scores into one-hop transition probabilities for the next iteration. The a new transition kernel is thus $\pi_{ij}' \propto u_{ij}$ with the constraint $\sum_j \pi_{ij}' = 1$. The matrix $M$ in the algorithm is a transposed stochastic transition matrix: each entry is the probability $\pi_{ij} = \mathbb{P}( i\rightarrow j )$.
Based on this description the following hypothese can be formulated with respect to the effects of the parameters $\alpha$, $p$ and $\theta$:
* $\theta\in(0,1)$ -- determines the level beyond which the path transition probabilities are considered important enough to contibute to a community. A high $\theta$ makes the detection less sensitive, and beyond some threshold may cause the algorithm to fail to detect anything. Lower significance thresholdmakes the MCL more likely to find larger communities;
* $p\geq1$ integer: determines how far the periphery of a community may spread. The higher the $p$ the the less sensitive is community detection to the separation of nodes, and the more likely it is to find larger communities;
* $\alpha>1$ real: governs the damping. With high $\alpha$ highly probable transitions stay probable in the final transition kernel $\pi_{ij}'$, while moderate become negligible. That is why with low $\alpha$ larger communities are detected, whereas high damping exponent tends to detect more concentrated, tighter communities, is any at all.
End of explanation
def mcl_test( A, a_grid, p_grid, theta = 1e-8 ) :
## Run the grid test
    res = [ mcl_iter( A, p = p, alpha = a, theta = theta )
        for p in p_grid for a in a_grid ]
## Extract the results
## Get the number of communities
    NC = np.array( [ len( extract_communities( C ) )
        for C,(s,i) in res ], dtype = np.int ).reshape( len( p_grid ), -1 )
## Extract the number of iterations
    NI = np.array( [ i for C,(s,i) in res ], dtype = np.int ).reshape( len( p_grid ), -1 )
    return NI, NC
## Not a good way of printing tables 
def show_table( S, r, c ): 
    print "    p\\a\t", "\t".join( len( c )*[ "%#4.3g" ] ) % tuple( c )
    for j in xrange( S.shape[0] ) :
        if np.all( S[j,:] == 0 ) :
            break
        row = [ "%#2d"%(v) if v > 0 else "  " for v in S[j,:] ]
        print "%#6d\t"%( r[ j ] ),  "\t".join( len( c )*[ "%s" ] ) % tuple( row )
    
## Produce a visually appealling picture of the adjacency
##  matrix and community detection results 
def show_network( A, C, title = "" ) :
    plt.spy( A, color = "gray", markersize = .5 )
    plt.spy( C, color = "magenta", markersize = 5 )
    if title : plt.title( title )
Explanation: Below are procedure for running the test, collecting the results and visualizing and printing them.
End of explanation
alpha_grid, p_grid = np.logspace( 4.5e-3, 1, num = 10, dtype = np.float ), np.arange( 2, 21, dtype = np.int )
Explanation: Create $10\times 7$ grid of $(\alpha, p)$ pairs, with $\alpha$ spaced evenly apart from $1.1$ to $10$ on the $\log$-scale, and $p$ running from $1$ to $20$ with step $3$.
End of explanation
files = [ 'network1.mat', 'network2.mat', 'network3.mat' ]
titles = [ "Noiseless network: "+files[0], "Noisier network: "+files[1], "Noisiest network: "+files[2] ]
matrices = [ np.array( loadmat( './data/hw4/'+f )[ 'A' ], np.float ) for f in files ]
Explanation: This grid will be used throughout this task.
Let's load the first netork and run the MCL community detection procedure with default parameters.
End of explanation
plt.figure( figsize=(18,6) )
plt.subplot( 131 )
C0, _ = mcl_iter( matrices[ 0 ] )
show_network( matrices[ 0 ], C0, titles[0] )
plt.subplot( 132 )
C1, _ = mcl_iter( matrices[ 1 ] )
show_network( matrices[ 1 ], C1, titles[1] )
plt.subplot( 133 )
C2, _ = mcl_iter( matrices[ 2 ] )
show_network( matrices[ 2 ], C2, titles[2] )
plt.show( )
Explanation: Below are the communities detected by the MCL with default parameters $\alpha=2$, $p=2$ and $\theta=10^{-8}$ for various leves of noise in the source data.
End of explanation
theta = 1.0E-8
Explanation: Case: $\theta = 10^{-8}$
Let's being with the prunnig threshold set to almost neglibile value.
End of explanation
I11, S11 = mcl_test( matrices[ 0 ], alpha_grid, p_grid, theta = theta )
Explanation: Run the mcl for the first network.
End of explanation
I12, S12 = mcl_test( matrices[ 1 ], alpha_grid, p_grid, theta = theta )
Explanation: Run the mcl for the second network.
End of explanation
I13, S13 = mcl_test( matrices[ 2 ], alpha_grid, p_grid, theta = theta )
Explanation: And finally for the last one.
End of explanation
print files[ 0 ]
show_table( S11, p_grid, alpha_grid )
print files[ 1 ]
show_table( S12, p_grid, alpha_grid )
print files[ 2 ]
show_table( S13, p_grid, alpha_grid )
Explanation: Analysis
Display the number of communities detected (empty cells represent failure to detect anything, empty rows are omitted).
End of explanation
theta = 1.0E-4
Explanation: As expected the higher the $\alpha$ the less perceptive the alogrithm is to the macrostructure. In contrast, high scouting range $p$ permits longer paths within a community, thereby making the iterative procedure absorb smaller communities into larger ones.  
Case: $\theta = 10^{-4}$
Let's set a higher threshold: $\theta = 10^{-4}$:
End of explanation
I21, S21 = mcl_test( matrices[ 0 ], alpha_grid, p_grid, theta = theta )
I22, S22 = mcl_test( matrices[ 1 ], alpha_grid, p_grid, theta = theta )
I23, S23 = mcl_test( matrices[ 2 ], alpha_grid, p_grid, theta = theta )
Explanation: Run the experiment
End of explanation
print files[ 0 ]
show_table( S21, p_grid, alpha_grid )
print files[ 1 ]
show_table( S22, p_grid, alpha_grid )
print files[ 2 ]
show_table( S23, p_grid, alpha_grid )
Explanation: Show the resulting tables
End of explanation
theta = 0.01
I31, S31 = mcl_test( matrices[ 0 ], alpha_grid, p_grid, theta = theta )
I32, S32 = mcl_test( matrices[ 1 ], alpha_grid, p_grid, theta = theta )
I33, S33 = mcl_test( matrices[ 2 ], alpha_grid, p_grid, theta = theta )
Explanation: The conclusions about the effects of $\alpha$ and $p$ match the $\theta=10^{-8}$ case. Higher $\theta$ narrowed the area, where the algorithm was able to detect any community.
Case: $\theta = 10^{-2}$
Finally consider the case of $1\%$-probability threshold
End of explanation
print files[ 0 ]
show_table( S31, p_grid, alpha_grid )
print files[ 1 ]
show_table( S32, p_grid, alpha_grid )
print files[ 2 ]
show_table( S33, p_grid, alpha_grid )
Explanation: Show the results
End of explanation
def cases( n, pairs, theta = 1.0E-8 ) :
    plt.figure( figsize=(18,6) )
    for i, (a, p) in enumerate( pairs, 1 ) :
        plt.subplot( int( "1" + str( len( pairs ) ) + str( i ) ) )
        C, _ = mcl_iter( matrices[n], alpha = alpha_grid[a], p = p_grid[p], theta = theta )
        show_network( matrices[n], C, titles[n] + ": (a=%.2f, p=%d)" % ( alpha_grid[a], p_grid[p] ) )
    plt.show( )
## The first network
cases( 0, [ (1, 0), (9, 1), (9, 15) ] )
## Second network
cases( 1, [ (8, 1), (9, 1), (9, 2) ] )
## Third network
cases( 2, [ (6, 1), (9, 1), (8, 2) ] )
Explanation: The significance threshold $\theta$ severely reduced the range of $\alpha$-$p$ pairs, for which MCT detects anything at all. Within that region, however, parameters affect the detection similarly to the previous cases. 
Intersting cases
Let's have a look at some "intersting" parameter pairs for the case of $\theta = 10^{-8}$.
End of explanation
## Compute the scout probabilites: an m-scout probability is the
##  chance of the first ever visit to j of a random walk from i
##  taking place no later than m-th step.
def scout( P, n = 3, cumulative = True ) :
    pi = P.copy() ; Sc = pi.copy() ; i = 1
    while i < n :
## Get the probability of a j-j loop with i hops
        pi_jj = pi.diagonal( ).copy( )
## The probability of the first ever visit to j from i being
##  on the m-th step: V^m_{ij} = \sum_{k\neq j} P_{ik} V^{m-1}_{kj}
        pi = P.dot( pi ) - np.multiply( P, pi_jj )
## The scout probability: the chance that the first visit is earlier
##  or exactly at the m-th step.
        Sc += pi
        i += 1
    return SC if cumulative else pi
Explanation: Summary
In conclusion:
* the size of communitites detected is inversely related to $\alpha$;
* $p$ regulates the scouting radius of the underlying random walk, and higher values make the procedure less sensitive to the micro-level topology of the conncetivity graph, which causes detection of coarser communities;
* The significance level $\theta$ is inversely proportional to the sensitivity of the procedure: higher values make MCL ignore more probable paths between nodes, which leads to breaking up of large communities into finer, more concentrated ones. The major effect of significance threshold is on the size of the region of $\alpha$ and $p$, where the procedure detect anything at all;
* the MCL is extrmely sensitive to noise or missing links in the connectivity data.
Appendix
Below is a failed attempt to alter the MCL algorithm.
The issue with the $p$-hop probability is that $\pi^p_{ij}$ reflects the total probability of arbitrary $p$-hop $i\rightarrow j$ paths, even those that use $j$ as a transitory vertex.
However one could also propose the following heuristic: if a randomly wandering agent, who sets off to purposefully scout the graph for the member $j$, reaches its goal within $m$ hops, then $i-j$ connection must be strong within the graph. Thus the air $i$ and $j$ should lie withi the same community.
Such scouting probability is computed recursively via the following formula:  
$$\sigma^p_{ij} = \sum_{k_1\neq j}\ldots \sum_{k_{p-1}\neq j} \pi_{ik_1}\pi_{k_1k_2}\ldots \pi_{k_{p-2}k_{p-1}} \pi_{k_{p-1}j} = \sum_{k\neq j} \pi_{ik} \sigma^{p-1}{kj} = \big[\pi \sigma^{p-1}\big]{ij} - \pi_{ij}\sigma^{p-1}{jj}$$
which implies that in matrix notation:
$$\sigma^p = \pi \sigma^{p-1} - \pi \odot \bigg(\begin{smallmatrix}(\sigma^{p-1}{ii})i\\ldots\(\sigma^{p-1}{ii})_i\end{smallmatrix}\bigg)$$
where $\odot$ is the element-wise matrix product.
End of explanation
if False :
    A = matrices[ 0 ]
    P = np.multiply( 1.0 / A.sum( axis = 1, dtype = np.float64 ).reshape(-1,1), A )
    S = scout(P, 100, False)
    l, v = np.linalg.eig( S )
    from scipy.cluster.vq import kmeans, vq
    K=4
    V = v[:,np.argsort( l )[ -K: ] ]
    O = vq(V,kmeans( V, K )[ 0 ])
    i = np.argsort( O[ 0 ] )
    plt.imshow( np.log(S[np.ix_(i,i)]))
    plt.colorbar( )
Explanation: Residual code from the failed attempt.
End of explanation
def floyd_warshall( A ) :
## Create a matrix object (not just a 2D array -- different broadcasting properties)
    pi = np.matrix( A, dtype = np.float, copy = True )
## Fill as of yet unreachable vertices
    pi[ pi != 1 ] = np.inf
## And show that the shortest path to oneself is staying.
    np.fill_diagonal( pi, 0 )
## For each transitory vertex
    for v in xrange( pi.shape[ 0 ] ) :
## Decide which is faster: to use a path without it, or to pass through it.
        np.minimum( pi, pi[:,v] + pi[v,:], pi )
## Return the shortest path matrix
    return pi
Explanation: Maybe the all-pars shortest path matrix could yield some insight into the effects of the different parameter settings.
End of explanation
import numpy as np
import networkx as nx
from matplotlib import pyplot as plt
%matplotlib inline
from scipy.io import loadmat
import warnings
warnings.filterwarnings( 'ignore' )
import scipy.io
import scipy.sparse as spma
import scipy.sparse.linalg as spla
Explanation: <hr />
Task 2
Criterion for partitioning
Lets consider a problem of partitioning a  similarity matrix $A=\big(a_{ij}\big)_{i,j\in I}$ in two clusters non overlapping clsuters $A, B\subseteq I$ with $A\uplus B = I$. The matrix $A$ is a symmetric matrix of nonnegative elements.
The cut, or in other words, the between-cluster similarity, is defined as the sum of all edges across the boundary of the partition in a weighted graph with adjacency matrix $A$:
$$\text{cut}(A, B) = \sum_{i\in A}\sum_{j\in B} a_{ij}$$
Let $u$ be a partition vector with $u_i = +1$ if $i\in A$ and $-1$ otherwise (if $i\notin A$, or $i\in B$). Then the cut can be represented as
$$\text{cut}(A,B) = \frac{1}{8} \sum_{i,j} a_{ij}(u_i-u_j)^2$$
since edges are counted twice, and $(u_i-u_j)^2 = 4$ whenever $i,j\in A$ or $i,j\in B$, and $0$ otherwise.
Furthermore one has
$$\sum_{i,j} a_{ij}(u_i-u_j)^2 = \sum_{i,j}a_{ij}(u_i^2+u_j^2) - 2\sum_{i,j}a_{ij}u_iu_j
= 2\Big( \sum_i u_i^2 \sum_j a_{ij} - u'Au \Big)
= 2\Big( \sum_i u_i \delta_i u_i - u'Au \Big) = 2 u'\big(D-A\big)u$$
where $D=\text{diag}\big(\delta_i\big)$ and $\delta_i = \sum_j a_{ij}$.
Now since $u'u = n$, the cut is equivalently given by
$$\text{cut}(A, B) = \frac{n}{4} \frac{u'Lu}{u'u}$$
where the matrix $L = D-A$ is also known as the Laplacian.
It can be shown that $L\mathbf{1} = 0$, and the number of zero eigenvalues of $L$ is equal to the number of connected components in a graph represented by $A$.
Consider a weighted cut criterion:
$$\mathcal{Q}(A,B) = \frac{\text{cut}(A,B)}{W_A} + \frac{\text{cut}(A,B)}{W_B} = \frac{W_A + W_B}{W_A W_B}\text{cut}(A,B)$$
where $W_A = \sum_{i\in A} \omega_i$ and $W = \text{diag}\big(\omega_i\big)$ -- a full rank diagonal matrix. In the case of normalised cut, the weighting matrix $W$ equal to $D$.
Define a new partition vector $q$ by $q_i = \sqrt{\frac{W_B}{W_A}}$ for $i\in A$ and $q_i = - \sqrt{\frac{W_A}{W_B}}$ otherwise. Then
$$q'W\mathbf{1} = \sum_i q_i \omega_i = \sqrt{\frac{W_B}{W_A}}\sum_{i\in A} \omega_i - \sqrt{\frac{W_A}{W_B}}\sum_j \omega_j = \sqrt{W_AW_B}-\sqrt{W_AW_B}=0$$
and
$$q'Wq = \sum_i q_i^2 \omega_i = \frac{W_B}{W_A}\sum_{i\in A} \omega_i^2 + \frac{W_A}{W_B}\sum_j \omega_j^2 = W_A + W_B$$
Now, let's express $q$ in terms of $u$:
$$q_i = \frac{W_B}{\sqrt{W_A W_B}} \text{ if } i\in A \text{ and } q_i = - \frac{W_A}{\sqrt{W_A W_B}} \text{ otherwise}$$
Whence the following expression follows:
$$q = \frac{1}{2 \sqrt{W_A W_B}} \Big( W_B (u+\mathbf{1}) + W_A (\mathbf{1}-u) \Big)
= \frac{1}{2 \sqrt{W_A W_B}} \Big( (W_B+W_A) u + (W_B-W_A) \mathbf{1} \Big)$$
because  $u_i=+1$ if $i\in A$ and $u_i=-1$ otherwise.
Since $L\mathbf{1} = D\mathbf{1} - A\mathbf{1} = 0$, one the following expresssion is true:
$$q'Lq
= \frac{1}{4 W_A W_B} \Big( (W_B+W_A) u +(W_B-W_A) \mathbf{1} \Big)' L \Big( (W_B+W_A) u +(W_B-W_A) \mathbf{1} \Big)
= \frac{(W_B+W_A)^2}{4 W_A W_B} u'L u$$
Rearranging:
$$\frac{q'Lq}{q'Wq} = \frac{W_B+W_A}{W_A W_B} \frac{u'L u}{4} = \frac{W_A + W_B}{W_A W_B}\text{cut}(A,B) = \mathcal{Q}(A,B)$$
This demonstrates that the problem of finding $A,B\subseteq I$ such that $\mathcal{Q}(A, B)$ is minimized is equivalent to finding such a partition vector $q$ that the ratio $\frac{q'Lq}{q'Wq}$ is minimal subject to $q'W\mathbf{1}=0$ constraint.
References:
Dhillon, I. S. (2001, August). "Co-clustering documents and words using bipartite spectral graph partitioning". In Proceedings of the seventh ACM SIGKDD international conference on Knowledge discovery and data mining (pp. 269-274). ACM.
Optimial solution
The problem
$$\frac{q'Lq}{q'Wq}\rightarrow \min_{q\neq 0} \text{ subject to } q'W\mathbf{1} = 0$$
for a real valued vector $q$ can easily be solved. Indeed, the Lagrangian of the problem is
$$\mathcal{L} = \frac{q'Lq}{q'Wq} - 2q'W\mathbf{1}\mu \rightarrow \min$$
yields the First Order Conditions 
$$\frac{\partial}{\partial q} \mathcal{L} = \frac{ 2Lq (q'Wq) - 2Wq (q'Lq) }{(q'Wq)^2} - W'\mathbf{1} (2\mu) = 0$$
and
$$\frac{\partial}{\partial \mu} \mathcal{L} = 2q'W\mathbf{1} = 0$$
The FOC for $q$ simplifies to $Lq (q'Wq) - Wq (q'Lq) = W'\mathbf{1} \mu (q'Wq)^2$. Left multiplying by $\mathbf{1}'$ results in
$$\mathbf{1}'Lq (q'Wq) - \mathbf{1}'Wq (q'Lq) = \mathbf{1}'W\mathbf{1} \mu (q'Wq)^2$$
However $\mathbf{1}'L = 0'$, whence $-\mathbf{1}'Wq (q'Lq) = \mathbf{1}'W\mathbf{1} \mu (q'Wq)^2$ and the FOC for $\mu$ finally gives $\mathbf{1}'W\mathbf{1} \mu (q'Wq)^2 = 0$. This implies $\mu=0$, as $W$ is diagonal with positive values, whence the first order conditions simplify to $Lq = Wq \lambda$ for $\lambda = \frac{q'Lq}{q'Wq}$.
Thus solution of the relaxed problem is equivalent to finding a generalsed eigenvector $q$ of $Lq=Wq\lambda$ corresponfing to the smallest generalised eigenvalue subject to the orthogonality constraint $q'W\mathbf{1} = 0$. 
Therefore the relaxed weighted cut minimization problem is equivalent to finding the smallest eigenvalues corresponding to an eigenvector in the subspace orthogonal to $W\mathbf{1}$.
Suppose $x$ is some solution to the generalised eigenvalue problem $Lx=Wx\lambda$. Since $L$ and $W$ are positive-semi deifinite, $x'Lx,x'Wx\geq0$, $x'Lx = x'Wx\lambda$ implies that $\lambda\geq0$.
Note that $q=\mathbf{1}$ is not a solution, even thought it is a generalised eigenvalue of $L$ and $W$ pair. Indeed, since $W$ is a sitriclty positive definite matrix, $\mathbf{1}'W\mathbf{1} \neq 0$ violating the orthogonality constraint.
In fact the problem needs to be reformulated so that the taks would become finding the largest or the second larges eigenvalue, because almsot all numerical methods for finding eigenvalue-egienvector pairs are optimised for computing the largest eigenvalues.
For the case of normalised cut this reformulation is trivial. Since the equation $Lx = Dx\lambda$ is eqivalent to
$Dx - Ax = Dx\lambda$, whence $Ax=Dx(1-\lambda)=Ax$. Thue the eqivalent formulation of the problem is to find the largest eienvalues and the corresponding generalised eigenvector of the pair $(A, D)$ orthogonal to $D\mathbf{1}$:
$$\text{ find largest } \mu \text{ with } Ax=Dx\mu \text{ subject to } x'D\mathbf{1} = 0 \text{ and } x\neq 0$$
The egienvalues of the reformulated problem are nonegative, since the both matrices are positive semidefinite. Therfore $\mu\geq 0$ and $\lambda\geq0$ with the ovious relationship $\mu = 1-\lambda$, imply that $\mu,\lambda\in [0,1]$.
Finally, since $D$ is invertible, the problem simplifies to finding the ordinary eigenvector of $D^{-1}A$ corresponding to the second largest eigenvalue.
Problem solution
End of explanation
data = scipy.io.loadmat('./data/hw4/music_data.mat')
A = spma.csc_matrix( data[ 'A' ], dtype = np.float )
Explanation: Load Yahoo Music network. Edges in this network appear if enough number of users have given ratings to both music bands. Note, that edges are weighted with similarity of the ratings.
Implement multilevel spectral recursive partitioning algorithm that was described during the lecture
Visualize community structure of the network and output some of the dense clusters (with interpretation, if you can)
You can load .mat files with the following commands:
End of explanation
def cluster( A, T = 100, Q = None, _index = None, mincut = False, depth = float( "inf" ), density_threshold = .05 ) :
## If the recursion depth is exceeded return
    if depth <= 0 :
        return np.arange( A.shape[ 0 ] )
## Create master indices if necessary
    if _index is None :
        _index = np.arange( A.shape[ 0 ] )
## Compute the global similarity of each vertex/element
    deg = A.sum( axis = 1 ).getA1( )
## Detect non-isolated items
    nz, zz = np.where( deg != 0 )[ 0 ], np.where( deg == 0 )[ 0 ]
    if len( nz ) < T :
        return np.arange( len( deg ) )
## This fiddling with tocsr() and tocsc() makes slicing faster, since format
##  conversions are extremely fast.
    S = A[:,nz].tocsr()[nz,:].tocsc()
    try :
        if mincut :
## Compute the unnormalised laplacian
            L = spma.diags( deg[ nz ], offsets = 0 ) - S
## Get the eigenvector of the second least eigenvalue
            l, e = spma.linalg.eigs( L, k = 2, which = 'SM',
                v0 = np.ones( L.shape[ 1 ], np.float ) )
            e = e[ :, np.argmax( l ) ].real
        else :
## Compute the stochastic transition kernel of non-isolated vertices for the
##  normalised cut problem.
            L = spma.diags( 1.0 / deg[ nz ], offsets = 0 ).dot( S )
## Find the eigenvector corresponding to the 2nd largest eigenvalue
            l, e = spma.linalg.eigs( L, k = 2, v0 = np.ones( L.shape[ 1 ], np.float ) )
## Get the real part of the second largest eigenvector
            e = e[ :, np.argmin( l ) ].real
    except Exception, e:
#             print "size = %d : \t%s" %( A.shape[ 0 ], e.message )
## Return the global cluster if the eigenvalue computation failed to converge
        return np.arange( len( deg ) )
## Set the threshold: use the zero threshold
    t = 0
## Separate the items in two sets: left(n) and right (p)
    n, p = np.where( e <= t )[ 0 ], np.where( e > t )[ 0 ]
    N, P = S[:,n].tocsr()[n,:].tocsc( ), S[:,p].tocsr()[p,:].tocsc( )
## Compute the densities of the halves
    nd, pd = N.nnz, P.nnz
    nw, pw = len( n ) * ( len( n ) - 1.0 ), len( p ) * ( len( p ) - 1.0 )
## If there is enough elements in a set, split it.
    if ( len( p ) > T ) and ( density_threshold * pw > pd ) :
        p = p[ cluster( P, T = T, Q = Q, _index = _index[ nz[ p ] ],
            mincut = mincut, depth = depth - 1, density_threshold = density_threshold ) ]
    if ( len( n ) > T ) and ( density_threshold * nw > nd ) :
        n = n[ cluster( N, T = T, Q = Q, _index = _index[ nz[ n ] ],
            mincut = mincut, depth = depth - 1, density_threshold = density_threshold ) ]
## Update the queue of clusters
    if Q is not None :
        Q.append( _index[ nz[ p ] ] )
        Q.append( _index[ nz[ n ] ] )
        Q.append( _index[ zz ] )
#         Q.append( _index[ np.concatenate( ( zz, nz[ n ], nz[ p ] ) ) ] )
## Reorder the clusters so that the denser groups are shifted to the left
    if nd*pw > pd*nw : p, n = n, p
    return np.concatenate( ( zz, nz[ n ], nz[ p ] ) )
Explanation: Consider a graph $G=(V,E)$ given by a similarity matrix $A$. The basic idea of recursive spectral clsutering is to continue splitting $V$ in two clusters until either a cluster is depleted or a clique is detected. To cut a long story short, using the equivalent reformulation of the spectral clustering problem, the goal is to find an eigenvector corresponding to the second largest eigenvalue.
End of explanation
fig = plt.figure( figsize = ( 16, 16 ) )
axs = fig.add_subplot(1, 1, 1, axisbg = 'black')
I = cluster( A, T = 5, depth = 3, density_threshold = .2 )
axs.spy( A[:,I].tocsr( )[I,:], marker = '.', markersize = 3, precision = 0, alpha = .35, color = 'magenta' )
I = cluster( A, T = 5, depth = 5, density_threshold = .2 )
axs.spy( A[:,I].tocsr( )[I,:], marker = '.', markersize = 2, precision = 0, alpha = .25, color = 'green' )
I = cluster( A, T = 5, depth = 7, density_threshold = .2 )
axs.spy( A[:,I].tocsr( )[I,:], marker = '.', markersize = 1, precision = 0, alpha = .25, color = 'cyan' )
I = cluster( A, T = 5, depth = np.inf, density_threshold = .2 )
axs.spy( A[:,I].tocsr( )[I,:], marker = '.', markersize = 1, precision = 0, alpha = .25, color = 'gold' )
Explanation: Results
Run the recursive spectral clustering on the Yahoo music similarity data to uncover the hidden cluster strucutre.
End of explanation
cluster_tree = list( )
I = cluster( A, T = 5, Q = cluster_tree, depth = np.inf, density_threshold = .2 )
## Compute the density
den = np.zeros( len( cluster_tree ), np.float )
dia = A.diagonal( )
for i, c in enumerate( cluster_tree ) :
## Omit residual clusters
    if len( c ) < 2 : continue
    within = A[:,c].tocsr()[c,:].nnz # sum( ) - dia[c].sum( )
    weight = len( c ) * ( len( c ) - 1.0 )
    den[ i ] = within / weight
## Reorder
order = np.argsort( den )
top = order[ ~np.isnan( den[ order ] ) ][::-1]
Explanation: Recursive spectral clustering does indeed reveal certain connectivity structure in the network even for a restriction on depth of the recusion tree to at most seven levels ($\leq$128 clusters).
Now compute the density of each cluster in the full recursive tree of clusters.
End of explanation
plt.loglog( den[ top ], [ len( cluster_tree[ i ] ) for i in top ] )
plt.xlabel( 'Density (log)' ) ; plt.ylabel( 'Size (log)' )
Explanation: Oviously denser cluster have less nodes:
End of explanation
fig = plt.figure( figsize = ( 16, 16 ) )
axs = fig.add_subplot(1, 1, 1, axisbg = 'black')
axs.spy( A[:,I].tocsr( )[I,:], marker = '.', markersize = 3, precision = 0, alpha = .50, color = 'cyan' )
H = spma.csc_matrix( A.shape, dtype = np.int8 )
for f in top[ np.where( np.logical_and( den[ top ] > .05, den[ top ] <= .2 ) ) ] :
    C = cluster_tree[ f ]
    H[ np.ix_( C, C ) ] = 1
axs.spy( H[:,I].tocsr( )[I,:], marker = '.', markersize = 1, precision = 0, alpha = .25, color = 'magenta' )
H = spma.csc_matrix( A.shape, dtype = np.int8 )
for f in top[ np.where( den[ top ] > .30 ) ] :
    C = cluster_tree[ f ]
    H[ np.ix_( C, C ) ] = 1
axs.spy( H[:,I].tocsr( )[I,:], marker = '.', markersize = 1, precision = 0, alpha = .15, color = 'yellow' )
Explanation: Lets visualize the densest clusters: the final spectral image is in vyan, clusters with density $\delta$ between $0.05$ and $0.2$ are in magenta, and the densest clusters ($\delta\geq0.3$) are in gold.
End of explanation
for n, i in enumerate( top, 1 ) :
    c = cluster_tree[ i ]
    if len( c ) and den[ i ] > .3 :
        print "%#4d (%#4d, %0.3f)\t"%(n,len( c ), den[ i ]), ", ".join(sorted( [ s.strip() for s in data['artists'][ c ] ] ) ), "\n"
Explanation: It is actaully quite difficult to analyse how well the clustering preformed, since we were not given the information on the artist's genre of music. However moderate googling revealed the genres of some of the densest clusters (#rank):
*  2 -- country;
* 48 -- Heavy metal cl&uuml;ster (strangely it includes the Kings of Metal Manowar, but not Iron Maiden);
*  3 -- metal/hard rock;
*  6 -- jazz;
* 10 -- blues;
* 42 -- rap;
*  7 -- R&B;
* 51 -- christian pop/rock/rap music.
End of explanation
<END_TASK>