1 Linear Algebra 1.1 Inverse Cholesky Factor of a positive definite Hermitian symmetric Toeplitz matrix using Durbin recursions 1.3 Determinant Form for KummerU Function of a Matrix Argument

1.2 Complex-valued Durbin Algorithm with MATLAB code

This post is a rehash of my earlier post in Section 1.1 where the complex-valued Durbin’s algorithm was used indirectly for matrix inversion. Here, the algorithm is written down separately as it may be useful for some to have it “out of the box.”

The complex-valued versions of Durbin’s and Levinson’s algorithms may be arrived upon by straightforward application of description in [7] sec. 4.7.2 “Solving the Yule-Walker Equations”, and 4.7.3 “The General Right Hand Side Problem.”

Let $r_{0}=1,r_{1},r_{2},\cdots,r_{n}\in\mathbb{C}$ and $\boldsymbol{T}\in\mathbb{C}^{\text{N}\times\text{N}}$ be a positive definite Hermitian symmetric Toeplitz matrix constructed with the scalars $r_{0},r_{1},\cdots,r_{n-1}$ (see example below), Durbin’s algorithm solves the system $\boldsymbol{T}\boldsymbol{y}=-\begin{bmatrix}r_{1}&r_{2}&\vdots&r_{n}\end{% bmatrix}^{\mathrm{T}}.$ The algorithm is as follows.

Algorithm 2 Algorithm for finding

\boldsymbol{y}.

z(1)=-r_{1}^{*};\beta=1;\alpha=-r_{1}^{*}

2:for

k=1:n-1

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

\alpha=\frac{-1}{\beta}(r^{*}(k+1)-r(k:-1:1)^{H}z(1:k))

z(1:k)=z(1:k)+\alpha z^{*}(k:-1:1)

z(k+1)=\alpha

7:end for

y=z^{*}

MATLAB Reference Code

⬇

function Y = durbin(R)

% Y = durbin(R)

% where R is a vector with r0, r1, r2, ..., rn.

% Solves RY = -K, where R is the NxN symmetric Toeplitz matrix

% with elements r0, ..., r(n-1) and K is a vector with elements

% r1, ..., rn.

R = R/R(1,1) ;

N = length(R)-1 ;

R = R(2:end) ;

Z = zeros(N,1) ;

Z(1) = -conj(R(1)) ;

A = -conj(R(1)) ;

B = 1 ;

for k=1:N-1

B = (1-abs(A)^2)*B ;

A = -(conj(R(k+1)) + R(k:-1:1)'*Z(1:k))/B ;

Z(1:k) = Z(1:k) + A*conj(Z(k:-1:1)) ;

Z(k+1) = A ;

end

Y = conj(Z) ;

MATLAB Test Code

⬇

>> r = complex(rand(10,1),rand(10,1)) ;

>> r(1) = 1 ;

>> T = toeplitz(r(1:9),conj(r(1:9))) ;

>> k = r(2:end) ;

>> y1 = inv(T)*-k ;

>> norm(y1-durbin(r))

ans = 1.0709e-10

Version History

1.

First published: 30th Jan. 2016 on aravindhk-math.blogspot.com
2.

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