1. Introduction
A B S T R A C T
Principal components analysis (PCA) is a popular descriptive multivariate method for handling quantitative data and it can be extended to deal with qualitative data and mixed measurement level data. The existing algorithms for extended PCA are PRINCIPALS of Young et al. (1978) and PRINCALS of Gifi (1989) in which the alternating least squares algorithm is utilized. These algorithms based on the least squares estimation may require many iterations in their application to very large data sets and variable selection problems and may take a long time to converge. In this paper, we derive a new iterative algorithm for accelerating the convergence of PRINCIPALS and PRINCALS by using the vector s algorithm of Wynn (1962). The proposed acceleration algorithm speeds up the convergence of the sequence of the parameter estimates obtained from PRINCIPALS or PRINCALS. Numerical experiments illustrate the potential of the proposed acceleration algorithm. © 2010 Elsevier B.V. All rights reserved. Principal components analysis (PCA) is a popular descriptive multivariate method for handling quantitative data and is extended to deal with qualitative data and mixed measurement level data. The existing algorithms for extended PCA are PRINCIPALS of Young et al. (1978) and PRINCALS of Gifi (1989) in which the alternating least squares (ALS) algorithm is utilized. These algorithms alternate between optimal scaling for quantifying qualitative data and the analysis of the optimally scaled data using the ordinary PCA approach. We will refer to PRINCIPALS and PRINCALS as PCA.ALS when not distinguishing between them. In the application of extended PCA for very large data sets and variable selection problems, many iterations and much computation time may be required for convergence of PCA.ALS. For example, for PCA based on a subset of variables for qualitative data, the PRINCIPALS approach taken by Mori et al. (1997) obtained estimates only after a large number of iterations. In this paper, we derive a new iterative algorithm for speeding up the convergence of PCA.ALS. The proposed algorithm uses the vector ε (νε) algorithm of Wynn (1962) for accelerating the convergence of PCA.ALS. We refer to νε accelerated PCA.ALS as v ε -PCA.ALS. During iterations of v ε -PCA.ALS, the v ε algorithm generates an accelerated sequence of optimal scaling data estimated by PCA.ALS. Then the v ε accelerated sequence converges faster than the original sequence of the estimated optimal scaled data. The paper is organized as follows. We briefly describe PCA for a mixture of quantitative and qualitative data in Section 2, and describe PRINCIPALS and PRINCALS for finding least squares estimates of the model and optimal scaling parameters in M. Kuroda et al. / Computational Statistics and Data Analysis 55 (2011) 143–153 Section 3. Section 4 presents the procedure of v -PCA.ALS that adds the v ε algorithm to PCA.ALS for speeding up convergence. Numerical experiments in Section 5 illustrate the performance and properties of the v ε acceleration of PRINCIPALS. In Section 6, we present our concluding remarks.
2.  Principal components analysis with variables measured with a variety of scale levels
PCA transforms linearly an original data set of variables into a substantially smaller set of uncorrelated variables that contains much of the information in the original data set. The original data matrix is then replaced by an estimate constructed by forming the product of matrices of component scores and eigenvectors. Let X = (X 1 X 2  . . . X g ) be an n × p matrix of n observations on p variables and be columnwise standardized. In PCA, we postulate that X is approximated by the following bilinear form:
X = ZA T ,

where Z = ( Z 1 Z 2  . . .  Z r ) is an n × r matrix of n component scores on r ( 1 ≤ r  ≤ p ) components, and A = ( λ 1 A 2  . .  A r ) is a p × r matrix consisting of the eigenvectors of X T X /n and A T A = I,. Then we determine model parameters Z and A such that θ = tr(X − Î) T (X − Î) = tr(X − ZA T ) T (X − ZA T )
is minimized for the prescribed r components. Ordinary PCA assumes that all variables are measured with interval and ratio scales and can be applied only to quantitative data. When the observed data contain several different types of variables with nominal, ordinal, interval and ratio scales, ordinary PCA cannot be directly applied to such data. In such situations, optimal scaling is used to quantify the observed qualitative data and then ordinary PCA can be applied. To quantify X i of qualitative variable j with K i categories, the vector is coded by using an n × K i indicator matrix G i with entries g ( ) ik = 1 if object i belongs to category k, and g ( ik = 0 if object i belongs to some other category k ( ≠ k ), i = 1, . . , n and k = 1, . . . , Kj. Then the optimally scaled vector X 2 of X 2 is given by X 2 =  G 2 j, where α j is a K 2 × 1 score vector for categories of X i . Let X * = (X * X *  . . . X * ) be an n × p matrix of optimally scaled observations to satisfy restrictions
where 1 n and O n are vectors of ones and zeros of length n and p respectively. In the presence of nominal and/or ordinal variables, the optimization criterion (2) is replaced by
To apply PCA to data with mixed measurement levels, we determine the optimal scaling parameter X*, in addition to estimating Z and A.
3.  Alternating least squares algorithm for principal components analysis
A possible computational algorithm for estimating simultaneously Z, A and X* is the ALS algorithm. The algorithm involves dividing an entire set of parameters of a model into the model parameters and the optimal scaling parameters, and finds the least squares estimates for these parameters. The model parameters are used to compute the predictive values of the model. The optimal scaling parameters are obtained by solving the least squares regression problem for the predictive values. Krijnen (2006) gave sufficient conditions for convergence of the ALS algorithm and discussed convergence properties in its application to several statistical models. Kiers (2002) described setting up the ALS and iterative majorization algorithms for solving various matrix optimization problems.
3.1. PRINCIPALS
PRINCIPALS proposed by Young et al. (1978) is a method for utilizing the ALS algorithm for PCA of data with mixed measurement levels of single discrete and continuous, single nominal, ordinal and numerical variables. PRINCIPALS alternates between ordinary PCA and optimal scaling, and minimizes θ* defined by Eq. (4) under the restriction (3). Then θ* is to be determined by model parameters Z and A and optimal scaling parameter X*, by updating each of the parameters in turn, keeping the others fixed. For the initialization of PRINCIPALS, we determine initial data X*(0) . The observed data X may be used as X*(0) after it is standardized to satisfy the restriction (3). For given initial data X*(0) with the restriction (3), PRINCIPALS iterates the following two steps: M. Kuroda et al. / Computational Statistics and Data Analysis 55 (2011) 143–153 ● Model parameter estimation step: Obtain A (t) by solving A = ADr .
where A T A = I r and D r is an r × r diagonal matrix of eigenvalues, and the superscript ( t ) indicates the t-th iteration. Compute Z (t) from Z (t) = X*(t)A (t) ● Optimal scaling step: Calculate î (t+1) = Z (t)A (t)T from Eq. (1). Find X*(t+1) such that = arg min tr(X* -  ¥(t+1) T (X* -  ¥(t+1) ) for fixed X ( '+') under measurement restrictions on each of the variables. Scale X *( '+') by columnwise centering and normalizing.
3.2. PRINCALS
PRINCALS by Gifi (1989) can handle multiple nominal variables in addition to the single nominal, ordinal and numerical variables accepted in PRINCIPALS. We denote the set of multiple variables by f M and the set of single variables with single nominal and ordinal scales and numerical measurements by 𝑠𝑠𝑗𝑗 . For X consisting of a mixture of multiple and single variables, the algorithm alternates between estimation of Z, A and X* subject to minimizing θ* = tr(Z − X*A) T (Z − X*A) under the restriction Z T 1 n = 0 r   and   Z T Z = n I p . (6) For the initialization of PRINCALS, we determine initial data Z ( 0 ) , A ( 0 ) and X *( 0 ) . The values of Z ( 0 ) are initialized with random numbers under the restriction (6). For j ∈ $ M , the initial value of Xj f is obtained by Xj + (0)  = Gj (Gj Gj) − 1Gj Z(0). For j ∈ g s, X t (0) is defined as the first K i successive integers under the normalization restriction, and the initial value of A i is calculated as the vector A ( 0  =  Z (0 T X * (0 ) . Given these initial values, PRINCALS as provided in Michailidis and de Leeuw ( 1998) iterates the following two steps: ● Model parameter estimation step: Calculate Z (t+1) by Columnwise center and orthonormalize Z ( t+ ) . Estimate A ( t+ )  for the single variable j by A ( + 1)  =  Z ( + 1) T X + ( ) / X + ( ) T X + ( ) ● Optimal scaling step: Estimate the optimally scaled vector for j ∈ J M by = G j (G j G j )- 1G j Z (t+1) and for j ∈ gs by X*(t+1) = G j (G j G j )- 1 G j Z(t+1)A(t+1)/A(t+1)TA(t+1) under measurement restrictions on each of the variables.
4.  The vε acceleration of the ALS algorithm
We briefly introduce the v ε algorithm of Wynn (1962) used in the acceleration of PCA.ALS. The v ε algorithm is utilized to speed up the convergence of a slowly convergent vector sequence and is very effective for linearly converging sequences. Kuroda and Sakakihara (2006) proposed the ε-accelerated EM algorithm that speeds up the convergence of the EM sequence via the νε algorithm and demonstrated that its speed of convergence is significantly faster than that of the EM algorithm. Wang et al. (2008) studied the convergence properties of the ε-accelerated EM algorithm. Let (Y(0) l >0 = {Y(0) , Y(1) , Y(2) , . . .} be a linear convergent sequence generated by an iterative computational procedure and let {Y(');t>o = {Y(0', Y(1', Y(2', . . .} be the accelerated sequence of {Y(');t>o . Then the v ε algorithm generates {Y(');t>o by using
where [Y]-'  =  Y/||Y|| 2  and ||Y|| is the Euclidean norm of Y. For the detailed derivation of Eq. (7). see the Appendix. When {Y(')},-,o converges to a limit point Y(o) of {Y(')},-,o, it is known that, in many cases, {Y(')},-,o generated by the v ε algorithm converges to Y(∞) faster than {Y(t)}t>0. M. Kuroda et al. / Computational Statistics and Data Analysis 55 (2011) 143–153
We assume that { X*( t ) t≥0 generated by PCA.ALS converges to a limit point X*( ∞ ) . Then v ε -PCA.ALS produces a faster convergent sequence {X*(')t≥0 of {X*(')t≥0 by using the v ε algorithm and enables the acceleration of convergence of PCA.ALS. The general procedure of vε-PCA.ALS iterates the following two steps: ● PCA.ALS step: Compute model parameters A ( ) and Z ( ) and determine optimal scaling parameter X*( +1 ) . ● Acceleration step: Calculate X*(t-1) using {X*(t-1) , X*(t) , X*(t+1)} from the v ε algorithm: [vec(X*(t −1) = vec X*(t) + = (X*T X*T . . . X*T) T , and check the convergence by where vec X*  = X*(t-1) _ X*(t-2)) where δ is a desired accuracy. Before starting the iteration, we determine initial data X*(0) satisfying the restriction (3) and execute the PCAALS step twice to generate {X*(0) , X*(2)}. vε-PCA.ALS is designed to generate ( X*( t ) t>o converging to X*( ∞ ). Thus the estimate of X* can be obtained from the final value of {X*(t)t>o when v ε -PCAALS terminates. The estimates of Z and A can then be calculated immediately from the estimate of X* in the Model parameter estimation step of PCA.ALS. Note that X*( − 1) obtained at the t-th iteration of the Acceleration step is not used as the estimate X*( +1) at the (t + 1)-th iteration of the PCA.ALS step. Thus vε-PCA.ALS speeds up the convergence of {χ*(')}t≥o without affecting the convergence properties of ordinary PCA.ALS.
5.  Numerical experiments
In this section, we examine the performance of ve-PCA.ALS by employing simulated and real data, and demonstrate the advantage of vε-PCA.ALS over PCA.ALS in terms of the number of iterations and CPU time (in seconds) required for convergence. All computations are performed with the statistical package R (R Development Core Team, 2008) executing on a Celeron 2.8 GHz computer with 1 GB of memory. CPU times taken are measured by the function proc.time. 1 Here we use PRINCIPALS as PCA.ALS. For all experiments, δ for convergence of vε-PRINCIPALS is set to 10~ 8 and PRINCIPALS terminates when |θ*(t+1) — θ*(t)| < 10~ 8 , where θ*(t) is the t-th update of θ* calculated from Eq. (4).
5.1. Simulated data: comparison of the number of iterations and CPU time
In this experiments, we study how much faster v ε -PRINCIPALS converges than PRINCIPALS. We apply these algorithms to a random data matrix of 60 observations on 40 variables with 10 levels and measure the number of iterations and CPU time taken for r = 1, . . . , 15. The procedure is replicated 50 times. M. Kuroda et al. / Computational Statistics and Data Analysis 55 (2011) 143–153 Table 1 shows the summaries of the number of iterations and CPU time taken by PRINCIPALS and vs-PRINCIPALS. We see from the table that PRINCIPALS requires more iterations and takes a longer time than νε-PRINCIPALS. Table 2 reports the iteration and CPU time speed-ups for comparing the speed of convergence of PRINCIPALS with that of vs-PRINCIPALS. The iteration speed-up is defined as the number of iterations required for PRINCIPALS divided by the number of iterations required for v ε -PRINCIPALS. The CPU time speed-up is calculated similarly to the iteration speed-up. We see from the second to fifth columns in Table 2 that v ε -PRINCIPALS converges about 3 times faster than PRINCIPALS in terms of the mean number of iterations and the mean CPU time. The boxplots of the iteration and CPU time speed-ups in Figs. 1 and 2 also indicate that v ε -PRINCIPALS converges 2 to 4 times faster than PRINCIPALS for both. The advantage of the νε acceleration is very obvious. The computation time per iteration of v s -PRINCIPALS is greater than that of PRINCIPALS due to computation of the Acceleration step. We calculate the CPU time ratio of PRINCIPALS to νε-PRINCIPALS. The ratio is defined by
CPU time of vε-PRINCIPALS the number of iterations of PRINCIPALS the number of iterations of vε-PRINCIPALS The values of the mean CPU time ratio reported in the sixth and seventh columns of Table 2 demonstrate that the computation time of the PRINCIPALS step accounts for 90–95% of that of v ε -PRINCIPALS except when r  =  1 and thus considerable time is required for computation of the PRINCIPALS step . This means that the computational cost of the νε acceleration in the Acceleration step is less expensive than that of PRINCIPALS. Therefore the v ε acceleration enables a significant reduction in the numbers of iterations and computation times with less computational effort.
5.2. Real data
5.2.1. Studies of convergence We use data obtained in teacher evaluation by students. The data are obtained from 56 students and consisted of 13 categorical variables with 5 levels each; the lowest evaluation level is 1 and the highest 5. The second and third columns in Table 3 and Fig. 3 show the numbers of iterations of PRINCIPALS and v ε -PRINCIPALS. The fourth and fifth columns in Table 3 give CPU times taken by PRINCIPALS and νε-PRINCIPALS. The values of the iteration speed-up in the sixth column of Table 3 also indicate that v ε -PRINCIPALS converges 3–4 times faster than PRINCIPALS. The CPU time speed-up in the seventh column of the table demonstrates that the computation times of νε-PRINCIPALS are 2.5–3.5 times shorter than those of PRINCIPALS. To check the convergence of both algorithms, we calculate the maximum absolute errors defined by Err(X*(t) = max |vec(X*(t) — X*(∞))|, Err(X*(t) = max |vec(X*(t) — X*(∞))|, where X*( ∞ ) is the final value of {X*( )} t≥o. Fig. 4 illustrates the traces of {Err(X*( ) } t≥o and {Err(X*( ) } t≥o for r = 2, 4, 6 and 8. The figures demonstrate that ( ˙ X*( ˙ ) t ≥ 0 converges to X*( ∞  significantly faster than { X*( ˙ ) t ≥ 0. For each r, Err( X*( ˙ ) M. Kuroda et al. / Computational Statistics and Data Analysis 55 (2011) 143–153
Table 3 The numbers of iterations and CPU times of PRINCIPALS and νε -PRINCIPALS, and the iteration and CPU time speed-ups.
after the convergence of νε-PRINCIPALS is smaller than 10~4; this therefore means that {X*(');:=o converges to the same M. Kuroda et al. / Computational Statistics and Data Analysis 55 (2011) 143–153 Table 4
Next we measure the rates of convergence of PRINCIPALS and νε-PRINCIPALS. The rates of convergence of these algorithms are assessed as
If the inequality 0 < ↑ < τ < 1 holds, we say that { x*( ′ ) t ≥ 0 converges faster than { x*( ′ ) t ≥ 0 . We also investigate the speed of convergence of νε-PRINCIPALS by
From Brezinski and Zaglia (1991) if ( X*( t ) t ≥ 0 converges to th same limit point X*( ∞ ) as { X*( t ) t ≥ 0 and Eq. (8) holds, we say that (X*(') t ≥ 0 accelerates the convergence of (X*(') t ≥ 0. The second and third columns of Table 4 provide the rates of sey that (t ' ) { 2, 6 and t , Fig. 5 shows the traces of (τ ' )}, ( τ ' ή ) and { ρ ( ' )} for r = 2, 4, 6 and 8. From these columns in Table 4 and Fig. 5, we see that (X*(t) t=0 converges faster than (X*(t) t=0 in comparison between τ and τ for each r and thus conclude that νε-PRINCIPALS significantly improves the rate of convergence of PRINCIPALS. The behavior of {م{م} for each r in Fig. 5 also indicates that v e -PRINCIPALS accelerates the convergence of {X*(t)}. For other values of r not drawn in the figure, ρ is reduced to zero. 5.2.2.  Application to variable selection problems We consider a situation where more computation is required and investigate how much the ε acceleration algorithm improves the computational efficiency. In order to do this, we apply the algorithm to the variable selection problem in Modified PCA (M.PCA) for qualitative data (Mori et al., 2007). Suppose that we wish to obtain the best subset of variables consisting of q variables among p original variables (1 ≤ q < p). Let X be an n × p original data set including categorical variables, Xv, an n × q submatrix of X and Xv, the n × (p − q) remaining submatrix of X. Variable selection using the criterion in M.PCA is to select the best subset of q variables which has the largest value of the proportion P =
[(S 21 + S 12 S 21 ) — D r S 11 ]A = 0,
(9) and S, S 11 , and S 12 (= S 1 ) are covariance matrices of X, of X 1 , and of X v 1 and X v 2 , respectively (Tanaka and Mori, 1997). Based on the above, variable selection in M.PCA for qualitative data (Mori et al., 2007) uses the PRINCIPAL technique, in which Eq. (5) is replaced by Eq. (9), to obtain λ i in the Model parameter estimation step and S  =  X*TX*/n in the Optimal scaling step for each q. M. Kuroda et al. / Computational Statistics and Data Analysis 55 (2011) 143–153 Here we employ Backward elimination and Forward selection as cost-saving stepwise selection procedures in which only one variable is removed or added sequentially (see http://mo161.soci.ous.ac.jp/vaspca/theory/procedureE.html for the details of procedures). In Backward elimination, to find the best subset of q — 1 variables, we perform M.PCA with the ALS algorithm (PRINCIPALS for M.PCA) for each of q possible subsets of the q — 1 variables among q variables selected in the previous selection step. The total number of estimations by PRINCIPALS for M.PCA from q = p − 1 to q = r is, therefore, p + (p − 1) + · · · + (r + 1) = (p − r)(p + r + 1)/2. Similarly, in Forward selection, the total number of estimations by PRINCIPALS for M.PCA from q = r to q = p − 1 is p Cr + (p − r) + (p − (r + 1)) +۰· · + 2 = p Cr + (p − r − 1)(p − r + 2)/2. In the experiments, we use the same data as in Section 5.2.1 and apply both variable selection procedures using PRINCIPALS and v ε -PRINCIPALS to the data. The selection criterion is the proportion P. For each r ( r = 1, . . . , 12), we measure the total number of iterations and total CPU time when all the best subsets of q variables (q = r, . . . , 12) are selected. The values in the second to fifth columns of Table 5 for Backward elimination and those of Table 6 for Forward selection indicate that the number of iterations of PRINCIPALS is very large and a long time is taken for convergence, while ve-PRINCIPALS converges considerably faster than PRINCIPALS. We can see from the sixth and seventh columns in these tables that PRINCIPALS requires the number of iterations 3–4 times greater and CPU time 2.5–3.4 times longer than those of M. Kuroda et al. / Computational Statistics and Data Analysis 55 (2011) 143–153
Table 6 The numbers of iterations and CPU times of PRINCIPALS and νε-PRINCIPALS, and their speed-ups in applying the variable selection problems using Forward
vs-PRINCIPALS. The values of the iteration and CPU time speed-ups demonstrate that the vs acceleration algorithm works well to accelerate the convergence of {X*(t)} t ≥ 0 and consequently results in greatly reduced computation times.
6. Conclusion
In this paper, we presented vε-PCA.ALS that accelerates the convergence of PCA.ALS by using the vε algorithm. The algorithm generates the νε accelerated sequence ( X*( t )} using ( X*( t ) t≥o but it does not modify the estimation equations of PCA.ALS. Therefore the algorithm enables an acceleration of the convergence of PCA.ALS while still preserving the stable convergence property of PCA.ALS. The ν ε algorithm in itself is a fairly simple computational procedure and, at each iteration, it requires only O(np) arithmetic operations. For each iteration, the computational complexity of the νε algorithm may be less expensive than that for computing a matrix inversion and for solving the eigenvalue problem in PCA.ALS. The most appealing points of vε-PCA.ALS are that, if {X*(')t≥o convergences to a limit point X*( ∞ ), then {X*(')t≥o converges to X*( ∞ ) of {X*( i ) t=0 and its speed of convergence is faster than that of PCA.ALS. In the application of v=- PCA.ALS, we do not have to take into account the structure of optimally scaled data sets. The numerical experiments employing simulated and real data demonstrated that νε accelerated PRINCIPALS improves the speed of convergence of ordinary PRINCIPALS and significantly speeds up the convergence of {X*(')t==0 in terms of the number of iterations and the computation time. Moreover, in the variable selection problems for finding a suitable variable set using Backward elimination and Forward selection procedures, the νε acceleration enables greatly reducing computation times. We intend to evaluate theoretically the convergence properties of v ε -PCA.ALS. As listed in Young (1981) and Krijnen (2006), there exist many other ALS types of algorithms. We attempt to speed up the convergence of their ALS algorithms by incorporating the νε acceleration.
Acknowledgements
The authors would like to thank the editor and two referees whose valuable comments and kind suggestions that led to an improvement of this paper. This research is supported by the Japan Society for the Promotion of Science (JSPS), Grant-in-Aid for Scientific Research (C), No 20500263. M. Kuroda et al. / Computational Statistics and Data Analysis 55 (2011) 143–153
Appendix.  The vε algorithm
Let Y(t) denote a vector of dimensionality d that converges to a vector Y( ∞ ) as t → ∞. Let the inverse [Y]−¹ of a vector Y be defined by where ||Y|| is the Euclidean norm of Y. In general, the νε algorithm for a sequence {Υ {t'}t≥0 starts with and then generates a vector ε(t,k+1)
For practical implementation, we apply the νε algorithm for k = 1 to accelerate the convergence of {Υ( 1 )t≥o. From Eq. (10), we have
for k = 0. Then the vector ε (t,2) becomes as follows: