q-Voter em networks

Junho-julho/2026, Mestrado, TF.

Implementação do modelo em redes com distribuição de graus $p(k)$ (Poisson, lei de potência, $\dots$).

Mean-field e q-Voter¶

O ponto de partida aqui, então, é o artigo:

📄 Mean-Field Analysis of the q-Voter Model on Networks
Paolo Moretti, Suyu Liu, Claudio Castellano, and Romualdo Pastor-Satorras
Journal of Statistical Physics 151, 113–130 (2013).

Neste artigo, os autores implementam a teoria de campo médio num grafo completo e com a teoria heterogênea, HMF. Nos dois casos, a teoria prevê duas fases, para (fase de fragmentação da opinião média) e ferromagnética (consenso da opinião média), bem como um região de biestabilidade, que, dependendo da condição inicial, leva ao consenso ou à fragmentação, como mostra a Fig Figure 1.

A curva que define a fronteira crítica entre a fase paramagnética e ferromagnética é dada por

onde $\phi$ é a magnetização (opinião) média da rede. Para a teoria HMF, sendo $\langle k\rangle$ o grau médio da rede, a fronteira é dada por

Algumas conclusões centrais deste artigo:

• As teorias de campo médio prevem biestabilidade para o q-voter para $q\geq4$, cuja fronteira é dada pelas equações (1) e (1);

• A taxa do modelo em redes heterogêneas fica definida para um nó de grau $k$ com $n$ vizinhos em um estado diferente por:

  $$f_k(n;q,\varepsilon) = \left(\frac{n}{k}\right)^q + \varepsilon\left[1 - \left(\frac{n}{k}\right)^q - \left(\frac{k-n}{k}\right)^q\right]\,;$$

  (3)

• Em RRNs com $K=4$, a biestabilidade não é observada, como previsto pelas teorias de campo médio, por conta da esparsidade da rede;

• Questões se RRNs muito densas e redes heterogêneas, como as em lei de potência, recuperam a biestabilidade permanecem abertas no artigo.

Nos interessa, aqui, o caso $q=4$, o menor valor de $q$ em que a biestabilidade é observada na teoria de campo médio. Procuramos fazer um diagrama de fases semelhante ao da Figura Figure 1, onde, para um dado valor $\varepsilon$, examinamos a dinâmica varrendo as condições iniciais da magnetização, $\phi_0\in[0,1]$.

Implementação computacional: benchmarks e resultados¶

Nas subseções que seguem, explico detalhadamente de que forma as redes complexas e a dinâmica (em tempo contínuo) foram implementadas. Também apresento benchmarks sobre a coerência e exatidão dos algoritmos, bem como as otimizações (lista dinâmica e paralelismo).

Redes com diferentes distribuições de graus¶

A implementação da rede utiliza de tipos derivados em Fortran. Em suma, uma rede corresponde a um objeto com algumas propriedades intrínsecas, a saber:

• adjacency_list (variável alocável): guarda o grau, $k_i$, e a lista de vizinhos de cada nó $i$ da rede (coleção de $N$ vetores alocáveis, cada um de tamanho $k_i$);

• gen_network (sub-rotina): gera a lista de graus da rede a partir de uma dada distribuição (power law, Poisson ou constante);

• build_network (sub-rotina): aplicação do modelo configuracional a partir de uma dada distribuição de graus e o tamanho da rede N;

• report (função): verifica os possíveis erros e inconsistências que podem aparecer com o modelo configuracional.

A grande vantagem desta abordagem é a capacidade de gerar redes totalmente independentes sem acúmulo de variáveis (ou uso de dummies), compacticidade do código e possibilidade de paralelização com OpenMP (basta passar a rede, e todas suas propriedades, como variável privada de uma thread), sem custo computacional extra.

Para gerar a distribuição de graus, dois algoritmos foram utilizados. No caso de uma distribuição power-law, o método mais interessante é a inversão de CDF com aproximação para dados discretos (acrescentar citação do Clauset); isto é, considere um número aleatório $u^{\text{r}}\in[0,1)$, então o número inteiro aleatório $x^{\text{r}}$, dado por:

é um número distribuído em lei de potência $p(x;\gamma)$ no intervalo $[x_{\min},\infty)$. Para controlar o surgimento de outliers e/ou o grau máximo da rede em lei de potência, pode-se atribuir um determinado cut-off: estrutural $\max(k) = \sqrt{N}$, natural $\max(k) = N^{\frac{1}{\gamma-1}}$ ou rígido $\max(k) = N^{\frac{1}{\gamma}}$

Para a de Poisson, o algoritmo de D. Knuth (ou sua variante, Junhao, acrescentar citações) é ideal. Ele parte da ideia que um número de Poisson conta o número de eventos que ocorrem numa janela temporal, em que o intervalo de tempo entre eventos sucessivos é dado por uma distribuição exponencial. Em outras palavras, é gerar números exponencialmente distribuídos com média $1/\lambda$ de tal forma que o número $N^{\text{r}}=m-1$, em que $m$ é o menor número tal que a sequência gerada é maior ou igual a um (isto é, $x_1^{\text{r}}+x_2^{\text{r}}+\cdots+x_m^{\text{r}}\geq 1$), é um número com distribuição de Poisson com média $\lambda$. Por razões de consistência, a distribuição gerada é truncada à esquerda, rejeitando graus inferiores a um dado $k_{\min}$.

O modelo configuracional segue o algoritmo abaixo:

• Dada a sequência de graus gerada, cria-se um vetor de “pontas soltas” (stubs) de tamanho $\sum k_i$ em que cada nó $i$ aparece (sequencialmente) $k_i$ vezes;

• Para evitar correlações na geração sequencial (seeds colision), o vetor de stubs é embaralhado com o Algoritmo P de Fisher-Yates (citar o livro do Knuth);

• O vetor é percorrido emparelhando stubs adjacentes (índices $j$ e $j+1$) para formar uma possível nova aresta entre os nós:

  • Caso a possível conexão for um auto-laço ou já existir, ela é rejeitada e o algoritmo continua:

    • Em caso de rejeição, um stub é sorteado (uniforme) dentre aqueles ainda não processados;

    • Caso falhe novamente, uma busca linear procura pelo primeiro stub válido disponível. Se nenhuma reatribuição for possível, o par é descartado;

  • Se a conexão for válida, a lista de adjacência de ambos os nós é atualizada, contabilizando a aresta final, e o processo avança até o esgotamento do vetor de stubs.

O método de sorteio e busca linear (two-stage cascade) serve para tentar ao máximo (para uma dada seed) encontrar conexões válidas para aquela distribuição. A Figure 2 mostra que a distribuição efetiva de graus (a distribuição de graus real da rede, após a aplicação do modelo configuracional) de fato segue as distribuições usadas para a amostragem dos graus, com boa convergência mesmo para redes com $N=10^4$ nós.

Finalmente, para verificar a autoconsistência das redes sintéticas e medir os desvios da lista de graus induzidos pelo método, a função de report faz uma grande checagem estrutural e estatística da topologia gerada. A verificação é dividida em duas etapas:

1. Autoconsistência Topológica:

   • Teorema do Grau: o número de arestas tem que ser igual à metade da soma dos graus efetivos dos nós;

   • Nós isolados: Contabiliza a existência de sítios com grau nulo ($k=0$);

   • Fragmentação da rede: Utilizando o algoritmo BFS, o código determina o tamanho da componente gigante (CG), de modo que, caso ela não equivalha à totalidade da da rede (sinalizando fragmentação), um alerta é emitido para que a componente principal seja isolada/extraída antes de iniciar a dinâmica.

2. Consistência Estatística (Teoria vs. Efetiva):

   • Arestas rejeitadas: Compara o número de arestas originais sorteadas com o número de arestas efetivamente conectadas, indicando a taxa de descarte exigida para manter a rede puramente simples;

   • Momentos da distribuição: Calcula o grau máximo ($k_{\max}$), o grau médio ($\langle k \rangle$) e o segundo momento matemático ($\langle k^2 \rangle$) da lista original e compara com a rede conectada. Uma tolerância estatística de $1/\sqrt{N}$ é adotada para classificar a confiabilidade da topologia gerada.

Estes relatórios on-the-fly servem para garantir robustez e validade ao algoritmo de redes, evitando que desvios induzidos pelo método enviesem os resultados.

Para várias redes geradas, nenhum desvio significativo foi encontrado. Assim, a robustez e validade deste algoritmo ficam comprovadas. A tabela abaixo reúne estes resultados comparativos para os casos de interesse.

Dinâmica em tempo contínuo¶

Para utilizar o algoritmo de Gillespie, o método de aceitação-rejeição (eventos reais e nulos) com lista dinâmica é a forma otimizada. A referência principal para implementação do algoritmo foi:

📄 Dynamic sampling from a discrete probability distribution with a known distribution of rates
Federico D’Ambrosio, Hans L. Bodlaender, and Gerard T. Barkema
Computational Statistics 37, 1203–1228 (2022).

A ideia é considerar apenas os nós que podem flipar. Basta definir uma lista que contém todos os $N_{\text{ativos}}$ nós ativos da rede (isto é, aqueles com $f_k>0$) e atualizá-la dinamicamente durante a evolução. Em suma:

• A um dado tempo $t$, sorteia-se aleatoriamente um nó $i$ exclusivamente dentre os $N_\text{ativos}$ sítios com taxa não-nula;

• Aceitação-rejeição padrão: dado um número aleatório uniforme $r^\text{r}\in[0,1)$, o spin flipa caso $f_{k_i}(n_i) \geq r^\text{r} \cdot \max_f$;

• Caso o estado seja alterado, a vizinhança local é atualizada. Sítios vizinhos que passam a ter 0 vizinhos discordantes são removidos da lista dinâmica (tornam-se inativos), enquanto sítios que adquirem discordância são adicionados à lista;

• O tempo físico do sistema avança independentemente do flip ter sido aceito ou rejeitado (evento nulo), considerando $u^\text{r}\in[0,1)$, por:

  $$t \leftarrow t - \frac{\ln(1-u^\text{r})}{N_\text{ativos} \cdot \max_f}\,.$$

  (5)

No pior dos casos, a lista dinâmica tem complexidade temporal igual ao método ingênuo de aceitação-rejeição (onde o sorteio do nó $i$ é feito sobre toda a rede). Mas em dinâmicas próximas ao estado de consenso, o tamanho da lista decresce rapidamente, de forma que o uso da lista dinâmica reduz a complexidade temporal. Outros métodos de otimização, a princípio, não geram melhorias no tempo computacional.

A Figura Figure 3 mostra a equivalência (estatística) das implementações em tempo discreto e tempo contínuo com lista dinâmica. As Figuras abaixo mostram o ganho computacional com o método em tempo contínuo.

Amostragem dos resultados¶

Para obtenção dos resultados, a condição de parada da simulação é dada por um tempo máximo numericamente igual ao tamanho da rede Gillespie_time = t_max=real(N) ou à obtenção do consenso, N_ativos = 0. Isto é feito para diferenciar claramente as fases de fragmentação e consenso. Na fase paramagnética o consenso é obtido em tempos muito grandes ($\tau\sim\exp(N)$), de modo que este “tempo máximo” serve para relevar o efeito de tamanho finito. De qualquer forma, a análise de tamanho finito ainda não foi realizada.

Para reproduzir a refbiestabilidade, foram 50 amostras (feitas com uma única rede) de cada ponto no espaço $(\phi_0,\varepsilon)$. Aqui, $\phi_0$ entra como condição inicial: $(\phi_0+1)/2$ é a proporção (probabilidade) de agentes no estado +1. Os diagramas abrangem 123 pontos distribuídos linearmente para $\varepsilon\in[0,0,0.3]$ e 101 para $\phi_0\in[-0.95,0.95]$ (totalizando 12423 pontos * 50 amostras = 621150 runs para cada rede).

Para adiantar a obtenção de amostras, a ferramenta OpenMP foi utilizada, utilizando somente núcleos de desempenho das máquinas do GISC. Para realizar o teste de benchmarking, 100 amostras (50 na fase paramagnética e 50 na fase ferromagnética) foram rodadas sequencialmente (em único núcleo) e com/sem hyper-threading na máquina Weiss. Como os processadores das máquinas Ising, Parisi, Fisher e Weiss seguem a mesma topologia (verificado com lstopo), o resultado percentual da tabela abaixo se mantém. O atraso do hyper-threading se explica porque a memória que a dinâmica ocupa é maior que a metade da memória disponível no cache L2: matriz de estados: 80 KB, matriz de taxas: 80 KB, lista de adjacência: 781 KB e lista dinâmica: 160KB, memória total de aproximadamente 1,1 MB, o que inviabiliza alocar duas threads simultâneas disputando o mesmo cache L2 de 2 MB sem causar cache misses.

Resultados preliminares¶