Adaptative method to detect convergence
L.-C. Canon
June 27, 2017
Let’s design a method that stop the number of iterations when a specific criterion is reached and provide an estimation of when the convergence is reached.
test_convergence <- function(values) {
n <- length(values)
if (n == 1)
return(0)
first_half <- head(values, n / 2)
second_half <- tail(values, n / 2)
sum(first_half >= min(second_half) & first_half <= max(second_half)) / (n / 2)
}
termination_criterion <- function(values) {
test_convergence(values) >= 0.9
}
convergence_func <- function(values) {
max(which(map_dbl(seq_along(values), ~ test_convergence(values[1:.])) <= 0.5)) / 2
}
detect_convergence <- function(object, modif_func, measures_func,
termination_criterion, convergence_func,
min_iteration = 0, max_iteration = 1e6) {
require(purrr)
values <- matrix(NA, nrow = max_iteration, ncol = length(measures_func))
iteration <- 1
conv <- rep(NA, length(measures_func))
max_iter <- conv
while (iteration <= max_iteration) {
object <- modif_func(object)
values[iteration,] <- map_dbl(measures_func, ~.(object))
stopifnot(all(!is.na(values[1:iteration,])))
term <- apply(matrix(values[1:iteration,], nrow = iteration),
2, termination_criterion)
while (any(term) && iteration >= min_iteration) {
meas <- which(term)[1]
res_idx <- which(map_int(seq_along(conv), ~ sum(is.na(conv[1:.]))) == meas)[1]
conv[res_idx] <- convergence_func(values[1:iteration,meas])
max_iter[res_idx] <- iteration
new_meas <- setdiff(seq_along(measures_func), meas)
values <- matrix(values[,new_meas], ncol = length(new_meas))
measures_func <- measures_func[new_meas]
term <- term[new_meas]
}
if (length(measures_func) == 0)
break
iteration <- iteration + 1
}
data.frame(conv = conv, max_iter = max_iter)
}Let’s test the method with some settings:
n <- 30
m <- 10
M <- 10
col <- sample(1:M, n, replace = TRUE)
row <- sample(1:M, m, replace = TRUE)
mat <- matrix(col) %*% t(matrix(row))
modif_mat <- function(mat) {
rows <- sample(nrow(mat), 2)
cols <- sample(ncol(mat), 2)
submat <- mat[rows,cols]
shift_min <- -min(submat[1,1], submat[2,2])
shift_max <- min(submat[1,2], submat[2,1])
shift <- sample(shift_min:shift_max, size = 1)
submat[1,1] <- submat[1,1] + shift
submat[2,2] <- submat[2,2] + shift
submat[1,2] <- submat[1,2] - shift
submat[2,1] <- submat[2,1] - shift
mat[rows,cols] <- submat
stopifnot(all(mat >= 0))
mat
}
measures <- list(sd_meas, chi2_meas, sd_meas, chi2_meas)
detect_convergence(mat, modif_mat, measures, termination_criterion,
convergence_func, min_iteration = n * m, max_iteration = 3e4)## conv max_iter
## 1 3678.5 20849
## 2 3680.0 18210
## 3 3678.5 20849
## 4 3680.0 18210This is consistent with previous results. Let’s repeat the test with more measures several times to see the robustness:
measures <- list(sd_meas, mean_CV_row, mean_CV_col,
mean_cor_row, mean_cor_col, chi2_meas)
detect_convergence(mat, modif_mat, measures, termination_criterion,
convergence_func, min_iteration = n * m, max_iteration = 3e4)## conv max_iter
## 1 1672.5 10260
## 2 2212.0 9140
## 3 2442.5 13023
## 4 862.5 6788
## 5 1098.5 8320
## 6 2130.5 10000detect_convergence(mat, modif_mat, measures, termination_criterion,
convergence_func, min_iteration = n * m, max_iteration = 3e4)## conv max_iter
## 1 3058.5 12807
## 2 2292.5 12565
## 3 4653.0 12563
## 4 1550.0 11620
## 5 2446.0 12836
## 6 2975.0 12661detect_convergence(mat, modif_mat, measures, termination_criterion,
convergence_func, min_iteration = n * m, max_iteration = 3e4)## conv max_iter
## 1 2841.0 14955
## 2 2752.0 7220
## 3 2795.0 22374
## 4 3358.0 11494
## 5 3147.0 20620
## 6 2803.5 14930There are some consistencies (convergence between 1000 and 4000 and analyzing 10000 to 20000 iterations is sufficient). Let’s gather more data with a smaller matrix.
n <- 5
m <- 5
M <- 10
col <- sample(1:M, n, replace = TRUE)
row <- sample(1:M, m, replace = TRUE)
mat <- matrix(col) %*% t(matrix(row))
results <- replicate(30, detect_convergence(mat, modif_mat, measures,
termination_criterion,
convergence_func,
min_iteration = n * m))Let’s plot the convergence values:
convs <- do.call(rbind, results[1,])
plot.ecdf(convs[,1])
for (i in 2:length(measures))
plot.ecdf(convs[,i], col = i, add = TRUE)
legend("bottomright", c("sd", "meanCVrow", "meanCVcol", "corrow", "corcol", "chi2"),
fill = seq_along(measures))
Most measures are consistent. The convergence occurs with some variability between 20 and 100 most of the times.
iters <- do.call(rbind, results[2,])
plot.ecdf(iters[,1])
for (i in 2:length(measures))
plot.ecdf(iters[,i], col = i, add = TRUE)
legend("bottomright", c("sd", "meanCVrow", "meanCVcol", "corrow", "corcol", "chi2"),
fill = seq_along(measures))
400 iterations is often what is necessary to identify the convergence value.
Let’s see what happen with a 10 by 5 matrix:
n <- 10
m <- 5
M <- 10
col <- sample(1:M, n, replace = TRUE)
row <- sample(1:M, m, replace = TRUE)
mat <- matrix(col) %*% t(matrix(row))
results <- replicate(30, detect_convergence(mat, modif_mat, measures,
termination_criterion,
convergence_func,
min_iteration = n * m))par(mfrow = c(1, 2))
convs <- do.call(rbind, results[1,])
plot.ecdf(convs[,1])
for (i in 2:length(measures))
plot.ecdf(convs[,i], col = i, add = TRUE)
iters <- do.call(rbind, results[2,])
plot.ecdf(iters[,1])
for (i in 2:length(measures))
plot.ecdf(iters[,i], col = i, add = TRUE)
legend("bottomright", c("sd", "meanCVrow", "meanCVcol", "corrow", "corcol", "chi2"),
fill = seq_along(measures))
We observe a time to converge of 300 with some variation between the measures for a total number of necessary iterations of 2000. It is around 4 times more than with a 5 by 5 matrix and 10 times less than with a 30 by 10 matrix. Let’s increase further with a 10 by 10 matrix:
n <- 10
m <- 10
M <- 10
col <- sample(1:M, n, replace = TRUE)
row <- sample(1:M, m, replace = TRUE)
mat <- matrix(col) %*% t(matrix(row))
results <- replicate(30, detect_convergence(mat, modif_mat, measures,
termination_criterion,
convergence_func,
min_iteration = n * m))par(mfrow = c(1, 2))
convs <- do.call(rbind, results[1,])
plot.ecdf(convs[,1])
for (i in 2:length(measures))
plot.ecdf(convs[,i], col = i, add = TRUE)
iters <- do.call(rbind, results[2,])
plot.ecdf(iters[,1])
for (i in 2:length(measures))
plot.ecdf(iters[,i], col = i, add = TRUE)
legend("bottomright", c("sd", "meanCVrow", "meanCVcol", "corrow", "corcol", "chi2"),
fill = seq_along(measures))
Convergence around 800 and max iterations around 5000.
The next step is to design a general method that repeats the markov process 30 times from the same initial matrix and select the 0.8-quantile.
## R version 3.4.0 (2017-04-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/openblas-base/libblas.so.3
## LAPACK: /usr/lib/libopenblasp-r0.2.18.so
##
## locale:
## [1] LC_CTYPE=fr_FR.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=fr_FR.UTF-8 LC_COLLATE=fr_FR.UTF-8
## [5] LC_MONETARY=fr_FR.UTF-8 LC_MESSAGES=fr_FR.UTF-8
## [7] LC_PAPER=fr_FR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] stringr_1.0.0 tidyr_0.2.0 ggplot2_2.2.1 dplyr_0.4.3 purrr_0.2.0
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.10 knitr_1.15.1 magrittr_1.5 munsell_0.4
## [5] colorspace_1.2-2 R6_2.1.0 rlang_0.1 plyr_1.8.1
## [9] tools_3.4.0 parallel_3.4.0 grid_3.4.0 gtable_0.1.2
## [13] pacman_0.4.6 DBI_0.3.1 htmltools_0.3.6 yaml_2.1.13
## [17] lazyeval_0.2.0 rprojroot_1.2 digest_0.6.3 assertthat_0.1
## [21] tibble_1.3.1 evaluate_0.10 rmarkdown_1.5 stringi_0.5-5
## [25] compiler_3.4.0 scales_0.4.1 backports_1.0.5