extfieldsofstudy:
Computer Science
Engineering
Mathematics

text:
Multi-Entity and Multi-Enrollment Key Agreement with Correlated Noise

A basic model for key agreement with a remote (or hidden) source is extended to a multi-user model with joint secrecy and privacy constraints over all entities that do not trust each other. Multiple entities using different measurements of the same source through broadcast channels (BCs) to agree on mutually-independent local secret keys are considered. Our model is the proper multi-user extension of the basic model since the encoder and decoder pairs are not assumed to trust other pairs, unlike assumed in the literature. Strong secrecy constraints imposed on all secret keys jointly, which is more stringent than separate secrecy leakage constraints for each secret key considered in the literature, are satisfied. Inner bounds for maximum key rate, and minimum privacy-leakage and database-storage rates are proposed for any finite number of entities. Inner and outer bounds for degraded and less-noisy BCs are given to illustrate cases with strong privacy. A multi-enrollment model that is used for common physical unclonable functions is considered to establish inner and outer bounds for key-leakage-storage regions that differ only in the Markov chains imposed. For this special case, the encoder and decoder measurement channels have the same channel transition matrix and secrecy leakage is measured for each secret key separately. We illustrate cases for which it is useful to have multiple enrollments as compared to a single enrollment and vice versa.

I. INTRODUCTION
A natural source of randomness is biometric identifiers such as fingerprints that are generally transformed into a frequency domain and quantized to obtain bit sequences that are unique to an individual [1]. Similarly, physical identifiers such as fine variations of ring oscillator (RO) outputs or random startup values of static random access memories (SRAMs) that are caused by uncontrollable manufacturing variations, are safer and cheaper alternatives to key storage in a non-volatile memory [2]. Physical identifiers for digital devices, such as Internet-of-Things (IoT) devices, can be implemented using physical unclonable functions (PUFs) [2]. One can use PUFs in various coding schemes as a source of local randomness [3, Chapter 1], e.g., in the randomized encoder of the wiretap channel [4] and of the strong coordination problem [5], [6].
We use the basic source model for key agreement from [7], [8] to find achievable rate regions for key agreement with PUFs and biometric identifiers. In this classic model, an encoder observes a source output to generate a secret key and O. Günlü is with the Information Theory and Applications Chair, Technische Universität Berlin, 10623 Berlin, Germany (e-mail: guenlue@tuberlin.de). sends public side information, i.e., helper data, to a decoder, so the decoder can reliably reconstruct the same secret key by observing another source output and the helper data. The main constraints are that the information leaked about the secret key, i.e., secrecy leakage, is negligible and the information leaked about the identifier output, i.e., privacy leakage, is small [9], [10]. Furthermore, the amount of storage should also be minimized to limit the hardware cost [11].
Suppose the encoder generates a key from a noisy measurement of a hidden (or remote) source output, and a decoder has access to another noisy measurement of the same source and the helper data to reconstruct the same key. We call this model the generated-secret (GS) model with a hidden source. This model is introduced in [12] as an extension of the visible (noiseless) source outputs observed by the encoder, considered in [9], [10]. Similarly, for the chosen-secret (CS) model, an embedded (or chosen) key and noisy identifier measurements are combined by the encoder to generate the public helper data. We consider both models to address different applications.

A. Related Work and Motivation
The same identifier is used by multiple encoder and decoder pairs in [13], where the identifier outputs observed by different encoders are the same because the encoder measurements are assumed to be noiseless. Therefore, the multiple use of the same noiseless source output allows all encoders to know the secret key of the other encoders. This model does not fit well to the practical key agreement with identifiers scenarios because there is noise in every identifier measurement.
Multiple enrollments of a hidden source using noisy measurements are considered in [14], where weakly secure secret keys are generated without privacy leakage and storage constraints. Furthermore, there is a causality assumption in [14] on the availability of the helper data, i.e., any decoder has access to all previously-generated helper data. This assumption is not necessarily realistic as a decoder of, e.g., an IoT device that embodies a PUF should be low complexity and the amount of data to process increases linearly with the number of enrollments. In addition, any manipulation in any of the helper data can cause the complete multi-enrollment system to fail.
A classic method used for key agreement, i.e., the fuzzy commitment scheme (FCS) [15], is used in [16] in combination with an SRAM PUF to enroll the noisy outputs of the same SRAM multiple times. The symmetry condition in [16,Eq. (16)] conditioned on a fixed SRAM cell state is entirely similar to the symmetry satisfied by binary-input symmetric output (BISO) channels; see e.g., [17, p. 613], [12,Eq. (14)]. For SRAM outputs that satisfy this symmetry, the normalized (weak) secrecy leakage about each separate secret key is shown to be zero. It is discussed in [18,Section 3.4] that any uniformly-distributed hidden identifier output with BISO measurement channels satisfies the results in [16]. In [18, Theorem 1] the secret-key capacity of the two-enrollment key agreement problem is established for measurement channels with the same channel transition matrix. However, these multienrollment models do not consider the privacy leakage and storage constraints, there is no constraint on the independence of the secret keys of different enrollments, and the secrecy leakage constraint is weak and is not applied jointly on all secret keys. Furthermore, optimal random linear code constructions that achieve the boundaries of the key-leakage-storage regions are given in [19], where the classic code constructions FCS and code-offset fuzzy extractors [20] are shown to be strictly suboptimal. Therefore, the multi-enrollment models and constructions in the literature are strictly suboptimal and not necessarily realistic. We therefore list stronger secrecy constraints jointly on all entities, which approximates the reality better in combination with storage rate and joint privacyleakage rate constraints. These constraints define the multientity key agreement problem, where the entities that use the same identifier do not have to trust other entities. Therefore, the multi-entity key agreement problem is a proper multi-user extension of single-enrollment models.
Every measurement of an identifier is considered to be noisy due to, e.g., local temperature and voltage changes in the hardware of the PUF circuit or a cut on the finger. Noise components at the encoder and decoder measurements of a hidden source can be also correlated due to, e.g., the surrounding logic in the hardware [21] or constant fingertip moisture. This correlation between the noise sequences is modeled in [22] as a broadcast channel (BC) [23] with an input that is the hidden source output and with outputs that are the noisy encoder and decoder measurements. We use this model for multi-entity key agreement with identifiers, where each entity (i.e., each encoder and decoder pair) observes noisy identifier outputs of the same hidden source through different BCs. We allow the BCs to be different as trusted entities generally use different hardware implementations of the encoder and decoder pairs, which results in different correlations between noise components.
We also consider physically-degraded (PD) and less-noisy (LN) BCs to give finer inner and outer bounds to the keyleakage-storage regions for the GS and CS models of the multi-entity key agreement problem. For the considered PD and LN BCs, we prove that strong privacy can be achieved. In [9], [10], [24], an extra common randomness that is available to the encoder and decoder and that is hidden from the eavesdropper is required to obtain strong privacy. This assumption is not realistic since such a common randomness requires hardware protection against invasive attacks, and if such a protection is feasible, then it is not necessary to use an identifier for key agreement.

B. Models for Identifier Outputs
We study physical and biometric identifier outputs that are independent and identically distributed (i.i.d.) according to a given probability distribution. These models are reasonable if one uses transform-coding algorithms from [25] that occupy a small hardware area to extract almost i.i.d. bits from PUFs under varying environmental conditions. Similar transformcoding based algorithms have been applied to biometric identifiers to obtain independent output symbols [26]. These transform-coding algorithms provide almost i.i.d. identifier outputs and noise sequences; however, the correlation between the noise components on the encoder and decoder components are not removed using these methods. Furthermore, PUFs are used for on-demand key reconstruction and physical attacks on PUFs permanently change the identifier outputs [27], so we assume that the eavesdropper cannot obtain information correlated with the PUF outputs, unlike biometric identifiers.

C. Summary of Contributions
We extend the key-leakage-storage rate tuple analysis of the single-enrollment model for hidden identifier outputs measured through general BCs in [22] to consider multi-entity key agreement with a set of stringent secrecy constraints. A summary of the main contributions is as follows.
• We derive achievable key-leakage-storage rate tuples for the GS model with strong secrecy for any finite number of entities using the same identifier's measurements through different BCs for key agreement. Separate identifier measurements considered in [12], [28] correspond to a PD BC and the visible source model in [9], [10] corresponds to a semi-deterministic BC. • For a set of PD and LN BCs, the privacy-leakage rates for the two-entity key agreement problem are calculated. These PD and LN BCs are shown to provide strong privacy without the need of a common randomness. An outer bound is given for the considered PD and LN BCs. • We consider a special case of the multi-enrollment key agreement problem, where all measurement channels are separate (i.e., PD BCs) and they have the same transition matrix. This is a common model used for SRAM PUFs [16]. Using a less stringent secrecy leakage constraint that bounds the information leakage for each secret key separately and without the mutual independence constraint on the secret keys, we establish inner and outer bounds for the strong-secrecy key-leakage-storage region for this two-enrollment key agreement problem. The bounds differ only in the Markov chains imposed. This result is a significant improvement to the two-enrollment secretkey rate region (without storage and privacy-leakage rate constraints) established in [18] for weak secrecy, which is recovered by eliminating auxiliary random variables in the proposed rate regions. • All inner and outer bounds for the GS model are extended to the CS model, which comprises secret-key binding methods that embed a chosen secret key to the encoder. • We give two scenarios to compare single-enrollment and two-enrollment models and illustrate that for different assumptions on measurement channels, either of the two models can perform better in terms of the privacy-leakage vs. secret-key rate boundary tuples.

D. Organization
This paper is organized as follows. In Section II, we describe the multi-entity key agreement problem with BC measurements. We give achievable key-leakage-storage regions for the GS and CS models with strong secrecy and BC measurements for any finite number of entities in Section III in addition to inner and outer bounds for PD and LN BCs that satisfy strong privacy. The proposed inner bounds for the two-entity key agreement problem in Section IV are shown to differ from the outer bounds only in the Markov chains imposed for a special case with less stringent secrecy constraints. In Sections V and VI, proofs of the given rate regions for the general multi-entity key agremeent problem and for the twoenrollment key agreement problem, respectively, are given. Section VII concludes the paper.

E. Notation
Upper case letters represent random variables and lower case letters their realizations. A superscript denotes a string of variables, e.g., X n = X 1 , X 2 , . . . , X i , . . . , X n , and a subscript i denotes the position of a variable in a string. A random variable X has probability distribution P X . Calligraphic letters such as X denote sets, set sizes are written as |X | and their complements as X c .
is the binary entropy function, where we take logarithms to the base 2, and H −1 b (·) denotes its inverse with range [0, 0.5]. X ∼ Bern(α) is a binary random variable with Pr[X = 1] = α. A binary symmetric channel (BSC) with crossover probability p is denoted by BSC(p). Q(·) is the Q-function that gives the tail probability for the standard normal distribution.

II. PROBLEM DEFINITIONS
Consider hidden identifier outputs X n that are i.i.d. according to a probability distribution P X . If there are a finite number J of trusted entities that use the same identifer, the j-th encoder and decoder pair observes noisy source measurements that are outputs of a BC P Xj Yj |X , with abuse of notation, for j ∈ [1 : J], where X j , Y j , and X are finite sets.
For the GS model shown in Fig. 1(a) for J = 2 trusted entities, the j-th encoder f GS,j (·) generates helper data W j and a secret key S j from X n j . All secret keys are stored in a secure database, whereas helper data are stored in a public database so that an eavesdropper has access only to the helper data. Using the helper data W j and the noisy measurements Y n j , the j-th decoder g j (·, ·) generates the key estimateŜ j . Similar steps are applied for the CS model in Fig. 1(b) also for J = 2 trusted entities, except that each S j should be embedded into the j-th encoder f CS,j (·, ·).
Denote a set of secret keys as and a set of helper data as Similarly, we denote a set of secret-key rates, for any K ⊆ [1 : J], as and a set of storage rates as We next define the key-leakage-storage regions for the GS and CS models.
is achievable for the GS and CS models with j-th encoder and decoder measurements through a BC P Xj Yj |X if, given any δ > 0, there is some n ≥ 1, and J encoder and decoder pairs for which R s,j = log |S j | n for all j ∈ [1 : J] and The key-leakage-storage regions C gs for the GS model and C cs for the CS model are the closures of the set of all achievable rate tuples (R s,[1:J] , R , R w, [1:J] ).
The mutual key independence constraint in (7) is not imposed in the multi-enrollment key agreement problem considered in [16]. Furthermore, a normalized (weak) version of this constraint is imposed in the multi-enrollment key agreement problem considered in [14], where the j-th decoder g j (·, ·) is assumed to have access to the set of helper data W [1:j] for all j ∈ [1 : J]. The lack of the mutual key independence constraint and the assumption of availability of all previous helper data require that different encoder and decoder pairs should trust each other. This can be the case, e.g., if all enrollments are made by the same entity. Therefore, the multientity key agreement problem imposes strictly more stringent constraints than the multi-enrollment key agreement problem.
The unnormalized privacy leakage I(X n ; W [1:J] ) cannot be bounded by a finite number in general. We illustrate special strong privacy cases in the next section.

III. INNER BOUNDS
We are interested in characterizing the optimal trade-off among the secret-key, privacy-leakage, and storage rates with strong secrecy for BC measurements at the encoders and decoders of any finite number J of entities that use the same hidden identifier outputs. We give achievable rate regions for the GS and CS models in Theorem 1. The proofs are given in Section V. Denote and define a function max{·, ·} that gives the maximum of the input values as its output.
Theorem 1 (Inner Bounds for Multi-entity GS and CS Models). An achievable rate region R gs for the GS model with J entities is the union over all P Uj | Xj for all j ∈ [1 : J] of the rate tuples such that R s,j ≥ 0 for all j ∈ [1 : J] and An achievable rate region R cs for the CS model with J entities is the union over all P Uj | Xj for all j ∈ [1 : J] of the rate tuples such that R s,j ≥ 0 for all j ∈ [1 : J], (12), (13), and For the achievable rate regions R gs and R cs , we have For these cases, strong privacy, i.e., can be achieved for the GS and CS models in combination with the other corresponding bounds given in Theorem 1.
The proof of Corollary 1 follows from Theorem 1 because I(U j ; X) − I(U j ; Y j ) ≤ 0 for all j ∈ [1 : J] for BCs considered in Corollary 1.
Corollary 1 illustrates that it is possible to obtain strong privacy, i.e., negligible unnormalized privacy leakage, without the requirement of a common randommness that is hidden from an eavesdropper assumed in [9], [10], [24]. This is the case because the observation Y n j of each decoder is "better" than the observation X n j of the corresponding encoder with respect to the hidden source X n for all entities. Remark 1. The rate regions for our problem depend on the joint conditional probability distributions P XYj | Xj rather than only the marginal conditional distributions. Thus, the keyleakage-storage regions for the stochastically-degraded BCs are not necessarily equal to the regions for the corresponding PD BCs, unlike in the classic BC problem. Furthermore, since P X [1:J] XY [1:J] is fixed, the distinction between the LN BCs and essentially-less noisy BCs [29], is not necessary.
We next give simple outer bounds for the key-leakagestorage regions C gs for the GS model and C cs for the CS model when the BCs P XYj | Xj for all j ∈ [1 : J] are PD BCs or LN BCs as defined in Corollary 1. These simple outer bounds give insights into the reason for different bounds on the secret-key rates. Based on these insights, we show a special case in the next section with a less stringent secrecy constraint, for which the inner and outer bounds differ only in the Markov chains imposed and we illustrate that they match for simpler models.
Lemma 1. Suppose one of the cases given in Corollary 1 is satisfied by the BCs P XYj | Xj for all j ∈ [1 : J]. An outer bound on the key-leakage-storage region C gs is the union over (14), (19), and An outer bound to the key-leakage-storage region C cs for the same BCs P XYj | Xj is the union over all (19), (20), and The proof of Lemma 1 follows straightforwadly by following the steps in [12, Section VI], defining the auxiliary random variables U j,i = (S j , W j , Y i−1 j ) for all j ∈ [1 : J] and i ∈ [1 : n], and by bounding I(X n ; W [1:J] ) ≥ 0; therefore, we omit the proof.
The outer bounds do not include the inequalities in (15) and (17). Furthermore, the secret-key rate achieved by the inner bound in (12) is smaller than the outer bound given in (20), where the difference is the term −I(U j ; U [1:J]\{j} ). This term is a result of the constraint in (44) that is imposed to satisfy the strong and mutual key independence constraint given in (7). Therefore, we next consider a model without the constraint in (7) and use a secrecy-leakage constraint that is less stringent than the one in (9), i.e., replace (9) by which is also a strong secrecy metric. Due to the lack of a mutual key independence constraint, the model in the next section is not a multi-entity model but rather a multienrollment model. For a special case of this multi-enrollment key agreement problem, we establish inner and outer bounds for the key-leakage-storage regions that comprise the same bounds but for different Markov chains.

IV. BOUNDS FOR A MULTI-ENROLLMENT MODEL
Consider J = 2 entities that measure the noisy outputs of the same hidden source X n through separate channels that have the same channel transition matrices, i.e., for all j ∈ [1 : 2],x j ∈ X , and y j ∈ X we have P Xj Yj |X (x j , y j |x) = P X|X (x j |x)P X|X (y j |x).
Using (23), we define a new multi-enrollment model.

Definition 2.
A key-leakage-storage rate tuple ( s R s,1 , s R s,2 , s R , s R w,1 , s R w,2 ) is achievable for the GS and CS models with measurements through a BC P XY |X (x, y|x) as in (23) if, given any δ > 0, there is some n ≥ 1, and two encoder and decoder pairs for which s R s,1 = log |S 1 | n , 1 n H(S j ) = s R s,j − δ, j = 1, 2 (key uniformity) (25) The key-leakage-storage regions s C gs,J=2 for the GS model and s C cs,J=2 for the CS model are the closures of the set of all achievable rate tuples ( s We characterize in Theorem 2 inner and outer bounds for key-leakage-storage regions s C gs,J=2 and s C cs,J=2 . The proofs of Theorem 2 are given in Section VI, where the reason for the necessity of the secrecy-leakage constraint in (27) that is less stringent than the joint secrecy-leakage constraint in (9) is given in Remark 2. Similarly, the reason for the necessity of the strong helper data (storage) independence constraint in (29) is discussed in Remark 4. Denote Theorem 2. (Inner Bounds): An achievable key-leakagestorage region Ď R gs,J=2 is the union over all P U1| X1 and P U2| X2 of the rate tuples such that s R s,j ≥ 0 for j = 1, 2 and An achievable key-leakage-storage region Ď R cs,J=2 is the union over all P U1| X1 and P U2| X2 of the rate tuples such that s R s,j ≥ 0 for j = 1, 2, (31)-(33), and For both achievable rate regions Ď R gs,J=2 and Ď R cs,J=2 , we have P U1U2 X1 X2XY1Y2 (u 1 , u 2 , x 1 , x 2 , x, y 1 , y 2 ) = P U1| X1 (u 1 |x 1 )P U2| X2 (u 2 |x 2 )P X|X (x 1 |x)P X|X (x 2 |x) × P X|X (y 1 |x)P X|X (y 2 |x)P X (x).
(Outer Bounds) An outer bound for the key-leakage-storage region s C gs,J=2 is the union over all P U1| X1 and P U2| X2 of the rate tuples such that s R s,j ≥ 0, (31) - (36), and U j − X j − X − Y j form a Markov chain for j = 1, 2.
An outer bound for the key-leakage-storage region s C cs,J=2 is the union over all P U1| X1 and P U2| X2 of the rate tuples such that s R s,j ≥ 0, (31) -(33), (37) -(39), and U j − X j − X − Y j form a Markov chain for j = 1, 2.
The inner and outer bounds differ because the outer bounds define rate regions for the Markov chains U 1 − X 1 − X − Y 1 and U 2 − X 2 − X − Y 2 , which are larger than the rate regions defined by the inner bounds that satisfy (40). For instance, in the achievability proof of Theorem 2, we apply the properties of the Markov chain U 2 − X 2 − U 1 in (87)(b), which does not form a Markov chain for the choice of U 1 and U 2 in the outer bounds. Therefore, inner and outer bounds do not match in general.
Corollary 2. Choosing U 1 = X 1 and U 2 = X 2 , it is straightforward to show that inner and outer bounds in Theorem 2 match if we do not impose any storage or privacy constraints, i.e., impose only (24), (25), and (27). This result improves on the secret-key capacity region given in [18, Theorem 1] for a weak secrecy constraint.

Example 1.
Consider the RO PUF model from [25, Section 4.1] where a transform-coding method is applied to conservatively model the measurement channels P Y |X = P X|X as independent BSCs with the same crossover probability of p A and where the hidden source output is Bern( 1 2 ). We therefore can apply the achievability results from Theorem 2 to this RO PUF model. Using [12,Theorem 3] to evaluate the boundary tuples of Ď R gs,J=2 , it suffices to consider probability distributions P Uj | Xj for j = 1, 2 such that P Xj |Uj are BSCs with crossover probabilities Consider the projection of the boundary tuples of Ď R gs,J=2 onto key-leakage plane, i.e., (31) and (32). We plot in Fig. 2 singleenrollment results where the privacy-leakage rate is measured with respect to single helper data and two-enrollment results for the sum rate of the two keys, both for p A = 0.06 [25]. To achieve a total secret-key rate of I( X 1 ; Y 1 ) = I( X 2 ; Y 2 ), the privacy-leakage rate for the two-enrollment model is approximately 13.5% less than the privacy-leakage rate for the single-enrollment model for RO PUFs. The reason for this gain is the information bottleneck problem that arises from (31) and (32) to find the boundary tuples.  Define the signal power as P S and the noise power as P N , so we have a signal-to-noise ratio (SNR) of SN R = P S P N . If a matched filter, which maximizes the SN R at the sampling instant for the AWGN channel, is applied at the encoder and decoder, the bit error probability P b is given by [31, pp. 96] The channel between the input binary symbols and the outputs of the matched filter is a BISO channel. Using [12, Theorem 3], we have that P Xj |Uj for j = 1, 2 that are BSCs with crossover probabilities given in (41) by replacing p A with P b , suffice to obtain the boundary tuples of Ď R gs,J=2 . We remark that p A = 0.06 used in Example 1 corresponds to an SNR of approximately 3.83dB.
In Fig. 3, the privacy-leakage rate vs. secret-key rate boundary tuples are depicted for two cases. First, a two-enrollment model at SN R = 3.83dB with a sum rate for secret keys is depicted. For comparison, we plot a single-enrollment model where the signal power P s is doubled, i.e., we have an SNR of approximately 6.84dB. Fig. 3 shows for the two cases with the same total signal power of 2P s that the single enrollment boundary tuple can result in a gain of approximately 228.55% at the corner point in terms of the secret-key rates achieved for a given privacyleakage rate. For such an AWGN channel with a fixed total signal power; therefore, the single-enrollment model can result in significant gains in terms of achieved secret-key rates as compared to the two-enrollment model for small s R values.
V. PROOF OF THEOREM 1 We provide a proof that follows from the output statistics of random binning (OSRB) method, proposed in [32] and further extended in [33], by applying the steps in [ (18). Assign three random bin indices (S j , W j , C j ) to each realization u n j for all j ∈ [1 : J], where S j represents the secret key, W j the helper data, and C j a public index. Assume S j ∈ [1 : 2 nRs,j ], W j ∈ [1 : 2 nRw,j ], and C j ∈ [1 : 2 nRc,j ] such that R s,j , R w,j , R c,j ≥ 0 for all j ∈ [1 : J].
Apply the union bound to the reliability constraint in (5) to obtain the sum of J error probabilities. This sum vanishes for any finite number J when n → ∞ by using a Slepian-Wolf (SW) [35] decoder to estimate U n j from (C j , W j , Y n j ) if [32, The key uniformity (6), mutual key independence (7), and strong secrecy (9)  Similarly, the public randomness C j is almost independent of X n j , so it is almost independent of ( X n Thus, the public indices C [1:J] can be fixed by generating them uniformly at random. The j-th encoder can generate U n j according to P U n j | X n j Cj obtained from the binning scheme above to compute the bins (S j , W j ) from U n j for all j ∈ [1 : J]. This procedure induces a joint probability distribution that is almost equal to P U [1:J] X [1:J] XY [1:J] fixed in (18) [34,Section 1.6].
To satisfy the constraints (46)-(48), we can fix the rates to where (a) follows from the fact that C j is almost independent of ( X n [1:J] , X n , Y n [1:J] ) for all j ∈ [1 : J] due to (45). Therefore, C j is almost independent of (X n , C [1:J]\{j} ) because U n k determines C k and U n k − X n k − X n j − U n j form a Markov chain for all k = j and j, k ∈ [1 : J], so we have I(X n ; C [1:J] ) = I(X n ; C 1 ) + I(X n ; C 2 |C 1 ) + . . . + I(X n ; C J |C [1: for some n > 0 such that n → 0 when n → ∞; i.e., H(U j |X) > H(U j |Y j ), so (W j , C j , X n ) are almost mutually independent [32, Theorem 1]. Therefore, we have for some n > 0 such that n → 0 when n → ∞. Combining (52) and (55) proves strong privacy.
where (a) follows because U n j determines (W j , C j ), (W j , C j , X n ) can realiably estimate U n for some n > 0 such that n → 0 when n → ∞, and (U n j , X n ) are i.i.d., and (b) follows by (50) and (51).
Combining (52) and (57), we obtain These prove the achievability of the rate region R gs .

B. Proof for the CS Model
We use the achievability proof for the GS model. Suppose the key S j , generated as in the GS model together with the helper data W j and public index C j , have the same cardinality as the corresponding embedded secret key S j , i.e., |S j | = |S j | for all j ∈ [1 : J]. The chosen key S j is uniformly distributed and independent of (X n , X n Consider the j-th encoder f cs,j (·, ·) with inputs ( X n j , S j ) and output W j = (S j + S j , W j ), and the j-th decoder g j (·, ·) with inputs (Y n j , W j ) and outputŜ j = S j + S j −Ŝ j . All addition and subtraction operations are modulo-|S j | for all j ∈ [1 : J]. The j-th decoder of the GS model is used to obtainŜ j for all j ∈ [1 : J].
We have the error probability which is small due to the proof of achievability for the GS model. Using (49) and (50), and from the one-time padding operation applied above, we can achieve a storage rate of are almost mutually independent and each S j is almost uniformly distributed due to (44) for some n > 0 such that n → 0 when n → ∞, and (e) follows because the GS model satisfies the strong secrecy constraint (9) due to (44) for some (4) n > 0 such that (4) n → 0 when n → ∞. Using the selection lemma, there exists a binning that achieves all rate tuples given in (49), (51), and (60). Consider the privacy leakage for this binning. This proves the achievability of the rate region R cs .

A. Achievability Proof for the GS Model
Fix Let (U n 1 , U n 2 , X n 1 , X n 2 , X n , Y n 1 , Y n 2 ) be i.i.d. according to (40). Assign three random bin indices (S j , W j , C j ) to each realization u n j for all j = 1, 2. Assume S j ∈ [1 : 2 n Ď Rs,j ], W j ∈ [1 : 2 n Ď Rw,j ], and C j ∈ [1 : 2 n Ď Rc,j ] such that s R s,j , s R w,j , s R c,j ≥ 0 for j = 1, 2. Apply the union bound to the reliability constraint in (24), which vanishes when n → ∞ by using an SW decoder to estimate U n j from (C j , W j , Y n j ) if [32, Lemma 1] The key uniformity (25) constraint is satisfied if [32, Theorem 1] since (65) ensures that the three random indices (S j , W j , C j ) are almost mutually independent and uniformly distributed. Suppose a virtual joint encoder assigns six indices (S 1 , W 1 , C 1 , S 2 , W 2 , C 2 ) to each realization pair (u n 1 , u n 2 ). This virtual encoder is an operational dual of the virtual decoder used in the proof of [18,Theorem 1]. Using the virtual joint encoder, the strong secrecy constraint in (27) and the strong helper data independence constraint in (29) are satisfied if [32, Theorem 1] and because (66) ensures that (S 1 , W 1 , C 1 , W 2 , C 2 ) are almost mutually independent; whereas, (67) ensures that (S 2 , W 2 , C 2 , W 1 , C 1 ) are almost mutually independent.
Remark 2. The set of equations considered in (65)-(67) cannot be imposed for the joint secrecy-leakage constraint in (9) for general probability distributions P X1 X2XY1Y2 , since to impose (9) one would replace (66) and (67) with which would also imply the mutual independence of secret keys in (7). However, the inequalities in (65) and (68) cannot be satisfied simultaneously in general as H(U 1 ) + H(U 2 ) ≥ H(U 1 , U 2 ). This problem is avoided in the proof of Theorem 1 by imposing the inequality in (44) rather than (65).
The public randomness C j is almost independent of X n j , so it is almost independent of ( X n 1 , X n 2 , X n , Y n 1 , Y n 2 ), if we have [32, Theorem 1] Thus, the public indices (C 1 , C 2 ) can be fixed by generating them uniformly at random. U n j can be generated according to P U n j | X n j Cj for j = 1, 2 obtained from the binning scheme above to compute the bins (S j , W j ) from U n j for j = 1, 2. This procedure induces a joint probability distribution that is almost equal to P U1U2 X1 X2XY1Y2 that is fixed in (40) [34,Section 1.6].
Applying the Fourier Motzkin elimination to (64)-(67) and (69), we obtain the inequalities Observe that we have due to (23) and (63). We therefore obtain where (a) follows by (84) for some > 0 such that → 0 when n → ∞ due to (84)-(87). Furthermore, using the selection lemma, there exists a binning that achieves all rate tuples given in (90)-(92). Consider the privacy leakage for a binning that satisfies (90)-(92). Since C 1 and C 2 are public, we can bound the privacy leakage as follows.

B. Achievability Proof for the CS Model
The achievability proof for the CS model follows by applying the one-time padding step used in Section V-B.

C. Outer Bound Proofs for the GS and CS Models
Suppose for some δ n > 0 and n, there is a pair of encoders and decoders such that (24)- (29) are satisfied by some keyleakage-storage tuple ( s R s,1 , s R s,2 , s R , s R w,1 , s R w,2 ). Using (24) and Fano's inequality, we obtain where (a) permits randomized decoding, n = δ n max{ s R s,1 , s R s,2 }+H b (δ n )/n such that n → 0 if δ n → 0. Let U j,i (S j , W j , X i−1 ), which satisfies the Markov chain U j,i − X j,i − X i − Y j,i for all i ∈ [1 : n] and j = 1, 2.
Remark 3. For the choice of U j,i = (S j , W j , X i−1 ) (and similarly for U j,i = (S j , W j , Y i−1 j )) for j = 1, 2, U 1,i − X 1,i − U 2,i do not form a Markov chain for all i ∈ [1 : n] although for the inner bound we use this Markov chain. This is the reason why inner and outer bounds do not match in general.
Proof for (31): We obtain for the GS and CS models for j = 1, 2 that (a) follows by (25) and (96), (b) follows by (27), (c) follows by applying the data-processing inequality to the Markov chain where (a) follows by (26) and from the Markov chain W 1 − X n − W 2 , (b) follows because I(W 1 ; Y n 1 ) = I(W 1 ; X n 2 ) due to (23), (c) follows from the Markov chain W 1 − X n 2 −W 2 , (d) follows by (96), (e) follows because the channel and source are memoryless and from the Markov chain in (98), and (f ) follows from the definition of U j,i .
Proof for (33): We have for the GS and CS models that where (a) follows by (26) and from the Markov chain W 1 − X n −W 2 , (b) follows by (96) and from the Markov chain S j − (W j , X n ) − Y n for j = 1, 2, (c) follows because the channel and source are memoryless, (d) follows from the Markov chain 101) and (e) follows from the definition of U j,i . Proof for (34): Observe for the GS model for j = 1, 2 that where (a) follows by (28), (b) follows from the encoding steps, (c) follows by (96), (d) follows because the source and channel are memoryless, (e) follows from the data-processing inequality applied to the Markov chains in (98) and (101), and (f ) follows from the definition of U j,i . Proof for (37): Observe for the CS model for j = 1, 2 that n( s R w,j +δ n ) (a) ≥ I(S j , W j ; X n j )−H(S j |W j )+H(S j , W j | X n j )

(b)
≥ I(S j , W j ; X n j ) + I(S j ; W j ) where (a) follows by (28), (b) follows because X n is independent of S j and from the encoding step, (c) follows because the source and channel are memoryless, (d) follows by applying the data-processing inequality to the Markov chain in (101), and (e) follows from the definition of U j,i . Proof for (35): We have for the GS model for j = 1, 2 that n( s R s,j + s R w,j ) where (a) follows by (25), (b) follows by (27), and (c) follows from the definition of U j,i .
Proof for (38): Similarly, we have for the CS model for j = 1, 2 that n s R w,j ≤ n i=1 H(S j , W j , X i−1 ) where (a) follows from the definition of U j,i . Proof for (36): We obtain for the GS model for j = 1, 2 and j as defined in (30)  = H(S j , W j , W j ) + I(S j ; W j , W j ) + I(W j ; W j ) + nδ n H(S j , W j , W j , S j , X i−1 ) + 2δ n n + δ n (106) where (a) follows by (25), (b) follows by (27) and (29), and (c) follows from the definitions of U j,i and U j ,i . Proof for (39): We have for the CS model for j = 1, 2 and j as defined in (30) that H(W j , W j , S j , S j , X i−1 ) + δ n n + s R s,j (108) H(U j,i , U j ,i ) + δ n n + s R s,j where (a) follows by (29) and (b) follows from the definitions of U j,i and U j ,i .
Remark 4. (106) and (108) are the only places we use the strong helper data independence constraint in (29) and it does not seem straightforward to obtain the inequalities in (106) and (108) without (29).

VII. CONCLUSION
We derived inner bounds for the multi-entity key-leakagestorage regions for GS and CS models with strong secrecy, a hidden identifier source, and correlated noise components at the encoder and decoder measurements that are modeled as BCs. The inner bounds are valid for any finite number of entities that use the same hidden source to agree on a secret key. We argued that the mutual key independence constraint we impose makes the proposed multi-entity key agreement problem a proper multi-user extension of the classic single-enrollment key agreement problem, unlike the multienrollment key agreement problem considered in the literature. A set of degraded and less-noisy BCs was shown to provide strong privacy without a need for a common randomness. We established inner and outer bounds for the key-lekagestorage regions for a two-enrollment model with measurement channels that are valid for SRAM and RO PUFs. Inner and outer bounds were shown to differ only in the Markov chains imposed and they match if the storage and privacy-leakage rate constraints are removed. Two examples illustrated that depending on the constraints of the practical scenario, single or multiple enrollment might perform better in terms of the secret-key vs. privacy-leakage rate ratio. In future work, we will find a set of probability distributions for which the strong helper data independence constraint in the two-enrollment model can be eliminated.
ACKNOWLEDGMENT O. Günlü thanks Rafael F. Schaefer for fruitful discussions.


extfieldsofstudy:
Computer Science
Mathematics

text:
QROSS: QUBO Relaxation Parameter Optimisation via Learning Solver Surrogates

An increasingly popular method for solving a constrained combinatorial optimisation problem is to first convert it into a quadratic unconstrained binary optimisation (QUBO) problem, and solve it using a standard QUBO solver. However, this relaxation introduces hyper-parameters that balance the objective and penalty terms for the constraints, and their chosen values significantly impact performance. Hence, tuning these parameters is an important problem. Existing generic hyper-parameter tuning methods require multiple expensive calls to a QUBO solver, making them impractical for performance critical applications when repeated solutions of similar combinatorial optimisation problems are required. In this paper, we propose the QROSS method, in which we build surrogate models of QUBO solvers via learning from solver data on a collection of instances of a problem. In this way, we are able capture the common structure of the instances and their interactions with the solver, and produce good choices of penalty parameters with fewer number of calls to the QUBO solver. We take the Traveling Salesman Problem (TSP) as a case study, where we demonstrate that our method can find better solutions with fewer calls to QUBO solver compared with conventional hyper-parameter tuning techniques. Moreover, with simple adaptation methods, QROSS is shown to generalise well to out-of-distribution datasets and different types of QUBO solvers.

Introduction
The Quadratic Unconstrained Binary Optimisation problem (QUBO) has become a unifying formulation for a wide range of combinatorial optimisation problems [1,2]. Its objective form is closely related to the Hamiltonian of the Ising model in statistical physics [3], where a variety of sampling methods based on Markov Chain Monte-Carlo have been developed to probe the energy landscape and find its ground states [4]. These form the basis of commercially available QUBO solvers, such as those implemented on quantum annealers [5] and quantum-inspired computers [6].
Many combinatorial optimization problems can be expressed in the form of min x∈{0,1} n x T Qx subject to Cx = d where decision variables are binary and the constraints are linear. The conversion of such problem into a QUBO form is straightforward where the problem can be rewritten as min x∈{0,1} n x T Qx + A Cx − d 2 [3], where A is referred as the relaxation parameter. The challenge is that the quality of the solutions found by a QUBO solver is sensitive to the choice of the parameter value. Inappropriate choice of relaxation parameter value could lead to a solution that is either infeasible or far from optimal. Hence, relaxation parameter optimisation is an important step in this problem-solving process.
In the context of emerging computing technology and new QUBO solvers, existing methods do not capture the common structure shared by instances of a problem. For industry applications like vehicle route planning [7] and resource allocation [8], instances of the same problem are solved repeatedly. We argue that there is valuable information in the history of solved instances that can be extracted as prior knowledge for a relaxation parameter optimisation method.
In this paper, we propose a QUBO Relaxation parameter Optimisation method based on QUBO Solver Surrogates (QROSS) to mitigate these issues. As indicated by its name, QROSS relies on a surrogate model that approximates important characteristics of a QUBO solver. We build the surrogate model using machine learning on data obtained from the QUBO solver when applied to solve instances of a problem. QROSS then proposes promising relaxation parameters using the surrogate model on new instances, effectively reducing the number of calls to QUBO solver.
The solver surrogate in QROSS only models certain aspects of the original QUBO solver that are necessary for relaxation parameter tuning. Given a problem instance and a relaxation parameter, the solver surrogate predicts the probability that the QUBO solver finds a feasible solution, as well as the objective energies. Note that no explicit solutions are predicted by the solver surrogate. An evaluation on the solver surrogate, which we will model as a carefully designed neural network, is much cheaper/faster than a call to a QUBO solver.
Based on the solver surrogate, we propose two search strategies that will provide the trial values of the relaxation parameters. Hence, given a new instance of a problem, QROSS can propose promising relaxation parameters that can be used to form the QUBO problem that is then solved by a QUBO solver.
QROSS has the following features that make it unique among other hyper-parameter tuning methods: • QROSS captures the common structure of a class of problems. This is achieved by learning a solver surrogate from a history of problems. Knowledge in previous problems help to solve new problems of the same type.
With QROSS, users are expected to find better relaxation parameters and solutions with fewer calls to the QUBO solver; • Given a new problem of the same class, QROSS is able to predict the landscape of the objective function and help users understand the expectations without resorting to the expensive QUBO solving step; • QROSS allows trade-off between optimality and number of calls to a QUBO solver. If an application only allows one call per problem, QROSS produces one parameter candidate, which has a good chance to satisfy the constraints of the problem. If a user can afford more calls for a problem, QROSS can change the strategy and propose more candidates to find better parameters and solutions.
Our contributions are as follows: 1. We introduce QROSS, a machine learning method to extract knowledge from QUBO problem instances solved in the past to facilitate the relaxation parameter optimisation process for a new instance of the same problem.
2. We propose two offline and one online parameter selection strategies. Offline strategies learn from instances in the past and propose promising parameters without calling QUBO solvers. The online strategy learn from the instance-to-solve and improves parameter by exploiting results from QUBO solvers.
3. We take Traveling Salesman Problem (TSP) as a case study to demonstrate that QROSS learns knowledge from history effectively and proposes promising parameters.
In experiments, we use Fujitsu Digital Annealer (DA) as a QUBO solver and compare QROSS with representative optimisation methods. Results show that QROSS outperforms these methods by a large margin. We train QROSS on synthetic dataset and find it is also adequate on out-of-distribution and real-world problem settings. We repeat all experiments with another QUBO solver, the hybrid solver Qbsolv from DWave Quantum Annealer. Results suggest that QROSS generalises well for the different QUBO solvers we used.

Related Work
Many parameter-sampling methods have been proposed in the literature, for example, SMAC [9] uses random forest; GPyOpt [10] and Spearmint [11] use Gaussian Process, and Hyperopt [12] uses tree-structured Parzen estimator. Optuna [13] allows the user to dynamically construct the search space.
While the above papers are generic in that they do not take problem specific features into consideration, SATzilla [14] is an automated approach for constructing per-instance algorithm portfolios for the Satisfiability Problem (SAT) that uses the so-called empirical hardness models to choose among their constituent solvers. Instance-Specific Algorithmic Configuration (ISAC) [15] predicts the impact of different parameters and the performance by learning from instances. In [16], ISAC was proposed to choose appropriate MAXSAT problem. In [17], a per-instance configurable portfolio which is able to adapt itself to every planning task was proposed.
Parameter tuning aside, generic constraint programming solvers have been proposed, for example [18] and [19], that perform local search by taking the constraints explicitly into the consideration.
All methods mentioned above are typically not concerned with feasibility of solutions, since feasibility is guaranteed through the algorithm itself. However feasibility is not guaranteed by a QUBO solver. Tuning QUBO relaxation parameters will need to take feasibility into consideration.

The QROSS Method
QROSS exploits knowledge learned from past instances of a problem when solving a new instance. In this section, explain what is the knowledge, why is it useful, how to learn it, and finally how to exploit it.

Motivation
As noted in the Introduction, the choice of the relaxation parameter value directly impacts feasibility. In other words, the concern is whether a solution for the QUBO model instantiated with a particular parameter value can be found that satisfies the constraints of the original constrained problem. We first investigate the relation between feasibility and parameter value. We denote the relaxation parameter as A, a problem instance as g, and define the probability of feasibility as the proportion of feasible solutions produced by the QUBO solver. Eq.1 gives an estimate of it when evaluated on a large number of solutions on a common parameter setting: The value of P f lies in the range [0, 1]. The relation between A and P f is shown in the left column of Fig.1. The upper part of the figure shows the results from a QUBO solver, which is the Digital Annealer from Fujitsu [6]. For each point, there are 128 solutions returned by the QUBO solver, given a problem instance g and a relaxation parameter A. Because of the stochastic nature of QUBO solvers, the result could be different even if the input is the same. The dashed line that we fitted to the points looks like a sigmoid shape. The QUBO solver for the lower part of the figure is Simulated Annealing on CPU. The same shape can be observed from the lower part. The relaxation parameter A is also closely related to the objective energy of a solution. The upper right part of Fig.1 is produced by Fujitsu DA. Each point represents the minimum of fitness of 128 solutions from the QUBO solver, given an instance g and a parameter A. The orange dashed line sketches the envelope of the points and presents as a dipper shape. The red star-shaped marker is the bottom of the dip, where the optimal solution appears. The parameter that leads to near optimal solution is called optimal parameter. The lower right part of the figure is for Simulated Annealing, the dipper shape is not as obvious as the upper part. This is because SA gets stuck in local minima more often for the current parameter setting and has an almost flat landscape.
We explain the intuition behind the observation as follows. As one increases the relaxation parameter value, the constraints, in the form of penalty, gradually dominates the QUBO objective. Thus, A QUBO solver will be more likely to find feasible solutions. Meanwhile, the part of the QUBO objective corresponding to the original objective becomes less prominent and sometimes even disappears. The QUBO solver will not be able to respond to the tiny difference presented by the objective part and fails to find (near-)optimal solutions. 1 On the other hand, if one decreases the relaxation parameter over a certain threshold, i.e., to the left part of the dip, there will be fewer feasible solutions, and therefore fewer chances for the QUBO solver to find a solution with a smaller fitness value. An analytic approximation of "expectation of minimum fitness" will be given in section Minimum Fitness Strategy, show that a promising relaxation parameter balances the weight of objective and feasibility.
Through these observations we arrive at the following Hypothesis: Optimal solutions appear within 0 < P f < 1, i.e., on the slope of the Sigmoid shape. With more experiments we confirm this hypothesis holds true for all instances in TSPLIB [20] with Fujitsu DA and QAPLIB [21] with SA on CPU. 2 Hence, if one can effectively predict P f and objective energy, the parameter optimisation process would be much easier. This is especially useful when the objective landscape is flat, because conventional parameter optimisation methods cannot effectively find clues from a flat landscape.
Fortunately, instances of a problems share common structure. For example, a car company has to do vehicle routing in a city many times a day. A logistic company has to manage [7] allocations in a warehouse repeatedly [8]. For instances of a problem, the relation between relaxation parameter A, P f and objective energy are learnable. With these knowledge, we can select promising parameter values before calling the QUBO solver. The selection of parameter value does not involve exact-solution-finding process, thus it is much cheaper and faster than using the QUBO solver.

Network Architecture
Next we describe the network architecture for the surrogate. We have two inputs for the solver surrogate: the problem instance, denoted as g, and relaxation parameter, denoted as A. A could be a scalar of vector, depending on the number of relaxation parameter of a problem. Learning from g is tricky because the size of an instance could vary largely. We use an feature extraction layer that handle problems of different sizes. After the feature extraction, instances of different sizes are converted into fixed-size feature vectors. The feature vectors, together with the relaxation parameter A, are passed to a fully connected (FC) layer for further calculation.
We use the surrogate to predict two types of outputs. The first one is the probability of feasibility P f , which is a function of g and A. P f (g, A) is in the range [0, 1]. We can use sigmoid function as the activation of the last layer of the FC network and use Binary Cross Entropy loss (BCE) to learn P f (g, A). The other target is to predict fitness, i.e., the objective energy of a solution that satisfies the constraints of the original problem. Predicting the fitness is tricky. When there is no feasible solutions, there will be no fitness available for the instance, and this is problematic for neural network training. Our workaround is to predict the statistics of QUBO objective energy, E avg (g, A) and E std (g, A). E avg (g, A) and E std (g, A) are functions of g and A. Any reasonable loss criteria for regression should work for the learning these functions. We use Huber loss as we are expecting many outliers in the dataset, due to the stochastic nature of a QUBO solver. Once we have P f (g, A), E avg (g, A) and E std (g, A), we can calculate the fitness using the analytical expression given in section Minimum Fitness Strategy. A detailed description of the architecture can be found in Appendix.
The training workflow of the QROSS is described in the upper part of Fig.2. Given an instance g and a parameter A, we find the prediction and the ground truth of P f , E avg and E std by QUBO solver and solver surrogate. The difference between the prediction and the ground truth, i.e., the loss, can be used to train solver surrogate.

Data Preparation
To train the solver surrogate, we prepare a dataset from past instances and solutions of a solver. Because heuristic QUBO solvers have a stochastic nature, it usually returns a batch of solutions and corresponding objective energy for an instance g and parameter A. One can check the feasibility of each solution and calculate P f for the batch. We observed that the objective energy in a batch usually follows a bell-shaped distribution. Therefore, we use Gaussian sufficient statistics E avg (g, A), E std (g, A), i.e., mean and standard deviation, to characterise the objective energy.
A good coverage of sampling in terms of A helps to improve training. Since the slope of the sigmoid shape is the region that we are interested in, we make sure that {A | 0 < P f (g, A) < 1} are well sampled for each problem. Including the plateau region of the sigmoid shape helps to prevent model overfitting and improves accruacy. Hence, we make sure that at least a sizable number of samples in {A | P f (g, A) = 0 or 1} for each problem are included in the dataset.
Data augmentation also helps the training. For example, pre-processing techniques, e.g. shifting or scaling, moves A of different problems to the same order of magnitude so that learning and prediction become easier. Normalisation helps the convergence of the training curve. These techniques are applied in our experiments.

Inference
Next, we describe how to exploit the knowledge in the surrogate. The diagram of the inference is shown in the lower part of Fig.2. Given an instance g and a parameter A, the solver surrogate predicts P f , E avg and E std . Parameter selection strategies try different A, until it finds a promisingÃ. We propose two offline strategies, which do not need to call a QUBO solver. Since now we focus on optimising A for one instance g, we drop g and use notation P f (A), E avg (A) and E std (A) for the sake of simplicity.

Minimum Fitness Strategy
First, we propose Minimum Fitness Strategy (MFS), which attempts to predict the optimal relaxation parameterÃ that leads to minimum fitness. Given P f (A), E avg (A) and E std (A), with a Gaussian assumption on the conditional distribution of E avg (A) and E std (A) given A, we are able to calculate the expectation of minimum fitness analytically: where E d represents the expectation of minimum fitness. Φ is the Cumulative Distribution Function of Gaussian function. B is the number of soltuions in a batch. (see appendix for the details of derivation and calculation). Therefore the expectation of minimum fitness D min can be seen as a function of P f (A), E avg (A) and E std (A). Then, the optimal relaxation parameterÃ can be found by conventional optimisation method.
We use shgo optimiser from scipy to search parameter search. The search is based on the solver surrogate and does not involve QUBO solver. The search can be done on the CPU within a few seconds.

P f -based Strategy
Since P f provides clue for the location of optimal solution, we can use P f -based Strategy (PBS) to find promising parameters. PBS allows customisable parameter search. The equation for searching the parameters is shown as eq.(3) Here, p is a user-defined parameter representing the desired feasibility probability and Eq.(3) finds parameter A that corresponds to the desired value. As before, it does not involve QUBO solvers. This strategy is useful for a variety of application scenarios. For example, if obtaining a feasible solution in one trial is of primary importance and its objective value is of secondary importance, then p = 90% Would be a reasonable choice. If the user can afford a few trials for each instance, say 5 trials, one could set p = 90%, 70%, 50%, 30%, 10%, and find the corresponding A for each instance.

Travelling Salesman Problem as A Case Study
The travelling salesman problem (TSP) is a classical combinatorial optimisation problem, where we are given a list of vertices and their pairwise distances, and we want to visit every vertex exactly once and return to the starting vertex.
The goal is to minimize the distance of the tour, i.e., we want to find the shortest Hamiltonian cycle.

QUBO Form and Relaxation Parameter
In [22], a QUBO formulation that only involves a quadratic number of terms in the number of cities is proposed. Without loss of generality, we can focus on the case where the graph is fully connected, as we can always introduce edges of infinite distances otherwise. We let d uv be the distance between city u and city v. We require n 2 variables for an n-city instance. The first subscript of x represents the city and the second indicates the order that the city is going to be visited at. Let x v,j be the indicator variable that the city v is the j-th city to be visited. Notice that the constraint implies that this satisfies the permutation condition. The formulation is as follows: where describes the total distance travelled and describes the constraints to be a feasible cycle.
To find the optimal relaxation parameter A for a TSP instance, conventional methods require a few calls to QUBO solvers. Meanwhile QROSS is designed to work without calling a QUBO solver. In order to have comparable experiments, we design an online strategy and evaluate it along with the two offline strategies we proposed in section Inference.

Online Fitting Strategy for Case Study
Online Fitting Strategy (OFS) exploits results of an instance from QUBO solvers to improves parameter search for the same instance. As the relation of A and P f resembles the sigmoid shape shown in Fig.1, We can approach promising parameter through the process of curve fitting. The ansatz function to fit is: where θ s and θ o represent the scaling and the offset of the sigmoid shapes in the direction of A, respectively. As more hyper-parameters are evaluated, the parameters θ s and θ o will be estimated more accurately and facilitates the hyper-parameter optimisation better. We use P f (A) and F(A) to represent the ground truth of probability of feasibility and fitness with respect to A found by a QUBO solver. The pseudo algorithm for the Online Fitting Strategy is listed in Algorithm 1. Evaluate P f (A next ) and F(A next ); 7 Update A left or A right with A next if applicable; 8 end 9 Return the best A among history of F; Line 1 and 2 find the left bound and right bound of {A | 0 < P f (A) < 1}. The bound does not have to be accurate, as a QUBO solver has stochastic nature in its results. The approximated left bound can be found by evaluating a sequence of {A, A 2 , A 4 , . . .} on the QUBO solver, until P f (A left ) = 0 is found. A similar way applies for the right bound. The offline strategies in section Inference is able to provide good guess of initial A without calling QUBO solvers.
The search for A left and A right forms a history of P f , i.e., a collection of points {(A 0 , P f (A 0 )), (A 1 , P f (A 1 )), . . .}. The history of P f can be used to fit S(A, θ s , θ o ) in Line 4. In Line 5 and 6, A randomly-picked candidate A next is evaluated on the QUBO solver to find P f (A next ) and F(A next ). Then P f (A next ) will join the history of P f for the Sigmoid fitting (Line 4) in the next iteration. The loop terminates when maximum iteration is reached.

Experiments
In the experiments we compare QROSS with baseline methods on synthetic and real-world datasets. We evaluate the generalisation of QROSS on Fujitsu DA and Qbsolv from DWave. 3 All results are summarised in Table 1.
Dataset We construct a synthetic TSP dataset for the experiments. The whole dataset consists of 300 TSP instances with number of cities ranging from 20-30. The detailed generation method of the dataset is given in the Appendix. We use 270 instances for training the solver surrogate. The remaining 30 instances are used for testing.
Strategy One can use any of the strategies we proposed depending on application requirements. For the purpose of benchmarking, we use a mixture of the three strategies in the experiments. The composed strategy is shown as the following list.
1. Use MFS to proposes the first candidates 2. Then use PBS to proposes the next two candidates {A | P f (A) = 80% and 20%} 3. Use OFS to proposes more candidates The candidates proposed by the composed strategy are evaluated on a QUBO solver in sequence. The {A | P f (A) = 80% and 20%} in the second step is just an example of the usage. One can use other probability setting that suits the application. The trials in the first two step can be used for curve fitting in the third step.

Comparison with Generic Methods
We compare QROSS with three popular methods on the synthetic dataset for testing. The baseline methods are Bayesian Optimisation (BO), Tree-structured Parzen Estimator (TPE), and Random Search. We choose BO and TPE because they are famous for finding optimal solutions with fewer trials. Random Search is also included as a representative of exhaustive methods.
The optimal relaxation parameters of the instances in the synthetic dataset are all within the range [1,100]. Therefore we restrict the exploration of relaxation parameter space to [1,100] for all baseline methods. BO requires some random samples before the actual exploration and exploitation. For each instance, we draw 5 random samples from a uniform distribution U (1, 100), as we are only comparing the quality of the first 20 trials of baseline methods. : Comparison on test instances of the synthetic dataset. X-axis is the number of trials a method has taken. Y-axis is normalised gap between the near-optimal fitness and the best fitness found so far by a method. The normalised gap is averaged across all test instances. Fig.3 shows the comparison between QROSS and baseline methods. The blue solid curve for QROSS outperforms other baseline methods by at least 5% at the first trial and 2.9% at the third trial. This means MFS and PBS find better relaxation parameter. We emphasise that the parameters proposed by QROSS so far require no calls to the QUBO solver, whereas TPE and BO improve the parameters based on the results from the QUBO solver. Then the blue curve continue to decline gradually and stay below all other curves. This suggests OFS are also working well on the synthetic dataset. The shade around the curve represents 95% confidence interval. QROSS has a narrower confidence interval, which suggests that it generalise better to the synthetic dataset than other methods do.

Performance on A Real-world Dataset
The real-world dataset contains eleven TSP instances from TSPLIB [20]. We exclude TSP instances with number of cities N >= 90 because decomposition technique is not the focus of this paper. TSP with N <= 14 is also excluded because they are not challenging enough for QUBO sovlers. We use techniques mentioned in section Data Preparation to pre-process these instances. Fig.4 shows the comparison. The plot settings are the same as Fig.3.
From Fig.4 we know that the blue curve for QROSS outperforms all other baselines by a large margin of at least 8% at the first trial and 2.9% at the third trial. We emphasise that the parameters QROSS proposes in the first three trials require no calls to the QUBO solver. for the rest of trials, QROSS with OFS retains its leading position all the way. The  95% confidence interval of QROSS suggests that its generalisation on the real-world dataset is as good as these baseline methods.
We emphasise that the surrogate is trained on a synthetic dataset, which has problem size of 20-30. The problem size in the real-world dataset is 14-90. QROSS outperforms all baseline methods on this out-of-distribution dataset. This suggests the knowledge in the previous instances generalise well to instances of different size. In other words, QROSS benefits from the learned knowledge about general TSP instances.

Generalisation
We repeat all previous experiments on Qbsolv. The last four rows of Table 1 are based on Qbsolv. We construct the training dataset using solutions generated by Qbsolv. The rest of the experiment settings are similar to that for DA. In this experiments QROSS also outperforms all baseline methods on in-distribution and out-of-distribution dataset. This suggest QROSS generalises well to different QUBO solvers.

Conclusion
Relaxation parameter optimisation is important for finding optimal solutions in a QUBO problem. Most existing parameter tuning methods do not exploit knowledge from previous problem instances. In this paper, we frame the parameter optimisation problem as surrogate model minimisation problem, in which the surrogate is learned from instances of a problem. We capture the common structure of the problem and produce optimal relaxation parameter with fewer calls to QUBO solver.

A Ablation Study
In our experiments, QROSS works well with DA and Qbsolv solver because it is trained on dataset generated by Fujitsu DA [8] and Qbsolv [23] from DWave respectively. Next we cross these experiments, i.e., train QROSS on DA-based dataset and test it on Qbsolv, to check if it is something else other than the knowledge in the dataset helps QROSS to outperform all baseline methods. The experiment is carried out on the 30 testing instances of the synthetic dataset. The comparison is shown in Fig.5. We know from Fig.5 that the optimality gap achieved by QROSS on Qbsolv is larger than TPE achieves on Qbsolv. This is because Fujitsu DA and Qbsolv work differently, the knowledge from DA would not be able to generalise to Qbsolv. The performance lag is what we expected for the ablation study.

B The choice of Penalty Weight
Due to the limitation of manufacturing and computation technique, the quantum computers and classical computers we have today are not oracle computers. Errors in these computers could spoils the solutions they produced. In the context of QUBO problem, the choice of penalty weight is critical to the quality of the solutions. We use Minimum Vertex Cover (MVC) problem to demonstrate the vulnerability of classical computers and quantum computers to large penalty weights.

B.1 Introduction to Minimum Vertex Cover (MVC)
The MVC problem is a classical NP-hard optimisation problem. Given an undirected graph with a vertex set V and an edge set E, a vertex cover is a set of vertices such that every edge of the graph has at least one endpoint in this set. A minimum vertex cover is the vertex cover of the smallest size. A standard optimisation model for MVC can be formulated as follows. Let v j = 1 if vertex j is in the cover (i.e., in the subset) and v j = 0 otherwise. Then the standard optimisation model for this problem is: . The constraints ensure that at least one of the endpoints of each edge will be in the cover and the objective function seeks to find the cover using the least number of vertices.
The constraints in the standard MVC model can be represented by a penalty of the form where σ represents a positive scalar penalty. A common extension of this problem allows a weight w i to be associated with each vertex i. The QUBO form for the weighted minimum vertex cover problem is The choice of σ is application specific. For example, if w i = 1, σ = 2, adding every additional vertex into a minimum vertex cover will increase the objective energy by three. Removing every vertex from a minimum vertex cover will increase the objective energy by one. Theoretically, any σ > max(w i ) would ensure that a solver can find feasible solutions to the weighted MVC problem.

B.2 Quantum and Simulated Annealer on MVC
We are using Quantum Annealer DW_2000Q [24] from DWave (qa) and Simulated Annealing on CPU (sa) in this experiment. The MVC problem instances are randomly generated graphs with 65 nodes and 50% probability of connections between any pair of nodes. 4 Fig. 6 shows the relation between penalty weight A and normalised objective energy for a MVC problem on DWave DW_2000Q and Simulated Annealing on a classical computer. 1.14 Energy normalised to optimal SA QA Figure 6: Penalty weight v.s. energy normalised. X axis is the penalty weight in log scale. Y axis is the objective energy normalised to the minimum energy state discovered in a run of the experiment. The result is an average over 4 runs for different random seeds. qa is DWave quantum annealer DW_2000Q. sa is simulated annealing on classical computer. The MVC problem instances are randomly generated graphs with 65 nodes and 50% probability of connections between any pair of nodes. The weights of nodes are random numbers following uniform distribution over [0, 1).
From Fig.6 we understand that the objective energy increases along with the increase of penalty weight. Although this is true for both QA and SA, the reason behind the phenomenon is different for them.
For classical computers, one of the reason is due to the limited range and precision of the numbers and error in floating point operations [25]. The widely used double precision floating point numbers has floating point error in the calculation could accumulates and spoils the results. The floating point error problem is especially prominent when the penalty term dominates the objective energy, as demonstrated in Fig.6. Fujitsu DA uses FP64. We have similar observation in the main body of the paper as well. The workaround is to tune the penalty weight to mitigate the issue, which one of main focus of this paper. The other is to use Unums representation [26] in the cost of additional hardware/software complexity.
The choice of penalty weight is critical to Quantum computers as well. The imperfections in the quantum hardware is one of the reasons to the degradation in the quality of solutions. More specifically, all physical implementations of quantum computers suffer from analog control errors, in which the coefficients of the Hamiltonian implemented differ from those intended. [27] If the penalty term dominates the Hamiltonian, the objective of the original problem could be overwhelmed by the analog control error. The analog control error threatens to spoil the results of computations due to the accumulation of small errors. This problem was recognized early on in the gate model [28] and later in the annealing based model [29].

C Feature Extraction
In order to learn the common structure of problem instances, we employ deep neural network as it is able to extract feature related to the task in an automatic way. The following literature focus on feature extraction for Traveling Salesman Problem, which is taken as a case study in this paper.
[30] uses deep Graph Convolutional Networks to build TSP graph representations and output tours in a nonautoregressive manner via beam search algorithm. The representations in [30] is the probability of each edge being in the route of the optimal solution. [31] approaches this edge-level features from a different route. [32] and [33] uses reinforcement learning technique to facilitate the heuristic search processing for TSP.
These methods aim to estimate solutions and/or objective energies, whereas QROSS is for relaxation parameter tuning. We apply the feature extraction technique in [30,31] to capture the common structure shared by TSP.

D TSP Problem Generation
To create augmented data set, We use uniform distribution and exponential distribution as our random number generators to create the coordinates of the cities. The parameter for the exponential distribution is generated from uniform distributions over a range. The uniform distribution is generated on a bounded domain. After we generated the coordinate data, we then compute the corresponding Euclidean distance. E TSP problem pre-processing and post-processing In [34], it is shown that by changing the distance matrix, , an optimal tour corresponding to the original distance matrix is also the optimal tour for the updated distance matrix. In [35], Minimizing the Variance of Distance Matrix (MVODM) is proposed where π i is chosen such that the resulting distance matrix has minimal variance. It is shown empirically that such procedure improve greedy search algorithm. We perform such pre-processing to our distance matrix.
After we solve the TSP instance using the solver, we use the original distance matrix to compute the original distance of the original Hamiltonian cycle.

F Calculate Expectation of Minimum Fitness
Suppose the number of solutions in a batch is B for a given problem instance g and a relaxation parameter A, the number of feasible solutions in the batch is: Let d i represent the ith feasible solution in a batch,d is the minimum fitness in batch, then We define the expectation of minimum distance D min (·) as a function of P f (A), E avg (A) and E std (A): where E stands for expectation operator. All the input of eq.(10) can be obtained from the output of the surrogate model. Next we are going to find the analytical expression to estimate the expectation of minimum fitness. We use Φ (z) to represent the probability of a fitness d i in a batch of replicas being smaller than a given value z: Then we have the analytical expression of P (d i > z) and P d < z : The expectation ofd can be expressed by: The fitness (d 0 , d 1 , . . . , d m ) are non-negative random variables, which means eq.(14) could be approximated as: The exact distribution of d i is unknown. We observed in our experiments that the distribution is a bell shape in most of cases. Therefore, we assume the d i follows Gaussian distribution, which can be expressed in the following equation: Then Φ could be expressed as the cumulative density function of the Gaussian distribution: One can plug eq.(17) into eq.(15) to calculate the expectation of minimum fitness. Please note that when P f is very close to zero, there is no feasible solutions in a batch of replicas. the calculation of P d < z would be meaningless.
We set lim P f →0 D min = +∞. Lastly we can solve the following equation using off-the-shelf optimiser to find the optimal hyper-parameterÃ:Ã

G Surrogate model
The structure of the neural network is shown in Fig.7. It takes two inputs: 1. a problem instance g; 2. a relaxation parameter A. It predict probability of feasibility P f , statistics of objective energy E avg and E avg , as functions of g and A.
g  [30,31], in which the authors take adjacency matrix of TSP as fully connected graph and apply graph convolutional neural network to extract features.
Our experiments are based on the pre-trained model for feature extraction in [30]. The edge-level features cannot be directly applied in our case. Still, we can aggregate the edge-level features into graph-level ones, which is what we needed for QROSS method.
Then the feature vectors of g, together with the A, goes into fully connected layer for prediction. Since the nature of P f is different from that of E avg and E std , we train these target separately. We use Binary Cross Entropy (BCE) loss for P f and use Huber loss for E avg and E std . the prediction is much cheaper than the solution-find process in terms of time and energy, QROSS is able to find promising parameters.