Abstract-Square-root (in particular, Cholesky) factorization of Toeplitz matrices and of their inverses is a classical area of research. The Schur algorithm yields directly the Cholesky factorization of a symmetric Toeplitz matrix, whereas the Levinson algorithm does the same for the inverse matrix. The objective of this letter is to use results from the theory of rational orthonormal functions to derive square-root factorizations of the inverse of an n × n positive definite Toeplitz matrix. The main result is a new factorization based on the Takenaka-Malmquist functions, that is parameterized by the roots of the corresponding auto-regressive polynomial of order n. We will also discuss briefly the connection between our analysis and some classical results such as Schur polynomials and the Gohberg-Semencul inversion formula

Fellow, IEEE Bo Wahlberg

Fellow, IEEE Petre Stoica

CiteSeerX

1
New Square-Root Factorization of Inverse Toeplitz
Matrices
Bo Wahlberg, Fellow, IEEE and Petre Stoica, Fellow, IEEE
Abstract—Square-root (in particular, Cholesky) factorization
of Toeplitz matrices and of their inverses is a classical area of
research. The Schur algorithm yields directly the Cholesky fac-
torization of a symmetric Toeplitz matrix, whereas the Levinson
algorithm does the same for the inverse matrix. The objective of
this letter is to use results from the theory of rational orthonormal
functions to derive square-root factorizations of the inverse of an
n×n positive definite Toeplitz matrix. The main result is a new
factorization based on the Takenaka-Malmquist functions, that is
parameterized by the roots of the corresponding auto-regressive
polynomial of order n. We will also discuss briefly the connection
between our analysis and some classical results such as Schur
polynomials and the Gohberg-Semencul inversion formula.
Index Terms—Toeplitz matrix, Square-root and Cholesky fac-
torization, AR processes, Rational orthonormal functions.
I. INTRODUCTION
Consider the symmetric positive-definite Toeplitz matrices
(for given covariances {rk ∈ R} and for j = 1, . . . , (n+ 1))
Rj =

r0 r1 . . . rj−1
r1 r0 . . . rj−2
...
. . . . . .
...
rj−1 . . . r1 r0
 > 0. (1)
Let [aj,1, . . . , aj,j ] and σ2j be the solution to
Rj+1

1
aj,1
...
aj,j
 =

σ2j
0
...
0
 , j = 0 . . . n. (2)
This implies that {aj,k, k = 1 . . . j} are the coefficients of
the auto-regressive (AR) process of order j, with innovation
variance σ2j , whose covariances at lags 0 to j exactly match
r0 to rj . Define
Lj+1 =

1 0 . . . 0
aj,1 1 . . . 0
aj,2 aj−1,1 . . . 0
...
...
. . . . . .
...
aj,j aj−1,j−1 . . . a1,1 1
 . (3)
Copyright (c) 2008 IEEE. Personal use of this material is permitted.
However, permission to use this material for any other purposes must be
obtained from the IEEE by sending a request to pubs-permissions@ieee.org
B. Wahlberg is with the Automatic Control Lab and ACCESS, School
of Electrical Engineering, KTH, SE-100 44 Stockholm, Sweden. (e-mail:
bo.wahlberg@ee.kth.se.)
P. Stoica is with the Department of Information Technology, Uppsala
University, SE-75105 Uppsala, Sweden (e-mail: ps@it.uu.se)
This work was partially supported by the Swedish Research Council and
the Linnaeus Center ACCESS at KTH
Then
R−1j+1 = Lj+1D
−1
j+1L
T
j+1, (4)
where
Dj+1 = diag(σ2j , σ
2
j−1, . . . , σ
2
0)
is a diagonal matrix, see e.g. [1]-[2], and the many references
therein. This well-known result, which forms the foundation
of the Levinson algorithm, shows that LnD
−1/2
n is a lower
triangular Cholesky factor of the inverse Toeplitz matrix R−1n .
The factorization in (4) is a function of the coefficients of all
AR models of orders 0 to n−1. In contrast to (4), the Gohberg-
Semencul (G-S) formula for R−1n+1 , see, e.g., [3],[2], and [4],
(which is not of the Cholesky or, more generally, square-root
factorization type) is a function only of the coefficients of the
AR model of order n, i.e. {an,1, . . . , an,n, σ2n}. There is also
a similar formula for R−1n , see e.g. [2], which again is an
explicit function of the coefficients of the AR model of order
n.
The objective of this letter is to show how the theory of
rational orthonormal functions can be used for finding a class
of square-root factorizations of R−1n or, if desired, of R
−1
n+1. In
particular, using the Takenaka-Malmquist functions, [5, page
18], we obtain an apparently new square-root factorization of
R−1n that is an explicit function of the roots of the nth order
AR polynomial. We should stress that, similarly to the G-S
theory, our main contribution consists of a theoretical inversion
formula, rather than a computational inversion algorithm.
II. PRELIMINARIES ON RATIONAL ORTHONORMAL
FUNCTIONS
Define the normalized nth order AR polynomial associated
with the given covariances {rk}nk=0 :
An(z) =
1
σn
(zn + an,1zn−1 + . . .+ an,n), (5)
and let
v(z) = [zn−1 zn−2 . . . 1]T , V (z) =
v(z)
An(z)
. (6)
Since the covariances are given by
rk =
1
2π
∫ π
−π
eiωk
|An(eiω)|2
dω, k = 0, 1, 2 . . . ,
we have
Rn =
1
2π
∫ π
−π
V (eiω)V ∗(eiω)dω,
2
where the superscript ∗ denotes the conjugate transpose. The
following function, where z1, z2 ∈ C, will play an important
role in the analysis of this letter:
∆(z1, z2) =
vT (1/z1)R−1n v(z2)
An(1/z1)An(z2)
= V T (1/z1)
[
1
2π
∫ π
−π
V (eiω)V ∗(eiω)dω
]−1
V (z2).
The function ∆(z1, z2) is called the reproducing kernel for the
space spanned by the elements of V (z). A key observation
is that ∆(z1, z2) is invariant under a change of coordinates
Ṽ (z) = TV (z), V (z) = T−1Ṽ (z), where T is an n × n
invertible matrix. Indeed
∆(z1, z2) = Ṽ T (1/z1)T−T
×
[
1
2π
∫ π
−π
T−1Ṽ (eiω)Ṽ ∗(eiω)T−T dω
]−1
T−1Ṽ (z2)
= Ṽ T (1/z1)
[
1
2π
∫ π
−π
Ṽ (eiω)Ṽ ∗(eiω)dω
]−1
Ṽ (z2).
A natural idea is then to use an ortho-normal basis for the
space spanned by {zk/An(z), k = 0, . . . , (n−1)}, for which
1
2π
∫ π
−π
Ṽ (eiω)Ṽ ∗(eiω)dω = I,
(I denotes the identity matrix) and, therefore,
∆(z1, z2) = Ṽ T (1/z1)Ṽ (z2).
The theory of rational ortho-normal (ON) functions is well
developed. Two books on this subject are [6] and [5], which
both contain extensive reference lists. Let {ξj}nj=1 be the roots
of An(z) = 0, i.e.
An(z) =
1
σn
n∏
j=1
(z − ξj).
Furthermore, assume that the roots satisfy |ξj | < 1, which is
always the case under (1), see e.g. [1]. Applying the Gram-
Schmidt procedure to the functions {1/(z−ξj)k}, k = 1 . . . nj
(where nj is the multiplicity of root j) over all roots (these
functions span the same space as the elements of V (z), (6))
leads to the Takenaka-Malmquist (T-M) ON functions, [5],
Bk(z) =
√
1− |ξk|2
z − ξk
k−1∏
i=1
[
1− ξ∗i z
z − ξi
]
, k = 1, 2, . . . , n. (7)
The T-M construction is valid for arbitrary root configurations
with |ξj | < 1. The special case of A(z) = (z − a)n,
i.e. repeated roots in ξj = a, -1 < a < 1, leads to the
well-known Laguerre functions, which further simplify to the
delay functions Bk(z) = z−k for ξj = 0. In the case
of complex conjugate roots, say ξj and ξj+1 = ξ∗j , define
(z− ξj)(z− ξj+1) = z2 + bj(cj − 1)z− cj with −1 < bj < 1
and −1 < cj < 0. Then take
Bj(z) =
√
1− c2j (z − bj)
z2 + bj(cj − 1)z − cj
j−1∏
i=1
[
1− ξ∗i z
z − ξi
]
Bj+1(z) =
√
1− c2j
√
1− b2j
z2 + bj(cj − 1)z − cj
j−1∏
i=1
[
1− ξ∗i z
z − ξi
]
to obtain ON functions with real coefficients, [5]. For repeated
complex roots, the functions obtained by this construction are
often refereed to as Kautz functions. The vector function made
from {Bk(z)} is Ṽ (z) = [B1(z), . . . , Bn(z)]T and thus
∆(z1, z2) =
n∑
k=1
Bk(1/z1)Bk(z2) (8)
Another set of ON functions can be constructed from the so-
called Schur polynomials Ak(z),
Ak(z) =
1
σk
(zk + ak,1zk−1 + . . .+ ak,k), k = 0, . . . n, (9)
as follows, [5, Section 2.6.2],
Bk(z) =
Ak−1(z)
An(z)
, k = 1, . . . , n,
where {aj,1, . . . , aj,j , σj} are the solutions to (2) for j = 0
to n. There are simple algorithms to find Ak(z) directly from
An(z) using Schur iterations.
In a somewhat different context, we remark on the fact that
an important result in the theory of rational ON functions is
the so-called Christoffel-Darboux formula, see [5]
Ṽ T (1/z1)Ṽ (z2) =
H(1/z1)H(z2)− 1
1− z2/z1
,
H(z) =
znAn(1/z)
An(z)
. (10)
We will make use of (10) to briefly show how to connect
the square-root factorization approach of this letter to the G-S
formula for R−1n .
III. NEW SQUARE-ROOT FACTORIZATION
Using expression (8), we have that
vT (1/z1)R−1n v(z2) = An(1/z1)An(z2)∆(z1, z2)
=
n∑
k=1
Ck(1/z1)Ck(z2), (11)
Ck(z) = An(z)Bk(z) (12)
It is possible to recover the matrix R−1n from (11) by noting
that the vectors v(eiωnk), ωn = 2π/n, for k = 0, . . . , (n− 1)
are orthogonal (a well known result from the Discrete Fourier
Transform (DFT) theory). Define the n× n DFT matrix
Wn =
1√
n
[
v(eiωn0) v(eiωn1) . . . v(eiωn(n−1))
]
,
for which W ∗nWn = I . We also define the row vectors
γk =
1√
n
[
Ck(eiωn0) Ck(eiωn1) . . . Ck(eiωn(n−1))
]
.
3
Using this notation we can write from (12):
W ∗nR
−1
n Wn =
n∑
k=1
γ∗kγk ⇒
R−1n =
n∑
k=1
[γkW ∗n ]
∗γkW
∗
n = [ΓnW
∗
n ]
∗ΓnW ∗,
Γn =

γn
γn−1
...
γ1
 . (13)
We have thus derived a square-root factorization1 of R−1n .Note
that the functions Ck(z) = An(z)Bk(z) are polynomials of
degree n− 1,
Ck(z) = ck,1zn−1 + ck,2zn−2 + . . .+ ck,n,
for all choices of minimal (of dimension n) ON bases that
span the same space as {zk/An(z)}(n−1)k=1 . This means that
we can write Γn = UTnWn with
Un =
 cn,1 . . . c1,1... . . . ...
cn,n . . . c1,n
 . (14)
and therefore the factorization (13) can be re-written
R−1n = UnU
T
n . (15)
The simplest case is that of Schur ON polynomials, see (9),
for which Ck(z) = Ak−1(z). For the Takenaka-Malmquist
functions we have
Ck(z) =
√
1− |ξk|2
σn
k−1∏
i=1
[1− ξ∗i z]
n∏
i=k+1
[z − ξi] ,
k = 1, . . . , n, (16)
with the modification described in the previous section to
handle complex-valued roots.
To summarize, let {Bk(z)} be a rational ON basis for
the space spanned by {zk/An(z)}(n−1)k=0 , and define Ck(z) =
An(z)Bk(z), k = 1 . . . n and Un as in (14). The key result
of this letter is the general square-root factorization formula
(15).
In the Schur polynomials case (15) reduces to the well-
known formula (3) with LnD
−1/2
n = Un. On the other
hand for the Takenaka-Malmquist functions we obtain a novel
factorization of R−1n that depends only on the zeros of An(z)
and on σn, c.f. (16).
Note that it is easy to modify the result (15) to factorize
R−1n+1, instead of R
−1
n , still using only An(z). Let
vn+1(z) = [zn zn−1 zn−2 . . . 1]T , S(z) =
vn+1(z)
An(z)
,
1By a matrix square-root factorization of M we mean M = NN∗. This
includes the Cholesky factorization of M , for which N has a lower or upper
triangular structure. It also includes the symmetric (Hermitian) matrix square
root R of M , for which M = RR.
for which
Rn+1 =
1
2π
∫ π
−π
S(eiω)S∗(eiω)dω.
By extending {B1(z) . . . Bn(z)} with Bn+1(z) = 1, we
obtain an ON basis for the n+ 1 dimensional space spanned
by {zk/An(z)}nk=0. This is so because Bn+1(z) = 1 is
orthogonal to all Bk(z) (which all have relative degree one).
Define as before Ck(z) = An(z)Bk(z), k = 1 . . . n, and in
addition Cn+1(z) = An(z), and let:
Cn+1(z) = cn+1,1zn + cn+1,2zn−2 + . . . . . . cn+1,n+1.
Define
Ūn =

cn+1,1 0 . . . 0
cn+1,2 cn,1 . . . c1,1
...
...
. . .
...
cn+1,n+1 cn,n . . . c1,n
 . (17)
Using the same arguments as for the case of Rn now gives
R−1n+1 = ŪnŪ
T
n .
This result can be viewed as a consequence of applying the
Levinson algorithm to go from dimension n to n+ 1.
Finally, we note that the use of the Christoffel-Darboux
formula leads to:
vT (1/z1)R−1n v(z2) = An(1/z1)An(z2)∆(z1, z2)
=
(z2/z1)nAn(z1)An(1/z2)−An(1/z1)An(z2)
1− z1z2
. (18)
Using the DFT technique of this letter to recover R−1n from
(18) will result in the Gohberg-Semencul formula. The connec-
tion between the Christoffel-Darboux and Gohberg-Semencul
formulas is discussed in, e.g. [4].
IV. APPLICATION OF THE NEW FACTORIZATION FORMULA
Let {y(t)}nt=1 be a sequence of a stochastic signal whose
power spectral density is Φ(eiω), ω ∈ [−π, π). Also let Rn
denote a generic Toeplitz (covariance) matrix associated with
an arbitrary mth order AR process with poles strictly inside
the unit circle and with m < n (m being a fixed integer that
does not depend on n). Under the Gaussian assumption, the
dominant term of the negative log-likelihood function that can
be used to fit an mth order AR model to {y(t)}nt=1 is given
by, see e.g., [1]:
Y TR−1n Y, Y = [y(1) y(2) . . . y(n)]
T . (19)
Because (19) is a complicated function of the AR parameters,
the minimization of (19) is often replaced by that of
n∑
t=m+1
[Ām(q)y(t−m)]2, (20)
Ām(q) = qm + a1qm−1 + . . .+ am,
where q is the forward shift operator. Notice that Ā(z) =
σmAm(z), due to the normalization of the coefficient of
zm. Minimizing (20) is a quadratic least-squares optimiza-
tion problem in the coefficients of Ām with a closed form
4
solution. In what follows, we make use of the new square-
root factorization of R−1n introduced in this letter, see (15), to
prove that2:
lim
n→∞
1
n
E{Y TR−1n Y } =
1
σ2m
E[Ām(q)y(t)]2, (21)
which provides support to the fact that (20) is a good approxi-
mation of (19) asymptotically (as n increases). Using (15) and
(16) yields:
E{Y TR−1n Y } = E{
n∑
k=1
[Ck(q)y(t− n)]2}
=
n∑
k=1
1
2π
∫ π
−π
|Ck(eiω)|2Φ(eiω)dω
=
1
2π
∫ π
−π
|An(eiω)|2Φ(eiω)
n∑
k=1
|Bk(eiω)|2dω, (22)
where, because we assume the process to be of order m,
An(z) = Am(z)zn−m, and thus An(z) has n − m roots at
z = 0. This implies that Bk = z−k, k = m + 1, . . . n both
for the Takenaka-Malmquist functions and for the Schur basis.
Hence,
n∑
k=1
|Bk(eiω)|2 =
m∑
k=1
|Bk(eiω)|2 + n−m, (23)
where an expression for the first term on the right hand side is
most easily obtained using the Takenaka-Malmquist functions,
|Bk(eiω)|2 =
1− |ξk|2
|eiω − ξk|2
Hence, we have
1
n
E{Y TR−1n Y } =
1
2π
∫ π
−π
|Am(eiω)|2Φ(eiω)
×
(
1
n
m∑
k=1
|Bk(eiω)|2 +
n−m
n
)
dω
→ 1
σ2m
E[Ām(q)y(t)]2, n→∞, (24)
which proves (21). Note that the function (23) is actually
a fundamental scale factor for the asymptotic variance of
the corresponding autoregressive spectral density estimate as
shown in [7]. The papers [8] and [9] contain further results in
this direction. In fact, these results inspired our derivation of
the square-root factorization (15). The Schur version of this
formula was used in [10] to establish convergence rates for
quadratic forms associated with inverse Toeplitz matrices as n
tends to infinity.
Another possible application of the new square root factor-
ization formula for R−1n , which is left for future research, is
to establishing a novel relationship between the Capon and the
AR method in the spectral analysis of stationary signals. See
e.g. [11] for an introduction to this subject and some known
relationships.
2E denotes the expectation operator
V. CONCLUSIONS
Toeplitz matrices and their inverses play a fundamental role
both in the implementation and in the analysis of parametric
estimation methods, and elsewhere. The Levinson algorithm
and the Schur algorithm (as well as the Gohberg- Semencul
formula) are classical results on how to find factorizations
of such matrices. In this letter we have derived a class of
square-root factorizations of inverse Toeplitz matrices, using
rational ON functions. In particular, we have shown that
the use of Takenaka-Malmquist functions in our framework
leads to an apparently new square-root factorization formula.
Interestingly, the magnitude of the roots of the autoregressive
polynomial appear to play a role for the existence of this
factorization, similar to that played by the magnitude of the
reflection coefficients in the Levinson case.
REFERENCES
[1] T. Söderström and P. Stoica, System Identification. Prentice Hall
International, Hemel Hempstead, 1994.
[2] M. H. Gutknecht and M. Hochbruck, “The stability of inversion formulas
for Toeplitz matrices,” Linear Algebra and its Applications, vol. 223-
224, pp. 307 – 324, 1995, honoring Miroslav Fiedler and Vlastimil Ptak.
[3] G. Heinig and K. Rost, Algebraic Methods for Toeplitz-like Matrices and
Operators. Akademie-Verlag, Berlin and Birkhuser, Basel/Stuttgart,
1984.
[4] T. Kailath and J. Chun, “Generalized Gohberg-Semencul formulas for
matrix inversion,” Oper. Theory Adv. Appl, vol. 40, p. 231246, 1989.
[5] P. Heuberger, B. Wahlberg, and P. Van den Hof, Eds., Modelling and
Identification with Rational Orthogonal Basis Functions. Springer
Verlag. London, 2005.
[6] A. Bultheel, P. González-Vera, E. Hendriksen, and O. Njåstad, Orthog-
onal Rational Functions. Cambridge Monographs on Applied and
Computational Mathematics, Cambridge University Press, UK, 1999.
[7] L.-L. Xie and L. Ljung, “Variance expressions for spectra estimated
using auto-regressions,” Journal of Econometrics, vol. 118, pp. 247 –
256, 2004.
[8] B. Ninness and H. Hjalmarsson, “Variance error quantifications that are
exact for finite model order,” IEEE Transactions on Automatic Control,
vol. 49, no. 8, pp. 1275–1291, Aug. 2004.
[9] H. Hjalmarsson and B. Ninness, “Least-squares estimation of a class
of frequency functions,” Automatica, vol. 42, no. 4, pp. 589–600, May
2006.
[10] E. Hannan and B. Wahlberg, “Convergence rates for inverse Toeplitz
matrix forms,” Journal of Multivariate Analysis, vol. 31, pp. 127–135,
1989.
[11] P. Stoica and R. L. Moses, Spectral Analysis of Signals. Prentice Hall,
New Jersey, 2005.


New Square-Root Factorization of Inverse Toeplitz Matrices

http://citeseerx.ist.psu.edu/viewdoc/download;jsessionid=370D0F7B39DBC7C41EAEC951BDCA7437?doi=10.1.1.1046.5242&rep=rep1&type=pdf

New Square-Root Factorization of Inverse Toeplitz Matrices

Abstract

Similar works

Full text

Available Versions

CiteSeerX