Embedding nonnegative definite Toeplitz matrices ... - Semantic Scholar
1206
IEEE TRANSACTIONS ON INFORMATION THEORY, VOL.
35, NO. 6, NOVEMBER 1989
Embedding Nonnegative Definite Toeplitz Matrices in Nonnegative Definite Circulant Matrices, with Application to Covariance Estimation AMIR DEMBO, COLIN L. MALLOWS, AND LARRY A. SHEPP
Ahtract -The class of nonnegative definite Toeplitz matrices that can be embedded in nonnegative definite circulant matrices of larger s u e is characterized. An equivalent characterization in terms of the spectrum of the underlying process is also presented, together with the corresponding extrema1 processes. It is shown that a given finite duration sequence p can be extended to be the covariance of a periodic stationary process whenever the Toeplib matrix R generated by this sequence is strictly positive definite. The sequence p =l,cosa,cos2a with (a/.) irrational, which has a unique nonperiodic extension as a covariance sequence, demonstrates that the strictness is needed. Our simple constructive proof supplies a bound on the abovementioned period in terms of the minimal eigenvalue of R, and also yields, under the same conditions, an extension of to covariances that eventually decay to zero. For the maximum likelihood estimate of the covariance of a stationary Gaussian process, we determine the extension length required for using the EM iterative algorithm.
I. INTRODUCTION
R
ECENTLY, the problem of estimating structured covariance matrices was addressed in [l], with particular emphasis on Toeplitz matrices (corresponding to stationary processes). Given a set of independent samples from a zero-mean multivariate Gaussian random process, the problem is to select a covariance matrix of specified structure that maximizes the probability of occurrence of these samples. Suppose that we have M independent samples x,, 1s m 5 M , from an N-variate Gaussian process. Their joint probability density is given by
transpose of vector or matrix), and
g ( S , R ) = -lnIRl-tr(R-'S).
(2)
Throughout, IAl denotes the determinant of the matrix, A and tr( A ) denotes the sum of the diagonal elements of the square matrix A . In a certain set I , of covariance matrices having a specified structure we look for the matrix R, which maximizes the joint probability density given in (1) or, alternatively, maximizes the function (2). This is a constrained optimization problem. denotes the class of nonnegative definite Hereafter symmetric matrices, with real entries. Unless stated explicitly, all matrices in this paper are in II. Let I I + C I I denotes the class of positive definite symmetric matrices; we shall use the notation R 2 0 for R E IX, and R > 0 for R E II+.The following existence theorem is due to Burg et al. [l].
n
Theorem 1: If S E II+ and I , taining at least one element of problem:
c II is a closed set con-
n+,then the optimization
max g ( S , R )
R
E IR
has positive definite solutions. Furthermore, if S then the unique solution of (3) is R = S .
(3) E I,,
Perhaps the most important estimation problem of this type, is that of a stationary Gaussian random process, for which I , is the class of nonnegative definite symmetric Toeplitz N X N dimensional matrices (denoted throughout where R is the covariance matrix, S A ( l / M ) Z ~ = l x m xis~ by T,). Usually, for M > N , the sample covariance matrix the sample covariance matrix (where the prime denotes the S is in Il+,and therefore (3) has positive definite solutions for I , = T,. However, as generally S T,,,, there appears to be no known analytical expression for these solutions Manuscript received June 19, 1987; revised February 28, 1989. and even the question of uniqueness remains open (cf., [2], A. Dembo was with AT&T Bell Laboratories, Murray Hill, NJ. He is 131). now with the Departments of Statistics and Mathematics, Stanford Let us denote by C, the class of N X N nonnegative University, Stanford. CA 94305-2099. C. L. Mallows and L. A. Shepp are with AT&T Bell Laboratories. definite symmetric circulant matrices, i.e., R E C, if Murray Hill. NJ 07974. Rkl=R(L-l,modN,O for every O s k , l s ( N - l ) , and R is IEEE Log Number 8931525. 0018-9448/89/llOO-1206$01.00 01989 IEEE
DLMBO f f U / . : EhZHl DIXNG NONNEGATIVE DEFINITE TOEPLITZ MATRICES
also a nonnegative definite symmetric matrix. For I , = C, there is a unique analytical solution of (3), as summarized by the following theorem from [3]. Theorem 2: If S E II+, then the unique solution of (3) for Z, = C,. is the matrix R* E C:, whose entries are RZ,=
c
1 -
N
(r
s,,,
O s k , I < (N-1)
(4)
1 ) = (/. p / ) n i o d N
where C s corresponds to the subset of positive definite matrices in C,. Likewise, TN+ denotes the subset of positive definite matrices in T N .For convenience we hereafter use the notation R* = circ(S) to denote the relation gwen in (4). Since C , T N . R* E TN and g ( S , R*) I maxRGTN g ( S , R ) . Thus R* given by (4) (whch is the periodogram estimator, e.g., [4]), is a simple but suboptimal solution of the original problem (solving (3) for I , = T,). The existence of these simple solutions to the estimation problem for I R = C l , (any L ) motivated us to use the EM iterative algorithm (cf., [5]) for solving our original problem (i.e., I , = T,). This is done by regarding the underlying covariance matrix R as the upper left N X N-dimen2ional block of an L x 1, dimensional covariance matrix R E C,, L 2 2 N - 1. Then, each iteration involves solving (3) for ZR = with s” (the sample covariance matrix corresponding to R ) being an extension of S , based on the current estimate of k . This algorithm was used in [lo], [ll]for the case of one sample path, i.e., rank S = l . While the case of one long sample path may be more common, positive definite solutions are not guaranteed (since S is singular). Furthermore, even if the true underlying stationary covariance matrix is positive definite, the maximum likelihood estimator for one sample (rank S = 1, and I , = T N )is almost surely singular (as shown in [13, theorem 31). As this estimation problem is ill-posed, it was suggested in [13] to restrict I , to the subset TNC, of N x N nonnegative symmetric Toeplitz matrices which can be embedded inside L X L nonnegative symmetric circulant matrices. Under this restriction the estimation problem is well-posed (cf., [13]). However, the dependence of L on the true covariance matrix is not explored there. For the sake of completeness we describe the EM algorithm in Section 111. We also show there that, for any sequence of estimates R(”) E TNC, generated by the EM algorithm, g ( S , R ( ” ) 2 ) g ( S , R(’lpl)), and when a limit point R* of the sequence R ( ” )belongs to TN+C,‘ (i.e., both R* > 0 and k* > 0), then R* is a stationary point of g(S9 R ) . These results (which are summarized in Proposition 4) are quite standard for the EM algorithm, and follow from the general results of [8]. Similar results were proven in [lo], [13], for the case of S singular. The main contributions of this work are presented in the next section. We characterize the set TNCi.in Proposition 3 and its corollaries and use these results to prove con-
cL
1207
structively that any positive definite Toeplitz matrix can be embedded in a positive definite circulant matrix of larger size. In particular, we have the following. Proposition 1: Let R E TN+; then R E TGC; provided that L > ( 2 N - 1 ) + N2R,/Xrnin(R). Thus when one has an a priori lower bound on the minimal eigenvalue of the true underlying covariance matrix, the value of L can be determined by (5) and the EM algorithm may then be used. Further, combining (5) with some of the results of [l]we obtain the following. Proposition 2: For any S E II’ and any R , E T;, let 0 < x1 < 1 < x, < 00. be the two solutions of g ( S , R , ) )- g ( S , S ) = l n x - x + l .
(6)
Then, the maximum likelihood (ML) estimator (the solution of (3) for I , = T,), belongs to T”C2, for every L such that
In particular, it is a limit point of an EM instance. Thus for S of full rank, eigendecomposition of S determines a value of L such that, at least for some initial conditions, the EM will converge to the ML estimator. We note that the sets T N ,C,, TNCL (as well as TN+,C; , TLC, and T;C; ) are convex sets. This is since R,, R , E rI (or II’) implies a , R , + a,R2 E rI(III’) for any a,, a 2 > 0, and both the Toeplitz and circulant structures are preserved under this operation. Also, symmetric Toeplitz and circulant matrices are determined by their first row, which we denote throughout by { r,, r,; . ., r N - , } .Thus T,, C,, TNCl-,and their strictly positive definite subsets are isomorphic to convex subsets of Iw., The convexity allows us to normalize r, to one without loss of generality. Let p A (p,,, p,: . . , p N p I ) be the normalized first row of R ( po 4 1); then the notation p E T,( T N C L ) , is equivalent to R E T,(C,, TNCi-,respectively). It is well-known (cf., [9] for example) that p E T, if and only if (iff) p is part of a covariance sequence, and that p E iff p is part of a periodic covariance sequence (with period N ) . From this point of view, Proposition 1 states that if p E T; then it has a periodic extension of some period L (as given in (5)). We note here that our constructive proof shows also the existence of an extension of p which is zero except maybe for the first L samples (with the same bound for L ) . However, there are singular matrices in T, that have no finite nonnegative definite circulant extension as, for example, the matrix
c,,
cN
J
1 cosa cos2a cosa 1 cosa cos2a cosa 1 with (a/..) irrational. We sketch the analogous results for complex-valued and continuous-time processes, as well as their interpretation
[
~
1208
IEEE TRANSACTIONS ON INFORMATIONTHEORY, VOL.
when the roles of time and frequency are interchanged. The last section is devoted to a short summary of our results.
11. THECHARACTERIZATION OF TNCLAND RELATED COVARIANCE EXTENSION PROBLEMS
35, NO. 6 , NOVEMBER 1989
Corollaries: C1) C, c TNC, iff N divides L.' C2) TNCL,c TNCL,if L, divides L2. T T:, NCL but U L > ( 2 N - 1 ) T N C L # c3) U L ~ ( ~ N - ~ ) 2 TN.
Proof: Note that all the geFrators of FN_arein FNCL iff N divides L, and those of TNCL,are in T,.CL2 iff L, We characterize below the convex subsets (of R ), E,, FN and TNCL.We then study their relationships and use divides L,. Since a convex cone is contained in another convex cone iff its generators are contained there, we get them to prove Propositions 1 and 2 stated above. Proposition 3: a) C , is the convex polygon in R gener- corollaries C1 and C2 from the isom_orphisms between C , ated by the 1N/2] 1 vertices { p('), p(l), . ., p [ N / 2 1 } ,with and ac,, and between T,C, and uTNCL, a > 0. The second part of C3 is an immediate consequence of p f ) A c 0 ~ ( 2 ~ i k / NO ) ,I k < ( N - 1 ) . b) TN is the convex hull of the one-dimensional mani- parts b) and c) of Proposition 3. Consider p(') for any irrational value of (t/lr). It is in f, but never in f N C L . fold p('), 0 I t I T, with p f ) = cos tk, 0 I k I ( N - 1). c) F,Cl. is the convex cone in R N , generated by Thus it generates a matrix in TN which is not in - . for the first part, without loss of { p(O),. . pLL/'l}, with pf) = cos ( 2 ~ i k / L ) , 0 I k 5 U L t ( 2 N - 1 ) T N C IAs generality we prove it for FNCLand !f;. Since T; is the ( N - 1). interior of the closed set FN, for any p E Fg there is a Proof: a) The eigenvalues of an N-dimensional circu- positive Euclidean diztance between p and the boundary lant matrix (denoted by A,; . A N - ' ) are exactly the of F,. The set FN - TNCLis the union of a finite number components of the N-dimensional discrete Fourier trans- of disjoint components, and for L + m , the diameters of form (DFT) of its first row, and it is well-known that CN these components approach zero uniformly. Thus the seis the subclass of circulant matrices with A, 2 0 , i = quence O;..,N-l and A , = A N - , , i = l , . - ., N - 1 (cf., [3]). The = tr R = N. constraint po = 1 implies that C~=~O'A, Thus EN= { p ( p = DFT-'(A), CLi'A, = N, A, = A N - , , I =I,. . . , N - 1, and A , 2 O}. The corresponding set of A's approaches zero as L+m. If there was a PE?;,' which is the convex polygon in R generated by {A('), A(1), . . . , was always in ( F, - TNCL),then d ( L ) would be bounded with A(/) = (N/2)(6,, + S ( N - , ) k ) for i 21, and A(): = away from zero yielding a contradiction. Thus Q.E.D. Nsok; 0 I k s N - 1. Therefore, is generated by U L 2 ( 2 N - ,,TNC, 2 DFT-'( A(,)), and the result follows from the DFT properNotes: 1) AS we have just proved, TN- U L 2 ( 2 N - 1 TNCL ties. contains the set of matrices having the first row ap('! with b) It is well-known that R E TN iff there exists a positive a > 0 and ( t / T ) irrational. They correspond to the atomic and even measure on the unit circle p ( e " ) , such that even measures p.(elB)= (a/2)6,, +(a/2)6,(-,,, where eJ' . . n is not a finite root of unity. 2) We essentially proved in C3 a weaker form of Propok = 0 ; . . , N 1, elk0 dp. (e"), rk = 2I.r - - n sition 1, namely, that R E T; implies R E T; C , for every L 2 Lo, with Lo sufficiently large. Since R E TG at least is the first row of R (cf., [9]). The constraint r, =1 implies whenever p ( e!') is not a finite combination of atomic even that p ( e @ )is a probability measure on the unit circle. The positive measures, we conclude that for every process with class of even probability measures on the unit circle is a a nonzero absolutely continuous component of its specconvex set, whose extremal points are the atomic measures trum, there is a value Lo(N), such that for every L 2 Lo(N ) p(eJO)=$Sot + $So(-t, for 0 s t I T (cf., [9]). Since the the N-dimensional covariance matrix R is the upper left elements of p are linear functionals-of the measure p.(*), block of a nonnegative definite L-dimensional circulant the extremal points that generate T,, are in correspon- matrix. dence with the extremal measures. Canonical representation of T,c,: Let us now interc) R E TNCldiff R E T,, and there exists I? E C,, such pret T,CI. in terms of the underlying process. First, we that r, = i = 0,. . ., N - 1. Since R is the upper left block shall construct a set of 1L/2] + 1 (scalar) processes correof I?, I? 2 0 implies R 2 0, and we can view TNCLas the sponding to the generators mentioned in part c) of Propoto R N obtained by omitting the last sition 3. The processes are X,") A cos(2mni/L)Y(') projection of ( L - N ) entries of the vector p. As both CL and its s i n ( 2 ~ n i / l ) Z ( ' ) with Y ( ' ) Z(') , two uncorrelated proprojection F,CL are convex polygons, the generators of cesses with zero mean and unit variance (as one can easily the latter set are merely the projections of the generators of check). Now, each p € F N C Lis part of a covariance sethe former set which are given in part a). Q.E.D.
+
e,
*,
cN
G-
/
e,
cl-
The following corollaries of Proposition 3, are of interest.
+
'Note. however. that even when N does not divide L , there may be element\ of C , that are alw in T , C,
DEMBO et
U/. :
1209
I,!v~BI.III)INGNONNEGATIVE DEFINITE TOEPLITZ .MATRICES
+
quence of a linear combination of these 1L / 2 ] 1 extremal processes, where Y ( ' ) Z , ( I ) are uncorrelated with Y ( ' ) ,Z ( ' ) for J f i. The proof of this statement is as follows: By Proposition 3 part c), p E TNCL means that ph = Cji/,21AlE(X ! ' ) X ! 2 x ) with A , 2 0. Define A, = a,' and use the absence of correlation between X i 1 ) and Xi') to conclude that p x = E [ X,,X,-,] with X , Li C ~ i ( 2 1 a l X i 1 ' . At this point we turn to the proof of Proposition 1 which is an immediate consequence of the following two lemmas. Lemmu 1: Let the sequence r1 -...
0
-1 :(
*
M
1 ' 1-- M
'
N-1 1--r N - lM
i
be the first row of R , E TL; then {r,, r l ; . ., r N P l }represents a matrix R E T$C:, for every L 2 ( 2 M - 1). Lemniu 2: If the minimal eigenvalue of the matrix R E T; I S greater than r,,N2/2(M - N ) , then R , E T;. Proof of Lemmu 1 . If R,,, E T;, then there is a positive even absolutely continuous measure fi( e''), such that FA (1/277)lTVe/" dfi(e/'), k E h, (cf., [ 9 ] ) and r, = (I - ( k / M ) ) F A , ( kI I ( N - 1).Furthermore, there is such a measure fi( el') with dfi( e / ' ) / d 8 > 0 (for example, the maximum entropy extension of the positive definite matrix R , defines one such measure). Define
where (M-1)
c
M,'(x)
(1- $)e-."20. x
=-(M-1)
Thus p( e'') is an absolutely continuous positive measure that is also even since both fi(e'?) and W,(x) are even measures. Furthermore, p( e'') is a positive measure whose first N moments are {r,, r 1 ; . - , r N P 1 } ,and ( d ~ / d 8 ) ( c / ( ' " " / ' -> ) )0 for every 0 I n I ( L - l ) , and L 2 ( 2 M - 1). As this is the DFT of size L of the zero extended sequence ((1- ( k / M ) ) F h } ~ I ! & l ) , we conclude that the circulant matrix k of which this is the first row is a positive definite extension of R . Q.E.D. Proof of Lemmu 2: roN
'
'
(N-1)
I
=1
as can bc easily verified. Since R E T N , lr,,l I r,, n = 0 , 1 ; . . , N - 1 , and thus
where A
=
R,w - R . Thus
and if min, A I ( R ) > r,,N2/2(M- N ) . then min, A I ( R ) > max, ( A , ( R , - R)I. Finally,
m i n A l ( R ) 2 - m a xI I A l ( R , - R ) I ,
minA,(R,)I
I
since min ( x ' R x + x'Ax) 2 min x'Rx \'\ =
1
\'\
=1
-
max (x'Axl. \'\
=I
Q.E.D. Remark: We have constructed an extension with length Lo(N ) such that
A
= - ((
-
of R ,
1)/2)
for every 8 E [ - 77, 771. This assures not only the periodic covariance extension of { r,;. ., r N P 1 }but , also the existence of an eventually zero covariance extension of this sequence (with the same value of Lo). We now use Proposition 1 to prove Proposition 2. Proof of Proposition 2: The existence theorem of Burg et al. (Theorem 1 in Section I) is proved in [1]based on the following upper bound on g(S, R ) : g(S,R ) - g ( S , S ) I l n x - x + l (8) whenever x 2 1 is an upper bound on Amin(S)/Amin(R), or x I 1 is a lower bound on A,, (S)/A,,( R ) . Assuming that the EM algorithm is initiated by R(O)= R , with R , E TN, let 0 < x1,as
where
and the 'matrices p ( " )of dimension ( L - N ) x N and fen) of dimension ( L - N ) x (L,- N ) are given by the following block partitioning of [ R ( " ) ] - ' :
Observing (12), it is easily verified that k(,)E II+ implies that F ( " )E n+,and then both Fen) and pen) are welldefined. Observing now (ll),it is easily verified that for and F ( ' I )E II+, iCn) E II+ as well. Therefore, SE applying Theorem 2 on (loa), we deduce that I?("+') E n + and is given by
n+
(13)
i ( n + l ) = circ(S(n))
where here circ( - ) is of course done modulo L, instead of modulo N . Thus (11)-(13) together with the definition of 2):, fully characterize the EM iteration. Fast implementation of (ll), (1 3) involves ( M 2 ) DFT's of size L (wh;_Chcan also give [ I?("+')] p 1 instead of i?("+')). Evaluating T ( " )involves the inversion of a positive definite Toeplitz matrix of dimension ( L - N ) , and evaluating 2):, m = 1; M involves solving M sets of ( L - N ) linear equations having the same positive definite Toeplitz matrix of coefficients. There are many efficient algorithms for this task, and-we shall omit the details here (cf., [6],[7]).We note that T'") is in general not circulant and so the numerical difficulties reported for example in [12] are irrelevant. The convergence properties of the EM algorithm which were discussed in [8] guarantee the following. Proposition 4: If the initial estimate Reo) E TSC; (and S E II+ ), then the sequence of upper left submatrices { R(O),R('),. . ., R'"), . . . } generated by the EM algorithm has the following properties:
+
-
ln,
where y,,, 1I m 5 M are vectors in R L missing data in the given samples, and
we obtain the following form for a typical iteration
a ,
R'") E T;C:, n = 0,1,. . . ; 2 g ( S , R(,-')), with equality iff R(")=
a) b)
g( S , I?'''))
c)
the sequence { R(O),R('),. . ., R("),. . . } is contained in a compact subset of TSC, and thus has at least
R(jI-1).
e,
DEMHO
et
U / .:
I.MHI3)DING NONNEGATIVE DEFINITE TOEPLITZ MATRICES
1211
trix R as the upper left block in a larger nonnegative definite circulant matrix I?. For N that divides L , our initial estimate has a positive definite extension to a circulant E(()).The EM algorithm then generates a sequence Proof: a) As R(O)E TGC‘L, kc’)E l l + . We observed of positive definite estimates { R(O); . ., R ( ” ) ., . . ) with earlier that then io’) E C: for n = 0,1,2, . . . and thus monotonically increasing likelihood. This sequence has R ( ” ) E T C C T ,n = 0 , 1 ; . . . limit points, all of which are either stationary points of the b) The monotone increase of g( S, R ) , is one of the basic likelihood, or on the boundary of TNC, the subclass of properties of the EM algorithm (cf., [5]). It follows from nonnegative definite Toeplitz matrices which can be emJensen’s inequality, that leads in our case to bedded in nonnegative definite Circulant matrices of dimension L. Furthermore, if the likelihood function is g ( S , R ) - g ( S, R ( ” - l ) ) unimodal on T N C , then convergence to the global maxi2 g ( i ( t 1 - 1 ) , i)- g ( i ( “ - l ) , j 0 1 - 1 ) 1 mum is guaranteed. (14) The characterization of the class of matrices T,C,. and since E ( ” ) is the unique maximum of g(S^(”-’),I?) in which has its own theoretical importance was also preC,, the desired inequality follows and is strict whenever sented, together with the extrema1 processes, and a canoniI?(’!)# (i.e., R ( ” )# R(“-’), since g depends only on cal representation in terms of these processes. It was also R!). shown that T , 3 3 U I . 2 ( Z N - I ) T N3C L 3 TG (where 3 3 c) We apply [8, Theorem 21 to our case. For any E(’) E denotes ‘strict inclusion’). Thus there are degenerate proII+, R(’) E n+,and therefore (S E II’ ) the level set cesses for which there is no circulant embedding, but for {RI g(S, R ) 2 g ( S , R“))) is a compact set contained in every process whose spectrum has a non-zero absolutely II+ (the proof of this claim appears in [l]and relies on continuous component, (or, whenever R > 0) there exist (8)). Let Q denote this set. Since 52 c II+,g(S, R ) is con- Circulant embeddings for large enough values of L. Using tinuous and differentiable in 52 and is bounded above our concrete bounds on this value, and the analysis of [l], there (since 52 is compact). Since TNC, is closed, TNCLn 52 we gave a bound on L , in terms of the eigendecomposition is a compact subset of TGC, ip which { R(O),-. ., R‘”),. . . } of S (the sample covariance matrix) and the initial estiis containe.. Since g(S(”),R ) is continuous both with mate R“), such that the ML estimate is in T$CCL+,and thus whenever i(“) E H+, respect to R and with respect to is a limit point of an EM iteration. it follows from [8] that any limit point of { R(O), Our motivation led to the presentation of TNCL as a . . . , R‘”),. . . ) in TN+CL is a stationary point of g ( S , R ) . class of matrices. However, another view of our results is Q.E.D. possible. Suppose we are given a finite duration sequence Remark: If we choose L according to Proposition 2 (see p = { p o , p l , . . . , p N P I } ; then p is the first row of (7)), then following the arguments we used in the proof of an element of TNCl., iff it can be extended to be the Propositions 1 and 2, we may conclude that 52 is a com- covariance of a periodic stationary process. Since pact subset of T;CT (i.e., the minimal eigenvalue of R is U I. (ZN-l)TNC,-3 3 TN+,any sequence p that generates a bounded away from zero for any R E a). Therefore, we strictly positive definite matrix, has such an extension. the strictness is can strengthen part c) of Proposition 4 and claim that all However, since T%3 3 Ur.2(ZN-1)TNC12, limit points of the sequence R(“)are stationary points of g indeed required for this property. In view of our remark following the proof of Proposition in TNCT. Convergence of the EM to the global maximum of 1, whenever p generates a strictly positive definite matrix, g ( S , R ) in TNCl, is guaranteed whenever g(S, R ) is uni- there also exists an extension of p which is an eventually modal in this region. However, the question whether or not zero covariance sequence. However, how to find the minimum periodic or minimum eventually zero extension reg(S, R ) is unimodal in TNCLis still open. We note that by restricting L to be of the form k N , k mains open. integer, Corollary C1 guarantees that the periodogram estimator R* given by (4) is in TGC;. Thus R* generates a ACKNOWLEDGMENT natural initial estimate R(’) for the EM. The authors thank H. J. Landau for many helpful discussions. They also thank D. J. Thomson for careful review IV. CONCLUSION of a preliminary version of this work. We have presented an EM iterative algorithm for the ML estimate of the covariance matrix of a zero-mean REFERENCES stationary Gaussian random process. The initial estimate is obtained by applying the periodogram method to the J. P. Burg. D. G. Luenberger. and D. L. Wegner. “Estimation of covariance estimation, getting a simple suboptimal estistructured covariance matrices.” Pro
IEEE TRANSACTIONS ON INFORMATION THEORY, VOL.
35, NO. 6, NOVEMBER 1989
Embedding Nonnegative Definite Toeplitz Matrices in Nonnegative Definite Circulant Matrices, with Application to Covariance Estimation AMIR DEMBO, COLIN L. MALLOWS, AND LARRY A. SHEPP
Ahtract -The class of nonnegative definite Toeplitz matrices that can be embedded in nonnegative definite circulant matrices of larger s u e is characterized. An equivalent characterization in terms of the spectrum of the underlying process is also presented, together with the corresponding extrema1 processes. It is shown that a given finite duration sequence p can be extended to be the covariance of a periodic stationary process whenever the Toeplib matrix R generated by this sequence is strictly positive definite. The sequence p =l,cosa,cos2a with (a/.) irrational, which has a unique nonperiodic extension as a covariance sequence, demonstrates that the strictness is needed. Our simple constructive proof supplies a bound on the abovementioned period in terms of the minimal eigenvalue of R, and also yields, under the same conditions, an extension of to covariances that eventually decay to zero. For the maximum likelihood estimate of the covariance of a stationary Gaussian process, we determine the extension length required for using the EM iterative algorithm.
I. INTRODUCTION
R
ECENTLY, the problem of estimating structured covariance matrices was addressed in [l], with particular emphasis on Toeplitz matrices (corresponding to stationary processes). Given a set of independent samples from a zero-mean multivariate Gaussian random process, the problem is to select a covariance matrix of specified structure that maximizes the probability of occurrence of these samples. Suppose that we have M independent samples x,, 1s m 5 M , from an N-variate Gaussian process. Their joint probability density is given by
transpose of vector or matrix), and
g ( S , R ) = -lnIRl-tr(R-'S).
(2)
Throughout, IAl denotes the determinant of the matrix, A and tr( A ) denotes the sum of the diagonal elements of the square matrix A . In a certain set I , of covariance matrices having a specified structure we look for the matrix R, which maximizes the joint probability density given in (1) or, alternatively, maximizes the function (2). This is a constrained optimization problem. denotes the class of nonnegative definite Hereafter symmetric matrices, with real entries. Unless stated explicitly, all matrices in this paper are in II. Let I I + C I I denotes the class of positive definite symmetric matrices; we shall use the notation R 2 0 for R E IX, and R > 0 for R E II+.The following existence theorem is due to Burg et al. [l].
n
Theorem 1: If S E II+ and I , taining at least one element of problem:
c II is a closed set con-
n+,then the optimization
max g ( S , R )
R
E IR
has positive definite solutions. Furthermore, if S then the unique solution of (3) is R = S .
(3) E I,,
Perhaps the most important estimation problem of this type, is that of a stationary Gaussian random process, for which I , is the class of nonnegative definite symmetric Toeplitz N X N dimensional matrices (denoted throughout where R is the covariance matrix, S A ( l / M ) Z ~ = l x m xis~ by T,). Usually, for M > N , the sample covariance matrix the sample covariance matrix (where the prime denotes the S is in Il+,and therefore (3) has positive definite solutions for I , = T,. However, as generally S T,,,, there appears to be no known analytical expression for these solutions Manuscript received June 19, 1987; revised February 28, 1989. and even the question of uniqueness remains open (cf., [2], A. Dembo was with AT&T Bell Laboratories, Murray Hill, NJ. He is 131). now with the Departments of Statistics and Mathematics, Stanford Let us denote by C, the class of N X N nonnegative University, Stanford. CA 94305-2099. C. L. Mallows and L. A. Shepp are with AT&T Bell Laboratories. definite symmetric circulant matrices, i.e., R E C, if Murray Hill. NJ 07974. Rkl=R(L-l,modN,O for every O s k , l s ( N - l ) , and R is IEEE Log Number 8931525. 0018-9448/89/llOO-1206$01.00 01989 IEEE
DLMBO f f U / . : EhZHl DIXNG NONNEGATIVE DEFINITE TOEPLITZ MATRICES
also a nonnegative definite symmetric matrix. For I , = C, there is a unique analytical solution of (3), as summarized by the following theorem from [3]. Theorem 2: If S E II+, then the unique solution of (3) for Z, = C,. is the matrix R* E C:, whose entries are RZ,=
c
1 -
N
(r
s,,,
O s k , I < (N-1)
(4)
1 ) = (/. p / ) n i o d N
where C s corresponds to the subset of positive definite matrices in C,. Likewise, TN+ denotes the subset of positive definite matrices in T N .For convenience we hereafter use the notation R* = circ(S) to denote the relation gwen in (4). Since C , T N . R* E TN and g ( S , R*) I maxRGTN g ( S , R ) . Thus R* given by (4) (whch is the periodogram estimator, e.g., [4]), is a simple but suboptimal solution of the original problem (solving (3) for I , = T,). The existence of these simple solutions to the estimation problem for I R = C l , (any L ) motivated us to use the EM iterative algorithm (cf., [5]) for solving our original problem (i.e., I , = T,). This is done by regarding the underlying covariance matrix R as the upper left N X N-dimen2ional block of an L x 1, dimensional covariance matrix R E C,, L 2 2 N - 1. Then, each iteration involves solving (3) for ZR = with s” (the sample covariance matrix corresponding to R ) being an extension of S , based on the current estimate of k . This algorithm was used in [lo], [ll]for the case of one sample path, i.e., rank S = l . While the case of one long sample path may be more common, positive definite solutions are not guaranteed (since S is singular). Furthermore, even if the true underlying stationary covariance matrix is positive definite, the maximum likelihood estimator for one sample (rank S = 1, and I , = T N )is almost surely singular (as shown in [13, theorem 31). As this estimation problem is ill-posed, it was suggested in [13] to restrict I , to the subset TNC, of N x N nonnegative symmetric Toeplitz matrices which can be embedded inside L X L nonnegative symmetric circulant matrices. Under this restriction the estimation problem is well-posed (cf., [13]). However, the dependence of L on the true covariance matrix is not explored there. For the sake of completeness we describe the EM algorithm in Section 111. We also show there that, for any sequence of estimates R(”) E TNC, generated by the EM algorithm, g ( S , R ( ” ) 2 ) g ( S , R(’lpl)), and when a limit point R* of the sequence R ( ” )belongs to TN+C,‘ (i.e., both R* > 0 and k* > 0), then R* is a stationary point of g(S9 R ) . These results (which are summarized in Proposition 4) are quite standard for the EM algorithm, and follow from the general results of [8]. Similar results were proven in [lo], [13], for the case of S singular. The main contributions of this work are presented in the next section. We characterize the set TNCi.in Proposition 3 and its corollaries and use these results to prove con-
cL
1207
structively that any positive definite Toeplitz matrix can be embedded in a positive definite circulant matrix of larger size. In particular, we have the following. Proposition 1: Let R E TN+; then R E TGC; provided that L > ( 2 N - 1 ) + N2R,/Xrnin(R). Thus when one has an a priori lower bound on the minimal eigenvalue of the true underlying covariance matrix, the value of L can be determined by (5) and the EM algorithm may then be used. Further, combining (5) with some of the results of [l]we obtain the following. Proposition 2: For any S E II’ and any R , E T;, let 0 < x1 < 1 < x, < 00. be the two solutions of g ( S , R , ) )- g ( S , S ) = l n x - x + l .
(6)
Then, the maximum likelihood (ML) estimator (the solution of (3) for I , = T,), belongs to T”C2, for every L such that
In particular, it is a limit point of an EM instance. Thus for S of full rank, eigendecomposition of S determines a value of L such that, at least for some initial conditions, the EM will converge to the ML estimator. We note that the sets T N ,C,, TNCL (as well as TN+,C; , TLC, and T;C; ) are convex sets. This is since R,, R , E rI (or II’) implies a , R , + a,R2 E rI(III’) for any a,, a 2 > 0, and both the Toeplitz and circulant structures are preserved under this operation. Also, symmetric Toeplitz and circulant matrices are determined by their first row, which we denote throughout by { r,, r,; . ., r N - , } .Thus T,, C,, TNCl-,and their strictly positive definite subsets are isomorphic to convex subsets of Iw., The convexity allows us to normalize r, to one without loss of generality. Let p A (p,,, p,: . . , p N p I ) be the normalized first row of R ( po 4 1); then the notation p E T,( T N C L ) , is equivalent to R E T,(C,, TNCi-,respectively). It is well-known (cf., [9] for example) that p E T, if and only if (iff) p is part of a covariance sequence, and that p E iff p is part of a periodic covariance sequence (with period N ) . From this point of view, Proposition 1 states that if p E T; then it has a periodic extension of some period L (as given in (5)). We note here that our constructive proof shows also the existence of an extension of p which is zero except maybe for the first L samples (with the same bound for L ) . However, there are singular matrices in T, that have no finite nonnegative definite circulant extension as, for example, the matrix
c,,
cN
J
1 cosa cos2a cosa 1 cosa cos2a cosa 1 with (a/..) irrational. We sketch the analogous results for complex-valued and continuous-time processes, as well as their interpretation
[
~
1208
IEEE TRANSACTIONS ON INFORMATIONTHEORY, VOL.
when the roles of time and frequency are interchanged. The last section is devoted to a short summary of our results.
11. THECHARACTERIZATION OF TNCLAND RELATED COVARIANCE EXTENSION PROBLEMS
35, NO. 6 , NOVEMBER 1989
Corollaries: C1) C, c TNC, iff N divides L.' C2) TNCL,c TNCL,if L, divides L2. T T:, NCL but U L > ( 2 N - 1 ) T N C L # c3) U L ~ ( ~ N - ~ ) 2 TN.
Proof: Note that all the geFrators of FN_arein FNCL iff N divides L, and those of TNCL,are in T,.CL2 iff L, We characterize below the convex subsets (of R ), E,, FN and TNCL.We then study their relationships and use divides L,. Since a convex cone is contained in another convex cone iff its generators are contained there, we get them to prove Propositions 1 and 2 stated above. Proposition 3: a) C , is the convex polygon in R gener- corollaries C1 and C2 from the isom_orphisms between C , ated by the 1N/2] 1 vertices { p('), p(l), . ., p [ N / 2 1 } ,with and ac,, and between T,C, and uTNCL, a > 0. The second part of C3 is an immediate consequence of p f ) A c 0 ~ ( 2 ~ i k / NO ) ,I k < ( N - 1 ) . b) TN is the convex hull of the one-dimensional mani- parts b) and c) of Proposition 3. Consider p(') for any irrational value of (t/lr). It is in f, but never in f N C L . fold p('), 0 I t I T, with p f ) = cos tk, 0 I k I ( N - 1). c) F,Cl. is the convex cone in R N , generated by Thus it generates a matrix in TN which is not in - . for the first part, without loss of { p(O),. . pLL/'l}, with pf) = cos ( 2 ~ i k / L ) , 0 I k 5 U L t ( 2 N - 1 ) T N C IAs generality we prove it for FNCLand !f;. Since T; is the ( N - 1). interior of the closed set FN, for any p E Fg there is a Proof: a) The eigenvalues of an N-dimensional circu- positive Euclidean diztance between p and the boundary lant matrix (denoted by A,; . A N - ' ) are exactly the of F,. The set FN - TNCLis the union of a finite number components of the N-dimensional discrete Fourier trans- of disjoint components, and for L + m , the diameters of form (DFT) of its first row, and it is well-known that CN these components approach zero uniformly. Thus the seis the subclass of circulant matrices with A, 2 0 , i = quence O;..,N-l and A , = A N - , , i = l , . - ., N - 1 (cf., [3]). The = tr R = N. constraint po = 1 implies that C~=~O'A, Thus EN= { p ( p = DFT-'(A), CLi'A, = N, A, = A N - , , I =I,. . . , N - 1, and A , 2 O}. The corresponding set of A's approaches zero as L+m. If there was a PE?;,' which is the convex polygon in R generated by {A('), A(1), . . . , was always in ( F, - TNCL),then d ( L ) would be bounded with A(/) = (N/2)(6,, + S ( N - , ) k ) for i 21, and A(): = away from zero yielding a contradiction. Thus Q.E.D. Nsok; 0 I k s N - 1. Therefore, is generated by U L 2 ( 2 N - ,,TNC, 2 DFT-'( A(,)), and the result follows from the DFT properNotes: 1) AS we have just proved, TN- U L 2 ( 2 N - 1 TNCL ties. contains the set of matrices having the first row ap('! with b) It is well-known that R E TN iff there exists a positive a > 0 and ( t / T ) irrational. They correspond to the atomic and even measure on the unit circle p ( e " ) , such that even measures p.(elB)= (a/2)6,, +(a/2)6,(-,,, where eJ' . . n is not a finite root of unity. 2) We essentially proved in C3 a weaker form of Propok = 0 ; . . , N 1, elk0 dp. (e"), rk = 2I.r - - n sition 1, namely, that R E T; implies R E T; C , for every L 2 Lo, with Lo sufficiently large. Since R E TG at least is the first row of R (cf., [9]). The constraint r, =1 implies whenever p ( e!') is not a finite combination of atomic even that p ( e @ )is a probability measure on the unit circle. The positive measures, we conclude that for every process with class of even probability measures on the unit circle is a a nonzero absolutely continuous component of its specconvex set, whose extremal points are the atomic measures trum, there is a value Lo(N), such that for every L 2 Lo(N ) p(eJO)=$Sot + $So(-t, for 0 s t I T (cf., [9]). Since the the N-dimensional covariance matrix R is the upper left elements of p are linear functionals-of the measure p.(*), block of a nonnegative definite L-dimensional circulant the extremal points that generate T,, are in correspon- matrix. dence with the extremal measures. Canonical representation of T,c,: Let us now interc) R E TNCldiff R E T,, and there exists I? E C,, such pret T,CI. in terms of the underlying process. First, we that r, = i = 0,. . ., N - 1. Since R is the upper left block shall construct a set of 1L/2] + 1 (scalar) processes correof I?, I? 2 0 implies R 2 0, and we can view TNCLas the sponding to the generators mentioned in part c) of Propoto R N obtained by omitting the last sition 3. The processes are X,") A cos(2mni/L)Y(') projection of ( L - N ) entries of the vector p. As both CL and its s i n ( 2 ~ n i / l ) Z ( ' ) with Y ( ' ) Z(') , two uncorrelated proprojection F,CL are convex polygons, the generators of cesses with zero mean and unit variance (as one can easily the latter set are merely the projections of the generators of check). Now, each p € F N C Lis part of a covariance sethe former set which are given in part a). Q.E.D.
+
e,
*,
cN
G-
/
e,
cl-
The following corollaries of Proposition 3, are of interest.
+
'Note. however. that even when N does not divide L , there may be element\ of C , that are alw in T , C,
DEMBO et
U/. :
1209
I,!v~BI.III)INGNONNEGATIVE DEFINITE TOEPLITZ .MATRICES
+
quence of a linear combination of these 1L / 2 ] 1 extremal processes, where Y ( ' ) Z , ( I ) are uncorrelated with Y ( ' ) ,Z ( ' ) for J f i. The proof of this statement is as follows: By Proposition 3 part c), p E TNCL means that ph = Cji/,21AlE(X ! ' ) X ! 2 x ) with A , 2 0. Define A, = a,' and use the absence of correlation between X i 1 ) and Xi') to conclude that p x = E [ X,,X,-,] with X , Li C ~ i ( 2 1 a l X i 1 ' . At this point we turn to the proof of Proposition 1 which is an immediate consequence of the following two lemmas. Lemmu 1: Let the sequence r1 -...
0
-1 :(
*
M
1 ' 1-- M
'
N-1 1--r N - lM
i
be the first row of R , E TL; then {r,, r l ; . ., r N P l }represents a matrix R E T$C:, for every L 2 ( 2 M - 1). Lemniu 2: If the minimal eigenvalue of the matrix R E T; I S greater than r,,N2/2(M - N ) , then R , E T;. Proof of Lemmu 1 . If R,,, E T;, then there is a positive even absolutely continuous measure fi( e''), such that FA (1/277)lTVe/" dfi(e/'), k E h, (cf., [ 9 ] ) and r, = (I - ( k / M ) ) F A , ( kI I ( N - 1).Furthermore, there is such a measure fi( el') with dfi( e / ' ) / d 8 > 0 (for example, the maximum entropy extension of the positive definite matrix R , defines one such measure). Define
where (M-1)
c
M,'(x)
(1- $)e-."20. x
=-(M-1)
Thus p( e'') is an absolutely continuous positive measure that is also even since both fi(e'?) and W,(x) are even measures. Furthermore, p( e'') is a positive measure whose first N moments are {r,, r 1 ; . - , r N P 1 } ,and ( d ~ / d 8 ) ( c / ( ' " " / ' -> ) )0 for every 0 I n I ( L - l ) , and L 2 ( 2 M - 1). As this is the DFT of size L of the zero extended sequence ((1- ( k / M ) ) F h } ~ I ! & l ) , we conclude that the circulant matrix k of which this is the first row is a positive definite extension of R . Q.E.D. Proof of Lemmu 2: roN
'
'
(N-1)
I
=1
as can bc easily verified. Since R E T N , lr,,l I r,, n = 0 , 1 ; . . , N - 1 , and thus
where A
=
R,w - R . Thus
and if min, A I ( R ) > r,,N2/2(M- N ) . then min, A I ( R ) > max, ( A , ( R , - R)I. Finally,
m i n A l ( R ) 2 - m a xI I A l ( R , - R ) I ,
minA,(R,)I
I
since min ( x ' R x + x'Ax) 2 min x'Rx \'\ =
1
\'\
=1
-
max (x'Axl. \'\
=I
Q.E.D. Remark: We have constructed an extension with length Lo(N ) such that
A
= - ((
-
of R ,
1)/2)
for every 8 E [ - 77, 771. This assures not only the periodic covariance extension of { r,;. ., r N P 1 }but , also the existence of an eventually zero covariance extension of this sequence (with the same value of Lo). We now use Proposition 1 to prove Proposition 2. Proof of Proposition 2: The existence theorem of Burg et al. (Theorem 1 in Section I) is proved in [1]based on the following upper bound on g(S, R ) : g(S,R ) - g ( S , S ) I l n x - x + l (8) whenever x 2 1 is an upper bound on Amin(S)/Amin(R), or x I 1 is a lower bound on A,, (S)/A,,( R ) . Assuming that the EM algorithm is initiated by R(O)= R , with R , E TN, let 0 < x1,as
where
and the 'matrices p ( " )of dimension ( L - N ) x N and fen) of dimension ( L - N ) x (L,- N ) are given by the following block partitioning of [ R ( " ) ] - ' :
Observing (12), it is easily verified that k(,)E II+ implies that F ( " )E n+,and then both Fen) and pen) are welldefined. Observing now (ll),it is easily verified that for and F ( ' I )E II+, iCn) E II+ as well. Therefore, SE applying Theorem 2 on (loa), we deduce that I?("+') E n + and is given by
n+
(13)
i ( n + l ) = circ(S(n))
where here circ( - ) is of course done modulo L, instead of modulo N . Thus (11)-(13) together with the definition of 2):, fully characterize the EM iteration. Fast implementation of (ll), (1 3) involves ( M 2 ) DFT's of size L (wh;_Chcan also give [ I?("+')] p 1 instead of i?("+')). Evaluating T ( " )involves the inversion of a positive definite Toeplitz matrix of dimension ( L - N ) , and evaluating 2):, m = 1; M involves solving M sets of ( L - N ) linear equations having the same positive definite Toeplitz matrix of coefficients. There are many efficient algorithms for this task, and-we shall omit the details here (cf., [6],[7]).We note that T'") is in general not circulant and so the numerical difficulties reported for example in [12] are irrelevant. The convergence properties of the EM algorithm which were discussed in [8] guarantee the following. Proposition 4: If the initial estimate Reo) E TSC; (and S E II+ ), then the sequence of upper left submatrices { R(O),R('),. . ., R'"), . . . } generated by the EM algorithm has the following properties:
+
-
ln,
where y,,, 1I m 5 M are vectors in R L missing data in the given samples, and
we obtain the following form for a typical iteration
a ,
R'") E T;C:, n = 0,1,. . . ; 2 g ( S , R(,-')), with equality iff R(")=
a) b)
g( S , I?'''))
c)
the sequence { R(O),R('),. . ., R("),. . . } is contained in a compact subset of TSC, and thus has at least
R(jI-1).
e,
DEMHO
et
U / .:
I.MHI3)DING NONNEGATIVE DEFINITE TOEPLITZ MATRICES
1211
trix R as the upper left block in a larger nonnegative definite circulant matrix I?. For N that divides L , our initial estimate has a positive definite extension to a circulant E(()).The EM algorithm then generates a sequence Proof: a) As R(O)E TGC‘L, kc’)E l l + . We observed of positive definite estimates { R(O); . ., R ( ” ) ., . . ) with earlier that then io’) E C: for n = 0,1,2, . . . and thus monotonically increasing likelihood. This sequence has R ( ” ) E T C C T ,n = 0 , 1 ; . . . limit points, all of which are either stationary points of the b) The monotone increase of g( S, R ) , is one of the basic likelihood, or on the boundary of TNC, the subclass of properties of the EM algorithm (cf., [5]). It follows from nonnegative definite Toeplitz matrices which can be emJensen’s inequality, that leads in our case to bedded in nonnegative definite Circulant matrices of dimension L. Furthermore, if the likelihood function is g ( S , R ) - g ( S, R ( ” - l ) ) unimodal on T N C , then convergence to the global maxi2 g ( i ( t 1 - 1 ) , i)- g ( i ( “ - l ) , j 0 1 - 1 ) 1 mum is guaranteed. (14) The characterization of the class of matrices T,C,. and since E ( ” ) is the unique maximum of g(S^(”-’),I?) in which has its own theoretical importance was also preC,, the desired inequality follows and is strict whenever sented, together with the extrema1 processes, and a canoniI?(’!)# (i.e., R ( ” )# R(“-’), since g depends only on cal representation in terms of these processes. It was also R!). shown that T , 3 3 U I . 2 ( Z N - I ) T N3C L 3 TG (where 3 3 c) We apply [8, Theorem 21 to our case. For any E(’) E denotes ‘strict inclusion’). Thus there are degenerate proII+, R(’) E n+,and therefore (S E II’ ) the level set cesses for which there is no circulant embedding, but for {RI g(S, R ) 2 g ( S , R“))) is a compact set contained in every process whose spectrum has a non-zero absolutely II+ (the proof of this claim appears in [l]and relies on continuous component, (or, whenever R > 0) there exist (8)). Let Q denote this set. Since 52 c II+,g(S, R ) is con- Circulant embeddings for large enough values of L. Using tinuous and differentiable in 52 and is bounded above our concrete bounds on this value, and the analysis of [l], there (since 52 is compact). Since TNC, is closed, TNCLn 52 we gave a bound on L , in terms of the eigendecomposition is a compact subset of TGC, ip which { R(O),-. ., R‘”),. . . } of S (the sample covariance matrix) and the initial estiis containe.. Since g(S(”),R ) is continuous both with mate R“), such that the ML estimate is in T$CCL+,and thus whenever i(“) E H+, respect to R and with respect to is a limit point of an EM iteration. it follows from [8] that any limit point of { R(O), Our motivation led to the presentation of TNCL as a . . . , R‘”),. . . ) in TN+CL is a stationary point of g ( S , R ) . class of matrices. However, another view of our results is Q.E.D. possible. Suppose we are given a finite duration sequence Remark: If we choose L according to Proposition 2 (see p = { p o , p l , . . . , p N P I } ; then p is the first row of (7)), then following the arguments we used in the proof of an element of TNCl., iff it can be extended to be the Propositions 1 and 2, we may conclude that 52 is a com- covariance of a periodic stationary process. Since pact subset of T;CT (i.e., the minimal eigenvalue of R is U I. (ZN-l)TNC,-3 3 TN+,any sequence p that generates a bounded away from zero for any R E a). Therefore, we strictly positive definite matrix, has such an extension. the strictness is can strengthen part c) of Proposition 4 and claim that all However, since T%3 3 Ur.2(ZN-1)TNC12, limit points of the sequence R(“)are stationary points of g indeed required for this property. In view of our remark following the proof of Proposition in TNCT. Convergence of the EM to the global maximum of 1, whenever p generates a strictly positive definite matrix, g ( S , R ) in TNCl, is guaranteed whenever g(S, R ) is uni- there also exists an extension of p which is an eventually modal in this region. However, the question whether or not zero covariance sequence. However, how to find the minimum periodic or minimum eventually zero extension reg(S, R ) is unimodal in TNCLis still open. We note that by restricting L to be of the form k N , k mains open. integer, Corollary C1 guarantees that the periodogram estimator R* given by (4) is in TGC;. Thus R* generates a ACKNOWLEDGMENT natural initial estimate R(’) for the EM. The authors thank H. J. Landau for many helpful discussions. They also thank D. J. Thomson for careful review IV. CONCLUSION of a preliminary version of this work. We have presented an EM iterative algorithm for the ML estimate of the covariance matrix of a zero-mean REFERENCES stationary Gaussian random process. The initial estimate is obtained by applying the periodogram method to the J. P. Burg. D. G. Luenberger. and D. L. Wegner. “Estimation of covariance estimation, getting a simple suboptimal estistructured covariance matrices.” Pro
