An unsupervised information-theoretic approach to identifying formulaic clusters in textual data

Abstract

Texts, whether literary or historical, exhibit structural and stylistic patterns shaped by their purpose, authorship and cultural context. Formulaic texts, which are characterized by repetition and constrained expression, tend to differ in their information content (as defined by Shannon) compared to more dynamic compositions. Identifying such patterns in historical documents, particularly multi-author texts like the Hebrew Bible, provides insights into their origins, purpose and transmission. This study aims to identify formulaic clusters: sections exhibiting systematic repetition and structural constraints, by analyzing recurring phrases, syntactic structures and stylistic markers. However, distinguishing formulaic from non-formulaic elements in an unsupervised manner presents a computational challenge, especially in high-dimensional and sample-poor data sets where patterns must be inferred without predefined labels.

To address this, we develop an information-theoretic algorithm leveraging weighted self-information distributions to detect structured patterns in text. Our approach directly models variations in sample-wise self-information to identify formulaicity. By extending classical discrete self-information measures with a continuous formulation based on differential self-information in multivariate Gaussian distributions, our method remains applicable across different types of textual representations, including neural embeddings under Gaussian priors.

Applied to hypothesized authorial divisions in the Hebrew Bible, our approach successfully isolates stylistic layers, providing a quantitative framework for textual stratification. This method enhances our ability to analyze compositional patterns, offering deeper insights into the literary and cultural evolution of texts shaped by complex authorship and editorial processes.

Plain Language Summary

This work presents a method for identifying formulaic clusters – segments of text that rely heavily on repetition, structural regularity or stylistic constraint. Such clusters can be found in a wide range of contexts, including religious texts, technical manuals, oral traditions and institutional writing, where language is shaped by convention, genre or the need for consistency. Detecting these patterns can offer insight into how texts were composed, transmitted and shaped by cultural forces over time.

Our approach is grounded in information theory and focuses on the concept of self-information, which measures how predictable or unpredictable a text unit is with respect to the corpus with which it is associated. Formulaic segments tend to be more predictable and, therefore, carry higher self-information. By measuring the distribution of self-information across a corpus, we identify areas where language use is more constrained, indicating a higher degree of formulaicity.

Benchmarking against standard clustering methods on synthetic and real data confirms the method’s robustness, particularly in high-dimensional, low-sample contexts. It outperforms conventional models in isolating internally consistent structures, without requiring strong distributional assumptions.

We apply our method to three books of the Hebrew Bible: Genesis, Exodus and Leviticus. In Genesis, we focus on the distinction between genealogical lists and surrounding narrative material. The genealogical sections are known for their rigid, repetitive structure, and our method effectively isolates them as highly formulaic. In Exodus, we examine the division between Priestly (P) and non-Priestly texts, with the Priestly material comprising legal and cultic instructions. The method reliably separates these layers in line with traditional source-critical scholarship. In Leviticus, we test whether the method can distinguish between the Priestly source and the Holiness code (H), two stylistically distinct strata. The results show that subtle differences in recurring expressions and structural regularity allow for a clear separation between them under appropriate conditions.

This work offers a flexible, information-theoretic framework for analyzing the internal structure of texts. It is well suited to the study of layered, edited or composite documents, where patterns of repetition and constraint signal deeper organizational principles. The method is broadly applicable to literary, historical, legal or technical corpora, wherever formulaic language plays a role in shaping textual form.

Introduction

An intriguing challenge in textual analysis is distinguishing formulaic layers – those characterized by recurring patterns, motifs or stylistic conventions – from non-formulaic ones. In this context, a formulaic element is defined as a word or a sequence thereof, syntactic or grammatical structures, or thematic patterns that recur systematically within a corpus, either due to conventionalized linguistic usage, genre-specific constraints or cultural transmission mechanisms. Formulaic expressions can function as structural scaffolds within texts, shaping their composition and reflecting underlying conventions of specific literary traditions, genres or modes of authorship (e.g., Paquot and Granger 2012; Read and Nation 2008; Wood 2019). Beyond their structural significance, analyzing formulaic patterns provides insights into the dynamic interplay between creativity and constraint in language, illuminating the socio-cultural and cognitive forces that have influenced textual transmission and stylistic evolution (e.g., Jensen 1980; Magoun Jr 1953).

The importance of formulaicity is particularly pronounced in historical texts, where composition layers often reflect complex authorship and redaction processes. Many historical documents, created and transmitted over centuries, are products of collective memory and institutional influence. As such, these texts often contain formulaic structures that encode cultural norms, narrative conventions, and theological or political ideologies (e.g., Knohl 2007; Polak 2017). Understanding these layers can reveal how texts were constructed, transmitted and adapted over time (e.g., Coleman 2019; Gitay 1980; Polak 2006).

Among historical texts, the Hebrew Bible provides a compelling case study. As a composite text shaped by centuries of authorship, redaction and transmission, the Hebrew Bible embodies a rich interplay of formulaic and non-formulaic elements. Its stylistic diversity reflects contributions from different schools of authorship, genres and literary motifs, ranging from tightly structured legal codes to more fluid poetic and narrative traditions (Shectman and Baden 2009; Stipp 2017; Wärters 1976). Detection and analysis of formulaic patterns within the Hebrew Bible offer a pathway to better understanding its compositional history, shedding light on how diverse voices and traditions coalesced into the canonical form.

Understanding the distinction between formulaic and non-formulaic layers requires rigorous methodologies to quantify stylistic and structural patterns across literary feature sets. One promising metric for addressing these challenges is entropy, a concept from information theory that measures the unpredictability or randomness of a dataset. Initially introduced by Rudolf Clausius in 1865 as a measure of energy dispersal in a system, and subsequently formalized statistically by Shannon (1948), entropy has found applications in diverse fields, including thermodynamics, cryptography and data compression. In the context of textual analysis, entropy provides a framework for quantifying the variability and regularity inherent in linguistic structures. Formally, given a discrete probability distribution $P = \{p_1, p_2, \dots , p_n\}$ over n possible elements, Shannon entropy is defined as

(1) $$ \begin{align} H(P) = -\sum_{i=1}^{n} p_i \log p_i ,\end{align} $$

where $p_i$ represents the probability of the ith element. Higher entropy values indicate greater variability, while lower entropy values correspond to constrained, repetitive structures, making entropy a powerful tool for distinguishing formulaic clusters from more unpredictable distributions.

The relevance of entropy to textual analysis lies in its ability to capture the degree of predictability of literary or linguistic features within a sequence, such as a sentence, paragraph or chapter. Texts with high formulaicity, such as legal codes or liturgical compositions, are likely to exhibit low entropy, such that text units within them resemble one another, for example, in vocabulary or grammatical structure. Conversely, more creative or extemporaneous texts often exhibit higher entropy, reflecting their greater unpredictability. Information-based metrics, therefore, serve as a natural tool for distinguishing between formulaic and non-formulaic elements within texts (e.g., Church and Hanks 1990).

In textual contexts, entropy is closely related to the concept of perplexity – the exponential of the average of self-information of all features in a sequence, which is widely used in language modeling and computational linguistics. Perplexity measures how well a probabilistic model predicts a sequence of words, with lower perplexity indicating a better fit. Formally, perplexity is defined as the exponential of the average negative log probability of the observed sequence:

(2) $$ \begin{align} \text{Perplexity} = \exp\left(-\frac{1}{N} \sum_{i=1}^N \log P(w_i \mid w_{i-1}, w_{i-2}, \dots)\right), \end{align} $$

where $w_i$ represents the ith word in the sequence, and $P(w_i \mid w_{i-1}, w_{i-2}, \dots )$ denotes the conditional probability of $w_i$ given its preceding context. Perplexity has proven instrumental in analyzing textual corpora over various linguistic angles (e.g., Gamallo, Campos, and Alegria 2017; Huang and Hansen 2007; Klakow and Peters 2002). In the context of formulaic analysis, perplexity can reveal how tightly a text adheres to specific linguistic conventions (e.g., Kurzynski 2023).

The task of identifying formulaic layers is inherently complex, as clustering formulaic elements within a corpus involves evaluating an exponentially large set of possible ways to divide the text, as each partition represents a different possible grouping of elements that exhibit formulaicity. This combinatorial challenge necessitates algorithmic approaches that efficiently approximate optimal solutions. One common strategy involves modeling linguistic feature distributions using probabilistic frameworks, with Gaussian-based methods being among the most widely applied.

For instance, Gaussian mixture models (GMMs) (Vlassis and Likas 2002) assume that the dataset can be represented as a weighted sum of Gaussian components, enabling efficient clustering through iterative optimization of parameters, such as means, covariance matrices and mixing weights. Similarly, cross-entropy clustering (CEC) (Tabor and Spurek 2014) explicitly incorporates entropy as a criterion to partition data, refining clusters based on their statistical coherence. These methods, while effective in many contexts, rely on assumptions about the underlying data distribution, which may not always align with the structural and stylistic non-Gaussian variability inherent in textual corpora.

Recent advances in representation learning, including the contrastive language–image pretraining model (CLIP; Radford et al. 2021) and bidirectional encoder representations from transformers (BERT; Devlin et al. (2019)), have introduced powerful mechanisms for unsupervised clustering and alignment in high-dimensional semantic spaces. These models often optimize mutual information across paired or contextualized inputs, effectively capturing latent structural relationships without requiring explicit probabilistic assumptions. Building on these architectures, recent work has proposed information-theoretic evaluation metrics computed over embedding spaces to quantify diversity and structural complexity. One such class of methods includes entropy measures based on Rényi kernel entropy (RKE), which estimate entropy non-parametrically using similarity kernels to capture the dispersion of representations in high-dimensional spaces (Giraldo, Rao, and Principe 2014; Principe 2010). These approaches are increasingly used to assess representational diversity and clustering compactness in generative models and learned embeddings. Many of these entropy-based methods derive from or build upon neural embedding models such as variational autoencoders (VAEs), which assume that latent features follow a Gaussian distribution in a continuous space (Dilokthanakul et al. 2016; Li et al. 2020; Yang et al. 2019). VAEs excel in capturing the underlying structure of data by encoding inputs into compact latent representations governed by a Gaussian prior (Casale et al. 2018), which enables smooth interpolations and the generation of new data points resembling the original distribution.

While these approaches provide useful frameworks for clustering and representation learning, their reliance on Gaussian assumptions and metrics that require the estimation of covariance presents notable challenges. In textual data, where features often exhibit categorical, skewed or multimodal distributions, the Gaussian prior may not accurately capture structural variability. Moreover, many of these methods are highly sensitive to the sample-to-dimension ratio (Altonji and Segal 1996; Yao, Zheng, and Bai 2015). When the sample size is small relative to the feature space, covariance estimation becomes unstable, leading to poorly conditioned matrices or singularities (e.g., Ashurbekova et al. 2021). Regularization techniques, such as adding small positive constants to the diagonal of covariance matrices, can partially mitigate these issues but introduce biases that require careful tuning (Bickel and Levina 2008). This bias arises because regularization alters the estimated eigenvalues of the covariance matrix, effectively shrinking them toward a predetermined value, which can distort the true structure of the data and influence clustering or classification outcomes.

Crucially, while these techniques are effective in high-resource, general-purpose NLP settings, their applicability to low-resource or singular historical corpora remains limited. In domains such as ancient languages, where pretrained embeddings are unreliable or unavailable, and where interpretability is essential for source-critical and literary-historical inquiry, alternative approaches are required.

These challenges are particularly pronounced in textual analysis, where statistical significance often necessitates dividing texts into relatively large units, such as paragraphs or entire chapters. While this improves the reliability of individual measurements, it also limits the number of available samples. At the same time, textual data often exhibit high dimensionality due to feature-rich representations, such as frequent words, n-grams or syntactic constructions, leading to a severe imbalance between the number of samples and features. This issue is especially prevalent in stylometric analysis, where subtle stylistic differences are captured through sparse, high-dimensional feature spaces (Stamatatos 2009). Given these constraints, clustering methods that rely on covariance estimation or strong parametric assumptions struggle to produce stable and interpretable results.

To address these limitations, we introduce an entropy-driven soft clustering framework that identifies formulaic structures based on the log probability (log p) of data points within a given distribution. The choice of − log p as our core metric stems from fundamental principles in information theory: the information content of an observation x, given a probability distribution $p(x)$ , is defined as $-\log p(x)$ , known as its self-information. This quantity directly captures the degree of predictability of x – highly predictable elements (i.e., those frequently occurring in structured, formulaic segments) have higher probabilities and thus lower self-information values, whereas less structured, unpredictable elements exhibit lower probabilities and correspondingly higher self-information values.

By leveraging this information-theoretic perspective, our method provides a natural and rigorous approach to distinguishing formulaic from non-formulaic structures. Our approach evaluates the information content of each data point within its feature space, allowing for a more granular and structurally sensitive identification of formulaic patterns. This circumvents the instability issues associated with covariance estimation in high-dimensional spaces while ensuring robustness to sparsity and multimodal distributions. By directly quantifying the predictability of textual elements, our method provides a principled, data-driven alternative to traditional clustering approaches that struggle with complex linguistic structures.

Our framework integrates soft clustering weights (Peters et al. 2013), which adaptively re-weight samples based on their self-information contributions, mitigating the influence of sparse or noisy data points. This probabilistic approach ensures robust partitioning of text segments by capturing both localized structural patterns and broader stylistic trends. Unlike traditional clustering methods that assume specific data distributions, such as GMMs, our information-based approach remains effective across a wide range of feature distributions, given an adequate choice of the score function, including multimodal and categorical datasets. Furthermore, by avoiding direct covariance inversion, our method enhances numerical stability and reduces susceptibility to the curse of dimensionality. Through this paradigm, we offer a flexible, data-driven solution for identifying formulaic structures in textual corpora while maintaining robustness in challenging high-dimensional settings.

To accommodate different types of textual feature representations, we consider two complementary formulations of the method: one for discrete data and one for continuous data. The discrete formulation operates directly on binary or count-based features, avoiding any form of covariance estimation. In contrast, the continuous formulation, used for real-valued vector representations, such as embeddings, does involve estimating a soft covariance matrix. While this introduces some of the challenges seen in classical approaches, our method confines the estimation to weighted subsets of the data associated with candidate clusters exhibiting internal coherence. This restriction improves numerical stability by avoiding cross-cluster covariance terms, which are often the most ill-conditioned. Moreover, in both formulations, the objective function remains grounded in self-information, prioritizing structural predictability over explicit density modeling. This unifying information-theoretic perspective allows the method to remain robust across data modalities while still benefiting from statistical sensitivity to underlying stylistic regularities.

The article is structured as follows: In the “Soft Clustering Algorithm” section, we introduce the mathematical foundations of the information-based soft clustering framework, defining the self-information distribution and formulating an optimization scheme for distinguishing formulaic from non-formulaic textual components in the “Benchmarking” section, we benchmark the method against existing clustering techniques on synthetic datasets, assessing its ability to classify structured and unstructured elements under controlled conditions, in the “Application to Textual Data: The Priestly Source(s) in the Pentateuch” section, we apply the framework to the analysis of the biblical corpus, evaluating its capacity to differentiate hypothesized literary strata based on stylistic and formulaic properties, and in the “Discussion” section, we conclude.

Soft clustering algorithm

The following sections establish the mathematical foundations of the proposed information-based soft clustering framework and formulate an optimization scheme for distinguishing formulaic from non-formulaic text layers. In this framework, “soft clustering” refers to the relaxation techniques employed during optimization to avoid indeterminate solutions (e.g., Royer 2017). This approach follows established clustering methodologies, where constraints are relaxed to transform discrete, hard assignment problems into continuous optimization problems, making them more tractable. An example of this is semidefinite relaxation in graph clustering, which reformulates discrete clustering problems as continuous ones, enabling more efficient optimization. By utilizing such relaxation techniques, our framework can more effectively identify structured patterns in text while maintaining flexibility in handling ambiguity in cluster assignments.

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 2014), 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.

Benchmarking

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 1979), 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. 1996), 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 2006), (2) Gaussian Mixture EM ¹ (GMM), (3) Regularized EM (Reg. EM) (Houdouin, Ollila, and Pascal 2023), and (4) CEC (Tabor and Spurek 2014). 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.

Application to textual data: The priestly source(s) in the Pentateuch

The composition of the Pentateuch, or Torah, has been a longstanding subject of critical scholarship, with various theories proposed to account for its linguistic diversity, stylistic inconsistencies, narrative discrepancies, and theological tensions. One of the most widely discussed models is the Documentary Hypothesis, which suggests that the Torah is a composite text formed through the combination of multiple sources, each reflecting different historical, ideological, or scribal traditions (Holzinger 1893; Kuenen and Wicksteed 1886; Wellhausen 1885). Among the sources reconstructed by this framework are the Yahwist (J), the Elohist (E), the Deuteronomist (D), and the Priestly (P) source. While the precise delineation of these sources remains debated, with alternative models proposing alternative groupings or rejecting the hypothesis altogether (Albertz 2018), the existence of the P material is still widely acknowledged, owing to its distinctive style and content and relative internal consistency (e.g., Albertz 2018; Haran 1981; Holzinger 1893; Kuenen and Wicksteed 1886). The texts typically classified as P exhibit certain recurring characteristics, including highly structured formulations, genealogical lists, legal material and a focus on ritual and cultic concerns (e.g., Albertz 2018; Vink 1969). Scholars have further proposed that P itself is not monolithic but contains sub-strata, most notably distinguishing between the preceding (or succeeding) Priestly Source (P) and a later Holiness School (H) (e.g., Ross 1997), which is thought to emphasize purity laws and ethical concerns (e.g., Klostermann 1907; Knohl 2007). However, the extent to which these strata represent distinct compositional layers rather than later redactional activity is an ongoing debate, and no consensus exists regarding their precise historical development. Given these complexities, our study does not assume the validity of any particular compositional model but instead seeks to evaluate whether quantitative, entropy-based methods can identify distinct stylistic patterns in the Pentateuch better than previous computerized characterization attempts (Dershowitz et al. 2015; Faigenbaum-Golovin et al. 2025; Yoffe et al. 2023). If such patterns emerge, this may lend support to certain divisions proposed in source-critical scholarship. If they do not, it may suggest that the distinction between the various traditions, particularly between P and non-P, or P and H, is more fluid than often assumed. We selected P as the main test case for our clustering algorithm due to its highly structured and repetitive nature, which contrasts with the more vivid and colorful style of other Pentateuchal traditions. Because our approach leverages entropy-based clustering to detect formulaic patterns, it provides an opportunity to evaluate whether the hypothesized Priestly texts exhibit distinct formulaic stylistic tendencies that set them apart from surrounding material. Additionally, if a similar distinction between P and H exists (as proposed by Knohl (2007)), our method may provide an independent and complementary means of quantifying and assessing it. Notably, there is a broad literary consensus that P constitutes a stylistically distinct body of text, often considered a literary outlier due to its rigid structure, formulaic expressions and characteristic lexical choices (e.g., Bühler et al. 2024; Smith 1996). Here, we consider three hypothesized partitions of the three books in the Pentateuch for evaluation:

• • Genesis: We test the distinction between genealogical lists and the surrounding narrative material. The great majority of the lists are attributed to P due to their molded expression (e.g., the Toledot formula) and emphasis on lineage continuity (e.g., Long, Tengström, and Tengstrom 1984; Plum 1989). However, a few texts belong to non-P. Hence, the focus in this case is on matters of genre and form in the light of entropy rather than source division. By isolating these passages, we assess whether their self-information properties diverge significantly from the surrounding material. In the second stage, the possibility of internal division between P and non-P within the genealogical lists will be examined/demands further inquiry.

• • Exodus: We analyze the division between Priestly (P) and non-Priestly (non-P) material. Priestly passages in Exodus largely consist of legal and cultic instructions, particularly concerning the construction of the Tabernacle, whereas non-Priestly sections retain a more narrative-driven style (Römer 2009; Schmid 2015). This partition has been previously explored (Bühler et al. 2024; Yoffe et al. 2023), making it a useful benchmark for gauging the performance of our method against previous results.

• • Leviticus: We investigate whether a distinction can be quantitatively detected between the Priestly (P) and Holiness (H) materials. While both exhibit high degrees of formulaicity, the Holiness texts are often argued to introduce distinctive stylistic and theological elements (Knohl 2007). If our algorithm can recover a meaningful distinction, this may provide additional empirical support for the hypothesized P/H division.

Each of these test cases represents a distinct type of contrast in structural predictability. In Genesis, the distinction between genealogical lists and surrounding narrative material constitutes the most pronounced formulaicity contrast, as the genealogies are highly repetitive and rigid in structure. In Exodus, the contrast between Priestly and non-Priestly material also exhibits a clear difference, particularly in sections related to cultic legislation and Tabernacle construction, which display strongly regularized patterns (Bühler et al. 2024; Yoffe et al. 2023). The P/H distinction in Leviticus had not previously been analyzed using computerized methods. Its inclusion here allows us to assess whether subtle stylistic differences, such as those associated with the Holiness corpus, manifest in statistically measurable differences in self-information. Together, these examples are chosen to gauge the method’s ability to quantify varying degrees of formulaicity across multiple strata of the biblical text.

Experimental setup

Digital biblical corpus and annotation by experts

We use a digital corpus of the Masoretic variant of the Hebrew Bible in (biblical) Hebrew, which is a version of the Leningrad codex freely available by STEPBible. ² This dataset comes parsed with full morphological and semantic tags for all words, prefixes and suffixes. In this work, we consider only the morphological (i.e., grammatical) representation of the text, as word-based representations may introduce additional variability due to synonymy and context-dependent meanings. By focusing on morphological features, we ensure greater consistency in capturing structural patterns relevant to authorship and textual stratification.

For each hypothesized partition, we obtain expert annotations provided by biblical scholars specializing in source-critical analysis. These annotations serve as reference classifications, delineating the hypothesized textual divisions. The expert-labeled datasets allow us to evaluate the extent to which our clustering algorithm aligns with established scholarly hypotheses and assess its ability to recover traditionally proposed literary strata in an unsupervised manner.

While this study primarily employs morphological features, due to their interpretive clarity and availability in the STEPBible corpus, the method is not restricted to such representations. In settings where morphological annotation is unavailable, surface word forms can also be used. However, surface-level representations present certain challenges when evaluating formulaicity, as variation in inflection, cliticization or syntactic context may obscure structural repetition. Morphological normalization helps reduce this variability, revealing underlying patterns of grammatical constraint. Nonetheless, there are textual domains where surface-level lexical repetition may be more meaningful.

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 2014). 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. 2022), 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 2023), 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. 2015; Radday and Shore 1985). 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. 2023), 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.

Results

The results of our clustering analysis for the Genesis, Exodus and Leviticus hypothesized partitions are presented in Figures 3–5, respectively. These figures compare the performance of our information-based clustering method with k-means in identifying the hypothesized partitions of the text. To assess clustering accuracy, we categorize MCC scores into discrete performance intervals: 50%–74%, 75%–84%, 85%–89%, 90%–95%, and 96%–100%, which provide a benchmark for evaluating the reliability of different parameter settings. In the case of the genealogical lists in Genesis and the P/non-P partition in Exodus, our approach resulted in a substantially larger number of parameter combinations yielding MCC scores within the two highest MCC intervals (85%–89%, 90%–95%) compared to k-means.

While each individual parameter combination in our grid (defined by n, $\ell $ and f) may yield results of interpretive interest, the primary objective of this study is to introduce and validate the clustering framework rather than to evaluate any specific configuration in detail. Accordingly, we assess method performance by analyzing the distribution of MCC scores across the entire parameter space. If a textual partition genuinely reflects differences in structural predictability, then this signal should emerge robustly across a wide range of embedding configurations, not just under fine-tuned settings.

Conversely, when only a limited subset of parameter configurations yields strong agreement with a hypothesized division, this suggests that certain levels of segmentation and feature granularity are particularly well suited to capturing certain formulaic patterns present in the data. Rather than treating this as a sign of feature selection artifacts, we interpret such behavior as evidence that distinct formulaic patterns may emerge more clearly at specific scales of representation. Some clusters, especially those linked to genre, register or discourse framing, may only become statistically salient when, for example, the n-gram size and running window align with their internal structure. In this view, variation in clustering performance is not a sign of instability. Rather, it reflects both methodological sensitivity and the heterogeneous nature of formulaicity across textual layers.

This perspective is especially relevant for corpora where no predefined labels or gold standards exist. In such cases, exploring the consistency or divergence of outputs across parameter settings can help identify stable signals and locate boundaries where stylistic coherence shifts. Disagreements between parameter configurations are not necessarily failures of the method. Instead, they offer interpretive entry points into the structure of the data and the scale at which formulaic patterns operate.

Figure 6. Distinctive n-grams extracted from the formulaic cluster for two parameter combinations, capturing 30% of the variance (see Section 2.8 in Yoffe et al. (2023)). Left panel: Clustering of morphologically-represented Leviticus using $\ell = 28$ , $n = 4$ and $f = \text {all}$ , achieving an MCC score of 94%. Right panel: Similar to the left panel but with $\ell = 20$ , $n = 2$ and $f = 300$ , achieving an MCC score of 93%. n-grams discussed in the “Formulaic Structure and Parameter Sensitivity in the P/H Partition of Leviticus” section as examples for H- and P-associated features are outlined in red. The insets display the self-information distributions of both clusters, with blue and orange representing the non-formulaic and formulaic clusters, respectively.

In the case of the P/H partition in Leviticus, our method yielded more parameter combinations in the highest MCC interval by roughly a factor of three, but k-means yielded slightly more MCC scores in the highest two intervals overall. This can be attributed to the fact that not all parameter combinations produce partitions with strong auto-entropy contrasts. Unlike our information-based approach, which depends on isolating structural predictability differences, k-means clusters based on direct feature-space separations, allowing it to identify partitions even when distinctions arise from lexical or thematic variation rather than underlying structural consistency. Because k-means optimizes intra-cluster similarity without requiring an information-theoretic distinction, it remains more robust across a wider range of parameter settings, whereas information-based clustering excels only when the chosen parameters successfully capture fundamental differences in textual predictability. In the “Formulaic Structure and Parameter Sensitivity in the P/H Partition of Leviticus” section, we conduct an in-depth analysis of the results for the P/H partition in Leviticus and discuss their implications.

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 2007). 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 2007). 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.

Discussion

In this study, we introduced an information-theoretic soft clustering framework designed to identify structured patterns within textual data. Our approach leverages self-information as a statistical indicator of structural regularities, enabling a systematic exploration of linguistic and stylistic consistency across different types of corpora. By treating clustering as a probabilistic process rather than enforcing hard boundaries, our method accommodates uncertainty in textual attributions, making it particularly suited for cases where transitions between sources or styles are gradual rather than abrupt.

We validated our method on both Gaussian-simulated and categorical datasets to assess its effectiveness across different data structures. Benchmarking against traditional clustering techniques, including k-means and various GMMs, demonstrated that our approach provides a more stable classification of structured and unstructured components, especially in cases where the sample-to-dimension ratio is small. The information-based constraint further improved cluster separability, ensuring that formulaic structures were identified with higher precision.

The Gaussian simulations allowed us to evaluate performance under controlled conditions, where we modeled clusters using multivariate Gaussian distributions with varying covariance structures. Our method proved effective in correctly classifying clusters despite high-dimensional covariance estimation challenges, particularly in settings where the number of observations was small relative to the number of features. Unlike conventional clustering approaches that rely on full covariance estimation, which can be unstable in low-sample regimes, our information-based approach remains robust by leveraging structural differences rather than density estimation alone.

The categorical data experiments, in turn, tested the robustness of our framework in realistic scenarios where feature dependencies, sparsity and small sample sizes of high-dimensional datasets introduce significant challenges to clustering. These experiments demonstrated that our method remains effective even when conventional clustering approaches struggle due to the high-dimensional and discrete nature of the data. Our model successfully captured co-occurrence structures within the data, identifying clusters with statistically significant variations in information-based distributions, making it particularly suitable for structured linguistic data where underlying dependencies are not easily captured by Euclidean distance-based methods.

Applying this method to the biblical corpus, we analyzed key textual partitions to determine whether the identified clusters align with traditional source-critical hypotheses. Our approach consistently yielded the highest classification accuracy across various parameter configurations, surpassing k-means in the number of high-confidence classifications. In Genesis and Exodus, our method produced a significantly larger number of classifications within the two highest MCC intervals (85%–89% and 90%–95%) compared to k-means, indicating its robustness in capturing the structural properties of these partitions. In Leviticus, while k-means exhibited more high-MCC classifications overall, our approach yielded approximately twice the number of classifications in the highest MCC interval (90%–95%), demonstrating its effectiveness in isolating the most formulaic structures within the text.

These findings suggest that information-theoretic clustering offers a promising avenue for unsupervised textual analysis, providing an independent measure of structural consistency that does not rely on predefined linguistic features. The stability of our method across different corpora highlights its potential applicability beyond biblical studies, particularly in fields where latent structures within textual datasets require tractable and unbiased identification.

More broadly, our findings demonstrate that structural predictability, quantified through sample-wise self-information, captures a core dimension of textual organization that aligns with, and in some cases clarifies, traditional source-critical hypotheses. In the case of computerized source criticism of the Hebrew Bible, prior computational studies have shown that blocks, such as P and non-P differ in feature distributions (Dershowitz et al. 2015; Faigenbaum-Golovin et al. 2025; Yoffe et al. 2023). Here, we show that they also differ in their degree of formulaicity. This distinction is especially significant when made explicit through statistical modeling, as it formalizes what has often been an implicit intuition in biblical exegesis: that certain textual strata are not only lexically distinct but also stylistically more rigid, repetitive or constrained. Crucially, our results further reveal that the expression of formulaicity is dependent on the feature set used; certain structures yield different interpretations of which passages are most internally consistent. This variability underscores the interpretive value of entropy-based methods: rather than producing a single fixed partition, they allow scholars to examine how different linguistic scales foreground different dimensions of textual structure, enabling a more nuanced and data-grounded exploration of compositional layers.

This makes entropy-based clustering particularly well suited for historical texts with compositional depth, where multiple stylistic logics may operate simultaneously. In such settings, soft clustering based on self-information variation provides a principled way to detect internal regularities without relying on predefined labels or explicit boundaries. Exploring the parameter space, across different n-gram granularities, segmentation schemes and feature types, becomes an interpretive tool for probing how distinct textual regimes express constraint, formulaicity or stylistic coherence.

For example, the distinction between Early and Late Biblical Hebrew is often expressed in systematic lexical and idiomatic shifts, making diachronic strata promising candidates for unsupervised clustering based on word-level features (Hurvitz 1968, 2014). In the Deuteronomistic history, recurring thematic formulas and narrative structures may yield internally coherent segments whose predictability aligns with redactional layers (Peckham 2019). Rabbinic legal corpora, such as the Mishnah, show compositional depth through repeated conditional constructions and formulaic attributions that vary across tractates or schools (Fraade 1991; Neusner 1988). In such cases, entropy-based clustering may help identify stylistic layers that are not organized by topic or subject matter, but by the underlying patterns of discourse structure and phrasing. Similarly, in oral-derived corpora like the Old Norse sagas, genealogical prologues, episodic transitions and set-piece scenes create localized structural regularities that could be identified by entropy-guided segmentation (Byock 1984; Clover 1982). In all these contexts, the method can be used to test source-critical hypotheses, uncover latent boundaries or evaluate competing assumptions about where and how regularity is embedded in the textual fabric.

In future work, we will extend our framework by incorporating additional embedding techniques, such as tf-idf representations, neural embeddings and transformer-based contextualized representations. These approaches will allow us to assess the model’s performance across a broader range of textual features and evaluate its adaptability to modern NLP methodologies. By integrating embeddings that capture both statistical word importance and deep semantic relationships, we aim to refine the clustering process and further improve the detection of latent structural patterns in diverse textual corpora.