1 Data

Part of the iris dataset

X <- iris[1:5,1:4] %>% as.matrix
X
#>   Sepal.Length Sepal.Width Petal.Length Petal.Width
#> 1          5.1         3.5          1.4         0.2
#> 2          4.9         3.0          1.4         0.2
#> 3          4.7         3.2          1.3         0.2
#> 4          4.6         3.1          1.5         0.2
#> 5          5.0         3.6          1.4         0.2

2 Centering

\[ \mathbf{H} = \mathbf{I} - \frac{1}{n} \mathbf{1}\mathbf{1}^T , \] \[ \mathbf{X}_c = \mathbf{H}\mathbf{X} \] \[ \mathbf{H}\mathbf{X}_c = \mathbf{X}_c \]

H <- diag(nrow(X)) - matrix(1/nrow(X),nrow=nrow(X),ncol=nrow(X))
Xc <- H%*%X
colMeans(Xc)
#>  Sepal.Length   Sepal.Width  Petal.Length   Petal.Width 
#>  2.220446e-16  0.000000e+00  1.576517e-17 -8.326673e-18
H%*%Xc-Xc
#>       Sepal.Length   Sepal.Width  Petal.Length  Petal.Width
#> [1,] -2.220446e-16  5.551115e-17 -1.521006e-17 8.326673e-18
#> [2,] -2.289835e-16  0.000000e+00 -1.521006e-17 8.326673e-18
#> [3,] -2.498002e-16 -1.387779e-17 -2.775558e-17 8.326673e-18
#> [4,] -2.220446e-16 -2.775558e-17 -1.387779e-17 8.326673e-18
#> [5,] -2.220446e-16  0.000000e+00 -1.756928e-17 8.326673e-18

3 Gram matrix

Gram matrix: \(\mathbf{G}=\mathbf{X}\mathbf{X}^T\) Here we work on centered data so \(\mathbf{G}=\mathbf{X}_c\mathbf{X}_c^T\)

G <- Xc%*%t(Xc)

4 Squared Distance matrix

\[ \mathbf{D}_X = \mathbf{N} - 2\mathbf{X}\mathbf{X}^T + \mathbf{N}^T , \]

N <- matrix(diag(G),nrow(Xc),nrow(Xc))
N
#>       [,1]  [,2]  [,3]  [,4]  [,5]
#> [1,] 0.106 0.106 0.106 0.106 0.106
#> [2,] 0.080 0.080 0.080 0.080 0.080
#> [3,] 0.042 0.042 0.042 0.042 0.042
#> [4,] 0.110 0.110 0.110 0.110 0.110
#> [5,] 0.122 0.122 0.122 0.122 0.122
dist2 <- N-2*G+t(N)
dist2
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,] 0.00 0.29 0.26 0.42 0.02
#> [2,] 0.29 0.00 0.09 0.11 0.37
#> [3,] 0.26 0.09 0.00 0.06 0.26
#> [4,] 0.42 0.11 0.06 0.00 0.42
#> [5,] 0.02 0.37 0.26 0.42 0.00
dist(Xc)^2
#>      1    2    3    4
#> 2 0.29               
#> 3 0.26 0.09          
#> 4 0.42 0.11 0.06     
#> 5 0.02 0.37 0.26 0.42

6 Show that N cancels out when multiplying with H

\[ \mathbf{H N H} = \mathbf{0}\]

\[ \mathbf{H N}^T \mathbf{H} = \mathbf{0}\]

N%*%H
#>               [,1]          [,2]          [,3]         [,4]          [,5]
#> [1,]  9.048318e-19  9.048318e-19 -6.034062e-18 9.048318e-19 -3.619327e-18
#> [2,] -8.520962e-18 -1.582068e-18  5.356826e-18 1.887379e-18  6.328271e-18
#> [3,] -4.662937e-19  3.003153e-18  3.003153e-18 1.268430e-18  1.865175e-18
#> [4,]  8.604228e-19  8.604228e-19  8.604228e-19 8.604228e-19 -3.441691e-18
#> [5,]  7.271961e-19  7.271961e-19  7.271961e-19 7.271961e-19 -2.908784e-18
H%*%t(N)
#>               [,1]          [,2]          [,3]          [,4]          [,5]
#> [1,]  9.048318e-19 -8.520962e-18 -4.662937e-19  8.604228e-19  7.271961e-19
#> [2,]  9.048318e-19 -1.582068e-18  3.003153e-18  8.604228e-19  7.271961e-19
#> [3,] -6.034062e-18  5.356826e-18  3.003153e-18  8.604228e-19  7.271961e-19
#> [4,]  9.048318e-19  1.887379e-18  1.268430e-18  8.604228e-19  7.271961e-19
#> [5,] -3.619327e-18  6.328271e-18  1.865175e-18 -3.441691e-18 -2.908784e-18

Session info

Session info
#> [1] "2024-10-07 12:47:20 CEST"
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value
#>  version  R version 4.4.0 RC (2024-04-16 r86468)
#>  os       macOS Big Sur 11.6
#>  system   aarch64, darwin20
#>  ui       X11
#>  language (EN)
#>  collate  en_US.UTF-8
#>  ctype    en_US.UTF-8
#>  tz       Europe/Brussels
#>  date     2024-10-07
#>  pandoc   3.1.1 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/ (via rmarkdown)
#> 
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package     * version date (UTC) lib source
#>  bookdown      0.40    2024-07-02 [1] CRAN (R 4.4.0)
#>  bslib         0.8.0   2024-07-29 [1] CRAN (R 4.4.0)
#>  cachem        1.1.0   2024-05-16 [1] CRAN (R 4.4.0)
#>  cli           3.6.3   2024-06-21 [1] CRAN (R 4.4.0)
#>  colorspace    2.1-1   2024-07-26 [1] CRAN (R 4.4.0)
#>  digest        0.6.37  2024-08-19 [1] CRAN (R 4.4.1)
#>  dplyr       * 1.1.4   2023-11-17 [1] CRAN (R 4.4.0)
#>  evaluate      1.0.0   2024-09-17 [1] CRAN (R 4.4.1)
#>  fansi         1.0.6   2023-12-08 [1] CRAN (R 4.4.0)
#>  fastmap       1.2.0   2024-05-15 [1] CRAN (R 4.4.0)
#>  forcats     * 1.0.0   2023-01-29 [1] CRAN (R 4.4.0)
#>  generics      0.1.3   2022-07-05 [1] CRAN (R 4.4.0)
#>  ggplot2     * 3.5.1   2024-04-23 [1] CRAN (R 4.4.0)
#>  glue          1.8.0   2024-09-30 [1] CRAN (R 4.4.1)
#>  gtable        0.3.5   2024-04-22 [1] CRAN (R 4.4.0)
#>  hms           1.1.3   2023-03-21 [1] CRAN (R 4.4.0)
#>  htmltools     0.5.8.1 2024-04-04 [1] CRAN (R 4.4.0)
#>  jquerylib     0.1.4   2021-04-26 [1] CRAN (R 4.4.0)
#>  jsonlite      1.8.9   2024-09-20 [1] CRAN (R 4.4.1)
#>  knitr         1.48    2024-07-07 [1] CRAN (R 4.4.0)
#>  lifecycle     1.0.4   2023-11-07 [1] CRAN (R 4.4.0)
#>  lubridate   * 1.9.3   2023-09-27 [1] CRAN (R 4.4.0)
#>  magrittr      2.0.3   2022-03-30 [1] CRAN (R 4.4.0)
#>  munsell       0.5.1   2024-04-01 [1] CRAN (R 4.4.0)
#>  pillar        1.9.0   2023-03-22 [1] CRAN (R 4.4.0)
#>  pkgconfig     2.0.3   2019-09-22 [1] CRAN (R 4.4.0)
#>  purrr       * 1.0.2   2023-08-10 [1] CRAN (R 4.4.0)
#>  R6            2.5.1   2021-08-19 [1] CRAN (R 4.4.0)
#>  readr       * 2.1.5   2024-01-10 [1] CRAN (R 4.4.0)
#>  rlang         1.1.4   2024-06-04 [1] CRAN (R 4.4.0)
#>  rmarkdown     2.28    2024-08-17 [1] CRAN (R 4.4.0)
#>  rstudioapi    0.16.0  2024-03-24 [1] CRAN (R 4.4.0)
#>  sass          0.4.9   2024-03-15 [1] CRAN (R 4.4.0)
#>  scales        1.3.0   2023-11-28 [1] CRAN (R 4.4.0)
#>  sessioninfo   1.2.2   2021-12-06 [1] CRAN (R 4.4.0)
#>  stringi       1.8.4   2024-05-06 [1] CRAN (R 4.4.0)
#>  stringr     * 1.5.1   2023-11-14 [1] CRAN (R 4.4.0)
#>  tibble      * 3.2.1   2023-03-20 [1] CRAN (R 4.4.0)
#>  tidyr       * 1.3.1   2024-01-24 [1] CRAN (R 4.4.0)
#>  tidyselect    1.2.1   2024-03-11 [1] CRAN (R 4.4.0)
#>  tidyverse   * 2.0.0   2023-02-22 [1] CRAN (R 4.4.0)
#>  timechange    0.3.0   2024-01-18 [1] CRAN (R 4.4.0)
#>  tzdb          0.4.0   2023-05-12 [1] CRAN (R 4.4.0)
#>  utf8          1.2.4   2023-10-22 [1] CRAN (R 4.4.0)
#>  vctrs         0.6.5   2023-12-01 [1] CRAN (R 4.4.0)
#>  withr         3.0.1   2024-07-31 [1] CRAN (R 4.4.0)
#>  xfun          0.47    2024-08-17 [1] CRAN (R 4.4.0)
#>  yaml          2.3.10  2024-07-26 [1] CRAN (R 4.4.0)
#> 
#>  [1] /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library
#> 
#> ──────────────────────────────────────────────────────────────────────────────
LS0tCnRpdGxlOiAiTURTOiBMaW5rIFNxdWFyZWQgRGlzdGFuY2UgTWF0cml4IGFuZCBHcmFtIE1hdHJpeCBpbiBQcmFjdGljZSIKYXV0aG9yOiAiTGlldmVuIENsZW1lbnQiCmRhdGU6ICJzdGF0T21pY3MsIEdoZW50IFVuaXZlcnNpdHkgKGh0dHBzOi8vc3RhdG9taWNzLmdpdGh1Yi5pbykiCi0tLQoKYGBge3IsIGNoaWxkPSJfc2V0dXAuUm1kIn0KYGBgCgpgYGB7ciBlY2hvPUZBTFNFfQpsaWJyYXJ5KHRpZHl2ZXJzZSkKYGBgCgojIERhdGEKClBhcnQgb2YgdGhlIGlyaXMgZGF0YXNldAoKYGBge3J9ClggPC0gaXJpc1sxOjUsMTo0XSAlPiUgYXMubWF0cml4ClgKYGBgCgojIENlbnRlcmluZwoKXFsKICBcbWF0aGJme0h9ID0gXG1hdGhiZntJfSAtIFxmcmFjezF9e259IFxtYXRoYmZ7MX1cbWF0aGJmezF9XlQgLApcXQpcWwogXG1hdGhiZntYfV9jID0gXG1hdGhiZntIfVxtYXRoYmZ7WH0KXF0KXFsKXG1hdGhiZntIfVxtYXRoYmZ7WH1fYyA9IFxtYXRoYmZ7WH1fYwpcXQoKYGBge3J9CkggPC0gZGlhZyhucm93KFgpKSAtIG1hdHJpeCgxL25yb3coWCksbnJvdz1ucm93KFgpLG5jb2w9bnJvdyhYKSkKWGMgPC0gSCUqJVgKY29sTWVhbnMoWGMpCkglKiVYYy1YYwpgYGAKCiMgR3JhbSBtYXRyaXgKCioqR3JhbSBtYXRyaXgqKjogJFxtYXRoYmZ7R309XG1hdGhiZntYfVxtYXRoYmZ7WH1eVCQKSGVyZSB3ZSB3b3JrIG9uIGNlbnRlcmVkIGRhdGEgc28KJFxtYXRoYmZ7R309XG1hdGhiZntYfV9jXG1hdGhiZntYfV9jXlQkCgpgYGB7cn0KRyA8LSBYYyUqJXQoWGMpCmBgYAoKIyBTcXVhcmVkIERpc3RhbmNlIG1hdHJpeAoKXFsKICBcbWF0aGJme0R9X1ggPSBcbWF0aGJme059IC0gMlxtYXRoYmZ7WH1cbWF0aGJme1h9XlQgKyBcbWF0aGJme059XlQgLApcXQoKYGBge3J9Ck4gPC0gbWF0cml4KGRpYWcoRyksbnJvdyhYYyksbnJvdyhYYykpCk4KZGlzdDIgPC0gTi0yKkcrdChOKQpkaXN0MgpkaXN0KFhjKV4yCmBgYAoKIyBMaW5rIEdyYW0gTWF0cml4IGFuZCBTcXVhcmVkIERpc3RhbmNlIE1hdHJpeAoKXFtcbWF0aGJme0d9ID1cbWF0aGJme1h9X2NcbWF0aGJme1h9X2NeVD0tXGZyYWN7MX17Mn1cbWF0aGJme0h9XG1hdGhiZntEfV9YXG1hdGhiZntIfVxdCgpgYGB7cn0KRwotMS8yKkglKiVkaXN0MiUqJUgKLTEvMipIJSolZGlzdDIlKiVIIC0gRwpgYGAKCiMgU2hvdyB0aGF0IE4gY2FuY2VscyBvdXQgd2hlbiBtdWx0aXBseWluZyB3aXRoIEgKClxbIFxtYXRoYmZ7SCBOIEh9ID0gXG1hdGhiZnswfVxdCgpcWyBcbWF0aGJme0ggTn1eVCBcbWF0aGJme0h9ID0gXG1hdGhiZnswfVxdCgpgYGB7cn0KTiUqJUgKSCUqJXQoTikKYGBgCgpgYGB7ciwgY2hpbGQ9Il9zZXNzaW9uLWluZm8uUm1kIn0KYGBgCg==