Fast Computation of Robust Subspace Estimators
Holger Cevallos-Valdiviezo Stefan Van Aelst

















mymainaddress,mysecondaryaddress]Holger Cevallos-Valdiviezo
holgceva@espol.edu.ec

myaddress]Stefan Van Aelstmycorrespondingauthor
[mycorrespondingauthor]Corresponding author: KU Leuven, Department of Mathematics, Celestijnenlaan 200B, 3001 Leuven, Belgium. Phone: +3216372383, fax: +3216327998.
Stefan.VanAelst@kuleuven.be

[mymainaddress]ESPOL Polytechnic University, Escuela Superior Politécnica del Litoral, ESPOL, (Facultad de Ciencias Naturales y Matemáticas, FCNM), Campus Gustavo Galindo Km. 30.5 Vía Perimetral, P.O. Box 09-01-5863, Guayaquil, Ecuador
[mysecondaryaddress]Ghent University, Department of Applied Mathematics, Computer Science and Statistics, Krijgslaan 281 S9, 9000 Gent, Belgium
[myaddress]KU Leuven, Department of Mathematics, Celestijnenlaan 200B, 3001 Leuven, Belgium


Dimension reduction is often an important step in the analysis of high-dimensional data.
PCA is a popular technique to find the best low-dimensional approximation of high-dimensional data. However, classical PCA is very sensitive to atypical data. 
Robust methods to estimate the low-dimensional subspace that best approximates the regular data have been proposed.
However, for high-dimensional data these algorithms become computationally expensive. Alternative algorithms for the robust subspace estimators are proposed that are better suited to compute the solution for high-dimensional problems. The main ingredients of the new algorithms are twofold. First, the  principal directions of the subspace are estimated directly by
iterating the first order solutions corresponding to the estimators. Second, to reduce the computation time even further five robust deterministic values are proposed to initialize the algorithms instead of using random starting values. 
It is shown that the new algorithms yield robust solutions and the computation time is largely reduced, especially for high-dimensional data. 



Deterministic algorithm High-dimensional data Least trimmed squares M-scale Principal Component Analysis










§ INTRODUCTION


Principal component analysis (PCA) is a popular exploratory tool for multivariate data.
In particular, PCA is extremely useful to find a low-dimensional representation of high-dimensional data that yields the best possible approximation to the original data. 
Classical PCA minimizes the squared euclidean distances between the original observations and their orthogonal projections onto the lower dimensional subspace. However, classical PCA is very sensitive to atypical data due to the use of quadratic loss. Therefore, several approaches to robustify PCA have been proposed. 

The earliest and easiest approach to robust PCA consists of taking the eigenvectors and eigenvalues of a robust scatter estimator instead of the standard sample covariance matrix <cit.>. However, this approach cannot be used for high-dimensional data because calculating high-dimensional robust scatter matrices is computationally complex or even infeasible if the sample size is small compared to the dimension. Moreover, while the efficiency of robust scatter estimators increases with dimension  this comes at the expense of a loss of robustness. Therefore, <cit.> introduced spherical PCA which uses the covariance matrix of the data projected onto the unit sphere and is fast to compute. 

A second approach to robust PCA sequentially looks for univariate directions that maximize a robust estimator of scale and are orthogonal to each other. This robust projection pursuit (PP) approach has been studied by e.g. <cit.> and <cit.>. 

Instead of looking for one direction at a time as in PP, one can robustly estimate a lower-dimensional subspace directly <cit.>. ROBPCA <cit.> seeks for a lower-dimensional subspace via a multiple-step procedure. Briefly, the first step aims to identify a subset of n/2≤ h<n observations based on the Stahel-Donoho outlyingness. Classical PCA is then applied on this h-subset and the optimal dimension of the subspace is determined. Next, the data points are projected onto the lower-dimensional subspace. Finally, the reweighted MCD estimator is computed on the projected data to obtain the estimates for the principal directions of the subspace. 

Principal Component Pursuit method (PCP) aims to decompose the data matrix into a low-rank component and a gross outlier component <cit.>. See e.g.  <cit.> and <cit.> for related work.
However, PCP may fail to detect outliers in the orthogonal complement of the subspace <cit.>. 
Therefore, <cit.> introduced a robust orthogonal complement
(ROCPCA) approach to deal with orthogonal complement outliers. The authors model the projected data onto the orthogonal complement subspace with a mean term, a sparse outlier matrix and a noise term. The sparse outlier matrix identifies orthogonal outliers. More recently, <cit.> combined ideas from PCP and ROCPCA to construct a reinforced robust PCP (RRPCP) method that considers outliers in both the observation space and in the orthogonal complement subspace. RRPCP decomposes the data matrix into a low-rank component, a sparse outlier matrix to represent outliers in the observation space, a noise term, a mean term and another sparse outlier matrix to represent outliers in the orthogonal complement space. The two latter components are represented in the orthogonal complement space and then transformed back to the observation space. However, according to their model formulation ROCPCA and RRPCP have to estimate the directions of the orthogonal complement space, which becomes computationally intensive for high-dimensional data.  

<cit.> proposed to robustly estimate the best lower-dimensional approximation by minimizing either an M-scale or a least trimmed squares (LTS) scale of the Euclidean distances corresponding to the observations. He also characterized the solutions by the orthogonal complement directions and showed that these directions correspond to the eigenvectors associated with the smallest eigenvalues of a weighted covariance matrix. Based on this characterization <cit.> proposed an iterative algorithm to compute the robust subspace estimators. 
The robustness of these subspace estimators has been widely investigated empirically. See e.g.   
<cit.> and <cit.> for the M-scale estimator, and <cit.> and <cit.> for the LTS-scale estimator. Moreover, <cit.> also provides a thorough theoretical study of the properties of the estimator based on the LTS scale.  

However, in case of a low-dimensional approximation for high-dimensional data Maronna's orthogonal complement algorithm requires to decompose a high dimensional covariance matrix and needs a large number of its eigenvectors to characterize the solution. This makes computing the subspace estimators time consuming or even infeasible in high dimensions.  Therefore, we propose an algorithm for the robust subspace estimators of <cit.> that directly calculates principal directions of the low-dimensional subspace. 

The main ingredients of our new algorithm are twofold. First, we use the first order conditions corresponding to the estimator to update the principal directions of the subspace iteratively. 
This approach only requires low-dimensional vector and matrix operations rather than manipulating high-dimensional covariance matrices. Second, instead of using random starting values, similarly to <cit.> we propose five robust deterministic values to initialize the algorithm. 
These starting values yield robust fits that are usually close to the sought after robust solution, so that convergence occurs quickly. 

The remainder of the manuscript is organized as follows. Section <ref> reviews the robust subspace estimators based on an M-scale or the LTS-scale.
Our definition is equivalent to the definition in <cit.> but characterizes the solution by the principal directions of the subspace which better serves our needs for the development of the new algorithm in Section <ref>. In Section <ref> we introduce the robust deterministic values to initialize the algorithm. 
In Section <ref> we compare the solutions calculated with the new algorithm to those obtained with the original algorithm by means of an extensive simulation study. We also include other robust subspace estimators in the simulations to compare their robustness properties. Since the deterministic starting values are not orthogonally equivariant, in Section <ref> we empirically investigate the effect of orthogonal transformations of the data on our algorithm with deterministic starting values. We compare all methods in terms of their computation time in Section <ref>, while in Section <ref> we extend the simulation study in Section <ref> to higher-dimensions.
Section <ref> contains a real data illustration and Section <ref> presents our final conclusions. 



§ ROBUST SUBSPACE ESTIMATORS
 
Consider a data matrix 𝐗=(𝐱_1…𝐱_n)^T∈ℝ^n× p which contains the measurements of p variables for n observations.
The goal is to approximate the n observations 𝐱_i by points 𝐱_i that lie in a q-dimensional subspace. That is, 
𝐱_i≡𝐱_i( 𝐁_q,𝐀_q, 𝐦)= 𝐦 + 𝐁_q𝐚_i
for some 𝐦∈ℝ^p, 𝐀_q=(𝐚_1,…,𝐚_n)^T∈ℝ^n × q and orthogonal matrix 𝐁_q∈ℝ^p × q, i.e. 𝐁_q^T𝐁_q = 𝐈_q.
Let 𝐛_1,…,𝐛_p denote the rows of 𝐁_q. 
The Euclidean distance between 𝐱_i and its approximation 𝐱_i is denoted by d_i(𝐁_q,𝐀_q, 𝐦)= d_i=‖𝐫_i‖, where 𝐫_i = 𝐱_i -  𝐱_i.
 
<cit.> proposed to robustly estimate the optimal subspace by
minimizing a robust scale estimator of the Euclidean distances d_i(𝐁_q,𝐀_q, 𝐦). 
Note that if the nonrobust standard deviation is used, then the estimator minimizes the sum of squared Euclidean distances and thus the classical PCA solution is retrieved. Maronna's estimators are thus robust extensions of classical PCA dimension reduction.   

Although the Euclidean distance between each observation 𝐱_i and its projection 𝐱_i onto the q-dimensional subspace is measured in the p-q dimensional orthogonal subspace in <cit.>, this is equivalent to our current formulation in the p-dimensional space. 




  §.§.§ Subspace S-estimator


The subspace S-estimator (𝐁_S,𝐀_S, 𝐦_S)  is obtained by minimizing an M-scale of the Euclidean distances d_i(𝐁_q,𝐀_q, 𝐦). That is, (𝐁_S,𝐀_S, 𝐦_S) is a solution  of 

    min_𝐁_q,𝐀_q, 𝐦σ_M(𝐝(𝐁_q,𝐀_q, 𝐦)),

over all 𝐦∈ℝ^p, 𝐀_q=(𝐚_1,…,𝐚_n)^T∈ℝ^n × q and orthogonal matrices 𝐁_q∈ℝ^p × q. For any 𝐝(𝐁_q,𝐀_q, 𝐦)
=(d_1,d_2,…,d_n), the corresponding M-scale σ_M(𝐝(𝐁_q,𝐀_q, 𝐦)) is defined as the solution in s of 

    1/n∑_i=1^nρ( d_i/s) = b,

where ρ is an even function that is differentiable and non-decreasing on the positive real line with ρ(0) = 0  (see e.g. ). 


Similarly as in  <cit.>  
first order conditions for the subspace S-estimator can be obtained by implicitly differentiating the M-scale in (<ref>). 
The solutions of these first order conditions can be written as








    𝐚_i   = (∑_j=1^p  𝐛_j 𝐛_j^T)^-1(∑_j=1^p    (x_ij - m_j)  𝐛_j), 
    
    𝐛_j   = (∑_i=1^n w_i 𝐚_i 𝐚_i^T)^-1(∑_i=1^n w_i  (x_ij - m_j)  𝐚_i),

 and

    m_j = ∑_i=1^n w_i  (x_ij - 𝐚_i^T𝐛_j)/∑_i=1^n w_i,

with weights 

    w_i≡ w(d_i) = ρ^'( d_i/σ_M)σ_M/d_i,

for the rows of 𝐀_q, i=1,…,n, and the rows of 𝐁_q and components of 𝐦, j=1,…,p, respectively. 
Expressions (<ref>)-(<ref>) suggest an iterative reweighted least squares procedure.
Initial values 𝐁_q^(0), 𝐀_q^(0) and 𝐦^(0) can be used to obtain starting weights w_i^(1) from expression (<ref>). 
Based on 𝐁_q^(0) and 𝐦^(0) one can obtain updated scores 𝐀_q^(1) from (<ref>). Next, based on 𝐀_q^(1), 𝐦^(0) and weights w_i^(1) one can obtain 𝐁_q^(1) from (<ref>). Finally, based on 𝐀_q^(1), 𝐁_q^(1) and weights w_i^(1) one can obtain 𝐦^(1) from (<ref>). This is the first iteration. 

This procedure is repeated until we obtain a stable solution. Note that the resulting estimate of the matrix 𝐁_q is not necessarily orthogonal. Therefore, we orthogonalize this estimate at the end of the procedure and update the scores accordingly. Using the orthogonality of 𝐁_q,  it can be seen from (<ref>) that 𝐚_i=𝐁_q^T(𝐱_i-𝐦).
Hence, once 𝐁_q and 𝐦 are known, the corresponding scores 𝐚_i of the observations are easily obtained. 
By combining this result  with (<ref>) we can also obtain that 𝐦= ∑_i=1^n w_i 𝐱_i/ (∑_i=1^n w_i ).
Note that if we put w_i=1 for all observations, then the procedure reduces to an  alternating least squares algorithm whose solution yields the classical PCA estimates <cit.>. 
Finally, by combining these expressions and using that 𝐁_S^T𝐁_S=𝐈_q
it can also be derived that the subspace S-estimators (𝐁_S, 
𝐦_S) satisfy the equation

    ∑_i=1^n w_i (𝐱_i - 𝐦)(𝐱_i - 𝐦)^T𝐁_q = 𝐁_q Λ.

where Λ=𝐁^T∑_i=1^n w_i (𝐱_i - 𝐦)(𝐱_i -  𝐦)^T𝐁.
From (<ref>) and (<ref>) it follows that the columns of 𝐁_S correspond to the first q eigenvectors of the weighted covariance matrix 

    𝐂(𝐦_S, 𝐁_S) = 1/n∑_i=1^n w_i(𝐱_i - 𝐦_S)(𝐱_i - 𝐦_S)^T,

which coincides with expression (9) of <cit.>. 



  §.§.§ Subspace LTS-estimator


The subspace LTS-estimator (𝐁_LTS,𝐀_LTS, 𝐦_LTS) is obtained by minimizing the LTS-scale of the Euclidean distances d_i(𝐁_q,𝐀_q, 𝐦). That is, (𝐁_LTS,𝐀_LTS, 𝐦_LTS) is a solution  of 

    min_𝐁_q,𝐀_q, 𝐦σ_LTS(𝐝(𝐁_q,𝐀_q, 𝐦)),

over all 𝐦∈ℝ^p, 𝐀_q∈ℝ^n × q and orthogonal matrices 𝐁_q∈ℝ^p × q. For any 𝐝(𝐁_q,𝐀_q, 𝐦)
=(d_1,d_2,…,d_n), the corresponding LTS scale is defined as

    σ^2_LTS(𝐝) = 1/h∑_i=1^h d_(i:n)^2,

where d_(1:n)≤…≤ d_(n:n) are the ordered Euclidean distances and h=n-⌊ nα⌋ with 0 ≤α < 0.5. 
A fraction α of the observations is not taken into account when calculating the LTS estimator and thus 
α determines the robustness of the estimator.

Similarly as in <cit.> it can be shown that the LTS solution satisfies the expressions (<ref>)-(<ref>) with weights now given by

    w_i≡ w(d_i) =
    1    if   d_i≤  d_(h:n),
    
    0    otherwise.

With these weights the LTS estimator (𝐁_LTS,𝐦_LTS) also satisfies (<ref>) and (<ref>) which again coincides with expression (9) in  <cit.>.



§ THE ALGORITHM


<cit.> characterizes a q-dimensional subspace by an equation 𝐁_p-q^T𝐱=𝐚 with  𝐚∈ℝ^p-q and 𝐁_p-q∈ℝ^p × (p-q) an orthogonal matrix. 
His algorithm can be summarized by the following steps. First, 50 random orthogonal matrices 𝐁_p-q are generated to initialize the algorithm. The corresponding optimal value of 𝐚 that minimizes the robust scale can then be calculated easily. A few iterative improvement steps are then applied to each of these starting values. In each iteration the weights w_i of the observations corresponding to the current solution are updated. Based on the updated weights, a new solution is then obtained by calculating the smallest p-q eigenvectors of the weighted covariance matrix in (<ref>). After two iterations, the 10 best solutions  are selected and these are iterated further until convergence (with a maximum of 10 iterations).

If p is large and the subspace dimension q is small, Maronna's algorithm requires the storage of a high-dimensional covariance matrix and the calculation of a large number (namely, p-q) of its eigenvectors which becomes computationally demanding. To improve the computation time and reduce the memory load, we instead propose to directly calculate the q basis directions of the subspace by iterating expressions (<ref>)-(<ref>). Extensive experiments have shown that iterating these expressions only a few times suffices to obtain close approximations to the first eigenvectors of the weighted covariance matrix in (<ref>). Note that our iterations of these expressions only require operations with q-dimensional vectors and matrices and thus will be more suitable for high-dimensional settings.

Algorithm <ref> contains a detailed description of the main part of our modified algorithm in pseudo-code.  The algorithm requires initial values of 𝐁_q and 𝐦 as input. It also depends on tuning parameters N_1, N_2, N_3 and tol. 
tol specifies the precision with which the solution is calculated. The tuning parameters N_1 and N_2 play the same role as in Maronna's original algorithm. That is, for each initial orthogonal matrix 𝐁_q first N_1 iterations are performed to improve the corresponding estimates of the location 𝐦 and the scores matrix 𝐀_q while keeping 𝐁_q fixed. In the next N_2 iterations, the estimates of all three quantities 𝐁_q, 𝐦 and 𝐀_q are updated. Since these updates are now calculated by iterating expressions (<ref>)-(<ref>), an additional tuning parameter N_3 specifies how often these expressions are iterated.  Using a similar proof as in <cit.>, it can easily be shown that either the M-scale σ_M or LTS scale σ_LTS decreases in each iteration of our algorithm. 

Moreover, the performance of our algorithm is similar to Maronna's algorithm if we use the same settings for the tuning parameters and start from the same initial values for 𝐁_q and 𝐦. 



For the subspace S-estimator we have used the popular Tukey biweight loss function ρ(y) = min(3y^2 - 3y^4 + y^6,1) to calculate the M-scales.
Two standard choices for the tuning parameters are 
b=0.5 which yields the maximal breakdown point (BDP) of 50% and 
b=0.2426 which yields a better compromise between efficiency and robustness (breakdown point ≈ 25%). 
For the subspace LTS estimator the  trimming fraction α determines the breakdown point. The most common choices are 
α=0.5 (50% BDP) and α=0.25 (25% BDP).
See <cit.> for more details on the breakdown point.  



§ STARTING VALUES


It is well-known that the objective functions in (<ref>) and (<ref>) are nonconvex and thus may have several local minima. Not all of these minima correspond to robust solutions.
A standard approach is to aim for the global minimum in (<ref>) or (<ref>) by generating many starting values. The solution that corresponds to the smallest value of the objective function that is reached by iterating these starting values is then the approximation for the global optimum <cit.>.
In Maronna's algorithm 50 random initial orthogonal matrices 𝐁_q are generated while the initial estimate for the location 𝐦 is the coordinatewise median of the data.
Each of these starting values is iterated N_1 = 3 times with 𝐁_q fixed, followed by N_2 = 2 iterations to improve the estimates of 𝐁_q, 𝐀 and 𝐦 together. Then, the 10 best solutions (with the lowest scale) are selected and these are iterated further until convergence (tol=0.001) with a maximum of N^'_2 = 10 iterations. Extensive experiments by <cit.> revealed that the algorithm shows good performance with these choices for the tuning parameters. Following <cit.> we denote the solutions obtained with his algorithm by S-M when the M-scale is minimized and by S-L when the LTS scale is minimized. 

Empirical comparisons have confirmed that if we use the same starting values for our new algorithm as well as the same settings for the tuning parameters, then we obtain a similar performance as Maronna's algorithm. Note that our algorithm requires the additional tuning parameter N_3 which determines the number of iterations to calculate the updated estimates by iterating expressions (<ref>)-(<ref>). Extensive experiments showed that it suffices to use N_3=3  iterations to obtain good approximations that yield stable results. Note that while Maronna's algorithm uses the fast-to-compute coordinatewise median as starting value for 𝐦, we prefer to use as default the orthogonal equivariant spatial median which can also be computed efficiently <cit.>. Moreover, we generate random orthogonal matrices 𝐁_q as starting values by the method of <cit.> which consists of orthogonalizing a matrix of normal random numbers while <cit.> instead orthogonalizes a matrix of uniform random numbers. We denote the solutions of our algorithm with random orthogonal starting values by RsubS for the subspace S-estimator and by RsubLTS for the subspace LTS-estimator.



While our adaptation of Maronna's algorithm is indeed faster for high-dimensional data (see results in Section 7), the computation time remains considerable because a sufficient number of random starting values is needed to obtain a robust solution. While the default setting is to use a fixed number of 50 random orthogonal matrices as starting values, it is clear that the search space increases dramatically with increasing dimension. Hence, it can be expected that for high-dimensional data (many) more random starting values may be needed to find a stable robust solution. Many algorithms for robust estimators encounter the same issue. To address this issue for the calculation of the minimum covariance determinant (MCD) estimator of multivariate location and scatter in high dimensions,  <cit.> introduced a deterministic MCD algorithm <cit.>. The main idea is to replace the random starting values by a few well-chosen robust starting values. 
These robust starting values should be fast to compute while at the same time they are expected to lie close to a robust minimum of the objective function, which in our case is given by (<ref>) or (<ref>). 
Hence, instead of exploring the whole parameter space to find the optimum, a  few robust starting values should point us to that part of the parameter space where robust solutions can be found. 
Since convergence from the robust starting values to their closest local minimum is generally fast as well, this results in an algorithm with a much lower computation time which allows us to handle problems in higher dimensions. 

For the deterministic version of our algorithm, we could use as starting values for 𝐁_q the q largest eigenvectors corresponding to the robust starting values for the p-dimensional scatter matrix proposed by <cit.>. However, we want to avoid having to calculate p-dimensional scatter matrices which may be unstable for large p and consumes a lot of memory.
Therefore, inspired by <cit.> we propose five robust starting values for 𝐁_q which are obtained by the following procedure:



  Step 1.
Robustly standardize each variable X_j by subtracting its median and dividing by its Q_n scale estimate <cit.>. 
Let 
𝐙 = (𝐳_1,…,𝐳_n)^T denote the standardized data matrix and Z_1,…, Z_p the standardized variables. Consider the following transformations of the data:


  * Hyperbolic tangent (sigmoid) transformation: let H_j = tanh(Z_j), j=1,…,p.
 
Set 𝐔_1^'=(H_1,… H_p).

  * Rank transformation: let R_j be the vector of ranks of column X_j, j=1,…,p.

Set 𝐔_2=(R_1,…,R_p).

  * Normal scores: compute normal scores from the ranks R_j, i.e. T_j = Φ^-1[(R_j - 1/3)/(n+1/3)  ], where Φ(.) is the normal cumulative distribution function. 
 
Set 𝐔_3=(T_1,…,T_p).

  * Spatial signs: let 𝐬_i= 𝐳_i/ 𝐳_i, i=1,…,n.

 Set 𝐔_4^'=(𝐬_1, …,𝐬_n)^T. 

Robustly standardize the columns of 𝐔_1^' and 𝐔_4^' obtained above by subtracting their median and dividing by their Q_n scale estimate. Denote these standardized transformed data matrices by 𝐔_1 and 𝐔_4  and take 𝐔_5=𝐙.

  Step 2. Calculate the classical q-dimensional PCA subspace using the alternating least squares algorithm (see Algorithm <ref> below) for 𝐘=𝐔_k.  Let 
𝐁_q^(k) be the resulting estimate of the q PC directions for k=1,…,5. 

  Step 3. 
For each of the score matrices 𝐀_q^(k) = 𝐙𝐁_q^(k), k=1,…,5, select the l=⌈ n/2  ⌉ rows with smallest Euclidean norm. Let ℐ_k ⊂{1,…,n} denote the indices of the selected rows. 

  Step 4. Set  𝐗^(k)=𝐗_ℐ_k for k=1,…,5.
For each of these five reduced data matrices, calculate the classical q-dimensional PCA subspace by applying Algorithm <ref> on 𝐘=𝐗^(k) - 𝐦^(k) where 𝐦^(k)= 1/l∑_i=1^l𝐱_i^(k). The resulting estimates of 𝐁_q together with 𝐦^(k) are the starting values for Algorithm <ref>.






The alternating least squares algorithm for classical PCA, summarized in Algorithm <ref>, calculates the q-dimensional PCA subspace without performing singular value decomposition of a high-dimensional matrix, and has the advantage that it only requires operations with q-dimensional vectors and matrices in the iterative updating of the estimates in step 4.i. 
Note that for simplicity Algorithm <ref> uses a centered data matrix as input data. 
Similarly as for Algorithm <ref>, we found that N_3 = 3 iterations in step 4 of Algorithm <ref> were enough to obtain close approximations to the first eigenvectors of the classical covariance matrix. The q-dimensional canonical basis or a random orthogonal matrix can be used to initialize 𝐁_q. In fact, we found that Algorithm <ref> converges to the global solution independent of the initial value that is chosen for 𝐁_q. 

The procedure to obtain deterministic starting values can be summarized in two parts. The first part  (first three steps) aims to identify subsets of the data with only regular observations. Transformations of the standardized data are used in Step 1 to reduce the effect of outliers. In Step 2 naive robust estimates of the subspace are calculated from the transformed data which are used to project the data in Step 3 before obtaining clean subsets. This projection makes the algorithm better suitable for correlated data. In Section <ref> we investigate by means of a simulation study the effect of orthogonal transformations of the data on the estimates obtained with these deterministic algorithms. The second part (last step) calculates the classical PCA subspace corresponding to each of the obtained subsets of the data matrix to obtain promising starting values for Algorithm <ref>. Similarly as for orthogonal starting values, each of these five starting values is iterated N_2 = 2 times to improve the estimates of 𝐁_q, 𝐀 and 𝐦. Then, the best solution (with the lowest scale) is selected and iterated further until convergence (tol=1 × 10^-6) with a maximum of N^'_2 = 10 iterations. We denote our algorithm with deterministic starting values by DsubS when the M-scale is minimized and by DsubLTS when the LTS scale is minimized. 



§ PERFORMANCE COMPARISON


To compare our new algorithms to Maronna's original algorithms, we repeat the simulations of <cit.>. We use our algorithms with random orthogonal matrices as starting values as well as with deterministic starting values. For the simulations we choose the 
tuning parameters for the M-scale and LTS scale that yield the maximal breakdown point.

We also consider some other computationally efficient methods for robust subspace estimation. These are robust projection pursuit (PP) <cit.>, Spherical PCA (SPC) <cit.> and ROBPCA <cit.>. For robust PP we use the algorithm of  <cit.> and maximize the LTS scale with α=0.5. Note that PP based on the M-scale was already considered in <cit.>. For ROBPCA we include the reweighting step as in <cit.>. The number of random directions through two data points to compute the measure of outlyingness was fixed to 250 in all experiments as in <cit.>. The parameter α in ROBPCA was set equal to α=0.5 for maximal robustness. 


Following <cit.> we generated M=200 samples of size n=100 and dimension p=10. The regular observations are generated according to N(0,Σ) with Σ = diag(λ_1,λ_2,…,λ_p) where the following two designs are considered for the diagonal elements.


  * An abrupt increase of the eigenvalues: λ_j=1+0.1j for 1 ≤ j ≤ (p-q) and λ_j = 20(1+0.5(j-p+q)) for (p-q+1) ≤ j ≤ p.

  * A smooth increase of the eigenvalues: λ_j=2^j-1 for 1 ≤ j ≤ p.

A fraction ϵ=0%, 10% or 20% of outlying observations is generated from N(k𝐱_0, 0.25Σ), where 𝐱_0 is a vector of length p  with x_0j = 1 for j ≤ (p-q) and 0 otherwise. The value of k runs between 0 and 20 with steps of 0.5. 
All methods are applied for q=2, so the best two-dimensional linear subspace approximation is estimated. 

To compare the methods we measure their prediction performance as in <cit.>. Hence, we measure the proportion of variance in independent regular data that remains unexplained by their approximation according to the estimated subspace. More formally, let 𝐱 be a N(0,Σ) vector independent of the random sample used to obtain 𝐁_q. Then, the variability of 𝐱 around the subspace generated by 𝐁_q equals

    E‖𝐱 - 𝐁_q𝐁_q^T𝐱‖^2 = tr[ Σ] - tr[𝐁_q^TΣ𝐁_q],

and the prediction proportion of unexplained variance is given by

    u_q^pred = E‖𝐱 - 𝐁_q𝐁_q^T𝐱‖^2/tr[ Σ]= 1 - tr[𝐁_q^TΣ𝐁_q]/tr[ Σ].

Since Maronna's algorithms characterize the subspace by an estimate 𝐁_p-q∈ℝ^p × (p-q) of its orthogonal complement, the corresponding prediction proportion of unexplained variability in this case becomes:

    u_q^pred = tr[𝐁_p-q^T Σ  𝐁_p-q]/tr[ Σ].

Comparing this prediction error with the optimal value, given by

    u_q^opt = ∑_j=1^p-qλ_j/∑_j=1^pλ_j,

yields the relative prediction error

    e_pred = u_q^pred/u_q^opt - 1.





















Table <ref> shows the mean of the relative prediction errors over M=200 samples. The values of k included in Table <ref> are those values at which some estimator attains its maximal error. Standard errors are not shown because they were small (below 0.08) in all cases.


For ϵ=20% of contamination, Figures <ref> and <ref> show the mean relative prediction errors as a function of k. From k=10 onwards these prediction errors stabilize so we only show results up to k=10. 

These results show that the estimates obtained by our new RsubS and RsubLTS algorithms are very similar to the results of Maronna's S-M and S-L algorithms, as expected. Moreover, the DsubS and DsubLTS algorithms often succeed in obtaining more robust estimates than their counterparts based on random starts. The advantage of using robust starting values is most pronounced when the outliers are at a small to moderate distance of the majority (k ≤ 6) as can be seen from Figures <ref> and <ref>. Looking at the robust competitors, we can see that PPLTS turns out to be inefficient for ϵ=0%, similarly as for the PP methods considered by <cit.>. Moreover, also for contaminated data PPLTS gives worse results in these settings. Spherical PCA is generally better than PPLTS but its performance decreases considerably for ϵ=20%. ROBPCA performs well in these settings, which is in accordance with findings in <cit.>. 

Overall, ROBPCA and DsubS show the best performance in these simulations.
















 




 

 
























§ ORTHOGONAL EQUIVARIANCE

In the previous section we showed that the RsubS/RsubLTS algorithms yield similar estimates as the original S-M/S-L algorithms while the DsubS/DsubLTS often performed even better. 

Since orthogonal equivariance is important in the context of PCA, we now investigate the effect of orthogonal transformations on our deterministic algorithms by means of a simulation study. Similarly to <cit.>, let 𝐗∈ℝ^n× p denote the data matrix and let 𝐏∈ℝ^p× p be a random orthogonal matrix, i.e. 𝐏^T𝐏=𝐏𝐏^T=𝐈_p. Let 𝐁_q(𝐗) denote an estimate of the first q principal component directions for the data matrix 𝐗. If the estimator is orthogonal equivariant, then it must hold that:

    𝐁_q(𝐗𝐏^T) = 𝐏𝐁_q(𝐗).

To investigate orthogonal equivariance, we therefore compare 𝐁_q(𝐗) to

    𝐁_q^'(𝐗) = 𝐏^T𝐁_q(𝐗𝐏^T).

To measure the effect of the transformation we calculate the standardized last principal angle between the subspaces generated by the columns of 𝐁_q and 𝐁_q^'(𝐗), using the algorithm of <cit.>. We divide this angle by π/2 to obtain values between 0 and 1. Values closer to 0 indicate that the subspaces are more comparable. 

We generated data matrices 𝐗 according to designs a) and b) in Section <ref> with ϵ=20% and estimate the best q=2 linear dimensional subspace. We generated M=200 replicates and for each data replicate we generated a random orthogonal matrix 𝐏 using the method of <cit.>. Figure <ref> displays the mean standardized angles as a function of k for both designs.

We can see that for both DsubS and DsubLTS the effect of the transformation is largest when the outliers are at a small to moderate distance of the regular data. In these settings the DsubS algorithm is also clearly less affected by orthogonal transformations than DsubLTS and the effect quickly disappears when the outlier distance increases.  




§ COMPUTATION TIME


We now compare the computation times of our algorithms to those of the original algorithms of Maronna  and the other methods considered in Section <ref>. For this purpose we generalize the data generating model in design a). More specifically, we generate clean data from the p-dimensional normal distribution with center zero and diagonal covariance matrix. The eigenvalues on the diagonal are set equal to

    λ_j =
    
    1+c_p,q j      for  1 ≤ j ≤ (p-q), 
     
    20(1+0.5(j-p+q))     for  (p-q+1) ≤ j ≤ p.


The constant c_p,q is chosen such that the q main directions explain about 80% of the total variability.
We consider sample sizes n=1000 or 5000 in dimension p=10, 500, 1000 with q=2, 5. The data contain ϵ=20% of contamination, generated as explained in Section <ref> with k=15. 

All algorithms are implemented in  <cit.>. 
For the S-M, S-L, RsubS, RsubLTS, DsubS, DsubLTS and PPLTS algorithms we used our own implementation. For S-M and S-L we use the  function  to calculate the smallest p-q eigenvectors of the weighted covariance matrix in (<ref>). This function uses LAPACK <cit.> routines which are written in FORTRAN. Alternatives such as the  function  for instance could be used, but the speed difference is generally small. For spherical PCA we use the function  and for ROBPCA we use the function , both in  package  <cit.>.  also uses  or  to compute classical PCA on the projected data onto the unit sphere while  uses a FORTRAN implementation to compute the MCD estimator of location and scatter on the projected data and  to compute its decomposition. To allow a fair comparison of computation times between these algorithms and our new algorithms, for our algorithms we implemented the iterative updating of the estimates based on expressions (<ref>)-(<ref>) in C++,  using  package  <cit.>. Note that the implementations of all methods in the comparison use plain , except for the calculation of covariance matrices and/or its decomposition. S-M, S-L, spherical PCA and ROBPCA use FORTRAN implementations for this purpose while our algorithms perform the iterative updates in C++. The algorithms were compared on a single Intel i7 CPU (3.4GHz) machine running Windows 7.  

Table <ref> shows the computation time of the algorithms in seconds, averaged over M=50 replications. 

Not surprisingly, we see that spherical PCA is the fastest to compute because it does not require any iterative 
process <cit.>.  
For n=1000 and q=2 projection-pursuit is relatively fast to compute regardless of the dimension of the data, as can be seen from Table <ref>. However, the computation time grows quickly with increasing subspace dimension (case q=5) and/or sample size (case n=5000).
 The computation time of Maronna's algorithms is reasonable for low dimensional data, but increases quickly when the dimension grows. This result was expected because these algorithms need to compute a large number (p-q) of eigenvectors of a high-dimensional covariance matrix which is very time-consuming. 
Clearly, our implementations of the new algorithms with random starting values are faster than Maronna's original algorithms, but computation time still increases quickly with growing dimension. 
The new algorithms with deterministic starting values are even faster and their computation time increases much slower when the dimension grows. 

When using 250 random directions, the computation time of ROBPCA is comparable to our deterministic algorithms for the different sample sizes and subspace dimensions considered. Increasing the number of directions considered to compute the measure of outlyingness makes  ROBPCA much slower (results not shown here), especially for the case n=5000 and q=5. Thus, computation time of ROBPCA critically depends on the number of directions considered. These results indicate that the algorithms with deterministic starting values do not only show good performance but are also computationally attractive. They make it possible to compute robust subspace estimates for a large scale of problems.




§ PERFORMANCE FOR HIGH DIMENSIONAL DATA


In Sections <ref> and <ref> we have shown that DsubS and DsubLTS show good performance for low-dimensional data and remain computationally attractive for high-dimensional data. We now investigate whether their good performance also carries over to the high-dimensional setting. 
To this end we consider the model of the previous section with eigenvalues according to (<ref>) for  q=2 or q=5, and ϵ=20% of contamination generated as in Section <ref>. We let k range between 0 and 20 as before and generated data of size n=1000 in dimension p=1000 or p=10000. 

 For p=1000 we did not include the S-L and S-M algorithms in the comparison because their performance was similar to the RsubS and RsubLTS algorithms, but computationally more demanding. Moreover, we did not consider the RsubS and RsubLTS algorithms in the case p=10000 because these algorithms require too much computation time in these settings. 
SPC implemented with the  function  also required excessive computation time for the case p=10000. Therefore, we wrote an alternative SPC  function that computes classical PCA on the projected data using Algorithm <ref>. 

Figures <ref> and <ref> show the mean relative prediction errors, averaged over M=200 samples, when p=1000 and p=10000, respectively. Relative prediction errors again stabilized from k=10 onwards, so we only present results up to k=10. From Figure <ref> we can see that in higher dimensions DsubS and DsubLTS behave similar as for low-dimensional data and clearly perform better than SPC and PPLTS. Figure <ref> shows that DsubS and DsubLTS keep their excellent performance when the dimension is increased further to p=10000. The performance of ROBPCA is mostly similar to our deterministic algorithms in this case, but it has slightly more difficulty to estimate the best q=5 dimensional subspace when the outliers are at a small distance (k=1.5) from the regular data. Similar conclusions were obtained for sample size n=100 so these results are not shown.




§ REAL DATA EXAMPLE

We illustrate our deterministic algorithms on real image data that have been vectorized to form a high-dimensional dataset with p>>n. We use the Extended Yale Face Database B <cit.> which contains aligned grayscale face images of 38 subjects under the same frontal pose and 64 different illumination conditions. This database contains cropped face images which have resolution 192 × 168 (=32256) pixels each. For the analysis we randomly sampled 11 subjects and selected for each subject the 6 images with the highest illumination contrasts, so that the face characteristics are clearly identifiable. In particular, we used the images with the following codes for light conditions: "P00A+000E+00.pgm", "P00A+010E-20.pgm", "P00A-010E-20.pgm", "P00A+000E-20.pgm",
"P00A-005E-10.pgm" and "P00A+005E-10.pgm". Similarly to 
<cit.>, next to these 66 face images we also sampled 9 non-face images from the Caltech 101 database <cit.> which constitute potential outliers. The Caltech 101 database contains 9144 color and grayscale images from 102 object categories that includes vehicles, plants, animals and cartoon characters. We only considered grayscale images with a height of at least 192 pixels and a width of at least 168 pixels. We randomly sampled 9 object images from this subset and cropped the sampled images to 192 pixel height and 168 pixel width when necessary. Combining the vectorized face and object images resulted in 75 observations in p=192× 168= 32256 dimensions. The last 9 observations (rows 67-75 in the data matrix) correspond to the object images. Hence, 12% of the observations in the dataset are contamination. Figure <ref> displays a random subset of the face images while Figure <ref> displays the 9 object images in the data matrix. The goal of our analysis is to recognize the faces and thus to flag the object images as outliers. We compare the results obtained by DsubLTS and DsubS to those of SPC, PPLTS and ROBPCA. We estimate the q=2 dimensional subspace which explains about 80% of the robustly estimated total variability. 




To identify outliers, we use the diagnostic plot introduced by <cit.>. It plots the (robust) orthogonal distances between the observations and the estimated subspace versus (robust) score distances of the projected observations in the subspace with respect to their center. To identify outlying observations cutoff values for both the robust orthogonal distances and the robust score distances were proposed by <cit.>. To compute the robust score distances robust estimates for the variances according to the basis directions within the estimated subspace are required. 
Since our algorithms do not yield such estimates, we estimate these variances robustly using univariate LTS or M-scales of the scores corresponding to the estimated basis directions. Also SPC does not yield such variance estimates, so we computed the MCD estimator on the scores to compute robust
distances in the PCA subspace, similarly to <cit.>. Note that as in Section <ref>, we used our SPC implementation in  based on Algorithm <ref> for these high-dimensional data. 

Figure <ref> shows the diagnostic plots of the solutions obtained by DsubLTS, DsubS, SPC, PPLTS and ROBPCA. All robust methods are able to flag the 9 object images as outliers with respect to the estimated 2-dimensional subspace. Object images are flagged as either orthogonal outliers or bad leverage points by the different methods. Note that SPC also detects four faces as outliers, but these false positives are merely borderline cases.

This example illustrates that also in very high-dimensional settings our deterministic algorithms are able to robustly estimate a low-dimensional subspace that best represents the regular observations and allows to identify outliers. 








§ CONCLUSIONS
 

We proposed new algorithms for the robust subspace estimation methods proposed by <cit.>. These algorithms directly estimate principal directions of the subspace, which makes them more suitable for high-dimensional problems. For the starting values of the algorithm we considered random orthogonal matrices, as well as five deterministic starting values.  
These well-chosen deterministic starting values can be computed fast in high-dimensional settings because they avoid the need to calculate a high-dimensional scatter matrix and only use simple manipulations of the data. 

Our experiments show that the deterministic algorithms yield results that are as good as or better than the results for the algorithms with random starting values, while having a much lower computation time.  
Moreover, the excellent performance of our deterministic algorithms carries over to high-dimensional settings. Maronna's algorithms and our algorithms with random orthogonal matrices attempt to find the global minimum of their objective function by using a sufficient number of starting values. On the other hand, 
our algorithms with deterministic starting values aim to find this minimum by using a few well-chosen robust but rough initial estimates of the subspace. While this may decrease the probability of obtaining the global minimum, it does often lead to a robust local minimum. 
This explains why these deterministic algorithms show a better performance. 

For the deterministic algorithms orthogonal equivariance is not guaranteed anymore. However, our experiments have shown that in most cases there is not a serious loss of equivariance. 
Note that the computation time of our deterministic algorithms can be reduced further on multi-core machines by calculating each of the five starting solutions in parallel on different cores. Finally, the subspace S-estimator usually shows a better compromise between robustness and efficiency than the subspace LTS-estimator.
Implementations of our algorithms in R <cit.> are available from the website <http://wis.kuleuven.be/stat/robust/software>.




§ ACKNOWLEDGMENTS

The research by Van Aelst was supported by the Internal Funds KU Leuven under Grant C16-15-068  and 
COST Action IC1408 CRoNoS. Their support is gratefully acknowledged.






