--- title: "admmDensestSubmatrix: ADMM for the densest submtarix problem " author: "Brendan Ames, Polina Bombina" date: "`r Sys.Date()`" output: rmarkdown::html_vignette code_folding: hide vignette: > %\VignetteIndexEntry{admmDensestSubmatrix: ADMM for the densest submtarix problem} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` \usepackage{amsmath} \usepackage{amssymb} \usepackage{amstext} \usepackage{mathtools} \newcommand{\tr}{\operatorname{Tr}} \newcommand{\rank}{\operatorname{rank}} \newcommand{\st}{\operatorname{s.t.}} \necommand{\P}{\operatorname{P}} # Introduction This is the R-package accompanying the paper ([Convex optimization for the densest subgraph and densest submatrix problems](https://arxiv.org/abs/1904.03272). The problem of identifying a dense submatrix is a fundamental problem in the analysis of matrix structure and complex networks. This package provides tools for identifying the densest submatrix of a given graph using first-order optimization methods. See the tutorials below to get started. # The densest submatrix problem Let $[M] = \{1,2,\dots, M\}$ for each positive integer $M$. Given a matrix $\mathbf{A} \in R^{M\times N}$, the densest $m\times n$-submatrix problem seeks subsets $\bar U \subseteq {[M]}$ and $\bar V \subseteq {[N]}$ of cardinality $|\bar U|=m$ and $|\bar V| = n$, respectively, such that the submatrix $\mathbf{A}{[\bar U, \bar V]}$ with rows index by $\bar U$ and columns indexed by $\bar V$ contains the maximum number of nonzero entries. That is, the densest $m\times n$-submatrix problem seeks the densest $m\times n$-submatrix of $\mathbf{A}$. The densest $m\times n$-submatrix problem can be formulated as: $$\begin{equation} \min_{\mathbf{X}, \mathbf{Y} \in { \{0,1\}}^{M\times N} } {\tr(\mathbf{Y} \mathbf{e} \mathbf{e}^T): \mathrm{P}_{\Omega}(\mathbf{X}-\mathbf{Y}) = \mathbf{0}, \tr(\mathbf{X} \mathbf{e} \mathbf{e}^T) = mn, \rank (\mathbf{X}) = 1 }, \end{equation}$$ where * $\mathrm{P}_{\Omega}$ is the projection onto the index set of zero entries of matrix $\mathbf A$; * $\tr$ is the matrix trace function; * $\Omega$ is the index set of zero entries of matrix $A$; * $\mathbf{X}$ is rank-one matrix with $mn$ nonzero entries; * $\mathbf{Y}$ is used to count the number of disagreements between $\mathbf A$ and $\mathbf X$; * $\mathbf{e}$ - all-ones vector. Unfortunately, optimization problems involving rank and binary constraints are intractable in general. Relaxing the rank constraint with a nuclear norm penalty term, $\|\mathbf{X} \|_* = \sum_{i=1}^N \sigma_i(\mathbf{X})$, i.e., the sum of the singular values of matrix, and the binary constraints with box constraints yields the convex problem: $$\begin{align} \min \; & \|\mathbf{X} \|_* + \gamma \tr(\mathbf{Y} \mathbf{e} \mathbf{e}^T) \\ \st \; & \tr(\mathbf{X} \mathbf{e} \mathbf{e}^T) = mn, \\ & \mathrm{P}_\Omega(\mathbf{X} - \mathbf{Y}) = \mathbf{0}, \\ & \mathbf{Y} \ge \mathbf{0}, \\ & \mathbf{0} \le \mathbf{X} \le \mathbf{e} \mathbf{e}^T, \end{align}$$ where $\gamma >0$ is a regularization parameter chosen to tune between the two objectives. It can be shown that the relaxation is exact when binary matrix $\mathbf{A}$ contains a single, relatively large dense $m\times n$ block. For more information, see ([Convex optimization for the densest subgraph and densest submatrix problems](https://arxiv.org/abs/1904.03272) # Alternating Direction Method of multipliers for densest submatrix problem The alternating direction method of multipliers (ADMM) has been succesfully used in a broad spectrum of applications. The ADMM solves convex optimization problems with composite objective functions subject to equality constraints. We direct the reader to Prof. Stephen Boyd’s website ([ADMM](http://stanford.edu/~boyd/papers/admm_distr_stats.html)) for a more thorough discussion of the ADMM. To apply ADMM to our problem, we introduce artificial variables $\mathbf{Q}$, $\mathbf{W}$ and $\mathbf{Z}$ to obtain the equivalent optimization problem: $$\begin{align} \min \; & \|\mathbf{X} \|_* + \gamma \|\mathbf{Y} \|_1 +{1}_{\Omega_Q}(\mathbf{Q})+{1}_{\Omega_W}(\mathbf{W})+{1}_{\Omega_Z}(\mathbf{Z})\\ \st \; & \mathbf{X}-\mathbf{Y}=\mathbf{Q},\mathbf{X}=\mathbf{W}, \mathbf{X}=\mathbf{Z} \end{align}$$ where * $\Omega_Q = \{\, \mathbf{Q}\in R^{M\times N} \mid P_{\tilde{N}}(Q)=0 \, \}$, * $\Omega_W =\{\, \mathbf{W}\in R^{M\times N} \mid \mathbf{e}^T\mathbf{W} \mathbf{e}=mn \, \}$, * $\Omega_Z =\{\, \mathbf{Z}\in R^{M\times N} \mid {Z}_{ij}\leq 1 \forall (i,j)\in M\times N \, \}$. Here ${1}_{S}: R^{M\times M} \rightarrow \left \{0,+\infty \right \}$ is the indicator function of the set $S \subseteq R^{M\times N}$, such that ${1}_S(\mathbf{X})=0$ if $\mathbf{X}\in S$, and $+\infty$ otherwise. Since our objective function is separable, we iteratively solve this optimization program using the ADMM. The basic idea is to rotate through $3$ steps: 1. minimize the augmented Lagrangian over primal variables, 2. update dual variables usng the updated primal variables, 3. calculate primal and dual residuals. Interested readers are referred to ([Convex optimization for the densest subgraph and densest submatrix problems](https://arxiv.org/abs/1904.03272). We include a summary of the algorithm below. ![](ALG.png){width=600px} # Examples We test this package on two different types of data: first, using random matrices sampled from the planted dense $m \times n$ submtarix model and, second, real-world collaboration and communication networks. ## Random matrices We first generate a random matrix with noise obscuring the planted submatrix using the function ``plantedsubmatrix``. and then call the function ``densub`` to recover the planted submatrix. ```{r, eval=FALSE} # Initialize problem size and densities # You can play around with these parameters M=100 #number of rows of sampled matrix N=200 #number of columns of sampled matrix m=50 #number of rows of dense submatrix n=40 #number of columns of dense submatrix p=0.25 # noise density q=0.85 #in-group density #Make binary matrix with mn-submatrix random<-plantedsubmatrix(M=M, N=N,m=m,n=n,p=p,q=q) ``` After generating the structure `random` containing the random matrix with desired planted structure, we can visually represent the matrix and planted submatrix as two-tone images, where dark pixels correspond to nonzero entries, and light pixels correspond to zero entries, using the following commands. ```{r, eval=FALSE} # Plot sampled G and matrix representations. image(random$sampled_matrix, useRaster=TRUE, axes=FALSE, main = "Matrix A") image(random$dense_submatrix, useRaster=TRUE, axes=FALSE, main = "Matrix X0") image(random$disagreements, useRaster=TRUE, axes=FALSE, main = "Matrix Y0") ``` Tne vizualization of the randomly generated matrix $\mathbf{A}$ helps us to understand its structure. It is clear that $\mathbf{A}$ contains a dense $50 \times 40$ block (in the bottom left corner). ![Visual representation of randomly generated $\mathbf{A}$](Rplot.jpeg){width=400px} We can remove all noise and isolate an image of a rank-one matrix $\mathbf{X0}$ with $mn$ nonzero entries. ![Visual representation of dense submatrix](Rplot01.jpeg){width=400px} Then we vizualize matrix $\mathbf{Y0}$ to see the number of disagreements between original matrix $\mathbf{A}$ and $\mathbf{X0}$. ![Disagreement between $\mathbf{A}$ and $\mathbf{X_0}$](Rplot02.jpeg){width=400px} We call the ADMM solver and visualize the output using the following commands. ```{r, eval=FALSE} #Call ADMM solver admm<-densub(G=random$sampled_matrix, m=m, n=n, tau = 0.35, gamma = 6/(sqrt(m*n)*(q-p)), opt_tol = 1.0e-4,maxiter=500, quiet = TRUE) #Plot results image(admm$X, useRaster=TRUE, axes=FALSE, main = "Matrix X") image(admm$Y, useRaster=TRUE, axes=FALSE, main = "Matrix Y") ``` The ADMM solver returns the optimal solutions $\mathbf{X}$ and $\mathbf{Y}$. It must be noted that matrices $\mathbf X$ and $\mathbf Y$ are identical to the actual structures of $\mathbf{X_0}$ and $\mathbf{Y_0}$. The planted submatrix is recovered. ![Optimal solution \mathbf{X}](Rplot03.jpeg){width=400px} ![Optimal Solution \mathbf{Y}](Rplot04.jpeg){width=400px} ## Collaboration Network The following is a simple example on how one could use the package to analyze the collaboration network found in the JAZZ dataset. It is known that this network contains a cluster of $100$ musicians which performed together. ![JAZZ Network](0001.jpg){width=400px} We have already prepared dataset to work with. More details can be found in the provided file `JAZZ_IN_R.R.` ```{r jazz, eval=FALSE} #Load dataset load(file="JAZZ.RData") #Initialize problem size and densities G=new #define matrix G equivalent to JAZZ dataset m=100 #clique size or the number of rows of the dense submatrix n=100 #clique size of the number of columns of the dense sumbatrix tau=0.85 #regularization parameter opt_tol=1.0e-2 #optimal tolerance verbose=1 maxiter=2000 #number of iterations gamma=8/n #regularization parameter #call ADMM solver admm <- densub(G = G, m = m, n = n, tau = tau, gamma = gamma, opt_tol = opt_tol, maxiter=maxiter, quiet = TRUE) # Planted solution X0. X0=matrix(0L, nrow=198, ncol=198) #construct rank-one matrix X0 X0[1:100,1:100]=matrix(1L, nrow=100, ncol=100)#define dense block # Planted solution Y0. Y0=matrix(0L, nrow=198, ncol=198) #construct matrix for counting disagreements between G and X0 Y0[1:100,1:100]=matrix(1L,nrow=100,ncol=1000)-G[1:100,1:100] #Check primal and dual residuals. C=admm$X-X0 a=norm(C, "F") #Frobenius norm of matrix C b=norm(X0,"F") #Frobenius norm of matrix X0 recovery = matrix(0L,nrow=1, ncol=1)#create recovery matrix if (a/b^2