%%%%%% March 31, 2004
%%%%% file name: higuchi04.tex
\documentclass[twoside,11pt]{article}

\usepackage{jmlr2e}

\newcommand{\bm}[1]{\mbox{\boldmath $#1$}}
\newcommand{\argmin}{\mathop{\rm argmin}}
\newcommand{\argmax}{\mathop{\rm argmax}}
\newcommand{\diag}{\mathop{\rm diag}}
\newcommand{\half}{\frac{1}{2}}

\def\labelenumi{(\theenumi)}
\def\theenumi{\roman{enumi}}

\jmlrheading{5}{2004}{453--471}{4/03; Revised 10/03}{5/04}{Isao Higuchi and Shinto Eguchi}

\ShortHeadings{Robust PCA with Adaptive Selection}{Higuchi and Eguchi}
\firstpageno{453}

\begin{document}

\title{Robust Principal Component Analysis with Adaptive Selection for Tuning Parameters}

\author{Isao Higuchi \email higuchi@amath.hiroshima-u.ac.jp \\
       \addr Department of Applied Mathematics\\
       Hiroshima University\\
       1-4-1, Kagamiyama\\
       Higashi-Hiroshima 739-8527, Japan
       \AND
       \name Shinto Eguchi \email eguchi@ism.ac.jp \\
       \addr Institute of Statistical Mathematics and \\
       Graduate University of Advanced Studies \\
       4-6-7, Minami-Azabu, Minato-ku\\
       Tokyo 106-8569, Japan}

\editor{David Madigan}

\maketitle

\begin{abstract}%
The present paper discusses robustness against outliers in a principal component analysis (PCA).
We propose a class of procedures for PCA based on the minimum psi principle, which unifies various approaches, including the classical procedure and recently proposed procedures.
The reweighted matrix algorithm for off-line data and the gradient algorithm for on-line data are both investigated with respect to robustness.
The reweighted matrix algorithm is shown to satisfy a desirable property with local convergence, and the on-line gradient algorithm is shown to satisfy an asymptotical stability of convergence.
Some procedures in the class involve tuning parameters, which control sensitivity to outliers.
We propose a shape-adaptive selection rule for tuning parameters using K-fold cross validation.
\end{abstract}

\begin{keywords}
  K-fold cross validation, on-line algorithm, reweighted matrix algorithm,
  influence function, data contamination
\end{keywords}

\section{Introduction}

In both neural network and statistical studies, PCA is one of the most fundamental tools of dimensionality reduction for extracting effective features from high-dimensional vectors of input data.
See Croux and Haesbroeck (2000) and De la Torre and Black (2001) for recent
discussions.
PCA is implemented by projecting input data onto the most informative subspace of lower dimension so that the hidden structure behind the input data may be clarified.
The procedure of detecting principal components from input data in an on-line manner is related to the mechanism of a single neuron by the Hebbian adaptation rule, for which the learning theory has been discussed (Amari, 1977; Oja, 1982).

One of the frequently occurring difficulties in PCA is that a few outliers give disturbance in finding the effective features in a bulk of input vectors.
The usual PCA satisfies the statistical optimality only under the assumption of a Gaussian distribution for all of the input data.
A small departure from the assumption produces a gross error in the performance of the principal component vector or subspace in the PCA. 
This motivates our study of robust PCA procedures.
In what follows, we discuss an $\varepsilon $-contamination model for a data distribution $F_\varepsilon $ defined by
\begin{equation}\label{Out1}
   F_\varepsilon=(1-\varepsilon)N(\bm{\mu},V)+\varepsilon H,
\end{equation}
where $N(\bm{\mu},V)$ is a Gaussian distribution with mean vector $\bm{\mu}$ and
covariance matrix $V$, and $H$ is a distribution of possible outliers.
In this modeling, it is assumed that $\varepsilon$, or the probability of outliers, is small and that the distribution $H$ is unspecified but qualitatively different from the supposed distribution $N(\bm{\mu},V)$.
Thus, the model (\ref{Out1}) lies in a kind of $\varepsilon $-neighborhood surrounding $N(\bm{\mu},V)$ taking into account all of the possible distributions $H$.
If $H$ has the mean vector $\bm{\mu}_H$ and the covariance matrix $V_H$, then the data distribution $F_\varepsilon$ has the covariance matrix
\begin{equation}\label{var}
   V_\varepsilon=(1-\varepsilon)V+\varepsilon V_H + \varepsilon
   (1-\varepsilon)(\bm{\mu} - \bm{\mu}_H)(\bm{\mu} - \bm{\mu}_H)^T.
\end{equation}
Thus, the classical procedure works properly as long as $V_H \simeq V$ and $\bm{\mu}_H \simeq \bm{\mu}$, because the procedure basically searches for the dominant eigenvectors of $V_\varepsilon$ by learning from input data generated by the assumed distribution $F_\varepsilon$.
However, even if the probability $\varepsilon$ is quite small, the classical procedure often breaks down when input data have a distribution $F_\varepsilon$, such that $\bm{\mu}_H$ or $V_H$ is far from the assumed $\bm{\mu}$ or $V$, as observed from (\ref{var}).

A variety of outlier distributions $H$ have infinite dimensionality, and the simplest candidate is a point-mass distribution $\delta_{\xi}(\bm{x})$ degenerated at $\bm{x}=\bm{\xi}$.
This choice corresponds to a situation in which the outliers occur deterministically in a singleton $\bm{\xi}$ with probability $\varepsilon $.
The influence function of a procedure in the PCA is defined as the derivative at $\varepsilon = 0$  of the procedure under $F_\varepsilon$ in (\ref{Out1}) with  $H=\delta_{\xi}$.
This concept will be more explicitly explored in a subsequent discussion. 
See Higuchi and Eguchi (1998) for the case of the PCA, and see Hampel (1974) and Hampel {\it et al.}~(1986) for the general case.

We discuss a class of principal component analyzers defined using generic functions which contain tuning parameters.  
For example if we adopt a log-sigmoidal function as a generic function,
the tuning parameters are the inverse temperature and saturation value parameters,
as will be discussed in detail.
In general the tuning parameter set makes a delicate trade 
between loss of information and degree of insensitivity to outliers.
The main objective in the present paper is to provide a reasonable selection of
tuning parameters of principal component analyzers.
The basic idea is to craft a loss function that reflects as appropriate trade off 
between loss of information and robustness to outliers.
We introduce K-fold cross validation for estimating the expected
loss based on a given data set.
As a result we build a method of data-adaptive selection of tuning parameters.
In a simulation study, we examine the performance of the adaptive selection under
three types of outlier distributions $H $ displaying deterministic, structural and distributional contaminations based on (\ref{var}).
The three types of outliers are simulated in a numerical experiment, and we test the performance in a few cases of principal component analyzers.
We provide an S implementation of the basic robust PCA at {\tt http://home.hiroshima-u.ac.jp/oxbow/RobustPCA/}.

The present paper is organized as follows.
Section 2 introduces a class of procedures in PCA derived by the minimum  psi method.
In Section 3 we discuss the robustness of the procedure in the class, and in Section 4 we present an adaptive method of selecting tuning parameters.
Finally, in Section 5 we provide the results of a simple simulation study to validate our theoretical discussion in previous sections and tests numerical behavior in three types of departures from the Gaussian distributional assumption.

\section{A Class of Principal Component Vectors}

In this section we propose a class of principal component vectors.
In general, PCA aims to extract the most informative $k$-dimensional output vector $\bm{y}$ from an input vector $\bm{x}$ of $p$-dimension.
This is achieved by learning the matrix $\Gamma $ which connects $\bm{x}$ to $\bm{y} = \Gamma^{\rm T} (\bm{x}-\bm{\mu})$ based on input data $\{\bm{x}_t;t=1,2,\cdots\}$, where $\bm{\mu}$ is a vector of center of the input data and $\Gamma $ is a $p\times k$ orthonormal matrix, or $\Gamma^{\rm T}\Gamma=I $ (the $k$-identity matrix).
In neural networks, $\Gamma $ is interpreted as the matrix of coefficients connecting $p$ neurons to $k$ neurons, where a learning process works by renewing $\Gamma $ according to a batch of inputs in an off-line manner or sequential input vectors in an on-line manner (Oja, 1982, 1989 and \S 8 in Haykin, 1999).
By combining these approaches, we propose a certain class of procedures for PCA.

We present a concise review of the classical PCA for detecting the principal $k$-subspace.
Let 
\begin{equation}\label{zxmG}
   z(\bm{x},\bm{\mu},\Gamma)=\frac {1}{2}
   \{\|\bm{x}-\bm{\mu}\|^2-\|\Gamma^{\rm T}(\bm{x}-\bm{\mu})\|^2\}
\end{equation}
be half the square of the residual distance of $\bm{x}-\bm{\mu}$ from the subspace spanned by the columns of $\Gamma $.
We note that $z(\bm{x},\bm{\mu},\Gamma)=\frac{1}{2}\min_{\beta \in R^k} \|\bm{x}-\bm{\mu}-\Gamma \bm{\beta}\|^2$.
See Hotelling (1933) for the original derivation.
The classical PCA is simply given by minimizing
\[
   \frac {1}{n}\sum_{t=1}^n z(\bm{x}_t,\bm{\mu},\Gamma)
\]
with respect to $\bm{\mu}$ and $\Gamma $, which reduces to solving $k$ dominant eigenvectors of the sample covariance matrix
\begin{equation}\label{S}
   S=\frac {1}{n}\sum_{t=1}^n
   (\bm{x}_t-\bar{\bm{\mu}})(\bm{x}_t-\bar{\bm{\mu}})^{\rm T},
\end{equation}
where the centralized vector $\bar{\bm{\mu}}$ is given by $\sum_{t=1}^n \bm{x}_t/n$.
Thus, we obtain a solution $\Gamma $ by stacking the $k$ dominant eigenvectors of $S$, which we write in the form
\[
   \Gamma={\rm eigen}(S).
\]


We propose a variant of this classical procedure for PCA obtained by minimizing an objective function
\begin{equation}\label{E1}
   {\cal E}(\bm{\mu},\Gamma)=\frac{1}{n} \sum_{t=1}^n \Psi(z(\bm{x}_t,\bm{\mu},\Gamma)), 
\end{equation}
where $\Psi(z)$ is assumed to be a monotonic increasing function of $z>0$.
Various $\Psi$ yield various procedures for PCA.
As typical examples, the identity function $\Psi_0(z)=z$ reduces to 
the classical PCA and  
\begin{equation}\label{psi1}
   \Psi_1(z)=\log \frac{1}{1+\exp\{-\beta(z-\eta)\}}
\end{equation}
defines Xu and Yuille's self-organizing rule, where $\beta $ and $\eta $ are tuning parameters, referred to as the inverse temperature and saturation value, respectively (Xu and Yuille, 1995).
Another possible function is
\begin{equation}\label{psi2}
   \Psi_2(z)=\frac{1- \exp(-\beta z)}{\beta}.
\end{equation}
In general, $\Psi$ is interpreted as the generic function which gives the total function ${\cal E}$, and we refer to the minimization of ${\cal E}$ in (\ref{E1}) as the ``minimum psi principle generated by $\Psi$".





Based on an argument similar to that of the classical PCA, we observe that the minimizer $(\tilde{\bm{\mu}},\tilde{\Gamma})$ of ${\cal E}(\bm{\mu},\Gamma)$ satisfies 
the stationary equations
\begin{eqnarray}\label{m1}
   \tilde{\bm{\mu}} &=& \sum_{t=1}^n p_t(\tilde{\bm{\mu}},\tilde{\Gamma})
   \bm{x}_t, \\ \label{G1}
   \tilde{\Gamma} &=& {\rm eigen}(S(\tilde{\bm{\mu}},\tilde{\Gamma})),
\end{eqnarray}
where
\[
   p_t(\bm{\mu},\Gamma)=\frac{\psi(z(\bm{x}_t,\bm{\mu},\Gamma))}
   {\sum_{s=1}^n \psi(z(\bm{x}_s,\bm{\mu},\Gamma))},
\]
\begin{equation}\label{SmG}
   S(\bm{\mu},\Gamma)=\sum_{t=1}^n
   p_t(\bm{\mu},\Gamma)(\bm{x}_t-\bm{\mu})(\bm{x}_t-\bm{\mu})^{\rm T},
\end{equation}
with $\psi(z)=(\partial/\partial z)\Psi(z)$.
Thus, the equilibrium point $(\tilde{\bm{\mu}},\tilde{\Gamma}) $ is expressed by
 the weighted mean and the covariance matrix, where the weight function $p_t$ depends upon $\tilde{\bm{\mu}}$ and $\tilde{\Gamma}$, except for the case of $\psi(z)=1$, which yields the classical procedure.
In effect, (\ref{m1}) is determined only up to the addition of a vector in the subspace associated with $\tilde{\Gamma} $.
Thus $\bm{\mu}^* =\tilde{\bm{\mu}} +\bm{\gamma}$ is also a solution of (\ref{m1}) if $\bm{\gamma} \in {\rm Im}(\tilde{\Gamma}) \equiv \{  \tilde{\Gamma}\bm{c} | \bm{c} \in {\bm R}^k\}$, because $\bm{\mu}^* +{\rm Im}(\tilde{\Gamma})=\tilde{\bm{\mu}}+{\rm Im}(\tilde{\Gamma})$.
However, (\ref{G1}) is independent of the choice of possible solutions, so we adopt the expression (\ref{m1}) for convenience.

In PCA research, centralization of input data by a vector other than a sample mean has not been considered. In the present paper, we investigate the problem of centralization and  explore the usefulness of a method using a vector such as (\ref{m1}).
In conventional robust statistics, the estimation of location has been
studied extensively. (See, for example, Huber, 1981.)
Our estimator $\tilde{\bm{\mu}}$ can be viewed as one of several variants for robust estimation.
However, our main objective concerning $\tilde{\bm{\mu}}$ is not the location estimation itself, but rather the data centralization for the extraction of principal components.
Thus, $\tilde{\bm{\mu}}$ is naturally linked to $\tilde{\Gamma}$ in the optimization of ${\cal E}(\bm{\mu},\Gamma)$.

For a batch of data $\{\bm{x}_t:1\le t\le n\}$, we propose a fixed-point algorithm.
See Hyvarinen and Oja (1997) for the related discussion on a fixed-point algorithm for ICA.
This algorithm alternates two steps associated with the stationary equations (\ref{m1}) and (\ref{G1}) in the following:

Step 1: Given $(\bm{\mu}_1,\Gamma_1)$, calculate
\[
   p_t^{(1)}=\frac {\psi(z(\bm{x}_t,\bm{\mu}_1,\Gamma_1))}
   {\sum_{s=1}^n \psi(z(\bm{x}_s,\bm{\mu}_1,\Gamma_1))}.
\]

Step 2: Using the estimated $\{p_t^{(1)}\}$ in step 1, perform the same task as in the classical PCA:
\begin{eqnarray}
%  \label{mu2}
   \bm{\mu}_2 &=& \sum_{t=1}^n p_t^{(1)} \bm{x}_t \quad {\rm and} \nonumber \\ 
   \label{G2}
   \Gamma_2 &=& {\rm eigen}(S^{(1)}),
\end{eqnarray}
where $S^{(1)}$ is a weighted matrix defined by $\sum p_t^{(1)}(\bm{x}_t-\bm{\mu}_2)(\bm{x}_t-\bm{\mu}_2)^{\rm T}$.
In this way, the algorithm alternates between two steps, and we refer to this as the reweighted matrix (RM) algorithm.

Assuming hereafter that the generic function $\Psi(z)$ is strictly concave in $z$, we have
\[
   {\cal E}(\bm{\mu}_2,\Gamma_2)-{\cal E}(\bm{\mu}_1,\Gamma_1) <
   \frac {\sum_{t=1}^n \psi(z(\bm{x}_t,\bm{\mu}_1,\Gamma_1))}{n}
   \{ \sum_{t=1}^n p_t^{(1)}z(\bm{x}_t,\bm{\mu}_2,\Gamma_2)
   -\sum_{t=1}^n p_t^{(1)}z(\bm{x}_t,\bm{\mu}_1,\Gamma_1) \}
\]
because, by assumption, $\Psi(z_2)-\Psi(z_1)<\psi(z_1)(z_2-z_1)$ for
$z_1 < z_2$.
In step 2, the procedure is equivalent to minimization of
\[
   \sum p_t^{(1)} z(\bm{x}_t,\bm{\mu},\Gamma)
\]
with respect to $(\bm{\mu},\Gamma)$.
Therefore, we conclude that the RM algorithm generated by $\{(\bm{\mu}_j,\Gamma_j):j\ge 1\}$ is responsible for the strict decrease of the objective function
\[
   {\cal E}(\bm{\mu}_1,\Gamma_1)>\cdots>{\cal E}(\bm{\mu}_j,\Gamma_j)>\cdots.
\]
This desirable property is mathematically the same as that of the EM algorithm.
The possible region of $(\bm{\mu},\Gamma)$ that the algorithm works in is ${\cal X} \times {\cal O}_{p,k}$, where $\cal X$ is the convex hull of data and ${\cal O}_{p,k}$ is the space of $p \times k$ orthonormal matrices.
We can easily check a condition for the convergence to the solution of the equations such that a set
\[
   \{(\bm{\mu},\Gamma)\in {\cal X} \times 
   {\cal O}_{p,k}:{\cal E}(\bm{\mu},\Gamma)\le c\}
\]
is compact for any fixed $c \le {\cal E}(\bm{\mu}_1,\Gamma_1)$.
This is referred to as the regularity condition of Wu (1983), found in \S 3.4.2 of McLachlan and Krishnan (1997), and this condition implies that the sequence $\{(\bm{\mu}_j,\Gamma_j):j \ge 1\}$ is convergent to a set of solutions of equations (\ref{m1}) and (\ref{G1}).
However, even when this regularity condition holds, the computational complexity with high dimensional data may be prohibitive since each iteration requires the solution of (\ref{G2}).

\subsection{Stability of On-line Gradient Algorithm}

We next discuss the on-line gradient algorithm.
The gradient vector of the objective function is the sum of
\[
   G(\bm{\mu},\Gamma)(\bm{x}_t)
   =\psi(z(\bm{x}_t,\bm{\mu},\Gamma))G_1(\bm{\mu},\Gamma)(\bm{x}_t)
\]
over $t=1,2,\cdots$, where
\[
   G_1(\bm{\mu},\Gamma)(\bm{x})=\left[
    \begin{array}{c}
    \bm{x}-\bm{\mu} \\
    (\bm{x}-\bm{\mu})(\bm{x}-\bm{\mu})^{\rm T}\Gamma
    -\Gamma\ {\cal LT}[\bm{y}\bm{y}^{\rm T}]
    \end{array}
   \right]
\]
with $\bm{y}=\Gamma^{\rm T}(\bm{x}-\bm{\mu})$, where the operator ${\cal LT}[\cdot]$ sets all of the elements above the diagonal of its matrix argument to zero.
Hence, the on-line gradient algorithm is given by
\begin{equation}\label{mG}
   \left[
    \begin{array}{c}
       \bm{\mu}_{t+1} \\
       \Gamma_{t+1}
    \end{array}
   \right]=
   \left[
    \begin{array}{c}
       \bm{\mu}_t \\
       \Gamma_t
    \end{array}
   \right]+r_t G(\bm{\mu}_t,\Gamma_t)(\bm{x}_t)
\end{equation}
for $t=1,2,\cdots$ with a learning rate $r_t$.
If we apply the classical procedure this algorithm reduces to the Oja algorithm (1982).
See \S 8 in Haykin (1999) for the related algorithmic developments in PCA.
The gradient algorithm (\ref{mG}) is different from the Oja algorithm only with respect to the factor $\psi(z)$, which depends on the $t$-step $(\bm{\mu}_t,\Gamma_t)$ and the $t$-th example $\bm{x}_t$ through $z=z(\bm{x}_t, \bm{\mu}_t, \Gamma_t)$ defined in (\ref{zxmG}).
In the classical procedure, the $\bm{\mu}$-part of the algorithm (\ref{mG}) reduces to the usual centralization and has no connection to $\Gamma $.
However, in the case of a non-constant weight function $\psi(z)$, the $\bm{\mu}$-part is essentially connected not only to $\bm{\mu}$ itself, but also to $\Gamma $ through the $z$-variable.
For the case of non-constant $\psi(z)$, we confine ourselves to the Xu-Yuille rule, which is generated by $\psi_1(z)=\beta/(1+e^{\beta(z-\eta)})$.
Xu and Yuille implemented the on-line algorithm in the same fashion as (\ref{mG}) for the $\Gamma $-part, but the centralizing mean $\bar{\bm{\mu}}$ was used for the $\bm{\mu}$-part.
We will make a simple comparison between the two methods for the $\bm{\mu}$-part in a subsequent discussion.

The on-line gradient algorithm does not satisfy the property of uniform decrease of the objective function possessed by the reweighted algorithm as shown above.
We first discuss the asymptotic convergence of (\ref{mG}) for the case of $k=1$, $\Gamma = \bm{\gamma}$.
See \S 8.4 in Haykin (1999) for the proof for the classical procedure.
The on-line gradient algorithm (\ref{mG}) is a special case of the generic stochastic approximation algorithm 
\begin{equation}\label{cla}
   \left[ \begin{array}{c}
      \bm{\mu}_{t+1} \\
      \bm{\gamma}_{t+1}
   \end{array}\right]
   =
   \left[ \begin{array}{c}
      \bm{\mu}_t \\
      \bm{\gamma}_{t}
   \end{array}\right]
   +r_t \bm{h}(\bm{\gamma}_t,\bm{\mu}_t,\bm{x}_t),
\end{equation}
where
\[
   \bm{h}(\bm{\gamma},\bm{\mu},\bm{x})=\psi(z(\bm{x},\bm{\mu},\bm{\gamma}))
   \left[\begin{array}{c}
      \bm{x}-\bm{\mu} \\
      (\bm{x}-\bm{\mu})(\bm{x}-\bm{\mu})^{\rm T}\bm{\gamma}
      -\{\bm{\gamma}^{\rm T}(\bm{x}-\bm{\mu})(\bm{x}-\bm{\mu})^{\rm T}
      \bm{\gamma}\}\bm{\gamma}
   \end{array}\right].
\]
If $\psi(z)\equiv 1$, then (\ref{cla}) leads to the classical PCA.
It is assumed that $\psi$ is a finite function, so the convergence is proved using an argument similar to that used in the case of classical PCA.


Take the expectation of $h(\bm{\gamma}_t,\bm{\mu}_t,\bm{x}_t)$ over $\bm{x}$, and then in the limit we have
\begin{eqnarray}
   \bar{\bm{h}}(\bm{\gamma}_\infty,\bm{\mu}_\infty)&=&\lim_{t \to \infty} 
   E[\bm{h}(\bm{\gamma}_t,\bm{\mu}_t,\bm{x}_t)] \nonumber \\
   &=&
   \left[
        \begin{array}{c}
        \bm{m}(\bm{\gamma}_\infty,\bm{\mu}_\infty)
        -\kappa(\bm{\gamma}_\infty,\bm{\mu}_\infty)\bm{\mu}_\infty \\
        R(\bm{\gamma}_\infty,\bm{\mu}_\infty)\bm{\gamma}_\infty
        -\{\bm{\gamma}_\infty^{\rm T}R(\bm{\gamma}_\infty,\bm{\mu}_\infty)
        \bm{\gamma}_\infty\}\bm{\gamma}_\infty
        \end{array}
   \right], \nonumber
\end{eqnarray}
where 
\[
\kappa(\bm{\gamma},\bm{\mu})=E\{\psi(z(\bm{x},\bm{\mu},\bm{\gamma}))\},
\ \ \ 
\bm{m}(\bm{\gamma},\bm{\mu})=E\{\psi(z(\bm{x},\bm{\mu},\bm{\gamma}))\bm{x}\}
\]
and
\[
R(\bm{\gamma},\bm{\mu})=E[\psi(z(\bm{x},\bm{\mu},\bm{\gamma})) (\bm{x}-\bm{\mu})
(\bm{x}-\bm{\mu})^{\rm T}].
\]
Thus, our differential equation is
\begin{eqnarray}
  \frac{d}{dt}
  \left[
        \begin{array}{c}
           \bm{\mu}_t \\
           \bm{\gamma}_t
        \end{array}
  \right]
  &=& \bar{\bm{h}}(\bm{\gamma}_t) \nonumber \\
  &=&
  \left[
        \begin{array}{c}
        \bm{m}(\bm{\gamma}_t,\bm{\mu}_t)
        -\kappa(\bm{\gamma}_t,\bm{\mu}_t)\bm{\mu}_t \\
         R(\bm{\gamma}_t,\bm{\mu}_\infty)\bm{\gamma}_t
        -\{\bm{\gamma}_t^{\rm T}R(\bm{\gamma}_t,\bm{\mu}_\infty)
        \bm{\gamma}_t\}\bm{\gamma}_t
        \end{array}
  \right].
  \label{eq2}
\end{eqnarray}
In the $\bm{\mu}$-part, we observe that $\bm{\mu}_t$  behaves asymptotically  as $e^{-\kappa_\infty t}\bm{a}+\bm{m}_\infty/\kappa_\infty$, which implies that $\bm{m}_\infty/\kappa_\infty$ is the stable-point, where $\kappa_\infty=\kappa(\bm{\gamma}_\infty,\bm{\mu}_\infty)$ and $\bm{m}_\infty=\bm{m}(\bm{\gamma}_\infty,\bm{\mu}_\infty)$.
Therefore, we consider only 
\[
   \frac{d}{dt}\bm{\gamma}_t= R(\bm{\gamma}_t,\bm{\mu}_\infty)\bm{\gamma}_t
   -\{\bm{\gamma}_tR(\bm{\gamma}_t,\bm{\mu}_\infty)\bm{\gamma}_t\}\bm{\gamma}_t.
\]
We expand $\bm{\gamma}_t$ in terms of the set of eigenvectors $\{\bm{q}_k(\infty):k=1,\cdots,p\}$ of $R(\bm{\gamma}_\infty,\bm{\mu}_\infty)$ with the dominant eigenvector $\bm{q}_1(\infty)$ as follows:
\[
   \bm{\gamma}_t=\sum \theta_k(t)\bm{q}_k(\infty),
\]
Let us decompose $\bar{\bm{h}}(\bm{\gamma})$ into $\bar{h}_1(\bm{\gamma},\bm{\gamma}_\infty)+\bar{h}_2(\bm{\gamma},\bm{\gamma}_\infty)$, where
\begin{eqnarray}
   \bar{h}_1(\bm{\gamma},\bm{\gamma}_\infty) &=&
   R(\bm{\gamma}_\infty,\bm{\mu}_\infty)\bm{\gamma}
   -\{\bm{\gamma}^{\rm T}R(\bm{\gamma}_\infty,\bm{\mu}_\infty)
   \bm{\gamma}\}\bm{\gamma} \nonumber \\
   \bar{h}_2(\bm{\gamma},\bm{\gamma}_\infty) &=&
   \{R(\bm{\gamma},\bm{\mu}_\infty)-R(\bm{\gamma}_\infty,\bm{\mu}_\infty)\}\bm{\gamma}\\
   & &-[\bm{\gamma}^{\rm T}\{R(\bm{\gamma},\bm{\mu}_\infty)
   -R(\bm{\gamma}_\infty,\bm{\mu}_\infty)\}\bm{\gamma}]\bm{\gamma}.\nonumber
\end{eqnarray}
Then the equilibrium condition $\bar{h}_1(\bm{\gamma}_\infty,\bm{\gamma}_\infty)=0$ implies that $\bm{\gamma}_\infty$ reduces to one of $k$ eigenvectors of $R(\bm{\gamma}_\infty,\bm{\mu}_\infty)$; $\bar{h}_2(\bm{\gamma}_\infty,\bm{\gamma}_\infty)=0$ holds identically.
The differential equation
\[
   \frac{d}{d t}\bm{\gamma}_t
   =\bar{h}_1(\bm{\gamma}_t,\bm{\gamma}_\infty)
\]
has an asymptotically stable point $\bm{q}_1(\infty)$ through the same discussion established in \S 8.4 in Haykin (1999).
Hence the differential equation (\ref{eq2}) leads to stable convergence to $\bm{q}_1(\infty)$.


Secondly, we observe that the case of $k$ principal component vectors also satisfies the stable convergence, noting that our differential equation is
\[
  \frac{d}{dt}\Gamma_t = R(\Gamma_t,\bm{\mu}_\infty)\Gamma_t
  -{\cal LT}[\Gamma_t^{\rm T}R(\Gamma_t,\bm{\mu}_\infty)\Gamma_t]\Gamma_t,
\]
where
\[
   R(\Gamma,\bm{\mu})=E\{\psi(z(\bm{x},\bm{\mu},\Gamma))
   (\bm{x}-\bm{\mu})(\bm{x}-\bm{\mu})^{\rm T}\}.
\]

In effect, the RM algorithm is applicable to on-line data by solving the eigen problem for batch data with a new observation incorporated in each step.
The computational burden is quite heavy relative to the on-line gradient algorithm, but we will pursue more rapid convergence property in a simulation study.

In the statistical literature another type of PCA methods has been proposed by minimizing
\[
   \frac {1}{n}\sum_{t=1}^n \Psi (d(\bm{x}_t,\bm{\mu},V))
\]
with respect to $(\mu, V)$, where $d$ is Mahalanobis squared distance, that is,
\[
        d(\bm{x}_t,\bm{\mu},V) = \frac{1}{2} (\bm{x}_t-\bm{\mu})^{\rm T} V^{-1}(\bm{x}_t-\bm{\mu}).
\]
See Campbell(1980), Devlin {\it et~al.}~(1981),  Caussinus and Ruiz
(1990), and Croux and Haesbroeck (2000).
The use of the nonlinear generic function $\Psi$ is the same, but the essential difference 
is that our method aims at estimating the
principal component vectors rather than estimating the scatter matrix $V$.
One advantage of our method is that it does not need
all the information of $V$.  In fact only the first $k$ dominant
eigenvalues and the corresponding eigenvectors are needed in the algorithm,
which is easily implemented by the singular-value decomposition algorithm
even if the data set is of high dimension.


\section{Robustness of the Proposed Principal Component Vectors}

Data analysts have frequently found that the classical PCA breaks down in the presence of outliers.
It can happen that a single outlier changes the principal component subspace into the orthogonal complement.
As a result, the PCA fails to capture an important feature of the bulk of the data, which will be observed from
a simple simulation study in Section 5.
In the statistical literature, the robustification of the classical likelihood-based procedures has been discussed and well established: see Huber (1981) for some notions on robustness.
Such contamination is typically expressed by the $\varepsilon $ -contamination model,
\[
   (1-\varepsilon)N(\bm{\mu},V)(\bm{x})+\varepsilon \delta_{\xi}(\bm{x}),
\]
with the point-mass distribution $\delta_{\xi} $ as discussed in the Introduction.
In the expression, $\varepsilon $ is undetectably small; nevertheless, the likelihood procedure based on the density function of $N(\bm{\mu},V) $ may sometimes break down for an extreme vector $\bm{\xi}$.
We next explore which first principal component vector or principal subspace is robust against outliers in our class.
First, we consider the case of the first principal component vector, or $k=1 $.
In Higuchi and Eguchi (1998), the influence function of the Xu and Yuille rule defined by $\Psi_1$ in (\ref{psi1}) is given by
\[
   {\rm IF}_{\Psi_1} (\bm{\xi})=-
   {\psi_1} (z(\bm{\xi},\bm{\mu},\bm{\gamma}_1))
   \tilde{\bm{\gamma}}_1^{\rm T}(\bm{\xi}-\bm{\mu})
   \sum_{j=2}^p \frac{\hat{\lambda}_j\tilde{\bm{\gamma}}_j^{\rm T}
   (\bm{\xi}-\bm{\mu})}
   {\tilde{\lambda}_j(\hat{\lambda}_j-\hat{\lambda}_1)}\tilde{\bm{\gamma}}_j,
\]
where $(\hat{\lambda}_j,\hat{\bm{\gamma}}_j)$  is the pair of the 
$j$-th dominant eigenvalue and its associated eigenvector of $S$ defined at (\ref{S}) and $(\tilde{\lambda}_j,\tilde{\bm{\gamma}}_j)$ is that of $S(\bm{\mu},\Gamma)$ defined at (\ref{SmG}),  and 
\begin{equation}\label{XY}
  {\psi_1} (z)= \frac{\partial \Psi_1(z)}
  {\partial z}= \frac{\beta}{1+\exp\{\beta(z-\eta)\}}.
\end{equation}

The influence function assesses the effect on the principal component subspace of the contamination of the data $\{\bm{x}_t|1 \le t \le n\}$ by the outlier $\bm{\xi}$.

Secondly, we discuss a general case of $k$ principal component vectors $\Gamma=(\bm{\gamma}_1,\ldots,\bm{\gamma}_k)$.
We consider the matrix
\[
   P=\Gamma\Gamma^{\rm T}.
\]
The matrix $P$ is the projection operator onto the subspace spanned by the eigenvectors 
${\bm{\gamma}}_I$, see Tanaka (1988).
Our estimator $\tilde{P}=\tilde{\Gamma}\tilde{\Gamma}^{\rm T}$ has the influence function
\begin{equation}\label{IFk}
   {\rm IF}_{\Psi_1} (\bm{\xi};\tilde{P})=-
   {\psi_1} (z(\bm{\xi},\tilde{\bm{\mu}},\tilde{\Gamma}))
      \sum \frac{\hat{\lambda}_j\tilde{\bm{\gamma}_I}^{\rm T}
      (\bm{\xi}-\tilde{\bm{\mu}})(\bm{\xi}-\tilde{\bm{\mu}})^{\rm T}
      \tilde{\bm{\gamma}_j}}
      {\tilde{\lambda}_j(\hat{\lambda}_j-\hat{\lambda}_I)}
      (\tilde{\bm{\gamma}}_I\tilde{\bm{\gamma}}_j^{\rm T}+
      \tilde{\bm{\gamma}}_j\tilde{\bm{\gamma}}_I^{\rm T}),
\end{equation}
where the summation is taken over $\{(I,j): I=1,\cdots,k, j=1,\cdots,p, I \neq j\}$, and we will also use this summation convention in a subsequent discussion.

The formula is valid for the minimum psi principle for a general $\Psi $, see Kamiya and Eguchi (2001) for a detailed discussion, as well as a discussion of relative efficiency under a Gaussian distribution.
The influence function, as a function of $\bm{\xi}$, assesses the smoothness of the principal component subspace around the supposed distribution $N(\bm{\mu},V)$.
The boundedness of the influence function in $\bm{\xi} $ qualitatively guarantees robustness for the target principal component subspace.

In PCA, $\bm{\mu}$ needs to be estimated in order to centralize the data before extracting the principal component subspace.
This is usually estimated as $\sum \bm{x}_t/n $, which can be expressed in the form of a functional $\int \bm{x} dF(\bm{x}) $ with $F = F_n$ (empirical distribution).
This usual estimation is Fisher consistent since $\int \bm{x} dF(\bm{x}) = \bm{\mu}$ when $F=N(\bm{\mu},V)$, but is quite sensitive to the outlier $\bm{\xi}$,
since the functional evaluated at $F$ is $(1-\varepsilon)\bm{\mu}+\varepsilon \bm{\xi}$, and so the influence function is $ \bm{\xi}-\bm{\mu}$.
In contrast, we will show that our PCA procedure automatically leads to a robust centralization of data.
We typically observe an unbounded case for the usual principal component subspace, and a bounded case for Xu and Yuille's principal component subspace.

The boundedness of ${\rm IF}_\Psi$ conditional on
\begin{equation}\label{d2}
  \|\tilde{\Gamma}^{\rm T}(\bm{\xi}-\tilde{\bm{\mu}})\|^2 \le d^2
\end{equation}
with a positive constant $d$ guarantees robustness against any outliers $\bm{\xi}$ satisfying (\ref{d2}) in the contamination model, cf. Higuchi and Eguchi (1998) for the justification for this robustness.
Using the formula (\ref{IFk}) with general $\Psi$, we find a sufficient condition for the robustness 
\begin{equation}\label{scon}
   \sup_{z>0} \sqrt{z} \psi(z) < \infty
\end{equation}
with $\psi(z)=\partial \Psi(z)/\partial z$.
The proof is given as follows.
First, we obtain
\begin{eqnarray}
   \|{\rm IF}_\Psi(\bm{\xi};\tilde{P})\| &\le& 
   \psi(z(\bm{\xi},\tilde{\bm{\mu}},\tilde{\Gamma}))
    \sum \left|\frac{\hat{\lambda}_j
   \tilde{\bm{\gamma}_I}^{\rm T}(\bm{\xi}-\tilde{\bm{\mu}})
   (\bm{\xi}-\tilde{\bm{\mu}})^{\rm T}\tilde{\bm{\gamma}_j}}
   {\tilde{\lambda}_j(\hat{\lambda}_j-\hat{\lambda}_I)}\right|
   \|\tilde{\bm{\gamma}}_I\tilde{\bm{\gamma}}_j^{\rm T}+
   \tilde{\bm{\gamma}}_j\tilde{\bm{\gamma}}_I^{\rm T}\| \nonumber
   \\
   &\le & d \sum \left|\frac{\hat{\lambda}_j
   \|\tilde{\bm{\gamma}}_I\tilde{\bm{\gamma}}_j^{\rm T}+
   \tilde{\bm{\gamma}}_j\tilde{\bm{\gamma}}_I^{\rm T}\|}
   {\tilde{\lambda}_j(\hat{\lambda}_j-\hat{\lambda}_I)}\right|
   |(\bm{\xi}-\tilde{\bm{\mu}})^{\rm T}\tilde{\bm{\gamma}_j}|
   \psi(z(\bm{\xi},\tilde{\bm{\mu}},\tilde{\Gamma})) \nonumber
\end{eqnarray}
from the assumption of (\ref{d2}).
Since
\[
   z(\bm{\xi},\tilde{\bm{\mu}},\tilde{\Gamma})=\frac{1}{2}
   \sum_{j=k+1}^p\{\tilde{\bm{\gamma}}_j^{\rm T}(\bm{\xi}-\tilde{\bm{\mu}})\}^2
   \ge \frac{1}{2}
   \{\tilde{\bm{\gamma}}_j^{\rm T}(\bm{\xi}-\tilde{\bm{\mu}})\}^2
\]
for any $j$, $k+1\le j \le p$ by the definition of $z$ at (\ref{zxmG}), we obtain 
\[
   \|{\rm IF}_\Psi(\bm{\xi};\tilde{P})\|
   \le \sum \left|\frac{\hat{\lambda}_j
   \|\tilde{\bm{\gamma}}_I\tilde{\bm{\gamma}}_j^{\rm T}+
   \tilde{\bm{\gamma}}_j\tilde{\bm{\gamma}}_I^{\rm T}\|}
   {\tilde{\lambda}_j(\hat{\lambda}_j-\hat{\lambda}_I)}\right|
   \left(d \sqrt{2} \sqrt{z(\bm{\xi},\tilde{\bm{\mu}},\tilde{\Gamma})}
   \psi(z(\bm{\xi},\tilde{\bm{\mu}},\tilde{\Gamma}))
   +d^2 \psi(z(\bm{\xi},\tilde{\bm{\mu}},\tilde{\Gamma}))\right).
\]
Finally, we conclude that
\[
   \|{\rm IF}_\Psi(\bm{\xi};\tilde{P})\|\le
   C\left(d \sqrt{2} \sup_{z>0} \sqrt{z}\psi(z)+d^2 \sup_{z>0}\psi(z)\right)
\]
for any outlier $\bm{\xi}$ in $\bm{R}^p$ satisfying (\ref{d2}), where
\[
   C=\sum \left| \frac{\hat{\lambda}_j
   \|\tilde{\bm{\gamma}}_I\tilde{\bm{\gamma}}_j^{\rm T}+
   \tilde{\bm{\gamma}}_j\tilde{\bm{\gamma}}_I^{\rm T}\|}
   {\tilde{\lambda}_j(\hat{\lambda}_j-\hat{\lambda}_I)}\right|.
\]
Therefore, the condition (\ref{scon}) for $\Psi$ leads to the boundedness of the influence function ${\rm IF}_\Psi$ conditional on the condition (\ref{d2}).


\section{Adaptive Selection for Tuning Parameters}

Let $\bm{\alpha}$ be a parameter in the generic function $\Psi$ which defines our objective function (\ref{E1}).
In this section, we focus on the role of the tuning parameter $\bm{\alpha}$.
The performance of PCA by the proposed method generated by $\Psi_\alpha$ depends on the value of the tuning parameter $\bm{\alpha}$.
As shown in (\ref{psi1}) and (\ref{psi2}), the generic function $\Psi_1$ involves tuning parameters $\beta$ and $\eta$, and $\Psi_2$ involves $\beta$.

Thus, generic functions control the sensitivity to outliers by these tuning parameters.
See Figure~1 for the graphs of the derivatives $\psi_1$ and $\psi_2$ of $\Psi_1$ and $\Psi_2$  for several tuning parameters.

\begin{figure}[htp]
\centerline{\psfig{figure=figure1.eps,width=12cm}}
\caption{Graphs of $\phi$ for several tuning parameters.}
\end{figure}

The generic function $\Psi_1$ where $\eta \rightarrow \infty $ or $\beta\rightarrow 0 $ yields the classical PCA.
If we exactly assume Gaussian distribution, or $\varepsilon=0$ in (\ref{Out1}), then the classical PCA or $\psi(z)\equiv 1$, is recommended as the standard method.
See Kamiya and Eguchi (2001) for a detailed discussion.
This suggests that under the situation in which $\varepsilon\neq 0$, there exists an optimal selection for tuning parameters giving a generic function other than $\psi(z) \equiv 1$.

We propose herein a method of determining tuning parameters as in (\ref{psi1}) and (\ref{psi2}), based on
K-fold cross-validation.
See Subsection 7.10 in Hastie {\it et al.}~(2001)
for the detailed discussion.
Throughout the section, we focus on batch data.

Let $F(\bm{x})$ be a data distribution which is assumed to generate input data $\bm{x}$.
Then, we adopt a generalization error function for assessing the performance of a given $\hat{\bm{\mu}}$ and $\hat\Gamma$ using
\[
   L(\hat{{\bm{\theta}}},F)=
   \int \cdots \int \Psi_0(z(\bm{x},\hat{{\bm{\mu}}},\hat{\Gamma}))dF(\bm{x})
\]
where 
\[
   \Psi_0(z)=\log \frac {1}{1+\exp\{-\beta_0(z\cdot\eta_0)\}}
\]
with $\bm{\alpha}_0=(\beta_0,\eta_0)=(50,10)$.
This choice of the error function is intended to achieve mild robustness to outliers.
This is because we cannot obtain a sensible result if the error function itself is sensitive to outliers.
The empirical error function is
\[
   L_{\rm emp}(\hat{\bm{\theta}})=\frac{1}{n}
   \sum_{t=1}^n \Psi_0(z(\bm{x}_t,\hat{\bm{\mu}},\hat{\Gamma}))
\]
for given data $\{\bm{x}_1,\ldots,\bm{x}_n\}$ and would be unchanged by data contamination if $(\hat{\bm{\mu}},\hat{\Gamma})$ is a robust estimator.
The choice of $\beta_0$ is universal, but that of $\eta_0$ should be adaptive. For example, the median of $\|\bm{x}_I-\bar{\bm{\mu}}\|$ with the usual centralizing vector $\bar{\bm{\mu}}$ as a default value.
Given a class of estimators $\hat{\bm{\theta}}_\alpha $ by the generic function $\Psi_\alpha $ with the tuning parameter $\alpha $, for example $\alpha=(\beta, \eta)$ in (\ref{psi1}) or (\ref{psi2}), we attempt to estimate the expected loss function, or the risk function associated with an estimator $\hat{\bm{\theta}}_\alpha $, which is essentially
\[
 R(\hat{\bm{\theta}}_\alpha,F) =
 {\bm E}_F\{ L(\hat{\bm{\theta}}_\alpha, F)\},
\]
where ${\bm  E}_F$ denotes the mathematical expectation when input data $\bm{x}_1,\cdots,\bm{x}_n$ follow from the underlying distribution $F$.

Here we provide a method of selecting $\bm{\alpha}^*$ which generates $\hat{\bm{\theta}}_{\alpha^*}$ with good performance.
We use K-fold validation to get a estimator of the generalization error $L(\hat{\bm{\theta}}_\alpha,F)$.
Here, we divide the data set $D=\{\bm{x}_1,\ldots,\bm{x}_n\}$ into K subsets $\{D_k=\{\bm{x}_1^{(k)},\ldots,\bm{x}_{n_k}^{(k)}\}:k=1,\cdots,K\}$ and so $D=\bigcup_{k=1}^K D_k$.
Define
\[
   CV(\bm{\alpha})=\frac{1}{K}
   \sum_{k=1}^K L_{\rm emp}^{(k)}(\hat{\bm{\theta}}_\alpha^{(-k)}),
\]
where $\hat{\bm{\theta}}_\alpha^{(-k)}$ is the estimator based on the data set $\bigcup_{k'\neq k} D_{k'}$ and
\[
   L_{\rm emp}^{(k)}(\bm{\theta})=\frac{1}{n_k}
   \sum_{I=1}^{n_k}L(\bm{x}_I^{(k)},\bm{\theta}).
\]
In this way, the estimator $\hat{\bm{\theta}}_\alpha$ and $D_k$ are statistically independent, which implies elimination of the bias in the empirical error $L_{\rm emp}(\hat{\bm{\theta}}_\alpha)$.
In this formulation, we now define the optimal $\bm{\alpha}^*$ by
\[
   \bm{\alpha}^*=\argmin_\alpha CV(\alpha).
\]
We will explore the performance of this selection method for synthetic data situations.

\section{Simulation Study}

We explore the robustness of our class of principal component subspaces in numerical experiments, focusing on the classical rule ($\psi(z)=1 $), the Xu-Yuille rule and the Gaussian kernel rule defined in (\ref{psi2}).
For our simulation study, we consider the following three types of outlier distributions $H$ in the $\varepsilon$-contamination model defined by (\ref{Out1}):
\begin{enumerate}
    \item Deterministic contamination: a sum of point-mass distributions at $\bm{x}=\bm{\xi}_j$ for $j=1,\cdots,M$.
    \item Structural contamination: the same Gaussian distribution $N({\bm{\mu}}_1,V_1)$, but with the structure in ${\bm{\mu}}_1$ and $V_1$ being quite 
   different from with that in ${\bm{\mu}}$ and $V$,
   that is to say, $\|\mu_1-\mu\|$ or ${\rm tr}(V_1-V)^2$ is
 substantially large.
    \item Distributional contamination: the same structure as in $({\bm{\mu}},V)$ but the distribution is 
    totally different from the Gaussian distribution $N(\bm{\mu},V)$.
\end{enumerate}
\begin{figure}[t]
\centerline{\psfig{figure=figure2.eps,width=10cm}}
\caption{The plot of inner product of the PC vectors with/without outliers against $\beta$ and $\eta$.}
\end{figure}
First, we investigate the case in which $\varepsilon$ is undetectably small, for example we take $\varepsilon=0.03$.
We have performed a numerical study for the behavior of our procedure for
seven-dimensional data in the following setting: $\bm{\mu} = (0,\cdots,0)^{\rm T}, V={\rm diag}(5,2,3,1,0.5,0.5,0.5)$.
As for (i), $\bm{\xi}_1 = (0,0,0,0,0,0,b)$, $\bm{\xi}_2 = (0,0,0,0,0,b,0)$ with a probability 0.5 for each.\\
In (ii) the distribution $H$ is a Gaussian distribution with
\[
\bm{\mu}_1 = (1,\cdots,1)^{\rm T}, V_1={\rm diag}(0.5,2,3,1,0.5,0.5,0.5).
\]
In (iii) the outlier, $\bm{\xi} $ has a distribution $H$ of Cauchy-type of which the location-scatter structure is the same as $(\bm{\mu},V)$, that is, all the components of $V^{-\frac{1}{2}}(\bm{\xi}-\bm{\mu}) $ are independently and identically distributed according to a standard Cauchy distribution with density function $1/(\pi(1+x^2))$.

We observe in a series of simulations in the above setting that the classical procedure $( \psi(z) = 1)$ breaks down for a batch of data with most observations from $N(\bm{\mu},V)$ and a few outliers from $H $.
The first principal component vector extracted only from the 98 simulated vectors is completely changed into the space of minor components after mixing with two outliers.
In our experience, PCA has never been weakly perturbed by outlier contamination under the situation of setting (i) with a small $b $.  
Whether the classical PCA resists or completely breaks down against outliers is determined by $b$.
If $b$ ranges from 14 to 16, the breakdown occurs with about 50 percent proportion.
If input data are simulated by setting (iii), then the classical PCA tends to break down with a higher frequency of occurrence than for setting (ii).
In almost all of the cases of setting (iii), the PCA breaks down.
Thus, we observe that the distributional contamination is more severe than the structural contamination for the classical PCA. 

We confine ourselves to a typical case of input vectors from the structural contamination.
We obtained 270 input vectors of 200-dimension from $N(\bm{\mu},V) $, where
\[
        \bm{\mu}=(0,\cdots,0)^{\rm T}, \quad V={\rm diag}(10,9,8,7,6,5,4,3,2,1,0.5, \cdots ,0.5),
\]
and 30 outliers from $N(\bm{\mu}_1,V_1) $, where
\[
        \bm{\mu}_1=(1,\cdots,1)^{\rm T}, \quad V_1={\rm diag}(1,9,8,7,6,5,4,3,2,1,1,1, \cdots,1),
\]
 and observed that the inner product of the first principal component vectors based on the data of 270 vectors and on the data added to 30 outliers is 0.678 when using the classical PCA.
On the other hand the inner product is 0.999 using procedure defined by (\ref{XY}) with $\beta = 0.5 $ and $\eta=130 $.
The RM algorithm is started with the initial vector $\bm{\gamma}=(1/\sqrt{200},\cdots,1/\sqrt{200})^{\rm T} $ and $\bm{\mu}=(0,\cdots,0)$, which assigns less weight to the 30 outliers after 10 iterations.

\begin{figure}[t]
\centerline{\psfig{figure=figure3.eps,width=10cm}}
\caption{The plot of 50 observations and 50 outliers in the minor subspace.}
\end{figure}

In this procedure, we heuristically choose the tuning parameters $\beta $ and $\eta $. 
We observe in Figure~2 that the performance is not so sensitive to the choice of $\beta $ if $\eta < 145$. In the subsequent simulation, we will investigate the data-adaptive selection for tuning parameters using the 10-fold CV method in Section 4.

Secondly, we investigate the case in which $\varepsilon$ is fairly large.
Hence, we take $\varepsilon=0.5$, which can be viewed as the worst case in our context.
We focus on the case of structural contamination with the same setting as that with $\varepsilon=0.03$.  
The situation is really an extreme case, and is beyond the usual context of outlier detection. 
Thus, we have several simulations involving this situation as the first step, as seen in Figure~3.
A typical result gives 0.101 as the inner product between the first principal component vectors with and without 50 outliers by the classical procedure. 
Alternatively our procedure gives 0.833, so we see that the procedure can detect a more sensible direction vector than the classical procedure.
The proposed procedure detects the heterogeneity of structure, as indicated in Figure~4.
If we have more information on the contamination or outlying structure, then we can build a shaper model for the outlier distribution.
For example, we might suggest a two component Gaussian mixture model and estimate the structure in a more complete situation via the EM algorithm.
However, to expect exact information on the outliers is often unrealistic in the present situation.

\begin{figure}[t]
\centerline{\psfig{figure=figure4.eps,width=10cm}}
\caption{The plot of the weight function of $\psi$ over 50 data with 50 outliers.}
\end{figure}



We apply the RM algorithm to on-line input vectors for comparison with the gradient learning algorithm.
For the simulation study, we assume a specific form of model (ii) with $\varepsilon=0.1$ and Gaussian density with $\bm{\mu}=\bm{0}$, $\bm{\mu}_1=(30,0,0,0,30)^{\rm T}$, $V={\rm diag}(9,7,5,3,1)$, and $V_1={\rm diag}(1,1,3,3,5)$.
Thus, the true principal component vector of $V$ is $(1,0,0,0,0)$ in the simulation design.
We observe that the RM algorithm is stable and attains rapid convergence from these on-line input vectors.
However, the computational burden is much heavier than the on-line gradient algorithm, as shown in Figure~5.

\begin{figure}[t]
\centerline{\psfig{figure=figure5.eps,width=12cm}}
\caption{The inner products of the true vector $(1,0,0,0)$ and the PC vectors by RM, Xu-Yuille gradient, the classical gradient and classical matrix algorithms.}
\end{figure}

We next investigate the numerical performance of the 10-fold CV method discussed in Section 4 under the model of the structural contamination as follows: 45 input vectors are simulated from a Gaussian density $N(0,{\rm diag}(9,7,5,3,1))$ and five outliers from $N((0,0,0,0,10),{\rm diag}(9,7,5,3,1))$.
We observe that the 10-fold CV method has detected the optimal tuning parameter $\eta=46$, as shown in Figures~6 (1) and (2), while the PCA is much less sensitive to $\beta$ than to $\eta$, so we fix the optimal tuning parameter as $\beta=1$.

\begin{figure}[t]
\centerline{\psfig{figure=figure6.eps,width=12cm}}
\caption{The plot of 10-fold CV for the Xu-Yuille rule for $\beta=1$ fixed.}
\end{figure}

We propose a robust procedure for centralization of the data in (\ref{mG}). 
In the neural networks literature, such a variant for centralizing data has been ignored until now.
Using the usual centralization is correct if all of the data are generated from Gaussian distribution.
This is because the mean vector and covariance matrix are orthonormal as parameters, so that any influence on the principal vectors is independent of that on the mean vector.
However, if the data in the mean vector are structured, this sometimes has a significant impact on the PCA.
Here we consider a simple simulation study for investigating the difference between two procedures defined by adopting the weighted sample mean vector $\hat{\mu}$ in the RM algorithm and the sample mean vector $\bar{\mu}$ in classical PCA as the centralizer.
We generate 180 observations from a Gaussian density $N(0,V)$ with $V={\rm diag}(9,7,5,3,1)$ and 20 outliers from $N((0,0,20,0,0)^{\rm T},V)$.
Figure~7 shows the two-dimensional score plot produced by the classical PCA based on only the 180 observations without any outliers, where the horizontal axis is taken as taken exactly as the first principal component vector.
We observe that the first principal component vector by $\hat{\mu}$-centralization yields a proper direction.

\begin{figure}[t]
\centerline{\psfig{figure=figure7.eps,width=11cm}}
\caption{The effect of centralization ways.}
\end{figure}

\

\section{Discussion}

We have discussed a class of procedures for  PCA based on the generic functions $\Psi$.
The derivative $\psi$ gives a weight expressing the degree of confidence 
for each input vector being an outlier.
The robust procedure for the PCA gives less weights to input vectors having long residual vectors when projected onto the principal component subspace. 
We emphasis that  the $\mu$-portion is defined to be a weighted mean with the same
weights as in the $\Gamma$-portion while the usual PCA
employs the naive centralization with constant weights.

Our major point is the adaptive selection of
 a set of tuning parameters which control the degree of robustification.
In empirical studies, we observe that the robustness performance 
is sensitive to the selection of tuning parameters.
K-fold cross validation properly gives the adaptive selection
for tuning parameters in accordance with data.
However the selection method is done only for batch data but it cannot be applied to on-line data, 
which we must post as a future research.
The RM algorithm needs the evaluation of eigenvalues and eigenvectors of the full matrix.
In this respect it requires  heavy computational burdens, whereas the convergence is 
stable and rapid relative to the gradient algorithm.
The RM algorithm must be improved when the dimension of the input vector is considerably high.
There is room for improvement in solving the $k$-dominant eigenvectors from a computational point
of view.


Another interesting issue would be the breakdown point of the method proposed in the present paper as a global measure of robustness.
In the previous literature the breakdown point has been 
considered as estimation of covariance (scatter) matrix
other than estimation of principal component vector.
However our method does not directly fit the theory since
the method is not only a function of covariance matrix.
We will need more discussion for this problem to be challenged as a future problem.

\

\bibliographystyle{apalike}
\begin{thebibliography}{10}
%\noindent
%{\bf References}

%\begin{list}{}{%
%    \setlength{\leftmargin}{0pt}
%    \setlength{\itemindent}{20pt}
%}

\item {Amari, S. -I. } Neural theory of association and concept formation. 
{\it Biological Cybernetics, \bf 26, \rm 175--185, 1977.}

\item {Campbell, N. A. } Robust procedures in multivariate analysis 1: Robust covariance estimation. 
{\it Applied Statistics \bf 29\rm, 231-237, 1980.}

\item {Caussinus, H. and Ruiz, A. } Interesting projections of multidimensional data by means of generalized principal component analysis. 
{\it COMPSTAT 90, \rm 121--126, 1990.}

\item {Croux, C. and Haesbroeck, G.} Principal component analysis based on robust estimators of the covariance or correlation matrix: Influence functions and efficiencies. 
{\it Biometrika, \bf 87\rm, 603-618, 2000.}

\item {De la Torre, F. and Black, M.}
 Robust principal component analysis for computer vision. 
 {\it International Conference on Computer Vision, \rm 2001.}

\item {Hampel, F. R.} The influence curve and its role in robust estimation.
{\it Journal of the American Statistical Association, \bf 69, \rm 383--393, 1974.}

\item {Hampel, F. R., Ronchetti, E. M., Rousseeuw, P. J. and Stahel, W. A. } 
{\it Robust Statistics: the Approach Based on Influence Functions.} Wiley, 1986.

\item {Hastie, T. J., Tibshirani, R. J. and Friedman, J. H. } 
{\it The Elements of Statistical Learning: Data Mining, Inference, and Prediction.} Springer-Verlag, 2001.

\item {Haykin, S. } {\it Neural Networks.} Prentice Hall, 1999.

\item {Higuchi, I. and Eguchi, S. } The influence function of principal component analysis by self-organizing rule.
{\it Neural Computation, \bf 10, \rm 1435--1444, 1998.}

\item {Hotelling, H. } Analysis of a complex of statistical variables into principal components. 
{\it Journal of Educational Psychology, \bf 24, \rm 417--441, 1933.}

\item {Huber, P. J. } {\it Robust Statistics.}  Wiley, 1981.

\item {Hyvarinen, A. and Oja, E. } A fast fixed-point algorithm for independent component analysis.
{\it Neural Computation, \bf 9, \rm 1483--1492, 1997}

\item {Kamiya, H. and Eguchi, S. } A class of robust principal component vectors. 
{\it Journal of Multivariate Analysis, \bf 76}, \rm 239--269, 2001. 

\item {McLachlan, G. J. and Krishnan, T. } 
{\it The EM Algorithm and Extensions.} Wiley, 1997.

\item {Oja, E. } A simplified neuron model as a principal component analyzer.
{\it Journal of Mathematical Biology, \bf 15, \rm 267--273, 1982.}

\item {Oja, E. } Neural networks, principal components and subspaces.
{\it International Journal of Neural Systems, \bf 1, \rm 61--68, 1989.}

\item {Tanaka, Y. } Sensitivity analysis in principal component analysis: influence on the subspace spanned by principal components. 
{\it Communications in Statistics -Theory and Methods, \bf 17, \rm 3157--3175, 1988.}

\item {Wu, C. F. J. } On the convergence properties of the EM algorithm. 
{\it Annals of Statistics, \bf 11, \rm 95--103, 1983.}

\item {Xu, L. and Yuille, A.} Robust principal component analysis by
self-organizing rules based on statistical physics approach.  
{\it IEEE Transactions on Neural Networks, \bf 6, \rm 131--143, 1995.}

%\end{list}
\end{thebibliography}


\end{document}
