Collection of datasets with DNS over HTTPS traffic

Recently, the Internet has adopted the DNS over HTTPS (DoH) resolution mechanism for privacy-aware network applications. As DoH becomes more disseminated, it has also become a network monitoring research topic. For comprehensive evaluation and comparison of developed classifiers, real-world datasets are needed, motivating this contribution. We created a new large-scale collection of datasets consisting of two classes of traffic: i) DoH HTTPS communication and ii) non-DoH HTTPS connections. The DoH traffic is captured for multiple DoH providers and clients to include nuances of various DoH implementations and configurations. The non-DoH HTTPS connections complement the DoH communication aiming to include a wide range of existing network applications. The dataset collection consists of network traffic generated in a controlled environment and traffic captured from a real ISP network. The resulting datasets thus provide real-world network traffic data suitable for evaluating existing classifiers and the development of new methods.

b s t r a c t
Recently, the Internet has adopted the DNS over HTTPS (DoH) resolution mechanism for privacy-aware network applications. As DoH becomes more disseminated, it has also become a network monitoring research topic. For comprehensive evaluation and comparison of developed classifiers, realworld datasets are needed, motivating this contribution. We created a new large-scale collection of datasets consisting of two classes of traffic: i) DoH HTTPS communication and ii) non-DoH HTTPS connections. The DoH traffic is captured for multiple DoH providers and clients to include nuances of various DoH implementations and configurations. The non-DoH HTTPS connections complement the DoH communication aiming to include a wide range of existing network applications. The dataset collection consists of network traffic generated in a controlled environment and traffic captured from a real ISP network. The resulting datasets thus provide real-world network traffic data suitable for evaluating existing classifiers and the development of new methods.

Value of the Data
• The datasets consist of captures of DNS over HTTPS (DoH) traffic and non-DoH HTTPS traffic that comes both from controlled and real-world environment. The real-world dataset was collected in an Internet Service Provider (ISP) backbone network servicing half-million users. DoH captures consist of traffic generated by several techniques and clients towards 16 selected DoH servers of various implementations and with different configurations. The generated part contains 64,0 0 0 web page accesses and related DoH communication. The datasets aim to provide a comprehensive sample of DoH traffic as observed in real networks and various available implementations and configurations.
• Researchers can use the provided collection as i) a benchmark for the DNS over HTTPS protocol classification and detection algorithms [1][2][3] , ii) for studying differences between the behavior and performance of various DNS resolution methods [4,5] , and because of the presence of both DoH and HTTPS communication also for other iii) encrypted traffic analysis research [6,7] . We redacted the real-world data by anonymization of IP addresses, preserving all other features of communication. Moreover, we also included domain names from TLS headers to allow TLS traffic classification [8,9] . • The datasets collection enables researchers to experiment with DNS over HTTPS traffic recognition and pattern analysis. However, since the data are provided in raw packet captures, it can be suitable for other network traffic analysis tasks, e.g., it can be also used as a realworld benign traffic sample in malware identification challenges [10,11] . • The datasets provide a large and unique combination of labeled traffic. The generated DoH and labeled real-world HTTPS traffic provide the ground truth data for training network classifiers and evaluating their performance. Moreover, packet captures can be used for IP Flowbased traffic analysis, detection, and classification to experiment with novel traffic features. • The datasets collection can also improve the understanding of DoH behavior and other HTTPS traffic phenomena in a large-scale network environment. It is possible to observe and analyze performance and other relevant metrics to improve the specification and implementations. • We are not aware of any other dataset of the comparable size and variety that provides real HTTPS traffic captured on a large-scale ISP in the form of raw packet captures that enable unbounded extraction of arbitrary available features.

Data Description
The datasets consist of labeled HTTPS traffic. The traffic is either the DoH or non-DoH HTTPS communication. The primary focus of this dataset collection is to provide the most variable and representative DoH traffic obtained from communication with many available DoH services. To allow comparison and discrimination of DoH characteristics with other HTTPS traffic types, we also provide general HTTPS traffic captures. Moreover, DoH and general HTTPS traffic comprehensively cover web traffic, allowing dataset utilization for other research tasks not limited to the DoH area.
The collection of datasets contains DoH and HTTPS packets in the total amount of 430 GB raw binary data. The datasets were captured in two enviroments -(i) generated DoH and HTTPS traffic from a controlled experimental system and (ii) real-world DoH and HTTPS traffic from a real large ISP network. The aim of generated data is to provide heterogeneous DoH traffic for various existing DoH implementations. Real-world captured data reflects the real network environment, e.g., packet timing, lost and retransmissions, and connection errors. Full-packet captures are in standardized PCAP format [12] , which is a default for libpcap library and broadly supported by network analysis software, e.g., Wireshark 1 , tcpdump 2 or IDS Suricata. 3 Full-packet captures can be considered the primary source of raw network data that can be used for further analysis or feature extraction. Fig. 1 exemplifies how various network features can be extracted from raw packets for network analysis.
In addition to raw packet captures, we also provide CSV representations of flow data enriched with information from TLS handshake. The flow data were generated for all captures (even the generated ones for dataset consistency). It contains server names (TLS SNI extension), used application protocol (TLS ALPN extension), and TLS JA3 fingerprints. 4 All of these values are commonly used in encrypted traffic analysis [7,13,14] . The description of provided flow data fields is in Table 1 .  Table 1 The description of column headers in CSV files with extended flow data.

Column Name
Column Description

Generated data
The datasets with generated traffic were created in a controlled environment to provide samples of DoH communication for different existing implementations in operating systems and web browsers. The traffic originates in DoH-capable browsers Firefox and Chrome. Browsers access a large collection of web pages and both DoH and HTTPS communication are captured. Table 2 lists of generated pcap files.
The captured files contain a mix of non-DoH HTTPS and DoH traffic. To discriminate the DoH traffic, we provide with each dataset a list of IP addresses of known DoH resolvers used in our experiments. This list can be used to identify and thus label the DoH flows in the PCAP files.
A sample of 16 different DoH servers was chosen to include a diversity of traffic characteristics related to domain name resolution. For each DoH server, we generated traffic by fetching Table 2 Properties of Generated datasets. Abbrevation stands for: Brws -Used Browser, F -Firefox, C -Chrome.  20 0 0 websites on each web browser. Totally, 64,0 0 0 websites requests were generated and corresponding HTTPS communication was captured in the dataset. Firefox web browser supports both GET and POST methods for resolving DoH, therefore 10 0 0 websites were visited by using the DoH GET method, and the other 10 0 0 websites access were performed with the DoH POST method. Note that traffic samples generated by the Firefox browser contain packets up to 64 KB because the browser was executed in a virtual environment with TCP offloading enabled [15] . We decided to include it in the datasets as this is one of the valid real-world scenarios. Datasets with Firefox data have captures located in the /data/generated/pcap/firefox/ directory. Each filename consists of the HTTP method involved in resolution (GET, POST) and the name of the DoH server used for resolution. Table 3 a shows summary metrics of all files created by Firefox across all generated datasets in the collection.
The datasets generated using by Google Chrome are all located at generated/pcap/chrome/ . Each subdirectory name determines used DoH server. Table 3 b depicts summary information about all files created with Chrome browser across all generated datasets in the collection. The difference in the number of connections can be caused by the generating process discussed later and browser-specific implementation. The ratio between used HTTP protocol version in generated datasets is shown in Fig. 2 .
The Chrome browser implements EDNS padding, slightly impacting the DoH flow characteristics when compared to the Firefox DoH flow which yields an increased diversity of DoH traffic.

Real-world data
The real-world data are provided in the DoH-Real-world dataset ( 10.5281/zenodo. 5956043 ) and contains DoH HTTPS and web-based HTTPS communication. The real-world traffic captures were anonymized to protect the privacy of real users (packets have no payload, clients IP, and MAC addresses were hashed). However, the anonymization does not pose a major limita-  Contrary to the generated traffic, the DoH communication was captured in the real network, which reflects ratios of most utilized DoH Internet services as seen in Fig. 3 . The most DoH   connections involve Google DoH resolvers, and the five most popular resolvers represent more than 93% of all DoH traffic. The importance of the real-world dataset is that it contains the HTTPS traffic of thousands of genuine users. The Table 6 shows overall statistics of captured dataset. It should be noticed that there is 116,263 unique client addresses in the network. Since the capture was created on an Internet Service Provider (ISP) backbone, each IP address may correspond to multiple physical/virtual devices due to Network Address Translation (NAT). The actual number of clients can thus be larger.

Supplementary files
In addition to the main dataset content, we also provide Supplementary files ( 10.5281/ zenodo.6024913 ) consisting of software scripts for data generation and post-processing. The list and short description of supplementary files, are provided in Table 7 . Table 7 The list and description of supplementary files.

doh-resolvers
The list of IP addresses of known DoH resolver, which was used during the real-world part creation. The list is in the CSV file format.

pcap-anonymizer
Python scripts used for the anonymization of real captures.

firefox|chrome/scripts/
Folder containing scripts and other files used for data generation.
firefox|chrome/domains/ < num > .csv Files containing the domains used for data generation.

doh_server_urls.txt
File containing DoH server urls used in generation process.

Experimental Design, Materials and Methods
In this section, we describe the data acquisition environments. Firstly, we describe the environment used for DoH traffic generation ( Section 2.1 ), then we provide an overview of Real-World Data acquisition procedure and anonymization ( Section 2.2 ).
The datasets contain all DoH and non-DoH captured HTTPS flows. The common attribute of the whole dataset collection is that the data are captured in a pure form, where the process and all necessary information are covered in this paper. Hence, the dataset may contain all sorts of DoH and non-DoH HTTPS flows, e.g., long, short, successful, and failed in the way they appeared during the capturing process. We have not performed any additional filtering that would lead to information loss except those mentioned in this section.

Generated data
Two DoH-enabled web browsers were used for making web requests. We selected Chrome and Firefox because they have the most advanced DoH implementations. Since we primarily target DoH and Non-DoH traffic generation, we have to cover all possible settings that can affect the DoH traffic characteristics (such as packet sizes) generated by those applications.
Chrome browser is based on Chromium code-base that shares DoH resolution implementation with other Chromium-based browsers such as Edge 5 , Brave. 6 Firefox browser provides possibility of choosing DoH HTTP method, either GET or POST. It provides more management options than Chrome, such as whether to force DoH usage or not and which server to use (Cloudflare's DoH service by default).
The generated data samples were created by visiting a collection of websites that are part of the Majestic million list [17] . The URLs were taken from the beginning of the list and each website was visited by one of the browsers. The list of visited websites is included in the supplementary files, see Section 1.3 .
DoH traffic characteristics are mostly because of client-server interaction. DoH servers can employ different implementations and be deployed with different configurations. Moreover, the HTTPS ecosystem on the server-side can consist of proxies, firewalls, load-balancers, and other systems that impact network communication. To cover various DoH traffic patterns, we included 16 different DoH servers in our datasets. The 16 DoH servers were chosen to cover various DoH traffic patterns as described in [18] , where we can find that the major difference is in the HTTP version and the number of HTTP headers used. The servers were taken from the same DoH monitoring tool as mentioned in the article [16] , supporting different features, and each server was deployed by a different provider. The CSV files with flow data enriched for TLS information were then created with ipfixprobe 7 open-source flow exporter. We used its TLS plugin to create flow data extended for TLS information.
The whole traffic generation process took place at the Brno University of Technology, Faculty of Information Technology. The machines were directly connected to the university network.

Firefox browser traffic generation
At first, the Firefox browser in version 76.0.1 was used in the Docker container. The use of the Docker container does not influence the traffic with specific properties; however, it simplified the generation process automation and provided an isolated application environment. The platform utilized in the data generation process using Firefox browser was Supermicro SuperTwin2 6026TT-TF server equipped with eight Intel (R) Xeon E5520 @ 2.26 GHz. The cluster consists of 4 nodes. The nodes were equipped with 48 GB RAM and 16 CPU cores. The platform was directly connected to the Brno University of Technology network. Network traffic was captured directly in the container using the tcpdump command-line tool.
Python and Selenium libraries were used to automate generation process. To be able to provide more realistic behavior of the browser, we used X-virtual frame buffer, which implements the X11 display server protocol where all graphical operations are made in virtual memory without graphical output [19] . The browser than works with a standard graphical output even within the container. In addition, we restricted the container to use 4 CPUs and 4 GB of RAM at maximum.
Firefox had disabled the local web cache during web page access to get more network traffic samples. The process consists of opening the browser, fetching the website, timeout, and waiting for the loading of the website, followed by browser closing. The schema of this generation is depicted in Fig. 4 . This process was repeated 10 0 0 times for DoH POST and GET methods and for each of the 16 DoH servers. Together, 32,0 0 0 different websites were fetched.

Chrome browser traffic generation
Since the used version of Chrome browser (94.0.4606.81) does not support DoH on Linux OS at the time of our data generation, we could not reuse the same generation environment as in the case of Firefox. The Chrome browser generation took place on a separate machine running Windows OS without Docker or any virtual environment, running directly on the bare metal. The machine was equipped with 3rd generation Intel Core i5, 8 GB RAM also connected directly to the same university network as in the case of Firefox.  The traffic generation process itself was similar to Firefox. It consisted of opening the browser with the given website, timeout, and waiting for loading the website, followed by browser closing. This generation is shown in supportive schema in Fig. 5 . The automation was handled by a python script with a manual setting of custom DoH servers. The same 16 DoH servers as in the case of Firefox were used. Chrome does not allow users to set HTTP GET or POST methods for querying DoH. Hence, the amount of fetched websites is similar to the Firefox case: 2 ×10 0 0 only with the default Chrome method used. In total, 32,0 0 0 different websites were fetched.

Real-World data
The Real-World data of the dataset were captured on the monitoring points of the CESNET organization, which is the Czech National Research and Education Network operator. CESNET operates a backbone network infrastructure called CESNET2, which is used by half-million users. CESNET2 provides internet connectivity to universities, campuses, research centers, schools, hospitals, and some government offices. The topology of the CESNET2 backbone network is depicted at Fig. 6 . The six monitoring points are located at the perimeter of the CESNET2 network in three locations in the Czech Republic -Prague, Brno, and Ostrava. Each of them is connected to one or multiple 100 Gbps lines via passive optical TAP. Since the packet capture is distributed across multiple points and locations, we need to deal with time shifts on each monitoring server resulting in inaccuracies in packet timestamps.
The points are synchronized via Network Time Protocol (NTP), with an accuracy of up to 10 ms to minimize the time inaccuracies. The time shift between monitoring probes does not affect connections that are routed symmetrically, meaning that all packets belonging to the one TCP connection pass through the same monitoring point.
As a result, the Real-World Data contains inaccurate interarrival time values between packets in the opposite directions -packets in the same direction are always routed via the same monitoring points. According to long-term measurement in the CESNET2 network, 40 % of all 443/TCP connections are routed asymmetrically, thus affected by the time shifts. However, the timestamp inaccuracies are small, and they do not cause any significant disturbances in packet order within the connection (HTTPS requests are always before HTTPS responses). Moreover, since the distributed monitoring infrastructure is common in large-scale network monitoring [20] , these data imperfections are viable for the dataset creating the detection algorithms prepared for such deployment.
As it can be seen in Fig. 7 , the creation of Real-World data can be divided into four steps: (i) Capturing filter preparation, (ii) Traffic Capture, (iii) Deduplication, (iv) Anonymization. Each step is described in the following sections.

Capturing filter preparation
The monitoring points can process 100 Gbps traffic; however, it is not possible to capture all ongoing traffic due to bandwidth limitation and writing speeds of the storage, which would result in missing packets in the dataset. Therefore, it was necessary to create a capturing filter with this limitation in mind. Since the use of DoH is still small, we could capture all DoH traffic by filtering based on the DoH resolver's list (which is included with the dataset as supplementary file -described in Section 1.3 ).
In the case of HTTPS traffic, we have limited the capturing to the selected /24, and /22 IPv4 address ranges assigned to the university campus. Additionally, we have monitored the whole capturing process for packet drops to ensure that all transmitted packets in filtered connections were saved on the disks and the points were not overloaded; thus, no artificial packet drops were added.
Overall, we created two capturing filters, which remained unchanged for all traffic captures during the dataset's creation. According to the monitoring of known DoH resolvers [16] , there was not any change during the capturing period.

Traffic capture
Since the CESNET2 network is a large infrastructure with multiple monitoring points, we have created custom scripts for concurrent full-packet capture on each point using parallel SSH 8 . Fig. 8 shows a scheme of the overall capturing process. As the first step, the capturing controller distributes the command to the monitoring points. The command contains a filter and maximal duration of capturing. When the packet capture exceeds the maximum duration, the pcap file is  distributed back to the controller. The controller then merges partial packet captures from the monitoring probes into a single one, using mergecap utility. 9

Packet deduplication
Since a backbone ISP network can also carry transit traffic, there is a chance of capturing duplicate packets on distributed monitoring infrastructure. The deduplication was performed with editcap 10 utility with a time window of 10 ms. Thus, if there are two, or more duplicate packets within 10 ms time window, only the first packet will be maintained.

Anonymization and flow export
Since we are using data from real ISP backbone lines, we needed to remove sensitive information. The anonymization of captured traffic was performed automatically; thus, nobody could view or analyze raw source data. The packet fields, that were anonymized are shown within the packet context in Fig. 9 .
At first, we anonymized all addresses (IPv4, IPv6, and MAC), except the IP addresses of DoH resolvers, which remained unmodified. The anonymization process replaced each address by a part of its SHA256 hash (the first N bytes, where N is the length of the corresponding address). The hashing algorithm ensures that a particular address is always mapped onto the same anonymized value.
PCAP files with anonymized IP addresses were then used to create TLS enriched flows with ipfixprobe flow exporter. The flows were created the same way as in the generated part.
After the flow export, we also anonymized the payload from all packets by substituting every byte with the letter 'X'. Removing packets' payload has only a negligible impact on the value of the data (since they are mostly encrypted anyway); however, it ensures the future privacy of the real users.

Ethics Statement
Privacy of users is our essential priority, so our whole research was done with extreme carefulness. The indisputable advantages of real traffic generated by hundreds of thousands of people come with the cost of potential privacy abuse of real users. Therefore, we used only automatic data processing with immediate data anonymization. With this, we declare that we did not analyze or manually process deanonymized data, and we did not perform any procedures that could lead us to the user's identity. Moreover, there is no possibility of revealing real users' identities from the provided dataset.
Brno University of Technology, Faculty of Information Technology internal grant FIT-S-20-6293, and also by Technology Agency of the Czech Republic, grant No. FW03010099: Context-based Encrypted Traffic Analysis Using Flow Data.


Rotation-inspired circuit cut optimization

Recent works have demonstrated that large quantum circuits can be cut and decomposed into smaller clusters of quantum circuits with fewer qubits that can be executed independently on a small quantum computer. Classical post-processing then combines the results from each cluster to reconstruct the output of the original quantum circuit. However, the runtime for such hybrid quantum-classical algorithms is exponential in the number of cuts on a circuit. We propose Rotation-Inspired Circuit Cut Optimization (RICCO), an alternative method which reduces the post-processing overhead of circuit cutting, at the cost of having to solve an optimization problem. RICCO introduces unitary rotations at cut locations to rotate the quantum state such that expectation values with respect to one set of observables are maximized and others are set to zero. We demonstrate practical application of RICCO to VQE by classically simulating a small instance of VQE and comparing it to one of the existing circuit-cutting methods.

I. INTRODUCTION
Simulating quantum circuits with a large number of qubits is intractable on a classical computer. Also, it is difficult to program these circuits on actual quantum hardware due to size constraints (insufficient qubits) and because they are also error-prone. To increase the use of current noisy small-scale quantum devices, methods have been developed to combine classical and quantum computers to simulate quantum circuits with a large number of qubits [1], [3]- [5].
For example, Peng et al. [1] demonstrated that a quantum circuit can be divided into subcircuits, each with fewer qubits, which can then be run separately on a quantum computer too small to run the original circuit. The quantum computer's results are sent to a classical computer, which combines them to replicate the expected output of the original quantum circuit. The number of quantum measurements and traditional post-processing involved in replicating the original quantum circuit increases with the number of cuts. Specifically, the cost of executing the original quantum circuit is exponential in the number of cuts. Tang et al. [4] proposed CutQC, a framework that employs heuristics to select a cut reducing classical post-processing overhead. To reduce classical postprocessing resources necessary to characterize circuit cutting, Perlin et al. [3] proposed maximum-likelihood fragment tomography (MLFT). More recently, Lowe et al. [5] employed a randomized measurements method to speedup quantum circuit cutting. We refer to these procedures generically as quantum circuit cutting.
In this paper we propose Rotation-Inspired Circuit Cut Optimization (RICCO), a method to reduce the cost of simulating a large quantum circuit on a small quantum computer. RICCO introduces a parameterized unitary operator and its adjoint at each cut location. After optimization, the unitary operator rotates a given quantum state such that it maximizes the expectation value of select Pauli observables. This is similar to rotating a quantum state vector such that it aligns with one of the Bloch sphere axes for a 1-qubit system. By this method, we effectively set the expectation value of the select Pauli observable(s) to maximum while the others are set to zero. This is relevant because we have eliminated the need to measure in multiple bases, as is the case in the original method, resulting in fewer circuit executions. We demonstrate through simulation that when compared to Peng et al. [1] RICCO reduces the total number of quantum circuit measurements required to reconstruct a cut circuit, at the expense of the optimization.

II. THEORY
Consider the quantum circuit in Fig. 1. A cut is applied on the qubits labeled B resulting in two separate subcircuits: upstream and downstream subcircuits. The ultimate goal is to combine the two subcircuits to reconstruct the output of the original uncut circuit. Our focus at the moment will be on the upstream subcircuit with the qubit register labeled A and the parts of qubit register B just before the cut. The downstream (a) Expectation value error of the observable Z ⊗n = Z ⊗(n A +n B +n C ) for random circuits. The error bar indicates a deviation of less than 0.3% for two qubit cuts (b) Circuit executions during optimization. The error bar is high because of the variance in the class of quantum circuits optimized using RICCO (c) The number of executions of quantum circuits for RICCO is significantly less than that for QCUT after optimization. subcircuit consists of the part of qubit register B just after the cut and beyond and the qubit register labeled C.
In prior circuit cutting methods, measurements are made in multiple bases at the cut point, and specific eigenstates prepared to subsequently run the downstream circuits. This results in an exponentially large number of measurements and circuits. The goal of our proposed method is to find a way to reduce the number of measurements.
Without loss of generality, we are interested in computing the expectation value of some observable, say, Z ⊗n = Z ⊗(n A +n B +n C ) (this also applies to observables with some Pauli X and Y , in which case we just apply the corresponding local basis rotations) for the circuit in Fig. 1, where n A , n B and n C are the number of qubits for register labels A, B and C respectively. To do this, we seek a unitary operator U that rotates the state ρ AB to the state ρ AB such that ρ AB aligns with one of the computational basis vectors of the observables in {{Z ⊗n A }⊗{I, Z} n B }. As a result, we only have to measure in the computational basis. To find U , we decompose U into parameterized rotations and optimize for its parameters in the quantum circuit.
To understand the intuition behind our method, consider that the density matrix ρ of an arbitrary 1-qubit quantum state can be expressed in terms of Pauli matrices, Here P i ∈ {I, X, Y, Z}, and Tr(P i ρ) is the expectation value, P i of the Pauli P i . By applying the unitary rotation U to ρ we transform ρ to the state ρ as follows: Since ρ is also a density matrix, we can express it as The first intuition behind our method is the Bloch sphere. The expectation values X , Y , and Z of the three Pauli matrices, allow one to identify the three coordinates with respect to X, Y , and Z axes. If we can find a unitary operator that rotates the Bloch vector of the density matrix to align with a particular state (e.g., the state with positive eigenvalue) of the Z-axis, we are by implication zeroing out the expectation values of Pauli X and Y observables. To avoid restricting ourselves to one eigenstate of Pauli Z, we can find a unitary, U such that when applied to ρ, its final state ρ ends up along either the |0 axis ( Z = +1) or |1 axis ( Z = −1).
The second intuition behind our method stems from the theory of optimal measurements and mutually unbiased bases (MUBs). MUBs are defined such that for any two vectors |µ and |ν from bases M and N , where d is the dimension of the system. MUBs for an n-qubit system have a standard construction based on the partitioning of the Pauli group into 2 n +1 disjoint sets of 2 n −1 commuting operators whose mutual eigenvectors form the measurement bases [6]. If we align our state with an eigenvector of one such basis, it will have equal overlap with all vectors in the other sets. This has further implications for the expectation values of Paulis: the state will have a non-zero expectation value for all Paulis in that set we measure with respect to, while for Paulis from other sets, they will be 0. This enables us to perform only a single measurement in the Z basis.

A. Cost Function to Determine the Unitary Operator
Let's consider the case where ρ is a single-qubit state. We want to rotate the state such that we point along either the |0 or |1 axis of the Bloch sphere; these are eigenstates of only Pauli Z, so we can ignore any observable, P i / ∈ {I, Z}. The identity operator is kept because |0 and |1 are also the eigenstates of I.
To design the right cost function, that optimizes the parameters of a unitary operator U such that it rotates ρ to a new state ρ , we can use a measure of distance such as fidelity [11]: to determine how close ρ is to our desired state |ψ . We decided to use fidelity in this case because it has a straightforward definition for the case of pure versus mixed state. We want to see how close our transformed mixed state ρ from the cut is to the desired pure state (a computational basis state). The best possible outcome is F = 1, when ρ is the same state as the desired state |ψ ψ|. The worst is F = 0 when the desired state is perfectly orthogonal to ρ . Another nice property of fidelity is that it can be expressed in terms of expectation values as [10], [12] These expectation values are what we measure in RICCO (see Algorithm 1). To determine U , we have to build a cost function that maximizes the fidelity with respect to either the |0 or |1 state. We propose to use the cost function This cost function is simply the squared difference between 1 and the product of the fidelities of |0 and |1 with ρ . We squared the cost function to ensure that it is a convex function for optimization purposes. Thus, for the one qubit case, our cost function is minimized when Z ρ = ±1. It is maximized when Z ρ = 0, which corresponds to states along the other two axes X and Y . Thus, by minimizing this cost function, we are effectively setting the expectation values of X and Y to zero.
For the general case, to zero out as many expectation values as possible, we need to rotate only the qubits being cut so that ρ aligns with any of the computational basis vectors. Then we can measure the upper subcircuit in only the computational basis. It is important to note that as a result of optimization error, the other uncut qubits of ρ may not align perfectly with the computational basis vectors of their observables. Our results show that this error is negligible in most cases and does not constitute a major limitation in reconstructing the original expectation value. For n qubits, we can express the density matrix of an arbitrary state ρ as in (8), The formula for the fidelity is the same, except now we must consider all observables that consist only of Zs and Is for the cut qubit. The uncut qubits still retain their respective Pauli observables. For n B cut qubits, the observables are in the set {I, Z} n B . The fidelity of the transformed state ρ and one of the computational basis states |ψ of the combined observables P i ∈ P (where P = {{Z ⊗n A } ⊗ {I, Z} n B }) of the cut qubits, n B and the uncut qubits, n A is Algorithm 1: RICCO algorithm Input: upstream and downstream subcircuit, observables in Hamiltonian and cost function equation 7 Output: Y , Expectation value ← 10 −7 set the tolerance of RICCO; t ← 1 set the initial tolerance; Initialize X (parameters of U ); Initialize cost previous ; for observable in optimization observables do Initialize params; while t > do params, cost new ← Optimize the cost function in equation 7; t ← |cost previous − cost new | ; end X ← params; end E ← measurement of expectation value of subcircuits with optimized X; Y ← combining results of subcircuits' expectation values using tensor product of corresponding observables; return Y Finally, we use these to construct our cost function. Let |x , forx ∈ {0, 1} n denote our n-qubit computational basis state. The function to minimize is Whenever we rotate the state ρ and align it with one of the computational basis states, the fidelities with respect to all others will be minimized. Furthermore, based on our choice of basis state, we can reconstruct the original expectation value for the observable of interest with only one measurement setting in the computational basis.

III. EXPERIMENTAL RESULTS
We evaluated RICCO by comparing the total number of quantum circuit executions after optimization against the method by Peng et al. [1], which we denote as QCUT. Because of the available computing resources, we simulated a total of 220 quantum circuits. The circuits were composed of two randomly-selected operations from unitary groups U (16) and U (32) for 1-qubit cut and 2-qubit cuts respectively. This is to ensure that we tested our method on the worst case possible. The total number of qubits for each circuit was 7 and 8 for 1-qubit cut and 2-qubit cuts respectively. Simulations and optimization (using Adam [7]) were performed using PennyLane [2], an open-source Python software framework for quantum differentiable programming. Fig. 2a plots the percentage error of the expectation values Z ⊗n = Z ⊗(n A +n B +n C ) of RICCO and QCUT circuit cutting methods compared with the ideal classical simulation of the uncut circuit for different number of cut qubits. This plot gives an insight into the accuracy of our method. We observed a reasonably high accuracy in RICCO with the worst case having less than 1.5% error. Fig. 2b plots the number of quantum circuit executions during optimization for the parameters of U in RICCO for 1-qubit and 2-qubits cut. On average, we require about 400 quantum circuit executions during optimization for the 15 parameters of U in the case of a 2-qubit cut. This number is very high partly because of the type of optimizer used. A better optimizer could reduce the number of executions required. We believe that more work needs to be done to explore how the optimization can be improved. quantum circuit executions for RICCO was 3 while QCUT was 7; and for 2 cuts, the total number of quantum circuit executions for RICCO was 5 while QCUT was 25. This result shows that, after optimization, RICCO improves runtime compared to QCUT.

IV. CIRCUIT CUTTING WITH VQE
After verifying correctness of the procedure and its behaviour in the worst case, it is of interest to evaluate its poten- indicates the cut location. Two subcircuits emerge after the cut: the upstream and downstream subcircuits. The upstream subcircuit consists of qubit A0, the part of qubit B0 just before the cut and A1. The downstream subcircuit consists of qubit C0 and the part of qubit B0 just after the cut location tial in practical applications. Variational quantum eigensolver (VQE) is one of the quantum algorithms in computational quantum chemistry used to predict the electronic structure and properties of molecules. In this section, we discuss the results of applying RICCO to perform numerical simulations of VQE to estimate the ground state energy of the hydrogen molecule.
PennyLane was used to build a representation of the electronic Hamiltonian. We used the circuit from the PennyLane demo by Li et al. [9] as our ansatz (see Fig. 5). The parametrized circuit prepares the trial quantum state of H 2 molecule which is trained using a cost function that measures the expectation value of the problem Hamiltonian in the trial state to compute the ground state energy. The paramters of VQE were randomly initialized at the start of optimization with Adam's optimizer.
To apply RICCO to VQE, one of the qubits of the original ansatz was cut which resulted into two subcircuits before training. During the training phase of the VQE for the cut circuit, RICCO was used to reconstruct the expectation value for the cost function (see Algorithm 2).
The Hamiltonian consists of 18 observables which contribute to its expectation value. After cutting the circuit, it is expected that each of these observables are treated as separate entities and optimized separately with RICCO (see Algorithm 1). The results of the expectation value of each observable are then recombined to reconstruct the original expectation value of the uncut circuit. To reduce the number of RICCO optimiza- At every step of VQE optimization, RICCO was deployed to reconstruct the original expectation value of the Hamiltonian. The set of 6 unique observables also consists of non-Z observables, {XIY, XZY } and {Y IX, Y ZX} that were also optimized separately with RICCO at each step of the VQE optimization. It is important to do this to ensure the accurate reconstruction of the original expectation value. This process is repeated until VQE converges to a preset tolerance of about 10 −6 . It took more than 70 iterations for RICCO to converge to this tolerance as shown in Fig. 3a. Fig. 3a compares VQE cost function versus the number of iterations for RICCO, QCUT and the conventional VQE optimization methods. Adam optimizer was used for the different methods. Results from the optimization showed that although RICCO coverged, it did not converge to the actual value with the other two methods. This is evident in Fig. 3b which shows that the difference in RICCO's energy compared with the true value was worse (about 10 −2 Hartree) than QCUT and conventional VQE with no cutting (about 10 −6 Hartree). RICCO's convergence threshold was set to 10 −7 . Increasing this threshold indicated some improvements in the energy difference at the expense of increased circuit executions and runtime. We noted that large difference in RICCO's energy could be as a result of the errors incurred during RICCO's optimization and the convergence threshold set for RICCO.
Notwithstanding, exploring other causes of this error and how to mitigate it would be the focus of our future work.

V. DISCUSSION
We want to understand how to reduce the number of measurements in cutting procedures, and see if this translates to improvements in actual applications such as VQE. Fig. 6a shows the number of executions required to implement VQE with RICCO and QCUT. RICCO requires less than half the number of quantum circuit executions compared to QCUT. Fig. 6b on the other hand, shows that number of executions required for RICCO optimization was at least 8X more than the number of executions for QCUT. This shows the tradeoff between the improved runtime for RICCO and number of executions required for its optimization. We observed that this increase in the number of executions was because for each iteration of VQE, RICCO was optimized. One way to solve this problem could be to combine RICCO optimization with VQE optimization for each iteration of VQE instead of running a separate optimization for each. We can also reseed RICCO parameters at each optimization step.
Finally, we observed that the order of optimization of the 6 unique sets of observables could have contributed to the increased number of executions for RICCO. This is because, the optimization of any set of the observables depends on the outcome of the optimized parameters for the previous set of observable that was optimized by RICCO. Alternative distance measures such as trace distance can be explored as an area for future work to determine how they affect optimization convergence


extfieldsofstudy:
Computer Science
Mathematics

text:
Improving IoT-over-Satellite Connectivity using Frame Repetition Technique

Through Non-Terrestrial Networks (NTN), a global coverage connecting areas with minimal or no terrestrial services is envisaged using satellite and air-borne platforms. However, one of the challenges to delivering NTN communications is the increased levels of path-loss, relative to typical terrestrial scenarios, due to the vast communication distance towards satellite platforms. Typically, cellular IoT technologies adopt frame repetitions where many redundant data transmissions are included to enhance its success rate. This can be also extended to satellite links to counter the effect of long distance communications on the link budget. In this paper, we put forth an analytic framework that captures the repetition behavior based on the probability of line-of-sight (LoS), as it heavily influences the propagation conditions. We analyze the coverage performance in terms of the frame success rate, which allows us to tune the repetitions behavior to suit a given satellite admittance region. Furthermore, we obtain an optimal global average success rate by considering the availability of satellites. The presented framework can facilitate coverage enhancements for IoT-over-Satellite networks.

I. INTRODUCTION
The current era is witnessing a massive surge in IoT adoptions arising from a wide range of applications such as, smart farming, logistics, environmental monitoring, and smart cities. The number of cellular IoT subscriptions are expected to reach a staggering number of approximately 5 billion by 2026 [1]. To support the increasingly ubiquitous IoT applications, it is required to have a network offering an extended worldwide coverage [2]. So far, the traditional terrestrial network has been limited to providing robust connectivity to urban, suburban and some rural areas. However, due to the lack of economical, and sometimes technical, feasibility such terrestrial networks fell short in providing such an envisioned global coverage. As such, Third Generation Partnership Project (3GPP) included Non-Terrestrial Networks (NTN) [3] as part of 5G road-map towards next generation hybrid networks [4], [5]. The NTN comprises space and aerial platforms, i.e., satellites and UAVs, that act either as relaying nodes or as stand-alone base-stations, i.e., gNB, and hence paving the road for global coverage services. However, the application of IoT over NTN comes with many challenges. The NTN channel has much higher path-loss when compared to terrestrial, where most of the propagation path is under free space condition and only interacts with This research is supported by the Australian Government Research Training Program (RTP) Scholarship. The authors are with the School of Engineering, RMIT University, Melbourne, VIC 3000, Australia a clutter layer, terrestrial scatterers, just before reaching the ground user. This is especially applicable for lower frequencies suitable for some IoT applications in the L and S bands. This is also quite different than terrestrial channels where all the propagation happens within the clutter layer. As such, LoS probability plays a major role in determining the radio channel conditions. In turn, the LoS probability is highly dependent on the elevation angle [6], which is the angle between the satellite and the horizon as seen by the ground user. The LoS dependency in an NTN channel impacts the signal power and coverage probability. One of the attributes adopted by cellular IoT access technologies, specifically by Narrow-band Internet-of-Things (NB-IoT), is to improve the coverage by means of frame repetition. In such mechanisms, devices are configured to repeat a message multiple times to improve the detection probability at the receiver. As a result, devices that operate under poor coverage conditions tend to transmit more frame repetitions than those under better conditions. Multiple works evaluate the performance of repetition mechanism on terrestrial network [7], [8]. Moreover, the evaluation of IoT over satellite access technologies has been proposed in several good works such as [9]- [11], where work in [9] analyzes the performance of NB-IoT and addresses the challenges related to delay, while the work in [10] presents the link budget analysis for IoT devices. However, many challenges remain unresolved such as the impact of interference induced by large numbers of repetition and the performance of frame repetition in satellite radio channels. In this paper, we leverage the concept of frame repetition, which is already deployed in several terrestrial IoT access technologies, and extend its application in IoT-over-satellite under realistic satellite-to-ground radio propagation conditions. We further evaluate the coverage probability of IoT devices in the admittance area of a satellite, while implementing repetition scheme. The contribution of this paper is the use of tools from stochastic geometry for analytically formulating a framework that models frame repetitions in terms of probability of lineof-sight. The paper further demonstrates the ability to tune the repetition factor and the admittance region such that an optimal repetition profile maximizes the global average success rate while considering the satellite availability.

II. SYSTEM MODEL A. Geometric Model
Consider a satellite in a circular orbit at an altitude h from the mean sea level. Ground IoT devices are assumed to be distributed by an isotropic and homogeneous Poisson Point Process (PPP) on the Earth's surface. An important parameter that impacts the performance is the angle formed between the satellite point, the Earth center, and the observer, this angle is termed as zenith angle, denoted ϕ.The elevation angle θ, which is the angle between satellite and horizon seen by ground user, can be obtained from the zenith angle using geometric reasoning as follows [4] The domain of the elevation angle is θ ∈ [0, π 2 ], where θ = 0 represents the case when the satellite is at the horizon and θ = π/2 is when the satellite is at the device zenith. Accordingly, , where ϕ h is the maximum zenith angle formed at the horizon. The scenario is depicted in Fig. 1. It can be shown that the maximal horizon zenith angle is given by [12], ϕ h = acos α. As such, the serving satellite can be geometrically seen from any point in a certain spherical cap region. However, note that IoT devices that have very low elevation angle are deemed to have an extremely poor radio link performance as a result of deep path-loss. Consequently, the probability of successful transmission in such cases is low. As such, we define a minimum elevation angle, denoted as θ min , which limits the admittance region of IoT devices. There are multiple ways to limit the admittance region practically; for example by controlling the beamwidth of the satellite or by programmaticaly preventing IoT devices from transmitting when the satellite is below the threshold elevation angle, i.e., θ ≤ θ min . Due to this impact, the zenith angle becomes lower than ϕ h and thus limited by the general threshold ϕ max , i.e. ϕ ∈ [0, ϕ max ] → θ ∈ [ π 2 , θ min ]. The relation between these two angles is given by (1) as follows Due to the near ground objects and foliage, i.e., also called clutter, geometric LoS is sometimes blocked by these obstacles. The probability of line-of-sight p LoS is an important parameter that heavily influences the path-loss and the probability of successful transmission. Therefore, we adopt the LoS model developed in [12]. However, we re-write the p LoS in terms of ϕ as follows where β is a controllable modeling parameter related to the underlaying propagation environment. Note that the probability of non-line-sight is given as p NLoS = 1 − p LoS .

B. Repetition Model
Since the purpose of repetitions to improve the coverage [13], the devices with lower line-of-sight probability are expected to rely on this feature more often in order to achieve successful transmission. This intuitively indicates that devices with larger zenith angle, ϕ, will need more repetitions. Hence, the developed model in [14] which primarily depends on the distance towards the terrestrial base station is adopted but with changing this dependency to the zenith angle. Assuming an initial duty cycle of D o 1 , the effective duty cycle, when, repetition is applied, is modeled as where g(ϕ) is a custom shaping function that satisfies the , since the duty cycle cannot exceed unity. A suitable function that satisfies the above requirement is a modified version of the probability of the nonline-of-sight as follows, where ϕ ∈ (0, ϕ max ) and a ∈ [0, 1] is a controllable tuning factor. Accordingly, the number of transmissions, i.e., repetitions is give by The impact of tuning factor a on the duty cycle and repetitions is depicted in Fig. 2 and Fig. 3, respectively. Smaller values of a correspond to lower overall repetitions and vice-versa. Furthermore, the figures also illustrate the effect of the admittance angle θ min . It is to be noted that both a and θ min can be tuned in order to optimize the overall success probability. Due to increased transmissions as a result of repetition, the effective device density in the admittance region becomes A thorough explanation of repetition model is given in [14].

C. Zenith angle distribution
For a given satellite admittance region defined by the angle domain, θ ∈ θ min , π 2 , the users will view the serving satellite within the zenith angle domain, ϕ ∈ (0, ϕ max ), where ϕ max ≤ ϕ h , as shown in Fig. 1. As introduced earlier, the elevation angle or the equivalent zenith angle of the satellite heavily impacts the repetition rate. Thus, the analysis of satellite's zenith angle distribution is essential to obtain the probability of successful transmission. Assuming PPP, the ground devices or IoT devices are randomly placed within the admittance region. Accordingly, the probability of finding devices within a given angle ϕ o is formulated as follows Fϕ(ϕo)=P(ϕ≤ϕo) = avg. no. of active devices in ϕ o avg. no. of active devices in ϕ max , which represents the cumulative distribution function (CDF) of the satellite zenith angle. The average number of points can be found by invoking Campbell theorem [15] for a PPP. Accordingly, the average number of points in a region bounded by an angle [0, ϕ] is given by because the spherical strip area is given by, 2πR 2 sin(ϕ)dϕ. Substituting in (8) and replacing λ(ϕ) = λ o D(ϕ) we obtain The probability density function (PDF) is thus obtained by differentiating (9) as follows An illustration of the CDF at θ min = 10 o is shown in Fig. 4 along with Monte-carlo simulations to verify the analytical results. Note that the distribution changes for different tuning factor a.

III. LINK PERFORMANCE
Majority of IoT technologies rely on simple random access techniques in order to reduce the complexity of the IoT devices. Examples can be found in LoRaWAN and in the initial random access stage of NB-IoT. Accordingly, collision of frames is bound to happen due to the lack of coordination.

A. Channel Model
For a given device with transmitting power p t , at an angle ϕ o , the received signal power in the uplink can be modeled is the excess gain caused by ground clutter, and its distribution follows simple Gaussian mixture, this path-gain heavily depends on the probability of LoS p LoS as given in eq. (7) in [17].

B. Interference
For a given device located at x o , the simultaneous transmissions imposed by other devices within the admittance region cause interference with a total aggregated power, denoted I, obtained as follows where Φ is the set of all devices in the admittance region, κ is a system parameter that captures the impact of interference mitigation caused by resource coordination in the access system. Note that κ = 1 indicates the worst case scenario when no coordination is implemented. Accordingly, the spatial average of the interference is obtained by invoking Campbell's theorem of sums and is given as whereζ is the average excess path-gain expressed as [6] ζ = p LoS exp −µ LoS + σ 2 where µ and σ are the mean and standard deviation parameters associated with respective LoS and non-LoS components. To validate the above analytic results, we depict in Fig. 5 Fig. 5. The average interference power given in (12) versus the tuning factor a with different admittance angles θ min .
comparison between the analytic results of (12) and Monte Carlo simulations at various values of a and θ min . Note the increase in the aggregated interference power when deploying the repetition scheme. Nevertheless, it is shown in the next subsection that the probability of success can be further increased with careful tuning of a and θ.

C. Probability of success
To measure the probability of success under influence of frame repetition, we utilize a vital metric of communication link which is the ratio of received signal power to the interference plus noise power (SINR). The probability of a frame having an SINR greater than a given design threshold γ, is equivalent to the frame success probability. For a given IoT device viewing the satellite at a zenith angle ϕ o , the success probability of a single transmission is formulated as follows where W is the average thermal noise power,Ī is given in (12) and F ζ is given as [4] F ζ =p LoS FX x,−ρµ LoS , where F X is the CDF of the log-normal distribution as and ρ = ln 10 10 . For a device attempting its transmission N number of times, the probability of at least one of the frames to be successfully received is formulated as We illustrate the success probability in Fig. 6 with respect to a range of elevation angles for three given examples of the admittance region. We note that for smaller admittance region, i.e., larger values of θ min , the probability of success increases. For small θ min the increase in admittance region causes devices to incur more repetition, hence, increasing the interference level, which overshadows the performance enhancement as shown for θ min = 5 • . Subsequently, to find the average success probability for a given admittance region bounded by ϕ max , we decondition (16) over the distribution of ϕ as follows It is illustrated in Fig. 7 where the success probability of the repetition scenario clearly outperforms the no repetition scenario. Note that this performance is for a single satellite spot, i.e. single admittance region. The following section extends results to a constellation of satellites.

D. Satellite Constellation Performance
Assuming most of the satellites are positioned in a non-GEO (Geosynchronous Equatorial Orbit), which are in a relative motion with respect to the ground IoT devices. As such, when reducing the size of the admittance region, the probability of having the coverage spot of the device is reduced. Accordingly there is a trade-off between satellite availability and its success rate within the spot. To capture this trade-off, we extend the performance by including the satellite availability, which is the ratio of satellite admittance area to the total surface area of Earth's sphere. For a small number of satellites, denoted k, and assuming non-overlapping spots, then the probability that a device is within the coverage is given as where 0 ≤ p spot ≤ 1. Accordingly, the probability of a successful frame delivery in a constellation of k satellite with non-overlapping spots is given by We depict in Fig. 8 the global average success probability evaluated for one tuning factor and various minimum elevation angles. We can clearly note the trade-off induced by adjusting the admittance region, i.e. by controlling θ min . Moreover, as the shape of the curve is impacted by the tuning factor a, an optimal point for a given admittance angle exists. In the following section we shed light on the joint optimization of these two parameters.
IV. DISCUSSION ON RESULTS In order to understate the joint impact of the admittance region size and the tuning factor on the global coverage probability, we depict in Fig. 9 the global success probability while varying factor a and admittance angle θ min . The solid line represents the optimal combination values for both parameters. Other factors that impact this optimal solution are related to the initial density of the ground devices λ o and the propagation conditions manifested in channel parameters since these parameters influence (19). As such, network designer can optimize the admittance region and tuning factor based on the operating parameters. To validate our analytical results we utilize Monte Carlo simulations that consider non-overlapping circular orbit satellites. Other simulation parameters include a satellite altitude of 550 km, noise average power of -138 dBm [6], center frequency at 2000 MHz, device p t at 23 dBm as utilized in many low power IoT devices, and a bit rate of 8 Kbps resulting in D o of 1 × 10 −6 .  This paper presented an analytic model that captures the performance of IoT-over-satellite in the uplink for non-GEO orbits using NTN radio channel. The model employs frame repetitions to enhance the success rate of transmissions. A trade-off between the satellite coverage region, admittance region, and satellite availability is formulated in order to evaluate the performance of a constellation of satellites. Accordingly, we found an optimal combination of parameters that maximizes the frame success rate. These results can facilitate network design and optimization for IoT-over-satellite links.


Graph Neural Networks as Multi-view Learning

Graph Neural Networks (GNNs) have demonstrated powerful representation capability in semi-supervised node classification. In this task, there are often three types of information -- graph structure, node features, and node labels. Existing GNNs usually leverage both node features and graph structure by feature transformation and aggregation, following end-to-end training via node labels. In this paper, we change our perspective by considering these three types of information as three views of nodes. This perspective motivates us to design a new GNN framework as multi-view learning which enables alternating optimization training instead of end-to-end training, resulting in significantly improved computation and memory efficiency. Extensive experiments with different settings demonstrate the effectiveness and efficiency of the proposed method.

INTRODUCTION
Graph is a fundamental data structure that denotes pairwise relationships between entities in a wide variety of domains (Wu et al., 2019b;. Semi-supervised node classification is one of the most crucial tasks on graphs. Given graph structure, node features, and labels on a part of nodes, this task aims to predict labels of the remaining nodes. In recent years, Graph Neural Networks (GNNs) have proven to be powerful in semi-supervised node classification (Gilmer et al., 2017;Kipf & Welling, 2016;Velickovic et al., 2017). Existing GNN models provide different architectures to leverage both graph structure and node features. Coupled GNNs, such as GCN (Kipf & Welling, 2016) and GAT (Velickovic et al., 2017), couple feature transformation and propagation to combine node feature and graph structure in each layer. Decoupled GNNs, such as APPNP (Klicpera et al., 2018), first transform node features and then propagate the transformed features with graph structure for multiple steps. Meanwhile, there are GNN models such as Graph-MLP (Hu et al., 2021) that extract graph structure as regularization when integrating with node features. Nevertheless, the majority of aforementioned GNNs utilize node labels via the loss function for end-to-end training.
In essence, existing GNNs have exploited three types of information to facilitate semi-supervised node classification. This understanding motivates us to change our perspective by considering these three types of information as three views of nodes. Then we can treat the design of GNN models as multi-view learning. The advantages of this new perspective are multi-fold. First, we can follow 2 THE PROPOSED FRAMEWORK We use bold upper-case letters such as X to denote matrices. X i denotes its i-th row and X ij indicates the i-th row and j-th column element. We use bold lower-case letters such as x to denote vectors. The Frobenius norm and trace of a matrix X are defined as X F = ij X 2 ij and tr(X) = i X ii . Let G = (V, E) be a graph, where V is the node set and E is the edge set. N i denotes the neighborhood node set for node v i . The graph can be represented by an adjacency matrix A ∈ R n×n , where A ij > 0 indices that there exists an edge between nodes v i and v j in G, or otherwise A ij = 0. Let D = diag(d 1 , d 2 , . . . , d n ) be the degree matrix, where d i = j A ij is the degree of node v i . The graph Laplacian matrix is defined as L = D − A. We define the normalized adjacency matrix asÃ = D − 1 2 AD − 1 2 and the normalized Laplacian matrix as L = I −Ã. Furthermore, suppose that each node is associated with a d-dimensional feature x and we use X = [x 1 , x 2 , . . . , x n ] ∈ R n×d to denote the feature matrix.
In this work, we focus on the node classification task on graphs. Given a graph G = {A, X} and a partial set of labels Y L = {y 1 , y 2 , . . . , y l } for node set V L = {v 1 , v 2 , . . . , v l }, where y i ∈ R C is a one-hot vector with C classes, our goal is to predict labels of unlabeled nodes. The labels of graph G can also be represented as a label matrix Y ∈ R n×C , where Y i = y i if v i ∈ V L and Y i = 0 if v i ∈ V U . The subscript U and L denote the sets of unlabeled and labeled nodes, respectively.

MULTI-VIEW LEARNING FOR GNNS
For the node classification task, we take a new perspective that considers node feature X, graph structure A, and node label Y as three views for nodes, and model graph neural networks as multiview learning. In particular, we need to jointly model each view and integrate three views. To achieve this goal, we introduce a latent variable F inspired by a traditional multi-view learning method (Xia et al., 2010). Then the loss function can be written as: where F is the introduced latent variable shared by three views, D X (·, ·), D A (·, ·) and D Y (·, ·) are functions to explore node feature, graph structure, and node label, respectively. These functions can contain parameters which we denote them as Θ. Hyper-parameters λ 1 and λ 2 are introduced to balance the contributions from these three views. One major advantage of the multi-view learning perspective is that it enables immense flexibility to design GNN models. Specifically, based on (1), there are numerous designs for D X (·, ·), D A (·, ·) and D Y (·, ·). Examples are shown below: • D X is to map node features X to F. In reality, we can first transform X before mapping. Thus, feature transformation methods can be applied including traditional methods such as PCA (Collins et al., 2001;Shen, 2009) and SVD (Godunov et al., 2021), and deep methods such as, MLP and self-attention (Vaswani et al., 2017). We also have various choices of the mapping functions such as Multi-Dimensional Scaling (MDS) (Hout et al., 2013) which preserves the pairwise distance between X and F and any distance measurements.
• D A aims to impose constraints on the latent variable F with the graph structure. Traditional graph regularization techniques can also be employed. For instance, the Laplacian regularization (Yin et al., 2016) is to guide a node i's feature F i to be similar to its neighbors; Locally Linear Embedding (LLE) (Roweis & Saul, 2000) is to force the F i be reconstructed from its neighbors. Moreover, modern deep graph learning methods can be applied, such as graph embedding methods (Perozzi et al., 2014;Grover & Leskovec, 2016) and Graph Contrastive Learning (Zhu et al., 2020;Hu et al., 2021), which implicitly encodes node similarity and dissimilarity. • D Y establishes the connection between the latent variable F L and the ground truth node label Y L for labeled nodes. It can be any classification loss function, such as the Mean Square Error and Cross Entropy Loss.
In this work, we set the dimensions of the latent variable F as R n×C , which can be considered as a soft pseudo-label matrix. Then the following designs are chosen for these functions: (i) for D X , we use an MLP with parameter Θ to encode the features of node i as MLP(X i ; Θ), and then adopt the Euclidean distance to map F i as MLP(X i ; Θ) − F i 2 2 ; (ii) for D A , Laplacian smoothness is imposed to constrain the distance between one node's pseudo labels F i and its neighbors as to constraint the pseudo label of a labeled node close to its ground truth.
These designs lead to our multi-view learning framework for graph neural networks (MULTIVIEW4GNN). Its loss function can be written in the matrix form as: where the first term maps node features into the label space, the second term indicates the pseudo labels should be smooth over the graph, and the last term constrains that the pseudo labels should be close to the ground-truth labels for labeled nodes.
Remark. There are recent works (Zhu et al., 2021;Yang et al., 2021) that aim to provide a unified optimization framework for understanding the message passing mechanism of different GNNs and designing new graph filter layers. However, they only focus on the forward process without taking the backward learning process into consideration, and they still follows the existing GNN architecture with end-to-end training. In this work, we do not aim to understand the message passing and design new GNN layers based on existing architectures. Instead, MULTIVIEW4GNN is a new graph deep learning framework as multi-view learning.

AN ALTERNATING OPTIMIZATION METHOD FOR MULTIVIEW4GNN
It is difficult to find the optimal solution of the loss function (2) for both F and Θ simultaneously due to the coupling between the latent variable F and model parameters Θ. The alternating optimization (Bezdek & Hathaway, 2002) based iterative algorithm can be a natural solution for this challenge. Specifically, for each iteration, we first fix the model parameters Θ and update the shared latent variable F on all three views. Then, we fix F and update the parameters Θ, which is effective in exploring the complementary characteristics of the three views. These two steps alternate until convergence. Next, we show the alternating optimization algorithm in detail.
Update F. Fixing MLP, we can minimize L with respect to the latent variable F using the gradient descent method. The gradient of L with respect to F (i.e., F U and F L ) is The gradient descent update of F with step sizes η L and η U is: According to the smoothness and strong convexity of the problem, we set η L = η U = 1 2(λ1+λ2+1) to ensure the decrease of loss value L (Nesterov et al., 2018), and the update becomes: Update Θ. Fixing F k+1 , we can minimize the loss function L with respect to MLP parameters: which equals training the MLP with the soft pseudo labels via the Mean Square Error Loss. Besides, we can also apply the Cross-Entropy Loss, and the details are in Appendix A.
Alternating Optimization and Scalability.
The multi-view learning perspective allows us to derive the alternating optimization solution which provides better flexibility in training than the end-to-end training in existing GNNs. The alternating optimization solution is a highly efficient and scalable training strategy, resulting in significantly improved computation and memory efficiency. Specifically, variable F and MLP model parameters Θ can be optimized separately. We can update F once and then train MLP multiple times. Meanwhile, there is no gradient backpropagation through the feature aggregation process so the aggregation steps do not need to store the activation and gradient values, which saves a significant amount of memory and computation. Moreover, due to the alternating updating scheme, we can use stochastic optimization to update F and train MLP by sampling the graph structure and node features. This can further improve the memory and computation efficiency as proved theoretically and empirically in stochastic optimization (Lan, 2020). In particular, only first-order neighbors are sampled in each update step of F, which avoids the neighborhood explosion problem in training large-scale GNNs (Hamilton et al., 2017;Fey et al., 2021). nd-to-end models, which can reduce the training epochs and correspondingly improve the training efficiency.

UNDERSTANDINGS OF MULTIVIEW4GNN
Another important advantage of alternating optimization is that it provides helpful insights to understand MULTIVIEW4GNN. In particular, based on the updating rules of F and Θ, we can naturally draw the following understandings on MULTIVIEW4GNN.
Understanding 1: Updating F is a feature-enhanced Label Propagation. Label Propagation (LP) (Zhou et al., 2003) is a well-known graph semi-supervised learning method based on the label smoothing assumption that connected nodes are likely to have the same label. LP can be written in an iteration form: k is the propagation step, and α is a hyper-parameter. Comparing LP with our update rule for F in Eq. (7) and Eq. (8), we can find that MLP(X) is involved in the propagation where MLP(X) can be regarded as labels generated by node features. In other words, our update rule for F takes advantage of node features, graph structure and labels while LP only uses graph structure and labels.
Understanding 2: Updating Θ is a pseudo-labeling approach. Pseudo-labeling (Lee et al., 2013;Arazo et al., 2020) is a popular method in semi-supervised learning that uses a small set of labeled data along with a large amount of unlabeled data to improve model performance. It usually generates pseudo labels for the unlabeled data and trains the deep models using both the true labels and pseudo labels with different weights. From this perspective, MULTIVIEW4GNN uses the pseudo labels F to train Θ such that it can take advantage of both labeled and unlabeled nodes. Figure 1: An overview of the proposed MULTIVIEW4GNN method, where grey color represents unlabeled nodes, and each node is associated with input features.

IMPLEMENTATION DETAILS OF MULTIVIEW4GNN
In this subsection, we detail the implementation of MULTIVIEW4GNN. As shown in Figure 1, we first preprocess the node feature through a diffusion, then alternatively update pseudo label F and MLP while taking into account the weight of pseudo labels and the class balancing problem, and finally get the prediction of unlabeled nodes. Next, we describe each step in detail.
Preprocessing. From Understanding 1, we use MLP to enhance the label propagation, so a good initialized MLP is needed. In real graphs, labeled data are usually scarce so that it is challenging to get a good initialization of MLP with a small number of labels. Therefore, we first diffuse the original node features with its neighbors to get smoothing and enhanced features. The new features are obtained from X = LP(X, α).
Then, we train MLP only using the labeled data for a few epochs to get an initialization, similar to pseudo-labeling methods (Iscen et al., 2019;Lee et al., 2013).
Update F. We initial F 0 = Y. Then we update F for labeled nodes and unlabeled nodes by Eq. (7) and Eq. (8), respectively. Since F acts as pseudo labels when training the MLP, we normalize F to be the distribution of classes by using the softmax function with temperature after the update: where τ is a hyperparameter to control the smoothness of pseudo labels. Pseudo-label certainty and class balancing. Directly using all pseudo labels to train MLP is not appropriate due to the following reasons. First, not all pseudo labels have the same certainty.
Second, pseudo-labels may not be balanced over classes, which will impede learning. To address the first issue, we assign a confidence weight to each pseudo-label. According to information theory, entropy can be used to quantify a distribution's uncertainty, so we define the weight for unlabeled nodes as To deal with the class imbalance problem, we adopt a simple method that chooses the same number of unlabeled nodes for each class with the highest weight to train MLP.
Update MLP. We train MLP using both labeled nodes set L and high confidence unlabeled nodes set U t . The loss function can be written as follows: where is a MSE loss and Θ is the parameters of MLP. Prediction. The inference of our method is based on the pseudo labels F, and the predicted class for the unlabeled node i can be obtained as c i = arg max j F ij . The overall algorithm, implementation details and code of MULTIVIEW4GNN are shown in Appendix B.
Suppose that p is the number of propagation layers, n is the number of nodes, m is the number of edges, and c is the number of classes. For simplicity, we assume that the hidden feature dimension is a fixed d for all transformation layers, and we have c d in most cases; all feature transformations are updated t epochs; Besides, the adjacent matrix A is a sparse matrix, and both forward and backward propagation have the same cost. Following , we only analyze the inherent differences across models by assuming that all models have the same transformation layers, allowing us to disregard the time required for feature transformation and the memory footprint of network parameters. The time and memory complexities are summarized in Table 1.
Time complexity. We first analyze the time complexity of feature aggregation. If we omit the difference in the dimension of the propagation features (d = c), the time and memory of GCN and APPNP are the same, as they require feature propagation in each epoch. We call the methods that need propagation every epoch as Persistent propagation methods. Similarly, the methods that only need to propagate once, such as SGC, SIGN (Rossi et al., 2020), and C&S (Huang et al., 2020) have the same time and memory complexity, namely One-time propagation methods. MULTIVIEW4GNN is a Lazy propagation method since the features are propagated k times during training with k being a small number. Thus, MULTIVIEW4GNN can be seen as a balance between these two groups of methods.

EXPERIMENT
In this section, we verify the effectiveness of the proposed method, MULTIVIEW4GNN, through the semi-supervised node classification tasks. In particular, we try to answer the following questions: • RQ1: How does MULTIVIEW4GNN perform when compared to other models? • RQ2: Is MULTIVIEW4GNN more efficient than state-of-the-art GNNs? • RQ3: How do different components affect MULTIVIEW4GNN?

EXPERIMENTAL SETTINGS
Datasets. For the transductive semi-supervised node classification task, we choose nine common used datasets including three citation datasets, i.e., Cora, Citeseer and Pubmed (Sen et al., 2008), two coauthors datasets, i.e., CS and Physics, two Amazon datasets, i.e., Computers and Photo (Shchur et al., 2018), and two OGB datasets, i.e., ogbn-arxiv and ogbn-products (Hu et al., 2020). For the inductive node classification task, we use Reddit and Flikcr datasets (Zeng et al., 2019). More details about these datasets are shown in Appendix C.
Following , we use 10 random data splits for the three citation datasets, and we run the experiments 3 times for each split. We report the average performance and standard deviation. Besides, we also test multiple labeling rates, i.e., 5, 10, 20, 60, 30% and 60% labeled nodes per class, to get a comprehensive comparison.
For other datasets, we use the fixed split and run 10 times.

PERFORMANCE COMPARISON ON BENCHMARK DATASETS
Transductive Node Classification. The transductive node classification results are partially shown in Table 2. We leave results on more datasets and methods in Appendix D due to the space limitation. From these results, we can make the following observations: • MULTIVIEW4GNN consistently outperforms other models at low label rates on all datasets. For example, in Cora and CiteSeer with label rate 5, our method can gain 1.2 % and 5.6 % relative improvement compared to the best baselines. This is because the pseudo labels generated by our framework are helpful for training models when there are few labels available. When the label rate is high, our method is also comparable to the best results. In addition, MULTIVIEW4GNN is alternately optimized but not end-to-end trained, which suggests that end-to-end training could not be necessary for node semi-supervised classification. • MULTIVIEW4GNN performs the best on OGB datasets. For example, in ogbn-products, it obtains 7.86% and 2.63% relative improvement compared to APPNP and SIGN, respectively. • Compared with the One-time propagation methods , Persistent propagation methods usually perform better when the labeling rate is low. In addition, the label propagation outperforms MLP in most cases, indicating the rationality of our proposed feature-enhanced label propagation. • The standard deviation of all models is not small across different data splits, especially when the label rate is very low. It demonstrates that splits can significantly affect a model's performance. A similar finding is also observed in the PyTorch-Geometric paper (Fey & Lenssen, 2019).
Inductive Node Classification. For inductive node classification, only training nodes can be observed in the graph during training, and all nodes can be used during the inference (Zeng et al., 2019). For MULTIVIEW4GNN, we first train an MLP with the training nodes' features and then do inference for the unlabeled node using the feature-enhanced label propagation in Eq (7) and Eq (8). Table 3, the MULTIVIEW4GNN outperforms other baselines on the inductive node classification task. The only difference between MULTIVIEW4GNN and MLP is the feature-enhanced label propagation, and the performance improvement can demonstrate its superiority.

EFFICIENCY COMPARISON
In this subsection, we compare the efficiency of our MULTIVIEW4GNN with other baselines, based on two large datasets, i.e., ogbn-arxiv and ogbn-products. To make a fair comparison, we choose the identical feature transformation layers for each method as there are no learnable parameters in feature propagation layers. Besides, we train model parameters with the same iterations in each method, i.e., 500 epochs for ogbn-arxiv and 1,000 epochs for ogbn-products. All the experiments are conducted on the same machine with a NVIDIA RTX A6000 GPU (48 GB memory).
For MULTIVIEW4GNN, we can update F with different frequency in training, i.e., 1, 2, 3, 4, 5, and "Full". "Full" means we update both F and MLP in each epoch. For MULTIVIEW4GNN-k, we only update the F for k times during the training procedure. The overall results are shown in Table 4.   MULTIVIEW4GNN with only a few update steps, such as MULTIVIEW4GNN-5, can achieve better accuracy with a minor increase in training time. For example, the whole training time of MULTIVIEW4GNN-5 is only 0.92s and 7s longer than SGC, but it has 8.63% and 11.24 % relative performance improvements in ogbn-arxiv and ogbn-products datasets, respectively. Meanwhile, we can observe that MULTIVIEW4GNN-5 has very similar performance with MULTIVIEW4GNN-Full, which suggests that there is no need to do propagation and train the model simultaneously for each epoch. This also suggests that the end-to-end training with propagation might not be necessary.
Memory Cost. Compared with the Persistent propagation methods, MULTIVIEW4GNN requires less memory with no requirement to store the hidden states in the propagation layers. Thus, MULTI-VIEW4GNN can keep a constant memory even with more propagation layers. Compared with MLP and SGC, MULTIVIEW4GNN only slightly increases memory as it needs to store the pseudo label matrix as analyzed in Section 2.5. In addition, MULTIVIEW4GNN is suitable for large-scale datasets that cannot fit in the GPU memory. First, it is easy to propagate the pseudo labels using CPUs as the label dimension is often lower than that of features. Then, MULTIVIEW4GNN is amenable to sampling training with mini-batch, which would significantly reduce the memory cost as discussed in Section 2.2. The only additional cost beyond One-time propagation methods is that we need to  transfer the result of MLP from GPU to CPU to do feature enhanced label propagation. Due to the small dimension of labels and the limited number of propagations during training, the cost is negligible.

ABLATION STUDY
In this subsection, we conduct ablation studies to gain a better understanding of how each component of our method works which correspondingly answers the third question. Feature Diffusion. It is expected that feature diffusion can improve the accuracy of the MLP in the pretraining procedure when the label rate is low and thus improve the quality of pseudo labels F during its following update steps. To validate this, we remove the feature diffusion step and also use the pseudo labels to train our method, which is called MULTIVIEW4GNN-w/o-diffusion. Experiments are conducted on both Cora and CiteSeer datasets. From Figure 2, we can see that at low label rates, MULTIVIEW4GNN is better than MULTIVIEW4GNN-w/o-diffusion which means that feature diffusion can boost the model's performance on low label rate setting. As the labeling rate increases, the performance gap becomes small, especially in the CiteSeer dataset. This shows that feature diffusion is not the key component in our method when the label rate is not very low.
Pseudo Labels. One of the most important advantages of MULTIVIEW4GNN is that we leverage pseudo labels to better train MLP. To study the contribution of pseudo labels in MULTIVIEW4GNN, we test the model variant MULTIVIEW4GNN-w/o-pseudo which only uses labeled data on Cora and CiteSeer datasets. Compared with MULTIVIEW4GNN, Figure 2 shows that pseudo labels have a large impact on model performance on both datasets, especially when the label rate is low.
Moreover, we choose the top K confidence pseudo labels per class after the first update of F to verify their accuracy. We adopt the same way to evaluate Label Propagation on Cora dataset with the label rate 20. As shown in Figure 3, after the first update of F, the accuracy of the top 180 nodes from each class can be 90%. So it is reasonable to use these pseudo labels to train MLP. Besides, the accuracy of our method at each K is much better than Label propagation, which suggests the effectiveness of the feature-enhanced label propagation update for F.
Hyperparameters Sensitivity. We test the parameter sensitivity of λ 1 and λ 2 in Eq (2) on Physics and Photo datasets by fixing one with the best parameters and tuning the other. From Figure 5, MULTIVIEW4GNN is not very sensitive to these two hyperparameters at the chosen regions.

RELATED WORK
Graph Neural Network (GNN) is an effective architecture to represent the graph-structure data. Two essential operations in the GNN are feature propagation and feature transformation. Considering how may times feature propagation in the training procedure, we categorize GNNs into Persistent propagation GNNs and One-time propagation GNNs. Persistent propagation GNNs (e.g., GCN (Kipf & Welling, 2016), GraphSAGE (Hamilton et al., 2017), GAT (Velickovic et al., 2017)) require feature propagation on each training step. GCN (Kipf & Welling, 2016) is the most commonly used method which performs feature transformation then feature aggregation in each layer. More recently, decoupled GNNs are proposed to alleviate the over-smoothness problem (Li et al., 2018;Oono & Suzuki, 2019). APPNP (Klicpera et al., 2018) is the first work to first apply multiple feature transformations then multiple aggregations. Similar architectures are also utilized in 2020;Zhou et al., 2021). However, these decoupled GNNs still require feature aggregation in each training step. One-time propagation GNNs are more efficient than the above Persistent propagation methods for they only propagate once despite the number of training steps. SGC (Wu et al., 2019a) is the first one to do feature aggregation and then transformation. SIGN (Rossi et al., 2020) adopts a similar strategy with a different aggregation scheme. Recently, more efficient neural networks are proposed which utilize graph structure for post-process. Huang et al. (2020) train a base predictor on labeled data and then apply a correct and a smooth step to post-process. Dong et al. (2021) further understand the relationship between decoupled GNNs and label propagation and utilizes soft pseudo labels for training. PPRGo (Bojchevski et al., 2020) precomputes the PageRank matrix but it lose much edge information due to aggressive sparsification of the PageRank matrix, which often leads to degraded performance and non-trivial tradeoff between efficiency and accuracy.
Multi-view learning on Graph. Our proposed multi-view learning framework for graph representation learning differs from existing works in literature. Specifically, multi-view graph cluster (Wang et al., 2019;Pan & Kang, 2021) consider a totally different setting where multi-view attributes and multiple structural graphs exist.
Multi-view graph contrastive learning algorithms (Hassani & Khasahmadi, 2020; use data augmentation to generate different graph views, then encourage the similarity between different views generated from the same graph while reduce the similarity in other view pairs. In contrast, our multi-view learning framework considers node features, graph structure, and node labels as three views of the nodes. Unified understanding on GNN. Recent works (Zhu et al., 2021;Yang et al., 2021) aim to provide a unified optimization framework for understanding the message passing mechanism of different GNNs and designing new graph filter layers. However, they only focus on the forward process without taking the backward learning process into consideration, and they are still following the existing GNN architecture with end-to-end training.
In this work, we do not aim to understand the message passing and design new layers based on existing architectures. Instead, MULTIVIEW4GNN is a new graph deep learning framework as multi-view learning. It provides a new perspective for graph representation learning with better flexibility, explanability, and efficiency.

CONCLUSION
In this work, we provide a new perspective to view the three types of information available for node classification (i.e., graph structure, node feature, and node label) as three views of nodes. This understanding inspires us to design GNN models as multi-view learning. The proposed MULTIVIEW4GNN framework can naturally be trained with the alternating optimization algorithm. Experimental results validate that MULTIVIEW4GNN is both computational and memory efficient with promising performance on the node classification task especially when the label rate is low.

A MULTIVIEW4GNN WITH CROSS ENTROPY LOSS
In section 2.2, we instance MULTIVIEW4GNN with a Mean Square Error Loss. In this section, we show it can be replaced by a Cross Entropy Loss. By replacing the first part of Eq (2) to be a cross entropy between MLP and F, the formulation becomes: where CE(·, ·) is the cross entropy function. Adopting the same gradient decent method, the update rule for F becomes: where the MLP is the output after the Softmax function and the step size can be set as η L = η U = 1 2(1+λ2) . Therefore, the update rule of F is: Then we can consider the hidden variable F as pseudo label and update the parameters of MLP based on the cross entropy loss CE(MLP, F). In practice, there are some situations that the cross entropy have a better performance than the original mean square error loss. Using two different losses would give similar results in most cases, and we report the best one.

B ALGORITHM OF MULTIVIEW4GNN
In this section, we provide the algorithm 1 and code of MULTIVIEW4GNN.
In line 1, we first initialize the pseudo label F as label matrix Y. In line 2, we preprocess data by feature diffusion. In line 3, we pretrain the MLP on the labeled data for a few epochs. From lines 4 to 9, we update F and MLP alternatively and iteratively.

C DATASET STATISTICS
In the experiments, the data statistics used in Section 3 are summarized in Table 5. For Cora, CiteSeer and PubMed dataset, we adopt different label rates, i.e., 5, 10, 20, 60, 30% and 60% labeled nodes per class, to get a more comprehensive comparison. For label rates 5, 10, 20, and 60, we use 500 nodes for validation and 1000 nodes for test. For label rates 30% and 60%, we use half of the rest nodes for validation and the remaining half for test. For Ogbn-arxiv dataset, we use the original fixed data split.