Brooks Paige

James Bell

Aurélien Bellet

Adrià Gascón

Daphne Ezer
email: daphne.ezer@york.ac.uk

Reconstructing Genotypes in Private Genomic Databases from Genetic Risk Scores

Keywords: Genomic privacy, Genetic risk scores, GWAS

Some organisations like 23andMe and the UK Biobank have large genomic databases that they re-use for multiple different genomewide association studies (GWAS). Even research studies that compile smaller genomic databases often utilise these databases to investigate many related traits. It is common for the study to report a genetic risk score (GRS) model for each trait within the publication. Here we show that under some circumstances, these GRS models can be used to recover the genetic variants of individuals in these genomic databases-a reconstruction attack. In particular, if two GRS models are trained using a largely overlapping set of participants, then it is often possible to determine the genotype for each of the individuals who were used to train one GRS model, but not the other. We demonstrate this theoretically and experimentally by analysing the Cornell Dog Genome database. The accuracy of our reconstruction attack depends on how accurately we can estimate the rate of co-occurrence of pairs of SNPs within the private database, so if this aggregate information is ever released, it would drastically reduce the security of a private genomic database. Caution should be applied when using the same database for multiple analysis, especially when a small number of individuals are included or excluded from one part of the study.

Introduction

In a survey of genomic privacy experts, the long-term privacy of genomic information was deemed both the most important and the most challenging problem to overcome [START_REF] Mittos | Systematizing genome privacy research: A privacy-enhancing technologies perspective[END_REF]. If an individual's password or ID number gets leaked, then it is always possible to change it. However, it is impossible for a person to change their genetic code and they will pass part of it onto their children, so any information leaks can have long-term impacts on both the individual and their descendants. While much of the research focus on long-term privacy of genomic databases rests on the longevity of the encryption scheme [START_REF] Huang | Genoguard: Protecting genomic data against brute-force attacks[END_REF], it is also important to remember that these genomic databases are not just sitting on a server somewhere, but are being continually utilised for making new scientific discoveries. Each time these databases are accessed and the scientific results are published, there is a risk that information will be leaked and that eventually this would enable an attacker to reconstruct private information held in the database.

Genomic researchers are already aware that some forms of aggregate data from their databases should not be released publicly, because there is a risk that an attacker may be able to determine whether a particular individual is a member of the database (a membership inference attack). For instance, such attacks have already been developed for summary statistics about the frequency of single nucleotide polymorphisms (SNPs) [START_REF] Cai | Deterministic identification of specific individuals from GWAS results[END_REF][START_REF] Dwork | Robust traceability from amounts[END_REF][START_REF] Simmons | One size doesn't fit all: Measuring individual privacy in aggregate genomic data[END_REF]. Membership inference attacks have also been developed for the case where a person is allowed to repeatedly query a database to learn if at least one individual contains a particular SNP [START_REF] Raisaro | Addressing beacon reidentification attacks: quantification and mitigation of privacy risks[END_REF][START_REF] Shringarpure | Privacy risks from genomic datasharing beacons[END_REF][START_REF] Von Thenen | Re-identification of individuals in genomic data-sharing beacons via allele inference[END_REF]. These kinds of aggregate statistics about the frequency or presence/absence of a particular SNP might be useful to release to the broader research community, but it is not an essential output of the research process. However, the main research findings -i.e. the SNPs associated with the trait of interest and their strength of association -are essential to publish since the entire purpose of these genomic research projects is to uncover the relationship between genetic variants and phenotypic traits. Moreover, knowledge of these SNPs can lead to new diagnosis procedures or new potential drug targets, so their release is important for the public interest [START_REF] Visscher | 10 years of gwas discovery: Biology, function, and translation[END_REF]. Yet, even this information can potentially leak private information about individuals in the database. For instance, [8] found that information about individuals in a genomic database is leaked when studies publish whether each SNP is correlated or anti-correlated to the trait of interest. It is important to quantify how much information is leaked by publishing these research findings, so that scientists can make informed decisions about when to publish their results and whether it is worth risking the privacy of the participants.

In this manuscript, we demonstrate that the kind of research output that is published from genome-wide association studies (GWAS) has the potential to leak enough information to recover the SNPs of individuals in the database (a reconstruction attack), under specific circumstances. In particular, we focus on the release of Genetic Risk Scores (GRS), a common research output for finding genetic associations with continuous traits [START_REF] Belsky | Development and evaluation of a genetic risk score for obesity[END_REF][START_REF] Chouraki | Evaluation of a Genetic Risk Score to Improve Risk Prediction for Alzheimer's Disease[END_REF]4,[START_REF] Knowles | Cardiovascular disease: The rise of the genetic risk score[END_REF][START_REF] Qi | Genetic risk score and risk of myocardial infarction in hispanics[END_REF][START_REF] Zhao | An obesity genetic risk score is associated with metabolic syndrome in Chinese children[END_REF]. We also focus on cases where a database is repeatedly used to perform a GWAS analysis, but not all the individuals are part of all the analyses. This could be the case because some individuals drop out of the study or skip specific survey questions. Alternatively, some databases, such as 23andMe, may grow in size over time and Fig. 1. We investigate the case where two GWAS studies are performed on two data sets that mostly contain the same individuals. We reconstruct the genotype of those individuals added to the second study, using the GRS from each study and an estimate of SNP frequencies.

( " # ! ) ( " # !"# ) (%) (& # )
allow several GWAS to be performed within a short period of time. Under these circumstances, we demonstrate that it is possible to completely reconstruct the SNPs of an individual using a custom Expectation-Maximisation (EM) algorithm. We also provide suggestions for avoiding this kind of attack.

To be clear, this manuscript focuses on the simpler case in which the exact same trait is investigated in multiple GWAS studies; however, we expect that some version of this attack may be developed in the near future for the case of multiple highly correlated traits.

Overview of scenarios that will be investigated

We demonstrate a series of reconstruction attacks that enable us to infer the genotypes of individuals in private genomic databases, based on publicly released GRS. These attacks will initially be deployed on a very favorable scenario, but the scope of the attack will be subsequently expanded, building up to the scenario shown in Figure 1. It is worth noting that the reconstruction attacks that we will describe do not depend on (i) how the SNPs were initially filtered or (ii) how strongly they associate with the trait of interest.

We will begin by investigating a simple scenario: two GWAS studies are performed to identify SNPs associated with the same trait, and the two studies use the same set of participants, except that the second study includes one extra individual. In addition, we will assume that we know the frequencies of each SNP and the frequencies that pairs of SNPs co-occur in the same individual. We assume that both studies publish the coefficients associated with the GRS models that they infer as part of the analysis.

Next, we will consider the case in which the second study includes more than one additional participant and demonstrate that in many circumstances this still allows us to easily reconstruct the individual genotypes of all the individuals that are found in the second study but not the first (see Section 3.2).

Afterwards, we will demonstrate that we do not need to know the precise frequencies of SNPs and frequencies of co-occurring SNPs, as long as we have a reasonable estimate of these values from public databases (see Section 3.3).

We also briefly discuss how loosening additional restrictions would impact our ability to predict individual genotypes. In particular, we analyse the case where the two sets of SNPs that are used by the two studies are not identical. These results imply that if two sets of GRS are released on two genetic data sets with largely overlapping populations, it may be possible to reconstruct the genotypes of those individuals who participated in one study but not the other (Figure 1).

Methods

Genetic risk score (GRS) models describe the relationship between a particular phenotype of interest and particular SNPs. These models are fit in a two-stage process: first, a reduced set of SNPs is selected from a potentially very large pool of candidates; then, this reduced set is used as the independent variables in a linear regression analysis. The set of SNPs is selected by first filtering for those that significantly correlate to the trait of interest, after controlling for other covariates. These SNPs are then further filtered to ensure that they are far apart from one another, in order to decrease the correlation between them.

In this setting, we suppose M individuals have taken part in a study, and N SNPs have passed the filtration steps to be used in a linear model. Let y M be the vector of M real-valued phenotypes, and X M be an M × N binary matrix, where X M [i, j] = 1 if individual i has SNP j. To include an intercept term in the linear model, we define a design matrix Φ M to be the M × (N + 1) matrix

Φ M = X M 1 M . (1)
The GRS model parameter β M is just the coefficient vector of the linear model

y M = Φ M β M + , (2)
where is independent Gaussian noise. Given Φ M and phenotypes y M , the maximum likelihood estimate of this parameter has a closed form

βM 1 M K -1 (Φ M y M ), (3)
where we have defined the symmetric (N + 1) × (N + 1) matrix K as

K = 1 M Φ M Φ M . (4)
Now, suppose a second study is run, targeting the same phenotype, which adds a single extra individual with SNPs represented by the N length vector x 0 . This corresponds to adding the row φ 0 = [ x 0 1 ] to the design matrix, and extending y with the additional phenotypic value y 0 for the new individual. The updated estimator (i.e. the GRS values for the second study) is given by βM+1

= (Φ M Φ M + φ 0 φ 0 ) -1 (Φ M y M + y 0 φ 0 ). ( 5
)
We assume that both GRS models βM and βM+1 are released publicly. An attacker aims to use this knowledge to reconstruct φ 0 (the genotype of the added individual). Through algebraic re-arrangement we find that:

φ 0 = 1 C K( βM+1 -βM ) ( 6
)
where C is a scalar, specifically C = 1 M (y 0 -φ 0 βM+1 ). This means that φ 0 is a scalar multiple of K( βM+1 -βM ).

Our approach thus centers on the use of the vector we define as d 1 ,

d 1 K( βM+1 -βM ) = Cφ 0 , (7)
corresponding to a rescaled copy of the input SNP data in the design matrix φ 0 , which can be easily computed from the two parameter vectors if the matrix K is known. As we will see in Section 3.1, we can use d 1 to exactly reconstruct the added individual with 100% accuracy. We additionally consider the case where m additional individuals have been included in the second study, yielding a new GRS model βM+m including these M + m participants. The extra rows of the design matrix now form a matrix Φ m of size m × (N + 1), where each row is an individual that was added to the second study and each column is a SNP (and the last column contains only 1). The corresponding analog to Eq. ( 7) for multiple individuals, which we derive in Appendix B, is

d m K( βM+m -βM ) = Φ m C m , (8)
where C m is a vector of length m. For sufficiently small m (relative to N ), exact reconstruction of all m added individual genomes is also possible in this setting, following the algorithm we will introduce in Section 3.2.

The previous examples have focused on cases in which the participants in the first study are a subset of the individuals in the second study. In Appendix C we consider the case in which the first study has some participants that are not found in the second study and vice versa. We show that the same strategies for reconstructing the genome can be used as in the previous scenario that we discussed, in which multiple participants are added to the second study.

Estimation of K

As it turns out, the entries of matrix K correspond to simple population-level statistics of the SNPs, which could either be inadvertently released (under the assumption they would be safe to share), or could be estimated from another sample from the same population. In fact, the entries of K depend only on the SNP frequencies and SNP co-occurrence frequencies in the dataset:

-For i = 1, . . . , N : K ii estimates the probability that SNP i has value 1 (i.e.

the frequency of the SNP in the population).

-For i = 1, . . . , N -1 and j > i: K ij = K ji estimates the probability that SNP i and SNP j are both 1 simultaneously (i.e. the frequency of SNP i and SNP j co-occurring in the same individual). -For i = 1, . . . , N and j = N + 1: K ij = K ji also estimates probability that SNP i has value 1, i.e.

K i,N +1 = K N +1,i = K ii . -Finally, K N +1,N +1 = 1.
Thus, knowledge of SNP frequencies and pairwise co-frequencies from the original study are all that is required in order to compute K. In the following Sections 3.1 and 3.2, we consider adding one and multiple individuals at once, respectively, in the setting where this matrix K can be estimated exactly.

However, while βM , βM+1 and M are likely to be published along with the study, an attacker would often need to estimate K from other publicly available data. Most studies will report some information about the study population (such as whether the study focused on individuals from a specific continent), which can help with estimating K. From this information, we can estimate the value of K in similar populations as those used in the study using publicly available data, e.g. from the HapMap project. Our additional experiments in Section 3.3 use a custom EM algorithm to find maximum likelihood estimates of φ 0 when the matrix K ≈ K is estimated from independent public data. The derivation of this EM algorithm is given in Appendix D.3, and a formal analysis of the reconstruction error of φ 0 given the error in K is found in Appendix D.1.

Results

The key observation from the previous section is that the vectors d 1 and d m , derived from the change in parameter vectors β from a first study to a second study, take only a finite number of values thanks to the fact that the design matrices Φ contain only zeros and ones. In particular, when m new individuals are added to the second study, each entry of the vector d m can only take at most 2 m values, and a zero value corresponds to the setting where all individuals have the most common variant for that SNP.

This section describes algorithmically how these vectors can be used to recover the genomes of the additional individuals, as well as empirical tests which use the Cornell Dog Genome dataset as a case study [6]. More details on the experimental setup can be found in Appendix A.

Complete reconstruction of one individual's genotype when SNP frequency information is known

The first, most straightforward case is when only one participant is added between the first and second studies, i.e. where βM is the GRS for the first study (containing M participants), and βM+1 is the GRS for the second study as described in Eq. ( 3) and ( 5). Both of these are vectors of length N + 1, where the first N indices correspond to the relationship between each SNP and the trait and the last element is the intercept of the linear model. For now, we also assume we are in the setting where the matrix K is known, e.g. because the SNP frequency information has been publicly released.

0 0 0 0 1 0 0 A 0 1 0 B 0 0 1 C 1 1 0 A+B 1 0 1 A+C 0 1 1 B+C 1 1 1 A+B+C !( # $!"# -# $!)[(]

Possible combination of SNPs at index

i A B !( # $!"# -# $!)[(]
Given K, βM+1 and βM , we can use d 1 (a vector of length N + 1) to precisely determine the genotype of the individual who was added to the database. For each i = 1, . . . , N , the i th entry of d 1 is either equal to 0 if φ 0 contains a 0 (i.e. the individual does not have the SNP at that index) or to C if φ 0 contains a 1 (i.e. the individual has the SNP at that index). In other words, it is possible to exactly read off the SNPs of the added individual in this setting. Indeed, we tested this strategy on the Cornell Dog Database and found that we were able to reconstruct the genotype of the dog that was added to the second study with 100% accuracy, both on common and uncommon SNPs (see Figure 2(A)).

Complete reconstruction of multiple individuals' genotype when SNP frequency information is known

We now consider the case where m additional individuals have been included in the second study, yielding a new GRS model βM+m including these M + m participants. Consider again Eq. ( 8) above. The i th row of Φ m is a binary vector that represents the combination of the m individuals who have SNP i. This means that, for a fixed value of C m , the value of the vector d m at index i is uniquely determined by the combination of individuals who have SNP i (Figure 2(B)). In other words, there will be at most 2 m unique values taken by entries of d m , each corresponding to a combination of the values in vector C m (see Figure 2(C)).

If we were to learn which values of d m are also found in C m , then we could infer the complete genotypes of all the m individuals added to the second study. We would be able to reconstruct m complete genotype vectors, although it would be impossible to know which of the genotypes corresponded to which of the m individuals. In fact, in many cases it is extremely straightforward to determine which values in d m correspond to values in C m . Here we describe a simple algorithm for finding C m when there are exactly 2 m unique values in d m . If this is not the case, please see the more complete algorithm in Section G.

1. First, extract all unique, non-zero values from d m . 2. Find the sum of all pairs of values in this vector. [START_REF] Chouraki | Evaluation of a Genetic Risk Score to Improve Risk Prediction for Alzheimer's Disease[END_REF]. Find all values that are in (1), but not in [START_REF] Cai | Deterministic identification of specific individuals from GWAS results[END_REF]. The values of C m appear in this list. There is no way to know which value of C m corresponds to which index, so for simplicity we can randomly assign them indices. 4. Each value in C m corresponds to a specific individual who was added to the second study. Each value in d m can be described as a sum of a unique combination of values in

C m . For instance, if d m [i] = C m [j] + C m [k],
this means that the SNP at position i is found in individual j and k, but no one else.

We tested this approach using the Cornell Dog Database, in a test scenario where the second study added three different dogs. We were able to uniquely identify the genotypes of all three dogs with 100% accuracy, both with common and uncommon SNPs (Figure 2(C)).

Accurate estimation of an individual's genotype when SNP

frequency information is estimated from a public database (right) adding three new participants. These figures are analogous to those in Figure 2, albeit in the case where K is not known and instead estimated from an independent public database.

Previously, we assumed that the attacker had access to the matrix K, which consists of population-level statistics on frequencies and co-occurrence frequencies of SNPs. While this could be released voluntarily by organisations which are not aware of the risk, we now consider the case where K is not directly available to the attacker but is instead estimated from a separate public database assumed to correspond to individuals from the same population. We simulated this scenario using the Cornell Dog Database by taking one random set of dogs for building the GRS model, and a second non-overlapping set of dogs for estimating K. We compared the value of d1 = K( βM+1 -βM ) with the known value of φ 0 . We observe that d1 has significantly different values at indices where φ 0 [i] = 0 and φ 0 [i] = 1; examples for the cases where one and three dogs are added can be seen in Figure 3. Public set: 800 dogs Accuracy Fig. 4. Accuracy at reconstruction of genomes x0 using EM estimation and a noisy estimate K, as compared to a natural baseline which always predicts the most common variant at each SNP locus. We use this as a baseline, because without any additional information about βM and βM+1, the most accurate prediction of the dog's genotype would be to predict the most common variant at each locus. Here we define accuracy as the proportion of SNPs that are correctly identified in the dog that was found in the second GWAS study, but not the first. Each distribution is constructed from 500 experimental test points, in which we (i) took 10 random splits of the full dog data set, assigning dogs to either the public and private data set (ii) for each split, we tested the reconstruction 50 times, each time adding a different randomly sampled dog to the second GWAS study. The private dataset always has 1000 individuals; the public test dataset is of increasing size, improving performance.

The main challenge is that the vector d1 now includes additional noise, so we cannot simply use its entry at index N + 1 to estimate C, nor do the entries i with φ 0 [i] = 0 also correspond directly to d1 [i] = 0. Instead, we develop a custom expectation-maximisation algorithm to find a maximum likelihood estimate of the constant C and recover φ 0 , i.e. to determine the probability that each φ 0 [i] = 0 or φ 0 [i] = 1, based on the value of d1 (see Section D.3 for details). We find that this method can successfully reconstruct the correct value of φ 0 [i] much better than a baseline which uses the public dataset to independently estimate the most common variant for each SNPs (see Figure 4). Crucially, we show that our approach is able to reconstruct, with relatively high accuracy, the genotypes of dogs even when they differ significantly from those in the public dataset (see Figure 5). This shows that our attack is able to extract information about the particular individuals that differ across the two studies, not merely about the general population as in the most-common-variant baseline. By definition, dogs that have genotypes that differ significantly from the general population have a higher proportion of uncommon SNPs, and the ability to recover these uncommon SNPs is particularly important from a privacy perspective. Indeed, uncommon SNPs can be used to identify a particular individual and are also more likely to be associated with disease phenotypes, which is sensitive information. In general, we find that the larger the public dataset available, and the more similar the dataset is to the unknown private dataset, the better we are able to reconstruct the genome of the added individual. Full details and description of the experimental setting are given in Section A. We also derive theoretical error bounds for our estimate of φ 0 based on the error in K in Section D.1. Here each point represents a dog and we define atypicality as the proportion of uncommon variants that the dog has compared to the public database-for instance, if 51% or more of dogs in the public database have a G in a specific locus, but this dog has a T, then this would count towards the dog's atypicality. In other words, dogs further to the right are less and less similar to average dog present in the public dataset (measured by percentage of different variants). In contrast to the most-common-variant baseline, our method generalizes well even to dogs which are highly dissimilar to those in the public dataset. Larger public databases (right) provide more accurate population estimates K, leading to more accurate reconstructions overall.

Accurate estimation of an individuals' genotype when different SNPs are used in each study

When GRS models are constructed, the first step is to filter the set of SNPs down to a small set of SNPs that are (i) significantly correlated to the trait after covariates are considered and (ii) far apart from one another along the genome.

If the two studies use two different sets of SNPs to construct the GRS model, it is still possible to recover whether or not each of the SNPs in the overlap is present in the new individual. This process is highly analogous to the previous cases and is detailed in Appendix F.

Discussion

In this manuscript, we demonstrate that private information is leaked when GRS models are published, specifically in the case where two sets of largely overlapping individuals are used for multiple studies. In particular, we show that we can recover SNPs from an individual in a private database-a reconstruction attack. Even though we would not have a name associated with this genotype, it may be possible to identify the individual once the genotypic data is available to the attacker. For instance, the attacker may have access to partial genotypic information of the individual and then be able to identify them. Alternatively, they could use the genotype information to predict ethnicity and other phenotypic traits that could then be used to uniquely identify the individual. We also note that even an incomplete reconstruction attack (in which only a proportion of the SNPs are correctly identified) is likely to be sufficient to perform a membership inference attack. Investigating the relationship between the reconstruction attack and the membership attack will be a subject of future research. Importantly, if the attackers were unable to link the genomic data with a particular individual, the reconstruction attack would still be a breach in privacy that could have serious consequences. For instance, the patient may have only consented to have their genomic data used in particular kinds of research studies, while the attacker may use the reconstructed genomic data for a different (potentially unethical) purpose.

Suggestions for good practice. We provide a number of simple suggestions for good practice that would help limit this attack.

1. Aggregate statistics about the frequency of SNPs in the database or the frequency of co-occurrence of SNPs should never be released. We have shown that this information, combined with GRS, allows to precisely reconstruct individual genomes in various settings. It may be possible to release noisy versions of SNP frequency data, but this would be equivalent to releasing K (our estimated K from the public database). With our EM algorithm, we have demonstrated that it is still possible to do some genotypic reconstruction with a noisy K, but this becomes harder as the noise in K increases. However, providing a very noisy K may be of limited utility to the scientific community.

2. If a genetic data set is intended to serve for multiple complementary analyses, it is important that all study participants are used in every analysis performed. If there is missing phenotypic data from a few individuals, they should not be included in any of the analyses performed, or their privacy may be compromised. 3. When multiple individuals are added in between two studies, then the ability to reconstruct the genomes depends on the number of SNPs being large relative to the number of individuals. In particular, if m new dogs are added, exact reconstruction is only possible using the approach in Section 3.2 if the number of SNPs N > 2 m . Thus, we suggest to avoid releasing multiple studies which differ by fewer than log 2 N individuals.

Extensions and future work. While we have analyzed the case where the genome is represented by binary values of 0 or 1, often studies instead count the number of times each allele is present, which would lead to a design matrix Φ containing values 0, 1, or 2. In this scenario, K no longer contains the frequencies of SNPs and their co-occurrences, but something slightly more complicated that we describe in Appendix H. This does not dramatically change the approach in this paper, except in that the vector d m can take 3 m possible values, rather than 2 m . In practice, then, studies which use allele counts are somewhat more robust to attacks; the multiple dog reconstruction attack would likely be ambiguous if 3 m > N , rather than 2 m > N .

A possible countermeasure to our reconstruction attack could consist in randomly perturbing the GRS models before releasing them, as done in differentially private linear regression [19]. However, a naive application of this strategy could destroy the utility of the models. A formal and empirical analysis of the effectiveness of such protection against reconstruction attacks, as well as of the usefulness of the resulting GRS models to genomic researchers, is beyond the scope of this paper and left for future work.

Another countermeasure could consist in refraining from releasing precise information about the population structure of the study population to prevent the attacker from estimating K effectively. This would however limit the utility of the research study, because the researchers would not know to what populations the research applies to. An investigation of how differences in population structures impact K est will be undertaken in public human genomic data sets in future work, but is also beyond the scope of this paper.

Our work has a number of limitations. For instance, we only test our EM algorithm on dog data. Dog populations may have different population structures than human populations due to selective breeding, so in the future we aim to test how properties of population structure will impact our ability to estimate K and the accuracy of our reconstruction attack. In addition, the expectationmaximisation algorithm described here is only explicitly described for the case in which one participant is added to the database, and the matrix K is estimated from a public database. We are currently extending the algorithm to the broader case in which multiple participants are added at once.

It may seem on the surface unlikely that two GWAS analyses will include nearly the same participants. One potentially common setting where this could arise is when a single study collects both genotype and phenotype data from a single set of participants, and releases multiple models to predict multiple traits. In this case, there may be a small number of individuals who are used in one analysis, but not the other; for instance, there may be a small subset of participants who skip a particular survey question that was used to collect phenotype information, and this is indeed evident in a recent study ([9]). In such settings, it could be very possible for multiple released GRS models to be computed on sets of individuals which differ by only a few participants. In future work, we aim to extend our analysis and attack to settings where multiple GRS models are released, each predicting different but highly correlated traits.

A Experimental details

Cornell Dog Database: To experimentally test the reconstruction attacks, we used data from the Cornell Dog Genome Database, which contains data about SNPs from a wide range of dog breeds and a number of associated phenotypic traits. The two traits we focused on were average breed weight and average breed height, because these two phenotypes had the fewest number of missing values. For the initial investigation, we binarised the genotype matrix-considering all heterogenous alleles to have a value of 1. (We also repeated the analysis with the original genotype matrix.) Only common SNPs (i.e. SNPs that were found in 25% to 75% of the dogs) were used, leaving 23,497 SNPs. For each linear model built, M = 1000 dogs were randomly sampled as the "private" dataset and N = 200 SNPs were randomly selected. To ensure that the SNPs that were sampled were spatially distributed, the SNPs were randomly sampled in a stratified way, so one SNP was selected in every 23,497 200 -sized bin.

Experiment with imprecise K: First, two linear models were constructed to predict average breed weights: one with the M = 1000 randomly sampled dogs and another that contained 1 additional randomly sampled dog. This gives βM and βM+1 . To mimic the process of estimating K from a public database, we randomly sampled an additional 200, 400, or 800 dogs that were not included as part of the original set and used this to estimate K, which we refer to as K. Now we could calculate K( βM+1 -βM ) and compare this to the known φ 0 for the additional dog from the second study. These additional dogs are taken from a third "test" dataset, disjoint from both the public and private data. The plots in Figures 4 and5 are produced by re-running the algorithms across 10 random public / private / test splits, where the "test" dataset has 50 dogs which are each individually considered as candidates for the (M + 1) th dog added to the private dataset.

B Adding multiple dogs

Here we explain Equations ( 7) and (8). Note that the former is a special case of the latter so we will only explain the latter in detail. First note that by definition βM

= (Φ M Φ M ) -1 Φ M y M = (M K M ) -1 Φ T M y M , βM+m = (Φ M +m Φ M +m ) -1 Φ M +m y M +m = (M K M + mK m ) -1 Φ M +m y M +m .
Substituting these into the left hand side of the following equation gives the right hand side:

(M K M + mK m ) βM+m -M K M βM = Φ T m y m . (9)
This equation can be rearranged to give

K M ( βM+m -βM ) = 1 M Φ T m y m - m M K m βM+m = 1 M Φ T m (y m -Φ m βM+m ).
Defining the length m vector C m = 1 M (y m -Φ m βM+m ) yields the form used in Equation (8). For the special case of m = 1, C m is a scalar and we recover Eq. [START_REF] Huang | Genoguard: Protecting genomic data against brute-force attacks[END_REF].

C Case in which each GWAS study adds two new sets of participants

This manuscript mostly explores the case in which one study's participants are a subset of the other study's participants. Here we demonstrate that this is equivalent to the case where each of the two studies contain a small number of participants that are not found in the other study.

In particular, let us say that the first study has M + a participants and the second study has M + b participants, where the first M participants are shared between the studies, but there are a participants that are found in the first study but not the second, and b participants that are found in the second study but not the first. Following on from Equation 9, we see that:

K M ( βM+a -βM+b ) = K M ( βM+a -βM ) -K M ( βM+b -βM ) = 1 M Φ T a (y a -Φ a βM+a ) -Φ T b (y b -Φ b βM+b ) .
Let us define the following (N +1)×(a+b) matrix obtained by concatenating the two genotype matrices:

Φ a+b = [Φ a , Φ b ] (10)
and the following a + b length vector:

r a+b = (y a -Φ a βM+a ), -(y b -Φ b βM+b ) (11)
Then this gives us:

K M ( βM+a -βM+b ) = 1 M Φ a+b r a+b (12)
This means that having two non-overlapping participant sets is equivalent to the setting in which the first study is a subset of the second (only m is now a + b).

D Estimating K

If the true matrix K is unknown, it can be estimated with public data. We denote this estimator by K. In order for K to be an accurate estimate the data that it is generated from must be drawn from the same (or a sufficiently similar) population as that used in the private study. We will model this assuming no discrepancy between population distributions, however when we discuss how to evaluate whether the estimate is good that assessment should account for this systematic error as well.

In the following we are primarily concerned with the error due to the subsampling in both the private and public data sets. For convenience we only consider the case of adding a single individual, though the generalization is quite straightforward.

D.1 Analytic bound on φ 0 -φ0

If K is substituted for K in our reconstruction equation [START_REF] Huang | Genoguard: Protecting genomic data against brute-force attacks[END_REF] we get an approximation of φ 0 which we denote φ0 . We would like to bound the (relative) error between φ 0 and φ0 . In the following, we ignore the constant factors and Ĉ for simplicity, noting that these scaling factors are estimated from the resulting φ 0 or φ0 anyway. We thus consider ϕ 0 = K( βM+1 -βM ) and φ0 = K( βM+1 -βM ). Using • on vectors, and also on matrices to denote the corresponding operator norm. The relative error between ϕ 0 and φ0 is given by:

ϕ 0 -φ0 ϕ 0 = ( K K-1 -KK -1 )ϕ 0 ϕ 0 ≤ K K-1 -KK -1 = K( K-1 -K -1 ) . Note that K-1 -K -1 = K-1 (K -K)K -1 and hence φ0 -ϕ 0 ϕ 0 ≤ K -1 K -K , (13)
This means we can bound the error by two quantities. The term K -1 is bounded above by 1/ min(eig(K)), which is finite as soon as K is non-singular. This is not a strong requirement as in the case of linear regression it is required for βM+1 and βM to exist. Note that in the case of L2-regularized linear regression (i.e., ridge regression), K is replaced by K + λI where λ is the regularization parameter, and we can directly bound this term by λ.

The key term in ( 13) is K -K , the error in estimating K by K. Let us assume that the public database used to obtain K follows the same distribution as the private database used to fit the GRS models. Denote by M the number of individuals used to estimate K. Then, under classic boundedness assumptions and leveraging matrix concentration inequalities such as matrix Bernstein [START_REF] Tropp | An introduction to matrix concentration inequalities[END_REF] we can show that E[ K -K ] = O(1/ min( M , M )). This shows that the error in estimating K is small as long as the private and public databases are large enough.

D.2 Modelling the error in K

In this section we define a model to capture the error in K, which leads to the expectation maximization algorithm for estimating φ 0 which is used in the experiments.

As our estimated K drifts from the true K, this expression K( βM+1 -βM ) would produce a wider range of values than just 0 and C. Let ij ∼ N (0, σ 2 ) be independent noise, which we assume corrupts each element of K ij ; i.e. given the estimated matrix K, suppose

K ij ∼ N ( Kij , σ 2 ), (14)
for some small σ 2 . This is clearly an oversimplification (as we know K is e.g. bounded and symmetric), but is a useful starting point that allows derivation of a simple estimation algorithm. For notational brevity, in this and the following section we define the vector

∆ = βM+1 -βM (15)
which corresponds to the difference between the two GRS model parameter vectors. Given the true value of K, the system of equations Cφ 0 = K∆ relates the known quantity ∆ and the Gaussian-distributed K with the unknown value of C and of the vector φ 0 . If K is Gaussian (following Eq. ( 14)), then the linear transformation K∆ is Gaussian as well. We denote each of the rows of K as a vector k i , i = 1, . . . , N ; then for each row, the scalar value

k i ∆ ∼ N ( k i ∆, σ 2 ∆ ∆),
meaning overall the vector K∆ is distributed N ( K∆, σ 2 ∆ ∆I). With some algebraic re-arrangement, and since for the true underlying value of K we have K∆ = Cφ 0 , we can write this as

K∆ ∼ N (Cφ 0 , σ 2 ∆ ∆I) (16)
where C and σ are parameters we need to estimate. The vector K∆ is observed "data", computed from the public SNP database and the two released parameter vectors. We can model each of the entries of φ 0 , which are zeros and ones, as Bernoulli distributions, whose prior probabilities correspond to the public dataset estimated frequencies. This suggests a model for K∆ which is akin to a constrained mixture of Gaussians.

D.3 Derivation of EM algorithm for reconstructing an individual's genotype

Given the model in Eq. ( 16), we can define an EM algorithm for maximum likelihood estimation of C and σ 2 , which then permits easy inference for each entry of φ 0 . Denote the entries of φ 0 as z 1 , . . . , z N +1 ; denote their prior probabilities as α 1 , . . . , α N +1 , where α 1 , . . . , α N are the (public) population frequencies for each SNP, and α N +1 = 1.0. Let x 1 , . . . , x M denote the entries of the fixed (observed) vector x = K∆, which is distributed

p(x|C, z, σ 2 ) = N (x|Cz, σ 2 ∆ ∆I).
Supposing we know values of C, σ 2 , to estimate the entries of φ 0 we want to find p(z|x, C, σ 2 ), p(z|x, C, σ 2 ) ∝ p(x|C, z, σ 2 )p(z)

For each z i , i = 1, . . . , N , we define the quantity

π i = p(z i = 1|x, C, σ 2 ) = α i N (x i |C, σ 2 ∆ ∆) α i N (x i |C, σ 2 ∆ ∆) + (1 -α i )N (x i |0, σ 2 ∆ ∆) , (17)
the posterior probability of each particular entry taking a value of 1, rather than 0. A maximum likelihood algorithm to estimate C, σ 2 proceeds by alternately:

1. Given C, σ 2 , estimate the posterior distribution π = p(z|x, C, σ 2 ); 2. Given the posterior π, maximize L = E π [log p(x|C, z, σ 2 )] with respect to C and σ 2 .

To maximize C and σ 2 , we first compute the derivatives of

L = i zi p(z i | . . . ) log p(x i |C, z i , σ 2 ) = N i=1 π i log N (x i |C, σ 2 ∆ ∆) + (1 -π i ) log N (x i |0, σ 2 ∆ ∆).
We take derivative w.r.t. C,

∂L ∂C = i π i σ 2 ∆ ∆ (x i -C)
and set equal to zero to find

Ĉ = i π i x i i π i . (18)
For the noise term σ 2 , we have

∂L ∂σ 2 = N i=1 π i ∂ ∂σ 2 log N (x i |C, σ 2 ∆ ∆) + (1 -π i ) ∂ ∂σ 2 log N (x i |0, σ 2 ∆ ∆)
which with some algebraic rearrangement becomes

σ2 = 1 N ∆ ∆ N i=1 π i (x i -C) 2 + (1 -π i )x 2 i . (19)
These updates taken together can be used to define an EM algorithm which optimizes the values of C and σ 2 , despite the fact that the entries of φ 0 are unknown; once C and σ 2 are then known, the vector π will give probability estimates for each entry of φ 0 . The overall EM algorithm can be summarized by the following iterative updates:

1. π i ≡ p(z i = 1|x, Ĉ, σ2 ) = α i N (x i | Ĉ, σ2 ∆ ∆) α i N (x i | Ĉ, σ2 ∆ ∆) + (1 -α i )N (x i |0, σ2 ∆ ∆) , 2. Ĉ ← i π i x i i π i , 3. σ2 ← 1 N ∆ ∆ N i=1 π i (x i -Ĉ) 2 + (1 -π i )x 2
i . To initialize the algorithm, we can set π i to some initial probabilities, and find initial values for Ĉ, σ2 ; we experimented with both setting to the prior probabilities per-SNP estimated from the public data, as well as to the vector of all zeros (corresponding to a "hard" initialization at the value of the baseline estimate), and found no qualitative difference in performance.

E Scaling of EM algorithm with size of private dataset

Figure 6 demonstrates the change in accuracy of the EM algorithm over a range of different private database sizes. For this test, a synthetic dataset with 100 SNPs and 1,000,000 individuals is generated; 10,000 are held out as a public database, and 30 individuals are taken as a fixed test dataset of new dogs to add and are used to estimate EM algorithm accuracy, across increasingly large private database sizes. The algorithm has stable performance for increasingly large private databases.

F Estimating φ 0 with different SNP sets

Here we analyse what can still be said in the event that the two studies do not use exactly the same set of SNPs. We will still assume the sets of SNPs considered to have a significant overlap.

For this purpose we will need a greater variety of notation. A primed variable denotes that it corresponds to the second set of SNPs, e.g. K is the co-occurrence matrix from the original M users for the second experiment. If a vector or matrix is surrounded by square brackets this denotes that object but with any rows or columns corresponding to SNPs not in the overlap removed, e.g. [K] denotes the co-occurrence matrix from the first experiment restricted to the overlapping SNPs.

As before, from the first experiment, we have

K β = Φ T y (20)
and now, from the second experiment, we have

(K + φ T 0 φ 0 ) β = Φ T y . (21)
Taking the difference between these expressions, as before, gives

K β -K β = Φ T y -Φ T y -φ T 0 φ 0 β . (22)
Restricting to the overlapping set gives that

[K β ] -[K β] = [Φ T y -Φ T y -φ T 0 φ 0 β ]. ( 23
) Noting that [K] = [K ] and that [Φ T y ] -[Φ T y] = [φ T 0 y 0 ] we get that [K]([ β ] -[ β]) = [φ T 0 ](y 0 -φ 0 β ). ( 24
)
Analogously to the previous cases (y 0 -φ 0 β ) is a scalar which we can label C and we get

[φ T 0 ] = 1 C [K]([ β ] -[ β]). ( 25
)
Thus if K is known it can be used to deduce whether the additional individual has each of the SNPs in the overlapping set. If K is not known exactly it can be estimated from public data just as in the same SNP case.

G Algorithm for identifying unique genotypes of multiple dogs, when a precise K is known While the simple approach described in the main manuscript will work in many cases, there are a few special circumstances where a more complex algorithm may be required. In particular, it would not work if there are combinations of SNPs that are not observed among the individuals added to the database. For instance, if there is not a single SNP location where the first individual has a SNP variant and the others do not, then we would miss the corresponding value in C m . However, it is still possible to identify all the values in C m through a more complex algorithm:

1. First, extract all unique, non-zero values from d m .

2. Find the sum of all pairs of values in (1).

3. Find all values that are in (1), but not in (2). 3) and (6) that sum to the last element d m . There may be more than one set of values for which this is true. 8. If this search is unsuccessful, repeat steps 6-7. Eventually, a set of values summing to d m should be found. 9. If more than one possible set of values is found for C m in [START_REF] Huang | Genoguard: Protecting genomic data against brute-force attacks[END_REF], it is still possible to compare these sets and identify which is the most likely to contain the true values of C m . For each possible C m vector, a set of genotypes can be constructed for the m additional individuals. Using the frequencies of each SNP, it is possible to calculate the probability of observing each genotype. The set of values that produces the most likely genotypes for the m individuals is most likely to be the correct one.

Additionally, this algorithm depends on the fact that it is extremely unlikely that if someone were to sample three random continuous numbers i, j and k, it would just so happen that i + j = k. There is an extremely small chance that a value of C m would be un-discoverable because of a coincidence of this nature.

H Description of K when the genotypes are non-binary

In many cases, GRS are calculated on genotype matrices that are non-binary. In particular, they may take on three discrete values 0, 1 and 2, where 0 indicates that the most common variant is homozygous, 1 indicates that the individual is heterozygous for the uncommon variant, and 2 indicates that the individual is homozygous for the uncommon variant. If this is the case, the description of K will change. However, it is still the case that the entries of K depend only on the SNP frequencies and SNP cooccurrence frequencies in the dataset, and that knowledge of SNP frequencies and pairwise co-frequencies from the original study, are all that is required in order to compute K.

-For i = 1, . . . , N : K ii = p Aa + 4p AA where p aa is the frequency of individuals being heterozygous for the uncommon variant and p AA is the frequency of individuals being homozygous for the uncommon variant.

-For i = 1, . . . , N -1 and j > i: K ij = K ji = p Aa/Bb + 2p AA/Bb + 4p AA/BB where p Aa/Bb is the frequency that both SNPs are simultaneously heterogygous, p AA/Bb is the frequency that one SNP is homozygous for the rare variant and the other is heterogygous simultaneously, and p AA/BB is the frequency that that uncommon variants are found to be homozygous simultaneously.

-For i = 1, . . . , N and j = N + 1: K ij = K ji = p Aa + 2p AA .

-Finally, K N +1,N +1 = 1.

CFig. 2 .

 2 Fig. 2. (A) We have perfect accuracy in reconstructing the genotype when K is known (using 200 random SNPs to estimate average breed weight in the Cornell Dog Database). (B) We can reconstruct all the genotypes of multiple dogs that are added to the second study and (C) this works in practice using the data from the Cornell Dog Database, as in (A).

Fig. 3 .

 3 Fig. 3. Example values taken by the noisy vector d, given the true value of the corresponding SNP in the genome. (Left) adding one new participant;(right) adding three new participants. These figures are analogous to those in Figure2, albeit in the case where K is not known and instead estimated from an independent public database.

Fig. 5 .

 5 Fig. 5. Results of Figure 4 broken down by individual dogs.Here each point represents a dog and we define atypicality as the proportion of uncommon variants that the dog has compared to the public database-for instance, if 51% or more of dogs in the public database have a G in a specific locus, but this dog has a T, then this would count towards the dog's atypicality. In other words, dogs further to the right are less and less similar to average dog present in the public dataset (measured by percentage of different variants). In contrast to the most-common-variant baseline, our method generalizes well even to dogs which are highly dissimilar to those in the public dataset. Larger public databases (right) provide more accurate population estimates K, leading to more accurate reconstructions overall.

Fig. 6 .

 6 Fig.6. Accuracy at reconstruction of the genome of one additional individual, using EM estimation and a noisy estimate K, measured as the size of the initial private increases. For very small private databases, accuracy is very high, as changes in entries of β are clearly attributable to the new individual. Beyond a certain threshold, overall accuracy is quite stable. Error bars show mean and two standard deviations.

4 .

 4 If there are exactly m values in (3) and the sum of these values equal the last value of d m (corresponding to the intercept term), then you have found the correct values of C m . 5. Otherwise, this suggests that there are one or more elements of C m that are missing from (3) and possibly a few values in (3) that are not in C m . 6. Begin by subtracting every pair of values in (3). These are now also potential values of C m 7. Search for a set of m values from (

This project was funded by the Alan Turing Institute Research Fellowship under EPSRC Research grant (TU/A/000017); EPSRC/BBSRC Innovation Fellowship (EP/S001360/1), and under the EPSRC grant EP/N510129/1. It was also partly funded by a grant from CPER Nord-Pas de Calais/FEDER DATA Advanced data science and technologies 2015-20 . CC-BY 4.review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made


Mariana Vargas-Vieyra

Aurélien Bellet

Pascal Denis

Joint Learning of the Graph and the Data Representation for Graph-Based Semi-Supervised Learning

Graph-based semi-supervised learning is appealing when labels are scarce but large amounts of unlabeled data are available. These methods typically use a heuristic strategy to construct the graph based on some fixed data representation, independently of the available labels. In this paper, we propose to jointly learn a data representation and a graph from both labeled and unlabeled data such that (i) the learned representation indirectly encodes the label information injected into the graph, and (ii) the graph provides a smooth topology with respect to the transformed data. Plugging the resulting graph and representation into existing graph-based semi-supervised learning algorithms like label spreading and graph convolutional networks, we show that our approach outperforms standard graph construction methods on both synthetic data and real datasets.

Introduction

An important bottleneck for the development of accurate Natural Language Processing (NLP) tools for many applications and languages is the lack of annotated data. A natural remedy to this issue lies in semi-supervised learning (SSL) methods, which are able to leverage smaller labeled datasets in combination with large amounts of unannotated text data (which are often more easily available). In particular, graph-based SSL algorithms [START_REF] Subramanya | Graph-Based Semi-Supervised Learning[END_REF], among which variants of label propagation [START_REF] Zhu | Learning from Labeled and Unlabeled Data with Label Propagation[END_REF][START_REF] Zhu | Semi-supervised learning using Gaussian fields and harmonic functions[END_REF][START_REF] Zhou | Learning with Local and Global Consistency[END_REF]) and more recently graph convolutional networks [START_REF] Thomas | Semi-supervised Classification with Graph-convolutional Neural Networks[END_REF][START_REF] Chen | FastGCN: Fast Learning with Graph Convolutional Networks via Importance Sampling[END_REF][START_REF] Wu | Simplifying Graph Convolutional Networks[END_REF], have attracted a lot of attention due to their interesting theoretical properties and good empirical performance. They have successfully been applied to several NLP problems, such as sentiment analysis [START_REF] Goldberg | Seeing stars when there aren't many stars: Graph-based semisupervised learning for sentiment categorization[END_REF], word sense disambiguation [START_REF] Alexandrescu | Data-Driven Graph Construction for Semi-Supervised Graph-Based Learning in NLP[END_REF], text categorization [START_REF] Subramanya | Soft-supervised Learning for Text Classification[END_REF], POS tagging [START_REF] Subramanya | Efficient Graph-based Semi-supervised Learning of Structured Tagging Models[END_REF], semantic parsing [START_REF] Das | Semi-supervised frame-semantic parsing for unknown predicates[END_REF], machine translation [START_REF] Saluja | Graph-based semi-supervised learning of translation models from monolingual data[END_REF], or lexicon induction [START_REF] Faruqui | Morpho-syntactic lexicon generation using graphbased semi-supervised learning[END_REF]. As their name suggests, graph-based SSL methods represent all data points (that is, labeled and unlabeled) as nodes in a graph with weighted edges encoding the similarity between pairs of points. This graph is then used as a propagation operator to transfer labels from labeled to unlabeled points. Despite differences in the way this propagation is achieved, graph-based SSL approaches all rely on two assumptions: (i) the graph representing the data provides a faithful approximation of the manifold on which the data actually live, and (ii) the underlying labels are smooth with respect to this manifold.

In many NLP problems, there is often no a priori-known graph, which raises the question of how to best construct this graph over the dataset given some data representation. Most existing work rely on classic graph construction heuristics such as k-nn graphs, -graphs or radial kernel graphs [START_REF] Subramanya | Graph-Based Semi-Supervised Learning[END_REF], which may poorly adapt to the intrinsic structure of the data manifold and hence violate assumption (i). More recently, in the machine learning and signal processing communities, some algorithms were proposed to learn more flexible graphs by enforcing the topology to be smooth with respect to the input data [START_REF] Samuel | Fitting a graph to vector data[END_REF][START_REF] Kalofolias | How to learn a graph from smooth signals[END_REF][START_REF] Dong | Learning Laplacian Matrix in Smooth Graph Signal Representations[END_REF]. All these approaches This work is licensed under a Creative Commons Attribution 4.0 International License. License details: http:// creativecommons.org/licenses/by/4.0/.

heavily depend on the choice of data representation and disregard the label information, making them unable to adapt to the prediction task and therefore potentially violating assumption (ii). While supervised representation learning techniques such as metric learning [START_REF] Bellet | Metric Learning[END_REF] could be used to adapt the representation to the task of interest, for instance by bringing closer points with the same label, the lack of labeled data in the semi-supervised learning scenario makes them prone to overfitting.

In this paper, we propose an original semi-supervised algorithm for graph construction that adapts to both the data and the predictive task. Specifically, our approach leverages the labeled and unlabeled data to jointly learn a graph and a data representation. On the one hand, the graph is learned to provide a smooth topology with respect to the learned representation. On the other hand, the representation should bring closer (labeled and unlabeled) points that are neighbors in the graph as well as similarly labeled points, while pulling away points of different labels. A key feature of our approach is that the learned representation indirectly encodes and injects label information into the graph beyond the labeled points alone. We formulate our problem as a joint optimization problem over the representation and the graph weights, with a hyperparameter to easily control the sparsity of the resulting graph and thereby obtain a good approximation of the underlying manifold. We discuss some appropriate parameterizations for learning the representation, which revolve around adapting pre-trained embeddings so as to avoid overfitting. We then propose to solve our joint problem by alternating optimization on the representation and the graph. We validate our approach through several graph-based SSL experiments using label spreading [START_REF] Zhou | Learning with Local and Global Consistency[END_REF] and graph convolutional networks (GCN) (Kipf and Welling, 2017), both on synthetic and real text classification datasets. Incidentally, note that our approach is generic and could in principle be used in combination with any existing graph-based SSL framework. The results show that our approach outperforms previous methods which rely on heuristic graphs, generally by a considerable margin. Interestingly, we also observe that our approach effectively bridges the accuracy gap between a simple method like label spreading and a richer neural-based approach like GCN.

The rest of this paper is organized as follows. We introduce some notations and discuss the related work in Section 2. We then describe our approach and algorithm in Section 3. Our experimental results are presented in Section 4, and we conclude with future work directions in Section 5.

Notations and Related Work

Our work lies in the intersection of two topics: graph-based semi-supervised learning and graph learning. In this section, we start by introducing some useful notations. We then briefly review some related work in both areas.

Notations

We consider a dataset consisting of l labeled points L = {(x i , y i )} l i=1 and u unlabeled points U = {x i } l+u i=l+1 , where data points x i lie in some space X (typically X ⊂ IR d ) and labels y i ∈ {1, . . . , C} are discrete. We denote the combined data by X ∈ IR n×d , where n = l+u. We place ourselves in the typical semi-supervised scenario where l u, and the goal is to predict y l+1 , . . . , y l+u . Let G = (V, W ) be a graph composed of a set of nodes V = {v 1 , . . . , v n } and a symmetric nonnegative weighted adjacency matrix W ∈ IR n×n . We denote the set W of admissible weighted adjacency matrices by

W = {W : W ≥ 0, diag(W ) = 0, W = W }.
(1)

Every observation x i (labeled or unlabeled) can be seen as a signal occurring at node v i , and W assigns a weight to each pair of nodes. We say that two nodes v i and v j are connected when W ij > 0.

Graph-based SSL

Graph-based SSL takes as input a graph G = (V, W ) over the labeled and unlabeled data, the known labels and optionally some feature representation X ∈ IR n×d of the data, and aims to predict the labels of the unlabeled points. Many approaches exist for graph-based SSL, see [START_REF] Subramanya | Graph-Based Semi-Supervised Learning[END_REF] for a review of standard approaches. Classic algorithms take as input only the graph and use it to propagate the known labels to the unlabeled points [START_REF] Zhu | Learning from Labeled and Unlabeled Data with Label Propagation[END_REF][START_REF] Zhou | Learning with Local and Global Consistency[END_REF][START_REF] Bengio | Label Propagation and Quadratic Criterion[END_REF]. A popular approach is label spreading [START_REF] Zhou | Learning with Local and Global Consistency[END_REF], which can be formulated as a convex optimization problem implementing a trade-off between two terms. The first term is a graph Laplacian regularization term which encourages similar predictions for strongly connected points in the graph, while the second one tries to keep the predictions accurate for points with known labels. There also exist algorithms formalized as graph partitioning problems inspired from spectral clustering [START_REF] Joachims | Transductive Learning via Spectral Graph Partitioning[END_REF].

Graph-based SSL methods can also leverage a feature representation of the data in addition to the graph. A generic strategy is to add the graph Laplacian regularization term to existing supervised learning algorithms such as SVM or neural networks [START_REF] Belkin | Manifold Regularization: A Geometric Framework for Learning from Labeled and Unlabeled Examples[END_REF][START_REF] Weston | Deep Learning via Semi-supervised Embedding[END_REF]. Following a different direction, [START_REF] Yang | Revisiting semi-supervised learning with graph embeddings[END_REF] builds upon graph embeddings approaches to propose a method to predict labels based on the input representation as well as some embeddings learned from the input graph and known labels. Finally, there has been some recent interest in graph convolutional networks (GCN) [START_REF] Thomas | Semi-supervised Classification with Graph-convolutional Neural Networks[END_REF][START_REF] Chen | FastGCN: Fast Learning with Graph Convolutional Networks via Importance Sampling[END_REF][START_REF] Wu | Simplifying Graph Convolutional Networks[END_REF]. Much like CNNs for images, they rely on an (approximation of) a notion of graph convolution, allowing them to learn a nonlinear transformation of the input representation while encoding the graph structure when propagating inputs from a layer to the next.

In all these approaches, the graph is fixed and given as input to the algorithm (for several methods, this is also true for the data representation). Their performance is thus very sensitive to the relevance of the graph for the task at hand: in particular, the underlying labeling should be smooth with respect to the graph. Our method is a principled approach to learn such a task-specific graph, and also infers a relevant data representation.

Graph Construction and Learning

As pointed out by [START_REF] Subramanya | Graph-Based Semi-Supervised Learning[END_REF], graph-based SSL methods typically rely on graphs constructed from the input data representation with a simple heuristic strategy. Popular choices include k-nearest neighbor graphs (connecting pairs of points that are among the k-closest to each other), -graph (connecting points that are within distance ), and radial kernel graphs (building a fully connected graph with exponentially decreasing weights W ij = e -γ x i -x j 2 ). Recently, more sophisticated methods that learn the graph weights as the solution of an optimization problem have been introduced [START_REF] Samuel | Fitting a graph to vector data[END_REF][START_REF] Kalofolias | How to learn a graph from smooth signals[END_REF][START_REF] Dong | Learning Laplacian Matrix in Smooth Graph Signal Representations[END_REF]. Essentially, the weights are learned to be smooth over the data representation (i.e., assigning large weights to nearby points) with some regularization to enforce or control some properties such as connectedness and sparsity. In any case, the graphs obtained with the above approaches are task-independent in the sense that they ignore the labels.

To the best of our knowledge, there have been very few attempts to learn task-specific graphs for SSL. [START_REF] Alexandrescu | Data-Driven Graph Construction for Semi-Supervised Graph-Based Learning in NLP[END_REF] propose to train a supervised classifier on labeled points and using the soft label predictions as the representation to build the graph. While this gives a way to incorporate label information, the supervised predictions are very dependent on the initial representation and the classifier itself can heavily overfit due to scarce labels.

Proposed Model

Our approach learns a graph and a data representation for use in downstream graph-based SSL algorithms. In this section, we start by introducing our formulation as a joint optimization problem over the representation and the graph. We then discuss some relevant choices for the parameterization of the learned representation, and finally present our alternating optimization scheme.

Problem Formulation

For the sake of generality, in this section we formulate our problem with respect to a generic representation function φ Θ : X → IR k , parameterized by Θ, which represents any data point x ∈ X as a k-dimensional vector φ Θ (x) ∈ IR k . We discuss some relevant choices of representation functions in Section 3.2.

We propose to learn a weighted adjacency matrix W * and a representation function φ Θ * by minimizing a joint objective function f that involves both the labeled and unlabeled data points:

W * , Θ * = arg min W ∈W,Θ f (W, Θ).
Once the above optimization problem has been solved, the learned graph W * (which is based on the learned representation function φ Θ * ) and possibly the representation φ Θ * can then be given as input to any graph-based SSL algorithm to obtain predictions for the unlabeled data.

Our objective function f (W, Θ) decomposes into three terms:

f (W, Θ) = f 1 (Θ) + α[f 2 (W ) + f 3 (W, Θ)] (2)
where f 1 (Θ) and f 2 (W ) are respectively the representation and graph specific terms, while f 3 (W, Θ) is the joint term. Hyperparameter α ≥ 0 controls the trade-off between the (supervised) representation learning term f 1 and the unsupervised part (f 2 and f 3 ). We now define these three terms. For notational convenience, let us denote by Z ∈ IR n×n the matrix whose entries are the normalized squared Euclidean distances between data points in the transformed space, i.e.

(Z Θ ) ij = ||φ Θ (x i )-φ Θ (x j )|| 2 i<j ||φ Θ (x i )-φ Θ (x j )|| 2 .
The normalization conveniently removes the dependency on the scale of the data and Θ. The representation term f 1 (Θ) is defined on the labeled data points only and takes the following form:

f 1 (Θ) = x i ,x j ,x k ∈L y i =y j ,y i =y k (Z Θ ) ij -(Z Θ ) ik + 1 + , (3)
where

[•] + = max(0, •)
. This is a large-margin triplet loss similar to those used in metric learning [START_REF] Bellet | Metric Learning[END_REF]: it attempts to learn a representation function φ Θ that brings each point x i closer to points x j with the same label than to differently labeled points x k , with a safety margin of 1. In practice, we can subsample instead of summing over all possible triplets. The graph term f 2 (W ) is inspired from the (unsupervised) graph learning approach proposed by [START_REF] Kalofolias | How to learn a graph from smooth signals[END_REF]:

f 2 (W ) = β||W || 2 F -1 log(1 W ), (4)
The log-barrier term on the degrees prevents any node from being isolated in the graph, while the Frobenius norm is a shrinkage term over the graph weights. Combined with our joint term (6) defined below, hyperparameter β ≥ 0 directly controls the sparsity of the learned graph: the smaller β, the more concentrated the weights of each point on its nearest neighbors in the learned representation (hence the sparser the graph). On the other hand, as β → +∞, the graph becomes complete with uniform weights. Sparsity allows to enforce the locality property (only close points are connected in the graph) which is necessary to obtain a good approximation of the data manifold. It also reduces the computational cost in downstream graph-based SSL algorithms, whose complexity typically depends on the number of edges in the graph.

Other options are possible for f 2 (W ) depending on the prior we want to have on the structure of the graph. For instance, one may use

f 2 (W ) = (1/γ) i,j W ij [log(W ij ) -1], (5)
where γ > 0 is a hyperparameter. This will force the graph to be fully connected. Finally, we introduce the joint term bringing together the graph and the representation:

f 3 (W, Θ) = tr(W Z Θ ) = i,j W ij (Z Θ ) ij . (6)
This can be seen as a weighted L 1 norm term on W (which is why it induces sparsity), and equivalently written as a quadratic form of the Laplacian matrix of the graph encoded by the symmetric matrix W .

It is also used in approaches based on graph Laplacian regularization (see Section 2), but in our case both the graph and the representation are learned in joint manner. This term makes the graph and the representation as smooth as possible with respect to each other on both labeled and unlabeled points.

Overall, our joint objective function ( 2) is designed to produce a sparse topology that tends to be smooth with respect to the data manifold and the underlying labeling function through an appropriate representation. We now discuss the choice of representation function φ Θ .

Choices of Representation Functions

Many options are possible for the representation function φ Θ depending on the nature of the data and task at hand. However, it is important to keep in mind that the amount of labeled information is scarce, hence learning complex text representations from scratch is likely to lead to severe overfitting. We argue that it is preferable to adapt pre-trained representations, which generally requires to optimize much fewer parameters. We give some concrete examples below.

Linear transformation. Pre-trained word embeddings [START_REF] Mikolov | Efficient Estimation of Word Representations in Vector Space[END_REF][START_REF] Pennington | Glove: Global vectors for word representation[END_REF] are commonly used to represent texts in a vectorial space, e.g. by averaging the embeddings of the words occurring in a document. In order to adapt the representation to the task, we can learn a simple linear mapping φ Θ (x) = Θx which transforms the initial d-dimensional representation into a k-dimensional one, with Θ ∈ IR k×d and k ≤ d. Such a strategy has been previously explored in the supervised setting to "re-embed" words in a task-specific manner [START_REF] Denis | Online Learning of Task-specific Word Representations with a JointBiconvex Passive-Aggressive Algorithm[END_REF]. This is the representation function that we use in our experiments (see Section 4).

Weighted combination. Recent work in learning deep contextualized word representations such as ELMo [START_REF] Peters | Deep contextualized word representations[END_REF] and BERT [START_REF] Devlin | BERT: Pre-training of Deep Bidirectional Transformers for Language Understanding[END_REF] allows to learn a task-specific combination of the token representations obtained at the K layers of the model, which typically capture different aspects of tokens (from syntax to semantics). In this case, we have K initial d-dimensional representations x ∈ IR K×d for each text x and we learn a weighted combination φ Θ (x) = Θx ∈ IR d where Θ ∈ IR K is simply a K-dimensional parameter vector.

Optimization

We propose to optimize the cost function f (W, Θ) by alternating minimization over W and Θ, which is guaranteed to converge to a local optimum. This is a natural approach: one step learns a smooth graph given the current representation Θ, while the other learns a smooth representation with respect to the current graph (this can be seen as a regularizer for Θ based on unlabeled data) and also tries to keep labeled points of the same class closer than points of different class.

As the joint problem is nonconvex, initialization plays an important role. We propose to initialize the graph weights to zero and to start by optimizing Θ so that the initial representation focuses only on the (scarce) labeled data. The graph learned on this representation will thus strongly connect together the labeled points as well as unlabeled points that are very close to the labeled points and are thus likely to share the same label. At the next iteration, these unlabeled points will then contribute in learning a better representation and in turn a graph which strongly connects new unlabeled points. This process can be seen as a principled version of self-training heuristics popular in traditional (non-graph-based) semi-supervised learning [START_REF] Triguero | Self-labeled techniques for semi-supervised learning: taxonomy, software and empirical study[END_REF].

The subproblem of optimizing W given Θ is convex regardless of whether we define f 2 (W ) as (4) or (5). Using ( 5) is computationally convenient as the subproblem has a closed-form solution: the weights are exponentially decreasing with the distance in the current representation φ Θ , as given by the radial kernel W ij = exp(-γ(Z Θ ) ij ) [START_REF] Kalofolias | How to learn a graph from smooth signals[END_REF]. Note that unlike the classic radial kernel baseline construction method mentioned in Section 2.3, our graph is computed based on the learned representation φ Θ by minimizing the joint objective function with respect to W . One drawback of using ( 5) is that the resulting graphs are always fully connected. Using (4) instead, we can obtain sparse graphs but the solution must be computed with an iterative algorithm. We found that the primal-dual algorithm introduced by [START_REF] Kalofolias | How to learn a graph from smooth signals[END_REF] converges slowly in practice -we instead optimize W by simple gradient descent over the "effective" n(n -1)/2 weights, adding a small positive constant inside the log term in (4) to make the objective function smooth.

As φ Θ is typically differentiable in Θ (as in the examples outlined in Section 3.2), we also solve the subproblem in Θ by (stochastic) gradient descent. Note that this subproblem is generally nonconvex due to the distance difference in f 1 (Θ).

Remark. Updating W requires to optimize over O(n 2 ) variables, which was manageable for the datasets used in our experiments. To scale to larger datasets, one can restrict the optimization to the weights corresponding to pairs of points that are close enough in the learned representation space 1 (other weights are kept to 0). This has a negligible impact on the solution in sparse regimes (small β).

Experiments

In this section, we study the practical behavior of our method by comparing the accuracy of downstream graph-based SSL algorithms when the graph (along with the underlying representation) is learned with our approach ( ours) rather than constructed with the following baseline strategies:

• radial: Complete graph with weights W ij = exp(-γ x i -x j 2 ).

• knn: W ij = 1 for x i in the k-neighborhood of x j (or vice versa), and W ij = 0 otherwise.

• kalo: Unsupervised graph learning with the method of [START_REF] Kalofolias | How to learn a graph from smooth signals[END_REF]. This corresponds to our approach when using the graph term (4) and keeping the original representation fixed.

In all cases the graph is constructed over the union of labeled (train set) and unlabeled data (validation and test sets). For experiments with our method, the learned representation is a linear transformation of the initial features as explained in Section 3.2.

We perform experiments with two graph-based SSL approaches: Label Spreading (LS) [START_REF] Zhou | Learning with Local and Global Consistency[END_REF] and the Graph Convolutional Network (GCN) method of (Kipf and Welling, 2017). We used the scikit-learn [START_REF] Pedregosa | Scikitlearn: Machine Learning in Python[END_REF] implementation of LS. For GCN, we used the TensorFlow implementation provided by the authors 2 and follow the recommended architecture:

Ŷ = softmax( L max(0, LXH 0 )H 1 ), (7)
where L = D1/2 W D1/2 is the normalized Laplacian corresponding to the graph W = W +λI (the input graph augmented with self-loops), D is the diagonal degree matrix ( Dii = j Dij ), and H 0 ∈ IR d×h , H 1 ∈ IR h×C are the parameters to be learned. We set the number of hidden units h to 16 and λ to 1 as done in [START_REF] Thomas | Semi-supervised Classification with Graph-convolutional Neural Networks[END_REF].

To illustrate the behavior of our approach, we first present some experiments on synthetic data. We then show some results on real text classification datasets.

Synthetic Data

We generated a 3-dimensional dataset consisting of 100 points evenly distributed in two classes (Figure 1). We have two clusters per class placed far from each other while keeping clusters from different classes closer. We randomly picked 60% of the points and removed their labels.

We compare the classification error of GCN and Label Spreading when the input graph is given by our approach instead of using baseline graph construction methods. For GCN, we also give as input the representation learned with our approach. For our approach, we use the graph term (4) and for each labeled point x i , we random sample 2 points x j of the same class and 3 points x k of different class and construct all combinations (x i , x j , x k ), leading to 6 triplets for each x i in the triplet loss (3). The results given in Table 1 show that our approach clearly and consistently outperforms all methods in both GCN and Label Spreading. 3 The improvements are especially large for Label Spreading, as LS makes 1 These can be identified in near-linear time using approximate nearest-neighbor techniques [START_REF] Muja | Scalable Nearest Neighbor Algorithms for High Dimensional Data[END_REF]).

2 https://github.com/tkipf/gcn 3 For this illustrating experiment, we picked the values of hyperparameters giving the best results for each method. predictions based on the graph only. In contrast, GCN learns its own (nonlinear) transformation of the representation given as input in an end-to-end manner. Still, our method is able to provide some gains for GCN as well, by providing it with a better graph. Note for instance the significant improvement compared to kalo, which learns the graph on the original representation.

To visualize this difference, Figure 2a shows the graph learned by kalo. Although the graph is learned to minimize the smoothness criterion with respect to the data, it fails to accurately capture the label distribution due to the limitations of the initial representation. Our alternating optimization approach overcomes this issue by learning a task-specific graph through an appropriate representation. In Figure 2b-2c-2d, we can see how label information is gradually injected at each step: after the first iteration, the graph is already significantly more smooth with respect to the underlying labeling and the graph is also sparser, but some edges between differently labeled points as well as an overly connected point remain. The following iterations further improve the graph quality. This explains the better performance obtained in downstream semi-supervised algorithms.

Real Data

We now evaluate our method on three text classification tasks derived from the 20NewsGroups dataset,4 a collection of documents categorized into 20 topics, each one of which is partitioned into sub-topics. We chose the topics of computers with classes IBM and Mac (n = 1945 documents), religion with classes atheism and Christian (n = 1796), and sports with classes baseball and hockey (n = 1993).

For all datasets, we represent data points using the average token embedding based on word2vec [START_REF] Mikolov | Efficient Estimation of Word Representations in Vector Space[END_REF]. These embeddings are of dimension d = 300 and were trained on a 100B word corpus of Google news data (vocabulary size is 3M). 5We experiment with different proportions of unlabeled points in the training set (90%, 75%, 60% and 40%), while the rest of the data is evenly split into a validation and a test set. As commonly done in semisupervised learning, we train on the union of the (labeled) training set and the (unlabeled) validation and test sets, select the values of hyperparameters based on the accuracy on the validation set, and report the corresponding accuracy on the test set.

To evaluate our approach we optimize the objective (2) as described in Section 3.3 with the graph term defined as in (4). To compute the representation term of our objective defined in (3), we construct triplets as follows: for each pair (x i , y i ) in the labeled set we obtain the closest points with labels other than y i ("imposters"), and the closest points with label y i ("targets"). We picked 8 imposters and 3 targets. We tune the hyperparameters α from {0.001, 0.01, 1}, β from {0.00001, 0.001, 0.1, 1}, the dimension k of the learned representation from {16, 32, 64}, and perform early stopping with respect to the number of alternating steps between learning the graph and learning the representation (up to 10 alternating steps). We also tuned the hyperparameters of each baseline method (γ for radial, k for knn and β for kalo) and the trade-off hyperparameter of Label Spreading. Finally, we computed the McNemar test of significance [START_REF] Quinn | Note on the sampling error of the difference between correlated proportions or percentages[END_REF] to compare the performance of our method against the best baseline. Results marked with a dagger symbol † yield a statistically significant test for a significance level of 0.05.

Label Spreading. Table 2a reports test classification accuracies obtained on the test set for each configuration of dataset and proportion of unlabeled data. Our approach clearly outperforms all baselines, most of the time by a large margin. Also, McNemar test indicates that we tend to be significantly better than the best baseline in the more challenging settings where labeled data is the most scarce. The results also show that learning the representation along with the graph makes a clear difference compared to learning the graph only (as seen by the superior performance of ours over kalo).

As LS only uses the graph to make predictions, these results provide strong evidence of the superior quality of the graphs learned with our method.

Graph Convolutional Networks. We now turn to the more complex GCN prediction model. We reuse the same setup as for LS and feed GCN with both the learned representation and the learned graph.

Table 2b summarizes the results. The gains obtained with our approach are smaller than those ob- tained in LS, which is to be expected since GCN has the ability to learn nonlinear transformations of the data. Nevertheless, we do observe some performance gains, as our approach generally improves upon or closely matches the performance of the best baseline. An interesting finding is that our method tends to close the gap of performance between LS and the richer neural-based GCN model. This suggests that simple propagation approaches may be sufficient in practice for many datasets, if provided with the right graph.

Visualization. We provide visualizations of the representation and the graph learned with our approach on the rel dataset. Figure 3 shows 3D PCA visualizations of the original representation and the representation learned with our approach. We see that the two classes are quite mixed up in the original representation while the learned representation is much smoother with respect to the underlying labeling (even in this crude low-dimensional summary). Figure 4 gives a snapshot of the graph learned with our approach by showing a subgraph of 50 randomly sampled nodes (subsampling helps to avoid clutter). The graph is very smooth with respect to the underlying labeling, and suggests that the learned high-dimensional representation has a nice manifold structure, with some regions of higher densities.

Discussion and Future Work

We presented a novel method bringing together graph learning, representation learning and SSL by jointly inferring the graph and the representation from semi-supervised data. The output of our approach can then be plugged into any graph-based SSL algorithm in place of using common graph constructions. Our experimental results suggest that the gains are especially significant for graph-based SSL algorithms that are unable to adapt the data representation (like label spreading and its variants), although we observe some gains also for GCN. To further improve the performance with richer models like GCNs, a promising direction is to extend our approach to learn the graph, the representation and the classifier in an end-toend manner. We note that there has been very recent attempts in this general direction [START_REF] Franceschi | Learning Discrete Structures for Graph Neural Networks[END_REF], though specific to GCN and with completely different modeling and assumptions.

The ideas underlying our approach could also be useful to tackle transfer learning settings and in particular domain adaptation [START_REF] Ben-David | Transductive Learning via Spectral Graph Partitioning[END_REF]. The latter can be seen as a SSL problem where the distribution of the target (unlabeled) data follows a different distribution that the source (training) data. Our objective function could be modified to encourage the learned representation and graph to serve as a "bridge" between the source and target distributions.

Figure 1 :

 1 Figure 1: Original 3-dimensional synthetic data. Semi-transparent points are unlabeled.

  (a) kalo on original data. (b) ours after 1st iteration. (c) ours after 2nd iteration. (d) ours at last iteration.

Figure 2 :

 2 Figure 2: Force-directed drawing (spring layout) of graphs learned with kalo, and with our method at several iterations of our alternating optimization algorithm. Semi-transparent points are unlabeled.

Figure 3 :

 3 Figure 3: 3D PCA visualization of the original representation (left) and the representation learned with our approach (right) on the rel dataset (%75 unlabeled). Transparent dots represent unlabeled documents.

Figure 4 :

 4 Figure4: Force-directed drawing (spring layout) of a random 50-node subgraph of the graph learned with our approach on the rel dataset (%75 labeled).

Table 1 :

 1 Classification accuracy on the synthetic dataset.

		Label Spreading GCN
	radial	50.7	93.3
	knn	81.3	93.3
	kalo	77.3	88.0
	ours	96.0	96.0

Table 2 :

 2 Classification accuracies of Label Spreading and GCN for different graph construction methods and proportions of unlabeled data. McNemar test to compare ours vs. the best baseline is statistically significant for those results marked with a dagger symbol †.

http://qwone.com/ ˜jason/20Newsgroups/

https://code.google.com/archive/p/word2vec/

Acknowledgements

The authors would like to thank the Agence Nationale de la Recherche for funding this project with the grant number ANR-16-CE33-0011, and the reviewers for their feedback and suggestions. We also thank Grid5000 team for providing the necessary resources and technical support for this project, and CPER Nord-Pas de Calais/FEDER DATA Advanced data science and technologies 2015-2020.