<table width = 900 border = 0><tr><td>
# <b>NMST539 | Cvičenie 6</b>

## <b>Metóda hlavných komponent | Teoretická časť</b>
### <font color = "#880000"><b>LS 2020/2021 |  06/04/21 | (online výuka)</b></font>
###### <a href = "cvicenie2021_6.Rmd">zdrojový Rmd súbor (kódovanie UTF8)</a>

**************

<b>Outline šiesteho cvičenia:</b> 

<ul>
<li>metóda hlavných komponent (PCA- Principal Component Analysis);</li>
<li>teoretické aspekty metódy hlavných komponent;</li>
<li>súvislosť hlavných komponent a lineárnej regresie;</li>
</ul><br>


<b>Odborná literatúra:</b><br>
<ul>
<li>Hardle, W. and Simar.L.: Applied Multivariate Statistical Analysis. Springer, 2015</li>
<li>Mardia, K., Kent, J., and Bibby, J.:Multivariate Analysis, Academic Press, 1979.</li>
</ul>

<br><br>

**************

### <font color = "#880000"><b>1. Metóda hlavných komponent / teoretické aspekty</b></font>

Metóda hlavných komponent (PCA) je mnohorozmerná (prevažne exploratívna) štatistická metóda založena na teoretickej variančnej-kovariančnej matici náhodného vektoru $\boldsymbol{X} \in \mathbb{R}^{p}$. V praxi sa samozrejme metóda hlavných komponent 
aplikuje na základe empricikého odhadu teoretickej variančnej-kovariančnej matice, ktorý získame na základe náhodného výberu $\boldsymbol{X}_1, \dots, \boldsymbol{X}_n$ z rovnakého rozdelenia, ako pôvodný náhodný vektor $\boldsymbol{X} \in \mathbb{R}^{p}$. 

Jednoducho môžeme hlavné komponenty definovať ako 'vhodné' lineárne kombinácie pôvodných zložiek náhodneho vektoru $\boldsymbol{X} = (X_1, \dots, X_p)^\top$. Hlavná myšlienka spočíva v definovaní/určení $\ell <  p$ hlavných komponent (lineárnych kombinácii zložiek $X_1, \dots, X_p$), ktoré využijeme namiesto $p \in \mathbb{N}$ pôvodných zložiek (tzv. redukcia dimenzionality). Nové zložky -- $\ell$ hlavných kompoment -- sú v určitom zmysle najlepšie možné. Hlavné kompomenty sú orgogonálne a navyše sú usporiadané vzhľadom k postupne klesajúcej variabilite jednotlivých kompoment. V praktických úlohach často požadujeme, aby $\ell \ll p$, čo predstavuje výraznú redukciu pôvodnej dimenzie. Na druhej strane sa ale snažíme zachovať čo najväčšiu mieru z pôvodnej informácie.  <br><br>

Z formálneho hľadiska je možné hlavne komponenty definovať rôzne: 

<ul>
<li>ortogonálne lineárne kombinácie zložiek náhodného vektoru s postupne klesajúcou variabilitou;</li>
<li>hlavné komponenty definujú najlepší lineárny podpriestor pre aproximáciu dat;</li>
<li>z geometrického hľadiska sa jedná o EIV rozklad teoretickej (alebo empirickej) variančnej-kovariančnej matice;</li>
<li>z pohľadu dat sa jedná o SVD rozklad $n \times p$ rozmernej datovej matice;</li>
</ul>

Z pohľadu štatistika sú hlavné komponenty definované ortogonálne lineárne kombinácie pôvodných premenných $(X_1, \dots, X_p)^\top$ v náhodnom vektore $\boldsymbol{X} \in \mathbb{R}^{p}$, tak aby tieto lineárne kombinácie boli vzájomne nekorelované, s nulovou strednou hodnotou a diagonálnou variančnou-kovariančnou maticou, ktorá ma na hlavnej diagonále vlastné čísla variančnej-kovariančnej matice náhodného vektoru $\boldsymbol{X} \in \mathbb{R}^{p}$. 

V tomto zmysle sa jedná o maximalizačnú úlohu, kde prvá hlavná komponenta $\boldsymbol{\delta}_1^\top \boldsymbol{X}$ (t.j. lineárná kombinácia zložiek vektora $\boldsymbol{X} \in \mathbb{R}^{p}$) je definovaná pre $\boldsymbol{\delta}_1 \in \mathbb{R}^p$ také, že
$$
\boldsymbol{\delta}_1 = Argmax_{\boldsymbol{\{\delta};~\|\boldsymbol{\delta}\|  = 1 \}} ~~Var (\boldsymbol{\delta}^\top \boldsymbol{X}). 
$$
Druhá hlavná komponenta $\boldsymbol{\delta}_2^\top \boldsymbol{X}$ je definovaná ako analogická maximalizačná uloha, ale pri dodatočnej podmienke, že sa jedná o ortogonálne lineárne kombinácie, t.j. $\boldsymbol{\delta}_1^\top\boldsymbol{\delta}_2  = 0$, resp. náhodné veličiny $\boldsymbol{\delta}_1^\top \boldsymbol{X}$ a $\boldsymbol{\delta}_2^\top \boldsymbol{X}$ sú nekorelované, t.j., že platí
$$
Cov(\boldsymbol{\delta}_1^\top \boldsymbol{X}, \boldsymbol{\delta}_2^\top \boldsymbol{X}) = 0.
$$
Na hlavné komponenty sa ale môžeme pozerať aj z geometrického pohľadu projekcií do lineárných podpriestorov a hlavne konceptu lineárnej regresie (prípadne tzv. totálne najmenších štvorocov -- TLS / Total Least Squares).

<br><br>

### <font color = "#880000"><b>2. Hlavné komponenty, ortogonálne projekcie a regresia</b></font>

Nech  $\boldsymbol{Y} \in \mathbb{R}^n$ je nejaký náhodný vektor  dĺžky $n \in \mathbb{N}$ a nech $\mathbb{X}$ predstavuje (náhodnú) maticu napr. vysvetlujúcich premenných, ktorá je typu $n \times p$. Vo všeobecnosti predpokládame, že $p \in \mathbb{N}$ (väčšinou počet rôznych premenných) je menšie, ako $n$ (väčšinou počet pozorovaní). Klasická úloha lineárnej regresie spočíva vo využití informácie obsiahnutej v matici $\mathbb{X}$ k vhodnej aproximácii $n$ rozmerného náhodného vektoru $\boldsymbol{Y}$ (tzv. vektor závislej premennej) 
v $p$ rozmernom lineárnom podpriestore $\mathcal{L}(\mathbb{X})$ -- lineárnom podpriestore $R^n$ generovanom stĺpcami matice $\mathbb{X}$. 

Ak označíme ortogonálnu projekciu náhodného vektoru $\boldsymbol{Y}$ do $p$ rozmerného lineárneho podpriestoru $\mathcal{L}(\mathbb{X})$ ako $P_{\mathbb{X}}(\boldsymbol{Y})$, tak potom klasický lineárny regresný model pre vektor závislej premennej $\boldsymbol{Y}$
a nezávisle premenné v stĺpcoch matice $\mathbb{X}$ vo forme $P_{\mathbb{X}}(\boldsymbol{Y}) = a_1 \boldsymbol{X}_1 + \dots + a_p \boldsymbol{X}_p$ (lineárna kombinácia stĺpcov matice $\mathbb{X}$) dáva najmenšiu možnú štvorcovú chybu v zmysle, že výraz 

$$
\| \boldsymbol{Y} - P_{\mathbb{X}}(\boldsymbol{Y}) \|_2^2
$$
je minimalizovaný (pripomeňte si Gauss-Markovovu vetu pre tzv. BLUE (Best Linear Unbiased Estimate)  odhad v lineárnej regresii). <br><br>

Príslušná (idempotentná) projekčná matica je definovaná predpisom
$$
\mathbb{H} = \mathbb{X}(\mathbb{X}^\top \mathbb{X})^{-1}\mathbb{X}^\top = \mathbb{Q}\mathbb{Q}^\top,
$$
kde matica $\mathbb{X}$ obsahuje vektory (v stĺpcoch) generujúce lineárný podpriestor $\mathcal{L}(\mathbb{X})$ 
resp. stĺce matice $\mathbb{Q}$ predstavujú ortogonálnu bázu lineárneho podpriestoru $\mathcal{L}(\mathbb{X})$ . 
Samozrejme platí, že 
$$
P_{\mathbb{X}}(\boldsymbol{Y}) = \mathbb{H}\boldsymbol{Y} = \mathbb{X}(\mathbb{X}^\top \mathbb{X})^{-1}\mathbb{X}^\top \boldsymbol{Y} = \mathbb{Q}\mathbb{Q}^\top \boldsymbol{Y} = P_{\mathbb{Q}}(\boldsymbol{Y})
$$

<br><br>

Koncept lineárnej regresie (t.j. ortogonálnej projekcie do lineárneho podpriestoru generovaného stĺpcami matice $\mathbb{X}$)
môže byť zobecnený na problém hľadania najlepšej ortogonálnej projekcie $\boldsymbol{Y}$ do lineárneho podpriestoru, ktorý je generovaný stĺpcami nejakej obecnej matice  $\mathbb{B}$ typu $n \times p$, tak že platí $\mathbb{B}^\top\mathbb{B} = \mathbb{I}$.

Označme obecne takúto projekciu ako $P_{\mathbb{B}}(\boldsymbol{Y})$. Potom platí, že $P_{\mathbb{B}}(\boldsymbol{Y}) = \mathbb{B}\mathbb{B}^\top\boldsymbol{Y}|$ a, navyše, lze ukázať, že zároveň platí aj 

$$
\|\boldsymbol{Y} - P_{\mathbb{\Gamma}}(\boldsymbol{Y})\|_2^2 \leq \|\boldsymbol{Y} - P_\mathbb{B}(\boldsymbol{Y})\|_2^2
$$
kde $\Gamma$ je matica typu $n\times p$ ktorá obsahuje (v stĺpcoch) prvých $p$ vlastných vektorov teoretickej variančnej-kovariančnej matice $\Sigma$ náhodného vektoru $\boldsymbol{Y}$ (tzv. Eckart a Young (1936) a Mirsky (1960) approximácia).
Pripomeňme, že v súvislosti s konceptom lineárnej regresie -- t.j. projekciou $P_{\mathbb{X}}(\boldsymbol{Y})$ -- 
predpokládame vlasnosť linearity (B<font color = "red">L</font>UE). <br><br>

V nasledujúcom príklade použijeme tri konkrétne ortogonálne projekcie náhodného vektoru $\boldsymbol{Y}$ do lineárneho podpriestoru, 
ktorý je generovaný stĺpcami vhodnej matice typu $n \times p$. Rovnaké data -- t.j. náhodný vektor $\boldsymbol{Y}$ -- použijeme pre všetky projekcie, ale lineárny podpriestor, do ktorého bude vektor projektovaný, bude zakaždým iný. Výsledné ortogonálne projekcie sa budu preto vájomne líšiť. Ktorá ortogonálna projekcia je najlepšia (a v akom zmysle)?


```{r, fig.width = 10, fig.height=4.5}
set.seed(1234)
x <- seq(0,1, length = 100)
y <- 0 + 2 * x + rnorm(100)

par(mfrow = c(1,3))
### regression of Y on X
plot(y ~ x, pch = 21, bg = "lightblue", xlab = "X variable", ylab = "Y variable", main = "Projection of Y onto L(X)", ylim = c(-1.5, 4))
m <- lm(y ~ x)
abline(m, col = "red", lwd = 2)

index <- c(10, 60, 80)

for (i in 1:100){
  lines(c(x[i], x[i]), c(y[i], m$fitted[i]), col = "blue", lty =2)
}


### regression of X on Y 
plot(y ~ x, pch = 21, bg = "lightblue", xlab = "X variable", ylab = "Y variable", main = "Projection of X onto L(Y)", ylim = c(-1.5, 4))
m <- lm(x ~ y)
lines((x - m$coeff[1])/m$coeff[2] ~ x, col = "red", lwd = 2)

index <- c(10, 60, 80)

for (i in 1:100){
  lines(c(x[i], m$coeff[1] + m$coeff[2] * y[i]), c(y[i], y[i]), col = "blue", lty =2)
}


### regression onto the first PC

### data standardization
means <- apply(cbind(x,y), 2, mean)
sdErs <- apply(cbind(x,y), 2, sd)
cd <- cbind((x - means[1])/sdErs[1], (y - means[2])/sdErs[2])

### first eigen vector
firstEigenVector <- eigen(cov(cd))$vectors[,1]

### projection onto the linear space generated by the first eigen vector
project2pc1 <- t(firstEigenVector) %*% t(cd)

### rotation back to the original coordinates
rotation <- as.numeric(project2pc1) %*% t(firstEigenVector)

### rescaling to the original mean/variance
rescaled <- cbind(rotation[,1] * sdErs[1] + means[1], rotation[,2] * sdErs[2] + means[2])

plot(y ~ x, pch = 21, bg = "lightblue", xlab = "X variable", ylab = "Y variable", main = "Projection onto 1.PC", ylim = c(-1.5,4), xlim = c(-0, 1))

### first PC
lines(rescaled[,2] ~ rescaled[,1], col = "red", lwd = 2)

### orthogonal projections 
for (i in 1:100){
  lines(c(x[i], rescaled[i,1]), c(y[i], rescaled[i,2]), col = "blue", lty =2)
}
```


Na základe predchádzajúceho príkladu môžeme obecne povedať, že pri hľadaní hlavných komponent sa dá namiesto spektrálnej dekompozície (EIV) teoretickej variančnej-kovariančnej matice náhodného vektoru $\boldsymbol{Y}$ postupovať vrámci konceptu klasickej lineárnej regresie
a iteratívnym postupom (alternáciou medzi dvoma podobnými lineárnymi regresnými úlohami) sa dopracovať k ekvivalentnému riešeniu -- hlavným komponentám. Pre háhodný vektor závislej premennej $\boldsymbol{Y} \in \mathbb{R}^n$  múžeme tieto dve úlohy formulovať následovne: 

$$
\|\boldsymbol{Y} - P_\mathbb{B}(\boldsymbol{Y})\|_2^2 =  \|\boldsymbol{Y} - \mathbb{B}\mathbb{B}^\top\boldsymbol{Y}\|_2^2 =   \|\boldsymbol{Y} - \mathbb{B} \boldsymbol{v}\|_2^2 = \sum_{i = 1}^{n} [Y_{i} - \boldsymbol{b}_i^\top\boldsymbol{v}]^2,  
$$

kde $\boldsymbol{v} = \mathbb{B}^\top\boldsymbol{Y}$, pre $i = 1, \dots, n$, a  $\boldsymbol{b}_i$ je príslušný řádek matice $\mathbb{B}$.
Musíme si uvedomiť, že ani matica $\mathbb{B}$ ani vektor $\boldsymbol{v}$  nie sú známe (preto potrebujeme riešiť dve vzálomne alternujúce regresné problémy).<br><br>

Výraz 
$$
\sum_{i = 1}^{n} [Y_{i} - \boldsymbol{b}_i^\top\boldsymbol{v}]^2
$$
budeme preto minimalizovať aj vzhľadom k neznámej ortogonálnej matici $\mathbb{B}$ (v jednom kroku), tzn., že platí $\mathbb{B}^\top\mathbb{B} = \mathbb{I}$ a tiež vzhľadom k neznámemu vektoru $\boldsymbol{v}$ (v následujúcom kroku).
Oba kroky dostatočne dlho alternujeme, až kým nie je dosiahnutá konvergencia. <br><br>

Hlavná myšlienka samozrejme spočíva v tom, že sa snazíme nájsť najlepšiu ortogonálnu projekciu do lineárneho podpriestoru, ktorý je generovaný orgogonálnou maticou $\mathbb{B}$, tak, aby súčet štvorcových chýb bol najmenší možný.  Z teórie lineárnej regresie vieme, 
že za platnosti dodatočného predpokladu linearity je najlepším odhadom práve model lineárnej regresie. 
V obecnom prípade (bez predpokladu linerity) ale lze získať lepšiu ortogonálnu projekciu. Ak bude matica $\mathbb{B}$ v jednotlivých stĺpcoch obsahovať prvých $p \in \mathbb{N}$ hlavných kompoment, tak získame najlepšiu možnú ortogonálnu projekciu náhodného vektoru $\boldsymbol{Y}$ do $p$ rozmerného lineárneho podpriestoru (samozrejme všetko pouze v zmysle najmenšej štvorcovej chyby). 

Jednoduchý príklad alternujúcej regresie od <a href = "https://www.stat.ubc.ca/~matias/">Matíasa Salibiána Barreru</a>. 

```{r}
alter.pca.k1 <- function(x, max.it = 500, eps=1e-10) {
  n2 <- function(a) sum(a^2)
  p <- dim(x)[2]
  x <- scale(x, scale=FALSE)
  it <- 0
  old.a <- c(1, rep(0, p-1))
  err <- 10*eps
  while( ((it <- it + 1) < max.it) & (abs(err) > eps) ) {
    b <- as.vector( x %*% old.a ) / n2(old.a)
    a <- as.vector( t(x) %*% b ) / n2(b)
    a <- a / sqrt(n2(a))
    err <- sqrt(n2(a - old.a))
    old.a <- a
  }
  conv <- (it < max.it)
  return(list(a=a, b=b, conv=conv))
}


set.seed(1234)
n <- 20
p <- 5
x <- matrix(rt(n*p, df=2), n, p)
```

Prvú hlavnú kompomentu pomocou alternujúcich lineárnych regresných modelov získame ako 

```{r}
alter.pca.k1(x)$a
```

Jedná sa o stĺpcový vektor matice $\mathbb{B}$, ktorý nam dá minimálnú štvorcovú chybu projekcie 
do jednorozmerného lineárneho podpriestoru, ktorý je generovaný prvým vlastnym vektorom variančnej-kovariančnej matice -- 
čo múžeme priamo porovnať s využitím EIV alebo SVD rozkladu: 


```{r}
svd(cov(x))$u[,1]
eigen(cov(x))$vectors[,1]
```

Je nutné si uvedomiť, že rozdiel v znamienku je v tomto prípade irelevantný. 

<br><br>

Z hľadiska praktických aplikacii môže byť užitočné porovnať fungovanie a časovú náročnosť oboch postupov napr. pomocou simulačnej štúdie. 
Využijeme výrazne vyšší počet dimenzii, aby bolo porovnanie viac relevantné a výpovedné. 

```{r}
n <- 2000
p <- 500
x <- matrix(rt(n*p, df=2), n, p)

system.time( tmp <- alter.pca.k1(x) )

system.time( e1 <- svd(cov(x))$u[,1] )
```

Výpočetná  časová náročnosť sa zdá byť v prospech alternujúcej regresie. Obecne platí, že v programe R je nutné dbať 
na maticové operácie a pristupovať k ním efektívne, čo do časovej náročnosti výpočtu. 

<br><br>

### <font color = "#880000"><b>3. Metóda hlavných komponent v programe R</b></font>

Metóda hlavných komponent je v programe R implementovaná pomocou príkazu `princomp()` 
ktorý využíva spektrálný rozklad (EIV) teoretickej (resp. empirickej) variančnej-kovariančnej matice. 
Pre konkrétne data $\mathcal{X} = (\boldsymbol{X}_1, \dots, \boldsymbol{X}_p)$ s odhadnutou variančnou-kovariančnou maticou 
$\mathcal{S}$ získame spektrálny rozklad -- t.j. $EIV(\mathcal{S})$ (analogicky, ak je známa teoretická matica, získame $EIV(\mathcal{\Sigma})$). 

Alternatívne a pri praktických (t.j. empirických) úlohach aj ekvivalentne je možné v programe R využiť aj príkaz `prcomp()`, 
ktorý funguje výhradne na empirickom princípe a využíva singulárný rozklad (SVD)
datovej matice $\mathcal{X}$ -- teda $SVD(\mathcal{X})$). 

Označme oba rozklady ako
$$
EIV(\mathcal{S}) = \Gamma \Lambda \Gamma^\top \quad \quad \textrm{a analogicky} \quad \quad SVD(\mathcal{X}) = UDV^\top.
$$

Ekvivalentnosť oboch rozkladov je viac-menej okamžitá: 
$$
n \mathcal{S} = \mathcal{X}^\top\mathcal{X} =  VDU^\top UD V^\top = V D^2 V^\top = \Gamma \widetilde{\Lambda} \Gamma^\top = EIV(n\mathcal{S}),
$$ 
pre $\Gamma \equiv V$ a $\widetilde{\Lambda} \equiv D^2$. Taktiež platí, že
$$EIV(n\mathcal{S}) = \Gamma \widetilde{\Lambda} \Gamma^\top = \Gamma (n \Lambda) \Gamma^\top = n \Gamma \Lambda \Gamma^\top = n EIV(\mathcal{S}).$$

<br><br>

Využitie metódy hlavných komponent z praktického hľadiska (prevažne ako exploratívny a dimenziu redukujúci nástroj)
bude predmetom následujúceho cvičenia. <br><br>


#### <font color = "#880000">Samostatne</font> 
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
Pre fungovanie hlavných komponent je podstatné porozumenie fungovania vlastných čísel a vlastných vektorov. 
Interaktívne grafické ilustrácie pre pripomenutie je možné nájsť napr. tu: 

<ul>
<li>Vlastné čísla a vlastné vektory matice:<br>
<a href = "http://setosa.io/ev/eigenvectors-and-eigenvalues/">http://setosa.io/ev/eigenvectors-and-eigenvalues/</a></li>

<li>Hlavné komponenty v 2D a v 3D:<br>
<a href = "http://setosa.io/ev/principal-component-analysis/">http://setosa.io/ev/principal-component-analysis/</a></li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>


<br><br>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<hr size="3" width="100%" noshade style="color:#440000" align="left" /><br><br> 

## <font color = "#880000"><b>Domáca (samostatná) úloha </b></font>
### (Deadline: Siedme cvičenie / 13.04.2021)

<hr size="3" width="100%" noshade style="color:#440000" align="left" />
Uvažujte náhodný výber $(X_1,Y_1)^\top, \dots (X_n,Y_n)^\top$ generovaný, resp. simulovaný z nejakého vhodného dvojrozmerného rozdelenia. Náhodný vektor $\boldsymbol{Y} = (Y_1, \dots, Y_n)^\top$ lze chápať aj ako element v $n$ rozmernom lineárnom priestore $\mathbb{R}^n$. 
Budu nás zaujímať tri rôzne ortogonálne projekcie -- zakaždým do iného dvojrozmerného lineárneho podpriestoru v $\mathbb{R}^n$. Pre každú projekciu explicitne spočítajte štvorcové chyby  (t.j. empiricky overte kvalitu danej projekcie vhľadom k teoretickým očakávaniam). 
Projekcie sa pokúste nejak vizualizovať. <br><br>

<ul>
<li>Zvoľte si ľubovolný dvojrozmerný lineárny podpriestor generovaný stĺpcovými vektormi vhodnej matice typu $n \times 2$ (ale uvažujte obecne iný lineárny podpierstor, ako v nasledujúcich dvoch úlohach). Pomocou projekčnej matice data ortogonálne projektujte do 
vami zvoleného lineárneho podpriestoru. </li>

<li>Využijte metodológiu lineárnej regresie a aplikujte ortogonálnu projekciu na regresnú přímku, t.j. do lineárneho podpriestoru 
ktorý je generovaný stĺpcami regresnej matice $\mathbb{X} = (\boldsymbol{1}, \boldsymbol{X})$, kde prvý stĺpec je 
stĺpec jednotiek a druhý stĺpec je $\boldsymbol{X} = (X_1, \dots, X_n)^\top$.</li>

<li>Nájdite ortogonálnu projekciu do dvojrozmerného lineárneho podpriestoru, ktorý je daný prvou a druhou principálnou kompomentou.
Pokúste sa túto úlohu riešiť bez pomoci príkazov `princomp()` a `prcomp()`.</li>

<li> Riešenie umiestnite na svoju webovú stránku, <font color = "red">najneskôr v utorok, 13.04.2021, do 14:00</font>.</li>
</ul>
<hr size="3" width="100%" noshade style="color:#440000" align="left" />
<br><br>



</td></tr></table>
<br><br><br>