extfieldsofstudy:
Computer Science
Mathematics
Medicine

text:
Classification from Triplet Comparison Data

Learning from triplet comparison data has been extensively studied in the context of metric learning, where we want to learn a distance metric between two instances, and ordinal embedding, where we want to learn an embedding in a Euclidean space of the given instances that preserve the comparison order as much as possible. Unlike fully labeled data, triplet comparison data can be collected in a more accurate and human-friendly way. Although learning from triplet comparison data has been considered in many applications, an important fundamental question of whether we can learn a classifier only from triplet comparison data without all the labels has remained unanswered. In this letter, we give a positive answer to this important question by proposing an unbiased estimator for the classification risk under the empirical risk minimization framework. Since the proposed method is based on the empirical risk minimization framework, it inherently has the advantage that any surrogate loss function and any model, including neural networks, can be easily applied. Furthermore, we theoretically establish an estimation error bound for the proposed empirical risk minimizer. Finally, we provide experimental results to show that our method empirically works well and outperforms various baseline methods.

Introduction
Recently, learning from comparison-feedback data has received increasing attention (Heim, 2016;Kleindeßner, 2017). It is usually argued that humans perform better in the task of evaluating which instances are similar, rather than identifying each individual instance (Stewart et al., 2005). It is also argued that humans can achieve much better and more reliable performance on assessing the similarity on a relative scale ("Instance A is more similar to instance B than to instance C") rather than on an absolute scale ("The similarity score between A and B is 0.9 while the one between A and C is 0.4") (Kleindeßner, 2017). Collecting data in this manner has the advantage of avoiding the problem caused by individuals' different assessment scales. On the other hand, the collected absolute similarity scores may only provide information on a comparison level in some applications, e.g., sensor localization (Liu et al., 2004). It was shown that keeping only the relative comparison information can help an algorithm be resilient against measurement errors and achieve high accuracy (Xiao et al., 2006).
In this paper, we focus on the problem of learning from triplet comparison data, which is a common form of comparison-feedback data. A triplet comparison (x a , x b , x c ) contains the information that instance x a is more similar to x b than to x c . As one example, search-engine query logs can readily provide feedback in the form of triplet comparisons (Schultz and Joachims, 2004). Given a list of website links {A, B, C} for a query, if links A and B are clicked and the link C is not clicked, we can formulate a triplet comparison as (A, B, C).
Learning from triplet comparison data was initially studied in the context of metric learning (Schultz and Joachims, 2004), in which a consistent distance metric between two instances is assumed to be learned from data. The well-known triplet loss for face recognition was proposed in this line of research (Schroff et al., 2015;Yu et al., 2018).
Using this loss function, an inductive mapping function can be efficiently learned from triplet comparison image data. At the same time, the problem of ordinal embedding has also been extensively studied (Agarwal et al., 2007;Van Der Maaten and Weinberger, 2012). It aims to learn an embedding of the given instances to the Euclidean space that preserves the order given by the data. Algorithms for large scale ordinal embedding have been developed (Anderton and Aslam, 2019). In addition, many other problem settings have been considered for the situation of using only triplet comparison data, such as nearest neighbor search (Haghiri et al., 2017), kernel function construction  and outlier identification .
However, learning a binary classifier from triplet comparison data remained untouched until recently. A random forest construction algorithm (Haghiri et al., 2018) was proposed for both classification and regression. However, it first requires a labeled dataset and needs to actively access a triplet comparison oracle many times. For passively collected triplet comparison data, a boosting based algorithm (Perrot and von Luxburg, 2018) was recently proposed without accessing a triplet comparison oracle. However, a set of labeled data is still indispensable to initiating the training process. To the best of our knowledge, this paper is the first to tackle the problem of learning a classifier only from passively obtained triplet comparison data, without accessing either a labeled dataset or an oracle.

Contributions:
We show that we can learn a binary classifier from only passively obtained triplet comparison data. We achieve this goal by developping a novel method for learning a binary classifier in this setting with theoretical justification. We use the direct risk minimization framework given for the classification problem. We then show that the classification risk can be empirically estimated in an unbiased way given only triplet comparison data. Theoretically, we establish an estimation error bound for the proposed empirical risk minimizer, showing that learning from triplet comparison data is consistent. Our method also returns an inductive model, which is different from clustering and ordinal embedding, and can be applied to unseen test data points. The test data would consist of single instances instead of triplet comparisons since our primitive goal is to perform a binary classification task on unseen data points.
In summary, for the problem of classification using only triplet comparison data, our contributions in this paper are three-fold: • We propose an empirical risk minimization method for binary classification us-ing only passively obtained triplet comparison data, which gives us an inductive classifier.
• We theoretically establish an estimation error bound for our method, showing that the learning is consistent.
• We experimentally demonstrate the practical usefulness of our method.

Related Work
Our problem setting of learning a binary classifier from passively obtained triplet comparison data can be considered as a type of a weakly-supervised classification problem, where we do not have access to ground-truth labels (Zhou, 2017).
An approach based on constructing an unbiased risk estimator of the true classification risk from weakly-supervised data has been explored in many problem settings; for example, positive-unlabeled classification (du Plessis et al., 2014;Niu et al., 2016) and similarity-unlabeled classification (Bao et al., 2018) can be handled by the framework of learning from two sets of unlabeled data (Lu et al., 2018). Nevertheless, our problem setting is not a special case addressed by Lu et al. (2018) since we have only one set of triplet comparison data. We later show that we can formulate three different distributions, which is fundamentally different from the framework used by Lu et al. (2018).
Moreover, our problem setting is also different from similarity-dissimilarity-unlabeled classification (Shimada et al., 2019) in the sense that we have no access to unlabeled data and similarity and dissimilarity pairs, but only triplet comparison information. Fur-thermore, it is important to note that our problem setting is also different from preference learning (Fürnkranz and Hüllermeier, 2010), since we do not want to learn a ranking function but construct a binary classifier. Although we can first learn a ranking function and then decide a proper threshold to construct a binary classifier (Narasimhan and Agarwal, 2013), it is not straightforward to choose a proper threshold. Therefore, instead of this two-stage method, we focus on a method that can directly learn a binary classifier from triplet comparison data.

Learning A Classifier from Triplet Comparison Data
In this section, we first review the ordinary fully supervised classification setting. Then we introduce the problem setting and assumption for the data generation process of triplet comparison data. Finally, we describe the proposed method for training a binary classifier from only passively obtained triplet comparison data.

Preliminary
We first briefly introduce the traditional binary classification problem. We denote X ⊂ R d as a d-dimensional sample space and Y = {+1, −1} as a binary label space. In the fully supervised setting, we usually assume the labeled data (x, y) ∈ X × Y are drawn from the joint probability distribution with density p(x, y) (Vapnik, 1995). The goal is to obtain a classifier f : X → R that minimizes the classification risk where the expectation is over the joint density p(x, y) and ℓ : R × Y → R + is a loss function that measures how well the classifier estimates the true class label.
In the traditional fully supervised classification setting, we are given both positive and negative training data collectively drawn from the joint density p(x, y). However, in our case, we still want to train a binary classifier that minimizes the classification risk, although we do not have fully labeled data.

Generation Process of Triplet Comparison Data
We formulate the underlying generation process of triplet comparison data in order to perform empirical risk minimization. Three samples in a triplet are first generated independently, then shown to a user. The user can mark the triplet to be proper or not. A proper triplet means that the similarity between the first and second samples is stronger or the same as the similarity between the first and third samples. Specifically, it means that three labels (y a , y b , y c ) in a triplet appear to be one of the following cases: Otherwise, it means the first sample is more similar to the third sample than to the second sample; thus, the user chooses to mark the triplet as not proper. Similarly, it means (y a , y b , y c ) appears to be one of the following cases First, three data samples are generated independently from the underlying joint density p(x, y), then D = {(x a , x b , x c )} are collected without knowing the underlying true labels (y a , y b , y c ). However, we can collect information about which case a triplet belongs to from user feedback. After receiving feedback from users, we can actually obtain two distinct datasets. The data the user chooses to keep the order is denoted as Similarly, the data the user chooses to flip the order is denoted as Note that the ratio of n 1 |D 1 | to n 2 |D 2 | is fixed because we assume the three samples in a triplet are generated independently from p(x, y); thus, the ratio n 1 n 2 is only dependent on the underlying class prior probabilities, which are fixed unknown values.
The two datasets can be considered to be generated from two underlying distributions as indicated by the following lemma.
We denote the pointwise data collected from D 1 and D 2 by ignoring the triplet , the marginal densities of which can be expressed by the following theorem.
Theorem 1. Samples in D 1,a , D 1,c , D 2,a and D 2,c are independently drawn from and samples in D 2,b are independently drawn from A proof is given in Appendix B.
Theorem 1 indicates that from triplet comparison data, we can essentially obtain samples that can be drawn independently from three different distributions. We denote the three aggregated datasets as

Unbiased Risk Estimator for Triplet Comparison Data
We now attempt to express the classification risk, on the basis of the three pointwise densities presented in Section 3.2.
The classification risk can be separately expressed as the expectations over p + (x) and p − (x). Although we do not have access to data drawn from these two distributions, we can obtain data from three related densitiesp 1 (x),p 2 (x), andp 3 (x) as indicated in Theorem 1. Letting we can express the relationship between these densities as Our goal is to solve the above equation so that we can express p + (x) and p − (x) in terms of the three densities from which we have i.i.d. data samples. To this end, we can rewrite the classification risk, which we want to minimize, in terms ofp 1 (x),p 2 (x) and p 3 (x). An answer to Eq. (8) is given by the following lemma.
Detailed derivation is given in Appendix C.
As a result of the above lemma, we can express the classification risk using only triplet comparison data. Letting ℓ + (x) ℓ(f (x), +1) and ℓ − (x) ℓ(f (x), −1), we have the following theorem.
Theorem 2. The classification risk can be equivalently expressed as where π test p test (y = +1) denotes the class prior of the test dataset.
A proof is given in Appendix D.
In this paper, we consider the common case in which π test = π + , which means the test dataset shares the same class prior as the training dataset. However, even when π test = π + , which means the class prior shift (Sugiyama, 2012) occurs, our method can still be used when π test is known.

Estimation Error Bound
In this section, we establish an estimation error bound for the proposed unbiased risk estimator. Let F ⊂ R X represent a function class specified by a model. First, let R(F ) be the (expected) Rademacher complexity of F which is defined as where n is a positive integer, Z 1 , · · · , Z n are i.i.d. random variables drawn from a probability distribution with density µ, and σ = (σ 1 , · · · , σ n ) are Rademacher variables, which are random variables that take the value of +1 or −1 with even probabilities.
We assume for any probability density µ, the specified model F satisfies R(F ) ≤ Then for any δ > 0, with probability at least 1 − δ: where A proof is given in Appendix E.
Since n appears in the denominator, it is obvious that when the class prior is fixed, the bound will get tighter as the amount of triplet comparison data increases. However, it is not clear how the bound will behave when we fix the amount of triplet comparison data and change the class prior. Thus in Figure 1, we show the behavior of the coefficient term C R |ac−b 2 | with respect to the same class prior of both training and test datasets.
From the illustration, we can capture the rough trend that the bound gets tighter when the class prior becomes further from 0.5. We will further investigate this behavior in experiments.

On Class Prior
In the previous sections, the class prior π + is assumed known. For this simple case, we can directly use the proposed algorithm to separate test data as well as identify correct classes. However, it may not be true for many real-world applications. There are two situations that can be considered. For the worst case, no information about the class prior is given. Although we still can estimate a result for the class prior from data and obtain a classifier that is able to separate data for different classes, we cannot identify the correct class without the information of which class has a higher class prior. A better situation is that we have the information of which class has a higher class prior.
By setting this class as the positive class, we can successfully train a classifier to identify the correct class. Thus, we assume that the positive class has a higher class prior, which means π + > 1 2 .

Experiments
In this section, we conducted experiments using real world datasets to evaluate and investigate the performance of the proposed method for triplet classification.
6.1 Baseline methods KMEANS: As a simple baseline, we used k-means clustering (Macqueen, 1967) with k = 2 on all the data instances of triplets while ignoring all the relation information.
ITML: Information-theoretic metric learning (Davis et al., 2007) is a metric learning method that requires pairwise the relationship between data instances. From a triplet (x a , x b , x c ), we constructed pairwise constraints as (x a , x b ) being similar and (x a , x c ) being dissimilar. Using the metric returned by the algorithm, we conducted k-means clustering on test data. We used the identity matrix for prior knowledge and fix the slack variable as γ = 1.
TL: Triplet loss (Schroff et al., 2015) is a loss function proposed in the context of deep metric learning which can learn a metric directly from triplet comparison data.
Using the metric returned by the algorithm, we conducted k-means clustering on test data.
SERAPH: Semi-supervised metric learning paradigm with hyper sparsity ) is a metric learning method based on entropy regularization. We formulated a pairwise relationship in the same manner as with ITML. Using the metric returned by ITML, we conducted k-means clustering on test data.
SU: SU learning (Bao et al., 2018) is a method for learning a binary classifier from similarity and unlabeled data. We used the same method for estimating the class prior, and considered the less similar sample in a triplet as unlabeled data.

Datasets
UCI datasets: We used six datasets from the UCI Machine Learning Repository (Asuncion and Newman, 2007). They are binary classification datasets and we use the given labels for further triplet comparison data generation.
Image datasets: We used the following three image datasets.

Proposed method
For the proposed method, we used a fully-connected neural network with only 1 hidden layer of width 100 and rectified linear units (ReLUs) (Nair and Hinton, 2010) for all the datasets except for CIFAR-10. The width of the hidden layer was set to be 100 through out all experiments. Adam (Kingma and Ba, 2014) was used for optimization.
The neural network architecture used for CIFAR-10 is specified in Appendix.

Results
The proposed method estimates the unknown class prior first. For baseline methods, performances are measured by the clustering accuracy 1 − min(r, 1 − r) where r is the error rate. The results of different triplet numbers are listed in Tables 1, 2,and 3.
The best and equivalent methods are shown in bold face on the one-sided t-test with a significance level of 5%. Also as shown in Figure 2, the performance of the proposed method with respect to the class prior and the size of training dataset followed the prediction by the theory in most of the cases.

Conclusion
In this paper, we proposed a novel method for learning a classifier from only passively obtained triplet comparison data. We established an estimation error bound for the proposed method, and confirmed that the estimation error decreases as the amount of triplet comparison data increases. We also empirically confirmed that the performance of the proposed method surpassed multiple baseline methods on various datasets. For future work, it would be interesting to investigate alternative methods that can handle a   Proof. From the data generation process, we can consider the generation distribution for data of D 1 as Note that the denominator in Eq. (15) can be rewritten as then we have Moreover, the distribution p(x a , x b , x c , (y a , y b , y c ) ∈ Y 1 ) at the numerator of Eq. (17) can be explicitly expressed as from the assumption that three instances in each triplet comparison is generated independently.
Similarly, the underlying density for data of D 2 can be expressed as

B Proof of Theorem 1
Proof. For simplicity, we give the proof of D 2,a and the other 5 cases follow the similar proof. Noticing p 2 (x a , x b , x c ) = π + p + (x a )p − (x b )p + (x c ) + π − p − (x a )p + (x b )p − (x c ).
In order to decompose the triplet comparison data distribution into pointwise distribution, we marginalize p 2 (x a , x b , x c ) with respect to x b and x c : C Proof of Lemma 2 Proof. Notice that the equation has an infinite number of solutions. Letting we resort to finding the Moore-Penrose pseudo inverse (Moore, 1920;Penrose and Todd, 1954), which provides the minimum Euclidean norm solution to the above system of linear equations.
Let T * denote the conjugate transpose. We have


extfieldsofstudy:
Engineering
Computer Science

text:
Compressive Initial Access and Beamforming Training for Millimeter-Wave Cellular Systems

Initial access (IA) is a fundamental physical layer procedure in cellular systems where user equipment (UE) detects nearby base station (BS) as well as acquire synchronization. Due to the necessity of using antenna array in millimeter-wave (mmW) IA, the channel spatial information can also be inferred. The state-of-the-art directional IA (DIA) uses sector sounding beams with limited angular resolution, and thus requires additional dedicated radio resources, access latency and overhead for refined beam training. To remedy the problem of access latency and overhead in DIA, this work proposes to use a quasi-omni pseudorandom sounding beam for IA, and develops a novel algorithm for joint initial access and fine resolution initial beam training without requiring extra radio resources. We provide the analysis of the proposed algorithm miss detection rate under synchronization error, and further derive Cram\'er-Rao lower bound of angular estimation under frequency offset. Using QuaDRiGa simulator with mmMAGIC model at 28 GHz, the numerical results show that the proposed approach is advantageous to DIA with hierarchical beam training. The proposed algorithm offers up to two order of magnitude access latency saving compared to DIA, when the same discovery, post training SNR, and overhead performance are targeted. This conclusion holds true in various propagation environments and 3D locations of a mmW pico-cell with up to 140m radius.

I. INTRODUCTION
T HE millimeter-wave (mmW) communication is a promising technology for the future cellular network including 5G New Radio (5G-NR) [2]. Due to abundant spectrum, it is expected that the mmW network will support ultra-fast data rate. As shown in both theory and prototypes, mmW system requires beamforming (BF) with large antenna arrays at both base station (BS) and user equipment (UE) to combat severe propagation loss [3]. Significant differences in propagation characteristics and hardware architectures for mmW band compared to microwave band require novel signal processing techniques [4] and physical layer procedures [5].
Initial access (IA) is the fundamental physical layer procedure that allows UE to discover and synchronize with nearby BS before further communication. However, IA for mmW networks brings new challenges and opportunities as compared to IA for sub-6GHz band networks. In mmW system, conventional omni-directional IA with single antenna can not be reliable, and as a result IA needs to leverage transmitter and receiver antenna array to exploit BF gain [6], [7]. A key design challenge in mmW IA is the design of sounding beams for reliable discovery. In addition, beam training is required to achieve high BF gain enabled by large arrays and establish communication link. However, beam training now introduces additional access latency and signaling overhead due to repeated channel probing.

A. Related works
A number of works investigated various sounding beam designs and signal processing algorithms for mmW IA and beam training. Directional beams for IA and beam training are the most popular and extensively investigated in recent literature [6]- [15]. Directional IA (DIA) is first studied in [6] where a Generalized Likelihood Ratio Test (GLRT) is proposed to the solve the cell discovery problem under unknown multipleinput multiple-output (MIMO) channel and synchronization parameters. The authors concluded that the directional IA signal improves discovery range as compared to omni-directional IA. The DIA is further investigated in [8] where overhead and access latency are analyzed. Works [9] and [10] study DIA and its access latency in large networks using stochastic geometry. Impact of beam-width of sounding beams in DIA is researched in [13]. The comparison between omni-directional and DIA is also discussed in [11]. IA using out-of-band information, e.g., location, sub-6GHz measurement, are discussed in [7], [12]. The aforementioned works mostly focused on the overhead and latency for the cell discovery, while beam training is either not discussed or assumed to have coarse resolution [8]. It is common that DIA is paired with directional beam training [14], [15] where hierarchical sounding beams are used in multiple stages to achieve fine angular resolution for each user individually. However, such user-specific hierarchical sounding beams introduce prohibitive latency when a BS is connected to large number of UEs.
The alternative approaches for beam training are based on parametric channel estimation [16]- [23]. Exploiting the mmW sparse scattering nature, compressive sensing (CS) approaches have been considered to effectively estimate channel parameters based on channel observations obtained via various sounding beams. Works [16], [17] proposed a CSbased narrowband BF training with pseudorandom sounding beamformers in the downlink, and [18] extended this approach for a wideband channel. Other related works include channel covariance estimation [19]- [21] which requires periodic channel observations, and UE centric uplink training [22], [23]. It is worth nothing that all recent works focus on channel estimation alone while assuming perfect cell discovery and synchronization. The 5G-NR frame structure that supports IA is rarely considered, and further the feasibility of joint initial access and CS-based beam training has not been investigated.
There are also recent works that consider some practical aspects of IA. For example, frequency offset robust algorithms in narrowband mmW beam training are reported in [1], [24], [25]. There are several hardware prototypes that consider a practical approach of using received signal strength (RSS) in CS-based beam training. Channel estimation problem without phase measurement is a challenging problem, which was solved via novel signal processing algorithms based on RSS matching pursuit [26], Hash table [27], and sparse phase retrieval [28]. Note that phase free measurements were associated with a particular testbed, and this constraint does not necessarily apply to mmW systems in general. In summary, while IA and beam training algorithms have been extensively studied in the literature, there is a lack of understanding about the theoretical limits and signal processing algorithms that jointly achieve cell discovery and accurate BF training using asynchronous IA signal in mmW frequency selective channel.

B. Contributions
In this work, we propose to use quasi-omni pseudorandom sounding beams and novel signal processing algorithm to jointly achieve initial cell discovery, synchronization, and fine resolution beam training. More specifically, we provide answers to the following questions.
How to use pseudorandom sounding beams for IA? We propose an energy detection algorithm for initial discovery tailored for pseudorandom sounding beams. We derive the optimal detection threshold, analyze the miss detection probability and the impact of synchronization errors, i.e., carrier frequency offset (CFO) and timing offset (TO).
How to reuse received IA signal for beam training? We propose a novel CS-based beam training algorithm that reprocesses the frequency asynchronous IA signals to provide well aligned beam pair. We derive the Cramér-Rao lower bound (CRLB) of asynchronous training in line-of-sight (LOS) channel. We show that proposed algorithm reaches CRLB in LOS and remains effective in non-LOS (NLOS).
What are the benefits of compressive IA? We compare the proposed approach with DIA followed by hierarchical directional beam training. Key performance indicators for both approaches are numerically compared, including discovery rate, post beam training SNR, overhead and access latency. The simulation study based on 5G-NR frame structure and measurement-endorsed 3D 28GHz channel shows that the proposed approach is advantageous to DIA for UEs across wide range of locations in a small cell.

C. Organizations and notations
The rest of the paper is organized as follows. We start with a brief introduction of 5G-NR frame structure, IA and beam training in Section II. In Section III, we present the system model and problem statement. Section IV includes the proposed algorithm for cell discovery and timing acquisition followed by associated performance analysis. In Section V we present the algorithm and analysis for initial beam training and (x) are the real and imaginary parts of x, respectively. Set S = [a, b] contains all integers between a and b.

II. PRELIMINARIES: INITIAL ACCESS AND BEAM TRAINING
In this section, we introduce the mmW physical layer initial access procedure in 5G-NR cellular network. We briefly review the frame structure, synchronization sequences, and directional IA scheme as well as beam training. The reader is referred to work [7] for a more detailed survey.
Frame Structure: Fig. 1 shows the frame structure of 5G-NR. We focus on two functional blocks, namely synchronization signal (SS) burst and channel state information reference signal (CSI-RS). 5G-NR uses orthogonal frequency division multiplexing (OFDM), and the subcarrier spacing is either 120 or 240 KHz for mmW band. The SS signal is transmitted by a BS with period T F , typically 20 ms. The SS consists of up to M = 64 burst blocks. In each one of the burst blocks of duration T B , a specific sounding beam pair is used by BS and UE. The CSI-RS block with duration T r is dedicated to specific UE(s) for beam training and tracking. CSI-RS can use all frequency resources, i.e., up to B tot , and it has periodicity of T R , an implementation dependent value.
Synchronization Signal: Referring to Fig. 1, each SS burst has 4 OFDM symbols, i.e., primary synchronization signal (PSS), physical broadcast channel (PBCH), and secondary synchronization signal (SSS), followed by another PBCH. PSS is used in cell detection and synchronization, and it is assigned to the middle P = 128 subcarriers of the first OFDM symbol. The PSS in 4G-LTE is based on Zadoff-Chu (ZC) sequences due to their perfect cyclic-autocorrelation property and their Fourier duals [29], while in 5G-NR PSS is replaced by Maximum Length Sequences (M-sequences) [30]. There are N PPS = 3 and 336 unique sequences of PSS and MIMO channel at d-th delay sample a T (θ), a R (φ) Spatial responses of BS and UE φ l , θ l , g l , τ l Gain/AoA/AoD/delay of l-th multipath α l , β l Real and imaginary parts of g l s,s,s[n] F/T domain PSS vector and sequence vm, wm RF precoder/combiner of the m-th burst z[n], zm, σ 2 n AWGN sequence, vector, and power Initial discovery (detection) P FA Target FA prob. in initial discovery P MD,PT , P MD,NT MD prob. w/ and w/o perfect timing γ PT , η PT Detection stat. and TH w/ perfect timing γ NT , η NT Detection stat. and TH w/ unknown timing Initial beamforming traning (estimation) ξ Unknown parameters in BF training ym Received OFDM symbols at m-th burst d, t, r Vectors with candidates delay/AoA/AoD G D , G T , G R Parameter grid in delay/AoA/AoD est. Q, F ICI matrix and DFT matrix SSS, respectively, and these 1008 combination define the cell identifier (ID) of BS. PBCH carries control information.
Beamformed Initial Access: The BS periodically transmits IA blocks and such signals are processed by UEs which desire to establish the initial access, reconnect after beam misalignment, and search for additional BSs for potential handover. The sounding beams in SS bursts are intended to facilitate multi-antenna processing in BS and UE when no a priori channel information is available. Referring to Fig. 1, BS and UE in the DIA scheme use M T and M R transmitter and receiver beams to cover angular space at both ends. One T/Rx beam is used at a time, for all M = M T M R SS bursts.
Beam Training: The purpose of beam training is to identify the best beam pairs between BS and UE. The sounding beams in DIA typically have large beam-width and flat response inside angular sectors [31]. Such design covers the angular space of BS and UE within M bursts, but achieves coarse propagation directions estimation [8]. Thus DIA relies on directional beam training to refine angular resolution where BS and UE steer narrow sounding beams within the sectors of interest during CSI-RS periods.
III. SYSTEM MODEL This section introduces the system model that adopts the 5G-NR frame structure and problem formulation. All important notations are summarized in Table I sparse multipath channel to UEs. We focus on the IA and BF training procedure for a single UE. It is straightforward to extend it to multiple UEs since there is no UE-dependent processing. The UE uses analog array architecture, i.e., phased array, with N R antennas. We assume that a single stream of IA signal is transmitted by the BS regardless of its architecture. We first consider the received signal model when a UE searches for BS to initialize the connection. In this procedure, UE follows a periodic SS burst structure and uses predefined receiver beamformers to capture the signal according to [7]. As illustrated in Fig. 2, when the signal is present, the received samples, sampled at T s , is denoted as In the above equation, T is the unknown integer sample TO within range 1 0 ≤ T ≤ T,max ≤ N B , where T,max is the largest offset known to the system and N B is the number of samples in one SS burst, i.e., N B = T B /T s . The phase measurement error e j( F n+ψ[n]) comes from two sources. F is the normalized initial CFO, i.e., F = 2πT s ∆f where ∆f is absolute CFO in Hz between BS and UE. ψ[n] is the phase noise process in the UE receiver. N c is the maximum excessive multipath delay in discrete time, based on which cyclic prefix (CP) N cp > N c for OFDM symbols is designed. s[n] is the time domain signals of SS bursts. Referring to Fig. 2, we focus on the PSS and treat other symbols as zero [6], i.e., = v m , if n/N B = m. Beamformer switching is assumed not to introduce latency or phase offset in the transmission and reception. In this work, we focus on the system where each element of v m and w m is randomly and independently chosen from a set Such sounding beams require only 4-level phase quantization when steered by phased array and have randomized quasi-omnidirectional beam pattern.
The discrete time MIMO channel at delay where L and R are the number of multipath clusters (typically small, L ≤ 4 [32]) and sub-paths (rays), respectively. Scalar g l,r , τ l,r , θ l,r and φ l,r are the complex gain, excessive delay, angle of departure (AoD) and angle of arrival (AoA) of the rth sub-path within the l-th cluster, respectively. Function p c (t) is the time domain response filter due to limited temporal resolution T s . With half wavelength antenna spacing, the angular response vectors at the BS and UE are denoted as a T (θ) Note that the above model aligns with measurementendorsed mmMAGIC channel model [33] and is used for the system performance evaluation in Section VII. However, for the sake of tractable algorithm design and analysis, the following assumptions and definitions are made.
Assumption 1: BS and UE have ULA with omni-directional element pattern in 2D environment. Intra-cluster AoA, AoD, and delay offsets are zero, i.e., R r=1 g l,r g l , φ l,r = φ l , θ l,r = θ l , τ l,r = τ l , ∀r. Index r is omitted in the rest of paper for clarity. The phase error process is solely from CFO i.e., phase noise process is ψ[n] = 0, ∀n in (1). The complex path gain g l is deterministic complex value, i.e., L l=1 |g l | 2 = σ 2 g . Definition 1: The pre-BF signal to noise ratio (SNR) is defined as SNR σ 2 g /σ 2 n .

B. Problem formulations
We intend to address the following two problems, and their connection to the existing works are remarked.
Problem 1 (Initial Discovery and Timing Acquisition): The UE needs to detect the SS burst from in-band received samples (1). This problem is a binary hypothesis testing with unknown channel H[d] and synchronization errors T and F .
In addition, the TO T is estimated at this stage.
Problem 2 (Initial BF Training): The BF training is triggered once UE has detected IA signals. In this stage, UE reuses the asynchronous signal samples (1) to estimate the AoD and AoA of a path with significant power, say θ and φ , and which are then used as the steering vectors, v = a T (θ ) and w = a R (φ ), in data communications phase.
Remark 1: The above problems can be solved by DIA and directional beam training with the help of CSI-RS, while our solution relies on processing IA block only. In additional, although Problem 2 has overlap with parametric channel estimation, approaches from this class are not directly comparable. In fact, [18], [22] estimate the entire wideband channel, which facilitates optimal MIMO processing, but the assumptions of perfect synchronization and capability of wideband probing, i.e., B tot , do not apply in our work. Our goal is to provide wellaligned beam pair within IA block, i.e. without requiring CSI-RS slots. Finally, the cell ID recognition and PBCH decoding are important tasks but are not studied in this work.

IV. INITIAL DISCOVERY AND TIMING SYNCHRONIZATION
This section presents the proposed initial discovery and timing synchronization algorithm followed by their performance analysis.

A. Initial discovery and timing synchronization algorithm
The UE processes the received signal using the correlation filter with s zc [n], and obtains the detection statistics: Intuitively, there are M correlation peaks across M SS bursts.
The magnitude of the m-th peaks depends on the array gain of the m-th sounding beamformer, CFO, and TO. Our proposed detector combines energy from all M SS bursts and compares it with the threshold. In contrast to previous works [6], [9], [10] where the detection threshold is a fixed constant, we propose to use the optimal detection threshold based on Neyman-Pearson criterion that meets target false alarm (FA) rate P FA .
To understand the impact of timing synchronization error, we first consider a Genie scenario where the UE has perfect timing (PT) information, i.e., T = 0. In this case, the proposed PSS detection scheme is an energy detector over all M bursts. In addition, a sample time window with N c is used to collect energy from all multipaths. Specifically, the proposed hypothesis testing scheme is expressed as where the detection threshold η PT is used to reach false alarm rate constraint such that Pr(γ PT > η PT |H 0 ) = P FA . In a practical scenario without initial timing information (NT), i.e., T = 0, we propose to use the following detector that searches all possible instances within TO window T ∈ [0, T,max ] and uses the highest energy collected for the hypothesis test. The sample index corresponding to the highest energy in (5) is the estimate of TO, namelŷ B. Performance of initial discovery and timing acquisition In this subsection, we analyze performance of the proposed discovery algorithm in terms of miss detection rate, and the impact of initial synchronization error F and T . The exact expression is challenging and tedious, if not impossible, and therefore we provide a tight closed-form approximation in the following proposition. To be concise, the subscripts of γ and η that indicate the timing information assumption are denoted as binary variable E ∈ {NT, PT}.
Proposition 1: The optimal threshold of (5) that reaches where ξ E ( T,max , P FA ) is the threshold adjustment factor dependent on synchronization computed as where Q(.) and Q −1 (.) are Q-function and inverse Q-function, respectively. The associated miss detection (MD) rate P MD,E Pr (γ E < η E |H 1 ) using the optimal threshold η E is where the SNR degradation factor κ( T , F ) is defined as where K( T ) is the number of samples during PSS reception that UE switches beamformer due to TO.
Proof: See Appendix A. Remark 2: 1 − P MD,NT is a close approximation of probability that UE detects IA and correctly estimates T .
We gain two main insights from MD expressions (9) corresponding to threshold adjustment factor ξ z ( T,max , P FA ) and SNR degradation factor κ( F , T ). Firstly, the CFO affects MD performance by effectively reducing SNR via term κ( F , T ). Under maximum CFO at UE of ±5ppm and typical frame parameters P, M, N c specified in Section VII, the SNR degradation is bounded by 4 dB, i.e., 10 log 10 [κ( F , T )] ≥ −4dB, ∀ T . Secondly, the TO has impact on both factors. As seen in 2 Approximation is tight when TO search window size T,max ≥ 100.
(10), the SNR in the detection problem degrades when severe TO exists. In fact, K( T ) in κ( F , T ) models phenomenon that receiver sounding beam switches during the reception of PSS, i.e., K( T ) = 0. In addition, the presence of TO forces system to use peak detection scheme (5) where system searches peak location over a sample window with length T,max , i.e.,the worst case in (5). Under H 0 , the algorithm picks strongest noise realization over T,max samples and thus system needs to use higher threshold than in PT scenario, as seen in (7) and (8). Note that such degradation does not depend on the value of T , and the degradation in (9) is not critical with practical maximum TO uncertainty T,max ≤ N B . In summary, synchronization offset does not severely affect discovery performance of the proposed scheme.

C. Benchmark approach: directional initial discovery
For completeness, we briefly introduce the benchmark approach using directional sounding beam in initial discovery [6]. The system model of DIA is similar to Section III, except that sounding beamformers W and V are codebooks that steer directional sector beams, e.g., [14], [34]. Adapting the approach in [6] for the wideband channel and known PSS in SS burst, the cell discovery in DIA uses the following detector where γ DIA and γ DIA are the detection statistic and threshold in DIA. Sequenceỹ DIA [n] is the correlation output in (3) that corresponds to directional sounding beams. Refer to Fig. 1, the UE detects the burst with maximum power and denotes the index as m DIA which is used in directional beam training.
V. COMPRESSIVE INITIAL BEAM TRAINING This section presents the proposed initial access based BF training. We start with signal rearrangement based on information obtained from successful cell discovery and timing acquisition. Then, we introduce the CS problem formulation followed by the proposed algorithm. Finally, we analyze the CRLB of AoA/AoD estimation in LOS.

A. Signal rearrangement after timing acquisition
The further processing requires correct detection and CP removal, and therefore we make a following assumption.
Assumption 2: In beam training, the received IA signal (1) is correctly detected and TO T is correctly estimated.
The UE first removes CPs of P PSS samples from y[n] corresponding to M bursts and rearranges them into vector For notation convenience, in the rest of subsection, we restate the received time domain signal after CP-removal at the m-th SS burst y m ∈ C P according to the model in Section III, In the above equation, deterministic vector x m (ξ) ∈ C P is observations model of unknown parameters ξ [ F , · · · , θ l , φ l , τ l , α l , β l , · · · ] T , where α l = (g l ) and β l = (g l ). z m ∈ C P is the vectorized random noise. We also define x(ξ) = [x T 1 (ξ), · · · , x T M (ξ)] T . Specifically, in (14) vector s ∈ C P contains PSS symbols assigned to P subcarriers. Vector f (τ l ) ∈ C P is the frequency response corresponding to the excessive delay τ l of a multipath, i.e., the contribution of τ l on the p-th subcarrier is Matrix F ∈ C P ×P is discrete Fourier transform (DFT) matrix 3 . The effective channel gain is defined asg m,l = e j F NB(m−1) g l w H m a R (φ l )a H T (θ l )v m , and it includes the contribution of phase rotation across SS bursts due to CFO and IA beamformers v m and w m . Matrix Q( F ) = diag 1, e j F , · · · , e j(P −1) F T contains phase rotations within an OFDM symbol.

B. Baseline CS formulation
Directly estimating ξ from (14) via maximum likelihood (ML) requires multi-dimensional search with prohibitive complexity. In the following subsections, we re-formulate Problem 2 to facilitate sequential parameter estimation. With straightforward extension of the derivation in [4, Sec. V], the vector [g l ] m =g m,l in (14) can be re-formulated as whereÃ ∈ C GTGR×M is defined by the Hermitian conjugate of its m-th column as . Note that the above equation is different from [4, Sec. V] which requires M 2 sounding beam pairs. The matrixQ( F ) ∈ C M ×M contains the phase rotation in each SS burst due to CFO.
In fact, matrices A T ∈ C NT×GT and A R ∈ C NR×GR are the dictionaries of angular responses with AoAs and AoDs from grids with G T and G R uniform steps from −π/2 to π/2, respectively. In order words, the k-th columns in A T and where {r} k are the vectors that contain angle candidates.
Also note that the steps ∆θ and ∆φ depend on the desired resolution. In this work, G T and G R are used as number of steps and namely ∆θ = 2π/G T and ∆φ = 2π/G R . Matrix H l ∈ C GR×GT contains the complex path gain of the l-th path, i.e., it has 1 non-zero element whose location depends on the AoA and AoD of the l-th cluster in the angular grids. Remark 3: Assuming noisy observation ofg l and zero CFO, (16) reduces to the baseline problem in [4, Sec. V]. However, (14) implies that the former assumption is non-trivial unless s = 1, τ l = 0, ∀l, e.g., [17]. Moreover, algorithm designed with latter assumption is sensitive to CFO [1]. Finally, the AoA/AoD estimators are commonly confined in r and t [17]. We address these challenges in the following three subsections.

C. Effective gain estimation
To address the challenge discussed in Remark 3, we propose the following approach. We treat Q( F ) in (14) as identity matrix and estimate delay of dominant path and gain, say τ l andg m,l , by ML approach. Actually, the proposed algorithm uses sparse impulse support [d] q = q∆τ to construct a dictionary, where ∆τ = N c T s /G D is the step-size of delay candidates. Based on the knowledge of the model (15) and PSS signal s, the delay estimation is implemented aŝ where (15). The estimated delay tapτ enables estimating effective gain of a significant path byĝ where I M is the M × M identity matrix.

D. Joint AoA and AoD estimation robust to CFO
The second step uses a modified matching pursuit to solve CS problem (16) fromĝ while incorporating the existence of CFO inQ. In the conventional matching pursuit step, say the k-th, the anticipated effective channel response corresponding to an AoA and AoD pair, i.e., [Ã] k , is used to evaluate inner product with g [16]. The proposed heuristic treats AoA and AoD as known in the k-th step, and uses the ML estimator of CFOˆ F,k which is available in closed form. The modified matching pursuit is expressed aŝ whereã k [Ã] k from (16). The matrixQ(ˆ F,k ) has structure as (17). The inputˆ F,k is the ML CFO estimator when treating AoA/AoD as they correspond to ones in [Ã] k . Specifically, the CFO estimator relies on the estimator in [35] by treatinḡ y k =ã * k •ĝ as a tone with frequency F .
Operation ∠(x) = tan −1 [ (x)/ (x)] evaluates angle based on complex samples. To get estimates of the AoA, AoD, and CFO, indexk is used to select candidates from grids (18) after the following adjustmentk R = (k − 1)/G T + 1 and E. Off-grid refinement The aformentioned heuristics provide estimates of delay, AoA, and AoD that are restricted to the grid, i.e., elements of d, r and t. Grid refinement is a technique to provide off-grid estimation accuracy. There are several approaches considered in the literature including multi-resolution refinement [36] and the Newtonized gradient refinement [37]. In this work, we propose to use first order descent approach. As initialization of refinement, the estimator from previous steps is saved intô ξ (k) for k = 1. In the k-th iteration, the following error vector is evaluated where y is the received signal after rearrangement as (13), ) is obtained by plugging in estimated parameters into parametric model (14). In other words, e (k) is the error vector between observed signal sequence and received signal model using current estimates, which is then used to update parameters. The complex gain in iteration k is computed aŝ where ∇x g = (∂x(ξ)/∂g)| ξ=ξ (k) is the partial derivative of x(ξ) over parameter g in (14) evaluated atξ (k) . The refinement steps for delay, CFO, AoA, and AoD are moving towards the gradient of their estimators in the previous iterations. For concise notation, in the following equation and paragraph we use x to denote the parameter to be refined, i.e., x = {τ, F , θ, φ}. The refinement steps arê (27) where µ x is the step-size, vector ∇x x = (∂x(ξ)/∂x)| ξ=ξ (k) is the the partial derivative of x(ξ) in (14) over parameter of interest. The above approach iteratively runs by appending updated parameter into x(ξ (k+1) ) for the next iteration until the error e (k) 2 converges or falls below threshold 0 .
It is worth noting that the proposed approach can be extended to support multi-path training which has been covered by a variety of works in CS-based approaches [16], [18], [20]- [22], [24], [25]. However, the main motivation of this work is to showcase and analyze pseudorandom sounding beams in the initial access and initial beam training. Thus the only metric directly comparable to its counterparts [6]- [13], namely single path training, is evaluated.
The entire compressive initial beam training algorithm is summarized in Algorithm 1.
Proof: See Appendix B.

G. Benchmark approach: hierarchical directional BF training
The directional beams in SS burst allow BS and UE to coarsely estimate the propagation directions [38]. Although approach in [38] is not tailored for wideband channel with synchronization offset, it relies on RSS measurement within burst and therefore it is robust to the model mismatch. Using the SS burst index that corresponds to the maximum received power, the system uses the knowledge of directional sounding beams to infer channel propagation angles. Specifically, as illustrated in Fig. 1, the estimated θ and φ are the centers of them T -th andm R -th sounding beams in BS and UE [38], respectively. Note that the estimated angle sector indiceŝ m T andm R are computed from the SS burst index m DIA in (12), i.e.,m R = (m DIA − 1)/M T + 1, andm T = m DIA − (m R − 1)M T . The large width of a sector beam results in poor angular resolution in DIA. In order to improve the resolution, hierarchical directional beam training scans narrower beams within the sector of interest. Such procedure occurs during CSI-RS bursts which are scheduled for individual UEs.

VI. ACCESS LATENCY, OVERHEAD AND DSP COMPLEXITY
In this section, we present a model for analyzing three system performance indicators, namely access latency, overhead, and computational complexity. Note that this unified model applies to both directional scheme and the proposed approach.
Based on [7], we propose to use the latency model for both SS burst and CSI-RS as shown in Fig. 3. In both IA schemes, the failure of cell discovery introduces penalty of T F for a new IA block. When cell discovery occurs, the additional latency is required for scheduled CSI-RS according to the required number N train . Thus the access latency is where the first term includes latency for cell discovery. In the second term,T R is the average time for the UE to get the scheduled CSI-RS for beam training and it is expressed as In the above equation, N U denotes the number of UEs in the network. They share available CSI-RS in a time division manner to combat channel dynamic. Due to the limited number of CSI-RS K R = (T F − M T B )/T R within one IA period, more than one frame duration is required to meet scheduling of large number of UE N U . Therefore, in (29) K F = (N U −1)/K R is the number of frames required to assign all CSI-RS to UEs and K res = N U − K cyc K R is the residual delay in the last frame. As shown in the next section, DIA and directional BF training typically require larger N train than the proposed approach. Following [7], the overhead (OH) ratio is modeled by counting the time-frequency resource in IA and CSI-RS where B IA = 1/T s is the bandwidth in IA and the channel usage is M T s every period T F . We focus on varying CSI-RS density K R . Note that with reduced K R (increased T R ), the OH reduces with a cost of additional latency. The required baseband operations of the proposed approach are summarized in Table II, where only the complex multiplications are taken into account. As shown in later sections, the off-grid refinement in the proposed algorithm provides CRLB reaching accuracy but is necessary to reach adequate beam alignment. As a consequence, we do not include its complexity here. It is worth noting that the above analysis assumes there is an off-line pre-computation of all required dictionaries for matching pursuit, i.e., p q in (19),ã k in (22). In addition, the directional IA requires the first two steps in Table II.

Function Block Equation Operations
Initial Discovery PSS FIR corr.
(3) P N B Detection and time sync.
. RESULTS This section presents the numerical comparison between the proposed approach and DIA with directional beam training.

A. Simulation settings
The simulations follow 5G-NR frame structure. We first evaluate performance in the simplified 2D S-V channel model. The maximum excess delay is set as N c = 4 samples. As for the benchmark DIA, we use two approaches to design directional sector beams, i.e., least-squares based sector beamforming (LS-Sec) codebook [14] and frequency sampling method based sector beamforming (FSM-Sec) codebook [34, C23.4]. Examples of beam patterns 4 are shown in Fig. 4. In each of the Monte Carlo simulations, we generate an independent random realization of pseudorandom sounding beam codebook and channel parameters, unless otherwise mentioned.
We also evaluate the efficacy of the proposed approach in a realistic 3D mmW environment. In Section VII-C we simulate the system with QuaDRiGa simulator [39] based on mmMAGIC model [33] in 28 GHz urban-micro (UMi). We remove Assumptions 1, 2 from Sections III-A and V-A. Uniform planar arrays (UPA) N T = 16 × 4, N R = 4 × 4 are used at BS and UE, respectively, to exploit the higher sparsity in the elevation plane. The proposed algorithm is extended to UPA and 3D mmW channel model accordingly. In the simulations, the transmit power is set to P out = 46 dBm. The large scale channel model includes pathloss and shadowing. The AWGN on the receiver with 4dB noise figure is added with power of −170+10 log 10 (BW) dBm, where BW = 1/T s and BW = B tot = 400 MHz [40] for IA and data stage, respectively. Moreover, the UE phase noise, ψ[n] in (1), is modeled as Weiner process [41] that corresponds to oscillator with phase noise spectrum −114 dBc/Hz at 1 MHz offset [42]. The frame structure remains as described previously and other detailed simulations setting in QuaDRiGa can be found in the supplementary material [43]. The benchmark approaches are also extended for UPA and 3D channel, i.e., FSM-Sec beams are extended in both azimuth and elevation plane. During each one of N train CSI-RS, BS and UE use 16 sounding beams pairs which bisect previous scanned azimuth and elevation angular regions. We use post-training SNR as performance indicator, which is evaluated by dividing channel gain P out Nc d=0 |(w ) H H[d]v | 2 over noise power in B tot . Unless otherwise mentioned, the simulation parameters are summarized in Table. III.   [14], [34] with M T = 16 transmit sectors and one realization of 16 pseudorandom beams. In the polar plot, the r-axis refers to the gain in dB and the angular axis refers to steering angle in degrees. All patterns correspond to N T = 128 ULA.

B. Performance in simplified S-V channel model
The miss detection rate 5 of the proposed approach for initial discovery is shown in Fig. 5, and it is verified against the theoretical expressions (9). We have the following findings. Firstly, the lack of perfect timing synchronization introduces around 3 dB sensitivity loss as shown between the blue circled curve and red solid curve. However, this issue is unavoidable in practical systems. Secondly, less than 3 dB sensitivity loss occur when ±5 ppm CFO is present in addition to STO, as shown by the light blue dashed and green dashedand-dotted curves. Finally, the practical STO (≤ 10µs) is noncritical as shown by red solid and blue dashed curves. But when STO is large enough to cause transmitter and receiver burst beamforming window mismatch, e.g., 17 µs STO which corresponds to large K( T ) in (11), severe sensitivity loss is introduced as shown in grey dotted curves. In summary, these simulations verified the findings from Section IV that practical initial synchronization error introduces up to few dB sensitivity loss as compared to perfect synchronization scenario.
The comparison among proposed approach and benchmark DIA based discovery approaches is also presented in Fig. 5. Although common sense may doubt the efficacy of the proposed approach since there is no significant angular gain for any beam pattern, as illustrated in Fig. 4  there is only a couple of dB difference among the proposed approach and benchmark. However, such gap is less than the performance fluctuation of DIA with difference codebooks. The rationale behind this result is that the proposed scheme collects signal energy spread over all M SS bursts which in fact gives equivalent energy measurement as directional approach where energy collection occurs only when a sector beam aligns with true propagation direction. The beam training performance of the proposed BF training algorithm in LOS is presented in Fig. 6. The performance metrics are the residual mean square error defined by RMSE AoA = E|φ 1 − φ 1 | 2 and RMSE AoD = E|θ 1 − θ 1 | 2 . The simulations are conducted with Assumption 2. The same pseudorandom setting is used in both simulation and theoretical CRLB evaluation. The refinement steps are forced to terminate in up to 100 iterations. We have the following findings. Firstly, when the off-grid refinement are used, the proposed algorithm reaches CRLB in high SNR regime. Secondly, the coarse estimation in high SNR has a compromised performance as compared to CRLB. However coarse estimation (without refinement) has adequate accuracy for beam steering since RMSE is order of magnitude lower than 3 dB beam-width in steering, i.e., 0.29π/N T and 0.29π/N R . Finally, Fig. 5 and 6 reveal that in SNR region between −15 dB and −7.5 dB reliable detection occurs but beam training performance is poor. Admittedly, this implies a compromised experience for UEs at the cell edge, which is worth further investigation.
C. Performance in QuaDRiGa channel simulator Fig. 7 (a) illustrates the network setting implemented in QuaDRiGa. We simulate the performance of typical UEs distributed in two planes, with different distance towards the pico cell mmW BS. We present the following findings based on Fig. 7 (b), which shows the cumulative distribution function (CDF) of post-training beam steering SNR. Firstly, the pro- posed approach provides comparable performance to DIA with N train = 2 CSI-RS. In fact, in LOS, both approaches closely achieve beam steering towards true LOS path. Although the SNR seems excessively high in LOS, this implies that the transmit power can be reduced to save power. Secondly, DIA with less than N train = 2 CSI-RS has compromised SNR performance. This drawback is intuitive because wide sounding sector beam fails to extract precise angle information. The SNR improvement of using higher N train is more significant in LOS. Thirdly, although the proposed approach is tailored for sparse channels and presence of phase measurement error due to CFO, it is robust in NLOS scenarios where channel sparsity is compromised and practical phase noise occurs. Admittedly, the algorithm has a certain chance to completely fail when NLOS UEs are distributed in the second plane. However, in these cases the counterparts based on DIA and CSI-RS training cannot do much better job either. In fact, they have lower probability to reach post-training beam steering SNR above 0 dB compared to the proposed approach. The overhead and initial access latency savings of the proposed approach are significant, since it does not require CSI-RS, as shown in Fig. 7 (c). As explained in Section VI, for DIA based approaches when number of UEs in the network increases, the latency increases dramatically due to CSI-RS scheduling. Increasing the density of CSI-RS effectively reduces latency, but it results in increased overhead. The proposed approach relies on advanced signal processing to digitally conduct beam training and avoids requesting CSI-RS after initial access. In summary, up to two order of magnitudes saving in initial access latency is reached as compared to DIA.

D. Baseband processing requirements
Using the simulation parameters in Table III to evaluate  required operations in Table II, the baseband resource of the proposed method are in the same order of magnitude with DIA, i.e., (P N B +P G d +3M G T G R )/(P N B ) ≈ 7.2. There are two reasons for this finding. Firstly, exhaustive PSS correlation filter (3) is extreme computational demanding in IA. This filter 6 is required by IA regardless of sounding beam design. Secondly, the proposed approach sequentially estimates parameters and avoids multi-dimensional grid search.

VIII. DISCUSSION ON OPEN ISSUES
In this section, we discuss relevant issues in practical implementation of compressive IA and beam training.
Required a priori knowledge: Firstly, this work assumes coarse timing is available. It would be also important to study the case when timing is completely unknown, i.e., there is no a priori information about the range of T in (1), which could cause SS burst index misalignment to occur. Secondly, the compressive approach requires precise information about the sounding beam patternã k in (22). As a results, array geometry and sounding codebooks of both BS and UE need to be known a priori. This raises new challenges in communication protocol design to effectively incorporate this information. It also requires an increase in baseband operations if all dictionaries need to be computed on-the-fly. Further, mmW testbed experiments in [45] showed that the measured beam patterns commonly have mismatches with patterns predicted by codebook and array geometry model. Future research should address these impairments.
Channel sparsity: The efficacy of compressive approach is affected by the sparsity level in AoAs, AoDs, and multipath delays. Sparsity is endorsed by various mmW channel measurement campaigns, and urban NLOS, which is known with infavorable sparsity, is tested in this work. However, severely rich scattering situation are modeled from standard perspective, e,g, there are up to L = 20 multipath clusters in the 3GPP specified mmW channel [46]. It is important for system that utilize CS-based approach to flexibly handle situation when channel sparsity disappears.
Array architecture: This work focuses on the scenario where UE uses a single RF-chain to process a single stream of IA signals. This allows other RF-chains, if available at BS or UE, to operate in the band of data communication during IA. Since [6] shows that the hybrid analog/digital array and fully digital array are advantageous for DIA, it would be interesting to investigate benefits of compressive IA and beam training algorithm when they are adapted to utilize multiple RF-chain.
MIMO Multiplexing: The proposed beam training is compatible with multiuser multiplexing for hybrid array architecture. In fact, multiplexing designs [23], [47] rely on each RF-chain and corresponding analog beamformer to provide adequate post-BF SNR, and use the digital baseband processing to handle multi-beam interference. However, as mentioned in Remark 1, the comparison with channel estimation based approaches, i.e., estimation of the entire wideband channel or its covariance during CSI-RS for optimal MIMO processing, is rarely investigated.   (c) Access latency versus overhead of both IA scheme (29) and (30) with different number of UEs that share the scheduled CSI-RS. P MD is 0.04. Fig. 7. Initial access and beam training of proposed and benchmark approaches evaluated in 3D outdoor UMi network using 28GHz mmMAGIC channel model [33]. The trade off between post-training SNR in the data phase, required overhead, and access latency are also studied.
notable prototype [48] operates at 8GHz with two phased arrays synchronized by cabled reference clock. In addition to CFO, as emphasized in this work, the phase noise can also severely degrade coherency among channel observations. The phase noise detrimental impact becomes more severe with increased carrier frequency. Proper phase noise compensation as well as non-coherent CS-based beam training [26]- [28] are naturally immune to phase error and are worth investigation.

IX. CONCLUSIONS
In this work, quasi-omni pseudorandom sounding beam is proposed for the mmW initial access, synchronization, and beam training. We design associated signal processing algorithm based on the proposed sounding beam structure that is compatible with 5G-NR frame format. We provide theoretical analysis of cell discovery rate and beam training performance, and evaluate them via simulations using the mmW hardware and urban channel models from the literature that are supported by measurements. The results showcase that the proposed approach provides comparable performance to the state-of-the-art directional cell search for initial discovery, but provides significantly more accurate angle estimation during initial beam training. This advantage holds true across different propagation condition (LOS/NLOS) and UE-BS distance at 28 GHz band. Due to the saving of additional radio resource (CSI-RS) for beam refinement, the proposed approach reduces up to two order of magnitude access latency compared to the directional initial access when the same signaling overhead and post-training beam steering SNR are targeted.
All numerical results are reproducible with scripts in [43].