Short description

In every iteration \(i\), an algorithm is in a permutation, say, \(\sigma_i\). Then a random transposition is drawn uniformly \(t_i = (j,k)\) and the value of \(f(\sigma_i \circ t_i)\) is computed.

• If new value is bigger than the previous one, (i.e. \(f(\sigma_i \circ t_i) \ge f(\sigma_i)\)), then we set \(\sigma_{i+1} = \sigma_i \circ t_i\).
• If new value is smaller (\(f(\sigma_i \circ t_i) < f(\sigma_i)\)), then we will choose \(\sigma_{i+1} = \sigma_i \circ t_i\) with probability \(\frac{f(\sigma_i \circ t_i)}{f(\sigma_i)}\). Otherwise, we set \(\sigma_{i+1} = \sigma_i\) with complementary probability \(1 - \frac{f(\sigma_i \circ t_i)}{f(\sigma_i)}\).

The final value is the best \(\sigma\) that was ever computed.

Notes

This algorithm was tested in multiple settings and turned out to be an outstanding optimizer. Especially given it does not need any hyperparameters tuned.

The only parameter it depends on is max_iter, which determines the number of steps described above. One should choose this number rationally. When decided too small, there is a missed opportunity to find potentially a much better permutation. When decided too big, there is a lost time and computational power that does not lead to growth. Our recommendation is to plot the convergence plot with a logarithmic OX scale: plot(g_map, type = "both", logarithmic_x = TRUE). Then decide if the line has flattened already. Keep in mind that the OY scale is also logarithmic. For example, a small change on the OY scale could mean \(10000\) times the change in A Posteriori.

This algorithm has been analyzed extensively by statisticians. Thanks to the ergodic theorem, the frequency of visits of a given state converges almost surely to the probability of that state. This is the approach explained in the [1, Sec.4.1.2] and shown in [1, Sec. 5.2]. One can obtain estimates of posterior probabilities by setting return_probabilities = TRUE.

Example

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

Let’s say we have the data Z from the unknown process:

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 <- find_MAP(g, max_iter = 10, optimizer = "Metropolis_Hastings")
#> ========================================================================
#> Warning: The found permutation has n0 = 67 which is bigger than the number_of_observations = 50.
#> ℹ The covariance matrix invariant under the found permutation does not have the likelihood properly defined.
#> ℹ For more in-depth explanation, see 'Project Matrix - Equation (6)' section in `vignette('Theory', package = 'gips')` or its pkgdown page: https://przechoj.github.io/gips/articles/Theory.html.
g_map
#> The permutation (5,25)(9,45,60)(37,70)
#>  - was found after 10 log_posteriori calculations
#>  - is 15208559032166170453082112 times more likely than the starting, () permutation.

Only after ten iterations the found permutation is unimaginably more likely than the original, () permutation.

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

Short description

In every iteration \(i\), an algorithm is in a permutation; call it \(\sigma_i\). Then all the values of \(f(\sigma_i \circ t)\) are computed for every possible transposition \(t = (j,k)\). Then the next \(\sigma_{i+1}\) will be the one with the biggest value:

\[\sigma_{i+1} = argmax_{\text{perm} \in \text{neighbors}(\sigma_{i})}\{\text{posteriori}(perm)\}\]

Where: \[\text{neighbors}(\sigma) = \{\sigma \circ (j,k) : 1 \le j < k \le \text{p}\}\]

The algorithm ends when all neighbors are less likely, or the max_iter was achieved. In the first case, the algorithm will end in a local maximum, but there is no guarantee that this is also the global maximum.

Pseudocode

hill_climb <- function(g, max_iter) {
  perm <- g[[1]]
  perm_posteriori <- 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_posteriori <- -Inf
  
  i <- 1
  
  while (best_neighbor_posteriori > perm_posteriori && i < max_iter) {
    best_neighbor <- NULL
    best_neighbor_posteriori <- -Inf
    
    for (j in 1:(perm_size-1)) {
      for (k in (j + 1):perm_size) {
        neighbor <- gips:::compose_with_transposition(perm, c(j, k))
        neighbor_posteriori <- log_posteriori_of_gips(gips(
          S,
          number_of_observations,
          perm = neighbor
        ))
        
        if (neighbor_posteriori > best_neighbor_posteriori) {
          best_neighbor <- neighbor
          best_neighbor_posteriori <- neighbor_posteriori
        } # end if
      } # end for k
    } # end for j
    i <- i + 1
  } # end while
  
  return(best_neighbor)
}

Notes

The max_iter parameter works differently for this optimizer and Metropolis-Hastings. The Metropolis-Hastings will compute a posteriori of max_iter permutations. The Brute Force optimizer, \({p\choose 2} \cdot\) max_iter.

Example

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

Let’s say we have the data Z from the unknown process:

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")

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 matrix invariant under the found permutation does not have the likelihood properly defined.
#> ℹ For more in-depth explanation, see 'Project Matrix - Equation (6)' section in `vignette('Theory', package = 'gips')` or its pkgdown page: https://przechoj.github.io/gips/articles/Theory.html.
g_map
#> The permutation (2,16)(14,20)
#>  - was found after 601 log_posteriori calculations
#>  - is 4414155601613037260524331466752 times more likely than the starting, () permutation.
plot(g_map, type = "both")

Example

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 <- MASS::mvrnorm(number_of_observations, mu = mu, Sigma = sigma_matrix)

Let’s say we have the data Z from the unknown process:

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 after 720 log_posteriori calculations
#>  - is 80063.152438581 times more likely than the starting, () permutation.