Projection Sparse Principal Component Analysis: an efficient least squares method
Giovanni Maria Merola \hspace{.2cm
We propose a new sparse principal component analysis (SPCA) method in which the solutions are obtained by projecting the full cardinality principal components onto subsets of variables. The resulting components are guaranteed to explain a given proportion of variance. The computation of these solutions is very efficient. The proposed method compares well with the optimal least squares sparse components. We show that other SPCA methods fail to identify the best sparse approximations of the principal components and explain less variance than our solutions. We illustrate and compare our method with the analysis of a real dataset containing socioeconomic data and the computational results for nine datasets of increasing dimension with up to 16,000 variables.
PSPCA_arXiv_master.tex
\documentclass[12pt, a4paper, authoryear]{article}%preprint,
\pdfoutput=1
\usepackage[english]{babel}
\usepackage{amsmath}
%\usepackage{graphicx}%,psfrag,epsf}
\usepackage{enumerate}
\usepackage{natbib}
\usepackage[]{graphicx}
%\usepackage[]{color}
%% mine
\usepackage{amssymb}%
\usepackage{mathrsfs}
\usepackage{amsbsy}
\usepackage{amsthm}
\newtheorem{proposition}{Proposition}
\usepackage{enumitem}
\usepackage{booktabs}%% for tables
\usepackage{relsize}%% for fonts mathrm
%\usepackage{multirow}
\usepackage{rotating}
\usepackage{algorithm}
\usepackage{algpseudocode}
\usepackage{algorithmicx}
\usepackage{url}
\usepackage{caption}
\usepackage{array}
\usepackage[flushleft]{threeparttable}
\usepackage{bookmark}

\usepackage{natbib}
\usepackage{float}
\usepackage{caption}
\usepackage{subcaption}
\captionsetup{font = scriptsize}%,labelfont={bf,sf}}
%\captionsetup[sub]{font=small}%,labelfont={bf,sf}}
%\usepackage{subfig} deprecated
%\renewcommand{\figurename}{\bf{Fig.}}

\IfFileExists{upquote.sty}{\usepackage{upquote}}{}

\usepackage{listings}
\usepackage{color}

%\usepackage[modulo]{lineno}
%\setpagewiselinenumbers
%%\modulolinenumbers[5]
%
%\def\linenumberfont{\normalfont\small\sffamily}

%%%%%%%%%%%%%%%%%%%%%% end knitr ==================================================
%\immediate\write18{doFigures.bat}
%\immediate\write18{erase /Q figure/*.*}
%\immediate\write18{copy R-jss/figure/* figure}
%% end of declarations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% macros
\newcommand{\bti}[1]{{\tilde{\mathbf{#1}}}}
\newcommand{\refp}[1]{(\ref{#1})}
%$$ \bti{a}$$
\newcommand{\sub}[1]{{_{#1}}}

\newcommand{\subd}{{[d]}}
\newcommand{\subnd}{{{[nd]}}}
\newcommand{\subnr}{{{[nr]}}}
\newcommand{\subnc}{{{[nc]}}}

\newcommand{\bht}{{\hat{\mathbf{t}}}}
\newcommand{\bhX}{{\hat{\mathbf{X}}}}

\newcommand{\xa}{\bX_a}


\newcommand{\bas}{{\mathbf{a^*}}}

\newcommand{\ba}{{\mathbf{a}}}
\newcommand{\bc}{{\mathbf{c}}}
\newcommand{\bd}{{\mathbf{d}}}
\newcommand{\bb}{{\mathbf{b}}}
\newcommand{\be}{{\mathbf{e}}}
\newcommand{\bt}{{\mathbf{t}}}
\newcommand{\bs}{{\mathbf{s}}}
\newcommand{\bp}{{\mathbf{p}}}
\newcommand{\bq}{{\mathbf{q}}}
\newcommand{\bu}{{\mathbf{u}}}
\newcommand{\bx}{{\mathbf{x}}}
\newcommand{\by}{{\mathbf{y}}}
\newcommand{\br}{{\mathbf{r}}}
\newcommand{\bbw}{{\mathbf{w}}}
\newcommand{\bz}{{\mathbf{z}}}

\newcommand{\bv}{{\mathbf{v}}}
\newcommand{\bX}{{\mathbf{X}}}
\newcommand{\bY}{{\mathbf{Y}}}
\newcommand{\bA}{{\mathbf{A}}}
\newcommand{\bB}{{\mathbf{B}}}
\newcommand{\bC}{{\mathbf{C}}}
\newcommand{\bD}{{\mathbf{D}}}
\newcommand{\bE}{{\mathbf{E}}}
\newcommand{\bF}{{\mathbf{F}}}
\newcommand{\bG}{{\mathbf{G}}}
\newcommand{\bK}{{\mathbf{K}}}
\newcommand{\bT}{{\mathbf{T}}}
\newcommand{\bL}{{\mathbf{L}}}
\newcommand{\bM}{{\mathbf{M}}}
\newcommand{\bV}{{\mathbf{V}}}
\newcommand{\bU}{{\mathbf{U}}}
\newcommand{\bZ}{{\mathbf{Z}}}
\newcommand{\bW}{{\mathbf{W}}}
\newcommand{\bR}{{\mathbf{R}}}
\newcommand{\bw}{{\mathbf{w}}}
\newcommand{\bS}{{\mathbf{S}}}
\newcommand{\bO}{{\mathbf{O}}}
\newcommand{\bI}{{\mathbf{I}}}
\newcommand{\bJ}{{\mathbf{J}}}
\newcommand{\bP}{{\mathbf{P}}}
\newcommand{\bQ}{{\mathbf{Q}}}

\newcommand{\bta}{{\tilde{\mathbf{a}}}}
\newcommand{\btt}{{\tilde{\mathbf{t}}}}
\newcommand{\btv}{{\tilde{\mathbf{v}}}}
\newcommand{\btw}{{\tilde{\mathbf{w}}}}
\newcommand{\btA}{{\tilde{\mathbf{A}}}}
\newcommand{\btV}{{\tilde{\mathbf{V}}}}
\newcommand{\btU}{{\tilde{\mathbf{U}}}}
\newcommand{\btT}{{\tilde{\mathbf{T}}}}
\newcommand{\btS}{{\tilde{\mathbf{S}}}}
\newcommand{\btR}{{\tilde{\mathbf{R}}}}
\newcommand{\btW}{{\tilde{\mathbf{W}}}}
\newcommand{\btQ}{{\tilde{\mathbf{Q}}}}
\newcommand{\btY}{{\tilde{\mathbf{Y}}}}

\newcommand{\bLa}{{\mathbf{\Lambda}}}
\newcommand{\bDe}{{\mathbf{\Delta}}}
\newcommand{\bSi}{{\mathbf{\Sigma}}}
\newcommand{\bla}{{\boldsymbol{\lambda}}}
\newcommand{\bOm}{{\boldsymbol{\Omega}}}
\newcommand{\bGa}{{\boldsymbol{\Gamma}}}
\newcommand{\bTh}{{\boldsymbol{\Theta}}}
\newcommand{\bPhi}{{\boldsymbol{\Phi}}}



%\newcommand{\cX}{\check{\mathbf{X}}}
\newcommand{\cX}{{\overset{\boldsymbol{\vee}}{\mathbf{X}}}}
\newcommand{\dX}{{\overset{\boldsymbol{.}}{\mathbf{X}}}}
\newcommand{\dQ}{{\overset{\boldsymbol{.}}{\mathbf{Q}}}}
\newcommand{\da}{{\overset{\boldsymbol{.}}{\mathbf{a}}}}

%\newcommand{\cX}{\vee{\mathbf{X}}}
%\newcommand{\dX}{\dot{\mathbf{X}}}
%\newcommand{\dQ}{\dot{\mathbf{Q}}}
%\newcommand{\da}{\dot{\mathbf{a}}}


\newcommand{\bdX}{{\overset{\boldsymbol{.}}{\mathbf{X}}}}
\newcommand{\bda}{{\overset{\boldsymbol{.}}{\mathbf{a}}}}
\newcommand{\bdt}{{\overset{\boldsymbol{.}}{\mathbf{t}}}}
\newcommand{\bdv}{{\overset{\boldsymbol{.}}{\mathbf{v}}}}
\newcommand{\bdw}{{\overset{\boldsymbol{.}}{\mathbf{w}}}}
\newcommand{\bdA}{{\overset{\boldsymbol{.}}{\mathbf{A}}}}
\newcommand{\bdV}{{\overset{\boldsymbol{.}}{\mathbf{V}}}}
\newcommand{\bdU}{{\overset{\boldsymbol{.}}{\mathbf{U}}}}
\newcommand{\bdT}{{\overset{\boldsymbol{.}}{\mathbf{T}}}}
\newcommand{\bdQ}{{\overset{\boldsymbol{.}}{\mathbf{Q}}}}
\newcommand{\bdR}{{\overset{\boldsymbol{.}}{\mathbf{R}}}}
\newcommand{\bdS}{{\overset{\boldsymbol{.}}{\mathbf{S}}}}
\newcommand{\bdY}{{\overset{\boldsymbol{.}}{\mathbf{Y}}}}
\newcommand{\bdW}{{\overset{\boldsymbol{.}}{\mathbf{W}}}}




\newcommand{\btLa}{{\tilde{\boldsymbol{\Lambda}}}}
\newcommand{\btOm}{\tilde{\boldsymbol{\Omega}}}
\newcommand{\btTh}{\tilde{\boldsymbol{\Theta}}}

\newcommand{\bdLa}{{{\overset{\boldsymbol{.}}}{\boldsymbol{\Lambda}}}}
\newcommand{\bdOm}{{\overset{\boldsymbol{.}}{\boldsymbol{\Omega}}}}
\newcommand{\bdTh}{{\overset{\boldsymbol{.}}{\boldsymbol{\Theta}}}}
\newcommand{\bdGa}{{\overset{\boldsymbol{.}}{\boldsymbol{\Gamma}}}}







\newcommand{\nc}{{\mathrm{nc}}}
\newcommand{\un}{{uncorrelated}}
\newcommand{\Un}{{Uncorrelated}}
\newcommand{\uns}{{uncorrelatedness}}
\newcommand{\Uns}{{Uncorrelatedness}}
%\newcommand{\trasp}{^{\prime}}
%\newcommand{\trasp}{^{\top}}

%^\mathsf{T}
\makeatletter
\newcommand*{\trasp}{%
  {\mathpalette\@trasp{}}%
}
\newcommand*{\@trasp}[2]{%
  % #1: math style
  % #2: unused
  ^{\raisebox{\depth}{$\m@th#1\scriptstyle{\intercal}$}}%\mathsf{T}
}
\makeatother






\newcommand{\eg}{{\emph{e.g.}}}
\newcommand{\vexp}{vexp}
\newcommand{\vexpq}{{vexp_Q}}
\newcommand{\evexp}{{evexp}}

\newcommand{\SPCA}{\emph{SPCA}}
\newcommand{\diag}{\text{diag}}


\DeclareMathOperator*{\argmin}{arg\,min}
\DeclareMathOperator*{\argmax}{arg\,max}
\DeclareMathOperator*{\Max}{max}


\newcommand{\Tr}{{\mathrm{tr}}}
\newcommand{\Lo}{{\text{L}_0}}
\newcommand{\Li}{\text{L}_1}
\newcommand{\Lt}{\text{L}_2}
\newcommand{\sq}{\frac{1}{2}}
\newcommand{\msq}{{-\frac{1}{2}}}

%%%


\def\spacingset#1{\renewcommand{\baselinestretch}{#1}}
%\small\normalsize} \spacingset{1}

\newcommand\bsln{1}
\renewcommand\baselinestretch{1.0} %for jasa{1.425}
\setlength{\parskip}{0\baselineskip}%
%\setlength{\parindent}{0pt}%

%\usepackage{draftwatermark}
%\SetWatermarkText{PREPRINT}
%\SetWatermarkScale{4}
%\SetWatermarkColor[gray]{0.85}
%\SetWatermarkLightness{0.15}

\begin{document}
%% almost as usual
  \title{\bf Projection Sparse Principal Component Analysis: an efficient least squares method}
  \author{Giovanni Maria Merola \hspace{.2cm}\\
    Department of mathematical Sciences\\ 
    Xi'an Jioatong Liverpool University\\
    giovanni.merola@xjtlu.edu.cn}
  \maketitle
\begin{abstract}
We propose a new sparse principal component analysis (SPCA) method in which the solutions are obtained by projecting the full cardinality principal components onto subsets of variables. The resulting components are guaranteed to explain a given proportion of variance. The computation of these solutions is very efficient. The proposed method compares well with the optimal least squares sparse components. We show that other SPCA methods fail to identify the best sparse approximations of the principal components and explain less variance than our solutions. We illustrate and compare our method with the analysis of a real dataset containing socioeconomic data and the computational results for nine datasets of increasing dimension with up to 16,000 variables.
\end{abstract}
%\linenumbers

%\lhead[\footnotesize\thepage\fancyplain{}\leftmark]{}\rhead[]{\fancyplain{}\rightmark\footnotesize\thepage}%Put
\section{Introduction}
Data dimensionality reduction methods approximate a set of variables on a lower
dimensional space. Principal components analysis (PCA) is the oldest and most popular
of these methods. Given the measurements of $p$ variables (or features) on $n$
subjects (or objects) the principal components (PCs) are the linear combinations of
the variables that best approximate the data on an orthogonal subspace with respect to
the least squares criterion (LS) \citep{pea}.

The PCs  are difficult to interpret and not informative on important features of the
dataset because they are combinations of all the observed variables, as already
pointed out by \citet{jef}. Sparse principal components analysis (SPCA) is a variant
of PCA in which the components computed are combinations of just a few of the observed variables. The aim of SPCA is to produce a reduction of the dimensionality that is easy to interpret.

Formally, let $\bX$ be an ($n\times p$) matrix containing $n$ observations on $p$ mean
centred variables.  The information contained in the dataset is summarised by its
\emph{total variance} defined by $||\bX||^2 = trace(\bS)$, where $\bS =
\bX\trasp\bX$ is the covariance matrix (omitting a scaling constant which is irrelevant for our analysis).

A \emph{component} is any linear combination of the columns of $\bX$, defined by $\bt
= \bX\ba$, where the elements of the vector $\ba$ are called \emph{loadings}.
The least squares estimates of $c<p$ components $\bT = (\bt_1, \ldots, \bt_c)$
are obtained by minimising the squared (Frobenius) norm of the difference of the data matrix from its projection onto the components, $\widehat\bX(\bT)$. That is, they are the
solutions to
\begin{equation}\label{eq:vexp}
\bT = \argmin\limits_{\bT}   ||\bX - \widehat\bX(\bT)||^2 = \argmax\limits_{\bT} ||\widehat\bX(\bT)||^2.
\end{equation}
Further constraints are necessary in order to identify the individual components. The \emph{variance explained}, $||\widehat\bX(\bT)||^2$, maximised by the least squares criterion is used to measure the performance of the components in approximating the data.

The PCs are the mutually orthogonal components that sequentially explain the largest proportion of variance of the $\bX$ matrix. The (vectors of) loadings of the PCs are the eigenvectors of the covariance matrix $\bS$ and the variance explained by each PC is equal to its norm. For this reason the PCs are commonly defined as the components with unit norm loadings that sequentially have largest norm \citep{hot}. %In passing, we note that the PCs are also solutions to the singular value decomposition problem \citep{eck}.

A sparse component is a linear combination of a subset of columns of the $\bX$ matrix,
or \emph{block} of variables, $\dX$, defined by the loadings $\da$ as $\bt = \dX\da$. The number of variables in the block is referred to as the \emph{cardinality} of the loadings. Sparse components are not eigenvectors of the covariance matrix and the variance that they explain is not equal to their norm. Since, in most cases they are not required to be orthogonal, new, more general measures of their performance are needed.

In recent years a great many SPCA methods have been proposed \citep[among others,][]{jol00, mog, zou, sri, she, wan}. These use different optimisation approaches to find sparse components that have largest norm.
Because of the adoption of this objective function the resulting components do not maximise the variance explained. Furthermore these methods cannot reproduce the sparsest expression of the PCs when the data is rank deficient and the components can be combinations of perfectly correlated variables, as we will show later.

In this paper we adopt the least squares sparse components analysis (LS SPCA, \citealp{mer}) approach to derive sparse components. In this approach the components are derived by minimising the least squares estimation criterion under sparsity constraints. LS SPCA components explain the most possible variance for a given block of variables and do not suffer from the same drawbacks mentioned above for other methods.

Identifying the blocks of variables from which to compute the components is an intractable NP-hard problem \citep{mog}. Therefore, greedy approaches are used to solve the SPCA problem. The techniques used to find components with the largest norm used in conventional SPCA are not appropriate for the maximisation of the variance explained. \citet{mer} proposed a backward elimination procedure which gives excellent results but does not scale well to larger datasets.

In this paper we propose computing sequentially the sparse components by projecting each PC on blocks of variables that explain a predetermined proportion of their variance. We show that the total variance explained by these projections is not less than the variance of the PC explained by the variables in the blocks. This means that we can select the variables in the blocks with an economical supervised forward selection and then compute the sparse component with a simple orthogonal projection. We will refer to this method as projection SPCA.

We also show that the lower bound for the variance explained by the projection SPCA components is also valid for the LS SPCA components computed from the same blocks of variables. Furthermore, we provide a novel interpretation of the LS SPCA components showing that these are the PCs of the projection of the data matrix onto the current block of variables.

We propose a simple but efficient algorithm for computing projection sparse components. We will show that sparse components for datasets with over 16,000 variables can be computed in a matter of seconds with this algorithm. Empirically we find that the time complexity of the algorithm is slightly larger than $O(p^{2})$.

The paper is organised as follows: in the next section we derive projection SPCA comparing it to LS SPCA and to conventional SPCA methods. In Section \ref{sec:comp} we give computational details. In Section \ref{sec:examples} we illustrate the usefulness of projection SPCA with an application on a real set of social data and show its performance on different large datasets also comparing it with that of conventional SPCA. In Section \ref{sec:conclusions} we give some final remarks.%
\section{Projection SPCA} In this section, for ease of exposition, we will refer to quantities computed in a generic step of an algorithm without the index identifying the order of the step; the index will be used only when needed.

We denote a generic vector of loadings as $\ba$. When sparsity is included in the model the symbol $\da$  denotes the vector of $c<p$ non-zero loadings corresponding to the  variables in the block $\dX$. In some cases we write the sparse vector $\da$ in its full cardinality expression as $\ba$ where missing values are set equal to zero. Principal components will be denoted as $\bu$ and their loadings as $\bv$.
%The symbol $\bJ$ denotes the $(p\times c)$ matrix formed by the columns of the $(p\times p)$ identity matrix with indices in $ind$. Then the full  vector of loadings (containing the zero elements) is equal to $\ba = \bJ\da$ and the sparse component can be written as $\bt = \bX\ba = \bX\bJ\da = \dX\da$ with $\dX = \bX\bJ$.
\subsection{Relationship between variance explained and norm of the components}\label{sec:vexppca}
When estimating a set of components, the minimisation of the least squares criterion is equivalent to the maximisation of the variance explained, as shown in Equation (\ref{eq:vexp}). For a set of components $\bT = (\bt_1,\ldots,\bt_c) = \bX(\ba_1, \ldots\ba_c)$ this is the total variance of the projection of $\bX$ onto $\bT$. If we denote the projector onto $\bT$ by $\bP_\bT = \bT(\bT\trasp\bT)^{-1}\bT\trasp$ so that this projection is $\widehat\bX(\bT) = \bT(\bT\trasp\bT)^{-1}\bT\trasp\bX$, the variance explained by the components in $\bT$ is equal to
\begin{equation*}
vexp(\bT) = ||\widehat\bX(\bT)||^2 = trace(\bX\trasp\bT(\bT\trasp\bT)^{-1}\bT\trasp\bX). %= trace(\bS\bA(\bA\trasp\bS\bA)^{-1}\bA\trasp\bS) \label{eq:vexpgenm}
\end{equation*}

In PCA the variance explained is maximised simultaneously for all components, $\bU = (\bu_1,\ldots,\bu_{c}) = \bX(\bv_1,\ldots,\bv_{c})$, under the constraint that these are mutually orthogonal. Because of these constraints, the variance explained by the $c$ components can be broken down into the sum of the variances of the complement of $\bX$ orthogonal to the previously computed components explained by each component. This complement is the difference of the data matrix from its projection onto the previous components $\bX - \widehat\bX(\bU)$ and is equal to $\bQ = \bX - \bU(\bU\trasp\bU)^{-1}\bU\trasp\bX$ . Therefore, the variance explained by the next principal component $\bu=\bX\bv$ is equal to
\begin{equation*}\label{eq:vexppcind} vexp(\bu) = \frac{\bu\trasp\bQ\bQ\trasp\bu}{\bu\trasp\bu} = \frac{\bv\trasp\bQ\trasp\bQ\bQ\trasp\bQ\bv}{\bv\trasp\bQ\bQ\trasp\bv},
\end{equation*}
where the second equality comes from the orthogonality constraints for which $\bu\trasp\bX = \bu\trasp\bQ = \bv\trasp\bQ\trasp\bQ$. Setting the derivatives of $vexp(\bu)$ with respect to $\bv$ equal to zero gives
\begin{equation*}
\frac{\mathop{d vexp(\bu)}}{\mathop{d\bv}} = \bQ\trasp\bQ\bQ\trasp\bQ\bv - \bQ\trasp\bQ\bv\lambda \equiv \bQ\trasp\bQ\bv - \bv\lambda = \mathbf{0},
\end{equation*}
with $\lambda= vexp(\bu) = \bv\trasp\bQ\trasp\bQ\bQ\trasp\bQ\bv/ \bv\trasp\bQ\trasp\bQ\bv$. For the $j$-th PC the solution that maximises the variance explained is the first eigenvector of $\bQ_j\trasp\bQ_j$, which is equal to the $j$-th eigenvector of $\bS$. Therefore, $vexp(\bu_j)$ simplifies to
\begin{equation}\label{eq:vexppc}
vexp(\bu_j) = \frac{\bv_j\trasp\bQ_j\trasp\bQ_j\bQ_j\trasp\bQ_j\bv_j}{\bv_j\trasp\bQ_j\trasp\bQ_j\bv_j} = \frac{\bv_j\trasp\bQ_j\trasp\bQ_j\bv_j\lambda_j}{\bv_j\trasp\bv_j\lambda_j} = \frac{\bv_j\trasp\bS\bv_j}{\bv_j\trasp\bv_j} = \lambda_j,
\end{equation}
because $\bQ_j\trasp\bQ_j\bv_j = \bS\bv_j = \bv_j\lambda_j$, where $\lambda_j$ is the $j$-th eigenvalue of $\bS$. This equation is the connection between Pearson's least squares definition of PCA and Hotelling's maximal norm definition.

Since the PCs loadings are eigenvectors of $\bS$, the orthogonality constraints $\bv_j\trasp\bS\bv_k=0, j\neq k$ are equivalent to the constraints $\bv_j\trasp\bv_k = 0, j\neq k$. By adding the requirement that $\bv_j\trasp\bv_j\ = 1$ Equation (\ref{eq:vexppc}) simplifies to Hotelling's definition of PCA:
\begin{eqnarray}\label{eq:pcahot}
\bv_j &=& \argmax_{\bv_j} \bv_j\trasp\bS\bv_j\\
&&\bv_j\trasp\bv_k = \delta_{jk},\nonumber
\end{eqnarray}
where $\bv_j\trasp\bS\bv_j=\lambda_j$ is the norm of the components with unit norm loadings (hereafter referred to simply as norm of the components) and  $\delta_{ij}$ is the Kronecker delta..


The simplified expression of the variance explained by a component in Equation (\ref{eq:vexppc}) is valid if and only if the loadings are eigenvectors of $\bS$. Hence it is not valid for a sparse component $\bt = \dX\da$ because in this case the first order conditions are
\begin{equation*}
\frac{\mathop{d vexp(\bt)}}{\mathop{d\da}} = \dQ\trasp\bQ\bQ\trasp\dQ\da - \dQ\trasp\dQ\da\gamma = \mathbf{0},
\end{equation*}
which cannot be simplified as above and the solutions are not eigenvectors of $\bS$. This should be expected because, in general, eigenvectors are not sparse. For this reason, the properties of the PCs are not valid for sparse components; in particular, the norm of the sparse components is not equal to the variance that they explain and the loadings are not required to be orthonormal \citep[see][for a more details]{mer}.
\subsection{Extra variance explained}\label{sec:vexp}
Since sparse loadings are not eigenvectors of the covariance matrix and sparse components are not necessarily orthogonal, we need to introduce new measures for the variance explained.

Let $\bT = (\bt_1,\ldots,\bt_c)$ be the matrix containing the first $c< p$ components computed. If a new component, $\bt$, that enters the model is correlated with the previously computed components --- that is, if $\bt_i\trasp\bt_j\not= 0,\ i < j$ --- the variance explained
\begin{equation*}\label{eq:vexpgen}
vexp(\bt) = \frac{\ba\trasp\bX\trasp\bX\bX\trasp\bX\ba}{\ba\trasp\bX\trasp\bX\ba}
\end{equation*}
includes part of the variance already explained by the other components. The net increment in total variance explained due to the new component is the variance explained by its residuals orthogonal to the components already in the model. Borrowing from the regression literature, we will refer to the net variance explained as \emph{extra} variance explained, and denote it by $evexp$.

%THIS ALREADY INTRODUCED\\
%Let $\bT = \{\bt_1,\ldots,\bt_c\}$ be the matrix containing the first $c< p$
%components computed, the orthogonal projector onto the columns of $\bT$ is equal to
%$\bP_{\bT}= \bT(\bT\trasp\bT)^{-1}\bT\trasp$; the residuals of $\bX$ orthogonal to
%$\bT$ are the difference of the data matrix from its projection onto the components
%$\bX - \widehat\bX(\bT)$ and are equal to
%\begin{equation*}
%\bQ = \bX - \bP_{\bT}\bX,
%\label{eq:q}
%\end{equation*}
%which we will refer to as \emph{deflated} $\bX$ matrix.
The orthogonal projector onto the columns of a set of components $\bT$ is equal to $\bP_{\bT}= \bT(\bT\trasp\bT)^{-1}\bT\trasp$; the residuals of $\bX$ orthogonal to $\bT$ are equal to
\begin{equation*}
\bQ = \bX - \bP_{\bT}\bX,
\label{eq:q}
\end{equation*}
which we will refer to as \emph{deflated} $\bX$ matrix. A new component $\bt = \bX\ba$ deflated of the previous components is equal to $\bt -\bP_{\bT} \bt = \bQ\ba$ and the extra variance that it explains is equal to
\begin{equation}
evexp(\bt) = ||\widehat\bX(\bQ\ba)||^2 =
\frac{\ba\trasp\bQ\trasp\bX\bX\trasp\bQ\ba}{\ba\trasp\bQ\trasp\bQ\ba} =
\frac{\ba\trasp\bX\trasp\bQ\bQ\trasp\bX\ba}{\ba\trasp\bQ\trasp\bQ\ba}.
%=\frac{\bt\trasp\bQ\bQ\trasp\bt}{\ba\trasp\bQ\trasp\bQ\ba}
\label{eq:evexp}
\end{equation}
The extra variance explained is equal to the variance explained when $\bt$ is the first component (in which case $\bQ = \bX$) and when it is uncorrelated with all the previously computed components. Since the PCs are mutually uncorrelated, the concept of extra variance explained does not need to be considered in their evaluation.

The extra variance explained is essential for the computation of the sparse components because it measures their contribution accounting for the components already computed. Furthermore, the total variance explained by a set of correlated components $(\bt_1, \ldots, \bt_c)$ can be expressed as the sum of the extra variances explained as
%\begin{equation*}
$
vexp(\bt_1, \ldots, \bt_c) = \sum_{j = 1}^c evexp(\bt_j).
$
%\end{equation*}
This formulation of the total variance explained is a generalisation of that commonly used in PCA.

A quantity that we will consider later is the variance of $\bQ$ explained by a component. For a component $\bt = \bX\ba$ this is equal to
\begin{equation*}
vexp_Q(\bt) =  ||\widehat\bQ(\bt)||^2 = \frac{\bt\trasp\bQ\bQ\trasp\bt}{\bt\trasp\bt}
=
\frac{\ba\trasp\bX\trasp\bQ\bQ\trasp\bX\ba}{\ba\trasp\bX\trasp\bX\ba} =
evexp(\bt) \frac{\ba\trasp\bQ\trasp\bQ\ba}{\ba\trasp\bX\trasp\bX\ba}.
\label{eq:vexpq}
\end{equation*}
Since ${\ba\trasp\bX\trasp\bX\ba}\geq {\ba\trasp\bQ\trasp\bQ\ba}$ it follows that
\begin{equation}
vexp_Q(\bt) \leq evexp(\bt) \leq vexp(\bt).
\label{eq:ineqvexp}
\end{equation}
$vexp_Q$ will be useful for deriving the sparse components because it is defined by the components expressed in terms of the observed variables (and not of their orthogonal residuals) and is a lower bound for $evexp$.
\subsection{Derivation of the projection sparse components}
The idea underpinning projection SPCA is to derive the sparse components at each step from a set of variables that explain a given percentage of the (full cardinality) first principal component of the deflated matrix $\bQ$, which is the component that explains the most extra variance in that step.

Let $\bu= \bQ\bv$ be the first PC of $\bQ$, assuming without loss of generality that $\bv\trasp\bv=1$. Since $\bu$ is orthogonal to the previous components $vexp(\bu) = evexp(\bu)$. By Equation (\ref{eq:vexppc}) the (extra) variance explained by the principal component $\bu$ simplifies to $vexp(\bu) = \bu\trasp\bu = \lambda_{\max}$, where $\lambda_{\max}$ is the largest eigenvalue of $\bQ\trasp\bQ$. The variance explained by the component $\bu$ is an upper bound for the extra variance explained by any component $\bt = \bX\ba$ added to the model, in symbols: $evexp(\bt)\leq \lambda_{\max}$.

Assume that the variables in a block $\dX$ explain a proportion of the variance of $\bu$ not less than $0< \alpha \leq 1$; that is, that $\hat\bu = \dX(\dX\trasp\dX)^{-1}\dX\trasp\bu$ is such that \begin{equation}\label{eq:rsqu} \frac{\hat\bu\trasp\hat\bu}{\bu\trasp\bu} \geq \alpha\; \mathrm{or, equivalently,} \; \hat\bu\trasp\hat\bu \geq \alpha\lambda_{\max}, \end{equation} where the right-hand-side inequality derives from Equation (\ref{eq:vexppc}). The projection $\hat\bu$ is a sparse component defined by $\hat\bu = \dX\hat\bv$, with
\begin{equation*}
    \hat\bv = (\dX\trasp\dX)^{-1}\dX\trasp\bu.
\end{equation*}
The variance of $\bQ$ explained by $\hat\bu$ is equal to:
\begin{equation}
vexp_Q(\hat\bu) = \frac{\hat\bu\trasp\bQ\bQ\trasp\hat\bu}{\hat\bu\trasp\hat\bu} =
\frac{\hat\bv\trasp\dX\trasp\bQ\bQ\trasp\dX\hat\bv}{\hat\bv\trasp\dX\trasp\dX\hat\bv}.
\label{eq:vexppc1}
\end{equation}

It is evident that if $\hat\bu$ explains more than a proportion $\alpha$ of the variance of the first PC of $\bQ$ it will explain a larger proportion of the variance of the whole matrix. This can be seen by writing in a singular value decomposition fashion $\bQ = \bu\bv\trasp + \bU_{/1}\bV_{/1}\trasp$ (with $\bu_j\trasp\bu_j=\lambda_j$), where the subscript `${/1}$' denotes that the first column of the matrix has been removed. Then, since $\bv\trasp\bV_{/1} =\mathbf{0}$,
\begin{equation}\label{eq:evexpuhat}
evexp(\hat\bu) \geq \vexpq(\hat\bu) = \hat\bu\trasp\hat\bu +
\frac{\hat\bu\trasp\bU_{/1} \bU_{/1}\trasp\hat\bu}{\hat\bu\trasp\hat\bu} \geq
\hat\bu\trasp\hat\bu \geq \alpha \lambda_{\max},
\end{equation}
because of Equations (\ref{eq:ineqvexp}) and (\ref{eq:rsqu}). This proves that the extra variance explained by a projection sparse component $\hat\bu$ is not less than $\alpha\lambda_{\max}$.
\subsection{Connection with least squares SPCA}
\citet{mer} derives the LS SPCA components from a block $\dX$, $\bt = \dX\da$, by maximising $vexp_Q(\bt)$. If the block explains at least a proportion $\alpha$ of the current PC, as in projection SPCA, from Equations (\ref{eq:vexppc1}) and
(\ref{eq:evexpuhat}) it follows that the extra variance explained by an LS SPCA component is not less than $\alpha\lambda_{\max}$, because
\begin{equation}
    evexp(\bt) \geq  \vexpq(\bt) =  \max\limits_a
    \frac{\da\trasp\dX\trasp\bQ\bQ\trasp\dX\da}{\da\trasp\dX\trasp\dX\da} \geq
    \frac{\hat\bv\trasp\dX\trasp\bQ\bQ\trasp\dX\hat\bv}{\hat\bv\trasp\dX\trasp\dX\hat\bv}
    \geq \alpha \lambda_{\max}.
\label{eq:evexplsspca}
\end{equation}
%This means that both projection SPCA and LS SPCA components computed from the same block of variables are guaranteed to explain more than an arbitrary proportion of the maximum variance explainable at each step.

The LS SPCA loadings are the generalised eigenvectors solution to
\begin{equation}
\dX\trasp\bQ\bQ\trasp\dX \da = \dX\trasp\dX\da\gamma_{\max},
\label{eq:lsspca}\end{equation}
where $\gamma_{\max}$ is the largest generalised eigenvalue and is equal to the variance of $\bQ$ explained by the component.

If we consider the orthogonal projection of $\bQ$ onto $\dX$, $\widehat\bQ = \dX(\dX\trasp\dX)^{-1}\dX\trasp\bQ$, the loadings of the first principal component of $\widehat\bQ$, $\bb$, say, must satisfy:
\begin{equation*}
\widehat\bQ\trasp\widehat\bQ\bb =  \bQ\trasp\dX(\dX\trasp\dX)^{-1}\dX\trasp\bQ \bb
= \bb\gamma_{\max}.
\end{equation*}
Then the vector $\da = (\dX\trasp\dX)^{-1}\dX\trasp\bQ\bb$ satisfies
\begin{equation}\label{eq:lspcapcsproj}
(\dX\trasp\dX)^{-1}\dX\trasp\bQ\bQ\trasp\dX \da = \da\gamma_{\max},
\end{equation}

which is equivalent to Equation (\ref{eq:lsspca}). Hence the LS SPCA components are the first PCs of the projection of the deflated $\bQ$ matrix onto the blocks $\dX$, while the projection SPCA components are the projections of the first PCs of $\bQ$ onto the blocks $\dX$.

By Equation (\ref{eq:evexplsspca}) the variance explained by LS SPCA solutions has a higher lower bound than that explained by projection SPCA components. However, due to the necessity of using greedy approaches to the identification of the blocks of variables, it is not always guaranteed that the LS SPCA components explain more variance than the corresponding projection SPCA components. We leave the comparison of these two SPCA approaches for future research. We can anticipate that in our studies we observed that the extra variance explained is very similar (to the third or fourth decimal figure) for all components.
%However, the cardinality and variables selected begin to diverge slightly as the number of components computed increases.
\subsection{Connection with conventional SPCA methods}
In conventional methods the SPCA problem is defined by adding sparsity constraints to Hotelling's formulation of PCA in Equation (\ref{eq:pcahot}). That is, the loadings are the solution to
\begin{eqnarray}\label{eq:spcaprob}
    \ba &=& \argmax\limits_{\ba} \ba\trasp\bS\ba\\
    &&card(\ba) = c,\;\ba\trasp\ba=1.\nonumber
\end{eqnarray}
%It is simple to see that the solution is the first eigenvector of the principal sub-matrix of $bS$ of size $c\times c$ with the maximal largest eigenvalue \citep{mog}.
In some cases the cardinality constraint is substituted with a lasso type penalty $\sum|a_i|$. It is well known that this penalty induces sparsity. Different constraints are added after the first component is computed to avoid trivial solutions.

For most SPCA methods the adoption of model \ref{eq:spcaprob} is justified by a not clearly defined \emph{statistical fidelity} criterion. This probably refers to the assumption that under sparsity constraints the norm of the components is equivalent to the variance that they explain. This is not correct, as shown in Section \ref{sec:vexppca}.

\citet{zou} reduce the minimization of the least squares criterion to the maximisation of the norm of the sparse components by adding an extra constraint to the PCA problem.
The least squares PCA problem for a generic component $\bt = \bX\ba$ can be defined as:
\begin{equation}\label{eq:spcalsprob}
(\ba, \bp) = \argmin_{\ba, \bp} || \bX - \bX\ba\bp\trasp||^2
\end{equation}
where $\bp$ is a vector of coefficients. From the first order conditions it follows that $\bp\trasp = (\ba\trasp\bX\trasp\bX\ba)^{-1}\ba\trasp\bX\trasp\bX$. \citet{zou} require that $||\bp||^2 = \bp\trasp\bp = 1$, which implies that $\bp = \ba$ (so that $\bp\trasp\ba = (\ba\trasp\bX\trasp\bX\ba)^{-1}\ba\trasp\bX\trasp\bX\ba = 1$) and that $||\ba|| = 1$. Since $||\bp\trasp\bp||^2 = (\ba\trasp\bS\bS\ba)/(\ba\trasp\bS\ba)^{-2}$, by the Cauchy-Schwarz inequality the condition $||\bp||^2 = 1$ is satisfied if and only if $\ba \propto \bS\ba$, that is, if and only if $\ba$ is an eigenvector of $\bS$. This condition is satisfied by the unconstrained PCs but not by the sparse components because when $\ba$ is sparse $\bS\ba \neq \ba\lambda$ and $\bp = \bS\ba(\ba\trasp\bS\ba)^{-1}\neq \ba$.

The consequence of adding the constraints $\bp = \ba$ and $\ba\trasp\ba = 1$ to the least squares criterion in Equation (\ref{eq:spcalsprob}) is that it becomes the maximisation of the norm of the sparse component $\dX\da$ in Equation (\ref{eq:spcaprob}). To see this, assume without loss of generality that the first $c$ variables form the block $\dX$ and that the block $\cX$ contains the remaining variables, so that $\bX = (\dX, \cX)$. Analogously we write $\ba = (\da\trasp, \mathbf{0}\trasp)\trasp$, where $\mathbf{0}$ is a vector of zeroes of dimension $p-c$. Then it follows that
\begin{eqnarray}\label{zou2}
    &&	\argmin_{\ba} || \bX - \bX\ba\ba\trasp||^2 = \argmin_{\da}|| (\dX, \cX) - \dX\da(\da\trasp, \mathbf{0}\trasp) ||^2  =\nonumber\\
%   && \argmin_{\da } tr(\cX\trasp\cX) + tr(\dX\trasp\dX) -2 tr(\dX\trasp\dX\da\da\trasp) + tr(\da\da\trasp\dX\trasp\dX\da\da\trasp) =\nonumber\\
       &&\argmin_{\da} || \dX - \dX\da\da\trasp||^2 = \argmax_{\da}\da\trasp\dX\trasp\dX\da\\
    &&subject\; to\;  \da\trasp\da = 1. \nonumber
\end{eqnarray}
The same argument is applicable to \citet{she}'s regularised SVD derivation of SPCA.

The maximisation of the norm of the sparse components does not maximise the variance explained by the components and also leads to results that are incompatible with the aim of finding the best sparse approximations of the PCs. %For example a component could be combination of two or more perfectly collinear variables or of a number of variables larger than the rank of the matrix.
This can be can be illustrated by considering the maximisation of the norm of the components on rank deficient data. Such data arise when some of the variables are collinear and when there are fewer observations than variables, which is often the case in gene expression and near infrared data sets.

When $rank(\bX) = r < min\{n, p\}$  $p - r$ columns of $\bX$ are linearly dependent and any linear combination of the $p$ columns can be expressed in terms of a linearly independent set of $r$ of its columns. Assume without loss of generality that the first $r$ columns of $\bX$ are linearly independent and denote them as $\dX$. Also let $\cX$ be the remaining columns. Then it is simple to prove that
\begin{equation}\label{eq:redpc}
\bX = [\dX, \cX] = \dX \left[\bI_r, \dX\trasp(\dX\dX\trasp)^{-}\cX\right]= \dX \left[\bI_r, (\dX\trasp\dX)^{-1}\dX\trasp\cX\right]= \bP_{\dX}\bX = \dX\bB,
\end{equation}
where the superscript `$^-$' denotes the Moore-Penrose generalised inverse and $\bB = \left[\bI_r, (\dX\trasp\dX)^{-1}\dX\trasp\cX\right]$. Therefore the PCs can be equivalently written as $\bX\bv \equiv \dX(\bB\bv) = \dX\dot\bv$, with $\dot\bv = \bB\bv$. Since the largest singular value of $\bX$ is not smaller than the largest singular value of $\dX$, $\dot\bv\trasp\dot\bv \geq 1$. Therefore, the PCs can be equivalently written as sparse components of cardinality $r$ defined by loadings with norm larger than one. However, conventional SPCA methods are prevented from identifying the sparse expression of the PCs because they consider only components with unit norm loadings and the only components with unit norm loadings that are equal to the PCs are the full cardinality PCs themselves. As a consequence, these methods would compute components with a number of variables larger than $r$. This fact is well documented by several examples available in the SPCA literature. It should be stressed that adding more than $r$ variables to the solution means including redundant variables that are perfectly correlated with the ones already selected, hence unnecessarily increasing its cardinality.

In LS SPCA a component of cardinality equal to the rank of the matrix will be recognised as optimal because by Equation (\ref{eq:lspcapcsproj}) the LS SPCA components are the PCs of $\bP_{\dX}\bX$ which, when $rank(\dX) = r$, is equal to $\bX$ by Equation (\ref{eq:redpc}). As a matter of fact, LS SPCA components of cardinality larger than the rank of the data matrix cannot be computed because in that case the matrix $\dX\trasp\dX$ would be singular.

As an example of inconsistent results given by conventional SPCA methods we consider a matrix with 100 observations on five perfectly collinear variables defined as
\begin{equation}
x_{ij} = (-1)^{i}\sqrt{j},\, i=1, \ldots, 100;\, j = 1,\ldots, 5.
\end{equation}
%The covariance matrix for these variable is equal to
%$$
%    \bS =
%    \begin{pmatrix}
%        100.00& 141.42& 173.21& 200.00& 223.61\\
%        141.42& 200.00& 244.95& 282.84& 316.23\\
%        173.21& 244.95& 300.00& 346.41& 387.30\\
%        200.00& 282.84& 346.41& 400.00& 447.21\\
%        223.61& 316.23& 387.30& 447.21& 500.00
%    \end{pmatrix}
%$$
The covariance matrix for these variables, $\bS = \bX\trasp\bX$, has rank one, its largest eigenvalues is equal to 1500 and the others are all zero. The first PC explains all the variance and can be written in terms of any of the variables as $\bx_j\sqrt{1500/s_{jj}},\, j = 1, \ldots, 5$, that is as a cardinality one component with loading larger than one. %by multiplying the variable by a weight larger than one. %It should be noted that the PC can equally be written as a component of cardinality less than five with loadings of norm larger than one.

The complete set of conventional SPCA solutions of different cardinalities are shown in Table \ref{tab:expl}. These indicate that $x_5$ "explains" just 33\% of the norm of the PC, $x_4$ and $x_5$ together only 60\% (with a net increase equal to 27\%) and so on. These results mislead the interpretation because they imply that some variables (the ones with larger variance) explain more variance than other collinear ones and that only the full cardinality component (the PC) explains the maximum variance. In reality, any variable explains the same variance as the PC but this cannot be revealed by conventional SPCA methods because the norm of a component with unit norm loadings is bounded by $trace(\dX\trasp\dX)$.
\renewcommand\baselinestretch{1} %for jasa{1.425}
\begin{table}[H]
 \centering
  \caption{Optimal conventional SPCA solutions.}
% Table generated by Excel2LaTeX from sheet 'Sheet1'
 \scriptsize{
\begin{tabular}{lrrrrr}
\toprule
      & \multicolumn{5}{c}{cardinality} \\
\cmidrule{2-6}
 variable     & 1     & 2     & 3     & 4     & 5 \\
\midrule
$x_1$ & 0     & 0     & 0     & 0     & 0.58 \\
$x_2$ & 0     & 0     & 0     & 0.60 & 0.52 \\
$x_3$ & 0     & 0     & 0.65 & 0.53 & 0.45 \\
$x_4$ & 0     & 0.75 & 0.58 & 0.46 & 0.37 \\
$x_5$ & 1     & 0.67 & 0.50 & 0.38 & 0.26 \\
\midrule
norm  & 500   & 900   & 1200  & 1400  & 1500 \\
rel norm & 0.33  & 0.60  & 0.80  & 0.93  & 1.0 \\
\bottomrule
\end{tabular}%
}
\label{tab:expl}%
\end{table}
\spacingset{\bsln}
%\renewcommand\baselinestretch{1.5} %for jasa{1.425}
The results of applying conventional SPCA to the correlation matrix of the variables in the above example follow a similar pattern.
The correlation matrix is a $5 \times 5$ matrix of ones with only one non-zero eigenvalue equal to 5. The results of conventional SPCA would be as shown in Table \ref{tab:expl2} irrespective of which variables are included because the variables are identical. These results lead to the meaningless conclusion that linear combinations of two or more identical variables explain more variance than a linear combination of fewer of them but less than the PC, which is collinear to them.
\renewcommand\baselinestretch{1} %for jasa{1.425}
\begin{table}[H]
 \centering
  \caption{Optimal conventional SPCA solutions for the correlation matrix.}
% Table generated by Excel2LaTeX from sheet 'Sheet1'
\scriptsize{
\begin{tabular}{lrrrrr}
\toprule
cardinality     & 1     & 2     & 3     & 4     & 5 \\
\cmidrule{2-6}
norm  & 1     & 2     & 3     & 4     & 5 \\
rel norm & 0.2  & 0.4  & 0.6  & 0.8  & 1.0 \\
\bottomrule
\end{tabular}%
}
\label{tab:expl2}%
\end{table}
%\renewcommand\baselinestretch{1.5} %for jasa{1.425}
\spacingset{\bsln}
Applying LS SPCA to this dataset would yield a single cardinality one component both for the covariance and the correlation matrices, as required. This can be seen by considering that the variance explained by any variable $x_j$ is equal to
$$
vexp(\bx_j) = \frac{\bx_j\trasp\bX\bX\trasp\bx_j}{\bx_j\trasp\bx_j} =
\frac{\sum_{i=1}^5 s_{ij}^2}{s_{jj}} =
\sum_{i=1}^5 s_{ii} = tr(S), %\frac{s_{jj}\sum_{i=1}^5 s_{ii}}{s_{jj}} =
$$
because $cor(x_i, x_j)^2 = 1 = s_{ij}^2/(s_{ii}s_{jj})$. Which variable is chosen depends on the algorithm used. An analogous example for blocks of collinear variables can be derived using the artificial data example introduced in \citet{zou}, as illustrated in \citet{mer}.

The case of collinear variables above is useful also to illustrate the fundamental difference between conventional SPCA and LS SPCA. The maximisation of the norm identifies blocks of highly correlated variables, which give similar information;  the maximisation of the variance explained selects less correlated variables which are important for explaining the variance also of the other variables in the set. For this reason, LS SPCA yield lower cardinality components that are better approximations of the PCS. \citet{mer} also discusses other unconvincing aspects of conventional SPCA methods.
\section{Computational details}\label{sec:comp}%
The basic projection SPCA algorithm is outlined in Algorithm \ref{algo:fspca}. The
algorithm is straightforward and simple to implement. Its complexity depends on how the variables for each block are selected (line \ref{algo:fwdselect}) and the eigenvectors are computed (line
\ref{algo:pwm}).

Determining the variables that explain a given percentage of the current PC is a problem analogous to model selection in multiple regression and feature selection in machine learning. These problems can be solved with different greedy approaches and algorithms. The question of which algorithm performs best is a long-standing problem which is difficult to answer because of the existence of competing desired objectives.
%Some of these algorithms are more computationally efficient than others, and some promise to yield more parsimonious solutions than others.
A recent literature review can be found in \citet{kum}. Many of these algorithms can be used for selecting the variables for computing projection SPCA.

In our implementation we chose the simple forward selection used in multiple
regression analysis for maximising the coefficient of determination because it is computationally very efficient.
%as it is well known that its computational complexity is $O(p)$.
In this algorithm the variables that sequentially explain the most extra variance of the
current PC enter the solution until a predetermined threshold for the coefficient of
determination is reached. This method can be seen as a QR decomposition with
supervised pivoting and can be implemented efficiently using updating formulae \citep[e.g.
see][]{bjo}.

Since for each component only the first eigenvector of the deflated covariance matrix
$\bQ_j\trasp\bQ_j$ is required, these can be computed economically using the power
method. The power method computes approximated results iteratively and is has complexity about $O(p^2)$, while algorithms that compute the whole set of eigenvectors accurately are about $O(p^3)$. The power method is used in extremely high dimensional problems such as, for example, Google page ranking \citep{bry}. The algorithm can be modified to compute LS SPCA components by changing line \ref{algo:pspca} accordingly and modifying the power method for the computation of generalised eigenvectors.

The theoretical time complexity of the algorithm cannot be computed exactly
%because the time taken by the forward selection and power method depends on the structure of
%the data. Since the complexity of all the operations is less than $O(p^3)$
but we expect the complexity of the whole algorithm to be less than $O(p^3)$ because it does not contain operations of this complexity or higher. In the next section we will analyse the time taken on the datasets empirically.
%\newcommand{\subki}{{_{ind_j[c]}}}
\renewcommand\baselinestretch{1.0} %for jasa{1.425}
%\algsetup{indent=2em}
\begin{algorithm}[H]\caption{Projection SPCA algorithm}
\footnotesize{
    \begin{algorithmic}[1]
    \Procedure {pspca}{$\bX,\, \alpha,\, stopRuleCompute$}
    \State \textbf{initialize}
        \State {\hspace{1em}$\bQ \gets \bX$}
        \State{\hspace{1em}$j \gets 0$}
        \State{\hspace{1em}stopCompute $\gets$ FALSE}
        %\Comment{stop rule for computation of components}
    \State \textbf{end initialize}
    %\Statex
    \While{stopCompute = FALSE}\Comment{\textbf{start components computation}}
        \State{$j \gets j + 1$}
        \State{$\bu_j :\, \bQ\bQ\trasp\bu_j = \bu_j\lambda_j$} \Comment{compute first
        PC of $\bQ$}\label{algo:pwm}
        \State{$ind_j \gets \{i_1,\ldots, i_{c_j}\}:\, ||\hat\bu_j||^2 \geq
        \alpha\lambda_j$} \Comment{variables selection output}\label{algo:fwdselect}
        %\Statex
%        \If{computeLSSPCA = FALSE}
        \State{$\ba_j \gets (\dX_j\trasp\dX_j)^{-1}\dX_j\trasp\bu_j$}
        \Comment{$j$-th sparse loadings, $\dX_j$ are columns of $\bX$ in
        $ind_j$}\label{algo:pspca}
%        \Else
%            \State{$\ba_j:\, \dX_j\trasp\bQ\bQ\trasp\dX_j\ba_j =
%\dX_j\trasp\dX_j\ba_j\gamma_{\max}$} \Comment{$j$-th sparse LS SPCA
%loadings}\label{algo:lsspca}
%        \EndIf
        \State{$ \bt_j \gets \dX_j\ba_j$} \Comment{$j$-th sparse component}
        \If{stopRuleCompute = FALSE}\quad\Comment{stop rule on total variance explained or
        number of components}
            \State{$ \bQ \gets \bQ - \frac{\bt_j\bt_j\trasp\bQ}{\bt_j\trasp\bt_j}$}
            \Comment{deflate $\bX$ of current component}\label{alg:deflX}
        \Else
            \State stopCompute $\gets$ TRUE \Comment{\textbf{terminate components
            computation}}
        \EndIf
    \EndWhile%
    \EndProcedure
    \end{algorithmic}}\label{algo:fspca}
\end{algorithm}%
\spacingset{\bsln}
%\renewcommand\baselinestretch{1.5} %for jasa{1.425}
\section{Numerical examples}\label{sec:examples}
In this section we will give an example of data analysis using projection SPCA and
illustrate the computational times and sparse reduction attained on large datasets.
For simplicity of exposition, we refer to the size of the loadings meaning the
absolute value of the non-zero ones.
\subsection{Sparse principal components analysis of crime data}
As an illustration of sparse principal component analysis we compare the results of
applying simple thresholding and projection SPCA to the correlation matrix of a dataset containing
socioeconomic measures collected in the 1990s for different US cities (available at
the UCI repository, \citealp{uci}).%\footnote{Communities and Crime Data Set,
%\url{https://archive.ics.uci.edu/ml/datasets/Communities+and+Crime}
The data were originally used to explain the rate of violent crime in the different
cities. We deleted 23 variables with missing values. The final set contains 1994
observations on 99 variables. The example is given only to illustrate an application
of sparse principal components analysis and does not claim to be a deep analysis of
these data, which is beyond the scope of this paper and our expertise.

The ordered loadings of the first and second PCs scaled to percentage
\emph{contributions} — that is, scaled to unit sum of absolute values — are shown in
Figure \ref{fig:crimepc1loads}. Clearly, it is difficult to choose a threshold value
to separate significant and nonsignificant loadings as they decrease at a regular
rate.
\begin{figure}[H]
  \centering
  % Requires \usepackage{graphicx}
  \includegraphics[width=0.95\textwidth, height =
  0.4\textwidth]{crime_ContrPC1_2_barplot.pdf}\\
  \caption{Ordered contributions of the first two PCs of the crime dataset.}
  \label{fig:crimepc1loads}
\end{figure}%
The summary results of projection SPCA with $\alpha = 95\%$ are shown in Table
\ref{tab:crimesum}. In this table the relative cumulative  variance explained (rCvexp)
denotes the ratio of the cumulative variance explained by the components over that
explained by the PCs. This is, it is computed as
\begin{equation*}
\mathrm{rCvexp(\hat\bu_1, \ldots, \hat\bu_c)} = \frac{vexp(\hat\bu_1, \ldots,
\hat\bu_c)} {vexp(\bu_1, \ldots, \bu_c)} =
\frac{\sum_{j = 1}^c evexp(\hat\bu_j)}{\sum_{j = 1}^c vexp(\bu_j)}.
\end{equation*}

As expected, rCvexp is always higher than 95\% and the cardinality of the sparse
components is quite low. The first two sparse components are highly correlated with
the corresponding PCs while the correlation decreases for successive components. The
$R^2$ statistics refer to the regression of the logarithm of crime rate on the
components computed thus far.
\renewcommand\baselinestretch{1} %for jasa{1.425}
\begin{table}[H]
 \centering
  \caption{Summary statistics for the first five projection SPCA components.}
 {\scriptsize
% Table generated by Excel2LaTeX from sheet 'Sheet3'
\begin{tabular}{lrrrrr}
\toprule
      & Comp1 & Comp2 & Comp3 & Comp4 & Comp5 \\
\midrule
\% Cumulative vexp & 24.4& 40.8& 49.8& 57.2& 62.7 \\
\% Relative cumulative vexp & 96.5  & 96.5  & 96.5  & 96.6  & 96.6 \\
Cardinality & 3     & 5     & 7     & 9     & 8 \\
Correlation with PC & 0.97  & 0.98  & 0.94 & 0.83 & 0.92 \\
\% $R^2$ log(crime rate) on PCs & 44.0 & 49.0 & 50.9 & 51.3 & 53.3 \\
\% $R^2$ log(crime rate) on sparse comp. & 45.7 & 51.9 & 56.4 & 57.2 & 57.2 \\
\bottomrule
\end{tabular}%
}
\label{tab:crimesum}%
\end{table}
\spacingset{\bsln}
Table \ref{tab:crimePC1} shows the ten largest contributions of the first PC and the
three non-zero contributions of the first projection sparse component. The last column
shows the coefficients of determination resulting from regressing each variable on all
other variables in the set (squared multiple correlations), $R_j^2$. %This is the same coefficient used to compute the variance inflation factor in regression; the closer the value to one, the more the variable is (multiply) correlated with the others.

The ten largest PC contributions assign about the same weight to different measures of
two indicators: income and the presence of both parents in a household. The
corresponding variables are highly correlated among themselves, as shown by the large
$\mathrm{R}^2_j$ values. The sparse contributions assign about 50\% of weight to
income, about 37\% to families with two parents, and 12\% to a new indicator: large
family size. This last variable is hardly correlated with the others and its
contribution to the PC ranks 46th in size. Combining just three variables out of 99,
the sparse component explains 96.5\% of the variance explained by the first PC,
indicates the key features of the PC and identifies another important variable which
would have been difficult to find by analysing the PC loadings.
\renewcommand\baselinestretch{1.1} %for jasa{1.425}
\begin{table}[H]
 \centering
  \caption{Largest ten contributions of the first full cardinality PC and those of the
  first projection SPCA component.}
 {\scriptsize
% Table generated by Excel2LaTeX from sheet 'Sheet2'
\begin{tabular}{rlr}
\toprule
Contribution & Variable & $R_j^2$ \\
\midrule
\multicolumn{1}{l}{PCA}   &       &  \\
\textbf{-2.2\%} & \textbf{median family income (differs from household income)} & 0.98
\\
-2.2\% & median household income  & 0.98 \\
\textbf{-2.1\%} & \textbf{percentage of children in family housing with two parents} &
0.98 \\
-2.1\% & percentage of households with investment/rent income in 1989  & 0.83 \\
2.1\% & percentage of people below the poverty level  & 0.84 \\
-2.1\% & percentage of families (with children) with two parents  & 0.98 \\
-2.0\% & percentage of children aged 4 and under in two-parent households  & 0.90 \\
-2.0\% & per capita income  & 0.92 \\
2.0\% & percentage of households with public assistance income in 1989  & 0.75 \\
2.0\% & percentage of occupied housing units without a phone (in 1990, this was rare!)
& 0.76 \\
\multicolumn{1}{l}{Proj SPCA} &       &  \\
-51\% & median family income (differs from household income)  & 0.51 \\
-37\% & percentage of children in family housing with two parents  & 0.52 \\
12\%  & \textbf{percentage of family households that are large (6 or more)} & 0.09 \\
\bottomrule
\end{tabular}%
}
\label{tab:crimePC1}%
\end{table}
\spacingset{\bsln}
%\renewcommand\baselinestretch{1.5} %for jasa{1.425}
The ten largest contributions of the second PC, shown in Table \ref{tab:crimePC2},
assign about 92\% of the weights (rescaled) to highly correlated variables regarding
immigrants: the percentage of immigrants by time of arrival and skills in English
language, and 8\% to living in dense housing conditions. Instead, the five sparse
contributions, also shown in Table \ref{tab:crimePC2}, assign  52\% of the weights to
one variable regarding immigrants, while the remaining weights are evenly divided into
four variables concerning different social aspects such as the percentage of older
residents and the percentage of people living in their own house. These variables have
low multiple correlation. The indicator for large families, which was also present in
the first component, has high correlation with the indicator for dense housing, which
is an important contributor to the PC. The positive contribution of the percentage of
people living in their own house is probably due to the cost of paying a mortgage over
household income and to the fact that older residents are more likely to own their
house.
\renewcommand\baselinestretch{1} %for jasa{1.425}
\begin{table}[H]
 \centering
  \caption{Largest ten contributions of the second full cardinality PC and those of
  the second sparse components computed with projection SPCA.}
 {\scriptsize
% Table generated by Excel2LaTeX from sheet 'Sheet2'
\begin{tabular}{rlr}
\toprule
Contribution & Variable  & $R_j^2$ \\
\midrule
\multicolumn{1}{l}{PCA}&       &  \\
\textbf{2.7\%} & \textbf{percentage of  population who immigrated within the last 10
years} & 1.00 \\
2.7\% & percentage of  population  who immigrated within the last 8 years  & 1.00 \\
2.6\% & percentage of  population  who immigrated within the last 5 years  & 1.00 \\
2.6\% & percentage of  population  who immigrated within the last 3 years  & 0.98 \\
2.6\% & percentage of foreign-born people & 0.96 \\
-2.3\% & percentage of people who speak only English  & 0.95 \\
2.3\% & percentage of people who do not speak English well  & 0.94 \\
\textbf{2.1\%} & \textbf{percentage of people in dense housing (more than 1 person per
room)} & \textbf{0.86} \\
2.0\% & percentage of population of Asian heritage  & 0.63 \\
2.0\% & percentage of population of Hispanic heritage  & 0.90 \\
\multicolumn{1}{l}{Proj SPCA} &       &  \\
42\%  & percentage of  population who immigrated within the last 10 years  & 0.11 \\
\textbf{-15\%} & \textbf{percentage of population aged 65 and over}  & 0.06 \\
\textbf{15\%}  & \textbf{owner-occupied housing - upper quartile value}  & 0.55 \\
\textbf{14\%}  & \textbf{percentage of family households that are large (6 or more)}
& 0.47 \\
\textbf{13\%} & \textbf{number of people living in areas classified as urban}  & 0.32
\\
\bottomrule
\end{tabular}%
}
\label{tab:crimePC2}%
\end{table}
\spacingset{\bsln}
%\renewcommand\baselinestretch{1.5} %for jasa{1.425}
In summary, in this example the sparse components combining few key variables explain
over 96\% of the variance explained by the PCs. In contrast, large PC loadings
correspond to clusters of mutually highly correlated variables which shadow other
important but less correlated variables.

The sparse components computed with LS SPCA are often better at predicting variables not
included in the dataset than the PCs. This probably depends on the sparse components
giving higher weight to variables that are less correlated with others in the set, but
which may be more correlated with an exogenous response. This is not guaranteed for
all applications but, in this example, the sparse components yield higher coefficients
of determination for the logarithms of crime rate than the PCs, as shown in Table
\ref{tab:crimesum}. These coefficients of determination are quite high considering
that the coefficient of determination for regressing the logarithm of crime rate on
all 99 variables is 69.1\%.
\subsection{Computational results for different datasets}
In this section we report the results of running 20 repetitions of projection SPCA on
datasets with different numbers of variables, from less than 100 to over 16,000. For
all datasets we required that each sparse component explained at least 95\% of the
corresponding PC. The computational times reported were measured using an
implementation of the algorithm in \texttt{C++} embedded in \texttt{R} using the
packages Rcpp and  RcppEigen.

It is well known that \texttt{R} is an inefficient language \citep{mor}, so this
choice makes the program run slower than if it were written only in \texttt{C++} but
has other advantages. Therefore, the times shown are only indicative of the
performance of the algorithms on large matrices, which can be improved if necessary.

The datasets that we consider in this section, listed in Table \ref{tab:descdata},
have been taken from various sources, mostly from the data distributed with the
Elements of Statistical Learning book (ESL) \citep{has}
website.
Other sets were taken from the UCI Machine Learning Repository \citep{uci}. The remaining sets
were taken from different sources, see Table \ref{tab:descdata} for details.
Most of the datasets contain gene expression data. The characteristic of these sets is
to have a large number of features and fewer objects.
The largest dataset \citep{rama} has been used to test other SPCA methods \citep{zou,
sri, wan}.
\renewcommand\baselinestretch{1} %for jasa{1.425}
\begin{table}[H]
  \caption{Description of the datasets used for numerical comparison.}
\begin{threeparttable}
 \centering
 {\scriptsize
% Table generated by Excel2LaTeX from sheet 'Sheet1'
\begin{tabular}{lrrll}
\toprule
Name  & Samples &Features & Description & Source \\
\midrule
Crime& 1994& 99    & social data & UCI Repository\tnote{a} \\
Isolet & 6238& 716   & character recognition & UCI Repository\tnote{b} \\
Ross NCI60 &60  & 1375  & gene expression & \texttt{R} package made4 \tnote{c}  \\
Khanh &88& 2318  & gene expression & ESL \\
Phoneme &256    & 4509  & speech recognition & ESL \\
NCI60 &60   & 6830  & gene expression & ESL \\
Protein &11    & 7466  & protein cytometry& ESL \\
Radiation &58    & 12625 & gene expression& ESL \\
Ramaswamy &144    & 16063 & gene expression & Broadinstitute repository\tnote{d} \\
\bottomrule
\end{tabular}%
\begin{tablenotes}
\item[a] \url{https://archive.ics.uci.edu/ml/datasets/Communities+and+Crime}
\item[b] \url{https://archive.ics.uci.edu/ml/datasets/ISOLET}
\item[c] \url{http://bioconductor.org/packages/release/bioc/html/made4.html}
\item[d]
    \url{http://software.broadinstitute.org/cancer/software/genepattern/datasets}
\end{tablenotes}
}
\end{threeparttable}
\label{tab:descdata}%
\end{table}
\spacingset{\bsln}
%\renewcommand\baselinestretch{1.5} %for jasa{1.425}

Tables \ref{tab:times1} and \ref{tab:times2} show the median computational time, over
20 repetitions, taken to compute ten sparse components for nine datasets considered,
together with the cardinalities. Table \ref{tab:times1} also shows the relative
cumulative variance explained, which is always larger than 95\% as required. The same
is true for larger datasets (not shown). The cardinalities are also very low for the
largest datasets. Such dramatic reduction may be partly due to the fact that these
sets have fewer observations than variables. The times taken to compute the first
component vary from four milliseconds to 2.7 seconds. Computing ten components for the
largest set takes less than one and a half minutes.

Computing the theoretical time complexity of the projection SPCA algorithm is not
possible because the convergence times for both the forward selection and power method
algorithms depend on the data structure. Since the power method is at least $O(n^2)$
and there are no operations $O(n^3)$, we expect the computational complexity to be
within these two limits. The time taken to compute additional components grows almost
linearly for fixed dataset size,  and the time taken to compute a fixed number of components increases non-linearly with the size of the datasets.%, as shown in Figure \ref{fig:timeby}.

%\begin{figure}[H]
\renewcommand\baselinestretch{1} %for jasa{1.425}
\begin{table}[H]
 \centering
  \caption{Median computational times, relative cumulative variance explained, and
  cardinalities of the first 10 projection SPCA components with $\alpha = 0.95$ for
  three medium-size datasets. Time is expressed in milliseconds.}
 {\scriptsize
% Table generated by Excel2LaTeX from sheet 'Tables'
\begin{tabular}{lrrrrrrrrrrrr}
\toprule
Name&&\multicolumn{3}{c}{Isolets}&     &\multicolumn{3}{c}{Ross NCI60}&
&\multicolumn{3}{c}{Khanh}\\
No. vars&&\multicolumn{3}{c}{716}&     &\multicolumn{3}{c}{1375}&
&\multicolumn{3}{c}{2318}\\
\midrule
%\cmidrule{2-4}\cmidrule{6-8}\cmidrule{10-12}
No. Components&&Time&rCvexp&card &     &Time&rCvexp&card &     &Time& rCvexp&card \\
\cmidrule{1-1}\cmidrule{3-5}\cmidrule{7-9}\cmidrule{11-13}
1 Comp&&4.0  &95.4\%&5&     &24.1&96.8\%&5&     &126.2&96.4\%&6 \\
2 Comps&&15.2 &95.4\%&17&     &65.1&96.6\%&6&     &260.6&96.2\%&6 \\
3 Comps&&26.9 &95.5\%&15&     &102.9&96.7\%&6&     &360.8&96.4\%&4 \\
4 Comps&&45.6 &95.5\%&26&     &155.0  &96.7\%&9&     &555.6&96.4\%&8 \\
5 Comps&&63.34&95.5\%&27&     &219.6&96.7\%&9&     &731.3&96.5\%&9 \\
6 Comps&&80.71&95.6\%&25&     &281.4&96.8\%&9&     &969.8&96.5\%&8 \\
7 Comps&&103.0  &95.6\%&32&     &343.1&96.8\%&11&     &1167.0&96.5\%&10 \\
8 Comps&&125.5&95.8\%&36&     &418.1&96.9\%&11&     &1447.0&96.6\%&12 \\
9 Comps&&148.7&95.9\%&32&     &499.5&96.9\%&12&     &1640.0&96.6\%&10 \\
10 Comps&&175.2&95.9\%&42&     &580.0&97.0\%&13&     &1841.0&96.7\%&12\\
\bottomrule
\end{tabular}%
}
\label{tab:times1}%
\end{table}
\newcolumntype{C}{@{\extracolsep{0.0cm}}l@{\extracolsep{0.25cm}}}%

\begin{table}[H]
 \centering
  \caption{Median computational times and cardinalities of the first 10 projection
  SPCA components
  with $\alpha = 0.95$ for five large-size datasets. Time is expressed in seconds.}
 {\scriptsize
% Table generated by Excel2LaTeX from sheet 'Tables'
\begin{tabular}{Crrrrrrrrrrrrrr}
\toprule
Name&\multicolumn{2}{c}{Phoneme}&     &\multicolumn{2}{c}{NCI60}&
&\multicolumn{2}{c}{Protein}&       &\multicolumn{2}{c}{Radiation}&
&\multicolumn{2}{c}{Ramaswamy}  \\
No. vars&\multicolumn{2}{c}{4509}&     &\multicolumn{2}{c}{6830}&
&\multicolumn{2}{c}{7466}&       &\multicolumn{2}{c}{12625}&
&\multicolumn{2}{c}{16063}  \\
\midrule
%\cmidrule{2-3}\cmidrule{5-6}\cmidrule{8-9}\cmidrule{11-12}\cmidrule{14-15}
No. comp.& Time  & card  &       & Time  & card  &       & Time  & card  &       &
Time  & card  &       & Time & card \\
\cmidrule{1-1}\cmidrule{2-3}\cmidrule{5-6}\cmidrule{8-9}\cmidrule{11-12}\cmidrule{14-15}
1 Comp & 0.1   &  1    &       & 0.6   & 4     &       & 0.4   & 1     &       & 2.7
& 5     &       & 2.3   & 3 \\
2 Comps & 0.4   & 4    &       & 1.9   & 5     &       & 1.1   & 2     &       & 6.1
& 5     &       & 6.9   & 4 \\
3 Comps & 1.1   & 13    &       & 2.9   & 6     &       & 1.8   & 1     &       & 12.8
& 7     &       & 12.5  & 6 \\
4 Comps & 2.0   & 16    &       & 4.5   & 8     &       & 2.4   & 2     &       & 17.7
& 7     &       & 19.5  & 8 \\
5 Comps & 3.3   & 20    &       & 6.7   & 10    &       & 3.3   & 2     &       & 23.7
& 7     &       & 28.0  & 9 \\
6 Comps & 4.4   & 23    &       & 8.2   & 9     &       & 4.2   & 3     &       & 30.3
& 8     &       & 37.6  & 8 \\
7 Comps & 5.8   & 27    &       & 10.0  & 8     &       & 5.6   & 3     &       & 37.3
& 8     &       & 46.2  & 10 \\
8 Comps & 7.0   & 28    &       & 12.0  & 10    &       & 6.5   & 3     &       & 44.6
& 9     &       & 60.4  & 13 \\
9 Comps & 8.7   & 27    &       & 14.2  & 10    &       & 7.4   & 4     &       & 52.0
& 9     &       & 73.0  & 17 \\
10 Comps & 10.6 & 26    &       & 16.4  & 10    &       & 8.2   & 4     &       & 59.0
& 8     &       & 85.6  & 16 \\
\bottomrule
\end{tabular}%
}
\label{tab:times2}%
\end{table}
\spacingset{\bsln}
To estimate the time complexity function for $c$ components and $p$ variables, $T(c,
p)$, we assume that it is polynomial in both parameters, which can be written as
\begin{equation*}
T(c, p) = kd^\alpha p^\beta\epsilon\; \mathrm{or}\; \log(T(c, p)) = \log(k) + \alpha
\log(c)  + \beta\log(p)+\eta.
\label{eq:polymod}\end{equation*}
The fit obtained with the data available is almost perfect with $\mathrm{R}^2 = 0.99$,
as shown in Table \ref{tab:timereg1}. As expected, the estimated coefficients show that computational time increases slightly faster than linearly with the number of components computed while it is approximately $O(p^2)$ with respect to the dataset size. It should be noted that the time taken for matrices with fewer objects than features can be reduced by using the inverse svd. We do not report these times because are peculiar to the matrices while the ones reported refer generically to matrices with a given number of variables.
% and Figure \ref{fig:timereg1}
\renewcommand\baselinestretch{1} %for jasa{1.425}
\begin{table}[H]
 \centering
  \caption{Output of the regression of $\log(T(c,p))$ on $\log(c)$ and $\log(p)$.}
 {\scriptsize
% Table generated by Excel2LaTeX from sheet 'Sheet5'
\begin{tabular}{lrrrrr}
\toprule
Coefficients & Estimate & Std. error & t-value & Pr($>|t|$) &  \\
$\log_{10}(k)$ & -5.06 & 0.08  & -61.44 & $<$ 2e-16 & *** \\
$\alpha$ & 1.46  & 0.05  & 30.81 & $<$ 2e-16 & *** \\
$\beta$ & 2.03  & 0.02  & 93.37 & $<$ 2e-16 & *** \\
%---   &       &       &       &       &  \\
%Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 &       &       &
&       &  \\
\midrule
\multicolumn{6}{l}{Residual standard error: 0.1362 on 87 degrees of freedom} \\
\multicolumn{6}{l}{Multiple R-squared:  0.9911,    Adjusted R-squared:  0.9909 } \\
\multicolumn{6}{l}{F-statistic:  4833 on 2 and 87 DF,  p-value: $<$ 2.2e-16} \\
\bottomrule
\end{tabular}%
}
\label{tab:timereg1}%
\end{table}
\spacingset{\bsln}
\subsection{Comparison with conventional SPCA}\label{sec:compspca}
In this section we compare the performance of the first sparse components computed with a conventional SPCA method and with projection SPCA. As conventional SPCA method we used SPCA-IE \citep{wan} with \emph{amvl} criterion. This method was shown to perform similarly to other SPCA methods, it does not require to choose arbitrary sparsity parameters and is simple to implement.

We compare the norm of the components, the relative cumulative variance that they explain and their correlation with the full PC for different cardinalities. Since the results for large matrices are quite similar we limit this comparison to four datasets: Crime, Isolets, Khanh and Ramaswamy. For the last dataset we computed the conventional sparse components using simple thresholding, stopping the computation to cardinality 200; details of the performance of components with larger cardinality computed with different conventional SPCA methods can be found in the papers cited above.

Figure \ref{fig:compspca} compares the performance of the components computed with different methods. The LS SPCA values for the rank deficient Khanh and Ramaswamy datasets are available only until the solutions reach full rank cardinality (87 and 143, respectively) at which the components explain the maximum possible variance. Clearly SPCA-IE outperforms projection SPCA in the norm of the components. However, the latter method guarantees higher variance explained and closer convergence to the PC with much lower cardinality.
\begin{figure}[H]
  \centering
  % Requires \usepackage{graphicx}
  \includegraphics[width=0.95\textwidth, height = 0.5\textwidth]{comp_all_3x2.pdf}\\
  \caption{Norm, rCvexp and correlation with the PCs versus cardinality for the first sparse components computed with SPCA-IE (spca) and projection SPCA (lsspca).}
  \label{fig:compspca}
\end{figure}
This phenomenon is more evident for large rank deficient datasets when conventional sparse components combination of hundreds of variables explain less variance than projection components which are combinations of much fewer variables. As explained in Section \ref{sec:compspca}, projection SPCA components reach the maximum when their cardinality is equal to the matrix rank while conventional SPCA components reach the maximum when their cardinality is equal to the number of variables. The plots show clearly that the components' norms are not related to the variance that they explain or to their correlation with the PC, confirming the theoretical conclusions given above.

Table \ref{fig:compspca} shows the cardinality with which the components computed with the two methods reached 99.9\% rCvexp (or 0.998 correlation with the PC). For the Ramaswamy dataset the threshold for the variance explained is not reached with cardinality 200. The last column of the table shows rCvexp for the conventional SPCA components with the same cardinality at which the projection SPCA components reach 99.9\% rCvexp. Again, it can be appreciated the much greater efficiency of projection SPCA at approximating the PCs.
\renewcommand\baselinestretch{1.0} %for jasa{1.425}
\begin{table}[H]
 \centering
  \caption{Cardinality needed to reach 99.9\% rCvexp by the components computed with projection SPCA, second column, and SPCA-IE, third column. The fourth column shows the rCvexp for the SPCA-IE components with the same cardinality as the corresponding projection SPCA component.}
% Table generated by Excel2LaTeX from sheet 'Sheet1'
 \scriptsize{
\begin{tabular}{lrrrr}
\toprule
 Dataset& Rank     & Proj. SPCA & SPCA  & rCvexp SPCA \\
\midrule
Crime & 99& 38    & 74    & 85.200\% \\
Isolet & 716& 123   & 439   & 59.600\% \\
Khanh & 87& 28    & 1338  & 0.060\% \\
Ramaswamy &143&  23    & $>$ 200    & 0.002\% \\
\bottomrule
\end{tabular}%
}
\label{tab:compspca}%
\end{table}
\spacingset{\bsln}
%\renewcommand\baselinestretch{1.5} %for jasa{1.425}
\section{Discussion}\label{sec:conclusions}
Conventional SPCA method have been shown to give results similar to simple thresholding if not worse \citep{zou, wan}. Thresholding has been proven to give misleading results \citep[e.g.][]{cad}. Since the loadings are proportional to the covariances of the variables with the PCs, the largest loadings correspond to variables that are highly correlated with the current PC and among themselves. These sets of variables are not very informative because they contain different measures of the same features.

\citet{zou} propose three properties of a good SPCA method:
\begin{itemize}
\item Without any sparsity constraint, the method should reduce to PCA.
\item It should be computationally efficient for both small p and big p data.
\item It should avoid misidentifying the important variables.
\end{itemize}
In the light of what explained above, the first property is not enough for a good method. The second is not necessary for the most commonly analysed datasets and the third one is vague because importance is not defined and variables known to be unimportant could be directly eliminated from the analysis.

We suggest the following properties for a good SPCA method:
\begin{itemize}
\item Without any sparsity constraint the method should reduce to PCA.
\item It should identify the sparsest expression of the principal components.
\item The addition of a variable perfectly correlated with one or more variables already in the solution should not improve the objective function.
\end{itemize}
The last property eliminates redundant variables from the solution and should deter the inclusion of highly correlated ones. Conventional SPCA methods do not have the last two properties while methods based on LS SPCA do. It is surely possible that different methods would have these properties.

Projection SPCA methodology is intuitive and can be understood by users who do not have a deep knowledge of numerical optimisation. The only parameter to be set for computing the solutions is the proportion of variance explained, the meaning of which is also easily understandable. The algorithm is simple to implement and scalable to large datasets. Most conventional SPCA methods are based on complicated numerical optimisations and require setting values of parameters the effect of which is difficult to understand.

Different variable selection methods could be used for projection SPCA. Future research could explore the effect of using different selection methods, such as least angle or lasso regression, for example.

In this work we have developed a framework for computing LS SPCA solutions. These solutions give sparse components that closely approximate the full PCs with low cardinality. These can be efficiently computed also for very large datasets. We also show that conventional SPCA methods suffer from a number a drawbacks which make the solution provided less attractive than the corresponding LS SPCA solutions.

LS SPCA is implemented in the R package \emph{spca} available on GitHub. Projection SPCA will be added to this package in the future.
%One aspect of projection SPCA that can be improved is the variable selection method used to identify the blocks for projecting the PCs. Many different algorithms can be used. Forward selection has the advantage of being extremely fast, while backward selection is more computationally demanding, especially when the number of variables is large; it is also not stable in the presence of highly correlated variables. Stepwise and stagewise procedures can produce lower cardinality loadings but are also more computationally demanding than other methods. Methods based on analytical constraints such as the Lasso or regularised regression could also be adopted.to verify in which context these give good sparse solutions.

%\input{projection}
%\input{partCorr}
%\input{bPca_algo}


%\input{apaperJSS_02_model}
%\input{apaperJSS_04_BE}
%\input{r-jss//baseball_from_knitr}
%\input{apaperJSS_06_examples}
%\bibliographystyle{anzjs}
%\bibliographystyle{anzjs}
%\nolinenumbers
\section*{References}
\bibliographystyle{elsarticle-harv}
\begin{thebibliography}{18}
\expandafter\ifx\csname natexlab\endcsname\relax\def\natexlab#1{#1}\fi
\expandafter\ifx\csname url\endcsname\relax
  \def\url#1{\texttt{#1}}\fi
\expandafter\ifx\csname urlprefix\endcsname\relax\def\urlprefix{URL }\fi

\bibitem[{Bj\"{o}rck(2015)}]{bjo}
Bj\"{o}rck, {\AA}., 2015. Numerical Methods in Matrix Computations. Vol.~59 of
  Texts in Applied Mathematics. Springer International Publishing.

\bibitem[{Bryan and Leise(2006)}]{bry}
Bryan, K., Leise, T., 2006. The \$25,000,000,000 eigenvector: The linear
  algebra behind {Google}. SIAM REVIEW 48~(3), 569--581.

\bibitem[{Cadima and Jolliffe(1995)}]{cad}
Cadima, J., Jolliffe, I.~T., 1995. Loadings and correlations in the
  interpretation of principal components. Journal of Applied Statistics 22~(2),
  203--214.

\bibitem[{Frank and Asuncion(2010)}]{uci}
Frank, A., Asuncion, A., 2010. {UCI} machine learning repository. University of
  California, Irvine, School of Information and Computer Sciences, isolet
  dataset available at {https://archive.ics.uci.edu/ml/datasets/ISOLET}.

\bibitem[{Hastie et~al.(2009)Hastie, Tibshirani, and Friedman}]{has}
Hastie, T., Tibshirani, R., Friedman, J., 2009. The elements of statistical
  learning: data mining, inference and prediction, 2nd Edition. Springer.

\bibitem[{Hotelling(1933)}]{hot}
Hotelling, H., 1933. {Analysis of a Complex of Statistical Variables with
  Principal Components}. Journal of Educational Psychology 24, 498--520.

\bibitem[{Jeffers(1967)}]{jef}
Jeffers, J., 1967. Two case studies in the application of principal component.
  Applied Statistics 16, 225--236.

\bibitem[{Jolliffe and Uddin(2000)}]{jol00}
Jolliffe, I., Uddin, M., 2000. The simplified component technique: An
  alternative to rotated principal components. Journal of Computational and
  Graphical Statistics 9~(4), 689--710.

\bibitem[{Kumar and Minz(2014)}]{kum}
Kumar, V., Minz, S., 2014. Feature selection: {A} literature review. Smart {CR}
  4~(3), 211--229.

\bibitem[{Merola(2015)}]{mer}
Merola, G., 2015. Least squares sparse principal component analysis: a backward
  elimination approach to attain large loadings. Australia \& New Zealand
  Journal of Statistics 57, 391--429.

\bibitem[{Moghaddam et~al.(2006)Moghaddam, Weiss, and Avidan}]{mog}
Moghaddam, B., Weiss, Y., Avidan, S., 2006. Spectral bounds for sparse pca:
  Exact and greedy algorithms. In: Advances in Neural Information Processing
  Systems. MIT Press, pp. 915--922.

\bibitem[{Morandat et~al.(2012)Morandat, Hill, Osvald, and Vitek}]{mor}
Morandat, F., Hill, B., Osvald, L., Vitek, J., 2012. Evaluating the design of
  the r language: Objects and functions for data analysis. In: Proceedings of
  the 26th European Conference on Object-Oriented Programming. ECOOP'12.
  Springer-Verlag, Berlin, Heidelberg, pp. 104--131.

\bibitem[{Pearson(1901)}]{pea}
Pearson, K., 1901. {On lines and planes of closest fit to systems of points in
  space}. Philosophical Magazine 2~(6), 559--572.

\bibitem[{Ramaswamy et~al.(2001)Ramaswamy, Tamayo, Rifkin, Mukherjee, Yeang,
  Angelo, Ladd, Reich, Latulippe, Mesirov, Poggio, Gerald, Loda, Lander, and
  Golub}]{rama}
Ramaswamy, S., Tamayo, P., Rifkin, R., Mukherjee, S., Yeang, C.~H., Angelo, M.,
  Ladd, C., Reich, M., Latulippe, E., Mesirov, J.~P., Poggio, T., Gerald, W.,
  Loda, M., Lander, E.~S., Golub, T.~R., dec 2001. {Multiclass cancer diagnosis
  using tumor gene expression signatures.} Proceedings of the National Academy
  of Sciences of the United States of America 98~(26), 15149--15154.

\bibitem[{Shen and Huang(2008)}]{she}
Shen, H., Huang, J., 2008. Sparse principal component analysis via regularized
  low rank matrix approximation. Journal of Multivariate Analysis 101,
  1015–1034.

\bibitem[{Sriperumbudur et~al.(2009)Sriperumbudur, Torres, and Lanckriet}]{sri}
Sriperumbudur, B., Torres, D., Lanckriet, G., 2009. A {D.C.} programming
  approach to the sparse generalized eigenvalue problem. Computer Engineering
  1~(1), 40.

\bibitem[{Wang and Wu(2012)}]{wan}
Wang, Y., Wu, Q., 2012. Sparse pca by iterative elimination algorithm. Advances
  in Computational Mathematics 36, 137--151.

\bibitem[{Zou et~al.(2006)Zou, Hastie, and Tibshirani}]{zou}
Zou, H., Hastie, T., Tibshirani, R., 2006. Sparse principal component analysis.
  Journal of Computational and Graphical Statistics 15~(2), 265--286.

\end{thebibliography}

\end{document}

1612.00939.extracted_bib
[{Bj\"{orck(2015)]{bjo Bj\"{orck, {\AA., 2015. Numerical Methods in Matrix Computations. Vol.~59 of Texts in Applied Mathematics. Springer International Publishing.
1612.00939.extracted_bib
[{Bryan and Leise(2006)]{bry Bryan, K., Leise, T., 2006. The \$25,000,000,000 eigenvector: The linear algebra behind {Google. SIAM REVIEW 48~(3), 569--581.
1612.00939.extracted_bib
[{Cadima and Jolliffe(1995)]{cad Cadima, J., Jolliffe, I.~T., 1995. Loadings and correlations in the interpretation of principal components. Journal of Applied Statistics 22~(2), 203--214.
1612.00939.extracted_bib
[{Frank and Asuncion(2010)]{uci Frank, A., Asuncion, A., 2010. {UCI machine learning repository. University of California, Irvine, School of Information and Computer Sciences, isolet dataset available at {https://archive.ics.uci.edu/ml/datasets/ISOLET.
1612.00939.extracted_bib
[{Hastie et~al.(2009)Hastie, Tibshirani, and Friedman]{has Hastie, T., Tibshirani, R., Friedman, J., 2009. The elements of statistical learning: data mining, inference and prediction, 2nd Edition. Springer.
1612.00939.extracted_bib
[{Hotelling(1933)]{hot Hotelling, H., 1933. {Analysis of a Complex of Statistical Variables with Principal Components. Journal of Educational Psychology 24, 498--520.
1612.00939.extracted_bib
[{Jeffers(1967)]{jef Jeffers, J., 1967. Two case studies in the application of principal component. Applied Statistics 16, 225--236.
1612.00939.extracted_bib
[{Jolliffe and Uddin(2000)]{jol00 Jolliffe, I., Uddin, M., 2000. The simplified component technique: An alternative to rotated principal components. Journal of Computational and Graphical Statistics 9~(4), 689--710.
1612.00939.extracted_bib
[{Kumar and Minz(2014)]{kum Kumar, V., Minz, S., 2014. Feature selection: {A literature review. Smart {CR 4~(3), 211--229.
1612.00939.extracted_bib
[{Merola(2015)]{mer Merola, G., 2015. Least squares sparse principal component analysis: a backward elimination approach to attain large loadings. Australia \& New Zealand Journal of Statistics 57, 391--429.
1612.00939.extracted_bib
[{Moghaddam et~al.(2006)Moghaddam, Weiss, and Avidan]{mog Moghaddam, B., Weiss, Y., Avidan, S., 2006. Spectral bounds for sparse pca: Exact and greedy algorithms. In: Advances in Neural Information Processing Systems. MIT Press, pp. 915--922.
1612.00939.extracted_bib
[{Morandat et~al.(2012)Morandat, Hill, Osvald, and Vitek]{mor Morandat, F., Hill, B., Osvald, L., Vitek, J., 2012. Evaluating the design of the r language: Objects and functions for data analysis. In: Proceedings of the 26th European Conference on Object-Oriented Programming. ECOOP'12. Springer-Verlag, Berlin, Heidelberg, pp. 104--131.
1612.00939.extracted_bib
[{Pearson(1901)]{pea Pearson, K., 1901. {On lines and planes of closest fit to systems of points in space. Philosophical Magazine 2~(6), 559--572.
1612.00939.extracted_bib
[{Ramaswamy et~al.(2001)Ramaswamy, Tamayo, Rifkin, Mukherjee, Yeang, Angelo, Ladd, Reich, Latulippe, Mesirov, Poggio, Gerald, Loda, Lander, and Golub]{rama Ramaswamy, S., Tamayo, P., Rifkin, R., Mukherjee, S., Yeang, C.~H., Angelo, M., Ladd, C., Reich, M., Latulippe, E., Mesirov, J.~P., Poggio, T., Gerald, W., Loda, M., Lander, E.~S., Golub, T.~R., dec 2001. {Multiclass cancer diagnosis using tumor gene expression signatures. Proceedings of the National Academy of Sciences of the United States of America 98~(26), 15149--15154.
1612.00939.extracted_bib
[{Shen and Huang(2008)]{she Shen, H., Huang, J., 2008. Sparse principal component analysis via regularized low rank matrix approximation. Journal of Multivariate Analysis 101, 1015–1034.
1612.00939.extracted_bib
[{Sriperumbudur et~al.(2009)Sriperumbudur, Torres, and Lanckriet]{sri Sriperumbudur, B., Torres, D., Lanckriet, G., 2009. A {D.C. programming approach to the sparse generalized eigenvalue problem. Computer Engineering 1~(1), 40.
1612.00939.extracted_bib
[{Wang and Wu(2012)]{wan Wang, Y., Wu, Q., 2012. Sparse pca by iterative elimination algorithm. Advances in Computational Mathematics 36, 137--151.
1612.00939.extracted_bib
[{Zou et~al.(2006)Zou, Hastie, and Tibshirani]{zou Zou, H., Hastie, T., Tibshirani, R., 2006. Sparse principal component analysis. Journal of Computational and Graphical Statistics 15~(2), 265--286.