On the inverse problem of fractional Brownian motion and the inverse of infinite Toeplitz matrices

The inverse problem of fractional Brownian motion and other Gaussian processes with stationary increments involves inverting an infinite hermitian positively definite Toeplitz matrix (a matrix that has equal elements along its diagonals). The problem of inverting Toeplitz matrices is interesting on its own and has various applications in physics, signal processing, statistics, etc. A large body of literature has emerged to study this question since the seminal work of Szeg\"o on Toeplitz forms in 1920's. In this paper we obtain, for the first time, an explicit general formula for the inverse of infinite hermitian positive definite Toeplitz matrices. Our formula is explicitly given in terms of the Szeg\"o function associated to the spectral density of the matrix. These results are applied to the fractional Brownian motion and to $m$-diagonal Toeplitz matrices and we provide explicit examples.


Introduction
The statistical inverse solution to the inverse problem involves quite often the inverse of a large or infinite matrix. Generally the measurement vector is modeled as a realisation of a stochastic process with unknown parameters and the problem is to recover these parameters from the measurements. In this paper we are mainly concerned with the fractional Brownian motion which is a Gaussian process { ( ) }  X t t : 0 defined by the covariance function where ( ) Î H 0, 1 is a parameter called the Hurst index and > a 0 is a scale parameter. The fractional Brownian motion has stationary increments but these increments are not independent. It is an important model that has been successfully used to model real data that appear in different contexts such as anomalous diffusion, singlefile dynamics, volatility in financial markets, hydrology, geophysics, queueing networks such as internet traffic, etc ( [1][2][3][4]). In what follows we shall assume that a=1 and in this case X(t) is called the standard fractional Brownian motion. The inverse problem of fractional Brownian motion deals with the recovery of the Hurst parameter H from a given vector of measurements. For this it is necessary to consider the increments The process { } D Î  n : n is called the (discrete) fraction Gaussian noise. The statistical inverse method for fractional Brownian motion requires the computation of the inverse of the matrix G. In fact assume that ( ) ( ) , , . Of course it is necessary to take n very large for a better approximation and ultimately  ¥ n for the exact value. To recover the value of the unknown parameter H, it is important to have an analytic expression of the elements of the inverse matrix ( ) -G n n 1 for large n and ultimately for -G 1 as functions of H. However such analytic forms are not known. Up to so far in the literature, because no closed form for -G 1 is known, the matrix G has to be approximated by some simpler matrices. One example is to consider the inverse of the matrix˜(˜) = G g n m , given bỹ (See for example [5].) A more widely used approximation in the literature for the matrix -G 1 is the asymptotic approximation introduced by Whittle [6] in 1953. It consists of the matrix Γ given by: is the spectral density function of the fractional Gaussian noise of index H (introduced in section 2). (We refer to Beran [7] for more details.) This matrix Γ corresponds to the inverse -G 1 only asymptotically in the sense that The objective of this paper is to complete these earlier efforts of Whittle and others by providing for the first time an analytic expression that gives the exact value of the inverse matrix of a general infinite, Hermitian and positive definite Toeplitz matrix.
Our main result can be summarized as follows. If G is an infinite, Hermitian and positive definite Toeplitz matrix, with inverse -G 1 , then the elements of -G 1 are given by: where ( ) a n are fixed complex numbers, then for all k j , with  j k, We have that the Taylor coefficients ( ) a n of ( ) y z are given by: a  n  n  k  u a  n  0  1   1  1  1  ,  1, 2 The upper left-hand n×n bloc of - These general results are then applied to the case of the covariance matrix of the fractional Gaussian noise and other matrices. Toeplitz matrices appear in various contexts in physics such as condensed matter physics, entanglement theory, electrical engineering and chemical physics and has many important applications. (See [8,9] and references therein for more details.) The problem of inversion of Toeplitz matrices has been studied for many years since the seminal algorithm of Trench [10]. It has important applications and it is interesting on its own (see for example [11][12][13][14]).
From the point of view of numerical applications, the problem can be considered as solved as there exist fast algorithms to compute the inverse of these matrices (see for example [15]). However as mentioned earlier, in terms of parameter estimation, it is desirable to obtain closed analytic expressions for inverses of Toeplitz matrices. An analytic expression can also be used to analyse the performance of numerical algorithms. Furthermore from the mathematical modelling point of view, it is important to see how the elements of the inverse matrix are directly related to elements of the initial matrix or other associated characteristics such as eigenvalues. Many Toeplitz matrices that arise in Physics such as the sine-kernel for and i j ii , are given in analytical forms and it is important to find the analytic form for the corresponding inverse (in the case where it exists).
The rest of the paper is organised as follows. Section 2 discusses the notion of spectral density function section 3, the notion of Szegö function of a Toeplitz matrix. The link between the inverse of a Toeplitz matrix and its Szegö function is given section 4 and the explicit formula for the elements of the inverse matrix is given in sections 5 and 6. The case of the inverse of the covariance matrix of the fractional Gaussian noise is given in section 7 and further illustrative examples are given in section 8. Concluding remarks are given in section 9.

Spectral density function
1 ,2 , be an infinite, Hermitian and positive definite Toeplitz matrix. We assume that G is such that where γ is a positive definite function defined on the set of integers  and taking values in the set of complex numbers  and satisfies This implies (by the classical Bochner theorem) that there exists a probability measure μ on the unit interval We shall assume that μ admits a probability density function j  0 that is integrable (more precisely j p is integrable for some number > p 1) so that its Fourier series The function j is the spectral density function of the matrix G. We shall assume that G is invertible in the sense that there exists another positive definite matrix - where I is the infinite identity matrix. Our goal is to compute the inverse matrix -G 1 . If G is the covariance matrix of fractional Gaussian noise given by relation (1), then clearly G is the Toeplitz matrix corresponding to the function It is an important result obtained by Sinai [16] that G admits a spectral density function j H given by Here C(H) is a normalising constant given by ,. is the classical Hurwitz zeta function. Moreover it is well-known that The Fourier series of j H (given by relation (7)) yields the following representation: where Li is the analytic continuation of the classical polylogarithm function in the complex plane (except at the single point z = 1). This representation is amenable to calculations (for example the function Li is implemented in Wolfram Mathematica).

Szegö function of a Toeplitz matrix
Assume that the matrix G admits a spectral density function j. Under the condition that the spectral density function j satisfies the Szegö condition The function S(z) is known as the Szegö function associated to the matrix G (or equivalently associated to j).
(We refer the reader to Grenander and Szegö [18] and Simon [19] for more details.) We shall mainly use the inverse ( ) given by (3). It is also known that if the Szegö function S(z) is written in the form, are complex numbers, then these coefficients satisfy the system: together with the condition that c 0 is a positive real number. Assume moreover that ( ( )) = -P z n degree 1 n for all = ¼ n 1, 2, and further that the leading coefficient ofz n 1 is a positive real number. These conditions uniquely determine the sequence ( ( )) P z n . Szegö proved that the function ψ is such that for all (See the book by Grenander and Szegö [18], pp. 37-51.) For the case of the fractional Gaussian noise, since as mentioned before, the spectral density function at the boundary of the interval, that is for t near 0 and t near 1, then it is clear that j H satisfies the Szegö condition (8) and hence the fractional Gaussian noise admits a Szegö function. We shall denote its inverse by y H .

Szegö function and the inverse of an infinite Toeplitz matrix
be arbitrary sequence of orthonormal polynomials on the unit circle with respect to the weight function j such that Here it is not necessary to impose any extra condition on the coefficient ofz n 1 in Q n (z) as it is the case for the polynomials P n (z). We have that any such general sequence of orthogonal polynomials satisfies: This is a consequence of an important result about the structure of the space H of holomorphic functions   f D : of the form: where ( ) a n is a sequence of complex numbers satisfying the following two conditions: In fact it is shown in [20] that the space H is a Hilbert space with respect to the inner product . .) Moreover H is a reproducing kernel Hilbert space associated to the kernel This means that for every function H Î f and for all | | < z 1, .
it it 0 1 2 2 Or equivalently, by the fact that G is hermitian. We shall write where A key property of a reproducing kernel Hilbert space is that the quantity is invariant for all possible choices of orthonormal basis f f f ¼ , , and in fact it is equal to the kernel of the space. Now it is again proven in [20] that the monomials ( ) =f z z n 1 for all = ¼ n 1, 2, are elements of H and indeed the sequence ¼ z z z 1, , , , 2 3 forms a basis in H (but in general it is not an orthonormal basis). This implies that any infinite sequence of polynomials ( ) such that ( ) = -Q z n deg 1 n for all = ¼ n 1, 2, 3 is a basis of H. In particular, the sequence ( ) ( ) ¼ P z P z , , 1 2 (of the previous section) is an orthonormal basis of H. Since this particular orthonormal basis satisfies (12), it follows that the kernel  of H is such that: That is, This is a key result that will now be used to compute the matrix -G 1 .

Explicit expression of the inverse matrix
for all z w , in the unit disc, then taking = = z w 0 yields that that is the element of the first row and first column of -G 1 . Hence To obtain ( ) -G 1 2,1 one can take differentiate (14) with respect to z at the point = = z w 0 and obtain In general for all = ¼ k 1, 2, 3, Since the matrix -G 1 is also hermitian, then clearly All other elements can also be obtained by taking successive partial derivatives with respect to z and w , that is, Since the Szegö function S(z) is analytic in the unit disc 1 and the Szegö condition implies that it does not vanish anywhere in D, then its inverse ( ) is also an analytic function in D. Assume then that Thus the -G 1 is fully determined by the coefficients of the function ( ) y z . Hence (this is the Whittle approximation of the matrix -G 1 ).

The upper left-hand bloc of -G 1
From the discussion above, it is clear that once the elements of the first row of -G 1 are determined, the other elements can easily be determined. Let ( ) be the upper left-hand n×n bloc of -G 1 (that is the bloc of -G 1 consisting of the first n rows and first n columns of -G 1 ). Then relation (18) implies that all the elements of ( ) -Ǵ n n 1 are explicitly determined by the firstn 1 derivatives of the function ( ) y z at the origin (or the first n coefficients of ( ) y z in its Taylor expansion at zero). In fact, a  a  a   a  a  a  a  a   a a a a  a  a a a  a  a  This is yields an LU decomposition of the matrix -G 1 . This decomposition also yields a system of orthonormal polynomials with respect to the weight function ( ) j t 1 . Indeed it is immediately seen that the system given by: is orthonormal on the unit circle with respect to the weight function ( ) j t 1 (the inverse of the spectral density of G).

Explicit expression for the inverse matrix
We can calculate the derivatives of ( ) y z at the origin with respect to z. Note that    Hence this yields using relation (19) that the upper left-hand 5×5 bloc of -G 1 (that is the bloc corresponding of the first 5 rows and 5 columns) is: One can compare that with the inverse of the bloc ( )Ǵ m m of the original covariance matrix G for m sufficiently large. Taking m=1000 yields the following: These two matrices are very close and the equality holds when  ¥ m .

Illustrating example: finite-diagonal Toeplitz matrices
In this section we now consider the particular case where the Toeplitx matrix G has only a finite number + m 2 1 of non-zero diagonals (such matrices are also called + m 2 1-diagonal Toeplitz matrices). It is also an interesting problem to compute the inverse of such matrices. Assume as previously that ( ) The corresponding spectral density if obviously given by The condition that G is positive definite implies that ( ) is integrable, then the corresponding Szegö function S(z) and its inverse ( ) y z exist as already discussed. The formulas of the previous section can be used to compute the elements ( ) a n and from which the inverse matrix -G 1 follows.

Case of infinite 3-diagonal Toeplitz matrix
This corresponds to m=1. Assume that the matrix G is defined by the function where q is a complex number such that | | < q 1 2.The corresponding spectral density is We can explicitly compute the corresponding function ( ) y z as follows. First consider the Szegö function S(z) associated to ( ) f t and write as in (10) and c k =0 for all = ¼ k 2, 3, Hence the Szegö function associated to j is: The corresponding function ψ is: That is: In the case where q is a real number, this function was obtained in [20] where it is expressed as: 1 where 4 and 2 sign . We can now explicitly compute the inverse of the matrix G. Using (21), it is clear that the Taylor expansion of ( ) y z is given by: Therefore the element of order ( ) k j , of the inverse matrix -G 1 (with  j k) is explicitly given by: This yields an explicit analytic expression of the inverse of the matrix G. Summarising we have the following result: Proposition 1. For any complex number q such that | | < q 1 2, the infinite matrix G given by is such that its inverse is the matrix -G 1 given explicitly by: To illustrate this result, take =q 1 5. Then the upper left-hand 5×5 bloc of -G 1 is: which is very close to the exact matrix.

Case of infinite 5-diagonal Toeplitz matrices
Assume G is given by the function: For example the first 5 coefficients in the Taylor expansion of ( ) y z at z=0 are" This yields that the the upper left-hand 5×5 bloc of -G 1 is: We can again compare with the inverse of the bloc ( )Ǵ m m of G for m sufficiently large. Taking m=100 yields the following: (This is exactly the same as the previous matrix and hence, the error of approximation is -10 6 .)

General case
In the general case of   From this explicit inverse Szegö function ( ) y z one can now compute every element of the inverse matrix -G 1 . For example using (15), we have that ( ) = -+ -G i 0.282433 0.183 806 . 1 10,9 We can also compute the first 10 coefficients ¼ a a a , , , The first 5×5 bloc of -G 1 is: Again this exact value is very close to the first 5×5 bloc of the inverse of ( )Ǵ m m for m large. In fact m=100 yields an approximation error -10 6 .