enh: add gle-mm, gle-admm, and cgl implementations

Author: mirca

NAMESPACE

@@ -4,6 +4,7 @@ export(A)
 export(L)
 export(block_diag)
 export(cluster_k_component_graph)
+export(fscore)
 export(learn_bipartite_graph)
 export(learn_bipartite_k_component_graph)
 export(learn_combinatorial_graph_laplacian)

R/combinatorialGraphLaplacian.R

@@ -22,9 +22,12 @@
 #'             Selected Topics in Signal Processing, vol. 11, no. 6, pp. 825-841, Sept. 2017.
 #'             Original MATLAB source code is available at: https://github.com/STAC-USC/Graph_Learning
 #' @export
-learn_combinatorial_graph_laplacian <- function(S, A_mask, alpha, prob_tol = 1e-4,
-                                                max_cycle = 1000, regtype = 1, verbose = TRUE) {
+learn_combinatorial_graph_laplacian <- function(S, A_mask = NULL, alpha = 0, prob_tol = 1e-5,
+                                                max_cycle = 10000, regtype = 1,
+                                                record_objective = FALSE, verbose = TRUE) {
   n <- nrow(S)
+  if (is.null(A_mask))
+    A_mask <- matrix(1, n, n) - diag(n)
   e_v <- rep(1, n) / sqrt(n)
   dc_var <- t(e_v) %*% S %*% e_v
   isshifting <- c(abs(dc_var) < prob_tol)
@@ -48,6 +51,8 @@ learn_combinatorial_graph_laplacian <- function(S, A_mask, alpha, prob_tol = 1e-
     pb <- progress::progress_bar$new(format = "<:bar> :current/:total  eta: :eta",
                                      total = max_cycle, clear = FALSE, width = 80)
   time_seq <- c(0)
+  if (record_objective)
+    fun <- vanilla.objective(O_best - (1/n), K)
   start_time <- proc.time()[3]
   for (i in c(1:max_cycle)) {
     O_old <- O
@@ -99,6 +104,8 @@ learn_combinatorial_graph_laplacian <- function(S, A_mask, alpha, prob_tol = 1e-
     }
     O_best <- O
     C_best <- C
+    if (record_objective)
+      fun <- c(fun, vanilla.objective(O_best - (1/n), K))
     if (verbose)
       pb$tick()
     time_seq <- c(time_seq, proc.time()[3] - start_time)
@@ -113,7 +120,10 @@ learn_combinatorial_graph_laplacian <- function(S, A_mask, alpha, prob_tol = 1e-
   }
   O <- O_best - (1 / n)
   C <- C_best - (1 / n)
-  return(list(Laplacian = O, frob_norm = frob_norm, elapsed_time = time_seq))
+  results <- list(Laplacian = O, frob_norm = frob_norm, elapsed_time = time_seq)
+  if (record_objective)
+    results$obj_fun <- fun
+  return(results)
 }
 
 update_sherman_morrison_diag <- function(O, C, shift, idx) {

R/graphLaplacianEstimation.R

@@ -40,27 +40,28 @@ get_incidence_from_adjacency <- function(A) {
 #'             IEEE Trans. on Signal Processing, vol. 67, no. 16, pp. 4231-4244, Aug. 2019
 
 #' @export
-learn_laplacian_gle_mm <- function(S, A, alpha = 0, maxiter = 10000, reltol = 1e-5,
-                                   abstol = 1e-5, record_objective = FALSE,
-                                   verbose = TRUE) {
+learn_laplacian_gle_mm <- function(S, A_mask = NULL, alpha = 0, maxiter = 10000, reltol = 1e-5,
+                                   record_objective = FALSE, verbose = TRUE) {
+  # number of nodes
+  p <- nrow(S)
   Sinv <- MASS::ginv(S)
-  mask <- Ainv(A) > 0
+  if (is.null(A_mask))
+    A_mask <- matrix(1, p, p) - diag(p)
+  mask <- Ainv(A_mask) > 0
   w <- w_init("naive", Sinv)[mask]
   wk <- w
-  # number of nodes
-  p <- nrow(S)
   # number of nonzero edges
-  m <- .5 * sum(A > 0)
+  m <- sum(mask)#.5 * sum(A_mask > 0)
   # l1-norm penalty factor
   J <- matrix(1, p, p) / p
   H <- 2 * diag(p) - p * J
   K <- S + alpha * H
-  E <- get_incidence_from_adjacency(A)
+  E <- get_incidence_from_adjacency(A_mask)
   R <- t(E) %*% K %*% E
   r <- nrow(R)
   G <- cbind(E, rep(1, p))
   if (record_objective)
-    fun <- obj_func(E, K, wk, J)
+    fun <- vanilla.objective(L(wk), K)
   if (verbose)
     pb <- progress::progress_bar$new(format = "<:bar> :current/:total  eta: :eta",
                                      total = maxiter, clear = FALSE, width = 80)
@@ -72,11 +73,10 @@ learn_laplacian_gle_mm <- function(S, A, alpha = 0, maxiter = 10000, reltol = 1e
     Q <- Q[1:m, 1:m]
     wk <- sqrt(diag(Q) / diag(R))
     if (record_objective)
-      fun <- c(fun, obj_func(E, K, wk, J))
+      fun <- c(fun, vanilla.objective(L(wk), K))
     if (verbose)
        pb$tick()
-    werr <- abs(w - wk)
-    has_converged <- all(werr <= .5 * reltol * (w + wk)) || all(werr <= abstol)
+    has_converged <- norm(w - wk, "2") / norm(w, "2") < reltol
     if (has_converged && k > 1) break
     w <- wk
   }
@@ -89,6 +89,7 @@ learn_laplacian_gle_mm <- function(S, A, alpha = 0, maxiter = 10000, reltol = 1e
   return(results)
 }
 
+
 obj_func <- function(E, K, w, J) {
   p <- ncol(J)
   EWEt <- E %*% diag(w) %*% t(E)
@@ -117,17 +118,20 @@ obj_func <- function(E, K, w, J) {
 #'             Optimization Algorithms for Graph Laplacian Estimation via ADMM and MM.
 #'             IEEE Trans. on Signal Processing, vol. 67, no. 16, pp. 4231-4244, Aug. 2019
 #' @export
-learn_laplacian_gle_admm <- function(S, A, alpha = 0, rho = 1, maxiter = 10000,
-                                     reltol = 1e-5, record_objective = FALSE,
-                                     verbose = TRUE) {
+learn_laplacian_gle_admm <- function(S, A_mask = NULL, alpha = 0, rho = 1, maxiter = 10000,
+                                     reltol = 1e-5, record_objective = FALSE, verbose = TRUE) {
   p <- nrow(S)
+  if (is.null(A_mask))
+    A_mask <- matrix(1, p, p) - diag(p)
   Sinv <- MASS::ginv(S)
   w <- w_init("naive", Sinv)
   Yk <- L(w)
   Theta <- Yk
-  P <- eigvec_sym(Yk)
-  P <- P[, 2:p]
   Ck <- Yk
+  C <- Ck
+  # ADMM constants
+  mu <- 2
+  tau <- 2
   # l1-norm penalty factor
   J <- matrix(1, p, p) / p
   H <- 2 * diag(p) - p * J
@@ -136,23 +140,35 @@ learn_laplacian_gle_admm <- function(S, A, alpha = 0, rho = 1, maxiter = 10000,
     pb <- progress::progress_bar$new(format = "<:bar> :current/:total  eta: :eta",
                                      total = maxiter, clear = FALSE, width = 80)
   if (record_objective)
-    fun <- c()
+    fun <- c(vanilla.objective(Theta, K))
   # ADMM loop
+  # P <- qr.Q(qr(rep(1, p)), complete=TRUE)[, 2:p]
+  P <- eigvec_sym(Yk)
+  P <- P[, 2:p]
   for (k in c(1:maxiter)) {
     Gamma <- t(P) %*% (K + Yk - rho * Ck) %*% P
     U <- eigvec_sym(Gamma)
     lambda <- eigval_sym(Gamma)
     d <- .5 * c(sqrt(lambda ^ 2 + 4 / rho) - lambda)
     Xik <- crossprod(sqrt(d) * t(U))
     Thetak <- P %*% Xik %*% t(P)
-    Ck <- (diag(pmax(0, diag(Yk / rho + Thetak))) +
-           A * pmin(0, Yk / rho + Thetak))
-    Yk <- Yk + rho * (Thetak - Ck)
+    Ck_tmp <- Yk / rho + Thetak
+    Ck <- (diag(pmax(0, diag(Ck_tmp))) +
+           A_mask * pmin(0, Ck_tmp))
+    Rk <- Thetak - Ck
+    Yk <- Yk + rho * Rk
     if (record_objective)
-      fun <- c(fun, aug_lag(K, P, Xik, Yk, Ck, d, rho))
-    has_converged <- norm(Theta - Thetak, "F") / norm(Thetak, "F") < reltol
+      fun <- c(fun, vanilla.objective(Thetak, K))
+    normF_Rk <- norm(Rk, "F")
+    has_converged <-  norm(Theta - Thetak) / norm(Theta, "F") < reltol
     if (has_converged && k > 1) break
+    #s <- rho * norm(C - Ck, "F")
+    #if (normF_Rk > mu * s)
+    #  rho <- rho * tau
+    #else if (s > mu * normF_Rk)
+    #  rho <- rho / tau
     Theta <- Thetak
+    C <- Ck
     if (verbose)
       pb$tick()
   }

R/learnGraphTopology.R

@@ -560,6 +560,7 @@ learn_bipartite_graph <- function(S, is_data_matrix = FALSE, z = 0, nu = 1e4, al
 #'      vertex.color = c("red","black")[V(estimated_bipartite)$type + 1],
 #'      vertex.shape = c("square", "circle")[V(estimated_bipartite)$type + 1],
 #'      vertex.label = NA, vertex.size = 5)
+
 #' @export
 learn_bipartite_k_component_graph <- function(S, is_data_matrix = FALSE, z = 0, k = 1,
                                               w0 = "naive", m = 7, alpha = 0., beta = 1e4,

R/objectiveFunction.R

@@ -46,3 +46,8 @@ joint.prior <- function(beta, nu, Lw, Aw, U, V, lambda, psi) {
   return(laplacian.prior(beta = beta, Lw = Lw, lambda = lambda, U = U) +
          bipartite.prior(nu = nu, Aw = Aw, psi = psi, V = V))
 }
+
+vanilla.objective <- function(Theta, K) {
+  p <- nrow(Theta)
+  return(sum(diag(Theta %*% K)) - sum(log(eigval_sym(Theta)[2:p])))
+}

R/utils.R

@@ -28,6 +28,22 @@ relative_error <- function(A, B) {
 }
 
 
+#' Computes the fscore between two matrices
+#'
+#' @param A first matrix
+#' @param B second matrix
+#' @param eps real number such that edges whose values are smaller than eps are
+#'            not considered in the computation of the fscore
+#' @examples
+#' library(spectralGraphTopology)
+#' X <- L(c(1, 0, 1))
+#' fscore(X, X)
+#' @export
+fscore <- function(A, B, eps = 1e-4) {
+  return(metrics(A, B, eps)[1])
+}
+
+
 # Compute the prial value between two matrices
 # @param Ltrue true Laplacian matrix
 # @param Lest estimated Laplacian matrix

benchmarks/cospectral/example.R

@@ -8,7 +8,7 @@ set.seed(123)
 p <- 15
 f <- .05
 w <- runif(p * (p - 1) / 2)
-Laplacian <- L(w)
+Laplacian <- block_diag(L(w), L(w), L(w))
 lambda <- spectralGraphTopology:::eigval_sym(Laplacian)
 Adjacency <- diag(diag(Laplacian)) - Laplacian
 # construct the network

papers/1909.11594.pdf

@@ -17,11 +17,12 @@ test_that("learn_k_component_graph with single component random graph", {
 
 
 with_parameters_test_that("we can recover a simple connected graph with the GLE-MM and GLE-ADMM methods", {
-    w <- sample(1:10, 6)
+    w <- sample(0:10, 6)
     Laplacian <- L(w)
     n <- ncol(Laplacian)
-    Y <- MASS::mvrnorm(n * 500, rep(0, n), MASS::ginv(Laplacian))
-    res <- func(cov(Y), A = A(rep(1, length(w))), record_objective = TRUE)
+    A_mask <- 1 * (Laplacian < 0)
+    Y <- MASS::mvrnorm(n * 250, rep(0, n), MASS::ginv(Laplacian))
+    res <- func(cov(Y), A_mask = A_mask, record_objective = TRUE, reltol = 1e-7)
     expect_true(res$convergence)
     expect_true(relative_error(Laplacian, res$Laplacian) < 1e-1)
     expect_true(metrics(Laplacian, res$Laplacian, 1e-1)[1] > .9)