Unveiling Molecular Secrets: An LLM-Augmented Linear Model for Explainable and Calibratable Molecular Property Prediction

## Abstract

Abstract Explainable molecular property prediction is essential for various scientific fields, such as drug discovery and material science.
Despite delivering intrinsic explainability, linear models struggle with capturing complex, non-linear patterns. Large language models (LLMs), on the other hand, yield accurate predictions through powerful inference capabilities yet fail to provide chemically meaningful explanations for their predictions. This work proposes a novel framework, called MoleX , which leverages LLM knowledge to build a simple yet powerful linear model for accurate molecular property prediction with faithful explanations. The core of MoleX is to model complicated molecular structure-property relationships using a simple linear model, augmented by LLM knowledge and a crafted calibration strategy. Specifically, to extract the maximum amount of task-relevant knowledge from LLM embeddings, we employ information bottleneck-inspired fine-tuning and sparsity-inducing dimensionality reduction. These informative embeddings are then used to fit a linear model for explainable inference. Moreover, we introduce residual calibration to address prediction errors stemming from linear models’ insufficient expressiveness of complex LLM embeddings, thus recovering the LLM’s predictive power and boosting overall accuracy. Theoretically, we provide a mathematical foundation to justify MoleX ’s explainability. Extensive experiments demonstrate that MoleX outperforms existing methods in molecular property prediction, establishing a new milestone in predictive performance, explainability, and efficiency. In particular, MoleX enables CPU inference and accelerates large-scale dataset processing, achieving comparable performance 300 × \times × faster with 100,000 fewer parameters than LLMs. Additionally, the calibration improves model performance by up to 12.7% without compromising explainability. The source code is available at https://github.com/MoleX2024/MoleX .

## 1 Introduction

Molecular property prediction, aiming to analyze the relationship between molecular structures and properties, is crucial in various scientific domains, such as computational chemistry and biology (Xia et al., ; Yang et al., ). Deep learning advancements have significantly improved this field, showcasing the success of AI-driven problem-solving in science. Representative deep models for predicting molecular properties include graph neural networks (GNNs) (Lin et al., ; Wu et al., ) and LLMs (Chithrananda et al., ; Ahmad et al., ). In particular, recently developed LLMs have exhibited remarkable performance by learning chemical semantics from text-based molecular representations, e.g., Simplified Molecular Input Line Entry Systems (SMILES) (Weininger, ). By capturing the chemical semantics and long-range dependencies in text-based molecules, LLMs show promising capabilities in providing accurate molecular property predictions (Ahmad et al., ). Nevertheless, the black-box nature of LLMs hinders the understanding of their decision-making mechanisms. Inevitably, this opacity prevents people from deriving reliable predictions and insights from these models (Wu et al., ).

Figure: Figure 1: The framework of MoleX is divided into the following stages: ($1$) given ChemBERTa-2 as the pre-trained LLM, ($2$) fine-tune it on Group SELFIES (functional group-based molecular representation) with an information bottleneck-inspired objective to produce embeddings with maximum task-relevant information, ($3$) extract high-dimensional LLM embeddings and apply sparsity-inducing dimensionality reduction to exclude redundant information, ($4$) train a linear model using the preserved task-relevant information, ($5$) integrate the linear model with a residual calibrator that recaptures prediction errors for explainable inference.
Refer to caption: x1.png

To narrow this gap, numerous explainable GNN and LLM methods have been proposed to identify molecular substructures that contribute to specific properties (Xiang et al., ; Proietti et al., ; Wang et al., ). Among these, Lamole (Wang et al., ) represents the state-of-the-art LLM-based approach attempting to provide both accurate predictions and chemically meaningful explanations—chemical concepts-aligned substructures along with their interactions. However, it still suffers from several flaws: first, the attention weights used for explanations do not correlate directly with feature importance (Jain & Wallace, ); second, it is model-specific due to varying implementations and interpretations of attention mechanisms across models (Voita et al., ); and third, the provided explanations are local, struggling to approximate global model decisions using established chemical concepts (Liu et al., ). Therefore, it is imperative to design a globally explainable method that provides accurate predictions as well as contributing substructures, along with their interactions, to explain molecular property predictions.

This work proposes a new framework (illustrated in Figure 1), dubbed MoleX, that leverages a linear model augmented with LLM knowledge for explaining complex, non-linear molecular structure-property relationships, motivated by its simplicity and global explainability. To capture these complex relationships, MoleX extracts informative knowledge/embeddings from the LLM, which serve as inputs to fit a linear model. Moreover, we design information bottleneck-inspired fine-tuning and sparsity-inducing dimensionality reduction to maximize task-relevant information in LLM embeddings. Following prior work (Wang et al., ), we use Group SELFIES (Cheng et al., )—a text-based molecular representation that partitions molecules into functional groups—as the LLM’s input (as shown in section A.14). Group SELFIES enables LLMs to tokenize molecules into units of functional groups, aligning with chemical concepts at the substructure level. To properly quantify individual functional groups’ contributions to properties, we extract n-gram functional groups from Group SELFIES. These n-grams are fed into the LLM, yielding embeddings where functional groups are semantically distinct and enables a more nuanced analysis of their isolated impacts. Notably, MoleX’s simplicity allows for the approximation of model behaviors over the entire input space during inference, offering global explanations rather than interpreting particular samples.

Despite being augmented with LLM knowledge, linear models still underfit complex, non-linear relationships. Thus, we propose a residual calibration strategy to address prediction errors arising from the gap between high-dimensional LLM embeddings and linear models’ limited expressiveness. Specifically, the residual calibrator is designed to learn the linear model’s residuals. Its outputs are integrated into the linear model, encouraging iterative adjustment of prediction errors. Through sequentially driving residuals toward target values, the residual calibrator recaptures samples missed by the linear model and recovers the original LLM’s predictive power. Consequently, the linear model, incorporated with LLM knowledge and a residual calibrator, provides outstanding predictive performance while enjoying the linear model’s explainability benefit. Our contributions are summarized as

- 1.
We propose MoleX, which extracts LLM knowledge to build a simple yet powerful linear model that identifies chemically meaningful substructures with their interactions for explainable molecular property predictions.
- 2.
We develop optimization-based methods to maximize and preserve task-relevant information in LLM embeddings and theoretically demonstrate their explainability and validity.
- 3.
We design a residual calibration strategy that refits samples missed by the linear model, recovering the LLM’s predictive performance while preserving explainability.
- 4.
We introduce n-gram coefficients, with a theoretical justification, to assess individual functional group contributions to molecular property predictions.

Experiments across 7 datasets demonstrate that MoleX achieves state-of-the-art classification and explanation accuracy while being 300$\times$ faster with 100,000 fewer parameters than alternative baselines, highlighting its superiority in predictive performance, explainability, and efficiency.

## 2 Related Work

Explainable Molecular Property Prediction. Given that molecules can be naturally represented as graphs, a collection of explainable GNNs have been proposed to explain the relationship between molecular structures and properties (Lin et al., ; Pope et al., ). However, these atom or bond-level explanations are not chemically meaningful to interpret their sophisticated relationships. Besides, through learning chemical semantics, the transformer-based LLMs can effectively capture interactions among substructures (Wang et al., ) and thus demonstrated their potential in understanding text-based molecules (Ross et al., ; Chithrananda et al., ). However, the opaque decision-making process of LLMs obscures their operating principles, risking unfaithful predictions with severe consequences, especially in high-stakes domains like drug discovery (Chen et al., ).

Explainability Methods for LLMs. To obtain trustworthy output, various techniques were introduced to unveil the LLM’s explainability. The gradient-based explanations analyze the feature importance by computing output partial derivatives with respect to input (Sundararajan et al., ). These methods, nevertheless, lack robustness in their explanations due to sensitivity to data perturbations (Kindermans et al., ; Adebayo et al., ). The attention-based explanations use attention weights to interpret outputs (Hoover et al., ). Yet, recent studies challenge their reliability as attention weights may not consistently reflect true feature importance (Jain & Wallace, ; Serrano & Smith, ). The perturbation-based explanations elucidate model behaviors by observing output changes in response to input alterations (Ribeiro et al., ). However, these explanations are unstable due to the randomness of the perturbations (Agarwal et al., ). To resolve these issues, we extract informative embeddings from the LLM to fit a linear model for inference. This approach leverages both the LLM’s knowledge and the linear model’s explainability, offering reliable substructure-level explanations.

## 3 Preliminaries

Let $g^{(i)}\in\mathcal{G}=\{(g^{(i)},y^{(i)})\}$ be a molecular graph, our goal is to train a model $f$ to map the molecular representation $g$ to its property $y$, denoted as $f:g\rightarrow y$. We first convert the molecular graph $g^{(i)}$ into a text-based molecular representation (i.e., Group SELFIES), denoted as $x^{(i)}=\{x_{1}^{(i)},\ldots,x_{j}^{(i)}\}$, where $x_{j}^{(i)}$ is the $j$-th functional group and $\mathcal{S}_{D}$ is the dataset we used. The LLM-augmented linear model consists of two modules: an explainable model and a residual calibrator. After the explainable model $h$ predicts, its residuals are fed into the residual calibrator $r$, which boosts the performance without incurring any explainability impairment. We denote $f_{H}$ and $f_{R}$ as features of explainable model and residual calibrator, $\mathcal{L}(\hat{y},y)$ as the training loss. To learn $h$ and $r$, we freeze parameters of the former and sequentially refit its residuals with the objective:

$$ $\min_{h,r}\,\mathbb{E}_{(x,y)\sim\mathcal{S}_{D}}\left[\mathcal{L}\left(h\left (f_{H}(x)\right)+r\left(f_{R}(x)\right),y\right)\right].$ (3.1) $$

Adapting the approach by Sebastiani (), we use n-gram coefficients in the linear model to measure the contributions of decoupled n-gram features/functional groups to molecular properties. Let the n-gram feature $x_{j}$ takes the coefficient $w_{j}$ in the linear model, then its contribution score is computed as $c_{j}=w_{j}\cdot\text{Embedding($x_{j}$)}$. This allows us to quantify the contribution of $j$-th functional group to property $y$. Our proof of the validity of using n-gram coefficients as contribution scores is provided in section A.1. In the following parts, we omit the superscript ^(i) for simplicity.

## 4 Our Framework: MoleX

MoleX does two things, i.e., ($1$) maximizing and preserving the task-relevant information in LLM embeddings via fine-tuning and dimensionality reduction and ($2$) extracting these embeddings to build an LLM-augmented linear model with residual calibration. We thus divide it into two stages: LLM knowledge extraction and LLM-augmented linear model fitting. This section details our framework and provides theoretical foundations for its explainability.

### 4.1 LLM Knowledge Extraction with Improved Informativeness

Fine-tuning. We fine-tune a pre-trained chemical LLM using Group SELFIES data to enhance its understanding of functional group-based molecules. To maximize the linear model’s effectiveness, we aim to extract embeddings that are as informative as possible, fully exploiting the LLM’s interior knowledge. However, empirically fine-tuning the LLM to produce embeddings with the desired informativeness poses a significant challenge. To address this, we apply the Variational Information Bottleneck (VIB) (Alemi et al., ) into the fine-tuning process by crafting a training loss to improve embedding informativeness. By optimizing the VIB-based loss function, we encourage the LLM to produce embeddings with the maximum amount of downstream task-relevant information, thus significantly enhancing their informativeness. Particularly, given Group SELFIES inputs $x$, properties $y$, and LLM embeddings $t$, we define $p_{0}(t)$ as the prior distribution over $t$ (e.g., standard normal distribution), and $q_{\theta}(y|t)$ as the variational approximation to the conditional distribution of the properties given the embeddings $t$. The mutual information between $t$ and $y$ is defined as:

$$ $I(t;y)=\mathbb{E}_{p(t,y)}\left[\log\frac{p(t,y)}{p(t)p(y)}\right]=\mathbb{E}_ {p(t,y)}\left[\log\frac{p(y|t)}{p(y)}\right],$ $$

and the mutual information between $t$ and $x$ is defined as:

$$ $I(t;x)=\mathbb{E}_{p(t,x)}\left[\log\frac{p(t|x)}{p(t)}\right]=\mathbb{E}_{p(x )}\left[D_{\mathrm{KL}}\left(p_{\theta}(t|x)\,\big{\|}\,p(t)\right)\right].$ $$

Since the marginal distribution $p(t)$ is intractable, we approximate it with the prior $p_{0}(t)$. Under this approximation, we use $D_{\mathrm{KL}}\left(p_{\theta}(t|x)\,\big{\|}\,p_{0}(t)\right)$ as a tractable surrogate for $I(t;x)$. Inspired by Kingma et al. (), we approximate the encoder $p_{\theta}(t|x)$ by a Gaussian distribution. Let $f_{e}^{\mu}(x)$ and $f_{e}^{\Sigma}(x)$ be neural networks that output the mean and covariance matrix of the latent variable $t$. Then, the encoder is given as:

$$ $p_{\theta}(t|x)=\mathcal{N}\left(t\,\big{|}\,f_{e}^{\mu}(x),f_{e}^{\Sigma}(x) \right).$ $$

Applying the reparameterization trick, we sample $t$ as:

$$ $t=f_{e}^{\mu}(x)+f_{e}^{\Sigma}(x)^{1/2}\cdot\epsilon,\quad\text{where } \epsilon\sim\mathcal{N}(0,I).$ $$

Putting all these together, we design our training loss as:

$$ $\mathcal{L}(\theta)=\sum_{(x_{i},y_{i})\in\mathcal{S}_{F}}\left(\mathbb{E}_{p_ {\theta}(t|x_{i})}\left[-\log q_{\theta}(y_{i}|t)\right]+\beta\cdot D_{\mathrm {KL}}\left(p_{\theta}(t|x_{i})\,\big{\|}\,p_{0}(t)\right)\right),$ (4.1) $$

where $\beta$ is the tuning parameter between compression and performance, $q_{\theta}$ is the decoder (predictive model), and $\mathcal{S}_{F}$ is the dataset used for fine-tuning. In particular, the first component, $\mathbb{E}_{p_{\theta}(t|x_{i})}\left[-\log q_{\theta}(y_{i}|t)\right]$, maximizes the task-relevant information about $y$ in the produced embeddings $t$. The second component, $\beta\cdot D_{\mathrm{KL}}\left(p_{\theta}(t|x_{i})\parallel p_{0}(t)\right)$, minimizes the redundant information from $x$ in the produced embeddings $t$. Empirically, we use ChemBERTa-2 (Ahmad et al., ) as the foundation LLM for fine-tuning.

In essence, this training loss ensures that the fine-tuned model produces embeddings $t$ that capture relevant information about the molecular property $y$ while compressing redundant information from $x$. Our design, grounded in the information bottleneck principle, demonstrably leads to informative embeddings with strong explainability. We claim this loss will converge to informative embeddings and provide the proof in section A.2.

###### Theorem 4.1 (Convergence to an Informative Representation in LLM Fine-tuning) .

Let $\mathcal{L}(\theta)$ be the loss defined in eq. 4.1. Under the assumptions of the reparameterization trick and the use of stochastic gradient descent, the optimization process converges to a local minimum that yields an informative representation $t$ while retaining only the most relevant information from the task.

Embedding Extraction. To capture individual functional group contributions and contextual information, we extract n-grams from Group SELFIES, with $n$ selected via cross-validation. To ensure explainability, each n-gram is processed separately by the fine-tuned LLM using a functional group-level tokenizer, encoding a fixed-size embedding vector. These vectors are then aggregated into a single fixed-size embedding that encompasses semantics of all individual n-grams. More precisely, this single embedding contains all chemical semantics at the functional group level and represents the knowledge LLM learned during its training and fine-tuning.

###### Theorem 4.1 (Convergence to an Informative Representation in LLM Fine-tuning) .

### 4.2 Dimensionality-reduced embeddings for linear model fitting

Dimensionality Reduction. As the aggregated n-gram embeddings are still high-dimensional and noisy, eliminating the redundancy in these informative embeddings becomes our new problem. Drawing inspiration from Lin et al. (), we design an explainable functional principal component analysis (EFPCA) that leads to effective dimensionality reduction. Accordingly, we can preserve a compact yet statistically significant set of features for the linear model. In this part, we demonstrate how to convert this dimensionality reduction into an optimization problem by introducing a sparsity-inducing penalty into the optimization objective function. We define our method as

###### Definition 4.1 (EFPCA) .

Let $X(t)$ be a stochastic process defined on a compact interval $[a,b]$ with mean function $\mu(t)=\mathbb{E}[X(t)]$. Assume that $X(t)$ has a covariance operator $\hat{\mathcal{C}}$ derived from the centered data $X(t)-\mu(t)$. The EFPCA seeks to find functions $\xi_{k}(t)$ that maximize the variance explained by the projections of $X(t)$ while promoting sparsity for explainability. Then, the EFPCA solves the following optimization objective:

$$ $\max_{\xi_{k}}\left\{\langle\xi_{k},\hat{\mathcal{C}}\xi_{k}\rangle-\rho_{k} \cdot\mathcal{S}(\xi_{k})\right\}$ $$

subject to $\|\xi_{k}\|_{\gamma}^{2}=\|\xi_{k}\|^{2}+\gamma\|\mathcal{D}^{2}\xi_{k}\|^{2}=1$ and $\langle\xi_{k},\xi_{j}\rangle_{\gamma}=0\quad\text{for all }j<k$,

where $\gamma>0$ is a tuning parameter that balances the fit and smoothness, $\mathcal{S}(\xi_{k})$ measures the length of the support of $\xi_{k}(t)$, promoting sparsity and explainability, and $\rho_{k}>0$ is a tuning parameter controlling the sparsity of $\xi_{k}(t)$.

The $\ell_{0}$ penalty term $\rho_{k}\|a_{k}\|_{0}$ promotes sparsity by encouraging many coefficients $a_{kj}$ to be exactly zero when $\rho_{k}$ is sufficiently large. This forces principal components $\xi_{k}(t)$ to be exactly zero over extensive portions of the domain $[a,b]$. Since $\xi_{k}(t)$ is a linear combination of basis functions, zero coefficients directly cause zero contributions from those basis functions across their supports. Specifically, the optimization balances maximizing the variance captured by $\xi_{k}(t)$ while minimizing the number of nonzero coefficients, thereby preserving only the most significant components. Furthermore, since the basis functions $\phi_{j}(t)$ have local support on subintervals $S_{j}\subset[a,b]$, the nonzero coefficients correspond to basis functions that are active only over those specific intervals. Therefore, the sparsity in the coefficients $a_{k}$ and the local support of the basis functions $\phi_{j}(t)$ result in principal components $\xi_{k}(t)$ being nonzero only over certain intervals. Mathematically, the support of $\xi_{k}(t)$ is the union of the supports $S_{j}$ corresponding to the nonzero coefficients $a_{kj}$. The EFPCA thus produces principal components that are sparse and explainable due to their localized structure, as they highlight the regions where the data exhibits significant variation.

In summary, EFPCA provides a framework for obtaining highly explainable principal components, enabling effective dimensionality reduction. We demonstrate that by combining the sparsity-inducing penalty with the local support of the basis functions, the resulting principal components are both sparse and localized, successfully capturing significant features of the data. In our implementation, we exclude irrelevant functional groups and detect principal ones from the high-dimensional embeddings. Based on this, we claim the following theorem (see our proof in section A.3):

###### Theorem 4.2 .

The EFPCA produces sparse functional principal components $\xi_{k}(t)$ that are exactly zero in intervals where the sample curves exhibit minimal variation. Consequently, the FPCs $\xi_{k}(t)$ are statistically significant and explanatory, facilitating effective dimensionality reduction.

Linear Model Fitting. Applying dimensionality-reduced n-gram embeddings as features, we train a logistic regression model for our classification tasks, which takes the form:

$$ $h(f_{H}(\mathbf{x}))=\sigma\left(\mathbf{w}^{\top}f_{H}(\mathbf{x})+b\right)= \frac{1}{1+e^{-(\mathbf{w}^{\top}f_{H}(\mathbf{x})+b)}},$ (4.2) $$

where $\sigma$ is the sigmoid function, $\mathbf{w}\in\mathbb{R}^{n}$ is the weight vector, $b\in\mathbb{R}$ is the bias term, and $f_{H}(\mathbf{x})$ are explainable features as defined in eq. A.11. In our settings, the logistic regression is explainable since the log-odds transformation establishes a linear relationship between features and the target variable such that
$\log\left(\frac{h(f_{H}(\mathbf{x}))}{1-h(f_{H}(\mathbf{x}))}\right)=\mathbf{w
}^{\top}f_{H}(\mathbf{x})+b$.
Differentiating with respect to a feature component $[f_{H}(\mathbf{x})]_{j}$ shows that each coefficient $w_{j}$ quantifies the impact of that feature on the log-odds such that $\frac{\partial}{\partial[f_{H}(\mathbf{x})]_{j}}\log\left(\frac{h(f_{H}(
\mathbf{x}))}{1-h(f_{H}(\mathbf{x}))}\right)=w_{j}$.
Moreover, if $f_{H}$ is a linear transformation, i.e., $f_{H}(\mathbf{x})=\mathbf{C}\mathbf{x}$, the chain rule relates changes in the original features to the log-odds which can be expressed as $\frac{\partial}{\partial x_{j}}\log\left(\frac{h(f_{H}(\mathbf{x}))}{1-h(f_{H}
(\mathbf{x}))}\right)=\sum_{k=1}^{n}w_{k}C_{kj}.$
Therefore, this linearity allows straightforward interpretation of each feature’s influence on the predicted probabilities, making logistic regression highly explainable (Hastie et al., ). By directly linking coefficients to features, it offers an intuitive understanding of each feature’s impact on the output.

Residual Calibration. The final step of MoleX involves learning a residual calibrator $r$ to improve predictive accuracy while preserving explainability. We freeze parameters of the explainable model $h$ and calibrate samples that $h$ fails to predict. By optimizing the objective function in eq. A.11, we iteratively fix prediction errors, sequentially driving the overall predictions closer to the target values. The residual calibrator is designed as a linear model, ensuring that the overall model remains explainable. Specifically, we define the residual calibrator $r$ with weights $w_{r}\in\mathbb{R}^{d_{r}}$ corresponding to each residual feature and bias $b_{r}$:

$$ $r(f_{R}(x))=w_{r}^{\top}f_{R}(x)+b_{r}.$ $$

Here, $f_{R}(x)$ represents the residual features obtained from the decomposition of the feature space $\mathbb{R}^{d}$ into orthogonal subspaces such that $f(x)=f_{H}(x)+f_{R}(x),$ where $f_{H}(x)\in\mathbb{R}^{d_{c}}$ contains the explainable features used by $h$, and $f_{R}(x)\in\mathbb{R}^{d_{r}}$ contains the residual features used by $r$, with the orthogonality condition given by $\langle f_{H}(x),f_{R}(x)\rangle=0.$ Then, the overall prediction combines the contributions from $h$ and $r$:

$$ $\hat{y}(x)=\underbrace{w_{h}^{\top}f_{H}(x)+b_{h}}_{\text{Explainable Model Contribution}}+\underbrace{w_{r}^{\top}f_{R}(x)+b_{r}}_{\text{Residual Calibrator Contribution}}.$ $$

The orthogonality and linearity between $f_{H}(x)$ and $f_{R}(x)$ guarantee that the contributions from $h$ and $r$ are additive and independent, making the residual calibrator $r$ theoretically explainable. Moreover, each feature’s impact on the prediction can be directly understood through the corresponding weights in $w_{h}$ and $w_{r}$. Since $f_{H}(x)$ and $f_{R}(x)$ are orthogonal, the inner products $w_{h}^{\top}f_{R}(x)=0$ and $w_{r}^{\top}f_{H}(x)=0$ vanish. This ensures that $h$ and $r$ do not influence each other’s feature contributions, thus preserving the explainability of both models in the combined prediction. We thus formalize the following theorem (see our proof in section A.4):

###### Theorem 4.3 .

Let $\mathcal{X}$ and $\mathcal{Y}$ be the input and output spaces, respectively. Let $f:\mathcal{X}\rightarrow\mathbb{R}^{d}$ be a pre-trained feature mapping, and let $h:\mathbb{R}^{d_{c}}\rightarrow\mathcal{Y}$ be an explainable linear model operating on the explainable features $f_{H}(x)$. The residual calibrator $r:\mathbb{R}^{d_{r}}\rightarrow\mathcal{Y}$, defined on the residual features $f_{R}(x)$, captures the variance not explained by $h$ in an explainable manner, thereby preserving the overall model’s explainability.

Quantifiable Functional Group Contributions. As described in section 3, we measure the functional group $x_{j}$’s contributions to molecular property $y$ using n-gram coefficients. The molecular property $y$ distributes its entire semantic information into individual functional groups $x_{j}$. Due to the linearity and additivity between $x_{j}$ and $y$, the scalar coefficient $w_{j}$ corresponding to $x_{j}$ in the linear model weighs $x_{j}$’s contributions to $y$ in terms of chemical semantics. By taking the dot product of $w_{j}$ and the embedding of $x_{j}$, we obtain a projection length of the functional group in the direction of weight vector, thus quantifying the impact of that functional group on the molecular property. Quantitatively, the larger the absolute value of an n-gram coefficient, the greater the contribution of the corresponding functional group to property. This metric provides a rigorous interpretation of feature contributions, ensuring unbiasedness and significance through OLS estimation (see our proof in section A.1). Using this method, we identify important functional groups from the LLM’s complex embedding space. Furthermore, by incorporating n-gram coefficients and identified functional groups into the molecular graph, we can determine whether identified functional groups bond with each other and infer interactions among them. Essentially, this metric, followed by a mathematical justification, demonstrably provides novel insights into quantifiable functional group contributions in molecular property prediction. Based on this, MoleX reveals chemically meaningful substructures along with their interactions to faithfully explain molecular property predictions.

###### Definition 4.1 (EFPCA) .

###### Theorem 4.2 .

###### Theorem 4.3 .

## 5 Experiments

### 5.1 Experimental Settings

Datasets. We empirically evaluate MoleX’s performance on six mutagenicity datasets and one hepatotoxicity dataset. The mutagenicity datasets include Mutag (Debnath et al., ), Mutagen (Morris et al., ), PTC family (i.e., PTC-FM, PTC-FR, PTC-MM, and PTC-MR) (Toivonen et al., ) and the hepatotoxicity dataset includes Liver (Liu et al., ). To demonstrate that MoleX can explain molecular properties using chemically meaningful substructures, we introduce the concept of ground truth: substructures verified by domain experts to have significant impacts on molecular properties. The ground truth substructures for six mutagenicity datasets are provided by Lin et al. (); Debnath et al. (), while those for the hepatotoxicity dataset are provided by Cheng et al. (). Further details are available in section A.5.

Evaluation Metrics. We show that MoleX achieves notable predictive performance, explainability performance, and computational efficiency. We apply a specific metric to evaluate each aspect of the model performance. For predictive performance, we define $\frac{1}{\mathrm{I}}\sum_{i=1}^{I}\mathbb{I}(y^{(i)}=\hat{y}^{(i)})$ to compute the classification accuracy. For explainability performance, we follow the settings in GNNExplainer (Ying et al., ) where explanations are treated as a binary classification of edges (the prediction that matches ground truth substructures as positive, and vice versa) and use the Area Under the Curve (AUC) to quantitatively measure the explanation accuracy. For computational efficiency, we measure the execution time in microseconds for each method.

Baselines. To extensively compare *MoleX* with different methods, we utilize (1) GNN baselines, including GCN (Kipf & Welling, ), DGCNN (Zhang et al., ), edGNN (Jaume et al., ), GIN (Xu et al., ), RW-GNN (Nikolentzos & Vazirgiannis, ), DropGNN (Papp et al., ), and IEGN (Maron et al., ); (2) LLM baselines, including Llama 3.1-8b (Dubey et al., ), GPT-4o (Achiam et al., ), and ChemBERTa-2 (Ahmad et al., ); (3) explainable model baselines, including logistic regression, decision tree (Quinlan, ), XGBoost (Chen & Guestrin, ), and random forest (Breiman, ).

Implementations. Our model is pre-trained on the full ZINC dataset (Irwin et al., ) using ChemBERTa-2, with $15\%$ of tokens in each input randomly masked. We then fine-tune this model on the Mutag, Mutagen, PTC-FM, PTC-FR, PTC-MM, PTC-MR, and Liver datasets (in Group SELFIES). To evaluate model performance, we compute the average and standard deviation of each metric for each method after 20 rounds of execution. Further details are provided in section A.6.

### 5.2 Results

Predictive Performance. Table 1 presents a comparison of predictive performance across different methods. MoleX consistently outperforms all baselines, indicating its robustness and generalizability. Specifically, as a combination of LLMs and explainable models, MoleX achieves better performance compared to either of them (i.e., $16.9\%$ and $23.1\%$ higher average classification accuracy than LLM and explainable model baselines, respectively) and demonstrates the effectiveness of augmenting explainable models with LLM knowledge. Moreover, by integrating residual calibration, MoleX raises the average classification accuracy by $7.0\%$ across seven datasets. Notably, the classification accuracy of our base model, logistic regression, improves by $27.8\%$ after LLM knowledge augmentation and then by an additional $5.5\%$ after residual calibration on the Mutag dataset. Therefore, by maximizing task-relevant semantic information in the LLM knowledge and employing a residual calibration strategy, we enable a simple linear model to achieve predictive performance even superior to that of GNNs and LLMs in molecular property predictions.

**Table 1: Classification accuracy over seven datasets (%). The best results are highlighted in bold.**
| Methods | Mutag | Mutagen | PTC-FM | PTC-FR | PTC-MM | PTC-MR | Liver |
| --- | --- | --- | --- | --- | --- | --- | --- |
| GCN (Kipf & Welling, ) | 83.4± 0.4 | 77.2± 0.7 | 56.5± 0.3 | 62.7± 0.5 | 58.3± 0.2 | 52.1± 0.6 | 40.6± 0.3 |
| DGCNN (Zhang et al., ) | 86.2± 0.2 | 73.7± 0.5 | 56.1± 0.4 | 64.0± 0.8 | 61.8± 0.7 | 57.1± 0.6 | 45.4± 0.9 |
| edGNN (Jaume et al., ) | 85.4± 0.6 | 76.5± 0.3 | 58.7± 0.4 | 66.3± 0.7 | 65.2± 0.6 | 55.1± 0.8 | 43.7± 0.4 |
| GIN (Xu et al., ) | 86.1± 0.3 | 81.0± 0.5 | 63.4± 0.8 | 67.8± 0.6 | 66.5± 0.4 | 65.5± 0.4 | 45.2± 0.9 |
| RW-GNN (Nikolentzos & Vazirgiannis, ) | 88.2± 0.6 | 79.6± 0.2 | 60.5± 0.7 | 63.2± 0.5 | 61.1± 0.4 | 58.2± 0.6 | 42.9± 0.3 |
| DropGNN (Papp et al., ) | 90.3± 0.5 | 82.2± 0.3 | 61.4± 0.8 | 65.3± 0.6 | 62.9± 0.2 | 63.5± 0.7 | 46.1± 0.6 |
| IEGN (Maron et al., ) | 83.9± 0.4 | 79.3± 0.5 | 61.9± 0.4 | 60.1± 0.3 | 62.1± 0.4 | 60.7± 0.5 | 44.8± 0.8 |
| LLAMA3.1-8b (Dubey et al., ) | 67.6± 3.4 | 50.7± 3.6 | 49.6± 2.6 | 46.2± 3.8 | 42.0± 2.8 | 47.5± 2.8 | 42.2± 2.2 |
| GPT-4o (Achiam et al., ) | 73.5± 3.6 | 51.2± 0.5 | 52.7± 2.3 | 53.8± 2.9 | 48.8± 2.4 | 53.7± 1.8 | 44.5± 2.5 |
| ChemBERTa-2 (Ahmad et al., ) | 87.3± 2.7 | 77.6± 2.2 | 59.2± 1.9 | 64.8± 2.2 | 59.7± 2.8 | 59.8± 2.4 | 46.3± 2.3 |
| Logistic Regression | 58.3± 1.2 | 55.4± 0.8 | 48.4± 1.1 | 48.3± 1.0 | 48.7± 1.1 | 44.9± 1.0 | 32.5± 0.5 |
| Decision Tree (Quinlan, ) | 60.8± 1.7 | 58.6± 1.5 | 43.3± 1.0 | 46.1± 0.7 | 47.2± 0.7 | 43.5± 0.5 | 36.9± 0.8 |
| Random Forest (Breiman, ) | 64.6± 1.9 | 60.6± 1.5 | 46.9± 1.2 | 51.4± 1.5 | 51.3± 1.8 | 46.4± 1.1 | 34.8± 1.9 |
| XGBoost (Chen & Guestrin, ) | 66.9± 1.2 | 67.6± 1.4 | 51.4± 1.3 | 53.1± 1.4 | 55.8± 1.2 | 49.3± 2.1 | 38.5± 1.8 |
| w/o Calibration | 86.1± 2.2 | 74.4± 1.0 | 59.7± 2.1 | 68.9± 1.9 | 69.3± 2.7 | 61.2± 2.4 | 45.0± 2.0 |
| w/ Calibration (Ours) | 91.6± 2.0 | 83.7± 0.9 | 64.2± 1.4 | 74.4± 1.9 | 76.4± 1.8 | 68.4± 2.3 | 54.9± 2.4 |

Explainability Performance. Table 2 reports the explanation accuracy across different methods. By encoding functional group-level molecular representation, MoleX offers significantly better explainability than baselines on six datasets. Additionally, the residual calibration improves average explanation accuracy by $8.8\%$, reflecting a remarkable enhancement of explainability. Similarly, on the Mutag dataset, the explanation accuracy of logistic regression is boosted by a total of $33.4\%$ via LLM knowledge augmentation and residual calibration. Interestingly, while most methods excel on simpler datasets like Mutag but falter on complex ones like Liver, our method maintains high explanation accuracy across both, showing adaptability and high explanation quality.

**Table 2: Explanation accuracy over seven datasets (%). The best results are highlighted in bold.**
| Methods | Mutag | Mutagen | PTC-FM | PTC-FR | PTC-MM | PTC-MR | Liver |
| --- | --- | --- | --- | --- | --- | --- | --- |
| GCN (Kipf & Welling, ) | 81.1± 0.2 | 76.4± 0.2 | 65.3± 0.4 | 67.8± 0.7 | 70.8± 0.8 | 65.1± 0.2 | 62.8± 0.2 |
| DGCNN (Zhang et al., ) | 86.3± 1.2 | 87.1± 0.5 | 63.0± 1.3 | 57.0± 1.2 | 63.0± 1.3 | 62.3± 0.8 | 67.5± 1.6 |
| edGNN (Jaume et al., ) | 94.7± 0.9 | 74.4± 0.7 | 65.9± 0.5 | 64.1± 0.5 | 66.6± 0.7 | 61.4± 0.7 | 63.2± 0.3 |
| GIN (Xu et al., ) | 92.1± 0.2 | 75.6± 0.3 | 67.5± 0.6 | 69.2± 0.5 | 68.5± 0.8 | 61.3± 0.5 | 68.3± 0.9 |
| RW-GNN (Nikolentzos & Vazirgiannis, ) | 89.9± 0.6 | 76.7± 0.2 | 65.8± 0.3 | 55.5± 0.3 | 66.9± 0.1 | 59.3± 0.2 | 64.7± 0.5 |
| DropGNN (Papp et al., ) | 83.4± 0.2 | 77.4± 0.3 | 68.4± 0.2 | 64.7± 0.4 | 63.2± 0.2 | 57.4± 0.7 | 64.5± 0.8 |
| IEGN (Maron et al., ) | 82.0± 0.2 | 77.5± 0.2 | 61.6± 0.6 | 62.6± 0.9 | 69.3± 0.7 | 59.1± 0.7 | 66.6± 0.6 |
| Logistic Regression | 59.2± 0.4 | 50.6± 0.9 | 54.4± 0.3 | 47.7± 0.8 | 49.9± 0.7 | 44.3± 0.7 | 53.8± 0.7 |
| Decision Tree (Quinlan, ) | 61.2± 0.2 | 55.7± 1.0 | 56.7± 0.8 | 46.4± 1.1 | 48.1± 0.9 | 39.9± 0.8 | 56.4± 1.0 |
| Random Forest (Breiman, ) | 66.7± 1.2 | 57.2± 1.2 | 59.9± 1.7 | 50.9± 1.2 | 55.0± 0.8 | 46.6± 1.1 | 60.7± 1.4 |
| XGBoost (Chen & Guestrin, ) | 65.2± 1.2 | 61.3± 1.1 | 58.5± 1.8 | 49.4± 1.8 | 51.6± 1.3 | 50.2± 0.8 | 69.0± 1.4 |
| w/o Calibration | 90.0± 0.9 | 77.7± 1.0 | 68.0± 1.7 | 66.6± 1.1 | 62.0± 1.5 | 67.5± 1.5 | 72.0± 2.0 |
| w/ Calibration (Ours) | 92.6± 1.7 | 89.0± 1.2 | 77.9± 1.5 | 79.3± 1.4 | 72.3± 1.7 | 73.4± 1.3 | 80.3± 1.4 |

Figure: Figure 2: Explanation visualization of a molecule from the Mutag dataset (left), and the contribution scores of the identified functional groups offered by MoleX (right).
Refer to caption: x2.png

Figure 2 visualizes the explanation for a randomly selected molecule from the Mutag dataset. The ground truth, verified by domain experts, shows that mutagenicity arises from an aromatic functional group (e.g., benzene ring) bonded with another group like nitro or carbonyl. MoleX precisely identifies this ground truth substructure, faithfully explaining molecular structure-property relationships. In contrast, other methods only identify a collection of individual atoms and bonds, failing to recognize chemically meaningful substructures as a whole. For instance, PGExplainer identifies single atoms from multiple benzene rings, whereas atoms alone are insufficient to explain overall molecular properties. Notably, MoleX without calibration identifies two additional elements beyond the ground truth, thus suggesting the significance of residual calibration to explanation accuracy. Moreover, contribution scores elucidate interactions among functional groups, with the benzene-nitro substructure on the upper left receiving a high score, showcasing its importance to mutagenicity as a bonded/interacting entity. More explanation visualizations are in Section A.11.

Computational Efficiency. Figure 3 displays the inference time of different methods. Unlike methods that rely on iterative optimization in neural networks, MoleX enables considerably faster inference. Generally, MoleX outperforms both GNNs (at least $15\times$ faster) and LLMs (at least $120\times$ faster) in speed while achieving higher classification and explanation accuracy. MoleX consistently costs the least inference times across all datasets, reinforcing its scalability for real-world applications and large-scale computations on molecular data. In addition to faster inference, MoleX also significantly reduces GPU memory usage compared to baselines by avoiding numerous iterative parameter updates and storage in optimization algorithms. Consequently, the inference power of the linear model is critically augmented by LLM knowledge and residual calibration while preserving the advantage of explainability and computational efficiency.

Figure: Figure 3: Execution time across different methods over seven datasets. Ours achieves the best inference efficiency.
Refer to caption: x3.png

### 5.3 Ablation Studies

In this section, we introduce ablation studies on the number of $n$ in n-gram, principal components in EFPCA, training iterations of the residual calibrator, and the selection of the base model.

Number of $n$ in N-grams. We empirically compare the choice of $n$ in n-grams. As shown in fig. 6, the overall model performance improves as $n$ increases from 1 to 3, then declines for $n$ from 4 to 9. Three of four datasets in our studies indicate the optimal performance at $n=3$. Increasing $n$ captures more contextual semantics, including functional group interactions and raises the model performance. However, overlarge $n$ values incorporate excessive or irrelevant contextual information and reduce model utility correspondingly. Further details are in section A.10.

Dimensionality Reduction via EFPCA. We use EFPCA to reduce the dimensionality of LLM embeddings, obtaining explainable and compact embeddings. As shown in fig. 5, cross-validation across four datasets determines the optimal number of principal components. Empirically, components beyond 20 contribute minimally to the molecular property prediction. Additional components yield diminishing returns while increasing model complexity and reducing explainability. Further details are in section A.8. Moreover, we also investigate the effect of our dimensionality reduction. As presented in table 5, we compare the model performance without dimensionality reduction. We find that models using only 20 principal components achieve performance within $5\%$ of models using all components. This means selected components effectively preserve task-relevant information while excluding redundancy. Further details are offered in section A.9.

Training Iterations of the Residual Calibrator. We apply the training objective in A.11 to learn a residual calibrator that iteratively refits prediction errors. As shown in fig. 4, we observe that the model performance improves substantially with increasing training iterations until reaching a threshold. Beyond this point, the model overfits the data, leading to a performance decline. This finding suggests the need for an appropriate stopping criterion to balance model performance and prevent overfitting. Empirically, the optimal number of training iterations is 5. Further details and a theoretical demonstration are offered in section A.7.

Selection of the Base Model. Aside from the logistic regression, we examine the effect of LLM augmentation using other statistical learning models as the base model. The classification and explanation accuracy are reported in table 6 and table 7, respectively. All statistical learning models augmented with LLM knowledge and residual calibration outperform GNNs and LLMs. Besides, more complicated models, like XGBoost and random forest, achieve better performance in both classification and explanation accuracy than simple models like LASSO. Therefore, LLM knowledge is capable of augmenting a model on top of its original predictive capabilities, evidencing the effectiveness and generalizability of our method. However, model complexity generally trades off with explainability. Considering this, we select the logistic regression as our base model for its optimal balance between explainability and performance. Further details are offered in section A.12.

## 6 Conclusion

This work develops MoleX, a novel framework utilizing LLM knowledge to build a powerful linear model for accurate molecular property predictions with chemically meaningful explanations. Specifically, MoleX extracts task-relevant knowledge from LLM embeddings using information bottleneck-inspired fine-tuning and sparsity-inducing dimensionality reduction to train a linear model for explainable inference. Additionally, a residual calibration module is designed to recover the original LLM’s performance and further enhance the linear model via recapturing prediction errors. During its inference, MoleX precisely reveals crucial substructures with their interactions as explanations. Notably, MoleX enjoys the advantage of LLM’s predictive power while preserving the linear model’s intrinsic explainability. Extensive theoretical and empirical analysis demonstrate MoleX’s exceptional predictive performance, explainability, and efficiency.

## Appendix A Appendix

### A.1 Proof of N-Gram Coefficients as Valid Contribution Scores for Decoupled n-gram Features

In this section, we demonstrate that n-gram coefficients in the linear model can be interpreted as feature contribution scores based on the statistical properties of the linear model.

###### Proof.

Suppose $\mathbf{E}\in\mathbb{R}^{n\times d}$ is the matrix of n-gram embeddings, where each row $\mathbf{e}_{i}^{\top}$ is the embedding of the $i$-th n-gram. Let $\mathbf{v}_{ij}\in\mathbb{R}^{d}$ be the embedding of the $j$-th feature in the $i$-th n-gram, and suppose that each n-gram consists of $m$ features ($m$ is a constant across all n-grams). Let $c_{ij}$ denote the contribution score of the $j$-th feature in the $i$-th n-gram. We formulate the following assumptions based on OLS properties to ensure the validity of using n-gram coefficients as contribution scores:

- 1.
Linearity. The relationship between the input embeddings and the output is linear. Namely, for all $i$,
$y_{i}=\mathbf{e}_{i}^{\top}\mathbf{w}^{*}+\epsilon_{i},$
where $\mathbf{w}^{*}\in\mathbb{R}^{d}$ is the true coefficient vector, and $\epsilon_{i}$ is the error term.
- 2.
N-gram Embedding Decomposition. Each n-gram embedding $\mathbf{e}_{i}$ is the average of its constituent feature embeddings:
$\mathbf{e}_{i}=\frac{1}{m}\sum_{j=1}^{m}\mathbf{v}_{ij}.$
- 3.
Ordinary Least Squares (OLS). The linear model is estimated using OLS by minimizing the residual sum of squares:
$\hat{\mathbf{w}}=\arg\min_{\mathbf{w}}\sum_{i=1}^{n}(y_{i}-\mathbf{e}_{i}^{
\top}\mathbf{w})^{2}.$
- 4.
Error Properties.
(a)
Zero Mean Errors. The errors $\epsilon_{i}$ have zero mean given the embeddings:
$\mathbb{E}[\epsilon_{i}\mid\mathbf{E}]=0.$
(b)
Homoscedasticity. The errors have constant variance given the embeddings:
$\operatorname{Var}[\epsilon_{i}\mid\mathbf{E}]=\sigma^{2},$
where $\sigma^{2}$ is a positive constant.
(c)
No Autocorrelation. The errors are uncorrelated with each other:
$\operatorname{Cov}[\epsilon_{i},\epsilon_{j}\mid\mathbf{E}]=0\quad\text{for }i
\neq j.$
- 5.
Full Rank. The matrix $\mathbf{E}^{\top}\mathbf{E}$ is invertible (i.e., $\mathbf{E}$ has full column rank).

We define the contribution score of each decoupled n-gram feature as follows:

###### Definition A.1 .

The *feature contribution score* $c_{ij}$ for the $j$-th feature in the $i$-th n-gram is defined as

$$ $c_{ij}=\mathbf{v}_{ij}^{\top}\hat{\mathbf{w}},$ $$

where $\hat{\mathbf{w}}$ is the estimated coefficient vector from the linear model.

###### Lemma A.1 (Prediction as Sum of Feature Contributions) .

Under Assumption 2, the predicted output for the $i$-th n-gram is

$$ $\hat{y}_{i}=\mathbf{e}_{i}^{\top}\hat{\mathbf{w}}=\frac{1}{m}\sum_{j=1}^{m}c_{ ij}.$ $$

###### Proof.

Using the embedding decomposition and the definition of the contribution scores, we have

$$ $\displaystyle\hat{y}_{i}$ $\displaystyle=\mathbf{e}_{i}^{\top}\hat{\mathbf{w}}$ $\displaystyle=\left(\frac{1}{m}\sum_{j=1}^{m}\mathbf{v}_{ij}\right)^{\top}\hat {\mathbf{w}}$ $\displaystyle=\frac{1}{m}\sum_{j=1}^{m}\mathbf{v}_{ij}^{\top}\hat{\mathbf{w}}$ $\displaystyle=\frac{1}{m}\sum_{j=1}^{m}c_{ij}.$ $$

This completes the proof.
∎

###### Theorem A.2 (Contribution Scores Quantify Individual Feature Contributions) .

Under the Linearity assumption (Assumption 1), the feature contribution scores $c_{ij}$ quantify the contributions of individual features to the prediction $\hat{y}_{i}$.

###### Proof.

From Lemma A.1, the predicted value $\hat{y}_{i}$ is given as the average of the feature contribution scores $c_{ij}$. Specifically,

$$ $\hat{y}_{i}=\frac{1}{m}\sum_{j=1}^{m}c_{ij}.$ $$

This equation shows that each feature’s contribution score $c_{ij}$ directly influences the prediction $\hat{y}_{i}$. Therefore, $c_{ij}$ quantifies the contribution of the $j$-th feature in the $i$-th n-gram to the prediction.

This completes the proof.
∎

Due to the statistical properties of the OLS estimator, we formulate the following theorem:

###### Theorem A.3 (Properties of the OLS Estimator) .

Under Assumptions 1–5, the OLS estimator $\hat{\mathbf{w}}$ satisfies:

- 1.
Unbiasedness. $\mathbb{E}[\hat{\mathbf{w}}\mid\mathbf{E}]=\mathbf{w}^{*}$.
- 2.
Variance-Covariance Matrix. $\operatorname{Var}[\hat{\mathbf{w}}\mid\mathbf{E}]=\sigma^{2}(\mathbf{E}^{\top
}\mathbf{E})^{-1}$.
- 3.
Consistency. As $n\to\infty$, $\hat{\mathbf{w}}\xrightarrow{P}\mathbf{w}^{*}$.

###### Proof.

We prove each property as follows.

(1) Unbiasedness: The OLS estimator is given by

$$ $\hat{\mathbf{w}}=(\mathbf{E}^{\top}\mathbf{E})^{-1}\mathbf{E}^{\top}\mathbf{y}.$ $$

Substituting $\mathbf{y}=\mathbf{E}\mathbf{w}^{*}+\bm{\epsilon}$, we have

$$ $\hat{\mathbf{w}}=\mathbf{w}^{*}+(\mathbf{E}^{\top}\mathbf{E})^{-1}\mathbf{E}^{ \top}\bm{\epsilon}.$ $$

Taking expectations conditional on $\mathbf{E}$ and using Assumption 4(a),

$$ $\mathbb{E}[\hat{\mathbf{w}}\mid\mathbf{E}]=\mathbf{w}^{*}+(\mathbf{E}^{\top} \mathbf{E})^{-1}\mathbf{E}^{\top}\mathbb{E}[\bm{\epsilon}\mid\mathbf{E}]= \mathbf{w}^{*}.$ $$

(2) Variance-Covariance Matrix: The variance conditional on $\mathbf{E}$ is

$$ $\displaystyle\operatorname{Var}[\hat{\mathbf{w}}\mid\mathbf{E}]$ $\displaystyle=\operatorname{Var}\left((\mathbf{E}^{\top}\mathbf{E})^{-1} \mathbf{E}^{\top}\bm{\epsilon}\mid\mathbf{E}\right)$ $\displaystyle=(\mathbf{E}^{\top}\mathbf{E})^{-1}\mathbf{E}^{\top}\operatorname {Var}[\bm{\epsilon}\mid\mathbf{E}]\mathbf{E}(\mathbf{E}^{\top}\mathbf{E})^{-1}$ $\displaystyle=\sigma^{2}(\mathbf{E}^{\top}\mathbf{E})^{-1},$ $$

using Assumptions 4(b) and (c).

(3) Consistency: As $n\to\infty$, under the Law of Large Numbers,

$$ $\frac{1}{n}\mathbf{E}^{\top}\mathbf{E}\xrightarrow{P}\mathbf{Q},$ $$

where $\mathbf{Q}$ is positive definite due to Assumption 5. Additionally,

$$ $\frac{1}{n}\mathbf{E}^{\top}\bm{\epsilon}\xrightarrow{P}\mathbf{0},$ $$

since $\bm{\epsilon}$ has zero mean and finite variance. Therefore,

$$ $\hat{\mathbf{w}}=\mathbf{w}^{*}+(\mathbf{E}^{\top}\mathbf{E})^{-1}\mathbf{E}^{ \top}\bm{\epsilon}\xrightarrow{P}\mathbf{w}^{*}.$ $$

This completes the proof.
∎

To validate the convergence of the contribution scores, we introduce the asymptotic normality of the OLS estimator.

###### Corollary A.1 (Asymptotic Normality) .

If the error terms $\bm{\epsilon}$ are independently and identically normally distributed with mean zero and variance $\sigma^{2}$, then we have

$$ $\sqrt{n}(\hat{\mathbf{w}}-\mathbf{w}^{*})\xrightarrow{d}\mathcal{N}\left( \mathbf{0},\sigma^{2}\mathbf{Q}^{-1}\right),$ $$

where $\mathbf{Q}=\lim_{n\to\infty}\frac{1}{n}\mathbf{E}^{\top}\mathbf{E}$.

###### Proof.

Under the given conditions, the Central Limit Theorem applies to the sum $\mathbf{E}^{\top}\bm{\epsilon}$. Specifically,

$$ $\sqrt{n}(\hat{\mathbf{w}}-\mathbf{w}^{*})=(\mathbf{E}^{\top}\mathbf{E})^{-1} \mathbf{E}^{\top}\bm{\epsilon}=\left(\frac{1}{n}\mathbf{E}^{\top}\mathbf{E} \right)^{-1}\left(\frac{1}{\sqrt{n}}\mathbf{E}^{\top}\bm{\epsilon}\right).$ $$

As $n\to\infty$, $\frac{1}{n}\mathbf{E}^{\top}\mathbf{E}\xrightarrow{P}\mathbf{Q}$ and $\frac{1}{\sqrt{n}}\mathbf{E}^{\top}\bm{\epsilon}\xrightarrow{d}\mathcal{N}(
\mathbf{0},\sigma^{2}\mathbf{Q})$. Therefore,

$$ $\sqrt{n}(\hat{\mathbf{w}}-\mathbf{w}^{*})\xrightarrow{d}\mathcal{N}\left( \mathbf{0},\sigma^{2}\mathbf{Q}^{-1}\right).$ $$

This completes the proof.
∎

###### Lemma A.4 (Variance of $\hat{c}_{ij}$ ) .

The variance of the estimated feature contribution score $\hat{c}_{ij}=\mathbf{v}_{ij}^{\top}\hat{\mathbf{w}}$ is

$$ $\operatorname{Var}[\hat{c}_{ij}\mid\mathbf{E}]=\sigma^{2}\mathbf{v}_{ij}^{\top }(\mathbf{E}^{\top}\mathbf{E})^{-1}\mathbf{v}_{ij}.$ $$

###### Proof.

Since $\hat{c}_{ij}$ is a linear function of $\hat{\mathbf{w}}$, its variance conditional on $\mathbf{E}$ is

$$ $\displaystyle\operatorname{Var}[\hat{c}_{ij}\mid\mathbf{E}]$ $\displaystyle=\operatorname{Var}\left(\mathbf{v}_{ij}^{\top}\hat{\mathbf{w}} \mid\mathbf{E}\right)$ $\displaystyle=\mathbf{v}_{ij}^{\top}\operatorname{Var}[\hat{\mathbf{w}}\mid \mathbf{E}]\mathbf{v}_{ij}$ $\displaystyle=\sigma^{2}\mathbf{v}_{ij}^{\top}(\mathbf{E}^{\top}\mathbf{E})^{- 1}\mathbf{v}_{ij},$ $$

using the result from Theorem A.3(2).

This completes the proof.
∎

Finally, we demonstrate the statistical significance of the feature contribution scores based on the n-gram coefficients.

###### Theorem A.5 (t-Statistic for Feature Contribution Scores) .

Under the above assumptions, the t-statistic for testing $H_{0}:c_{ij}=0$ is given as

$$ $t_{ij}=\frac{\hat{c}_{ij}}{\operatorname{SE}[\hat{c}_{ij}]}=\frac{\mathbf{v}_{ ij}^{\top}\hat{\mathbf{w}}}{\sigma\sqrt{\mathbf{v}_{ij}^{\top}(\mathbf{E}^{ \top}\mathbf{E})^{-1}\mathbf{v}_{ij}}}.$ $$

###### Proof.

The standard error of $\hat{c}_{ij}$ is

$$ $\operatorname{SE}[\hat{c}_{ij}]=\sqrt{\operatorname{Var}[\hat{c}_{ij}\mid \mathbf{E}]}=\sigma\sqrt{\mathbf{v}_{ij}^{\top}(\mathbf{E}^{\top}\mathbf{E})^{ -1}\mathbf{v}_{ij}}.$ $$

Therefore, the t-statistic is

$$ $t_{ij}=\frac{\hat{c}_{ij}}{\operatorname{SE}[\hat{c}_{ij}]}.$ $$

Under the null hypothesis $H_{0}:c_{ij}=0$ and the assumption of normality, $t_{ij}$ follows a t-distribution with $n-d$ degrees of freedom.

This completes the proof.
∎

From Theorem A.2, we have shown that the feature contribution scores $c_{ij}$ represent the contributions of individual features to the predictions $\hat{y}_{i}$. The statistical properties outlined in Theorem A.3 and Lemma A.4 guarantee that these estimates are reliable and their statistical significance can be assessed.

Therefore, we prove that each feature’s contribution to the prediction can be quantified by its corresponding coefficient in the linear model, enable us to assess the importance of individual features. By mathematically linking the model coefficients to the feature contributions, we validate the use of these coefficients as measures of feature importance. We also conclude that using n-gram coefficients derived from feature embeddings and model coefficients as contribution scores for input features is valid and grounded in the statistical properties of the linear model. The n-gram coefficients used in our implementation faithfully reflect the importance and contribution of each input feature in our linear model. By expressing the predicted output as the sum of individual feature contributions, we effectively decouple the influence of each feature/functional group on the output/molecular property. This decoupling allows us to isolate the effect of each n-gram feature or functional group $x$ on the molecular property $y$. Consequently, the contribution scores $c_{ij}$ provide a quantitative measure of how each functional group impacts the molecular property.

This completes the proof.
∎

###### Proof.

###### Definition A.1 .

###### Lemma A.1 (Prediction as Sum of Feature Contributions) .

###### Proof.

###### Theorem A.2 (Contribution Scores Quantify Individual Feature Contributions) .

###### Proof.

###### Theorem A.3 (Properties of the OLS Estimator) .

###### Proof.

###### Corollary A.1 (Asymptotic Normality) .

###### Proof.

###### Lemma A.4 (Variance of c ^ i ⁢ j subscript ^ 𝑐 𝑖 𝑗 \hat{c}_{ij} over^ start_ARG italic_c end_ARG start_POSTSUBSCRIPT italic_i italic_j end_POSTSUBSCRIPT ) .

###### Proof.

###### Theorem A.5 (t-Statistic for Feature Contribution Scores) .

###### Proof.

### A.2 Proof of theorem 4.1 (Explainability of VIB-based Training Objectives)

###### Proof.

We demonstrate the Variational Information Bottleneck (VIB) framework, which aims to learn a compressed representation $Z$ of input variable $X$ that preserves maximal information about the target variable $Y$ while being minimally informative about $X$ itself. This is achieved by optimizing the objective function as follows:

$$ $\mathcal{L}_{\text{IB}}(\theta)=I(Z;X)-\beta I(Z;Y)$ $$

where $I(\cdot;\cdot)$ is mutual information, $\beta\geq 0$ is a tuning parameter, and $\theta$ is the parameters of the encoder. Our goal is to derive a tractable variational lower bound of this objective function that can be optimized using stochastic gradient descent.

###### Definition A.2 (Mutual Information) .

For random variables $X$ and $Z$ with joint distribution $p(X,Z)$, the mutual information $I(X;Z)$ is defined as

$$ $I(X;Z)=\mathbb{E}_{p(X,Z)}\left[\log\frac{p(X,Z)}{p(X)p(Z)}\right]$ $$

Alternatively, it can be expressed as

$$ $I(X;Z)=\mathbb{E}_{p(X)}\left[D_{\mathrm{KL}}(p(Z|X)\|p(Z))\right]$ $$

###### Definition A.3 (Kullback-Leibler Divergence) .

For probability distributions $P$ and $Q$ over the same probability space, the KL divergence from $Q$ to $P$ is defined as

$$ $D_{\mathrm{KL}}(P\|Q)=\int p(x)\log\frac{p(x)}{q(x)}\,dx=\mathbb{E}_{p(x)} \left[\log\frac{p(x)}{q(x)}\right]$ $$

###### Definition A.4 (Conditional Entropy) .

The conditional entropy $H(Y|Z)$ is defined as

$$ $H(Y|Z)=-\mathbb{E}_{p(Z,Y)}\left[\log p(Y|Z)\right]$ $$

We then formulate the problem. Let $\mathcal{D}=\{(X_{i},Y_{i})\}_{i=1}^{N}$ be a dataset of input-output pairs sampled from an unknown distribution $p(X,Y)$. The encoder $p_{\theta}(Z|X)$ parameterizes the conditional distribution of $Z$ given $X$, and the decoder $q_{\phi}(Y|Z)$ parameterizes the conditional distribution of $Y$ given $Z$. Our objective is to optimize the parameters $\theta$ and $\phi$ by maximizing the Information Bottleneck Lagrangian as follows:

$$ $\mathcal{L}_{\text{IB}}(\theta,\phi)=I(Z;Y)-\beta I(Z;X)$ $$

However, direct computation of $I(Z;Y)$ and $I(Z;X)$ is intractable. Therefore, we derive variational bounds to make the optimization objective tractable. We start by applying the lemma as

###### Lemma A.6 (Variational Upper Bound on $I(Z;X)$ ) .

The mutual information $I(Z;X)$ can be upper-bounded as

$$ $I(Z;X)\leq\mathbb{E}_{p(X)}\left[D_{\mathrm{KL}}(p_{\theta}(Z|X)\|r(Z))\right]$ $$

where $r(Z)$ is an arbitrary prior distribution over $Z$.

###### Proof.

We start by rewriting $I(Z;X)$ as

$$ $I(Z;X)=\mathbb{E}_{p(X)}\left[D_{\mathrm{KL}}(p_{\theta}(Z|X)\|p(Z))\right]$ $$

Since $p(Z)=\int p_{\theta}(Z|X)p(X)\,dX$, and $p(Z)$ is intractable, we introduce an approximate prior $r(Z)$ and apply the following decomposition:

$$ $\displaystyle D_{\mathrm{KL}}(p_{\theta}(Z|X)\|p(Z))$ $\displaystyle=D_{\mathrm{KL}}(p_{\theta}(Z|X)\|r(Z))-D_{\mathrm{KL}}(p_{\theta }(Z|X)\|r(Z))+D_{\mathrm{KL}}(p_{\theta}(Z|X)\|p(Z))$ $\displaystyle=D_{\mathrm{KL}}(p_{\theta}(Z|X)\|r(Z))-D_{\mathrm{KL}}(p(Z)\|r(Z))$ $$

Thus, we have:

$$ $I(Z;X)=\mathbb{E}_{p(X)}\left[D_{\mathrm{KL}}(p_{\theta}(Z|X)\|r(Z))\right]-D_ {\mathrm{KL}}(p(Z)\|r(Z))$ $$

Since $D_{\mathrm{KL}}(p(Z)\|r(Z))\geq 0$, it follows that:

$$ $I(Z;X)\leq\mathbb{E}_{p(X)}\left[D_{\mathrm{KL}}(p_{\theta}(Z|X)\|r(Z))\right]$ $$

This completes the proof.
∎

###### Lemma A.7 (Variational Lower Bound on $I(Z;Y)$ ) .

The mutual information $I(Z;Y)$ can be lower-bounded as

$$ $I(Z;Y)\geq\mathbb{E}_{p(X,Y)}\left[\mathbb{E}_{p_{\theta}(Z|X)}\left[\log q_{ \phi}(Y|Z)\right]\right]-H(Y)$ $$

###### Proof.

By the definition of mutual information:

$$ $I(Z;Y)=H(Y)-H(Y|Z)=H(Y)+\mathbb{E}_{p(Z,Y)}\left[\log p(Y|Z)\right]$ $$

Since $p(Y|Z)$ is intractable, we introduce a variational approximation $q_{\phi}(Y|Z)$ to get:

$$ $\mathbb{E}_{p(Z,Y)}\left[\log p(Y|Z)\right]=\mathbb{E}_{p(X,Y)}\left[\mathbb{E }_{p_{\theta}(Z|X)}\left[\log p(Y|Z)\right]\right]$ $$

Using Jensen’s inequality and the non-negativity of KL divergence, we have:

$$ $\mathbb{E}_{p_{\theta}(Z|X)}\left[\log p(Y|Z)\right]\geq\mathbb{E}_{p_{\theta} (Z|X)}\left[\log q_{\phi}(Y|Z)\right]-D_{\mathrm{KL}}(p_{\theta}(Z|X)\|q_{\phi }(Y|Z))$ $$

Neglecting the KL divergence term (assuming $q_{\phi}(Y|Z)$ approximates $p(Y|Z)$ well), we have:

$$ $I(Z;Y)\geq H(Y)+\mathbb{E}_{p(X,Y)}\left[\mathbb{E}_{p_{\theta}(Z|X)}\left[ \log q_{\phi}(Y|Z)\right]\right]$ $$

Thus:

$$ $I(Z;Y)\geq\mathbb{E}_{p(X,Y)}\left[\mathbb{E}_{p_{\theta}(Z|X)}\left[\log q_{ \phi}(Y|Z)\right]\right]-H(Y)$ $$

This completes the proof.
∎

Now we can formulate the Variational Information Bottleneck (VIB) Objective. Particularly, combining Lemmas lemma A.6 and lemma A.7, we obtain a tractable objective function.

###### Proposition A.8 (Variational Lower Bound on the Information Bottleneck Objective) .

The Information Bottleneck Lagrangian can be lower-bounded by the variational objective function:

$$ $\mathcal{L}(\theta,\phi)=\mathbb{E}_{p(X,Y)}\left[\mathbb{E}_{p_{\theta}(Z|X)} \left[-\log q_{\phi}(Y|Z)\right]+\beta D_{\mathrm{KL}}(p_{\theta}(Z|X)\|r(Z))\right]$ $$

###### Proof.

Starting from the original objective:

$$ $\mathcal{L}_{\text{IB}}(\theta,\phi)=I(Z;X)-\beta I(Z;Y)$ $$

Applying the upper bound of $I(Z;X)$ from Lemma lemma A.6 and the lower bound of $I(Z;Y)$ from Lemma lemma A.7, we get:

$$ $\displaystyle\mathcal{L}_{\text{IB}}(\theta,\phi)$ $\displaystyle\leq\mathbb{E}_{p(X)}\left[D_{\mathrm{KL}}(p_{\theta}(Z|X)\|r(Z)) \right]-\beta\left(\mathbb{E}_{p(X,Y)}\left[\mathbb{E}_{p_{\theta}(Z|X)}\left[ \log q_{\phi}(Y|Z)\right]\right]-H(Y)\right)$ $\displaystyle=\mathbb{E}_{p(X)}\left[D_{\mathrm{KL}}(p_{\theta}(Z|X)\|r(Z)) \right]+\beta H(Y)-\beta\mathbb{E}_{p(X,Y)}\left[\mathbb{E}_{p_{\theta}(Z|X)} \left[\log q_{\phi}(Y|Z)\right]\right]$ $$

Since $H(Y)$ is constant with respect to $\theta$ and $\phi$, we can ignore it for optimization purposes. We define the variational objective function as

$$ $\mathcal{L}(\theta,\phi)=\mathbb{E}_{p(X,Y)}\left[\mathbb{E}_{p_{\theta}(Z|X)} \left[-\log q_{\phi}(Y|Z)\right]+\beta D_{\mathrm{KL}}(p_{\theta}(Z|X)\|r(Z))\right]$ $$

which serves as an upper bound on $\mathcal{L}_{\text{IB}}(\theta,\phi)$. Minimizing $\mathcal{L}$ will therefore minimize $\mathcal{L}_{\text{IB}}$, satisfying our optimization goal.

This completes the proof.
∎

In our fine-tuning stage, since the expectation of $p(X,Y)$ is approximated by empirical samples from the dataset $\mathcal{D}$, and the expectations of $p_{\theta}(Z|X)$ are approximated by Monte Carlo sampling with reparameterization trick. Thus, the loss function is expressed as (note that this is a generalized form of our designed loss function shown in eq. 4.1):

$$ $\hat{\mathcal{L}}(\theta,\phi)=\frac{1}{N}\sum_{i=1}^{N}\left(-\mathbb{E}_{p_{ \theta}(Z|X_{i})}\left[\log q_{\phi}(Y_{i}|Z)\right]+\beta D_{\mathrm{KL}}(p_{ \theta}(Z|X_{i})\|r(Z))\right)$ $$

To demonstrate the convergence, we have:

###### Theorem A.9 (Convergence of Stochastic Gradient Descent) .

Under standard assumptions of stochastic optimization (bounded gradients, appropriate learning rates, etc.), stochastic gradient descent (SGD) converges to a local minimum of $\hat{\mathcal{L}}(\theta,\phi)$.

###### Proof.

While the neural network training is non-convex, empirical and theoretical results in optimization, SGD can converge to critical points (which may be local minima, maxima, or saddle points) provided the loss function is smooth (i.e., continuously differentiable) and the gradients are Lipschitz continuous.

Given that $\hat{\mathcal{L}}(\theta,\phi)$ is composed of differentiable functions, and the gradients with respect to $\theta$ and $\phi$ can be computed via back-propagation, the convergence to a local minimum is attainable under proper value of learning rate and optimization parameters.

This completes the proof.
∎

We express the corollary of our learned molecular representation after fine-tuning as

###### Corollary A.2 (Informative and Compressed Molecular Representation) .

At convergence, the learned representation $Z$ satisfies:

$$ $I(Z;Y)\text{ is maximized, and }I(Z;X)\text{ is minimized (subject to the tuning parameter }\beta)$ $$

###### Proof.

By optimizing the variational objective function $\hat{\mathcal{L}}(\theta,\phi)$, we are effectively minimizing an upper bound on $I(Z;X)$ (Lemma lemma A.6) and maximizing a lower bound on $I(Z;Y)$ (Lemma lemma A.7). The trade-off between the two objectives is controlled by $\beta$.

As $\beta$ increases, more emphasis is placed on minimizing $I(Z;X)$, leading to a more compressed representation $Z$ that preserves only the most task-relevant information about $Y$.

This completes the proof.
∎

Specifically, as the first term in the loss function encourages the embeddings $t$ to be highly predictive of $y$, it intrinsically captures the task-relevant information. Meanwhile, the second term penalizes the complexity of $t$ by forcing it to be close to the prior $p_{0}(t)$, thereby excluding unnecessary information from $x$. These objectives ensure that the embeddings are both task-relevant and compact, containing minimal spurious data. Additionally, through the derivation of variational bounds and the construction of a tractable objective function, we have shown that minimizing $\mathcal{L}(\theta,\phi)$ allows us to learn a molecular representation $Z$ that captures maximal information about $Y$ while being minimally informative about $X$, in accordance with the Information Bottleneck principle. The optimization of $\mathcal{L}$ via SGD converges to a local minimum under standard optimization assumptions. Therefore, we learn an informative embedding after fine-tuning the pre-trained LLM, and we thus can extract the embedding with improved informativeness.

In conclusion, by framing the fine-tuning within the VIB framework, we derive this approach that balances the essential information for property prediction $y$ with the elimination of irrelevant details from the input molecular representation $x$. This theoretical foundation ensures that MoleX effectively focuses on extracting the most relevant features needed for accurate predictions.

This completes the proof.
∎

###### Proof.

###### Definition A.2 (Mutual Information) .

###### Definition A.3 (Kullback-Leibler Divergence) .

###### Definition A.4 (Conditional Entropy) .

###### Lemma A.6 (Variational Upper Bound on I ⁢ ( Z ; X ) 𝐼 𝑍 𝑋 I(Z;X) italic_I ( italic_Z ; italic_X ) ) .

###### Proof.

###### Lemma A.7 (Variational Lower Bound on I ⁢ ( Z ; Y ) 𝐼 𝑍 𝑌 I(Z;Y) italic_I ( italic_Z ; italic_Y ) ) .

###### Proof.

###### Proposition A.8 (Variational Lower Bound on the Information Bottleneck Objective) .

###### Proof.

###### Theorem A.9 (Convergence of Stochastic Gradient Descent) .

###### Proof.

###### Corollary A.2 (Informative and Compressed Molecular Representation) .

###### Proof.

### A.3 Proof of theorem 4.2 (Explainability of EFPCA)

###### Proof.

To demonstrate the explainability of the EFPCA, we will showcase how the incorporation of a sparsity-inducing penalty and the use of basis functions with local support lead to explainable FPCs.

We first formulated the EFPCA as an optimization problem. The EFPCA seeks to find FPCs $\xi_{k}(t)$ by solving the following optimization

$$ $\max_{\xi_{k}}\left\{\langle\xi_{k},\hat{\mathcal{C}}\xi_{k}\rangle-\rho_{k} \mathcal{S}(\xi_{k})\right\}$ (A.1) $$

subject to

$$ $\|\xi_{k}\|_{\gamma}^{2}=\|\xi_{k}\|^{2}+\gamma\|\mathcal{D}^{2}\xi_{k}\|^{2}=1$ (A.2) $$

and

$$ $\langle\xi_{k},\hat{\xi}_{j}\rangle_{\gamma}=0\quad\text{for all }j<k$ (A.3) $$

where

- •
$\hat{\mathcal{C}}$ is the empirical covariance operator of the centered stochastic process $X(t)-\mu(t)$.
- •
$\langle f,g\rangle=\int_{a}^{b}f(t)g(t)\,dt$ is the standard $L^{2}$ inner product.
- •
$\|\xi_{k}\|^{2}=\langle\xi_{k},\xi_{k}\rangle$ is the squared $L^{2}$ norm.
- •
$\mathcal{D}^{2}\xi_{k}$ denotes the second derivative of $\xi_{k}(t)$.
- •
$\|\mathcal{D}^{2}\xi_{k}\|^{2}=\langle\mathcal{D}^{2}\xi_{k},\mathcal{D}^{2}
\xi_{k}\rangle$ penalizes the roughness of $\xi_{k}(t)$.
- •
$\gamma>0$ is a smoothing parameter balancing the trade-off between variance explanation and smoothness.
- •
$\langle f,g\rangle_{\gamma}=\langle f,g\rangle+\gamma\langle\mathcal{D}^{2}f,
\mathcal{D}^{2}g\rangle$ is the roughness-penalized inner product.
- •
$\mathcal{S}(\xi_{k})=\int_{a}^{b}\mathbf{1}_{\{\xi_{k}(t)\neq 0\}}\,dt$ measures the support length of $\xi_{k}(t)$.
- •
$\rho_{k}>0$ is a tuning parameter controlling the sparsity of $\xi_{k}(t)$.

We then construct the expansion in basis functions with local Support. In particularly, we choose a set of basis functions $\{\phi_{j}(t)\}_{j=1}^{p}$ that have local support on the interval $[a,b]$, such as B-spline basis functions. Each $\phi_{j}(t)$ is nonzero only over a subinterval $S_{j}\subset[a,b]$. We expand $\xi_{k}(t)$ in terms of these basis functions:

$$ $\xi_{k}(t)=\sum_{j=1}^{p}a_{kj}\phi_{j}(t)$ (A.4) $$

where $a_{k}=(a_{k1},a_{k2},\dots,a_{kp})^{\top}$ is the coefficient vector for the $k$-th principal component. Therefore, the reformulation of the optimization problem in terms of coefficients can be expressed. We substitute the expansion (A.4) into the optimization problem(A.1), and express the objective function and constraints in terms of $a_{k}$:

$$ $\max_{a_{k}}\left\{a_{k}^{\top}\mathbf{Q}a_{k}-\rho_{k}\|a_{k}\|_{0}\right\}$ (A.5) $$

subject to

$$ $a_{k}^{\top}\mathbf{G}a_{k}=1$ (A.6) $$

and

$$ $a_{k}^{\top}\mathbf{G}\hat{a}_{j}=0\quad\text{for all }j<k$ (A.7) $$

where $\mathbf{Q}\in\mathbb{R}^{p\times p}$ is the matrix with entries $Q_{ij}=\langle\phi_{i},\hat{\mathcal{C}}\phi_{j}\rangle$, $\mathbf{G}\in\mathbb{R}^{p\times p}$ is the roughness-penalized Gram matrix with entries $G_{ij}=\langle\phi_{i},\phi_{j}\rangle+\gamma\langle\mathcal{D}^{2}\phi_{i},
\mathcal{D}^{2}\phi_{j}\rangle$, $\|a_{k}\|_{0}=\sum_{j=1}^{p}\mathbf{1}_{\{a_{kj}\neq 0\}}$ counts the number of nonzero coefficients in $a_{k}$, and $\hat{a}_{j}$ are the coefficient vectors of previously computed principal components.

Through the introduction of sparsity-inducing $\ell_{0}$ penalty, we can observe that the term $\rho_{k}\|a_{k}\|_{0}$ in the objective function (A.5) is a sparsity-inducing penalty that encourages many coefficients $a_{kj}$ to be exactly zero. When $\rho_{k}$ is sufficiently large, the optimization process favors solutions where $\|a_{k}\|_{0}$ is minimized, effectively shrinking less significant coefficients to zero.

Furthermore, define the index set of nonzero coefficients:

$$ $\mathcal{I}_{k}=\{j\mid a_{kj}\neq 0\}$ (A.8) $$

The principal component $\xi_{k}(t)$ can then be expressed as:

$$ $\xi_{k}(t)=\sum_{j\in\mathcal{I}_{k}}a_{kj}\phi_{j}(t)$ (A.9) $$

Since each $\phi_{j}(t)$ has support only on $S_{j}$, the support of $\xi_{k}(t)$ is:

$$ $\operatorname{supp}(\xi_{k})=\bigcup_{j\in\mathcal{I}_{k}}S_{j}$ (A.10) $$

This means $\xi_{k}(t)$ is exactly zero outside the union of the supports of the basis functions corresponding to nonzero coefficients. The localization of $\xi_{k}(t)$ to specific intervals $S_{j}$ directly indicates regions where the functional data exhibits significant variation. The sparsity in the coefficients $a_{k}$ leads to sparsity in $\xi_{k}(t)$, making the principal components explainable:

- •
Identification of Significant Intervals. The nonzero coefficients correspond to basis functions whose supports cover intervals where the data carries significant information.
- •
Reduction of Complexity. By focusing on a subset of intervals, $\xi_{k}(t)$ simplifies the understanding of functional variability.
- •
Improved Explainability. We can easily interpret $\xi_{k}(t)$ by examining the intervals where it is nonzero, relating these regions back to the original context of the data.

Besides, we introduce the roughness penalty for smoothness. The roughness penalty $\gamma\|\mathcal{D}^{2}\xi_{k}\|^{2}$ in the normalization constraint (A.2) and the inner product $\langle\cdot,\cdot\rangle_{\gamma}$ ensures that $\xi_{k}(t)$ remains smooth within its support. This avoids overfitting and maintains the functional integrity of the PCs.

In conclusion, by incorporating a sparsity-inducing $\ell_{0}$ penalty and utilizing basis functions with local support, the EFPCA produces FPCs that are both sparse and localized. This leads to principal components $\xi_{k}(t)$ that are nonzero only over intervals where the data contains significant variation, making them intrinsically explainable. By effectively encouraging sparsity, the EFPCA filters out task-irrelevant information in the high-dimensional LLM embeddings, preserving only the most informative features associated with the task. Consequently, the dimensionality reduction achieved is efficient and designed to exclude noisy and redundant information in the embeddings. Therefore, the EFPCA achieves explainability through its optimization framework that deliver sparsity, assisting dimensionality reduction while providing insights into the underlying functional data. The EFPCA can effectively reduce the dimensionality of the extracted LLM embeddings as the solution of the objective function can be optimized to become sparser. In our implementation, the statistically significant features are kept in an explainable manner.

This completes the proof.
∎

###### Proof.

### A.4 Proof of theorem 4.3 (Explainability of Residual Calibration)

###### Proof.

We demonstrate that the residual calibrator $r$ is explainable when combined with the explainable linear model $h$, under certain conditions.

Let $\mathcal{X}$ be the input space, and $\mathcal{Y}$ be the output space. Let $f:\mathcal{X}\rightarrow\mathbb{H}$ be a feature extraction function mapping inputs $x\in\mathcal{X}$ to a Hilbert space $\mathbb{H}$ of dimension $d$, equipped with inner product $\langle\cdot,\cdot\rangle_{\mathbb{H}}$ and norm $\|\cdot\|_{\mathbb{H}}$. Assume that $f(x)$ can be decomposed into two orthogonal components via the direct sum decomposition of $\mathbb{H}$:

$$ $\mathbb{H}=\mathbb{H}_{C}\oplus\mathbb{H}_{R},$ $$

where $\mathbb{H}_{C}$ and $\mathbb{H}_{R}$ are closed linear subspaces of $\mathbb{H}$, satisfying the following conditions:

$$ $\mathbb{H}_{C}\cap\mathbb{H}_{R}=\{0\},\quad\text{and}\quad\forall u\in\mathbb {H}_{C},\ \forall v\in\mathbb{H}_{R},\ \langle u,v\rangle_{\mathbb{H}}=0.$ $$

Then, for each $x\in\mathcal{X}$, we have:

$$ $f(x)=f_{H}(x)+f_{R}(x),$ $$

where $f_{H}(x)\in\mathbb{H}_{C}$ is the projection of $f(x)$ onto $\mathbb{H}_{C}$, and $f_{R}(x)\in\mathbb{H}_{R}$ is the projection onto $\mathbb{H}_{R}$. The projections are given by the orthogonal projection operators $P_{C}:\mathbb{H}\rightarrow\mathbb{H}_{C}$ and $P_{R}:\mathbb{H}\rightarrow\mathbb{H}_{R}$, such that:

$$ $f_{H}(x)=P_{C}f(x),\quad f_{R}(x)=P_{R}f(x),\quad\text{and}\quad P_{C}+P_{R}=I _{\mathbb{H}},$ $$

where $I_{\mathbb{H}}$ is the identity operator on $\mathbb{H}$. Let $h:\mathbb{H}_{C}\rightarrow\mathbb{R}$ be an explainable model, assumed to be linear:

$$ $h(f_{H}(x))=\langle w_{h},f_{H}(x)\rangle_{\mathbb{H}}+b_{h},$ $$

where $w_{h}\in\mathbb{H}_{C}$ is the weight vector, and $b_{h}\in\mathbb{R}$ is the bias term. Similarly, let $r:\mathbb{H}_{R}\rightarrow\mathbb{R}$ be the residual calibrator, also linear:

$$ $r(f_{R}(x))=\langle w_{r},f_{R}(x)\rangle_{\mathbb{H}}+b_{r},$ $$

where $w_{r}\in\mathbb{H}_{R}$ and $b_{r}\in\mathbb{R}$. Then the combined model’s prediction is:

$$ $\hat{y}(x)=h(f_{H}(x))+r(f_{R}(x))=\langle w_{h},f_{H}(x)\rangle_{\mathbb{H}}+ \langle w_{r},f_{R}(x)\rangle_{\mathbb{H}}+b_{h}+b_{r}.$ $$

Define $w=w_{h}+w_{r}\in\mathbb{H}$ and $b=b_{h}+b_{r}$, so that we have:

$$ $\hat{y}(x)=\langle w_{h},f_{H}(x)\rangle_{\mathbb{H}}+\langle w_{r},f_{R}(x) \rangle_{\mathbb{H}}+b=\langle w,f(x)\rangle_{\mathbb{H}}+b.$ $$

Since $w_{h}\in\mathbb{H}_{C}$, $w_{r}\in\mathbb{H}_{R}$, and $\mathbb{H}_{C}\perp\mathbb{H}_{R}$, it follows that:

$$ $\langle w_{h},f_{R}(x)\rangle_{\mathbb{H}}=0,\quad\langle w_{r},f_{H}(x) \rangle_{\mathbb{H}}=0.$ $$

To demonstrate the explainability of $r$, we consider the following statistical perspectives to verify its explainability (Hastie et al., ):

Linearity. Since $r$ is linear and operates on $f_{R}(x)$, its influence on the prediction is straightforward. Each component $[f_{R}(x)]_{i}$ contributes to $\hat{y}(x)$ proportionally to $[w_{r}]_{i}$, which is:

$$ $\frac{\partial\hat{y}(x)}{\partial[f_{R}(x)]_{i}}=[w_{r}]_{i}.$ $$

This allows us to quantify the effect of each residual feature on the prediction.

Orthogonality. The orthogonality of $\mathbb{H}_{C}$ and $\mathbb{H}_{R}$ ensures that the contributions from $h$ and $r$ are independent. Specifically, the covariance between $h(f_{H}(x))$ and $r(f_{R}(x))$ is zero:

$$ $\operatorname{Cov}\left(h(f_{H}(x)),\,r(f_{R}(x))\right)=\mathbb{E}_{x}\left[ \left(h(f_{H}(x))-\mu_{h}\right)\left(r(f_{R}(x))-\mu_{r}\right)\right]=0,$ $$

where $\mu_{h}=\mathbb{E}_{x}[h(f_{H}(x))]$ and $\mu_{r}=\mathbb{E}_{x}[r(f_{R}(x))]$, assuming $h(f_{H}(x))$ and $r(f_{R}(x))$ have finite second moments.

Variance Decomposition. The total variance of the prediction $\hat{y}(x)$ can be decomposed as

$$ $\operatorname{Var}\left(\hat{y}(x)\right)=\operatorname{Var}\left(h(f_{H}(x)) \right)+\operatorname{Var}\left(r(f_{R}(x))\right),$ $$

due to the independence of $h$ and $r$. This decomposition assists in understanding how much variability each component contributes.

Optimization Decoupling. Consider the expected loss over the data distribution $\mathcal{D}$:

$$ $\mathcal{L}(w_{h},w_{r},b_{h},b_{r})=\mathbb{E}_{(x,y)\sim\mathcal{D}}\left[ \ell\left(\hat{y}(x),y\right)\right],$ $$

where $\ell:\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}_{\geq 0}$ is a convex loss function, such as squared loss $\ell(\hat{y},y)=\frac{1}{2}(\hat{y}-y)^{2}$.

Due to the orthogonality, the gradients with respect to $w_{h}$ and $w_{r}$ are:

$$ $\nabla_{w_{h}}\mathcal{L}=\mathbb{E}_{(x,y)}\left[\ell^{\prime}\left(\hat{y}(x ),y\right)f_{H}(x)\right],$ $$

$$ $\nabla_{w_{r}}\mathcal{L}=\mathbb{E}_{(x,y)}\left[\ell^{\prime}\left(\hat{y}(x ),y\right)f_{R}(x)\right].$ $$

where $\ell^{\prime}$ denotes the derivative of $\ell$ with respect to its first argument. Since $f_{H}(x)$ and $f_{R}(x)$ are orthogonal and belong to different subspaces, the optimization of $w_{h}$ and $w_{r}$ can be treated separately, or decoupled. Specifically, we can minimize $\mathcal{L}$ by independently minimizing:

$$ $\mathcal{L}_{h}(w_{h},b_{h})=\mathbb{E}_{(x,y)}\left[\ell\left(h(f_{H}(x)),y-r (f_{R}(x))\right)\right],$ $$

$$ $\mathcal{L}_{r}(w_{r},b_{r})=\mathbb{E}_{(x,y)}\left[\ell\left(r(f_{R}(x)),y-h (f_{H}(x))\right)\right].$ $$

Preservation of Explainability. The explainability of $h$ is preserved because $r$ operates on features orthogonal to those used by $h$, and their contributions to $\hat{y}(x)$ are additive:

$$ $\hat{y}(x)=h(f_{H}(x))+r(f_{R}(x))=\underbrace{\langle w_{h},f_{H}(x)\rangle_{ \mathbb{H}}+b_{h}}_{\text{Explainable Component}}+\underbrace{\langle w_{r},f_ {R}(x)\rangle_{\mathbb{H}}+b_{r}}_{\text{Residual Component}}.$ $$

Each term can be examined independently. The explainable component aligns with the original explainable model, while the residual component adds adjustments based on the residual features.

Feature Contributions Interpretation. For any given feature $[f(x)]_{i}$, its contribution to the prediction can be determined based on whether it belongs to $\mathbb{H}_{C}$ or $\mathbb{H}_{R}$. Specifically:

$$ $\text{If }[f(x)]_{i}\in\mathbb{H}_{C}\implies\frac{\partial\hat{y}(x)}{ \partial[f(x)]_{i}}=[w_{h}]_{i},$ $$

$$ $\text{If }[f(x)]_{i}\in\mathbb{H}_{R}\implies\frac{\partial\hat{y}(x)}{ \partial[f(x)]_{i}}=[w_{r}]_{i}.$ $$

This allows for a granular interpretation of the impact of each feature on the prediction.

Mathematical Justification. From a functional analysis perspective, since both $h$ and $r$ are bounded linear operators on $\mathbb{H}_{C}$ and $\mathbb{H}_{R}$, respectively, and these subspaces are closed in $\mathbb{H}$, their sum $h+r$ is a bounded linear operator on $\mathbb{H}$. The boundedness ensures stability and well-posedness of the model. Moreover, the spectral properties of the operators can be analyzed separately. For example, if we consider the operator norms:

$$ $\|h\|_{\text{op}}=\sup_{f_{H}\in\mathbb{H}_{C},\|f_{H}\|_{\mathbb{H}}=1}|h(f_{ H})|,$ $$

$$ $\|r\|_{\text{op}}=\sup_{f_{R}\in\mathbb{H}_{R},\|f_{R}\|_{\mathbb{H}}=1}|r(f_{ R})|.$ $$

These norms provide insight into the maximum amplification of inputs by the respective operators.

The orthogonality and linearity between $f_{H}(x)$ and $f_{R}(x)$ ensure that the contributions from $h$ and $r$ are additive and independent, making the residual calibrator $r$ explainable when combined with $h$. Each feature’s impact on the prediction can be directly understood through the corresponding weights $w_{h}$ and $w_{r}$, enabling straightforward model interpretation. Specifically, since $f_{H}(x)$ and $f_{R}(x)$ are orthogonal, the cross terms $w_{h}^{\top}f_{R}(x)=0$ and $w_{r}^{\top}f_{H}(x)=0$ vanish, so $h$ and $r$ do not influence each other’s feature contributions.

Therefore, under the conditions of linearity and orthogonality, the residual calibrator $r$ preserves explainability when combined with $h$. The linearity ensures transparency in how input features affect the output, as each feature’s contribution is directly proportional to its weight. The orthogonality guarantees that the explainability of $h$ is not compromised by $r$ since they operate on separate feature spaces. In our implementation, these conditions are satisfied by utilizing logistic regression models for both $h$ and $r$, making our method provably explainable.

This completes the proof.
∎

###### Proof.

### A.5 Dataset Details

We use six mutagenicity datasets and one hepatotoxicity dataset. The mutagenicity datasets are: Mutag (Debnath et al., ), Mutagen (Morris et al., ), PTC-FM (Toivonen et al., ), PTC-FR (Toivonen et al., ), PTC-MM (Toivonen et al., ), PTC-MR (Toivonen et al., ), and the hepatotoxicity dataset is the Liver (Liu et al., ). Followed by Morris et al. (), we list the summary statistics of these datasets as

**Table 3: Summary statistics of seven datasets**
| Dataset | Mutag | Mutagen | PTC-FM | PTC-FR | PTC-MM | PTC-MR | Liver |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Samples | 188 | 4337 | 349 | 351 | 336 | 344 | 587 |
| Classes | 2 | 2 | 2 | 2 | 2 | 2 | 3 |
| Ground truth | 120 | 724 | 58 | 49 | 51 | 61 | 187 |

The ground truth indicates the true molecular substructures that impact molecular properties. As verified by Lin et al. (); Debnath et al. (), the ground truth substructures for six mutagenicity datasets consist of an aromatic group, such as a benzene ring, bonded with another functional group, such as methoxy, oxhydryl, nitro, or carboxyl groups (note that ground truth exists only for the mutagenic class). For the Liver dataset, the ground truth substructures annotated by chemists are: fused tricyclic saturated hydrocarbon moiety, hydrazines, arylacetic acid, sulfonamide moiety, aniline moiety, a class of proton pump inhibitor drugs, acyclic bivalent sulfur moiety, acyclic di-aryl ketone moiety, para oxygen and nitrogen di-substituted benzene ring, a relatively small number of com- pounds in the expanded LiverTox dataset, halogen atom bonded to a $sp^{3}$ carbon, and fused tricyclic structural moiety. A detailed illustration of Liver’s ground truth substructures are provided by Liu et al. ().

### A.6 Implementation Details

Our model is pre-trained on all data in the ZINC dataset (over $230$ million compounds) using ChemBERTa-2, with $15\%$ (default setting) of tokens in each input randomly masked. We extract all functional groups in the ZINC dataset as the vocabulary to expand the LLM’s tokenizer so that the fine-tuned LLM can better encode functional group-level inputs. We then fine-tune this model on Mutag, Mutagen, PTC-FM, PTC-FR, PTC-MM, PTC-MR, and Liver datasets. The fine-tuning is conducted on 1$\times$ NVIDIA RTX3090 GPU for about 3 hours. The detailed hyperparameters with their values are given in table 4. For experiments on model performance, we employ chain-of-thought prompting for the molecular property prediction tasks on LLMs.

**Table 4: Hyperparameters and their values we used for fine-tuning**
| Hyperparameter | Value |
| --- | --- |
| learning rate | 1e-5 |
| batch size | 128 |
| epochs | 30 |
| weight decay | 0.01 |
| gradient clipping | 1.0 |
| warmup proportion | 0.06 |
| max sequence length | 1024 |
| optimizer | AdamW |
| dropout rate | 0.1 |
| gradient accumulation steps | 1 |
| mixed precision training | True |

### A.7 Does the Residual Calibrator Improves Model Performance by training with more iterations?

We employ the training objective in A.11 to learn a residual calibrator that iteratively recaptures samples the linear model fails to predict accurately. We empirically study how training iterations influence the overall model predictions. As shown in fig. 4, we visualize the model performance on the Mutag, Mutagen, PTC-MR, and Liver datasets under different numbers of training iterations. As training iterations increase, model performance improves significantly until reaching a threshold. This suggests that more iterations on our designed loss lead to better performance. After the threshold, the model overfits the data, resulting in performance degradation. Therefore, increasing the number of training iterations helps improve model performance. Empirically, we found that 5 iterations yield optimal performance. A theoretical demonstration shows that training with multiple iterations increases model performance until a threshold, after which it declines, as follows.

Figure: (a) Mutag
Refer to caption: x4.png

Problem Setup. Given the objective the residual calibrator minimized during training:

$$ $\min_{h,r}\;\mathbb{E}_{(x,y)\sim\mathcal{S}_{\text{train}}}\left[\mathcal{L} \left(h(f_{H}(x))+r(f_{R}(x)),\;y\right)\right],$ (A.11) $$

where $\mathcal{S}_{\text{train}}$ is the empirical distribution of the training data and $\mathcal{L}:\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}_{\geq 0}$ is a convex, differentiable loss function, e.g., the squared loss $\mathcal{L}(\hat{y},y)=\frac{1}{2}(\hat{y}-y)^{2}$. We demonstrate that: initially, as the residual calibrator $r$ is trained, the model’s performance on unseen data improves, i.e., the generalization loss decreases. Beyond a certain threshold, further minimization of the training loss leads to overfitting, where the generalization loss starts to increase, and prediction accuracy on unseen data degrades.

###### Proof.

We aim to demonstrate that learning the residual calibrator $r$ with multiple training iterations initially improves the model accuracy, but after a certain training threshold, continued minimization of the training loss leads to overfitting, leading to the predictive accuracy on unseen data decline.

Let $\mathcal{X}$ and $\mathcal{Y}$ be the input and output spaces, respectively. Consider a feature extraction function $f:\mathcal{X}\rightarrow\mathbb{R}^{d}$ that maps inputs to a $d$-dimensional feature space. We assume that $f$ can be decomposed into two components:

$$ $f(x)=f_{H}(x)+f_{R}(x),$ $$

where $f_{H}(x)\in\mathbb{R}^{d_{c}}$ represents the explainable features used by the explainable model $h$, and $f_{R}(x)\in\mathbb{R}^{d_{r}}$ represents the residual features used by the residual calibrator $r$, with $d=d_{c}+d_{r}$. We assume that the feature components $f_{H}(x)$ and $f_{R}(x)$ are orthogonal, which means:

$$ $\langle f_{H}(x),f_{R}(x)\rangle=0\quad\text{for all }x\in\mathcal{X}.$ $$

The explainable model $h:\mathbb{R}^{d_{c}}\rightarrow\mathbb{R}$ is defined as a linear model:

$$ $h(f_{H}(x))=W_{h}^{\top}f_{H}(x)+b_{h},$ $$

where $W_{h}\in\mathbb{R}^{d_{c}}$ and $b_{h}\in\mathbb{R}$ are the weights and bias of $h$. The residual calibrator $r:\mathbb{R}^{d_{r}}\rightarrow\mathbb{R}$ is also defined as a linear model:

$$ $r(f_{R}(x))=W_{r}^{\top}f_{R}(x)+b_{r},$ $$

where $W_{r}\in\mathbb{R}^{d_{r}}$ and $b_{r}\in\mathbb{R}$ are the weights and bias of $r$. Due to the orthogonality of $f_{H}(x)$ and $f_{R}(x)$, the overall prediction model becomes:

$$ $\hat{y}(x)=h(f_{H}(x))+r(f_{R}(x))=W_{h}^{\top}f_{H}(x)+W_{r}^{\top}f_{R}(x)+b _{h}+b_{r}.$ $$

Our objective is to minimize the expected loss:

$$ $\mathcal{L}(W_{h},W_{r},b_{h},b_{r})=\mathbb{E}_{(x,y)\sim\mathcal{D}}\left[ \ell\left(\hat{y}(x),y\right)\right],$ $$

where $\ell\left(\hat{y}(x),y\right)$ is a convex and differentiable loss function, such as the squared loss $\ell(\hat{y},y)=\frac{1}{2}(\hat{y}-y)^{2}$, and $\mathcal{D}$ is the data distribution. We begin by considering the training loss over a finite training dataset $\{(x_{i},y_{i})\}_{i=1}^{n}$:

$$ $\mathcal{L}_{\text{train}}(W_{h},W_{r},b_{h},b_{r})=\frac{1}{n}\sum_{i=1}^{n} \ell\left(\hat{y}(x_{i}),y_{i}\right).$ $$

Initially, when $r$ is untrained or minimally trained, the model may be underfitting, and both the training loss $\mathcal{L}_{\text{train}}$ and generalization loss $\mathcal{L}_{\text{gen}}$ are high. By updating $W_{r}$ and $b_{r}$ via gradient descent to minimize $\mathcal{L}_{\text{train}}$, we have the updates:

$$ $W_{r}^{(t+1)}=W_{r}^{(t)}-\eta\nabla_{W_{r}}\mathcal{L}_{\text{train}}(W_{h},W _{r}^{(t)},b_{h},b_{r}^{(t)}),$ $$

$$ $b_{r}^{(t+1)}=b_{r}^{(t)}-\eta\nabla_{b_{r}}\mathcal{L}_{\text{train}}(W_{h},W _{r}^{(t)},b_{h},b_{r}^{(t)}),$ $$

where $\eta>0$ is the learning rate, and $t$ denotes the iteration number. Since $\ell$ is convex and differentiable, these updates ensure that the training loss decreases:

$$ $\mathcal{L}_{\text{train}}^{(t+1)}\leq\mathcal{L}_{\text{train}}^{(t)}.$ $$

During this phase, $r$ captures genuine patterns in the residual features $f_{R}(x)$ that are not explained by $h$. Consequently, the generalization loss decreases as well:

$$ $\mathcal{L}_{\text{gen}}^{(t+1)}\leq\mathcal{L}_{\text{gen}}^{(t)},$ $$

where

$$ $\mathcal{L}_{\text{gen}}(W_{h},W_{r},b_{h},b_{r})=\mathbb{E}_{(x,y)\sim \mathcal{D}}\left[\ell\left(\hat{y}(x),y\right)\right].$ $$

However, as training continues, $W_{r}$ and $b_{r}$ may begin to fit the noise or idiosyncrasies specific to the training data, especially if the model has a high capacity (i.e., $d_{r}$ is large relative to $n$). The fitting capacity of $r$ allows it to minimize $\mathcal{L}_{\text{train}}$ further, but this comes at the cost of increasing model complexity.

To formalize this, we consider the concept of Rademacher complexity $\mathfrak{R}_{n}(\mathcal{H})$ for the hypothesis class $\mathcal{H}$ associated with $r$. The Rademacher complexity provides a measure of the model’s ability to fit random noise in the data. The generalization error can be bounded as:

$$ $\mathcal{L}_{\text{gen}}(W_{h},W_{r},b_{h},b_{r})\leq\mathcal{L}_{\text{train} }(W_{h},W_{r},b_{h},b_{r})+2\mathfrak{R}_{n}(\mathcal{H})+\delta,$ $$

where $\delta$ is a constant dependent on the loss function and confidence level. As $\|W_{r}\|$ increases due to continued training, $\mathfrak{R}_{n}(\mathcal{H})$ increases, reflecting the higher complexity of $r$. This leads to circumstances that:

$$ $\mathcal{L}_{\text{train}}^{(t+1)}<\mathcal{L}_{\text{train}}^{(t)}\quad\text{ but}\quad\mathcal{L}_{\text{gen}}^{(t+1)}>\mathcal{L}_{\text{gen}}^{(t)}\quad \text{for }t\geq t^{*},$ $$

where $t^{*}$ is the iteration threshold beyond which overfitting occurs.

For linear models, the Rademacher complexity can be bounded by:

$$ $\mathfrak{R}_{n}(\mathcal{H})\leq\frac{B\|W_{r}\|}{\sqrt{n}},$ $$

where $B=\sup_{x\in\mathcal{X}}\|f_{R}(x)\|$. As $\|W_{r}\|$ increases, $\mathfrak{R}_{n}(\mathcal{H})$ increases, leading to a wider generalization gap. This increase in model complexity without a corresponding increase in true predictive power causes the model to generalize poorly on unseen data, despite the training loss decreasing. This phenomenon is a bias-variance trade-off: the variance increases significantly due to overfitting, outweighing any small reductions in bias achieved by further minimizing the training loss.

In conclusion, while initial training of the residual calibrator $r$ improves model accuracy by reducing both the training loss and the generalization loss, continued training beyond a certain threshold leads to overfitting. The residual calibrator begins to model noise in the training data, increasing its complexity and causing the generalization loss to increase. This results in a decline in prediction accuracy on unseen data, suggesting the importance of strategies such as early stopping or regularization to prevent overfitting.

This completes the proof.
∎

###### Proof.

### A.8 How to choose the optimal number of principal components?

We exclude task-irrelevant embeddings using a dimensionality reduction method, namely EFPCA. To empirically determine the optimal number of principal components for our implementation, we compare model performance metrics (classification accuracy and explanation accuracy) across four datasets under different numbers of principal components. As shown in fig. 5, both metrics tend to converge as the number of principal components exceeds 20. This indicates that when the number of components surpasses 20, the contribution of additional components to molecular property prediction becomes trivial. In this scenario, adding more components produces diminishing marginal benefits while significantly increasing model complexity, which in turn reduces explainability. Therefore, we choose the top 20 principal components to explain the variance in molecular properties, seeking for a balance between performance and explainability.

Figure: (a) Mutag
Refer to caption: x8.png

### A.9 Does EFPCA effectively works?

In addition to the analysis in section A.8, we demonstrate that the dimensionality reduction by EFPCA effectively preserves the most explanatory components. We compare the model performance (classification accuracy and explanation accuracy) across seven datasets with and without dimensionality reduction. As shown in table 5, when using only 20 principal components, the model performance improves by no more than $5\%$ compared to using all 384 components (i.e., no dimensionality reduction). This indicates that EFPCA effectively preserves the most task-relevant and important information in LLM embeddings while excluding noisy components. These preserved components achieve comparable performance to the models with all components while being significantly simpler and more explainable. This showcases the success of our dimensionality reduction in maintaining model performance while enhancing explainability.

**Table 5: Model performance without EFPCA over seven datasets**
| Dataset | Classification Accuracy ($\%$) | Explanation Accuracy ($\%$) |
| --- | --- | --- |
| Mutag | 94.9±1.6 | 96.1±3.0 |
| Mutagen | 86.4±1.4 | 91.2±1.6 |
| PTC-FR | 78.7±1.2 | 82.7±1.7 |
| PTC-FM | 68.1±1.5 | 81.1±2.0 |
| PTC-MR | 70.5±1.7 | 76.5±2.6 |
| PTC-MM | 80.9±2.7 | 75.3±2.2 |
| Liver | 57.3±1.6 | 83.8±1.9 |

### A.10 Does The Choice of n 𝑛 n italic_n in N-gram makes a difference?

We compare the different values of $n$ in n-gram via cross-validation based on our two evaluation metrics, classification accuracy and explanation accuracy. The results in fig. 6 suggest an overall trend that as $n$ goes from 1 to 3, both classification accuracy and explanation accuracy improve; as $n$ goes from 4 to 9, both classification accuracy and explanation accuracy drop. On the four datasets we used for experiments, three of them show that good model performance can be achieved when $n$ is taken to be 3. As $n$ grows from small to large, it encourages the model to capture more contextual semantics, including interactions between functional groups, which allows for a significant improvement in prediction. When $n$ exceeds a certain threshold, irrelevant or even toxic information emerges from the captured contextual information (i.e., irrelevant long-range dependencies), making the overall model utility gradually decreases.

Figure: (a) Mutag
Refer to caption: x12.png

### A.11 More explanation visualizations

We randomly select one sample from each of the six remaining datasets and provide explanation visualizations based on MoleX. Specifically, fig. 7, fig. 8, fig. 9, fig. 10, fig. 11, and fig. 12 display the samples selected from the Mutagen, PTC-FM, PTC-MM, PTC-FR, PTC-MR, and Liver datasets, respectively. On the left, we visualize the identified molecular substructures based on different methods, with the ground truth indicating the true substructures influencing molecular properties as given by domain experts. Red marks on the molecular graph highlight important components impacting molecular properties revealed by different methods. We compare with three baselines: OrphicX (Lin et al., ), GNNExplainer (Ying et al., ), and PGExplainer (Luo et al., ), as well as MoleX with and without residual calibration (w/ denotes with and w/o denotes without). On the right, we showcase the n-gram coefficients-based contribution scores of important functional groups revealed by MoleX, ranging from 0 0 to $100$, where higher scores suggest a more significant contribution to molecular properties.

Taking fig. 7 as an example, MoleX precisely identifies the ground truth substructures for the sample from the Mutagen dataset. Specifically, MoleX highlights the benzene ring bonded with an amino group on the upper left as vital substructures to explain the molecule’s mutagenicity. The contribution scores computed by MoleX indicate that the benzene ring has the highest contribution to molecular properties, followed by the amino group. This aligns with the ground truth that a benzene ring bonded with an amino group leads to mutagenicity (Lin et al., ; Debnath et al., ). Therefore, MoleX accurately captures the important functional groups (i.e., the benzene ring and the amino group) and the interaction between them, revealing their precise bonding. As the ground truth indicates, only the bonded benzene and amino group together impact the molecular properties. In contrast, other methods provide only atom or bond-level explanations and fail to discover important functional groups as a whole. They identify only a few atoms and bonds in the benzene or amino group and fail to capture the interaction between these two functional groups. Consequently, these atom or bond-level explanations are insufficiently faithful in explaining molecular properties, as individual atoms or bonds have limited impact on overall molecular properties (Mirghaffari et al., ). The explanation visualizations for samples from other datasets also demonstrate MoleX’s effectiveness in identifying important substructures and their interactions, aligning with chemical concepts to explain molecular property predictions.

Figure: Figure 7: Explanation visualization of a molecule from the Mutagen dataset (left), and contribution scores of the identified functional groups offered by MoleX (right).
Refer to caption: x16.png

Figure: Figure 8: Explanation visualization of a molecule from the PTC-FM dataset (left), and contribution scores of the identified functional groups offered by MoleX (right).
Refer to caption: x17.png

Figure: Figure 9: Explanation visualization of a molecule from the PTC-MM dataset (left), and contribution scores of the identified functional groups offered by MoleX (right).
Refer to caption: x18.png

Figure: Figure 10: Explanation visualization of a molecule from the PTC-FR dataset (left), and contribution scores of the identified functional groups offered by MoleX (right).
Refer to caption: x19.png

Figure: Figure 11: Explanation visualization of a molecule from the PTC-MR dataset (left), and contribution scores of the identified functional groups offered by MoleX (right).
Refer to caption: x20.png

Figure: Figure 12: Explanation visualization of a molecule from the Liver dataset (left), and contribution scores of the identified functional groups offered by MoleX (right).
Refer to caption: x21.png

### A.12 Can other statistical learning models be augmented with the LLM knowledge?

In addition to the linear model, we augment various statistical learning models with the LLM knowledge and test them on seven datasets. The classification accuracy and explanation accuracy are shown in table 6 and table 7, respectively. Other linear models, such as ridge regression, LASSO, and linear discriminant analysis, achieve comparable performance to MoleX and showcase the generalizability of LLM knowledge augmentation on linear models. Additionally, the polynomial regression, as a more complicated linear model, achieves better performance compared to the simpler ones shown above. For more complex models, such as tree-based and ensemble learning models, the performance is even better, achieving incredible results across all seven datasets. These empirical studies suggest that augmenting statistical machine learning models with LLM knowledge significantly improves performance. Moreover, compared to simple models, the models exhibit more powerful data fitting capabilities become more predictive after the LLM augmentation. However, model complexity generally trades off with explainability. Considering this, we select the logistic regression as our base model due to its optimal balance between explainability and performance.

**Table 6: Classification Accuracy across different machine learning models over seven datasets (%)**
| Method | Mutag | Mutagen | PTC-FR | PTC-FM | PTC-MR | PTC-MM | Liver |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Ridge Regression | 90.7±1.2 | 84.1±1.3 | 72.4±2.0 | 65.2±2.0 | 69.8±1.4 | 77.5±1.5 | 58.1±1.6 |
| LASSO | 91.9±1.7 | 84.4±0.7 | 75.1±2.1 | 65.8±1.7 | 65.2±0.9 | 74.2±1.2 | 58.7±1.8 |
| Linear Discriminant Analysis | 89.9±1.9 | 83.6±1.2 | 75.2±1.9 | 65.7±1.9 | 69.3±1.8 | 76.8±2.0 | 57.7±1.3 |
| Polynomial Regression | 93.9±2.4 | 87.2±2.0 | 77.1±2.1 | 67.3±1.8 | 70.2±2.3 | 79.5±1.8 | 60.2±2.4 |
| Support Vector Machine | 93.9±1.6 | 86.6±1.5 | 73.4±1.9 | 69.3±2.6 | 69.5±2.0 | 78.6±1.3 | 61.5±2.9 |
| Decision Tree | 89.7±2.1 | 79.5±1.2 | 72.4±1.8 | 64.3±2.1 | 68.5±1.5 | 74.4±1.4 | 59.5±2.2 |
| Random Forest | 92.8±2.7 | 84.4±1.7 | 77.3±2.1 | 68.6±2.5 | 71.0±2.2 | 77.2±2.1 | 62.7±2.7 |
| Gradient Boosting Machine | 94.8±2.1 | 85.3±1.9 | 78.9±1.9 | 69.4±2.8 | 72.2±2.1 | 79.2±1.9 | 63.9±2.6 |
| XGBoost | 94.6±2.3 | 85.0±2.0 | 78.7±2.2 | 70.1±2.3 | 73.4±2.9 | 78.1±2.1 | 63.0±2.3 |
| MoleX (Ours) | 91.6±2.0 | 83.7±0.9 | 74.4±1.9 | 64.2±1.4 | 68.4±2.3 | 76.4±1.8 | 54.9±2.4 |

**Table 7: Explanation Accuracy across different machine learning models over seven datasets (%)**
| Method | Mutag | Mutagen | PTC-FR | PTC-FM | PTC-MR | PTC-MM | Liver |
| --- | --- | --- | --- | --- | --- | --- | --- |
| Ridge Regression | 92.8±1.1 | 89.5±1.3 | 79.0±1.2 | 78.1±1.6 | 72.5±2.5 | 69.7±2.3 | 82.4±1.7 |
| LASSO | 92.3±1.5 | 89.6±0.9 | 76.9±1.8 | 81.2±1.9 | 70.4±2.3 | 70.7±2.1 | 81.3±1.8 |
| Linear Discriminant Analysis | 92.9±1.8 | 88.5±1.9 | 80.7±2.3 | 80.1±2.2 | 71.7±2.8 | 71.3±1.6 | 87.8±1.6 |
| Polynomial Regression | 94.3±2.1 | 91.9±1.6 | 80.1±1.9 | 82.9±1.9 | 79.3±2.3 | 75.4±1.7 | 81.0±2.2 |
| Support Vector Machine | 92.0±1.7 | 92.0±1.6 | 84.7±2.2 | 86.3±2.0 | 80.1±2.3 | 76.0±2.3 | 81.9±2.1 |
| Decision Tree | 87.6±1.9 | 89.1±1.5 | 78.6±2.0 | 80.7±1.6 | 73.1±2.1 | 74.2±1.8 | 76.0±1.8 |
| Random Forest | 93.2±1.9 | 90.5±1.8 | 82.1±2.1 | 84.2±2.2 | 74.2±2.0 | 74.5±2.1 | 81.2±2.0 |
| Gradient Boosting Machine | 92.7±2.2 | 92.4±1.5 | 82.9±2.3 | 85.2±2.4 | 73.9±2.9 | 77.7±2.6 | 84.5±2.4 |
| XGBoost | 95.6±1.8 | 90.7±1.7 | 84.0±2.2 | 82.0±2.3 | 74.4±2.7 | 77.4±2.2 | 86.2±2.5 |
| MoleX (Ours) | 92.6±1.7 | 89.0±0.9 | 79.3±2.6 | 77.9±2.6 | 73.4±2.8 | 72.3±3.0 | 80.3±2.5 |

### A.13 Classification analysis via Confusion Matrix

As shown in fig. 13, we visualize the classification result via confusion matrix at a random round on the Mutag and PTC-MR datasets. For Mutag, we achieve high precision in predicting the positive class due to fewer false positives and high recall for the positive class, reflecting the model’s effectiveness in identifying positive instances. Furthermore, the model shows a good balance between precision and recall, with a low number of false positives and false negatives. For PTC-MR, the model achieves lower precision compared to the Mutag due to a higher number of false positives. The confusion matrix also suggests that the model struggles with false negatives and false positives, indicating areas for improvement. This analysis highlights the strengths and weaknesses of the model, providing insight for further model refinement.

Figure: (a) Mutag
Refer to caption: x22.png

### A.14 An illustration of Group SELFIES

As illustrated in fig. 14, the 4-Nitroanisole (\ceC7H7NO3) can be represented by Group SELFIES with three functional groups separated by square brackets: a benzene ring, a nitro group, and a methoxy group (different functional groups are displayed in different colors).

Figure: Figure 14: Molecular representation of 4-Nitroanisole (\ceC7H7NO3)
Refer to caption: x24.png

### A.15 Broader Impact

This study on explainable molecular property prediction using the LLM-augmented linear model has significant real-world application. The efficiency of linear models allows for fast inference of large-scale molecular data, potentially accelerating drug discovery and materials design processes. This efficiency, combined with the predictive accuracy enhanced by LLM-derived features, and lead to more cost-effective and computationally efficient development in the pharmaceutical and materials industries. The high explanation accuracy of our approach addresses an imperative need in scientific contexts. It provides explainable insights into structure-property relationships, offering faithfulness in predictions and facilitating their adoption in high-stake fields like healthcare. This explainability supports scientific insights, potentially leading to new discoveries in molecular science. Moreover, our method’s balance of predictive accuracy, explainability, and efficiency could serve as a template for developing explainable AI models in other domains, which contributes to the broader goal of trustworthy AI. This could have far-reaching impacts on fields ranging from personalized medicine to sustainable chemistry where understanding the decision-making mechanisms behind predictions is crucial. However, we must also consider potential risks, such as over-reliance on model predictions without appropriate expert validation or the potential misuse in designing harmful molecules. Responsible implementation and research into model limitations are essential.

### A.16 Limitations and future works

The proposed explainable molecular property prediction method has some limitations and needs further studies.

- •
Generalizability: Enhancing the generalizability of explainable models to deal with different molecular datasets across various chemical domains while preserving explainability to structure-property relationships remains a persistent challenge.
- •
Impact of LLM quality: Though our empirical studies discuss the model performance of Llama3.1 and GPT-4o on molecular property prediction, LLM quality is still a topic that deserves to be explored in-depth. Future studies may discuss how LLM quality impacts the augmented linear model, e.g., model performance change using weak LLMs or LLMs without fine-tuning.
- •
Trade-off between complexity and performance: In pursuit of explainability, we employ a linear model, which inherently risks underfitting when faced with complex data patterns. Our preliminary experiments comparing MoleX with more sophisticated statistical learning models show marginally better performance from these complex models. Future research could explore the trade-off between model complexity and performance in the context of LLM knowledge augmentation and investigate optimal balances between explainability and performance.

## Peer Reviews

confidence: 2
contribution: 3
initial_rating: 8
presentation: 3
questions: Can you discuss the potential limitations of the linearity assumption in MoleX? How might non-linear relationships between features impact the accuracy of the framework?

How do you interpret the n-gram features in the context of molecular properties? Can you provide examples of how specific n-grams might be related to particular molecular properties or mechanisms?

How robust is MoleX to different types and levels of noise? Are there any mechanisms in place to mitigate the effects of noise or outliers on the computed feature contributions?

You mention that "Theoretically, we provide a mathematical foundation to justify MoleX's explainability". Could you clarify what you mean by "theoretically" in this context? Are you implying that the mathematical foundation is theoretical in nature, or is this a statement about the theoretical soundness of the foundation?
soundness: 3
strengths: This paper exhibits several notable strengths that contribute to its overall quality and impact.
- Well-Written and Clear Structure
  The author does an good job of discussing the main points in the article, providing a clear and concise overview of the framework, its theoretical background, and its advantages. The text is well-organized, making it easy to follow the author's reasoning and understand the key takeaways.
- Relevant Problem
  The problem addressed in the paper is relevant and timely, tackling a significant challenge in the field. The author provides a clear motivation for the research, highlighting the importance of understanding feature contributions in molecular properties.
- Advantages Over State-of-the-Art
  The proposed solution, MoleX, offers several advantages over existing state-of-the-art methods. By employing n-gram coefficients in a linear model, MoleX provides a more nuanced understanding of feature contributions, enabling researchers to better interpret the results. The author effectively highlights the benefits of MoleX, demonstrating its potential to improve the field.
summary: The article discusses a framework called MoleX, which aims to provide a granular interpretation of the impact of each feature on predictions. The framework uses a linear model to measure the contributions of decoupled n-gram features to molecular properties. It discusses the explainable component and residual component of a model and how they can be used to determine the contribution of each feature to the prediction.
weaknesses: - Limited Experimental Evaluation
  The paper's experimental evaluation is limited, with a narrow focus on a specific dataset and task. To further demonstrate the effectiveness of MoleX, it would be beneficial to conduct more extensive experiments on diverse datasets and tasks.
- Lack of Comparison to Other Methods
  The paper could benefit from a more comprehensive comparison to other state-of-the-art methods. While the author mentions some existing approaches, a more detailed analysis of their strengths and weaknesses would provide a clearer understanding of MoleX's advantages.
- Future Work
  The paper concludes somewhat abruptly, without discussing potential avenues for future research. Outlining possible extensions or applications of MoleX would provide a clearer direction for future studies and encourage further exploration of the framework.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: 1. How is “faithfulness” in MoleX’s explanations quantified or evaluated? Could you provide more details on the metrics or assessment criteria used to determine that the model’s explanations are accurate and reliable from a chemical standpoint?

2. How does MoleX differ from or improve upon prototype-based approaches, which are also commonly used to enhance interpretability by referencing representative examples?

3. Could you provide more details about the residual calibration strategy? Specifically, how is it implemented, what types of prediction errors does it address, and how does it affect the interpretability of the final predictions?

4. Since MoleX relies on a linear model, are there cases where it struggles to represent the complexity of molecular data, even with LLM embeddings and residual calibration? How do you manage this limitation?

5. Which metrics were used to evaluate the explainability of MoleX compared to other interpretable models?
soundness: 3
strengths: The paper stands out for its innovative approach to bridging the gap between prediction accuracy and interpretability in molecular property prediction. While explainability and accuracy have traditionally been seen as trade-offs, MoleX proposes a hybrid framework that leverages the strengths of both linear models and LLMs. This combination—using LLM embeddings to build a more interpretable linear model, while employing information bottleneck fine-tuning and sparsity to adapt these embeddings—is a creative solution to a long-standing problem. Furthermore, the introduction of a residual calibration strategy to enhance a linear model’s capacity in representing complex LLM-derived features demonstrates originality by addressing a specific limitation of linear methods. Overall, this paper is a meaningful contribution that combines originality with technical rigor, clarity, and significant practical implications.
summary: The paper introduces MoleX, a framework designed to achieve accurate and interpretable molecular property predictions by leveraging the knowledge embedded in large language model (LLM) embeddings while maintaining a simple, interpretable linear model structure. MoleX combines the predictive power of LLM embeddings with the transparency of linear models, aiming to bridge the gap between complex models that excel in prediction and simpler models that offer more interpretability.
weaknesses: While MoleX offers a creative use of LLM embeddings for interpretable molecular property prediction, the novelty might be limited as the work relies heavily on existing LLM representations and linear models. The use of embeddings from pre-trained LLMs combined with linear regression isn’t entirely new, as similar approaches have been applied in other fields for simplifying complex models while preserving interpretability. The authors could enhance the novelty by providing a clearer comparison of MoleX with similar hybrid approaches in other domains, or by further distinguishing their method’s unique aspects in the context of molecular property prediction. The authors can compare their method with prototype-based approaches as ProtoPNet [1] and xDNN [2].

[1] Chen, Chaofan, et al. "This looks like that: deep learning for interpretable image recognition." Advances in neural information processing systems 32 (2019).

[2] Angelov, Plamen, and Eduardo Soares. "Towards explainable deep neural networks (xDNN)." Neural Networks 130 (2020): 185-194.

---

confidence: 5
contribution: 2
initial_rating: 3
presentation: 1
questions: 1. **Line 261**: How are the $\ell_0$ penalty, $a_{kj}$, and $\phi_j(t)$ derived? These symbols should be explained in the main text rather than in the appendix to ensure readers understand their significance without having to refer to supplementary materials.

2. **Input Consistency for Explainable Model and Residual Calibrator**: Are the inputs for the explainable model and residual calibrator identical, specifically the reduced LLM embeddings after EFPCA? Clarifying this point would help readers better understand the workflow and the rationale behind the model architecture.

3. **Line 308**: How is $f$ defined, and are $f_H$ and $f_R$ the outputs of EFPCA? If so, why are they different, and how does EFPCA produce distinct outputs from the same input? Detailed clarification on this process would enhance the comprehension of the model’s internal mechanisms.

4. **Inconsistent Notation**: The notation varies across different sections and should be made consistent to avoid confusion.
soundness: 3
strengths: 1. **Reduced Complexity via Linear Model**: The use of a linear model effectively lowers overall computational complexity, making the framework suitable for large-scale applications.

2. **Efficient Dimensionality Reduction**: The dimensionality reduction technique minimizes computational demands, optimizing the handling of high-dimensional LLM embeddings.

3. **High Interpretability**: The model offers strong interpretability, providing reliable, chemically meaningful explanations for molecular predictions.
summary: This paper introduces MoleX, a framework for explainable molecular property prediction that combines large language model (LLM) embeddings with a linear model to achieve both accuracy and interpretability. MoleX leverages fine-tuned LLM knowledge and residual calibration to enhance prediction accuracy while enabling CPU-efficient inference, showing up to 12.7% improved performance and significant parameter reduction compared to standard LLMs.
weaknesses: This paper, though innovative, exhibits certain weaknesses:

1. **Lack of Comparison with Relevant Explainability Studies**: The authors did not compare their approach with other information bottleneck based methods, for example, SGIB and VGIB [1,2], which also employ IB techniques to identify explainable subgraphs most relevant to molecular properties. Including these comparisons would provide a more comprehensive evaluation of the proposed method.

2. **Poor Presentation**: The presentation of the paper is lacking in clarity, with several areas insufficiently explained. This detracts from the paper’s readability and hinders a full understanding of the methodology. This is my primary reason for recommending rejection, as the paper requires major revision to improve clarity. Specific issues are detailed in the "Questions" section.

3. **Omission of LLM-Based Method Comparisons in Table 2**: In Table 2, the authors do not compare LLM-based methods evaluated in Table 1. The reasons for this omission should be explicitly addressed within the paper to clarify the experimental choices and maintain consistency in the evaluation.

[1] Yu, Junchi, et al. "Graph Information Bottleneck for Subgraph Recognition." International Conference on Learning Representations.

[2] Yu, Junchi, Jie Cao, and Ran He. "Improving subgraph recognition with variational graph information bottleneck." Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. 2022.

---

confidence: 4
contribution: 1
initial_rating: 3
presentation: 1
questions: - I think the authors would want to reformulate the problem and reorganize the paper to clearly and precisely state their contributions and claims, and describe their method framework, better following academic integrity.
soundness: 2
strengths: - Experimental results across seven datasets demonstrate that MoleX outperforms existing baselines (both GNN and LLM models) in predictive and explanation accuracy.
- The evaluation metrics, including classification accuracy, explanation AUC, and computational efficiency, present a holistic view of the model's advantages. The ablation studies reinforce the robustness of the method by exploring different parameter choices and model components.
summary: The paper introduces MoleX, a framework combining linear models with large language model (LLM) embeddings for explainable and accurate molecular property prediction. The framework's unique features include:
- Information bottleneck-inspired fine-tuning and sparsity-inducing dimensionality reduction to enhance the informativeness of LLM embeddings.
- Residual calibration to address prediction errors, thereby boosting the performance of the linear model without sacrificing explainability.
- Quantification of functional group contributions to molecular properties through n-gram coefficients.
weaknesses: - The methodology section, while comprehensive, largely restates existing techniques without significant novel contributions. Many of the core ideas, such as the information bottleneck fine-tuning and dimensionality reduction via EFPCA, are established approaches adapted with minimal modification. This raises concerns regarding the originality of the proposed framework.
- The claims surrounding the methodological innovations appear overstated, given that the fundamental components, such as fine-tuning with variational information bottleneck and principal component analysis, are widely known and previously employed in similar contexts. The application of these techniques, while effective, may not sufficiently differentiate this work as a significant leap forward in the field.

---

confidence: 2
contribution: 2
initial_rating: 5
presentation: 1
questions: - Please add an algorithm describing the whole procedure of MoleX including both its training and inference.
- I still don't understand why the explanation accuracy of the proposed method improves when with calibration. Could it be because the dimensionality reduction was too aggressive that a part of the features necessary for explanation was removed? I would appreciate it if the authors could point out any misunderstanding of mine if any / provide explanation on it.
soundness: 2
strengths: - The problem setting to seek for a predictive model for molecules with explainability is a high-demand setting for practitioners, and I would appreciate the authors for achieving the state-of-the-art performance on this problem setting.
- The experiments are exhaustive in terms of the methods compared and experiments conducted, although many of the details are not in the main part of the paper.
summary: The authors are interested in learning a molecular property prediction model using a large language model (LLM) with explainability. The proposed framework called MoleX consists of four components.

First, it extracts a useful feature from LLM embedding using the variational information bottleneck (VIB) principle (Section 4.1).
Second, it applies dimensionality reduction to the feature extracted in the first step based on the explainable functional principal component analysis, which roughly corresponds to PCA with sparsity-inducing penalty to each component (Section 4.2).
Third, using the explainable feature obtained in the second component, it fits a linear model that maps the explainable feature to the target variable (Section 4.2).
Fourth, with the predictive model fixed, it fits another linear model using the features not used in the third component for residual calibration (Section 4.2).

The authors evaluate the effectiveness of MoleX in terms of its predictive power as well as its explainability and computational efficiency.
The experimental results suggest MoleX is better than baseline methods in every aspect.
Sensitivity analyses on hyperparameters are conducted in Section 5.3, but most of the details are sent to the appendix.
weaknesses: - The overview of MoleX is not clear, especially, how we train it and use it for obtaining prediction with explanations in a step-by-step manner. Although Figure 1 suggests a vague overview of it, as far as I am aware of, there is no rigorous explanation of it. Since it is very essential to describe MoleX, I would like to request the authors to explain it in an algorithmic way in the main part of the paper (not in the appendix).
- There are many mathematically incorrect/unclear descriptions, which pose challenges in understanding the paper.  A few examples are as follows:
  - Line 152: Since $g$ and $y$ are not sets, $f\colon g\to y$ should be $f\colon g\mapsto y$.
  - Definition 4.1: There are multiple concepts that are not explained, namely, an operator $\mathcal{D}^2$ and an inner product $\langle \cdot, \cdot \rangle_{\gamma}$. $k$ appears both in the objective function and constraints, but there is no explanation about it.
  - Line 260: $\ell_0$ penalty does not directly appear in Def. 4.1 (in a finite dimensional case, the regularization term boils down to the $\ell_0$ penalty, but the definition is given for the infinite dimensional case).
  - Lines 309-311: Given that $f_H(x)\in\mathbb{R}^{d_c}$ and $f_R(x)\in\mathbb{R}^{d_r}$, $\langle f_H(x), f_R(x) \rangle$ cannot be defined in a natural way.

Meta-review: **Summary:** The authors use LLMs to extract the maximum amount of task-relevant knowledge from their dimensionality-reduced embeddings and fit a linear model for explainable inference. This is used to determine the contribution of each feature for prediction. They also introduce residual calibration to address prediction errors stemming from linear models' insufficient expressiveness of complex LLM embeddings. The faster performance of MoleX is demonstrated experimentally.  

**Strengths:** The paper tackles an interesting problem that holds practical importance. It also provides an accurate yet highly interpretable approach. The main advantage seems to be using an LLM for feature extraction, which captures the higher-order interaction of the components. The method appears to outperform other techniques experimentally.  

**Weaknesses:** The main shortcomings of the paper are its novelty and experimental evaluation. Specifically, the paper combines well-established techniques to present a new framework. However, it does not add much as a technical contribution. In terms of experiments, several reviewers point out the lack of comparison to previously known interpretable techniques. There are also some issues with the writing and notation.  

**Decision:** Although the paper provides some interesting ideas, I believe the novelty is still the main shortcoming of the work. The authors should also try to justify the rationale for each design choice and perform separate ablations. Moreover, additional comparisons to previous techniques would also improve the work. Finally, the authors should verify that the notation and the presentation are clear and consistent. At this stage, I recommend rejection.

Decision: Reject