---
title: 'Generalized PLS-SVD: Explicit Whitening Reference'
output:
  rmarkdown::html_vignette:
    toc: yes
    toc_depth: 2.0
    css: albers.css
    includes:
      in_header: albers-header.html
params:
  family: red
  preset: interaction
resource_files:
- albers.css
- albers.js
- albers-header.html
- fonts

vignette: |
  %\VignetteIndexEntry{Generalized PLS-SVD: Explicit Whitening Reference}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
if (requireNamespace("ragg", quietly = TRUE)) knitr::opts_chunk$set(dev = "ragg_png")
if (requireNamespace("systemfonts", quietly = TRUE) && requireNamespace("albersdown", quietly = TRUE)) albersdown::albers_register_fonts()
if (requireNamespace("ggplot2", quietly = TRUE) && requireNamespace("albersdown", quietly = TRUE)) ggplot2::theme_set(albersdown::theme_albers(family = params$family, preset = params$preset))
knitr::opts_chunk$set(
  collapse   = TRUE,
  comment    = "#>",
  message    = FALSE,
  warning    = TRUE,
  fig.width  = 6,
  fig.height = 4,
  out.width  = "85%"
)
library(genpca)
library(Matrix)
```

```{r albers-classes, echo=FALSE, results='asis'}
cat(sprintf(
  paste0(
    '<script>document.addEventListener("DOMContentLoaded",function(){',
    'document.body.classList.remove("palette-red","palette-lapis","palette-ochre","palette-teal","palette-green","palette-violet","preset-homage","preset-interaction","preset-study","preset-structural","preset-adobe","preset-midnight");',
    'document.body.classList.add("palette-%s","preset-%s");',
    '});</script>'
  ),
  params$family,
  params$preset
))
```




Start here to relate two data blocks measured on the same observations.
The first example fits `genpls()` and projects both blocks into latent
coordinates. The later sections provide a dense reference for contributors
checking the whitening identities behind `gplssvd_op()`.

## Quick practical use

Rows must refer to the same observations in the same order. Here `X` has
150 rows and 8 variables, `Y` has 150 rows and 5 variables, and both contain
one shared simulated signal. The example metrics are positive diagonal
weights.

```{r quick-use}
set.seed(123)
N <- 150
shared <- rnorm(N)
X <- outer(shared, seq(0.5, 1.2, length.out = 8)) + matrix(rnorm(N * 8), N, 8)
Y <- outer(shared, seq(0.5, 1.2, length.out = 5)) + matrix(rnorm(N * 5), N, 5)
row_wt   <- diag(runif(N, 0.5, 1.5))
col_wt_x <- diag(runif(8, 0.8, 1.2))
col_wt_y <- diag(runif(5, 0.8, 1.2))

fit <- genpls(X, Y, ncomp = 2,
              preproc_x = multivarious::center(),
              preproc_y = multivarious::center(),
              Mx = row_wt, My = row_wt,
              Ax = col_wt_x, Ay = col_wt_y)
round(fit$d, 3)
```

Project each block using its fitted preprocessing and projection weights:

```{r quick-project}
Sx <- multivarious::project(fit, X)
Sy <- multivarious::project(fit, Y, source = "Y")
cor(Sx[, 1], Sy[, 1])
```

```{r quick-plot, echo = FALSE, fig.cap = "Leading projected coordinates of the two blocks, which share a simulated signal.", fig.height = 3.5}
plot(Sx[, 1], Sy[, 1], pch = 19, col = "steelblue",
     xlab = "X latent coordinate 1", ylab = "Y latent coordinate 1")
```

```{r quick-project-check, include = FALSE}
stopifnot(identical(dim(Sx), c(150L, 2L)), identical(dim(Sy), c(150L, 2L)),
          all(is.finite(Sx)), all(is.finite(Sy)),
          abs(cor(Sx[, 1], Sy[, 1])) > 0.6)
```

The leading coordinates track the shared variation. This is a training-data
association, not a held-out predictive assessment. For new observations,
pass matrices with the same columns to `project()` in the same way. With
nonidentity row metrics these ambient coordinates differ from the stored
row-whitened `lx` and `ly` used in the identities below.

## Reference for contributors

### Notation

GPLSSVD decomposes the relationship between two data blocks `X` (`N x I`) and `Y` (`N x J`) with optional row and column metrics:

- `MX`, `MY`: row metrics (`N x N`) -- weight observations differently for the two blocks
- `WX`, `WY`: column metrics (`I x I` and `J x J`) -- encode within-block variable relationships
- `p`, `q`: generalized singular vectors (saliences) satisfying `p' WX p = I`, `q' WY q = I`
- `Fi`, `Fj`: factor scores (loadings scaled by singular values)
- `Lx`, `Ly`: latent variables (data projections onto components)
- `d`: singular values of the whitened cross-product matrix

### Reference implementation

The function below builds the whitened cross-product `S = (M_X^{1/2} X W_X^{1/2})' (M_Y^{1/2} Y W_Y^{1/2})` explicitly and runs a dense SVD. It is deliberately the most literal transcription of the algebra above, with no attention to speed or to input types beyond what the check below needs; the package's operator path can apply the corresponding products without explicitly assembling these matrices, with dense fallbacks where needed.

Everything rests on one helper. Each metric enters through its symmetric PSD
square root, and singular metrics need the *pseudo*-inverse of that root --
zero eigenvalues stay zero rather than blowing up. This small reference
assumes valid PSD inputs; it is not a replacement for package validation or
numerical rank handling:

```{r ref-sqrt}
psd_sqrt <- function(W, n) {
  if (is.null(W)) return(list(h = diag(n), hinv = diag(n), full = diag(n)))
  W    <- as.matrix(W)
  e    <- eigen(W, symmetric = TRUE)
  lam  <- pmax(e$values, 0)                    # PSD assumed; clip numerical negatives
  half <- function(f) e$vectors %*% (f * t(e$vectors))
  list(h    = half(sqrt(lam)),                 # W^{1/2}
       hinv = half(ifelse(lam > 0, 1 / sqrt(lam), 0)),  # W^{-1/2}, pseudo
       full = W)
}
```

The decomposition itself is then a whitening, one SVD, and an unwhitening:

```{r ref-impl}
dense_gplssvd_ref <- function(X, Y, MX = NULL, MY = NULL,
                              WX = NULL, WY = NULL, k = NULL,
                              center = FALSE, scale = FALSE) {
  X <- scale(as.matrix(X), center = center, scale = scale)
  Y <- scale(as.matrix(Y), center = center, scale = scale)
  stopifnot(nrow(X) == nrow(Y))

  mx <- psd_sqrt(MX, nrow(X)); wx <- psd_sqrt(WX, ncol(X))
  my <- psd_sqrt(MY, nrow(Y)); wy <- psd_sqrt(WY, ncol(Y))

  # whiten both blocks, then SVD their cross-product
  S  <- crossprod(mx$h %*% X %*% wx$h, my$h %*% Y %*% wy$h)
  sv <- svd(S)
  keep <- seq_len(if (is.null(k)) length(sv$d) else min(k, length(sv$d)))

  # unwhiten: saliences are W^{-1/2} u, so that p' WX p = I
  p <- wx$hinv %*% sv$u[, keep, drop = FALSE]
  q <- wy$hinv %*% sv$v[, keep, drop = FALSE]
  D <- diag(sv$d[keep], nrow = length(keep))

  list(d  = sv$d[keep], p = p, q = q,
       fi = wx$full %*% p %*% D,             # factor scores
       fj = wy$full %*% q %*% D,
       lx = mx$h %*% X %*% wx$full %*% p,    # latent variables
       ly = my$h %*% Y %*% wy$full %*% q)
}
```

### Cross-checking the operator

Run the reference on a small block, run `gplssvd_op()` with the same metrics, and compare:

```{r example}
set.seed(1)
N <- 20; I <- 8; J <- 6
X  <- matrix(rnorm(N * I), N, I)
Y  <- matrix(rnorm(N * J), N, J)
MX <- diag(runif(N, .5, 1.5))
MY <- diag(runif(N, .5, 1.5))
WX <- diag(runif(I, .5, 1.5))
WY <- diag(runif(J, .5, 1.5))

ref <- dense_gplssvd_ref(X, Y, MX, MY, WX, WY,
                         k = 3, center = TRUE, scale = FALSE)
op  <- gplssvd_op(X, Y,
                  XLW = MX, YLW = MY,
                  XRW = WX, YRW = WY,
                  k = 3, center = TRUE, scale = FALSE)

all.equal(ref$d, op$d, tolerance = 1e-6)
all.equal(diag(crossprod(op$lx, op$ly)), op$d, tolerance = 1e-6)
round(op$d, 4)
```

```{r ref-vs-op-plot, echo = FALSE, fig.cap = "Reference vs operator singular values agree to plotting precision (left). Latent variables show the expected diagonal cross-product structure (right).", fig.width = 7, fig.height = 3.5}
op_par <- par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))
plot(ref$d, type = "b", pch = 19, col = "grey30",
     xlab = "Component", ylab = "Singular value",
     main = "Singular values")
lines(op$d, type = "b", pch = 21, col = "steelblue")
legend("topright", legend = c("reference", "operator"),
       col = c("grey30", "steelblue"), pch = c(19, 21),
       bty = "n", cex = 0.85)

cp <- as.matrix(crossprod(op$lx, op$ly))
image(t(cp)[, ncol(cp):1], axes = FALSE,
      main = "t(Lx) %*% Ly",
      col = grey.colors(20, start = 0.95, end = 0.2))
par(op_par)
```

The diagonal of `t(Lx) %*% Ly` recovers the singular values, as the GPLSSVD identity guarantees.

## Where next

See [Getting Started](genpca.html) for a getting-started walkthrough and [GPCA Metrics](gpca-metrics.html) for metric recipes that apply to both GPCA and GPLSSVD.