Efficient Rank-Adaptive Least-Square Estimation and Multiple-Parameter Linear Regression Using Novel Dyadically Recursive Hermitian Matrix Inversion

Abstract

Least-square estimation (LSE) and multiple-parameter linear regression (MLR) are the important estimation techniques for engineering and science, especially in the mobile communications and signal processing applications. The majority of computational complexity incurred in LSE and MLR arises from a Hermitian matrix inversion. In practice, the Yule-Walker equations are not valid, and hence the Levinson-Durbin algorithm cannot be employed for general LSE and MLR problems. Therefore, the most efficient Hermitian matrix inversion method is based on the Cholesky factorization. In this paper, we derive a new dyadic recursion algorithm for sequential rank-adaptive Hermitian matrix inversions. In addition, we provide the theoretical computational complexity analyses to compare our new dyadic recursion scheme and the conventional Cholesky factorization. We can design a variable model-order LSE (MLR) using this proposed dyadic recursion approach thereupon. Through our complexity analyses and the Monte Carlo simulations, we show that our new dyadic recursion algorithm is more efficient than the conventional Cholesky factorization for the sequential rank-adaptive LSE (MLR) and the associated variable model-order LSE (MLR) can seek the trade-off between the targeted estimation performance and the required computational complexity. Our proposed new scheme can benefit future portable and mobile signal processing or communications devices.

1. Introduction

Least-square estimation (LSE) or multiple-parameter linear regression (MLR) is a crucial mathematical technique for science and engineering ubiquitously. A wide variety of applications in signal processing devices, mobile communications, computational physics, and statistics arise from the LSE or MLR. It has been adopted as a powerful tool in nuclear science [1, 2], system identification [3], data mining [4], linear system modeling [5], medical imaging analysis [6], and so forth. In particular, for signal processing and communications, least-square approach dominates a large majority of applications, such as the equalizations [7–16] as well as the interference estimation [17–19]. An overview of linear regression can be found in [20].

However, all of the aforementioned applications of LSE or MLR require the predetermined model order (the number of parameters to be estimated), and therefore the associated computational complexities are fixed and not flexible. In practice, the predetermined model order should be large enough to assure that the square error can be minimized due to the sufficient degree of freedom. Therefore, the model-order selection is essential to all LSE and MLR problems. Intuitively, we may adjust the model order such that an appropriate model order can be achieved to satisfy a preset square-error requirement. In the meantime, the computational complexity can be variable and depends on the error tolerance using such an adjustable model order. Nevertheless, the computational complexity would augment tremendously during the model order adjustment. How to make the most use of the calculated estimates from the previous model order and estimate the new parameters thereupon draws the researchers’ interest always.

Rank-one update scheme or iterative Levinson-Durbin algorithm is the well-known approach for computing the linear prediction coefficients with adjustable model orders [21]. The Levinson-Durbin algorithm is based on the Yule-Walker equations [22] and the strict assumption of Hermitian Toeplitz correlation matrix [21], but such an assumption is valid only when the multiple time series (the observed data) possess the exact temporal relationship (one is the delayed version of another) and the sample size approaches infinity under the ergodicity [23]. Besides, the Levinson-Durbin algorithm can only increase the model order by one at each iteration. A new class of extended Levinson-Durbin algorithms has been developed to achieve computational efficiency: in speech or audio coders, an adaptive multi-stage Levinson-Durbin algorithm was investigated in [24], and it was shown to be more robust than the conventional Levinson-Durbin algorithm when the input signal has large spectral dynamics; low-complexity approaches were also provided to determine the optimal delay and calculate the linear prediction filter coefficients in [7, 25]; iterative Levinson-Durbin algorithms were investigated for the frequency-domain decision feedback and MMSE equalizers in [8–12]. However, the Levinson-Durbin family of algorithms are still restricted by their limited model-order scalability due to their common rank-one update characteristics and Hermitian Toeplitz correlation matrix assumption.

For a general linear regression problem, the temporal relationship among multiple observed time series does not necessarily sustain. Hence, the Hermitian Toeplitz correlation matrix can exist only for some circumstances. To be more precise, the correlation matrix for a general linear regression problem should be Hermitian only [22, 26]. If the correlation matrix is Hermitian, the most efficient technique for the least-square solutions involves Cholesky factorization [21]. Recursion-based factorization methods were proposed to reduce the memory usage and computational complexity [27, 28]. However, they are not related to the direct Hermitian matrix inversion, which is the backbone of the linear estimation and the multiple-parameter linear regression. The existing Cholesky-factorization-based Hermitian matrix inversion procedure has to recalculate all the estimates when the model order is enlarged. In this paper, we present a novel rank-adaptive Cholesky factorization scheme, which can be employed to establish a low-complexity variable model-order LSE or MLR algorithm. In our newly derived LSE or MLR algorithm, the model order can be adjusted dyadically (model order is 2^𝑚, 𝑚=1,2,3,…) while that in the Levinson-Durbin algorithm can be incremented one by one. To the best of our knowledge, our new algorithm is the only estimation or regression method which can adjust the model order dyadically. The magnificent dyadic-scalability property of our proposed LSE and MLR algorithm can well serve for the future signal processing or communications applications which demand the adaptability in both memory storage and computational complexity.

The rest of this paper is organized as follows. In Section 2, the problem formulation for the LSE and MLR is presented where the Wiener-Hopf equations will result in a Hermitian correlation matrix inversion. In Section 3, we derive the dyadic recursion formula for such a Hermitian matrix inversion. In Section 3.2, we design a novel dyadic rank-adaptive Hermitian matrix inversion algorithm based on the results in Section 3.1. The computational complexity analyses for the fixed model-order and rank-adaptive Hermitian matrix inversion using the conventional Cholesky factorization and our new dyadic recursion are manifested in Section 4. The frequent communications and signal processing applications using the sequential rank-adaptive LSE thereupon for channel estimation and equalization are introduced in Section 5. The numerical evaluations for the comparisons of the computational complexities and the least mean-square errors between our new dyadic recursion scheme and the conventional method are delineated in Section 6. Ultimate concluding remarks will be drawn in Section 7.

Notations: 𝒩, 𝒞 denote the sets of positive integers and complex numbers, respectively. ^→𝑎 denotes a column vector and 𝐴 denotes a matrix; 𝐴^𝑇, 𝐴^𝐻 denote the transpose and the Hermitian adjoint of 𝐴, respectively. The symbols 𝑖, 𝑗, and 𝑘 are reserved as the integer indices. The symbol 𝐴_𝑘(𝑖,𝑗) represents a 𝑘×𝑘 square matrix with the dimensionality 𝑘, and this matrix is in a larger matrix which contains 𝐴_𝑘(𝑖,𝑗) as its (𝑖,𝑗)-block. ≡ represents the mathematical definition and ∗ denotes the complex-conjugate operator. The computational complexity function ℂ_•(𝑁) specifies the number of complex multiplications incurred when • of dimension 𝑁 is calculated.

2. Least-Square Estimation and Multiple-Parameter Linear Regression

The LSE and MLR can be formulated as follows: given the degree of freedom or model order 𝑁∈𝒩 and the 𝑁 observed time series 𝑥_𝑖(𝑛), 𝑛=0,1,…,𝑀−1 and the target time series (desired response) 𝑑(𝑛), 𝑛=0,1,…,𝑀−1 where 𝑀 is the sample size for all time series, we want to search for the best-fitting coefficients 𝑤_𝑖, for 𝑖=1,2,…,𝑁, such that ∑^𝑀−1_𝑛=0∑^𝑁_𝑖=1|𝑑(𝑛)−𝑤^∗_𝑖𝑥_𝑖(𝑛)|² is minimized. In general, 𝑥_𝑖(𝑛), 𝑑(𝑛), and 𝑤_𝑖 are complex-valued, ∀𝑛, ∀𝑖. Thereby, in matrix form, the LSE or MLR problem can be formulated as 𝑤^→_opt=argmin^→_𝑤^𝑀−1_𝑛=0𝑑(𝑛)−𝑤^→𝑋^𝐻→(𝑛)², (1) where 𝑤^→≡𝑤₁𝑤₂⋯𝑤_𝑁^𝑇,_→𝑋(𝑛)≡𝑥₁(𝑛)𝑥₂(𝑛)⋯𝑥_𝑁(𝑛)^𝑇, (2)(3) and 𝑤^→_opt is the collection of the best-fitting (optimal) coefficients. Equation (1) can be rewritten as 𝑤^→_opt=argmin𝑤_→𝑤^→𝐻𝑅𝑤_𝑁^→−^→𝑃𝑤^𝐻_𝑁^→−𝑤^→𝐻→𝑃_𝑁, (4) where 𝑅_𝑁≡1𝑀^𝑀−1→_𝑛=0𝑋(𝑛)𝑋^→𝐻(𝑛),_→𝑃_𝑁≡1𝑀^𝑀−1_𝑛=0𝑑^∗(𝑛)𝑋^→(𝑛). (5)(6) According to [22], the least-square solution to the LSE or MLR problem stated in (4) is given by the Wiener-Hopf equation as follows: 𝑤^→_opt=𝑅_𝑁^→−1𝑃_𝑁, (7) and the least mean-squared error 𝜎²_opt is given by 𝜎²_opt≡1𝑀^𝑀−1_𝑛=0𝑑(𝑛)−𝑤^→𝑋^𝐻_opt^→(𝑛)²=1𝑀^𝑀−1_𝑛=0|𝑑(𝑛)|²−^→𝑃𝑤^𝐻_𝑁^→_opt. (8) It is noted from (7) that the major portion of the computational complexity for solving 𝑤^→_opt involves the calculation of 𝑅_𝑁⁻¹. Therefore, the computational complexity can be further reduced if the special mathematical properties of the correlation matrix 𝑅_𝑁 are exploited.

In particular, when the LSE or MLR is adopted for filter or equalizer design, the temporal relationship 𝑥_𝑖(𝑛−1)=𝑥_𝑖+1(𝑛), for 𝑖=1,2,…,𝑁−1, for all 𝑛, is assumed. Thus, as the sample size is large (𝑀→∞), the correlation matrix 𝑅_𝑁 is both Hermitian and Toeplitz under the stationary and ergodic assumption [23]. However, very often, the multiple time series, 𝑥_𝑖(𝑛), 𝑖=1,2,…,𝑁, specify the different observations generated from independent sources, or the sample size 𝑀 is not large. Consequently, ^𝑀−1_𝑛=0𝑥_𝑖(𝑛)𝑥^∗_𝑖+𝑞(𝑛)≠^𝑀−1_𝑛=0𝑥_𝑗(𝑛)𝑥^∗_𝑗+𝑞(𝑛),∀𝑖≠𝑗,1≤𝑖,𝑗≤𝑁,−𝑁+1≤𝑞≤𝑁−1,1≤𝑖+𝑞,𝑗+𝑞≤𝑁. (9)(10) According to (9), the correlation matrix 𝑅_𝑁 is not necessarily Toeplitz. Based on (5), the only ever-sustaining property is that 𝑅_𝑁 is always Hermitian. Instead of employing the QR factorization to compute 𝑅_𝑁⁻¹ with high complexity [21], we can benefit from this Hermitian property to carry out. Since 𝑅_𝑁 is Hermitian, we can use the Cholesky factorization such that 𝑅_𝑁=𝐿_𝑁𝐿^𝐻_𝑁, (11) where 𝐿_𝑁 and 𝐿^𝐻_𝑁 are the lower- and upper-triangular matrices, respectively. If 𝑅_𝑁⁻¹ exists, according to (11), we can calculate 𝑅_𝑁⁻¹ as 𝑅_𝑁⁻¹=𝐿_𝑁^−𝐻𝐿_𝑁⁻¹, (12) where 𝐿_𝑁^−𝐻≡(𝐿_𝑁⁻¹)^𝐻=(𝐿^𝐻_𝑁)⁻¹, which is a more efficient procedure than the inversion based on the QR factorization [21].

3. Novel Dyadically Recursive Hermitian Matrix Inversion

3.1. Recursion of Hermitian Matrix Inversion Using Cholesky Outer Product

From the previous studies in Section 2, we cannot adopt the Levinson-Durbin algorithm for LSE or MLR when the correlation matrix 𝑅_𝑁 is not Toeplitz. To the best of our knowledge, no exiting literature had ever addressed any form of recursive formula for the procedure of Hermitian matrix inversion. In this section, we would like to introduce such a new procedure, which can dyadically extend the inverse of any Hermitian matrix from those of its submatrices with a half dimension. This novel dyadic recursion for Hermitian matrix inversion can reduce the computational complexity as we will analyze later on.

Any lower-triangular matrix 𝐿_𝑁 of dimension 𝑁 can be decomposed as 𝐿_𝑁=𝐿_𝑁/2(1,1)𝟎_𝑁/2𝜓_𝑁/2𝐿_𝑁/2(2,2), (13) where 𝟎_𝑁/2 is the (𝑁/2)×(𝑁/2) matrix containing all zeros, 𝐿_𝑁/2(1,1), 𝐿_𝑁/2(2,2) are the two (𝑁/2)×(𝑁/2) submatrices both of which are lower-triangular, and 𝜓_𝑁/2 is the (𝑁/2)×(𝑁/2) submatrix. According to (11) and (13), we have 𝑅_𝑁≡𝑅_𝑁/2(1,1)𝑅^𝐻_𝑁/2()𝑅2,1_𝑁/2(2,1)𝑅_𝑁/2()=𝐿2,2_𝑁/2(1,1)𝐿^𝐻_𝑁/2(1,1)𝐿^𝐻_𝑁/2(1,1)𝜓^𝐻_𝑁/2𝜓_𝑁/2𝐿^𝐻_𝑁/2(1,1)𝜓_𝑁/2𝜓^𝐻_𝑁/2+𝐿_𝑁/2(2,2)𝐿^𝐻_𝑁/2(2,2), (14)(15) where 𝑅_𝑁/2()≡11,1𝑀^𝑀−1_𝑛=0𝑥₁(𝑛)𝑥₂(𝑛)⋯𝑥_𝑁/2(𝑛)^𝑇×𝑥^∗₁(𝑛)𝑥^∗₂(𝑛)⋯𝑥^∗_𝑁/2(𝑛=𝐿_𝑁/2)(1,1)𝐿^𝐻_𝑁/2(),𝑅1,1_𝑁/2()≡12,1𝑀^𝑀−1_𝑛=0𝑥_(𝑁/2)+1(𝑛)𝑥_(𝑁/2)+2(𝑛)⋯𝑥_𝑁(𝑛)^𝑇×𝑥^∗₁(𝑛)𝑥^∗₂(𝑛)⋯𝑥^∗_𝑁/2(𝑛=𝜓_𝑁/2)𝐿^𝐻_𝑁/2(),𝑅1,1_𝑁/2()≡12,2𝑀^𝑀−1_𝑛=0𝑥_(𝑁/2)+1(𝑛)𝑥_(𝑁/2)+2(𝑛)⋯𝑥_𝑁(𝑛)^𝑇×𝑥^∗_(𝑁/2)+1(𝑛)𝑥^∗_(𝑁/2)+2(𝑛)⋯𝑥^∗_𝑁(𝑛=𝜓_𝑁/2)𝜓^𝐻_𝑁/2+𝐿_𝑁/2(2,2)𝐿^𝐻_𝑁/2(2,2). (16)(17)(18) According to (11)–(14), we obtain 𝐿_𝑁⁻¹=𝐿⁻¹_𝑁/2(1,1)𝟎_𝑁/2𝜙_𝑁/2𝐿⁻¹_𝑁/2(2,2), (19) where 𝜙_𝑁/2=−𝐿⁻¹_𝑁/2(2,2)𝜓_𝑁/2𝐿⁻¹_𝑁/2(1,1). (20)

3.2. New Efficient Dyadic Recursion for Hermitian Matrix Inversion

In this paper, we assume that the ill-posed problem can be neglected during the matrix inversion and no additional pivoting technique is required. Therefore, the computational complexity for the pivoting is not considered here. From the results in the previous subsection, we can design the following efficient dyadic recursion algorithm for any Hermitian matrix inversion.

Step 1 (Initialization). Set 𝑁=1 and 𝑅_𝑁(1,1)≡𝑟₀=(1/𝑀)∑^𝑀−1_𝑛=0|𝑥₁(𝑛)|². 𝐿_𝑁(1,1)≡𝑙₁=√𝑟₀. Hence, 𝐿_𝑁⁻¹(1,1)=1/𝑙₁.

Step 2. Construct the correlation matrix of dimension 2𝑁 such that 𝑅_2𝑁=𝑅_𝑁(1,1)𝑅^𝐻_𝑁()𝑅2,1_𝑁(2,1)𝑅_𝑁(2,2). (21)

Step 3 (Dyadic Expansion). Compute 𝜓_𝑁=𝑅_𝑁(2,1)𝐿_𝑁^−𝐻(1,1). (22) Use the Cholesky factorization to calculate 𝐿_𝑁(2,2)𝐿^𝐻_𝑁(2,2)=𝑅_𝑁(2,2)−𝜓_𝑁𝜓^𝐻_𝑁. (23) Then, calculate 𝐿_𝑁⁻¹(2,2) accordingly.
Compute 𝜙_𝑁=−𝐿_𝑁⁻¹(2,2)𝜓_𝑁𝐿_𝑁⁻¹(1,1). (24)

Step 4 (Dyadically Recursive Hermitian Matrix Inversion). Construct 𝐿_2𝑁=𝐿_𝑁(1,1)𝟎_𝑁𝜓_𝑁𝐿_𝑁(),𝐿2,2_−12𝑁=𝐿_𝑁⁻¹(1,1)𝟎_𝑁𝜙_𝑁𝐿_𝑁⁻¹().2,2 (25)

Step 5. Repeat Steps 2–4 for 𝑁=2^𝑝, 𝑝=1,2,3,… and stop at 𝑁=2^𝑝, which is determined by the preset square-error tolerance. Use the ultimate 𝐿⁻¹_2𝑁 to calculate 𝑅⁻¹_2𝑁=𝐿^−𝐻_2𝑁𝐿⁻¹_2𝑁.

4. Computational Complexity Analyses for Hermitian Matrix Inversion

In this section, we will study the computational complexities of the two algorithms for the Hermitian matrix inversion, namely (i) conventional Cholesky factorization and (ii) our new dyadic recursion presented in Section 3.2. Since the correlation matrix assumption for the Levinson-Durbin method is different and it also has to be based on the Yule-Walker equations as previously discussed in Section 2, we cannot compare the Levinson-Durbin method with the two underlying schemes here on a fair ground. Without loss of generality, we assume that all of the observed time series are normalized with respect to their energies such that ∑^𝑀−1_𝑛=0|𝑥_𝑖(𝑛)|²=1, for 𝑖=1,2,…,𝑁 in our subsequent analyses. Thus, the diagonal elements in the correlation matrices 𝑅_𝑁 of any dimension 𝑁 are always 1.

4.1. Computational Complexity for Hermitian Matrix Inversion Using Conventional Cholesky Factorization

The Hermitian matrix inversion using the conventional Cholesky factorization has to be based on the predetermined model order 𝑁. Since the diagonal elements of 𝑅_𝑁 are 1, it yields 𝐿_𝑁=𝑙10⋯⋯⋯0₂₁𝑙10⋯⋯0₃₁𝑙₃₂𝑙10⋯0⋮⋮⋯⋱⋱⋮⋮⋮⋯⋯⋱0_𝑁1𝑙_𝑁2𝑙_𝑁3⋯𝑙_{𝑁,(𝑁−1)}1,𝑅_𝑁=𝐿_𝑁𝐿^𝐻_𝑁. (26) According to (26), the corresponding computational complexity ℂ_𝐿𝐿^𝐻(𝑁) in terms of the total complex multiplications is ℂ_𝐿𝐿^𝐻(𝑁)=𝑁(𝑁−1)(𝑁−2)6. (27) To compute 𝐿_𝑁⁻¹, we can write 𝐿_𝑁⁻¹=𝜆10⋯⋯⋯0₂₁𝜆10⋯⋯0₃₁𝜆₃₂𝜆10⋯0⋮⋮⋯⋱⋱⋮⋮⋮⋯⋯⋱0_𝑁1𝜆_𝑁2𝜆_𝑁3⋯𝜆_{𝑁,(𝑁−1)}1, (28) where 𝜆_𝑖𝑗∈𝒞, for 1≤𝑗<𝑖≤𝑁, are the values to be determined from 𝐿_𝑁. Since 𝐿_𝑁⁻¹𝐿_𝑁=𝐼_𝑁̃ where 𝐼_𝑁̃ is the 𝑁×𝑁 identity matrix, the computational complexity ℂ𝐿⁻¹(𝑁) for calculating 𝐿_𝑁⁻¹ is also ℂ𝐿⁻¹(𝑁)=𝑁(𝑁−1)(𝑁−2)6. (29) Once 𝐿_𝑁⁻¹ is available, the inverse matrix 𝑅_𝑁⁻¹=𝐿^−𝐻(𝑁)𝐿⁻¹(𝑁) can be calculated accordingly and the corresponding computational complexity for this multiplication is ℂ𝐿^−𝐻𝐿⁻¹(𝑁)=𝑁(𝑁−1)(2𝑁−1)6. (30) Thus, from (27), (29), and (30), the total computational complexity ℂ𝑅⁻¹_Chol(𝑁) for calculating 𝑅_𝑁⁻¹ using the conventional Cholesky factorization is ℂ𝑅⁻¹_Chol(𝑁)=ℂ_𝐿𝐿^𝐻(𝑁)𝐿+ℂ⁻¹(𝑁)𝐿+ℂ^−𝐻𝐿⁻¹(𝑁)=4𝑁³−9𝑁²+5𝑁6. (31)

4.2. Computational Complexity for Hermitian Matrix Inversion Using New Efficient Dyadic Recursion

From the previous discussion in Section 3.2, we can calculate the computational complexity for our new scheme thereupon. First, let us focus on the incurred computational complexity in terms of complex multiplications when we apply the recursion in Section 3.2 to solve 𝑅₂^−1𝑝 from the given 𝐿₂^𝑝−1(1,1)=𝐿₂^𝑝−1 and 𝐿₂^−1𝑝−1(1,1)=𝐿₂^−1𝑝−1, for 𝑝∈𝒩. According to (22), the computational complexity ℂ_𝜓(2^𝑝−1) for computing 𝜓_2𝑝−1=𝑅₂^𝑝−1(2,1)𝐿₂^{−𝐻𝑝−1}(1,1) is ℂ_𝜓2^𝑝−1=4^𝑝−12^𝑝−1−12. (32) Next, the calculation of 𝜓_2𝑝−1𝜓^𝐻_2𝑝−1 involves the computational complexity ℂ_𝜓𝜓^𝐻2^𝑝−1=8^𝑝−1. (33) Then, we can compute 𝐿_2𝑝−1(2,2)𝐿^𝐻₂^𝑝−1(2,2) according to (23); the corresponding computational complexity, similar to (27), is ℂ𝐿_2,2𝐿^𝐻_2,22^𝑝−1=2^𝑝−12^𝑝−1−12^𝑝−1−26. (34) In addition, similar to (29), the computational complexity for calculating 𝐿⁻¹_2𝑝−1(2,2) is also ℂ𝐿⁻¹_2,22^𝑝−1=2^𝑝−12^𝑝−1−12^𝑝−1−26. (35) According to (24), we can calculate the computational complexity ℂ_𝜙(2^𝑝−1) for 𝜙_2𝑝−1=−𝐿₂^−1𝑝−1(2,2)𝜓_2𝑝−1𝐿₂^−1𝑝−1(1,1) such that ℂ_𝜙2^𝑝−1=4^𝑝−12^𝑝−1−1. (36) Thus, the computational complexity involved in the calculation of 𝐿₂^−1𝑝 is equivalent to ℂ𝐿⁻¹2^𝑝=ℂ_𝜓2^𝑝−1+ℂ_𝜓𝜓^𝐻2^𝑝−1𝐿+ℂ_2,2𝐿^𝐻_2,22^𝑝−1𝐿+ℂ⁻¹_2,22^𝑝−1+ℂ_𝜙2^𝑝−1. (37) Consequently, according to the steps in Section 3.2 and (30), as well as (32)–(37), the total computational complexity ℂ𝑅⁻¹_new(2^𝑝) for our new scheme in one dyadic recursion to solve 𝑅₂^−1𝑝 can be calculated as ℂ𝑅⁻¹_new2^𝑝𝐿=ℂ⁻¹2^𝑝+2^𝑝2^𝑝−12^𝑝+1−16=118^𝑝−184^𝑝+82^𝑝.16 (38)

4.3. Computational Complexities for Sequential Rank-Adaptive Hermitian Matrix Inversions

From the studies in Sections 4.1 and 4.2, we have found that the asymptotical complexity ratio lim_𝑝→∞ℂ𝑅⁻¹_Chol2^𝑝ℂ𝑅⁻¹_new2^𝑝=lim_𝑝→∞8348^𝑝−94^𝑝+52^𝑝118^𝑝−18(4^𝑝)+82^𝑝=3233≈96.97%, (39) where 𝑁=2^𝑝. It is noted that from (39), our new dyadic recursion scheme does not possess any advantage over the conventional Cholesky method. The reason is quite obvious. The ratio presented in (39) involves the fixed model order for the conventional Cholesky factorization. If the model order 𝑁=2^𝑚 is predetermined as a constant, our new scheme described in Section 3.2 still needs to calculate all the intermediate matrices 𝐿₂^𝑝 and 𝐿₂^−1𝑝, for all 𝑝=1,2,…,𝑚 while the conventional Cholesky factorization method does not.

However, in practice, the model order 𝑁=2^𝑝 should be variable and the sequential rank-adaptive LSE or MLR should take place. The least-square error will decrease when the model order 𝑁=2^𝑝 increases, that is, 𝑝=1,2,3,…. A stop criterion for the error tolerance (the maximum squared error which is allowed) can be introduced to terminate the model order enlargement. If such a sequential rank-adaptive LSE or MLR is considered, we need to recalculate the corresponding computational complexities for the sequential computations of 𝑅₂^−1𝑝, 𝑝=1,2,3,…. Assume that the stop criterion terminates the sequential computations of 𝑅₂^−1𝑝 at 𝑁=2^𝑚. We can compute the total computational complexity ℂ𝑅⁻¹_Chol,seq(2^𝑚) for such a sequential procedure using the conventional Cholesky factorization as ℂ𝑅⁻¹_Chol,seq2^𝑚=^𝑚_𝑝=1ℂ𝑅⁻¹_Chol2^𝑝=16821^𝑚−24^𝑚+532^𝑚−37. (40) On the other hand, if our new dyadic recursion scheme is used for the sequential calculations of 𝑅₂^−1𝑝 up to 𝑁=2^𝑚, the corresponding computational complexity is ℂ𝑅⁻¹_new,seq2^𝑚𝐿=ℂ⁻¹2^𝑚+^𝑚_𝑝=1138^𝑝−124^𝑝+162^𝑝=2478336^𝑚−31424^𝑚+232^𝑚−1,21 (41) where ℂ𝐿⁻¹(2^𝑚) is given by (37). Similarly, we may rederive the complexity ratio for the sequential calculations such that lim_𝑚→∞ℂ𝑅⁻¹_Chol,seq2^𝑚ℂ𝑅⁻¹_new,seq(2^𝑚)=lim_𝑚→∞(16/21)8^𝑚−24^𝑚+(5/3)2^𝑚−(3/7)(247/336)(8^𝑚)−(31/24)(4^𝑚)+(2/3)(2^𝑚)−()=1/21256247≈103.64%. (42)

5. Applications of Sequential Rank-Adaptive Hermitian Matrix Inversions for Channel Estimation and Equalization

Based on the sequential rank-adaptive Hermitian matrix inversion approach in Section 4.3, we can design a new communication channel estimation and equalization scheme with model order adaptability. Other similar applications can be easily extended using our efficient dyadic Hermitian matrix inversion algorithm in Section 3.2. Here, we introduce this proposed new channel estimation and equalization method using our dyadic recursion for adjusting the model order.

A basic transmission model for communication systems can be formulated as 𝑟(𝑛)=𝑠(𝑛)⊗ℎ(𝑛)+𝜉(𝑛),𝑛=0,1,…,𝑀^′−1, (43) where ⊗ denotes the linear convolution; 𝑠(𝑛), 𝑟(𝑛), ℎ(𝑛), and 𝜉(𝑛) represent the transmitted training sequence, the received signal, the channel impulse response, and the background noise, respectively; 𝑀^′ is the sample size. Without loss of generality, the windowing of observed input data is assumed to comply with the covariance method [22]. For the channel estimation problem, the multiple time series can be set as 𝑥_𝑖(𝑛)≡𝑠(𝑛−𝑖+1),𝑛=𝑁−𝑖,𝑁−𝑖+1,…,𝑀^′−𝑖,𝑖=1,2,…,𝑁, (44) where 𝑥₁(𝑛)≡𝑠(𝑛) and 𝑁 is the estimated channel filter length. Thereby, the desired response is the received signal such that 𝑑(𝑛)≡𝑟(𝑛+𝑁−1),𝑛=0,1,…,𝑀^′−𝑁. (45) On the other hand, for the equalization problem, the multiple time series are set as 𝑥_𝑖(𝑛)≡𝑟(𝑛−𝑖+1),𝑛=𝑁−𝑖,𝑁−𝑖+1,…,𝑀^′−𝑖,𝑖=1,2,…,𝑁, (46) where 𝑥₁(𝑛)≡𝑟(𝑛) and 𝑁 is the estimated equalizer length. The desired response is therefore the transmitted training signal instead such that 𝑑(𝑛)≡𝑠(𝑛+𝑁−1),𝑛=0,1,…,𝑀^′−𝑁. (47)

According to (43), we can acquire the data specified by (44) and (45) for channel estimation and that specified by (46) and (47) for equalization. Thus, setting 𝑀≡𝑀^′−𝑁+1, we can carry out 𝑅_𝑁, 𝑃_𝑁, 𝑤^→_opt, 𝜎²_opt according to (5)–(8). Note that the least mean-squared error 𝜎²_opt given by (8) is a function of both 𝑀 and 𝑁, and hence we can write 𝜎²_opt=𝜎²_opt(𝑀,𝑁). (48) For a fixed sample size 𝑀, we can increase the model order 𝑁. For a sufficiently large 𝑁, we may achieve the error floor such that 𝜎²_opt(𝑀,𝑁)≈𝜎²_opt(𝑀,𝑁+𝑘),∀𝑘∈𝒩. (49) We can design a simple threshold criterion to determine the minimum model order 𝑁_min for satisfying the condition in (49). For the dyadic model order adaptation, we propose a stop criterion 𝐽(𝑁)=𝜎²_opt(𝑀,𝑁)−𝜎²_opt(𝑀,2𝑁)𝜎²_opt().𝑀,𝑁 (50) Therefore, given a predetermined threshold 𝜏_𝐽, the minimum model order 𝑁_min for the dyadic adaptation is obtained as 𝑁_min=inf_{𝑁∈𝑆𝐽}𝑁, (51) where 𝑆_𝐽≡𝑁∈𝒩∶𝐽(𝑁)≤𝜏_𝐽, (52) and inf denotes the infimum operation [29]. Ultimately, we can apply (5)–(8) and (50)-(51) to build a sequential rank-adaptive LSE or MLR procedure which is capable of adapting the model order dyadically.

6. Numerical Evaluations

In this section, we would like to provide the numerical evaluations for the illustrations of the computational complexities associated with the two Hermitian matrix inversion algorithms (conventional Cholesky factorization and our new efficient dyadic recursion) and the LSE applications using the adaptive model order for channel equalization. Table 1 lists the comparison of the two methods for a few model orders. It is noted that when a Hermitian matrix inversion is carried out for a fixed model order, the conventional Cholesky factorization approach is more computationally efficient than our new dyadic recursion. However, in comparison of the total computational complexity for the sequential calculation of variable model order Hermitian matrix inversion, our new dyadic recursion is more efficient than the conventional Cholesky factorization when the same terminal model order is set, since the latter does not utilize any information from the previous calculation of the Hermitian inversion with a smaller model order.

Table 1 

Computational complexities of hermitian matrix inversion algorithms in comparison.

6.1. Comparison of Computational Complexities

Figures 1 and 2 illustrate the total computational complexities in terms of complex multiplications for the Hermitian matrix inversion using the conventional Cholesky factorization and our new dyadic recursion method; they depict the complexities versus the model order for the fixed model order (one-time LSE or MLR) and the variable model order (sequential rank-adaptive LSE or MLR). According to Figures 1 and 2, it is clear that our dyadic recursion algorithm is more efficient than the conventional Cholesky factorization for sequent LSE (MLR) and vice versa. The computational complexity margin ratios in between for the two situations are also depicted in Figure 3. The asymptotical complexity margin ratios of around ±3% derived in (39) and (42) can be observed therein.

Figure 1 

Computational complexity comparison (model order 𝑁=2𝑝 is fixed for one-time LSE or MLR) for Hermitian matrix inversion using conventional Cholesky factorization 𝑅(ℂ−1Chol(𝑁)) and our proposed new dyadic recursion 𝑅(ℂ−1new(𝑁)).

Figure 2 

Computational complexity comparison (accumulated complexity for the sequential rank-adaptive LSE or MLR with terminal model order 𝑁=2𝑚) for Hermitian matrix inversion using conventional Cholesky factorization 𝑅(ℂ−1Chol,seq(𝑁)) and our proposed new dyadic recursion 𝑅(ℂ−1new,seq(𝑁)).

Figure 3 

Computational complexity margin ratio versus model order 𝑁 (fixed model order: 𝑅(ℂ−1new𝑅(𝑁)−ℂ−1Chol𝑅(𝑁))/(ℂ−1Chol(𝑁))%, variable model order: 𝑅(ℂ−1new,seq𝑅(𝑁)−ℂ−1Chol,seq𝑅(𝑁))/(ℂ−1Chol,seq(𝑁))%, equality line is for the margin ratio = 0%; 𝑁 is the fixed or the terminal model order therein).

6.2. Sequential Rank-Adaptive LSE for Channel Equalization

Next, we would like to present some simulations for illustrating the application of our proposed sequential rank-adaptive Hermitian matrix inversions for channel equalization as discussed in Section 5. Assume that the modulation type of the transmitted signal is BPSK and the channel noise is additive white Gaussian noise (AWGN). The channel transfer function is arbitrarily chosen as 𝐻(𝑧)=0.58+0.56𝑧⁻¹+0.55𝑧⁻²−0.18𝑧⁻³−0.12𝑧⁻⁴, where the channel gain is normalized as unity such that ∑_𝑛|ℎ(𝑛)|²=1, and it complies with the LOS (light-of-sight) channel model such that the leading coefficient is the largest among all in magnitude. The received signals are generated by the computer using (43). The signal sample size is 𝑀=500.

We carry out 200 Monte Carlo trials with randomly generated AWGN (various signal-to-noise ratios from 0 to 30 dB) and the aforementioned transmission model. The least mean-squared errors 𝜎²_opt on average are calculated using (8). Here, we compare three schemes, namely (i) one-time Hermitian matrix inversion using the Cholesky factorization with the fixed model order 𝑁=12 (specified as “fixed 𝑁=12”), (ii) one-time Hermitian matrix inversion using the Cholesky factorization with the fixed model order 𝑁=25 (specified as “fixed 𝑁=25”), and (iii) sequential rank-adaptive Hermitian matrix inversion using our proposed dyadic recursion with variable model order (specified as “Variable Model Order”). Two different threshold values are employed, namely 𝜏_𝐽=1/200,1/30. It is noted that the allowed maximum model order for our dyadic recursion scheme is defaulted as ⌊(𝑁−𝑀+1)/10⌋, where ⌊⌋ is the flooring operator.

Figures 4–7 show the simulation results. According to Figures 4 and 5, our proposed new dyadic recursion method achieves the smallest least mean-squared errors most of time especially for low threshold values. The trade-off of our scheme is the increased average computational complexities, which are shown in Figures 6 and 7. From Figures 6 and 7, it can be observed that we may save a huge computational burden for a bad channel condition (low signal-to-noise ratio) for a very slight equalization quality deterioration. According to Figures 1–7, we can justify that our dyadic recursion algorithm is more efficient than the conventional Cholesky factorization in the Hermitian matrix inversion for the sequential rank-adaptive LSE or MLR. Besides, the sequential rank-adaptive LSE or MLR with variable model orders using our recursion scheme can also seek the trade-off between the estimation performance and the required computational complexity.

Figure 4 

Various least mean-squared errors 𝜎2opt versus signal-to-noise ratio (the threshold for our variable model scheme is 𝜏𝐽=1/200).

Figure 5 

Various least mean-squared errors 𝜎2opt versus signal-to-noise ratio (the threshold for our variable model scheme is 𝜏𝐽=1/30).

Figure 6 

The corresponding computational complexities to Figure 4 versus signal-to-noise ratio (the threshold for our variable model scheme is 𝜏𝐽=1/200).

Figure 7 

The corresponding computational complexities to Figure 5 versus signal-to-noise ratio (the threshold for our variable model scheme is 𝜏𝐽=1/30).

7. Conclusion

Hermitian matrix inversion is the pivotal computation for least-square estimation and multiple-parameter linear regression. Conventional Cholesky factorization can be applied for one-time LSE and MLR with a predetermined fixed model order. In this paper, we decompose any Hermitian matrix into submatrices with a half of the original dimension and derive the new dyadic recursion algorithm for the Hermitian matrix inversion accordingly. Moreover, for the theoretical comparison, we elaborate the computational complexity analyses to derive the asymptotical complexity ratios for the Hermitian matrix inversion using the conventional Cholesky factorization and our new dyadic recursion scheme. We show that our new method can achieve the complexity margin ratio of −3.03% for the sequential rank-adaptive Hermitian matrix inversions and that of +3.64% for the one-time Hermitian matrix inversion over the conventional Cholesky factorization method. We also present the applications of the sequential rank-adaptive LSE with variable model order using our new dyadic recursion procedure for channel estimation and equalization in telecommunications. Our proposed new variable model order LSE scheme can seek the trade-off between the estimation performance and the computational complexity.

Acknowledgments

This work was supported by Information Technology Research Award for National Priorities (NSF-ECS 0426644) from National Science Foundation, Faculty Research Initiation Award from the Southeastern Center for Electrical Engineering Education, Research Enhancement Award from the Louisiana-NASA Space Consortium and Faculty Research Award from Louisiana State University.