The problem of matrix completion has been deployed in numerous applications of computer science and bioinformatics. One can mention, for instance, recommender systems [13] , wireless sensor networks [12] and MIMO channel estimation [73] . Matrix completion also forms the foundation of many biological interaction problems including gene expression imputation [34] , drug-target interaction prediction [24, 16, 85] , miRNA disease association [11, 40] and gene disease association [54] . As seen in the aforementioned studies, matrix completion is turning out to one of the most promising approaches for modelling biological interactions.
Matrix completion for in silico drug re-positioning is an established topic [24, 84] . In such context, the matrix rows relates to drugs while the columns correspond to the diseases. The information about a drug being effective for a disease is encoded with a one into the matrix and zero, otherwise. The matrix is only partially filled since information about all drug disease associations are not known. Matrix completion amounts to filling the missing entries by relying and interpolating over the available knowledge.
Our recent work [52] introduced a dataset of antiviral drug virus association (DVA). It shows how the matrix completion framework can be used to computationally predict the drugs that could be effective against SARS-CoV-2 (severe acute respiratory syndrome coronavirus 2), the virus responsible for the ongoing pandemic, COVID-19 (COrona VIrus Disease-2019). In the aforesaid dataset, we manually curated drug-virus associations from drugbank.ca. We collected all the anti-viral drugs proved to be effective against viruses infecting humans along with the similarity information between the viruses (genomic structure-based similarity) and drugs (chemical structurebased similarity). The dataset contains both DNA and RNA viruses known to infect human beings. State-of-the-art off-the-shelf matrix completion techniques were then employed to solve the association prediction problem. It was found that graph regularised techniques, including similarities based on genomic structure of viruses and chemical structure of antiviral drugs, performed better than the non-regularised competitors. The most successful techniques were graph regularized versions of matrix completion [51] , binary matrix completion [48] and (shallow) matrix factorization [25] .
In this present work, we propose to solve the problem of drug-virus association prediction via graph regularized deep matrix factorization. We propose a novel theoretically sounded algorithm for this task. Moreover, since our final goal is to repurpose drugs for COVID-19, an RNA virus, we have pruned our DVA dataset to consists only RNA viruses. DNA viruses are not meaningful in this context [55] , hence we will not be considering them in this work. Furthermore, in addition to the graph based on chemical structure similarity of drugs [52] , one more similarity graph has been built on the mode of action of the drugs. Similarly, an additional similarity graph on viruses is built, based on the human symptoms caused by the virus' infection. The current work uses the symptom based similarities along with the genomic similarities used in the prior work [52] .
A prior study on graph regularised deep matrix factorization [49] was based on the alternating direction method of multipliers (ADMM) approach [76] . A major issue with ADMM is that the convergence guarantees are rather mild in such challenging non convex scenario. Furthermore, ADMM method requires the resolution of costly inner steps as well as the introduction of several intermediary variables that may be detrimental to the practical efficiency of the method. Owing to the aforesaid reasons, we propose a novel resolution scheme called HyPALM (Hybrid Proximal Alternating Linearized Minimization), relying on the recent techniques developed in [14, 7, 1] . To summarize the main contributions of our proposed work are as follows:
• An original graph structure introducing (i) symptoms based viral similarity in the virus-virus graph, (ii) mechanism-of-action based drug similarity in the drug-drug graph, leading to the so-called multi graph regularized deep matrix factorization (GRDMF) model; We call the approch GRDMF throughout this paper.
• A novel optimization method, named HyPALM, for addressing the resulting matrix factorization problem with sounded convergence guarantees on its iterates.
In this subsection, we study the choice of technique to solve the GRDMF problem (10) . In the past work [50, 49] , the proposed formulation has been solved using ADMM (alternating direction method of multipliers) [8, 37] . We introduce here a novel resolution strategy named HyPALM. In contrast with ADMM, the iterates of HyPALM are guaranteed to converge to a stationnary point of the non-convex problem. The HyPALM algorithm steps are depicted in the Methods section of the paper. We randomly drop 10 % of association values in the drug-virus association matrix and report the AUC (Area under the ROC curve) and AUPR (Area under the precision recall curve) on the test entries (averaged over 10 runs of cross validation) after solving 2-layer and 3-layer GRDMF problem with ADMM (first two columns) and HyPALM (last two columns). As can be seen, HyPALM solution gives much better prediction results, for both 2 and 3 layers, and both evaluation metrics . Hence, owing to the better performance and the assessed convergence guarantees, we propose to retain HyPALM in the remaining of the paper for our simulations with the GRDMF model. 
We compare the performance of our method with other recent graph regularized matrix completion techniques such has graph regularized matrix completion (GRMC) [51] , graph regularized binary matrix completion (GRBMC) [52] and graph regularized techniques like graph regularized matrix factorization (GRMF) [25] . The authors of [52] show that GRMF and GRBMC perform well in predicting drug-virus association and hence estimating drugs for SARS-CoV-2. We evaluate the performance of 2-layers and 3-layers graph regularized deep matrix factorization model (GRDMF) with the similarities provided by the method from [52] and by the newly proposed models GRDMF-2L and GRDMF-3L.
To quantitatively measure and compare the performance of the proposed and other benchmark algorithms, we carry out 10-fold cross validation on the drug-virus association dataset in 3 settings and observe the standard metrics such as mean AUC and AUPR. The first setting (CV1) corresponds to randomly hiding 10 % association entries as the test set in the association matrix, while the other two settings CV2 and CV3 correspond to randomly hiding 10 % viruses and drugs in the data as test set respectively. Such a validation technique is a standard to evaluate the accuracy in association/interaction prediction problems [16, 74] . As can be observed from Table 2 , both version of our algorithm (GRDMF-2L and GRDMF-3L) perform much better than the state-of-the-art in terms of AUPR and gives comparable results in terms of AUC. This metric is more important than AUC because it penalizes more highly ranked false positives so that only highly ranked drug-virus pairs in prediction would be recommended for biological or chemical tests. It can be seen that the prediction results do not improve much in the 3-layers version. Hence, we use the 2-layers GRDMF for further experiments, as it presents lower complexity. Let us analyse deeper the results in Table 2 . The dataset contains more drugs (86) than viruses (23) . In CV1 settings, when entries are randomly omitted, without skipping drugs and viruses our method performs the best. There is enough data and hence our method does not overfit. The results improve from one (GRMF) to two and three layers. In the CV2 setting, when viruses are omitted, the size of the dataset reduces. The proposed deeper method overfits and hence performs poorly compared to shallower techniques. CV3 also reduces the size of the dataset by omitting drugs, but since the number of drugs are significantly larger than the number of viruses, the effect is not pronounced on our proposed algorithm and we see an improved prediction.
We then measure the precision and recall at top k (Pre@k and Rec@k, k = 3, k = 5 and 7) drug recommendations for viruses while implementing Leave-one-out-cross validation (LOOCV) by hiding (i.e., leaving out) the association profile of every virus. Performing LOOCV tells us how well is the algorithm doing in terms for predicting drugs for current known viruses. This is the most realistic setting, where all the drugs and all the viruses (barring the test virus) are considered. One can see that our algorithm gives comparable precision values and exhibits the best recall at all values of k. This appears as the most realistic case since, on average, there are 4 to 5 drugs for each virus. Table 3 : Precision@k and Recall@k for association prediction for k=5 and 10 with multiple similarities used for drugs and viruses.
The values for the GRDMF model parameters µ, k 1 , k 2 , k 3 and for the HyPALM parameters α, ϑ have been set by performing cross-validation on the training set (see Supplementary Material for the list of retained values). Moreover, we set the number of iterations K = 10 in HyPALM, which appears sufficient to reach stability of the iterates. The time taken by the proposed approach is comparable to the state-of-the art techniques and has been reported in Table 4 . All experiments have been performed on a system with Intel(R) Core(TM) i7 x64-based processor and 8 GB RAM. We also illustrate in Figure 1 the practical convergence profile of our algorithm HyPALM, on an example. Table 4 : Execution times (in seconds) for association prediction for all techniques.
In order to quantitatively assess the necessity/importance of each of the similarity models at use, we conducted an ablation study where we used various combinations of similarity (graph regularizations) matrices within the model to predict 10 % of the test association matrix, which were randomly hided as test set (cross validation setting CV1).
The results have been reported in Table 5 .
Here, S1 d corresponds to chemical structure based similarity for drugs and S2 d corresponds to the similarity based on mechanism of action of drugs. Similarly, S1 v corresponds to the genomic structure based similarity and S2 v corresponds to symptomatic similarity between the viruses. Let us make a few observations from this study. Firstly, taking into account the new similarities gives the best prediction results. Secondly, when we use the chemical structure based drug similarity (S1 d ) and vary the viral similarities from S1 v to S2 v (columns 1 and 2), we observe that symptomatic similarity (S2 v ) yields better results than genomic similarity (S1 v ) and hence contributes more towards the prediction task. Similar behavior appears when using the genomic similarity (S1 v ) for viruses and varying the kind of drug similarity used (columns 1 and 4 of Table 5 ). In other words, our results assess that the new similarities contribute more towards the drug recommendation and greatly improve the prediction results when taken into account in combination with the previous ones. Lastly, we notice that simply replacing the older similarities with the new ones (column 3 of Table 5 ) improves the results. Hence, this ablation study establishes the benefits of the new similarity models that we introduced.
In this subsection, we predict top-10 drugs using the best performing solution GRDMF-2L and the baseline algorithms for the novel Coronavirus (SARS-CoV-2). We display the results in Table 6 Laninamivir Taribavirin  Simeprevir  Taribavirin  10 Dolutegravir Chloroquine Taribavirin Paritaprevir as c Table 6 : Top-10 ranked recommendations/drugs predicted for SARS-Cov-2 by the DVA computational methods.
As can be observed, Table 6 reports the drug recommendations for the novel Corona virus (nCov/SARS-CoV-2) using the proposed algorithm and other baselines. The results from our proposed algorithm are sensible. Ribavarin is an FDA approved antiviral being considered as a treatment for COVID19 either by itself [35] or as a part of combination [30] . Remdesivir is a promising candidate for treating COVID19 having shown efficacy against SARS-CoV-1. The antiviral has been approved for emergency use by the US FDA [20, 60] ; it has been approved by the drug controller general of India (DGCI) to treat mild to moderate COVID-19 patients 1 1 . The science behind the mechanism of Remdesivir's action for inhibiting COVID19 has been thoroughly studied [82] . [80] showed that Umifenovir (brand name Arbidol) is effective as a prophyalactic against COVID19 infection. It has also found to be effective in the treatment of COVID affected subjects [75, 75, 72, 79] . DGCI has approved it for phase 3 trials in India 2 . Favipiravir, a Russian drug for influenza, has also been undergoing trials in Japan; it has been sent to 40+ countries for COVID19 trials 3 . [63] explains the scientific reason behind the prospective success of this antiviral. Baloxavir Marboxil has undergone trial in China for treating SARS-COV2 4 . Interferon alfa-2a in combination with Umifenovir have been found to be effective against COVID19 infected pneumonia [79] . In silico docking studies shows the potential of Baloxavir marboxil [59] ; [42] suggests that this antiviral can be effective if administered within 24 hours of the onset of infection. Several in silico docking studies have shown Dolutegravir may be effective in treating COVID19 infection [4, 39] . Trials for Chloroquine and alternately hydroxychloroquine was once suspended by WHO based on a controversial article, but they were resumed within a few days of suspension 5 . This antimalarial drug has been found to be effective against corona in multiple studies [19] ; the scientific reason behind the efficacy of this drug against COVID19 has been explained in [28] . Geldanamycin and Laninamivir, predicted by our algorithm has not been considered for treating COVID19 patients.
The benchmark techniques GRMC, GRBMC and GRMF have three drugs common to ours -Ribavarin, Remdesivir and Umifenovir. These three techniques also predict an anti-hepatitis C drug -Sofosbuvir; this drug was known to be effective against the coronavirus family [22] ; it's mechanism of action against COVID19 has also been studied [21] . Evidence from clinical studies have corroborated the same [31] . GRMC and GRBMC also predicts chloroquine (as do we); but GRMF is unable to predict it.
Note that the results shown in Table 6 are ranked. The ranking is important, since it gives us a faith in the efficacy of the drugs. Our method predicts 6 drugs under trial among the top 7. The benchmarks only recommend 3 drugs in correct order. In fact, had we looked at the top 5 recommendations, only our algorithm would have recommended all of them correctly.
The benchmark algorithms predict three common drugs -Telaprevir [62] , Boceprevir [26] and Simeprevir [44, 70] . These drugs have the potential to be considered for COVID19 clinical trials. All the three benchmarks predict the antibody Palivizumab; there is no evidence in medical literature to support its usage for treating COVID19 patients or using it as a prophylactic. They also predict Taribavirin, which is similar to Ribavirin, and is used for treating Hepatitis-C; it has not been considered for clinical trials for SARS-COV2.
GRMF predicts a unique drug -Paritaprevir. An in silico docking study [62] showed that it has the potential for treating COVID19 much like Telaprevir, Boceprevir and Simeprevir.
The summary of the findings are in the following Only our proposed algorithm predicts 6 drugs that are under trial (Remdesivir, Ribavarin, Umifenovir, Chloroquine, Favipiravir and Interferon alfa-2a). Two of the benchmarks predict 5 drugs that are under trial (Remdesivir, Ribavarin, Umifenovir, Sofosbuvir and Chloroquine); one of them (GRMF) misses out on Chloroquine.
It is encouraging to observe that all the techniques (proposed and benchmarks) are predicting common drugs -Remdesivir, Ribavarin, Umifenovir.
AI researchers all across the world are contributing towards fighting COVID19. This is our humble contribution towards the same goal [45, 29] . In this paper, we have proposed a novel optimization approach to solve the drug-virus association prediction problem. The contribution of this work is two-fold. Biologically, we have extended the phenomena of solving drug-virus association prediction for drug re-positioning using multiple similarities (graphs) and introduced new similarity measures for both drugs (based on mechanism of action) and viruses (based on symptoms). Technically, the proposed novel technique brings together the benefits of deep learning (using deep matrix factorization) and multiple graph learning (by regularization with multiple graphs). The algorithmic solution to the problem, relying on the novel HyPALM approach, does not require the setting of any extra hyper-parameter and has sounded convergence guarantees. We have made the software publicly available at https://github.com/aanchalMongia/DeepVir.
A significant amount of research has been done in formulating drug re-positioning as a machine learning task. Various approaches such as network diffusion, supervised classification, neighborhood based prediction, clustering and matrix completion have been used in the past to address the said problem. Of these matrix completion techniques have turned out to be the most promising approach [24] . Various biological applications which model the problem of drug-re-positioning using matrix completion include:
• Drug target interaction prediction: The problem is to predict interaction score between a drug and target protein in a matrix assumed to have partial interaction data with drugs on rows and proteins on columns [43, 25, 47, 74, 51, 85].
• Drug Disease association prediction: One predicts the probability that a certain drug will interact with a disease or not in a partially filled drug-disease association matrix. [84, 46, 81, 48] .
• Drug-Drug interaction prediction: The drug-drug relationship is detected through a symmetrical drugdrug matrix/network [83] . This helps predict drugs similar to the ones known to be effective against a pathogen/disease. Attempts have also been made to propose a machine learning approach for predicting drug-likeness of a molecule to distinguish potential drugs from small molecules that lack drug-like features [5] .
In addition to drug re-positioning, other biological interaction problems have also be modeled as matrix completion with the common goal of discovering or improving treatment of human diseases, for example we can cite:
• Gene disease association prediction: Disease and gene features are used to learn genes for diseases as a gene-disease association prediction [54, 78, 36] . This assists deciphering genetic basis of human diseases and better understanding of gene functions, interactions, and pathways.
• Protein-protein interaction prediction: Interactions are estimated in a protein-protein network [77, 38] .
Kshirsagar et al. [38] deploy a matrix completion variant to model interactions between host (human here) proteins and pathogen (viruses causing infectious diseases here) proteins. This helps identifying the interaction between viral proteins and the human proteins, enabling deeper understanding of infectious diseases (which may involve biologically similar pathogens).
Motivated by such works, the problem of Drug virus association prediction was considered in [52] . Although solvable by any computational approach, the authors propose to tackle it using matrix completion motivated by the success of matrix completion in the above-mentioned applications.
Among all the methodologies compared in [52] , graph regularized matrix factorization based technique (GRMF) provided the best results for the validation setting where drugs are predicted for novel viruses. This impels us to come up with an improved algorithm built upon GRMF and leverage the benefit of deep learning and various kinds of metadata associated with drugs/viruses into the current method to aid drug re-positioning for viruses (including SARS-CoV-2).
The goal of matrix completion is to recover all the entries of a matrix given a subset of known entries. Assume that Y is a masked / undersampled version of the complete matrix X. Then, the problem can be expressed as:
Here, M denotes the mask which is element-wise multiplied to X and has ones at positions where the values are known and zeros elsewhere. In general, the problem is under-determined with infinitely many solutions. However, when the sought matrix X is of low-rank, several efficient resolution methods exist. The simplest of them is the matrix factorization approach [15, 25] . Although known to work empirically, a theoretical understanding of factorization based completion is relatively recent [68] . In this approach to matrix completion, X is factored into a product of two low-rank matrices U and V . Formally, this is expressed as:
The factor matrices (U and V ) are obtained by minimizing the following objective function:
with · F the Frobenius norm. The setting of the size of the latent factors U and V allow to impose rank constraints on their product. Motivated by the success of deep learning in numerous fields [64, 67, 61, 53] , the shallow (two factor) models have been extended to deeper versions. The general factorization problem when the matrix is completely observed have been proposed in [41, 71] . The solution to matrix completion via deep factorization has been proposed very recently, in [53] . This can be formally expressed as follows, in the 2-layers case:
Hereagain, the size of the latent factors U 1 , U 2 and V control the rank of the resulting product. One can extent to deeper layer formulations (3-layers and 4-layers) by factorizing X into 4 and 5 matrices respectively. The first factor matrix U 1 (of size m × k 1 ) contains the representation of the row elements while the last factor matrix V (of size k 2 × n) is associated to its column representations in the latent space. The values for k 1 and k 2 here represent the number of latent factors, and are typically set using cross-validation.
Matrix factorization yields non-convex optimization problems which may be difficult to solve. A more direct approach to low rank matrix completion is to directly estimate matrix X under low-rank constraints. Since the minimization of the rank leads to NP hard problems, one usually replaces it by its convex proxy, the nuclear norm [9, 10] , defined as the sum of singular values. The problem formulation is then given as follows:
In regular matrix completion, the entries in the unknown matrix are assumed to be real-valued. In communications and signal processing, this is actually never the case as those values are almost always quantized. This led to the quantized matrix completion problem [23, 6] . One extreme case is binary matrix completion [17] where the matrix entries are represented by only one bit. Binary matrix completion is particularly appropriate for modeling biological interactions, so unsurprisingly studies in the past have used it for microbe disease association [65] .
In all the aforesaid approaches, there was no scope for incorporating associated metadata, to constrain the completion problem and reduce the indeterminacy. In problems of biological imputations, there are often several sources of associated metadata. For example, in drug target interaction, one can compute the SIMCOMP score between the drugs based on their chemical structure as a drug-drug similarity measure [27] . Similarly, similarities between target protein can be found a priori by computing the Smith-waterman score between the amino acid sequences of the target proteins [66] . Studies have shown that the incorporation of associated metadata can improve the prediction results considerably [18] . There can be several approaches to model the associated metadata, in the framework of matrix completion. One of the most promising approaches is to model the associations as graphs. This led to weighted graph regularised matrix factorization [25] . The formulation was expressed as:
Hereabove, µ ≥ 0 is the parameter penalizing the graph regularization terms, and tr(·) denotes the trace of a matrix. Furthermore, the cost function in (10) incorporates Laplacian weight matrices L r and L c imposing prior correlation between row entities and column entities, respectively. Those matrices are obtained from known similarity factors S r and S c , according to
D r and D c are diagonal degree matrices for the corresponding similarity matrices S r or S c , that is, their i-th diagonal entry equals D ii r = j S ij r and D ii c = j S ij c . We recently proposed a deep formulation for (10) in [50] which amounts to solving:
Graph regularised version of nuclear norm minimization was also proposed in [33] and was employed for drug target interactions [74] . The formulation is given by:
All the aforementioned graph regularised techniques have been compared in our prior work [52] in the context of drug-virus association prediction.
In this work, we assume the unknown matrix X to be a complete drug-virus association matrix, with drugs along the rows (m drugs) and viruses along the columns (n viruses). It is a binary matrix with each entry x ij either equal to zero, representing no proven association between the corresponding drug and virus or to one, denoting that the i th drug is known to be effective against the j th virus. We set Y as the observed matrix, corresponding to a sampled version from X. The information of the known association positions is encoded in the masking binary operator M . In our proposed technique, multi-graph regularization is incorporated in the deep matrix factorization formulation with the aim to incorporate the metadata associated with the drugs and viruses in the form of similarity information as shown below:
Hereabove,
with n d and n v the number of available similarity matrices (S ℓ d ) ℓ and (S ℓ v ) ℓ , for row entities (i.e. drugs) and column entities (i.e., viruses). The sizes of matrices U 1 , U 2 and V are m × k 1 , k 1 × k 2 and k 2 × n respectively, where k 1 ≥ 1 and k 2 ≥ 1 are the numbers of latent factors assumed to be involved. Note that our formulation relies on the observation that computing the Laplacians associated to the sums n d ℓ=1 S ℓ d and nv ℓ=1 S ℓ v is actually equals to the sum of all the individual Laplacians involved [51] , given by (11) .
In our previous work on drug-virus association prediction [52] , we only consider one kind of similarity for drugs (based on the chemical structure of the drugs) and viruses (based on the genomic sequence of the viruses), so that n d and n v were equal to one. In this work, we propose to integrate a second type of similarities for drugs and for viruses by taking into account the mechanism of action of drugs and symptoms of viruses, so that n d = n v = 2.
The sought drug-virus association matrix X to be recovered is assumed to be a positive entry matrix with a low-rank structure as similar drugs are known to affect biological systems (target pathways) in a similar fashion by having a similar mechanism of action [32] , which motivates its formulation as the product of rectangular matrices, X = U 1 U 2 V and the addition of the positivity constraint on it.
We deploy the sparsification method to process the similarity matrices (as done in [51, 25] ) using p-nearest neighbor graph which is obtained by taking into account only the similarity values which correspond to the nearest neighbors for each drug/virus. p here is set by performing cross-validation on the training set.
The minimization problem (10) is non-convex. The multi-linear structure of the involved operators motivates us for using the alternating minimization algorithm from [14] . This algorithm has been shown to be very powerful, in several applications such as seismic signal recovery [58] , phase retrieval [56] and hyperspectral imaging [57] . It can be seen as an alternating Majorization-Minimization (MM) technique [69] , including the powerful proximity operator [3] which allows to stabilize and speed up the convergence. In order to apply this method, one needs to introduce explicitly the variable X in the problem formulation, as follows:
with
The advantage is to remove the complicated nonlinear constraint U 1 U 2 V ≥ 0 while not deteriorating the solution (for ϑ > 0 sufficiently large). Function F is then the sum of convex terms involving only a single parameter, and a quadratic, thus differentiable, non convex coupling term weighted by ϑ. We are thus in the framework depicted in [14] (see also [7] for a particular case with two variables), so that we can apply the block coordinate variable metric forward-backward to solve Problem (12) . This amounts to updating the four variables X, U 1 , U 2 and V sequentially, combining gradient (possibly with preconditioning) and proximal updates (see [3] for more details about the proximity operator). As emphasized in [1] , the two first quadratic terms can actually be processed either through proximal or gradient descent updates, without altering the convergence of the method. This nice flexibility allows us to propose the following alternating hybrid scheme, that we called HyPALM for the resolution of (12):
• (update of X) Here, we perform a gradient step over the first quadratic term of (12), followed by a proximal step over the remaining terms:
Hereabove, ς > 0 is a parameter of the algorithm, and α is a stepsize that must be chosen within the interval ]0, 2[ to secure the convergence of the method. The maximum operation in (16) is taken entrywise, so it is equivalent to a simple cropping of the negative values.
• (update of U 1 ) We solve the following minimization problem, equivalent to evaluate the proximity operator of F restricted to the variable U 1 , at the previous version of U 1 denoted U 1 :
By taking the gradient of the above expression with respect to U 1 equals to zero, we obtain:
with I the identity matrix. System (19) takes the form of a Sylvester equation for which efficient solvers are available.
• (update of U 2 ) Similarly to the latter update, we compute:
with U 2 the previous version of U 2 . This leads to the resolution of the Sylvester equation:
• (update of V ) Hereagain, we solve the minimization problem:
with V the previous version of V . This leads to solving the Sylvester equation:
Iterating over the above four steps, with the simplified setting ς = 1, leads to Algorithm 1, where we denote sylv{A, B, C} the solution of the Sylvester equation AX + XB = C. The convergence of the iterates to a stationary point of Problem (12) is guaranteed as a consequence of the theoretical analysis in [14, 1] . It should be noted that the solution for deeper layers (3-layers, 4-layers, etc) of HyPALM method for the GRDMF problem can easily be obtained in a similar fashion to the 2-layers framework shown above.
V (k−1) (U (k−1) 2 V (k−1) ) ⊤ , ϑX (k) (U (k−1) 2 V (k−1) ) ⊤ + U (k−1) 1 } 8: 9: U (k) 2 = sylv{((U (k) 1 ) ⊤ U (k) 1 ) −1 , ϑV (k−1) (V (k−1) ) ⊤ , ((U (k) 1 ) ⊤ U (k) 1 ) −1 U (k−1) 2 + ϑ((U (k) 1 ) ⊤ U (k) 1 ) −1 (U (k) 1 ) ⊤ X (k) (V (k−1) ) ⊤ }
In this subsection, we describe in detail how each of the new similarities between drugs and viruses is obtained.
For drugs, we relied on the chemical structure-based SIMCOMP scores [27] to represent drug-drug similarity as in [52] . We also integrate an additional type of similarity by finding cosine similarities between one hot encoded representation of drug class based on mechanism of action of the drug. Similarly, we used the genomic structurebased d2 * distance based on ONF (Oligonucleotide frequency) measure [2] as similarities between the viruses. To obtain the second similarity measure, we calculate the cosine similarity between one hot encoded representations of symptomatic profiles of viruses i.e. symptoms caused by the virus. Both the drug class information and the file encoding which symptoms are caused by the virus (metadata used to generate new similarity measures) are available as supplementary. Hence, the drug similarity matrix (of size 86 × 86) and virus similarity matrix (of size 23 × 23) are fixed and encode the metadata available.