Asymmetrical Overlap Matrix

[luedsch]

Good morning,

I am currently updating an old workflow for solvation free energy simulations with GROMACS that uses the old alchemical_analysis.py to alchemlyb. While testing alchemlyb (version 2.5.0) I realised that I get different overlap matrices, eventough I am using the same dhdl files from gmx energy. Using alchemical_analysis.py the overlap matrix is symmetric (O_MBAR_alchemical-analysis.pdf) , with alchemlyb on the other hand not (O_MBAR_alchemlyb.pdf). The code for alchemlyb is pasted below and a zipped folder containing the dhdl.*.xvg files can be accessd here: dhdl-files.zip.

from alchemlyb.parsing.gmx import extract_u_nk
from alchemlyb.preprocessing.subsampling import decorrelate_u_nk#
from alchemlyb.estimators import MBAR
from alchemlyb.visualisation import plot_mbar_overlap_matrix

xvg_files = sorted(glob.glob("dhdl-files/dhdl.*.xvg"))
u_nk_list = [extract_u_nk(outfile=f, T=298.15, filter=True) for f in xvg_files]
decorrelated_u_nk_list = [decorrelate_u_nk(u_nk) for u_nk in u_nk_list]
mbar = MBAR()
mbar.fit(alchemlyb.concat(decorrelated_u_nk_list)) 
mbar.delta_f_    # -> this output is still a symmetric matrix
ax = plot_mbar_overlap_matrix(mbar.overlap_matrix)

I do get some information about the statistical inefficiency when decorrelating the data, if that matters.

2026-02-05 08:45:49.800 | DEBUG    | alchemlyb.preprocessing.subsampling:statistical_inefficiency:557 - Statistical inefficiency: 14.78.
2026-02-05 08:45:49.801 | DEBUG    | alchemlyb.preprocessing.subsampling:statistical_inefficiency:563 - Number of uncorrelated samples: 1334.

Do you have any idea why the matrix computed with alchemlyb isn't symmetric? From my current understanding I would expect it to be. Would you revise my code snippet, eventough it is based on the tutorial on Read the Docs.

[xiki-tempula]

Why the Overlap Matrix is Not Symmetric

The overlap matrix in pymbar is asymmetric because of how it scales the raw overlap weights by the number of samples in each state.

1. The Code Implementation

In pymbar/mbar.py, the overlap matrix $O$ is calculated as:

608:        W = self.weights()
609:        O = self.N_k * (W.T @ W)

2. Mathematical Explanation

• Symmetric Part: The term W.T @ W computes the inner product of the weight matrix. This resulting matrix is mathematically symmetric: $(W^T W){ij} = (W^T W){ji}$

• Asymmetric Factor: The code multiplies this symmetric matrix by self.N_k (the array of sample counts). In NumPy, multiplying a 1D array (K,) with a 2D matrix (K, K) broadcasts across the last dimension (columns).

• Result:

$$ O_{ij} = (W^T W)_{ij} \times N_j $$

$$ O_{ji} = (W^T W)_{ji} \times N_i $$

Since $(W^T W){ij} = (W^T W){ji}$ but generally $N_i \neq N_j$, the final matrix elements $O_{ij}$ and $O_{ji}$ are not equal.

3. Physical Interpretation

The matrix is designed to be row-stochastic, meaning each row sums to 1.

• $O_{ij}$ represents the estimated probability that a sample generated from state $i$ would be observed in state $j$.
• To ensure the probabilities for a given source state $i$ sum to 1 ($\sum_j O_{ij} = 1$), the weights must be scaled by the number of samples in the destination state ($N_j$). This normalization makes the matrix useful for analyzing phase space overlap but breaks symmetry.