Thierry Arrabal

Marija Stojanova

Isabelle Guérin-Lassous

Joris Picot

Isabelle Guerin-Lassous

Experiment-driven platform for link quality estimation in IEEE 802.11 WLANs

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

The Wi-Fi technology, based on the IEEE 802.11 standard, is constantly undergoing many changes (including Wi-Fi 4, Wi-Fi 5, and Wi-Fi 6), and is now used in a variety of contexts ranging from local networks (home, universities, shopping malls, etc.) to fleets of UAVs (Unmanned Aerial Vehicles), vehicular networks, or IoT (Internet of Things). Wi-Fi performance is impacted by many factors including its environment and potential interferences (from Wi-Fi or other radio networks), or the distance and antenna orientation between the communicating devices. The performance can vary greatly, both on short and long time scales. Some applications rely on low-level performance metrics in order to properly operate, a.k.a., the performance obtained at the physical and data link layers. In this study, we focus on the Frame Delivery Ratio (FDR) of a Wi-Fi link consisting of one source node and one destination node. The FDR is also an indicator of throughput and packet delay, as frame losses imply frame retransmissions and thus delays in the segment transmission times.

Monitoring the FDR is not an easy task. As far as we know, there are currently no simple tools available for the FDR monitoring during application run-time and it requires either to modify the Wi-Fi driver to collect statistics from the Wi-Fi interfaces or to deploy at least one extra Wi-Fi interface configured in monitor mode, in addition to the Wi-Fi interfaces of the source and destination nodes. In monitor mode, a Wi-Fi interface can collect all the frames it detects, regardless of their destination, from which the FDR can be computed. Both methods have the disadvantage of requiring the modification of the system (at the hardware or software level).

Overcoming these difficulties calls for developing an FDR estimation tool that requires only on-line measurements available at the source and/or destination nodes. But, very few such tools exist for Wi-Fi devices and they are based on datasets collected using old Wi-Fi technologies not using all the features implemented since Wi-Fi 4 which can have an impact on the FDR. In addition, the authors are not aware of the existence of any publicly-available datasets that use recent Wi-Fi versions and provide the full physical and link layer headers. The lack of such datasets makes it difficult to design the needed estimation tools. To bridge this gap, the main contributions presented in this paper can be summarized as follows:

• Deployment of an experimental Wi-Fi platform, in various environments, for the collection of physical and data link layer headers as experimental data, and providing an open-source FDR computation method using that data (see Sec. III). • Development of an open-source framework for FDR analysis and estimation using Machine Learning (ML) methods in Python (see Sec. IV). • Publication of the obtained experimental datasets for any future uses [START_REF] Arrabal | Wi-Fi 5 packets capture[END_REF]. We provide a state of the art on the problem studied in this paper in Sec. II and we conclude in Sec. V.

II. RELATED WORK

Network monitoring tools: Linux commands like ip and iw, provide networking statistics (such as the number of transmitted and received packets) on a network interface since the interface was set up. It is difficult to use these commands to obtain statistics on a very short time scale. Wavemon is a Linux tool that provides lowlevel Wi-Fi parameters, e.g., the MCS index (Modulation and Coding Scheme). However, these statistics are shown in an ncurses window, making it difficult to extract raw data. Monitoring tools, e.g., Wireshark, can capture all the frames received and transmitted by a network interface and provide the data included in all the headers of the different TCP/IP layers. Unfortunately, to get data from the physical and link layers of Wi-Fi frames, the Wi-Fi interface has to been configured in monitor mode preventing it from communicating in the monitored network. This implies using a new interface solely for monitoring purposes, and not for communication. Finally, some Wi-Fi drivers may be modified in order to get statistics on the FDR, however this approach requires a modification of the driver, which is not always straightforward or desirable.

FDR estimation tools:

The survey of Gregor et al. [START_REF] Cerar | Machine learning for wireless link quality estimation: A survey[END_REF] lists the different Machine Learning (ML) solutions that estimate a wireless link's quality. The survey includes only a few works relating to the Wi-Fi technology [3]- [START_REF] Cerar | Learning to fairly classify the quality of wireless links[END_REF] and, among them, only two estimate the FDR [START_REF] Srinivasan | PRR is not enough[END_REF], [START_REF] Cerar | Learning to fairly classify the quality of wireless links[END_REF]. Moreover, these works are based on old Wi-Fi versions (the IEEE 802.11a/b/g standard amendments) which are significantly different from Wi-Fi 4/5/6. Scanzio et al. [START_REF] Scanzio | Predicting Wi-Fi link quality through artificial neural networks[END_REF] study a Wi-Fi 4-based network, however their FDR estimation method requires active-probing packets and the addition of a specialized Software-Defined MAC for managing the Wi-Fi interface.

Publicly available datasets and testbeds: Gregor et al. [START_REF] Cerar | Machine learning for wireless link quality estimation: A survey[END_REF] provide a list of open-source datasets that can be used to estimate link quality metrics. Several datasets use Wi-Fi networks [START_REF] Aguayo | Link-Level Measurements from an 802.11b Mesh Network[END_REF]- [START_REF] Anderson | The impact of directional antenna models on simulation accuracy[END_REF], however none of them use recent Wi-Fi devices. We found datasets using recent Wi-Fi devices, but they do not provide the entire content of the physical and link headers of the collected frames 123 [START_REF] Barrachina-Muñoz | Wi-Fi All-Channel Analyzer[END_REF]. Concerning the publicly-available testbeds, some of them implement the Wi-Fi technology, e.g., NICTA, w-iLabt, NITOS, NTUA, federated within the Fed4FIRE, FIT wireless, and ORBIT testbeds. Among them, only OR-BIT implements Wi-Fi versions 4 or higher. Moreover, these testbeds offer little control in terms of node location and antenna orientation, and only a limited diversity of environments.

1 https://zenodo.org/record/5616432#.Y34opTOZOV4 2 https://zenodo.org/record/814898#.Y34saTOZOV4 3 https://zenodo.org/record/3952557#.Y3ID-76ZOV4

III. EXPERIMENTAL DATA COLLECTION

Our experiments are focused on estimating the quality of a Wi-Fi link between a source and a destination. The source emits a stream of data which the destination receives over an imperfect radio medium whose quality we wish to assess.

A. Platform

The experiments are performed using a homemade platform, allowing us to have full flexibility in its configuration and location. We use an IEEE 802.11ac (Wi-Fi 5) network operating in the 5 GHz frequency band with a maximum channel bandwidth of 80 MHz. All the technical details on the used Wi-Fi interfaces and the system can be found in [START_REF] Arrabal | Wi-Fi 5 packets capture[END_REF].

In the interest of keeping a high degree of fidelity to an operating IEEE 802.11 WLAN, the choice of transmission parameters is left up to the default rate adaptation algorithm of the Wi-Fi interface. We use iPerf3 to generate UDP/CBR data flows, tshark to capture the exchanged frames, and Bash and awk scripts to automate experiments, and to extract and format data.

We opted for the approach using a Wi-Fi interface configured in monitor mode 4 to collect the Wi-Fi physical and link layer headers. Therefore, our platform consists of three nodes: a source, a destination, and a monitor node. The monitor node acts as a sniffer that simply collects all the frames exchanged in its vicinity, but does not participate in any of the communications. The monitor is placed close to the destination in the effort to closely replicate its wireless environment and the frames it receives. Since the monitor is simply an intermediary for the destination, there is of course an inherent measurement error.

B. Experimental setup and environment

One of the main goals in our data collection campaign was to have a diverse set of measurements reflecting the impact of various environments on the FDR. In this section, we describe the different experimental setups, as well as some of the boundaries of those setups. A single experimental setup on our platform is described by: 1) the physical environment and Line of Sight (LoS); 2) the distance (in meters) between the source and the destination and the position of the nodes' antennas. Each experimental setup with given values for the different parameters (physical environment, LoS, distance, antennas' position) is called a scenario. Every scenario is replicated in three independent runs of 120 seconds.

1) Physical environment and LoS:

The small size of our experimental platform makes it easy to transport in different environments. As a result, we performed experiments in three different floors of our laboratory. In the three environments, we ran experiments with and without the presence of an obstacle obstructing the LoS between the source and destination. Changing the LoS creates diversity in the experiments at a minimal cost.

It should be noted that during all the experiments, traffic from surrounding networks may exist. The 5GHz band is seldom used in our laboratory, however a small amount of traffic may come from other access points in the building and it impacts the collected measurements. In the interest of system fidelity, we decided to keep the environment in its existing conditions and not try to limit this traffic or the other perturbations created by the movement of our colleagues.

The experiments are always conducted in a corridor, however each floor has different characteristics:

• 1st floor: The corridor is slightly curved, and the nodes are placed in a way that guarantees an LoS between them when an obstacle is not present. • 2nd floor: The corridor is straight and empty.

• 3rd floor: The corridor is also straight, but different elements are present (e.g., tables along the walls).

2) Distance and antenna position: The distance between the source and the destination varies between zero and 26 meters. The maximum distance that permits successful communication depends on the chosen environment and the position of the antennas. The antennas can be in one of the four positions depicted in Figure 1. It was our experience that different floors of our building result in different propagation environments due to their shape and construction.

The segment length is set to 1470 Bytes and the CBR datarate is 1 Gbps, corresponding to a saturated flow.

C. Data extraction

At the end of every experimental run, the monitor has collected a trace containing the traffic exchanged between the source and destination nodes 5 . We extract the following transmission parameters from every data frame captured by the monitor:

• Modulation and Coding Scheme (MCS) index:

Higher MCS indexes offer higher datarates, however, they are less robust to a lossy channel. • The link layer throughput computed based on the average frame payload length, denoted as throughput monitor. These measured parameters and computed metrics constitute our database. In our final database used for building the FDR estimation model, each observation is the average of 100ms of experimental data. Typically, over 600 frames are exchanged during each interval, ensuring that this averaged value is representative of the behaviour during the interval.

D. Link quality metric

The performance metrics of interest in this paper is the Frame Delivery Rate (FDR) defined as:

However, calculating the FDR from experimental data is far from a straightforward task as both the number of received and sent frames are difficult to accurately compute from traces. Since the monitor acts as an intermediary for the destination, we regularly observed several undesired events:

• Frames are well received by the monitor, however marked as missing in the BACKs sent by the destination. • Frames are not at all received by the monitor, but they are acknowledged by the destination in the BACKs. • Frames in retransmission do not consistently have the "retry" flag set to True. In an effort to minimize the effect of these events, we decided to calculate two different FDRs: a first one based on the frames that the monitor received from the source node, and a second one based on the BACKs received from the destination node.

Both FDR versions rely on a mechanism called the Sequence Number (SN) 6 . Every data frame has a single SN that is constant during all retransmissions. A stream of frames from the same flow will have increasing SNs. The SN enables the destination to 1) properly reassemble a stream of frames in case they arrive out of order and 2) acknowledge only the well-received frames and indicate the missing frames in the BACKs sent back to the source. The source can then consult each individual BACK to know which frames were not properly received and only retransmit those frames.

Thus, we compute the FDR using two different versions of Eq. 1, one based on the BACK information:

FDR ack = 1 -
Nb. missing frames in BACKs Nb. appeared frames in BACKs ,

and a second one based on data frames' Sequence Numbers (SN):

FDRseq = 1 - Nb. missing SNs for monitor Max received SN -Min received SN . (3)
Equation 3 calculates the FDR based on the number of missing sequence numbers in a series of data frames captured by the monitor and the difference between the maximum and minimum SNs received in that series. Equation 2 calculated the FDR using the number of frames missing and the total number of frames appearing in the BACKs captured by the monitor.

IV. DATA ANALYSIS AND PROPOSED SOLUTION

This section introduces a novel framework for a tracedriven performance evaluation of IEEE 802.11 networks, developed in Python, making it fully open-ource and available online [START_REF] Stojanova | marija-stojanova/hpsr-wifi-link-quality-estimation: HPSR 2023 Link Quality Estimation[END_REF].

We focus on 1) identifying the highest-impact features, i.e., transmission parameters and computed metrics, 2) building an FDR estimation model based on these features, 3) studying its accuracy, and 4) discussing a posteriori insight we can obtain using the model. The collected traces can be divided into six datasets: experiments with LoS named First, Second, and Third floor, and experiments with no LoS named Obstacle 1st, Obstacle 2nd, and Obstacle 3rd.

A. Correlation study

We first study the Pearson's correlation coefficient between the linear and quadratic features and the FDR. Fig. 2 shows these correlations for both versions of FDR for the full dataset containing all six scenarios. FDR ack does not correlate strongly with any of the other features, and this is true even when considering the different scenarios separately. In the interest of brevity, for the remainder of the paper we focus solely on FDR seq .

FDR seq is strongly correlated with throughput monitor and frames monitor. All the Pearson's coefficients higher than 0.4 are a combination of at least one of these features. Surprisingly, the FDR seq is not as highly affected by the distance between the nodes and the RSSI as one would expect. Many existing works have considered the FDR to be a function of these two features, and it is the authors' opinion that the complexity of MIMO has made the relationship between these features highly unpredictable. We remove the BW and NSS features, as they have almost constant values and weak correlations to the FDR.

B. FDR estimation model

The estimation model we aim to develop needs to be computationally simple so that it can be executed on small Wi-Fi-enabled devices. To this effort, we first limit the number of input features used to build the estimation model. We then choose a model that reflects the configuration of our dataset, i.e., a small number of features and a large number of samples.

The model presented has the following characteristics:

• We minimize the number of input features by selecting only the features that have a Pearson's correlation coefficient of at least 0.4 in at least one of the datasets. These are the MCS index, RSSI, throughput monitor, frames monitor, datarate, and distance 7 . • We use Ridge linear regression that has a low complexity in the number of samples, and (potentially) high complexity in the number of features.

Ridge regressions also allow us to circumvent any multicollinearity issues that are potentially present in experimental Wi-Fi traces. • We use five-fold cross validation on 80% of the data, and use the remaining 20% as a test set. • We use grid search to find the optimal α parameter for the Ridge linear regression. We tested all values in the interval [0, 1] with a step of 0.01 and found an optimal α = 0.04. • We use polynomial features, i.e., our input features are a collection of linear features and all their possible quadratic combinations. • We make sure that each scenario is represented with an equal number of samples in the complete dataset 8 . 7 We refer the reader to our Jupyter notebook for a complete correlations heatmap [START_REF] Stojanova | marija-stojanova/hpsr-wifi-link-quality-estimation: HPSR 2023 Link Quality Estimation[END_REF]. 8 Our complete dataset contains 91494 samples with 13082 samples per scenario. A scatter-plot of the estimation vs. ground truth is shown in Fig. 3. We notice that most of the points lie closely to the ideal estimation, with a Mean Absolute Error (MAE) of 0.06 across the datasets, and an R2 score of 0.89. Given the simplicity of the proposed model and its input features, we believe this accuracy level makes it relevant for many Wi-Fi devices and applications.

C. A posteriori data knowledge

We look at the knowledge we can extract from the models built using the different datasets separately. We train one model per dataset and we only keep the features that have a Pearson correlation coefficient larger than 0.4 in that particular dataset 9 . The MAE is shown in Fig. 4, where the y-axis is the name of the dataset used in training, and the x-axis is the name of the test dataset.

There are several key observations we can make:

1) The obstacle datasets are both harder to estimate and less suited for the estimation of other datasets.

2) The datasets First floor and Third floor are the best choice when it comes to estimating the FDR of other datasets. Interestingly, these are also the two simplest models, as they only contain only frames monitor and throughput monitor. 3) As expected, the best choice when it comes to overall estimation is to use the complete dataset named here All Data.

Fig. 4: Mean Absolute Error (MAE) when using the yaxis train dataset to estimate the x-axis test dataset.

V. CONCLUSION

In this work, we have set up a custom-built Wi-Fi 5 experimental platform that collects frames exchanged between a source node and a destination node, with all their headers including the physical and link headers. We have carried out experiments in various environments, distances between nodes, and antennas' positions. We made available a database of all the parameters accessible in the physical and link layer headers, as well as several computed metrics. Using this database, we have provided two methods for computing the FDR and a separate correlation study, from which we have built a simple FDR estimation model. We have shown that 9 Separate training and testing datasets are always used even when working on a single scenario. using the complete dataset and the proposed estimation model, we can estimate the FDR obtained in different environments with high accuracy.

As a future work, we plan to conduct additional experiments in other environments (anechoic chamber, in a fleet of UAVs) and with other traffic profiles (in terms of stream datarate and segment size, and also with TCP) in order to complement the database.

  FDR =Number of well received frames Number of sent frames .

Fig. 1 :

 1 Fig. 1: Possible antenna positions.

Fig. 2 :

 2 Fig. 2: Pearson's feature correlations for the FDR ack and FDR seq .

Fig. 3 :

 3 Fig. 3: Estimated FDR seq vs. ground truth FDR seq .

https://wiki.wireshark.org/CaptureSetup/WLAN

We recall that the monitor's trace is merely an approximation of the real stream of traffic exchanged between the source and destination.

We only present a high-level explanation of this mechanism, for a detailed review we refer the reader to the IEEE 802.11 standard[START_REF]IEEE Standard for IT Part 11: Wireless LAN MAC and PHY Specifications[END_REF]

VI. ACKNOWLEDGEMENTS

The authors would like to thank the ANR and the AID for their funding support via the CONCERTO project (ANR-20-ASTR-0003) and the CNRS for funding a research engineer position.


Stuti Jain

Senior Member, IEEE Emilie Chouzenoux

Member, IEEE Kriti Kumar

Senior Member, IEEE Angshul Majumdar

Graph Regularized Probabilistic Matrix Factorization for Drug-Drug Interactions Prediction

Keywords: Matrix factorization, Probabilistic matrix factorization, Graph regularization, Drug-drug interaction prediction

Co-administration of two or more drugs simultaneously can result in adverse drug reactions. Identifying drug-drug interactions (DDIs) is necessary, especially for drug development and for repurposing old drugs. DDI prediction can be viewed as a matrix completion task, for which matrix factorization (MF) appears as a suitable solution. This paper presents a novel Graph Regularized Probabilistic Matrix Factorization (GRPMF) method, which incorporates expert knowledge through a novel graph-based regularization strategy within an MF framework. An efficient and sounded optimization algorithm is proposed to solve the resulting non-convex problem in an alternating fashion. The performance of the proposed method is evaluated through the DrugBank dataset, and comparisons are provided against state-of-the-art techniques. The results demonstrate the superior performance of GRPMF when compared to its counterparts.

I. INTRODUCTION

Drug-Drug Interaction (DDI) refers to the effects of a given drug when taken together with another drug, at the same time. Co-administration of two or more drugs simultaneously can affect the pharmokinetics and/or pharmacodynamics of one or more drugs, which can cause unexpected and even adverse drug reactions [START_REF] Vilar | Detection of drug-drug interactions through data mining studies using clinical sources, scientific literature and social media[END_REF]. These effects can cause severe injuries to the patients and even be responsible for deaths. Thus, it is necessary to know the DDI for the drugs used in the market, for clinical safety. Knowledge of DDI is also vital for developing new drugs, as well as for repurposing old drugs from clinical and public health perspectives. Preclinical identification of DDIs is an ill-posed problem as clinical testing (e.g., in vitro, in vivo, and in populo) is usually conducted on a small group of drugs. Such a process is time and cost intensive. Thus, DDI prediction employing computational approaches that are implementable on a large scale has become a popular research topic in recent years [START_REF] Roblek | Drug-drug interaction software in clinical practice: a systematic review[END_REF].

Computational approaches for DDI prediction can be broadly divided into two categories: (i) Similarity-based approaches, that are based on the similarity of drug information (e.g., chemical structure [START_REF] Vilar | Similarity-based modeling in large-scale prediction of drug-drug interactions[END_REF], targets [START_REF] Luo | Ddi-cpi, a server that predicts drug-drug interactions through implementing the chemical-protein interactome[END_REF], side-effects [START_REF] Shi | Predicting comprehensive drug-drug interactions for new drugs via triple matrix factorization[END_REF], and (ii) knowledge-based approaches, that employ text S. Jain and E. Chouzenoux are with CVN, Inria Saclay, Univ. Paris Saclay, 91190 Gif-sur-Yvette, France. e-mail: (emilie.chouzenoux@centralesupelec.fr, stuti.jain@inria.fr).

A. Majumdar and K. Kumar is with Dept. of ECE, IIIT -Delhi, India, 110020. e-mail: (angshul@iiitd.ac.in, kritik@iiitd.ac.in).

This work received support from the Associate Team COMPASS between Inria and IIIT Delhi. E.C. and S.J. acknowledge support from the European Research Council Starting Grant MAJORIS ERC-2019-STG-850925. mining from scientific literature [START_REF] Zhang | Leveraging syntactic and semantic graph kernels to extract pharmacokinetic drug drug interactions from biomedical literature[END_REF], electronic medical record database [START_REF] Duke | Literature based drug interaction prediction with clinical assessment using electronic medical records: novel myopathy associated drug interactions[END_REF] and the FDA Adverse Event Reporting System [START_REF]U.s. food and drug administration[END_REF] to predict DDI. Note that the latter group of approaches do not perform learning per se, as it is mostly a tool to mine clinical findings. Although different machine learning and artificial intelligence models have been used to address the problem [START_REF] Han | A review of approaches for predicting drug-drug interactions based on machine learning[END_REF], DDI prediction remains a challenging problem to address.

In this work, the DDI prediction task is formulated as a matrix completion problem, involving a symmetric matrix with rows/columns corresponding to drugs. We adopt a matrix factorization paradigm, when the sought matrix is defined as the product of two latent factors satisfying some prior knowledge. Our contribution lies in the construction of an original prior well suited to DDI prediction. It incorporates expert knowledge on the DDIs within a graph-based regularization term similar to [START_REF] Strahl | Scalable probabilistic matrix factorization with graph-based priors[END_REF], with the aim to favor expected (dis)similarities between drug pairs. This formulation results in the so-called Graph Regularized Probabilistic Matrix Factorization (GRPMF) method, for which we also propose a sounded optimization algorithm relying on modern proximal methods. We evaluate the performance of our method on the DrugBank dataset, and present comparisons against state-ofthe-art techniques.

The paper is organized as follows. Section II discusses related works on DDI prediction. Section III provides a brief overview of various MF-based methods. Section IV presents our main contribution, that is the proposed GRPMF formulation and the optimization algorithm to resolve it. Section V presents our experimental results and comparisons with the state-of-the-art methods, and finally, Section VI concludes the work.

II. RELATED WORKS

Clinical trials are time and money consuming. Many DDIs thus remain unknown, mostly because they are not tested during trials. Machine learning-based methods have been widely investigated to predict the unobserved DDIs. A heterogeneous network-assisted inference (HNAI) method is proposed in [START_REF] Cheng | Machine learning-based prediction of drug-drug interactions by integrating drug phenotypic, therapeutic, chemical, and genomic properties[END_REF], gathering five prediction models (naive Bayes, decision tree, k-nearest neighbor, logistic regression, and support vector machine) to perform DDIs prediction from drug phenotypic, therapeutic, structural, and genomic similarities. The work in [START_REF] Rohani | Drug-drug interaction predicting by neural network using integrated similarity[END_REF] presents a neural network (NN) based method that proposes a heuristic selection of several drug similarity scores integrated with a nonlinear similarity fusion strategy to obtain abstract features for DDI prediction. Another work [START_REF] Yan | Predicting drug-drug interactions based on integrated similarity and semisupervised learning[END_REF] proposes a semi-supervised learning method for DDI prediction, that uses drug chemical, biological, and phenotype data to calculate the feature similarity of drugs using a regularized least square score minimization. Some works address the DDI prediction problem as an edge detection (i.e., link prediction) problem where the edges to infer represent connections between drugs. The work in [START_REF] Kastrin | Predicting potential drug-drug interactions on topological and semantic similarity features using statistical learning[END_REF] presents both unsupervised and supervised techniques for link prediction using binary classifiers such as tree, knearest neighbors, support vector machine, random forest, and gradient boosting machine based on topological and semantic similarity features to estimate the drug interactions. Another work [START_REF] Shtar | Detecting drug-drug interactions using artificial neural networks and classic graph similarity measures[END_REF] proposes two methods based on NNs and factor propagation over graph nodes, namely, adjacency matrix factorization (AMF) and adjacency matrix factorization with propagation (AMFP) for link prediction for discovering DDIs.

The superior performance of deep learning (DL) techniques across different domains has triggered the interest in such techniques to estimate drug interactions. The work in [START_REF] Xu | Leveraging biomedical resources in bi-lstm for drug-drug interaction extraction[END_REF] presents a biomedical resource LSTM (BR-LSTM) that combines biomedical resources with lexical information and entity position information together to extract DDI from the biomedical literature. Note that this model is not DDI prediction per se, but only an automatic tool for mining of information from clinical literature. The work in [START_REF] Shukla | Efficient prediction of drug-drug interaction using deep learning models[END_REF] proposes a convolutional mixture density recurrent NN model that integrates convolutional neural networks, recurrent NNs, and mixture density networks for DDI prediction. An autoencoder-based semi-supervised learning algorithm for feature extraction from FDA adverse event reports to identify potential high priority DDIs for medication alerts is presented in [START_REF] Liu | Semi-supervised learning algorithm for identifying high-priority drug-drug interactions through adverse event reports[END_REF]. Another work [START_REF] Lee | Novel deep learning model for more accurate prediction of drug-drug interaction effects[END_REF] employs autoencoders and a deep feed-forward network trained on structural similarity profiles (SSP), Gene Ontology (GO) term similarity profiles (GSP), and target gene similarity profiles (TSP) of known drug pairs to predict the effects of DDIs. Due to the black-box nature of the DL models, some work has been done on seeking for explainable DL-based DDI techniques. A comprehensive review of the explainable AIbased techniques to promote the trust of AI models for the critical task of DDI prediction is presented in [START_REF] Vo | On the road to explainable ai in drug-drug interactions prediction: A systematic review[END_REF].

Recently, attention-based deep learning models have been explored for DDI prediction task that can predict novel drugdrug interactions and drug-drug interaction-associated events [START_REF] Schwarz | Attentionddi: Siamese attention-based deep learning method for drug-drug interaction predictions[END_REF], [START_REF] Pan | Multisource attention-mechanism-based encoder-decoder model for predicting drug-drug interaction events[END_REF], [START_REF] Liu | Enhancing drug-drug interaction prediction using deep attention neural networks[END_REF]. The work in [START_REF] Schwarz | Attentionddi: Siamese attention-based deep learning method for drug-drug interaction predictions[END_REF] presented a Siamese selfattention multi-modal neural network for DDI prediction that utilizes multiple drug similarity measures derived from drug characteristics, including drug targets, pathways, and gene expression profiles. Here, model explainability is offered via an Attention mechanism for identifying salient input features. Another work in [START_REF] Pan | Multisource attention-mechanism-based encoder-decoder model for predicting drug-drug interaction events[END_REF] proposed two attention-mechanism-based encoder-decoder models that incorporate multisource information derived from drugs. They considered three encoding methods to encode the drugs that consisted of (i) integrating the self-attention mechanism, (ii) cross-attention mechanism, and (iii) graph attention network to construct a multisource feature fusion network. Different encoders are constructed based on the different features to learn the interrelationship information between drug pairs, and that information is used to predict DDI. The work in [START_REF] Liu | Enhancing drug-drug interaction prediction using deep attention neural networks[END_REF] presented a Deep Attention Neural Network based DDI prediction framework, abbreviated as DANN-DDI that uses a graph embedding method to learn drug representations from multiple drug feature networks and adopts an attention neural network to learn representations of drug-drug pairs that are finally used to predict DDI.

Another direction of work formulates the DDI prediction problem as a Matrix Completion task. Here, given the partially observed DDI matrix, the task is to compute the unobserved interactions between the drugs. Some of the popular generic (not tailored for DDI) MF techniques are (i) singular value decomposition (SVD) [START_REF] Sarwar | Application of dimensionality reduction in recommender systems: A case study[END_REF], (ii) non-negative matrix factorization (NMF) [START_REF] Lee | Learning the parts of objects by nonnegative matrix factorization[END_REF] and, (iii) probabilistic matrix factorization (PMF) [START_REF] Salakhutdinov | Probabilistic matrix factorization[END_REF]. We will present the two later approaches in detail in our next section, as MC/MF constitutes the core of our contribution.

In addition to the conventional binary DDI prediction, the work in [START_REF] Yu | Predicting and understanding comprehensive drug-drug interactions via semi-nonnegative matrix factorization[END_REF] presents an NMF-based approach utilizing drug features for comprehensive DDI prediction. Here, the comprehensive DDI matrix is a signed binary matrix with +1 for enhancive drugs, -1 for degressive drugs, and 0 for no drug interactions, respectively, which is rather useful to predict the (positive/negative) behaviors of the interacting drugs. The work in [START_REF] Zhu | Attribute supervised probabilistic dependent matrix tri-factorization model for the prediction of adverse drug-drug interaction[END_REF] presents an attribute supervised learning model probabilistic dependent matrix tri-factorization (PDMTF) approach for adverse DDI prediction. They utilized two drug attributes, molecular structure, side effects, and their correlation to compute the adverse interactions among drugs. The work in [START_REF] Zhang | Manifold regularized matrix factorization for drug-drug interaction prediction[END_REF] introduces a manifold regularized MF (MRMF) technique to predict DDIs using drug similarities based on drug features like substructures, targets, enzymes, transporters, pathways, indications, side effects, and off side effects. Different drug features and similarity measures are considered to calculate different manifold regularization terms to construct the MRMF models. The inclusion of drug featurebased manifold regularization in the matrix factorization formulation is reported to provide more accurate and robust predictions.

The publicly available large structured biomedical databases has enabled the use of knowledge graph (KG) based approaches for different applications in the biomedical domain. KGs are used to synthesize large biomedical graphs that map similar drug-related entities in the drug database. The work in [START_REF] Celebi | Evaluation of knowledge graph embedding approaches for drug-drug interaction prediction in realistic settings[END_REF] uses KGs embeddings, namely, RDF2Vec, TransE, TransD, and machine learning algorithms for DDI prediction. A KG NN method (KGNN) that captures the drug and its potential neighborhoods by mining their associated relations in KG for DDI prediction is proposed in [START_REF] Lin | Kgnn: Knowledge graph neural network for drug-drug interaction prediction[END_REF]. This method utilizes the drugs' topological structures in KG for potential DDI prediction. Another work [START_REF] Karim | Drug-drug interaction prediction based on knowledge graph embeddings and convolutional-lstm network[END_REF] utilizes KGs combined with DL techniques for estimating DDIs. This work considers the DDI matrix and KG in the form of learned embeddings (like ComplEx, TransE, RDF2Vec, etc.) as input to the Convolutional Neural Networks (CNN) and Long-Short Term Memory (LSTM) model to predict DDIs.

In this work, we focus on the MC/MF based framework, as it presents the advantage of being non supervised and highly interpretable. Our contribution is to incorporate expert knowledge within this family of approach, so as to take advantage of the aforementioned progressed in database availability.

III. BACKGROUND

This section presents an overview of MF techniques for MC. We choose here to remain in a generic setting where the matrix to complete is real-valued and rectangular. Note that, for the DDI task, the sought matrix is square symmetric and, in most cases, binary valued, which might lead to simplified formulations.

A. Matrix Completion Problem

Let us consider the problem of a full matrix R ∈ R N ×M to recover from partially known matrix

Y ∈ R N ×M . Let D = {i ∈ {1, . . . , N }, j ∈ {1, . . . , M } s.t. (i, j) is observed}.
Non observed entries are typically set to zero. The masking of the indexes outside the set D is modeled through a Hadamard product ⊙ with a matrix B ∈ {0, 1} N ×M , such that B ij = 1 if (i, j) ∈ D, and B ij = 0 otherwise. The partially known matrix Y can be expressed as:

Y = B ⊙ R. (1)
The task of matrix completion amounts to recovering the entries of R that do not belong to the set of observed indexes D.

B. Matrix Factorization (MF)

MF [START_REF] Koren | Matrix factorization techniques for recommender systems[END_REF] consists of recovering missing entries in matrix R by minimizing the simple least-squares function minimize

R ∥Y -B ⊙ R∥ 2
F under some specific structural prior constraints on R. Under MF prior, R is recast as a product of two matrices U ∈ R N ×Z and V ∈ R Z×M , where Z ≥ 1 defines a latent space dimension, typically low compared to (N, M ). The matrices U and V are inferred by solving:

minimize U,V ∥Y -B ⊙ (U V )∥ 2 F . (2)
Subsequently, the complete matrix R is simply recovered by R = U V . Problem ( 2) is however highly under-determined and extra priors are typically introduced to obtain meaningful solutions. The most widely used being probably the positivity of the entries of the latent factors (U, V ), yielding the NMF (nonnegative MF) formulation [START_REF] Naik | Non-negative Matrix Factorization Techniques: Advances in Theory and Applications[END_REF].

Let us discuss related formulations for MC. First, another formulation strategy to impose low rank is to resort to nuclear norm minimization [START_REF] Recht | A simpler approach to matrix completion[END_REF]. Regularization strategies, based on graph modeling, have been considered in [START_REF] Gu | Collaborative filtering: Weighted nonnegative matrix factorization incorporating user and item graphs[END_REF] for the MF formulation and in [START_REF] Mongia | Matrix completion on multiple graphs: Application in collaborative filtering[END_REF] for the nuclear norm formulation. MF models with more than two factors lead to the so-called deep MF approach, investigated for instance in [START_REF] Mongia | Deep latent factor model for predicting drug target interactions[END_REF] in the context of MC. Graph regularized version for deep MC has been proposed in the recent work [START_REF] Mongia | Deepvir: Graphical deep matrix factorization for in silico antiviral repositioningapplication to covid-19[END_REF].

C. Probabilistic Matrix Factorization (PMF)

PMF introduces probabilistic models on the latent factors U and V in the MF formulation [START_REF] Salakhutdinov | Probabilistic matrix factorization[END_REF]. More precisely:

• Each observed entry (Y ij ) (i,j)∈D is assumed to follow a Gaussian distribution, with mean [U V ] i,j and variance σ 2 (positive scalar assumed to be known). • Each entry (U iz ) 1≤i≤N,1≤z≤Z is assumed to follow a Gaussian distribution with zero mean and variance σ 2 U (positive scalar assumed to be known).

• Each entry (V zj ) 1≤z≤Z,1≤j≤M is assumed to follow a Gaussian distribution with zero mean and variance σ 2 V (positive scalar assumed to be known). The maximum a posteriori (MAP) estimator of (U, V ) given Y , associated with the above model can be obtained by solving: minimize

U,V 1 2σ 2 ∥Y -B⊙(U V )∥ 2 F + 1 2σ 2 U ∥U ∥ 2 F + 1 2σ 2 V ∥V ∥ 2 F . (3)
The minimization with respect to U (resp. V ) in this formulation amounts to invert a linear system, which can be performed using a conjugate gradient solver. The PMF formulation can be enhanced by incorporating correlated Gaussian distributions, which results in the PMFG formulation described hereafter.

D. Probabilistic Matrix Factorization with Graph regularization (PMFG)

In PMFG, the prior distributions of U and V now include a graph regularization strategy [START_REF] Strahl | Scalable probabilistic matrix factorization with graph-based priors[END_REF], which amounts to inferring correlations along the rows (resp. columns) of U (resp. V ) jointly with the factors U and V . These correlations are modeled through two precision matrices, Γ U ∈ S ++ N and Γ V ∈ S ++ M , where S ++ N ∈ R N ×N and S ++ M ∈ R M ×M denote symmetric positive definite matrices while S N , S M denote symmetric matrices. The prior is the following:

• The columns u z ∈ R N , z ∈ {1, . . . , Z}, of U are independent realizations of a multivariate Gaussian distribution with zero mean and covariance

C U = Γ -1 U ; • The lines v z ∈ R M , z ∈ {1, . . . , Z} are independent
realizations of a multivariate Gaussian distribution with zero mean and covariance C V = Γ -1 V . The precision matrices Γ U and Γ V are related to Gaussian graphical models associated to the two underlying Gaussian distributions [cite the book Elements-Statistical-Learning-Inference-Prediction], which justifies the name for "graph" regularization. Specifically, matrix Γ U (resp. Γ V ) can be understood as the adjacency matrix of an undirected graph where each edge identifies with two entries of u z (resp. v z ) being correlated, given all the others.

The MAP estimate can now be obtained by solving the following:

minimize U, V, ΓU , ΓV 1 2σ 2 ∥Y -B ⊙ (U V )∥ 2 F + 1 2 tr(U ⊤ Γ U U ) + 1 2 L(Γ U ) + 1 2 tr(V Γ V V ⊤ ) + 1 2 L(Γ V ), (4)
where tr(•) denotes the trace operation and

(∀Γ U ∈ S N ) L(Γ U ) = -ln det(Γ U ) if Γ U ∈ S ++ N +∞ otherwise (5) with det(•) the determinant operation.
PMFG approach provides promising results in [START_REF] Strahl | Scalable probabilistic matrix factorization with graph-based priors[END_REF]. However, it does not incorporate any physical-oriented knowledge in the sought factors (U, V ). Actually, in many applications, such as DDI, expert knowledge is available, that dictates more or less likely correlations among the variables. The aim of this present work is to propose a novel formulation to account for such prior knowledge, within the PMFG paradigm.

IV. PROPOSED GRAPH REGULARIZED PROBABILISTIC MATRIX FACTORIZATION (GRPMF)

Let us specify our targeted application. We focus on solving the matrix completion problem arising when predicting interactions between the different drugs, which is the socalled DDI problem. Let Y ∈ R N ×N a partially known drug interaction matrix (with unobserved entries set to 0) for N different drugs. The aim is to recover the full drug interaction matrix R ∈ R N ×N1 . The proposed method includes expert knowledge within the formulation of PMFG, with the aim of estimating DDIs. Note that unlike general matrix completion problems discussed in previous section, in DDI the matrix to infer is square and symmetric. The sought interaction matrix R can thus be factored as:

R = U U ⊤ . (6)
Under this new setting, the graphical model is reduced to a single graph with adjacency matrix Γ U , and the previously presented PMFG formulation can simply be modified as:

minimize U, ΓU 1 2σ 2 ∥Y -B⊙(U U ⊤ )∥ 2 F + 1 2 tr(U ⊤ Γ U U )+ 1 2 L(Γ U ). (7)

A. Integrating Expert Knowledge

As already discussed, in DDI application, one might have prior knowledge about the position of the graph edges (i.e. non zero elements in the precision matrix Γ U ), thanks to some expert analysis of the database. This prior knowledge is available through an extra symmetric matrix with positive real entries, namely A U ∈ [0, +∞[ N ×N . For instance, for the DDI prediction task, A U could result from precomputing the similarity between the N drugs of the dataset in terms of the SIMCOMP (SIMilar COMPound) scores [START_REF] Hattori | Simcomp/subcomp: chemical structure search servers for network analyses[END_REF]. Such matrix can then be used as a structural prior on the sought matrix Γ U , so as to remove spurious edges with no physical meaning, and to promote expected ones in the restored graph. Otherwise stated, if for some (i, j)

with i ̸ = j, [A U ] ij is large, then [Γ U ] ij = [Γ U ]
ji should be encouraged to be high as well. In contrast, an entry [A U ] ij close or equal to zero should promote the removal of the edge between nodes i and j in the sought graph.

In order to build a suitable regularization function associated to this new prior, let us introduce the following sets:

E U = {(i, j) ∈ {1, . . . , N } 2 , i ̸ = j and [A U ] ij > τ } (8)
E U = {(i, j) ∈ {1, . . . , N } 2 , i ̸ = j and [A U ] ij ≤ τ } (9)
with τ ≥ 0 a given detection threshold. In order to avoid a tedious non-convex coupling term in the resulting minimization problem, we furthermore introduce the proxy variable X ∈ R N ×N that we penalize so as to be close to the sought product U U ⊤ . We thus propose to solve: minimize

X, U, ΓU 1 2σ 2 ∥Y -B ⊙ X∥ 2 F + 1 2 tr(U ⊤ Γ U U ) + λ U (i,j)∈E U |[Γ U ] ij | -λ U (i,j)∈E U ln(|[Γ U ] ij | + δ) + λ R 2 ∥X -U U ⊤ ∥ 2 F + 1 2 L(Γ U ) + λ U 2 ∥Γ U ∥ 2 F . ( 10
)
Hereabove, parameter λ R > 0 controls the fulfillment of the equality constraint X = U U ⊤ while parameter λ U > 0 controls the regularization imposed on the precision matrix Γ U . Our proposed regularization term is made of two parts.

We introduce an ℓ 1 term to promote sparsity on the regions where edges should not appear (i.e. E U ), and a log-barrier term, smoothed by δ > 0, in regions where edges (i.e., nonzero entries) should be promoted (i.e. E U ). Finally, a quadratic term, that can be viewed as an elastic-net penalty, is added in order to avoid too large values in the entries of the sought covariance matrix.

B. Optimization Algorithm

Problem ( 10) is highly non-convex, as it is commonly the case in MF formulations. We propose to use an alternating optimization strategy to solve it. Let F (X, U, Γ U ) denote the loss function presented in [START_REF] Strahl | Scalable probabilistic matrix factorization with graph-based priors[END_REF]. Starting from a given initialization (X 0 , U 0 , Γ 0 U ), for every iteration k ∈ N, the parameters of the algorithm are updated as:

           X k+1 = argmin X∈R N ×N F (X, U k , Γ k U ), U k+1 = argmin U ∈R Z×Z F (X k+1 , U, Γ k U ), Γ k+1 U = argmin Γ U ∈S N F (X k+1 , U k+1 , Γ U ). (11)
Such procedure ensures the monotonical decrease of function F . We explicit hereafter each update. Let us remark that, for variables U and Γ U , the subproblems remain non-convex.

The proposed subroutines only amounts to finding a stationary point for each, which might not be a global minimum. No numerical instabilities were observed in our experiments though. 1) Update of X: The matrix X is updated as:

X k+1 = argmin X∈R N ×N 1 2σ 2 ∥Y -B⊙X∥ 2 F + λ R 2 ∥X-U k (U k ) ⊤ ∥ 2 F . ( 12
)
This is a strictly convex quadratic problem, whose solution satisfies the following optimality condition:

1 σ 2 Y + λ R U k (U k ) ⊤ = 1 σ 2 B ⊙ X k+1 + λ R X k+1 . (13)
The above equation is a linear system that can be solved efficiently by conjugate gradient [START_REF] Shewchuk | An introduction to the conjugate gradient method without the agonizing pain[END_REF].

2) Update of U : The matrix U is updated as:

U k+1 = argmin U ∈R Z×Z 1 2 tr(U ⊤ Γ k U U ) + λ R 2 ∥X k+1 -U U ⊤ ∥ 2 F . (14)
The above optimization problem is non-convex and differentiable. An efficient nonlinear conjugate gradient method is proposed in [START_REF] Duan | Low rank approximation of the symmetric positive semidefinite matrix[END_REF] to address it.

3) Update of Γ U : The matrix Γ U is updated using:

Γ k+1 U = argmin Γ U ∈S N tr((U k+1 ) ⊤ Γ U U k+1 )+L(Γ U )+λ U ∥Γ U ∥ 2 F + 2λ U (i,j)∈E U |[Γ U ] ij | -2λ U (i,j)∈E U ln(|[Γ U ] ij | + δ) (15)
or, equivalently,

Γ k+1 U = argmin Γ U ∈S N f (Γ U ) + g(Γ U ) (16)
with

(∀Γ U ∈ S N ) f (Γ U ) = tr((U k+1 ) ⊤ Γ U U k+1 )+L(Γ U ) (17)
and

(∀Γ U ∈ S N ) g(Γ U ) = 2λ U (i,j)∈E U |[Γ U ] ij | -2λ U (i,j)∈E U ln(|[Γ U ] ij | + δ) + λ U ∥Γ U ∥ 2 F . ( 18
)
The minimization of f +g does not have a close form solution and hereagain, an inner solver is required. Function f is convex, differentiable on its domain S ++ N , while g is nonconvex, non-differentiable. Luckily, the latter is separable over each of the entries of Γ U , that is:

(∀Γ U ∈ S N ) g(Γ U ) = 1≤i,j≤N g ij ([Γ U ] ij ) (19)
with, for every (i, j) ∈ {1, . . . , N } 2 ,

(∀ω ∈ R) g ij (ω) = 2λ U |ω| + λ U ω 2 if (i, j) ∈ E U , -2λ U ln(|ω| + δ) + λ U ω 2 if (i, j) ∈ E U . (20)
We thus opt for running L ≥ 1 iterations of a proximal gradient algorithm [START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized Gauss-Seidel methods[END_REF], initialized using the previous value Γ k U ∈ S ++ N (by construction). This reads as follows:

Γ (0) U = Γ k U For ℓ = 1, 2, . . . , L Γ (ℓ) U = Γ (ℓ) U -θ (ℓ) ∇f (Γ (ℓ) U ) Γ (ℓ+1) U = prox θ (ℓ) g ( Γ (ℓ) U ) Γ k+1 U = Γ (L) U . (21)
Hereabove, (θ (ℓ) ) 1≤ℓ≤L is a sequence of positive stepsizes obtained through a suitable backtracking strategy, so that all iterates remain in the (open) domain of f . Moreover, prox denotes the proximity operator, for which a definition in the non convex case can be found in [START_REF] Attouch | Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized Gauss-Seidel methods[END_REF].

Hereafter, we provide the expression for the gradient of f and the proximity operator of g. First,

(∀Γ U ∈ S ++ N ) ∇f (Γ U ) = U k+1 (U k+1 ) ⊤ -Γ -1 U , ( 22
)
while it is not defined for non definite positive matrices of S N . Second, due to the separability of function g, we have, for any θ > 0 [START_REF] Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF],

(∀Γ U ∈ S N ) prox θg (Γ U ) = prox θgij ([Γ U ] ij ) 1≤i,j≤N . (23)
The expression for the proximity operator of each term g ij , defined in [START_REF] Vo | On the road to explainable ai in drug-drug interactions prediction: A systematic review[END_REF], depends if (i, j) ∈ E U or not.

For (i, j) ∈ E U , g ij is a convex, proper, lower semicontinuous function on R. Its proximity operator is thus uniquely defined, and it reads:

(∀ω ∈ R) prox θgij (ω) = argmin ξ∈R 1 2 (ξ -ω) 2 +θ(2λ U |ξ| + λ U ξ 2 ) = prox (1+2λ U θ) -1 |•| (1 + 2λ U θ) -1 ω , = S (1+2λ U θ) -1 (1 + 2λ U θ) -1 ω , (24)
with S τ the soft thresholding operator with parameter τ > 0:

(∀u ∈ R) S τ (u) = sign(u) max(0, |u| -τ ). (25)
For (i, j) ∈ E U , g ij is non-convex. Its proximity operator reads:

(∀ω ∈ R) prox θgij (ω) = argmin ξ∈R 1 2 (ξ -ω) 2 +θ(-2λ U ln(|ξ| + δ) + λ U ξ 2 ) , (26)
= prox 2θλ U 1+2θλ U ζ 1 1 + 2θλ U ω , (27)
with ζ : ω → -ln(|ω| + δ). Function ζ is non-convex. It is lower semi-continuous, and lower bounded by polynomial ω → -(ω 2 + δ), thus its proximity operator exists though it is not uniquely defined (i.e., it is set-valued). Let τ > 0, and

ω ∈ R. Set φ : ω → τ ζ(ω) + 1 2 (ω -ω) 2 .
The operator prox τ ζ evaluated at ω is defined as the set of (global) minimizers of φ. Hereafter, we study φ, for positive or negative input values, so as to deduce the proximity set of τ ζ.

First, for every ω > 0, the first and second order derivatives of φ read:

φ ′ (ω) = - τ ω + δ + ω -ω, (28)
and

φ ′′ (ω) = τ (ω + δ) 2 + 1. ( 29
)
Function φ ′ is strictly increasing on ]0, +∞[. Moreover, 28) is equivalent to search the roots for the second order polynomial:

lim ω→0 + φ ′ (ω) = -τ δ -ω. Thus, φ ′ cancels on ]0, +∞[ if ω > -τ δ . Canceling (
-τ + (ω -ω)(ω + δ) = ω 2 -(ω -δ)ω -ωδ -τ. (30)
The discriminant of such polynomial is

∆ + = (ω -δ) 2 + 4(ωδ + τ ) (31) = (ω + δ) 2 + 4τ, (32)
which is strictly positive under the condition ω > -τ δ . Two roots exist and read:

ω + 1,2 = (ω -δ) ± √ ∆ + 2 . ( 33
)
There is a sole positive root, equals to:

ω + 2 = (ω -δ) + (ω + δ) 2 + 4τ 2 . ( 34
)
Second, when ω ∈] -∞, 0[, we have

(∀ω < 0) φ ′ (ω) = τ -ω + δ + ω -ω, (35)
and

(∀ω < 0) φ ′′ (ω) = τ (ω -δ) 2 + 1. (36)
Hereagain,

φ ′ is strictly increasing. Moreover, lim ω→0 -φ ′ (ω) = τ δ -ω. Therefore, φ ′ cancels on ] -∞, 0[ if ω < τ δ .
Canceling [START_REF] Recht | A simpler approach to matrix completion[END_REF] is equivalent to search the roots for the second order polynomial:

-τ + (ω -ω)(ω -δ) = ω 2 -(ω + δ)ω + ωδ -τ. ( 37
)
The discriminant of such polynomial is

∆ -= (ω + δ) 2 -4(ωδ -τ ) = (ω -δ) 2 + 4τ, (38)
which is strictly positive under the condition ω < τ δ . Two roots exist and are given by:

ω - 1,2 = (ω + δ) ± √ ∆ - 2 . ( 39
)
There is a sole negative root, equals to

ω - 1 = (ω + δ) -(ω + δ) 2 -4(ωδ -τ ) 2 . ( 40
)
Finally, studying the variations of φ on all R allows us to distinguish three cases:

• Case 1: ω ∈] -τ δ , τ δ [. Then, prox τ ζ (ω) = argmin p∈{ω - 1 ,ω + 2 } φ(p). (41)
• Case 2: ω ≤ -τ δ . Then,

prox τ ζ (ω) = ω - 1 . ( 42
)
• Case 3: ω ≥ τ δ . Then,

prox τ ζ (ω) = ω + 2 . ( 43
)
4) Summarized algorithm: The proposed Graph Regularized Probabilistic Matrix Factorization (GRPMF) method for DDI prediction task is summarized in Algorithm 1. Given the partially observed drug interaction matrix Y , we recall that we aim to recover the full drug interaction matrix R by exploiting the drug similarity information captured in A U . Matrices X and Γ U are initialized using identity matrices scaled by a positive value s 0 . The first Z left-singular vectors of Y obtained using singular value decomposition (SVD) are considered for initialization of U . The three unknowns, namely X, U and Γ U , are updated using the updates defined in the previous subsection, in an iterative manner for K iterations. The final matrix R is recovered using R = U U ⊤ .

Algorithm 1 GRPMF for DDI Prediction 1: Input: Y , AU 2: Parameters: Z, σ, s 0 , λR, λU , K, L. 3: Initialization: U 0 (using svd(Y )), X 0 = Γ 0 U = s 0 × IN . 4: Compute EU and EU using ( 8) and ( 9). 5: for k = 1, 2, . . . , K iterations 6: Update X k+1 using (13); Update ΓU k+1 using (21). 9: end 10: Return: R = U K (U K ) ⊤ .

V. EXPERIMENTAL RESULTS

This section presents our experimental results illustrating the validity of the proposed method. We first introduce the dataset considered for DDI prediction along with the necessary data pre-processing steps carried out. Subsequently, comparison of the proposed GRPMF method against the benchmark algorithms, and an ablation study for GRPMF, are presented.

A. Dataset

We rely on the DDI data from Stanford University [START_REF]Drug-drug interaction network[END_REF] that is approved by the U.S. Food and Drug Administration. It contains 48, 514 interactions from 1, 514 drugs extracted from drug labels and scientific publications [START_REF] Wishart | Drugbank 5.0: a major update to the drugbank database for 2018[END_REF], [START_REF]Miner: Gigascale multimodal biological network[END_REF]. Our experiments rely on the publicly available drug dataset from DrugBank that gathers 14, 315 drugs along with their KEGG ID or compound ID. In particular, this latter dataset is used to build the expert knowledge for our regularization term. To do so, we rely on the SIMCOMP score that amounts to comparing the chemical structures of the drugs [START_REF] Hattori | Simcomp/subcomp: chemical structure search servers for network analyses[END_REF]. SIMCOMP computes the similarity of two chemical compounds by counting the number of matched atoms in those atom alignments, using the KEGG ID reference system. We computed the SIMCOMP score using the SIMCOMP search tool with a cutoff of 0.01 [START_REF]Simcomp search[END_REF]. Due to the unavailability of KEGG ID, and thus SIMCOMP scores, for few drugs, and to the non overlap between [START_REF]Drug-drug interaction network[END_REF] and [START_REF] Wishart | Drugbank 5.0: a major update to the drugbank database for 2018[END_REF] drug lists, the final data reduced to N = 927 drugs. For this subset of drugs, we made sure that they each have at least 10 known interactions with other drugs.

From our dataset curating, we end up with a binary matrix Y ∈ {0, 1} N ×N , where y i,j = 1 if drugs i and j are known to interact. An entry y i,j = 0 means that no interaction has been reported between drugs i and j so far. Our SIMCOMP calculation yields a drug similarity data represented as a symmetric matrix A U ∈ R N ×N . The computation for matrix A U is only required once. On our computational device, it required about 12 hours for 927 drugs on this dataset. Note that the initial SIMCOMP matrix has been sparsified so as to retain only the p-nearest neighbours of each drug with the aim to preserve the local geometry of the original data. Such operation promotes that drugs that are close to one another in the original (chemical structure) space should also be close to one another in the learnt (drug interaction) manifold (i.e., local invariance assumption) [START_REF] Ezzat | Drug-target interaction prediction with graph regularized matrix factorization[END_REF]. We then compute sets E U , and E U using ( 8) and ( 9) with τ = 0.

We then adopt a supervised learning paradigm, as in [START_REF] Celebi | Evaluation of knowledge graph embedding approaches for drug-drug interaction prediction in realistic settings[END_REF], [START_REF] Lin | Kgnn: Knowledge graph neural network for drug-drug interaction prediction[END_REF], [START_REF] Karim | Drug-drug interaction prediction based on knowledge graph embeddings and convolutional-lstm network[END_REF]. 20% drug pairs are considered for training (i.e., assumed to be observed), and the remaining 80% are taken for testing (i.e., masked and thus must be estimated). This means that the cardinality of the observed set D equals 20%N 2 . Since the number of known interactions (i.e., 1 entries) are much less than unknown interactions (i.e., 0 entries) in matrix Y , we made sure that the training data contains at most 60% samples of known interactions for each drug, and the remaining are the samples from unknown interactions. The remaining known and unknown interactions for each drug are considered in the test data. The training and test data samples from each of the interaction categories per drug are selected randomly. Also, only one pair of interactions are kept that belong either from the upper triangle or the lower triangle of the DDI matrix. The test data comprises of 3,32,716 samples or interactions. Out of this, 3,18,714 samples correspond to class 0 (i.e., no observed drug interaction) and 14,002 correspond to class 1 (i.e., observed drug interaction) that makes the problem very challenging.

B. GRPMF and benchmarks settings

In all our experiments, GRPMF is applied with the iteration numbers K = 10 and L = 5, that appeared enough to reach stability for the inner and outer loops of our algorithm. As a post-processing, the interaction matrix R recovered using GRPMF is processed by removing the diagonal entries that indicate self interactions (interaction of the drug with itself), and are not of interest in this application. The class labels are assigned using a simple thresholding operation, such that if The performance of the proposed method is compared against the following benchmark algorithms for DDI prediction:

|R ij | ≤ 0.
• Graph DDI [START_REF] Celebi | Evaluation of knowledge graph embedding approaches for drug-drug interaction prediction in realistic settings[END_REF]: In this work, the KGs are constructed using different types of embeddings namely, RDF2Vec, TransE and TransD. These embeddings, along with the DDI matrix, are fed one by one to different machine learning techniques like Random Decision Forest (RF), Gaussian Naive Bayes (GNB), and Logistic Regression (LogR) for DDI prediction. Here, only the RDF2Vec embedding vectors with Skip-Gram is retained for comparison as it achieved the best performance in the study [START_REF] Celebi | Evaluation of knowledge graph embedding approaches for drug-drug interaction prediction in realistic settings[END_REF].

• KGNN [START_REF] Lin | Kgnn: Knowledge graph neural network for drug-drug interaction prediction[END_REF]: Here, the KG results from the drugs' topological structures for DDI prediction. The KG and DDI matrix is fed to Graph Neural Network (GNN). This method focuses on drug neighborhood sampling and aggregates the entities to represent the drugs' potential neighbors in three different ways: (i) Sum, (ii) Concat, and (iii) Neighbor. • Conv-LSTM [START_REF] Karim | Drug-drug interaction prediction based on knowledge graph embeddings and convolutional-lstm network[END_REF]: In this method, KG is learned using different embeddings (e.g., ComplEx, TransE, RDF2Vec, etc.) and fed as input to the CNN and LSTM along with the input DDI matrix, to perform DDI prediction. Only ComplEx embedding is considered here for comparison as it gave the best results in [START_REF] Karim | Drug-drug interaction prediction based on knowledge graph embeddings and convolutional-lstm network[END_REF]. • Additionally, comparisons with Graph Regularized Matrix Factorization GRMF [START_REF] Ezzat | Drug-target interaction prediction with graph regularized matrix factorization[END_REF] and PMFG (7) techniques are presented. The former uses the SIMCOMP drug similarity as an expert knowledge term but, in contrast with our method, it does not infer any graph (i.e., the graph is imposed from the beginning by the user). The method PMFG has been described in the beginning of the paper. Its results are generated using [START_REF] Strahl | Scalable probabilistic matrix factorization with graph-based priors[END_REF] with λ U = 0, using the same post-processing procedure as in our method.

The hyperparameters for GRMF, PMFG and the proposed GRPMF method are tuned using grid search. The other techniques make use of the hyperparameter values mentioned in their respective works. Two main performance metrics are used to evaluate the models, namely the Area under the ROC curve (AUC) and the Area under the PR curve (AUPR). In addition, weighted Precision, weighted Recall, weighted F1 score and Accuracy are calculated for the prediction results.

C. Comparison with benchmarks

Table I summarizes the average values of the performance metrics obtained with the benchmark techniques and the proposed method on the test data obtained using 5 randomly generated train-test splits. Since we have considered ≈ 20% data for training and the rest for testing, the AUPR of all the models, especially the benchmark models employing deep learning and machine learning models are low. It can be seen from the table that among the learning based models, KGNN (Concat) performs best in terms of AUPR and AUC. The remaining learning based methods have AUC ≈ 0.5, indicating no discrimination capability between the classes. They may require more training data for improved performance. The methods based on matrix factorization work much better in this case. Among them, the proposed GRPMF method achieves the best performance across all metrics, with a considerable increase in the AUPR and AUC metrics, respectively, over KGNN (Concat). Additionally, it provides ≈ 28% and 40% improvement in AUPR, and ≈ 5% and 7% improvement in AUC, with respect to GRMF and PMFG respectively.

The high AUPR score of GRPMF indicates that the extra terms in the formulation (10) compared to PMFG [START_REF] Duke | Literature based drug interaction prediction with clinical assessment using electronic medical records: novel myopathy associated drug interactions[END_REF] are able to learn the drug-drug interactions in an effective manner.

Due to the important class imbalance, the values of accuracy and the weighted metrics, namely, Precision, Recall, and F1 score (given in Table I), are more biased towards the majority class (i.e., class 0 in this case). This is the reason for the difference between the AUPR and AUC metrics compared to all other metrics. Confusion matrices for one test set of all the methods are presented in Fig. 1 to provide additional insights on the prediction results. It can be observed that,

D. Ablation Study

Three main hyperparameters govern the performance of the proposed GRPMF method: (i) latent dimension Z, (ii) expert prior weight λ U , and (iii) data fidelity weight λ R . We now study the effect of different values of these hyperparameters on the performance of the proposed GRPMF by varying one of the hyperparameters while keeping the others fixed. Figs. 2 (a) and (b) present the plots of all the performance metrics for different values of the Z. We observe that AUPR metric is high for Z = 20 while the other metrics are comparable for other values of Z. Similarly, from Figs. 2 (c) and (d), we observe that λ U = 0.5 seems to be an optimal choice across all metrics. Figs. 2 (e) and (f) present the performance variation with the change in λ R . Since λ R controls the equality constraint X = U U ⊤ , and as a consequence, the fidelity to the observed data Y , the performance is observed to increase with λ R , with best value λ R = 1. The best values for these hyperparameters were considered for generating the results for GRPMF method in Table I. Finally, Fig. 3 displays the plot of the evolution of the GRPMF loss function F along the iterations for the proposed GRPMF method with the retained settings. It can be seen that F decreases monotonically as expected, and that it reaches stability within a few iterations.

VI. CONCLUSION

This paper presents GRPMF, a novel matrix factorization framework for DDI prediction. Through the introduction of an original regularization strategy, our approach encodes efficiently prior expert knowledge so as to perform jointly DDI prediction and graph similarity inference. Our experimental results obtained with the DrugBank dataset demonstrates the superior performance of our approach compared to deep learning and machine learning-based methods. Although the expert knowledge prior in our experiments relied on the SIMCOMP score for drug similarity, other drug similarity features could be considered in future work. Another extension would be to incorporate an explicit constraint favoring binary entries in the sought matrix, which appeared to be beneficial in DDI, and more generally bioinformatics studies [START_REF] Davenport | 1-bit matrix completion[END_REF]. A possible avenue would be to peruse our recent formulation [START_REF] Talwar | Binary matrix completion on graphs: Application to collaborative filtering[END_REF], by introducing probabilistic graph modeling terms into it.

7 :

 7 Update U k+1 using (14); 8:

Fig. 1 :

 1 Fig. 1: Confusion matrices for all compared methods.

Fig. 2 :

 2 Fig. 2: Performance of GRPMF with different values of Z, λ U and λ R

Fig. 3 :

 3 Fig. 3: Evolution of F along iterations, for GRPMF method.

TABLE I :

 I Quality metrics for all compared methods.

	Method	AUPR	AUC	Precision Recall	F1	Accuracy
	Graph DDI (GNB)	0.0395 0.5	0.9349	0.9559 0.9406 0.9559
	Graph DDI (LogR)	0.0385 0.5008 0.9176	0.9579 0.9373 0.9579
	Graph DDI (RF)	0.0385 0.5004 0.9176	0.9579 0.9373 0.9579
	KGNN (Sum)	0.1830 0.8255 0.9482	0.8867 0.9117 0.8867
	KGNN (Concat)	0.1894 0.8954 0.9502	0.8995 0.9200 0.8995
	KGNN (Neighbor)	0.1123 0.758	0.9344	0.8086 0.8613 0.8086
	Conv-LSTM	0.0402 0.5006 0.9226	0.8932 0.9073 0.8932
	GRMF

  [START_REF] Shi | Predicting comprehensive drug-drug interactions for new drugs via triple matrix factorization[END_REF], decision is class 0, otherwise this is class 1.

			1															1
			0.9							0.962					F1 score Accuracy		0.9
			0.8							0.9615							0.8
		Metrics	0.6 0.7					AUPR AUC	Metrics	0.9605 0.961 0.96						Metrics	0.6 0.7	AUPR AUC
			0.5							0.9595							0.5
			0.4	0	50		100		150	0.959	0	50		100		150	0.4	0	0.2	0.4	0.6	0.8	1
						Z								Z					U
					(a) AUPR, AUC vs. Z					(b) F1 score, Accuracy vs. Z			(c) AUPR, AUC vs. λU
											1							0.97
		0.9615 0.962					F1 score Accuracy			0.8							0.96
	Metrics	0.9605 0.961							Metrics	0.4 0.6						Metrics	0.95
			0.96															0.94
		0.9595								0.2					AUPR AUC			F1 score Accuracy
		0.959								0							0.93
				0	0.2	0.4	0.6	0.8	1			0	0.2	0.4	0.6	0.8	1		0	0.2	0.4	0.6	0.8	1
						U								R					R
				(d) F1 score, Accuracy vs. λU					(e) AUPR, AUC vs. λR			(f) F1 score, Accuracy vs. λR

In our experimental part, due to specificity of the retained dataset, R is a binary valued matrix. However, our framework holds for any type of realvalued symmetric DDI matrix R.