Self-information as a measure for identifying formulaic structures

The use of self-information stems from its foundational role in information theory, where the information content of an observation x is defined as

(3) $$ \begin{align} I(x) = -\log p(x), \end{align} $$

which quantifies how predictable or surprising an observation is. More predictable elements (those that occur frequently in structured, formulaic segments) have higher probabilities and thus lower − log p(x) values, while less predictable elements exhibit lower probabilities and higher − log p(x) values.

In the context of textual analysis, formulaic structures, by their nature, are repetitive and highly predictable, leading to higher − log p(x) values for formulaic segments. Our method directly measures the degree of structural predictability of each text segment, providing a granular, data-driven approach to distinguishing formulaic from non-formulaic patterns. Unlike existing information-based clustering methods, which rely on global variability measures, our approach evaluates the local predictability of individual elements, ensuring sensitivity to complex and varied linguistic structures.

Layman explanation: Self-information provides a measure of how surprising or expected a given unit of text is, based on how likely it is to appear. Frequent and structurally constrained expressions are more predictable and therefore contribute higher self-information values. By measuring this quantity locally for each segment, our method identifies areas of the text that exhibit heightened regularity or constraint, which are often associated with formulaic language.

Clustering approach

For a two-cluster problem, clustering methods aim to partition a dataset $X = \{x_1, x_2, \dots , x_n\}$ into two disjoint subsets, or clusters, based on a given similarity measure. Traditional hard clustering assigns each data point $x_i$ to exactly one of the two clusters $\mathcal {C}_1$ or $\mathcal {C}_2$ such that

(4) $$ \begin{align} x_i \in \mathcal{C}_j \quad \text{where} \quad \mathcal{C}_1 \cup \mathcal{C}_2 = X, \quad \mathcal{C}_1 \cap \mathcal{C}_2 = \emptyset. \end{align} $$

Soft clustering relaxes this constraint by allowing each data point $x_i$ to have a continuous membership weight $s_i \in [0,1]$ , representing its degree of association with one of the two clusters. Instead of assigning $x_i$ to a single cluster, we define

(5) $$ \begin{align} s_i \in [0,1], \quad s_i &= 0 \text{ for full membership in } \mathcal{C}_1,\nonumber\\s_i &= 1 \text{ for full membership in } \mathcal{C}_2. \end{align} $$

Thus, each data point has a soft assignment, where

$$\begin{align*}\begin{cases} x_i \in \mathcal{C}_1, & \text{if } s_i = 0, \\ x_i \in \mathcal{C}_2, & \text{if } s_i = 1, \\ x_i \in \mathcal{C}_1 \text{ and } \mathcal{C}_2 \text{ with weight } (1 - s_i, s_i), & \text{if } 0 < s_i < 1. \end{cases} \end{align*}$$

Our clustering approach employs a relaxation strategy within the optimization process, allowing for continuous rather than binary sample assignments. At each step of optimization, we iteratively update the probabilistic membership of each sample in the formulaic class, refining its classification based on structural predictability. This dynamic adjustment enhances the stability of the optimization, preventing indeterminate solutions while maintaining flexibility in high-dimensional or noisy datasets. However, this relaxation is not a final conclusion; rather, it is an intermediate step that ensures a more robust and well-conditioned optimization process, at the end of which each sample is determined by an integer label through a thresholding process.

Our method utilizes self-information to quantify structural predictability in text. Unlike hard clustering, which enforces discrete partitioning, we assign probabilistic memberships based on information content. Since formulaic and non-formulaic text layers differ in their information-theoretic properties – formulaic sections exhibiting higher predictability and lower self-information – our optimization scheme ensures a stable, information-guided classification.

Layman explanation: Rather than forcing each text unit of text into one cluster or the other, our method allows for partial or uncertain membership between the two. This flexibility is important in cases where text units may contain a mix of formulaic and non-formulaic features. During optimization, each segment is evaluated based on how structurally predictable it is. Highly predictable text units lean toward the formulaic cluster, while less predictable ones lean toward the non-formulaic cluster. This “soft” classification makes the method more stable and better suited for analyzing noisy or stylistically mixed texts.

Clustering formalism for discrete categorical data

To determine optimal cluster assignments, we minimize the cross-entropy between the empirical data distribution $\mathcal {P}_{\text {data}}(x)$ and the model distribution $\mathcal {P}_{\text {model}}(x)$ . The empirical distribution $\mathcal {P}_{\text {data}}(x)$ represents the observed frequencies of different textual elements in the dataset, while $\mathcal {P}_{\text {model}}(x)$ is derived from the soft clustering assignments, representing the likelihood of each sample belonging to a given cluster. Minimizing cross-entropy ensures that $\mathcal {P}_{\text {model}}(x)$ closely approximates the observed data, allowing the model to generalize the underlying structure while maintaining probabilistic flexibility.

Layman explanation: The goal is to ensure that the feature distributions within each cluster align well with the actual data. Cross-entropy serves as a standard measure of divergence between observed and modeled probabilities. Minimizing it improves the model’s ability to describe how features are distributed in the data.

Formally, cross-entropy is defined as

(6) $$ \begin{align} H(\mathcal{P}_{\text{data}}, \mathcal{P}_{\text{model}}) &= - \mathbb{E}_{\mathcal{P}_{\text{data}}} \left[ \log \mathcal{P}_{\text{model}}(x) \right]\nonumber\\&= - \sum_{x} \mathcal{P}_{\text{data}}(x) \log \mathcal{P}_{\text{model}}(x). \end{align} $$

We adopt a Bernoulli distribution-based model to characterize cluster assignments, as it provides a simple yet effective representation of feature occurrences in categorical data. This assumption is particularly relevant for textual analysis, where features – such as word occurrences or syntactic markers – tend to follow distinct distributions in formulaic and non-formulaic clusters. Given a soft clustering assignment, we estimate the model probability of a given word empirically by computing its weighted (i.e., soft) frequency within a cluster.

Layman explanation: The Bernoulli model treats each feature as either present or absent in a given segment of text. This is well-suited to binary or count-based feature representations. We estimate the frequency of each feature within a cluster by aggregating across samples, weighted by how strongly they belong to the cluster.

Specifically, the soft probability of the jth word appearing in a given cluster is estimated as

(7) $$ \begin{align} f(w_j | \vec{s}) = \frac{\sum_{i=1}^{n} s_i x_{ij}}{\sum_{i=1}^{n} s_i}, \quad p(w_j | \vec{s}) = \frac{f(w_j | \vec{s})}{\sum_{w' \in \chi} f(w' | \vec{s})}, \end{align} $$

where $\chi $ is the set of unique categorical features satisfying $w \in \chi $ , and $p(w | \vec {s})$ represents the probability of category w given a weight sequence $\vec {s}$ . The variable $x_i$ denotes the feature vector of the ith sample. Under the Bernoulli assumption, $x_{ij} = 1$ if the jth feature is present in the sample and $x_{ij} = 0$ otherwise. In the cumulative case, $x_{ij}$ instead represents the number of times the jth feature appears in the sample. This formulation allows us to compute empirical word distributions within each cluster while incorporating soft clustering assignments as weighting factors.

Layman explanation: This step produces a probability distribution over features for each cluster. These distributions capture what makes each cluster statistically distinct, based on how often specific features appear in text segments that are likely to belong to that cluster.

With these estimated word probabilities, we can compute the likelihood of a complete message (e.g., a verse) given a cluster. Under the Bernoulli assumption, the probability of a sample $x_i$ conditioned on its cluster assignment is

(8) $$ \begin{align}\kern-14pt P(x_i | \vec{s}) = \prod_{j=1}^{m} p(w_j | \vec{s})^{x_{ij}} (1 - p(w_j | \vec{s}))^{1 - x_{ij}}, \end{align} $$

(9) $$ \begin{align} P(x_i | 1 - \vec{s}) = \prod_{j=1}^{m} p(w_j | 1 - \vec{s})^{x_{ij}} (1 - p(w_j | 1 - \vec{s}))^{1 - x_{ij}}. \end{align} $$

Here, $p(w_j | \vec {s})$ and $p(w_j | 1 - \vec {s})$ represent the soft Bernoulli probabilities of the jth feature being present in each cluster. These probabilities are derived directly from the empirical word distributions calculated in Eq. (7).

Layman explanation: These expressions compute how likely a full set of features is for a given sample under each cluster model. If a text segment contains many features that are common in one cluster but rare in the other, the value will reflect that asymmetry.

Given this formulation, our score function naturally takes the form of the negative log-likelihood under a Bernoulli mixture model:

(10) $$ \begin{align} S(\vec{s}, X) = -\sum_{i=1}^{n} \left( s_i \cdot \log P(x_i | \vec{s}) + (1 - s_i) \cdot \log P(x_i | 1 - \vec{s}) \right), \end{align} $$

where $\vec {s} \in \mathbb {R}^n$ is the weight vector representing the soft assignment of samples to clusters, and $X \in \mathbb {R}^{n \times d}$ is the binary feature matrix of n samples with d features.

Layman explanation: The score function evaluates how well the current cluster assignments explain the data. The algorithm updates these assignments to minimize the score, which corresponds to improving the model’s alignment with the observed feature patterns.

To accommodate non-binary data, we extend the Bernoulli model to a Binomial formulation, where the number of trials per sample is determined by the total feature occurrences in that sample:

(11) $$ \begin{align}\kern-30pt P(x_i | \vec{s}) = \prod_{j=1}^{m} \binom{n_i}{x_{ij}} p(w_j | \vec{s})^{x_{ij}} (1 - p(w_j | \vec{s}))^{n_i - x_{ij}}, \end{align} $$

(12) $$ \begin{align} P(x_i | 1 - \vec{s}) = \prod_{j=1}^{m} \binom{n_i}{x_{ij}} p(w_j | 1 - \vec{s})^{x_{ij}} (1 - p(w_j | 1 - \vec{s}))^{n_i - x_{ij}}, \end{align} $$

where $n_i = \sum _{j=1}^{m} X_{ij}$ represents the total count of present features in sample $X_i$ . This ensures that the likelihood formulation remains valid when dealing with count-based or real-valued feature representations.

Layman explanation: When features can appear multiple times in a sample, such as when counting words as many times as they appear in a text unit, rather than whether the word appears at all or not, as in the binary case, the Binomial model accounts for these counts. This extension preserves the logic of the model while improving its fit to data with elements that are counted repeatedly.

Finally, the soft weighted log-likelihood of the dataset under the Binomial mixture model, which serves as the basis for our clustering optimization, is given by

(13) $$ \begin{align} \log P(x_i | \vec{s}) = \sum_{i=1}^{n} \sum_{j=1}^{m} \left[ \log \binom{k_j}{x_{ij}} + x_{ij} \log p(w_j | \vec{s}) + (k_j - x_{ij}) \log (1 - p(w_j | \vec{s})) \right]. \end{align} $$

This formulation connects back to the cross-entropy minimization framework introduced in Eq. (6), ensuring that cluster assignments reflect the probabilistic structure of the data.

Layman explanation: This final expression combines all the elements described above into a unified scoring objective. The clustering process then proceeds by adjusting the sample weights to find the configuration that most efficiently captures the structural differences in the data.

While this formulation provides a principled probabilistic basis for clustering, it does not explicitly account for competing structural signals that may emerge from different generative processes within the data. However, assessing the presence and dominance of such alternative signals is beyond the scope of this work, as our primary focus is on the probabilistic clustering framework itself. That said, the relative influence of entropy-based clustering versus clustering driven by feature distributions is inherently dependent on dataset-specific parameters, such as the baseline activation probability of features and the presence of structured formulaic dimensions. A heuristic analysis of this interplay is provided in Appendix C, where we illustrate how these factors affect clustering outcomes and discuss the need for future work to develop a predictive framework for determining when each signal will dominate.

Clustering formalism for continuous multivariate Gaussian data

In the case of continuous data, applying a Bernoulli model is not feasible, as such data does not inherently exhibit discrete binary properties. Instead, we adopt an information-based formalization, which provides a natural and mathematically elegant means of characterizing the structure of high-dimensional continuous distributions. Specifically, for data modeled as multivariate Gaussians, entropy offers a well-parameterized measure of dispersion and structure, making it a suitable alternative for defining soft cluster assignments.

Layman explanation: In the discrete case, each text segment is represented by binary indicators, capturing the presence or absence of individual features (e.g., words or n-grams) in a text unit. However, this representation does not extend to settings where features are continuous-valued, such as neural embeddings or normalized frequency vectors. In these cases, one must model how the values themselves vary across samples rather than simply whether a feature occurs. The multivariate Gaussian distribution provides a natural and mathematically grounded way to describe such continuous variation. Within this framework, entropy reflects how concentrated or dispersed the distribution is, offering a principled measure of structural regularity in continuous feature space.

Formally, in a d-dimensional space, the entropy of a multivariate Gaussian distribution is given by

(14) $$ \begin{align} H(X) = \frac{d}{2}(1 + \ln(2\pi)) + \frac{1}{2} \ln |\Sigma|, \end{align} $$

where $\Sigma $ is the covariance matrix, and $|\Sigma |$ denotes its determinant.

To incorporate weights in our soft clustering framework, we define a weighted covariance matrix:

(15) $$ \begin{align} \mu_w = \frac{\sum_{i=1}^n s_i x_i}{\sum_{i=1}^n s_i}, \quad \Sigma_w = \frac{\sum_{i=1}^n s_i (x_i - \mu_w)(x_i - \mu_w)^\top}{\sum_{i=1}^n s_i}, \end{align} $$

where $\{s_i\}$ are the weights assigned to samples. To enhance numerical stability in cases where the covariance matrix may be ill-conditioned, a regularization term is applied:

(16) $$ \begin{align} \Sigma_w \gets \Sigma_w + \epsilon I, \end{align} $$

where $\epsilon $ is a small positive constant, and I is the identity matrix of rank d. Using this regularized covariance, the total entropy of the weighted dataset is computed as

(17) $$ \begin{align} H(X) = \frac{d}{2}(1 + \ln(2\pi)) + \frac{1}{2} \ln |\Sigma_w|. \end{align} $$

This expression provides a global measure of the dataset’s uncertainty based on its covariance structure. However, in the context of soft clustering, we require a measure that quantifies the contribution of individual samples to the overall uncertainty. Since entropy characterizes unpredictability, a natural way to approximate the uncertainty associated with a specific sample $x_{i'}$ is by considering its likelihood under the dataset’s distribution.

For a multivariate Gaussian, the log-likelihood of a sample is

(18) $$ \begin{align} \log P(x_{i'}|\vec{s}) = -\frac{d}{2} \log(2\pi) - \frac{1}{2} \log |\Sigma_w| - \frac{1}{2} (x_{i'} - \mu_w)^\top \Sigma_w^{-1} (x_{i'} - \mu_w). \end{align} $$

This decomposition separates the likelihood into three terms: (i) a normalization constant, (ii) a term encoding the dataset’s covariance structure, and (iii) a term measuring a sample’s deviation from the mean under the Mahalanobis metric. Subsequently, the self-information of a sample is

(19) $$ \begin{align} I(x_{i'}|\vec{s}) = -\log p(x_{i'}) = \frac{d}{2} \log(2\pi) + \frac{1}{2} \log |\Sigma_w| + \frac{1}{2} D_M(x_{i'}|\vec{s}), \end{align} $$

where $D_M(x_{i'}|\vec {s}) = (x_{i'} - \mu _w)^\top \Sigma _w^{-1} (x_{i'} - \mu _w)$ is the Mahalanobis distance given the weighting by $\vec {s}$ (see Appendix A).

Layman explanation: A sample’s self-information increases the further it lies from the mean of the data, especially if it does so along directions where the data varies less. This distance, measured by the Mahalanobis metric, provides a principled way of quantifying how typical or atypical a given point is relative to the distribution of other points.

This equation directly relates a sample’s entropy contribution to its likelihood under the dataset’s distribution.

The Mahalanobis distance quantifies how far a sample deviates from the mean, making it a natural proxy for self-information. Samples close to $\mu _w$ have higher likelihoods and lower entropy, while those further away contribute greater uncertainty. This formulation provides a probabilistic approach to clustering, where samples are distinguished based on their deviation from the global covariance structure.

Score function for cluster identification: To identify structured (formulaic) clusters, we define a score function that optimizes cluster coherence while ensuring a sufficient number of samples:

(20) $$ \begin{align} S(I(X|\vec{s})) = \text{std}(I(X|\vec{s})) \cdot \text{avg}(I(X|\vec{s})) \cdot \frac{1}{||\vec{s}||}, \end{align} $$

where $I(X|\vec {s})$ is the soft distribution of self-information values (negative log-likelihood), and $\vec {s}$ is the weight vector. The terms $\text {std}$ and $\text {avg}$ represent the empirical standard deviation and mean of $I(X|\vec {s})$ , respectively.

Layman explanation: The score function favors clusters that are internally consistent (low variability), strongly predictable (low average self-information), and include enough samples to be meaningful. These three criteria together help isolate regions of the data that exhibit formulaic regularity.

This score combines the standard deviation, the mean and the inverse of the sample weights’ norm to quantify cluster coherence by measuring variability and predictability. The optimization minimizes variance while maintaining sufficient samples, distinguishing the formulaic cluster. Here, several differences from cross-entropy-based clustering (see the “Clustering Formalism for Discrete Categorical Data” section) are noteworthy: (1) Gaussian-based cross-entropy (e.g., Tabor and Spurek Reference Tabor and Spurek2014), relies on full covariance estimation, which becomes unstable when the sample size is small relative to the feature dimensionality, leading to singular covariance matrices. Our approach mitigates this issue by refining the score function to enhance the separability of the formulaic cluster. By focusing on isolating this cluster, we effectively reduce reliance on cross-cluster covariance terms, which are typically the most ill-conditioned. This ensures that covariance estimation is primarily constrained to a more compact and homogeneous subset of the data, where the conditioning of the covariance matrix remains more stable. (2) In Gaussian data, formulaic clusters have low self-information due to high predictability. Low entropy in these clusters corresponds to high likelihood, yielding lower negative log-likelihood. In contrast, in discrete categorical data, formulaic clusters are associated with higher self-information due to sparse but predictable occurrences of categories. Thus, Gaussian formulaic clusters exhibit low self-information, while discrete ones exhibit high self-information.

Optimization scheme

The algorithm optimizes weights $\vec {s} \in \mathbb {R}^n$ to isolate the self-information distribution of the structured or formulaic layer within the data. The optimization scheme is summarized in Algorithm 1. While our method is flexible and conceptually transparent, it is also more computationally intensive than conventional iterative clustering algorithms, such as k-means or expectation–maximization (EM). Those methods benefit from closed-form updates for centroids or covariance matrices, leading to fast convergence in low-dimensional, well-behaved data. In contrast, our approach involves nonlinear optimization over a continuous weight vector. This added cost reflects the complexity of modeling self-information and structural predictability, rather than geometric proximity or global variance. In scenarios where sample sizes are large and feature dimensionality is low, iterative methods may indeed be more efficient and comparably effective. However, this is rarely the case in historical or literary corpora, which tend to exhibit sparse, high-dimensional feature spaces and relatively small sample counts. It is precisely these conditions that motivated the design of an information-theoretic clustering approach tailored to such data.

Clustering benchmarking on categorical data

To evaluate the resolving power of our algorithm on one-hot encoded data, we design an experiment that generates synthetic datasets with features reflective of textual embeddings. These datasets simulate challenges, such as sparsity, feature dependencies and formulaic activation patterns.

The dataset is generated as follows: For each cluster, we create n samples with dimensionality d, where features are binary, representing the activation (1) or absence (0) of specific elements. One cluster, referred to as the uniform cluster, has an equal small activation probability p across all dimensions, resulting in sparse, randomly distributed features. The other cluster, referred to as the formulaic cluster, increases the activation probability for a subset of $d_{\mathrm {{form}}}$ dimensions by a bias factor f, while the remaining dimensions maintain the uniform activation probability, such that $p_{\text {form}} = p + f$ .

To introduce feature interdependencies within the formulaic cluster, we enforce correlations between pairs of dimensions. Specifically, for each sample in the formulaic cluster, m dimensions are randomly selected from the formulaic subset, and the activation state of one dimension is set to match the other. This process ensures that certain features are not independent but instead exhibit a structured co-occurrence pattern. Such interdependencies mimic the stylistic or structural dependencies commonly found in textual data, where the presence of one linguistic feature often implies the presence of another.

However, it is important to acknowledge the inherent difficulty of modeling textual data. Texts are complex, with nuanced structures and dependencies that cannot be fully captured by simplistic binary features or synthetic generation processes. Despite these limitations, the purpose of this benchmark is not to replicate textual data exactly but to demonstrate the ability of our algorithm to solve a difficult nonlinear problem – optimizing the distinction between sparse, discrete datasets that differ in entropy. This serves as proof of the algorithm’s robustness in handling challenging clustering problems.

We consider three nominal clustering algorithms that represent different clustering paradigms: (1) GMM, which assumes data is generated from a mixture of Gaussian distributions, making it suitable for probabilistic clustering in continuous spaces, (2) k-means (Hartigan and Wong Reference Hartigan and Wong1979), which partitions data into clusters based on minimizing intra-cluster variance, making it effective for Euclidean distance-based separation, and (3) DBSCAN (Ester et al. Reference Ester, Kriegel, Sander, Xiaowei, Simoudis, Han and Fayyad1996), which identifies clusters as dense regions of points separated by low-density areas, making it robust to arbitrary-shaped clusters but sensitive to parameter tuning. The results of this experiment, shown in Figure 1, demonstrate the superior performance of our algorithm in resolving formulaic clusters in sparse categorical encoded datasets whose sample-to-dimension ratio is small.

Figure 1.

Classification results for the benchmarking experiment of discrete categorical one-hot encoded discrete data described in the “Clustering Benchmarking on Categorical Data” section. The test datasets included 100 samples of (equally-sized) formulaic and non-formulaic classes, of 200 (top panel), 50 (middle panel), and 20 dimensions (bottom panel), with varying degrees of the probability of the base- (p) and formulaic- feature activation ( $p_{\text {form}}$ ), respectively, and the fraction of formulaic dimensions in the formulaic class. The colored areas represent one-standard-deviation intervals derived from 100 simulations.

The classification accuracy is measured using the Matthews correlation coefficient (MCC), given by

(21) $$ \begin{align} \text{MCC} = \frac{\text{TN} \times \text{TP} - \text{FN} \times \text{FP}}{\sqrt{(\text{TP} + \text{FP})(\text{TP} + \text{FN})(\text{FN} + \text{FP})(\text{TN} + \text{FN})}}, \end{align} $$

where T (F) stands for true (false) and P (N) stands for positive (negative), which are evaluations of some sequence against another. We subsequently normalize the original metric, spanning $\text {MCC}\in [-1, 1]$ to percent within the range of 50%–100%, where 50% suggests an arbitrary overlap between the output and real label sequences, and 100% suggests perfect overlap.

As shown in Figure 1, our method consistently outperforms traditional clustering approaches in the sparse low sample-size-to-dimensions ratio case by effectively distinguishing structured formulaic patterns even in highly sparse discrete encoded data.

Clustering benchmarking on multivariate Gaussian data

We perform a series of classification experiments to demonstrate the classification power of our algorithm compared to well-established classification routines designed for multivariate Gaussian data. These methods are: (1) Differential-Entropy (Davis and Dhillon Reference Davis, Dhillon, Schölkopf, Platt and Hofmann2006), (2) Gaussian Mixture EMFootnote ¹ (GMM), (3) Regularized EM (Reg. EM) (Houdouin, Ollila, and Pascal Reference Houdouin, Ollila and Pascal2023), and (4) CEC (Tabor and Spurek Reference Tabor and Spurek2014). Additionally, we test the resolving power of the $L_2$ norm distribution, incorporating it into a similar optimization scheme as described for the self-information distribution due to its relevance in high-dimensional settings where the concentration of measure ensures the norm becomes a distinguishing feature (but is independent of the number of samples).

We consider a baseline number of dimensions of $d = 50$ and evaluate performance across varying sample sizes per class. Each class is drawn from a multivariate Gaussian distribution, with both distributions centered at zero and defined by covariance matrices $\Sigma _1$ and $\Sigma _2$ . The eigenvalues of these covariance matrices are set to be uniformly $\lambda _1 = 10$ and $\lambda _2 = 30$ , respectively.

To introduce noise, we modify several eigenvalues of $\Sigma _1$ to match those of $\Sigma _2$ . The noise level is defined as the relative proportion of eigenvalues in $\Sigma _1$ replaced with corresponding eigenvalues of $\Sigma _2$ . In Figure 2, we present the results of these experiments, showcasing the classification performance of each method under varying sample sizes, numbers of dimensions, and noise levels.

Figure 2.

Classification results of the experiment described in the “Clustering Benchmarking on Multivariate Gaussian Data” section, for a varying number of sample sizes and numbers of dimensions of multivariate Gaussian classes of varying entropy. Upper panel: Varying sample sizes for $d = 50$ . Bottom panel: Varying sample sizes for $d = 10$ . The colored areas represent one-standard-deviation intervals, derived from 100 simulations.

We find that our classification approach benefits both from the sample size and the number of dimensions and exhibits a considerably weaker sensitivity to a small sample-size-to-dimensions ratio than other classification routines that leverage the multivariate Gaussian assumption. This is because estimating covariance matrices in high-dimensional, small-sample settings is ill-conditioned, leading to unstable classification. Our method circumvents this issue by leveraging auto-entropy distributions rather than direct covariance estimation. Specifically, as shown in Appendix A, self-information-based classification avoids reliance on explicit covariance inversion by computing the soft auto-entropy of individual samples with respect to the weighted dataset entropy.

This formulation improves numerical stability and mitigates singularity issues that arise when the sample size is small relative to the number of dimensions, as illustrated in Figure 2.

Experimental setup

Embedding

We consider a cumulative-one-hot-encoded representation of the corpus D, such that D ∈ ℝ ^n×χ , where n is the number of text units, χ is the set of all unique n-grams in the corpus and d _ij represents the number of occurrences of the jth feature in the ith text unit.

We consider a parameter space of embedding possibilities whose permutations we fully explore. Specifically, each book is embedded according to a combination of several parameters across a grid of all permutations thereof: n and ℓ, spanning the following ranges: n ∈{1, 2, 3, 4, 5} for word n-grams, and ℓ ∈{2, 3, 4, 6, 8, 10, 12, 14, 18, 22, 24, 26, 28}, which defines the number of (overlapping) consecutive verses over which features are aggregated. This windowing approach allows the model to capture local textual dependencies by smoothing over short-range fluctuations, helping to account for stylistic continuity while preserving meaningful segment-level distinctions. To further balance feature richness and sparsity, we restrict our analysis to the $f\in \{100, 300,500, \text {all} \}$ most frequent features in each embedding configuration. This ensures that selected features are statistically meaningful while preventing noise from rare or idiosyncratic terms that may introduce spurious clustering patterns. By systematically varying the feature set size, we evaluate the robustness of our method across different levels of vocabulary complexity.

This scope is motivated by the desire to exhaustively cover potential feature spaces while ensuring that the extracted features remain both meaningful and sufficiently frequent within the biblical corpus, thereby avoiding excessive sparsity due to the limited length of biblical texts (e.g., Antonia, Craig, and Elliott Reference Antonia, Craig and Elliott2014). Given the relatively small size of the corpus compared to modern datasets, this balance is particularly crucial in detecting meaningful linguistic patterns while maintaining statistical robustness.

We opt for traditional, count-based embeddings, operating under the Binomial score function (Eqs. 11 and 12), to ensure interpretability, particularly given the highly structured and constrained nature of biblical Hebrew. While modern deep learning techniques, such as transformer-based models, have revolutionized text representation in many languages (e.g., Shmidman et al. Reference Shmidman, Guedalia, Shmidman, Shmidman, Handel and Koppel2022), no robust pre-trained language model for biblical Hebrew currently exists. As a result, neural embeddings would require extensive domain-specific pertaining (e.g., Huertas-Tato, Martín, and Camacho Reference Huertas-Tato, Martín and Camacho2023), making their applicability to this study uncertain. By employing interpretable feature representations, we ensure that clustering results remain transparent and analyzable within the context of existing linguistic and source-critical scholarship.

Figure 3.

Clustering results for the book of Genesis across different parameter combinations, evaluated against expert annotations distinguishing between the main textual body and genealogical lists traditionally attributed to P. Results are shown for our cross-information-based clustering method (left) and k-means (right). Top panel: The 20 feature combinations that yield the highest MCC scores, indicating the strongest agreement with expert annotations. Bottom panel: Distribution of MCC scores across all parameter combinations, sorted into discrete performance intervals.

To further clarify the effect of parameter settings, we note that varying the n-gram size directly affects both sparsity and entropy in the feature space. Larger n-grams tend to emphasize rigid phraseology and reduce overlap across segments, making them useful for identifying formulaic repetitions. Smaller n-grams increase coverage but may conflate stylistically distinct phrases. Likewise, broader running windows ( $\ell $ ) smooth over local variation and can help reveal segmental regularities, while narrower windows are more sensitive to fine-grained shifts. Feature count thresholds (f) balance interpretability and statistical stability by filtering out rare or idiosyncratic features. These parameters jointly determine which aspects of formulaicity are foregrounded and how they align with the internal segmentation of the text.

Comparison with previous work

There is very limited prior work on unsupervised clustering approaches for exploring hypothesized divisions of biblical texts. Most previous studies have relied on supervised classification methods, where textual divisions were predefined, and models were trained using manually curated feature sets (Dershowitz et al. Reference Dershowitz, Akiva, Koppel and Dershowitz2015; Radday and Shore Reference Radday and Shore1985). These approaches often rely on cherry-picked linguistic or stylistic markers that align with existing scholarly expectations, making the claim of statistical significance somewhat circular. Since the feature selection process is influenced by prior assumptions about textual divisions, such methods do not provide independent verification of whether distinct literary strata emerge naturally from the data.

Figure 4.

Clustering results for the P/non-P partition in the book of Exodus, similar to Figure 3.

Figure 5.

Clustering results for the P/H partition in the book of Leviticus, similar to Figure 3.

By avoiding reliance on predefined features and instead allowing clusters to emerge based on intrinsic textual properties, we aim to assess whether computational methods can independently recover traditionally hypothesized partitions. Due to the lack of prior unsupervised studies, we resort to comparing our results with previous supervised approaches, evaluating whether our method aligns with existing classifications or suggests alternative structuring. In addition, we employ k-means clustering, previously utilized in unsupervised clustering of biblical corpora (e.g., Yoffe et al. Reference Yoffe, Bühler, Dershowitz, Römer, Piasetzky, Finkelstein and Sober2023), as a baseline method to provide a comparative reference for our information-based approach. We also performed this analysis using GMM clustering for completeness, whose results are nearly identical to those of k-means (similarly to what was observed on synthetic data; see Figure 1) and can be found in Appendix D.

Formulaic structure and parameter sensitivity in the P/H partition of Leviticus

Among the three case studies, the Leviticus P/H division shows the most variation in performance across parameter settings. In Genesis and Exodus, the contrasts between genealogical and narrative passages in Genesis, and between priestly and non-priestly material in Exodus, lead to many strong results clustered in the 85%–90% and 90%–95% MCC ranges. This suggests that the stylistic differences in those cases are captured consistently across a wide range of feature combinations. In contrast, the P/H division in Leviticus produces more extreme outcomes. It yields more high-performing results in the 90%–95% range than k-means-based partitions, but also many more low-performing results below 75%. This behavior stands in contrast to the k-means baseline, which produces fewer strong results but also fewer weak ones. These patterns suggest that the formulaic differences between P and H are more sensitive to how the text is segmented and represented. In the analysis that follows, we explore how different parameter settings lead to shifts in which passages are classified as more internally regular.

In Figure 6, we illustrate how varying parameter choices, specifically, using 2-grams versus 4-grams with different running window widths, affect the identification of formulaic patterns in the text. These differences in parameterization lead to both hypothesized classes (P and H) being classified as the formulaic cluster under different settings. In the first case, we apply the parameter combination that yields the highest MCC score of 94% (see Figure 5). Here, the formulaic cluster is characterized by 4-grams, such as (“I am the LORD thy God”) and { } (“[Speak to] the Israelites and say [to them]”), which are prosaic refrains strongly associated with H (Knohl Reference Knohl2007). Notably, in this setting, H is classified as the formulaic cluster, while P is identified as the non-formulaic one. There, the first 4-gram serves as a theological and moral refrain, reinforcing divine authority and the covenantal obligations of Israel. It often punctuates commandments, emphasizing that obedience is not merely a legal or social requirement but a direct response to God’s sovereignty (e.g., Lev 19:2, 19:4, 19:10), and the second introduces divine laws and instructions, repeatedly marking the transition between legislative units in Leviticus (e.g., Lev 17:2, 18:1, 19:1), which constitute the P-associated corpus. Notably, the identification of a running-window width of 28 is not coincidental; it corresponds to the average length of the priestly legislative units (Knohl Reference Knohl2007). This choice of parameters (n = 4, $\ell $ = 28, $f = \text {all}$ ) distinguishes the more formulaic nature of H (only in the case of 4-grams), from the P material, whose structural organization aligns with the larger 28-verse window.

In the case of 2-grams-based clustering results (third highest MCC score with $n = 2$ , $\ell = 20$ and $f = 300$ ), we demonstrate that the distinctive features of the formulaic cluster are explicitly associated with the priestly legislative units. These include recurring phrases such as {} (“The Priest”), which frequently appear in ritual contexts, delineating the priestly role in sacrificial procedures, purity regulations and atonement rites (e.g., Lev 4:5, 6:7, 16:32). Similarly, {} (“The Altar”) is central to the sacrificial system, often marking legal prescriptions concerning offerings and the sanctity of cultic spaces (e.g., Lev 1:5, 6:9, 8:15). These formulaic elements reflect the institutional and hierarchical focus of P, where priestly duties and sacrificial procedures are codified with precise and repetitive language. The strong association of these 2-grams with priestly legislative discourse underscores the structured, procedural nature of P, distinguishing it from the moral and theological refrains characteristic of H.

These results underscore that the identification of formulaic structure is not fixed to a single representation but emerges differently depending on the linguistic and structural scales chosen. In the case discussed here, distinct sets of formulaic features become salient under different n-gram and segmentation settings, yielding competing clusterings where either P or H is classified as more formulaic. This divergence does not undermine the method, but rather illustrates its strength: it allows researchers to test how formulaicity is distributed under different structural assumptions, exposing layered regularities that vary in grammatical, semantic, and discourse scale. The ability to generate high-confidence classifications under some parameterizations, and not others, suggests that formulaicity itself is not an absolute property, but a perspective-sensitive feature of textual organization.