Ultra-Fast Accurate AoA Estimation via Automotive Massive-MIMO Radar

1 Ultra-F ast Accurate AoA Estimation via Automoti v e Massi v e-MIMO Radar Bin Li 1 , 2 , Shusen W ang 3 , Jun Zhang 2 , Xianbin Cao 4 , Chenglin Zhao 2 Abstract —Massive multiple-input multiple-output (MIMO) radar , enabled by millimeter -wav e virtual MIMO techniques, pro vides gr eat pr omises to the high-resolution automotive sens- ing and target detection in unmanned ground/aerial vehicles (UGA/U A V). As a long-established problem, howev er , existing subspace methods suffer from either high complexity or low accuracy . In this work, we propose two efficient methods, to accomplish fast subspace computation and accurate angle of arrival (AoA) acquisition. By leveraging randomized lo w-rank approximation, our fast multiple signal classification (MUSIC) methods, relying on random sampling and projection techniques, substantially accelerate the subspace estimation by orders of magnitude. Moreo ver , we establish the theoretical bounds of our proposed methods, which ensure the accuracy of the approxi- mated pseudo-spectrum. As demonstrated, the pseudo-spectrum acquired by our fast-MUSIC would be highly precise; and the estimated AoA is almost as accurate as standard MUSIC. In contrast, our new methods are tr emendously faster than standard MUSIC. Thus, our fast-MUSIC enables the high-resolution real- time en vironmental sensing with massiv e MIMO radars, which has gr eat potential in the emerging unmanned systems. Index T erms —Millimeter -wa ve radar , massive MIMO, auto- motive sensing, AoA estimation, subspace method, real-time. I . I N T RO D U C T I O N Unmanned aircrafts and automoti ve vehicles are recei ving more and more attention from industry and academia [1], [2], owing to its great potential in the widespread applications [3]. The success of unmanned systems critically depends on the fusion of Global Positioning System (GPS), automotive radar , and other en vironment sensors (e.g. lidar , ultrasound, and cam- era) [4], [5]. In comparison to others techniques, millimeter- wa ve (mm-wav e) automoti ve radar is extremely attractiv e to future unmanned systems, due to its two inherent merits. First, it is immune to adverse environment such as dust, fog, and smoke; and second, it is also robust to dazzling and no light conditions [6]. Furthermore, recent advancement in mm- wa ve semiconductor circuit (e.g. 24 and 77GHz) [7] makes it possible to deploy the large-scale arrays economically onto unmanned systems. As such, assisted by virtual Multiple-Input Multiple-Output (MIMO) or co-located MIMO techniques [8], hundreds of equivalent receiving channels are made av ailable to achieve the high-resolution environment sensing in dynam- ical automotiv e scenarios [9], [10]. 1 School of Information and Electronics, Beijing Institute of T echnology , Beijing, 100081, China. 2 School of Information and Communication Engineering, Beijing Uni ver - sity of Posts and T elecommunications, Beijing, 100876, China. 3 Department of Computer Science, Stev ens Institute of T echnology , Hoboken, NJ 07030, USA. 4 School of Electronic and Information Engineering, Beihang University , Beijing, 100191, China. In practice, unmanned vehicles use three types of mm-wav e radars [7]: (1) long-range radar (LRR) for automotiv e cruise control, in the range ∼ 200 m and with field of views (FO Vs) of azimuth ± 9 ◦ ; (2) medium-range radar (MRR) for cross- traffic alert ( ∼ 70 m, FO Vs ± 15 ◦ ); and (3) short-range radar (SRR) for park assist ( ∼ 10 m, FO Vs ± 75 ◦ ); as illustrated in Figure 1(a). T o meet div erse requirements, various modulation wa veforms are also designed for automotive radars [11], [12]. Popular solutions include frequency-modulated continuous- wa veform (FMCW) [13]–[15] and pulsed continuous wa ve (CW) radar, etc. Among these, FMCW sweeps a broad Radio Frequency (RF) bandwidth (GHz) while maintaining a small Intermediate Frequency (IF) bandwidth (MHz), which permits the high-resolution sensing via lo w-cost circuit [9], thus pro- viding great potential to the practical deployment. Despite the great advances in hardware integration and system design in automotiv e massi ve-MIMO radars [11], [16], [17], high-resolution and low-comple xity signal processing remains one substantial challenge. For massiv e-MIMO radars, the high-resolution en vironment sensing requires a large num- ber of channels [9], which incurs a high computation cost. Meanwhile, computational data processing causes an intolera- ble latency , which may lead to disastrous consequences. Thus, the deployment of massi ve MIMO automotive radar calls for high-resolution yet real-time processing algorithms. Theoretically , the maximum-lik elihood (ML) method achiev es the optimal accuracy [18], with an exhaustiv e search in two-dimensional large space (i.e. range and AoA). Provided a large number of samples, subspace methods, represented by Multiple Signal Classification (MUSIC) [19] and Estimation of Signal Parameters via Rational Inv ariance T echniques (ES- PRIT) [20], are able to attain the near-optimal accuracy [21], relying on signal or noise subspaces extracted by singular value decomposition (SVD) on a cov ariance matrix. Despite the high-resolution sensing, the high complexity , i.e. O ( M 3 ) for 1-D estimation and O ( N 3 M 3 ) for 2-D estimation ( M is the number of antennas and N is the number of snapshots), as well as the intolerable latency seriously limit the practical use [22], especially in massiv e-MIMO radars ( M > 200 ) 1 . In practice, the dimension of interests of MUSIC can be reduced by pre-estimation [23], which is still computational (e.g. an extended matrix in 2-D MUSIC is very large). T o simplify such subspace methods, a block-Lanczos scheme was applied [24], which focuses only on the first K singular v ectors and thus reduces the complexity to 1 For example, when the number of elements is 1000 as in emerging MIMO automotiv e radars, the processing latency would be around 500 ms even on high-performance CPUs. 2 O ( p M K M 2 ) , whereby p M ∼ O (log M ) is the number of blocks in the Krylov subspace and K is the number of targets. The so-called beamspace method, i.e. transforming element space to beam-space via a beamforming weight v ector [25], [26], can be used to simplify the subspace computation. Al- though it attains a promising performance in uniform circular array (UCA) [26], it may seriously degrade the AoA accuracy in the case of uniform linear array (ULA), as for automoti ve MIMO radars. Compressiv e sensing (CS) can be also applied to AoA estimation, which may be computational for massive MIMO radar [17], [27]–[29]. In ref. [30], the complex SVD was replaced with matrix in v erse [22]. Although the latter may be realized more ef ficiently , it has the same complexity O ( M 3 ) . When the signal-to-noise-ratio (SNR) is low , the estimated noise subspace would be e ven inaccurate, making the approximated pseudo-spectrum largely de viated. In ref. [31], the Propagator scheme was designed. By introducing a Propagator operator, it simplifies the complexity [32], yet scarifying the spatial resolution and the estimation reliability (e.g. in the range of [70 , 90] degree) [33]. As one popular so- lution, unkno wn targets can be quickly estimated via FFT [18], [34], [35]. Due to its low complexity , O ( M log 2 M ) , FFT has been widely used in automotive radars [10], [35], [36], which, howe v er , is insuf ficient for the high-resolution sensing and high-quality point cloud acquisition [11]. In this work, we break this major bottleneck in massiv e MIMO automotive radar . That is to say , we aim to accom- plish real-time automoti ve sensing at a scalable computational complexity , b ut achie ve the high-resolution target estimation as the standard MUSIC. T o achieve this, we apply randomized low-rank approximation (RLAR) to a cov ariance matrix, and dev elop two fast-MUSIC methods to approximately compute the signal subspace of it – the general stumbling block in all subspace-based methods. As such, we are able to acquire the highly accurate AoA estimation at a linear complexity – O ( K 2 M ) . Our fast-MUSIC would boost the widespread use of massive MIMO radars in the emerging unmanned systems. In summary , we offers the follo wing contributions. • W e introduce a novel concept, i.e. appr oximated rather than traditional exact matrix computation, to subspace methods. W e first approximate a large cov ariance S via three small matrices, in one special form S ' CW C H . Here, the matrix sketch C ∈ C M × p ( K < p  M ) is abstracted by uniform column-sampling on S ; while W ∈ C p × p is a weight matrix minimizing the approx- imation residue. Relying on such a randomized low- rank approximation, our fast-MUSIC reduces the time complexity of subspace computation to O ( p 2 M ) . • W e sho w our fast-MUSIC method enables the highly accurate pseudo-spectrum approximation and the precise AoA estimation. As the theoretical analysis suggests, when the signal subspace is approximately computed as in our fast-MUSIC, the estimated pseudo-spectrum q ˜ P music ( θ ) has a relativ e error O  σ K +1 ( S ) σ K ( S ) q M 2 p  , compared to the standard MUSIC pseudo-spectrum p P music ( θ ) , when the user -specify parameter meets p ∼ O ( K log K ) . Here, σ k ( S ) is the k -th largest singular value of S . If the signal-to-noise ratio (SNR) is relati vely high, the spectral gap σ K +1 ( S ) σ K ( S ) is very small, making our estimation highly accurate. • T o further improv e the approximation accuracy , we resort to another random projection technique, by iterati vely finding a good orthogonal projection matrix V ∈ C M × p . In contrast to a direct uniform sampling, we obtain an improv ed sketch C = SV ∈ C M × p and a low-rank approximation S ' CW C H . By incorporating more information into C , the accuracy of the approximated pseudo-spectrum is further enhanced. • W e establish the theoretical bound of our refined fast- MUSIC method. As shown, after t iterations, the approx- imated pseudo-spectrum q ˜ P music ( θ ) has a relativ e error O   σ K +1 ( S ) σ K ( S )  t +1 √ M 2 K  , compared to the standard MUSIC pseudo-spectrum p P music ( θ ) . In this case, even if the SNR is not very high, sev eral power iterations (e.g. t = 2 ) suffices for the near-optimal estimation. • W e ev aluate our fast-MUSIC methods with comprehen- siv e simulations, and compare them with its subspace counterparts (e.g. MUSIC, Lanczos, Propagator , etc). Our numerical studies v alidate the efficienc y and the accuracy of fast-MUSIC, which provides the great potential to high-resolution and real-time massive-MIMO radars in the emerging automotiv e applications. The remaining of this work is structured as follows. In Section II, we briefly introduce the system model and the subspace methods for automoti ve sensing. In Section III, we present two fast-MUSIC methods inspired by randomized matrix sketching. In Section IV , we deri ve some theoretical bounds of interests. In Section V , the numerical simulations are provided. Finally , we conclude this work in Section VI. The related notations are summarized as follows: A denotes an matrix with a rank r = rank( A ) ; its SVD is A = UΣV H = P r i =1 σ i ( A ) u i v H i , where σ i ( A ) ( > 0 ) is the i -th singular value; u i and v i are the i -th left and right singular vector . The Moore-Penrose pseudo-in v erse of A is A † . The Frobenius norm of A is k A k F = p tr ( A H A ) . The best rank- k ( k < r ) approximation to A is A k , P k i =1 σ i ( A ) u i v H i . I I . S U B S PAC E M E T H O D S I N M I M O R A DA R In this section, we briefly introduce FMCW MIMO radar of autonomous systems. On this basis, some popular subspace methods for high-resolution estimation are analysed. A. System Model – FMCW Radar A massi ve-MIMO radar system (e.g. virtual MIMO) con- sists of M recei ving antennas, i.e. typically the ULA. For 3 each element, the emitted FMCW wav eform within a symbol duration T sym reads [13], [35]: s ( t ) = exp  j  w s t + µ 2 t 2   , 0 ≤ t < T sym , (1) where w s is the initial frequency and µ is the changing rate of instantaneous frequency of the chirp signal. Given the bandwidth of FMCW signals w B , as well as the emission duration, T sym , the changing rate is µ = w B /T sym . After the signal calibration (i.e. mitigating non-ideal antenna effects based on the measured response), the recei ved signal at the m -th antenna reads: y m ( t ) = K − 1 X k =0 α k,m s ( t − τ k ) × exp  j 2 π λ s md sin θ k  + n m ( t ) . Here, r k ( t ) , α k,m s ( t − τ k ) is the signal reflected from the k -th target; τ k and θ k are times of arriv al and angles of arriv al; α k,m is the complex channel gain between the k -th target and the m -th element, and in the far -filed cases 2 we hav e α k,m = α k ; λ s denotes the wa velength of carrier; d is the spacing distance of two adjacent elements (we assume d = λ s / 2 ); n m ( t ) = 1 √ 2 [ N (0 , σ 2 n ) + j N (0 , σ 2 n )] denotes the complex additiv e Gaussian noise with a variance σ 2 n . A significant adv antage of FMCW MIMO radar is that it is capable of demodulating signals via one simple mixer [13], which greatly facilitates the low-cost RF implementation. T o be specific, at the reception-end the mixer-based de-chirping process is used to extract the beat signal ˜ y m ( t ) , i.e. ˇ y m ( t ) = s ( t ) ×   K − 1 X k =0 r k ( t ) exp  j 2 π λ s md sin θ k  + n m ( t )   . By substituting s ( t ) into eq. (1), and after a lo w-pass filter with an impulse response h ( t ) , the target signal can be obtained, which is a composition of target-modulated signals: ˜ y m ( t ) = K − 1 X k =0 α k e j ( µτ k t + ω s τ k − µ 2 τ 2 m ) · e j ( πk sin θ k ) + ˜ n m ( t ) . Here, ˜ n m ( t ) = h ( t ) ⊗ [ n m ( t ) s ( t )] is the residual noise after the de-chirping procedure; ⊗ denotes the conv olution process. Then, the analog to digital con vertor (ADC) with a sampling frequency f s = 1 /T s outputs the discrete signal: ˜ y m ( n ) = ˜ y m ( nT s ) , for m = 0 to M − 1 and n = 0 to N − 1 , with N , T sym /T s . On this basis, the estimation of AoA, range and v elocity of targets can be estimated from the discrete signal ˜ y m ( n ) . Combined with the massi ve MIMO techniques, the reso- lution of automoti ve radars can be ef fectiv ely improv ed. In practice, this can be implemented via virtual MIMO techniques [8]–[10], [12]. As in Figure 1-b, each subsystem in volv es 3 TX and 4 RX elements – the recei ving channels is then 12. When two sub-systems are cascaded ( N c = 2 ), the equiv alent 2 When the near-field ef fect occur , the subspace method would be still applicable, e.g. with the proper modification [37]. In such cases, our proposed methods would be applied directly to reduce the time complexity of subspace computation. L R R MRR SRR Automotive r adar Master Slave FMCW Clock Local Oscillator Clock 2 RX RX M M    2 TX TX M M    4 RX TX RX TX M M M M M      (b) (d) (c) TX RX Virtual MIMO Rada r Virtual MIMO Rada r RX TX RX (a) Fig. 1. Massive MIMO automotive radar with cascaded sub-arrays . (a) T ypical radars for autonomous driving systems, including LRR, MRR and SRR. (b) A schematic framework on cascaded MIMO radar , whereby the master generates the local oscillator signal and controls the slave via two synchronized clocks, i.e. FMCW clock and LO clock. Here, each co-located MIMO radar is equipped with M T X = 3 transmitting antennas and M RX = 4 recei ving antennas. (c) The equiv alent MIMO radar system of two cascaded sub-systems. (d) The cascaded-based massiv e MIMO radar system with the equiv alent channels of M = 4 × M RX M T X . receiving channels would be (3 × 2) × (4 × 2) = 48 , see Figure 1-c. In the general case, for the N c cascading system, the equiv alent receiving channels would be 12 N 2 c , as illustrated in Figure 1-d. Thus, hundreds of elements can be integrated to enhance the spatial resolution [9], [38]. Moreover , this leads to one compact MIMO radar (e.g. ∼ 10 cm array), thus alle viating the near-field effects to some extent. B. Subspace Detection & Estimation For automotiv e sensing, unkno wn AoA (i.e. θ k ) can be estimated by various methods [11], [12], [16]. When the high resolution estimation was emphasised as in automotiv e sensing, subspace methods tend to be more attractiv e [11], [17]. Popular subspace methods, such as MUSIC [19] and ESPRIT [20], all start from a cov ariance matrix, i.e. S = 1 / N × YY H , (2) whereby the signal matrix Y ∈ C M × N is arranged as: Y ,       ˜ y 0 (0) ˜ y 0 (1) · · · ˜ y 0 ( N − 1) ˜ y 1 (0) ˜ y 1 (1) · · · ˜ y 1 ( N − 1) . . . . . . . . . . . . ˜ y M − 1 (0) ˜ y M − 1 (1) · · · ˜ y M − 1 ( N − 1)       . (3) There are two important things to be noted for the above co- variance matrix. First, the matrix multiplication in computing S incurs also a high complexity , which yet can be ef ficiently calculated in parallel. Second, for massiv e-MIMO radars (with a large M ), the required snapshot number N will be very large, in order to estimate an unbiased cov ariance. For a small 4 length, the standard MUSIC pseudo-spectrum may be biased [39], e.g. due to the so-called subspace swapping problem [40]. Such theoretical limitations of the standard MUSIC may be ov ercame by some post-processing algorithms, e.g. correcting the estimated pseudo-spectrum via well-designed weights [41]. Recently , a single-snapshot MUSIC was dev eloped [42]. As shown, even in the extreme case N = 1 , the high-resolution AoA estimation can be attained. Note that, the emphasis of this work is on the subspace computation which is one prerequisite to the accurate estima- tion of AoA. Assume there are K unknown target points, we first compute SVD of the large cov ariance S ∈ C M × M , i.e. S = UΣV H = U K Σ K V H K +  2 n U − K V H − K , (4) where U K = [ u 1 , · · · , u K ] ∈ C M × K corresponds to the signal subspace, and U − K = [ u K +1 , · · · , u M ] ∈ C M × ( M − K ) corresponds to the noise subspace. For a Hermitian matrix S H = S , we further hav e V = U . As common, we assume σ 1 ( S ) > σ 2 ( S ) > · · · > σ M ( S ) , and σ m ( S ) =  2 n ( m > K ) . Then, we uniformly partition the spatial angle range [0 , π ] into a grid of L v alues, { θ 0 , θ 1 , · · · , θ L − 1 } . For each angle θ l , we define the M -dimensional steering vector as: a ( θ ) =  a 0 ( θ ) , a 1 ( θ ) , · · · , a M − 1 ( θ )  H where a m ( θ ) , exp  j 2 π d m · sin ( θ ) λ  , m = 0 , 1 , · · · , M − 1 . Finally , we are able to e valuate the standard MUSIC pseudo- spectrum, based on the computed subspaces, i.e. P music ( θ ) = 1 a ( θ ) H U − K U H − K a ( θ ) = 1 a ( θ ) H ( I M − U K U H K ) a ( θ ) . (5) And then, each peak in the pseudo-spectrum P music ( θ ) in- dicates one target point. So, unknown AoAs, θ k ( k = 0 , 1 , · · · , K − 1) , are attained by ˆ θ θ θ =  θ k , θ k = peak { P music ( θ ) }  . (6) In practice, the pseudo-spectrum P music ( θ ) needs to be computed for ev ery θ ∈ { θ 0 , · · · , θ L − 1 } . The computation of the whole pseudo-spectrum P music ( θ ) requires a time com- plexity O ( M K L ) . Fortunately , P music ( θ ) can be e valuated in a highly efficient manner [20], [43], which is hence out of the scope of this work. In the following, we are interested in the fast computation of a signal subspace U K . As seen, the in volv ed complexity mainly comes from the decomposition of large covariance S , which is measured by O ( M 3 ) ; and more importantly , there is no parallel algorithm to accelerate it. Thus, despite its high-resolution in automotive sensing, the processing latency of subspace methods would be intolerable for massiv e-MIMO radars, e.g. > 100 ms for M > 500 . C. Considerations on Subspace Methods As mentioned, there are many works which target at reducing the time complexity of subspace methods, which hav e dif ferent adv antages and disadv antages, as summarized in T able I. From this comparativ e result, one long-standing challenge in MIMO radar signal processing is that the high T ABLE I T I ME C O M P LE X I T IE S O F C L AS S I C AL D O A E S T I MATI O N M E T HO D S Methods Complexity Resolution Conditions MUSIC O ( M 3 ) high all Lanczos O ( p M K M 2 ) high K ≤ M Matrix-in verse O ( M 3 ) high  n → 0 Propagator O ( K M N ) medium | θ | ≤ 70 deg FFT Method O ( M log M ) low all F ast MUSIC 1 O ( M p 2 ) high K  M F ast MUSIC 2 O ( tpM 2 ) high K  M 1. Here, p ( K ≤ p  M ) is an over-sampling parameter; t ( t ≥ 1 ) is the number of iterations in updating the projection matrix. 2. Note that, the real complexity of our f ast method 2 is largely lo wer than O ( tpM 2 ) , as the polynimial complexity on M comes from the matrix multiplication which is yet more efficient than SVD or matrix inv erse. 3. p M ∼ O (log M ) is the number of blocks in constructing a Krylov subspace in the Lanczos method. resolution and the lo w complexity hav e the inherent conflict. Meanwhile, most existing algorithms manage to compromise such two contradictory performance aspects. I I I . F A S T M U S I C M E T H O D S In this section, we dev elop two highly ef ficient methods for subspace estimation. Note that, various low-rank approx- imation techniques, e.g. SVD, ha ve been widely applied to the subspace-based AoA estimation [44]–[46]. In contrast to previous works mainly balancing tw o contradictory aspects, here we design two fast-MUSIC algorithms which permits the highly accurate AoA estimation at a linearly scalable complexity . In principle, our new methods are inspired by the random sampling/projection techniques [47]–[50]. Despite their great interests in linear algebra or scientific computing, such randomized sketch methods have been rarely considered in the field of MIMO signal processing [38], [51]. In sharp contrast to our pre vious work [38] that adopts a less efficient sketch form and enjoys no theoretical bound, this would be the first attempt to develop ultra-fast subspace methods with the close-form error bound in the approximated pseudo-spectrum. A. F ast MUSIC Method – 1 In our fast-MUSIC algorithm, we firstly approximate the large covariance S ∈ C M × M with three small matrices (re- ferred as matrix sketches), in the special form: S ' CWC H (as in Figure 2). Meanwhile, in order to minimize the com- putational complexity , we resort to a simple random sampling technique to abstract the in v olved matrix sketch C ∈ C M × p ( K ≤ p  M ) from S . T o be specific, we first configure one over -sampling param- eter 3 , which is denoted as p ; and thus we uniformly sample p items from { 0 , 1 , · · · , M − 1 } to form the indexing set I , i.e. |I | = p . Then, the sketching matrix is abstracted as C = S ( : , I ) = SΠ , whereby Π ∈ R M × p is the equi v alent sampling matrix. As seen, C ∈ R M × p contains p columns of 3 Here, the term ov er-sampling indicates that the sampling parameter p should be always greater than the number of target points K , i.e. p > K . 5 =  M M   S  M p   C Π S p H H M    Π S C   p p  W { 1 , 2, 4, 9}, =4  s  Fig. 2. Illustration of randomized low-rank approximation of a large covari- ance matrix S ∈ R N × N , with rank ( S ) = 4 . Here, the randomly selected column indexes set is I = { 1 , 4 , 2 , 9 } , with |I | = 4 . Thus, the low- dimensional sketch C consists of p = |I | columns of S , while C H in volves 4 rows. The resulting approximation error is less than 10 − 12 . S indexed by I . In the equi valent sampling matrix Π , there are p columns contains only one non-zero entry p M /p (the remaining ( M − p ) columns are all zeros). By minimizing the approximation error , we further obtain another weight matrix W opt , i.e. W opt = arg min W k S − CW C H k 2 F = C † SC † H . (7) Since the in v olved pseudo-inv erse and matrix multiplication incur a high complexity , we tend to the other scheme to compute W ' W opt . T o be specific, we solve the above ov er-determined system by sketching both S and its random approximation. This leads to the so-called Nystr ¨ om approxi- mation [47], [52]; and the weight matrix is determined by: W = arg min W k Π H ( S − CW C H ) Π k 2 F = ( Π H SΠ ) † . (8) Finally , a large cov ariance matrix S is approximated by: S ' ˜ S , CWC H . (9) Relying on this low-rank representation of S , its SVD can be ef ficiently computed via multiple SVDs of these small sketches. Assume C = U c Σ Σ Σ c V H c , the cov ariance matrix S would be further approximated by: S ' U c Σ Σ Σ c V H c W ( U c Σ Σ Σ c V H c ) H . After defining B , Σ Σ Σ c V H c WV c Σ Σ Σ H c and computing its SVD, i.e. B = U B Σ Σ Σ B V H B ∈ C p × p , we would obtain the rank-restricted approximate SVD of S , i.e. S ' U c U B Σ Σ Σ B V H B U H c = ¯ U K Σ Σ Σ B ¯ U H K , whereby ¯ U K , U c U B is one unitary matrix. Relying on it, we finally obtain an approximated pseudo-spectrum, i.e. ˜ P music ( θ ) = 1 a ( θ ) H ( I M − ˜ U K ˜ U H K ) a ( θ ) . (10) The schematic flo w of our fast-MUSIC method 1 is sum- marized in Algorithm 1 , which actually obtains the rank- restricted SVD of a large cov ariance and thus reduces the time complexity significantly . Based on the above analysis, we easily see that: Step 1 (Lines 2 to 5) costs O ( pM ) time. Step 2 (Lines 6 to 9) costs O ( p 2 M ) time. Combined together , the overall complexity of our fast method 1 in estimating the signal subspace is O ( p 2 M ) . Algorithm 1 Fast MUSIC Method – 1. 1: Input : a cov ariance matrix S , the number of target points K , and over -sampling parameter p ( ≥ K ). 2: // Step 1 : Random Sampling Method 3: I ← − randomly sampling p indices from { 0 , 1 , · · · , M } ; 4: C ← − the column of S indexed by I ; 5: W ← − the pseudo-in verse of S ( I , I ) ; 6: // Step 2 : Rank Restriction 7: Compute the SVD of C = U c Σ Σ Σ c V H c ; 8: Compute the SVD: Σ Σ Σ c V H c WV c Σ Σ Σ T c = U B Σ Σ Σ B V H B ; 9: Compute the signal subspace: ˜ U K ← − U c U B ; 10: // Step 3 : Approximate MUSIC 11: Compute ˜ P music for all θ . B. F ast MUSIC Method – 2 The above random sampling method enables a highly ef fi- cient approximation of lar ge cov ariance as well as the pseudo- spectrum. In practice, this approximated pseudo-spectrum in eq. (10) is almost accurate, which suffices to attain the high- resolution DoAs estimation. Ne vertheless, one potential risk is that its error bound would be relativ ely weak in the worst case; see Section IV -C for the theoretical analysis. If the precise matrix approximation was emphasized, we may resort to the other random projection technique [53]. Here, the main difference with regard to our fast method 1 is that, when abstracting a sk etch C , a more informati ve projection matrix Π ∈ C M × p is adopted, i.e. C = SΠ , which will be computed in an iterative manner . Compared with the random sampling based matrix sketch as in Algorithm 1 , this random projection based sketch incorporates more information of S and permits the more accurate approximation. T o begin with, we first initialize one M × p random matrix Π , whose entries are i.i.d. values from the Normal distrib ution N (0 , 1) . Then, we obtain one randomly projected sketch C (1) = SΠ ∈ C M × p . T o improv e the quality of this sketching matrix, we further compute the orthonormal basis of C (1) , [ Q c , R c ] = qr ( SΠ ) , V (1) , Q c , where qr ( · ) denotes the QR decomposition with the time complexity O ( p 2 M ) ; and therefore, the M × p matrix V gi ves another impro ved projection matrix. Subsequently , we obtain one updated sk etch C (2) = SV (1) ∈ C M × p , and further obtain its orthonormal basis V (2) , i.e. [ Q c , R c ] = qr  SV (1)  , V (2) , Q c . (11) The abov e process would be repeated for t times, and finally we achiev e a sufficiently good matrix sketch C ( t ) = SV ( t ) ∈ C M × p . Similar to eq. (8), the p × p weight matrix is: W = ( V ( t ) H SV ( t ) ) † . W ith C = C ( t ) and W , we are able to obtain the projection- based matrix approximation, i.e. S ' CWC H . In a similar way , the SVD of S is approximately attained by multiple SVDs on the small sketches, i.e. S ' ˜ U K Σ Σ Σ B ˜ U H K . Finally , the pseudo-spectrum ˜ P music ( θ ) is estimated, as in eq. (10). The time complexity of our fast method 2 is measured by O ( tpM 2 ) , see Algorithm 2 . Note that, despite the similar random projection as in the Lanczos method, here our method 6 targets at deriving small sketches and a matrix approximation in eq. (9), which is faster than the block-Lanczos method, as sho wn in T able I. In our fast method 2, the user -specific parameter t trades of f the complexity and the accuracy . E.g., a large t impro ves the projection matrix and the matrix approximation accuracy , but also increases the complexity . One practical concern is that the noise may become non- Gaussian. In this case, a generalized eigen-decomposition problem can be similarly formulated [54], by resorting to high- order cumulates techniques. As a result, EVD/SVD can be still used to acquire unknown AoA. In this regard, our ne w methods are readily applicable to reduce the time complexity of subspace computation in the case of non-Gaussian noises. I V . P E R F O R M A N C E B O U N D S F O R F A S T - M U S I C Our primary concern is then the accuracy of the approxi- mated pseudo-spectrum ˜ P music ( θ ) relative to its exact version P music ( θ ) . In the following, we analyze the theoretical perfor- mance of our fast-MUSIC methods, giving the relative error between the approximated pseudo-spectrum in eq. (10) and the exact pseudo-spectrum in eq. (5) as standard MUSIC. It should be noted that, from a theoretical perspective, the existing theories or error bounds dev eloped for randomized matrix sketching are inapplicable to our concerned problem. This is to say , the ultimate aim of our fast-MUSIC method is to approximate an exact pseudo-spectrum P music ( θ ) , and then to acquire the high-resolution AoAs estimation. Unfortunately , most existing theories only bound the matrix approximation errors [47], [48], [53], e.g. k S − ˜ S k F or k S − ˜ S k 2 . Accordingly , such matrix norm bounds should not lend any support to our fast-MUSIC methods. In this section, we show the relativ e error of our estimated pseudo-spectrum ˜ P music ( θ ) would be also bounded, after ap- proximating S with ˜ S , which is more challenging than to deriv e the classical matrix norm bounds. Our Theorems 1 and 2 establish the theoretical lower bounds for ˜ P music ( θ ) . Theorem 3 establishes the theoretical upper bound for it. A. Lower/Upper Bounds of ˜ P music ( θ ) First, we expect to have the lower bound in the follo wing form: there is a bounded positiv e number α l such that ˜ P music ( θ ) ≥ 1 1+ α l P music ( θ ) , (12) for all θ . Such a lo wer bound is of great importance, as it guarantees the approximated pseudo-spectrum ˜ P music ( θ ) shall not miss the target peaks (see the theoretical analysis in Section IV -E). Theorems 1 and 2 are such lo wer bounds for our fast methods 1 and 2, respectively . Second, we also expect to have an upper bound in the following form: there is a positiv e number α u such that ˜ P music ( θ ) ≤ 1 1 − α u P music ( θ ) , (13) for all θ . Such an upper bound guarantees our approximated pseudo-spectrum ˜ P music ( θ ) shall not interpret the baseline fluctuation in the non-target region as one false target peak (see the analysis Section IV -E). Theorem 3 is such an upper bound for the fast method 2. For the fast method 1, we do not Algorithm 2 Fast MUSIC Method 2. 1: Input : spatial covariance matrix S ( M × M ) and target rank K . 2: // Step 1 : Random Pr ojection 3: Π ← − M × K matrix with entries i.i.d. drawn from N (0 , 1) ; 4: Random projection based matrix sketch: C ← − SΠ ; 5: Repeat V ← − orth ( C ) and C ← − SV for t times; 6: // Step 2 : Randomized Matrix Appr oximation 7: C ← − SV and W ← − ( V H SV ) † ; 8: Compute the SVD of C = U c Σ Σ Σ c V H c ; 9: Compute the SVD: Σ Σ Σ c V H c WV c Σ Σ Σ T c = U B Σ Σ Σ B V H B ; 10: Compute the signal subspace: ˜ U K ← − U c U B ; 11: // Step 3 : Approximate MUSIC 12: Compute ˜ P music ( θ ) for all θ . hav e such an upper bound at the current stage, which would be left to the future study . B. Lower Bounds for F ast Method 1 Let S K = U K Σ K U H K be the best rank- K approximation to S ; ˜ S = CW C H be randomized matrix approximation to S via uniform column sampling (i.e. C = SΠ and Π is the sampling matrix); and ˜ U be the orthonormal basis of ˜ S . Theorem 1. Let the notation be defined in the above . Let δ ∈ (0 , 1) be any user-specified constant; the r ow coher ence of U K be defined as µ ( U K ) , M K max i   ( U K ) i :   2 2 ∈  1 , M K  [53], [55], [56]. F or a user-specific sampling par ameter p ≥ 4 . 5 µ ( U K ) K · log K δ , the following r elation holds with pr obability at least (1 − δ ) : r P music ( θ ) ˜ P music ( θ ) ≤ 1 + 2 q M 2 p σ K +1 ( S ) σ K ( S ) . (14) Thus, Theorem 1 establishes the lower bound on ˜ P music ( θ ) . As found, the r elative approximation error is basically O  2 q M 2 p · σ K +1 ( S ) σ K ( S )  . If the SNR is high, we hav e σ K ( S )  σ K +1 ( S ) =  2 n ; and therefore; the approximation error could be very small. The proof of Theorem 1 can be found in Appendix C. C. Lower Bound for F ast Method 2 Let Π be an M × p initial projection matrix with i.i.d. entries drawn from N (0 , 1) . Let the M × p matrix V ( t ) be an improv ed projection matrix after t iterations; in other words, V ( t ) is the orthonormal basis of S t Π . By definition, C = SV ( t ) and W = ( V ( t ) H SV ( t ) ) † , and ˜ S = CWC H is an approximation to S . Let ˜ U be the orthonormal basis of ˜ S . Theorem 2. Let the notation be defined in the above . Repeat the iteration pr ojection for t ( t ≥ 0 ) times. Let δ ∈ (0 , 1) be any constant. The following r elation holds with probability at least 1 − δ : r P music ( θ ) ˜ P music ( θ ) ≤ 1 + √ M 2 K δ  σ K +1 ( S ) σ K ( S )  t +1 . (15) 7 Theorem 2 establishes a lo wer bound for ˜ P music ( θ ) which is obtained via random projection in Section (III-B). As seen, the relativ e approximation error is basically O  √ M 2 K ·  σ K +1 ( S ) σ K ( S )  t +1  . That means, the approximation errors exponentially con ver ge to zero with t . Even if the SNR is not very high, our fast method 2 con v erges to the precise estimation in se veral iterations (i.e. t ≥ 2 ). The proof can be found in Appendix D. D. Upper Bound for F ast Method 2 Theorem 3. Let the notation be defined in Section IV -C. Repeat the iterative projection for t ( t ≥ 0 ) times. Then, with pr obability at least 1 − δ , r P music ( θ ) ˜ P music ( θ ) ≥ 1 − √ M 2 K δ  σ K +1 ( S ) σ K ( S )  t +1 . (16) Theorem 3 establishes an upper bound for ˜ P music ( θ ) , which is also estimated via the random projection. Similarly , the relativ e approximation error is basically O  √ M 2 K ·  σ K +1 ( S ) σ K ( S )  t +1  . It also con verges to zero exponentially . Again, ev en if the SNR is not very high, few iterations suffice for the high precision. The proof of Theorem 3 is structured into Appendix E. E. Accurate DoA Estimation As shown latterly , in the context of massive-MIMO radars, the target peaks of pseudo-spectrum, which are located exactly at target AoAs, are significantly higher than noise baseline, P music ( θ k )  P 0 , where P 0 , max  P music ( θ 0 )  denotes the maximum ampli- tude of non-target regions, i.e. θ 0 / ∈ θ θ θ K × 1 . Without loss of generality , we assume it almost satisfies P music ( θ k ) /P 0 > γ , when the SNR is not very small (e.g. γ ∼ M ; see the following Figure 3). If the approximated pseudo-spectrum misses one target peak located at θ k , we should hav e q ˜ P music ( θ k ) ≤ p P 0 . (17) At the same time, according to our theoretical lower bounds in eq. (13), we are supposed to hav e 1 1+ α l p P music ( θ k ) ≤ q ˜ P music ( θ k ) ≤ p P 0 , and therefore, √ γ ≤ p P music ( θ k ) /P 0 ≤ 1 + α l . (18) On the other hand, according to the follo wing numerical analysis (e.g. in Figure 5), we would have 1 + α l → 1 and 1 + α l < √ γ . (19) By checking eq. (18) and (19), we immediately conclude this contradiction result indicates eq. (17) should be inv alid. I.e., the approximated pseudo-spectrum ˜ P music ( θ ) should not miss the true peaks of targets. Similarly , from the upper 25 30 35 40 45 50 AoA /deg 10 0 10 1 Pseudo-spectrum Fast MUSIC: p=K Fast MUSIC: p=2K Fast MUSIC: p=3K MUSIC Ground-truth 25 35 0.15 0.5 Fast MUSIC: p=K Fast MUSIC: p=2K Fast MUSIC: p=3K MUSIC Ground-truth (a) Pseudo-spectrum of the Nystr ¨ om method with differ- ent p . 25 30 35 40 45 50 AoA /deg 10 0 10 1 Pseudo-spectrum Fast MUSIC: t=0 Fast MUSIC: t=1 MUSIC Ground-truth 25 35 0.15 0.5 Fast MUSIC: t=0 Fast MUSIC: t=1 MUSIC Ground-truth (b) Pseudo-spectrum of the power method with different t . Fig. 3. W e plot the pseudo-spectra of standard MUSIC and the fast-MUSIC methods with different tuning parameters. bound in eq. (14), we show the approximated pseudo-spectrum ˜ P music ( θ ) should not mistake non-target AoA as one target peak, thus prev enting the false alarm in automotive sensing. V . N U M E R I C A L S I M U L A T I O N & A N A L Y S I S In the section, we examine our fast-MUSIC methods via nu- merical simulations; and compare them with popular subspace algorithms in high-resolution massiv e MIMO radars. A. Experiment Settings Since we focus on the subspace computation as well as sub- sequent AoA estimation, we directly start from a signal matrix Y of M receiving elements. In the following simulations, we configure the operation frequency to be 77GHz, as for mm- wa ve MIMO radar systems. The IF bandwidth is assumed to be 80MHz. The element spacing distance is equal to the half of wa velength. The time comple xity can be fairly measured by the CPU runtime (CPU 2.59GHz, 32GB memory), which is proportional to the required number of multiplication flops. The a veraged runtime is reported based on 100 independent realizations. B. T uning P arameters Our fast-MUSIC method 1 (Algorithm 1) in volv es one key user-specific parameter p ( p ≥ K ), which is used for trading 8 100 200 300 400 500 600 700 800 900 1000 Number of elements M 10 -4 10 -3 10 -2 10 -1 10 0 CPU time /sec MUSIC Matrix inverse Block Lanczos Randomized MUSIC Propagator Fast method 2 Fast method 1 Fig. 4. W all-clock runtime (second) against the number of antennas, M . W e calculate only the runtime for computing the signal subspace which is the major bottleneck in automotive radar . Here we set N = M and SNR = 0 dB. off the accuracy and the complexity . I.e., a lar ge p is necessary to attain the tighter error bound (as in Theorem 1), which in turns increases the time complexity . Like wise, the fast-MUSIC method 2 (Algorithm 2) in v olves a tuning parameter t ( ≥ 1 ). In the following, we discuss their practical settings, whereby we fix M = 200 , N = 2 M , K = 4 , and SNR = 0 dB. In Figure 3(a), we plot the approximated pseudo-spectrum of our fast MUSIC method 1 under different p . W e observ ed that, under three different settings ( p = K , 2 K , and 3 K ), the fast algorithm successfully obtains the accurate AoA estimation of K targets (e.g. with overlapped peaks). Define the approximation error relativ e to standard pseudo-spectrum as P L − 1 l =0  ˜ P music ( θ l ) − P music ( θ l )  2 . When p is increased from K to 2 K , the error drops from 0 . 59 to 0 . 11 ; while p is further increased to 3 K , the approximation error would be marginal. In Figure 3(b), we show the approximated pseudo-spectrum of our fast method 2 under different settings of t ; here we assume P = K . When t = 2 the approximation error is 0 . 03 , which is suf ficiently small for the accurate AoA estimation. In practical automotiv e sensing, we may directly set t = 2 . C. Comple xity vs Accuracy W e now ev aluate the runtime of v arious popular subspace methods. Figure 4 shows the computation latency (seconds) against the number of antenna elements, i.e. M . W e assume N = M and K = 10 in the numerical analysis. For our fast- MUSIC method 1, we assume p = 12 (i.e. p = d 1 . 2 × K e ); and in fast-MUSIC method 2, we have t = 2 and p = 12 . Despite the near-optimal accuracy in high-resolution au- tomotiv e sensing, the standard MUSIC has O ( M 3 ) time complexity , due to the computational SVD on lar ge co variance. As M increase, it easily becomes impractical for real-time automotiv e applications, by producing an intolerable latency . When M = 1000 , the latency of SVD itself consumes more than 500 milliseconds (ms), which is formidable to real- time sensing of unmanned aircrafts/vehicles; similar problems persist in matrix in verse and Lanczos methods. Here, we do not compare the ESPRIT method, as it relies on SVD and has exactly the same runtime as MUSIC in subspace computation. As shown in Figure 4, our f ast-MUSIC method 1 is the most efficient, with its time complexity linearly scalable with M , 0.2 0.4 0.6 0.8 1 1.2 0.2 0.4 0.6 0.8 1 1.2 (a) 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 0.5 1 1.5 (b) Fig. 5. Comparing the theoretical bounds (Theorems 1, 2, and 3) with the numerical results. (a) Theoretical bound of fast method 1. (b) Theoretical bounds of fast method 2, green cross – upper bound, blue cr oss – lower bound. The constants κ , κ l , and κ u are defined in Section V -D. i.e. O ( p 2 M ) . For a lar ge array M = 1000 , its required compu- tation latency is only 0.5 ms; while MUSIC consumes 500 ms. So, it would be 1000 times faster than the standard MUSIC. It is noteworthy that, although our fast method 2 theoretically has a time complexity O ( tpM 2 ) , its time cost, as discussed, is dramatically lower than the block-Lanczos algorithm and the other Randomized MUSIC [38]. As discussed, the iteration number in f ast method 2 is set to 2 , while p ∼ O ( K ) . As validated by Figure 4, its time complexity is much lower than the Lanczos method, whereby the parameter p M grows as O (log M ) (e.g. p M = 7 for M = 1024 ). Its processing latency is only 4 ms when M = 1000 , only slightly higher than the Propagator method (2.7 ms). As such, our fast method 2 is also applicable to real-time automotiv e sensing. D. Err or Bound vs Actual Err or W e then study the accuracy of our theoretical error bounds (Theorems 1, 2, and 3), by comparing the pseudo-spectrum ˜ P music ( θ ) approximated by our fast-MUSIC with the exact one P music ( θ ) . W e will show how much the theoretical bounds under-estimate or ov er-estimate an exact P music ( θ ) . In this simulation, we hav e M = N = 400 , p = K = 12 and t = 3 . 1) F ast MUSIC – Method 1 : W e first study the lo wer bound for our fast-MUSIC method 1, which is premised on 9 30 40 50 60 70 80 AoA /deg 10 -6 10 -4 10 -2 10 0 Spatial spectrum Fig. 6. Robustness of our fast MUSIC algorithm to inexact guess of K . Here, the actual number of targets is K = 10 . The red curve plots MUSIC with the exact K . The other four curves plot our fast MUSIC algorithm with different guesses of K . random sampling. W e denote the righthand side of the bound in Theorem 1 by κ = 1 + α l . Thus, we hav e q ˜ P music ( θ ) > p P music ( θ ) /κ. (20) W e repeat our numerical simulations to calculate a set of ˜ P music ( θ ) and P music ( θ ) . In the simulation, we assume K = 9 , M = 200 , N = 2 M and SNR=1dB. W e set p = 10 K log ( K /δ ) ( δ = 0 . 01 ) when deri ving the theoretical lower bound. Interestingly , our fast method 1 attains the accurate pseudo-spectrum, even if the user-specific parameter is p ∼ O ( K ) . In this analysis, K unknown DoAs are randomly generated in the angle range [0, 80] deg. In Figure 5(a), we plot p P music ( θ ) /κ (the y -axis) against q ˜ P music ( θ ) (the x -axis) as the blue crosses. Ideally , if the lower bound is tight ( α l → 0 , κ l → 1 ), then the blue crosses should fall on the line y = x . Since the lower bound (righthand side of eq. (20)) would ov er-estimate q ˜ P music ( θ ) , the blue crosses can fall below the line y = x with the probability approaching 1. That means, the empirical results match the deriv ed theoretical result, i.e. no blue cross is above y = x . If θ 0 falls in the non-target region ( θ 0 / ∈ θ θ θ K × 1 ), the normalized pseudo-spectrum is very small, e.g. P music ( θ 0 ) ∼ O (1 / M ) . Figure 5(a) shows that our lower bound is very tight in the non-target region. As seen, the blue crosses are almost on the line y = x and hence α l → 0 , e.g. P music ( θ 0 ) < 0 . 2 . Follo wing this, the approximated pseudo-spectrum ˜ P music ( θ 0 ) should not interpret non-target peaks as targets. If θ k is in the target region ( θ k ∈ θ θ θ K × 1 ), the normal- ized pseudo-spectrum is relatively large, P music ( θ k ) ∼ 1 . Figure 5(a) shows that in the target region, the empirical results (blue cross) deviate slightly from the lower bound, indicating our lower bound is less tight in this re gion (e.g. P music ( θ k ) ∼ 1 ). Howe ver , this would be fine for realistic applications. I.e. ˜ P music ( θ k ) for target AoAs is sufficiently large to be detected via its peak amplitude, i.e. 1 + α l  √ γ . 2) F ast MUSIC – Method 2 : W e study the theoretical error bounds of our fast method 2, as in Theorems 2 and 3. W e denote the righthand side of the lower bound in Theorem 2 by κ l = 1 + α l ; and the righthand side of the upper bound in Theorem 3 by κ u = 1 + α u . Then, we expect q ˜ P music ( θ ) ≥ p P music ( θ ) /κ l , q ˜ P music ( θ ) ≤ p P music ( θ ) /κ u . In the numerical analysis, we assume K = 9 , M = 200 , N = 2 M , δ = 0 . 1 , SNR=1dB, and t = 2 . In Figure 5(b), we show p P music ( θ ) /κ l against q ˜ P music ( θ ) as blue crosses, and p P music ( θ ) /κ u against q ˜ P music ( θ ) as green crosses. As shown, the blue crosses and green crosses lie almost on the line y = x , indicating that our theoretical lower and upper bounds are very tight, i.e. α l → 0 and α u → 0 . Combined with the theoretical results in Section IV -E, we conclude there will be no false/missed peaks in the approximated pseudo-spectrum. E. Rob ustness to Inaccurate K Similar to standard MUSIC or ESPRIT , our fast-MUSIC requires the prior kno wledge of the number of targets, K . Unfortunately , in practice K would become unkno wn. Thus, we hav e to use an over -estimate of K as the input. In Figure 6, we study the effects on the approximated pseudo-spectrum from an inexact guess of K . In this sim- ulation, we assume M = 200 , N = 400 , K = 10 , SNR = 0 dB and p = d 1 . 2 × ˜ K e . Here, the standard MUSIC is assumed with the e xact knowledge, i.e. K = 10 . As in Figure 6, our fast MUSIC method 1 is practically robust to an over -estimate of K . In other words, once the guess value ˜ K is larger than K , our fast-MUSIC produces almost the same pseudo-spectrum (e.g. with the ov erlapped peaks). Another more promising solution to combat this prior uncertainty of target numbers is to apply the pr e-estimation scheme [17], by directly acquiring a priori estimation of K . For example, both the Akaike information criterion (AIC) [57] or minimum description length (MDL) [58] could be used to estimate the number of targets K . On this basis, we may configure the sampling length p more reliably , which introduces a slight increase on time complexity b ut excludes the potential risk of underestimating K . F . Accuracy in AoA Estimation W e now ev aluate the AoA estimation accuracy and compare fast-MUSIC (e.g. Fast Method 1) with the existing subspace algorithms. In Figure 7, the normalized pseudo-spectrum is ˇ P ( θ ) = P ( θ ) − min θ { P ( θ ) } max θ { P ( θ ) }− min θ { P ( θ ) } where P ( θ ) can be the spectrum of MUSIC, ESPRIT , Propaga- tor , etc. In the numerical simulation, we assume M = 200 and K = 9 and then ev aluate the relativ e accuracy under different N . Here, we focus on our fast method 1 with p = 12 . As illustrated by Figure 7 only showing AoAs that fall into [65 86] deg), among all counterpart methods, MUSIC and block-Lanczos could attain the near-optimal pseudo-spectrum, by exploiting the exact subspace information, which thus acquire the highly accurate AoAs estimation. Our fast-MUSIC obtains an approximated pseudo-spectrum that closely tracks the standard MUSIC. Moreov er , the pseudo-spectrum of both 10 65 70 75 80 85 AoA /deg 10 -4 10 -3 10 -2 10 -1 10 0 Pseudo-spectrum MUSIC / Lanczos Inversion-based Propagator ESPRIT Fast MUSIC Ground-truth (a) M = 200 , N = 800 , and SNR = 0 dB. 65 70 75 80 85 AoA /deg 10 -4 10 -3 10 -2 10 -1 10 0 Pseudo-spectrum MUSIC / Lanczos Inversion-based Propagator ESPRIT Fast MUSIC Ground-truth (b) M = N = 200 and SNR = 0 dB. Fig. 7. Estimated pseudo-spectrum ˇ P ( θ ) of v arious algorithms. Here, the proposed algorithm refers to the fast-MUSIC method with p = 12 . MUSIC and our fast-MUSIC will not change, when the snapshot number N drops from N = 800 to N = 200 . In comparison, another ESPRIT method attains less accurate AoAs. The estimation error (i.e. false targets nearby 65 degrees) seems to be ine vitable, especially when the exact number of sources (i.e. K ) remains unkno wn. Despite its low complexity , the Propagator method produces the less accurate pseudo-spectrum. For one thing, the fluc- tuated baseline is relativ ely high, which may produce false alarms in low SNR cases. For another, it gets the less accurate estimation when the target DoAs surpass 70 deg. For example, it may misinterpret two tar gets (between 80 and 85 degrees) as one fake object, resulting in the significant estimation error . Note that, a modified Propagator method was proposed to improve the accuracy [33], which, howe ver , requires one special array structure (e.g. L -shape array) and tends to be less attractiv e for compact automotive MIMO radars. For another matrix-in version method, ev en through the time complexity in approximating noise subspace is reduced to some extent, it still has drawbacks in realistic applications. Its estimation accuracy may be guaranteed only in high SNRs, subject to the sufficiently accurate approximation of subspace. In the case of N ≤ M or strong noise ( σ K ( S ) ' σ 2 n ), it may be unstable and produces false objects or inaccurate AoAs. Except for the linear complexity and real-time implementa- tion, our fast-MUSIC is highly accurate and stable. Unlike its counterparts that hav e to scarify accuracy for complexity , our -8 -6 -4 -2 0 2 SNR /dB 10 -3 10 -2 10 -1 AoA Estimation MSE Inversion-based Propagator ESPRIT Fast MUSIC MUSIC / Lanczos Fig. 8. Plot of AoA estimation MSE against SNRs. new method resolves the long-standing dilemma between high resolution and low complexity in MIMO radar signal process- ing. It thus achiev es the high-resolution AoA estimation at a scalable complexity , which provides also the great potential to massiv e MIMO signal processing. One practical challenge for our fast-MUSIC is that, when multiple targets of dif ferent SNRs are taken into accounts, the detection accuracy may be af fected slightly . From Figure 3, we observe the targets with low SNRs may be probably missed, when the number of elements is relativ ely small. Fortunately , owing to the substantial processing gain of massiv e MIMO array and the subspace methods, this baseline fluctuation is extremely lo w , i.e. O (1 / M ) of the maximum peak. In other words, those lo w SNR targets can be still detected and their AoAs can be estimated, provided the smallest SNR was larger than O (1 / M ) of the maximum SNR. G. MSE in AoA Estimation Finnaly , we e valuate the AoA estimation error of various methods in the conte xt of massiv e-MIMO radar . The mean square error (MSE) of estimated AoAs is measured by: MSE , E θ k ∈U (0 ,π / 2)  1 K X K k =1 k θ k − ˆ θ k k 2 2  . (21) Here, θ k denotes the ground truth AoA of the k th target point, while ˆ θ k denoted its estimation. In Figure 8, we show the MSE performance of various subspace methods. W e set M = 200 , K = 10 and N = 220 ; our fast method 1 is used with p = 11 . The AoAs estimation MSE of all subspace methods would be reduced, as SNR increases. As noted, SNR affects the spectral gap σ K +1 ( S ) σ K ( S ) , and hence the approximation accuracy of the MUSIC pseudo-spectrum. Ho wev er , it is note worthy that the MSE of estimated AoA has a non-linear relation with the approximated pseudo-spectrum. I.e. only when the error of estimated pseudo-spectrum surpasses one threshold, can the attained AoA suffer a serious deviation [41], [59]. In this regard, our fast-MUSIC attains the same MSE as standard MUSIC 4 , e ven with the pseudo-spectrum approximation error . That is to say , it scarifies no accuracy in estimated AoA, while substantially reducing the time complexity . 4 Note that, when the number of elements is relatively small (e.g. M ≤ 100), the MSE of our fast-MUSIC will be still comparable to standard MUSIC. 11 V I . C O N C L U S I O N W e inv estigate the high-resolution target detection and AoA estimation problem in the automotiv e massive-MIMO radar . Existing subspace methods (e.g. MUSIC and ESPRIT) are computationally expensi v e, due to the complex decomposition of a large co variance matrix. Other simplified schemes, e.g. matrix-in verse and propagator methods, lead to the degraded accuracy . As one long-standing dilemma in MIMO radar signal processing, the high-resolution estimation and the real-time computation remain mutually exclusi ve. T o enable real-time and high-resolution en vironmental sens- ing, we le verage on randomized matrix approximation and dev elop two fast-MUSIC methods, whereby a large cov ariance was approximated by three small sketches that are abstracted via random sampling/projection. W e theoretically show that, in high SNRs, our fast-MUSIC is almost as good as the near- optimal MUSIC. The numerical study also demonstrates that, whilst our method is nearly as accurate as MUSIC in acquiring unknown AoA, it is dramatically faster than standard MUSIC by orders-of-magnitude. As such, our fast-MUSIC enjoys both the linear complexity and the high accuracy , breaking the theoretical bottleneck in MIMO radar signal processing. It would provide the great potential to high-resolution and real- time massi ve-MIMO radars in the emer ging automotiv e appli- cations. In the future study , we may consider the extension of this new concept of randomized matrix sketching to the other methods, e.g. subspace methods (e.g. Capton method) and the optimal ML method. A C K N O W L E D G M E N T This work was supported by Natural Science Foundation of China (NSFC) under Grants U1805262. A P P E N D I X Here, we show the detailed proofs of Theorems 1, 2, and 3. For simplicity , we leav e out the θ in P music ( θ ) , ˜ P music ( θ ) , and a ( θ ) . A. Analysis of Uniform Sampling An M × p ( p < M ) matrix Π is called a uniform sampling matrix if its columns are sampled from the columns of √ M √ p I M uniformly at random. For any N × M matrix X , the multiplication XΠ ( N × p ) contains p uniformly sampled columns of X . Lemma 1. Let Π be an M × p uniform sampling matrix. Then   Π   2 2 = M /p. Lemma 2 ( [49], [56]) . Let U K be an M × K matrix with orthonormal columns. Let the r ow coherence of U K be defined as µ ( U K ) = M K max i   ( U K ) i :   2 2 ∈  1 , M K  . Let Π be an M × p uniform sampling matrix. F or p ≥ (6 + 2 η ) µ ( U K ) K 3 η 2 log K δ the K singular values of U H K ΠΠ H U K ar e at least 1 − η with pr obability at least 1 − δ . B. Analysis of Gaussian Pr ojection Lemma 3 ( [60]) . Let Π be an M × K matrix whose entries ar e i.i.d. drawn fr om N (0 , 1) . If M is substantially lar ger than K , the spectral norm of Π is bounded by   Π   2 ≤ √ M + √ K + O (1) , almost sur ely . Lemma 4 ( [61], [62]) . Let G be a K × K matrix whose entries ar e i.i.d. drawn fr om N (0 , 1) . Then the smallest singualr value of G satisfies σ K ( G ) ≥ δ √ K with pr obability at least 1 − δ − o (1) . C. Pr oof of Theor em 1 Lemma 5. Let the notation be defined in Section IV -B. Let Π be any matrix with rank( U H K Π ) ≥ K . Then, for all a ,   U K U H K a − ˜ U ˜ U H U K U H K a   2 ≤    ( S − S K ) Π  U H K Π  † Σ − 1 K U H K a    2 , Pr oof. Let ˜ U be the orthonormal bases of ˜ S . Lemma 12 of [56] shows that  CW † C H  CW † C H  † C = C . Thus C = SΠ is in the subspace spanned by the columns of ˜ U , and therefore, ˜ U ˜ U H SΠ = SΠ . (22) It is not hard to prove that ˜ U H U K U H K a = argmin x   U K U H K a − ˜ Ux   2 . Thus, the inequality holds for all ˜ x :   U K U H K a − ˜ U ˜ U H U K U H K a   2 = min x   U K U H K a − ˜ Ux   2 ≤   U K U H K a − ˜ U ˜ x   2 . (23) W e artifically construct ˜ x = ˜ U H SΠ  U H K Π  † Σ − 1 K U H K a . It follows from (22) that ˜ U ˜ x = ˜ U ˜ U H SΠ  U H K Π  † Σ − 1 K U H K a = SΠ  U H K Π  † Σ − 1 K U H K a =  S K + S − S K  Π  U H K Π  † Σ − 1 K U H K a . (24) Since rank( U H K Π ) ≥ K , we hav e ( U H K Π )( U H K Π ) † = I K . Thus S K Π  U H K Π  † Σ − 1 K U K a = U K Σ K U H K Π  U H K Π  † Σ − 1 K U H K a = U K Σ K Σ − 1 K U H K a = U K U H K a . (25) It follows from (24) and (25) that ˜ U ˜ x =  S K + ( S − S K )  Π  U H K Π  † Σ − 1 K U H K a = U K U H K a + ( S − S K ) Π  U H K Π  † Σ − 1 K U H K a . (26) 12 It follows from (23) and (26) that   U K U H K a − ˜ U ˜ U H U K U H K a   2 ≤   U K U H K a − ˜ U ˜ x   2 =    ( S − S K ) Π  U H K Π  † Σ − 1 K U H K a    2 , by which the lemma follows. Complete the proof of Theorem 1: Pr oof. The M × p matrix Π is a uniform sampling matrix. Then C = SΠ and W = Π T SΠ . Lemma 1 sho ws the spectral norm of Π is bounded by   Π   2 ≤ p M /p. Because U K has orthonormal columns, Lemma 2 shows that for p ≥ (6+2 η ) µ ( U K ) K 3 η 2 log K δ , the K -th singular v alue of U H K Π is bounded by σ K ( U H K Π ) ≥ p 1 − η with probability at least 1 − δ . Let follows from Lemma 5 that    I M − ˜ U ˜ U H  U K U H K a   2 ≤   U K U H K a − ˜ U ˜ U H U K U H K a   2 ≤    ( S − S K ) Π  U H K Π  † Σ − 1 K U H K a    2 ≤   ( S − S K )   2   Π   2    U H K Π  †   2   Σ − 1 K   2   U H K a   2 ≤ σ K +1 ( S ) σ K ( S ) p M /p √ 1 − η   U H K a   2 . W e set η = 0 . 75 . Then, for p ≥ 4 . 5 µ ( U K ) K · log K δ , it holds with probability at least 1 − δ that    I M − ˜ U ˜ U H  U K U H K a   2 ≤ 2 q M p σ K +1 ( S ) σ K ( S )   U H K a   2 . It follows that    I M − ˜ U ˜ U H  a   2 =    I M − ˜ U ˜ U H  U K U H K + U − K U H − K  a   2 ≤    I M − ˜ U ˜ U H  U K U H K a   2 +   U − K U H − K a   2 . Hence,    I M − ˜ U ˜ U H  a   2 −    I M − U K U H K  a   2 ≤    I M − ˜ U ˜ U H  U K U H K a   2 ≤ 2 q M p σ K +1 ( S ) σ K ( S )   U H K a   2 , where the last inequality holds with probability at least 1 − δ . W e leave out θ in P music ( θ ) , ˜ P music ( θ ) , and a ( θ ) . It follows from (5) and (10) that q P music ˜ P music =   ( I − ˜ U ˜ U H ) a   2   ( I − U K U H K ) a   2 ≤ 1 + 2 q M p σ K +1 ( S ) σ K ( S )   U H K a   2   ( I − U K U H K ) a   2 . Furthermore, when the number of element is suf ficiently large, for an variant pseudo-spectrum spectrum || U H K a || 2 and   ( I − U K U H K ) a   2 , the following relations always hold: || U H K a || 2 2 ≤ || U H K a ( θ k ) || 2 2 → M , (27) and   ( I − U K U H K ) a   2 ≥   ( I − U K U H K ) a ( θ k )   2 → 1 , (28) where θ k denotes the k -th AoA of target point. As a result, we further hav e:   U H K a   2   ( I − U K U H K ) a   2 ≤ || U H K a || 2 || a || 2 ≤ √ M . (29) On this basis, we can finally obtain the Theorem 1: q P music ˜ P music ≤ 1 + 2 q M 2 p σ K +1 ( S ) σ K ( S ) . D. Pr oof of Theor em 2 Lemma 6. Let the notation be defined in Section IV -C. Let Π be any matrix satisfying that rank( U H K Π ) ≥ K . Then, for all a ,   U K U H K a − ˜ U ˜ U H U K U H K a   2 ≤    ( S − S K ) t +1 Π  U H K Π  † Σ − t − 1 K U H K a    2 . Pr oof. Since ˜ S = CW † C H , it follows from [56, Lemma 12] that ˜ S ˜ S † C = C . Since ˜ U is an orthonormal bases of ˜ S , ˜ U ˜ U H C = ˜ S ˜ S † C = C . By the definition C = SV and that V is the orthonormal basis of S t Π , we ha ve that C and S t +1 Π have the same column space. ˜ U ˜ U H S t +1 Π = S t +1 Π . (30) In the same way as the proof of Lemma 5, we can prove that for all ˜ x :   U K U H K a − ˜ U ˜ U H U K U H K a   2 ≤   U K U H K a − ˜ U ˜ x   2 . (31) W e artificially construct ˜ x = ˜ U H S t +1 Π  U H K Π  † Σ − t − 1 K U H K a . It follows from (30) that ˜ U ˜ x = ˜ U ˜ U H S t +1 Π  U H K Π  † Σ − t − 1 K U H K a =  S t +1 K + S t +1 − S t +1 K  Π  U H K Π  † Σ − t − 1 K U H K a . (32) In the same way as the proof of Lemma 5, we can use rank( U H K Π ) ≥ K to show that S t +1 K Π  U H K Π  † Σ − t − 1 K U H K a = U K Σ t +1 K U H K Π  U H K Π  † Σ − t − 1 K U H K a = U K Σ t +1 K Σ − t − 1 K U H K a = U K U H K a . 13 It follows from (32) that ˜ U ˜ x = U K U H K a +  S t +1 − S t +1 K  Π  U H K Π  † Σ − t − 1 K U H K a . It follows from (31) that   U K U H K a − ˜ U ˜ U H U K U H K a   2 ≤     S t +1 − S t +1 K  Π  U H K Π  † Σ − t − 1 K U H K a    2 , by which the lemma follows. Complete the proof of Theorem 2: Pr oof. Lemma 3 shows that the spectral norm of the M × K standard Gaussian matrix Π satisfies   Π   2 ≤ √ M + √ K + O (1) , almost surely . Since U K has orthonormal columns, the K × K matrix U H K Π is a standard Gaussian matrix. It follows from Lemma 4 that σ − 1 K ( U H K Π ) ≤ √ K δ with probability at least 1 − δ − o (1) . It follows from Lemma 6 that   U K U H K a − ˜ U ˜ U H U K U H K a   2 ≤     S t +1 − S t +1 K  Π  U H K Π  † Σ − t − 1 K U H K a    2 ≤  σ K +1 ( S ) σ K ( S )  t +1   Π   2   ( U H K Π ) †   2   U H K a   2 ≤  σ K +1 ( S ) σ K ( S )  t +1 √ M K + K + O ( √ K ) δ   U H K a   2 . The proof follows from the proof of Theorem 1 that q P music ˜ P music =   ( I − ˜ U ˜ U H ) a   2   ( I − U K U H K ) a   2 ≤ 1 +    I M − ˜ U ˜ U H  U K U H K a   2   ( I − U K U H K ) a   2 . Thus, it holds with probability at least 1 − δ − o (1) that q P music ˜ P music ≤ 1 +  σ K +1 ( S ) σ K ( S )  t +1 √ M K  1+ o (1)  δ   U H K a   2   ( I − U K U H K ) a   2 , by incorporating eq. (29), then the theorem can be prov ed, i.e. q P music ˜ P music ≤ 1 +  σ K +1 ( S ) σ K ( S )  t +1 √ M K  1+ o (1)  δ √ M , = 1 + √ M 2 K δ  σ K +1 ( S ) σ K ( S )  t +1 . Proof of Theorem 3: Pr oof. It can be shown that for all ˜ X :   U K U H K − ˜ U ˜ U H U K U H K   2 ≤   U K U H K − ˜ U ˜ X   2 . W e artifically construct ˜ X = ˜ U H S t +1 Π  U H K Π  † Σ − t − 1 K U H K . In the same way as the proof of Lemma 6, we can show that ˜ U ˜ X = U K U H K +  S t +1 − S t +1 K  Π  U H K Π  † Σ − t − 1 K U H K . Thus   U K U H K − ˜ U ˜ U H U K U H K   2 ≤    S t +1 − S t +1 K  Π  U H K Π  † Σ − t − 1 K U H K   2 ≤  σ K +1 ( S ) σ K ( S )  t +1   Π   2   ( U H K Π ) †   2 . It follows from the proof of Theorem 2 that   U K U H K − ˜ U ˜ U H U K U H K   2 ≤  1 + o (1)  √ M K δ  σ K +1 ( S ) σ K ( S )  t +1 holds with probability at least 1 − δ . Since U K and ˜ U are both M × K , [63, Eqn 2.54] shows that   ( I M − ˜ U ˜ U H ) U K   2 =   ( I M − U K U H K ) ˜ U   2 . It follows that with probability at least 1 − δ ,   ( I M − U K U H K ) ˜ U   2 =   ( I M − ˜ U ˜ U H ) U K   2 ≤  1 + o (1)  √ M K δ  σ K +1 ( S ) σ K ( S )  t +1 . (33) W e hav e that   ( I M − U K U H K ) a   2 ≤   ( I M − U K U H K ) ˜ U ˜ U H a   2 +   ( I M − U K U H K )( I M − ˜ U ˜ U H ) a   2 ≤   ( I M − U K U H K ) ˜ U ˜ U H a   2 +   ( I M − ˜ U ˜ U H ) a   2 . (34) It follows from (33) and (34) that   ( I M − U K U H K ) a   2 −   ( I M − ˜ U ˜ U H ) a   2 ≤   ( I M − U K U H K ) ˜ U   2   a   2 ≤  1 + o (1)  √ M K δ  σ K +1 ( S ) σ K ( S )  t +1   a   2 , where the latter inequality holds with probability at least 1 − δ . Equiv alently ,   ( I M − ˜ U ˜ U H ) a   2 −   ( I M − U K U H K ) a   2 ≥ −  1 + o (1)  √ M K δ  σ K +1 ( S ) σ K ( S )  t +1   a   2 , Thus, with probability at least 1 − δ , q P music ˜ P music =   ( I − ˜ U ˜ U H ) a   2   ( I − U K U H K ) a   2 ≥ 1 −  1 + o (1)  √ M K δ  σ K +1 ( S ) σ K ( S )  t +1 k a k 2 k ( I − U K U H K ) a k 2 . By further applying the inequality in eq. (28), we hav e k a k 2 k ( I − U K U H K ) a k 2 ≤ √ M . by which the theorem can be proved, i.e. q P music ˜ P music =   ( I − ˜ U ˜ U H ) a   2   ( I − U K U H K ) a   2 ≥ 1 − √ M 2 K δ  σ K +1 ( S ) σ K ( S )  t +1 . R E F E R E N C E S [1] W . D. Jones, “Keeping cars from crashing, ” IEEE Spectrum , vol. 38, no. 9, pp. 40–45, 2001. [2] D. Guizzo, “Ho w google’ s self-driving car works, ” IEEE Spectrum Online , vol. 18, 2011. 14 [3] E. P . Blasch, A. Lakhotia, and G. Seetharaman, “Unmanned vehicles come of age: The darpa grand challenge, ” Computer , vol. 39, no. 12, pp. 26–29, 2006. [4] K. Alonzo, A. Stentz, O. Amidi, M. Bode, D. Bradley , A. Diazcalderon, M. Happold, H. Herman, R. Mandelbaum, and T . Pilarski, “T oward reliable off road autonomous vehicles operating in challenging environ- ments, ” in Pr oc of International Symposium on Exprimental Robotics , 2004. [5] M. Murad, I. Bilik, M. Friesen, J. Nickolaou, J. Salinger, K. Geary , and J. S. Colburn, “Requirements for next generation automotive radars, ” in Pr oc. of 2013 IEEE Radar Conference , pp. 1–6, 2013. [6] C. W aldschmidt and H. H. Meinel, “Future trends and directions in radar concerning the application for autonomous driving, ” in Proc. of 2014 European Micr owave Conference , pp. 1719–1722, 2014. [7] J. Hasch, E. T opak, R. Schnabel, T . Zwick, R. W eigel, and C. W ald- schmidt, “Millimeter-wa ve technology for automotive radar sensors in the 77 ghz frequency band, ” IEEE Tr ansactions on Micr owave Theory and T ec hniques , vol. 60, no. 3, pp. 845–860, 2012. [8] J. Li and P . Stoica, “Mimo radar with colocated antennas, ” IEEE Signal Pr ocessing Magazine , vol. 24, no. 5, pp. 106–114, 2007. [9] I. Bilik, S. V illev al, D. Brodeski, H. Ringel, O. Longman, P . Goswami, C. Y . B. Kumar , S. Rao, P . Swami, A. Jain, et al. , “ Automotive multi- mode cascaded radar data processing embedded system, ” in Proceedings of 2018 IEEE Radar Confer ence , pp. 1–6, 2018. [10] I. Bilik, O. Longman, S. V illev al, and J. T abrikian, “The Rise of Radar for Autonomous V ehicles: Signal Processing Solutions and Future Research Directions, ” IEEE Signal Pr ocessing Magazine , vol. 36, no. 5, pp. 20–31, 2019. [11] S. P atole, M. T orlak, D. W ang, and M. Ali, “ Automoti ve radars: A review of signal processing techniques, ” IEEE Signal Pr ocessing Magazine , vol. 34, no. 2, pp. 22–35, 2017. [12] G. Hakobyan and B. Y ang, “High-Performance Automotive Radar: A Revie w of Signal Processing Algorithms and Modulation Schemes, ” IEEE Signal Pr ocessing Magazine , vol. 36, no. 5, pp. 32–44, 2019. [13] A. G. Stove, “Linear FMCW radar techniques, ” Pr oceedings of the IEEE , 1992. [14] C. H. Lin, Y . S. W u, Y . L. Y eh, S. H. W eng, G. Y . Chen, C. H. Shen, and H. Y . Chang, “A 24-GHz highly integrated transceiv er in 0.5- µ m E/D-PHEMT process for FMCW automotive radar applications, ” in Micr owave Conference Pr oceedings , pp. 512–515, 2011. [15] G. Babur , O. A. Krasnov , A. Y aro voy , and P . Aubry , “Nearly Orthogonal W aveforms for MIMO FMCW Radar, ” IEEE Tr ansactions on Aerospace & Electronic Systems , vol. 49, no. 3, pp. 1426–1437, 2013. [16] D. Kok and J. S. Fu, “Signal processing for automotiv e radar , ” in Radar Confer ence, 2005 IEEE International , pp. 842–846, 2005. [17] S. Sun, A. P . Petropulu, and H. V . Poor , “MIMO Radar for Advanced Driv er-Assistance Systems and Autonomous Driving: Advantages and Challenges, ” IEEE Signal Processing Magazine , vol. 37, no. 4, pp. 98– 117, 2020. [18] J. W enger, “ Automoti ve radar - status and perspectiv es, ” in Compound Semiconductor Integr ated Circuit Symposium, 2005 (CSIC ’05) , p. 4 pp., 2005. [19] R. O. Schmidt, “Multiple emitter location and signal parameter estima- tion, ” IEEE T ransactions on Antennas and Propa gation , vol. 34, no. 3, pp. 276–280, 1986. [20] R. H. Roy and T . Kailath, “ESPRIT-estimation of signal parameters via rotational inv ariance techniques, ” IEEE Tr ansactions on Acoustics, Speech, and Signal Pr ocessing , vol. 37, no. 7, pp. 984–995, 1989. [21] H. Krim and M. V iberg, “T wo decades of array signal processing research: the parametric approach, ” IEEE Signal Processing Magazine , vol. 13, no. 4, pp. 67–94, 1996. [22] D. Oh and J. Lee, “Low-complexity range-azimuth fmcw radar sensor using joint angle and delay estimation without svd and evd, ” IEEE Sensors Journal , vol. 15, no. 9, pp. 4799–4811, 2015. [23] M. D. Zolto wski, G. M. Kautz, and S. D. Silverstein, “Beamspace Root-MUSIC, ” IEEE T ransactions on Signal Pr ocessing , vol. 41, no. 1, pp. 344–364, 1993. [24] H. D. Simon, “The lanczos algorithm with partial reorthogonalization, ” Mathematics of Computation , vol. 42, no. 165, pp. 115–142, 1984. [25] Mathe ws, P . C., Zoltowski, and M. D ., “Eigenstructure techniques for 2-D angle estimation with uniform circular arrays, ” Signal Pr ocessing, IEEE T ransactions on , vol. 42, no. 9, pp. 2395–2407, 1994. [26] Z. Lin, T . Lv , and P . T . Mathiopoulos, “3-D Indoor Positioning for Millimeter-W ave Massi ve MIMO Systems, ” IEEE Tr ansactions on Com- munications , vol. 66, no. 6, pp. 2472–2486, 2018. [27] M. Rossi, A. M. Haimovich, and Y . C. Eldar , “Spatial Compressive Sensing for MIMO Radar, ” IEEE T ransactions on Signal Pr ocessing , vol. 62, no. 2, pp. 419–430, 2013. [28] Stoeckle, Christoph, Munir, Jawad, Mezghani, Nossek, Munich, and Germany), “DoA Estimation Performance and Computational Complex- ity of Subspace- and Compressed Sensing-based Methods, ” in in Proc. of 19th International ITG W orkshop on Smart Antennas, , 2015. [29] B. Li and A. P . Petropulu, “Distributed MIMO radar based on sparse sensing: Analysis and efficient implementation, ” IEEE T ransactions on Aer ospace & Electr onic Systems , vol. 51, no. 4, pp. 3055–3070, 2015. [30] J. Benesty , J. Chen, and Y . Huang, “ A generalized MVDR spectrum, ” IEEE Signal Pr ocessing Letters , vol. 12, no. 12, pp. 827–830, 2005. [31] S. Marcos, A. Marsal, and M. Benidir, “The propagator method for source bearing estimation, ” Signal Processing , vol. 42, no. 2, pp. 121– 138, 1995. [32] S. Marcos, A. Marsal, and M. Benidir, “Performances analysis of the propagator method for source bearing estimation, ” in Acoustics, Speech, & Signal Pr ocessing, on IEEE International Confer ence , 1994. [33] N. T ayem and H. M. Kwon, “L-shape 2-dimensional arriv al angle estimation with propagator method, ” IEEE T rans Antennas & Pr opagat , vol. 1, no. 5, pp. 1622–1630, 2005. [34] F . Engels, P . Heidenreich, A. M. Zoubir, F . K. Jondral, and M. W in- termantel, “ Advances in automotive radar: A framework on compu- tationally efficient high-resolution frequency estimation, ” IEEE Signal Pr ocessing Magazine , vol. 34, no. 2, pp. 36–46, 2017. [35] M. Y . Keegan Garcia and A. Purkovic, “Robust traffic and intersection monitoring using millimeter wave sensors, ” tech. rep., T exas Instruments Incorporated (TI), May 2018. [36] K. Ramasubramanian, “Using a complex-baseband architecture in FMCW radar systems, ” tech. rep., T exas Instruments Incorporated (TI), May 2017. [37] J. He, M. N. S. Swamy , and M. O. Ahmad, “Ef ficient Application of MUSIC Algorithm Under the Coexistence of Far -Field and Near- Field Sources, ” IEEE T ransactions on Signal Processing , vol. 60, no. 4, pp. 2066–2070, 2012. [38] B. Li, S. W ang, J. Zhang, X. Cao, and C. Zhao, “Fast Randomized- MUSIC for mm-W av e Massiv e MIMO Radars, ” IEEE T ransactions on V ehicular T echnology , vol. 70, no. 2, pp. 1952–1956. [39] P . Stoica and N. Arye, “MUSIC, maximum likelihood, and Cramer-Rao bound, ” IEEE T ransactions on Acoustics, Speech, and Signal Pr ocessing , vol. 37, no. 5, pp. 720–741, 1989. [40] S. Fortunati, L. Sanguinetti, M. S. Greco, and F . Gini, “Scaling up MIMO radar for target detection, ” in 2019 IEEE International Con- fer ence on Acoustics, Speech and Signal Pr ocessing (ICASSP 2019) , 2019. [41] P . V allet, X. Mestre, and P . Loubaton, “Performance Analysis of an Improved MUSIC DoA Estimator, ” IEEE Tr ansactions on Signal Pro- cessing , vol. 63, no. 23, pp. 6407–6422, 2015. [42] W . Liao and A. Fannjiang, “MUSIC for Single-Snapshot Spectral Estimation: Stability and Super-resolution, ” Applied & Computational Harmonic Analysis , vol. 40, no. 1, pp. 33–67, 2014. [43] B. Li, S. W ang, Z. Feng, J. Zhang, X. Cao, and C. Zhao, “Fast Pseudo- spectrum Estimation for Automotiv e Massi ve MIMO Radar, ” IEEE Internet of Things Journal , vol. 8, no. 20, pp. 15303–15316, 2021. [44] H. Cao, Q. Liu, and Y . Wu, “Transform Domain: Design of Closed- Form Joint 2-D DOA Estimation Based on QR Decomposition, ” Circuits Systems and Signal Processing , pp. 1–12, 2020. [45] Q. Liu, Y . Gu, and H. C. So, “DO A Estimation in Impulsiv e Noise via Low-Rank Matrix Approximation and W eakly Con vex Optimization, ” IEEE Tr ansactions on Aer ospace and Electronic Systems , vol. 55, no. 6, 2019. [46] Q. Liu, J. Xu, Z. Ding, and H. C. So, “T ar get Localization With Jammer Remov al Using Frequency Diverse Array , ” IEEE T ransactions on V ehicular T echnology , vol. 69, no. 10, pp. 11685–11696, 2020. [47] S. W ang, A. Gittens, and M. W . Mahoney , “Scalable kernel k-means clustering with Nystrom approximation: Relative-error bounds, ” Journal of Machine Learning Researc h , vol. 20, no. 12, pp. 1–49, 2019. [48] P . Drineas and M. W . Mahoney , “On the Nystr ¨ om method for approx- imating a gram matrix for improved kernel-based learning, ” J ournal of Machine Learning Research , vol. 6, pp. 2153–2175, 2005. [49] D. P . W oodruff, “Sketching as a tool for numerical linear algebra, ” F oundations and T r ends® in Theoretical Computer Science , vol. 10, no. 1–2, pp. 1–157, 2014. [50] B. Li, P . Chen, H. Liu, W . S. Guo, X. B. Cao, J. Z. Du, C. L. Zhao, and J. Zhang, “Random Sketch Learning for Deep Neural Networks in Edge Computing, ” Nature Computational Science , vol. 1, no. 3, pp. 221–228, 2021. 15 [51] B. Li, S. W ang, J. Zhang, X. Cao, and C. Zhao, “Randomized Ap- proximate Channel Estimator in Massi ve-MIMO Communication, ” IEEE Communications Letters , vol. 24, no. 10, pp. 2314–2318, 2020. [52] J. A. T ropp, A. Y urtsever , M. Udell, and V . Cevher , “Fixed-rank approximation of a positive-semidefinite matrix from streaming data, ” in Advances in Neural Information Pr ocessing Systems (NIPS) , 2017. [53] A. Gittens and M. W . Mahoney , “Revisiting the Nystrom method for improved large-scale machine learning, ” Journal of Machine Learning Resear ch , vol. 17, no. 1, pp. 3977–4041, 2016. [54] B. M. Sadler, G. B. Giannakis, and S. Shamsunder , “Noise subspace techniques in non-gaussian noise using cumulants, ” IEEE T ransactions on Aer ospace & Electr onic Systems , v ol. 31, no. 3, pp. 1009–1018, 1995. [55] E. J. Candes and B. Recht, “Exact matrix completion via con vex optimization, ” F oundations of Computational mathematics , vol. 9, no. 6, p. 717, 2009. [56] S. W ang, L. Luo, and Z. Zhang, “SPSD matrix approximation vis column selection: Theories, algorithms, and extensions, ” J ournal of Machine Learning Research , vol. 17, no. 49, pp. 1–49, 2016. [57] H. Akaike, “ A ne w look at statistical model identification, ” IEEE T ransactions on Automatic Contr ol , no. 6, pp. 716–723, 1974. [58] A. Barron, J. Rissanen, and B. Y u, “The minimum description length principle in coding and modeling, ” IEEE T ransactions on Information Theory , vol. 44, no. 6, pp. 2743–2760, 1998. [59] J. K. Thomas and L. L. Scharf, “The probability of a subspace swap in the SVD, ” IEEE T ransactions on Signal Processing , vol. 43, no. 3, pp. 730–736, 1995. [60] R. V ershynin, “Introduction to the non-asymptotic analysis of random matrices, ” arXiv preprint , 2010. [61] M. Rudelson and R. V ershynin, “The Littlew ood–Offord problem and in vertibility of random matrices, ” Advances in Mathematics , vol. 218, no. 2, pp. 600–633, 2008. [62] T . T ao and V . V u, “Random matrices: The distribution of the smallest singular v alues, ” Geometric And Functional Analysis , vol. 20, no. 1, pp. 260–297, 2010. [63] P . Arbenz, “Lecture notes on solving large scale eigenv alue problems, ” D-MA TH, EHT Zurich , vol. 2, 2012.