Principal components analysis

Suppose we conduct an experiment $n$ times and collect $p$ pieces of information from each experiment. Then we can suppose that we collected a dataset $\mathcal{D} = \{ (x_{11},x_{12},\cdots, x_{1p} ), (x_{21},x_{22},\cdots, x_{2p} ), \cdots, (x_{n1},x_{n2},\cdots, x_{np} ) \}$ and can organize this dataset into the following dataframe

\begin{align} \mathcal{F} = \begin{bmatrix} x_{11} & x_{12} & \cdots & x_{1p}\\ x_{21} & x_{22} & \cdots & x_{2p}\\ \vdots & \vdots & \vdots & \ddots \\ x_{n1} & x_{n2} & \cdots & x_{np}\\ \end{bmatrix} \end{align}

(1)

It is important to acknowledge that a single observation is the tuple $(x_{i1} \cdots x_{1p})$ and so we assume that a tuple of random variables jointly generated this single observation. In other words, it may not be appropriate to treat each piece of information (often called a covariate) as being generated independently of other covariates: they are generated together.

The goal of principal components analysis is to find combinations of covariates that “best” describe our original dataset $\mathcal{D}$ .

In this lesson we will aim to understand:

How principal components defines “best”
The interplay between the (linear) algebra and geometry of principal components
How to interpret the outputs of out principal components analysis
How principal components analysis differs from linear regression

Simulated data¶

We will first build a simulated dataset to illustrate the main concepts in principal components analysis. Let the random variables $[X_{1}, X_{2}]$ have a multivariate normal distribution (MVN)

\begin{align} &\begin{bmatrix} X_{1} \\ X_{2} \end{bmatrix} \sim \text{MVN}( \mu, \Sigma ) \\ & \mu = \begin{bmatrix} -1 \\ 1 \end{bmatrix} ; \; \Sigma = \begin{bmatrix} 1 & 0.3 \\ 0.3 & 2 \end{bmatrix} \end{align}

(2)

As a recap, the vector $\mu$ decribes the expected value for $X_{1}$ (the value -1) and the expected value for $X_{2}$ (the value 1). The covariance matrix $\Sigma$ describes the variance for $X_{1}$ (the value 1), the variance for $X_{2}$ (the value 2), and the covariance between $X_{1}$ and $X_{2}$ (the value 0.3).

We can use numpy to help use generate $n=100$ observations from the above MVN where every observation contains two pieces of information (i.e. two covariates). In addition to storing our dataset, we will visualize this information by scatterplotting values of the first covariate versus values of the second covariate.

Plotting and centering our data¶

Below, you will see two plots side-by-side. The first plot includes the scatter (blue circles) as well as a vertical and horizontal line at the values 0 (dashed blue lines). Though we typically draw the horizontal and vertical axes for a plot at 0, the value zero is arbitrary. A more reasonable place to draw our axes may be the expected values for the first (-1) and second (1) covariate. We drew those axes as solid blue lines.

The solid blue lines (expected values) look to be better centered over our dataset than the dashed lines (zeros).

The above is the motivation behind centering. Instead of plotting our raw data $\mathcal{D} = \{ (x_{11}, x_{12}),(x_{21}, x_{22}),\cdots, (x_{n1}, x_{n2}) \}$ we can plot a centered dataset that describes the distance of each observation from their expected value or $(x_{i1} - \mu_{1}, x_{i2}-\mu_{2})$ for all $i$ observations. Another way to look at this is with our dataframe

\begin{align} \underset{-}{\mathcal{F}} = \begin{bmatrix} x_{11} - \mu_{1} & x_{12} - \mu_{2} \\ x_{21} - \mu_{1} & x_{22} - \mu_{2} \\ \vdots & \vdots \\ x_{n1}- \mu_{1} & x_{n2} - \mu_{2} \\ \end{bmatrix} = \mathcal{F}_{n \times 2} - \mu_{1 \times 2} \end{align}

(3)

Importantly, the original dataset describes the distance of each covariate from the (arbitrary) value zero. In other words,

\begin{align} \mathcal{F} = \begin{bmatrix} x_{11} & x_{12} \\ x_{21} & x_{22} \\ \vdots & \vdots \\ x_{n1} & x_{n2} \\ \end{bmatrix} = \begin{bmatrix} x_{11} - 0 & x_{12} - 0 \\ x_{21} - 0 & x_{22} - 0 \\ \vdots & \vdots \\ x_{n1}- 0 & x_{n2} - 0 \\ \end{bmatrix} = \mathcal{F}_{n \times 2} - 0_{1 \times 2} \end{align}

(4)

The rightmost panel in our figure below plots the centered dataset (solid blue circles) as well as vertical and horizontal axes at zero. The reason that these horizonal/vertical lines at zero are centered over the data now is because the value zero for our centered dataset represents the expected value. For example, the $i^{\text{th}}$ observation is $(x_{1}= x_{i1}-\mu_{1}, x_{2} = x_{i1}-\mu_{2} )$ . When $x_{1}=0$ then $x_{i1}-\mu_{1}=0$ or $x_{i1}=\mu_{1}$ . The same argument could be used for the second covariate as well. The centered datset describes the distance from each covariate’s expected value.

import numpy as np 
import pandas as pd 
import seaborn as sns
import matplotlib.pyplot as plt 

#--Specify the MVN expected value M and the covariance COV
#--Generate 100 observations from this vector of random variables  
M   = np.array([-1,1]) 
COV = np.array([[ 1, 0.6 ], [0.6 , 2]])

#--Note: these three lines are not important to know. 
#--You can just assume that this is the process by which we generate MVN observations. 
L   = np.linalg.cholesky(COV)
z   = np.random.normal( 0,1 ,size = (100,2) )
d   = M + z@L


#--Visualize this dataset 
fig,axs = plt.subplots(1,2)

#--non-centered dataset
ax=axs[0]
ax.scatter(d[:,0], d[:,1])

ax.axhline(0, alpha=0.5, ls="--")
ax.axvline(0, alpha=0.5, ls="--")

ax.axhline(1,lw=2)
ax.axvline(-1,lw=2)

ax.set_xlabel("Covariate 1 values")
ax.set_ylabel("Covariate 2 values")

#--centered dataset
d_centered = d - M #<--subtract the expected value from each covariate

ax=axs[1]

ax.scatter(d_centered[:,0], d_centered[:,1])

ax.axhline(0)
ax.axvline(0)

ax.set_xlabel("Covariate 1 values centered")
ax.set_ylabel("Covariate 2 values centered")

fig.set_size_inches(9,4)

plt.show()

Geometric intuition behind Principal Components Analysis¶

Centering our dataset (right panel above) helped us more naturally descirbe our dataset. The zero horizonal/vertical values now correspond to the expected value for the first/second covariate. However, centering only allows us to shift our axes up/down and left/right to better decribe our data. What about rotation?

Geometrically, Principal Components Analysis asks “Can we rotate our axes to better characterize our dataset?”. Well, if we are willing to consider rotating our axes then we need to decide on what we mean when say better.

Lets start with just the horizontal axis and no vertical axis.

Below is a short video that shows several iterations of the following steps:

(Top left) Produce a scatter plot of $x_{1}$ on the horizonal axis and $x_{2}$ on the vertical axis.
(Top right) Include only the horizontal axis and project all of the points onto this axis. To be clear, this projection maps every data point $(x_{1},x_{2}) \to x_{1}$ .
(Bottom left) Plot a histogram of the projected points. In other words, produce a histogram of the $x_{1}$ values and include a red shaded region that charazterizing 95% probability. From this sample of $x_{1}$ points we will compute the variance or $\text{Var}(x_{1},x_{2},\cdots,x_{n})$ . This value is titled at the top of the plot.
(Bottom right) Append this variance value to a list of all the variances, one for each iteration, and plot a histogram of these variances.

Its important here to note what objects are random and which are not. We can see in this video os simultions that randomness does not come from the projection (our choise of the horizontal axis). Over each iteration the horizontal axis never changes. Instead, the randomness comes from our sample of data points.

import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

from matplotlib.animation import FuncAnimation
from IPython.display import HTML

# --------------------------------------------------
# Parameters
# --------------------------------------------------

M   = np.array([-1,1])
COV = np.array([[1,0.6],[0.6,2]])

variances = []

fig,axs = plt.subplots(2,2, figsize=(6,6))
fig.tight_layout()

# --------------------------------------------------
# Animation frame
# --------------------------------------------------

def animate(frame):

    global variances

    # clear axes
    for ax in axs.ravel():
        ax.clear()

    # -------------------------
    # sample new dataset
    # -------------------------

    L = np.linalg.cholesky(COV)

    z = np.random.normal(
        0,1,
        size=(100,2)
    )

    d = M + z @ L.T

    d_centered = d - M

    var = np.var(d_centered[:,0])
    variances.append(var)

    # -------------------------
    # panel 1
    # -------------------------

    ax = axs[0,0]

    ax.scatter(
        d_centered[:,0],
        d_centered[:,1]
    )

    ax.axhline(0)
    ax.axvline(0,ls="--",alpha=0.4)

    ax.set_xlim(-4,4)
    ax.set_ylim(-4,4)

    ax.set_title(
        f"Dataset #{frame+1}"
    )

    # -------------------------
    # panel 2
    # -------------------------

    ax = axs[0,1]

    ax.scatter(
        d_centered[:,0],
        d_centered[:,1]
    )

    ax.axhline(0,color="black")

    ax.vlines(
        x=d_centered[:,0],
        ymin=0,
        ymax=d_centered[:,1],
        alpha=0.2
    )

    ax.set_xlim(-4,4)
    ax.set_ylim(-4,4)

    # -------------------------
    # panel 3
    # -------------------------

    ax = axs[1,0]

    x = d_centered[:,0]
    var = np.var(x, ddof=1)
    sd  = np.sqrt(var)

    ci_low  = -1.96 * sd
    ci_high =  1.96 * sd

    sns.histplot(
        x=x,
        bins=15,
        ax=ax
    )

    ax.axvspan(
        ci_low,
        ci_high,
        color="red",
        alpha=0.20
    )

    ax.axvline(ci_low, color="red", linestyle="--")
    ax.axvline(ci_high, color="red", linestyle="--")

    ax.set_xlim(-4,4)
    ax.set_title(
        f"Sample variance={var:.2f}"
    )

    # -------------------------
    # panel 4
    # -------------------------

    ax = axs[1,1]

    sns.histplot(
        x=variances,
        bins=20,
        ax=ax
    )

    # red dot showing the newest variance estimate
    ax.scatter(
        var,
        0,
        color="red",
        s=50,
        zorder=5
    )
    ax.scatter(
        np.mean(variances),
        0.25,
        color="red",
        s=50,
        marker="x",
        zorder=5
    )

    ax.set_title(
        f"{len(variances)} variance estimates"
    )
# --------------------------------------------------
# Create movie
# --------------------------------------------------

ani = FuncAnimation(
    fig,
    animate,
    frames=100,
    interval=500
)

plt.close(fig)

HTML(ani.to_jshtml())

Loading...

Now that we see how to map a projection of $x$ onto the horizontal axis, it is intuitive that we might want to consider projections onto more complicated directions than just the horizontal (or vertical for that matter) projection.

A projection of one vector $x$ onto another direction $u$ ¶

Suppose we have a vector $x$ (in red) and want to project this vector onto another called $u$ (in blue). See the picture below. Then we know that this projection, called $x_{u}$ must be in the same direction as $u$ . We can write this down mathemtically as

\begin{align} x_{u} = \alpha u \end{align}

(5)

where $\alpha$ defines how far in the direction of $u$ we should move. But, we are not sure how far in the direction of $u$ we should extend this vector to be considered a “projection” of $x$ onto $u$ .

One reasonable property we might suggest for our projection $x_{u}$ , is that this vector will be in the direction of $u$ and as close as possible to our original $x$ . The vector that is closest to $x$ but in the direction of $u$ is the on such that $x-u$ and $u$ are orthogonal to one another (instead or orthogonal, perpendicular is ok to use here as well).

To see why, intuitively, we plotted the difference between two projections onto u. The first difference we plotted as a black dashed line that is perpendicular to $u$ . The second difference, drawn as a purple dashed-dot line is not perpendicular to $u$ . You can see that the line that is perpendicular to $u$ is shorter in length than the other.

We can also see why this is the case algebraically by looking at the distance between $x$ and the projection of $x$ onto $u$ or $x_{u}$ . The typical way to define the length of a vector is

\begin{align} ||x|| & = (x^{T}x)^{1/2}. \\ \end{align}

(6)

It is important to note that if the length of $x$ is smallest than the squared length is smallest too. The squared length is almost always better to work with because

\begin{align} ||x||^{2} & = (x^{T}x) \end{align}

(7)

In other words, we want to find the value $\alpha$ that minimizes

\begin{align} ||x - x_{u}||^{2} = ||x - \alpha u||^{2} & = (x-\alpha u)^{T} (x-\alpha u) \\ &= x^{T}x - 2\alpha x^{T} u + \alpha^{2} u^{T}u \end{align}

(8)

Lets take the derivative and set this to zero to find that

\begin{align} f(\alpha) &= x^{T}x - 2\alpha x^{t} u + \alpha^{2} u^{T}u\\ \frac{df}{d \alpha} = - 2 x^{T} u + 2 \alpha u^{T}u = 0 \alpha = \frac{x^{T} u }{ u^{T} u } \end{align}

(9)

Now when we look at the vector $u$ and the residual $(x - \alpha u)$ we see that

\begin{align} u^{T} \left ( x - \frac{x^{T} u }{ u^{T} u } u \right) &= u^{T} x - \frac{x^{T} u}{u^{T} u} u^{T} u \\ &= u^{T} x - x^{T} u = 0 \end{align}

(10)

There is an important lesson here. When we want to choose a vector in the direction of $u$ such that it is closest to $x$ then the vector that we need to choose is one such that $u$ and the residual, $(x-\alpha u)$ , are orthogonal.

import matplotlib.pyplot as plt 
import numpy as np 

fig, ax = plt.subplots()

ax.axhline(0,alpha=0.10,color="k")
ax.axvline(0,alpha=0.10,color="k")
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['bottom'].set_visible(False)
plt.grid()

x = np.array([-1,2])
u = np.array([2,3])

ax.annotate('', 
            xy=( x[0], x[1]),             # Point the arrow points to
            xytext=(0,0),         # Where the arrow starts (tail)
            arrowprops=dict(arrowstyle="->", color="red", lw=1.,ls="-"))
ax.text(x[0],x[1],s= r"$x$",fontsize=12,ha="right",va="top",color="red")


ax.annotate('', 
            xy=( u[0], u[1]),             # Point the arrow points to
            xytext=(0,0),         # Where the arrow starts (tail)
            arrowprops=dict(arrowstyle="->", color="blue", lw=1.,ls="-"))
ax.text(u[0],u[1],s= r"$u$",fontsize=12,ha="left",va="top",color="blue")


ax.text(-0.5,.75, rotation=-25,s= r"$x_{u}-x$",fontsize=12,ha="left",va="bottom",color="black")


alpha = (x.T @ u) / (u.T@u)
proj = alpha*u

ax.annotate('', 
            xy=( proj[0], proj[1]),             # Point the arrow points to
            xytext=(x[0], x[1]),         # Where the arrow starts (tail)
            arrowprops=dict(arrowstyle="->", color="black", lw=1.,ls="--"))

eps=0.05
ax.annotate('', 
            xy=( proj[0]+eps, proj[1]-eps),             # Point the arrow points to
            xytext=(0+eps, 0-eps),         # Where the arrow starts (tail)
            arrowprops=dict(arrowstyle="->", color="purple", lw=1.5))
ax.text(proj[0]+eps, proj[1]-eps,s= r"$x_{u}$",fontsize=12,ha="left",va="top")

#--alternative projection 
alpha = 1.5*alpha
alt_u = alpha*u

ax.annotate('', 
            xy=( alt_u[0], alt_u[1]),             # Point the arrow points to
            xytext=( x[0], x[1]),         # Where the arrow starts (tail)
            arrowprops=dict(arrowstyle="->", color="black", lw=1.5,ls="-."))


ax.set_xlim(-2,3)
ax.set_ylim(-2,3)

plt.show()

Recap If i want to project the vector $x$ onto (or in the direction of) the vector $u$ then that projection is

\begin{align} x_{u} = \frac{ x^{T}u }{u^{T}u} u = \frac{ u^{T}x }{u^{T}u} u \end{align}

(11)

We can also extend this idea to project many vectors onto the vector $u$ . For example, suppose that we have three vectors: $x_{1}$ , $x_{2}$ , and $x_{3}$ and wish to project them all onto $u$ . Then we can form the following matrix multiplication to do this in one swoop.

its important to note below that we need to divide the vector $u$ by $u^{T}u$

Let $\tilde{u} = \frac{u}{u^{T}u}$ then

\begin{align} X \tilde{u} = \frac{Xu}{u^{T}u} = \begin{bmatrix} x_{1}^{T} u \Big / u^{T}u \\ x_{2}^{T} u \Big / u^{T}u \\ \vdots \\ x_{n}^{T} u \Big / u^{T}u \\ \end{bmatrix} = \begin{bmatrix} \alpha_{1} \\ \alpha_{2} \\ \vdots\\ \alpha_{n} \\ \end{bmatrix} \end{align}

(12)

Projecting onto different directions¶

Back to our data as a stacked matrix of vectors or

\begin{align} X = \begin{bmatrix} \rule[0.5ex]{1.2em}{0.4pt} x_{1} \rule[0.5ex]{1.2em}{0.4pt} \\ \rule[0.5ex]{1.2em}{0.4pt} x_{2} \rule[0.5ex]{1.2em}{0.4pt}\\ \vdots \\ \rule[0.5ex]{1.2em}{0.4pt} x_{n} \rule[0.5ex]{1.2em}{0.4pt}\\ \end{bmatrix} \end{align}

(13)

Given a direction vector, $u$ , we can project our data onto this $u$ by first computing $\tilde{u} = u / (u^{T}u)$ and and then computing $X_{u} = X u$ . The notation $X_{u}$ stands for the vectors $x_{1}, x_{2}, \cdots, x_{n}$ projected in the direction of $u$ .

Lets see how this works for the horizontal and vertical axes. The horizontal axis is in the direction $u = [1,0]$ , and we see that $u^{T}u = [1,0]^{T}[1,0] = 1^{2} + 0^{2} = 1$ and so we can project our data onto the horizontal axis by computing

\begin{align} X_{ [1,0] } &= X \begin{bmatrix} 1 \\ 0 \end{bmatrix} \\ &= \begin{bmatrix} x_{11} & x_{12} \\ x_{21} & x_{12} \\ \vdots & \vdots \\ x_{n1} & x_{n2} \\ \end{bmatrix}\begin{bmatrix} 1 \\ 0 \end{bmatrix} = \begin{bmatrix} x_{11} \\ x_{21} \\ \vdots\\ x_{n1} \\ \end{bmatrix} \end{align}

(14)

In other words “projecting” onto the horizontal axis is just selecting the first entry in each vector. In this 2-dimensional example, we can also project onto the vertical axis by choosing the direction vector $u = [0,1]^{T}$ , and so then

\begin{align} X_{ [0,1] } = \begin{bmatrix} x_{12} \\ x_{22} \\ \vdots\\ x_{n2} \\ \end{bmatrix} \end{align}

(15)

Below we will project our data, $X$ , onto three different $u$ vectors: (Top right) $u=[1,0]^{T}$ , (Bottom Left) $u = [1,-1]^{T}$ , and (Bottom Right) $u=[1,1]^{T}$ . The panel on the Top Left shows the typical visualization of our data not projected onto a $u$ vector.


import numpy as np 
import matplotlib.pyplot as plt 
fig,axs = plt.subplots(2,2)

ax=axs[0,0]
ax.scatter(d_centered[:,0], d_centered[:,1])

ax.axhline(0)
ax.axvline(0,ls="--",alpha=0.4)

ax.set_xlabel("Covariate 1 values centered")
ax.set_ylabel("Covariate 2 values centered")
ax.set_xlim(-4,4)
ax.set_ylim(-4,4)


ax=axs[0,1]

ax.scatter(d_centered[:,0], d_centered[:,1])
ax.axhline(0,color="black")
ax.vlines(x=d_centered[:,0], ymin=0, ymax=d_centered[:,1],alpha=0.2)

u = np.array([1,0]).reshape(2,1)
alphas = d_centered @ u
var    = np.var(alphas)

dist   = d_centered - alphas@u.T
dist   = np.linalg.norm(dist,axis=1)
avg_dist = np.mean(dist)

ax.text(0.01,0.99,s="V = {:.1f}".format(var),ha="left",va="top",transform=ax.transAxes)
ax.text(0.99,0.01,s="D = {:.1f}".format(avg_dist),ha="right",va="bottom",transform=ax.transAxes)



ax.set_xlabel("Covariate 1 values centered")
ax.set_ylabel("Covariate 2 values centered")
ax.set_xlim(-4,4)
ax.set_ylim(-4,4)


ax=axs[1,0]
u      = np.array([-1,1]).reshape(2,1)
u      = u/np.linalg.norm(u)

alphas = (d_centered @ u) #<-- N X 1

var    = np.var(alphas)

dist   = d_centered - alphas@u.T
dist   = np.linalg.norm(dist,axis=1)
avg_dist = np.mean(dist)


for point,alpha in zip(d_centered, alphas):
    alpha = float(alpha[0])
    u0    = float(u[0,0])
    u1    = float(u[1,0])
    ax.plot( [alpha*u0, point[0] ] ,  [alpha*u1, point[1]  ], color="blue", alpha=0.1  )
    ax.scatter( point[0], point[1], color ="purple" )
    
ax.plot( [-4,4], [4,-4], color="black" )  
ax.text(0.01,0.99,s="V = {:.1f}".format(var),ha="left",va="top",transform=ax.transAxes)
ax.text(0.99,0.01,s="D = {:.1f}".format(avg_dist),ha="right",va="bottom",transform=ax.transAxes)

ax.set_xlabel("Covariate 1 values centered")
ax.set_ylabel("Covariate 2 values centered")
ax.set_xlim(-4,4)
ax.set_ylim(-4,4)

ax=axs[1,1]
u      = np.array([1,1]).reshape(2,1)
u      = u/np.linalg.norm(u)

alphas = (d_centered @ u) #<-- N X 1

var    = np.var(alphas)

dist   = d_centered - alphas@u.T
dist   = np.linalg.norm(dist,axis=1)
avg_dist = np.mean(dist)

for point,alpha in zip(d_centered, alphas):
    alpha = float(alpha[0])
    u0    = float(u[0,0])
    u1    = float(u[1,0])
    ax.plot( [alpha*u0, point[0] ] ,  [alpha*u1, point[1]  ], color="blue", alpha=0.1  )
    ax.scatter( point[0], point[1], color ="red" )
    
ax.plot( [-4,4], [-4,4], color="black" )  
ax.text(0.01,0.99,s="V = {:.1f}".format(var),ha="left",va="top",transform=ax.transAxes)
ax.text(0.99,0.01,s="D = {:.1f}".format(avg_dist),ha="right",va="bottom",transform=ax.transAxes)

ax.set_xlim(-4,4)
ax.set_ylim(-4,4)

ax.set_xlabel("Covariate 1 values centered")
ax.set_ylabel("Covariate 2 values centered")

fig.set_tight_layout(True)
fig.set_size_inches(6,6)

plt.show()

\begin{align} \mathbb{E}( || X_{n \times p} - \alpha_{n \times 1} u_{1 \times p} || ) \end{align}

(16)

\begin{align} \mathbb{E} \left[ ( X-\alpha u )( X-\alpha u )^{T} \right] = \mathbb{E} \left[ XX^{T} - 2X\alpha u^{T} + uu^{T} \right] \end{align}

(17)

\begin{align} x^{u}_{i} &= \alpha \cdot u \\ (x^{u}_{i} - u)^{T}x &= 0 \\ (x^{u}_{i})^{T}u - u^{T}x &= 0 \\ \alpha \cdot u^{T}u - u^{T}x &= 0 \\ \alpha &= \frac{u^{T}x}{u^{T}u} = \frac{x^{T}u}{u^{T}u} \\ \end{align}

(18)

\begin{align} Var(\alpha u ) &= Var\left[ (Xu)u \right]\\ &=\mathbb{E}\left( \left[ (Xu) u\right]^{T} \left[ (Xu) u\right]\right) - \left[\mathbb{E}( (Xu) u) \right]^{2} \\ &=\mathbb{E}\left( (Xu)^{2} \left[u\right]^{T} \left[u\right]\right) - 0^{2}\\ &=\mathbb{E}\left[(Xu)^{2}\right] \\ &=\mathbb{E}\left[\alpha^{2}\right] \\ \end{align}

(19)

\begin{align} X = \alpha u + (x - \alpha u)\\ || X || = || \alpha u || + || (X - \alpha u) ||\\ X^{T} X = (\alpha u)^{T}(\alpha u) + || (X - \alpha u) ||\\ X^{T} X = \alpha^{2} u^{T}u + || (X - \alpha u) ||\\ X^{T} X = \alpha^{2} + || (x_{i} - \alpha_{i}u) ||\\ \mathbb{E}( X^{T} X ) = \mathbb{E}\left[ \alpha^{2} \right] + \mathbb{E}\left[ || (x_{i} - \alpha_{i}u) ||\right] \\ \text{Var}(X) = \text{Var}( \alpha u ) + \mathbb{E}\left[ || (x_{i} - \alpha_{i}u) ||\right] \\ \end{align}

(20)

\begin{align} \text{Var}( \alpha u )\\ \text{Var}( (Xu) u )\\ \mathbb{E}\left[ ( (Xu) u )^{2} \right] \\ \mathbb{E}\left[ u^{T}u ( (Xu) )^{2} \right] \\ \mathbb{E}\left[ ( (Xu) )^{2} \right] \\ \mathbb{E}\left[ ( (Xu) )^{T}(Xu) \right] \\ \mathbb{E}\left[ u^{T} X^{T}Xu \right] \\ u^{T} \Sigma u \\ \end{align}

(21)


l,u         = np.linalg.eigh(COV)
l_est,u_est = np.linalg.eigh( d_centered.T @ d_centered/N )

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[11], line 2
      1 l,u         = np.linalg.eigh(COV)
----> 2 l_est,u_est = np.linalg.eigh( d_centered.T @ d_centered/N )

NameError: name 'N' is not defined


u      = np.array([1,1]).reshape(2,1)
u      = u/np.linalg.norm(u)

alphas = (d_centered @ u) #<-- N X 1

fig, ax = plt.subplots()

for point,alpha in zip(d_centered, alphas):
    alpha = float(alpha)
    u0    = float(u[0])
    u1    = float(u[1])
    ax.plot( [alpha*u0, point[0] ] ,  [alpha*u1, point[1]  ], color="blue", alpha=0.1  )
    ax.scatter( point[0], point[1] )
    
ax.plot( [-3,3], [-3,3], color="black" )  

ax.set_aspect("equal")
ax.set_xlim(-3,3)
ax.set_ylim(-3,3)
plt.show()

fig = plt.figure(figsize=(10, 5))
gs  = fig.add_gridspec(2, 3)

ax00 = fig.add_subplot(gs[0, 0])
ax10 = fig.add_subplot(gs[1, 0])

ax_proj = fig.add_subplot(gs[:, 1])  # joined second column
ax_box  = fig.add_subplot(gs[:, 2])  # joined third column

# ------------------------------------------------------------
# First row, first column: original data + projection to x-axis
# ------------------------------------------------------------
ax = ax00
ax.scatter(d_centered[:, 0], d_centered[:, 1])
ax.axhline(0, color="black", linewidth=1)
ax.vlines(x=d_centered[:, 0], ymin=0, ymax=d_centered[:, 1],colors="blue",alpha=0.2)
ax.set_xlim(-3, 3)
ax.set_ylim(-3, 3)
ax.set_aspect("equal")
ax.set_title("Projection onto x-axis")

# ------------------------------------------------------------
# Second row, first column: original data + projection to u
# ------------------------------------------------------------
u = np.array([1, 1], dtype=float).reshape(2, 1)
u = u / np.linalg.norm(u)

alphas = d_centered @ u          # N x 1
projs  = alphas @ u.T            # N x 2

ax = ax10
for point, proj in zip(d_centered, projs):
    ax.plot([proj[0], point[0]],[proj[1], point[1]],color="blue",alpha=0.1,)
    ax.scatter(point[0], point[1])

ax.plot([-3, 3], [-3, 3], color="black")
ax.set_xlim(-3, 3)
ax.set_ylim(-3, 3)
ax.set_aspect("equal")
ax.set_title(r"Projection onto u")

# ------------------------------------------------------------
# Joined second column: compare 1D projected coordinates
# ------------------------------------------------------------
x_proj = d_centered[:, 0]
u_proj = alphas.reshape(-1)

ax = ax_proj

ax.scatter(x_proj
           ,np.zeros_like(x_proj)
           ,facecolor="white",edgecolor="blue",alpha=0.5,label="x-axis")

ax.scatter(u_proj ,np.ones_like(u_proj) ,facecolor="white",edgecolor="darkorange",alpha=0.5,label="u")

ax.set_xlim(-4, 4)
ax.set_ylim(-1, 2)
ax.set_yticks([0, 1])
ax.set_yticklabels(["x-axis", r"$u$"])
ax.axvline(0, color="black", linewidth=1, alpha=0.4)
ax.set_xlabel("Projection coordinate")
ax.set_title("Projected points")

# ------------------------------------------------------------
projected_data = pd.DataFrame({
    "p": np.concatenate([x_proj, u_proj]),
    "projection": ( ["x-axis"] * len(x_proj) + ["u"] * len(u_proj) )})

ax = ax_box

sns.boxplot(x="p"
    ,y="projection"
    ,data=projected_data
    ,ax=ax
    ,width=0.35
    ,showcaps=True
    ,showfliers=False
    #,boxprops={"facecolor": "none", "edgecolor": "black", "linewidth": 1.5}
    #,whiskerprops={"color": "black", "linewidth": 1.2}
    #,capprops={"color": "black", "linewidth": 1.2}
    #,medianprops={"color": "black", "linewidth": 2}
)

sns.stripplot(x="p",y="projection",data=projected_data,ax=ax,jitter=0.18,size=4,alpha=0.55,color="black")

var_xaxis = np.var( x_proj )
var_proj  = np.var( u_proj )

ax.text( 0,1.25, s= "Var = {:.2f}".format(var_proj), va="bottom", ha="center" )
ax.text( 0,0.25, s= "Var = {:.2f}".format(var_xaxis), va="bottom", ha="center" )

ax.set_xlim(-4, 4)
ax.set_ylim(-1,2)
ax.set_xlabel("Projection coordinate")
ax.set_ylabel("")
ax.set_title("Summary of spread")

sns.despine(fig=fig)

plt.tight_layout()
plt.show()

np.var( d_centered[:,0] )

np.var(alphas)