add matrix_X argument

.gitignore

@@ -49,3 +49,7 @@ po/*~
 # RStudio Connect folder
 rsconnect/
 
+
+/.quarto/
+**/*.quarto_ipynb
+_freeze/

R/singular_value_plot.R

@@ -3,6 +3,7 @@ source(here::here("R", "singular_values.R"))
 source(here::here("R", "graphon_distribution.R"))
 
 
+# expr_to_label ----------------------------------------------------------------
 # Convert a call or character to a nicely formatted character string.
 #   * If the user supplied a character, we keep it unchanged.
 #   * If the user supplied a call (e.g. quote(20 / sqrt(x))) we deparse it
@@ -17,6 +18,7 @@ expr_to_label <- function(expr) {
 }
 
 
+# smallest_sv_sequence ---------------------------------------------------------
 #' Compute the smallest singular value of a sequence of matrices Q(K)
 #'
 #' @title Smallest singular values for a family of matrices Q(K)
@@ -149,10 +151,10 @@ smallest_sv_sequence <- function(
       sample_X_fn = sampler_fn,
       fv = fv,
       Fv = Fv,
-      guard = guard
+      guard = guard,
+      scaled = FALSE
     )
     
-    Q <- 1 /sqrt(n) * Q
     
     sv_res <- compute_minmax_sv(Q)
     if (!is.list(sv_res) || is.null(sv_res$smallest_singular_value)) {

R/singular_values.R

@@ -44,6 +44,8 @@ source(here::here("R", "graphon_distribution.R"))
 #' @param Fv     Cumulative distribution function of the latent variable
 #'               \eqn{v}.  Also has to be vectorised.  Typical examples are
 #'               `pnorm`, `pexp`, ….
+#' @param matrix_X matrix with the covariates at each node. Each row corresponds
+#' to a single node with p attributes.
 #' @param guard Positive numeric guard value.  Default is `sqrt(.Machine$double.eps)`,
 #'   which is about `1.5e‑8` on most platforms – small enough to be negligible
 #'   for most computations. If it is null, then it is not used. 
@@ -107,6 +109,7 @@ compute_matrix <- function(
     sample_X_fn,
     fv,
     Fv,
+    matrix_X = NULL,
     guard = sqrt(.Machine$double.eps),
     scaled = FALSE
 ) {
@@ -118,14 +121,21 @@ compute_matrix <- function(
   if (!is.function(sample_X_fn)) stop("'sample_X_fn' must be a function")
   if (!is.function(fv)) stop("'f_v' must be a function")
   if (!is.function(Fv)) stop("'F_v' must be a function")
+  if (!is.null(matrix_X) && !is.matrix(matrix_X)) stop("matrix_X must be either null or a matrix")
   
   ## 1.2 Generate the Matrix X of covariates ===================================
-  # The withr environment allows us to capsulate the global state like the seed
-  # and enables a better reproduction
-  X <- withr::with_seed(seed, {
-    as.matrix(sample_X_fn(n))
-  })
-  if (nrow(X) != n) stop("`sample_X_fn` must return exactly `n` rows")
+  # If the argument matrix_X is present, use this matrix, otherwise generate one
+  # with sample_X_fn.
+  if (!is.null(matrix_X)) {
+    X <- matrix_X
+  } else {
+    # The withr environment allows us to encapsulate the global state like the seed
+    # and enables a better reproduction
+    X <- withr::with_seed(seed, {
+      as.matrix(sample_X_fn(n))
+    })
+  }
+  if (nrow(X) != n) stop(" the covariate matrix `X` must have exactly `n` rows")
   if (ncol(X) != length(a)) {
     stop("Number of columns of X (", ncol(X), ") must equal length(a) (", length(a), ")")
   }

_quarto.yml

@@ -0,0 +1,6 @@
+project:
+  type: default
+
+execute:
+  freeze: auto
+  cache: false

scripts/.gitignore

@@ -0,0 +1,2 @@
+*.html
+*plots_dimensions_files/

scripts/plots_dimensions.qmd

@@ -0,0 +1,122 @@
+---
+title: "Plots of n vs. k"
+author: "Niclas"
+format: html
+editor: visual
+execute:
+  echo: true
+  working-directory: ../
+---
+
+# Plots of the dimensions
+
+## Setup 
+We consider the matrix $QQ^\top$ and look at the smallest eigenvalue, i.e. the
+smallest non-zero singular value of $Q$.
+
+The matrix $Q$ is given by
+$$
+Q_{ik} = \int_{\frac{k}{K}}^{\frac{k+1}{K}} p_a(u| X_i) \, du
+$$
+with
+$$
+p_a(u|X) = \frac{f_v(F_a^{-1}(u) - a^\top X)}{f_a(F_a^{-1}(u))}
+$$
+
+## Plots of n vs. k
+
+- The $v$'s are normally distributed with $v \sim \mathcal N(0,1)$
+- Plot $n = 100, 200, 300, 400$ and $k = 1, \dots, K$ with $K = \sqrt n$.
+
+```{r Load Libraries}
+# load local files
+source(here::here("R", "singular_values.R"))
+source(here::here("R", "graphon_distribution.R"))
+source(here::here("R","singular_value_plot.R"))
+
+# load libaries for data handling
+library(ggplot2)
+library(dplyr)
+library(latex2exp)
+```
+
+```{r Compute the data}
+#| cache: true
+#| echo: false
+#| collapse: true
+ns <- c(100, 200, 300, 400, 500)
+Ks <- floor(sqrt(ns))
+as <- c(0.5, 1.0, 1.5, 2.0)
+
+# set a global seed
+set.seed(42)
+results <- data.frame(dim_n = integer(),
+                      dim_k = integer(),
+                      param_a = double(),
+                      ssv = double())
+for (a in as) {
+  for (i in 1:length(ns)) {
+    n <- ns[i]
+    K <- Ks[i]
+    # use the default seed 1L
+    out <- smallest_sv_sequence(
+      a = a,
+      n = n,
+      maxK = K,
+      sampler_fn =function(n) matrix(rnorm(n), ncol = 1L),
+      guard=1e-12,
+      plot=FALSE,
+      fv = function(x) {dnorm(x, mean=0, sd=1)},
+      Fv = function(x) {pnorm(x, mean=0, sd=1)}
+      )
+    
+    current_res <- data.frame(dim_n = rep(n, K), dim_k = out$K, param_a = rep(a, K), ssv = out$sv)
+    results <- rbind(results, current_res)
+  }
+}
+
+```
+
+```{r plot the results}
+#| cache: true
+#| echo: false
+#| collapse: true
+#| fig-cap: "Simulation of the smallest singular values w.r.t. a, n and k"
+results |>
+  mutate(param_a =  as.factor(param_a),
+         dim_n = as.factor(dim_n)) |>
+  group_by(param_a, dim_n) |>
+  ggplot(aes(dim_k, ssv, col=dim_n, shape=param_a, interaction(dim_n, param_a))) +
+  geom_point(size=1.5) +
+  geom_line() +
+  scale_y_log10() +
+  theme_bw() +
+  labs(x=latex2exp::TeX("$k$"),
+       y=latex2exp::TeX("Smallest singular value of $Q$"),
+       title=latex2exp::TeX("Smallest singular value of $Q$ with respect to $n$, $k$, and $a$."),
+       colour=latex2exp::TeX("$n$"),
+       shape=latex2exp::TeX("$a$"))
+```
+The data for $n = 100$ is covered by the data for $n = 200$.
+
+## Analysis of the convergence
+
+We assume that the smallest singular value $\sigma$ can be approximated by:
+$$
+\sigma = C \cdot n^\eta \cdot k^\kappa \cdot a^\alpha
+$$
+to estimate the coefficients we make a log-transform and perform a linear regression, i.e.
+$$
+\log(\sigma) = \log (C) + \eta \log(n) + \kappa\log(k) + \alpha \log(a).
+$$
+
+```{r estimate coeffs}
+model1 <- results |>
+  filter(ssv > 1e-15) |> # exclude to small values
+  lm(formula = log(ssv) ~ log(dim_n) + log(dim_k) + log(param_a))
+
+summary(model1)
+plot(model1)
+  
+```
+