{'fraction_non_alphanumeric': 0.05379058664733671, 'fraction_numerical': 0.01779688554502962, 'mean_word_length': 4.16298877274487, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Leveraging the Restricted Isometry Property: Improved Low-Rank Subspace Decomposition for Hybrid Millimeter-Wave Systems

Communication at millimeter wave frequencies will be one of the essential new technologies in 5G. Acquiring an accurate channel estimate is the key to facilitate advanced millimeter wave hybrid multiple-input multiple-output (MIMO) precoding techniques. Millimeter wave MIMO channel estimation, however, suffers from a considerably increased channel use overhead. This happens due to the limited number of radio frequency (RF) chains that prevent the digital baseband from directly accessing the signal at each antenna. To address this issue, recent research has focused on adaptive closed-loop and two-way channel estimation techniques. In this paper, unlike the prior approaches, we study a non-adaptive, hence rather simple, open-loop millimeter wave MIMO channel estimation technique. We present a random phase rotation design of channel subspace sampling signals and show that they obey the restricted isometry property (RIP) with high probability. We then formulate the channel estimation as a low-rank subspace decomposition problem and, based on the RIP, show that the proposed framework reveals resilience to a low signal-to-noise ratio. It is revealed that the required number of channel uses ensuring a bounded estimation error is linearly proportional to the degrees of freedom of the channel, whereas it converges to a constant value if the number of RF chains can grow proportionally to the channel dimension while keeping the channel rank fixed. In particular, we show that the tighter the RIP characterization the lower the channel estimation error is. We also devise an iterative technique that effectively finds a suboptimal, but stationary, solution to the formulated problem. The proposed technique is shown to have improved channel estimation accuracy with a substantially low channel use overhead as compared to that of previous closed-loop and two-way adaptation techniques.


I. INTRODUCTION
With the growth in the demand for wireless broadband, the operating frequency of modern wireless systems is steadily shifting upward to the millimeter wave band to provide a much wider bandwidth [2]. In the millimeter wave bands, the radio channel experiences severe pathloss that is compensated for by using large-scale antenna arrays [3]- [5]. To make the system available at low cost, large-sized analog arrays that are fed by a limited number of radio frequency (RF) chains are popularly discussed. This system is often referred to as a millimeter wave hybrid multiple-input multiple-output (MIMO) system. Despite the limited RF chains, a remarkable throughput boost is possible if the system employs advanced hybrid MIMO precoding techniques [5]- [9]. To enable such techniques, accurate channel state information (CSI) is the key.
Useful CSI has to be attained via channel estimation. However, millimeter wave hybrid MIMO channel estimation is limited by an insufficient number of RF chains. This limitation increases the channel use overhead required to extract a useful channel estimate of the system. To address this difficulty, recent research has focused on adaptive channel estimation frameworks, including adaptive closed-loop [5], [9]- [12] and two-way [8] channel estimation methods. Closed-loop channel estimation includes techniques such as the multi-level adaptive subspace search techniques [5], [9]- [12], where the receiver at each search level conveys a feedback to the transmitter to guide the next level subspace sampling signal design [13]- [16]. The two-way technique [8] exploits the uplink-downlink channel reciprocity to iteratively refine the channel estimate via uplink-downlink channel echoing [17]- [20]. Improved estimation accuracy compared to other non-adaptive techniques, e.g., open-loop techniques, has been evidenced in [5], [8]- [12].
Nevertheless, these processes still require a relatively large number of channel uses because of the repeated channel sampling [5], [8]- [12]. Achieving lower channel use overhead is particularly important to the millimeter wave communications that are limited by a shortage of channel coherence resources.
Open-loop MIMO channel estimation has been thoroughly studied in the literature [21], [22]. It is nonadaptive and simple. However, it is impractical for estimating large-dimensional channels, such as sub-6 GHz massive MIMO channels [13]- [15], due to its inaccuracy. Meanwhile, there have been convincing experimental studies showing that the practical millimeter wave MIMO channel exhibits substantial rank sparsity [23], [24]. A general belief is that the sparsity of a signal, either in its support or subspace, provides tremendous clues when estimating it from distorted observations. Useful directions for open-loop millimeter wave MIMO channel estimation are recently suggested by the low-rank matrix reconstruction techniques [25]- [33].
Methods for low-rank matrix reconstruction have recently spurred considerable interest due to their exploitation of rank sparsity [26]- [33]. The reliability of low-rank matrix reconstruction is best described by the restricted isometry property (RIP) [26], [31]. The RIP models, in probability, the convergence of the subspace-sampled matrix to its original matrix. It has also been used to analyze the estimation error performances of various algorithms [29], [31], [34]. The RIP concept is useful and can be exploited in devising new design criteria for the open-loop channel estimation. However, low-rank matrix reconstruction techniques [26]- [33], in their current form, face difficulties in estimating the millimeter wave MIMO channels. Reliable performance of low-rank matrix reconstruction techniques [26]- [33] is only ensured in the no noise or asymptotically high signal-to-noise ratio (SNR) regime. Reliable operation at low SNR is critical for millimeter wave systems that are limited by heavy mixed signal processing with an excessive power consumption [4]. Moreover, the probabilistic analysis of when the RIP holds has been studied for the subspace samples with independently and identically distributed (i.i.d.) Gaussian or Bernoulli entries [26], [31]. Yet, the subspace samples that the millimeter wave arrays generate are not well modeled by these distributions. Thus far, its RIP characterization has not been reported in the literature.
In this paper, we study an open-loop channel estimation technique by leveraging low-rank matrix reconstruction for millimeter wave hybrid MIMO systems. This is in contrast with prior approaches like [5], [8]- [12], that focus on closed-loop and two-way techniques. We begin by introducing a random subspace sampling signal design, which does not rely on any prior knowledge of channel subspace information. We adopt the precoder and combiner gain maximization criterion [35] to formulate the channel estimation as a low-rank subspace decomposition problem. We show through an analysis that the random ensembles of the subspace sampling signals satisfy the RIP with high probability. On the basis of the established RIP, we carry out channel estimation error analysis and show that the proposed low-rank subspace decomposition reveals resilience to low SNRs.
In particular, we reveal that a tighter RIP bound characterization results in a lower estimation error for the devised low-rank subspace decomposition problem. It is also shown that the required number of channel uses ensuring a bounded estimation error performance scales in proportional to the degrees of freedom of the channel when the number of RF chains is fixed, whereas it can converge to a constant value if the number of RF chains grows proportionally to the channel dimension while keeping the channel rank fixed. It is worthwhile to note that the proposed technique does not assume any explicit channel model and statistics, making it useful in a wide range of sparse scenarios. The subspace decomposition is a non-convex operation. To effectively solve the problem, we devise an alternating optimization technique that finds a suboptimal but stationary solution to the problem. The proposed technique results in improved channel estimation accuracy compared with prior approaches while consuming a considerably reduced channel use overhead.
This result is interesting since it is achieved without advanced adaptation techniques such as closedloop or two-way adaptations [5], [8]- [12]. We are not aware of similar work establishing the RIP of a specific wireless signaling and estimation problem such as the low-rank millimeter wave MIMO channel estimation problem by drawing a connection between the RIP and estimation performance. Though we focus exclusively on millimeter wave MIMO channel estimation in this work, we hope that our approach will motivate the use of the RIP in other related communication problems.
The rest of the paper is organized as follows. Section II provides the system model and a brief review of the low-rank matrix reconstruction. Section III presents the low-rank subspace decomposition problem and the subspace sampling signal design. In Section IV, the RIP of the subspace sampling signals is analyzed. Section V presents the performance analysis and the algorithm to effectively solve the low-rank subspace decomposition problem. Finally, the numerical simulations and concluding remarks are presented in Sections VI and VII, respectively.
Notations: A bold lower case letter a is a vector, a bold capital letter A is a matrix. A T , A H , A −1 , tr(A), |A|, A F , A * , and a 2 are, respectively, the transpose, conjuagate transpose, inverse, trace, determinant, Frobenius norm, nuclear norm (i.e., the sum of the singular values of A) of A, and l 2 -norm of a.

II. SYSTEM MODEL
We provide the signal model and briefly review the low-rank matrix reconstruction techniques under consideration.

A. Signal Model
Consider a point-to-point millimeter wave hybrid MIMO system. An independent block fading channel with a coherence block length T c (channel uses) is assumed. The open-loop channel estimation is based on sounding a set of predefined channel subspace sampling signals from the transmitter to the receiver by using K (K T c ) channel uses. We denote by N r , N t , and N , respectively, the numbers of receive antennas, transmit antennas, and RF chains (N min(N r , N t )). Fig. 1 illustrates the open-loop channel sounding at the kth channel use, k = 1, . . . , K. The received signal y k is described as where W k ∈ C Nr×N , H ∈ C Nr×Nt , f k ∈ C Nt×1 , and n k ∈ C Nr×1 are, respectively, the receive subspace sampling signal, channel, transmit subspace sampling signal, and the noise.
Seen from Fig. 1, the subspace sampling signals W k and f k in (1) are further divided into the analog and digital components, i.e., where W A,k ∈ C Nr×N (F A,k ∈ C Nt×N ) and W D,k ∈ C N ×N (F D,k ∈ C N ×N ) are the analog and digital receive (transmit) sampling signals, respectively. s k ∈ C N ×1 is the transmitted symbol. The entries of F A,k and W A,k obey the constant modulus property such that |[F A,k ] m,n |=1/ √ N t and |[W A,k ] m,n |=1/ √ N r , ∀m, n, representing the analog phased-array elements. We assume the entries of n k are independent and identically distributed (i.i.d.) Gaussian with zero mean and σ 2 variance. The low-rank channel H has rank(H) = L with L ≤ N min(N r , N t ). With the power constraint E[ f k 2 2 ] = 1, the SNR per channel use is given by 1 σ 2 . For the generation of the subspace sampling signal pair (W k , f k ), we take a simple approach. The transmitter only activates one RF chain per channel use for channel sounding. This is to say the transmit subspace sampling signal is in the form f k = f A,k f D,k s k where f A,k ∈ C Nt×1 , f D,k ∈ C, and s k ∈ C for simplicity. Moreover, we admit the following random generation, where random variables {η k,m,n } and {φ k,n } are mutually independent and i.i.d. according to the uni- 6 form distribution over [0, 2π), i.e., η k,m,n , φ k,n ∈ U[0, 2π), ∀k, m, n. The f k = f A,k f D,k s k in (3) satisfies Although simple, this random design, as we will show in the following sections, permits theoretical insights to conduct channel estimation performance analysis based on the RIP.
The digital baseband signal in (1) cannot directly access H, but it instead observes a noisy version of the analog distorted effective channel W H k Hf k ∈ C N ×1 per channel use. Hence, the channel sounding has to rely on the subspace sampling signal set S = {(W k , f k )} K k=1 . In general, observing the entire H using S requires K = O(N r N t ) channel uses. This overhead could overwhelm the channel coherence time T c in millimeter wave bands, which is over ten-fold shorter than that in sub-6 GHz bands. On the other hand, rank sparsity of the channel matrix is an inherent characteristic of radio channels in the millimeter wave bands [23], [24].

B. Low-Rank Matrix Reconstruction
Low-rank matrix reconstruction has received much recent attention in the field of large-scale sparse signal processing [25]- [33]. In its problem formulation, distorted observations are modeled by using an affine map A. For the sounding model in (1), we can define the affine map A S : Stacking {y k } K k=1 in (1) into a long column vector gives where y = [y T 1 , . . . , y T K ] T ∈ C KN ×1 and the effective noisẽ We will use the notation M KN to indicate the total number of observations. In low-rank matrix reconstruction, popularly adopted approaches include the nuclear norm minimization (NNM) [26] and matrix factorization (MF) [32]- [34]. These approaches can directly be extended to estimate H in (5).

1) Nuclear Norm Minimization (NNM):
The NNM problem in [26] is formulated as This problem is convex and can be solved reliably. However, as seen from (7), it does not guarantee the low-rank solution. One critical drawback is that it does not scale to large-dimensional matrices due to prohibitively high computational complexity [1], [26].

2) Matrix Factorization (MF):
The development of MF [32]- [34] has revealed a new path towards handling large-scale low-rank matrices. It achieves very accurate reconstruction for various classes of matrices at high SNR [25]. Given (5), the MF problem can be formulated as Solving (8) yields the estimate H = L (R ) H . The problem in (8) is non-convex. A widely adopted suboptimal treatment is to use the alternating minimization with power factorization [33]. The complexity of the latter approach is much lower than that of NNM. Interestingly, MF outperforms NNM when the SNR is high and the matrix is substantially low-rank [1], [33], [34]. A critical drawback of MF is that it has very poor performance at low SNR [1]. This is because the objective function in (8) is highly susceptible to noise. For example, detrimental over-matching can frequently happen when the noise overpowers the .

3) Restricted Isometry Property (RIP):
A reliability guarantee of the low-rank matrix reconstruction is best described by the RIP [26], [31]. For any H whose rank is at most L and an arbitrary affine map A, the RIP can be expressed as the bound, where 0 ≤ δ L < 1 is the L-RIP constant [26]. The bound in (9) is justified when A(H) is as significant as if H were directly observed. The reliable reconstruction of NNM [26], [29] and MF [34] is guaranteed with high probability when the A(H) obeys the RIP. The RIP also provides an operational characterization that any A obeying (9) ensures a bounded estimation error for K = O(L · max(N r , N t )) [31].

C. General Statement of Technique
The desired open-loop millimeter wave MIMO channel estimation must leverage the sparse property of the channel while consuming far less channel use overhead than K = O(N r N t ). The RIP offers a useful statistical characterization that ensures a bounded estimation error performance while consuming a reduced amount of channel uses. However, the RIP has only been modeled for the S with i.i.d. Gaussian or Bernoulli entries [26]. The S that hybrid arrays generate based on the random design in (3), however, is in different random ensembles. Moreover, the conventional low-rank matrix reconstructions are critically limited at low SNR. For instance, the rank constraint in (8), i.e., rank(Ĥ) ≤ L, is insufficient to capture the sparsity under the noise. The original MF [32], [33] was originally devised for the noise-free scenario where the rank constraint was sufficient. Note that with the noise, the matrix becomes full rank. The desired open-loop millimeter wave MIMO channel estimation will require further stringent approaches to incorporate the sparsity under the noise.

III. LOW-RANK SUBSPACE DECOMPOSITION
We begin by addressing the questions related to how to enhance the low SNR performance. In MIMO transmission, high resilience to noise is conventionally achieved by matching the precoder and combiner to the channel subspace [35]. This points that when recovering a sparse matrix under the noise, constraints on the channel subspace would be more essential than the rank constraint in (8).

A. Problem Formulation
Consider the unitary precoding and combining gain maximization problem, where W ∈ C Nr×d and F ∈ C Nt×d are semi-unitary matrices, and d ≤ L denotes the dimension of W and F. Equation (10) can be written in an alternative form [1] U, Σ, V = argmin where U ∈ C Nr×d and V ∈ C Nt×d . In (11), Σ ∈ C d×d is a diagonal matrix. The relation between (10) and (11) is that the optimal solution U and V in (11) will also be optimal for (10). However, W and F in (10) are not necessarily optimal for (11). An important implication of (10) and (11) is that maximizing the precoding gain is equivalent to approximating H to a subspace decomposition UΣV H . (4), and including the effective noise in (6) to (11) yields the low-rank subspace decomposition problem, The significance of (12) compared to (8) is that now (12) explicitly contains the subspace constraints and Σ can be exploited to improve robustness at low SNR.
When the noise overpowers the signal, the objective in (12) is likely misled by ñ−A S (UΣV H ) 2 2 . This results in the critical over-matching problem [1]. As a remedy, we propose to control Σ to prevent over-matching by adding a condition to (12), The channel estimate is given by H = U Σ V H . For the selection of β, we first notice that the true channel H is feasible to the problem in (13). This being said, choosing the smallest β satisfying the condition However, we will see in the simulation study in Section VI, that a rough guess of β is sufficient to provide robust estimation performance. Notice that the subspace decomposition in (13) is non-convex due to the coupled variables and quadratic equality constraints. An algorithm for effective minimization of (13) will be discussed in Section V.

B. Random Subspace Sampling Signal
Based on the design in (3), inserting (3) in (1) gives The ith entry of W H A,k Hf A,k ∈ C N ×1 can be rewritten as where (4) can now be rewritten as The following lemma characterizes the distribution of each entry of X k,i ∈ C Nr×Nt . (16) are i.i.d. and each entry follows where d = denotes equality in distribution and θ ∼ U[0, 2π).
Proof: From (16) and (3), one can write where η k,m,i and φ k,n are independent and η k, Since e jθt d = e j mod(θt,2π) , where mod(θ t , 2π) is the modulo 2π of θ t , defining θ mod(θ t , 2π) yields Since η k,m,i and φ k, The distribution in Lemma 1 will be used in the next section to establish the RIP.

IV. RESTRICTED ISOMETRY PROPERTY
In this section, we show that the affine map A S in (16) satisfies the RIP in (9) with high probability.

A. Preliminaries
The following theorem states a sufficient condition for the RIP of an arbitrary affine map A.

Theorem 1 ( [31]):
Let an affine map A : C Nr×Nt → C M ×1 be a random ensemble obeying the following condition: for any H ∈ C Nr×Nt , whose rank is at most L, and fixed 0 < α < 1, for C, c > 0 which depends on α. Then, if M ≥ LC(N t + N r + 1), A satisfies the RIP in (9) with In (21), the probability is taken over the affine map A. Here, A is not restricted to the sounder set S in (16) but it is modeled for any arbitrary S. Reindexing X k,i in (16) as X m with m being a function of (k, i) such that m = (k−1)N +i, k = 1, . . . , K and i = 1, . . . , N , we have It is tractable to transform the sufficient condition in (21) in a vector form where A ∈ C M ×NrNt and The condition in (23) specifies that the random ensembles of A vec(H) 2 2 should be concentrated around vec(H) 2 2 in order to meet the RIP. One such example of A is the matrix with i.i.d. Gaussian entries [36], i.e., [A] m,n ∼ N (0, 1 M ). Another example is when A has entries sampled from an i.i.d. symmetric Bernoulli distribution [37], i.e., Specifically, in [37], it has been shown that the random construction in (24) satisfies (23) with C = 2 and c = α 2 /4−α 3 /6.
Accoridng to Lemma 1, the matrix A associated with S = {W k , f k } K k=1 , which we will denote as A S , has i.i.d. entries following with θ ∈ U[0, 2π). Each entry of (25) can be viewed as a random phase rotation of [A] m,n in (24). In what follows, we will establish the RIP of A S in (25) by showing (23) (i.e., (21) in Theorem 1).

B. RIP Proof
The following lemma establishes a connection between the vector drawn from the distribution in (25) and the vector consisting of Bernoulli random variables in (24).
tively. Then, for any fixed unit norm vector x ∈ C NrNt×1 , x 2 2 = 1, where u is a non-negative integer,x is the element-wise absolute of x, and the expectations are taken √ NrNt e jθ , ∀m, n, where θ ∈ U[0, 2π). Then, for any fixed x ∈ C NrNt×1 with x 2 2 = 1 and 0 < α < 1, the following holds Proof: See Appendix B.
Based on Theorem 1 and Theorem 2, we provide our main theorem of this section as follows. C Nr×Nt , whose rank is at most L, and for any 0 < α < 1, the following holds Moreover, if M ≥ 2L(N t + N r + 1), thenÂ S satisfies the RIP in (9) i.e., the sufficient condition in Theorem 1. Then, by Theorem 1, if (29) holds and A S satisfies the RIP for L-RIP constant δ L with probability exceeding 1 − 2e −qM and fixed constant . This concludes the proof.  N r , N t )). The RIP will further be consolidated in Section V by analyzing the estimation error performance of the low-rank subspace decomposition problem in (13).

C. Asymptotic Analysis
To gain intuition on the bound in (28), we model the distribution of Â S (H) 2 2 in the asymptotically large M , N r , and N t regime. The central limit theorem (CLT) [38] provides a tractable means for this.
Based on (32), the l.h.s. of (28) is upper bounded by where the last step follows from the Chernoff bound of the normal Gaussian random variable [40].
Remark 2: It is interesting to examine the two-sided tail bounds in (29) and (33). The comparison can be made by looking at the two exponents. It can be seen that α 2 2 − α 3 3 ≤ α 2 2 and α 2 2 − α 3 3 → α 2 2 as α → 0, meaning for asymptotically large M , N r , and N t , the bound in (29) converges to the CLT bound in (33) as α tends to zero. This is to say, in the asymptotic sense, the RIP characterization is equivalent to the CLT whose convergence is in probability. A fundamental difference is that the RIP in Theorem 3 gives a tighter characterization than the convergence in probability since it does not require that M , N r , and N t tend to infinity, while the bound in (33) only holds for asymptotically large M , N r , and N t .

V. ESTIMATION ERROR ANALYSIS AND AN ALGORITHM FOR LOW-RANK SUBSPACE DECOMPOSITION
On the basis of the RIP, the estimation error performance of the low-rank subspace decomposition problem in (13) is analyzed. We also provide an algorithm for effective minimization of the low-rank subspace decomposition problem.   (13) where H = H d +H L\d and A S is generated according to (16).

Remark 4:
The fact that the bound (34) is proportional to δ 2d suggests that optimizing the subspace sampling signal set S by minimizing the RIP bound in (9) would result in better execution of the channel estimation because the tighter the RIP bound, the lower δ 2d will be. In this work, we do not aim to optimize S, but designing S to minimize the RIP constant in (9) (i.e., reducing the error bound in (34)) is an interesting future research problem.
Remark 5: When L = d, H L\d = 0 in (34), and this results in the least error bound Now, taking the expectation yields the mean squared error (MSE) bound where we use the fact that E[ ñ 2 2 ] = M σ 2 . It is noteworthy that, in the low SNR regime (i.e., when σ 2 is large), the MSE bounded by 2β, preventing the estimation error from growing along with the noise power σ 2 . This can sharpen the MSE performance at low SNR. Overall, Theorem 4 formulates a required number of channel uses, guaranteeing a bounded estimation error performance with high probability.

B. Alternating Optimization
It is difficult to directly find the solution to (13) due to the semi-unitary constraints, which are nonconvex. We turn to a suboptimal approach that effectively solves the problem. To this end, we first consider the convex relaxation of (13) as where Σ ∈ C d×d is the same diagonal matrix as in (13). The equality constraints are relaxed and the feasible set of (37) becomes convex. However, the problem is still non-convex since the optimization parameters are coupled with each other. The coupling among U, V, and Σ in the objective function makes block coordinate decent [41] an ideal approach for (37) to iteratively optimize U, V, and Σ. In particular, optimizing one parameter by fixing the other two parameters in (37) is convex. Denoting as the interation index, we can iteratively optimize U , V , and Σ for = 1, 2, . . . by solving the following subproblems: (S1) Fix the row subspace V and power allocation Σ , optimize the column subspace U ; (S2) Fix the column subspace U and power allocation Σ , optimize the row subspace V ; (S3) Fix the row subspace V and column subspace U , optimize the power allocation Σ .
A formal description of the alternating optimization is provided in Algorithm 1. In Step 2 of Algorithm 1, Since (S1), (S2), and (S3) are all convex, we can obtain the optimal solution of each subproblem in closed-form. We limit our discussion to the column subspace optimization (S1) and the power matrix optimization (S3), keeping in mind that the same approach for solving (S1) applies to the row subspace optimization (S2). For simplicity, we omit the iteration index attached to the variables. The following lemma provides the solution to (S1). Update the column subspace U ( ) to U ( +1) by solving (S1) subject to tr(U H U) ≤ d, 5: Update the row subspace V ( ) to V ( +1) by solving subject to tr(V H V) ≤ d, 6: Update the subspace power allocation matrix Σ ( ) to Σ ( +1) by solving where Then, the optimal solution U is given by either vec where µ > 0 is the unique solution to the fixed point equation g(µ) = d and g(µ) is monotonically decreasing in µ > 0.
Proof: Suppose the affine map A S in (4). By the Kronecker product equality vec(W H k UΣV H f k ) = ((ΣV H f k ) T ⊗W H k ) vec(U), k = 1, . . . , K and collecting them as a long column vector yields an equivalent problem to (38) as where G is given in (39). This problem is a standard differentiable convex optimization problem, which can be solved by applying Karush-Kuhn-Tucker (KKT) conditions [42]. Since it is a standard procedure, we omit the details here.
We now turn our attention to the subspace power allocation problem (S3). To obtain a closed-form solution to (S3), we have the following lemma.
Lemma 4: Consider the following quadratic program, where y ∈ C M ×1 , U ∈ C Nr×d , V ∈ C Nt×d , and Σ ∈ C d×d is diagonal. Define the matrix P ∈ C M ×d with the ith column The solution Σ is given by either The optimal solution to (43) is found by applying the KKT conditions [42].
Notice that since the subsequent minimizations in (S1), (S2), and (S3) are optimally solved, a sequence is non-increasing and converges to a stationary point over the iterations in Algorithm 1.

A. Simulation Setup
In this section, we evaluate the performance of low-rank subspace decomposition (SD) in Algorithm 1.
We adopt the physical representation of sparse millimeter wave MIMO channels [4], [9], where L scatters are assumed to constitute the propagation paths between the transmitter and receiver, In (44), h l ∈ C is the complex gain of the lth path, i.i.d. according to h l ∼ CN (0, 1), and ξ t,l and ξ r,l are the angles of departures (AoDs) and arrivals (AoAs) at the transmitter and receiver, respectively. Both ξ t,l and ξ r,l are i.i.d. and uniform over [−π/2, π/2]. The a t (ξ t,l ) ∈ C Nt×1 and a r (ξ r,l ) ∈ C Nr×1 are the array response vectors at the transmitter and receiver, respectively. We assume, for simplicity, a t (ξ t,l ) and a r (ξ r,l ) are uniform linear arrays (ULAs) with the antenna spacing set to half of the wavelength. The channel model in (44) Throughout the simulation, we assume N r = 36 and N t = 144, and the parameter β in (13) for the SD method is set to β = E[ H 2 F ] = N r N t . We consider four other estimation techniques as our benchmarks, two open-loop channel estimation technique, i.e., the MF in [34] and CS-based method in [43], and two adaptive channel estimation techniques, i.e., the two-way Arnoldi approach in [8] and closed-loop adaptive channel estimation (ACE) in [9]. The work in [43] provides an open-loop compressed sensing (CS) approach. Unlike the proposed SD, the CS-based method in [43] leverages the sparse representation of millimeter wave channels in angular domain. Hence, in [43], the orthogonal matching pursuit (OMP) is employed to recover the angle supports of the channel matrix. To enhance the sensing capability, it designs the measurement matrix to meet the incoherency property [44]. For the comparison, we admit two performance metrics, i.e., normalized MSE (NMSE) and spectrum efficiency. Specifically, the NMSE is defined as The spectrum efficiency is defined as , and W D ∈ C N ×d , respectively, being the analog precoder, digital precoder, analog combiner, and digital combiner. Here, we consider d-stream multiplexing. For a fair comparison, we admit the near-optimal precoding technique proposed in [7]. This is to say, given the channel estimate H attained by the proposed SD, MF [34], CS-based method [43], ACE [9], and Arnoldi [8] approaches, the precoders and combiners are designed based on the technique in [7] and then R in (46) for each approach is evaluated. For the NMSE evaluation, we involve the proposed SD, MF, CS-based method, and ACE 1 , while we consider all the benchmarks for the spectrum efficiency.

B. Numerical Results
In Fig. 2, the NMSE of the proposed SD is compared with MF, CS-based method, and ACE. We assume the number of paths L = 4, the rank of the desired channel estimate d = 4, and the number of RF chains N = 4. According to Theorem 3, channel use overhead K satisfies K ≥ 2 L N (N r + N t + 1) to meet the RIP. We set the channel uses for SD, MF and ACE at K = 192 and K = 384, which are approximately half and equal to the bound defined in (30), respectively. Because of the assumed angle grid channel model [43], the number of channel uses for the CS-based method is reduced to K = 2N , where N is the number of RF chains. Seen from Fig. 2, the proposed SD results in the most accurate channel estimation performance. Moreover, the NMSE of the proposed SD shows sharp improvement at low SNR and it further decreases as the SNR grows, which is consistent with the analysis in (36).
Specifically, when the SNR is high, i.e., 0 dB ∼ 20 dB, and the number of channel uses is sufficiently large, i.e., K = 384, MF and SD can both achieve accurate estimation. The latter scenario represents an ideal case at high SNR with a sufficiently large number of observations. On the other hand, when K = 192, the MF experiences large performance degradation, while the proposed SD is still compatible with the case of K = 384, demonstrating its robustness. For ACE and CS-based method, even if the SNR and/or the number of channel uses are high, they can not provide an accurate channel estimate.
This happens because the channel estimation of these two methods is based on the angle grids [9] [43].
When the true AoDs and AoAs do not fall exactly on the grid points, these two methods always suffer from residual NMSE even at very high SNR, thus the NMSE flooring in Fig. 2. When the SNR is less than 0 dB in Fig. 2, the performance of SD does not improve by increasing the number of channel uses K. This phenomenon is consistent with the mean square error bound analysis in (36). When the SNR is low, the NMSE in (36) is noise-limited and dominated by 2β, which is irrespective of K. As a result, the proposed SD significantly outperforms MF, especially at low SNR, due to its exploitation of sparse subspace and power regularization.
In Fig. 3, the spectrum efficiency of the proposed SD, MF, ACE, CS-based method, and Arnoldi approaches are illustrated under the same simulation settings as in Fig. 2. The Perfect curve in Fig. 3 is based on perfect CSI and fully digital precoding to evaluate the spectrum efficiency. Given the system dimension, the minimum channel use overhead required for Arnoldi [8] is given by K = 360. When the K for ACE increases to K = 384, its performance becomes compatible with Arnoldi. Seen from Fig. 3, the proposed SD achieves the best spectrum efficiency and both MF and SD closely approach perfect CSI near optimal performance. Note that in Fig. 3, though CS-based method requires only K = 8 channel uses, the spectrum efficiency achieved by it is very low. This is because the CS-based method, assuming the discrete AoDs and AoAs as their channel priors, can not be extended to the practical channel model assumed in (44). Thus, in what follows, we shall only evaluate ACE, Arnoldi, MF, and the proposed SD, and exclude the CS-based method in the rest of the simulations. In Fig. 2, the observation was that when the SNR is in between −10 dB and 5 dB for K = 384, the NMSE performance gap between ACE and SD is minor. However, this is not the case when comparing the spectrum efficiency in Fig. 3. In Fig. 3, the proposed SD achieves considerably improved spectrum efficiency than the ACE. This is the direct consequence of the improved subspace estimation of SD.
In Fig. 4, we evaluate the impact of the number of channel uses K on the spectrum efficiency under the same simulation setting as before. The curves are displayed at SNR 10 dB across various K values.
There is a clear performance improvement for ACE when K = 384 and above. However, ACE requires considerably larger K to further enhance its performance. In the range of K in Fig. 4, the Arnoldi approach [8] can only set K = 360 and K = 720. It should be noted that the Arnoldi approach becomes ideal if the system can be equipped with a larger number of RF chains (e.g., N = N t /8 in [8]). The MF and SD In this set of simulations, the effects of the number of RF chains N to the NMSE and spectrum efficiency are evaluated. In Fig. 6, the NMSE curves for N ∈ {4, 6} are displayed. As illustrated in Fig. 6, for fixed K = 192, the more RF chains the better NMSE performance is expected. It can also be seen that the proposed SD is less sensitive to the change of N at low SNR, while achieving the lowest NMSE. The spectrum efficiency under the different numbers of RF chains is demonstrated in Fig. 7. In Fig. 7, we set  shows the best performance.

C. Complexity Analysis
In this subsection, we compare the computational complexities of the proposed SD, ACE [9], MF [33], and Arnoldi [8] in terms of the number of complex scalar multiplications. Provided the equal channel uses K, the complexity introduced by the precoding and combining is the same for all the approaches. Therefore, we only need to compare the complexities involved in the alternating minimizations for SD and MF, the two-way echoing vector adaptation for Arnoldi [8], and the training precoder and combiner adaptation for ACE [9]. Specifically, for the proposed SD, the computational complexity is determined by the maximum number of iterations max for (S1), (S2), and (S3) and the associated computations in Algorithm 1. Since these three subproblems are all least squares over a sphere, the SVD constitutes the dominant computational overhead. Hence, for a fixed max , the complexity of SD is given by O(KN d 2 (N 2 t + N 2 r )). The same overhead O(KN d 2 (N 2 t + N 2 r )) results for MF. The number of AoD and AoA search levels of ACE in [9] is given by O (KN/d 3 ). At each search level of ACE, the complexity is dominated by the computations of training precoders and combiners, which is given by [9], where D ( max(N t , N r )) denotes the cardinality of an over-complete dictionary. Therefore, the complexity of the ACE is given by O (KN D(N r + N t )/d 2 ). The main computational complexity of Arnoldi is the design of channel echoing vectors, which is given by O (K 2 N 2 /(N r + N t )) [8]. In summary, the Arnoldi method has the lowest complexity O (K 2 N 2 /(N r + N t )). The proposed SD has the similar complexity with MF and ACE, which is about O(KN d 2 (N 2 t + N 2 r )).

VII. CONCLUSIONS
In this paper, we have proposed an open-loop low-rank subspace decomposition technique for millimeter wave MIMO channel estimation. The RIP that ensures reliability of the low-rank matrix reconstruction was established for the channel subspace sampling signals randomly generated by the hybrid analog and digital arrays. When the number of channel uses is around L N (N t + N r ), the considered random subspace sampling signals satisfy the RIP with high probability. The established RIP was applied to analyze the channel estimation error performance. This analysis showed the resilience of the proposed technique at low SNR. Moreover, we devised an alternating optimization algorithm that effectively finds a suboptimal but stationary solution to the problem. The simulation studies corroborate our analysis and showed that the proposed low-rank subspace decomposition technique can achieve near optimal performance at high SNR, while offering robust performance at low SNR with channel use overhead less  Both sides of (26) can be expanded as Since a part of the proof is based on a similar technique in [37], we first introduce the following lemma, which will be used in the proof. and We first focus on the upper tail probability in (27), i.e., Pr( N t N r /M A S x 2 2 > 1 + α). For any h > 0, applying Markov's inequality to the upper tail probability yields Since the inequality holds for any h > 0, insearting h = NrNtα 2(1+α) in (53) gives where the last inequality follows from (1 + α) M 2 e −α M 2 = e M 2 (ln(1+α)−α) and ln(1 + α) ≤ α − α 2 /2 + α 3 /3, ∀α ∈ (0, 1). This proves the upper tail bound in (27).
Using the same approaches in (51), the lower tail bound of (27) is upper bounded by Then, plugging in (55) into (54) with the substitution h = NrNtα 2(1+α) results in where the last inequality follows from 1− α(4+α) We start by realizing that the desired channel H d = d l=1 λ l ψ l φ H l is feasible to (13). Thus, one should have Inserting y = A S (H d )+A S (H L\d )+ñ in (57) yields Applying the triangular inequality to the l.h.s. of (58) gives Thus, combining (58) and (59) results in Note that rank( H − H d ) ≤ 2d. Then, based on the RIP of NrNt M A S (·) by Theorem 3, we have where δ 2d is 2d-RIP constant. Hence, combining (60) and (61) results in the error bound The inequality constraint in (13) is equivalent to Σ 2 F ≤ β. Moreover, since H d is feasible to (13), one should have by the triangular inequality Finally, combining (62) and (63) leads to This bound holds, by Theorem 3, with probability greater than or equal to 1 − exp(−qM ) when M ≥ 4d(N t + N r + 1) for q = α 2 4 − α 3 6 − ln(36 and 0 < α < 1.