# Principal Component Analysis


```python
from IPython.display import Image 
Image('../../../python_for_probability_statistics_and_machine_learning.jpg')
```




    

    



The features from a particular dataset that will ultimately prove important for
machine learning can be difficult to know ahead of time.  This is especially
true for problems that do not have a strong physical underpinning. The
row-dimension of the input matrix ($X$) for fitting data in Scikit-learn is the
number of samples and the column dimension is the number of features. There may
be a large number of column dimensions in this matrix, and the purpose of
dimensionality reduction is to somehow reduce these to only those columns that
are important for the machine learning task.

Fortunately, Scikit-learn provides some powerful tools to help uncover the most
relevant features.  Principal Component Analysis (PCA) consists of taking
the input $X$ matrix and (1) subtracting the mean, (2) computing the covariance
matrix, and (3) computing the eigenvalue decomposition of the covariance
matrix. For example, if $X$ has more columns than is practicable for a
particular learning method, then PCA can reduce the number of columns to a more
manageable number.  PCA is widely used in statistics and other areas beyond
machine learning, so it is worth examining what it does in some detail. First,
we need the decomposition module from Scikit-learn.


```python
from sklearn import decomposition
import numpy as np
pca = decomposition.PCA()
```

Let's create some very simple data and apply PCA.


```python
x = np.linspace(-1,1,30)
X = np.c_[x,x+1,x+2] # stack as columns
pca.fit(X)
print(pca.explained_variance_ratio_)
```

    [  1.00000000e+00   4.44789028e-32   6.56191018e-33]


**Programming Tip.**

The `np.c_` is a shorcut method for creating stacked column-wise arrays.




 In this example, the columns are just constant offsets of the first
column. The *explained variance ratio* is the percentage of the variance
attributable to the transformed columns of `X`. You can think of this as the
information that is relatively concentrated in each column of the transformed
matrix `X`. [Figure](#fig:pca_001) shows the graph of this dominant
transformed column in the bottom panel.


```python
%matplotlib inline

from matplotlib.pylab import subplots
fig,axs = subplots(2,1,sharex=True,sharey=True)
ax = axs[0]
_=ax.plot(x,X[:,0],'-k',lw=3)
_=ax.plot(x,X[:,1],'--k',lw=3)
_=ax.plot(x,X[:,2],':k',lw=3)
ax=axs[1]
_=ax.plot(x,pca.fit_transform(X)[:,0],'-k',lw=3)
#ax.tick_params(labelsize='x-large')
```

<!-- dom:FIGURE: [fig-machine_learning/pca_001.png, width=500 frac=0.75]  The
top panel shows the columns of the feature matrix and the bottom panel shows the
dominant component that PCA has extracted.  <div id="fig:pca_001"></div> -->
<!-- begin figure -->
<div id="fig:pca_001"></div>

<p>The top panel shows the columns of the feature matrix and the bottom panel
shows the dominant component that PCA has extracted.</p>


<!-- end figure -->


 To make this more interesting, let's change the slope of each of the
columns as
in the following,


```python
X = np.c_[x,2*x+1,3*x+2,x] # change slopes of columns
pca.fit(X)
print(pca.explained_variance_ratio_)
```

    [  1.00000000e+00   1.36443189e-32   1.71524927e-33   3.68262901e-34]


However, changing the slope did not impact the explained variance
ratio. Again, there is still only one dominant column. This means that PCA is
invariant to both constant offsets and scale changes. This works for functions
as well as simple lines,


```python
x = np.linspace(-1,1,30)
X = np.c_[np.sin(2*np.pi*x),
          2*np.sin(2*np.pi*x)+1,
          3*np.sin(2*np.pi*x)+2] 
pca.fit(X)
print(pca.explained_variance_ratio_)
```

    [  1.00000000e+00   3.63254364e-32   3.47777102e-33]


Once again, there is only one dominant column, which is shown in the
bottom panel of [Figure](#fig:pca_002). The top panel shows the individual
columns of the feature matrix. To sum up, PCA is able to identify and eliminate
features that are merely linear transformations of existing features. This also
works when there is additive noise in the features, although more samples are
needed to separate the uncorrelated noise from between features.

<!-- dom:FIGURE: [fig-machine_learning/pca_002.png, width=500 frac=0.85] The top
panel shows the columns of the feature matrix and the bottom panel shows the
dominant component that PCA has computed. <div id="fig:pca_002"></div>  -->
<!-- begin figure -->
<div id="fig:pca_002"></div>

<p>The top panel shows the columns of the feature matrix and the bottom panel
shows the dominant component that PCA has computed.</p>


<!-- end figure -->


```python
fig,axs = subplots(2,1,sharex=True,sharey=True)
ax = axs[0]
_=ax.plot(x,X[:,0],'-k')
_=ax.plot(x,X[:,1],'--k')
_=ax.plot(x,X[:,2],':k')
ax=axs[1]
_=ax.axis(xmin=-1,xmax=1)
_=ax.plot(x,pca.fit_transform(X)[:,0],'-k')
# ax.tick_params(labelsize='x-large')
```


```python
fig,axs=subplots(1,2,sharey=True)
fig.set_size_inches((8,5))
ax=axs[0]
ax.set_aspect(1/1.6)
_=ax.axis(xmin=-2,xmax=12)
x1 = np.arange(0, 10, .01/1.2)
x2 = x1+ np.random.normal(loc=0, scale=1, size=len(x1))
X = np.c_[(x1, x2)]
good = (x1>5) | (x2>5) 
bad = ~good
_=ax.plot(x1[good],x2[good],'ow',alpha=.3)
_=ax.plot(x1[bad],x2[bad],'ok',alpha=.3)
_=ax.set_title("original data space")
_=ax.set_xlabel("x")
_=ax.set_ylabel("y")

_=pca.fit(X)
Xx=pca.fit_transform(X)
ax=axs[1]
ax.set_aspect(1/1.6)
_=ax.plot(Xx[good,0],Xx[good,1]*0,'ow',alpha=.3)
_=ax.plot(Xx[bad,0],Xx[bad,1]*0,'ok',alpha=.3)
_=ax.set_title("PCA-reduced data space")
_=ax.set_xlabel(r"$\hat{x}$")
_=ax.set_ylabel(r"$\hat{y}$")
```

To see how PCA can simplify machine learning tasks, consider
[Figure](#fig:pca_003) wherein the two classes are separated along the diagonal.
After PCA, the transformed data lie along a single axis where the two classes
can be split using a one-dimensional interval, which greatly simplifies the
classification task. The class identities are preserved under PCA because the
principal component is along the same direction that the classes are separated.
On the other hand, if the classes are separated along the direction
*orthogonal* to the principal component, then the two classes become mixed
under PCA and the classification task becomes much harder. Note that in both
cases, the `explained_variance_ratio_` is the same because the explained
variance ratio does not account for class membership.

<!-- dom:FIGURE: [fig-machine_learning/pca_003.png, width=500 frac=0.85]  The
left panel shows the original two-dimensional data space of two easily
distinguishable classes and the right panel shows the reduced the data space
transformed using PCA. Because the two classes are separated along the principal
component discovered by PCA, the classes are  preserved under the
transformation.  <div id="fig:pca_003"></div>  -->
<!-- begin figure -->
<div id="fig:pca_003"></div>

<p>The left panel shows the original two-dimensional data space of two easily
distinguishable classes and the right panel shows the reduced the data space
transformed using PCA. Because the two classes are separated along the principal
component discovered by PCA, the classes are  preserved under the
transformation.</p>


<!-- end figure -->


```python

fig,axs=subplots(1,2,sharey=True)
ax=axs[0]
x1 = np.arange(0, 10, .01/1.2)
x2 = x1+np.random.normal(loc=0, scale=1, size=len(x1))
X = np.c_[(x1, x2)]
good = x1>x2
bad = ~good 
_=ax.plot(x1[good],x2[good],'ow',alpha=.3)
_=ax.plot(x1[bad],x2[bad],'ok',alpha=.3)
_=ax.set_title("original data space")
_=ax.set_xlabel("x")
_=ax.set_ylabel("y")

_=pca.fit(X)
Xx=pca.fit_transform(X)
ax=axs[1]
_=ax.plot(Xx[good,0],Xx[good,1]*0,'ow',alpha=.3)
_=ax.plot(Xx[bad,0],Xx[bad,1]*0,'ok',alpha=.3)
_=ax.set_title("PCA-reduced data space")
_=ax.set_xlabel(r"\hat{x}")
_=ax.set_ylabel(r"\hat{y}")

```

<!-- dom:FIGURE: [fig-machine_learning/pca_004.png, width=500 frac=0.85]  As
compared with [Figure](#fig:pca_003), the two classes differ along the
coordinate direction that is orthogonal to the principal component. As a result,
the two classes are no longer distinguishable after transformation. <div
id="fig:pca_004"></div>  -->
<!-- begin figure -->
<div id="fig:pca_004"></div>

<p>As compared with [Figure](#fig:pca_003), the two classes differ along the
coordinate direction that is orthogonal to the principal component. As a result,
the two classes are no longer distinguishable after transformation.</p>


<!-- end figure -->


PCA works by decomposing the covariance matrix of the data using the Singular
Value Decomposition (SVD). This decomposition exists for all matrices and
returns the following factorization for an arbitrary matrix $\mathbf{A}$,

$$
\mathbf{A} =  \mathbf{U} \mathbf{S}  \mathbf{V}^T
$$

 Because of the symmetry of the covariance matrix, $\mathbf{U} =
\mathbf{V}$. The elements of the diagonal matrix $\mathbf{S}$ are the singular
values of $\mathbf{A}$ whose squares are the eigenvalues of $\mathbf{A}^T
\mathbf{A}$. The eigenvector matrix $\mathbf{U}$ is orthogonal: $\mathbf{U}^T
\mathbf{U} =\mathbf{I}$. The singular values are in decreasing order so that
the first column of $\mathbf{U}$ is the axis corresponding to the largest
singular value. This is the first dominant column that PCA identifies. The
entries of the covariance matrix are of the form $\mathbb{E}(x_i x_j)$ where
$x_i$ and $x_j$ are different features [^covariance]. This means that the
covariance matrix is filled with entries that attempt to uncover mutually
correlated relationships between all pairs of columns of the feature matrix.
Once these have been tabulated in the covariance matrix, the SVD finds optimal
orthogonal transformations to align the components along the directions most
strongly associated with these correlated relationships. Simultaneously,
because orthogonal matrices have columns of unit-length, the SVD collects the
absolute squared lengths of these components into the $\mathbf{S}$ matrix.  In
our example above in [Figure](#fig:pca_003), the two feature vectors were
obviously correlated along the
diagonal, meaning that PCA selected that diagonal direction as the principal
component.

[^covariance]: Note that these entries are constructed from the data
using an estimator of the covariance matrix because we do not have
the full probability densities at hand.

We have seen that PCA is a powerful dimension reduction method that is
invariant to linear transformations of the original feature space. However,
this method performs poorly with transformations that are nonlinear. In that
case, there are a wide range of extensions to PCA, such as Kernel PCA, that are
available in Scikit-learn, which allow for embedding parameterized
non-linearities into the PCA at the risk of overfitting.

## Independent Component Analysis

Independent Component Analysis (ICA) via the `FastICA` algorithm is also
available in Scikit-learn. This method is fundamentally different from PCA
in that it is the small differences between components that are emphasized,
not the large principal components.  This method is adopted from signal
processing.  Consider a matrix of signals ($\mathbf{X}$) where the rows are
the samples and the columns are the different signals. For example, these
could be EKG signals from multiple leads on a single patient. The analysis
starts with the following model,

<!-- Equation labels as ordinary links -->
<div id="eq:ICA"></div>

$$
\begin{equation}
\mathbf{X} =  \mathbf{S}\mathbf{A}^T
\label{eq:ICA} \tag{1}
\end{equation}
$$

 In other words, the observed signal matrix is an unknown mixture
($\mathbf{A}$) of some set of conformable, independent random sources
$\mathbf{S}$,

$$
\mathbf{S}=\left[ \mathbf{s}_1(t),\mathbf{s}_2(t),\ldots,\mathbf{s}_n(t)\right]
$$

 The distribution on the random sources is otherwise unknown, except
there can be at most one Gaussian source, otherwise, the mixing matrix
$\mathbf{A}$ cannot be identified because of technical reasons.  The problem in
ICA is to find $\mathbf{A}$ in Equation [eq:ICA](#eq:ICA) and thereby un-mix the
$s_i(t)$ signals, but this cannot be solved without a strategy to reduce the
inherent arbitrariness in this formulation.

To make this concrete, let us simulate the situation with the following code,


```python
np.random.seed(123456)
```


```python
from numpy import matrix, c_, sin, cos, pi
t = np.linspace(0,1,250)
s1 = sin(2*pi*t*6)
s2 =np.maximum(cos(2*pi*t*3),0.3)
s2 = s2 - s2.mean()
s3 = np.random.randn(len(t))*.1

# normalize columns
s1=s1/np.linalg.norm(s1)
s2=s2/np.linalg.norm(s2)
s3=s3/np.linalg.norm(s3)
S =c_[s1,s2,s3] # stack as columns

# mixing matrix
A = matrix([[  1,  1,1],
            [0.5, -1,3],
            [0.1, -2,8]])
X= S*A.T # do mixing
```


```python
fig,axs=subplots(3,2,sharex=True)
fig.set_size_inches((8,8))
X = np.array(X)

_=axs[0,1].plot(t,-X[:,0],'k-')
_=axs[1,1].plot(t,-X[:,1],'k-')
_=axs[2,1].plot(t,-X[:,2],'k-')
_=axs[0,0].plot(t,s1,'k-')
_=axs[1,0].plot(t,s2,'k-')
_=axs[2,0].plot(t,s3,'k-')

_=axs[2,0].set_xlabel('$t$',fontsize=18)
_=axs[2,1].set_xlabel('$t$',fontsize=18)
_=axs[0,0].set_ylabel('$s_1(t)$        ',fontsize=18,rotation='horizontal')
_=axs[1,0].set_ylabel('$s_2(t)$        ',fontsize=18,rotation='horizontal')
_=axs[2,0].set_ylabel('$s_3(t)$        ',fontsize=18,rotation='horizontal')
for ax in axs.flatten():
    _=ax.yaxis.set_ticklabels('')

_=axs[0,1].set_ylabel('        $X_1(t)$',fontsize=18,rotation='horizontal')
_=axs[1,1].set_ylabel('        $X_2(t)$',fontsize=18,rotation='horizontal')
_=axs[2,1].set_ylabel('        $X_3(t)$',fontsize=18,rotation='horizontal')
_=axs[0,1].yaxis.set_label_position("right")
_=axs[1,1].yaxis.set_label_position("right")
_=axs[2,1].yaxis.set_label_position("right")

```

<!-- dom:FIGURE: [fig-machine_learning/pca_008.png, width=500 frac=0.85] The
left column shows the original signals and the right column shows the mixed
signals. The object of ICA is to recover the left column from the right. <div
id="fig:pca_008"></div>  -->
<!-- begin figure -->
<div id="fig:pca_008"></div>

<p>The left column shows the original signals and the right column shows the
mixed signals. The object of ICA is to recover the left column from the
right.</p>


<!-- end figure -->


 The individual signals ($s_i(t)$) and their mixtures ($X_i(t)$) are
shown in [Figure](#fig:pca_008). To recover the individual signals using ICA,
we use the `FastICA` object and fit the parameters on the `X` matrix,


```python
from sklearn.decomposition import FastICA
ica = FastICA()
# estimate unknown S matrix
S_=ica.fit_transform(X)
```

The results of this estimation are shown in [Figure](#fig:pca_009),
showing that ICA is able to recover the original signals from the observed
mixture. Note that ICA is unable to distinguish the signs of the recovered
signals or preserve the order of the input signals.


```python
fig,axs=subplots(3,2,sharex=True)
fig.set_size_inches((8,8))
X = np.array(X)

_=axs[0,1].plot(t,-S_[:,2],'k-')
_=axs[1,1].plot(t,-S_[:,1],'k-')
_=axs[2,1].plot(t,-S_[:,0],'k-')
_=axs[0,0].plot(t,s1,'k-')
_=axs[1,0].plot(t,s2,'k-')
_=axs[2,0].plot(t,s3,'k-')

_=axs[2,0].set_xlabel('$t$',fontsize=18)
_=axs[2,1].set_xlabel('$t$',fontsize=18)
_=axs[0,0].set_ylabel('$s_1(t)$        ',fontsize=18,rotation='horizontal')
_=axs[1,0].set_ylabel('$s_2(t)$        ',fontsize=18,rotation='horizontal')
_=axs[2,0].set_ylabel('$s_3(t)$        ',fontsize=18,rotation='horizontal')
for ax in axs.flatten():
    _=ax.yaxis.set_ticklabels('')

_=axs[0,1].set_ylabel('        $s_1^\prime(t)$',fontsize=18,rotation='horizontal')
_=axs[1,1].set_ylabel('        $s_2^\prime(t)$',fontsize=18,rotation='horizontal')
_=axs[2,1].set_ylabel('        $s_3^\prime(t)$',fontsize=18,rotation='horizontal')
_=axs[0,1].yaxis.set_label_position("right")
_=axs[1,1].yaxis.set_label_position("right")
_=axs[2,1].yaxis.set_label_position("right")

```

<!-- dom:FIGURE: [fig-machine_learning/pca_009.png, width=500 frac=0.85] The
left column shows the original signals and the right column shows the signals
that ICA was able to recover. They match exactly, outside of a possible sign
change. <div id="fig:pca_009"></div>  -->
<!-- begin figure -->
<div id="fig:pca_009"></div>

<p>The left column shows the original signals and the right column shows the
signals that ICA was able to recover. They match exactly, outside of a possible
sign change.</p>


<!-- end figure -->


To develop some intuition as to how ICA accomplishes this feat, consider the
following two-dimensional situation with two uniformly distributed independent
variables, $u_x,u_y \sim \mathcal{U}[0,1]$. Suppose we apply the
following orthogonal rotation matrix to these variables,

$$
\begin{bmatrix}
u_x^\prime \\\
u_y^\prime
\end{bmatrix}=
\begin{bmatrix}
\cos(\phi) & -\sin(\phi) \\\
\sin(\phi) & \cos(\phi)
\end{bmatrix}
\begin{bmatrix}
u_x \\\
u_y
\end{bmatrix}
$$

<!-- dom:FIGURE: [fig-machine_learning/pca_005.png, width=500 frac=0.85] The
left panel shows two classes labeled on the $u_x,u_y$ uniformly independent
random variables. The right panel shows these random variables after a rotation,
which removes their mutual independence and makes it hard to separate the two
classes along the coordinate directions. <div id="fig:pca_005"></div>  -->
<!-- begin figure -->
<div id="fig:pca_005"></div>

<p>The left panel shows two classes labeled on the $u_x,u_y$ uniformly
independent random variables. The right panel shows these random variables after
a rotation, which removes their mutual independence and makes it hard to
separate the two classes along the coordinate directions.</p>


<!-- end figure -->


 The so-rotated variables $u_x^\prime,u_y^\prime$ are no longer
independent, as shown in [Figure](#fig:pca_005). Thus, one way to think about
ICA is as a search through orthogonal matrices so that the independence is
restored. This is where the prohibition against Gaussian distributions arises.
The two dimensional Gaussian distribution of independent variables is
proportional the following,

$$
f(\mathbf{x})\propto\exp(-\frac{1}{2}\mathbf{x}^T \mathbf{x})
$$

 Now, if we similarly rotated the $\mathbf{x}$ vector as,

$$
\mathbf{y} = \mathbf{Q} \mathbf{x}
$$

 the resulting density for $\mathbf{y}$ is obtained by plugging in
the following,

$$
\mathbf{x} = \mathbf{Q}^T \mathbf{y}
$$

 because the inverse of an orthogonal matrix is its transpose, we
obtain

$$
f(\mathbf{y})\propto\exp(-\frac{1}{2}\mathbf{y}^T \mathbf{Q}\mathbf{Q}^T
\mathbf{y})=\exp(-\frac{1}{2}\mathbf{y}^T \mathbf{y})
$$

 In other words, the transformation is lost on the $\mathbf{y}$
 variable. This means that ICA cannot search over orthogonal transformations if
 it is blind to them, which explains the restriction of Gaussian random
 variables. Thus, ICA is a method that seeks to maximize the non-Gaussian-ness
 of the transformed random variables.  There are many methods to doing this,
 some of which involve cumulants and others that use the
*negentropy*,

$$
\mathcal{J}(Y) = \mathcal{H}(Z)-\mathcal{H}(Y)
$$

 where $\mathcal{H}(Z)$ is the information entropy of the
Gaussian random variable $Z$ that has the same variance as $Y$. Further
details would take us beyond our scope, but that is the outline of how
the FastICA algorithm works.

The implementation of this method in Scikit-learn includes two different ways
of extracting more than one independent source component. The *deflation*
method iteratively extracts one component at a time using a incremental
normalization step. The *parallel* method also uses the single-component method
but carries out normalization of all the components simultaneously, instead of
for just the newly computed component.  Because ICA extracts
independent components, a whitening step is used beforehand to balance the
correlated components from the data matrix. Whereas PCA returns
uncorrelated components along dimensions optimal for Gaussian random
variables, ICA returns components that are as far from the Gaussian density
as possible.

The left panel on [Figure](#fig:pca_005) shows the orignal uniform random
sources. The white and black colors distinguish between two classes. The right
panel shows the mixture of these sources, which is what we observe as input
features. The top row of [Figure](#fig:pca_006) shows the PCA (left) and ICA
(right) transformed data spaces.  Notice that ICA is able to un-mix the two
random sources whereas PCA transforms along the dominant diagonal. Because ICA
is able to preserve the class membership, the data space can be reduced to two
non-overlapping sections, as shown. However, PCA cannot achieve a similiar
separation because the classes are mixed along the dominant diagonal that PCA
favors as the main component in the decomposition.

<!-- dom:FIGURE: [fig-machine_learning/pca_006.png, width=500 frac=0.85]  The
panel on the top left shows two classes in a plane after a  rotation. The bottom
left panel shows the result of dimensionality reduction using PCA, which causes
mixing between the two classes. The top right panel shows the ICA transformed
output and the lower right panel shows that, because ICA was able to un-rotate
the data, the lower dimensional data maintains the separation between the
classes. <div id="fig:pca_006"></div>  -->
<!-- begin figure -->
<div id="fig:pca_006"></div>

<p>The panel on the top left shows two classes in a plane after a  rotation. The
bottom left panel shows the result of dimensionality reduction using PCA, which
causes mixing between the two classes. The top right panel shows the ICA
transformed output and the lower right panel shows that, because ICA was able to
un-rotate the data, the lower dimensional data maintains the separation between
the classes.</p>


<!-- end figure -->


For a good principal component analysis treatment, see
[[richert2013building]](#richert2013building),
[[alpaydin2014introduction]](#alpaydin2014introduction),
[[cuesta2013practical]](#cuesta2013practical), and
[[izenman2008modern]](#izenman2008modern). Independent Component Analysis is
discussed in more detail
in [[hyvarinen2004independent]](#hyvarinen2004independent).

<!-- #  *Learning from Data*, p. 263 -->
<!-- #  *Building machine learning systems*, p.231 -->
<!-- #  *Introduction to machine learning by Alpaydin*, p148 -->
<!-- #  *Practical Data Analysis Cuesta*,p.143 -->
<!-- #  *Modern Multivariate statistical techniques Izenman*, p.553, 558 -->
<!-- #  *Independent Component Analysis Hyvarinen.pdf*, p.147 -->

9778d96513cd8bb69beeda7f4feaded7ac1e5d73
ipynb
Jupyter Notebook
chapters/machine_learning/notebooks/pca.ipynb
nsydn/Python-for-Probability-Statistics-and-Machine-Learning
d3e0f8ea475525a694a975dbfd2bf80bc2967cc6
MIT
2016-05-05T19:08:27.000Z
2022-03-31T05:09:19.000Z
chapters/machine_learning/notebooks/pca.ipynb
crlsmcl/https-github.com-unpingco-Python-for-Probability-Statistics-and-Machine-Learning
6fd69459a28c0b76b37fad79b7e8e430d09a86a5
MIT
2016-05-12T22:18:58.000Z
2019-11-06T14:37:06.000Z
chapters/machine_learning/notebooks/pca.ipynb
crlsmcl/https-github.com-unpingco-Python-for-Probability-Statistics-and-Machine-Learning
6fd69459a28c0b76b37fad79b7e8e430d09a86a5
MIT
2016-05-27T01:42:05.000Z
2022-03-27T11:20:27.000Z
Qwen/Qwen-72B
1. YES
2. YES