Bootstrap analiza stabilnosti klastera

🏠 dambamath.com

1 Učitavanje paketa

library(mclust)
library(tidyverse)
library(cluster)
library(factoextra)
library(furrr)
library(glue)
library(kableExtra)

2 Učitavanje podataka

Učitava se kvadratna matrica reda 369 koja predstavlja udaljenosti nizova stanja logova pri čemu je korištena optimal matching metrika.

mat2.om <- readRDS("MAT2_distOM.rds")

3 Definicije funkcija

3.1 Originalna implementacija (Tidyverse / Map)

Ova verzija koristi standardni tidyverse pristup s ugniježđenim mutate() i map() funkcijama unutar jednog tibble objekta.

• Prednosti: Kod je čitljiv, modularan i prati logiku “jedan redak podataka - jedna iteracija”. Odlična je za debugiranje i rad s manjim skupovima podataka.
• Mane: Izrazito spora na velikom broju iteracija (\(B > 1000\)). Zbog načina na koji R obrađuje ugniježđene map pozive unutar tibble-a, stvara se veliki memorijski overhead koji značajno usporava izvođenje.

calculate_bootstrap_stability <- function(dist_matrix, k_clusters, B_iterations = 200, seed = 123) {
  set.seed(seed)
  n <- nrow(dist_matrix)

  ref_pam   <- pam(dist_matrix, k = k_clusters, diss = TRUE)$clustering
  ref_agnes <- cutree(agnes(dist_matrix, diss = TRUE, method = "ward"), k = k_clusters)

  results <- tibble(iter = 1:B_iterations) %>%
    mutate(
      u_idx = map(iter, ~ unique(sample(n, replace = TRUE))),

      ari_agnes = map_dbl(u_idx, ~ {
        d_sub <- dist_matrix[.x, .x]
        cl_boot <- cutree(agnes(d_sub, diss = TRUE, method = "ward"), k = k_clusters)
        adjustedRandIndex(ref_agnes[.x], cl_boot)
      }),

      ari_pam = map_dbl(u_idx, ~ {
        d_sub <- dist_matrix[.x, .x]
        cl_boot <- pam(d_sub, k = k_clusters, diss = TRUE)$clustering
        adjustedRandIndex(ref_pam[.x], cl_boot)
      }),

      sil_agnes = map_dbl(u_idx, ~ {
        d_sub <- dist_matrix[.x, .x]
        cl_boot <- cutree(agnes(d_sub, diss = TRUE, method = "ward"), k = k_clusters)
        ss <- silhouette(cl_boot, d_sub, diss = TRUE)
        mean(ss[, "sil_width"])
      }),

      sil_pam = map_dbl(u_idx, ~ {
        d_sub <- dist_matrix[.x, .x]
        cl_boot <- pam(d_sub, k = k_clusters, diss = TRUE)$clustering
        ss <- silhouette(cl_boot, d_sub, diss = TRUE)
        mean(ss[, "sil_width"])
      })
    )

  return(results)
}

3.2 Fast implementacija (List / Map)

Optimizirana verzija koja izbjegava kompleksnu strukturu tibble objekta tijekom računanja. Funkcije se izvršavaju unutar obične liste, a tek se na kraju spajaju u tablicu pomoću bind_rows().

• Prednosti: Značajno brža od originalne verzije bez ikakvih promjena u hardverskom pristupu. Smanjuje broj internih operacija koje R mora obaviti pri svakom koraku iteracije.
• Mane: Dalje je ograničena na jednu procesorsku jezgru (single-threaded), što znači da snaga modernih višejezgrenih procesora ostaje neiskorištena.

calculate_bootstrap_stability_fast <- function(dist_matrix, k_clusters, B_iterations = 200, seed = 123) {
  set.seed(seed)
  n <- nrow(dist_matrix)

  ref_pam <- pam(dist_matrix, k = k_clusters, diss = TRUE)$clustering
  ref_agnes <- cutree(agnes(dist_matrix, diss = TRUE, method = "ward"), k = k_clusters)

  results_list <- map(1:B_iterations, function(i) {
    u_idx <- unique(sample(n, replace = TRUE))
    d_sub <- dist_matrix[u_idx, u_idx]

    cl_agnes <- cutree(agnes(d_sub, diss = TRUE, method = "ward"), k = k_clusters)
    ari_a <- adjustedRandIndex(ref_agnes[u_idx], cl_agnes)
    sil_a <- mean(silhouette(cl_agnes, d_sub, diss = TRUE)[, "sil_width"])

    cl_pam <- pam(d_sub, k = k_clusters, diss = TRUE)$clustering
    ari_p <- adjustedRandIndex(ref_pam[u_idx], cl_pam)
    sil_p <- mean(silhouette(cl_pam, d_sub, diss = TRUE)[, "sil_width"])

    list(ari_agnes = ari_a, ari_pam = ari_p, sil_agnes = sil_a, sil_pam = sil_p)
  })

  bind_rows(results_list)
}

3.3 Nuclear implementacija (Parallel / Future)

Vrhunac optimizacije koji koristi furrr paket i future backend za distribuciju zadataka na sve dostupne procesorske jezgre (u ovom slučaju 20 threadova na Intel Ultra Core 7).

• Prednosti: Drastično ubrzanje koje omogućuje izvođenje ekstremno velikog broja iteracija (\(B = 100\,000\)) u svega par minuta, što bi na standardan način trajalo satima.
• Mane: Zahtijeva više RAM memorije (iako 64GB to rješava bez problema) i pažljivo upravljanje paralelnim planom pomoću plan(multisession). Postoji mali inicijalni trošak (overhead) za pripremu radnih procesora pa je neisplativa za jako mali broj iteracija.

calculate_bootstrap_stability_nuclear <- function(dist_matrix, k_clusters, B_iterations = 200, seed = 123) {
  set.seed(seed)
  n <- nrow(dist_matrix)

  ref_pam <- pam(dist_matrix, k = k_clusters, diss = TRUE)$clustering
  ref_agnes <- cutree(agnes(dist_matrix, diss = TRUE, method = "ward"), k = k_clusters)

  results_list <- future_map(1:B_iterations, function(i) {
    u_idx <- unique(sample(n, replace = TRUE))
    d_sub <- dist_matrix[u_idx, u_idx]

    cl_agnes <- cutree(agnes(d_sub, diss = TRUE, method = "ward"), k = k_clusters)
    ari_a <- mclust::adjustedRandIndex(ref_agnes[u_idx], cl_agnes)
    sil_a <- mean(cluster::silhouette(cl_agnes, d_sub, diss = TRUE)[, "sil_width"])

    cl_pam <- cluster::pam(d_sub, k = k_clusters, diss = TRUE)$clustering
    ari_p <- mclust::adjustedRandIndex(ref_pam[u_idx], cl_pam)
    sil_p <- mean(cluster::silhouette(cl_pam, d_sub, diss = TRUE)[, "sil_width"])

    list(ari_agnes = ari_a, ari_pam = ari_p, sil_agnes = sil_a, sil_pam = sil_p)
  }, .options = furrr_options(seed = TRUE))

  bind_rows(results_list)
}

4 Benchmark funkcija

run_stability_benchmark <- function(dist_matrix, test_ks = c(3, 4), B = 1000) {  
  vremena_lista <- list()
  
  for(k in test_ks) {
    
    message(glue("\n{str_dup('=', 30)}"))
    message(glue("POKREĆEM BENCHMARK ZA k = {k} (B = {B})"))
    message(glue("{str_dup('=', 30)}"))
    
    # 1. ORIGINALNA FUNKCIJA
    message("--> Izvodi se: Originalna metoda...")
    t_orig <- system.time({
      calculate_bootstrap_stability(dist_matrix, k_clusters = k, B_iterations = B)
    })
    vremena_lista[[length(vremena_lista) + 1]] <- tibble(k = k, metoda = "Original", vrijeme = t_orig["elapsed"])
    
    # 2. FAST FUNKCIJA
    message("--> Izvodi se: Fast (Map) metoda...")
    t_fast <- system.time({
      calculate_bootstrap_stability_fast(dist_matrix, k_clusters = k, B_iterations = B)
    })
    vremena_lista[[length(vremena_lista) + 1]] <- tibble(k = k, metoda = "Fast", vrijeme = t_fast["elapsed"])
    
    # 3. NUCLEAR FUNKCIJA
    message(glue("--> Izvodi se: Nuclear (Parallel) na {parallel::detectCores() - 2} jezgre..."))
    plan(multisession, workers = parallel::detectCores() - 2)
    t_nuke <- system.time({
      calculate_bootstrap_stability_nuclear(dist_matrix, k_clusters = k, B_iterations = B)
    })
    plan(sequential) 
    vremena_lista[[length(vremena_lista) + 1]] <- tibble(k = k, metoda = "Nuclear", vrijeme = t_nuke["elapsed"])
  }
  
  # Spajanje rezultata i izračun ubrzanja
  res_df <- bind_rows(vremena_lista) %>%
    group_by(k) %>%
    mutate(speedup = first(vrijeme) / vrijeme) %>%
    ungroup()
  
  return(res_df)
}

5 Pokretanje benchmarka

Testiranje svih implementacija na \(1000\) bootstrap uzoraka. U ovom slučaju je testirano za \(k=3\) i \(k=4\) pri čemu je \(k\) zadani broj klastera. Uočite da svaka implementacija na svakom uzorku radi hijerarhijsko i PAM klasteriranje i za svako takvo klasteriranje određuje silhouette svih elemenata u uzorku i prosječnu vrijednost i računa adjusted random index za svako takvo klasteriranje u odnosu na klasteriranje na početnom skupu podataka.

moj_benchmark <- run_stability_benchmark(mat2.om, test_ks = c(3,4), B = 1000)
moj_benchmark %>% kbl(format = "html") %>%
  kable_styling(
    bootstrap_options = c("hover", "striped", "condensed"), 
    full_width = FALSE
  )

k	metoda	vrijeme	speedup
3	Original	20.245	1.000000
3	Fast	10.678	1.895954
3	Nuclear	2.144	9.442631
4	Original	22.741	1.000000
4	Fast	11.939	1.904766
4	Nuclear	2.279	9.978499

6 Vizualizacija

# Filtriramo samo stupac 'vrijeme' za prikaz na grafu
moj_benchmark_long <- moj_benchmark %>%
  select(k, metoda, vrijeme) %>% 
  pivot_longer(cols = vrijeme, names_to = "tip_mjerenja", values_to = "vrijednost")

ggplot(moj_benchmark_long, aes(x = metoda, y = vrijednost, fill = metoda)) +
  geom_col(position = "dodge") +
  facet_wrap(~k) +
  labs(title = "Benchmark performansi", 
       x = "Metoda", 
       y = "Vrijeme (s)") +
  theme_minimal()

7 Nuclear stress test (100 000 iteracija)

U donjem videu možete vidjeti testiranje nuclear implementacije na \(100\,000\) bootstrap uzoraka. U videu se vidi dodatno GPU opterećenje, ali to je zbog snimanja ekrana. Također, zbog snimanja ekrana, izvršavanje je trajalo malo više, oko 256 sekundi. U normalnim okolnostima izvršavanje je trajalo 193 sekunde. U svakom slučaju impresivna brzina pri čemu se istovremeno računalo može normalno koristiti za neke druge zadatke. Osim dobrog hardvera, zaslužan je i cachyOS.

Status: Nuclear benchmark uspješno završen.
Živio Linux! 🐧 ❄️ (CachyOS)