What are we optimizing?

The goal of the find_MAP() is to find the permutation \(\sigma\) that maximizes the a posteriori probability (MAP - Maximum A Posteriori). Such a permutation represents the most plausible symmetry given the data.

This a posteriori probability function is described in-depth in the Bayesian model selection section of the vignette("Theory", package="gips"), also available as a pkgdown page. gips can calculate the logarithm of it by the log_posteriori_of_gips() function. In the following paragraphs, we will refer to this a posteriori probability function as \(f(\sigma)\). We have \(f(\sigma) > 0\).

Available optimizers

The space of permutations is enormous - for the permutation of size \(p\), the space of all permutations is of size \(p!\) (\(p\) factorial). Even for \(p=12\), this space is practically impossible to browse. This is why find_MAP() implements multiple (3) optimizers to choose from:

• "brute_force", "BF", "full" | recommended for \(p\le 10\).
• "Metropolis_Hastings", "MH" | recommended for \(p\ge 11\).
• "hill_climbing", "HC"

library(gips)
# install.packages("pbmcapply")
library(pbmcapply)

get_gips_for_p <- function(p, n = 30) {
  S_inv <- diag(1, p, p)
  S_inv[1, 2:p] <- -1/sqrt(p)
  S_inv[2:p, 1] <- -1/sqrt(p)
  
  Sigma <- chol2inv(chol(S_inv))
  Z <- matrix(rnorm(n * p), nrow = n) %*% chol(Sigma)
  S <- cov(Z)
  S <- (S + t(S)) / 2 # Ensure exact symmetry for gips().
  
  gips(S, n)
}

get_time <- function(p, optimizer, max_iter = NA, seed = 1234) {
  set.seed(seed)
  g <- get_gips_for_p(p)
  g_MAP <- suppressWarnings(
    find_MAP(g, optimizer = optimizer, max_iter = max_iter, show_progress_bar = FALSE)
  )
  
  as.numeric(attr(g_MAP, "optimization_info")$optimization_time, units = "secs")
}

benchmark_optimization <- function(
    n_reps = 1, n_cores = 1, seed = 1234) {
  make_label <- function(optimizer, p, max_iter = NA) {
    if (is.na(max_iter)) {
      paste0(optimizer, " p=", p)
    } else {
      paste0(optimizer, " p=", p, ", ", max_iter)
    }
  }
  
  cases <- data.frame(
    p = c(2, 3, 4, 5, 6, 7, 8, 9, 10),
    max_iter = c(NA, NA, NA, NA, NA, NA, NA, NA, NA),
    stringsAsFactors = FALSE
  )
  
  cases$optimizer <- ifelse(is.na(cases$max_iter), "BF", "MH")
  
  cases$label <- mapply(
    make_label,
    optimizer = cases$optimizer,
    p = cases$p,
    max_iter = cases$max_iter
  )
  label_levels <- unique(cases$label)
  
  jobs <- cases[rep(seq_len(nrow(cases)), each = n_reps), ]
  jobs$rep <- rep(seq_len(n_reps), times = nrow(cases))
  jobs$seed <- seed + seq_len(nrow(jobs))
  
  rownames(jobs) <- NULL
  
  run_one_job <- function(j) {
    job <- jobs[j, ]

    # Warm up each worker on a small problem, outside the recorded timings.
    invisible(get_time(
      p = 6L,
      optimizer = job$optimizer,
      max_iter = if (is.na(job$max_iter)) NA else min(job$max_iter, 5L),
      seed = 1L
    ))

    # get_time() resets the seed before generating the measured sample.
    
    data.frame(
      label = job$label,
      p = job$p,
      optimizer = job$optimizer,
      max_iter = job$max_iter,
      rep = job$rep,
      seed = job$seed,
      time = get_time(
        p = job$p,
        optimizer = job$optimizer,
        max_iter = job$max_iter,
        seed = job$seed
      ),
      stringsAsFactors = FALSE
    )
  }
  
  results_list <- pbmcapply::pbmclapply(
    X = seq_len(nrow(jobs)),
    FUN = run_one_job,
    mc.cores = n_cores
  )
  
  results <- do.call(rbind, results_list)
  
  results$label <- factor(results$label, levels = label_levels)
  rownames(results) <- NULL
  
  results
}

start_time <- Sys.time()
# For 1.4.0, consider increasing n_reps while keeping runs sequential.
benchmark_results <- benchmark_optimization(n_reps = 1, n_cores = 1, seed = 1234)
end_time <- Sys.time()
end_time - start_time

aggregate(
  time ~ label,
  data = benchmark_results,
  FUN = function(x) {
    c(
      median = median(x),
      mean = mean(x),
      min = min(x),
      max = max(x),
      sd = sd(x)
    )
  }
)

perm_size <- 6
mu <- runif(perm_size, -10, 10) # Assume we don't know the mean
sigma_matrix <- matrix(
  data = c(
    1.0, 0.8, 0.6, 0.4, 0.6, 0.8,
    0.8, 1.0, 0.8, 0.6, 0.4, 0.6,
    0.6, 0.8, 1.0, 0.8, 0.6, 0.4,
    0.4, 0.6, 0.8, 1.0, 0.8, 0.6,
    0.6, 0.4, 0.6, 0.8, 1.0, 0.8,
    0.8, 0.6, 0.4, 0.6, 0.8, 1.0
  ),
  nrow = perm_size, byrow = TRUE
) # the real covariance matrix, that we want to estimate, is invariant under permutation (1,2,3,4,5,6)
number_of_observations <- 13
Z <- withr::with_seed(2022,
  code = MASS::mvrnorm(number_of_observations, mu = mu, Sigma = sigma_matrix)
)

dim(Z)
#> [1] 13  6
number_of_observations <- nrow(Z) # 13
perm_size <- ncol(Z) # 6
S <- cov(Z) # Assume we have to estimate the mean

g <- gips(S, number_of_observations)

g_map <- find_MAP(g, optimizer = "brute_force")
#> ================================================================================
g_map
#> The permutation (1,2,3,4,5,6):
#>  - was found in a search using 362 posterior calculations;
#>  - is 28979.967 times more likely than the () permutation.

perm_size <- 70
number_of_observations <- 50
Z <- withr::with_seed(2022, {
  mu <- runif(perm_size, -10, 10) # Assume we don't know the mean
  sigma_matrix <- crossprod(matrix(rnorm(perm_size * perm_size), nrow = perm_size))
  # sigma_matrix is the real covariance matrix, that we want to estimate based on Z
  MASS::mvrnorm(number_of_observations, mu = mu, Sigma = sigma_matrix)
})

dim(Z)
#> [1] 50 70
number_of_observations <- nrow(Z) # 50
perm_size <- ncol(Z) # 70
S <- cov(Z) # Assume we have to estimate the mean

g <- gips(S, number_of_observations)
suppressMessages( # message from ggplot2
  plot(g, type = "heatmap") +
    ggplot2::scale_x_continuous(breaks = c(1, 10, 20, 30, 40, 50, 60, 70)) +
    ggplot2::scale_y_reverse(breaks = c(1, 10, 20, 30, 40, 50, 60, 70))
)

g_map <- withr::with_seed(2023,
  code = find_MAP(g, max_iter = 150, optimizer = "Metropolis_Hastings")
)
#> ===============================================================================
g_map
#> The permutation (1,30,6,8,21,64,24,61,7,62,31,25,52,14,35,33,40,51,68,4,11,50,42,16,34,54,27,60,46,20,45,38,2,65,55,59,15,67,5,12,63,22,28,70,23)(3,37,66,18,56,36,57,29,58,9,17,10,41,43,53,69,39,48,47,44,19,32,26):
#>  - was found in a search using 150 posterior calculations;
#>  - is 5.915e+657 times more likely than the () permutation.

plot(g_map, type = "best", logarithmic_x = TRUE)

hill_climb <- function(g, max_iter) {
  perm <- g[[1]]
  perm_log_f <- log_posteriori_of_gips(g)
  perm_size <- attr(perm, "size")
  S <- attr(g, "S")
  number_of_observations <- attr(g, "number_of_observations")

  best_neighbor <- NULL
  best_neighbor_log_f <- -Inf

  i <- 1
  perm_i <- perm
  perm_i_log_f <- perm_log_f
  perm_i_minus_1_log_f <- -Inf

  while (i < max_iter) {
    best_neighbor <- NULL
    best_neighbor_log_f <- -Inf

    for (j in 1:(perm_size - 1)) {
      for (k in (j + 1):perm_size) {
        t <- c(j, k)
        neighbor <- gips:::compose_with_transposition(perm_i, t)
        neighbor_log_f <- log_posteriori_of_gips(gips(
          S, number_of_observations,
          perm = neighbor
        ))

        if (neighbor_log_f > best_neighbor_log_f) {
          best_neighbor <- neighbor
          best_neighbor_log_f <- neighbor_log_f
        } # end if
      } # end for k
    } # end for j
    i <- i + 1

    perm_i_minus_1_log_f <- perm_i_log_f
    
    if (best_neighbor_log_f <= perm_i_log_f) {
      break
    }

    perm_i <- best_neighbor
    perm_i_log_f <- best_neighbor_log_f
  } # end while

  return(perm_i)
}

perm_size <- 25
number_of_observations <- 20
Z <- withr::with_seed(2022, {
  mu <- runif(perm_size, -10, 10) # Assume we don't know the mean
  sigma_matrix <- crossprod(matrix(rnorm(perm_size * perm_size), nrow = perm_size))
  # sigma_matrix is the real covariance matrix, that we want to estimate based on Z
  MASS::mvrnorm(number_of_observations, mu = mu, Sigma = sigma_matrix)
})

dim(Z)
#> [1] 20 25
number_of_observations <- nrow(Z) # 20
perm_size <- ncol(Z) # 25
S <- cov(Z) # Assume we have to estimate the mean

g <- gips(S, number_of_observations)
plot(g, type = "heatmap")

# we will see a warning:
g_map <- find_MAP(g, max_iter = 2, optimizer = "hill_climbing")
#> ================================================================================
#> Warning: Hill Climbing algorithm did not converge in 2 iterations!
#> ℹ We recommend to run the `find_MAP(optimizer = 'continue')` on the acquired output.
#> Warning: The found permutation has n0 = 24, which is bigger than the number_of_observations = 20.
#> ℹ The covariance MLE is unavailable because the sample-size requirement is not met.
#> ℹ For a more in-depth explanation, see the 'Project Matrix - Equation (6)' section in the `vignette('Theory', package = 'gips')` or its pkgdown page: https://przechoj.github.io/gips/articles/Theory.html.
g_map
#> The permutation (1,13)(5,21):
#>  - was found in a search using 601 posterior calculations;
#>  - is 2.234e+21 times more likely than the () permutation.
plot(g_map, type = "best")

Continuing the optimization

When max_iter is reached during Metropolis-Hastings or hill climbing, the optimization stops and returns the result. Plot the best result found so far to assess whether it is still improving. A flat segment does not guarantee that the global maximum has been reached. Users can continue the optimization, as shown below.

# Generate data as in the Metropolis-Hastings example above,
# with a fixed seed for all random draws.

perm_size <- 70
number_of_observations <- 50
Z <- withr::with_seed(2022, {
  mu <- runif(perm_size, -10, 10) # Assume we don't know the mean
  sigma_matrix <- crossprod(matrix(rnorm(perm_size * perm_size), nrow = perm_size))
  # sigma_matrix is the real covariance matrix, that we want to estimate based on Z
  MASS::mvrnorm(number_of_observations, mu = mu, Sigma = sigma_matrix)
})

dim(Z)
#> [1] 50 70
number_of_observations <- nrow(Z) # 50
perm_size <- ncol(Z) # 70
S <- cov(Z) # Assume we have to estimate the mean

g <- gips(S, number_of_observations)

g_map <- withr::with_seed(2023,
  code = find_MAP(g, max_iter = 50, optimizer = "Metropolis_Hastings")
)
#> ==============================================================================
plot(g_map, type = "best")

We see the value is still increasing steadily. We now run 100 additional iterations and compare the best results before and after the continuation:

g_map2 <- withr::with_seed(2024,
  code = find_MAP(g_map, max_iter = 100, optimizer = "continue")
)
#> ===============================================================================
plot(g_map2, type = "best")

log_gain <- compare_log_posteriories_of_perms(
  g_map2, g_map, print_output = FALSE
)

The posterior probability of the best group found after 150 iterations is approximately \(10^{275}\) times that of the best group found after 50 iterations. Keep in mind that the y scale is logarithmic, so a visually small increase can represent a large posterior ratio. Use the plot to decide whether the best result is stabilizing and whether to run further iterations.

Additional parameters

The find_MAP() function has three additional logical parameters: show_progress_bar, save_all_perms, and return_probabilities.

When show_progress_bar = TRUE, gips prints “=” characters on the console during optimization. With return_probabilities = TRUE, a second progress bar shows probability calculation. For BF with \(p \le 10\), this bar tracks preparation of probability names from the canonical group representatives.

The save_all_perms = TRUE will save all visited permutations in the outputted object, which significantly increases the required RAM. For instance, with \(p=150\) and max_iter = 150000, we needed 400 MB to store it, whereas save_all_perms = FALSE only required 2 MB. However, save_all_perms = TRUE is necessary for return_probabilities = TRUE or more complex path analysis.

With return_probabilities = TRUE, probabilities are calculated from unnormalized posterior weights \(w(\Gamma)\), conditional on the set \(V\) of distinct visited cyclic groups:

\[P(\Gamma \mid Z,\Gamma\in V)=\frac{w(\Gamma)}{\sum_{H\in V}w(H)}.\]

Each group contributes one weight, regardless of repeated visits or different generators. BF uses all cyclic groups and gives the full posterior. A high conditional probability after MH does not establish that all important groups were found.

get_probabilities_from_gips() only reads the stored probabilities; saving permutations alone is not enough.

Probabilities after MH continuation or HC followed by MH use the combined saved visited groups. Keep save_all_perms = TRUE in every stage, and request probabilities in the final MH call:

g_small <- gips(diag(c(1, 1.4, 2)), 8, was_mean_estimated = FALSE)
g_hc <- find_MAP(g_small, optimizer = "HC", max_iter = Inf,
  save_all_perms = TRUE, show_progress_bar = FALSE)
set.seed(1)
g_mh <- find_MAP(g_hc, optimizer = "MH", max_iter = 10,
  save_all_perms = TRUE, return_probabilities = TRUE,
  show_progress_bar = FALSE)
get_probabilities_from_gips(g_mh)
#>     (1,2,3)       (2,3)          () 
#> 0.905590568 0.089759534 0.004649898

Standalone HC does not support return_probabilities = TRUE.

Discussion

We encourage everyone to discuss available and potential new optimizers on ISSUE#21. One can also see why some optimizers were implemented but not yet added to gips there.

References

[1] Piotr Graczyk, Hideyuki Ishi, Bartosz Kołodziejek, Hélène Massam. “Model selection in the space of Gaussian models invariant by symmetry.” The Annals of Statistics, 50(3) 1747-1774 June 2022. arXiv link; DOI: 10.1214/22-AOS2174

[2] “Learning Permutation Symmetry of a Gaussian Vector with gips in R” by gips developers Adam Chojecki, Paweł Morgen, and Bartosz Kołodziejek, Journal of Statistical Software.