Paper deep dive
Interior interpretability with attention rollout: contraction and propagation profiles in Transformers
Umberto Biccari, Qian Huang, Enrique Zuazua
Intelligence
Status: not_run | Model: - | Prompt: - | Confidence: 0%
Entities (0)
Relation Signals (0)
No relation signals yet.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Feature-attribution methods assign scores relating input variables to a model's output, but do not by themselves characterize how explicitly defined interaction operators compose across its intermediate layers. We introduce \emph{interior interpretability}, a propagation-based perspective on internal model organization, and instantiate it for tabular Transformers using attention rollout. We interpret rollout as a row-stochastic operator encoding attention-mediated propagation between feature tokens. By applying classical Doeblin--Dobrushin contraction theory, we show that a rollout operator with a small Dobrushin coefficient is quantitatively close to a rank-one stochastic matrix whose common row is determined by its normalized column sums. This result gives a structural interpretation to the corresponding rollout propagation profile. In Transformers trained for metabolomic age prediction, the measured rollout contraction strengthens with depth. Trained and randomly initialized models also exhibit different propagation profiles, although the present experiments do not establish the predictive relevance of individual rollout-ranked variables. Exploratory comparisons with PCA and GradientExplainer approximations to SHAP reveal localized agreement among highly ranked variables but weak agreement across complete rankings. Attention rollout is therefore used here as a diagnostic of attention-mediated propagation, not as a causal explanation or faithful attribution of the complete Transformer.
Tags
Links
- Source: https://arxiv.org/abs/2607.22367v1
- Canonical: https://arxiv.org/abs/2607.22367v1
Trouble viewing inline? Open PDF directly →
Full Text
97,132 characters extracted from source content.
Expand or collapse full text
Interior interpretability with attention rollout: contraction and propagation profiles in Transformers Umberto Biccari † † Chair of Computational Mathematics, DeustoTech, University of Deusto, Avenida de las Universidades 24, 48007 Bilbao, Basque Country, Spain. umberto.biccari@deusto.es , Qian Huang ∗ ‡ ∗ Universidad Carlos I de Madrid, Gregorio Millán Institute for Fluid Dynamics, Nanoscience and Industrial Mathematics. ‡ Universidad Carlos I de Madrid, and Universidad Carlos I de Madrid, Department of Mathematics, Leganés, 28911, Spain. qhuang@inst.uc3m.es and Enrique Zuazua † ¶ § † Chair of Computational Mathematics, DeustoTech, University of Deusto, Avenida de las Universidades 24, 48007 Bilbao, Basque Country, Spain. ¶ Chair for Dynamics, Control, Machine Learning, and Numerics (Alexander von Humboldt-Professorship), Department of Mathematics, Friedrich-Alexander-Universität Erlangen-Nürnberg, 91058 Erlangen, Germany § Universidad Autónoma de Madrid, Departamento de Matemáticas, Ciudad Universitaria de Cantoblanco, 28049 Madrid, Spain. enrique.zuazua@deusto.es, enrique.zuazua@fau.de, enrique.zuazua@uam.es Abstract. Feature-attribution methods assign scores relating input variables to a model’s output, but do not by themselves characterize how explicitly defined interaction operators compose across its intermediate layers. We introduce interior interpretability, a propagation-based perspective on internal model organization, and instantiate it for tabular Transformers using attention rollout. We interpret rollout as a row-stochastic operator encoding attention-mediated propagation between feature tokens. By applying classical Doeblin–Dobrushin contraction theory, we show that a rollout operator with a small Dobrushin coefficient is quantitatively close to a rank-one stochastic matrix whose common row is determined by its normalized column sums. This result gives a structural interpretation to the corresponding rollout propagation profile. In Transformers trained for metabolomic age prediction, the measured rollout contraction strengthens with depth. Trained and randomly initialized models also exhibit different propagation profiles, although the present experiments do not establish the predictive relevance of individual rollout-ranked variables. Exploratory comparisons with PCA and GradientExplainer approximations to SHAP reveal localized agreement among highly ranked variables but weak agreement across complete rankings. Attention rollout is therefore used here as a diagnostic of attention-mediated propagation, not as a causal explanation or faithful attribution of the complete Transformer. Key words and phrases: Interior interpretability, Attention rollout, Transformer architectures, Propagation analysis 2020 Mathematics Subject Classification: 68T07,68T01,62R07 This project has received funding from the European Research Council (ERC) under the European Union’s Horizon Europe research and innovation programme (grant agreement No. 101096251, CoDeFeL). This material is based upon work supported by the Air Force Office of Scientific Research under award number FA8655-22-1-7012. EZ was partially supported by the Alexander von Humboldt Professorship program; the European Union’s Horizon Europe MSCA project ModConFlex (HORIZON-MSCA-2021-DN-01, project 101073558); the Transregio 154 Project Mathematical Modelling, Simulation and Optimization Using the Example of Gas Networks of the DFG; and SURE-AI: The Norwegian Centre for Sustainable, Risk-Averse, and Ethical AI, grant 357482, Research Council of Norway. UB and EZ were partially supported by the Grant PID2023-146872OB-I00-DyCMaMod of MICIU (Spain) and by the COST Actions CA24122-Multiscale Stochastics, Patterns, and Analysis of Combinatorial Environments and CA24136-Interactions between Control Theory and Machine Learning. QH was supported by the RYC2024-048848-I grant funded by MICIU/AEI/10.13039/501100011033 and the FSE+. 1. Introduction and motivation Methods for feature attribution address a central interpretability question: how do the input variables of a learned predictor contribute to its output? For a model f:ℝd→ℝf:R^d , methods such as SHAP, Integrated Gradients, and saliency maps analyze the map x↦f(x)x f(x) and assign a score to each coordinate of x [5, 27, 28, 29, 36, 37, 39]. These methods provide output-oriented explanations, but they do not by themselves describe how explicitly identifiable interaction operators accumulate across the intermediate layers of a deep architecture. This paper studies that complementary question. We use the term interior interpretability for the analysis of propagation operators associated with intermediate model interactions. The objective is deliberately narrower than reconstructing the complete hidden-state dynamics: a propagation operator captures only the part of the computation represented by its construction. In the Transformer setting considered here, the relevant operators are obtained from self-attention and describe attention-mediated exchange between feature tokens. They do not include value and output projections, normalization, MLP sublayers, or all nonlinear transformations of the complete model. Attention rollout [1] provides a natural object for this analysis. It composes the attention operators of successive layers, including the standard residual averaging, into a row-stochastic matrix. For tabular Transformers, where individual variables are represented by dedicated tokens, the resulting matrix describes how attention-mediated propagation between source and destination features accumulates with depth. This stochastic-operator viewpoint leads to two distinct questions: (1) What structural form does the rollout operator acquire as depth increases? (2) Which properties of its propagation profile follow from the rollout construction, and how do the observed profiles differ between trained and randomly initialized models? The first question is mathematical. Products of stochastic matrices are governed by classical contraction theory, suggesting that the rows of a sufficiently contractive rollout operator should approach a common profile. The rollout product supplies a layerwise contraction bound, but this bound does not by itself guarantee decay with depth. The second question is empirical. The common profile need not be uniform and depends on the realized attention matrices. Separating the mathematical constraints of rollout composition from profile differences observed after optimization is therefore essential: rank-one organization alone is not evidence of learning or task relevance. Figure 1.1 places this propagation viewpoint alongside two reference analyses used in the empirical study. PCA describes variance structure in the input data, GradientExplainer provides an expected-gradients approximation to prediction-attribution scores, and rollout describes the propagation geometry encoded by the attention operators. These are not three interchangeable notions of explanation; they answer different questions and are compared to determine where their rankings agree or diverge. Input data x=(x1,…,xd)x=(x_1,…,x_d) Internal representations Z(0)→⋯→Z(L)Z^(0)→·s→ Z^(L) Prediction f(x)f(x) PCA What structures variance in the data? Attention rollout How does attention-mediated propagation accumulate? GradientExplainer Which variables receive prediction-attribution scores? Figure 1.1. Three complementary analytical viewpoints used in this work. PCA describes variance structure in the input data, attention rollout describes attention-mediated propagation across intermediate layers, and GradientExplainer provides an expected-gradients approximation to prediction-attribution scores. These quantities are not interchangeable and need not induce the same feature rankings. The paper makes four contributions: (1) We formulate interior interpretability as the analysis of propagation operators associated with intermediate model interactions and instantiate this perspective using attention rollout in tabular Transformers. (2) Using a direct finite-dimensional consequence of classical Doeblin–Dobrushin theory [11], we characterize the rank-one regime of rollout operators and identify their normalized column sums with the corresponding mean row profile. (3) We place the layerwise contraction bound for rollout composition alongside empirical differences between trained and randomly initialized propagation profiles. (4) We compare rollout propagation rankings with PCA and GradientExplainer approximations to SHAP in a metabolomic age-prediction case study. These comparisons are descriptive and exploratory: rollout is treated as a diagnostic of attention-mediated propagation rather than as a causal or faithful explanation of the complete predictor. The empirical setting is metabolomic age prediction, an application of current interest in precision medicine [21, 23, 24]. The real-world study is complemented by a synthetic example used only as a second qualitative illustration. Across the trained architectures, rollout contraction strengthens with depth. The trained–random-initialization comparisons address a different question: the two groups exhibit different propagation profiles and rankings in the reported runs. Agreement with GradientExplainer scores is localized among some highly ranked variables, while complete rankings remain weakly correlated. The present paper is organized as follows. Section 2 introduces the analytical framework underlying the proposed interpretability approach. In particular, we present the attention rollout operator and its structural characterization. In Section 3, we describe the dataset and models used in our empirical study, whose results are presented in Section 4. We conclude with a discussion of the implications of our findings and potential directions for future research in Section 5. Appendix A contains the proof of the structural result and additional mathematical details supporting our analysis. Finally, Appendix B places selected highly ranked variables in biomedical context. 1.1. Positioning with respect to the existing literature The perspective developed here belongs to the broader effort to understand internal model mechanisms beyond input–output behavior. Related work includes mechanistic interpretability and neural-circuit analysis [14, 30, 31, 32], as well as methods that probe learned representations [2, 22]. These approaches study different objects and levels of explanation. Our narrower focus is the structure of an explicitly defined propagation operator derived from attention matrices. Within this broader perspective, attention rollout is interpreted as a global stochastic operator describing how attention-mediated interactions compose across successive Transformer layers. This operator does not represent the complete evolution of the hidden states. Our work is also related to the extensive literature discussing whether attention mechanisms can be interpreted as explanations. Previous studies [6, 20, 25, 34] have shown that attention distributions may exhibit weak correlation with attribution measures and that substantially different attention patterns can lead to similar predictions. More generally, several attention-based interpretability methods have been proposed for Transformers, including relevance-propagation approaches that combine attention with gradient information [10]. Consequently, attention weights should not generally be regarded as faithful feature attributions or causal explanations. The viewpoint adopted here is fully consistent with these observations. We do not interpret attention as explaining model predictions. Instead, rollout is viewed as an attention-derived propagation surrogate describing how attention mass is redistributed across tokens when layerwise attention matrices are composed. From this perspective, attribution-based methods and propagation-based analyses address complementary aspects of interpretability. 2. Attention-based interpretability Attention rollout defines feature-to-feature propagation operators summarizing attention-mediated exchange across Transformer layers. Understanding the structural properties of these operators allows us to describe the organization encoded by rollout and motivates the rollout propagation scores introduced later in the paper. 2.1. Transformers and the self-attention mechanism In this work, we consider a single-output Transformer-based regression model of the form x∈ℝd↦ξ(x)∈ℝp↦f(x)=w⊤ξ(x)+b∈ℝ, x ^d ξ(x) ^p f(x)=w ξ(x)+b , where p indicates the dimension of the latent representation space. At layer ℓ , the hidden state is Z(ℓ)∈ℝd×pZ^( ) ^d× p, whose i-th row Zi(ℓ)∈ℝpZ_i^( ) ^p corresponds to the latent representation of feature i. The latent representation ξ(x)∈ℝpξ(x) ^p is obtained by mean pooling over the feature representations of the last Transformer layer. The central mechanism governing feature interactions is multi-head self-attention [38]. For each layer ℓ∈⟦L⟧ ∈ L and attention head h∈⟦H⟧h∈ H , the attention matrix is A(ℓ,h)=softmax(Q(ℓ,h)(K(ℓ,h))⊤ph)∈ℝd×d, A^( ,h)=softmax ( Q^( ,h)(K^( ,h)) p_h ) ^d× d, where php_h denotes the query/key dimension of head h. The matrices are constructed from queries: Q(ℓ,h)=Z(ℓ−1)WQ(ℓ,h),keys: K(ℓ,h)=Z(ℓ−1)WK(ℓ,h). : Q^( ,h)=Z^( -1)W_Q^( ,h), : K^( ,h)=Z^( -1)W_K^( ,h). Here, the notation ⟦q⟧≔1,…,q for all q∈ℕ∗, q \1,…,q\ for all q , indicates the set containing the first q strictly positive natural numbers. Notice that A(ℓ,h)A^( ,h), as all other attention-based quantities considered in this work, is actually dependent on the input datum x∈ℝdx ^d, that is A(ℓ,h)=A(ℓ,h)(x). A^( ,h)=A^( ,h)(x). However, for simplicity of notation, we will avoid writing explicitly this dependence unless strictly necessary. Notice also that A(ℓ,h)A^( ,h) is row-stochastic as Aij(ℓ,h)≥0 for all i,j∈⟦d⟧ and ∑j=1dAij(ℓ,h)=1 for all i∈⟦d⟧. A^( ,h)_ij≥ 0\;\; for all i,j∈ d and _j=1^dA^( ,h)_ij=1\;\; for all i∈ d . Moreover, to describe the average propagation geometry induced by multi-head attention, we introduce the head-averaged attention matrix A(ℓ)=1H∑h=1HA(ℓ,h), A^( )= 1H _h=1^HA^( ,h), that by construction remains row-stochastic. Finally, motivated by the presence of residual connections, we shall view the hidden-state evolution schematically through the effective propagation operators I+A(ℓ)I+A^( ), so that successive Transformer layers induce an iterative process of the form Z(L)↝(I+A(L))⋯(I+A(1))Z(0). Z^(L) (I+A^(L) )·s (I+A^(1) )Z^(0). This effective description provides a convenient representation of how local feature interactions accumulate through depth and organize the global propagation dynamics inside the network. Remark 2.1. While self-attention provides the principal mechanism through which information is exchanged across features, the standard Transformer block contains more components: the attention layer is followed by residual connections, normalization layers, and multilayer perceptrons (MLPs), all of which introduce additional nonlinear transformations of the hidden representations. The normalization layers rescale and stabilize the representations, whereas the MLP sublayers enrich the feature representations through nonlinear mixing in the latent space. Consequently, the exact evolution of the hidden states cannot be described solely in terms of attention matrices. Nevertheless, because self-attention is the component that explicitly governs feature-to-feature interactions, it provides a natural and tractable framework for interior interpretability. 2.2. Attention rollout Attention rollout [1] provides a global description of how local attention interactions accumulate across the sequence of attention layers. Given an input datum x∈ℝdx ^d, we define the local rollout operator as the ordered composition R(x)=(L)(x)(L−1)(x)⋯(1)(x)∈ℝd×d, R(x)= A^(L)(x) A^(L-1)(x)\,·s\, A^(1)(x) ^d× d, (2.1) with (ℓ)(x)=12(I+A(ℓ)(x)). A^( )(x)= 12 (I+A^( )(x) ). Since each (ℓ) A^( ) is row-stochastic, the rollout matrix R(x)R(x) is also row-stochastic for all x∈ℝdx ^d. Its coefficient Rij(x)R_ij(x) quantifies the cumulative attention-mediated propagation from source feature j to destination feature i under the rollout approximation. To obtain a cumulative description of these feature interactions at the level of a dataset xnn=1N⊂ℝd\x_n\_n=1^N ^d, we also define the corresponding global rollout operator =1N∑n=1NR(xn)∈ℝd×d, R= 1N _n=1^NR(x_n) ^d× d, (2.2) which provides a coarse-grained description of the dataset-averaged geometry encoded by rollout. Notice that also R is row-stochastic by construction. Throughout this paper, attention rollout is interpreted as an attention-mediated propagation operator. Accordingly, expressions such as propagation, information flow, and influence always refer to the propagation encoded by the rollout operator, rather than to the exact evolution of the Transformer hidden states. As discussed in Remark 2.1, the latter also depends on value projections, output projections, normalization layers, and MLP transformations, which are not represented by the rollout operator. Beyond providing a compact description of attention-mediated propagation, rollout raises a natural mathematical question. For each input x, the local rollout R(x)R(x) is an ordered product of row-stochastic propagation matrices, whereas the global rollout R is the dataset average of these local products. The structural result below applies to either operator because both are row-stochastic. The depth-dependent product estimate developed in Appendix A is first established for the local rollout and then transferred to the global rollout by convexity. To state the structural result, we recall the Dobrushin contraction coefficient (see [11] or [33, Chapter 4, Definition 4.6]), which quantifies the contraction of row-wise oscillations induced by a row-stochastic matrix. Definition 2.2. Given a row-stochastic matrix P∈ℝd×dP ^d× d, the Dobrushin contraction coefficient is defined by κ(P)≔1−mini,k∈⟦d⟧∑j=1dminPij,Pkj. κ(P) 1- _i,k∈ d _j=1^d \P_ij,P_kj\. Equivalently, κ(P)=12maxi,k∈⟦d⟧‖pi−pk‖1, κ(P)= 12 _i,k∈ d \|p_i-p_k\|_1, (2.3) where pip_i denotes the i-th row of P. Theorem 2.3. Let ℛ=(ℛij)i,j=1d∈ℝd×dR=(R_ij)_i,j=1^d ^d× d be a row-stochastic rollout operator, and let κ(ℛ)κ(R) denote its Dobrushin coefficient. Define v≔d−1ℛ⊤v d^-1R 1, where 1 denotes the all-ones vector in ℝdR^d. Then v is a probability vector, the matrix Πℛ≔v⊤ _R 1v is row-stochastic and has rank one and, if we define ‖ℛ−Πℛ‖∞,1≔maxi∈⟦d⟧∑j=1d|(ℛ−Πℛ)ij|, \|R- _R\|_∞,1 _i∈ d _j=1^d|(R- _R)_ij|, we have κ(ℛ)≤‖ℛ−Πℛ‖∞,1≤2(d−1d)κ(ℛ). κ(R)≤\|R- _R\|_∞,1≤ 2 ( d-1d )κ(R). Notice that Theorem 2.3 is stated for a generic rollout operator ℛR, allowing it to be applied both to local rollout matrices R(x)R(x) and to the global rollout R defined in (2.2). Theorem 2.3 shows that, when the Dobrushin coefficient is small, the rollout operator ℛR is well approximated by Πℛ=v⊤ _R= 1v , with v=d−1ℛ⊤v=d^-1R 1. Consequently, every row of ℛR is close to the same probability vector v. More precisely, for every destination index i∈⟦d⟧i∈ d , ∑j=1d|ℛij−vj|≤2(d−1d)κ(ℛ). _j=1^d |R_ij-v_j |≤ 2 ( d-1d )κ(R). Thus, when κ(ℛ)κ(R) is small, the distribution of source-feature contributions received by the different destination features depends only weakly on the destination index. Since ∑i=1dℛij=dvj, _i=1^dR_ij=dv_j, the coefficient vjv_j is exactly the normalized column mass of source feature j. Thus the common row profile v identifies how propagation is distributed among the source features, while the dependence on the destination feature is approximately lost. The contraction estimate used in Theorem 2.3 is a direct finite-dimensional consequence of classical results on stochastic matrices and ergodicity coefficients developed in the theory of nonhomogeneous Markov chains; see [11] and [33, Section 4.3]. Related contraction phenomena have also appeared in mathematical analyses of deep attention mechanisms [3, 4, 12, 18, 19]. Its role here is to provide an analytical lens for attention rollout: it identifies the normalized column sums as the mean row profile and quantifies when the destination-dependent variation around that profile is small. Remark 2.4. No novelty is claimed for the Doeblin–Dobrushin contraction estimate itself. The contribution made here is its application and interpretation in the attention-rollout setting, including the distinction between contraction of the rows and heterogeneity of the resulting propagation profile. This distinction motivates the propagation scores used in the empirical analysis. 2.3. Rollout propagation scores Theorem 2.3 identifies the normalized column-sum vector v as the common row profile appearing in the rank-one approximation of the rollout operator. This gives the column sums a precise structural interpretation: when the Dobrushin coefficient is small, they encode the remaining dependence of the common propagation profile on the features after the row-wise variability has contracted. Motivated by this observation, we use the column sums as propagation scores. Given the global rollout matrix ∈ℝd×d R ^d× d introduced in (2.2), we therefore define the rollout column-score sj=∑i=1dij for all j∈⟦d⟧, s_j= _i=1^d R_ij for all j∈ d , (2.4) which measures the total propagated mass associated with source feature j, aggregated over all destination features. Features with larger sjs_j therefore contribute more strongly to the global propagation profile encoded by the rollout operator. Consequently, sorting the features according to sjs_j yields a propagation-based ranking describing their relative contribution to that profile. Notice that, since sj=dvjs_j=dv_j, where v is the probability profile introduced in Theorem 2.3, the scores sjs_j differ from vjv_j only by the constant normalization factor d and therefore induce exactly the same ranking. Theorem 2.3 provides the structural motivation for these scores. Because rollout omits value and output projections, normalization, and MLP transformations, the scores should not be interpreted as causal effects or as faithful attributions of the complete predictor. 3. Experimental setup and interpretability pipeline In this section, we describe the experimental framework used to examine the proposed propagation analysis, including the datasets, the models, and the analysis pipeline. 3.1. Dataset For our experiments, we considered one biomedical dataset and one language-model-generated tabular dataset for age prediction. Each sample is described by a set of metabolic, biochemical, clinical, and possibly environmental variables. Their main characteristics are summarized in Table 3.1. Dataset 1 consists of anonymized metabolomic measurements obtained from an existing biomedical cohort under a data-sharing agreement with the data owners. Each participant contributed a single observation. The data are not publicly available because of ethical and legal restrictions. The original cohort design, recruitment procedures, laboratory protocols, and ethical approvals are described in [24], which should be consulted for clinical details beyond the scope of the present methodological study. Dataset 2, instead, consists of synthetic tabular data generated using ChatGPT. The generation prompt requested plausible ranges and qualitative relationships among common metabolic, biochemical, and clinical variables used in age-prediction studies. The dataset was not intended to reproduce the statistical distribution of any specific cohort, provide known ground-truth feature effects, or serve as a benchmark dataset. It is used only as a complementary qualitative example. The complete dataset, including the original csv dataframe, is publicly available in the project repository [7] and is sufficient to reproduce the reported experiments. Prior to model training, observations with missing values in any of the input variables were excluded, resulting in a complete-case dataset. The remaining observations were randomly divided into training, validation, and test subsets in proportions 80%80\%, 10%10\%, and 10%10\%, respectively. All subsequent data-dependent preprocessing decisions were made exclusively on the training subset. To reduce redundancy among the input variables, we computed pairwise Pearson correlations on the training data. For every pair of variables whose absolute correlation was at least 0.850.85, one variable was removed. The resulting feature set and column order were then fixed and applied unchanged to the validation and test subsets. Finally, the mean and standard deviation used to normalize each input variable were estimated from the training subset only and applied without refitting to the validation and test subsets. The complete preprocessing pipeline is available on GitHub [7]. Table 3.1. Summary of the datasets used in the experimental study. For each dataset, we report the number of samples in the train/validation/test split and the number of input features before and after the preprocessing pipeline. Dataset 1 Dataset 2 # train data 12744 24000 # validation data 1593 3000 # test data 1594 3000 # features before processing 75 20 # features after processing 72 18 3.2. Model and training setup All the models considered in our experiments share the same general structure: an embedding layer followed by a stack of Transformer encoder blocks and a final linear regression head to generate the predictions. Each scalar input feature xjx_j is mapped to a token through a feature-specific affine embedding xjwj+ejx_jw_j+e_j, where wj,ej∈ℝpw_j,e_j ^p are learned parameters associated with feature j. The embedding dimension is fixed to p=64p=64 and each attention block employs 44 attention heads of dimension ph=16p_h=16. The feedforward sublayers have hidden dimension 128128 and use 0.10.1 dropout. No positional encoding is used because the metabolomic variables have no canonical sequential ordering. Feature identity is nevertheless encoded by the feature-specific parameters wjw_j and eje_j; thus, the architecture does not impose a sequence order but is not permutation-invariant with respect to arbitrary reassignment of feature identities. The main parameter varying across architectures is the number of self-attention layers. More precisely, we train Transformer models with progressively increasing depth, ranging from shallow to deeper configurations. This family of models allows us to study how the attention mechanisms evolve as feature representations propagate through longer sequences of interaction operators. As shown in Section 4, deeper models in the present experiments exhibit increasing row homogenization of the rollout operator. Training is performed using the Adam optimizer with learning rate 10−310^-3, weight decay 10−510^-5, and batch size 128128, using early stopping and ReduceLROnPlateau. The loss function employed is SmoothL1Loss. Predictive performance is evaluated on the test split using standard regression metrics: Root Mean Squared Error (RMSE), Mean Absolute Error (MAE), and R2R^2. Since the target variable is normalized during training, RMSE and MAE are computed after denormalizing the predictions and ground-truth labels and are therefore reported in units of years. Simulations were performed on an Acer TravelMate P414-53 laptop with OS Ubuntu 24.04.3 LTS, 13th Gen Intel® Core™ i5-1335U × 12 processor and 32 GiB of RAM memory. Moreover, all experiments were repeated using three independent random seeds. Unless otherwise stated, quantitative performance metrics are reported as mean and standard deviation across the three independent training runs. 3.3. Interpretability methods The objective of our empirical analysis is to compare the three viewpoints introduced above. PCA supplies a variance-based score, GradientExplainer supplies an approximate prediction-attribution score, and rollout supplies a propagation score. Variance viewpoint: PCA-based feature importance. As an unsupervised reference, we perform PCA on the normalized data. This provides a description of the intrinsic variance structure of the metabolomic feature space independently of the predictor. Feature importance is quantified through the singular-value-weighted loading score rj=∑k=1dωk|ζkj|, for all j∈⟦d⟧, r_j= _k=1^d _k| _kj|, for all j∈ d , (3.1) which measures the overall contribution of feature j to the principal-component representation of the data, weighting the contribution of each principal component by its singular value. Here, ωk∈ℝζk=(ζkj)j=1d∈ℝd, for all k∈⟦d⟧ array[]l _k \\ _k=( _kj)_j=1^d ^d array, for all k∈ d are the singular values and right singular vectors in the SVD decomposition of the data matrix. Propagation viewpoint: rollout propagation scores. We analyze feature interactions through the global rollout matrix R constructed as in (2.2). Besides computing the rollout propagation scores (2.4), we also quantify the structural organization of the propagation operator by evaluating its Dobrushin coefficient κ()κ( R) and the approximation error ‖−Π‖∞,1\| R- _ R\|_∞,1 appearing in Theorem 2.3. These quantities measure the extent to which the rollout operator is close to a rank-one propagation profile at each Transformer depth. Attribution viewpoint: GradientExplainer scores. We estimate feature-attribution scores with the GradientExplainer implementation of the Python SHAP package [27]. This implementation computes expected-gradients attributions relative to an empirical background distribution and is used as an approximation to SHAP values for differentiable models. These numerical scores need not coincide exactly with either interventional or conditional Shapley values. The exact interventional Shapley operator used in Proposition A.6 is defined separately in Appendix A; that idealized proposition is not a validation theorem for the empirical GradientExplainer scores. The background supplied to GradientExplainer consists of a random subset of 256256 samples from the test set. Methodologically, using the test set as the GradientExplainer background does not introduce information leakage, because the background distribution is used exclusively after the models have been fully trained to compute post hoc SHAP values. It is not involved in training, hyperparameter tuning, model selection, or performance evaluation, and therefore cannot influence the learned models or their predictive performance. Moreover, for our dataset, the choice of training or test background has negligible impact. We verified that the training and test sets are statistically indistinguishable by training a Random Forest classifier to discriminate between them. The classifier achieved an AUC of 0.495±0.0200.495± 0.020, i.e., no better than random guessing, indicating that the two sets are effectively drawn from the same distribution. 4. Experimental results We first analyze the real metabolomic dataset (Dataset 1), which serves as the main case study throughout the paper. We then repeat selected analyses on the synthetic Dataset 2 as a complementary qualitative example. All experiments presented in this section, together with the corresponding code, are available on GitHub [7]. 4.1. Experimental results on Dataset 1 We present here the results of our experiments performed on Dataset 1. 4.1.1. Training and predictive performance To investigate the effect of depth, we trained Transformer architectures with varying numbers of self-attention layers ranging from 33 to 2020. The predictive performance obtained on the test set is reported in Table 4.1. Depth R2R^2 RMSE MAE Best validation loss 33 0.6942±0.00480.6942± 0.0048 7.5414±0.05977.5414± 0.0597 5.8941±0.05295.8941± 0.0529 0.1574±0.00170.1574± 0.0017 55 0.6980±0.00830.6980± 0.0083 7.4942±0.10357.4942± 0.1035 5.8787±0.06235.8787± 0.0623 0.1584±0.00090.1584± 0.0009 77 0.6867±0.01210.6867± 0.0121 7.6322±0.14747.6322± 0.1474 5.9504±0.05395.9504± 0.0539 0.1568±0.00160.1568± 0.0016 99 0.6866±0.00190.6866± 0.0019 7.6354±0.02307.6354± 0.0230 6.0094±0.02276.0094± 0.0227 0.1577±0.00120.1577± 0.0012 1111 0.6919±0.00560.6919± 0.0056 7.5705±0.06867.5705± 0.0686 5.9114±0.04965.9114± 0.0496 0.1582±0.00070.1582± 0.0007 1313 0.6872±0.00310.6872± 0.0031 7.6277±0.03807.6277± 0.0380 5.9542±0.02385.9542± 0.0238 0.1567±0.00210.1567± 0.0021 1515 0.6861±0.00590.6861± 0.0059 7.6424±0.07147.6424± 0.0714 5.9406±0.01505.9406± 0.0150 0.1569±0.00320.1569± 0.0032 1717 0.6910±0.00640.6910± 0.0064 7.5816±0.07877.5816± 0.0787 5.9594±0.05115.9594± 0.0511 0.1596±0.00250.1596± 0.0025 2020 0.6891±0.00310.6891± 0.0031 7.6043±0.03777.6043± 0.0377 5.9613±0.03815.9613± 0.0381 0.1618±0.00150.1618± 0.0015 Table 4.1. Training and predictive performance over Dataset 1 of Transformer architectures with increasing depth. The table reports predictive accuracy on the test split for Transformer models with different numbers of self-attention layers. Results are presented as the mean ± standard deviation over independent runs using different random seeds, summarizing the variability observed across initializations. The reported metrics include the coefficient of determination R2R^2, Root Mean Squared Error (RMSE), Mean Absolute Error (MAE), and the best validation loss achieved during training. Overall, performance remains stable across the considered depths. The coefficient of determination stays in the approximate range R2∈[0.68,0.70]R^2∈[0.68,0.70], while RMSE and MAE vary only moderately. Figure 4.1 complements this analysis by reporting the training and validation losses for representative shallow, intermediate, and deep Transformer architectures. All displayed curves converge rapidly during the first epochs and do not show a large sustained separation between training and validation loss. Figure 4.1. Evolution of training and validation losses during optimization over Dataset 1 for representative shallow (3-layer), intermediate (11-layer), and deep (20-layer) Transformer architectures. Shaded regions indicate ±1± 1 standard deviation over three runs. Across the depths considered, the observed differences in predictive performance are modest relative to the overall performance level, and no monotone improvement with depth is apparent. These results are descriptive and are not intended to establish statistical equivalence among architectures. This observation motivates the analysis developed in Sections 4.1.2–4.1.4: while the reported predictive performance remains broadly stable, the measured attention-propagation patterns continue to change with depth. 4.1.2. Attention-based analysis Motivated by the structural characterization provided by Theorem 2.3 and by the depth-dependent contraction estimate (A.8) proved in Appendix A, we now examine how the corresponding column-wise organization is manifested in the trained Transformers. Figure 4.2 displays the global rollout operator together with its rank-one approximation and the associated error. Experimentally, as the depth increases, the rollout operator becomes increasingly well described by a common propagation profile. For each depth, the corresponding approximation is quantitatively characterized by Theorem 2.3 through the measured Dobrushin coefficient. The learned propagation profile remains nonuniform, showing empirically that some source features contribute more strongly to the propagation dynamics than others. Notice that Theorem 2.3 does not predict this nonuniformity; it is an empirical property of the learned rollout operator. Figure 4.2. Global rollout matrix (2.2) computed over the test split of Dataset 1 for one representative trained 11-layer Transformer realization, together with its rank-one approximation and the corresponding approximation error. The small residual error illustrates that the rollout operator is well approximated by a rank-one stochastic matrix, in accordance with Theorem 2.3. Figure 4.3 reports the Dobrushin contraction coefficient κ()κ( R) together with the approximation error ‖−Π‖∞,1\| R- _ R\|_∞,1. By Theorem 2.3, these two quantities are analytically equivalent up to fixed multiplicative constants for every rollout operator. Moreover, estimate (A.8) controls the Dobrushin coefficient of the global rollout, making explicit the factors associated with residual averaging and the layerwise attention matrices. In the architectures considered here, both the Dobrushin coefficient and the approximation error are observed experimentally to decrease rapidly with depth. Depth Dobrushin coefficient κ(R)κ( R) Error ‖−Π‖∞,1\| R- _ R\|_∞,1 33 1.7×10−1±5.5×10−51.7× 10^-1± 5.5× 10^-5 2.8×10−1±5.9×10−32.8× 10^-1± 5.9× 10^-3 55 4.7×10−2±3.5×10−34.7× 10^-2± 3.5× 10^-3 7.3×10−2±1.4×10−37.3× 10^-2± 1.4× 10^-3 77 1.2×10−2±1.2×10−31.2× 10^-2± 1.2× 10^-3 1.8×10−2±1.6×10−31.8× 10^-2± 1.6× 10^-3 99 2.9×10−3±2.5×10−42.9× 10^-3± 2.5× 10^-4 4.3×10−3±1.9×10−44.3× 10^-3± 1.9× 10^-4 1111 7.6×10−4±1.2×10−57.6× 10^-4± 1.2× 10^-5 1.1×10−3±1.1×10−51.1× 10^-3± 1.1× 10^-5 1313 1.9×10−4±1.1×10−51.9× 10^-4± 1.1× 10^-5 2.7×10−4±3.8×10−62.7× 10^-4± 3.8× 10^-6 1515 5.4×10−5±3.1×10−65.4× 10^-5± 3.1× 10^-6 7.8×10−5±4.6×10−77.8× 10^-5± 4.6× 10^-7 1717 1.4×10−5±7.1×10−71.4× 10^-5± 7.1× 10^-7 1.8×10−5±2.3×10−71.8× 10^-5± 2.3× 10^-7 2020 1.5×10−6±2.2×10−71.5× 10^-6± 2.2× 10^-7 2.3×10−6±1.2×10−72.3× 10^-6± 1.2× 10^-7 Figure 4.3. Decay of the Dobrushin contraction coefficient κ()κ( R) and the approximation error ‖−Π‖∞,1\| R- _ R\|_∞,1 as functions of the Transformer depth on Dataset 1. The table reports the corresponding mean ± standard deviation over independent training runs. Theorem 2.3 guarantees that these two quantities remain quantitatively equivalent up to fixed multiplicative constants for every rollout operator. The observed decrease in the trained models is compatible with the product mechanism identified in Appendix A, but the present measurements do not quantify the separate contribution of training to contraction. Section 4.1.3 therefore addresses a narrower empirical question by comparing the propagation profiles of trained and randomly initialized Transformers. 4.1.3. Effect of learning on the propagation profile The previous section documented depth-dependent contraction for the trained Transformers. We now address a separate question: whether trained and randomly initialized models exhibit different propagation profiles. The comparison below concerns the profiles and their rankings; it does not compare their Dobrushin coefficients. Our analysis is based on two complementary experiments. First, we compare the propagation profiles v obtained from trained and randomly initialized (untrained) models using their ℓ1 ^1 distance. Second, we measure the ℓ1 ^1 distance of each profile from the uniform distribution u=d−1u=d^-1 1. Together, these descriptive quantities assess how the distribution of rollout mass differs between the two groups. The results are reported in Tables 4.2 and 4.3. Table 4.2 shows that the trained and randomly initialized profiles differ across all considered depths, with an ℓ1 ^1 distance of approximately 0.150.15 in these experiments. This comparison does not by itself establish the cause of the difference or whether the resulting ranking is predictively relevant. Table 4.3 compares each propagation profile with the uniform distribution. The reported randomly initialized profiles remain close to uniform, whereas the trained profiles have larger distances from uniformity. This descriptive comparison indicates more heterogeneous rollout profiles in the trained models considered here. Depth ℓ1 ^1 distance 3 0.1438±0.00060.1438± 0.0006 5 0.1644±0.02800.1644± 0.0280 7 0.1529±0.02120.1529± 0.0212 9 0.1452±0.01810.1452± 0.0181 11 0.1511±0.00380.1511± 0.0038 13 0.1489±0.00620.1489± 0.0062 15 0.1618±0.01640.1618± 0.0164 17 0.1573±0.01140.1573± 0.0114 20 0.1318±0.00870.1318± 0.0087 Table 4.2. ℓ1 ^1 distance between the propagation profiles v of trained and randomly initialized Transformer models as a function of the depth. Results are reported as mean ± standard deviation over independent random seeds. Depth Trained Random init. 3 0.1434±0.00150.1434± 0.0015 0.018273±0.0015440.018273± 0.001544 5 0.1616±0.03080.1616± 0.0308 0.018319±0.0014590.018319± 0.001459 7 0.1506±0.02010.1506± 0.0201 0.018258±0.0014790.018258± 0.001479 9 0.1428±0.01880.1428± 0.0188 0.018272±0.0014840.018272± 0.001484 11 0.1511±0.00220.1511± 0.0022 0.018272±0.0014870.018272± 0.001487 13 0.1469±0.00630.1469± 0.0063 0.018272±0.0014870.018272± 0.001487 15 0.1606±0.01680.1606± 0.0168 0.018272±0.0014870.018272± 0.001487 17 0.1553±0.00970.1553± 0.0097 0.018272±0.0014870.018272± 0.001487 20 0.1314±0.00730.1314± 0.0073 0.018272±0.0014870.018272± 0.001487 Table 4.3. ℓ1 ^1 distance between the propagation profile and the uniform distribution for trained and randomly initialized Transformer models as a function of the depth. Results are reported as mean ± standard deviation over independent random seeds. To illustrate how the profile differences appear in the rollout ranking, Figure 4.4 compares the six highest-ranked variables obtained from trained and randomly initialized Transformers for the representative 11-layer architecture. We stress that this comparison is intended only as a visual illustration of the quantitative differences reported in Tables 4.2 and 4.3. Figure 4.4. Comparison of rollout propagation scores over Dataset 1 between trained and randomly initialized Transformers with 11 attention layers. Consistent with the profile comparison, the trained and randomly initialized models produce different rollout rankings. Some variables highly ranked after training, including albumin and LDL cholesterol, are established clinical variables discussed in Appendix B. This provides a qualitative plausibility check only: it does not demonstrate that rollout scores are faithful, causal, or predictively relevant. These results show that the propagation profile and its ranking differ between trained and randomly initialized models. Together with the product analysis of Section A.1.1, they motivate a distinction between the structural contraction mechanism of rollout composition and the profile realized after optimization. Establishing a causal relation between training, the profile, and predictive behavior requires additional intervention-based experiments. 4.1.4. Comparison of feature-ranking scores Section 4.1.3 documented differences between the propagation profiles and rankings of the trained and randomly initialized models. We now compare three complementary scoring systems for the trained Transformer: • the PCA scores defined in (3.1), obtained from the singular-value-weighted principal-component representation of the data; • the rollout column-sum scores defined in (2.4), measuring the rollout mass associated with each source feature; • GradientExplainer scores, providing approximate prediction attributions relative to the selected background distribution. Before comparing the different interpretability methods, we describe the cross-seed stability of the corresponding feature rankings. Figure 4.5 reports the mean pairwise Spearman correlation of the rollout and GradientExplainer rankings across independent training runs. Figure 4.5. Cross-seed stability of rollout- and GradientExplainer-based feature rankings over Dataset 1. For each Transformer depth, the figure reports the mean pairwise Spearman correlation between rankings obtained from independent training runs with different random seeds. The solid curves represent the mean across runs, and the shaded bands indicate one standard deviation. Within the three reported initializations, the GradientExplainer rankings have higher mean pairwise Spearman correlations than the rollout rankings. The latter range from approximately 0.650.65 to 0.800.80. These values describe cross-seed stability in the present experiment; they do not establish population-level reproducibility. Figure 4.6 reports the six highest-ranked variables according to these scores for one representative trained Transformer of depth 1111. This intermediate depth is used throughout the qualitative illustrations because its measured Dobrushin coefficient is already small and its rollout operator already exhibits the rank-one propagation structure characterized by Theorem 2.3. Table 4.4 reports the meaning of the feature abbreviations. Figure 4.6. Comparison of top-66 rankings for Dataset 1 obtained with PCA scores, rollout propagation scores, and GradientExplainer feature-attribution scores. Feature Abbreviation Feature type 1,5-Anhydrosorbitol 1,5-AG Metabolite 3-Hydroxyisobutyric acid 3-HIB Metabolite Acetone Acetone Metabolite Albumin Alb Protein Dimethylsulfone DMSO2 Metabolite Erythrocyte Sedimentation Rate ESR Inflammatory marker Glucose Gluc Metabolite Glutamine to Glutamate ratio Gln/Glu Derived ratio LDL Cholesterol LDL-C Clinical biomarker Triglycerides TG Clinical biomarker Urea Urea Clinical biomarker Table 4.4. Correspondence between feature abbreviations and full variable names in Figure 4.6. Variables are grouped by category (metabolites, clinical biomarkers, proteins, inflammatory markers and derived ratios). The PCA ranking differs substantially from both model-dependent rankings. The rollout and attribution rankings, by contrast, share several features. This difference is unsurprising because PCA ranks variables according to the variance structure of the input data alone and therefore has no access to the prediction task. This is also reflected in the nearly uniform PCA scores shown in Figure 4.6, where the highest-ranked variables receive very similar scores. To compare the three interpretability methods quantitatively, Figure 4.7 summarizes their agreement across Transformer depths. The upper panel reports the normalized overlap between the six highest-ranked features identified by each pair of methods, averaged over three random initializations. The lower panel reports the Spearman correlation between the complete rankings, likewise averaged over the three initializations. Both panels are descriptive summaries of one cohort and a small number of training runs. Figure 4.7. Descriptive comparison of rankings obtained with PCA scores, rollout propagation scores, and GradientExplainer feature-attribution scores over Dataset 1. The left panel reports normalized top-66 overlap, and the right panel reports Spearman correlation across all 7272 variables. Values are averaged over three random initializations. The comparisons involving PCA exhibit weak descriptive agreement. Rollout and PCA share few top-ranked variables, while GradientExplainer and PCA share, on average, approximately two of their six highest-ranked variables. These observations are summaries of the present cohort and training runs, not evidence of generalizable agreement. Rollout and GradientExplainer identify between two and five common variables among their six highest-ranked features across the individual runs, with average normalized overlaps ranging from 0.500.50 to 0.780.78 across depths. Thus, in this dataset, their agreement is most visible at the top of the rankings. This descriptive comparison is sensitive to the selected top-66 threshold. When the comparison is extended to the complete rankings, the Spearman correlations remain close to zero for all method pairs and depths. The top-66 overlap and the full-ranking comparison therefore describe different aspects of the same results: rollout and GradientExplainer share some highly ranked variables in this cohort while assigning substantially different orders to many variables of intermediate or low rank. 4.2. Experimental results on Dataset 2 As an auxiliary illustration, we repeat the main propagation analyses on the synthetic Dataset 2. This analysis is not intended to test generalization beyond Dataset 1; it only shows how the same diagnostics behave in a second tabular example. The predictive performance obtained on Dataset 2 is reported in Table 4.5. As in Dataset 1, the point estimates vary little across the considered depths; no monotone association between depth and predictive performance is apparent. Depth R2R^2 RMSE MAE Best validation loss 33 0.9035±0.00050.9035± 0.0005 6.4480±0.01526.4480± 0.0152 5.1568±0.02645.1568± 0.0264 0.0504±0.00010.0504± 0.0001 55 0.9039±0.00010.9039± 0.0001 6.4334±0.00456.4334± 0.0045 5.1136±0.00815.1136± 0.0081 0.0501±0.00010.0501± 0.0001 77 0.9039±0.00010.9039± 0.0001 6.4352±0.00136.4352± 0.0013 5.1140±0.00435.1140± 0.0043 0.0501±0.00010.0501± 0.0001 99 0.9041±0.00050.9041± 0.0005 6.4295±0.01596.4295± 0.0159 5.1076±0.00615.1076± 0.0061 0.0511±0.00040.0511± 0.0004 1111 0.9042±0.00040.9042± 0.0004 6.4247±0.01226.4247± 0.0122 5.1229±0.03395.1229± 0.0339 0.0501±0.00020.0501± 0.0002 1313 0.9037±0.00080.9037± 0.0008 6.4408±0.02716.4408± 0.0271 5.1403±0.04235.1403± 0.0423 0.0500±0.00020.0500± 0.0002 1515 0.9032±0.00070.9032± 0.0007 6.4572±0.02286.4572± 0.0228 5.1402±0.02975.1402± 0.0297 0.0503±0.00050.0503± 0.0005 1717 0.9034±0.00070.9034± 0.0007 6.4526±0.02516.4526± 0.0251 5.1540±0.02465.1540± 0.0246 0.0505±0.00020.0505± 0.0002 2020 0.9034±0.00010.9034± 0.0001 6.4502±0.00426.4502± 0.0042 5.1386±0.01175.1386± 0.0117 0.0502±0.00010.0502± 0.0001 Table 4.5. Training and predictive performance over Dataset 2 of Transformer architectures with increasing depth. The table reports predictive accuracy on the test split for Transformer models with different numbers of self-attention layers. Results are presented as the mean ± standard deviation over independent runs using different random seeds, summarizing the variability observed across initializations. The reported metrics include the coefficient of determination R2R^2, Root Mean Squared Error (RMSE), Mean Absolute Error (MAE), and the best validation loss achieved during training. While the predictive point estimates vary little, Figures 4.8 and 4.9 show that the measured rollout matrices become more contractive with depth. This is a similar qualitative depth-dependent pattern to that observed for Dataset 1. Figure 4.8. Global rollout matrices (2.2) computed over the test split of Dataset 2 for shallow (3-layer), intermediate (11-layer), and deep (20-layer) Transformers. Figure 4.9. Evolution of the Dobrushin contraction coefficient κ()κ( R) and the approximation error ‖−Π‖∞,1\| R- _ R\|_∞,1 of the global rollout operator as functions of the Transformer depth on Dataset 2. Both quantities decrease rapidly with depth, showing that the rollout operator becomes progressively closer to its rank-one approximation, consistently with Theorem 2.3 and the depth-dependent estimate in Appendix A. After documenting the contraction and rank-one approximation pattern in this auxiliary example, we compare the three analytical viewpoints descriptively. Figure 4.10 reports the Spearman correlations between the rollout, GradientExplainer, and PCA rankings for Transformer architectures of increasing depth. Figure 4.10. Comparison of rankings obtained with PCA scores, rollout propagation scores, and GradientExplainer feature-attribution scores over Dataset 2 for Transformer architectures of varying depth. The heatmap reports the Spearman correlation computed over the complete rankings of all 1818 features. In this synthetic example, correlations involving PCA are weak, while rollout and GradientExplainer rankings have moderate positive agreement at several Transformer depths and remain far from identical. These observations describe this generated dataset only. Overall, Dataset 2 shows two patterns similar to the point estimates in Dataset 1: progressive contraction of the trained rollout operator and a non-uniform propagation profile. Because the dataset was generated by a language model rather than a fully specified stochastic generator with known ground truth, this experiment is an illustration rather than an independent validation or evidence of generality. 4.3. Limitations of the empirical analysis The empirical results should be interpreted within several limitations. First, attention rollout represents only propagation mediated by attention. It omits value and output projections, normalization layers, MLP sublayers, and other nonlinear transformations in the complete Transformer. Consequently, rollout propagation scores are neither causal effects nor guaranteed faithful attributions of the prediction. Second, the product analysis provides a layerwise upper bound on rollout contraction but does not by itself prove decay with depth; moreover, the experiments report trained contraction only, while the trained–random-initialization comparison concerns propagation profiles rather than contraction coefficients. Intervention-based controls would be required to establish the predictive relevance of individual rollout-ranked variables. Third, the real-data analysis uses a single biomedical cohort and three random initializations, so the ranking comparisons are descriptive and exploratory. Fourth, Dataset 2 was generated with a language model rather than a fully specified stochastic data-generating process and does not provide ground-truth feature effects. Finally, the GradientExplainer background used in the reported experiments was sampled from the test set. It was employed exclusively for post-hoc interpretability after model training and was not involved in model fitting, hyperparameter selection, or checkpoint selection. Accordingly, the conclusions of this work concern the mathematical characterization and empirical illustration of the propagation encoded by attention rollout. 5. Conclusions and perspectives We introduced interior interpretability as the analysis of explicitly defined operators that summarize selected interactions inside a learned architecture, and instantiated this perspective with attention rollout in tabular Transformers. The scope of the resulting explanation is operator-specific: rollout describes attention-mediated propagation, not the complete hidden-state dynamics and not a causal decomposition of the prediction. Using a direct finite-dimensional consequence of classical Doeblin–Dobrushin theory, we characterized the rank-one regime of any row-stochastic rollout operator. Its Dobrushin coefficient controls the distance to a canonical rank-one matrix, and the normalized column masses are exactly the mean row profile. This identifies two mathematically distinct aspects of rollout: contraction determines how much destination-dependent variation is lost, whereas the profile records how the remaining propagated mass is distributed among source features. Neither rank-one structure nor a nonuniform profile, by itself, establishes learning, predictive relevance, or attribution faithfulness. In the trained tabular Transformers studied here, the measured Dobrushin coefficient and rank-one approximation error decreased with depth. The trained and randomly initialized models also displayed different propagation profiles in the reported runs, but this experiment does not isolate the cause or predictive consequence of those differences. On Dataset 1, exploratory comparisons with PCA and GradientExplainer scores revealed localized overlap among some highly ranked variables and weak agreement across complete rankings. These observations motivate using rollout as a complementary diagnostic of propagation geometry while underscoring that variance, propagation, and prediction attribution answer different questions. The broader methodological lesson is that a defensible interpretation of an internal quantity requires its underlying operator, assumptions, invariances, scope, and failure modes to be made explicit. This principle suggests a research program with four priorities: (1) construct computation-aware propagation operators that incorporate value and output projections, residual paths, normalization, and MLP transformations, for example through local Jacobian or transport formulations; (2) formulate propagation–attribution relationships as explicit, falsifiable hypotheses, supported by interventions and counterexamples rather than ranking agreement alone; (3) quantify the formation, stability, and uncertainty of propagation profiles across initializations, subjects, cohorts, architectures, and appropriate random-model baselines; (4) extend the contraction–profile decomposition to local, subgroup-specific, continuous-depth, and mean-field operators, where population heterogeneity and depth scaling can be studied directly. Interior interpretability is therefore proposed not as a universal feature-importance score, but as an operator-first framework for deriving mathematically scoped claims and empirical tests that can be falsified, refined, and compared across models. Acknowledgments The authors would like to express their sincere gratitude to Alain Ibáñez de Opakua, José María Mato and Óscar Millet (ATLAS Molecular Pharma, CIC bioGUNE, and CIBERehd) for providing the biomedical Dataset 1 used in the experimental part of this work. The experimental investigations presented in this paper would not have been possible without their contribution in collecting, curating, and organizing the dataset, as well as their guidance concerning the biomedical context of the selected variables. The authors also wish to thank Ziqian Li (Chair for Dynamics, Control, Machine Learning, and Numerics at the Friedrich-Alexander-Universität Erlangen-Nürnberg) for his valuable assistance during the final revisions of this work. Data availability The data contained in Dataset 1 were collected through individual medical examinations and are not publicly available due to privacy and ethical restrictions. Access to anonymized data may be granted upon reasonable request and subject to the approval of the relevant data-governance and ethical procedures. Dataset 2 instead is available on the GitHub repository associated with this paper [7]. Appendix A Mathematical results This appendix collects the mathematical results supporting our discussion. A.1. Structural properties of attention rollout We first establish structural statements that apply to any row-stochastic propagation operator. They therefore apply both to each local rollout matrix R(x)R(x) defined in (2.1) and to the dataset-average global rollout R defined in (2.2). Only the later depth-dependent estimate uses the product representation of the local rollout. Throughout the first part of the appendix, we consider a generic matrix ℛ∈ℝd×d, ^d× d, ℛij≥0for every i,j∈⟦d⟧, _ij≥ 0 every i,j∈ d , ∑j=1dℛij=1for every i∈⟦d⟧. _j=1^dR_ij=1 every i∈ d . (A.1) We start by giving the following definition of oscillation in a vector. Definition A.1. For every x=(x1,…,xd)∈ℝdx=(x_1,…,x_d) ^d, we define its oscillation by osc(x)≔maxj∈⟦d⟧xj−minj∈⟦d⟧xj. (x) _j∈ d x_j- _j∈ d x_j. Notice that osc(x)osc(x) defines a seminorm on ℝdR^d. The following result is the finite-dimensional specialization applied to the rollout matrix ℛR of the Doeblin–Dobrushin contraction theorem (see [11] or [33, Theorem 4.14]). We include it for completeness, since the proof is short and because the notation adopted here is tailored to attention rollout operators. Proposition A.2. Let ℛ=(ℛij)i,j=1d∈ℝd×dR=(R_ij)_i,j=1^d ^d× d be a rollout matrix as in (A.1), and let κ(ℛ)κ(R) denote its Dobrushin coefficient given by Definition 2.2. Then 0≤κ(ℛ)≤10≤κ(R)≤ 1 and, for every x∈ℝdx ^d, osc(ℛx)≤κ(ℛ)osc(x). (Rx)≤κ(R)osc(x). In particular, if there exists δ>0δ>0 such that ℛij≥δR_ij≥δ for all i,j∈⟦d⟧i,j∈ d , then δ≤1/dδ≤ 1/d, κ(ℛ)≤1−dδ<1, κ(R)≤ 1-dδ<1, and osc(ℛx)≤κ(ℛ)osc(x)≤(1−dδ)osc(x). (Rx)≤κ(R)osc(x)≤ (1-dδ )osc(x). Proof. Fix x∈ℝdx ^d, and set m≔minj∈⟦d⟧xj,M≔maxj∈⟦d⟧xj. m _j∈ d x_j, M _j∈ d x_j. Then osc(x)=M−mosc(x)=M-m. If M=mM=m, then osc(x)=0osc(x)=0 meaning that x is constant. In this case, since ℛR is row-stochastic, we necessarily have ℛx=xRx=x, and therefore osc(ℛx)=0osc(Rx)=0. Due to the above argument, in what follows we may assume that M>mM>m. Define the vector y∈ℝdy ^d by yj≔xj−mM−m, for all j∈⟦d⟧. y_j x_j-mM-m, for all j∈ d . Since m≤xj≤Mm≤ x_j≤ M, we have 0≤yj≤1, for all j∈⟦d⟧. 0≤ y_j≤ 1, for all j∈ d . Moreover, x=m+(M−m)yx=m1+(M-m)y. Because ℛR is row-stochastic, ℛ=R1=1. Consequently, ℛx=mℛ+(M−m)ℛy=m+(M−m)ℛy. =mR1+(M-m)Ry=m1+(M-m)Ry. Adding a constant vector does not change oscillation, while multiplication by a nonnegative scalar scales the oscillation. Therefore, osc(ℛx)=(M−m)osc(ℛy). (Rx)=(M-m)osc(Ry). (A.2) It remains to estimate osc(ℛy)osc(Ry). For arbitrary i,k∈⟦d⟧i,k∈ d , we have (ℛy)i−(ℛy)k=∑j=1d(ℛij−ℛkj)yj. (Ry)_i-(Ry)_k= _j=1^d (R_ij-R_kj )y_j. (A.3) For each j∈⟦d⟧j∈ d , define cj≔minℛij,ℛkjc_j \R_ij,R_kj\. Then we may write ℛij=cj+pj and ℛkj=cj+qj, _ij=c_j+p_j and _kj=c_j+q_j, where pj≔ℛij−cj≥0 and qj≔ℛkj−cj≥0. p_j _ij-c_j≥ 0 and q_j _kj-c_j≥ 0. Notice that, for every j, at least one of pjp_j and qjq_j is zero. Using this decomposition in (A.3), we obtain (ℛy)i−(ℛy)k=∑j=1dpjyj−∑j=1dqjyj. (Ry)_i-(Ry)_k= _j=1^dp_jy_j- _j=1^dq_jy_j. Moreover, since 0≤yj≤10≤ y_j≤ 1, it follows that 0≤∑j=1dpjyj≤∑j=1dpj and 0≤∑j=1dqjyj≤∑j=1dqj. 0≤ _j=1^dp_jy_j≤ _j=1^dp_j and 0≤ _j=1^dq_jy_j≤ _j=1^dq_j. Therefore, using the row-stochasticity of ℛR, (ℛy)i−(ℛy)k (Ry)_i-(Ry)_k ≤∑j=1dpj=∑j=1d(ℛij−cj)=1−∑j=1dcj=1−∑j=1dminℛij,ℛkj ≤ _j=1^dp_j= _j=1^d (R_ij-c_j )=1- _j=1^dc_j=1- _j=1^d \R_ij,R_kj\ ≤1−minr,q∈⟦d⟧∑j=1dminℛrj,ℛqj=κ(ℛ). ≤ 1- _r,q∈ d _j=1^d \R_rj,R_qj\=κ(R). (A.4) Since i and k are arbitrary, we can choose i∗i_ and k∗k_ such that (ℛy)i∗=maxi∈⟦d⟧(ℛy)i and (ℛy)k∗=mink∈⟦d⟧(ℛy)k. (Ry)_i_ = _i∈ d (Ry)_i and (Ry)_k_ = _k∈ d (Ry)_k. Applying (A.1) to this pair yields osc(ℛy)=(ℛy)i∗−(ℛy)k∗≤κ(ℛ). (Ry)=(Ry)_i_ -(Ry)_k_ ≤κ(R). (A.5) Combining (A.2) and (A.5), we conclude that osc(ℛx)≤(M−m)κ(ℛ)=κ(ℛ)osc(x). (Rx)≤(M-m)κ(R)=κ(R)osc(x). It remains to verify the properties of κ(ℛ)κ(R). For any pair (i,k)(i,k), we have 0≤∑j=1dminℛij,ℛkj≤∑j=1dℛij=1. 0≤ _j=1^d \R_ij,R_kj\≤ _j=1^dR_ij=1. This immediately gives 0≤κ(ℛ)≤10≤κ(R)≤ 1. If, moreover, ℛij>0R_ij>0 for all (i,j)(i,j), then, for every (i,k,j)(i,k,j), minℛij,ℛkj>0. \R_ij,R_kj\>0. Since there are finitely many rows and columns, mini,k∈⟦d⟧∑j=1dminℛij,ℛkj>0. _i,k∈ d _j=1^d \R_ij,R_kj\>0. It follows that κ(ℛ)<1κ(R)<1. Finally, suppose that ℛij≥δ>0,for all i,j∈⟦d⟧R_ij≥δ>0, all i,j∈ d . Since ℛR is row stochastic, we have δ≤1/dδ≤ 1/d. Moreover, for every pair (i,k)(i,k), minℛij,ℛkj≥δ, \R_ij,R_kj\≥δ, and therefore ∑j=1dminℛij,ℛkj≥dδ. _j=1^d \R_ij,R_kj\≥ dδ. Thus, κ(ℛ)≤1−dδκ(R)≤ 1-dδ and we can conclude that osc(ℛx)≤(1−dδ)osc(x)osc(Rx)≤ (1-dδ )osc(x). ∎ We now present the proof of Theorem 2.3. Proof of Theorem 2.3. For all i∈⟦d⟧i∈ d , let ri≔(ℛi1,…,ℛid)∈ℝd r_i (R_i1,…,R_id ) ^d denote the i-th row of ℛR. Since ℛR is row-stochastic, every rir_i belongs to the probability simplex Δd−1≔z∈ℝd:zj≥0for all j∈⟦d⟧,∑j=1dzj=1. ^d-1 \z ^d:\;z_j≥ 0\ for all j∈ d , _j=1^dz_j=1 \. Define v≔1dℛ⊤=1d∑i=1dri∈ℝd. v 1dR 1= 1d _i=1^dr_i ^d. Since v is a convex combination of the rows ri∈Δd−1r_i∈ ^d-1, and since the probability simplex is convex, we have that v∈Δd−1v∈ ^d-1. Now define Πℛ≔v⊤ _R 1v . Every row of Πℛ _R is equal to v⊤v . Therefore, (Πℛ)ij=vj≥0 and ∑j=1d(Πℛ)ij=∑j=1dvj=1 ( _R)_ij=v_j≥ 0 and _j=1^d( _R)_ij= _j=1^dv_j=1 Thus Πℛ _R is row-stochastic. Moreover, since it is the outer product of the nonzero vectors 1 and v, we also have that rank(Πℛ)=1rank( _R)=1. It remains to estimate the distance between ℛR and Πℛ _R. By (2.3), the Dobrushin coefficient admits the representation κ(ℛ)=12maxi,k∈⟦d⟧‖ri−rk‖1. κ(R)= 12 _i,k∈ d \|r_i-r_k\|_1. Equivalently, for every i,k∈⟦d⟧i,k∈ d , we have ‖ri−rk‖1≤2κ(ℛ)\|r_i-r_k\|_1≤ 2κ(R). Fix i∈⟦d⟧i∈ d . Since v is the arithmetic mean of the rows, ri−v r_i-v =ri−1d∑k=1drk=1d∑k=1d(ri−rk). =r_i- 1d _k=1^dr_k= 1d _k=1^d(r_i-r_k). Using the triangle inequality, then we obtain ‖ri−v‖1≤1d∑k=1d‖ri−rk‖1=1d∑k=1k≠id‖ri−rk‖1≤1d∑k=1k≠id2κ(ℛ)=2(d−1d)κ(ℛ). \|r_i-v\|_1≤ 1d _k=1^d\|r_i-r_k\|_1= 1d _ subarrayck=1\\ k≠ i subarray^d\|r_i-r_k\|_1≤ 1d _ subarrayck=1\\ k≠ i subarray^d2κ(R)=2 ( d-1d )κ(R). Since the i-th row of Πℛ _R is v⊤v , this gives ‖ℛi⋅−(Πℛ)i⋅‖1=‖ri−v‖1≤2(d−1d)κ(ℛ). \|R_i·-( _R)_i·\|_1=\|r_i-v\|_1≤ 2 ( d-1d )κ(R). Taking the maximum over i∈⟦d⟧i∈ d , we conclude that ‖ℛ−Πℛ‖∞,1≤2(d−1d)κ(ℛ). \|R- _R\|_∞,1≤ 2 ( d-1d )κ(R). Finally, for the lower bound, notice that we can write v as v=1d∑k=1drk. v= 1d _k=1^dr_k. Therefore, for every i,k∈⟦d⟧i,k∈ d , ‖ri−rk‖1≤‖ri−v‖1+‖rk−v‖1≤2‖ℛ−Πℛ‖∞,1. \|r_i-r_k\|_1≤\|r_i-v\|_1+\|r_k-v\|_1≤ 2\|R- _R\|_∞,1. Using again (2.3), we conclude that κ(ℛ)≤‖ℛ−Πℛ‖∞,1κ(R)≤\|R- _R\|_∞,1. ∎ Theorem 2.3 gives the rollout column scores a structural interpretation. When the Dobrushin coefficient is small, the rows of ℛR are close to their mean profile v=d−1ℛ⊤v=d^-1R 1, whose components satisfy vj=sj/dv_j=s_j/d. Thus, in the approximately rank-one regime, the column scores encode the common distribution toward which the rows of the rollout operator homogenize. This statement concerns the representation of the rollout operator and does not imply that the profile v is concentrated or that dominant coordinates necessarily emerge. Any heterogeneity of v must be assessed separately from the contraction estimate. Remark A.3. We stress that small Dobrushin coefficient does not imply concentration of the propagation profile. For example, both ℛ1=d−1⊤R_1=d^-1 1 1 and ℛ2=e1⊤R_2= 1e_1 , where e1e_1 is the first vector in the canonical basis of ℝdR^d, satisfy κ(ℛi)=0κ(R_i)=0 and ℛi=ΠℛiR_i= _R_i. However, the associated profiles are respectively uniform and maximally concentrated. Theorem 2.3 therefore characterizes row homogenization, but not the heterogeneity of the common profile. A.1.1. Depth-dependent rank-one approximation through the product structure Theorem 2.3 characterizes the geometric structure of rollout operators in terms of their Dobrushin coefficient. We now examine how this coefficient is controlled by the product structure defining attention rollout. Throughout this subsection, fix x∈ℝdx ^d and write ℛ=R(x)=(L)(x)(L−1)(x)⋯(1)(x), =R(x)=A^(L)(x)A^(L-1)(x)·sA^(1)(x), where (ℓ)(x)≔(ℓ)(x)=12(I+A(ℓ)(x)). ^( )(x) A^( )(x)= 12 (I+A^( )(x) ). Here A(ℓ)(x)A^( )(x) is the attention matrix at layer ℓ . To simplify notation, the dependence on x is suppressed in the remainder of this subsection. Let ai(ℓ)a_i^( ) denote the i-th row of A(ℓ)A^( ). Then, the i-th row of (ℓ)A^( ) is i(ℓ)=12(ei+ai(ℓ)), a_i^( )= 12 (e_i+a_i^( ) ), where eie_i is the i-th canonical basis vector. Therefore, for every i≠ki≠ k, 12‖i(ℓ)−k(ℓ)‖1 12 \| a_i^( )- a_k^( ) \|_1 =14‖(ei−ek)+(ai(ℓ)−ak(ℓ))‖1 = 14 \|(e_i-e_k)+ (a_i^( )-a_k^( ) ) \|_1 ≤14‖ei−ek‖1+14‖ai(ℓ)−ak(ℓ)‖1=12+14‖ai(ℓ)−ak(ℓ)‖1. ≤ 14\|e_i-e_k\|_1+ 14 \|a_i^( )-a_k^( ) \|_1= 12+ 14 \|a_i^( )-a_k^( ) \|_1. Taking the maximum over all pairs and using (2.3) gives κ((ℓ))≤12(1+κ(A(ℓ))). κ (A^( ) )≤ 12 (1+κ (A^( ) ) ). (A.6) Moreover, in the standard self-attention mechanism considered in this work, the attention matrices A(ℓ)A^( ) are obtained by applying a row-wise softmax to finite attention logits. Consequently, they have strictly positive entries. This implies κ(A(ℓ))<1κ (A^( ) )<1 and, consequently, κ((ℓ))<1κ (A^( ) )<1. Thus, every rollout factor is strictly contractive. The residual identity component nevertheless limits the strength of this contraction, as shown by the following lemma. Lemma A.4. Let d≥2d≥ 2, let A∈ℝd×dA ^d× d be row-stochastic, and define P=12(I+A). P= 12(I+A). Then κ(P)≥d−22(d−1). κ(P)≥ d-22(d-1). Proof. Using Definition 2.2, we estimate the overlap between two distinct rows of P. For distinct i,k∈⟦d⟧i,k∈ d , we have Pii=12(1+Aii) and Pki=12Aki. P_i= 12(1+A_i) and P_ki= 12A_ki. Since Aki≤1≤1+AiiA_ki≤ 1≤ 1+A_i, it follows that minPii,Pki=12Aki and minPik,Pkk=12Aik. \P_i,P_ki\= 12A_ki and \P_ik,P_k\= 12A_ik. For j∉i,kj∉\i,k\, one therefore has minPij,Pkj=12minAij,Akj≤14(Aij+Akj). \P_ij,P_kj\= 12 \A_ij,A_kj\≤ 14(A_ij+A_kj). Consequently, ∑j=1dminPij,Pkj≤12(Aki+Aik)+14∑j∉i,k(Aij+Akj)=12+14(Aki+Aik−Aii−Akk). _j=1^d \P_ij,P_kj\≤ 12(A_ki+A_ik)+ 14 _j∉\i,k\(A_ij+A_kj)= 12+ 14 (A_ki+A_ik-A_i-A_k ). Averaging the last quantity over unordered pairs i<ki<k, we obtain (d2)−1∑i<k(Aki+Aik−Aii−Akk)=2(1−tr(A))d−1≤2d−1, d2^-1 _i<k (A_ki+A_ik-A_i-A_k )= 2 (1-tr(A) )d-1≤ 2d-1, where we used the row-stochasticity of A and tr(A)≥0tr(A)≥ 0. Hence there exists at least one pair i≠ki≠ k such that ∑j=1dminPij,Pkj≤12+12(d−1)=d2(d−1). _j=1^d \P_ij,P_kj\≤ 12+ 12(d-1)= d2(d-1). It follows that κ(P)≥1−d2(d−1)=d−22(d−1). κ(P)≥ 1- d2(d-1)= d-22(d-1). ∎ Lemma A.4 provides a dimension-dependent hard floor for the Dobrushin coefficient of an individual rollout factor. In particular, d−22(d−1)=12−12(d−1)→12 as d→+∞. d-22(d-1)= 12- 12(d-1)→ 12 as d→+∞. Thus, because one half of each rollout factor is reserved for the identity, a single layer cannot contract row oscillations by substantially more than a factor of one half. Finally, applying submultiplicativity to the local product above and using (A.6) (see, for instance, [33, Section 4.3]), we obtain κ(ℛ)≤∏ℓ=1Lκ((ℓ))≤∏ℓ=1L1+κ(A(ℓ))2=2−L∏ℓ=1L(1+κ(A(ℓ))). κ(R)≤ _ =1^Lκ (A^( ) )≤ _ =1^L 1+κ(A^( ))2=2^-L _ =1^L (1+κ(A^( )) ). (A.7) Estimate (A.7) separates the respective roles of residual averaging and of the layerwise attention patterns. The factor 2−L2^-L originates from the residual averaging in each rollout factor, while the terms 1+κ(A(ℓ))1+κ(A^( )) quantify how the attention matrices modulate the resulting bound. Since strict positivity only guarantees κ(A(ℓ))<1κ(A^( ))<1 separately at each layer, (A.7) does not by itself imply that κ(ℛ)→0κ(R)→ 0 as L→+∞L→+∞. Such a conclusion would follow, for instance, under the uniform condition supℓ≥1κ(A(ℓ))≤δ<1 _ ≥ 1κ(A^( ))≤δ<1, in which case κ(ℛ)≤(1+δ2)L→0. κ(R)≤ ( 1+δ2 )^L→ 0. In the absence of such a condition, the depth-dependent decrease of the Dobrushin coefficient remains an empirical property of the architectures studied in Section 4. Whenever the small-coefficient regime is reached, Theorem 2.3 guarantees that the rollout operator is close to a rank-one stochastic matrix whose propagation profile is determined by the realized attention matrices. Remark A.5. The product estimate (A.7) applies to the local rollout operator (2.1) which is by definition a product of attention propagation matrices. The global rollout operator R in (2.2) is instead defined as the dataset average and is generally not representable as a product of layerwise rollout matrices. Nevertheless, R is row-stochastic because it is the average of row-stochastic matrices. Consequently, Theorem 2.3 applies directly to R, independently of any product representation. Furthermore, by [17, Theorem 5.1], the Dobrushin coefficient coincides with the operator norm induced by Hopf’s oscillation seminorm. Hence, by the convexity of operator norms, κ()=κ(1N∑n=1NR(xn))≤1N∑n=1Nκ(R(xn)). κ( R)=κ ( 1N _n=1^NR(x_n) )≤ 1N _n=1^Nκ(R(x_n)). Combining this inequality with (A.7) yields κ()≤2−LN∑n=1N∏ℓ=1L(1+κ(A(ℓ)(xn))). κ( R)≤ 2^-LN _n=1^N _ =1^L (1+κ(A^( )(x_n)) ). (A.8) A.2. A simplified analytical relation between rollout and SHAP To make the analytical statement unambiguous, we work with exact interventional Shapley values. Let D≔⟦d⟧D d , let Z∼μZ μ be a reference random vector, and, for S⊆DS D and x∈ℝdx ^d, define the hybrid vector (xS⊕ZD∖S)k≔xk,k∈S,Zk,k∉S. (x_S Z_D S )_k casesx_k,&k∈ S,\\ Z_k,&k∉ S. cases The interventional coalition game is Vfint(S;x)≔[f(xS⊕ZD∖S)]. V_f^int(S;x) \! [f (x_S Z_D S ) ]. The corresponding interventional Shapley value [35] is ϕf,jint(x)≔∑S⊆D∖j|S|!(d−|S|−1)!d!(Vfint(S∪j;x)−Vfint(S;x)). _f,j^int(x) _S D \j\ |S|!\,(d-|S|-1)!d! (V_f^int(S∪\j\;x)-V_f^int(S;x) ). (A.9) Whenever the relevant expectations exist, efficiency gives f(x)=[f(Z)]+∑j=1dϕf,jint(x). f(x)=E[f(Z)]+ _j=1^d _f,j^int(x). This game differs in general from conditional SHAP, which is based on Vfcond(S;x)≔[f(Z)∣ZS=xS], V_f^cond(S;x) \! [f(Z) Z_S=x_S ], whose associated conditional Shapley values ϕf,jcond(x) _f,j^cond(x) are obtained by applying to it the Shapley-value construction (A.9). If the components of Z are independent, the interventional and conditional games coincide, up to the usual almost-everywhere qualification. The independence assumption in Proposition A.6 is imposed for this reason and should be regarded as a restrictive sufficient condition. To study how propagation- and attribution-based quantities can become related in a simplified setting, we fix a rollout matrix ℛ∈ℝd×dR ^d× d that is independent of the input x, and consider f(x)=g(x)+e(x),g(x)=β⊤ℛx,β∈ℝd. f(x)=g(x)+e(x), g(x)=β Rx, β ^d. (A.10) This decomposition is an analytical device, not a surrogate representation of the complete Transformer. Expanding g gives g(x)=∑j=1dγjxj,γj=∑i=1dβiℛij. g(x)= _j=1^d _jx_j, _j= _i=1^d _iR_ij. The following proposition supplies only a one-sided bound. It does not imply agreement between rollout and SHAP rankings. Proposition A.6. Let ℛR be a fixed, input-independent, rollout operator, let β∈ℝdβ ^d, and let f=g+ef=g+e be as in (A.10). Define β+≔maxi∈⟦d⟧|βi|,sj≔∑i=1dℛij. β^+ _i∈ d | _i|, s_j _i=1^dR_ij. Let X∼μX μ have independent components satisfying [Xj]=0E[X_j]=0 and [|Xj|]<+∞E[|X_j|]<+∞ for all j∈⟦d⟧j∈ d , and assume that e is globally bounded on ℝdR^d. Write ‖e‖∞≔supx∈ℝd|e(x)|<+∞. \|e\|_∞ _x ^d|e(x)|<+∞. Define the exact interventional SHAP importance by If,jint≔[|ϕf,jint(X)|]. I_f,j^int \! [ | _f,j^int(X) | ]. Then If,jint≤β+sj[|Xj|]+2‖e‖∞. I_f,j^int≤β^+s_j\,E[|X_j|]+2\|e\|_∞. (A.11) Under the stated independence assumption, the same bound holds for exact conditional SHAP values. Proof. The interventional Shapley operator is linear, so ϕf,jint(x)=ϕg,jint(x)+ϕe,jint(x). _f,j^int(x)= _g,j^int(x)+ _e,j^int(x). For every coalition S and every x∈ℝdx ^d, |Veint(S;x)|≤‖e‖∞ |V_e^int(S;x) |≤\|e\|_∞. The Shapley weights in (A.9) are nonnegative and sum to one. Hence |ϕe,jint(x)|≤2‖e‖∞ | _e,j^int(x) |≤ 2\|e\|_∞ for every x∈ℝdx ^d. For the linear function g, Vgint(S;x)=∑k∈Sγkxk+∑k∉Sγk[Xk]. V_g^int(S;x)= _k∈ S _kx_k+ _k∉ S _kE[X_k]. Every marginal contribution of feature j is therefore γj(xj−[Xj])=γjxj, _j (x_j-E[X_j] )= _jx_j, and consequently ϕg,jint(x)=γjxj _g,j^int(x)= _jx_j. Since rollout matrices are nonnegative, |γj| | _j| =|∑i=1dβiℛij|≤∑i=1d|βi|ℛij≤β+sj. = | _i=1^d _iR_ij |≤ _i=1^d| _i|R_ij≤β^+s_j. Thus, for every x∈ℝdx ^d, |ϕf,jint(x)|≤β+sj|xj|+2‖e‖∞. | _f,j^int(x) |≤β^+s_j|x_j|+2\|e\|_∞. Taking expectations proves (A.11). ∎ Remark A.7. Nonnegativity of ℛR is essential for the bound in terms of the unsigned column mass sjs_j. Row-stochasticity is not otherwise used in this proposition. For an arbitrary signed matrix B, the same argument requires the absolute column mass tj≔∑i=1d|Bij|, t_j _i=1^d|B_ij|, and yields a bound with tjt_j in place of sjs_j. The proposition therefore does not extend verbatim to arbitrary signed matrices while retaining the original score sjs_j. The independence assumption is substantive and is generally not satisfied by correlated biomedical covariates. Centering can enforce [Xj]=0E[X_j]=0, but centering, scaling, whitening, PCA, and decorrelation do not in general guarantee independence. Proposition A.6 should therefore be read only as an illustrative sufficient-condition result; it does not theoretically explain or validate the empirical agreement between rollout rankings and GradientExplainer rankings. Finally, because d is finite, the first-moment assumptions imply that M≔maxj∈⟦d⟧[|Xj|]<+∞, M _j∈ d E[|X_j|]<+∞, and hence If,jint≤β+Msj+2‖e‖∞. I_f,j^int≤β^+Ms_j+2\|e\|_∞. Appendix B Biomedical context for selected highly ranked variables This appendix places selected variables from the representative 11-layer ranking in biomedical context. It is a literature cross-check, not a validation of the ranking, a biomarker-discovery analysis, or evidence of causal physiological mechanisms. The discussion is based on the representative realization presented in Figure 4.6; the rankings vary across training realizations (Section 4.1.4). Several selected variables can be related to established biomedical literature: • Serum albumin has been studied as a marker of nutritional and inflammatory status and is associated with health outcomes in older adults [8]. • LDL cholesterol is an established causal risk factor for atherosclerotic cardiovascular disease [15]. • Glucose is central to glycemic status. The metabolite 1,5-AG is used as a marker of short-term glycemic excursions [9, 13]. • 3-hydroxyisobutyric acid has been reported to promote vascular fatty-acid transport and to contribute to insulin resistance in experimental models [26]. • Erythrocyte sedimentation rate is a nonspecific clinical marker of inflammation. Its appearance is compatible with, but does not establish, a link to the literature on chronic low-grade inflammation in aging [16]. The appearance of these highly ranked variables is compatible with prior biomedical literature, but it does not rule out confounding, redundancy, or spurious association. The partial overlap between attribution-based and propagation-based rankings is therefore treated only as a hypothesis-generating observation. The rollout operator should not be interpreted as reconstructing biochemical pathways or causal physiological mechanisms. The literature comparison supplies a qualitative plausibility check for this particular realization, not an independent validation of the framework or of any biomarker effect. References [1] S. Abnar and W. Zuidema (2020) Quantifying attention flow in transformers. In Proceedings of the 58th annual meeting of the association for computational linguistics, p. 4190–4197. Cited by: §1, §2.2. [2] G. Alain and Y. Bengio (2016) Understanding intermediate layers using linear classifier probes. arXiv preprint arXiv:1610.01644. Cited by: §1.1. [3] A. Alcalde, G. Fantuzzi, and E. Zuazua (2025) Clustering in pure-attention hardmax transformers and its role in sentiment analysis. SIAM J. Math. Data Sci. 7 (3), p. 1367–1393. Cited by: §2.2. [4] A. Álvarez-López, B. Geshkovski, and D. Ruiz-Balet (2026) Perceptrons and localization of attention’s mean-field landscape. arXiv preprint arXiv:2601.21366. Cited by: §2.2. [5] M. Arunika, S. Saranya, S. Charulekha, S. Kabilarajan, and G. Kesavan (2024) A survey on explainable AI using machine learning algorithms SHAP and LIME. In 2024 15th International Conference on Computing Communication and Networking Technologies (ICCCNT), p. 1–6. Cited by: §1. [6] A. Bibal, R. Cardon, D. Alfter, R. Wilkens, X. Wang, T. François, and P. Watrin (2022) Is attention explanation? An introduction to the debate. In Proceedings of the 60th Annual Meeting of the Association for Computational Linguistics (volume 1: long papers), p. 3889–3900. Cited by: §1.1. [7] U. Biccari (2026) Interior interpretability with attention rollout: contraction and propagation profiles in transformers. Note: https://github.com/umbertoBiccari/Attn_interpretability Cited by: §3.1, §3.1, §4, Data availability. [8] S. Cabrerizo, D. Cuadras, F. Gomez-Busto, I. Artaza-Artabe, F. Marín-Ciancas, and V. Malafarina (2015) Serum albumin and health in older people: review and meta analysis. Maturitas 81 (1), p. 17–27. Cited by: 1st item. [9] A. Ceriello, L. Monnier, and D. Owens (2019) Glycaemic variability in diabetes: clinical and therapeutic implications. Lancet Diabetes Endocrinol. 7 (3), p. 221–230. Cited by: 3rd item. [10] H. Chefer, S. G. Gur, and L. Wolf (2021) Transformer interpretability beyond attention visualization. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, p. 782–791. Cited by: §1.1. [11] R. L. Dobrushin (1956) Central limit theorem for nonstationary Markov chains. I. Theory Probab. Appl. 1 (1), p. 65–80. Cited by: §A.1, item 2, §2.2, §2.2. [12] Y. Dong, J. Cordonnier, and A. Loukas (2021) Attention is not all you need: Pure attention loses rank doubly exponentially with depth. In International Conference on Machine Learning, p. 2793–2803. Cited by: §2.2. [13] K. M. Dungan, J. B. Buse, J. Largay, M. M. Kelly, E. A. Button, S. Kato, and S. Wittlin (2006) 1,5-anhydroglucitol and postprandial hyperglycemia as measured by continuous glucose monitoring system in moderately controlled patients with diabetes. Diabetes Care 29 (6), p. 1214–1219. Cited by: 3rd item. [14] N. Elhage, T. Hume, C. Olsson, N. Schiefer, T. Henighan, S. Kravec, Z. Hatfield-Dodds, R. Lasenby, D. Drain, C. Chen, et al. (2022) Toy models of superposition. arXiv preprint arXiv:2209.10652. Cited by: §1.1. [15] B. A. Ference, H. N. Ginsberg, I. Graham, K. K. Ray, C. J. Packard, E. Bruckert, R. A. Hegele, R. M. Krauss, F. J. Raal, H. Schunkert, et al. (2017) Low-density lipoproteins cause atherosclerotic cardiovascular disease. 1. Evidence from genetic, epidemiologic, and clinical studies. A consensus statement from the European Atherosclerosis Society Consensus Panel. Europ. Heart J. 38 (32), p. 2459–2472. Cited by: 2nd item. [16] L. Ferrucci and E. Fabbri (2018) Inflammageing: chronic inflammation in ageing, cardiovascular disease, and frailty. Nat. Rev. Cardiol. 15 (9), p. 505–522. Cited by: 5th item. [17] S. Gaubert and Z. Qu (2015) Dobrushin’s ergodicity coefficient for Markov operators on cones. Int. Eq. Operator Theo. 81 (1), p. 127–150. Cited by: Remark A.5. [18] B. Geshkovski, C. Letrouit, Y. Polyanskiy, and P. Rigollet (2023) The emergence of clusters in self-attention dynamics. Adv. Neur. Inf. Process. Syst. 36, p. 57026–57037. Cited by: §2.2. [19] B. Geshkovski, C. Letrouit, Y. Polyanskiy, and P. Rigollet (2025) A mathematical perspective on transformers. Bull. Amer. Math. Soc. 62 (3), p. 427–479. Cited by: §2.2. [20] C. Grimsley, E. Mayfield, and J. R. Bursten (2020) Why attention is not explanation: surgical intervention and causal reasoning about neural models. In Proceedings of the Twelfth Language Resources and Evaluation Conference, p. 1780–1790. Cited by: §1.1. [21] G. Hannum, J. Guinney, L. Zhao, L. Zhang, G. Hughes, S. Sadda, B. Klotzle, M. Bibikova, J. Fan, Y. Gao, et al. (2013) Genome-wide methylation profiles reveal quantitative views of human aging rates. Mol. Cell 49 (2), p. 359–367. Cited by: §1. [22] J. Hewitt and P. Liang (2019) Designing and interpreting probes with control tasks. In Proceedings of the 2019 conference on empirical methods in natural language processing and the 9th international joint conference on natural language processing (emnlp-ijcnlp), p. 2733–2743. Cited by: §1.1. [23] S. Horvath (2013) DNA methylation age of human tissues and cell types. Genome Biol. 14 (10), p. 3156. Cited by: §1. [24] A. Ibáñez de Opakua, R. Conde, A. de Diego, M. Bizkarguenaga, N. Embade, S. Lu, J. Mato, and O. Millet (2025) Metabolomic-based aging clocks. npj Metabolic Health and Disease 3 (1), p. 35. Cited by: §1, §3.1. [25] S. Jain and B. C. Wallace (2019) Attention is not explanation. In Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers), p. 3543–3556. Cited by: §1.1. [26] C. Jang, S. F. Oh, S. Wada, G. C. Rowe, L. Liu, M. C. Chan, J. Rhee, A. Hoshino, B. Kim, A. Ibrahim, et al. (2016) A branched-chain amino acid metabolite drives vascular fatty acid transport and causes insulin resistance. Nature Med. 22 (4), p. 421–426. Cited by: 4th item. [27] S. M. Lundberg and S. Lee (2017) A unified approach to interpreting model predictions. Adv. Neur. Inf. Process. Syst. 30. Cited by: §1, §3.3. [28] C. Molnar (2020) Interpretable machine learning. Lulu. com. Cited by: §1. [29] E. Mosca, F. Szigeti, S. Tragianni, D. Gallagher, and G. Groh (2022) SHAP-based explanation methods: a review for NLP interpretability. In Proceedings of the 29th international conference on computational linguistics, p. 4593–4603. Cited by: §1. [30] N. Nanda, L. Chan, T. Lieberum, J. Smith, and J. Steinhardt (2023) Progress measures for grokking via mechanistic interpretability. arXiv preprint arXiv:2301.05217. Cited by: §1.1. [31] C. Olah, N. Cammarata, L. Schubert, G. Goh, M. Petrov, and S. Carter (2020) Zoom in: An introduction to circuits. Distill 5 (3), p. e00024–001. Cited by: §1.1. [32] T. Räuker, A. Ho, S. Casper, and D. Hadfield-Menell (2023) Toward transparent ai: A survey on interpreting the inner structures of deep neural networks. In 2023 ieee conference on secure and trustworthy machine learning (satml), p. 464–483. Cited by: §1.1. [33] E. Seneta (2006) Non-negative matrices and Markov chains. Springer Science & Business Media. Cited by: §A.1.1, §A.1, §2.2, §2.2. [34] S. Serrano and N. A. Smith (2019) Is attention interpretable?. In Proceedings of the 57th annual meeting of the association for computational linguistics, p. 2931–2951. Cited by: §1.1. [35] L. S. Shapley (1953) A value for n-person games. In Contributions to the Theory of Games, Volume I, H. W. Kuhn and A. W. Tucker (Eds.), Annals of Mathematics Studies, Vol. 28, p. 307–317. Cited by: §A.2. [36] K. Simonyan, A. Vedaldi, and A. Zisserman (2013) Deep inside convolutional networks: Visualising image classification models and saliency maps. arXiv preprint arXiv:1312.6034. Cited by: §1. [37] M. Sundararajan, A. Taly, and Q. Yan (2017) Axiomatic attribution for deep networks. In International conference on machine learning, p. 3319–3328. Cited by: §1. [38] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin (2017) Attention is all you need. Adv. Neur. Inf. Process. Syst. 30. Cited by: §2.1. [39] Y. Zhang, P. Tiňo, A. Leonardis, and K. Tang (2021) A survey on neural network interpretability. IEEE Trans. Emerg. Topics Comput. Intell. 5 (5), p. 726–742. Cited by: §1.