extfieldsofstudy:
Computer Science
Physics

text:
Scalable Experimental Bounds for Entangled Quantum State Fidelities

Estimating the state preparation fidelity of highly entangled states on noisy intermediate-scale quantum (NISQ) devices is an important task for benchmarking and application considerations. Unfortunately, exact fidelity measurements quickly become prohibitively expensive, as they scale exponentially as O(3^N) for N-qubit states, using full state tomography with measurements in all Pauli bases combinations. However, it is known [Somma et.al. 2006] that the complexity can be drastically reduced when looking at fidelity lower bounds for states that exhibit symmetries, such as Dicke States and GHZ States. For larger states, these bounds have so far not been tight enough to provide reasonable estimations on today's (2022) NISQ devices. In this work, for the first time and more than 15 years after the theoretical introduction, we report meaningful lower bounds for the state preparation fidelity of all Dicke States up to N=10 and all GHZ states up to N=20 on Quantinuum H1 ion-trap systems using efficient implementations of recently proposed scalable circuits for these states. For example, we give state preparation fidelity lower bounds of (i) 0.46 for the Dicke State |D10,5>and (ii) 0.73 for the GHZ State |G20>. These match or exceed exact fidelity records recently achieved on superconducting systems for the much smaller states |D6,3>and |G5>, respectively. Furthermore, we provide evidence that for large Dicke States |DN,N/2>, we can resort to a GHZ-based approximate state preparation to achieve better fidelity.

Abstract-Estimating the state preparation fidelity of highly entangled states on noisy intermediate-scale quantum (NISQ) devices is an important task for benchmarking and application considerations. Unfortunately, exact fidelity measurements quickly become prohibitively expensive, as they scale exponentially as O(3 N ) for N -qubit states, using full state tomography with measurements in all Pauli bases combinations.
However, it is known [Somma et.al. 2006] that the complexity can be drastically reduced when looking at fidelity lower bounds for states that exhibit symmetries, such as Dicke States and GHZ States. For larger states, these bounds have so far not been tight enough to provide reasonable estimations on today's (2022) NISQ devices. In this work, for the first time and more than 15 years after the theoretical introduction, we report meaningful lower bounds for the state preparation fidelity of all Dicke States up to N = 10 and all GHZ states up to N = 20 on Quantinuum H1 ion-trap systems using efficient implementations of recently proposed scalable circuits for these states.
For example, we give state preparation fidelity lower bounds of (i) 0.46 for the Dicke State |D 10 5 and (ii) 0.73 for the GHZ State |G20 . These match or exceed exact fidelity records recently achieved on superconducting systems for the much smaller states |D 6 3 and |G5 , respectively. Furthermore, we provide evidence that for large Dicke States |D N N/2 , we can resort to a GHZ-based approximate state preparation to achieve better fidelity.

I. INTRODUCTION
Any advantage that quantum computing may have over classic computing usually relies on the principle of superposition, informally defined as being in multiple computational basis states at the same time. Quantum states that are in a superposition such that individual parts cannot be described independently of the state of others are called entangled states. Superposition and entanglement are essential prerequisites for the successful execution of additional gate operations per the overall structure of the quantum algorithm executed to ultimately lead to a desired quantum end state. Thus, proving that a NISQ device actually produces entangled states is a key necessary condition for any quantum computing success. Such low-level testing of quantum mechanical properties is necessary on quantum devices because environmental noise quickly modifies or partially destroys a pure quantum state (as it would be theoretically prepared); this noise process, usually Research  called decoherence, turns the actual state in a computing device into a so-called mixed state, which is a probabilistic mixture of our target pure quantum state and other pure states. The fidelity of a (mixed) quantum state measures to what extent the intended target (pure) quantum state has been realized; fidelity is defined to be 1 if we measure a correct pure quantum state, and it is close to zero, i.e., 1/2 n for a maximally mixed n-qubit quantum state that consist purely of noise.
Let us use the state 1 √ 2 (|00 + |11 ), i.e., the Bell state as the prototypical entangled state to illustrate the concept of fidelity. From a classical vantage point, it may appear tempting to test whether a NISQ device has actually produced a Bell state as follows: we execute a large number of runs of the generating circuit each with measurement at the end, and count how often we measure 00, 01, 10, and 11; if we measure 00 and 11 about 50 percent of the time each, we declare that the devices actually produces the Bell state. However, this logic is flawed because any classical device that simply returns 00 and 11 with probability 0.5 would have also passed our test, but such a cheat device would be useless and unable to properly execute additional gate operations that propagate the quantum state and its entanglement properly to execute any full quantum algorithm. The appropriate way to test entanglement is to calculate the fidelity measure (formally defined later), which measures the outcome in all possible combination bases of X, Y and Z axes, thus requiring 3 N tests for N -qubit states. This is usually called full state tomography. A fidelity measure of 1 indicates a perfect quantum state (with entanglement), where as a value of 0 indicates a state without entanglement.
Showing that entanglement exists on NISQ devices across more than 10 qubits is computationally prohibitive due to the exponential number of tests, but it nevertheless remains vital to demonstrate that the device actually leverages quantum mechanical principles.
A solution to this dilemma is to approximate or compute bounds on the fidelity measure that require fewer runs to calculate. [20] discovered such bounds for certain symmetric entangled states based on angular momentum operators. In this paper, we adopt and modify these generic bounds for the Dicke states, GHZ states, and approximate Dicke states -three well-known classes of entangled states; Dicke states are equal amplitude superpositions of all computational states of the same Hamming weight (i.e., number of ones) and GHZ states are Bell states that are generalized to larger qubit counts. We then show through a series of NISQ experiments that the bound from [20] has become useful in practice to provide evidence for entanglement across up to 20 qubits without resorting to exponentially expensive full state tomography. This is a testament to both advances in state preparation algorithms as well as as NISQ hardware improvements: we show that the original bounds are not useful for NISQ devices from just two years ago for Dicke states on as few as 5 qubits as they give negative values; they become slightly positive with our modifications without providing strong evidence of entanglement; however, newer technology (by Quantinuum) produces results good enough for these bounds to become meaningful lower bounds for the fidelity.
For example, we give state preparation fidelity lower bounds of (i) 0.46 for the 10-qubit Dicke State |D 10 5 and (ii) 0.73 for the 20-qubit GHZ State |G 20 . These match or exceed exact fidelity records recently achieved on superconducting systems for the much smaller states |D 6 3 and |G 5 , respectively. The original as well as our modified fidelity bounds require only three different measurement runs, thus making fidelity based entanglement verification easily accessible on NISQ devices even at high qubit counts.

II. RELATED WORK
We use the quantum fidelity F as a similarity measure between a prepared mixed state ρ, expressed as a density matrix, and a target pure quantum state ρ ψ = |ψ ψ| [13]: where √ ρ denotes the unique positive semi-definite square root of ρ such that √ ρ √ ρ = √ ρ † √ ρ = ρ, and where the first equality would also hold for a mixed state ρ ψ . In this paper, we consider pure target states, notably Dicke states |D N K and GHZ states |G N . We remark that quantum fidelity is sometimes also defined as F = √ F [18], [20]. In a straightforward way, the fidelity of a prepared N -qubit state can be computed by sampling ρ in 3 N different Pauli basis (N -fold tensor products of Pauli operators σ x ,σ y and σ z ) to reconstruct its density matrix, from which we can compute the fidelity with the target state.
If one wants to upper bound the quantum fidelity, then one can consider the measurement success probability, i.e. the overall probability of sampling non-zero amplitude states when measuring in the computational Z-basis: Using a double application of the Cauchy-Schwarz inequality, one gets F(ρ ψ , ρ) ≤ MSP (ρ ψ , ρ) [2]. Several works have shown that it is possible to estimate a lower bound on the quantum fidelity using only few measurement settings avoiding the full state tomography. Earlier works [11], [20] focused on special states that have unique symmetry which require only a small random subset of Pauli operators to estimate fidelity. Somma et al. [20] proposed expressions for estimating quantum fidelity for highly symmetric classes of multi-qubit state preparation i.e. rotational invariant states, stabilizer states and generalized coherent states. To estimate lower bound of fidelity, they only used the measurements that have relation with the symmetry operators of the density matrix. Also, [11] derived fidelity estimation for GHZ and W states using 2N − 1 measurements. A generalized approach by Flammia et al. [9] showed estimation techniques for all possible state preparation by measuring lesser number of Pauli operators that are close to the desired state and of greater importance. Mentioning the scalabilty issue in generalized N qubit state fidelity estimation, Elben et al. [8] showed how randomized measurement can enable fidelity estimation by comparing two states. Recently, Zhang et al. [25] proposed machine learning based approach for direct fidelity estimation as a classification problem using a constant number of expectations.

III. FIDELITY LOWER BOUNDS FOR DICKE STATES
We provide experimental lower bound estimations on the state preparation fidelity of Dicke states |D N K (up to N = 10). We first discuss the circuit preparation strategies, then provide  lower bound expressions for estimating quantum fidelity and present experimental result for fidelity estimation in Quantinuum's H1-2 quantum processor.

A. State Preparation
A Dicke state D N K is the equal weight superposition of all N K N -qubit basis states x with K Ones and N − K Zeroes, i.e., with Hamming weight wt(x) = K: We optimize and generalize existing work on Dicke state preparation by defining Dicke state unitaries U N k,K which prepare Dicke states |D N upon any input of a unary encoded Hamming weight k ≤ ≤ K: Dicke state unitaries were first defined for k := 0 [3], and there exists a linear-depth circuit construction for k := 0, K := N with 5 N −1 2 + 2(N − 1) two-qubit CNOT gates between neighboring qubits in a 1D Linear Nearest Neighbor (LNN) connectivity only [2]. We start from this construction and further reduce the CNOT count by 3 N −K−1 2 and 5 k 2 using ideas for upper [3] and lower [17] bounds on the Hamming weight , respectively. This effectively generalizes and also improves all of the former constructions by constant factors while preserving LNN connectivity.
To prepare a fixed-K Dicke state |D N K , we further parallelize our construction using a divide-and-conquer strategy [2], [4]: The main idea is to split the Hamming weight K across qubit sets of size N 1 = N 2 and N 2 = N 2 using 2K − 1 CNOT gates on an all-to-all connectivity before applying the Dicke state unitaries U N 1 0,K and U N 2 0,K in parallel. Additional gains can be made by implementing these unitaries with layers of bit-flipping Pauli-X gates, U Ni 0,K = X ⊗Ni · U Ni Ni−K,Ni · X ⊗Ni , resulting in the circuit structure shown in Figure 2, which prepares |D 9 3 with the following steps: Note that to prepare the Dicke state |D 9 6 we can simply omit the last layer of Pauli-X gates. Thus in our experiments we only prepare Dicke states |D N K up to K ≤ N 2 , for which our two-qubit gates amount to 5K(N −K)−3(N +K)+5 CNOTs, giving linear depth O(N ) and total gate count O(KN ).

B. Lower bounds for fidelity estimation
Following [20], we observe that Dicke States are rotationally invariant and thus are the unique simultaneous eigenstates of the (squared) total angular momentum J 2 = J 2 acting on qubit i. This means we can write the Dicke state |D N K as |j, j z in terms of its quantum numbers j = N and j z = N − 2K with: In general, the quantum numbers j and j z satisfy 0 ≤ j ≤ N and |j z | ≤ j and the difference between any two eigenvalues is at least 2 (because we use Pauli operators instead of spin-1/2 operators).  and |D 9 1 -|D 9 2 may be due to these experiments taking place right after a recalibration downtime of the QPU.
Now we can write a target state ρ j,jz := |j, j z j, j z | in terms of the operators J 2 and J z : . (7) We could compute the fidelity F(ρ j,jz , ρ) by measuring every Pauli-String that appears when we expand the product, i.e., all correlations between the operators J 2 and J z in in Equation (7). This process still requires measurement of an exponentially large number of observables but less than 3 N measurements of performing full state tomography. Instead we turn towards the construction of operators that allow us to compute a lower bound on the fidelity. We will define operators S Jz and S J 2 inspired by the density operator of Equation (7) such that they satisfy [S Jz + S J 2 ] |j , j z = e j ,j z |j , j z with e j ,j z = 1 for the correct Dicke state and e j ,j z ≤ 0 for (j , j z ) = (j, j z ). We define the operator S Jz as the left-hand product in Equation (7), for which |j , j z has eigenvalue 1 if j z = j z and 0 otherwise. 1 Expanding the product results in Pauli-Strings consisting of identities and σ z operators only. These commute pairwise, hence we can measure S Jz in the computational basis and thus in a single measurement setting. In fact, we get the measurement success probability Tr [S Jz , ρ] = MSP (ρ j,jz , ρ).
The same is not true if we expanded the right-hand product of Equation (7): we would get a large number of noncommuting Pauli-Strings. Instead we use the following definition for S J 2 [20], which is based on the observation that Dicke states have maximum angular momentum, and all other eigenvalues for J 2 deviate by at least N (N +2)−(N −2)N = 4N : for which |j , j z has eigenvalue 1 if j = j = N and integer eigenvalues ≤ 0 otherwise. Similar to S Jz , we can measure the operator S J 2 by computing Tr(J 2 τ , ρ) through measurement in the τ -basis, i.e., measuring in the three X-basis, Y -basis and Z-basis settings is enough.
If the sum of the operators S Jz and S J 2 are applied to a target Dicke state ρ = |j , j z , we get eigenvalue 1 for the correct Dicke state and eigenvalues ≤ 0 for all other rotational states. n general, if |φ = j ,j z c j ,j z |j , j z is a pure state, then we get Here |c j,jz | 2 is the probability of projecting |φ onto the state |j = N, j z = N − 2K . Our experimentally prepared mixed state ρ is a convex combination of pure states, yielding our lower bound

C. Experimental Result
We construct Dicke state D N K circuits using the preparation scheme described in Section III-A. We use IBM Qiskit [1]  on Quantinuum H1-2 with increasing shot count. Experiments were run for 280 and 840 shots for |D 8 4 and |D 10 4 , respectively, to achieve comparable confidence intervals. Plots are scaled (1:3) accordingly.
to prepare the circuits and we did not apply the Qiskit tanspilation and optimization options. To obtain the bounds on fidelity, first we generate the untranspiled QASM circuits for Dicke State |D N K where 2 ≤ N ≤ 10 and 1 ≤ K ≤ N/2. Next, we compute the bounds on the quantum fidelity using only three measurements (X, Y and Z basis measurement) instead of measuring 3 N state tomography circuits to measure the closeness of prepared Dicke state ρ to the expected Dicke state |D N K . We use the formulas of Section III-B to estimate the lower bound S Jz + S J 2 ρ on quantum fidelity D N K | ρ |D N K . We also show the probability of finding expected basis states among a total number of circuit executions (informally called measured success probability). The definition of S Jz in 8 shows that it is the measured success probability of the prepared Dicke state.
1) Overall Result: Figure 3 shows the bounds for Dicke State |D N K where N = 2 to 10 and K = 1 to N/2. The Dicke states from left to right are sorted according to the number of required CNOT gates in the untranspiled circuit. In X axis, Dicke states |D N K are represented as D − N − K : (C) where N is number qubits, K is hamming weight and C is number of CNOT counts. The plot shows the lower bound S Jz + S J 2 ρ and measurement success probability along with two-sided 68% confidence interval below and above the mean of the distribution. We observe statistical noise in the lower bound estimation for the smaller Dicke states (with few CNOTs). We also find that the fidelity lower bound decreases with the increasing number of CNOT gates on both simulator and real hardware. We also observe that the lower bound estimation of Dicke states is better in hardware than the ones in the simulator. For our largest Dicke state |D 10 5 , the simulator lower bound estimation on fidelity is 0.22 while the hardware estimation is 0.46. Quantinuum specification mentions that although the simulators provide a high fidelity representation of the hardware device's output, some variances between the results may occur as the noise models in the simulator cannot fully capture the real hardware behavior. This explains the difference between simulator and hardware estimations of lower bounds on quantum fidelities. The sudden inconsistency of lower bound estimation from |D 8 2 could be because the experiments before and after |D 8 2 were executed in a different time span. Moreover, the experiments after |D 7 3 were run just after the H1-2 machine became online after calibration. The plot also shows that the measured success probability represents the upper bound of quantum fidelity and the gap between lower bound & upper bound starts to increase for higher Dicke states.
2) Number of Shots: We ran the Dicke state circuits on Quantinuum H1-2E simulator using max(150, 4 * num states) shots where num states = N C K for N = 2 to 10 and K = 1 to N/2. Figure 4 shows cumulative measurements of bounds on quantum fidelity D N K | ρ |D N K . The left plot shows bounds on quantum fidelity for Dicke state |D 8 4 & the right plot shows the bounds for Dicke state |D 10 4 on H1-2 quantum device. In each plot, the X-axis shows the number of shots and Y -axis shows the cumulative bound estimations. Each plot shows three cumulative lines: lower bound S Jz + S J 2 , upper bound measured success probability and S j 2 . The plots also show two-sided 68% confidence interval below mean and above mean for lower bound and measured success probability S Jz . The plots are scaled according to the maximum number of shots used in the execution of the circuit. In each plot, we find that initially there are a lot of fluctuation in the lower bound and upper bound estimation but the estimations start to get balanced with the increasing number of shots. Similarly, the confidence interval region follows the lower bound estimation and shows how the estimation starts to be stable at the end. The cumulative measurement success probability of the experiments represents the upper bound of quantum fidelity which stabilizes faster  than the lower bounds.

IV. FIDELITY LOWER BOUNDS FOR GHZ STATES
We provide experimental lower bound estimations on quantum fidelity for GHZ states |G N (upto N = 20). First, we discuss the circuit preparation strategies, then we provide lower bound expressions for estimating quantum fidelity in state preparation and finally present experimental result for fidelity estimation in Quantinuum's H1-1 quantum processor.

A. State Preparation
An N-qubit GHZ state is equal weight superposition of all zeros and all one's defined as The straightforward way to prepare |G N is first initializing all qubits to |0 N . The next step is adding a Hadamard gate H to the first qubit and then consecutively adding CNOT gate between i th and (i + 1) th qubit where i = 0 to N − 1. This circuit generates a linear time complexity of N . Recent work showed that some CNOT gates can be rearranged in GHZ state preparation without affecting the output state that will allow execution of the CNOT gates in parallel [7]. Such rearrangement reduces the execution step of the GHZ circuit and generates logarithmic depth GHZ circuits with time complexity 1 + log 2 N . Figure 5 shows N = 10 qubit linear (left) and logarithmic depth (right) GHZ state preparation from initial state |0 N . Both circuits contain N gates (1 Hadamard and (N − 1) CNOTs). The dotted lines in the right plot indicate slices within which the CNOT gates can be performed in parallel. The linear depth circuit (left) requires serial execution of the CNOT gates after the first Hadamard gate that induces linear time complexity of N = 10. In the logarithmic depth |G N circuit, the CNOT gates in each time slice executes in parallel that generates logarithmic complexity of 1 + log 2 N = 1 + log 2 (10) ≈ 4.

B. Lower bounds for fidelity estimation
GHZ states lies in the family of stabilizer states that are used in quantum error correction procedures. From [20], we find that stablilizer states can be defined by symmetry operators, Here the stablilizer operatorÔ s are products of Pauli operators that take +1 or -1 as possible eigen values. An stabilizer state |ψ can be uniquely defined by 13 where G s = {Ô 1 , ..,Ô s } forms the stabilizer group. G s can also be defined using its L linear independent generators [10], G s ≡ (ĝ 1 , ...,ĝ L ) wherê Thus, the state |ψ can be expressed as |ψ = |g 1 = 1, ..., g L = 1 while the density operator ρ ψ can be defined as ρ ψ = |g 1 = 1, ..., g L = 1 g 1 = 1, ..., g L = 1| The fidelity of a prepared state |ρ can be measured by computing the expectations listed in 15. From [20], we find that a lower bound on the fidelity can be estimated using the When the operator S Gs is applied to the state |ψ = |g 1 , ..., g L , we get, S Gs |g 1 , ..., g L = e g1,...,g L |g 1 , ..., g L , (g i = ±1) (16) where e 1,...,1 = 1 and e g1,...g L ≤0 . The lower bound on fidelity between prepared state ρ and expected state |G N becomes Using the formulas above for an N-qubit GHZ state |G N , the stabilizer group can be defined as G s ≡ (σ 1 lower bounds on GHZ state fidelities using only the Pauli X and Pauli Z measurement and S Gs will be defined as where S Gx represents product of Pauli X operator on i = 1 to N qubits and S Gz represents sum of the pairwise products of Pauli Z operator on i th and (i + 1) th qubit upto N th qubit. For a 4-qubit GHZ state, we get N = 4 and L = 4, Basically we need to compute the expectation of the product of Pauli X measurement on each of the four qubit by taking +1 for even parity and −1 for odd parity. Similarly, we compute pairwise consecutive product of Pauli Z measurement on qubit 1 & 2, 2 & 3 and 3 & 4 by checking the pairwise parity. Then we get an estimation of the fidelity lower bound for GHZ state by taking the average of those products and subtracting N − 2 from the sum.

C. Experimental Result
We prepare linear depth and logarithmic depth GHZ states |G N circuits mentioned in Section IV-A. Similar to the Dicke state circuits, we use IBM Qiskit [1] to construct the circuits with no compilation and optimizations. In this experiment, we generate untranspiled QASM circuits for GHZ states |G N , where 2 ≤ N ≤ 20. We compute the bounds on quantum fidelity using only two measurements (X and Z basis measurement) instead of measuring 3 N state tomography circuits to measure the closeness of prepared GHZ state ρ to the expected state |G N . To estimate lower bound S Gs ρ , we followed the idea described in Section IV-B. Similarly, we compute 68% confidence interval below and above mean for the lower bound and measurement success probability of our prepared GHZ state. Figure 6 shows the bounds for GHZ state |G N where N = 2 to 20. The GHZ states are sorted left to right according to CNOT gate counts in untranspiled GHZ circuits. The left plot shows bounds for logarithmic depth GHZ state circuits, while the right plot shows bounds for the linear depth GHZ circuits. Both plots show lower bound on quantum fidelity and measurement success probability along with 68% confidence interval below and above the mean of the distribution. For both log and linear depth circuits, the lower bound on fidelity decreases almost linearly, albeit few variations with increasing CNOT counts for higher N . The gap between lower bound and measured success probability starts to increase with increasing CNOT count. The linear circuit shows faster drop for higher GHZ states than log depth circuit which is expected as logarithmic depth circuits require less time steps for circuit execution. In both plot, we observe few irregular ups and downs in the lower bound estimation that could be because of the statistical noise. The gap between lower bound and confidence interval region stays almost close in all GHZ states. Our investigation find that expectation S Gs (defined as product of pauli X operations on each qubit) starts to drop from |G 11 and becomes more significant from |G 17 . There are two flat decrease from linear depth GHZ state |G 17 to |G 18 and |G 19 to |G 29 which could be caused by the qubit rearrangement in Quantinuum H1-1 device.

1) Overall Result:
2) Number of Shots:: To see the effect of the number of shots on quantum fidelity of GHZ states, we execute the GHZ states |G N where 2 ≤ N ≤ 20 for 200 shots in Quantinuum H1-1 quantum device. We use 200 number of shots for all GHZ states. We estimate the cumulative lower bound distribution on quantum fidelity G N | ρ |G N from initial to up to i th shot where i = 2 to 200. two plots (left and right) show linear depth |G 10 , |G 15 & |G 20 states. Each plot shows lower bound estimation, measured success probability and expectations of Pauli X & Z measurements for GHZ state |G N preparation. The plots also include 68% confidence interval region below and above the mean of lower bound and and measured success probability. In each plot, the cumulative measurements starts with a lot of fluctuations that becomes balanced with increasing number of shots. We observe that all the measurements are more steady in |G 10 plot for both linear and logarithmic depth circuits than |G 15 and |G 20 plots. We find that the expectation S G § (purple line) drops significantly for |G 20 . It could be possible that larger GHZ states require more shots to become stable. Additionally we observe that the confidence interval region for log depth |G 20 circuit is less stable than linear depth circuit. Besides, the upper bound measured success probability stays almost stable in both linear and log depth circuits.

V. FIDELITY LOWER BOUNDS FOR APPROXIMATE DICKE STATES
Section III-C shows the linear relationship between Dicke state fidelites with increasing number of qubits (N ). For higher N , direct fidelity computation becomes exponentially expensive in terms of time and resource. Also the lower bound estimation becomes infeasible because of the statistical noise in the system. As alternative, preparing an "approximate" Dicke state might be another option that yields Dicke states with weaker theoretical fidelities but less noise.

A. State Preparation
We prepare two versions of "Approximate Dicke states" with (polynomially) vanishing fidelity to an actual Dicke state |D N K : Product states |P S N K for small K and Even/Odd Hamming weight states |E N , |O n for large K with even/odd parity, respectively. These are defined as:   Fig. 8. Preparations of Approximate Dicke States, N = 9: (left) Product state |P S 9 4 = ( √ 5 9 |0 + √ 4 9 |1 ) ⊗9 , using parallel Ry(2 cos −1 ( 5/9)) rotations. (center) Schematic circuit for Even and Odd Hamming weight state preparation for K even or odd, respectively. The circuit adds the Hamming weight parity of a + 8 state into the first qubit q 0 . Moving a column of H-gates to the right (blue arrow) swaps controls and targets of the CNOT gates. This leads to (right) a logarithmic-depth |G 9 -based state preparation circuit for Even/Odd Hamming weight states |E 9 , |O 9 , depending on the parity of K.
Theoretical product state from [5] is a probabilistic approach of preparing N − qubit symmetric product state ( 1 − K/N |0 + K/N |1 ) N . The state is prepared by adding rotation R y (2wcos −1 ( (n − k)/n) on each qubit. We then measure them in X, Y , and Z basis. This yields approximate Dicke state with success probability . To find approximate Dicke state fidelities, we generate Product state |P S N K for N = 2 to 10 qubits and K = 1 to N/2.
Another way of generating approximate Dicke states could be by preparing N -qubit odd or even hamming weight states based on the GHZ circuits discussed in Section IV-A. Even hamming weight N-qubit states are close to Dicke states |D N K with K even and N-qubit odd hamming weight states give approximate Dicke state |D N K with K odd. Odd/Even hamming weight states provide Dicke state fidelities of N K /2 N −1 . N -qubit even states are constructed by adding H N at the end of N -qubit logarithmic depth GHZ circuits. For Nqubit odd states, we just flip one qubit of the even hamming state. For example, 3-qubit even hamming weight state 1/4(|000 +|011 +|101 +|110 ) are generated from 3-qubit GHZ states 1/2(|000 + |111 ). If the first qubit is flipped from 3-qubit even hamming weight state, we get odd hamming weight state 1/4(|001 + |010 + |100 + |111 ). We get a closer Dicke state with higher fidelity when hamming weight is higher as Dicke state |D N K has N K states and odd/even hamming weight states has 2 N −1 states. When N is higher but K is low, this approximate Dicke state has more non-Dicke states than Dicke states, therefore, the fidelity gets lower. Thus, we expect higher fidelity with increasing K/N . We construct odd hamming weight states |O N for N = 2 to 10 qubits and even hamming weight states |E N for N = 4 to 10 qubits.

B. Lower bounds for Fidelity Estimation
For N-qubit product state |P S N K , we compute lower bound estimations on Dicke state fidelity using the expressions from Section III-B. Similarly, we estimate Dicke state fidelity for odd and even hamming weight states. When K is odd, we estimate a lower bound on fidelity for odd hamming weight state |O N , and for even K, we compute a lower bound on fidelity for even hamming weight state |E N .

C. Experimental Result
We prepare untranspiled circuits for product states |P S N K where 1 ≤ 2K ≤ N ≤ 10, odd hamming weight states |O N where 2 ≤ N ≤ 10 and even hamming weight states |E N where 4 ≤ N ≤ 10. Figure 8 shows preparation of approximate Dicke States for N = 9. We execute the circuits in Quantinuum H1-1 quantum processor to get X,Y , and Z basis measurements similar to Dicke state circuits in Section III-C. We then compute experimental lower bound on Dicke state fidelities for product states and odd/even states. Figure 9 shows the experimental bounds on Dicke state fidelities using the approximate Dicke state circuits. The plot shows bounds on quantum fidelity D N K | ρ |D N K for N = 2 to N = 10 qubits. In each plot, the X axis shows values for K/N and the Y axis shows bounds on fidelity. Each plot contains two lines showing theoretical fidelity of odd/even state and product state preparation with Dicke state. In each plot, the bounds (theoretical and experimental) increase for odd/even state with higher K/N but decreases for product states. Experimental lower bounds are also shown for odd/even state, product state, and Dicke state (from Section III-C) preparation. We find that the experimental lower bounds of the product state closely follow theoretical fidelities in all the plots (N = 2 to 10). However, the experimental lower bounds of odd/even states are not close but follow theoretical fidelities. We also observe that the approximate lower bounds of product state and odd/even state are getting close to Dicke state lower bounds with higher N and for N = 10 and K = 5, the odd/even state lower bounds match with Dicke state lower bound estimations. Thus, we can get a good approximation of lower bounds on Dicke state fidelities from just preparing those approximate states for higher Dicke states.

VI. METHODOLOGY
Our access to the Quantinuum backends was granted through OLCF program of Oak Ridge National Laboratory. In this work, we ran the experiments on both hardware and emulator devices of Quantinuum H-systems using the Quantinuum Python API. Running jobs on Quantinuum backend costs requires H-System Quantum Credits (HQC) that is linear with both number of shots per circuit and number of CNOTs in the circuit. In general, the total cost scales with the product of number of shots and CNOT counts. Our initial goal was to determine number of experiments and number of shots for different state preparation permissible within the HQCs allocated to us. As the Quantinuum machine's are available throughout a calender month and HQCs are allocated monthly, we target to complete experiments over the duration of a couple of months with available hardware credits. First we ran the experiments on Quantinuum H1-2E emulator to get an estimation the confidence interval sizes we can expect from the quantum hardware. To see the effect of the number of circuit execution, we compute the cumulative bounds for fidelity and measured success probability. Similarly, we compute twosided 68% confidence interval below and above the mean of the distribution to check the stability of the bounds with the increasing number of shots. Analyzing the results of emulator experiments, we get an approximation on the meaningful number of different experiments than can be performed on quantum hardware for various states. We limit Dicke state preparation upto N = 10 as Dicke States above N = 11 becomes too costly. For GHZ states, we went upto N = 20 as we ran the experiement on 20-qubit H1-1 device. Additionally, we observe roughly equal confidence interval widths for the bounds across the experiments. Next, we execute the experiments of different state preparation on Quantinuum hardware backends depending on their availability. Dicke state circuits were executed in the 12 qubit H1-2 quantum processor while GHZ state and approximate Dicke state circuits were executed in the 20 qubit H1-1 quantum processor. For all experiments, we directly submit the QASM [6] code using the Python API.
To execute the circuits in the same time span, we submit the jobs in batch mode where we only have to wait the queue time once. However, Quantinuum queue only allows 500 HQC execution in a single batch. Thus, we had to execute subexperiments over a couple of months. Dicke states experiments were run in this way, |D 2 1 to |D 7 2 on May 24th, |D 8 1 to |D 9 2 on June 6th, |D 9 3 to |D 10 3 on June 10th, |D 10 4 on June 17th and |D 10 5 on July 25th of 2022. Both logarithmic depth and linear depth GHZ circuits were executed on August 19th of 2022. Approximate Dicke states (odd/even states and product states) were executed on July 12th of 2022.

VII. CONCLUSION
We have demonstrated good fidelity values for highly entangled quantum states such as Dicke and GHZ states, prepared on Quantinuum H1 ion trap systems. Since direct fidelity estimation through full state tomography does not scale, we have applied and improved existing lower bound techniques, achieving reasonable estimations on NISQ devices for the first time in 15 years.


extfieldsofstudy:
Computer Science
Mathematics

text:
Scheduling of UAV-assisted Millimeter Wave Communications for High-Speed Railway

To exploit richer spectrum resources for even better service quality, millimeter wave (mmWave) communication has been considered for high-speed railway (HSR) communication systems. In this paper, we focus on scheduling as many flows as possible while satisfying their QoS requirements. Due to interference, eavesdropping, or other problems, some flows may not be directly transmitted from the track-side BS. In this paper, we propose an UAV-assisted scheduling scheme which utilizes a UAV to serve as relay for such flows. The proposed scheme also utilize two mmWave bands, one for the BS links and the other for the UAV links. The proposed algorithm aims to maximize the number of flows with their QoS requirements satisfied. Simulations demonstrate that the proposed scheme achieves a superior performance on the number of completed flows and the system throughput over two baseline schemes.

I. INTRODUCTION
With the rapid development of high-speed railway (HSR) in recent years, there is a compelling need for high-quality broadband wireless access services for HSR. There is increasing demand to support data-and bandwidth intensive applications (e.g., multimedia applications) for HSR passengers. The current transmission rate requirement of each train cabin is about 37.5 Mbps according to a recent study, while may continue to grow and reach 0.5 ∼ 5 Gbps in the future [1]. Obviously, these current transmission schemes cannot meet such increasing demands.
To this end, the millimeter wave (mmWave) band from 30 GHz to 300 GHz has abundant spectrum resources and can provide multi-gigabit transmission rates for bandwidthintensive services, such as high-speed data transmission between devices (such as cameras, smartphones, tablets and laptops), compressed and uncompressed high-definition TV, real-time playback, wireless gaming, wireless personal area network, cellular access and wireless backhaul, etc. [2]- [5].
In this paper, we consider the mmWave communications between the HSR and the ground base station (BS). To overcome the severe penetration loss of the train body shell, multiple mobile relays (MRs) are deployed on the rooftop of the train to provide better connections to the track-side BS. The passengers in the train communicates with the BS through radio access links, using the MRs as relays. MRs serve passengers via the APs which are installed in each carriage to avoid the penetration loss and frequent handover. Note that, the communications inside the train can use multiple radio access technologies, such as LTE or WiFi. In fact, the links between the MRs and the track-side BS are the main cause of the capacity bottleneck [6]. Thus, we focus our work on the mmWave communications between the BS and MRs.
In the typical scenario of mmWave HSR systems, there are a certain number of traffic flows to be transmitted. Each flow has its own minimum throughput requirement (i.e., the lowest transmission rate requirement), which is also referred to as its quality of service (QoS) requirement in this paper. Due to the diversity of applications, the QoS requirements of different flows are diverse. Furthermore, the mutual interference from concurrent link transmissions and the presence of eavesdroppers may cause some flows fail to be scheduled directly [7]. Therefore, we propose to deploy a UAV to relay the flows that cannot be directly transmitted from the BS to MRs under certain conditions. We present an UAV-assisted scheduling algorithm that exploits two mmWave bands (including one lower band for links from the BS and a higher band for links from the UAV) to achieve high system throughput. We also consider protecting the mmWave communications from potential eavesdroppers by enforcing the security capacity constraint. The main contributions of this paper are summarized as follows.
• To increase the network capacity of the HSR communication system and satisfy the needs of HSR passengers, we exploit the available spectrum resources of two mmWave frequency bands for the scheduled flows with diverse QoS requirements. The lower one of the two mmWave frequency bands is used for links from the BS, and the higher one is used for links from the UAV. Besides, we consider the existence of potential eavesdroppers, we determine whether a flow can be scheduled based on the security capacity to ensure its communication security. • Due to interference, eavesdropping, or other problems, some flows cannot be directly scheduled from the BS to the MRs. To address this issue, we propose to use a UAV as a relay. With the assistance of the UAV, there is more flexibility for scheduling the flows. We formulate the problem of optimal scheduling as a nonlinear integer programming problem. Moreover, we propose a heuristic UAV-assisted scheduling algorithm to maximize the number of flows with their QoS requirements satisfied within a fixed time. • We evaluate the proposed UAV-assisted scheduling scheme for the mmWave HSR network with simulations. Compared with two baseline schemes, the simulation results demonstrate that the proposed scheme achieves a superior performance on both the number of completed flows and the system throughput. The remainder of this paper is organized as follows. In Section II, we review the related work. Section III presents the system model and assumptions. In Section IV, we formulate the optimal scheduling problem as a nonlinear integer programming problem. In Section V, we propose the relay decision algorithm and UAV-assisted scheduling algorithm. Section VI presents our simulation study of the proposed scheme. We conclude this paper in Section VII.

II. RELATED WORK
There are significant efforts on exploiting mmWave communications for augmented capacity for HSR systems. For the problem of flow forced interruption, utilizing other receivers/nodes as a one-hop or multihop relay to change the original transmission path is the most common solution [8]- [10]. In this paper, we do not use any original nodes to relay blocked flows, but use the UAV to create a new, high-altitude relay node. Due to the superiority of the UAV's open location, superior relay effects can usually be achieved. Many prior works use alternative APs, which are densely deployed, to substitute direct paths from transmitter to receiver [11]- [14]. In this work, we adopt alternative APs in the HSR system, i.e., the MRs are deployed on top of the train to assist transmissions from the BS to passengers in the train.
There have been considerable prior works on scheduling in mmWave networks. Time division multiple access (TDMA) based schemes are commonly used in mmWave networks. Such schemes consider different types of link transmissions, including radio access, backhaul, and so on [15]- [18]. At present, many existing schemes aim to satisfy the QoS requirements of flows, where the scheduling scheme is designed under the premise of meeting the QoS requirement of each flow, and aiming to finally achieve different optimization goals [19], [20]. In this work, the goal is to develop a QoSaware solution to schedule flows in HSR systems, which is assisted by an UAV relay. We consider the cooperation of two mmWave bands to achieve greater transmission capacity. We also consider the corresponding security measures against the communication security risks from potential eavesdroppers.
Recently, there have been a few studies on UAV-assisted relay networks [21]. The scheme of the optimal UAV placement is proposed for the single user or multiple users in the maritime communication, which can increase the capacity of wireless backhaul link between BS and UAV [22]. A cooperative secrecy transmission mechanism is presented for multiple UAV assisted relay selection to maximize the secrecy capacity [23].
In this paper, the transmission research is in HSR systems with the specific communication environment, and the completed link number and system throughput are maximized by scheduling link transmissions. UAV can serve as an energy source and provide RF energy for low power device-to-device (D2D) pairs [24]. Moreover, UAV enables the simultaneous wireless information and power transfer (SWIPT) system by jointly optimizing time switching-based relaying (TSR) and the optimization of its hovering location for better channel quality to maximize the network throughput [25]. These works usually assume the Rician fading channel as well as the lineof-sight (LoS) path between the UAV and each receiver [26]. To better model the actual channel state in HSR systems and to provide a better fit to real measurements observed in LoS scenarios, the κ-µ distribution was proposed in [27]- [30] as a general fading channel model where the Rician, Nakagamim, and Rayleigh distributions are its special cases. In this paper, we adopt the α-β path loss model which is more in line with the path loss measurements in the HSR scenario. Besides, there are works which have studied the application of UAVassisted communications in HSR system. Two UAV-assisted relaying transmission schemes are presented to allocate limited power for moving source and UAV, which are to save signaling overhead during the entire transfer process [31]. In this work, we proposed QoS-aware solution to maximize the number of scheduled flows. In addition, we allocate the time resource for two types of links and consider the cooperation of two mmWave bands for the links from BS to UAV and the links from UAV to MRs.
It is worth noting that these existing related works have not considered utilizing UAV to relay traffic flows which cannot be directly transmitted from the BS to MRs. In this paper, we integrate the dual-band cooperation of mmWave communications into the scheduling problem, aiming to maximize the number of scheduled flows. We also aim to satisfy the QoS requirements of each flow and avoid the potential risks caused by eavesdroppers. We will show through analysis and simulation that the proposed scheme can achieve a superior performance compared to the baseline schemes.

A. System Model
In this paper, we consider a single cell scenario of the mmWave HSR communication system, as shown in Fig. 1. A fixed BS is located beside the railway track, and a UAV is hovering above at a fixed position within the coverage of the track-side BS.
which are within the coverage range of the BS and UAV. The BS links are on frequency band f 1 , indicated by the black arrows in Fig. 1. The communications between the UAV and MRs are on frequency band f 2 , which are shown as red arrows in the figure. In addition, there is an eavesdropper, termed Eva, which is trying to overhear the communications. Without loss of generality, we assume Eva can only eavesdrop on the BS through frequency band f 1 .
There are F traffic flows between the BS and MRs, which need to be scheduled in the system. The UAV provides relay service in the air, thus each flow may be received from the BS directly or be relayed by the UAV for the destination MR. MRs can communicate with the track-side BS or the aerial UAV through radio access links. According to the selection of the transmission methods of the flows, some links are scheduled for transmission each time. In this paper, a link represents the actual connection of direct communications between two endpoints in the system.
In this paper, we consider the scenario where a train runs at a constant velocity v = 300km/h. In this case, the speed of the train is in line with reality. Although the train is moving at a high speed, the large bandwidth of mmWave allows the train to receive a large amount of data in a very short time. Therefore, the pressure of handover problem in high-speed railway communication can be greatly alleviated. We consider the downlink transmissions. The auxiliary transmission of UAV can effectively alleviate the link blockage problem caused by obstacles on the ground. All the devices in the system are assumed to operate in the half-duplex mode, i.e., they can only transmit or receive data at a time. Each MR is equipped with one steerable directional antenna, so it can only receive content from the BS or UAV at a time. The UAV is equipped with one directional antenna for receiving and N Ut antennas for transmitting. The BS is equipped with N Bt directional antennas for transmitting. In the proposed system, time is divided into a series of non-overlapping super frames. Each super frame consists of a scheduling phase and a transmission phase. The scheduling phase is the time duration to collect link requests and their QoS requirements. The transmission phase consists of M equal time slots (TS). For clarity of exposition, we use TSs to measure the transmission time. For reliable transmissions, the control messages and BS requests can be transmitted with 4G technologies, including the source/destination information and the QoS requirements of each link [32]. The proposed scheduler is to allow multiple links be transmitted concurrently using spatial-time division multiple access (STDMA) [33].

B. Antenna Model
As in prior works, directional antenna is necessary for mmWave communications [34]. We assume that the BS and UAV are both equipped with multiple directional antennas, which means they can send to multiple links simultaneously. However, each MR can receive from only the BS or UAV at a time. We adopt multiple antennas to form multiple directional beams at the BS and UAV for transmitting multiple data flows simultaneously [35]. The maximum numbers of beams formed simultaneously at the BS and UAV are limited by the numbers of theirs equipped transmitting antennas.
In this paper, we adopt the realistic directional antenna model which is mentioned in [36]. All directional antennas can steer their beams for the maximum directivity gain [37], [38]. The perfect beam alignment for the desired link has the maximum directivity antenna gain of G 0 . The antenna model can be expressed as where θ is an angle within the range [0 • , 180 • ]. G 0 is the maximum antenna gain, G 0 = 10 log(1.6162 sin(θ −3dB /2)) 2 . θ −3dB is the angle of the half-power beamwidth. θ ml is the main lobe width, θ ml = 2.6 · θ −3dB . G sl is the sidelobe gain, The directional antenna model is for the antenna gain calculations of the links.

C. Transmission Model
The directional path loss in the target environment is studied via the impacts of antenna beam and propagation conditions. After evaluating the different existing models, the α-β model is selected as the path loss model for quantitative analysis in this work [39]. For a link i, the path loss between its transmitter t i and its receiver r i is where α and β are constants, d tiri is the distance between t i and r i , and X σ is a zero mean Gaussian random variable with a standard deviation σ.
Compared with line-of-sight (LOS) transmissions, non-lineof-sight (NLOS) transmissions suffer from much higher attenuation in the mmWave band [40]. Therefore we only consider the directional LOS links. For a link i, the received signal power at receiver r i from transmitter t i is where G ti is the transmitter antenna gain of link i, G ri is the receiver antenna gain of link i, P t i is the transmit power at the transmitter of link i, and P L(d tiri ) is calculated by (2). Since we allow multiple links transmit concurrently by the STDMA transmission scheme, there is mutual interference between concurrent links using the same frequency band. For co-band concurrent links i and j, the interference at r i from t j can be calculated as where G tj is the antenna gain at the transmitter of link j, G ri is the antenna gain at the receiver of link j, which is calculated as in (1), P t j is the transmit power at the transmitter of link j, and P L(d tj ri ) is calculated by (2). According to Shannon's formula, the transmission rate of link i is given by where η describes the efficiency of the transceiver design, which in the range of (0, 1); W i is the transmission bandwidth of link i; and N 0 is the one-sided power spectra density of the white Gaussian noise. We denote R t i as the actual transmit rate of link i in time slot t, which is calculated by (5). Then the achieved throughput of link i in this time slot can be obtained as where T s is the duration of the scheduling phase, and ∆T is the duration of each TS. The denominator in (6) is the duration of a super frame.

D. Secrecy Capacity Model
In this paper, we consider the existence of an eavesdropper, which can overhear the BS transmissions on frequency band f 1 . It is important to ensure secure communications and prevent information leakage. The secrecy capacity, an important physical layer security performance metric, is adopted to determine the feasibility of the schedule scheme in this work. Let C M and C E denote the Shannon capacities of the main channel (C M = R i if link i is attacked) and the eavesdropping channel, respectively. The secrecy capacity is given by [41] where C M and C E are both calculated by (5). The position of the eavesdropper is constantly changing with the train moving. But the eavesdropper is fixed in a certain position of the train and is stationary relative to the train. Hence, compared with other randomly moving eavesdroppers, the location of the eavesdropper in this case is easier to estimate. The eavesdropper is not completely passive, and sends eavesdropping messages to the unknown receiver. The train is slow when it starts to run and can be regarded as stationary. At this time, we adopt the angle of arrival (AOA) positioning method to determine the specific location of the eavesdropper on the train [42]. The BS and MRs detect the signals sent by the eavesdropper to obtain the angles of incidences. Then the location of the eavesdropper can be determined by these angles of incidences. However, determination of the eavesdropper location facilitates the calculation of the security capacity.

IV. PROBLEM FORMULATION AND ANALYSIS A. Problem Formulation
In this paper, the goal is to accommodate as many flows as possible through effective scheduling. In fact, different flows carry different types of service. Thus each flow f has its QoS requirement q f in the form of its required minimum throughput. We define a binary variable δ f to indicate whether flow f is scheduled successfully. If that is the case, we have δ f = 1; otherwise, δ f = 0. Let there be F flows. Then the objective function of the scheduler is given by For link i, a binary variable a t i is defined to indicate whether the link is scheduled for transmission in time slot t. If so, a t i = 1; otherwise, a t i = 0. For the links generated from relayed flows, they are all transmitted to or received from the UAV. The following constraint enforces the half-duplex transmission condition.
The constraint means the link to the UAV and the link from the UAV can't be scheduled concurrently. Each MR and the UAV can only receive a single link at a time. In other words, the receivers of any two concurrent links i and j must be different. Thus, we have the following constraint.
Each link i has its QoS requirement q i . The actual achieved throughput of each flow must satisfy its QoS requirement. For link i, the QoS constraint is as follows.
Regardless of whether it be the BS or UAV, the numbers of links can be simultaneously transmitted are limited by the numbers of its transmit antennas. In time slot t, the set of links transmitted from the BS is defined as S t B , and the set of links transmitted from the UAV is defined as S t U . Then we have the following hardware constraint.
Some flows are relayed by the UAV and then forwarded to the target MRs. The total throughput of the links transmitted from the BS is larger than the sum throughput of the links transmitted from the UAV.
where q a i and q a j are calculated as in (6). Due to the existence of an eavesdropper, the confidentiality of transmissions must be guaranteed. The eavesdropping channel between a legitimate pair is a degraded version of the main channel [41]. In time slot t, the main channel capacity and the eavesdropping channel capacity of link i are denoted by C i M and C i E , with the following constraint [43].
The problem of maximizing the amount of scheduled flows in in each time slot can be formulated as follows.
It can be seen that Problem P1 is a nonlinear integer programming problem. There are also complex nonlinear terms in the constraints. The problem is actually NP-hard [19], [44]. Thus, we aim to develop heuristic algorithms to obtain competitive solutions with low computational complexity.

V. THE UAV-ASSISTED SCHEDULING SCHEME
In this section, we propose an UAV-assisted scheduling scheme. With this scheme, each flow can choose from two transmission methods: one is to directly transmit from the BS to the target MR, and the other is to use the UAV as a relay. The proposed scheme first determines the transmission method for each flow and then schedules the links to transmit, aiming to maximize the number of scheduled flows.

A. The Relay Decision Algorithm
The relay decisions of the flows are based on their QoS requirements and the transmission efficiency. The flows that need to be scheduled are determined based on the transmitterreceiver ranges and their QoS requirements. Then we schedule the flows in each time slot for transmission.
We first consider whether the QoS requirements of the flows can be satisfied. Then we make a selection for each flow between being directly transmitted and being relayed. The set of the links which select to be directly transmitted from BS to MR is defined as S 1 , and the set of the links relayed by the UAV is defined as S 2 .
For each flow f , if it is directly transmitted from the BS to MR in the f 1 band, then this decision will allow one link l 1 from BS to MR to be scheduled. The achieved throughput q a l1 of this link is computed as in (6). If flow f is to be relayed by the UAV in the f 2 band, then there are two links l 2 and l 2 to be scheduled: one from the BS to UAV and the other from the UAV to MR. The achieved throughput of these two links, q a l2 and q a l 2 , can be calculated by (6). Since it is uncertain which links are transmitted simultaneously, the above throughput calculated by (6) does not consider the interference from other links. We can judge whether the two decisions (relay and no relay) can successfully transmit the flows by comparing the calculated throughputs and the QoS requirements. For flow f , the choice between direct transmission and relayed transmission can be decided as follows.
1) If min{q a l1 , q a l2 , q a l 2 } ≥ q f , flow f can be scheduled by both choices.
2) If q a l1 ≥ q f and min{q a l2 , q a l 2 } < q f , flow f can only be directly transmitted from the BS.
3) If q a l1 < q f and min{q a l2 , q a l 2 } ≥ q f , flow f need to be relayed by the UAV. 4) If q a l1 < q f and min{q a l2 , q a l 2 } < q f , flowf should be abandoned. It should not be considered in the later scheduling process. Since the eavesdropper is stealthily listening on the BS transmissions in the f 1 frequency band, the confidentiality of transmission should be considered. The location of the UAV ensures that the BS-UAV links can satisfy constraint (14). Thus the flows relayed by the UAV are less likely of being eavesdropped. For flows which are directly transmitted from the BS in the f 1 frequency band, we judge whether a direct transmission is feasible by evaluating (14).
To further compare the two transmission choices, we make a preliminary estimation of the transmission time of all the flows in S all = S 1 ∪ S 2 . The number of estimated time slots of flow f can be obtained as follows.
Here R i is the estimate rate without considering other links' interference at this time. The remaining flows select a better transmission choice by comparing the number of TSs required by the schedule. The relay decision algorithm is summarized in Algorithm 1. First, we remove the flows whose QoS requirements cannot be satisfied by either of the two transmission choices in Line 2. In addition to the flows in set D , the other flows are divided into three categories based on whether their QoS requirements can be satisfied under different transmission choices in Lines 5-24. In Lines 5-6, the flows with their QoS requirements satisfied by being relayed by the UAV are added to S 2 . In Lines 9-13, the transmission methods for the other flows are determined, with their QoS requirements satisfied by direct transmission from the BS and satisfying the secrecy capacity condition. And these flows will be put in S 1 , the set of flows to be directly transmitted from the BS. In Lines 14-24, flows whose QoS

Algorithm 1: The Transmission Method Selection Algorithm
Input: S all , S1 = ∅, S2 = ∅, QoS requirements of flows, locations of the BS, UAV, eavesdropper and MRs Output: S1, S2 1 Calculate q a l 1 , q a l 2 and q a l 2 for each flow; 2 Remove D = {f |q a l 1 < q f and min{q a l 2 , q a l 2 } < q f then 9 if C l 1 E can satisfy (14) then 10 S1 = S1 ∪ f ; requirements can be satisfied by both transmission methods will be divided into three categories according to the secrecy capacity condition and transmission efficiency. If the secrecy capacity condition cannot be satisfied, it will be put in the set of flows relayed by the UAV, as in Lines 15-16. Finally, the remaining flows select their transmission methods by comparing the required transmission time of different choices as in  For the complexity of the transmission method selection algorithm, the number of the iterations of the outer for loop in Line 3 is |S all |, where |S all | in the worst case is O(F ). Hence, the computational complexity of the algorithm is O(F ).

B. The UAV-assisted Scheduling Algorithm
Next we present a transmission scheduling algorithm. This scheduling algorithm is developed on the basis of the above transmission method selection algorithm.
We select the appropriate links for transmission in each TS based on the principle of maximizing the number of transmitted flows. For ease of exposition, we first introduce the simultaneous transmission requirements of links and the priority of link scheduling. Regarding the contention among links, we consider the three cases as follows.
1) Due to half-duplex transmissions, the link from the BS to the UAV and the link from the UAV to an MR cannot be scheduled concurrently. 2) If the simultaneous transmissions of two links will cause at least one of them to fail to satisfy the QoS requirement, there is a conflict between the two links. (14), the links transmitted from the BS in the f 1 frequency band that do not meet the safe capacity threshold will not be scheduled simultaneously with other links. A flow f relayed by the UAV will require two links l 2 and l 2 to be scheduled, i.e., l 2 from the BS to UAV and l 2 from the UAV to MR. We define the set of links scheduled in the f 1 band as S f1 , and the set of links scheduled in the f 2 band as S f2 . Link l 2 should be scheduled before l 2 . Thus we first schedule link l 2 ∈ S f1 , and then link l 2 is generated in S f2 to wait for being scheduled.

3) As stated in
Then we consider the transmit order of the links in S f1 or S f2 . A parameter is defined to prioritize the scheduling order of the links, which is the inverse of a link's number of required TSs in a frame to satisfy its QoS requirement. For link i which is generated from flow f , the priority value is given by .
In fact, the less time the transmissions consume, the more flows can be completed in a fixed time interval. Prioritizing the flows that need fewer TSs can transmit more flows in a frame. Therefore, we calculate the priority values of links in S f1 and S f2 as (17), and then transmit the links in the descending order in the scheduling process. Since there are always new links joining S f2 , the scheduling order of the links in the set needs to be frequently updated.
We define the set of the links which are to be transmitted from the BS to UAV as S BU , and a parameter N f is used to indicate the number of completed flows. The QoS-aware scheduling algorithm is summarized in Algorithm 2. First, the input and output of the algorithm are specified. Then, we configure sets S f1 and S BU in Lines 1-2. In Lines 7-14, we find the links in S f1 that can be scheduled simultaneously with other selected transmitted links. In Lines 15-22, the simultaneous scheduled links in S f2 are determined. In Lines 23-40, we obtain the remaining throughput demands of all the scheduled links after this time slot, and examine whether there are flows which have finished transmission in this TS.
The outer for loop of the scheduling algorithm has O(M ) iterations. The for loop in Line 8 has |S f1 | iterations, and |S f1 | in the worst case is O(F ). The for loop in Line 18 has |S f2 | iterations with the worst case value O(F ). Moreover, the for loop in Line 23 has |SS| iterations, and its worst case value is also O(F ). The for loops in Line 8, Line 18 and Line 23 are parallel. Therefore, the complexity of the QoS-based scheduling algorithm is O(M F ). Such low complexity of the algorithm makes it suitable for practical implementation.

VI. PERFORMANCE EVALUATION
In this section, we evaluate the performance of the UAVassisted scheduling algorithm. The flow transmissions of the Algorithm 2: The UAV-assisted Scheduling Algorithm Input: S1, S2, SBU = ∅, S f 1 = S1, S f 2 = ∅, N f = 0, Calculate the priority values of the links in S f 1 ; 5 Sort the links in S f 1 in increasing order by priority value; scheme involve two bands: f 1 in the 28GHz band and f 2 in the 60GHz band. We also compare the performance of the proposed scheme with several existing baseline schemes.

A. Simulation Setup
We consider the single cell scenario of the mmWave HSR system. There are 24 MRs which are evenly distributed on top of the train. The flows which need to be scheduled are all transmitted from the BS. Then the flows are generated by randomly selecting the receivers among all the MRs. The train has a total of eight cars with a total length of 200m. Three MRs are deployed in each car. One eavesdropper is randomly placed on top of the train. The BS and UAV can reach a communication range of hundreds of meters, and all the links in the system are transmitted within the allowable communication ranges of their transceivers. The number of flows which need to be transmitted between the BS and MRs is smaller than the number of MRs.
Since the MRs move at the high speed with the train, we need to update theirs location information periodically. We collect and update the location information of MRs every 2000 time slots in the process. During this period of time, the change distances of MRs do not exceed 3m, which does not affect the data transmissions.
Besides, each carriage of the train has 50-100 seats, then each MR has to deal with the data traffic requirements of 10-20 passengers at most. The flow of each MR request may be a combination of the traffic flows of multiple passengers. Thus, the QoS requirement values of scheduled flows fluctuates in a wide range. We assume that the QoS requirement of each flow is uniformly distributed between 10Mbps and 500Mbps. Other parameters are presented in Table I [14], [35], [39]. For the path loss model, the specific expression is given in (2) and the parameter settings are from [39], which are summarized in Table II. For the performance evaluation study, the following performance metrics are used.
• Number of completed flows: the number of scheduled flows with their QoS requirements satisfied. If the QoS requirement cannot be met, the flow that has been scheduled will not be counted as a completed flow. • System throughput: the total throughput of the entire network system. This metric is the sum of the average throughputs of all the flows carried in the network. • Total amount of time slots: the total number of time slots used to transmit all the scheduled flows. The actual value of this metric must be smaller than the total number of TSs in the frame. For comparison purpose, we implement the state-of-theart QoS-aware Concurrent Scheme [17] and Maximum QoS-aware Independent Set (MQIS) based scheduling algorithm [19] as baseline schemes.
• QoS-aware Concurrent Scheme [17]: this is a concurrent scheduling scheme aiming to satisfy the QoS requirements of flows and to maximize the number of flows. If a flow cannot meet the current scheduling conditions, the scheme cannot relay the flow and has to give up scheduling it. • MQIS [19]: this is the scheduling algorithm based on the concept of maximum QoS-aware independent set proposed in [19]. The algorithm divides flows into multiple independent sets. Only after all the flows in one set are scheduled transmitted, the flows in another set can begin to be scheduled.

B. Comparison with Existing Schemes
In Fig. 2, we plot the number of completed flows under different numbers of requested flows. The horizontal distance between the UAV and BS is set to 150m, and the number of time slots is set to 8,000. With the increased number of requested flows, we can see that the curves for the three schemes are all rising. The more requested flows in the HSR network, the more flows can complete their scheduling. For the proposed UAV-assisted scheduling scheme, some flows which cannot be directly transmitted from the BS can choose to be scheduled by UAV relaying. Therefore, the proposed scheme achieves a superior performance with respect to the number of completed flows compared with the other two schemes. When the number of flows is 18, the proposed scheme outperforms the QoS-aware concurrent scheme with a gain of 35.0% and MQIS with a gain of 32.4%.
In Fig. 3, we plot the system throughput under different numbers of requested flows. The horizontal distance between the UAV and BS and the number of time slots are the same as in Fig. 2. From the figure, we can see that the system throughput of the three schemes are all increasing with the increased number of requested flows. For each flow, our proposed scheme can decide whether to use UAV relaying according to the principle of transmission rate maximization. The proposed scheme achieves the largest system throughput among the three scheme. It can also relay some flows which cannot be directly transmitted from the BS. The QoS-aware scheme and MQIS have similar performance on the system throughput. When the number of requested flows is 18, the  proposed UAV-assisted scheduling scheme outperforms the QoS-aware concurrent scheme with a gain of 64.6% and MQIS with a gain of 61.2%.
In Fig. 4, we plot the system throughput under different numbers of requested flows. The horizontal distance between the UAV and BS and the number of time slots remain the same as previous simulations. From this figure, we can see that the total number of slots of the three schemes are all increasing as more flows are requested. Although the proposed scheme can schedule the most flows among the three schemes, as shown in Fig. 2, it can be seen from Fig. 4 that the total numbers of time slots needed by the three schemes are very similar. This proves that the proposed UAV-assisted scheme has a superior performance on transmission efficiency compared with the QoS-aware concurrent scheme and MQIS.
In Fig. 5, we plot the number of completed flows under different numbers of time slots in each frame. The horizontal distance between the UAV and BS is set to 150 m, and the  number of requested flows is set to 18. When the number of time slots begin to increases, the curves of the three schemes are all rising. When the number of time slots increases to a certain value, the numbers of completed flows of three schemes are all becoming steady. Therefore, the number of time slots does not need to be increased to a very large value; an appropriate value of time slot number in each scheme can ensure a satisfactory performance.
In Fig. 6, we plot the system throughput under different numbers of time slots. The horizontal distance between UAV and BS and the number of requested flows are the same as that in Fig. 5. From the figure, we can see similar trends for the three schemes as in Fig. 5. With the increased number of time slots, we can see that the curves for the three schemes are all rising. The more time slots, the more system throughput can achieve during scheduling. The three curves start to stabilize when the number of slots is increased over 8,000. Hence, the number of TSs is set to 8,000 in the other simulations.  In Fig. 7, we plot the number of completed flows under different horizontal distances between UAV and BS, to demonstrate the impact of the UAV location on performance. The number of requested flows is set to 18, and the number of time slots is set to 8,000. With the increased horizontal UAV-BS distance, the number of completed flows increases first and then declines. Thus a suitable location of the UAV has a significant impact on the system performance. When UAV-BS distance is 150 m, the number of completed flows reaches the maximum value. Thus, the horizontal distance between the UAV and BS is set to 150 m in all the other simulations.

VII. CONCLUSIONS
In this paper, we considered the problem of scheduling flows with diverse QoS requirements in the scenario of mmWave HSR communications. To maximize the number of flows while satisfying their QoS requirements, we proposed the UAVassisted scheduling scheme to deploy a UAV relay to assist the scheduling of flows. Our simulation study showed that the proposed UAV-assisted scheme outperformed two baseline schemes on the number of completed flows, achievable system throughput, and transmission efficiency.


Feature Selection for Cross-Scene Hyperspectral Image Classification via Improved Ant Colony Optimization Algorithm

Hyperspectral images (HSIs) generally contain a large amount of spectral bands (features), and the redundant information in them will cause the Hughes phenomenon in the classification process. And feature extraction and feature selection are the two main existing methods to effectively reduce the redundancy of spectral information in the field of HSIs classification. Compared with feature extraction methods, feature selection methods can preserve most of the features of the original HSIs data without losing their valuable details. However, most existing feature selection methods based on single scene (domain) perform poorly in some scenes (domains) with insufficient labeled samples. Therefore, how to adopt an efficient feature selection method to select the optimal feature subsets of source scene and target scene and use the sample information of source scene to assist in the classification of target scene so as to improve the classification accuracy of images in the target scene as much as possible is still very challenging. In order to solve the above problem, this paper proposes a new cross-scene algorithm: Improved Ant Colony Optimization Algorithm-Based Cross-Scene Feature Selection Algorithm (IMACO-CSFS). In order to obtain more accurate feature subsets of the two scenes, IMACO-CSFS proposes a priority sorting-based ant colony strategy to make the subsequent search process focus on the global optimal solution (optimal feature subset) found in the previous iteration. In addition, in order to further accelerate the convergence speed of the global optimal solution, an ant colony strategy based on elite ants is proposed in IMACO-CSFS to more efficiently obtain the optimal feature subsets of the two scenes for training the classifier. Furthermore, this paper simultaneously considers overall classification accuracies of the optimal feature subsets for both scenes and dynamically adjusts their scale to ensure the consistency of the selected features between the two scenes, attenuating the effect of spectral shift and achieving the higher image classification accuracy in the target scene. Experimental results on three cross-scene HSI data pairs demonstrate that IMACO-CSFS is superior in cross-scene feature selection and cross-scene HSIs classification.

I. INTRODUCTION
HYPERSPECTRAL images (HSIs) typically comprise 26 dozens or hundreds of spectral bands, which contain a lot 27 of spatial and spectral information about ground objects [1], 28 so they have been widely used in agriculture, forestry, 29 The associate editor coordinating the review of this manuscript and approving it for publication was Gerardo Di Martino . environment, ecology and other fields [2], [3]. However, 30 problems of redundant features and high data dimensions 31 still exist in HSIs [4], which may cause Hughes phenomenon 32 in the process of HSIs classification (Hughes phenomenon 33 refers to the phenomenon that the performance of the classi-34 fier ''increases first and then decreases'' with the increase 35 of the feature/band dimension involved in the operation, 36 during the analysis and processing of HSIs). Therefore, 37 tion of the target scene, the problem of spectral offset will 92 basically be encountered because of the different material 93 composition of the same type of ground object in different 94 time and size and the different nonlinearity of sensors to 95 obtain the data of the two scenes. Therefore, it is particu-96 larly important to reduce spectral offset. At present, some 97 methods such as cross-perspective learning [27] and transfer 98 learning [28] have been proposed to solve such problems. 99 For example, Literature [29] calculates the distance of condi-100 tional distribution between two scenes in a regenerated kernel 101 hilbert space based on a new dataset shift measure. Based 102 on deep adaptive network (DAN), multi-kernel maximum 103 mean discrepancies (MK-MMD) is used in [30] to align the 104 deep features of source and target scenes. In [31], I-ReliefF 105 (Iterative-ReliefF) was applied in the field of cross-scene 106 classification, and the experiments on the newly proposed 107 CDIRF (Cross-Domain Iterative ReliefF) method on three 108 cross-scene datasets validate its superiority in cross-scene 109 feature selection and reduction of spectral shift between two 110 scenes.

111
Ant colony optimization algorithm (ACO) is a classi-112 cal heuristic random search algorithm inspired by nature, 113 because its positive feedback mechanism can help ant colony 114 find the advantage of optimal solution in a short time, ACO 115 has been widely used in image processing, feature/band selec-116 tion for HSIs data and other fields [32], [33]. However, 117 it is not suitable to be directly used to solve the problem 118 of feature/band selection due to its shortcomings such as 119 slow convergence rate and significant decrease of population 120 diversity in the process of cyclic iteration [34], [35]. Just 121 as the Improved Ant Colony Optimization Algorithm-Based 122 Band Selection Algorithm (IMACO-BS) newly proposed in 123 literature [34] is based on ACO and uses a new pre-filter to 124 reduce the number of candidate bands(n), accelerating the 125 convergence of the standard ant colony optimization algo-126 rithm. In addition, the adaptive information updating strat-127 egy was introduced in IMACO-BS to avoid the ant colony 128 falling into local optimum. Although it has the advantages of 129 quick ant colony convergence speed and difficulty in falling 130 into local optimal solution, this method is only suitable for 131 the selection of the optimal band subset of a single scene. 132 With the development of HSIs classification technology, it is 133 becoming more and more important to use the information of 134 two scenes to deal with the problem of high-dimension cross-135 scene HSIs classification. As Zhang et al. [31] extended the 136 iterative ReliefF method to a cross-scene HSIs classifica-137 tion method (CDIRF), CDIRF (Cross-Domain Iterative Reli-138 efF) extracts effective features from redundant spectral bands 139 for subset evaluation, improving the the classification accu-140 racy of the target scene. Therefore, IMACO-BS is applied 141 and improved in this paper to solve the problem of fea-142 ture selection and image classification of cross-scene HSIs, 143 and the new algorithm is called cross-scene IMACO-FS 144 (IMACO-CSFS). 145 At present, there are also some related works that apply 146 ACO to the HSIs classification and feature selection, and 147 some achievements have been achieved. For example: Lit-148 erature [36] proposes a remote sensing image classification 149 VOLUME 10,2022 technology based on the optimal support vector machine and  search process can focus on the global optimal solution found 205 so far. The strategy can obtain more accurate feature subsets 206 of source and target scenes than other cross-scene feature 207 selection methods; 208 2) In view of the shortcomings of the original IMACO-FS 209 that the optimal solution converges too slowly in each itera-210 tion process, this paper proposes an ant colony strategy based 211 on elite ants: in order to make the current optimal solution 212 more attractive to the ants in the next iteration, this strategy 213 gives the optimal solution an extra pheromone after each iter-214 ation. The ant colony strategy based on elite ants can further 215 accelerate the convergence speed of the global optimal solu-216 tion of IMACO-FS, so as to obtain the optimal feature subsets 217 in the source scene and the target scene more effectively for 218 training the classifier; 3) In this paper, the feature selection method IMACO-FS 220 based only on a single scene is successfully applied to the 221 fields of cross-scene HSIs feature selection and cross-scene 222 image classification, which further improves the update strat-223 egy of pheromone in the ant colony algorithm, weakens the 224 influence of spectral drift on cross-scene image classifica-225 tion, and greatly improves the feature selection accuracy and 226 image classification accuracy of the target scene.

227
The rest of this paper is as follows: Section II intro-228 duces basic knowledge of ACO and IMACO-FS. Section III 229 introduces a cross-scene feature optimal selection algo-230 rithm, IMACO-CSFS, which extends the original single-231 scene-based IMACO-FS algorithm to make it suitable for 232 cross-scene feature selection. In Section IV, several algo-233 rithms are compared on three cross-scene HSIs data pairs and 234 the experimental results are analyzed, proving the superiority 235 of the proposed IMACO-CSFS algorithm. Finally, Section V 236 summarizes the work of this paper.

238
A. ANT COLONY OPTIMIZATION ALGORITHM (ACO) 239 ACO is a heuristic global optimization algorithm used to 240 find the optimal path. It was first proposed by Marco Dorigo 241 in his doctoral dissertation in 1992. This algorithm has the 242 characteristics of distributed computing, positive feedback 243 of information and heuristic search [49]. ACO imitates the 244 behavior of ants looking for the optimal path in the process 245 of finding food. This principle of finding the optimal path is 246 shown in Fig. 1. While ants travel, they leave behind a volatile 247 secretion called pheromone. Ants can sense the presence of 248 this substance during foraging and walk along paths with 249 higher pheromone concentrations. And each passing ant con-250 tinues to leave pheromones on the way, forming a mechanism 251 similar to positive feedback. Through this mechanism, ants 252 can finally find the best action path [50]. In IMACO-FS, the set A = {A 1 , A 2 , . . . , A n } repre-257 sents n graph nodes (a node represents a HSI feature), 258 represents a subset of features selected by the ants. Each    [40], which is set to n/2 in this paper).

299
The calculation method of the heuristic expectation η ij (t) 300 is shown in formula (1), where O ij is the OA calculated by the 301 SVM with feature i and feature j.
The pheromone calculation formula between feature i and 304 other n-1 features is shown in formula (2).
where O ij is the OA achieved by SVM with feature i and 307 feature j, O max and O min are the maximum and minimum OA 308 achieved by SVM between feature i and the remaining n-1 309 features. The transformation mode of ants' states is as follows: where q 0 = 1 − e −1/s (s = 1, 2, . . . , I M , s is the number 314 of iterations), q is a random decimal between 0 and 1, p k ij (t) 315 represents the probability that ant k selects the next feature 316 node j at time t, which can be calculated as follows: where α is the pheromone, β is the heuristic information, T 319 is the sum of I M iteration times, τ ij (t) and η ij (t) represent 320 the pheromone concentration and heuristic expectation on the 321 edge (A i , A j ) at time t, respectively, and tabu table k stores 322 the feature nodes that ant k has visited.

324
After each ant completes one cycle (iteration), the pheromone 325 concentration on the path is updated as follows: where τ ij (t + 1) and τ ij (t) respectively represent the 328 pheromone concentration on the edge (A i , A j ) at time t +1 and 329 time t + 1, and τ ij k is the pheromone newly left by ant k on 330 the edge (A i , A j ), the calculation method is as follows: where Q is the initiavalue of the pheromone secreted by the 333 ants at the beginning of an iteration, which is a constant, O max 334 is the maximum OA achieved by the SVM between feature i 335 and other o − 1 features in feature subset F k , t represents 336 time t, and T is the sum of I M iteration times. 337 VOLUME 10, 2022

338
The common HSIs feature selection method generally only 339 selects the optimal feature based on a single scene. The

365
The optimal path found by ant k so far will obtain additional 366 pheromone τ ij * , and the edges passed by the remaining 367 P − 1 ants can also obtain a certain amount of pheromone 368 τ ij k (k = 1, 2, . . . , P − 1). The ant colony strategy based on 369 priority sorting can further search around the current optimal 370 solution δ, thus obtaining the more accurate optimal feature 371 subsets M s and M t in the source and target scenes.

372
In addition, since original IMACO-FS has the disadvantage 373 that the convergence speed of optimal solution δ is too slow in 374 each iteration process, this paper proposes an ant colony strat-375 egy based on elite ants: in order to make the current optimal 376 solution δ more attractive to ants in the next cycle, this strat-377 egy will give the optimal solution an additional pheromone 378 τ ij * * after each iteration, and the optimal solution δ * is 379 completed after I M iterations is the global optimal solution, 380 and the ant that finds this solution is an elite ant. The elite 381 ants-based ant colony strategy further accelerates the conver-382 gence speed of the global optimal solution δ * in the original 383 IMACO-FS, thereby obtaining the optimal feature subsets M s 384 and M t in the source and target scenes more efficiently.

385
The newly proposed pheromone update formula can be 386 calculated as follows: τ ij k is the update of the pheromone content on the 390 path by the P − 1 ants according to their rank k, the calculation method of τ ij k is shown in formula (8), the calculation method of τ ij * is shown in formula (9), and the calculation 393 method of τ ij * * is shown in formula (10).
Among them, Q is the initial value of pheromone secreted   (12): 435

Input:
The labeled samples of source scene T s , The labeled samples of target scene T t .

Output:
Objective function δ * . 1: Train a classifier using randomly selected 1/10 of the target scene labeled samples and calculate D t by using equation (12) (see Table 1). 2: Train a classifier using randomly selected 1/10 the number of labeled samples from the source scene and calculate O s (see Table 1). 3: Train a classifier using randomly selected 1/10 the number of labeled samples from the target scene and calculate O t (see Table 1). 4: Use (11) to compute δ * . 5: return δ *

Input:
Parameter µ, Pheromone α, The labeled samples of target scene T t , The labeled samples of source scene T s .

Output:
Objective function δ * . 1: Initialize pheromone α using (1).   Table 2, which is very suitable for the experiment for 468 classification of cross-scene HSIs. For better comparison, 469 we select Huston 2018 images and Huston 2013 images with 470 a size of 935 × 209 pixels in our experiments. Fig. 4 shows 471 their false color map and ground truth map (the color of each 472 type of ground object is marked next to the ground truth map). 473 The second data pair is University of Pavia (source scene) 474 and Center of Pavia (target scene). The Pavia dataset was 475 acquired by the DAIS hyperspectral sensor. The size of the 476 PaviaU (University of Pavia) image is 243 × 243 pixels and 477 the number of bands is 72, while the size of the PaviaC (Cen-478 ter of Pavia) image is 400 × 400 pixels and the number of 479 bands is also 72. Although they have different pixels, they 480 have 6 identical land cover classes, and the number distribu-481 tion of their labeled samples is shown in Table 3    It performs I-ReliefF on the source scene and the target 512 scene respectively. CDIRF proposes a cross-scene feature 513 weight update rule, which uses two distance measures in 514 experiments. They are CDIRF1 using the absolute dis-515 tance and CDIRF2 using the squared Euclidean distance, 516 and the latter CDIRF2 performs better in the experi-517 ments as a whole. To verify the performance of fea-518 ture selection, CDIRF chooses SVM with radial basis 519 function (RBF) as the classifier. The parameters σ ∈ 520  It performs the hybrid whale optimization algorithm with 527 simulated annealing (WOASA) on the source and target 528 scenes, respectively. CDWOASA selects feature subsets that 529 are discriminative and scene-invariant by using information 530 VOLUME 10, 2022    it uses a combination of particle swarm algorithm and 565 minimum distance classifier to select the optimal OBFs. 566 PSO-OBFS also uses the chaotic random inertia weight 567 which was regarded as one of the best strategies [51]. Further-568 more, PSO-OBFS linearly reduces the acceleration coeffi-569 cient c 1 from 2.5 to 0.5 and linearly increases the acceleration 570 coefficient c 2 from 0.5 to 2.5 to avoid being trapped in local 571 optima [52]. It performs DRB-RESNET on a single scene. DRB-RESNET 575 uses an improved deep residual network to extract the spa-576 tial and spectral information of the hyperspectral image 577 after dimensionality reduction, and forms a fully connected 578 two-branch feature extraction network. It mainly uses two 579 networks of different depths: the upper network arranges the 580 processed images of the same size into the DRB-ResNet net-581 work of a specified depth (11 layers) for training, and sum 582 the feature extraction data of each region; the lower network 583 first reduces the dimension of the spectral information of each 584 pixel to be classified, and then expands it to 11 × 11 pictures, 585 and then inputs it into the DRB-ResNet network (7 layers) for 586 feature extraction. Finally, combined with the feature extrac-587 tion results of the upper and lower layers, the fully connected 588 network is used to obtain the final classification result [53]. 589 7) IMACO-FS (Feature Select Based on Improved ACO): 590 It selects the optimal feature only based on the labeled 591 samples of a single scene (target scene) to achieve the 592 highest possible overall classification accuracy. A new pre-593 filting method is introduced into IMACO-FS to optimize the 594 pheromone initialization of the ant colony system, thereby 595 speeding up the convergence speed of the ant colony sys-596 tem, and IMACO-FS also adopts pseudo-random rules and 597 adaptive information update strategy to maintain the diver-598 sity of the ant colony. In IMACO-FS, P (user-defined num-599 ber) features with the highest OA value including feature i 600 are selected, and P is set to n/2 (n is the total number of 601 bands) [40]. IMACO-FS initializes the control parameters 602   [45], [46], and the initial values of β and ρ are set to 605 1 and 0.2, respectively. β is an integer between 1 and 10, 606 and the step size is 1 each time it's updated, while ρ is a 607 decimal between 0.1 and 0.9, and the step size of each change 608 is 0.1 [34]. 609

8) Cross-Scene IMACO-FS (IMACO-CSFS): IMACO-610
CSFS is a newly proposed algorithm in this paper, which is 611 executed on the source scene and the target scene respectively. 612 IMACO-CSFS simultaneously considers the consistency of 613 the optimal feature subsets of the two scenarios and the dis-614 criminability of different kinds of ground objects in the target 615 scene. In the experiment, we use SVM with RBF and squared 616 Euclidean distance as the classifier. The maximum number of 617 iterations I M in equation (12) is set to 100; the ant colony 618 size n on each data pair is set to 5; the ant quantity P of 619 the ant colony is set as about 1.5 times of the number of 620 bands (The initial value of P is set to 0.1 times of the number 621 of bands, and it is gradually increased by 0.1 times until it 622 VOLUME 10, 2022   Fig. 7 and Fig. 8. In Fig. 7, (a)   are shown in Table 5. Table 5 shows that the D t in the first    Table 1 for details).  finally makes the features selected in the right part of Fig-805 ure 13 show a shape similar to ''pile''. Therefore, the features 806 selected by IMACO-CSFS around the optimal solution are 807 more evenly distributed than those selected by other methods, 808 IMACO-CSFS tends to select features that are more dispersed 809 and uniform, which benefits from the priority sorting-based 810 ant colony strategy. Fig. 14 contains the classification graphs 811 of several algorithms on PaviaC dataset when F n =40. Com-812 bining Fig. 5 and Fig. 14 Fig. 15 shows the classification accuracy of several algo-817 rithms under different feature dimensions on Hangzhou 818 dataset. As shown in Fig. 15, the OA, AA and κ of 819 IMACO-FS and IMACO-CSFS almost all reached their max-820 imum at F n =8. On the whole, OA and κ of IMACO-CSFS 821 is almost always larger than that of other several algorithms, 822 however, AA of IMACO-CSFS is inferior to that of them 823 when F n is big, its disadvantages is revealed with the increase 824 of F n . Note that when F n ∈ {8, 10, 12}, the κ of CDWOASA 825 is slightly higher than that of CDIRF, which corresponds 826 to the slightly higher κ of CDWOASA than CDIRF in the 827 Hangzhou dataset in Table 7 below. It may be because 828 CDWOASA provides slightly more mutual information about 829 the strong agreement between the ground truth map and the 830 classification map of the Hangzhou dataset than CDIRF. The 831 results of the features selected by the eight algorithms when 832 F n =20 are shown in Figure 16. It can be seen from Fig. 16 that 833 IMACO-CSFS tends to choose more uniform and dispersed 834 features. Fig. 17 shows the classification graphs of the eight 835 algorithms on the Hangzhou dataset when F n =8 and its orig-836 inal ground truth. Combining Fig. 6 and Fig. 17, we can find 837 that IMACO-CSFS has advantages in the land cover classes 838 Water and Land/Building, but the classification effect in the 839 class Plant is not as good as the other seven methods, the 840 reason may be that IMACO-CSFS has a weaker recognition 841 ability for this class, resulting in lower classification perfor-842 mance on this class.

844
Firstly, this paper compares IMACO-CSFS with the origi-845 nal IMACO-FS. Table 6      of IMACO-CSFS on the Pavia data pair is 7.9 seconds 926 longer than that of DRB-RESNET; the computation time 927 of IMACO-CSFS on the Shanghai data pair is 4.7 seconds 928 longer than that of DRB-RESNET. On the three data pairs, the 929 computation time of DRB-RESNET is always the shortest, 930 the second is IMACO-CSFS (the method proposed in this 931 paper), and the last one is the original IMACO-FS method: 932 the computational complexity of all eight methods is shown 933 above the column in Fig. 18. Furthermore, we can find 934 that the execution time of these algorithms on the Pavia 935 data pair is shorter than that on the Huston data pair and 936 the Shanghai-Hangzhou data pair on the whole. It may be 937 because the number of labeled samples of the former is 938 smaller than that of the latter. And the reason why the running 939 time of the algorithms on Huston data pair is greater than that 940 VOLUME 10, 2022   In this paper, a new cross-scene feature selection algorithm 954 IMACO-CSFS is proposed. In order to make the subsequent 955 search process focus on the global optimal solution (optimal 956 feature subset) found in the previous iteration, IMACO-CSFS 957 proposes a priority sorting-based ant colony strategy to obtain 958 more accurate feature subsets of the two scenes than other 959 cross-scene feature selection methods. Furthermore, in order 960 to further accelerate the convergence speed of the global opti-961 mal solution, IMACO-CSFS introduces an ant colony strat-962 egy based on elite ants to obtain the optimal feature subsets of 963 two scenes more efficiently. Finally, the original IMACO-FS 964 method based only on a single scene is successfully applied 965 to the fields of cross-scene HSIs feature selection and image 966 classification in this paper, which further improves the update 967 strategy of pheromone in the ant colony algorithm and greatly 968 improves the feature selection accuracy and image classifica-969 tion accuracy of the target scene. IMACO-CSFS's superiority 970 is demonstrated on three public HSI data pairs-Huston data 971 pair, Pavia data pair and Shanghai data pair. However, this 972 method is not perfect, as the improvement of AA in this 973 research was small, and we might incorporate improving their 974 AA into the scope of the follow-up work, such as selecting 975 some feature subsets M from the candidate feature subset F 976 to pick the best combination of AA, OA and µ. 977 supervised classification of remotely sensed hyperspectral images,'' IEEE