svdGeometricInterpretation.Rmd

---
title: "2.3. Singular Value Decomposition - Geometric interpretation"
author: "Lieven Clement"
date: "statOmics, Ghent University (https://statomics.github.io)"
always_allow_html: true
---

<iframe width="560" height="315" src="https://www.youtube.com/embed/srUDBwP9IW8
" frameborder="0" style="display: block; margin: auto;" allow="autoplay; encrypted-media" allowfullscreen></iframe>

We introduce the geometric interpretation of the svd by using a toy example.

# Iris dataset

The iris dataset is a dataset on iris flowers.

- Three species (setosa, virginica and versicolor)
- Length and width of Sepal leafs
- Length and width of Petal Leafs

For didactical purposes we will use a subset of the data.

- Virginica Species
- 3 Variables: Sepal Length, Sepal Width, Petal Length
- This allows us to visualise the data in 3D plots
- Illustrate the data compression of the SVD from 3 to two dimensions.

## Subset the data

```{r}
library(tidyverse)
library(plotly)
irisSub <- iris %>%
  filter(Species == "virginica") %>%
  dplyr::select("Sepal.Length","Sepal.Width","Petal.Length")
```

```{r echo=FALSE}
p1 <- plot_ly(
    irisSub,
    x = ~Sepal.Width,
    y = ~Sepal.Length,
    z= ~Petal.Length) %>%
  add_markers(type="scatter3d") %>%
  layout(
    scene = list(
      aspectmode="cube",
      xaxis = list(range=c(-1,1)*max(abs(irisSub))), yaxis = list(range=c(-1,1)*max(abs(irisSub))), zaxis = list(range=c(-1,1)*max(abs(irisSub)))
      )
    )
p1
```

## Center the data

```{r}
X <- irisSub %>% scale(scale=FALSE)
```

The data is translated to a mean of [0, 0, 0].

```{r echo = FALSE}
p2 <- plot_ly(X %>% as.data.frame, x = ~Sepal.Length, y = ~Sepal.Width, z= ~Petal.Length) %>%
  add_markers(type="scatter3d") %>%
  layout(
    scene = list(
      aspectmode="cube",
      xaxis =  list(range=c(-1,1)*max(abs(irisSub))), yaxis = list(range=c(-1,1)*max(abs(irisSub))), zaxis = list(range=c(-1,1)*max(abs(irisSub)))
      )
    )
p2
```

We zoom in and add the original axis in grey in the origin.

```{r echo=FALSE}
p3 <- plot_ly(X %>% as.data.frame,
    x = ~Sepal.Length,
    y = ~Sepal.Width,
    z= ~Petal.Length) %>%
  layout(
    scene = list(
      aspectmode="cube",
      xaxis =  list(range=c(-1,1)*max(abs(irisSub))), yaxis = list(range=c(-1,1)*max(abs(irisSub))), zaxis = list(range=c(-1,1)*max(abs(irisSub)))
      )
    ) %>% layout(
    scene = list(aspectmode="cube",
    xaxis = list(range=c(-1,1)*max(abs(X))),
    yaxis = list(range=c(-1,1)*max(abs(X))),
    zaxis = list(range=c(-1,1)*max(abs(X))))) %>%
  add_trace(
    x = c(0,1),
    y = c(0,0),
    z = c(0,0),
    mode = "lines",
    line = list(width = 5, color = "grey"),
    type="scatter3d") %>%
  add_trace(
    x = c(0, 0),
    y = c(0, 1), z = c(0, 0),
    mode = "lines",
    line = list(
      width = 5,
      color = "grey"),
    type="scatter3d") %>%
  add_trace(
    x = c(0, 0),
    y = c(0, 0),
    z = c(0, 1),
    mode = "lines",
    line = list(width = 5, color  = "grey"),
    type = "scatter3d") %>%
hide_legend()
p3 %>%
  add_markers(
    type = "scatter3d",
    type = "markers",
    marker = list(color="#1f77b4")
    )
```

# SVD

1. We adopt the SVD on the centered data

```{r}
irisSvd <- svd(X)
```

2. We extract

  - the right singular vectors $\mathbf{V}$ and
  - the projections $\mathbf{Z}$

```{r}
V <- irisSvd$v
Z <- irisSvd$u %*% diag(irisSvd$d)
```

Note, that

- the SVD is essentially a rotation to a new coordinate system.
- we plotted $\mathbf{V}_3$ with dots because we will use the SVD for dimension reduction \[\text{3D} \rightarrow \text{2D}\]

```{r echo=FALSE}
p4 <- p3 %>%
  add_trace(
    x = c(0, V[1,1]),
    y = c(0, V[2,1]),
    z = c(0, V[3,1]),
    mode = "lines",
    line = list(width = 5,
      color = "red"),
    type="scatter3d") %>%
  add_trace(
    x = c(0, V[1,2]),
    y = c(0, V[2,2]),
    z = c(0, V[3,2]),
    mode = "lines",
    line = list(width = 5, color = "red"),
    type = "scatter3d") %>%
  add_trace(
    x = c(0, V[1,3]),
    y = c(0, V[2,3]),
    z = c(0, V[3,3]),
    mode = "lines",
    line = list(
      width = 5,
      color = "red",
      dash = "dash"),
    type="scatter3d") %>%
    hide_legend() %>%
    layout(
      scene = list(
        aspectmode="cube",
        annotations = list(
          list(
            x = V[1,1],
            y = V[2,1],
            z = V[3,1],
            text = "V1",
            textangle= 0,
            ax = 0,
            ay = 0,
            az = 0,
            xref='x',
            axref='x',
            yref = 'y',
            ayref = 'y',
            zref = 'z',
            azref='z',
            showarrow=FALSE
            ),
          list(
            x = V[1,2],
            y = V[2,2],
            z = V[3,2],
            text = "V2",
            textangle= 0,
            ax = 0,
            ay = 0,
            az = 0,
            xref='x',
            axref='x',
            yref = 'y',
            ayref = 'y',
            zref = 'z',
            azref='z',
            showarrow=FALSE
            ),
          list(
            x = V[1,3],
            y = V[2,3],
            z = V[3,3],
            text = "V3",
            textangle= 0,
            ax = 0,
            ay = 0,
            az = 0,
            xref='x',
            axref='x',
            yref = 'y',
            ayref = 'y',
            zref = 'z',
            azref='z',
            showarrow=FALSE)
          )
        )
      )
p4
p4 <- p4 %>%
  add_markers(
    type = "scatter3d",
    type = "markers",
    marker = list(color="#1f77b4")
    )
p4
```

Rotate the plot

- Note, that

  - V1 points in the direction of the largest variability in the data
  - V2 points in a direction orthogal on V1 pointing in the direction of the second largest variability in the data.


# Geometric Interpretation?

Write the **truncated SVD** as
\[
  \mathbf{X}_k = \mathbf{U}_k \boldsymbol{\Delta}_k \mathbf{V}_k^T = \mathbf{Z}_k \mathbf{V}_k^T
\]
with
\[
  \mathbf{Z}_k = \mathbf{U}_k \boldsymbol{\Delta}_k
\]
an $n \times k$ matrix.

Each of the $n$ rows of $\mathbf{Z}_k$, say $\mathbf{z}^T_{k,i}$, represents a point in a $k$-dimensional space.

```{r}
V2 <- V[,1:2]
Z2 <- Z[,1:2]
X2 <- Z2 %*% t(V2)
```

Because of the orthonormality of the singular vectors, we also have
\begin{eqnarray*}
  \mathbf{X}_k\mathbf{V}_k &=& \mathbf{Z}_k \mathbf{V}_k^T\mathbf{V}_k \\
  \mathbf{X}_k\mathbf{V}_k &=& \mathbf{Z}_k.
\end{eqnarray*}

Thus the matrix $\mathbf{V}_k$ is a **transformation matrix** that may be used to transform $\mathbf{X}_k$ into $\mathbf{Z}_k$, and $\mathbf{Z}_k$ into $\mathbf{X}_k$.

---

More importantly, it can be shown that (thanks to orthonormality of $\mathbf{V}$)
\[
\mathbf{X}\mathbf{V}_k = \mathbf{Z}_k.
\]

This follows from (w.l.g. rank($\mathbf{X}$)=$r$)
\begin{eqnarray*}
\mathbf{X}\mathbf{V}_k
&=& \mathbf{UDV}^T\mathbf{V}_k = \mathbf{UD}\begin{pmatrix}
\mathbf{v}_1^T \\
\vdots \\
\mathbf{v}_r^T
\end{pmatrix}
\begin{pmatrix}
\mathbf{v}_1 \ldots \mathbf{v}_k
\end{pmatrix} \\
&=& \mathbf{UDV}^T\mathbf{V}_k = \mathbf{UD}\begin{pmatrix}
1 & 0 & \ldots & 0 \\
0 & 1 & \ldots & 0 \\
\vdots & \vdots & \ddots & 0 \\
0 & 0 & \ldots & 1 \\
0 & 0 & \ldots & 0 \\
\vdots & \vdots & \vdots & \vdots \\
0 & 0 & \ldots & 0
\end{pmatrix} \
= \mathbf{U}_k\boldsymbol{\Delta}_k = \mathbf{Z}_k \end{eqnarray*}

The $p \times k$ matrix $\mathbf{V}_k$ acts as a transformation matrix: transforming $n$ points in a $p$ dimensional space to $n$ points in a $k$ dimensional space.

```{r}
Z2proj <- X %*% V2
range(Z2 - Z2proj)
```


## Projection of a single data point

- Zoom in to see projection.
- The projection is indicated for the blue point $X_{44}$ to the red point $X_{2,44}$ in the plane spaned by V2.

```{r echo=FALSE}
i <- 44
colnames(X2) <- colnames(X)
p7 <- p2 %>% layout(
  scene = list(
    xaxis = list(range=c(-1,1)*max(abs(X))),
    yaxis = list(range=c(-1,1)*max(abs(X))),
    zaxis = list(range=c(-1,1)*max(abs(X))))
    ) %>%
  add_trace(
    x = c(X[i, 1], Z[i, 1]*V[1, 1]),
    y = c(X[i, 2], Z[i, 1]*V[2, 1]),
    z = c(X[i, 3], Z[i, 1]*V[3, 1]),
    mode = "lines",
    line = list(
      width = 5,
      color = "black",
      dash = "dash"),
    type = "scatter3d") %>%
  add_trace(
    x = c(X[i, 1], Z[i, 2] * V[1, 2]),
    y = c(X[i, 2], Z[i, 2] * V[2, 2]),
    z = c(X[i, 3], Z[i, 2] * V[3, 2]),
    mode = "lines",
    line = list(
      width = 5,
      color = "black",
      dash = "dash"),
    type="scatter3d") %>%
  add_markers(type="scatter3d",mode="markers") %>%
  add_trace(
    x = c(0, V[1, 1]),
    y = c(0, V[2, 1]),
    z = c(0, V[3, 1]),
    mode = "lines",
    line = list(
      width = 5,
      color = "red"),
    type="scatter3d") %>%
    add_trace(
      x = c(0, V[1, 2]),
      y = c(0, V[2, 2]),
      z = c(0, V[3,2]),
      mode = "lines",
      line = list(
        width = 5,
        color = "red"),
      type="scatter3d") %>%
    add_trace(
      x = c(0, V[1, 3]),
      y = c(0, V[2, 3]),
      z = c(0, V[3, 3]),
      mode = "lines",
      line = list(
        width = 5,
        color = "red",
        dash = "dash"),
        type="scatter3d") %>%
    hide_legend() %>%
    layout(
      scene = list(
        aspectmode="cube",
        annotations = list(list(
          x = V[1,1],
          y = V[2,1],
          z = V[3,1],
          text = "V1",
          textangle= 0,
          ax = 0,
          ay = 0,
          az = 0,
          xref = 'x',
          axref = 'x',
          yref = 'y',
          ayref = 'y',
          zref = 'z',
          azref='z',
          showarrow=FALSE
          ),
        list(
          x = V[1,2],
          y = V[2,2],
          z = V[3,2],
          text = "V2",
          textangle= 0,
          ax = 0,
          ay = 0,
          az = 0,
          xref = 'x',
          axref = 'x',
          yref = 'y',
          ayref = 'y',
          zref = 'z',
          azref='z',
          showarrow=FALSE
          ),
        list(
          x = V[1,3],
          y = V[2,3],
          z = V[3,3],
          text = "V3",
          textangle= 0,
          ax = 0,
          ay = 0,
          az = 0,
          xref = 'x',
          axref = 'x',
          yref = 'y',
          ayref = 'y',
          zref = 'z',
          azref='z',
          showarrow=FALSE
          )
        )
      )
    ) %>%
  add_markers(type="scatter3d",mode="markers",marker=list(color="grey"),hoverinfo="none") %>%
add_trace(
  data = (X%>%as.data.frame)[i,],
  x = ~Sepal.Length,
  y = ~Sepal.Width,
  z = ~Petal.Length,
  type = "scatter3d",
  mode = "markers",
  marker = list(color="#1f77b4")
  )%>%
add_trace(
  data = (X2 %>% as.data.frame)[i,],
  x = ~Sepal.Length,
  y = ~Sepal.Width,
  z = ~Petal.Length,
  type = "scatter3d",
  mode = "markers",
  marker = list(color="red"))%>%
add_trace(
  x = c(X2[i, 1], Z[i, 1] * V[1, 1]),
  y = c(X2[i, 2], Z[i, 1] * V[2, 1]),
  z = c(X2[i, 3], Z[i, 1] * V[3, 1]),
  mode = "lines",
  line = list(
    width = 5,
    color = "red",
    dash = "dash"),
  type="scatter3d") %>%
add_trace(
  x = c(X2[i, 1], Z[i, 2] * V[1, 2]),
  y = c(X2[i, 2], Z[i, 2] * V[2, 2]),
  z = c(X2[i, 3], Z[i, 2] * V[3,2]),
  mode = "lines",
  line = list(
    width = 5,
    color = "red",
    dash = "dash"),
  type="scatter3d") %>%
add_trace(
  x = c(X2[i, 1], X[i, 1]),
  y = c(X2[i, 2], X[i, 2]),
  z = c(X2[i, 3], X[i, 3]),
  mode = "lines",
  line = list(
    width = 5,
    color = "orange",
    dash = "dash"),
  type="scatter3d")
p7
```

## Projection of all datapoints: project all rows of X on V2

- Zoom in first look orthonal via direction V2 (rotate until text V2 is viewed in the origin)
- Zoom in first look orthonal via direction V1 (rotate until text V1 is viewed in the origin)
- Note, that
  - V1 points in the direction of the largest variability in the data
  - V2 points in a direction orthogal on V1 pointing in the direction of the second largest variability in the data.

```{r echo=FALSE}
p8 <- p2 %>% layout(
  scene = list(
    xaxis = list(range=c(-1,1)*max(abs(X))),
    yaxis = list(range=c(-1,1)*max(abs(X))),
    zaxis = list(range=c(-1,1)*max(abs(X))))
    ) %>%
  add_trace(
    x = c(0, V[1, 1]),
    y = c(0, V[2, 1]),
    z = c(0, V[3, 1]),
    mode = "lines",
    line = list(
      width = 5,
      color = "red"),
    type="scatter3d") %>%
    add_trace(
      x = c(0, V[1, 2]),
      y = c(0, V[2, 2]),
      z = c(0, V[3,2]),
      mode = "lines",
      line = list(
        width = 5,
        color = "red"),
      type="scatter3d") %>%
    add_trace(
      x = c(0, V[1, 3]),
      y = c(0, V[2, 3]),
      z = c(0, V[3, 3]),
      mode = "lines",
      line = list(
        width = 5,
        color = "red",
        dash = "dash"),
        type="scatter3d") %>%
    hide_legend() %>%
    layout(
      scene = list(
        aspectmode="cube",
        annotations = list(list(
          x = V[1,1],
          y = V[2,1],
          z = V[3,1],
          text = "V1",
          textangle= 0,
          ax = 0,
          ay = 0,
          az = 0,
          xref = 'x',
          axref = 'x',
          yref = 'y',
          ayref = 'y',
          zref = 'z',
          azref='z',
          showarrow=FALSE
          ),
        list(
          x = V[1,2],
          y = V[2,2],
          z = V[3,2],
          text = "V2",
          textangle= 0,
          ax = 0,
          ay = 0,
          az = 0,
          xref = 'x',
          axref = 'x',
          yref = 'y',
          ayref = 'y',
          zref = 'z',
          azref='z',
          showarrow=FALSE
          ),
        list(
          x = V[1,3],
          y = V[2,3],
          z = V[3,3],
          text = "V3",
          textangle= 0,
          ax = 0,
          ay = 0,
          az = 0,
          xref = 'x',
          axref = 'x',
          yref = 'y',
          ayref = 'y',
          zref = 'z',
          azref='z',
          showarrow=FALSE
          )
        )
      )
    ) %>%
  add_trace(
    data = X2 %>% as.data.frame,
    x = ~Sepal.Length,
    y = ~Sepal.Width,
    z = ~Petal.Length,
    type = "scatter3d",
    mode = "markers",
    marker = list(color = "red"))

for (i in 1:nrow(X2))
p8 <- p8 %>%  add_trace(
  x = c(X2[i, 1], X[i, 1]),
  y = c(X2[i, 2], X[i, 2]),
  z = c(X2[i, 3], X[i, 3]),
  mode = "lines",
  line = list(
    width = 5,
    color = "orange",
    dash = "dash"),
  type="scatter3d")

p8
```

- Projection only.
- This clearly shows that the projected points X2 (X projected on V2) live in a two dimensional space

```{r echo = FALSE}
p9 <- plot_ly(
    X %>% as.data.frame,
    x = ~Sepal.Length,
    y = ~Sepal.Width,
    z= ~Petal.Length) %>%
  add_trace(
    data = X2 %>% as.data.frame,
    x = ~Sepal.Length,
    y = ~Sepal.Width,
    z = ~Petal.Length,
    type = "scatter3d",
    mode = "markers",
    marker = list(color="red")) %>%
  add_trace(
    x = c(0, 1),
    y = c(0, 0),
    z = c(0, 0),
    mode = "lines",
    line = list(
      width = 5,
      color = "grey",
      dash = "dash"),
    type="scatter3d") %>%
  add_trace(
    x = c(0, 0),
    y = c(0, 1),
    z = c(0, 0),
    mode = "lines",
    line = list(
      width = 5,
      color = "grey",
      dash = "dash"),
    type="scatter3d") %>%
  add_trace(
    x = c(0, 0),
    y = c(0, 0),
    z = c(0, 1),
    mode = "lines",
    line = list(
      width = 5,
      color = "grey",
      dash = "dash"),
    type="scatter3d") %>%
  add_trace(
    x = c(0, V[1, 1]),
    y = c(0, V[2, 1]),
    z = c(0, V[3, 1]),
    mode = "lines",
    line = list(
      width = 5,
      color = "red"),
    type="scatter3d") %>%
  add_trace(
    x = c(0, V[1, 2]),
    y = c(0, V[2, 2]),
    z = c(0, V[3, 2]),
    mode = "lines",
    line = list(
      width = 5,
      color = "red"),
    type="scatter3d") %>%
  hide_legend() %>%
  layout(
    scene = list(
      aspectmode="cube",
      annotations = list(list(
        x = V[1,1],
        y = V[2,1],
        z = V[3,1],
        text = "V1",
        textangle= 0,
        ax = 0,
        ay = 0,
        az = 0,
        xref = 'x',
        axref = 'x',
        yref = 'y',
        ayref = 'y',
        zref = 'z',
        azref='z',
        showarrow=FALSE
      ),
      list(
        x = V[1,2],
        y = V[2,2],
        z = V[3,2],
        text = "V2",
        textangle = 0,
        ax = 0,
        ay = 0,
        az = 0,
        xref = 'x',
        axref = 'x',
        yref = 'y',
        ayref = 'y',
        zref = 'z',
        azref='z',
        showarrow=FALSE
        )
      )
    )
  ) %>%
  layout(
    scene = list(
      aspectmode = "cube",
      xaxis = list(range = c(-1, 1) * max(abs(X))), yaxis = list(range=c(-1, 1) * max(abs(X))),
      zaxis = list(range=c(-1, 1)*max(abs(X))))
      )
p9
```