IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 34, NO.

7, JULY 1989 783

conditioned on Y,, Dl is a Gaussian random variable with mean El and REFERENCES

variance a,,
which are calculated sequentially u$ng the Kalman filtering D. L. Alspach and H. W. Sorenson, “Nonlinear Bayesian estimation using
equations over one step, with = 0 and @,I, = R, i.e., Gaussian sum approximations,” IEEE Trans. Automat. Contr., vol. AC-17, pp.
439-448, Aug. 1972.
a) P ; = R ( R + Q ) - ‘ Y , B. D. 0. Anderson and J. B. Moore, Optimal Filtering. Englewood Cliffs, NJ:
Prentice-Hall, 1979.
b) B, = R - R ( R + Q ) - ‘ R . (e.20) P. Bremand and J. Van Schuppen, “Discrete time stochastic systems, Parts I and
11,” Washington Univ., St. Louis, MO, SSM Rep. 7603a and b, 1977.
Therefore, C. D. de Benito, “Nonlinear filtering for discrete time systems,” Ph.D. thesis
proposal, Dep. Syst. Eng., Case Western Reserve Univ., Cleveland, OH, Dec.
1983.
-, “Nonlinear filtering for discrete time systems,” Ph.D. dissertation, Case
Westem Reserve Univ., Cleveland, OH, 1985.
C. D. de Benito and K. A. Loparo, “Nonlinear filtering for discrete time
systems,” submitted for publication.
G. B. Di Massi and W. J. Runggaldier, “On measure transformation for combined
fitering and parameter estimation in discrete time,’’ Syst. Contr. Lett., vol. 2,
pp. 51-62, July 1982.
(e.21)
I. V. Girsanov, “On transforming a certain class of stochastic processes by
absolutely continuous substitution of measure,” Theory of Probability and its
Then the conditional characteristic function is given by the following: Applications, vol. V, no. 3, pp. 285-301, 1960.
R. E. Kalman, “A new approach to linear fitering and prediction problems,” J .
Basic Eng., Trans. ASME, pp. 3 5 4 5 , Mar. 1960.
B. D. Koopman, “On distributions admitting a sufficient statistic,” Trans. Amer.
Math. Soc., vol. 39, pp. 3 9 9 4 0 9 , 1936.
a) & o ~ Y ~ U =
) D. G. Lahiotis, S. K. Park, and R. Krkhnaiah, “Optimal state vector estimation
for non-Gaussian initial state vector,” IEEE Trans. Automat. Contr., vol. AC-
16, pp. 197-198, Apr. 1971.
D. G. Lainiotis, “Partitioned linear estimation algorithms: Discrete case,” IEEE
Trans. Automat. Contr., vol. AC-20, pp. 255-257, Apr. 1975.
-, “Discrete Riccati equation solutions: Partitioned algorithms,” IEEE
(e.22) Trans. Automat. Contr., vol. AC-20, pp. 555-556, Aug. 1975.
-, “Discrete Riccati equation solutions: Generalized partitioned algorithms,”
Inform. Sci., vol. 15, no. 3, pp. 169-185, 1978.
For instance, for a one-dimensional problem, if the probability density M. Loeve, Probability Theory, vol. 11. New York: Springer-Verlag, 4th ed.,
function of the initial data is exponential, i.e., 1978.
[16] A. M. Makowski, “Results on the filtering problem for linear systems with non-
d p ( x o ) = e e o d ~ o , X>O. (e.23) Gaussian initial conditions,” in Proc. IEEE Conf. Decision Contr., Dec. 1982,
pp. 201-204.
Then the conditional characteristic function is given by [171 A. Makowski, W. S. Levine, and M. Ashev, “The nonlinear MMSE fiter for
partially observed systems driven by non-Gaussian white noise, with application to
failure estimation,” in Proc. 23rd IEEE Con5 Decision Contr., Las Vegas,
(e.24) NV, Dec. 1984, pp. 644-650.
[IS] V. S. Pugachev, “Recursive estimation of variables and parameters in stochastic
where systems described by difference equations,” Soviet Math. Dokl., vol. 19, no. 6,
pp. 1495-1497, 1978.
a) j?o=yo-XQ [I91 A. Sawitzki, “Finite dimensional filter systems in discrete time,” Slochastics,
vol. 5, pp. 107-114, 1981.
[201 E. T . Shapiro, “The random distribution method and its applications to the
and problem of nonlinear filtering in discrete time,” Radio Eng. Electron. Phys.,
vol. 26, pp. 48-54. June 1981.
b)$,=R(R+Q)-ly,: t2l (e.25)

and

a)po=Q

b) P , = R - R ( R + Q ) - ’ R for t 2 l (e.26) Cumulant Based Identification of Multichannel

Moving-Average Models
which is a finite-dimensional filtering realization for the problem.
CONCLUDING
REMARKS GEORGIOS B. GIANNAKIS, YUJIRO INOUYE, AND
JERRY M. MENDEL
Throughout the literature, a filtering structure is referred to as finite
dimensional if a finite set of sufficient statistics contains all the Abstmcl-Given cumulants of a stationary, perhaps noisy, non-
information necessary to determine the conditional probability law. Gaussian r-variate moving-average, MA(q) process, we study identifia-
Evidently, if a filtering system is such that the a posteriori density is in bility conditions, under which the MA coefficient matrices, the input
the same class as the initial density, then the filter is self-reproducing and statistics, and the order q can be uniquely determined. The selection of a
is finite dimensional if the initial density can be characterized by a finite unique representative from the equivalence class corresponding t o a given
set of statistics. In this context, Sawitzki [19] concluded that a necessary
and sufficient condition for the existence of a finite-dimensional filtering Manuscript received April 2 I , 1988; revised June 29, 1988. This work was supported
system is that the conditional distribution is of the “exponential class.” by the National Science Foundation under Grant ECS-8602531 and by NOSC Contract
N66001-85-D-0203.
(For analytic density functions, the “exponential class” describes a class G. B. Giannakis is with the Department of Electrical Engineering, University of
of densities that admit a set of sufficient statistics [lo].) The example Virginia, Charlottesville, VA 22901.
illustrates this idea since the a priori and a posteriori distribution belong Y. Inouye is with the Department of Control Engineering, Osaka University, Osaka,
to the exponential class. Theorem 2 states that for the linear filtering Japan.
J. M. Mendel is with the Department of Electrical Engineering-Systems, University of
problem, if the initial distribution can be characterized by a finite set of Southern California Signal and Image Processing Institute, Los Angeles, CA 90089-
statistics, the conditional distribution is also finite dimensional in this 0781.
context; however, not necessarily within the same class. IEEE Log Number 8927778.

0018-9286/89/0700-0783$01.oO 0 1989 IEEE

7 84 IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 34, NO. 7, JULY 1989

cumulant structure involves less restrictions than that corresponding to a Substituting (1) into (2) and using ASl-AS4, we obtain
given covariance structure. We derive two algorithms for estimating the r 7

(possibly) nonminimum-phase MA coefficient matrices. 4

C,(ml, m 2 ) = x H ( i + m l )
,=o
I. INTRODUCTION
Osm,, m l s q , (4a)
One of the basic problems in multivariate (or multichannel) time-series
model fitting is that of identifying parametric ARMA models given output =O m l > q andlor m2>q (4b)
statistics. Parametric modeling of multichannel processes is important in
where h,(i) is the ( m , I ) entry of H ( i ) . In deriving (4a) and (4b), we
applications such as multivariate linear prediction and spectral estimation
have also used the fact that third-order cumulants of the Gaussian process
[ 111, [ 181, multivariate control, econometrics [lo], array, and geophysi-
u ( i ) vanish for k 2 3; hence,
cal data processing [17].
Although multivariate models bear a superficial resemblance to the E { U( i + m , ) u T( i ) u , ( i + m 2 ) }= 0 . (4c)
corresponding univariate ones, their structure is more complicated, and
gives rise to quite deep identifiability and parameter estimation problems The latter makes cumulant based identification insensitive to additive
(see [l], 191, [15], and references therein). Up to now results are mostly (even colored) Gaussian noise. Similarly, for fourth-order cumulants we
restricted to identifying, under certain conditions, a unique model have
corresponding to a given covariance structure. Especially for MA models,
existing multichannel parameter estimation approaches are limited to
Gaussian processes (because they usually involve the maximum likelihood 1

method), and/or assume that the determinant of the MA transfer function

matrix has all of its zeros inside the unit circle [9], [IO], [13].
The objective of this note is to obviate the "multiminimum phase" and
Gaussianity assumptions, and derive theoretical identifiability and O s m , , m2, m l s q , (5a)
parameter estimation results of non-Gaussian multivariate MA processes =O m l > q andlor m 2 > q andlor m , > q (5b)
using third- or fourth-order statistics called cumulants. In [16], third-
order cumulants are employed for estimating multichannel AR models of where we have replaced AS2 by its fourth-order counterpart.
known orders. In the univariate case cumulants have received interest for From (4b) and (5b) we observe that determination of the MA order q
ARMA parameter estimation and other signal processing tasks (see [2]- amounts to checking the lag points for which C , ( m l , m2) or C m n ( m l ,
[8], 1121, [14], [19], and references therein). m2, m 3 ) become zero. The same order determination approach is
discussed in detail in 171 for univariate non-Gaussian MA processes.
II. IDENTIFIABILITY-PARAMETER
ESTIMATIONALGORITHMS Before addressing the identification of { H ( l ~ ) } i =and~ , { rm}k= I,
from third-order cumulants, we observe that y ( i ) in (1) can also come
Let y ( n ) be a stationary r-variate, non-Gaussian, MA(q) vector fromf?(k) = H(k)U, $(k) = U-lw(k),k = 1, ..., q, where U i s a n
process given by (s x s) arbitrary nonsingular matrix. Also, with w ( k ) = U $ ( k ) in AS2,
we find r m= U[C;=I u,,fl] U T ,where U,/ is the (m, I ) entry of U,and
4 f, is as in AS2 with w = $. To remove this ambiguity we need to fix one
y(n)= H ( k ) w ( n - k ) + u(n) (1) MA coefficient matrix, e.g., H ( 0 ) . As long as the set {r,};, is allowed
k=O
to be of general form we may assume, without loss of generality, that
r
and satisfying the following assumptions.
ASI: w ( i ) is an (s x 1) zero-mean, stationary, and non-Gaussian
process with components w,,,(i).
A S 2 : E { w ( i ) w r ( j ) w , ( k ) } = rm6,,P,
with {I'm};=l denotings x s With this assumption we establish the following identifiability theorem.
nonsingular input cumulant matrices. Notice that AS2 refers only to third- Theorem I: If ASl-AS6 hold true and the order q is known, then
order temporal independence, and that it does not include any spatial { H ( k ) }I=,, and { I', }; =,
can be uniquely identified from third-order
independence constraint. output cumulants.
AS3: u ( i ) is an r-variate, Gaussian, perhaps colored, zero-mean Proof: When H ( 0 ) = ZrX,, considering (4a) with ml = q, m2 = k ,
process with components u,,,(i). and m, = q, m2 = 0, we obtain
AS4: u ( i ) is independent of w ( i ) .
AS? H ( 0 ) and H ( q ) are r x s matrices of full rank s ( r 2 s).
To determine the MA coefficient matrices { H ( k ) } ; = o , the input
cumulants and the order q, we use third-order output
cumulants, which for the stationary r-variant process y ( i ) are defined at
the r x r matrices

(6b)
where y , ( i ) denotes the mth component of the vector process y ( i ) . A
similar definition to (2) is also used in [16] for third-order cumulants of From (6a) and (6b) we observe that the last r - s columns of C,(q, k)
vector AR processes. For notational simplicity we confine ourselves to and C,(q, 0) are zero. If we let C:(q, k) and C:(q, 0) denote the firsts
third-order cumulants although our results are expandable to fourth- and columns of C,(q, k ) and C,(q, 0 ) , respectively, and substitute (6b) into
higher order cumulant statistics (useful when r min AS2 is 0). In the (6a), we find
multichannel case, fourth-order cumulants are defined as s
Cz(q, k ) = x CF(q, O)hmi(k), m=l, ..., r. (W
Cmn(mi,
m2, md E E(~(i+m~)~~(i)ym(i+m,)y,(j+m,)J I= I

- E { Y ( i + m, )Y 7 ( i ) 1 E i y m ( i +m d y A i + m d ) If c1(1)(q,k ) denotes the ith column of C,*(q,k), then (6c) for m = 1,

- E { y ( i + m d y , ( i + ~ z ) I E {7Y( W n ( i + m 3 ) ) . . ., r yields
- E { Y ( i + ml ) y n ( i +4 1E { Y ' ( i ) y m ( i +m2)1. (3) Lc')(q,O)H'(k)=L'"(q, k) (7a)
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 34, NO. 7, JULY 1989 785

where the matrices L(;?,(q, 0) and L(;Lr(q,k) are given by

L'"(q, 0) = [c',"(q,0) *. . cj"(q, O)],

L'"(q, k ) [c',"(q, k ) . . . c;"(q, k ) ] . (7b)

Using the least-squares solution of (7a), the MA coefficient matrices

{ H ( k)} 1= I are given by the closed-form expression

H ( k ) = L ' ( ' ) ( q , k ) L ( ' ) ( g 0, ) [ L r ( " ( q O)L")(g,

, O)] (8)
where L ( ' ) ( q ,0) is of full rank because { c { ' ) ( q , O ) } ; = , are linearly
independent. To show the latter, observe from (6b) that cjf)(q, 0) =
H(q)y,('),where H ( q ) is of full rank (cf. A S ) , and y j r )stands for the ith
column of rr.But from the definition of rl(see AS2) we have that yj') =
7:') . Therefore, the columns {yjf)};=,,or equivalently the columns
{y!l)];=l are linearly independent because l', is nonsingular. Thus, the
,
columns {c,C')(q,0));- are linearly independent, and L ( ' ) ( q 0) , is of full
rank.
Notice that although nonsingularity of r, suffices for the uniqueness of
the solution in (8), if additional input cumulant matrices are nonsingular,
we may incorporate them in (8), after concatenating equations similar to
(7a) for different i's. When H(q) is known (cf. (8) with k = q), the input
cumulant matrices can be obtained from (6b) as

which in the square case (r = s) reduces to

where t denotes the generalized inverse.
r , = H - l ( q ) C m ( q ,0) m=l, ..., s. (9b) Notice that if Cm(O,q ) is nonsingular for some m = 1, 2, . , r , then
(13c) holds true if we use Cm(k,q ) instead of L ( k , 4). The least-squares
This completes the proof of the identifiability result and provides us with a solution in (13c) is employed because L(0, q ) is of rank s. The latter can
cumulant based parameter estimation algorithm for multichannel MA be shown if we substitute ( l l b ) into (13a) and rewrite L(0, q ) as L(0, q )
models. 0 = H(O)M(q)diag [HT(0)],where M ( q ) = [diag [ h t l ( q ) . . . h,,(q)l
If the true full rank matrix H(0) is not Zrx,, but known, we first operate .. . diag [ h , , ( q ) . . h,,(q)]],and diag [HT(0)]is a block diagonal rs x
on its rows to transform it to [3]= SH(O), where H, is an (s x s) r2 matrix with diagonal elements HT(0).To show that M ( q ) is of full
nonsingular matrix and S is an appropriate r x r elementary matrix. rank, let m,(q)denote the ith row of M ( q ) ,and { X I } : = , be a nontrivial
Then, using (8) with C,(q, k) E E { g ( i + q)yT(i)y,(i + k ) } , with set of constants. If
y ( i ) = S y ( i ) , we find the normalized MA model (with A(0)= Z,,, and
input cumulant matrix r,). Finally, we compare the true MA coefficient
matrices using H ( k ) = S-IH(k)H,, and the true input statistics using r m
= H;l[C;=, pl(Hs-l)ml]Hs-r = H ; l r , H ; r , where ( H ~ ] ) ~ ,the
i s (m,
I ) entry of HS-' (cf. discussion preceding Theorem 1 with U = H,-I). If
H(0) is not of full rank, H, will be singular, and although the H ( k ) and ( m i ( q ) } ; =are
, linearly dependent, then the right-hand side of (13d)
matrices can be identified, only part of the input cumulant matrices r mcan must be zero; i.e., at least one column of H ( q ) has to be the zero vector,
be determined via pseudoinverses because they involve H,-'. something that contradicts the full rank property of H ( q ) . Therefore,
If H ( q ) is not of full rank, the multivariate MA process lacks either M ( q ) and M ( q ) diag [HT(0)] are of full rank s. Because H ( 0 ) is of full
controllability or observability [101 and cannot be "block-identifiable" rank, we infer that L(0, q ) is also of rank s.
[I]. In (8) and (9a), if H ( q ) , and hence L(')(q,0) are not of full rank, To determine H(0) we consider (13c) for k = q and substitute H ( q )
pseudoinverses can be employed to identify the observable and controlla- into (1 IC) to obtain
ble part of the multichannel MA model.
As an alternative to the coefficient constraint H(0) = I,,, (or known), H(0)diag [h,,(O) ... h,(O)]H'(O)=E,, m=l, ..., r (14a)
we can impose a constraint on the input cumulant matrices { r,} =; , , where
namely we may substitute AS2 with

AS2': E{w,(l)w,(m)wk(n)}=l i=j=k, I=m=n,

Em L ( 0 , q)L'(O, 4 ) [ L ( 4 ,q)LT(O,4)1'Cm(q, 0). (14b)

=O otherwise Without any further assumptions (14a) is not guaranteed to have a

unique solution. To study the equivalence class of MA(q) models that
Notice that AS2 ' is equivalent to AS2, if we require r mto have all-zero satisfy (lo), and consequently (13c) and (14a), we use the following
entries, except the ( m , m ) entry which must be unity. Moreover, AS2' lemma.
refers to both spatial and temporal third-order independence, and it is Lemma I: An (s x s) matrix P = (p,,) is a permutation matrix, i.e..
meaningful when we have access to the input process (e.g., when we by definition
design inputs for stochastic identification). Using AS2 ' we can rewrite
(4a) as p,,= 1 ifj=j,,
=O ifj#j,, i , j = l , 2 , ' " , s

if and only if it satisfies

O s m 2 , m , i q (IO) Pdiag [ p , ...

~ p,]Pr=E,, i=l, ... , s
where diag [al . * . a,] stands for an (s X s) diagonal matrix with diagonal where E, = diag [ e l l .. . es], with e,, = 1 and ek, = 0 for k # i.
elements {a,} =; I . Proof: See Appendix A.
786 IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 34, NO. I, JULY 1989

Using Lemma 1 ye may show that if H ( 0 ) satisfies (14a), then any after estimating the AR matrices as in [16], and applying the results of this
matrix of the form H(0) = H(0)P also satisfies (14a), provided that P is note to the residual MA vector process. Detailed derivation of the
a permutation matrix. But because the nonuniqueness of H ( k ) in (13c) is multivariate ARMA case, a unified Kronecker-product formulation of the
due to the nonuniqueness in determining H ( 0 ) from (14a), we conclude ARMA identification approach, and simulation examples will be reported
that the following theorem holds true. in the near future.
Theorem 2: If ASI, AS2', and AS3-AS5 hold true, and { H ( k ) } f = ,
satisfy (lo), then { H ( k ) = H(k)P}f=,also satisfies (IO), where P is a APPENDIXA
permutation matrix. Conversely, if { H ( k ) } f = , ,and { H ( k ) } f = satisfy
,
(lo), then there exists a permutation matrix P such that { H ( k ) = Proof of Lemma I: For the only if part let T,, be the s x s matrix
H ( k ) P }f =O' obtained by interchanging the i th and j th columns of the identity matrix
I,,,. If P is a permutation matrix, then by definition it has a single
Proof: See Appendix B.
nonzero entry (which equals one) in each row and column, i.e.,
Note that, when second-order statistics are prescribed with input
covariance matrix R , = I,,,, the corresponding equivalence class is p,,= 1 if j = j , ,
{ H ( k ) = H ( k ) U } f = , ,where U is nonsingular and satisfies UUT =
I,,,. Moreover, autocorrelation based identification is largely restricted =O ifjtj,, i, j = l , 2, ..., s.
to multivariate Gaussian processes, and requires the zeros of det [ H ( w ) ]
to lie inside the unit circle [9], [IO], [13]. Therefore,
Although H(0) is in generalnot uniquely identifiablefrom (14a), in
the special case of triangular (s x s)H(O) matrices, unique solution of
diag[p,, ... p n ] = E J 8 , i = l , 2 , "',s (AI)
(14a) is possible. If a unique H(0) can be found from (14a), then and
{ H ( k ) }f = I can be obtained from ( 1 3 ~ ) Note
. that in the square case L (0,
q ) is of full rank, and that (13c) and (14b) can be simplified. In particular, P =El T I , ,+E2TzJ,+ ' . . +ESTsJs. (A2)
the matrix B, in (14b) will not contain H ( 0 ) . Let us consider the lower
triangular case, and write (14a) elementwise as Because the matrix E,T,,/ has all zero components except the (I, jl)th
component, which equals unity, we obtain
5 mm ( m J , J )

bm(i,j ) = hmk(O)h,k(O)hjk(O)= hmk(o)htk(o)hJk(o) (I5) E,T,,,E,,=E,T,,, if I = ;

k=l k=l
=O if I t i
where, for establishing the second equality in (15) we used the fact that
H ( 0 ) has nonzero entries hmk(0),h,k(O), /I,~(O) for m 2 k, i 2 k, and; and
2 k, respectively.
Notice that for m = 1, (15) yields b l ( l , 1) = h:,(O), b l ( l , 2) = El Tb,(E,,, TIJ,)
'= E, if I = i
h;,(O)h2,(0), ..., bl(1, s) = h~,(O)h,,(O).Because h,,(O) must be =O iflti.
nonzero, the entries of the first column of H(0) can be found using h l l(0)
= by3(1, I), h21(0)= b , ( i , 2)/b;2/3(1, I), ..., h,,(o) = b l ( i , SI/ It follows from (AI)-(A3) that
bC2I3(l, 1). The second column of H ( 0 ) can be determined if we
substitute m = 2 in (15) to obtain b2(2, 2) = h:,(O) + h:,(O), b2(2, 3) P diag I P , ~ . . . PJ = E ,T,,,
= h&(O)h31(0) + h:,(O)hn(o), " ' , b2(2, s) h;,(O)h~i(O) + which after using (A2) and (A4) yields
hi2(0)hs2(O).Similarly, the nth-column of H(0) can be computed if we
use m = n in (15), and observe that the equation corresponding to b,(n, T ,I , ,+ . . . +EsTSJJr=E , .
P diag [p,l . . . p s ] PT = E , T z J , ( E
I ) can be solved for h,,,(O), I = n , n + 1, . . * ,s. Note that the rest of the
quantities involved in the equation corresponding to b,(n, I ) are already For the if part, let P = (p,,) be an (s X s) matrix that satisfies
available from the equations corresponding to b , ( i , j ) with i In a n d j 5
1. An analogous procedure can be followed to identify upper triangular
H(0) matrices. The case of diagonal H(0) matrices is straightforward P d i a g [ p , , . . . p n ]P r = E i , ;=I, ..., s. (-45)
since from (15) we have that b,(m, m ) = hi,(O). Given { B , } ; = , our But from (A5) we have that
algorithm uniquely determines H(0) if it is triangular. The uniqueness in
the triangular case is also verified if we view (14a) as the L-U
decomposition of B, . p diag IPtl '..PlSlPT=~sxs 646)
Similar equations to (8) and (13c) can be derived for fourth-order ,=I
cumulants. From (5a), we infer that using fourth-order cumulants one can
first compute the products { h,,, , ( q ) h , , 2 ( k } f l , , 2 =from
I , which the which implies that P is nonsingular. Considering (A5) for i = 1, we
h,/(k) coefficients can be subsequently obtained. obtain
In practice, the identification algorithms described in this note require
substitution of the theoretical cumulant matrices by their estimates
P [ P ; , ...p;,]T=[l 0 ... 01T . (A7)
obtained through sample averaging. In this case, strong consistency of the Also, for i = 2, . . ., s (A5) yields
MA estimates is implied by the strong consistency of the sampled
cumulant estimates [8], [14].
PII[P21 . . . PSI1
111. CONCLUSIONS (A@
P1rIP2r . . . PSI
Cumulants of multichannel processes can be defined and used to study
identifiability conditions for non-Gaussian and nonminimum-phase MA From (A7) we infer that there exists a nonzero component in the first row
processes. Closed-form expressions relate the MA coefficient matrices of P. Moreover, because P is nonsingular, (A8) implies that if plk # 0,
with cumulants, and can be used, under certain assumptions, for then p,k = 0 for j = 2, . . ., s. Assume that there are two nonzero
parameter estimation. As in the univariate case, we have demonstrated components in the first row, say p l k # 0 and p I l # 0. In this case, pJk =
that cumulants play the role of cross-correlation in identifying nonmini- 0 and p,l = 0 for j = 2, * . * , s, something that contradicts the
mum-phase multichannel systems of finite impulse response. Cumulant nonsingularity of P. Thus, in the first row (and similarly in all rows) mere
based estimation of vector ARMA processes (useful for identification of is a single nonzero component, which must be unity because of (A7).
infinite impulse response multichannel stochastic systems) is possible, Therefore, P is a permutation matrix. 0
IEEE TRANSACTIONS ON AUTOMATIC CONTROL, VOL. 34. NO. 7. JULY 1989 787
APPENDIXB REFERENCES

Proof of Theorem 2: From (13c) we infer that the nonuniqueness of H. Akaike, “Canonical correlations analysis of time series and the use of an
information criterion.’’ Advances and Case Studies in System Identification, R.
H ( k ) is due to the nonuniqueness in determining H ( 0 ) from (14a). Mehrd and D. Lainiotis. Eds. New York: Academic, 1976.
Hence, it suffices to prove Theorem 2 for H ( 0 ) . For simplicity, consider r D. R. Brillinger and M. Rosenblatt, ‘Computation and interpretation of kth-order
= s. spectra.” Spectral Analysis of Time Series, B. Haria, Ed. New York. Wiley.
1967, pp. 189-232.
For the direct statement of the theorem we wish to show that if H ( 0 ) G. B. Giannakis. “Cumulants: A powerful tool in signal proceshg.” Proc.
satisfies (14a), then H ( 0 ) = H(O)P also satisfies (14a)provided that P i s IFFF y 1717-1714. 1987.
a permutation matrix. Indeed, if H ( 0 ) = H(O)P, and B , as in (14b) with G . R Giannahl? .tiid I \1 \ l ~ ~ i d ~“Identification
~l. n~mniiiiiliut:~ i’h.isc
A(0)used instead of H(O), then it is not hard to show that E,, = B,. system\ using higher-order statistics.‘ / & € E Tran.5. A c o u ~ t . ,Speech. I,gnal
Moreover, Processing. vol. 37. pp. 360-377. 198’). .
ci. B. Giaiiirakis d i d A. SHIIIIII.“New results on statc-apace and input-output
identification of non-Gaussian processes using cumulants,” in Proc. Soc. of
H(0) diag [fi,,,,(O) ... h;,,(O)]f?’(O) Phoro-Opt. Inslrumenl. Eng. (SPIE’87), San Diego, CA, Aug. 1987, vol. 826.
pp. 19-205.
Y. Inouye, G. Giannakia, and J. Mendel. “Cumulant based parameter estimation
of multichannel moving-average processes.” in Proc. In/. Con$ .4cousr.
Speech, Signal Processing, Apr. 1988, pp. I?i-1255.
G. B. Giannakis and J . M. Mendel. “ARMA order determination via higher-order
statistics.” presented at the Int. Conf. Math Theory of Networks and Sqst.,
Phoenix, AZ. June 1987.
G. B. Giannakis and A Swami. “On estimating non-causal non-minimum phase
ARMA models of non-Gaussian processes,” in Proc. 4rh A S S P Workshop on
Spectrum Estimalion and Modeling, 1988; also submitted for publication.
E. J. Hannan, Mulriple Time Series. New York: Wiley, 1970.
~ ~~, ”The identification and parameterisation of ARMAX and state space
forms.” Economerrica, pp. 713-723. 1976.
=H(O)diag [ h , , ( O ) ... h,,(O)lH‘(O) (B1) S. L. Marple Jr. and A. H. Nutall. “Experimental comparison of three
multichannel linear prediction spectral estimators.” Proc. /€E€, pp. 2 18-229.
where for establishing the third equality in (Bt) we have used Lemma I . 1983.
C . L. Nikias and M. R. Raghuveer. “Bispectrum estimation: A digital signal
To prove the cznverse statement, let us assume that (14a) is satisfied by processing framework,”Proc. I€€E, pp. 869-891, 1987.
both H ( 0 ) and H(O), i.e., M. S . Phadke, and G. Kedem. “Computation of the exact likelihood function of
multivariate MA models.” Biometrica, pp. 51 1-520, 1978.
H(O)diag [hml(O) . . . hms(O)lH’(O)-Bm B. Porat and B. Friedlander. “Performance analyais of parameter estimation
algorithms based on high-order moments.“ in Proc. Int. ConJ Acoust., Speech,
= H ( O ) diag [&,,,,(O). . . fims(0)]Hr(O)-&. (B2) Signal Processing, Apr. 1988, pp. 2412-2415, also submitted for publication.
M. Priestley. Spectral Analysis and Time Series. New York: Academic, I981
If we define P ) , B,,
[ H r ( 0 ) H ( 0 ) ] - i H 7 ( 0 ) f i ( Othen = B,, p(0)= M. R. Raghuveer, “Multichannel bispectrum estimation,” in Proc. 3rd A S S P
Workshop on Spectrum Estimation and Modeling, 1986, pp. 21-24.
H(O)P, and E. A . Robinson Multichannel Time Series Analysis with Digital Computer
Programs. San Francisco, CA: Hclden Day, 1967.
h;,m= h“p,, . (B3) 0 . N. Strand, “Multichannel complex maximum entropy (autoregressive) spectral
I= I analysis,’’ IEEE Trans Auromat. Contr., vol. AC-22, pp. 634-640, 1977.
J. K . Tugnait, “~dentificationof nonminimum phase linear stochastic systems,”
Automalica, vol. 22, pp. 457-464, 1986.
Substituting (B3) into (B2), we obtain

or, with E, defined as in Lemma 1,

Adaptive Control of a Class of Multivariable Nonlinear
Systems and Convergence Analysis

JINGXIN ZHANG AND SHIJUN LANG

But (B4) can also be written as
Abstract-This note is concerned with adaptive control of a class of
multivariable nonlinear systems which can be characterized by a stochas-
tic multivariable Hammerstein model whose linear part possesses an
arbitrary interactor matrix. A simple suboptimal control law is derived
which provides an efficient way to control a multivariable Hammerstein
where model whose linear part is not necessarily minimum phase. A direct
adaptation scheme is presented to implement the control law and the
A, = E , - P diag [ p , , . . . p,$]Pr. 036) global convergence of the algorithm is established.
Because H(0) is of full rank, (€35) implies that A , = 0, which after using
(B6) and Lemma 1 shows that P is a permutation matrix. c! I. INTRODUCTION

In a few specific situations, adaptive control alzorithms based on SISO

ACKNOWLEDGMENT
nonlinear models have recently been studied, e.g., Ill-[h], and some
The authors wish to thank A. Swami whose helpful comments analyses have appeared in the literature [3]-[61. cor MIMO nonlinear
improved the content of this technical note. Part of the work in this note
was presented at the International Conference on Acoustics, Speech, and Manuscript received March 1 I . 1988.
The authors are with the Department of Automatic Control, Northeast University of
Signal Processing ( I C A S S P ’ 8 8 ) , and was performed when the first two Technology, Shenyang, Limning, People’s Republic of China.
authors were at the University of Southern California. IEEE Log Number 8927771

0018-9286/891’0700-0787$01.OO O 1989 IEEE