tests/testthat/test_GramMatrix.R

In DistributionalDDA: Distributional Data Analysis in Bayes Spaces

library(testthat)
library(DistributionalDDA)

## Regression tests for the Gram matrix (metric tensor) S of ZB-spline bases.
##
## Background: until commit e3a7015, BayesSpaceBasis.basisfd stored S = I
## (orthonormality assumption; Issue 1 in TODO_orthonormality_issue.md).
## ZB-splines are NOT orthonormal, so S must be the true Gram matrix
## S = A %*% inprod(B, B) %*% t(A). These tests pin down:
##   (a) the closed form of S,
##   (b) its structural properties (symmetric, positive definite, != I),
##   (c) the CONVENTION: consumers evaluate <c,d>_S = c^T S d, i.e. S is the
##       Gram matrix itself, NOT its inverse. The inverse Gram is applied only in
##       the L2 projection (.SplineCoefsFromLogDensity, solve(S, .) -- formerly
##       Issue 3, now fixed).

## trapezoidal rule on an equispaced grid
.trapz <- function(y, dx) dx * (sum(y) - (y[1] + y[length(y)]) / 2)

test_that("BayesSpaceBasis.basisfd stores S = A %*% inprod(B,B) %*% t(A)", {
  bsb <- fda::create.bspline.basis(rangeval = c(0, 1), nbasis = 8, norder = 4)
  zb  <- BayesSpaceBasis(bsb)                    # dispatches to .basisfd
  A   <- zb$A
  expect_equal(unname(getMetricTensor(zb)),
               unname(A %*% fda::inprod(bsb, bsb) %*% t(A)),
               tolerance = 1e-10)
})

test_that("metric tensor of default ZB-spline bases is symmetric, PD, and clearly not I", {
  for (basis in list(B2Unif01Basis5, B2Unif01Basis7)) {
    S <- getMetricTensor(basis)
    expect_equal(unname(S), unname(t(S)), tolerance = 1e-10)
    expect_true(all(eigen(S, symmetric = TRUE, only.values = TRUE)$values > 0))
    ## ZB-splines are not orthonormal: pre-fix this was exactly 0
    expect_gt(sqrt(sum((S - diag(nrow(S)))^2)), 1)
  }
})

test_that("S agrees with grid quadrature of the evaluated ZB-spline products", {
  basis <- B2Unif01Basis5
  tt <- seq(0, 1, length.out = 20001)
  dx <- tt[2] - tt[1]
  Z  <- basisEval(basis, tt)                     # ZB-splines on the grid
  S  <- getMetricTensor(basis)
  p  <- ncol(Z)
  S_grid <- matrix(0, p, p)
  for (i in seq_len(p))
    for (j in seq_len(p))
      S_grid[i, j] <- .trapz(Z[, i] * Z[, j], dx)
  expect_equal(unname(S), unname(S_grid), tolerance = 1e-3)
})

test_that("scalarProduct on BDataSplines reproduces the L2 inner product (Gram, not inverse)", {
  basis <- B2Unif01Basis7
  p  <- getBayesSpaceBasisDimension(basis)
  c1 <- sin(seq_len(p))                          # deterministic, mixed-sign coefficients
  c2 <- cos(seq_len(p))
  x  <- BDataSplines(t(c1), basis)
  y  <- BDataSplines(t(c2), basis)

  tt <- seq(0, 1, length.out = 20001)
  dx <- tt[2] - tt[1]
  Z  <- basisEval(basis, tt)
  ip_quad <- .trapz(as.vector(Z %*% c1) * as.vector(Z %*% c2), dx)

  expect_equal(as.numeric(scalarProduct(x, y)), ip_quad, tolerance = 1e-3)

  ## guard against accidentally storing the inverse Gram matrix:
  ## c^T S^{-1} d must NOT reproduce the inner product
  S <- getMetricTensor(basis)
  ip_inverse <- as.numeric(t(c1) %*% solve(S) %*% c2)
  expect_false(isTRUE(all.equal(ip_inverse, ip_quad, tolerance = 0.05)))
})

test_that("princomp detects the non-orthonormal default basis and stores its Gram matrix", {
  set.seed(2)
  X <- BDataSplines(matrix(rnorm(10 * 7), ncol = 7), B2Unif01Basis7)
  pca <- princomp(X)
  expect_false(pca$orthonormal)                  # pre-fix: TRUE (correction never activated)
  expect_equal(unname(pca$GramMatrix),
               unname(getMetricTensor(B2Unif01Basis7)),
               tolerance = 1e-10)
})