Design of Robust Adaptive Beamforming Algorithms Based on Low-Rank and Cross-Correlation Techniques

This work presents cost-effective low-rank techniques for designing robust adaptive beamforming (RAB) algorithms. The proposed algorithms are based on the exploitation of the cross-correlation between the array observation data and the output of the beamformer. Firstly, we construct a general linear equation considered in large dimensions whose solution yields the steering vector mismatch. Then, we employ the idea of the full orthogonalization method (FOM), an orthogonal Krylov subspace based method, to iteratively estimate the steering vector mismatch in a reduced-dimensional subspace, resulting in the proposed orthogonal Krylov subspace projection mismatch estimation (OKSPME) method. We also devise adaptive algorithms based on stochastic gradient (SG) and conjugate gradient (CG) techniques to update the beamforming weights with low complexity and avoid any costly matrix inversion. The main advantages of the proposed low-rank and mismatch estimation techniques are their cost-effectiveness when dealing with high dimension subspaces or large sensor arrays. Simulations results show excellent performance in terms of the output signal-to-interference-plus-noise ratio (SINR) of the beamformer among all the compared RAB methods.

I. INTRODUCTION Adaptive beamforming has been one of the most important research areas in sensor array signal processing. It has also been recognized that traditional adaptive beamformers are extremely sensitive to environmental uncertainties or steering vector mismatches, which may be caused by many different factors (e.g., imprecise antenna size calibration, signal pointing errors or local scattering). Furthermore, some radar systems in advanced applications require antenna arrays with a very large number of sensor elements in highly dynamic environments, which leads to the increase of computational complexity and the decrease of the convergence rate for computing the parameters of the beamformer.

A. Prior and Related Work
In order to mitigate the effects of uncertainties on adaptive beamformers, robust adaptive beamforming (RAB) techniques have been developed. Popular approaches include worst-case optimization [2], [9], diagonal loading [4], and eigen-subspace decomposition and projection techniques [6], [8], [11]. However, these RAB approaches have some limitations such as their ad hoc nature, high probability of subspace swap at low signal-to-noise ratio (SNR) [3] and high computational cost due to online optimization or subspace decomposition techniques.
Furthermore, in the case of large sensor arrays the above mentioned RAB methods may encounter problems for their application. This is because in these RAB algorithms, a cubic or greater computational cost is required to compute the beamforming parameters. Therefore, dimensionality reduction (or rank-reduction) methods ( [15]- [31], [33]- [36]) have been employed and developed to reduce the complexity and improve the convergence rate.

B. Contributions
In this work, we propose and study novel RAB algorithms that are based on low-rank and cross-correlation techniques.
In the proposed techniques, we exploit the prior knowledge that the steering vector mismatch of the desired signal is located within an angular sector which is assumed known. The proposed algorithms are based on the exploitation of the cross-correlation between the array observation data and the output of the beamformer, which avoids costly optimization procedures. We firstly construct a linear system (considered in high dimension) involving the mismatched steering vector and the statistics of the sampled data. Then we employ an iterative full orthogonalization method (FOM) [12], [13] to compute an orthogonal Krylov subspace whose model order is determined by both the minimum sufficient rank [17], which ensures no information loss when capturing the signal of interest (SoI) with interferers, and the execute-and-stop criterion of FOM [12], [13], which automatically avoids overestimating the number of bases of the computed subspace. The estimated vector that contains the cross-correlation between the array observation data and the beamformer output is projected onto the Krylov subspace, in order to update the steering vector mismatch, resulting in the proposed orthogonal Krylov subspace projection mismatch estimation (OKSPME) method.
Moreover, we develop an analysis of the mean squared error (MSE) between the estimated and the actual steering vectors for the general approach of using a presumed angular sector associated with subspace projections. This analysis mathematically describes how precise the steering vector mismatch can be estimated. Upper and lower bounds are derived and compared with the approach in [6]. Another analysis on the computational complexity of the proposed and existing algorithms is also provided.
In the simulations, we consider local scattering scenarios (both coherent and incoherent) to model the mismatch effects. We also study the performance of the proposed algorithms by testing the output signal-to-interference-plus-noise ratio (SINR) of the beamformer with respect to training snapshots and different input SNRs. The number of sensor elements and interferers is also varied and compared in each scenario to provide a comprehensive performance study. In summary, the contributions of this work are: • The proposed OKSPME RAB method.
• The development of the modified SG and CG type OKSPME RAB algorithms. • An analysis of the computational complexity and the MSE performance of the proposed and existing RAB algorithms. The remaining sections of this paper are organized as follows: The system model and problem statement are described in Section II. Section III introduces the proposed OKSPME method, whereas Section IV introduces the proposed robust adaptive algorithms. Section V provides the MSE analysis of the steering vector estimation and the complexity analysis. Section VI presents and discusses the simulation results. Section VII gives the conclusion.

II. SYSTEM MODEL AND PROBLEM STATEMENT
Let us consider a linear antenna array of M sensors and K narrowband signals which impinge on the array. The data received at the ith snapshot can be modeled as where s(i) ∈ C K×1 are uncorrelated source signals, θ = [θ 1 , · · · , θ K ] T ∈ R K is a vector containing the directions of arrival (DoAs) and [.] T denotes the transpose, A(θ) = [a(θ 1 )+e, · · · , a(θ K )] = [a 1 , · · · , a K ] ∈ C M×K is the matrix which contains the steering vector for each DoA and e is the steering vector mismatch of the desired signal, n(i) ∈ C M×1 is assumed to be complex circular Gaussian noise with zero mean and variance σ 2 n . The beamformer output is given by where w = [w 1 , · · · , w M ] T ∈ C M×1 is the beamformer weight vector, where (·) H denotes the Hermitian transpose. The optimum beamformer is computed by maximizing the SINR and is given by where σ 2 1 is the desired signal power, R I+N is the interference-plus-noise covariance (INC) matrix. The problem of maximizing the SINR in (3) can be cast as the following optimization problem: which is known as the MVDR beamformer or Capon beamformer [1], [4]. The optimum weight vector is given by Since R I+N is usually unknown in practice, it can be estimated by the sample covariance matrix (SCM) of the received data asR Using the SCM for directly computing the weights will lead to the sample matrix inversion (SMI) beamformer w SMI = R −1 a1 a H 1R −1 a1 . However, the SMI beamformer requires a large number of snapshots to converge and is sensitive to steering vector mismatches [2], [3]. As previously mentioned, most of the conventional and existing RAB algorithms are computationally costly especially when encountering arrays with a very large number of sensors. Therefore, the RAB design problem we are interested in solving includes the following aspects: • To design cost-efficient algorithms that are robust against uncertainties and values of SNRs and interferers in the presence of uncertainties in the steering vector of a desired signal. • The proposed algorithms must preserve their robustness and low-complexity features for large sensor arrays.

III. PROPOSED OKSPME METHOD
In this section, the proposed OKSPME method is introduced. This method aims to construct a linear system involving only known or estimated statistics and then projects an estimated cross-correlation vector between the array observation data and the beamformer output onto an orthogonal Krylov subspace, in order to update the steering vector mismatch with reduced complexity. The SCM of the array observation data is estimated by (5). The cross-correlation vector between the array observation data and the beamformer output can be expressed as d = E[xy * ] (where [.] * denotes complex conjugation) or equivalently as Assuming that |a H k w| ≪ |a H 1 w| for k = 2, · · · , K and all signals have zero mean, the cross-correlation vector d can be rewritten as Note that we also assume that the desired signal is statistically independent from the interferers and the noise, i.e., E[s k s * 1 ] = 0 and E[s k a k s * 1 a H 1 w] = 0 for k = 2, · · · , K. With this assumption the desired signal power is not statistically affected by the interference and (7) can be rewritten as which can be estimated by the sample cross-correlation vector (SCV) given byd

A. Desired Signal Power Estimation
In this subsection, we describe an iterative method for the desired signal power (σ 2 1 ) estimation based on our prior work in [32], which can be accomplished by directly using the desired signal steering vector. Alternatively, a designer can employ a maximum likelihood (ML) or a minimum variance (MV) estimator for computing the desired signal power. However, the approach described here has lower complexity than the ML and the MV estimators. In the adopted method, we need to choose an initial guess for the steering vector mismatch within the presumed angular sector, sayâ 1 (0) and setâ 1 (1) =â 1 (0). By adding the snapshot index i, we can rewrite the array observation data as whereâ 1 (0) andâ 1 (i) (i = 1, 2, · · · ) designate the initial guess of the steering vector and its estimate at the ith snapshot, respectively. Pre-multiplying the above equation byâ H 1 (i) we havê (11) Here we assume that each of the interferers is orthogonal or approximately orthogonal to the desired signal. Specifically, the steering vector of each of the interferers is orthogonal (â H 1 (i)a k (i) = 0, k = 2, 3, · · · , K), or approximately orthogonal (â H 1 (i)a k (i) ≪â H 1 (i)â 1 (i), k = 2, 3, · · · , K) to the desired signal steering vector (i.e.,â 1 (i)), so that a H 1 (i)a k (i) (k = 2, 3, · · · , K) approaches zero and the term (11) can be neglected, resulting in Taking the expectation of |â H 1 (i)x(i)| 2 , we obtain . (13) Assuming that the noise is statistically independent from the desired signal, then we have where E[n(i)n H (i)] represents the noise covariance matrix R n (i) that can be replaced by σ 2 n I M , where the noise variance σ 2 n can be easily estimated by a specific estimation method. A possible approach is to use a Maximum Likelihood (ML) based method as in [14]. Replacing the desired signal power E[|s 1 (i)| 2 ] and the noise variance σ 2 n by their estimatesσ 2 1 (i) andσ 2 n (i), respectively, we obtain The expression in (15) has a low complexity (O(M )) and can be directly implemented if the desired signal steering vector and the noise level are accurately estimated.

B. Orthogonal Krylov Subspace Approach for Steering Vector Mismatch Estimation
An orthogonal Krylov subspace strategy is proposed in order to estimate the mismatch with reduced cost and deal with situations in which the model order is time-varying. Our idea is based on constructing a linear system, which considers the steering vector mismatch as the solution, and solving it by using an iterative Krylov subspace projection method. To this end, consider a general high-dimensional linear system model given by where B ∈ C M×M and b ∈ C M×1 . Then we need to express B and b only using available information (known statistics or estimated parameters), so that we can solve the linear system with the Krylov subspace of order m (m ≪ M ) described by Taking the complex conjugate of (12), we have Pre-multiplying both sides of (18) by the terms of (10), then adding an extra term δIâ 1 (i) (where δ is a small positive number defined by the user) and simplifying the terms, we obtain in which by further defining the expression on the right-hand side asb(i), we can rewrite (20) aŝ As can be seen (21) shares the same form as the linear system of equations in (16) andb(i) can be expressed in terms of a 1 (i),σ 2 1 (i) andσ 2 n (i) whereasR(i) can be estimated by (5). In the following step, we employ the Arnoldi-modified Gram-Schmidt algorithm from the FOM method [12], [13] associated with the minimum sufficient rank criterion discussed in [17] to compute an orthogonal Krylov subspace. We define a residue vector to represent the estimation error in the ith snapshot aŝ and let Then the Krylov subspace bases can be computed using the modified Arnoldi-modified Gram-Schmidt algorithm as in Table I (the snapshot index i is omitted here for simplicity).
In Table I, <, > denotes the inner product and the parameters h l,j (l, j = 1, 2, · · · , m) are real-valued coefficients, the model order is determined once if one of the following situations is satisfied: • The execute-and-stop criterion of the original Arnoldimodified Gram-Schmidt algorithm is satisfied (i.e., h j,j+1 = 0). • The minimum sufficient rank for addressing the SoI and interferers is achieved (i.e., j ≥ K + 1, where K is Compute u j =Rt j For l = 1, 2, · · · , j, do: the number of signal sources), so that no more subspace components are necessary for capturing the SoI from all the existing signal sources. Now by inserting the snapshot index, we havê and the Krylov subspace projection matrix is computed bŷ It should be emphasized that the Krylov subspace matrix T(i) obtained here is constructed by starting with the residue vectorr(i). In other words,T(i) is constructed with the estimation error of the steering vector. In order to extract the estimation error information and use it to update the steering vector mismatch, we can project the SCVd(i) in (9) ontô P(i) and add the estimation error to the current estimate of a 1 (i) asâ

C. INC Matrix and Beamformer Weight Vector Computation
Since we have estimated both the desired signal power σ 2 1 (i) and the mismatched steering vector in the previous subsections, the INC matrix can be obtained by subtracting the desired signal covariance matrix out from the SCM aŝ The beamformer weight vector is computed bŷ which has a computationally costly matrix inversionR −1 I+N (i). The proposed OKSPME method is summarized in Table II. In the next section, we will introduce adaptive algorithms to avoid matrix inversions and reduce the complexity.

IV. PROPOSED ADAPTIVE ALGORITHMS
This section presents adaptive strategies based on the OK-SPME robust beamforming method, resulting in the proposed OKSPME-SG, OKSPME-CCG and OKSPME-MCG algorithms, which are especially suitable for dynamic scenarios. In the proposed adaptive algorithms, we estimate the desired Choose an initial guessâ1(0) within the sector and setâ1(1) =â1(0); For each snapshot i = 1, 2, · · · : Step 1. Compute the desired signal power Step 2. Determine the Krylov subspacê Apply the algorithm in Table I to determine m and t1(i),· · · ,tm(i) Step 3. Update the steering vector End snapshot signal power and its steering vector with the same recursions as in OKSPME, whereas the estimation procedure of the beamforming weights is different. In particular, we start from a reformulated optimization problem and use SG and CG-based adaptive recursions to derive the weight update equations, which reduce the complexity by an order of magnitude as compared to that of OKSPME.

A. OKSPME-SG Adaptive Algorithm
We resort to an SG adaptive strategy and consider the following optimization problem: whereR(i) can be written as x(i)x H (i) andR 1 (i) represents the desired signal covariance matrix and can be written aŝ Then we can express the SG recursion as where and µ is the step size.
By substituting L into the SG equation (30) and letting w H (i + 1)â 1 (i + 1) = 1, λ L is obtained as By substituting λ L in (30) again, the weight update equation for OKSPME-SG is obtained as ).
(32) The adaptive SG recursion circumvents a matrix inversion when computing the weights using (28), which is unavoidable in OKSPME. Therefore, the computational complexity is (32) is ensured converging to a solution. To implement OKSPME-SG, Step 1, Step 2 and Step 3 from Table II and (32) are required.

B. OKSPME-CCG Adaptive Algorithm
In this subsection, the proposed OKSPME-CCG algorithm is introduced. In CG-based approaches, we usually employ a forgetting factor (e.g. λ) to estimate the second-order statistics of the data or the SCM [1], [10], which can be expressed bŷ whereas the SCVd(i) can be estimated with the same forgetting factor as described bŷ The proposed optimization problem that leads to the OKSPME-CCG algorithm is described by where v(i) is the CG-based weight vector. In OKSPME-CCG, we require N iterations for each snapshot. In the nth iteration, a 1,n (i) and v n (i) are updated as followŝ where pâ 1,n (i) and p v,n (i) are direction vectors updated by where gâ 1,n (i) and g v,n (i) are the negative gradients of the cost function in terms ofâ 1 (i) and v(i), respectively, which are expressed as The scaling parameters αâ 1,n (i), α v,n (i) can be obtained by substituting (36) and (37) into (35) and minimizing the cost function with respect to αâ 1,n (i) and α v,n (i), respectively. The solutions are given by .
(43) The parameters βâ 1,n (i) and β v,n (i) should be chosen to provide conjugacy for direction vectors [10], which results in Afterâ 1,n (i) and v n (i) are updated for N iterations, the beamforming weight vector w(i) can be computed by The computational cost of OKSPME-CCG algorithm is O(N M 2 ), which is higher than the cost required in OKSPME-SG due to the inner iterations at every snapshot. The proposed OKSPME-CCG is summarized in Table III.

C. OKSPME-MCG Adaptive Algorithm
In OKSPME-MCG, we let only one iteration be performed per snapshot, which further reduces the complexity compared to OKSPME-CCG. Here we denote the CG-based weights and steering vector updated by snapshots rather than inner iterations asâ As can be seen, the subscripts of all the quantities for inner iterations are eliminated. Then, we employ the degenerated scheme to ensure αâ 1 (i) and α v (i) satisfy the convergence bound [10] given by Instead of updating the negative gradient vectors gâ 1 (i) and g v (i) in iterations, now we utilize the forgetting factor to reexpress them in one snapshot as Choose an initial guessâ1(0) within the sector and setâ1(1) =â1(0); For each snapshot i = 1, 2, · · · : Step 1 from Table II Step 2 from Table II Step 3 from Table II Steering Vector and Weight Vector Estimationŝ For each iteration index n = 1, 2, · · · , N :

End snapshot
Pre-multiplying (51) and (52) by p Ĥ a1 (i) and p H v (i), respectively, and taking expectations we obtain where in (54) After substituting (54) in (50) we obtain the bounds for α v (i) as follows Then we can introduce a constant parameter η v ∈ [0, 0.5] to restrict α v (i) within the bounds in (55) as Similarly, we can also obtain the bounds for αâ 1 (i).

For simplicity let us define
in which we can introduce another constant parameter ηâ 1 ∈ [0, 0.5] to restrict αâ 1 (i) within the bounds in (57) as or (59) Then we can update the direction vectors pâ 1 (i) and p v (i) by pâ 1 (i + 1) = gâ 1 (i) + βâ 1 (i)pâ 1 (i), where βâ 1 (i) and β v (i) are updated by Finally we can update the beamforming weights by The MCG approach employs the forgetting factor λ and constant η for estimating α(i), which means its performance may depend on a suitable choice of these parameters. The proposed OKSPME-MCG algorithm requires a complexity of O(M 2 ). However, the cost is usually much lower compared to CCG approach for the elimination of inner recursions and it presents a similar performance in most studied scenarios. From an implementation point of view, the choice of using the CCG and MCG algorithms is based on the stationarity of the system: the CCG algorithm is more suitable for scenarios in which the system is stationary and we can compute the beamformer with K iterations while the MCG algorithm is suggested for non-stationary scenarios as we only run one iteration per snapshot and can track variations in the environment. Table  IV summarizes the OKSPME-MCG algorithm. w(1) = v(0) = 1; λ; ηv = ηâ 1 ; Choose an initial guessâ1(0) within the sector and setâ1(1) =â1(0); gv(0) = pv(1) =â1(1); gâ 1 (0) = pâ 1 (1) = v(0); For each snapshot i = 1, 2, · · · : Step 1 from Table II Step 2 from Table II Step 3 from Table II Steering Vector and Weight Vector Estimations In this section, we present an analysis of the following aspects of the proposed and existing algorithms: • An analysis of the MSE between the estimated and actual steering vectors for the general approach that employs a presumed angular sector. • MSE analysis results of the proposed OKSPME method and the SQP method in [6] and their relationships and differences. • A complexity analysis for the proposed and existing algorithms.

A. MSE analysis
Firstly, we present the MSE analysis of the general approach that employs a presumed angular sector. Since we have the steering vector estimateâ 1 (i) in the ith snapshot, by denoting the true steering vector as a 1 , we can express the MSE of the estimateâ 1 (i) as In the approach that employs an angular sector, we usually choose an initial guess (i.e.,â 1 (0)) from the presumed sector. Let us express the accumulated estimation error aŝ then (65) can be rewritten as (67) The initial guessâ 1 (0) can be described as the true steering vector plus a guess error vector (i.e., ǫ): Taking expectation of both sides of the above, we have Substituting (68) into (67), taking into account that the accumulated estimation error is uncorrelated with the initial guess error and simplifying the expression, we obtain Furthermore, it should be emphasized that both ǫ andê(i) are in vector forms, which means that their second-order statistics can be re-expressed in terms of their first-order statistics of their Euclidean norms. Then we can re-express (70) as Since both ǫ and ê(i) are scalars we have At this stage, we can employ Popoviciu's inequality [37] to obtain the upper bounds for the variances of the norms of the random vectors ǫ andê(i), which are given by However, the last term in (71) is not analytical when conducting a norm analysis. Actually, E[ǫ] depends on how the presumed sector is chosen: if the sector is chosen in an unbiased manner (i.e., the true steering vector lies in the centre of the sector), then we have E[ǫ] = 0 by symmetry criterion, in which case we can omit the last terms of (71). For convenience of carrying out the norm analysis as the next step, we focus on the unbiased case only, so that the MSE only depends on the expectation, the infimum and the supremum of ǫ and ê(i) . In Fig. 1, we utilize Euclidean geometry to illustrate the relationships among the norms of the errors and the norm of the steering vector, which is a fixed parameter due to the re-normalization procedure after it is estimated each time.
According to Fig. 1, we can use θ (i.e., half of the angular sector, assumed less than π/4) and a 1 to obtain E[ ǫ ] by the following (any angular parameter appeared in the equations should be measured in radians rather than degrees): ǫ is equivalent to the chord length which corresponds to the arc of a variable τ , which can be any value from 0 to θ with equal probability, in other words, the choice of τ is uniformly distributed within [0, θ]. If the sample size of the selected ǫ is large enough, we can approximately describe its probability density function (pdf) as a continuous function given by Meanwhile, we are also able to calculate the chord length ǫ from a simple geometric criterion as Then the expectation of ǫ can be computed by from which after simplification we obtain At this point, we can also compute the variance of ǫ by using (79) as from which after simplification we obtain In addition, it is clear that we have inf ǫ = 0 and sup ǫ = 2 a 1 sin θ 2 , which can be substituted in (74) and result in We can see that the right-hand side of (81) satisfies the inequality in (82). After substituting (79) and (81) in (72), we obtain   Regarding the computation of the norm of the accumulated estimation error ê(i) , we need to emphasize that even though the steering vector is always re-normalized each time after it is updated, the piecewise estimation error in each snapshot does not directly update the steering vector to its normalized version, which means it is inappropriate to calculate the estimation error by geometric methods directly from Fig. 1 because the accumulated estimation error partially comes from the unnormalized steering vectors. However, we can obtain the infimum and supremum values for ê(i) if we assume the update scheme is unidirectional (i.e., the steering vector is updated fromâ 1 (0) toâ 1 (i) in a single direction within the sector), with the unnormalized steering vectors considered.
We firstly look at the SQP method scenario in [6]. The steering vector update scheme is shown in Fig. 2. It is necessary to emphasize that now we focus on the angular sector range of θ i (i.e., the angle difference between the initially guessed steering vector and its estimate in the ith snapshot) rather than θ. In [6], an online-optimization program was used to iteratively solve for the piecewise estimation error in every snapshot, which was always orthogonal to the current steering vector estimate. Let us consider that at each time instant the steering vector is updated, its direction changes by θ i,k , where i is the snapshot index and k (1 ≤ k ≤ i) is the index for the kth update. Since the total direction change in a snapshot is θ i , then we have and the norm of the accumulated estimation error is no greater than the sum of the norms of all the piecewise estimation errors, which is given by the inequality If we assume θ i is less than π/2, then the right-hand side of (85) achieves its maximum value when θ i,k = tan θ i , which is also the supremum of ê(i) and equals On the other hand, we notice that the piecewise estimation error vector can never enter into the angular sector, but at most move along with the arc if the number of iterations is large enough. In this case, we can approximately and geometrically illustrate the arc length corresponding with θ i as the lower bound by taking the limit i → ∞, i.e., which is actually the infimum of ê(i) and cannot be achieved since the number of snapshots or iterations are always limited in practical situations. By combining (86) and (87), Different from the SQP method, the proposed OKSPME method utilizes the Krylov subspace and the cross-correlation vector projection approach to extract the error information then use it to update the steering vector. From (9) we havê Note that an initialization for vectord or matrixR should be considered to ensureR is full-rank and invertible, which can be done by either settingd(0) = δIw(0) orR(0) = δI. We also know that Pre-multiplying (90) byR(k) on both sides we obtain which is then substituted in (89) and results in whereσ 2 1 (k) is a scalar, which means the SCV contains the direction of the desired signal steering vector. Projectingd(i) onto the Krylov subspace represented byP(i) is therefore similar to projectingâ 1 (i). In our method, the estimation ofd(i) is separate from the update ofâ 1 (i), which means the steering vector estimation error used for the updates is obtained fromd(i), so that in the kth (1≤k < i) snapshot, the error does not have to be orthogonal toâ 1 (k), but should be orthogonal to another potentially better estimateâ 1 (j) (1≤k < j≤i), resulting in a situation where the error is located inside the sector (see Fig. 3). There are two benefits in the case which the error is inside the sector: faster convergence rate and smaller estimation error. We can obtain the infimum and supremum values in a similar way. By applying the inequality that the norm of the accumulated estimation error is no greater than the sum of the norms of all the piecewise estimation errors, we have where the parameters θ i,k (k = 1, 2, · · · , i) satisfy the constraint in (84). However, the right-hand side of (93) achieves its maximum value when all these parameters are equal (i.e., θ i,1 = θ i,2 = · · · = θ i,i = θi i ) and it is given by The right-hand side of (94) can be treated as a function of i which is an increasing function on i = 1, 2, · · · , ∞. Therefore, we can take the limit of it to obtain the upper bound of ê(i) max , and so as to ê(i) . In fact, when i → ∞, the piecewise estimation error moves along the arc corresponding with θ i , resulting in the upper bound obtained is the same as the lower bound of the SQP method case, which is given by the right-hand side expression of (87) and defines the supremum of ê(i) in this case. Since we have already assumed that θ i is less than π/2 so thatê(i) must be inside of the angular sector but its Euclidean norm cannot be smaller than the orthogonal distance betweenâ 1 (0) toâ 1 (i), so this orthogonal distance can define the lower bound of ê(i) , which is actually the infimum and calculated by a 1 sin θ i . Then, in the OKSPME method, ê(i) is bounded by inf ê(i) = a 1 sin θ i ≤ ê(i) < θ i a 1 = sup ê(i) .
(95) By taking expectations of both (88) and (95), we obtain On the other side, by substituting (88) and (95) in (75), we obtain Substituting (96), (98) and (97), (99) in (73), respectively, we have However, E[θ i ] also has its lower and upper bounds. Since our analysis focuses on the unbiased case only as mentioned, the true steering vector is located in the center of the angular sector and the estimateâ 1 (i) is always closer to the center thanâ 1 (0). Let us assume that even if the estimateâ 1 (i) always happens to be very close to either edge of the sector, no matter howâ 1 (0) is chosen within the sector, θ i will vary from 0 to 2θ with equal probability, or equivalently, uniformly distributed within  (100) and (101) respectively, resulting in Finally, by combining the expectation of the mean-squared initial guess error E[ ǫ 2 ] in (83) with (102) and (103), we obtain the bounds for the MSE of the steering vector estimate MSE{â 1 (i)} as From (104) and (105), we can see that the MSEs now only depend on two parameters: the norm of the true steering vector and the angular sector spread. The lower and upper bounds of the proposed OKSPME method are lower than those of the SQP method. As mentioned before, it is important that the presumed angular sector spread 2θ must be less than π/2 (i.e., 90 • ) to ensure the previous assumption of θ i < π/2 is always valid.

B. Complexity Analysis
The computational complexity analysis is discussed in this subsection. We measure the total number of additions and multiplications (i.e., flops) in terms of the number of sensors M performed for each snapshot for the proposed algorithms and the existing ones and list them in Table V. Note that the SQP method in [6] has a highly-variant computational complexity in different snapshots, due to the online optimization program based on random choices of the presumed steering vector. However, it is usually in O(M 3.5 ). The complexity of the LCWC algorithm in [9] often requires a much larger n than that in the proposed LOCSME-CCG algorithm. It is obvious that all of the proposed algorithms have their complexity depending on the Krylov subspace model order m, which is determined from Table I and is no larger than K + 1. For the convenience of comparison, we eliminate all parameters except M by setting them to common values (the values of n in LCWC and OKSPME-CCG is set to 50 and 5 respectively, m = K + 1 where K = 3) and illustrate their complexity with M varying from 10 to 100 as shown in Fig. 4. As can be seen that the proposed OKSPME-SG and OKSPME-MCG algorithms have lower complexity than the other algorithms.

VI. SIMULATIONS
In this section, we present and discuss the simulation results of the proposed RAB algorithms by comparing them to some of the existing RAB techniques. We consider uniform linear arrays (ULA) of omnidirectional sensors with half wavelength spacing. To produce all the figures (if unspecified in a few scenario), 100 repetitions are executed to obtain each point of the curves and a maximum of i = 300 snapshots are observed. The desired signal is assumed to arrive at θ 1 = 10 • . The signal-to-interference ratio (SIR) is fixed at 0dB. As the prior knowledge, the angular sector in which the desired signal is assumed to be located is chosen as [θ 1 −5 • , θ 1 +5 • ]. The results focus on the beamformer output SINR performance versus the number of snapshots, or a variation of input SNR (−10dB to 30dB) and both coherent and incoherent local scattering mismatch [5] scenarios are considered. RCB [4] SQP [6] LOCME [7] LCWC [9] OKSPME OKSPME−SG OKSPME−CCG OKSPME−MCG

A. Mismatch due to Coherent Local Scattering
All simulations in this subsection consider coherent local scattering. With time-invariant coherent local scattering, the steering vector of the desired signal is modeled as where p corresponds to the direct path while b(θ k )(k = 1, 2, 3, 4) corresponds to the scattered paths. The angles θ k (k = 1, 2, 3, 4) are randomly and independently drawn in each simulation run from a uniform generator with mean 10 • and standard deviation 2 • . The angles ϕ k (k = 1, 2, 3, 4) are independently and uniformly taken from the interval [0, 2π] in each simulation run. Both θ k and ϕ k change from trials while remaining constant over snapshots. We firstly compare our proposed methods with some classical RAB methods (i.e., standard diagonal loading method with a fixed loading factor equal to 10 times the noise variance, the RCB method in [4] which estimates the loading factor iteratively, and the method that solves an online quadratic optimization programming, which refers to the SQP method [6]). The numbers of sensors and signal sources (including the desired signal) are set to M = 10 and K = 3, respectively. For this case only, we set the INR (interferences-to-noise ratio) to 20dB and illustrate the SINR performance versus snapshots within 100 snapshots in Fig. 5. The two interferers are arranged to be in the directions of θ 2 = 30 • and θ 3 = 50 • , respectively. The other user-defined parameters, if unspecified, (e.g. the step size µ and the forgetting factor λ) are manually optimized to give the best algorithm performance, which is also applied for the other simulation scenarios.
We then set the number of sensors to M = 12, the number of signal sources as (including the desired signal) K = 3 and illustrate the SINR versus snapshots and the SINR versus input SNR performance in Fig. 6 and Fig. 7 respectively.  The two interferers are arranged to be in the directions of θ 2 = 30 • and θ 3 = 50 • , respectively. In either Fig. 6 or Fig. 7, we can see that the proposed OKSPME method has a very similar or slightly better performance compared to the LOCSME algorithm of [32] and both of them have the best performance. Furthermore, the proposed OKSPME-CCG and OKSPME-MCG algorithms also achieve very close performance to OKSPME. In Fig. 8, we assess the SINR performance versus snapshots of those selected algorithms in a specific time-varying scenario which encounters a halfway redistribution of the interferers at a certain snapshot. In this case, the number of sensors is kept at M = 12, whereas the details of the interferers are given in Table VI.
In Figs. 9 and 10, we set the number of signal sources to K = 3, but increase the number of sensors from M = 12 to M = 40 and study the SINR versus snapshots and the SINR versus input SNR performance of the selected and proposed dimensionality reduction RAB algorithms, respectively. We set the reduced-dimension as D = 4 for the beamspace based algorithm [23] in all simulations. This time, it is clear that the proposed OKSPME, OKSPME-SG, OKSPME-CCG and OKSPME-MCG algorithms all have a certain level of performance degradation compared to the scenario where M = 12. The proposed OKSPME based algorithms achieve better performances than the beamspace approach.

B. Mismatch due to Incoherent Local Scattering
In this case, the desired signal affected by incoherent local scattering has a time-varying signature and its steering vector Optimum SINR LOCSME [32] LCWC [9] SQP [6] OKSPME OKSPME−SG OKSPME−CCG OKSPME−MCG where s k (i)(k = 0, 1, 2, 3, 4) are i.i.d zero mean complex Gaussian random variables independently drawn from a random generator. The angles θ k (k = 0, 1, 2, 3, 4) are drawn independently in each simulation run from a uniform generator with mean 10 • and standard deviation 2 • . At this time, s k (i) changes both from run to run and from snapshot to snapshot. In order to show the effects caused by incoherent scattering only, we set the parameters M = 40 and K = 3, study the SINR versus SNR performance of the selected algorithms in Fig. 11 and compare the results with Fig. 10. As a result, a performance degradation is observed for all the studied algorithms. This is because the time-varying nature of incoherent scattering results in more dynamic and environmental Optimum SINR LOCSME [32] LCWC [9] SQP [6] OKSPME OKSPME−SG OKSPME−CCG OKSPME−MCG Optimum SINR OKSPME OKSPME−SG OKSPME−CCG OKSPME−MCG Beamspace [23] Fig. 9. Coherent local scattering, SINR versus snapshots, M = 40, K = 3 uncertainties in the system, which increases the steering vector mismatch.

VII. CONCLUSION
We have proposed the OKSPME algorithm based on the exploitation of cross-correlation mismatch estimation and the orthogonal Krylov subspace. In addition, low complexity RAB algorithms, OKSPME-SG, OKSPME-CCG and OKSPME-MCG have been developed to enable the beamforming weights to be updated recursively without matrix inversions. A detailed steering vector estimation MSE analysis for the general RAB design approach that relies on a presumed angular sector as prior knowledge has been provided. The computational complexity of the proposed and some of the existing algorithms have been compared and discussed. Simulation results have shown that the proposed algorithms have robustness against Optimum SINR OKSPME OKSPME−SG OKSPME−CCG OKSPME−MCG Beamspace [23]