Modified Recursive Cholesky (Rchol) Algorithm: An Explicit Estimation and Pseudo-inverse of Correlation Matrices

TL;DR

This paper introduces a modified recursive Cholesky (RChol) algorithm for explicit estimation and pseudo-inversion of correlation matrices, offering a potentially faster and stable alternative to traditional methods like SVD and LU.

Contribution

It presents a new recursive algorithm for Cholesky decomposition that improves estimation of correlation matrices compared to conventional methods.

Findings

RChol algorithm provides accurate correlation matrix estimation.

It demonstrates computational efficiency over traditional methods.

The method enhances numerical stability in matrix inversion.

Abstract

The Cholesky decomposition plays an important role in finding the inverse of the correlation matrices. As it is a fast and numerically stable for linear system solving, inversion, and factorization compared to singular valued decomposition (SVD), QR factorization and LU decomposition. As different methods exist to find the Cholesky decomposition of a given matrix. This paper presents the comparative study of a proposed RChol algorithm with the conventional methods. The RChol algorithm is an explicit way to estimate the modified Cholesky factors of a dynamic correlation matrix.

Equations

$$\mathbf{y}(n):=\sum _{k=1}^K\sum _{l=0}^{L-1}{\mathbf{h}}^k(n;l)s^k(n-l)+\mathbf{v}(n),n=0,1,..T-1$$

$${\tiny\left[\begin{array}[]{ccccc}{\bf r}^{n}_{00}&{\bf r}^{n}_{01}&\ldots\ldots&{\bf r}^{n}_{0(N-2)}&{\bf r}^{n}_{0(N-1)}\\ {\bf r}^{n}_{10}&{\color[rgb]{0,1,1}{\bf r}^{n}_{11}}&{\color[rgb]{0,1,1}\ldots\ldots}&{\color[rgb]{0,1,1}{\bf r}^{n}_{1(N-2)}}&{\color[rgb]{0,1,1}{\bf r}^{n}_{1(N-1)}}\\ \vdots&{\color[rgb]{0,1,1}\vdots}&{\color[rgb]{0,1,1}\ddots}&{\color[rgb]{0,1,1}\vdots}&{\color[rgb]{0,1,1}\vdots}\\ {\bf r}^{n}_{(N-1)0}&{\color[rgb]{0,1,1}{\bf r}^{n}_{(N-1)1}}&{\color[rgb]{0,1,1}\ldots\ldots}&{\color[rgb]{0,1,1}{\bf r}^{n}_{(N-1)(N-2)}}&{\color[rgb]{0,1,1}{\bf r}^{n}_{(N-1)(N-1)}}\\ \end{array}\right]}$$

$${\tiny\left[\begin{array}[]{ccccc}{\color[rgb]{0,1,1}{\bf r}^{n}_{11}}&{\color[rgb]{0,1,1}{\bf r}^{n}_{12}}&{\color[rgb]{0,1,1}\ldots\ldots}&{\color[rgb]{0,1,1}{\bf r}^{n}_{1(N-1)}}&{\bf r}^{n}_{1N}\\ {\color[rgb]{0,1,1}\vdots}&{\color[rgb]{0,1,1}\vdots}&{\color[rgb]{0,1,1}\ddots}&{\color[rgb]{0,1,1}\vdots}&\vdots\\ {\color[rgb]{0,1,1}{\bf r}^{n}_{(N-1)1}}&{\color[rgb]{0,1,1}{\bf r}^{n}_{(N-1)2}}&{\color[rgb]{0,1,1}\ldots\ldots}&{\color[rgb]{0,1,1}{\bf r}^{n}_{(N-1)(N-1)}}&{\bf r}^{n}_{(N-1)N}\\ {\bf r}^{n}_{N1}&{\bf r}^{n}_{N2}&\ldots\ldots&{\bf r}^{n}_{N(N-1)}&{\bf r}^{n}_{NN}\\ \end{array}\right]}$$

$$[{\bf R}]_{1:N,1}=\frac{[{\bf R}]_{1:N,1}}{\color[rgb]{0,1,1}\sqrt{[{\bf R}]_{1,1}}}$$

$$[\mathbf{R}{]}_{k:N,k}=[\mathbf{R}{]}_{k:N,k}-[\mathbf{R}{]}_{k:N,1:k-1}[\mathbf{R}{]}_{k,1:k-1}^H$$

$$[{\bf R}]_{k:N,k}=\frac{[{\bf R}]_{k:N,k}}{\color[rgb]{0,1,1}\sqrt{[{\bf R}]_{k,k}}}$$

$$\mathbf{L}=\text{tril}(\mathbf{R})$$

$$[\mathbf{R}{]}_{2:N,1}=\frac{[\mathbf{R}{]}_{2:N,1}}{[\mathbf{R}{]}_{1,1}}$$

$$[\text{v}{]}_k=[\mathbf{R}{]}_{k,k}-[\mathbf{R}{]}_{k,1:k-1}[\text{v}{]}_{1:k-1}$$

$$[\mathbf{R}{]}_{k,k}=[\text{v}{]}_k$$

$$[\mathbf{R}{]}_{k+1:N,k}=\frac{[\mathbf{R}{]}_{k+1:N,k}-[\mathbf{R}{]}_{k+1:N,1:k-1}[\text{v}{]}_{1:k-1}}{[\text{v}{]}_k}$$

$$\mathbf{D}=diag(daig(\mathbf{R}))$$

$$\mathbf{L}=\text{tril}(\mathbf{R})$$

$$H_1(n)=[r_{00}^n,r_{10}^n,.....r_{(N-1)0}^n{]}^T$$

$${\tilde{H}}_1(n)=[0,r_{10}^n,.....r_{(N-1)0}^n{]}^T$$

$${\sigma }_k{H}_k={\tilde{k}}_{re{f}_k}[k_{re{f}_k}({\sigma }_{k-1}{\tilde{H}}_{k-1})+Z_M({\sigma }_{k-1}{H}_{k-1})]$$

$${\sigma }_k{\tilde{H}}_k={\tilde{k}}_{re{f}_k}[({\sigma }_n{\tilde{H}}_k)+k_{re{f}_k}{Z}_M({\sigma }_k{H}_k)]$$

$$k_{re{f}_k}=-\frac{({\sigma }_k{\tilde{H}}_k{)}_k}{({\sigma }_k{H}_k{)}_{k-1}}$$

$${\tilde{k}}_{re{f}_k}=\frac{({\sigma }_k{H}_k{)}_{k-1}}{({\sigma }_{k+1}{H}_{k+1}{)}_k}$$

$${\mathbf{A}}_1(n)=[{\mathbf{r}}_{00}^n,{\mathbf{r}}_{10}^n,\dots {\mathbf{r}}_{(N-1)0}^n{]}^T$$

$${\tilde{\mathbf{A}}}_1(n)=[0,{\mathbf{r}}_{10}^n,\dots {\mathbf{r}}_{(N-1)0}^n{]}^T$$

$${\mathbf{A}}_k(n):={\mathbf{Z}}_M{\mathbf{A}}_{k-1}(n-1)-{\tilde{\mathbf{A}}}_{k-1}(n){\tilde{\mathbf{k}}}_{\mathbf{ref}}(n)$$

$${\mathbf{D}}_k(n)={\mathbf{D}}_k(n-1)[{\mathbf{I}}_M-{\mathbf{k}}_{ref}(n){\tilde{\mathbf{k}}}_{ref}(n)]$$

$${\mathbf{A}}_k(n)={{\mathbf{Z}}_{\mathbf{M}}\mathbf{A}}_k(n-1)$$

$${\mathbf{D}}_k(n)={{\mathbf{Z}}_{\mathbf{M}}\mathbf{D}}_k(n-1)$$

Full text

Cholesky decomposition is a fast and numerically stable for linear system solving, inversion, and factorization compared to singular valued decomposition (SVD), QR factorization and LU decomposition $[1]$. The wireless communication system is highly dependent on matrix inversion of the correlation matrix. Such system consists of a huge matrix inversion. An outdoor wireless communication has a time-varying channel which changes dynamically for mobile user. In case of narrowband channel, the channel is considered constant for a symbol duration, whereas for broadband, it is changing within a symbol period. Such time-varying channel forms the special structure of channel matrix and correlation matrix. To exploit such special structure, a novel modified recessive Cholesky (RChol) algorithm is introduced in $[2]$. Our proposed (RChol) algorithm is a computational efficient algorithm to compute the modified Cholesky factors of known as well as an unknown covariance matrix.

In this paper, we present the comparative study of conventional Cholesky algorithm and the RChol algorithm to manifest the importance of the proposed algorithm in a highly dynamic wireless communication.

I System Model

In wireless communication system, number of transmit and or received antennas are used to improve the diversity of the system. The channel $‘h’$ between transmitter and receiver has the different form and depends on the number of antennas used at the transmitter and the receiver side. The channel for Single-input-single-output (SISO) as $h=\{h^{0}_{n},h^{1}_{n},\ldots h^{L-1}_{n}\}$, for Single-input-multiple-output (SIMO) as ${\bf h}=\{{\bf h}^{0}_{n},{\bf h}^{1}_{n},\ldots{\bf h}^{L-1}_{n}\}$ and for Multiple-input-multiple-output (MIMO): ${\bf H}=\{{\bf H}^{0}_{n},{\bf H}^{1}_{n},\ldots{\bf H}^{L-1}_{n}\}$.

Let y(n) be received signal with the number of transmit antennas ${}^{\prime}K=1^{\prime}$, multipath $‘L-1’$ and channel noise $‘v’$, represented as

[TABLE]

Let ${\bf y}_{N}(n)$ be the received vector by stacking ${}^{\prime}N^{\prime}$ successive received vectors. Where ${\bf y}_{N}(n)=[{\bf y}(n),{\bf y}(n-1),\ldots{\bf y}(n-N+1)]^{T}$ and the transmitted symbol vector is ${\bf s}_{N}=[s(n),s(n-1),….s(n-N+1)]^{T}$. Then ${\bf y}_{N}(n)$ can be represented in matrix form as ${\bf y}_{N}(n)={\bf H}_{N}{\bf s_{N}(n)+{\bf v}_{N}(n)}$ and the correlation matrix for ${\bf y}_{N}$ can be written as ${\bf R}_{N}(n)=E[{\bf y}_{N}(n){\bf y}_{N}^{H}(n)]$. Let ${\bf r}^{n}_{00}=E[{\bf y}(n){\bf y}^{H}(n)]$ and ${\bf r}^{n}_{ij}=E[{\bf y}(n-i){\bf y}^{H}(n-j)]$ then the correlation matrix ${\bf R}_{N}(n)$ and ${\bf R}_{N}(n-1)$ at time instant ${}^{\prime}n^{\prime}$ and ${}^{\prime}n-1^{\prime}$ can be represented as equation (2) and equation (3) respectively.

[TABLE]

[TABLE]

II Cholesky Decomposition

The correlation matrix is complex matrix and the pseudo-inverse of R can be computed from Cholesky factors, such that if lower triangular matrix $‘{\bf L}’$ is Cholesky factors of the correlation matrix $‘{\bf R}’$ and can be represented as ${\bf R}={\bf LL}^{H}$ then pseudo-inverse of R can be computed as ${\hat{\bf R}}^{\dagger}:={\bf L}^{-H}{\bf L}^{-1}$. The section below details the conventional Cholesky algorithms and the RChol algorithm.

II-A Cholesky Decomposition (Gaxpy version)

The Cholesky Decomposition $[3]$ factorizes a complex (or real-valued) positive-definite Hermitian symmetric matrix into a product of a lower triangular matrix and its Hermitian transpose. ${\bf R}={\bf LL}^{H}$ where, L is a lower triangular matrix and ${\bf L}^{H}$ is Hermitian of L. The matrix R must be a positive definite and this method needs square root operation.

II-A1 Algorithm steps

Compute R at each time instant n 2. 2.

Find the square root of diagonal element of R 3. 3.

Modify each column of R 4. 4.

Equate lower triangular part of R to L 5. 5.

Repeat steps $(1)$ to $(4)$ for each time instant

II-B Modified Cholesky Algorithm ${\bf R}={\bf LDL}^{H}$

To avoid square root operation, a modified Cholesky algorithm $[3]$ is used, which avoids square root operation by introducing a diagonal matrix D in between Cholesky factors. The modified Cholesky algorithm does not require R to be a positive definite matrix but it’s determinant must be nonzero. R may be rank deficient to a certain degree i.e. D may contain negative main diagonal entries if R is not positive semi-definite.

II-B1 Algorithm steps

Compute R at each time instant n 2. 2.

Modify each column of R 3. 3.

Equate the strictly lower part of matrix R to ${\bf L}_{1}$ with ones on the main diagonal 4. 4.

Equate main diagonal of R with the main diagonal of D 5. 5.

Repeat step $(1)$ to $(4)$ for each time instant.

II-C Recursive Cholesky Algorithm (The Shcur Algorithm) ${\bf R}_{Schur}:={\bf HH}^{H}$

The Schur algorithm recursively compute the columns of the lower triangular matrix H form matrix R. It is shown in $[4]$ that Levinson recursion may be used to derive the Lattice recursion for computing QR factors of data matrices and Lattice recursion can be used to derive the Schur recursion for computing Cholesky factors of a Toeplitz correlation matrix. The detail algorithm is given in algorithm $3$. The Schur algorithm like previously mentioned algorithm computes all $‘N’$ inner product to compute matrix R for initialization.

II-C1 Algorithm steps

Compute R at each time instant n 2. 2.

Initialize first column of R to the first column of Cholesky factor H 3. 3.

Compute rest column recursively from columns of R 4. 4.

repeat step $(1)$ to $(3)$ for each time instant

II-D * The RChol Algorithm ${\hat{\bf L}}{\hat{\bf D}}{\hat{\bf L}}^{H}:={\bf R}$*

It is clear from above equation (2) and equation (3) that ${\bf R}_{N}(n)$ can be represented from submatrix of ${\bf R}_{N}(n-1)$. To utilize such special structure of correlation matrices, we propose a modified recursive Cholesky algorithm to compute the Cholesky factors recursively. This algorithm is modification of Schur algorithm mentioned above. The more general approach consists of using the Schur algorithm to induce recursion for columns of dynamic L. This algorithm does not need N inner products to compute the correlation matrix R. The Cholesky factors are computed explicitly such that Let ${\bf L}_{1}={\bf LD}^{1/2}$then pseudo-inverse can be computed as ${\hat{\bf R}}^{\dagger}={\bf L}_{1}^{-H}{\bf L}_{1}^{-1}$

II-D1 Algorithm steps

Initialize first the first column of Cholesky factor A as ${\bf A}_{1}$ 2. 2.

Compute second column recursively from ${\bf A}_{1}(n)$ and ${\bf A}_{1}(n-1)$ 3. 3.

Substitute sub-matrix ${\bf A}_{2:N-1,2:N-1}(n-1)$ to ${\bf A}_{3:N,3:N}(n)$ 4. 4.

Repeat step $(1)$ to $(3)$ for each time instant

In the Schur algorithm, columns of Cholesky factors at time instant $‘n’$ are computed recursively from the correlation matrix at that instant. Whereas in the RChol algorithm first two columns of Cholesky factors at time instant $‘n’$ is computed recursively from previous Cholesky factor and submatrix of that Cholesky factors are updated recursively from previous Cholesky factor i.e. at time instant $‘n-1’$. Conventional Cholesky algorithm mentioned here are introduced for normal matrices whereas proposed matrix is well suited for block matrices and simulations are shown for that only.

III Simulation results

To compare proposed the RChol algorithm with Schur algorithm, we compared the result of both the algorithm with theoretical results. Fig. $1$. Show the ratio and difference of matrices ${\hat{\bf R}}_{N}$, ${\hat{\bf R}}_{RChol}$ and ${\hat{\bf R}}_{Schur}$ , when the correlation matrix is unknown. That has the application in blind channel and or data estimation. Fig. 1 (a) and (b) shows the maximum error for the RChol algorithm, $[{\hat{\bf R}}_{N}-{\hat{\bf R}}_{RChol}]$ is $‘0.6’$ while for the Schur algorithm, $[{\hat{\bf R}}_{N}-{\hat{\bf R}}_{Schur}]$ is $‘4’$ i.e. nearly 6 times the RChol algorithm. In case of ratio Fig. 1 (a) and (b) shows the maximum ratio for the RChol algorithm, $[{\hat{\bf R}}_{N}./{\hat{\bf R}}_{RChol}]$ is $‘45’$ while for the Schur algorithm, $[{\hat{\bf R}}_{N}./{\hat{\bf R}}_{Schur}]$ is $‘1500’$.

Fig. 1 Show the ratio and difference of matrices ${\hat{\bf R}}_{N}$, ${\hat{\bf R}}_{RChol}$ and ${\hat{\bf R}}_{Schur}$ , when the correlation matrix is known. Fig. 1 (a) and (b) shows that the maximum error for the RChol algorithm, $[{\hat{\bf R}}_{N}-{\hat{\bf R}}_{RChol}]$ is $‘2.5’$ while for the Schur algorithm, $[{\hat{\bf R}}_{N}-{\hat{\bf R}}_{Schur}]$ is $‘0.03’$ i.e. nearly 6 times the RChol algorithm. In case of ratio Fig. 1 (e) and (f) shows that the maximum ratio for the RChol algorithm, $[{\hat{\bf R}}_{N}./{\hat{\bf R}}_{RChol}]$ is $‘1.15’$ while for the Schur algorithm, $[{\hat{\bf R}}_{N}./{\hat{\bf R}}_{Schur}]$ is $‘1’$.

From Fig. 1 it can be concluded that the Schur algorithm is best suited when the correlation matrix is known, but leads to huge error propagation through the column when R is unknown and cannot be applied for blind channel estimation. In converse, the RChol algorithm is best suited for blind channel estimation and reduces error propagation through the column.

IV Conclusion

Convention methods of Cholesky factorization requires the correlation matrix which needs inner product. While the recursive modified Cholesky algorithm (RChol) algorithm is an explicit way to recursively calculating the pseudo-inverse of the matrices without estimating the correlation matrix. It requires less number of iteration which avoids error propagation through column updates. The RChol algorithm has most of the use in calculating the pseudo-inverse of the of a time-varying matrix which is applicable to SIMO/MIMO, CDMA, OFDM, etc. wireless communication systems.

V. Pawar and K. Naik (DIAT, Pune, India)

E-mail: vanietaapawar@gmail.com