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

Authors: Hsiao-Chun Wu; Shih Yu Chang; Tho Le-Ngoc; Yiyan Wu
Journal: International Journal of Antennas and Propagation
 (2012)

---

## 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.


---

## Body

## 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 predeterminedmodel 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 2m, m=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 Section2, 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. a→ denotes a column vector and A~ denotes a matrix; A~T, A~H denote the transpose and the Hermitian adjoint of A~, respectively. The symbols i, j, and k are reserved as the integer indices. The symbol A~k(i,j) represents a k×k square matrix with the dimensionality k, and this matrix is in a larger matrix which contains A~k(i,j) as its (i,j)-block. ≡ represents the mathematical definition and * denotes the complex-conjugate operator. The computational complexity function ℂ•(N) specifies the number of complex multiplications incurred when • of dimension N is calculated.


## 2. Least-Square Estimation andMultiple-Parameter Linear Regression

The LSE and MLR can be formulated as follows: given the degree of freedom or model orderN∈𝒩 and the N observed time series xi(n), n=0,1,…,M-1 and the target time series (desired response) d(n), n=0,1,…,M-1 where M is the sample size for all time series, we want to search for the best-fitting coefficients wi, for i=1,2,…,N, such that ∑n=0M-1∑i=1N|d(n)-wi*xi(n)|2 is minimized. In general, xi(n), d(n), and wi are complex-valued, ∀n, ∀i. Thereby, in matrix form, the LSE or MLR problem can be formulated as
(1)w→opt=argminw→(∑n=0M-1‍|d(n)-w→HX→(n)|2),
where
(2)w→≡[w1w2⋯wN]T,(3)X→(n)≡[x1(n)x2(n)⋯xN(n)]T,
and w→opt is the collection of the best-fitting (optimal) coefficients. Equation (1) can be rewritten as
(4)w→opt=argminw→(w→HR~Nw→-P→NHw→-w→HP→N),
where
(5)R~N≡1M∑n=0M-1‍X→(n)X→H(n),(6)P→N≡1M∑n=0M-1‍d*(n)X→(n).
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:
(7)w→opt=R~N-1P→N,
and the least mean-squared error σopt2 is given by
(8)σopt2≡1M∑n=0M-1‍|d(n)-w→optHX→(n)|2=1M∑n=0M-1‍|d(n)|2-P→NHw→opt.
It is noted from  (7) that the major portion of the computational complexity for solving w→opt involves the calculation of R~N-1. Therefore, the computational complexity can be further reduced if the special mathematical properties of the correlation matrixR~N are exploited.In particular, when the LSE or MLR is adopted for filter or equalizer design, the temporal relationshipxi(n-1)=xi+1(n), for i=1,2,…,N-1, for all n, is assumed. Thus, as the sample size is large (M→∞), the correlation matrix R~N is both Hermitian and Toeplitz under the stationary and ergodic assumption [23]. However, very often, the multiple time series, xi(n), i=1,2,…,N, specify the different observations generated from independent sources, or the sample size M is not large. Consequently,
(9)∑n=0M-1‍xi(n)xi+q*(n)≠∑n=0M-1‍xj(n)xj+q*(n),∀i≠j,1≤i,j≤N,-N+1≤q≤N-1,(10)1≤i+q,j+q≤N.
According to (9), the correlation matrix R~N is not necessarily Toeplitz. Based on (5), the only ever-sustaining property is that R~N is always Hermitian. Instead of employing the QR factorization to compute R~N-1 with high complexity [21], we can benefit from this Hermitian property to carry out. Since R~N is Hermitian, we can use the Cholesky factorization such that
(11)R~N=L~NL~NH,
where L~N and L~NH are the lower- and upper-triangular matrices, respectively. If R~N-1 exists, according to  (11), we can calculate R~N-1 as
(12)R~N-1=L~N-HL~N-1,
where L~N-H≡(L~N-1)H=(L~NH)-1, which is a more efficient procedure than the inversion based on the QR factorization [21].


## 3. Novel Dyadically RecursiveHermitian Matrix Inversion

### 3.1. Recursion of Hermitian Matrix Inversion Using Cholesky Outer Product

From the previous studies in Section2, we cannot adopt the Levinson-Durbin algorithm for LSE or MLR when the correlation matrix R~N 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 matrixL~N of dimension N can be decomposed as
(13)L~N=[L~N/2(1,1)0~N/2ψ~N/2L~N/2(2,2)],
where 0~N/2 is the (N/2)×(N/2) matrix containing all zeros, L~N/2(1,1), L~N/2(2,2) are the two (N/2)×(N/2) submatrices both of which are lower-triangular, and ψ~N/2 is the (N/2)×(N/2) submatrix. According to (11) and (13), we have
(14)R~N≡[R~N/2(1,1)R~N/2H(2,1)R~N/2(2,1)R~N/2(2,2)](15)=[L~N/2(1,1)L~N/2H(1,1)L~N/2H(1,1)ψ~N/2Hψ~N/2L~N/2H(1,1)ψ~N/2ψ~N/2H+L~N/2(2,2)L~N/2H(2,2)],
where
(16)R~N/2(1,1)≡1M∑n=0M-1‍([x1(n)x2(n)⋯xN/2(n)]T×[x1*(n)x2*(n)⋯xN/2*(n)])=L~N/2(1,1)L~N/2H(1,1),(17)R~N/2(2,1)≡1M∑n=0M-1‍([x(N/2)+1(n)x(N/2)+2(n)⋯xN(n)]T×[x1*(n)x2*(n)⋯xN/2*(n)])=ψ~N/2L~N/2H(1,1),(18)R~N/2(2,2)≡1M∑n=0M-1‍([x(N/2)+1(n)x(N/2)+2(n)⋯xN(n)]T×[x(N/2)+1*(n)x(N/2)+2*(n)⋯xN*(n)])=ψ~N/2ψ~N/2H+L~N/2(2,2)L~N/2H(2,2).
According to (11)–(14), we obtain
(19)L~N-1=[L~N/2-1(1,1)0~N/2ϕ~N/2L~N/2-1(2,2)],
where
(20)ϕ~N/2=-L~N/2-1(2,2)ψ~N/2L~N/2-1(1,1).


### 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 followingefficient dyadic recursion algorithm for any Hermitian matrix inversion.Step 1 (Initialization).
SetN=1 and R~N(1,1)≡r0=(1/M)∑n=0M-1|x1(n)|2. L~N(1,1)≡l1=r0. Hence, L~N-1(1,1)=1/l1.Step 2.
Construct the correlation matrix of dimension2N such that
(21)R~2N=[R~N(1,1)R~NH(2,1)R~N(2,1)R~N(2,2)].Step 3 (Dyadic Expansion).
Compute(22)ψ~N=R~N(2,1)L~N-H(1,1).
Use the Cholesky factorization to calculate
(23)L~N(2,2)L~NH(2,2)=R~N(2,2)-ψ~Nψ~NH.
Then, calculate L~N-1(2,2) accordingly.
Compute(24)ϕ~N=-L~N-1(2,2)ψ~NL~N-1(1,1).Step 4 (Dyadically Recursive Hermitian Matrix Inversion).
Construct(25)L~2N=[L~N(1,1)0~Nψ~NL~N(2,2)],L~2N-1=[L~N-1(1,1)0~Nϕ~NL~N-1(2,2)].Step 5.
Repeat Steps2–4 for N=2p, p=1,2,3,… and stop at N=2p, which is determined by the preset square-error tolerance. Use the ultimate L~2N-1 to calculate R~2N-1=L~2N-HL~2N-1.


## 3.1. Recursion of Hermitian Matrix Inversion Using Cholesky Outer Product

From the previous studies in Section2, we cannot adopt the Levinson-Durbin algorithm for LSE or MLR when the correlation matrix R~N 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 matrixL~N of dimension N can be decomposed as
(13)L~N=[L~N/2(1,1)0~N/2ψ~N/2L~N/2(2,2)],
where 0~N/2 is the (N/2)×(N/2) matrix containing all zeros, L~N/2(1,1), L~N/2(2,2) are the two (N/2)×(N/2) submatrices both of which are lower-triangular, and ψ~N/2 is the (N/2)×(N/2) submatrix. According to (11) and (13), we have
(14)R~N≡[R~N/2(1,1)R~N/2H(2,1)R~N/2(2,1)R~N/2(2,2)](15)=[L~N/2(1,1)L~N/2H(1,1)L~N/2H(1,1)ψ~N/2Hψ~N/2L~N/2H(1,1)ψ~N/2ψ~N/2H+L~N/2(2,2)L~N/2H(2,2)],
where
(16)R~N/2(1,1)≡1M∑n=0M-1‍([x1(n)x2(n)⋯xN/2(n)]T×[x1*(n)x2*(n)⋯xN/2*(n)])=L~N/2(1,1)L~N/2H(1,1),(17)R~N/2(2,1)≡1M∑n=0M-1‍([x(N/2)+1(n)x(N/2)+2(n)⋯xN(n)]T×[x1*(n)x2*(n)⋯xN/2*(n)])=ψ~N/2L~N/2H(1,1),(18)R~N/2(2,2)≡1M∑n=0M-1‍([x(N/2)+1(n)x(N/2)+2(n)⋯xN(n)]T×[x(N/2)+1*(n)x(N/2)+2*(n)⋯xN*(n)])=ψ~N/2ψ~N/2H+L~N/2(2,2)L~N/2H(2,2).
According to (11)–(14), we obtain
(19)L~N-1=[L~N/2-1(1,1)0~N/2ϕ~N/2L~N/2-1(2,2)],
where
(20)ϕ~N/2=-L~N/2-1(2,2)ψ~N/2L~N/2-1(1,1).


## 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 followingefficient dyadic recursion algorithm for any Hermitian matrix inversion.Step 1 (Initialization).
SetN=1 and R~N(1,1)≡r0=(1/M)∑n=0M-1|x1(n)|2. L~N(1,1)≡l1=r0. Hence, L~N-1(1,1)=1/l1.Step 2.
Construct the correlation matrix of dimension2N such that
(21)R~2N=[R~N(1,1)R~NH(2,1)R~N(2,1)R~N(2,2)].Step 3 (Dyadic Expansion).
Compute(22)ψ~N=R~N(2,1)L~N-H(1,1).
Use the Cholesky factorization to calculate
(23)L~N(2,2)L~NH(2,2)=R~N(2,2)-ψ~Nψ~NH.
Then, calculate L~N-1(2,2) accordingly.
Compute(24)ϕ~N=-L~N-1(2,2)ψ~NL~N-1(1,1).Step 4 (Dyadically Recursive Hermitian Matrix Inversion).
Construct(25)L~2N=[L~N(1,1)0~Nψ~NL~N(2,2)],L~2N-1=[L~N-1(1,1)0~Nϕ~NL~N-1(2,2)].Step 5.
Repeat Steps2–4 for N=2p, p=1,2,3,… and stop at N=2p, which is determined by the preset square-error tolerance. Use the ultimate L~2N-1 to calculate R~2N-1=L~2N-HL~2N-1.


## 4. Computational Complexity Analyses forHermitian 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 Section3.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 ∑n=0M-1‍|xi(n)|2=1, for i=1,2,…,N in our subsequent analyses. Thus, the diagonal elements in the correlation matrices R~N of any dimension N 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 orderN. Since the diagonal elements of R~N are 1, it yields
(26)L~N=[10⋯⋯⋯0l2110⋯⋯0l31l3210⋯0⋮⋮⋯⋱⋱⋮⋮⋮⋯⋯⋱0lN1lN2lN3⋯lN,(N-1)1],R~N=L~NL~NH.
According to (26), the corresponding computational complexity ℂL~L~H(N) in terms of the total complex multiplications is
(27)ℂL~L~H(N)=N(N-1)(N-2)6.
To compute L~N-1, we can write
(28)L~N-1=[10⋯⋯⋯0λ2110⋯⋯0λ31λ3210⋯0⋮⋮⋯⋱⋱⋮⋮⋮⋯⋯⋱0λN1λN2λN3⋯λN,(N-1)1],
where λij∈𝒞, for 1≤j<i≤N, are the values to be determined from L~N. Since L~N-1L~N=I~N where I~N is the N×N identity matrix, the computational complexity ℂL~-1(N) for calculating L~N-1 is also
(29)ℂL~-1(N)=N(N-1)(N-2)6.
Once L~N-1 is available, the inverse matrix R~N-1=L~-H(N)L~-1(N) can be calculated accordingly and the corresponding computational complexity for this multiplication is
(30)ℂL~-HL~-1(N)=N(N-1)(2N-1)6.
Thus, from (27), (29), and (30), the total computational complexityℂR~Chol-1(N) for calculating R~N-1 using the conventional Cholesky factorization is
(31)ℂR~Chol-1(N)=ℂL~L~H(N)+ℂL~-1(N)+ℂL~-HL~-1(N)=4N3-9N2+5N6.


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

From the previous discussion in Section3.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 R~2p-1 from the given L~2p-1(1,1)=L~2p-1 and L~2p-1-1(1,1)=L~2p-1-1, for p∈𝒩. According to  (22), the computational complexity ℂψ~(2p-1) for computing ψ~2p-1=R~2p-1(2,1)L~2p-1-H(1,1) is
(32)ℂψ~(2p-1)=4p-1(2p-1-1)2.
Next, the calculation of ψ~2p-1ψ~2p-1H involves the computational complexity
(33)ℂψ~ψ~H(2p-1)=8p-1.
Then, we can compute L~2p-1(2,2)L~2p-1H(2,2) according to  (23); the corresponding computational complexity, similar to  (27), is
(34)ℂL~2,2L~2,2H(2p-1)=2p-1(2p-1-1)(2p-1-2)6.
In addition, similar to  (29), the computational complexity for calculating L~2p-1-1(2,2) is also
(35)ℂL~2,2-1(2p-1)=2p-1(2p-1-1)(2p-1-2)6.
According to  (24), we can calculate the computational complexity ℂϕ~(2p-1) for ϕ~2p-1=-L~2p-1-1(2,2)ψ~2p-1L~2p-1-1(1,1) such that
(36)ℂϕ~(2p-1)=4p-1(2p-1-1).
Thus, the computational complexity involved in the calculation of L~2p-1 is equivalent to
(37)ℂL~-1(2p)=ℂψ~(2p-1)+ℂψ~ψ~H(2p-1)+ℂL~2,2L~2,2H(2p-1)+ℂL~2,2-1(2p-1)+ℂϕ~(2p-1).
Consequently, according to the steps in Section 3.2 and (30), as well as (32)–(37), the total computational complexityℂR~new-1(2p) for our new scheme in one dyadic recursion to solve R~2p-1 can be calculated as
(38)ℂR~new-1(2p)=ℂL~-1(2p)+2p(2p-1)(2p+1-1)6=11(8p)-18(4p)+8(2p)16.


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

From the studies in Sections4.1 and 4.2, we have found that the asymptotical complexity ratio
(39)limp→∞ℂR~Chol-1(2p)ℂR~new-1(2p)=limp→∞[(83)4(8p)-9(4p)+5(2p)11(8p)-18(4p)+8(2p)]=3233≈96.97%,
where N=2p. 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 N=2m is predetermined as a constant, our new scheme described in Section 3.2 still needs to calculate all the intermediate matrices L~2p and L~2p-1, for all p=1,2,…,m while the conventional Cholesky factorization method does not.However, in practice, the model orderN=2p should be variable and the sequential rank-adaptive LSE or MLR should take place. The least-square error will decrease when the model order N=2p increases, that is, p=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 R~2p-1, p=1,2,3,…. Assume that the stop criterion terminates the sequential computations of R~2p-1 at N=2m. We can compute the total computational complexity ℂR~Chol,seq-1(2m) for such a sequential procedure using the conventional Cholesky factorization as
(40)ℂR~Chol,seq-1(2m)=∑p=1m‍ℂR~Chol-1(2p)=1621(8m)-2(4m)+53(2m)-37.
On the other hand, if our new dyadic recursion scheme is used for the sequential calculations of R~2p-1 up to N=2m, the corresponding computational complexity is
(41)ℂR~new,seq-1(2m)=ℂL~-1(2m)+∑p=1m‍13(8p)-12(4p)+16(2p)=247336(8m)-3124(4m)+23(2m)-121,
where ℂL~-1(2m) is given by  (37). Similarly, we may rederive the complexity ratio for the sequential calculations such that
(42)limm→∞ℂR~Chol,seq-1(2m)ℂR~new,seq-1(2m)=limm→∞(16/21)(8m)-2(4m)+(5/3)(2m)-(3/7)(247/336)(8m)-(31/24)(4m)+(2/3)(2m)-(1/21)=256247≈103.64%.


## 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 orderN. Since the diagonal elements of R~N are 1, it yields
(26)L~N=[10⋯⋯⋯0l2110⋯⋯0l31l3210⋯0⋮⋮⋯⋱⋱⋮⋮⋮⋯⋯⋱0lN1lN2lN3⋯lN,(N-1)1],R~N=L~NL~NH.
According to (26), the corresponding computational complexity ℂL~L~H(N) in terms of the total complex multiplications is
(27)ℂL~L~H(N)=N(N-1)(N-2)6.
To compute L~N-1, we can write
(28)L~N-1=[10⋯⋯⋯0λ2110⋯⋯0λ31λ3210⋯0⋮⋮⋯⋱⋱⋮⋮⋮⋯⋯⋱0λN1λN2λN3⋯λN,(N-1)1],
where λij∈𝒞, for 1≤j<i≤N, are the values to be determined from L~N. Since L~N-1L~N=I~N where I~N is the N×N identity matrix, the computational complexity ℂL~-1(N) for calculating L~N-1 is also
(29)ℂL~-1(N)=N(N-1)(N-2)6.
Once L~N-1 is available, the inverse matrix R~N-1=L~-H(N)L~-1(N) can be calculated accordingly and the corresponding computational complexity for this multiplication is
(30)ℂL~-HL~-1(N)=N(N-1)(2N-1)6.
Thus, from (27), (29), and (30), the total computational complexityℂR~Chol-1(N) for calculating R~N-1 using the conventional Cholesky factorization is
(31)ℂR~Chol-1(N)=ℂL~L~H(N)+ℂL~-1(N)+ℂL~-HL~-1(N)=4N3-9N2+5N6.


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

From the previous discussion in Section3.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 R~2p-1 from the given L~2p-1(1,1)=L~2p-1 and L~2p-1-1(1,1)=L~2p-1-1, for p∈𝒩. According to  (22), the computational complexity ℂψ~(2p-1) for computing ψ~2p-1=R~2p-1(2,1)L~2p-1-H(1,1) is
(32)ℂψ~(2p-1)=4p-1(2p-1-1)2.
Next, the calculation of ψ~2p-1ψ~2p-1H involves the computational complexity
(33)ℂψ~ψ~H(2p-1)=8p-1.
Then, we can compute L~2p-1(2,2)L~2p-1H(2,2) according to  (23); the corresponding computational complexity, similar to  (27), is
(34)ℂL~2,2L~2,2H(2p-1)=2p-1(2p-1-1)(2p-1-2)6.
In addition, similar to  (29), the computational complexity for calculating L~2p-1-1(2,2) is also
(35)ℂL~2,2-1(2p-1)=2p-1(2p-1-1)(2p-1-2)6.
According to  (24), we can calculate the computational complexity ℂϕ~(2p-1) for ϕ~2p-1=-L~2p-1-1(2,2)ψ~2p-1L~2p-1-1(1,1) such that
(36)ℂϕ~(2p-1)=4p-1(2p-1-1).
Thus, the computational complexity involved in the calculation of L~2p-1 is equivalent to
(37)ℂL~-1(2p)=ℂψ~(2p-1)+ℂψ~ψ~H(2p-1)+ℂL~2,2L~2,2H(2p-1)+ℂL~2,2-1(2p-1)+ℂϕ~(2p-1).
Consequently, according to the steps in Section 3.2 and (30), as well as (32)–(37), the total computational complexityℂR~new-1(2p) for our new scheme in one dyadic recursion to solve R~2p-1 can be calculated as
(38)ℂR~new-1(2p)=ℂL~-1(2p)+2p(2p-1)(2p+1-1)6=11(8p)-18(4p)+8(2p)16.


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

From the studies in Sections4.1 and 4.2, we have found that the asymptotical complexity ratio
(39)limp→∞ℂR~Chol-1(2p)ℂR~new-1(2p)=limp→∞[(83)4(8p)-9(4p)+5(2p)11(8p)-18(4p)+8(2p)]=3233≈96.97%,
where N=2p. 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 N=2m is predetermined as a constant, our new scheme described in Section 3.2 still needs to calculate all the intermediate matrices L~2p and L~2p-1, for all p=1,2,…,m while the conventional Cholesky factorization method does not.However, in practice, the model orderN=2p should be variable and the sequential rank-adaptive LSE or MLR should take place. The least-square error will decrease when the model order N=2p increases, that is, p=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 R~2p-1, p=1,2,3,…. Assume that the stop criterion terminates the sequential computations of R~2p-1 at N=2m. We can compute the total computational complexity ℂR~Chol,seq-1(2m) for such a sequential procedure using the conventional Cholesky factorization as
(40)ℂR~Chol,seq-1(2m)=∑p=1m‍ℂR~Chol-1(2p)=1621(8m)-2(4m)+53(2m)-37.
On the other hand, if our new dyadic recursion scheme is used for the sequential calculations of R~2p-1 up to N=2m, the corresponding computational complexity is
(41)ℂR~new,seq-1(2m)=ℂL~-1(2m)+∑p=1m‍13(8p)-12(4p)+16(2p)=247336(8m)-3124(4m)+23(2m)-121,
where ℂL~-1(2m) is given by  (37). Similarly, we may rederive the complexity ratio for the sequential calculations such that
(42)limm→∞ℂR~Chol,seq-1(2m)ℂR~new,seq-1(2m)=limm→∞(16/21)(8m)-2(4m)+(5/3)(2m)-(3/7)(247/336)(8m)-(31/24)(4m)+(2/3)(2m)-(1/21)=256247≈103.64%.


## 5. Applications of Sequential Rank-Adaptive Hermitian Matrix Inversions forChannel Estimation and Equalization

Based on the sequential rank-adaptive Hermitian matrix inversion approach in Section4.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(43)r(n)=s(n)⊗h(n)+ξ(n),n=0,1,…,M'-1,
where ⊗ denotes the linear convolution; s(n), r(n), h(n), and ξ(n) represent the transmitted training sequence, the received signal, the channel impulse response, and the background noise, respectively; M′ 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
(44)xi(n)≡s(n-i+1),n=N-i,N-i+1,…,M′-i,i=1,2,…,N,
where x1(n)≡s(n) and N is the estimated channel filter length. Thereby, the desired response is the received signal such that
(45)d(n)≡r(n+N-1),n=0,1,…,M′-N.
On the other hand, for the equalization problem, the multiple time series are set as
(46)xi(n)≡r(n-i+1),n=N-i,N-i+1,…,M′-i,i=1,2,…,N,
where x1(n)≡r(n) and N is the estimated equalizer length. The desired response is therefore the transmitted training signal instead such that
(47)d(n)≡s(n+N-1),n=0,1,…,M′-N.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 M≡M′-N+1, we can carry out R~N, P~N, w→opt, σopt2 according to (5)–(8). Note that the least mean-squared error σopt2 given by  (8) is a function of both M and N, and hence we can write
(48)σopt2=σopt2(M,N).
For a fixed sample size M, we can increase the model order N. For a sufficiently large N, we may achieve the error floor such that
(49)σopt2(M,N)≈σopt2(M,N+k),∀k∈𝒩.
We can design a simple threshold criterion to determine the minimum model order Nmin for satisfying the condition in (49). For the dyadic model order adaptation, we propose a stop criterion
(50)J(N)=|σopt2(M,N)-σopt2(M,2N)|σopt2(M,N).
Therefore, given a predetermined threshold τJ, the minimum model order Nmin for the dyadic adaptation is obtained as
(51)Nmin=infN∈SJN,
where
(52)SJ≡{N∈𝒩:J(N)≤τJ},
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. Table1 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.
ℂ
R
~
chol
-
1
(
N
)
,
ℂ
R
~
new
-
1
(
N
)
,
ℂ
R
~
chol, seq
-
1
(
N
)
,
ℂ
R
~
new, seq
-
1
(
N
)
N
=
2
1
N
=
2
2
N
=
2
4
N
=
2
6
Conventional cholesky factorization (fixed model order)
1
22
2360
168672
Dyadic recursion (fixed model order)
2
28
2536
175648
Conventional cholesky factorization (variable model order)
1
23
2635
191643
Dyadic recursion (variable model order)
2
30
2691
187459


### 6.1. Comparison of Computational Complexities

Figures1 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 orderN=2p is fixed for one-time LSE or MLR) for Hermitian matrix inversion using conventional Cholesky factorization (ℂR~Chol-1(N)) and our proposed new dyadic recursion (ℂR~new-1(N)).Figure 2
Computational complexity comparison (accumulated complexity for the sequential rank-adaptive LSE or MLR with terminal model orderN=2m) for Hermitian matrix inversion using conventional Cholesky factorization (ℂR~Chol,seq-1(N)) and our proposed new dyadic recursion (ℂR~new,seq-1(N)).Figure 3
Computational complexity margin ratio versus model orderN (fixed model order: (ℂR~new-1(N)-ℂR~Chol-1(N))/(ℂR~Chol-1(N))%, variable model order: (ℂR~new,seq-1(N)-ℂR~Chol,seq-1(N))/(ℂR~Chol,seq-1(N))%, equality line is for the margin ratio = 0%; N 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 Section5. 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 H(z)=0.58+0.56z-1+0.55z-2-0.18z-3-0.12z-4, where the channel gain is normalized as unity such that ∑n‍|h(n)|2=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 M=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σopt2 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 N=12 (specified as “fixed N=12”), (ii) one-time Hermitian matrix inversion using the Cholesky factorization with the fixed model order N=25 (specified as “fixed N=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 τJ=1/200,1/30. It is noted that the allowed maximum model order for our dyadic recursion scheme is defaulted as ⌊(N-M+1)/10⌋, where ⌊⌋ is the flooring operator.Figures4–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σopt2 versus signal-to-noise ratio (the threshold for our variable model scheme is τJ=1/200).Figure 5
Various least mean-squared errorsσopt2 versus signal-to-noise ratio (the threshold for our variable model scheme is τJ=1/30).Figure 6
The corresponding computational complexities to Figure4 versus signal-to-noise ratio (the threshold for our variable model scheme is τJ=1/200).Figure 7
The corresponding computational complexities to Figure5 versus signal-to-noise ratio (the threshold for our variable model scheme is τJ=1/30).


## 6.1. Comparison of Computational Complexities

Figures1 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 orderN=2p is fixed for one-time LSE or MLR) for Hermitian matrix inversion using conventional Cholesky factorization (ℂR~Chol-1(N)) and our proposed new dyadic recursion (ℂR~new-1(N)).Figure 2
Computational complexity comparison (accumulated complexity for the sequential rank-adaptive LSE or MLR with terminal model orderN=2m) for Hermitian matrix inversion using conventional Cholesky factorization (ℂR~Chol,seq-1(N)) and our proposed new dyadic recursion (ℂR~new,seq-1(N)).Figure 3
Computational complexity margin ratio versus model orderN (fixed model order: (ℂR~new-1(N)-ℂR~Chol-1(N))/(ℂR~Chol-1(N))%, variable model order: (ℂR~new,seq-1(N)-ℂR~Chol,seq-1(N))/(ℂR~Chol,seq-1(N))%, equality line is for the margin ratio = 0%; N 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 Section5. 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 H(z)=0.58+0.56z-1+0.55z-2-0.18z-3-0.12z-4, where the channel gain is normalized as unity such that ∑n‍|h(n)|2=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 M=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σopt2 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 N=12 (specified as “fixed N=12”), (ii) one-time Hermitian matrix inversion using the Cholesky factorization with the fixed model order N=25 (specified as “fixed N=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 τJ=1/200,1/30. It is noted that the allowed maximum model order for our dyadic recursion scheme is defaulted as ⌊(N-M+1)/10⌋, where ⌊⌋ is the flooring operator.Figures4–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σopt2 versus signal-to-noise ratio (the threshold for our variable model scheme is τJ=1/200).Figure 5
Various least mean-squared errorsσopt2 versus signal-to-noise ratio (the threshold for our variable model scheme is τJ=1/30).Figure 6
The corresponding computational complexities to Figure4 versus signal-to-noise ratio (the threshold for our variable model scheme is τJ=1/200).Figure 7
The corresponding computational complexities to Figure5 versus signal-to-noise ratio (the threshold for our variable model scheme is τJ=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.



---

*Source: 891932-2012-11-05.xml*