extfieldsofstudy:
Computer Science
Biology

text:
Privacy-Preserving Identification of Target Patients from Outsourced Patient Data

With the increasing affordability and availability of patient data, hospitals tend to outsource their data to cloud service providers (CSPs) for the purpose of storage and analytics. However, the concern of data privacy significantly limits the data owners' choice. In this work, we propose the first solution, to the best of our knowledge, that allows a CSP to perform efficient identification of target patients (e.g., pre-processing for a genome-wide association study - GWAS) over multi-tenant encrypted phenotype data (owned by multiple hospitals or data owners). We first propose an encryption mechanism for phenotype data, where each data owner is allowed to encrypt its data with a unique secret key. Moreover, the ciphertext supports privacy-preserving search and, consequently, enables the selection of the target group of patients (e.g., case and control groups). In addition, we provide a per-query based authorization mechanism for a client to access and operate on the data stored at the CSP. Based on the identified patients, the proposed scheme can either (i) directly conduct GWAS (i.e., computation of statistics about genomic variants) at the CSP or (ii) provide the identified groups to the client to directly query the corresponding data owners and conduct GWAS using existing distributed solutions. We implement the proposed scheme and run experiments over a real-life genomic dataset to show its effectiveness. The result shows that the proposed solution is capable to efficiently identify the case/control groups in a privacy-preserving way.

INTRODUCTION
With the advent of precision medicine, healthcare providers are empowered with the ability to select treatments based on a genetic understanding of the patient's disease. For instance, for cancer patients, a potential treatment may be a combination of surgery, chemotherapy, radiation, and immunotherapy, depending on the type of cancer and its stage. Precision medicine can help decide on specific personalized treatment plans with certain drugs proving more effective treatment for specific genetic profiles.
However, to pave the way to such precision medicine, research over large volumes of patient data is required. This leads to the popularity of large-scale patient data sharing, such as Patient-Centered Clinical Research Network (PCORNet) [25] in the US, TranSMART [6] in the EU, and the Global Alliance for Genomics and Health (GA4GH) [12]. These data-sharing systems are required to comply with the increasingly stringent privacy regulations (e.g., HIPAA [16] or GDPR [20]).
Both centralized and decentralized approaches have been explored in this context. In the centralized schemes [18,23], all data owners encrypt their data and outsource them to a common repository, where collective computation becomes possible. In the decentralized schemes [11,21], all data owners store their data locally and run a distributed protocol to conduct the computation over all databases. Although both of these approaches provide the ability to conduct computation over combined data, both approaches assume that the client (who sends a query for the computation) already knows which database owners to query for the analysis. However, in real-life, there is a pre-computation step to identify which databases are utilized in the query processing. For instance, if a physician wants to identify similar patients to a given one based on the symptoms and then conduct some statistical tests on the identified patients, they should first contact all database owners to identify the useful databases. Similarly, to conduct a genome-wide association study (GWAS), a client first needs to identify which database owners have genomes belonging to the desired case and control group specifications.
In this work we aim to fill this gap. In general, we propose to develop efficient, privacy-preserving algorithms that can conduct this pre-computation (i.e., identification of target patients) at a centralized setting, at a CSP, over data from multiple database owners while also providing inherent access control. After the precomputation step, we support either (i) computation of statistics over the identified records at the CSP (e.g., for simple statistical operations, such as calculation of minor allele frequencies and chisquare values) or (ii) utilization of other distributed solutions (e.g., [11,21]) between the data owners of the identified records (for more complex operations).
In the following, we discuss the main challenges of conducting large-scale privacy-preserving identification of target patients over multi-tenant patient data at the CSP. Encryption of data using unique cryptographic keys. To protect data from a potentially curious CSP, each data owner (e.g., hospital) should encrypt its dataset before outsourcing. Moreover, each hospital should use its unique secret key for encryption. By doing so, hospitals can maintain data privacy even if some of the hospitals are corrupted. On the other hand, when each hospital encrypts its dataset with its own unique key, computation over a combined dataset becomes challenging. Efficient and privacy-preserving search over encrypted data across hospitals. All the patient data should be encrypted against a potentially curious CSP. However, the ciphertext should support secure search operation to facilitate selection of target patients (e.g., case and control groups for the association study) in a privacypreserving and efficient way. Moreover, to avoid issuing a separate query for each hospital's dataset, searchability of the encrypted data should be supported across hospitals. Fine-grained access control of data. Since the data is stored at the CSP, the data owner loses the direct control of the data. In order to guarantee that the data is accessed properly, the data owner needs to enforce an access policy over the outsourced data.
In this work, we propose a privacy-preserving framework for the identification of target patients over outsourced patient data. To the best of our knowledge, this is the first framework that tackles all the aforementioned challenges. In the proposed scheme, each hospital encrypts its own dataset with its unique key and outsources the storage and processing of the search operation over encrypted data. We mainly focus on identifying the target patients based on their encrypted phenotype data (which is typically the case in GWAS). For privacy of phenotype data, we encrypt them with a novel pairing-based encryption algorithm. The proposed encryption algorithm also provides the ability to efficiently search over encrypted phenotype data from several hospitals. This enables efficient identification of target group of patients (e.g., case and control groups) across hospitals. Furthermore, to let the phenotype data be properly accessible by legitimate/authorized clients, we introduce a fine-grained authorization scheme. For a single authorization request, each phenotype is considered as a unit of authorization. We implement and evaluate the proposed scheme using a real-life genomic dataset. Our result shows that the proposed scheme requires less than 2 seconds to identify the case and control groups (each of size 50) over a dataset with 1052 patients.
The rest of the paper is organized as follows. In the next section, we discuss the state of the art. In Section 3, we introduce the primitives we used in this paper. In Section 4, we present the models, including system model, data model, query model, and security model. Then, in Section 5 we describe the proposed scheme in detail. In Section 6, we provide the privacy analysis of the proposed solution. In Section 7, we evaluate the performance. Finally, we conclude the paper in Section 9.

RELATED WORK
With the increasing volume of patient data, there are two lines of research work to process the data. The first approach is the centralized solution, which requires large amount of data to be stored in a single repository. In order to match the requirement of storage and computation power, the CSP is the most popular central repository. However, concerns about data privacy in cloud storage arise due to malicious attacks from inside and outside of the CSPs [14]. To resolve the concern, Kim et al. [18] applied the BGV scheme [13] and YASHE scheme [8] to encrypt the patient data and conduct secure evaluation of 2 distribution over the encrypted data. Sadat et al. [23] proposed a hybrid system called SAFETY, which combines Paillier encryption scheme and Intel Software Guard Extensions (Intel SGX) to improve the efficiency of the chi-square test. The second approach is the decentralized solution, where each participant stores its data locally. In order to inter-operate the patient data over multiple data providers, Kamm et al. [17] proposed to secretly share the sensitive data among several parties and compute GWAS over the distributed data. Similarly, Bogdandov et al. [7] adopted secret-sharing based techniques to implement a privacypreserving framework for statistical analysis on federated datasets. Raisaro et al [21] combines homomorphic encryption and obfuscation techniques to achieve privacy-preserving medical data sharing among many clinical sites. Froelicher et al. [11] proposed a multiparty homomorphic encryption algorithm. Based on the algorithm, each party can store the data locally and be able to run analysis algorithms over all the participants' data without privacy violation.
In the previous work [26], we proposed a scheme to find similar patients based on genomic makeup. In that scheme, with an assumption that all the data owners share a secret key, all the data owners build an index for their data using the secret key to support similar patient search. In this paper, we attempt to remove not only the assumption, but also the index. Contribution. Compared with the existing work, our main contributions are as follows: (1) Our scheme supports multiple hospitals to outsource their data to the CSP and each of them encrypts its data with its own key while computation can be properly conducted over all the data.
(2) Previous work assumes that the case and control groups for the association study are already known. This cannot meet the requirement of dynamic selection of case/control groups. In our paper, we propose a scheme that supports to dynamically select the case and control groups, which can be easily integrated to the existing solutions of association study.
(3) The authorization mechanism in the proposed scheme is designed based on per query request and it supports fine-grained authorization.

BACKGROUND
In this section, we first briefly introduce the background of genomics. Then, we introduce symmetric bilinear groups applied in the proposed pairing-based encryption algorithm.

Genomics Background
The most common mutation in human population is called single nucleotide polymorphism (SNP). It is the variation in a single nucleotide at a particular position of the genome [22]. There are about 5 million SNPs observed per individual and sensitive information about individuals (such as disease predispositions) are typically inferred by analyzing the SNPs. Two kinds of nucleotides (or alleles) are observed for each SNP: (i) major allele is the one that is observed with a high frequency and (ii) minor allele is the one that is observed with low frequency. The frequency of the minor allele in a given population is denoted as the minor allele frequency (MAF). Each SNP includes two nucleotides, one inherited from the father and the other one from mother. For simplicity, we represent the value of a SNP as the number of its minor alleles, and hence ∈ {0, 1, 2}. A SNP is represented by an (ID, value) pair, where the ID is taken from a large standardized set of strings and the value is in {0, 1, 2}. In the following sections, if we mention a SNP (or SNPs) without mentioning the ID or value, we mean both parts.

Symmetric Bilinear Groups
Let be an additive group of prime order and 1 ∈ be the generators of and 2 ∈ is a random element from . Let : × → be a function which maps two elements from to an element in the target group having prime order . The tuple ( , , , ) is an symmetric bilinear group if following properties hold: (a) the group operations in , are efficiently computable. (b) the mapping from to is efficiently computable. (c) the mapping is non-degenerate: ( 1 , 2 ) ≠ 1. (d) the mapping is bilinear: for all , ∈ Z , ( 1 , 2 ) = ( 1 , 2 ) . Based on the symmetric bilinear group, we design an encryption algorithm to encrypt phenotypes of patients.

MODELS
In this section, we describe the system, data, query, and security models for our proposed scheme. We present the frequently used notation in Table 1.

System Model
As shown in Figure 1, the system consists of three types of entities: clients, hospitals, and a cloud service provider (CSP). The hospital collects biological samples from patients and sequences them with patients' consent. In parallel, the hospital records various phenotypes of patients (e.g., height, eye color, and blood type). Instead of storing all genotype and phenotype data locally, the hospital encrypts the data and outsources them to the CSP. The client (e.g., a medical researcher) queries the CSP for different association studies on datasets of several hospitals. Before sending a query to the CSP, the client needs authorization from the involved hospital(s). If the client gets such an authorization, the corresponding computation is allowed to be conducted over all hospitals' datasets. Upon receiving a query from a client, the CSP first constructs the identification of target group of patients (e.g., case and control groups) by running the query over the encrypted phenotype data, and then executes the other algorithms (e.g., GWAS) over the encrypted genotype data of individuals in the identified groups (e.g., case and control groups).

Data Model
In the proposed scheme, we make an assumptions to make sure that the data model is uniform across hospitals. We assume there exists a common set of terms applied to describe all the phenotypes across hospitals (referred as "phenotype attributes"). That is, all the hospitals use the same terms to describe the same phenotypes. For efficiency, we represent the value of a phenotype attribute in a binary format (as 0 or 1). 1 means that the patient matches the phenotype attribute while 0 denotes the opposite. For instance, Table 2 illustrates a partial taxonomy of phenotypes that includes different height ranges in centimeters, as well as presence of breast cancer.
The value of each SNP is set to represent its number of minor alleles (0, 1, or 2) (see the details in Section 3.1). Based on these settings, the phenotype and SNP dataset for hospital are represented as in Tables 2 and 3, respectively. In the following sections, when we mention phenotype, it means both phenotype attribute and value, and when we mention SNP, it means both SNP identity and value.

Query Model
The query mechanism is designed to find target groups of patients for specified phenotypes. In the proposed scheme, the query includes two parameters: (i) a list of phenotypes of interest and (ii) a A client sends an authorization request to a hospital for accessing its data. If the request is approved, the client gets authorization and is able to send request to the CSP to access the hospital's data.
parameter to set the size of matching groups (e.g., case and control groups). As an instance, we use case and groups to illustrate the query mechanism. For simplicity, we assume the size of case and control groups to be equal, but it can be customized to support different sizes. Prior to outsourcing the data to the CSP, a hospital first encrypts its phenotype dataset by using the proposed pairing-based encryption algorithm that supports efficient search over encrypted phenotypes (as discussed in Section 5). A query is generated by the client with input phenotype attributes and sent to the CSP. The proposed scheme allows a client to form the case and control groups based on multiple phenotype criterion (e.g., an association study for a particular type of cancer can be conducted only on males within a specific age range). The CSP, using the properties of the proposed pairing-based encryption, can check whether an encrypted patient's record (in a hospital's dataset) contains the phenotype attributes inside the query. If all the phenotype attributes in the query are included in a patient's record, the record is added into case group. In contrast, if a patient's record does not include any of the phenotype attributes in the query, the patient's record is added into the control group. Note that the query is encrypted without disclosing any phenotype information to the CSP and the CSP cannot learn any information from the search process. Thus, the CSP identifies the individuals in the case and control groups in a privacy-preserving way. The case and control groups may possibly contain patients from multiple hospitals. Specifically, given a query from a client and transform keys (detailed in Section 5.6) from hospitals that authorize the client to access their data, the CSP can search multiple hospitals' phenotype datasets (encrypted by using different secret keys). Thus, the search result is not limited to one hospital's dataset.

Threat Model
Our threat model is consistent with previous work [19,24]. Client, hospital(s), and the CSP are assumed to be semi-honest, that is, they honestly follow the protocol while trying to learn extra information during the protocol. A hospital may try to use its own data and knowledge to infer another hospital's data either via collusion with the CSP or exploring the common patient records in different hospitals. The CSP may analyze the stored ciphertext and observe the encrypted queries. Based on this information, the CSP may try to extract sensitive information. Also, a client may try to infer patients' data without having proper access authorization. Specifically, we consider following attacks: Ciphertext analysis: Since all the data are outsourced to the CSP, the CSP can run different algorithms to analyze ciphertext and try to extract meaningful information. Query analysis: Since all the queries are sent to the CSP, the CSP may analyze the received queries and their frequency. Consequently, the CSP may try to infer the query pattern and content. Operation analysis: The CSP conducts search computation and is able to obtain all the transcripts of the operation. Based on this information, the CSP may try to infer the content of query and stored ciphertext. Unauthorized access: Since all the hospitals' data is outsourced to the CSP, a client may try to access a hospital's data without authorization from that hospital. Collusion between hospitals: To infer a target hospital's data, several hospital may collude with each other and combine their knowledge. Collusion between hospitals and the CSP: If some hospitals and the CSP reach consensus on common interest to learn another hospital's (victim's) data, all the related parties combine their knowledge and try to infer the sought information. For instance, if the search (at the CSP) includes two hospitals and if one of these hospital collude with the CSP, the CSP can learn which patients of the other (victim) hospital has the considered phenotype as a result of the search operation. However, to provide the common search functionality, this attack is unavoidable, and hence we do not consider it in this work.
We thoroughly analyze all these attacks and robustness of the proposed scheme against them in Section 6.

PROPOSED SCHEME
In this section, we first give an overview of the proposed scheme, and then describe its details.

Overview
In the proposed scheme, hospitals independently encrypt and outsource their phenotype datasets to a CSP and the CSP conduct search operation over the outsourced federated data (from multiple hospitals) to identify target group of patients. To process phenotype data in an efficient and privacy-preserving way, we propose a novel encryption mechanism to encrypt the phenotype data. The proposed encryption scheme allows different hospitals to encrypt their phenotype data independently with their own secret keys. Moreover, the encryption algorithm supports identification of case and control groups efficiently, without information leakage. Furthermore, the identification process requires less communication compared with secure multi-party computation-based approaches [24] and less computation compared to homomorphic encryption-based approaches [4]. In general, the execution of privacy-preserving identification of target group of patients can be divided into seven phases: data preprocessing, initialization, key generation, data encryption, client authorization, query generation, identification of case and control groups. We now present a high-level overview of these seven phases while the rest of this section provides in-depth details. Data Preprocessing. Phenotype and genotype data need to be properly encoded in the required data format so that further processing can be conducted. Initialization. In this phase, the Setup is performed to initialize the parameters and functions. Key generation. In this phase, the KeyGen function is executed to generate secret keys for hospitals. Data encryption. Hospitals recursively call SinglePhenotypeEncrypt to encrypt their phenotype data. Client authorization. Before requesting access to a hospital's dataset (from the CSP), a client obtains an authorization from the hospital. Query generation. To run query over outsourced data with specified phenotypes, a query is generated by recursively running query generation algorithm SinglePhenotypeQueryGen with the input of target phenotypes. Identification of the case and control groups. Given a query, the CSP identifies whether a patient record is in case group or control group by running the Search algorithm.

Data Preprocessing
In this section, we present the procedures we use to preprocess the phenotype data, so that the processed data matches the data format requirements.

Preprocessing Phenotype Data.
The representations of phenotype attributes should be uniform for all the hospitals, and hence we assume that all the hospitals share a common ordered set of phenotype attributes, denoted as A. We also assume that the value corresponding to a phenotype attribute is binary: 1 represents a patient having such phenotype attribute, while 0 means the opposite. In Table 2, we illustrate the processed phenotype data for height and breast cancer. As seen in the table, in the dataset (Ψ ) of hospital , a patient record can be represented as , = { , , (P.name 1 , P.val ,1 ), · · · , (P.name , P.val , )}, where P.name ∈ A, P.val , ∈ {0, 1}, ∈ [1, ] and ∈ [1, ].

System Initialization
In this section, we present the procedures during the initialization so that all algorithms use the same initial parameters.

Initialization for Phenotype Data Encryption.
With the input of the security parameter , the system first generates a symmetric bilinear mapping : × → , where the multiplicative cyclic group is generated by generator and has the prime order ( > 2 ). Then, a cryptographic hash function is selected. Above procedures are detailed in Setup algorithm that is shown in Algorithm 1.

Key Generation
In this section, we present the procedures to generate the required keys in the system.

Key Generation for Phenotype Data Encryption.
Hospital randomly selects a master key ∈ , and sets its private key as = ∈ . In addition, hospital selects a secret key from {0, 1} . This is detailed in the Keygen procedure presented in Algorithm 2.

Key
Generation for client. The system randomly selects a number ∈ Z and generates a key = for each client. Once a client joins the system, it is assigned and .

Data Encryption
In this section, we describe the proposed phenotype data encryption algorithm, which supports efficient privacy-preserving search and update.
5.5.1 Phenotype Data Encryption. The set of phenotype attributes is denoted by A, as described before. Phenotype data of a patient in hospital is represented as , = { , , (P.name 1 , P.val ,1 ), · · · , (P.name , P.val , )}. For each phenotype of the patient (phenotype name-value pair, i.e., P.name , P.val , ), the hospital first selects a random number , , from Z , where is a large prime that is larger than 2 . The pair of phenotype attribute and value is first hashed and then encrypted into a group value , , . Afterwards, a symmetric encryption algorithm SE (e.g., AES) is invoked with the input of a secret key and the pair of phenotype attribute P.name and value P.val , . Finally, the algorithm outputs the ciphertextˆ, , consisting of a random group element , , , an encoded group element , , , and a symmetrically encrypted ciphertext¯, , .
Algorithm 3 shows how hospital encrypts a single phenotype of patient (i.e., P.name , P.val , ). To encrypt all patients' phenotype data in a hospital, we iteratively encrypt each patient's phenotypes. In detail, for each patient in hospital , the hospital first reads its phenotype data , , and then invokes the SinglePheno-typeEncrypt to encrypt each phenotype of the patient. Once all the phenotypes of the patient is encrypted, the result denoted as C , is stored in the dataset Dict with the patient pseudonym , . This process is repeated for each patient, detailed in Algorithm 4. for all (P.name , P.val , ) ∈ , do 5:ˆ, , ← SinglePhenotypeEncrypt(P.name , P.val , , ,

Client Authorization
In the proposed scheme, a client needs to get authorization from a hospital before it can access to (operate on) a hospital's data. Meanwhile, the hospital is capable to authorise the client in finegranularity and the CSP should not learn the data that the hospital authorises the client to access. In this section, we first present a simple client authorization mechanism, which requires a client to generate a query for each hospital that authorizes the client to access its data. In addition, this mechanism allows a client to access a hospital's data without any limitations once it is authorized. Then, to achieve flexible authorization and let the client generate a single query that can be used to operate on all authorized hospitals' datasets, we further present an improved mechanism. The improved mechanism supports per-query based fine-grained authorization and it allows a single query to access multiple hospitals' data.
5.6.1 Simple Authorization. In the simple authorization mechanism, if a client gets authorization from a hospital, the client can access all the data of the corresponding hospital for all future queries.
To present this authorization mechanism, we assume that there exists a client with private key = who wants to access (operate on) data from hospital . The client first sends an authorization request to hospital with its private key . If the hospital approves the request, the hospital signs the private key , = and sends it back to the client. The client recovers = 1/ , . Upon obtaining , the client can generate a legitimate query. The details of this authorization mechanism are also shown in Algorithm 5.

Algorithm 5 SimpleAuthorization
Input: , , Output: 1: Client: send = to hospital 2: Hospital : compute , ← and send , to the client 3: Client: compute ← 1/ 5.6.2 Improved Authorization. In the simple authorization mechanism, a client can access hospital's data indefinitely once the hospital authorizes the client. Also, the hospital is unable to control the granularity of the authorization. In other words, the hospital authorizes either all its data to a client or none of them. Furthermore, simple authorization mechanism results in a query that cannot be executed across multiple hospitals' data (i.e., the client needs to generate separate queries for each hospital's dataset).
In order to overcome these concerns, we propose an improved authorization mechanism. The new authorization mechanism supports fine-grained authorization on a per-query basis. Moreover, using the transform key function, a single query sent to the CSP can be executed across hospitals. As a result, the complexity of a client's query generation is reduced from O( ) to (1), where is the number of accessible hospitals. The steps for the improved authorization mechanism are as follows. A client first generates an authorization request and sends it to a hospital. If the hospital approves the request, the hospital first generates a transform key for the client (per query) and sends it to the CSP, such that the CSP can apply the transform key to transform the ciphertext into the format that supports the client's query. After the transformation, the query can be executed by the CSP over the phenotype data of each hospital that authorizes the client. The construction of the transform key needs to consider the privacy of the transform key, transformation process, and the query. With these privacy requirements in mind, we set the transform key at the granularity of phenotypes and compute it as , , = − , ( − ) − , , where denotes the target hospital, represents the associated data from a client, and represents the phenotype attribute in phenotype attribute set A. The details of the transform key construction are as follows.
(1) A client randomly selects , ∈ Z (1 ≤ ≤ ) and sends it to a hospital with and , (1 ≤ ≤ ) to request corresponding data access.
(2) Upon receiving the query request, the hospital decides which phenotype attributes can be accessed and generates , , (1 ≤ ≤ ) for each approved phenotype attribute. Specifically, if the hospital approves the requested , , the corresponding , , is properly constructed based on , and , . Otherwise, a random group element is selected, such that the access structure can be protected from the CSP.
(3) For each phenotype attribute, a transform key , , (1 ≤ ≤ )) is generated. The hospital sends all of them to the CSP and replies to the client with a success message.
Given the transform key, the CSP can transform the ciphertext into the format that supports corresponding client's query. Using this technique, the CSP can process all the ciphertext before a query is launched. The details of the transformation are discussed in Section 5.8.

Query Generation
The query is generated based on client's input that specifies phenotype data of interest. Then, the query is executed over encrypted phenotype data. Algorithm 6 shows query generation based on the input of a single phenotype. The client first applies the hash function on the input phenotype attribute and value, obtaining , ( ∈ A) as the output. Then, the client computes a group element , , where , is selected during the generation of authorization request. , is included in the query. Afterwards, the client computes the other part of the query: − , , , . To support multiple phenotypes in a query, the client can iteratively invoke the SinglePhenotypeQueryGen algorithm. In detail, given a set Δ of phenotypes, for each phenotype inside Δ, the client calls Sin-glePhenotypeQueryGen to encrypt it. The result , (1 ≤ ≤ ) is stored in a vector v , which has the same dimension and the same order of phenotype attributes as A. For the phenotypes whose phenotype attributes are in A but not in Δ, the corresponding value is set to 0 in vector v .
is a set and it includes the vector v . Finally, in addition to v , the size of case and control groups is also included in before is sent to the CSP. Algorithm 7 shows how to extends the input from a single phenotype to multiple phenotypes. Here, for each phenotype inside the input phenotype data set Δ, the client invokes Algorithm 6. Once all the input data are processed, the algorithm outputs the final query . Note that all the positions, where phenotype attributes are not included in the input phenotype data are filled with zeros. For example, a client can input a set of pairs of phenotype attribute and value as {(breast cancer, 1), (lung cancer, 0), (blue eye color, 1)} to generate a query. Applying above algorithms, each target phenotype attribute is corresponding to a component , in the query, the remaining that is not included in the input phenotype attribute set is set to 0. The final query also includes the size of case and control groups.

Identification of Case and Control groups
Without loss of generality, we present the identification process at the CSP over a single hospital 's dataset. Multiple hospitals' datasets are processed separately and in parallel. We describe the identification process in two steps. First, we show the search operation over a single phenotype. Then, we extend the algorithm to multiple phenotypes. The details of search over a single phenotype are shown in Algorithm 8. Specifically, given the ciphertextˆ, , of phenotype of patient in hospital , we first extract the components , , , , , ,¯, , fromˆ, , , denoted asˆ, , ,0 ,ˆ, , ,1 , andˆ, , ,2 . Based on the input of the query , from client querying for phenotype , we extract , , − , , , from , , denoted as , ,0 and , ,1 . With the three pairs (ˆ, , ,0 , , ,0 ), (ˆ, , ,1 , , ,1 ), and (ˆ, , ,0 , , ), we compute the bilinear mapping over each pair and multiply all the results one by one. If the computed result is 1, the current phenotype matches the query, otherwise, it does not match. The correctness of the algorithm is shown in Eq. 1. To support a query that contains multiple phenotypes, we extend Algorithm 8 to Algorithm 9. In detail, we first initialize two lists Υ case and Υ control . Then, for each patient with pseudonym , in hospital 's encrypted phenotype dataset Dict , the corresponding encrypted phenotype data , is read. Given the encrypted phenotype data , , query , and transform key ( , ,1 , · · · , , , ), the per-hospital components of the given data are extracted and passed as the input to the SinglePhenotypeSearch algorithm. For example, both the first element of phenotype data query and transform key are extracted, and then input into the SinglePhenotypeSearch algorithm. If a patient contains all the phenotype data inside a query, the patient's pseudonym , is added to the case group Υ case . If a patient does not match query phenotypes, its pseudonym is added to the control group Υ control . This process is executed until the size of both case and control groups reaches or until all the data is processed.

Computation over Identified Target Patients
Here, we briefly discuss how our scheme can compute GWAS (or other statistics) over identified patients in centralized and decentralized approaches. In a centralized setting, the patient genomic data is also stored at the CSP, and the genomic data of each patient can be encrypted using multi-key fully homomorphic encryption mechanism (e.g., [9]) for storage at the CSP. Then, after identifying the case and control groups (as in Section 5.8), the CSP can directly execute the GWAS algorithm over the identified patient data using the homomorphic properties of multi-key fully homomorphic encryption mechanism. In a decentralized setting, the patient genomic data is not stored at the CSP. In this case, the CSP can return the identified patients to the client. Based on the identified patients, the client can send requests to all the involved hospitals to get access to their data for computing GWAS in a distributed way [11,21].

Managing Dynamic Phenotype Data
Here, we show how the proposed scheme supports efficient update of patient phenotype data stored at the CSP. Assume the hospital wants to update phenotype of patient . Given the patient pseudonym , and the phenotype (P.name , P.val ), the phenotype encryption algorithm (in Section 5.5) is called to encrypt the phenotype. Once the encryption is completed, the vector of the update query is constructed by inserting the ciphertext at the position where P.name is located in A and setting the remaining values to 0. Additionally, the command "update" is added to the query. The update query is sent to the CSP. Upon receiving the update query, the CSP first identifies the entry of the patient , , then identifies the location corresponding to the non-zero element inside the vector of the update query, and finally replaces the old cipher with the new cipher from the update query. To insert a new phenotype or to delete an existing phenotype from a patient record, the proposed update algorithm can also be used by directly appending or removing a record from stored ciphertext.

PRIVACY ANALYSIS
In this section, we analyze the privacy of phenotype data. We first provide a high level discussion on how the proposed scheme achieves robustness in the presence of attacks presented in Section 4.4. Then we present formal privacy definition and proof.

Privacy Against Ciphertext Analysis
Phenotype data is encrypted and stored at the CSP. This allows the CSP to analyze the stored ciphertext. The encrypted phenotype data can be split into two parts (as in Section 5.5.1). The first part of the encrypted phenotype data is constructed in two steps. Hospital first computes the hash of the phenotype attribute and value. Then, hospital randomly selects a number to mask above hash result and raised the result to the power of a group element. The privacy of this part relies on the hardness of discrete logarithm problem, one-wayness of the hash function, and randomness of the selected number. The second part of the encrypted phenotype data results from directly encrypting phenotype attribute and value. Hospital uses its secret key and invokes symmetric encryption algorithm (e.g., AES) to encrypt the concatenation of phenotype attribute and value. The privacy of the second part relies on the robustness of the utilized symmetric encryption algorithm (e.g., AES). From the above description, we can conclude that both parts of the encrypted phenotype are semantically secure. Thus, the CSP is unable to learn significant information from the ciphertext analysis.

Privacy Against Query Analysis
The input of the query includes a set of phenotypes. Each phenotype includes a pair of phenotype attribute and value, which is first hashed and the hash result is lifted as a power of a group element (as in Section 5.7). Additionally, a random mask is selected to hide this result, which enables the encrypted query to be semantically secure. Due to the semantic security of the query, the CSP is unable to learn the query content from the query analysis.

Privacy Against Operation Analysis
Since the ciphertext of phenotype genotype data is stored at the CSP, the CSP is responsible for conducting search over the ciphertext. For the search operation, the CSP applies the query to search over the ciphertext of phenotype data. If a phenotype is included in the query, the corresponding search result is 1, otherwise, it is 0. Through the execution of the search process, the CSP can learn the number of matching phenotypes of each patient. However, for each matching phenotype, its value can either be 0 or 1, and the CSP cannot distinguish between the two. Thus, the CSP cannot learn any information about the query and ciphertext from the search operation.

Robustness Against Unauthorized Access
Access control is designed through the transform key (as discussed in Section 5.6). A client selects random numbers , (1 ≤ ≤ ) and sends them to the target hospital. The hospital generates the transform key ( , , = − , ( − ) − , ) by lifting these numbers into the power of a group element. Due to the randomness of , and the hardness of the discrete logarithm problem, the CSP is unable to extract any information from the transform key. Analyzing the search operation, the output of using incorrect transform key is 0, which does not disclose any information to the CSP, as described in Section 6.3. Therefore, the proposed authorization scheme is robust against the unauthorized access.

Robustness Against Colluding Hospitals
Each hospital independently encrypts its phenotype data and genotype data with its own secret key. Even if several hospitals collude with each other, they cannot get any advantage to infer another hospital's data.

Robustness against Collusion between a Hospital and CSP
If one or more hospitals collude with the CSP, the CSP cannot obtain any advantage to infer the remaining hospitals' data since each hospital's data is encrypted independently. However, as we clarified in Section 4.4, we do not consider following case. For instance, if the search (at the CSP) includes two hospitals and if one of these hospital collude with the CSP, the CSP can learn which patients of the other (victim) hospital has the considered phenotype as a result of the search operation (but not the genomic data of the identified patients). To provide the common search functionality, this attack is unavoidable, and hence we do not consider it in this work.

Privacy Analysis
In the following, we provide a formal privacy analysis of the proposed scheme. Following previous works [9], the allowed leakage includes (i) size pattern and (ii) access pattern. The size pattern discloses the size of the ciphertext, while the access pattern reveals the access frequency of matching patient records. The allowed leakage is not considered violation of our privacy goal. The privacy of the proposed scheme is based on the hardness of discrete logarithm problem, the randomness of selected random mask, and the robustness of applied symmetric encryption (e.g., AES). The privacy of the proposed scheme can be divided into two elements: phenotype data privacy, and query privacy. The privacy of phenotype data can be further divided into two parts. One part is encrypted by using symmetric encryption algorithm (the third element in a ciphertext, detailed in Section 5.5.1) and the other part is not (the first and second elements in a ciphertext, detailed in Section 5.5.1). Due to the robustness of symmetric encryption algorithm (e.g., AES), the part with the symmetric encryption is also semantically secure. The other part is first protected by a hash function and then masked with a random value. Both the hash result and random mask are put into the power of a group element. Due to the random mask and the hardness of discrete logarithm problem, the ciphertext is semantically secure in the presence of chosen plaintext attack. The query privacy is achieved similarly, relying on the random mask and discrete logarithm problem.
Formally, we define the privacy experiments as follows. Let Π be the scheme, the advantage of the adversary is defined as ADV Π A ( ) = ( * = ) − 1/2, where and * are defined below. In the following, we detail the game. Init: The adversary A selects two datasets 0 and 1 with same size and sends them to the challenger. Setup: With the input of security parameter , the challenger runs Setup to initialize the parameters. Then, the challenger calls KeyGen to generate the keys. Phase 1: A is allowed to ask the following request: Phenotype data encryption request: A is allowed to send a dataset with phenotype data to ask for encryption. The challenger calls Encrypt algorithm to encrypt the dataset and sends the result back to A. Challenge: The challenger randomly selects from {0, 1}, encrypts the dataset , and sends it to the adversary A. Phase 2: A is allowed to send requests as in Phase 1. Additionally, A is allowed to send a query request. The challenger only authorizes a query containing the phenotype attributes, where two datasets have the same value. Then, it generates a transform key for A, where the mask ( , ) in transform key is randomly selected by the challenger (see details in Section 5.6 ). Guess: A outputs * as the guess for .
We say the scheme Π is privacy-preserving if the advantage of the adversary is negligible, i.e, Π A ≤ ( ), where ( ) is a negligible function in .
From above defined privacy game, the adversary A is only allowed to learn the information from Phase 1 and Phase 2. The difference of Phase 2 from Phase 1 is that A holds the challenged ciphertext and is allowed to ask a query request. Thus, we only need to prove that what the adversary A can learn from ciphertext request and query request is negligible as follows.
Phenotype data encryption request: A submits the dataset of phenotype data to ask for encryption from hospital . If the PhenotypeEnc algorithm is semantically secure, A is unable to distinguish ciphertext from a random string. The ciphertext of each pair of phenotype attribute and value includes three components. The first component , , is randomly selected from Z , which does not reveal any information. The second component , , is − , , /( − , , ) . Due to the hardness of discrete logarithm problem, A is unable to extract the power of a group element. Similarly, A is unable to distinguish the second component from a random element in G. The third component is encrypted by using a symmetric encryption algorithm (e.g., AES), which is semantically secure. Therefore, the ciphertext obtained through Encrypt algorithm is semantically secure.
Query request: First, the transform key is indistinguishable from a random element of group . Second, for each pair of phenotype attribute and value, the query is ( , , − , , , ), where is an element from . Based on the hardness of discrete logarithm problem and the randomness of , , the query is indistinguishable from two random elements from . Third, given the query and transform key, the adversary A is capable to run search operation over the ciphertext of phenotype data. Furthermore, A is also capable to run analysis algorithms on ciphertext of genotype data. However, due to the constraint of issuing client authorization, two datasets should output the same search result. Thus, A cannot learn any significant information via executing search operation.
Based on above analysis, we can conclude Π A ( ) is negligible and the proposed scheme is privacy-preserving.

EVALUATION
In this section, we evaluate the performance of the proposed scheme. We run the experiments on a commodity machine with 7 CPU and 16GB RAM. The proposed phenotype encryption algorithm is implemented by Charm [5] written in Python while the SNP encryption algorithm is implemented by HEAAN [1] written in C++. Each experiment is run 10 times; we report the average results.

Data Model
We use the rsnps tool [2] to obtain all the raw patient files from the publicly available OpenSNP dataset [3]. Then, we converted the raw patient files into the VCF format using an open source software named personal-genome-analysis [15]. Eventually, we ended up with 1052 valid VCF files. For the phenotype data, we also used the OpenSNP dataset. In total, we collected 1052 records and extracted 1052 phenotype attributes.

Experimental Results
In this section, we first show the efficiency and scalability of the phenotype and genotype data encryption algorithms. After that, we present the scalability and efficiency of client authorization and query generation.   We adopt symmetric pairing group SS512 to construct the bilinear mapping and AES to implement the symmetric encryption. Accordingly, the time required to encrypt the phenotype data can be divided into two constituents. The first constituent is due to the pairing group operation. The second constituent is due to using AES to encrypt phenotype attribute and value. The secret key of AES is set to 256 bits. Table 4 shows the performance of the phenotype data encryption for different number of phenotypes. We observed that with the linear increase in the number of patients, the time cost of phenotype encryption for both AES and non-AES (pairing based encryption) constituents increases linearly. Similarly, when the number of phenotypes increases, the required encryption time also increases linearly.

Client
Authorization. The process of client authorization can be divided into two phases. In the first phase, the client generates an authorization request while in the second phase, the hospital generates the transform key. Table 5 shows the efficiency of authorization request generation for different number of requested phenotypes. We observe that the required time of authorization request generation increases linearly upon increasing the number of requested phenotype attributes. Table 6 shows the performance of the transform key generation for different number of phenotype attributes. Here, we observe that the time required for transform key generation increases linearly as a function of the number of phenotype attributes.

Query Generation.
The performance of the query generation algorithm is affected by the number of input phenotypes. The experimental results are shown in Table 7. From the table, we observe that with the linear increase in the number of input phenotypes, the time required for the query generation also increases linearly.

Phenotype Data
Search. The search process involves a number of patient records to be processed to construct the case and control groups. To access each hospital's data, the query needs to be transformed using the transform key. In addition, the desired size of case and control groups can be a factor to stop the search process earlier once the required number of individuals are identified in case and control groups. Table 8 shows the effect of the number of hospitals, number of input phenotypes, and the size of case and control groups on the efficiency. We observed that if the number of patients is fixed, the number of hospitals almost does not affect the efficiency, while the performance is sensitive to the size of case and control groups and the number of input phenotypes. When the size of case/control groups are set to 10 and 50, the search time is reduced to 0.29s and 1.78s, respectively. The reason is that once the required number of patients are identified for the case and control groups, the search process is terminated. We also evaluate the performance for reduced number of input phenotypes. The observed search times are 32.2s, 166.7s, and 327.1s for 10, 50, and 100 phenotypes, respectively. From these results, we can say that the search time increases linearly as the number of input phenotypes grows.
To show the efficiency of the proposed algorithm for phenotype data search, we also designed a following fully homomorphic encryption (FHE)-based version of it for comparison. In detail, the alternative FHE-based approach includes four steps. First, all the phenotype data is encrypted by using CKKS [10]. Second, the client sends a query to the CSP in order to determine the case/control groups. Third, the CSP sends the computed result to different hospitals. Fourth, hospitals decrypt the result in parallel and send the result (identified case/control groups) to the client. Without considering the time cost of communication, Table 9 shows the required time to complete the search operation using this FHE-based scheme for different number of patients and phenotypes. Comparing Table  8 with Table 9, we observe that the proposed algorithm is more than 20 times faster than the FHE-based algorithm when the number of queried phenotypes is at most 10. Notably, when the size of case/control groups is smaller than 50, the proposed algorithm is more than 300 times faster. Furthermore, in Table 10 we show the comparison between the FHE-based solution and proposed algorithm for different number of patients. From the table, we observe that (i) the proposed scheme consistently exhibits similar performance advantage over the FHE-based scheme and (ii) the time cost of both algorithms increases linearly with the number of patients.

DISCUSSION
The proposed pairing-based encryption scheme for phenotype data (in Section 5.8) is not limited to efficient identification of case and control groups. The scheme can be extended to support additional functions, such as similar patient search, target record retrieval, and statistical analysis of phenotype data. In the following, we discuss some of these potential extensions. Similar patient search. In the proposed scheme, a client is allowed to input several phenotypes of interest into a query and send it to the CSP. Upon receiving the query, the CSP searches the stored datasets of several hospitals. Based on the number of matching phenotypes, the CSP can order the search results and return similar patient records. Similar functionality can also be provided in the genome level by encrypting genome data with the proposed pairingbased encryption. Phenotype data retrieval. In the proposed scheme, we show how to use the Search algorithm (in Section 5.8) to find target pseudonyms based on the query. One can extend this function to support a client to retrieve phenotype data of interest. For instance, a client (e.g., a physician) may be interested to know the phenotypes of patients having lung cancer. Then, the client can generate a query and send it to the CSP to retrieve the phenotype data from patients that are diagnosed with lung cancer. Once the client receives the phenotype data from the CSP, the SinglePhenotypeDecrypt algorithm 10 is called to decrypt the phenotype data.

CONCLUSION
In this paper, we have designed a privacy-preserving framework for identification of a target group of patients across multi-tenant data. To achieve this, we have proposed a novel phenotype encryption algorithm. To support search and computation over multi-tenant data by a cloud service provider (CSP), we have introduced a transform key to enable the CSP to transform a single query and execute it over different hospitals' datasets without privacy violation. To manage the authorization of clients, we have proposed a per-query based authorization mechanism supporting selective phenotype data authorization. Via simulations on real genomic data, we have shown the practicality and efficiency of the proposed scheme. We believe that the proposed scheme will further facilitate the use of genomic data in clinical settings and pave the way for personalized medicine. In future work, we will focus on improving the search efficiency of genomic data and batch queries.


extfieldsofstudy:
Mathematics
Computer Science

text:
On optimal convergence rates for discrete minimizers of the Gross-Pitaevskii energy in LOD spaces

In this paper we revisit a two-level discretization based on the Localized Orthogonal Decomposition (LOD). It was originally proposed in [P.Henning, A.M{\aa}lqvist, D.Peterseim. SIAM J. Numer. Anal.52-4:1525-1550, 2014] to compute ground states of Bose-Einstein condensates by finding discrete minimizers of the Gross-Pitaevskii energy functional. The established convergence rates for the method appeared however suboptimal compared to numerical observations and a proof of optimal rates in this setting remained open. In this paper we shall close this gap by proving optimal order error estimates for the $L^2$- and $H^1$-error between the exact ground state and discrete minimizers, as well as error estimates for the ground state energy and the ground state eigenvalue. In particular, the achieved convergence rates for the energy and the eigenvalue are of $6$th order with respect to the mesh size on which the discrete LOD space is based, without making any additional regularity assumptions. These high rates justify the use of very coarse meshes, which significantly reduces the computational effort for finding accurate approximations of ground states. In addition, we include numerical experiments that confirm the optimality of the new theoretical convergence rates, both for smooth and discontinuous potentials.

Introduction
In this paper we consider the stationary Gross-Piteavskii equation (GPE), which can be written as an eigenvalue problem with eigenvector nonlinearity seeking u ∈ H 1 0 (Ω) with u L 2 (Ω) = 1 and λ ∈ R such that −∆u + V u + β|u| 2 u = λu.
Here, Ω ⊂ R d (with d = 1, 2, 3) is a computational domain, V (x) models a trapping potential and β ∈ R ≥0 is a constant. The perhaps most prominent application of the GPE is the modeling of so-called Bose-Einstein condensates [11,19,23,42]. Bose-Einstein condensates (BECs) are formed when a dilute gas of bosons is cooled to ultra-low temperatures so that almost all particles occupy the same quantum state, i.e. they become indistinguishable from each other and behave like one giant "macro particle". The first experimental realization of a BEC goes back to the Nobel prize-winning works by Cornell, Ketterle and Wieman [5,17]. Bose-Einstein condensates are of great interest in the study of macroscopic quantum effects such as superfluidity [2,40] and the closely related phenomenon of superconductivity [24].
In the context of BECs, β ≥ 0 models a repulsion parameter that depends on the number of bosons, as well as their type and mass. In this setting, the GPE has a positive spectrum, where the smallest eigenvalue is simple. This eigenvalue λ is called the ground state eigenvalue of the problem and it describes the chemical potential of a BEC in the lowest energy state. The corresponding ground state eigenfunction u is unique up to its sign. An equivalent characterization of the ground state u is to find a minimizer of an energy functional E (defined in equation (2.1) below) over the space of functions in H 1 0 (Ω) with the normalization constraint u L 2 (Ω) = 1.
Due to the complexity of physical experiments, computing the ground sate numerically is a very important task. Starting from the energy minimizing perspective, a numerical method has two essential components. The first component is an iterative method that allows to find minimizers of the energy E on the constrained manifold. Here there exists a variety of well-established schemes, where we exemplarily refer to methods based on self consistent field iterations (SCF) [13,14,18,43], methods based on discrete Sobolev gradient flows / Riemannian optimization [7,6,4,9,10,15,16,30,25,33,45] or the "J-method" [32,3]. We shall however not further discuss this aspect in our paper. Instead, we shall focus on the second important component of any numerical method for the computation of ground states: the choice of a suitable discrete space in which the energy functional is minimized. Standard choices are Lagrange finite element spaces (which can be often associated with finite difference discretizations) or Fourier spaces. These standard choices were analytically explored in [12]. In theses discrete spaces, each iteration of a scheme for finding a discrete minimizer of E requires the solution of a (typically elliptic) problem. Hence, the complexity of the final method depends naturally on the choice and the dimension of the discrete spaces.
In addition to the aforementioned standard choices, the usage of an alternative discrete space was proposed in [27]. Besides achieving high convergence orders, the proposed spaces have the advantage that they achieve their full potential already under minimal regularity assumptions on u and the potential V . Their construction is based on techniques from numerical homogenization, which build a low-dimensional generalized finite element space using a localized orthogonal decomposition (LOD). The LOD splits an ideal solution space into a low-dimensional approximation space and a high-dimensional remainder space. The low-dimensional space (to which we refer to as the LOD-space) is enriched with problem-specific features and admits a set of quasi-local basis functions that can be efficiently used in practical computations. In particular, due to the very good approximation properties and the quasi-locality of a basis, LOD-spaces can be easily used in a traditional Galerkin approach. Historically, the LOD was first introduced in [38] for elliptic multiscale problems with a rough coefficient. Since then, it has been further developed to suit a range of different problems, e.g., parabolic equations [36,35,37], equations describing wave phenomena, see, e.g., [1,41,28,21,34], non-linear problems [26,44], and more. In the context of nonlinear Schrödinger equations and GPEs, LOD techniques have been suggested in [27,31]. For a review of the LOD we refer to the textbook by Målqvist and Peterseim [39] and the recent survey article on numerical homogenization [3].
In this paper we revisit the LOD for the computation of ground states of the Gross-Pitaevskii equation (as introduced in [27]) in order to prove the higher order convergence rates that were numerically observed by the authors, but which have not yet been theoretically established. To be precise, the a priori error estimates proved in [27] for the GPE predict convergence rates of order O(H 2 ) for the H 1error, of order O(H 3 ) for the L 2 -error and the eigenvalue error and of order O(H 4 ) for the error in energy, where H is the "coarse" mesh size of the LOD-space. As already noted in [27], these rates seemed suboptimal when compared to numerical experiments. In fact, in this contribution we prove the optimal rates, which are of order O(H 3 ) for the H 1 -error, of order O(H 4 ) for the L 2 -error and (with a dramatic improvement) of order O(H 6 ) for the eigenvalue-error and the error in energy. For that we will not require any additional regularity assumptions and the results still hold for rough potentials V ∈ L 2+σ (Ω) with σ > 0.
Besides the obvious advantage of efficiently approximating ground states in the LOD space, our new results have direct implications for the simulation of the dynamics of Bose-Einstein condensates. In such a setting, the dynamics are described by nonlinear Schrödinger equations (NLS), where the initial values are typically ground states of the GPE (with respect to some modified configuration). In [31] it was recently proved that when the NLS is discretized with LOD spaces, then the energy is approximated and conserved with an accuracy of order O(H 6 ), provided that this accuracy can be already guaranteed for the initial value. Hence, this is exactly what we are establishing in this paper, as it justifies that the ground states computed in LOD-spaces can be straightforwardly used as initial values in an NLS without sacrificing the accuracy with which the energy is approximated over time.
Finally we note that in [27], an additional post-processing step on a finer mesh is suggested, which improves the quality of the approximation even further. We shall not focus on the post-processing in this paper. However, we note that the convergence orders can be further improved by such a technique.
The outline of the paper is as follows; in Section 2 the Gross-Pitaevskii equation is presented in more detail, in Section 3 we describe the LOD discretization technique, in Section 4 we present proofs of the optimal convergence rates in the H 1 -and L 2norm as well as for the energy and eigenvalue. Finally, in Section 5 we provide two numerical experiments that confirm the predicted convergence rates, both for smooth and discontinuous potentials.
We consider the Gross-Pitaevskii energy functional E : Here we make the following assumptions on the domain Ω, the trapping potential V and the interaction parameter β: (A1) Ω is a convex domain with polygonal boundary (for d = 1, 2, 3).
Observe that E is two times Fréchet-differentiable, where we denote the corresponding first and second derivatives by E and E respectively.
Under the above assumptions, we can seek the ground state of the energy E. The ground state is defined as a function u ∈ H 1 0 (Ω) that minimizes E(·) under the normalization constraint that Ω |u| 2 dx = 1, i.e.
The ground state is unique up to a sign and it can be equivalently characterized by an eigenvalue problem with eigenvector nonlinearity. The following result summarizes these analytical properties and a corresponding proof can be found in [12]. Furthermore, the unique positive minimizer u can be equivalently characterized by the Gross-Pitaevskii eigenvalue problem (GPE), which seeks the nonnegative eigenfunction u ∈ H 1 0 (Ω) with u = 1 and corresponding smallest eigenvalue λ > 0, such that The smallest eigenvalue λ and the ground state energy E(u) are connected through Note the nontrivial observation in Theorem 2.1 that if λ is the smallest eigenvalue of problem (2.3), then the corresponding eigenfunction u is a ground state of E. This means that minimizing the energy and minimizing the eigenvalue is equivalent.
We also note that (2.2) is a weak formulation that can be written using the Fréchet derivative of E. More precisely, with we observe that any eigenpair (λ, u) of (2.3) satisfies This formulation expresses that the eigenvalue problem should be seen as finding the critical points of the (even, positiv and convex) functional E on the L 2 -sphere S := {φ ∈ H 1 0 (Ω)| φ = 1}. Consequently, we also know that the spectrum of the differential operator is real and unbounded.
In the more general form of the Gross-Pitaevskii eigenvalue problem, where a rotational term is included, the eigenfunctions u are typically complex-valued, but the spectrum remains real. We write |u| 2 u instead of u 3 (which is the same for real functions), to work with the standard notation for the Gross-Pitaevskii equation that remains valid also in complex-valued settings.
Using the equivalence of the energy minimization perspective (2.2) and the GPE (2.3) (respectively (2.5)), we can easily prove additional regularity for the ground state u that will be useful in the convergence analysis. Lemma 2.2. Assume (A1)-(A3). Let u be a minimizer of (2.2) and let λ > 0 denote the corresponding ground state eigenvalue in the sense of (2.3), then u ∈ C 0 (Ω) ∩ H 2 (Ω) and |u| 2 u ∈ H 1 0 (Ω) ∩ H 2 (Ω) and we have the following bounds where C is a generic constant that only depends on Ω and d. Note that u H 1 (Ω) can be bounded again by √ λ.
Proof. In the following we write a b to abbreviate a ≤ C b for a constant C that only depends on Ω and d. In the first step we write If g ∈ L q/2 (Ω), for q > d, it follows that u ∈ C 0 (Ω), see e.g. [22,Theorem 8.30]. To verify that this condition is fulfilled, we define q such that 6q/(12 − q) = 2 + σ, where σ is as in (A2). For d = 1 and d = 2, this gives us q = 3 and consequently also q > d as desired. For d = 3, the condition q > d is equivalent to which is consistent with the assumption σ > 0 in (A2). Without loss of generality, we further assume σ ≤ 1 so that q ≤ 4 for d = 3. For q as introduced above, we need to prove now that g ∈ L q/2 (Ω). For that we use Hölder's inequality and the Sobolev embedding where we exploited u = 1. The term V L 2+σ (Ω) is bounded by (A2) and for the last term we get Consequently, we verified g ∈ L q/2 (Ω) for q > d and we have u ∈ C 0 (Ω). It is worth to mention that u ∈ C 0 (Ω) was already proved in [12], however, we could not verify the result under the stated regularity assumptions which slightly differ from ours. Therefore we decided to explicitly elaborate the argument here. Now that we have u ∈ C 0 (Ω), we can use this result to show that we indeed even have g ∈ L 2 (Ω), as . Using standard elliptic regularity theory for equation (2.8) we conclude that u ∈ H 2 (Ω) with the estimate . Hence, the bound in (2.6) follows.
For (2.7) note that where the multi-index α denotes the partial derivatives. For |α| = 0 the term reduces to (2.9). For the first order derivatives, i.e., for |α| = 1, we get where Hölder's inequality and the Sobolev embedding gives . For the second order derivatives, i.e. for |α| = 2, we have where α 1 , α 2 denotes first order derivatives such that α = α 1 +α 2 . For the first term (2.11) we have using Hölder's inequality and the Sobolev embedding once again . For the second term in (2.11) we get . This gives . Combining all the terms and using u H 1 (Ω) ≤ u H 2 (Ω) we achieve (2.7).

Discretization
Approximations of the ground state can be found by minimizing the energy E over a discrete space V δ subject to the L 2 -normalization constraint. The following abstract approximation result was established in [12, Theorem 1 and Remark 2] and will be a central ingredient for our proofs.
Theorem 3.1. Assume (A1)-(A3) and let u be the ground state of (2.2) (up to sign). Furthermore, let (V δ ) δ>0 be a family of finite-dimensional subspaces of H 1 0 (Ω) with the property that with the property (u, u δ ) ≥ 0. Then there exist (problem-dependent) generic constants C, c 1 , c 2 > 0 such that i.e., u δ is a H 1 -quasi best approximation to u; and it also holds Finally, for the L 2 -error it holds where for ω ∈ H −1 (Ω), the function The results established in [12] also contain an abstract approximation result for the eigenvalue, which is however not optimal for all discrete spaces and needs to be therefore revisited individually depending on the considered discrete setting.
Following [27], our next goal is the construction of a suitable discrete space that allows for very high convergence orders under minimal regularity assumptions. For that we start from a conventional conforming finite element discretization based on a (coarse) quasi-uniform triangulation T H of the domain Ω. Here H denotes the mesh size parameter, i.e. H = diam(T ), for T ∈ T H . With this, let V H be a classical P 1 finite element space on T H with The classical finite element method seeks an approximation u H ∈ V H with u H = 1, and where the corresponding eigenvalue is given by . Since E is weakly lower semi-continuous on H 1 0 (Ω) and bounded from below by zero, discrete minimizers always exist on finite dimensional spaces. However, it is typically not clear if the solution to (3.3) is also unique up to sign, see e.g. the discussion in [12]. Independent of that uniqueness, it can be shown that for any minimizer u H with (u H , u) ≥ 0 (i.e. consistent with the sign of u) it holds where u is the unique (up to sign) ground state and C is a constant independent of H. The space that we construct in the next subsection will have the same dimension as V H , but it will be able to boost the convergence rates to 3rd order for the H 1 -error, 4th order for the L 2 -error and 6th order for eigenvalue and energy error.
3.1. Localized Orthogonal Decomposition. In the following we summarize the technique proposed in [27] to make the paper self-contained. Let P H : H 1 0 (Ω) → V H denote the classical L 2 -projection onto V H and let V f := ker(P H ) = {v ∈ H 1 0 (Ω)| P H (v) = 0}, be the kernel of the projection. Due to H 1 -stability of the L 2 -projection on quasiuniform meshes (cf. [8]), the space V f is a closed subspace of H 1 0 (Ω).
We expect V f to contain fine scale details from H 1 0 (Ω) that are not captured by the coarse space V H . With this definition we obtain the following L 2 -orthogonal splitting Next, we define an inner product a(·, ·) on H 1 0 (Ω), that is based on the linear terms in the GPE (2.3), by a(v, w) := (∇v, ∇w) + (V v, w).

This inner product induces another orthogonal splitting, namely
respect to a(·, ·). Note that, by construction, the V LOD space has the same dimension as V H . This is precisely the space in which we will look for discrete minimizers. This is fixed in the following definition.
Definition 3.2 (LOD ground state approximation). The LOD ground state approximation is obtained by using the V LOD space in (2.2). That is, we seek u LOD ∈ V LOD such that (u LOD , u) ≥ 0, u LOD = 1, and where the corresponding discrete eigenvalue is given by Remark 3.3 (Practical realization). The practical realization of (3.4) requires some remarks. First of all, the condition (u LOD , u) ≥ 0 is never enforced in the numerical method, as it does not matter if u LOD converges to the positive or the negative ground state, the respective approximation properties are identical. Second, for the construction of V LOD it is necessary to solve a set of saddle point problems that determine the basis functions of V LOD . For an efficient realization, this includes a localization step which is analytically justified by an exponential decay of the basis function in units of H (cf. [29,38]). Furthermore, the space V f needs a discrete representation on a finer mesh of size h < H. Comprehensive details on all these implementation details are given in [20].
The influence of these additional approximations (i.e. truncation/localization and fine scale representation of the LOD basis functions) has been studied extensively in the literature and we note that they do not lead to a reduction of the convergence orders of the final method when chosen appropriately. For the Gross-Pitaevskii eigenvalue problem, the influence has been analytically studied in [27].
Further practical tricks for an efficient treatment of the nonlinear term are discussed in [31].
Finally we also stress that aside from the construction of the LOD space, a practical realization also requires an iterative method for finding a discrete minimizer. Here various techniques are possible, where we refer to the introduction for a corresponding literature survey.
In the light of Remark 3.3, all the proofs below are based on the ideal setting of Definition 3.2, i.e. the setting without localization and fine scale representation. This simplifies both the arguments and the notation. By considering localization, an additional term depending on the size of the localization patches would appear in the error estimates. The implications of truncation and fine scale discretization are already thoroughly discussed in [27] for the GPE, which is why we omit this part in our paper. This allows us to keep the presentation short and to focus on the main novelty of this contribution, namely the arguments that actually improve the convergence orders in the proofs.
We have the following main result, which we shall prove in Section 4.
where C depends on V , λ, β, Ω, u L ∞ (Ω) and the mesh regularity of T H .
Compared to the original result proved in [27], the rates for the L 2 and H 1 -error are both improved by one order each, the rate for the energy by two orders and the rate for the eigenvalue is dramatically improved by three orders. However, we also stress that the equation considered in [27] is slightly more general, involving another (potentially discontinuous) coefficient in the kinetic term. Therefore the rates in [27] are optimal for this generalized setting, but only suboptimal for the GPE. Finally, we also note that the rates in (3.5) can be even further improved significantly by using a post-processing technique as in [27].

Proof of Optimal Convergence Rates
In this section we prove the optimal convergence rates of order O(H 3 ) in the H 1norm, O(H 4 ) in the L 2 -norm, and O(H 6 ) for both the energy and the eigenvalue. The error bounds presented below all depend on some constant C, which may depend on the ground state eigenvalue λ, the potential V , the interaction constant β, and the domain Ω, but not on the mesh size H.
In the analysis below we shall make use of the a(·, ·)-orthogonal projection P LOD : The following approximation properties of the LOD space are well known and are explicitly stated and proved in [31, Section 2.1].
The constant C only depends on Ω and the mesh regularity of T H .

4.1.
Convergence in H 1 -norm. We start with the H 1 -error estimate which follows from the abstract H 1 -best approximation result in Theorem 3.1, combined with the approximation properties in Lemma 4.1 and the regularity result in Lemma 2.2. We have: Proof. We note that

4.2.
Convergence in L 2 -norm. Next, we turn our attention to proving the L 2error estimate which reads as follows. The proof is similar to [27,Section 6.3] with several simplifying modifications that keep the proof short.
We analyze the terms on the right hand side one by one. Hölder's inequality and the Sobolev inequality u L 6 (Ω) ≤ C u H 1 (Ω) (for d ≤ 3) give us for the first term For the second term we readily have λψ eLOD ≤ Cλ e LOD . For the third term we use (2.9) and u = 1 to get With u = 1, the last term is estimated as Combining the estimates for all terms yields

4.4.
Convergence of Eigenvalues. The optimal convergence rates for the eigenvalue involve the most substantial changes compared to [27]. In the first step, we prove the following auxiliary result, which established optimal convergence rates in the H −2 -norm.
We conclude with (4.3) that For II 2 we use Lemma 2.2 For the last terms we use Lemma 2.2 again together with the regularity bound ψ ω H 2 (Ω) ≤ C ω to conclude We deduce II ≤ CH 6 ω H 2 (Ω) .
For the third and fourth term in (4.9) we use the error bound in the L 2 -norm from Proposition 4.3, the regularity of ψ ω , and Lemma 2.2 to derive Combining everything, we obtain (ω, u − u LOD ) ≤ CH 6 ω H 2 (Ω) , which gives With this we are ready to prove the final estimate for the eigenvalue, which also finishes the proof of Theorem 3.4. Proposition 4.6. Let λ and λ LOD be the corresponding eigenvalues to u and u LOD in (2.4) and (3.4) respectively. It holds that |λ − λ LOD | ≤ CH 6 .
Proof. We use the definition of the eigenvalues to estimate where the first term is bounded by Conclusion 4.4. It remains to bound the second term. We have For the first term we use Proposition 4.2 to achieve The second term splits into . We note that by Hölder's inequality and Sobolev embeddings we have the bound This relies on the fact that u LOD H 1 (Ω) ≤ C, which can simply be deduced by using the error bound in H 1 -norm Thus, using Lemma 4.5 and Proposition 4.2 we deduce II ≤ CH 6 , which completes the proof.

Numerical Experiments
In this section we perform two numerical experiments to validate the theoretical convergence rates for both a smooth and a discontinuous potential V . The numerical experiments are performed for the ideal method, without localization of the corrections, to eliminate any potential reduction in convergence rate coming from the size of the localization patches. We emphasize that these experiments are only to verify that the derived rates are indeed optimal for the Gross-Pitaevskii equation and in a practical application the localization should always be implemented. For further details regarding localization we refer to [27] where the method was first suggested.

Smooth potential.
For the smooth potential we consider a setup similar to [9]. The computational domain is given by the square Ω = [−6, 6] × [−6, 6], the potential is set to the smooth (harmonic) function V (x, y) = 0.5(x 2 + y 2 ), and the interaction parameter to β = 100.
The resulting minimizing problem (3.4) is solved by using the iterative scheme, referred to as normalized gradient flow, presented in [9].
A reference solution is computed on mesh a fine mesh of size h = 2 −4 . The LOD approximations u LOD are computed for meshes of decreasing size H = 2 1 , 2 0 , 2 −1 , 2 −2 and then compared to the reference solution. The resulting convergence rates for the H 1 -norm, the L 2 -norm, the energy, and the eigenvalue are plotted in Figure 1. The plots clearly confirm the convergence orders predicted by Theorem 3.4.
The potential is set to a checkerboard with squares of size 2 −2 which results in 24 × 24 squares on the domain. The values at the squares alternate between 0 and 1, see Figure 2.
A reference solution is computed on a mesh of size h = 2 −5 . The triangles in the mesh are aligned with the discontinuities in the potential, so that the potential can be computed exactly on each triangle to avoid additional numerical errors. The LOD approximations u LOD are computed for meshes of decreasing size    H = 2 0 , 2 −1 , 2 −2 , 2 −3 and then compared to the reference solution. The resulting convergence rates for the H 1 -norm, the L 2 -norm, the energy, and the eigenvalue are plotted in Figure 3. Again, the plots clearly confirm the predicted convergence orders.