1 Linear Algebra 1 Linear Algebra 1.2 Determinant Form for KummerU Function of a Matrix Argument

1.1 Novel Numerical Linear Algebra Algorithms

This is a consolidated post containing nifty numerical algorithms.

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

1.1.2 Complex-valued Durbin Algorithm with MATLAB code

This post is a rehash of my earlier post in Section 1.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