1 Linear Algebra 1 Linear Algebra 1.2 Complex-valued Durbin Algorithm with MATLAB code

1.1 Inverse Cholesky Factor of a positive definite Hermitian symmetric Toeplitz matrix using Durbin recursions

Let $\boldsymbol{T}\in\mathbb{C}^{n\times n}$ be a positive-definite Hermitian symmetric Toeplitz matrix defined with elements $t_{ij}=r_{i-j}$ for a set of scalars $r_{-n+1},\cdots,r_{0},\cdots,r_{n-1}$ , where $r_{-t}=r_{t}^{*}$ . Without loss of generality, assume $r_{0}=1$ .

We have the Cholesky factor $\boldsymbol{R}$ of the matrix $\boldsymbol{T}$ such that $\boldsymbol{R}^{\mathrm{H}}\boldsymbol{R}=\boldsymbol{T}$ . The inverse Cholesky factor is then given by $\boldsymbol{R}^{-1}$ . A closely related decomposition known as the LDL decomposition is given as follows: $\boldsymbol{\Lambda}=\text{diag}(\boldsymbol{R})$ , $\boldsymbol{L}=\boldsymbol{R}^{\mathrm{H}}\boldsymbol{\Lambda}^{-1}$ , and therefore, $\boldsymbol{T}=\boldsymbol{L}\boldsymbol{\Lambda}^{2}\boldsymbol{L}^{\mathrm{H% }}=\boldsymbol{L}\boldsymbol{D}\boldsymbol{L}^{\mathrm{H}}$ . Matrix $\boldsymbol{L}$ has ones along the principal diagonal. Let $\boldsymbol{\Sigma}=\boldsymbol{L}\boldsymbol{D}$ .

We have $\boldsymbol{T}\boldsymbol{L}^{-\mathrm{H}}=\boldsymbol{\Sigma}$ , where $\boldsymbol{L}^{-\mathrm{H}}$ is the Hermitian transpose of the inverse with ones along the diagonal. As $\boldsymbol{\Sigma}$ is lower triangular, we can solve for the strictly upper triangular parts of $\boldsymbol{L}^{-\mathrm{H}}$ as the R.H.S. is known to be zero.

For the $k$ -th column of $\boldsymbol{L}^{-\mathrm{H}}$ we have the first $k-1$ elements $\boldsymbol{z}_{k-1}$ , upon simplification, as:

\boldsymbol{T}_{k-1}\boldsymbol{z}_{k-1}=-\begin{bmatrix}r_{k+1}\\ \vdots\\ r_{1}\end{bmatrix}

For Hermitian symmetric Toeplitz matrix $\boldsymbol{T}$ we have $\boldsymbol{T}=\boldsymbol{E}_{n}\boldsymbol{T}^{*}\boldsymbol{E}_{n}$ , where $\boldsymbol{T}^{*}$ is complex conjugate of the matrix and $\boldsymbol{E}_{n}$ is the $n\times n$ exchange matrix. Using this property, and setting $\boldsymbol{y}_{k-1}=\boldsymbol{E}_{k}\boldsymbol{z}^{*}_{k-1}$ we have:

\boldsymbol{T}_{k-1}\boldsymbol{y}_{k-1}=-\begin{bmatrix}r^{*}_{1}\\ \vdots\\ r^{*}_{k+1}\end{bmatrix}

which may be solved efficiently using the Durbin iterations [7]. The algorithm for finding $\boldsymbol{W}=\boldsymbol{R}^{-1}=[w_{i,j}]$ is given as follows:

Algorithm 1 Algorithm for finding

\boldsymbol{W}.

y_{1}=-r^{*}_{1};\enspace\alpha=-r^{*}_{1};\enspace\boldsymbol{W}=\boldsymbol{% I}_{n}

\beta=1-|\alpha|^{2}

w_{1,2}=y^{*}_{1}

\boldsymbol{w}_{:,2}=\frac{\boldsymbol{w}_{:,2}}{\sqrt{\beta}}

5:for

k=1:n-2

\alpha=-(r^{*}_{k+1}+\frac{1}{\beta}\boldsymbol{r}^{\mathrm{H}}_{k:-1:1}% \boldsymbol{y}_{1:k})

\boldsymbol{y}_{1:k}=\boldsymbol{y}_{1:k}+\alpha\boldsymbol{y}^{*}_{k:-1:1}

y_{k+1}=\alpha

\beta=(1-|\alpha|^{2})\beta

10:

\boldsymbol{w}_{1:k+1,k+2}=\boldsymbol{y}^{*}_{k+1:-1:1}

11:

\boldsymbol{w}_{:,k+2}=\frac{\boldsymbol{w}_{:,k+2}}{\sqrt{\beta}}

12:end for

The algorithm requires $2n^{2}$ flops [7, Algorithm 4.7.1] for the Durbin iterations, $n^{2}/2$ flops for scaling with $1/\sqrt{\beta}$ and $n$ inverse square root computations; and may be easily adapted to find $\boldsymbol{L}^{-1}$ and $\boldsymbol{D}$ of the LDL decomposition. A MATLAB reference code is available in [14] as well as reproduced below.

MATLAB Reference Code

⬇

function R = invchol_durbin(T)

% R = invchol_durbin(T)

% Finds the inverse of the Cholesky factor of the positive definite

% Hermitian symmetric Toeplitz matrix T (N>=2) using Durbin recursions [1].

% Aravindh Krishnamoorthy, aravindh.krishnamoorthy@fau.de, 04-Sep-2015.

% Released under the 2-clause BSD license.

% [1] Gene H. Golub, Charles F. Van Loan, Matrix Computations, Third

% Edition, Algorithm 4.7.1 (Durbin).

% Exchange matrix of order k

E = @(k) rot90(eye(k)) ;

% Normalise matrix

ts = T(1,1) ;

T = T/ts ;

t = transpose(T(1,:)) ;

%% Durbin recursions [1]

% Setup variables

N = size(T,1) ;

r = t(2:end) ;

y = zeros(N,1) ;

R = eye(N) ;

% Initial values

y(1) = -conj(r(1)) ;

b = 1 ;

a = -conj(r(1)) ;

% Update 'b'

b = (1-abs(a)^2)*b ;

% Column 2

R(1,2) = conj(y(1)) ;

R(:,2) = R(:,2)/sqrt(b) ;

% Rest of the columns

for k=1:N-2

% Update 'a' and the result

a = -(conj(r(k+1)) + ctranspose(r(k:-1:1))*y(1:k))/b ;

y(1:k) = y(1:k) + a*conj(y(k:-1:1)) ;

y(k+1) = a ;

% Update 'b'

b = (1-abs(a)^2)*b ;

% Write out the result

R(1:k+1,k+2) = E(k+1)*conj(y(1:k+1)) ;

R(:,k+2) = R(:,k+2)/sqrt(b) ;

end

% Scale back the result

R = R/sqrt(ts) ;

Version History

1.

First published: 4th Sep. 2015 on aravindhk-math.blogspot.com
2.

Modified: 16th Dec. 2023 – Style updates for LaTeX