Kapteyn, A.J.

Neudecker, H.

Wansbeek, T.J.

English

Research Papers in Economics

PSYCHOMETRIKA--VOL, 51, NO. 2, 269--275 
JUNE 1986 
AN  APPROACH  TO  n-MODE  COMPONENTS  ANALYSIS 
ARIE KAPTEYN 
TILBURG UNIVERSITY 
HEINZ NEUDECKER 
UNIVERSITY  OF AMSTERDAM 
TOM WANSBEEK 
GRONINGEN  UNIVERSITY 
As  an  extension  of  Lastovicka's four-mode  components  analysis  an  n-mode components 
analysis is developed. Using a convenient notation, both a canonical and a least squares solution 
are derived. The relation between both solutions and their computational aspects are discussed. 
Key words: principal components, multidimensional matrices. 
Introduction 
In  a  recent  article,  Lastovicka (1981) extends Tucker's  three-mode factor analysis 
model to four modes and  presents  a  "canonical" solution. In this paper we extend his 
research in three respects: 
(i)  We extend his results to an arbitrary number of modes; 
(ii) We  introduce some new  notation, fit for handling multidimensional matrices, 
and show how the calculus involved works; 
(iii) We derive a least squares (LS) solution and compare it to Lastovicka's solution. 
The problem and the notation are introduced first. We then present the LS solution and 
compare it with Lastovicka's solution, which we will denote by the term "canonical." An 
iteration method to compute LS is given next.  We compare the canonical and the LS 
solution by means of an empirical example. 
Apart from Lastovicka's work cited above, various other related references  can be 
mentioned. Carroll and Chang (1970) present  methods for analyzing multidimensional 
matrices. This work is restrictive in the sense that the component dimensionality of each 
mode is taken equal. Carroll, Pruzansky and Kruskal (1980) extend these results to in- 
clude linear constraints on the parameters. The present paper can be considered to be an 
n-mode extension of the work of Kroonenberg and de Leeuw (1980) on three-mode data, 
which itself is  a  further development of Tucker's (1966) work on three-mode principal 
components. They also give an algorithm (TUCKALS3) based on what they call alter- 
nating least squares.  Hence, Kroonenberg (1983), in his comprehensive monograph on 
Requests for reprints should be sent to Tom Wansbeek, Econometrics Institute, Groningen University, PO 
Box 800, 9700 AV Groningen, THE NETHERLANDS. 
The first draft was  written while Wansbeek was with the Netherlands Central  Bureau of Statistics. We 
thank Jaap Verhees for performing the computations arid for many discussions on the subject, John Lastovicka 
for kindly making available his data to us, and the Editor, the referees,  Jeroen Weesie and Pieter Kroonenberg 
for their useful comments. 
0033-3123/86/0600-5022500.75/0 
©  1986 The Psychometric Society 
269 270  PSYCHOMETRIKA 
three-mode data, describes our model as a "Tuckern model", and our iteration method as 
a "TUCKALSn algorithm." 
We agree with the observation by Kroonenberg (1983): "In principle, the extension 
to n-mode data is straightforward, but it becomes increasingly complex to keep track of 
the  summations  over the  proper indices.  Moreover, the description of such  an  n-mode 
procedure  becomes exceedingly cumbersome without  new  notation"  (p.  73).  These are 
precisely the points where the present paper might offer a useful contribution. 
The Problem and Some Notation 
Let A i be an unknown :i x  mi matrix, :i >  mi, satisfying A~ A i =  Iml, i =  1  ..... n, and 
define 
A  =  A I  ®--" ®A..  (i) 
The symbol "®"  refers to  the (right) Kronecker product, with, for example, C ® D  = 
[cij D]. We do not introduce special notation for the multiple Kronecker product like the 
one appearing in (1). Let 
n 
t =  1-It,  (2) 
n 
m  =-  ]-Imi .  (3) 
i=1 
So,  t  _> m  and  A  is  of order t  × m.  Let y  be  a  (known) t-vector,  whose  elements  are 
identified by n indices, with the i-th index (i =  1  ..... n) assuming values from 1 to t~. The 
elements are arranged in such a  way that  the first index runs  slowly and  the last index 
runs fast. Analogously, let b be an (unknown) m-vector, also with elements identified by n 
indices, with the i-th index running from 1 to m~ (i =  1  ..... n). The elements are arranged 
in the same way as the elements of y. (The vectors y and b can be considered as "stacked" 
versions of n-dimensional matrices.) It is our aim to choose At, i =  1, ..., n and b in such 
a way that Ab represents y "as well as possible." 
Let Ci be the t  ×  t  commutation matrix that changes the running order of the indi- 
ces in such a  way that when applied to y the i-th index runs fastest and the running order 
of the other indices remains unaffected. Let further D~ be the m  ×  m commutation matrix 
that reshuffles the elements of b  in the same way as does Ci with the elements of y. The 
sizes of Cj and D~ differ. As C~ and D~ are commutation matrices they have the property 
C~ C~ =  C~ C, =  Ie,  (4) 
Di D~  =  D~D  i =  Ira.  (5) 
The C:s  and  D[s are generalizations of the well-known commutation matrix for n =  2, 
which has been extensively studied by, among many others, Balestra (1976), Magnus and 
Neudecker (1979), and Henderson and Searle (1981). 
It  will  be  convenient  to  introduce  the  matrix  A i -  A 1 ®--- ®Ai_ 1 ®Ai+ 1 ®  --. 
®  A,. The following definitional relationship then emerges: 
Ai® A i =  C~(A t ®'"  ®  A,)D[ =  Ci AD'i.  (6) 
For reasons to become clear later on we define the matrices Yi and Bi by 
vec Y'i =  C~ y  and  vec B'~ =  D~ b,  (7) 
where Yi is of order (:  i × :i) and Bi is of order (m  i x  m3, with : ~ =- :/:i  and m ~ =- m/mi. ARIE KAPTEYN,  HEINZ  NEUDECKER  AND  TOM  WANSBEEK  271 
It will be clear that in spite of the normalization A~ A~ =  Ira,, the representation (A 1 
®--- ®  A~)b is not  unique, because if T~, i =  1  .....  n,  are  orthogonal m~, rn~ matrices 
there holds 
(ArT x ®...  ®  AnT,,)(T  1, ®...  ® T',,)b =  (A t ®...  ®  A,,)b.  (8) 
This indeterminacy will be used to pick a convenient solution when discussing the repre- 
sentation problem. 
The Least Squares Solution 
The LS solution amounts to finding At ..... An and b such that 
Q -  (y -  Ab)'(y -  Ab)  (9) 
is minimal, subject to A'~ A~ -- Imt and (1). The solution to this problem is readily obtained. 
First, it follows from standard  LS theory that for any choice of At .....  A~, the solution 
for b is 
Notice that 
1~ =  A'y =  (A~ ®---  ®  A,,)'y.  (lO) 
y'(At ®  .. ® A,,)I~ =  y'(AxA'  1 ®...  ®  A,, A',,)y 
=  y'C~  C~(A tA'~ ®.. "An A'n)C~ C~ y 
=  (vec Y'i)'(AiA i' ~  A i A~)vec Y~ 
t  i  i,  =  tr AiYiAA  Y~AI,  (11) 
for any one index i. It is convenient to define 
1  i  if  Si -  Yi A A  Yi-  (12) 
Then 
Q =  y'y -  tr A~S~A~.  (13) 
Minimization of Q  subject to A~ Ai =  I,,~, i =  1  .....  n, can be done by differentiating the 
Lagrangean function: 
n 
L =  y'y -  tr A'ISiAi +  tr ~  Fi(A~Ai- Ira,),  (14) 
i=1 
with F~ a symmetric ms  x  m~ matrix of Lagrange multipliers, with respect to A~ and setting 
the result equal to zero. This yields: 
SiAi-  AiFi =  0.  (15) 
So a  solution is to choose for F i the diagonal matrix Fi containing the m~ largest eigen- 
values of Si, and for A~ the corresponding orthonormal eigcnvectors. 
SinceSi(-  '  ~  ~'  Y~ A A  Y~) depends on the unknown parameters in A ~ the solution to the 
minimization problem has to be obtained iterativcly. An iteration process and its conver- 
gence properties are discussed below. 
The fit of the LS solution is assessed as follows. (Carets indicate LS solutions.) Define 
the f-vector of residuals: 
~.._~ y -  ./~fi,  (16) 
then, analogous to (13) we have: 272  PSYCHOMETRIKA 
d'd  =  y'y  -  ~'~  =  y'y  -  tr A'~S,.~, =  y'y  -  tr F,,  (17) 
for any one i. So the sum of the m~ largest eigenvalues in the LS solution is the same for 
all i. It is obvious to define 
R 2 =  fi'----~- tr ~  (18) 
--  y'y  -  y,y 
as the coefficient of determination. 
The Canonical Solution 
Lastovicka (1981)  proposes a canonical solution Z,~, i =  1  ..... n, where the columns 
of J~ are the rn~ orthonormal eigenvectors corresponding to the m~ largest eigenvalues of 
Y~ Y~.  The  solution for b  is  ~ =  (.~ ®..'  ® ~'n)y. The  LS  and  the canonical solutions 
differ only with respect to the A~. Given A~, b is the same in both cases. With LS, we have 
vec B', =  D,~ =  Di(,~ ®'"  ® A'n)y =  Di(.~.~ ®""  ® An)CiC, y 
=  (Ai'®  .~'~)Ciy  =  (,~l,®  .~'~)vec Y'i =  vec A~ V~,'.  (19) 
So 
i~i =  ~i'Yi -~i.  (20) 
By the same procedure the canonical solution is 
~i  =  ~i'Yi~ki.  (21) 
As Y'~ ~,~.~'Y~  .~ =  ~i.~  ~  (cf. (15)), ~  =  $~'Y~  i.~ contains the first m~ principal components 
of .i,~'Y~, provided  that  mi  <-m i.  This  (necessary) condition  follows from  the  fact  that 
•  ~'Y~ A,i  is  an  m  ~ x  m~ matrix.  (It  is  clear  that  this  condition cannot  be  satisfied when 
n =  2, unless ml =  m2  .) 
Let us write the canonical equations as: 
Y~ Y~'~i =  ~  ~i,  where  -~'i ~i =  Ira,"  (22) 
Now, Ai does not contain the eigenvectors of Y~.~.~'¥i, but of Y~ Y~.  Hence if for all i 
except one A~ =  le~, the two solutions coincide, as A  ~ =  Ie~ in that case. Thus, when con- 
sidering index i, the canonical solution does not involve the data reduction with respect to 
any other index. 
Still another interpretation of the LS solution is obtained by observing that A~A  ~' is 
idempotent and hence a  projection matrix. It projects onto the space spanned by A ~. So 
the LS solution amounts to a  search for principal components, not of Y~ itself, but of its 
projection onto the space spanned by A ~. 
An Iteration Method 
An obvious iterative procedure to obtain the LS solution is as follows. Take ~2 ..... 
.Z,, as starting values for A2 ..... A,. Use these values to form a  first estimate of A ~, .~ 
say. Compute the first estimate of A1, say -~1tl), as the m~ eigenvectors of Y~,~).~  Y1 
corresponding to the ml largest eigenvalues. Using -~10) and ~a ..... .Z.  n, we can form ~-~1) 
and estimate -~2t1~ in a  similar manner.  Having computed A1~11 through ~,~), we start 
with A1 again and form estimates ~-~t2j  ..... Ant2). The process is continued until conver- 
gence. 
Convergence of the iterations follows from the following considerations. Q defined in 
(9) is quadratic and consequently nonnegative. For each i, the solution for A~ correspond- ARIE KAPTEYN, HEINZ NEUDECKER AND TOM WANSBEEK  273 
ing to (15) minimizes Q, for any value of A ~. So, in the above iterative process, each newly 
computed Ai lowers the values of Q or leaves it unaffected. Thus, we obtain a nonincreas- 
ing  sequence of values  of Q  which is  bounded from below by zero. As  a  result,  the  se- 
quence converges. 
Some Empirical Results 
In order to get some practical experience with the LS approach we applied it to the 
same  data  set  as  was  used  by  Lastovicka  (1981).  This  data  set  consists  of a  four- 
dimensional  matrix  of order  27 x  6  x  5 x  16,  containing  the  scores  (ranging  from  1: 
"strongly disagree" to 6: "strongly agree") given by fl =  27 individuals to each of ~z =  6 
TV commercials after each of ~a =  5 exposures for each item from a  list of E4 =  16. This 
data set was standardized such that the 16 x  16 matrix Y~, Y4 is a correlation matrix; see 
Lastovicka (1981) for details. (The factor 1/(27) t/2 mentioned there on page 51 should read 
1/(810)1/2.) 
We first recomputed Lastovicka's results, which are based on m x =  4, m 2 =  2, m 3 =  3 
and m4 =  3. In general, we were able to reproduce his results, with one notable exception 
(as  was  brought to our attention  by Jaap  Verhees, see also  Verhees,  1985,  the columns 
headed  "Ip." and  "IIp." in  his  Table  3,  containing the  varimax  rotated eigenvectors of 
Y3 Y3 corresponding to the largest and second largest eigenvalues, respectively, should be 
interchanged).  As  a  result,  most  of the  entries  in  his  Table  5,  which  gives  the  "core" 
matrix, a  two-dimensional display of b, are incorrect. (But we were also unable to repro- 
duce the entries that  are unaffected by the error, i.e.,  those in  the  two columns headed 
"Emotive Response" in his Table 5. The overall magnitudes of the entries in this table are 
too low to be correct. A rough check is provided by computing the sum of squares of the 
entries, i.e., b'b in our notation. This yields .069 or, cf. (18), R 2 =  .004, an improbably low 
value. In fact, ~'~ =  4.634 according to our computations.) 
We next analyzed the same data by the LS method. We used the same values for the 
re{s,  and followed the iteration method discussed above. The computer programming is 
straightforward, with the exception of the procedure to build up the matrices of the type 
Si.  In order to avoid storage capacity and time limits, one has  to break down the con- 
struction of Si into a set of nested "do loops" in a scalar fashion rather than to construct 
all intermediate matrices in full. 
As expected, convergence towards the optimum was rapid. The iteration was termi- 
nated as soon as tr F~, for some i, showed a  relative difference from tr F~ in the previous 
round of less than  5  x  10 -6.  This occurred with tr F 4  at the end of the tenth  round of 
iterations. (Notice that  tr Fi depends on all Ai.  In general  it will  therefore only remain 
constant between iterations if all A~ remain constant.) 
The most  interesting  question  is,  of course,  how  the  LS  method  compares  to  the 
canonical one with respect to the fit; see Table  1, where we show how y'y (= 16, as the 
16 x  16 matrix Y'4 Y4 is a correlation matrix, and y'y =  tr Y~Y~ for all i) is split up in an 
"unexplained" part (6'~) and an "explained" part (6'~), see (17). It appears that R 2 goes up 
from  .290  to  .323  when  replacing  the  canonical  method  by  ours.  The overall  level  of 
explanation may not seem to be very high, but one has to keep in mind that the job of 
reproducing  the  27 ×  6  x  5 x  16 =  t2,960  elements  of  y  has  to  be  borne  by  the 
4  x  2  x  3  x  3 =  72 elements of the "core matrix" t~. 
As  to  the  results  on  the  varimax-rotated  A~-matrices, there  were some  differences 
between  both methods.  The results  on  A2  (ads)  and  A 3  (exposures) did  not differ very 
much, but those on A t  (individuals) and  A4 (items) did. In these two cases, there was a 
marked difference in  that (a) the principal loadings per row were often in different col- 
umns (for 21  of the 27 rows of At, for 4  of the  16 rows of A4), and (b) the LS solution 274  PSYCHOMETRIKA 
TABLE  1. 
LS  versus  the  canonical  method:  some results 
Method  y'y  e 'e  b' b  R  2 
LS  16.000  10.835  5.165  .323 
Canonical  16.000  11.366  4.634  .290 
shows  a  tendency to concentrate the principal loadings (i.e., the largest factor scores) in 
the columns corresponding with the largest eigenvalues. For example, A1 has  11 of the 27 
row maxima (absolute values) in the first column and ~  only 7; for -&4 these figures are 9 
out of 16 and for .Z,  4 6 out of 16. 
As can  be expected from these differences, there were also considerable differences 
between G and ~. 
Conclusion 
In this paper we presented a least squares approach to multidimensional component 
analysis.  A  possibly  fruitful  byproduct  of the  paper  is  the  introduction  of a  notation, 
based  on  permutation  matrices,  that  allows  for easy handling  of the  algebra  involved, 
making derivations simple and transparent. 
As  we  have  LS as  our loss  function, our method scores, of course,  higher on  this 
criterion than the canonical one.  A  troublesome aspect of multidimensional component 
analysis appears  to be the apparent sensitivity of at least part  of the results to the cri- 
terion used. A possible way to study this phenomenon is to consider the individuals in the 
data set as a  random sample from some population and to set up a  statistical model to 
explain the observations. A model which would justify the LS procedure is the following 
one: 
y  =  Ab +  u,  (23) 
where u is an E-vector of i.i.d, errors with mean zero. If we assume the elements of u to be 
normally  distributed,  the  LS  procedure provides maximum  likelihood  estimators  of A 
and b. 
Given this model, we can carry out an asymptotic F-test to see whether Lastovicka's 
results are consistent with the data (e.g., Goldfeld &  Quandt,  1972). The only ingredients 
we need for this test are the  R2-values for the canonical and  the  LS solutions. We find 
that  F  (145,  12743)= 4.28,  which  is  highly  significant.  Thus,  assuming  the  statistical 
model, (23), the canonical solution is rejected by the data. 
References 
Balestra, P. (1976). La d#rivation matricielte [Matrix derivatives].  Paris: Sirey. 
Carroll, J. D., & Chang, J. J. (1970). Analysis of individual differences  in multidimensional scaling via an N-way 
generalization of "Eckart-Young"  decomposition.  Psychometrika,  35, 283-319. ARIE KAPTEYN,  HEINZ NEUDECKER  AND TOM WANSBEEK  275 
Carroll, J. D., Pruzansky, S., & Kruskal, J. B. (1980). CANDELINC: A general approach to multidimensional 
analysis  of many-way arrays with linear constraints on parameters. Psychometrika, 45, 3-23. 
Goldfeld, S. M., & Quandt, R. E. (1972). Nonlinear methods in econometrics. Amsterdam: North-Holland. 
Henderson, H. V., & Searle,  S. R. (1981). The Vec-permutation matrix, the Vec-operator and Kronecker prod- 
ucts: A review. Linear and Multilinear Algebra, 9, 271-288. 
K roonenberg, P. M. (1983). Three-mode principal component  analysis. Leyden: DSWO-Press. 
Kroonenberg,  P.  M., &  de  Leeuw,  J. (1980), Principal  component  analysis of three-mode data  by means of 
alternating least squares algorithms. Psychometrika, 45, 69-97. 
Lastovicka, J. L. (198 l). The extension of component analysis  to four-mode matrices.  Psychometrika,  46, 47-57. 
Magnus,  J.  R.,  &  Neudecker,  H.  (1979). The  commutation  matrix:  Some properties  and  applications.  The 
Annals of Statistics, 7, 381-394. 
Tucker, L. R. (1966). Some mathematical notes on three-mode factor analysis. Psychometrika, 31,279-311. 
Verhees, J, (1985). Four-mode principal component  analysis;  evaluation  of solutions from two alternative algorithms 
(Report).  Groningen: Groningen University, Econometrics Institute. 
Manuscript received 9/14/82 
Final version received 9/26/85 

Four-mode principal component analysis; evaluation of solutions from two alternative algorithms

La d#rivation matricielte [Matrix derivatives].

Nonlinear methods in econometrics.

Principal component analysis of three-mode data by means of alternating least squares algorithms.

Some mathematical notes on three-mode factor analysis.

The commutation matrix: Some properties and applications.

The Vec-permutation matrix, the Vec-operator and Kronecker products: A review. Linear and Multilinear Algebra,

Three-mode principal component analysis.

An approach to n-mode components analysis.

An approach to n-mode components analysis.

Abstract

Similar works

Full text

Available Versions

Research Papers in Economics