experiments with rho_n

R/build_network.R

@@ -135,10 +135,3 @@ compute_adj_matrix <- function(
   adj_mat
 }
 
-set.seed(1)
-X <- matrix(c(-1, -0.5, 0, 0.5, 1), nrow = 5, ncol=1)
-a <- 2.0
-v <- c(0, 0.2, 0.4, 0.6, 0.8)
-adj <- compute_adj_matrix(X, v, a, phi = function(x,y) {x * y}, 0.5, Fv = punif)
-adj
-

R/estimators.R

@@ -0,0 +1,246 @@
+# 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"))
+source(here::here("R", "build_network.R"))
+
+# Helper functions -------------------------------------------------------------
+# helper function for wrapping the parameters of the Q_a creation function
+# TODO rename this function
+make_matrix_creation <- function(seed, n, K, matrix_X, fv, Fv, guard) {
+  function(a) {
+    compute_matrix(seed=seed, a, n=n, K=K, matrix_X = matrix_X, fv=fv, Fv=Fv, guard=guard)
+  }
+}
+
+# estimators -------------------------------------------------------------------
+
+## Moore-Penrose Inverse -------------------------------------------------------
+#' Moore‑Penrose pseudoinverse of a matrix
+#'
+#' Computes the Moore‑Penrose generalized inverse of a numeric matrix using a
+#' singular‑value decomposition (SVD). Singular values smaller than the
+#' tolerance are treated as zero to improve numerical stability.
+#'
+#' @param A A numeric matrix (or an object coercible to a matrix) whose
+#'   pseudoinverse is required.
+#' @param tol Tolerance for treating singular values as zero.  By default it
+#'   is set to `max(dim(A)) * max(svd(A)$d) * .Machine$double.eps`, which
+#'   scales with the size of the matrix and machine precision.
+#' @return A matrix representing the Moore‑Penrose inverse of `A`.  The
+#'   dimensions of the result are `ncol(A) × nrow(A)`.
+#' @examples
+#' set.seed(123)
+#' A <- matrix(rnorm(12), nrow = 3)
+#' A_pinv <- pinv(A)
+#' # Verify the Moore‑Penrose properties (A %*% A_pinv %*% A ≈ A)
+#' all.equal(A %*% A_pinv %*% A, A)
+#' @importFrom stats svd
+#' @export
+pinv <- function(A, tol = NULL) {
+  # Coerce to matrix and check type
+  if (!is.matrix(A)) {
+    A <- as.matrix(A)
+  }
+  if (!is.numeric(A)) {
+    stop("`A` must be a numeric matrix.", call. = FALSE)
+  }
+  
+  # Singular value decomposition
+  s <- svd(A)
+  
+  # Determine tolerance if not supplied
+  if (is.null(tol)) {
+    tol <- max(dim(A)) * max(s$d) * .Machine$double.eps
+  }
+  
+  # Invert non‑zero singular values
+  d_inv <- ifelse(s$d > tol, 1 / s$d, 0)
+  
+  # Construct diagonal matrix of inverted singular values
+  D_plus <- diag(d_inv, nrow = length(d_inv))
+  
+  # Moore‑Penrose inverse: V %*% D⁺ %*% t(U)
+  s$v %*% D_plus %*% t(s$u)
+}
+
+
+
+## Estimate Matrix B -----------------------------------------------------------
+#' Estimate the matrix \$B\$ for a graphon model
+#'
+#' For a given graphon scaling parameter `rho_n`, a square matrix `Q_a`, and an
+#' adjacency matrix `A`, this function computes  
+
+#' $ B = \\rho_n \, Q_a^{+\\,T} \, A \, Q_a^{+} $
+
+#' where `Q_a^{+}` denotes the Moore‑Penrose pseudoinverse of `Q_a`. 
+#'
+#' @param rho_n Numeric scalar.  The graphon scaling parameter.
+#' @param Q_a   Square numeric matrix. TODO: write description of the Matrix
+#' @param A     Square numeric adjacency matrix (same dimension as `Q_a`).
+#' @return A numeric matrix of the same dimension as `Q_a` representing the
+#'   estimated \$B\$.
+#' @examples
+#' set.seed(42)
+#' n <- 5
+#' Q_a <- diag(n)               # simple identity basis
+#' A   <- matrix(rbinom(n^2, 1, 0.3), n, n)
+#' rho_n <- 0.5
+#' B_est <- estimate_B_matrix(rho_n, Q_a, A)
+#' str(B_est)
+#' @importFrom stats svd
+#' @export
+estimate_B_matrix <- function(rho_n, Q_a, A) {
+  # ---- Input checks ---------------------------------------------------------
+  if (!is.numeric(rho_n) || length(rho_n) != 1) {
+    stop("`rho_n` must be a single numeric value.", call. = FALSE)
+  }
+  if (!is.matrix(Q_a) || !is.numeric(Q_a)) {
+    stop("`Q_a` must be a numeric matrix.", call. = FALSE)
+  }
+  if (!is.matrix(A) || !is.numeric(A)) {
+    stop("`A` must be a numeric matrix.", call. = FALSE)
+  }
+  if (nrow(A) != ncol(A)) {
+    stop("`A` must be square.", call. = FALSE)
+  }
+  if (ncol(Q_a) != nrow(A)) {
+    stop("Dimensions of `Q_a` and `A` must agree for matrix multiplication.", call. = FALSE)
+  }
+  
+  # ---- Compute pseudoinverse ------------------------------------------------
+  pinv_Qa <- pinv(Q_a)   # assumes `pinv()` is available in the namespace
+  
+  # ---- Estimate B -----------------------------------------------------------
+  B <- rho_n * (t(pinv_Qa) %*% A %*% pinv_Qa)
+  
+  B
+}
+
+# TODO rename this function
+# TODO test the convergence of the function estimate_B
+# with given graphon (block function, Hölder continuous functions)
+# and given K with growing N
+# and other options
+# plot the loss function with respect to a
+estimate_a <- function(A, # adjacency matrix
+                       a0, # start value
+                       n,
+                       K,
+                       sample_X_fn,
+                       fv,
+                       Fv, 
+                       guard
+) {
+  calc_Q_a <- make_matrix_creation(seed, n, K, sample_X_fn, fv, Fv, guard)
+  
+  loss_func <- function(a) {
+    Q_a <- calc_Q_a(a)
+    pinv_Qa <- pinv(Q_a)
+    
+    norm(pinv_Qa %*% Q_a %*% A %*% pinv_Qa %*% Q_a - A, type="F")^2
+  }
+  
+  return(loss_func)
+}
+
+#' Calculate the edge density of an undirected graph
+#'
+#' @title Edge density from an adjacency matrix
+#' @description Computes the proportion of possible edges that are present in an
+#'   undirected, unweighted graph represented by a square adjacency matrix.  The
+#'   density is defined as \eqn{2E / (V(V-1))} where \eqn{E} is the number of edges
+#'   and \eqn{V} is the number of vertices.  This corresponds to the parameter
+#'   \eqn{rho_n}.
+#'   The function checks that the input is a square matrix and treats any
+#'   non‑numeric or missing entries as absent edges.
+#'
+#' @param adj_matrix A square numeric matrix representing the adjacency matrix
+#'   of an undirected graph.  Entries should be 0/1 (or any truthy numeric) with
+#'   `adj_matrix[i, j] == adj_matrix[j, i]`.  Missing values are ignored.
+#'
+#' @return A single numeric value giving the edge density \eqn{rho}.
+#'
+#' @examples
+#' # Simple triangle graph (3 nodes, 3 edges)
+#' A_tri <- matrix(c(0,1,1,
+#'                   1,0,1,
+#'                   1,1,0), nrow = 3, byrow = TRUE)
+#' calculate_edge_density(A_tri)
+#'
+#' # Empty graph (no edges)
+#' A_empty <- matrix(0, nrow = 4, ncol = 4)
+#' calculate_edge_density(A_empty)
+#'
+#' @export
+calculate_edge_density <- function(adj_matrix) {
+  # -------------------------------------------------------------------------
+  # Validate input
+  # -------------------------------------------------------------------------
+  if (!is.matrix(adj_matrix)) {
+    stop("`adj_matrix` must be a matrix.")
+  }
+  if (nrow(adj_matrix) != ncol(adj_matrix)) {
+    stop("`adj_matrix` must be square (same number of rows and columns).")
+  }
+  
+  # -------------------------------------------------------------------------
+  # Count edges and nodes
+  # -------------------------------------------------------------------------
+  edge_idx   <- which(upper.tri(adj_matrix, diag = FALSE), arr.ind = TRUE)
+  edge_count <- sum(adj_matrix[edge_idx], na.rm = TRUE)
+  node_count <- nrow(adj_matrix)
+  
+  # -------------------------------------------------------------------------
+  # Compute density
+  # -------------------------------------------------------------------------
+  rho <- (2 * edge_count) / (node_count * (node_count-1))
+  
+  return(rho)
+}
+# test the estimator routines
+seed <- 17L #  121L this seed works exceptionally well
+set.seed(seed)
+#X  <- matrix(seq(-1, 1, length.out = 5), ncol = 1)
+a  <- 20
+n <- 400
+K <- 4
+rho_n <- log(n) / n
+sample_X_fn <- function(n) {matrix(rnorm(n), ncol = 1L)}
+X <- sample_X_fn(n)
+fv <- function(x) {dnorm(x, mean=0, sd=1)}
+Fv <- function(x) {pnorm(x, mean=0, sd=1)}
+guard <- 1e-12
+v  <- rnorm(n) #seq(0, 0.8, length.out = n)
+phi_fun <-Vectorize(function(x, y) ifelse(((x > 0.5 && y <= 0.5) || (x <= 0.5 && y > 0.5)), 1.6, 0.4))        # multiplicative kernel
+adj <- compute_adj_matrix(
+   X_matrix = X,
+   v        = v,
+   a        = a,
+   phi      = phi_fun,
+   rho_n    = rho_n,
+   Fv       = Fv
+ )
+adj
+
+# Q_a matrix
+Qa <- compute_matrix(seed, a=a, n=n, K=K, fv=fv, Fv=Fv, guard=guard, matrix_X=X)
+
+calc_Q_a <- make_matrix_creation(seed, n=n, K=K, matrix_X = X, fv=fv, Fv=Fv, guard=guard)
+
+loss_func <- function(a) {
+  Q_a <- calc_Q_a(a)
+  pinv_Qa <- pinv(Q_a)
+  
+  norm(pinv_Qa %*% Q_a %*% adj %*% pinv_Qa %*% Q_a - adj, type="F")^2
+}
+
+
+plot_as <- seq(-10, 100, length.out=500)
+loss_vals <- sapply(plot_as, loss_func)
+plot(plot_as , loss_vals, type="b")
+abline(v=a)
+title("Test plot of the loss function")
+
+optim(25, loss_func, lower=10, upper=100, method="Brent", control=list(fnscale=1))

R/singular_values.R

@@ -33,10 +33,6 @@ source(here::here("R", "graphon_distribution.R"))
 #'               generated.
 #' @param K      Positive integer.  Number of divisions of the unit interval;
 #'               the resulting grid has length `K+1`.
-#' @param sample_X_fn
-#'               Function with a single argument `n`.  It must return an
-#'               \eqn{n \times p} matrix (or an object coercible to a matrix) of
-#'               covariate samples.
 #' @param fv     Density function of the latent variable \eqn{v}.  Must be
 #'               vectorised (i.e. accept a numeric vector and return a numeric
 #'               vector of the same length).  Typical examples are
@@ -44,8 +40,15 @@ 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 sample_X_fn
+#'               Function with a single argument `n`.  It must return an
+#'               \eqn{n \times p} matrix (or an object coercible to a matrix) of
+#'               covariate samples. It can be NULL, but then `matrix_X` must be
+#'               given.
 #' @param matrix_X matrix with the covariates at each node. Each row corresponds
-#' to a single node with p attributes.
+#' to a single node with p attributes. The default value is `NULL`. If it
+#' is `NULL` then `sample_X_fun` must be given. If both parameters are provided,
+#' then `matrix_X` is used.
 #' @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. 
@@ -106,9 +109,9 @@ compute_matrix <- function(
     a,
     n,
     K,
-    sample_X_fn,
     fv,
     Fv,
+    sample_X_fn=NULL,
     matrix_X = NULL,
     guard = sqrt(.Machine$double.eps),
     scaled = FALSE
@@ -118,10 +121,12 @@ compute_matrix <- function(
   if (!is.numeric(a) || !is.vector(a)) stop("'a' must be a numeric vector")
   if (!is.numeric(n) || length(n) != 1 || n <= 0) stop("'n' must be a positive integer")
   if (!is.numeric(K) || length(K) != 1 || K <= 0) stop("'K' must be a positive integer")
-  if (!is.function(sample_X_fn)) stop("'sample_X_fn' must be a function")
+  if (!(is.function(sample_X_fn) || is.null(sample_X_fn))) stop("'sample_X_fn' must be a function or Null and matrix_X must be given")
   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")
+  if (is.null(matrix_X) && is.null(sample_X_fn)) stop("Either 'matrix_X' or 'sample_X_fn' must be supplied!")
+  if (!is.null(matrix_X) && !is.null(sample_X_fn)) warning("Both arguments 'matrix_X' and `sample_X_fn` is given. Priority is given by to the first!")
   
   ## 1.2 Generate the Matrix X of covariates ===================================
   # If the argument matrix_X is present, use this matrix, otherwise generate one

notebooks/two-dimensional-calculation-example.qmd

@@ -0,0 +1,48 @@
+---
+title: "Two-dimensional-example"
+format: html
+editor: visual
+---
+
+# Zwei dimensionales Rechenbeispiel für die Matrix $Q_a$.
+
+Mit den Eingabedaten:
+
+- $a = 2.0$
+- $X = (-0.6264, 0.1836)^\top$
+- $n = K = 2$
+und $v \sim \mathcal{N}(0, 1)$
+
+Daraus ergibt sich eine Matrix $Q_a$ mit
+
+$$
+Q_a = \begin{pmatrix} 0.7911 & 0.2089 \\
+ 0.2089 & 0.7911 \\
+ \end{pmatrix}
+$$
+
+
+Anscheinend ist es so, dass für verschiedene Eingabewerte der Matrix $X$, es immer
+wieder zu verschiedenen, aber immer noch diagonaldominanten Einträgen kommt. 
+Wieso passiert dies?
+
+Falls allerdings $X = (0.2167549 -0.5424926)^\top$, dann ist 
+$$
+Q_a = \begin{pmatrix}
+0.2239 & 0.7761 \\
+0.7761 & 0.2239 \\
+\end{pmatrix}
+$$
+wo jetzt die Nebendiagonale dominant ist. Wenn man verschiedene Seeds ausprobiert,
+so sieht man ein Muster. Doch es gibt auch das Beispiel mit $X = (0.7667960 -0.8164583)^\top$
+und dann erhalten wir:
+$$
+Q_a = \begin{pmatrix}
+0.4802 & 0.5198 \\
+0.5198 & 0.4802 \\
+\end{pmatrix}
+$$
+
+Diese Matrix hat immer noch eine dominante Nebendiagonale, aber nicht so stark wie
+alle anderen bekannten Beispiele. Bei den anderen Berechnungen zeigte sich meistens
+ein Unterschied von $| a_{11} - a_{12}| > 0.5$.

scripts/plot_Qa_wrt_a.R

@@ -13,7 +13,7 @@ library(dplyr)
 a_grid <- seq(-20, 20, length.out = 200)
 
 # function which takes only a to compute Q_c
-make_matrix <- function(a) { compute_matrix(seed=4L,
+make_matrix <- function(a) { compute_matrix(seed=11513215L,
                               a= a,
                               n = 2,
                               K = 2,
@@ -46,17 +46,18 @@ ggplot(df_entries, aes(x = a, y = value, colour = entry, linetype = entry)) +
   theme_minimal()
 
 # Heat map for a single larger matrix ------------------------------------------
-
+# TODO Daten für 2x2 und 3x3 an Michael schicken
 # Choose a value of a
-a0 <- -10
-M0 <-  compute_matrix(seed=1L,
+a0 <- 2
+M0 <-  compute_matrix(seed=9L,
                       a= a0,
-                      n = 50,
-                      K = 50,
+                      n = 2,
+                      K = 2,
                       sample_X_fn = function(n) {matrix(rnorm(n), ncol = 1L)},
                       fv = function(x) {dnorm(x, mean=0, sd=1)},
                       Fv = function(x) {pnorm(x, mean=0, sd=1)},
                       guard = 1e-12)
+M0
 
 # Convert to a tidy data frame for ggplot
 df_heat <- as.data.frame(M0) %>%

scripts/plot_non_normal_X.qmd

@@ -76,14 +76,16 @@ for (a in as) {
 ```{r k = n^alpha plotting, rate = 1}
 # plot the results
 results01 |>
-  filter(param_a %in% c(0, 10, 20)) |>
+  filter(param_a %in% c(2, 6, 12) & param_alpha <= 0.12) |>
   mutate(param_a = as.factor(param_a),
          param_alpha = as.factor(param_alpha)) |>
   group_by(param_a, param_alpha) |>
   ggplot(aes(dim_n, ssv, col=param_a, shape=param_alpha, interaction(param_a, param_alpha))) +
   geom_point(size=1.5) +
   geom_line() +
-  geom_function(fun = function(x) {sqrt(x)}, colour="black") +
+  geom_line(aes(dim_n, sqrt(dim_n) / dim_k, shape=param_alpha), linetype = 2) +
+  geom_point(size=1.5) +
+  #geom_function(fun = function(x) {sqrt(x)}, colour="black") +
   #scale_y_log10() +
   theme_bw() +
   labs(x=latex2exp::TeX("$n$"),
@@ -135,14 +137,16 @@ for (a in as) {
 
 ```{r k = n^alpha plotting, rate = 3}
 results02 |>
-  filter(param_a %in% c(0, 10, 20)) |>
+  filter(param_a %in% c(0, 10, 20) & param_alpha < 0.12) |>
   mutate(param_a = as.factor(param_a),
          param_alpha = as.factor(param_alpha)) |>
   group_by(param_a, param_alpha) |>
   ggplot(aes(dim_n, ssv, col=param_a, shape=param_alpha, interaction(param_a, param_alpha))) +
   geom_point(size=1.5) +
   geom_line() +
-  geom_function(fun = function(x) {sqrt(x)}, colour="black") +
+  geom_line(aes(dim_n, sqrt(dim_n) / dim_k, shape=param_alpha), linetype = 2) +
+  geom_point(size=1.5) +
+  #geom_function(fun = function(x) {sqrt(x)}, colour="black") +
   #scale_y_log10() +
   theme_bw() +
   labs(x=latex2exp::TeX("$n$"),
@@ -203,7 +207,9 @@ results03 |>
   ggplot(aes(dim_n, ssv, col=param_a, shape=param_alpha, interaction(param_a, param_alpha))) +
   geom_point(size=1.5) +
   geom_line() +
-  geom_function(fun = function(x) {sqrt(x)}, colour="black") +
+  geom_line(aes(dim_n, sqrt(dim_n) / dim_k, shape=param_alpha), linetype = 2) +
+  geom_point(size=1.5) +
+  #geom_function(fun = function(x) {sqrt(x)}, colour="black") +
   #scale_y_log10() +
   theme_bw() +
   labs(x=latex2exp::TeX("$n$"),
@@ -259,14 +265,16 @@ for (a in as) {
 
 ```{r k = n^alpha plotting, U[0,1]}
 results04 |>
-  filter(param_a %in% c(0, 10, 20)) |>
+  filter(param_a %in% c(2, 10, 20) & param_alpha < 0.22 & param_alpha > 0.12) |>
   mutate(param_a = as.factor(param_a),
          param_alpha = as.factor(param_alpha)) |>
   group_by(param_a, param_alpha) |>
   ggplot(aes(dim_n, ssv, col=param_a, shape=param_alpha, interaction(param_a, param_alpha))) +
   geom_point(size=1.5) +
   geom_line() +
-  geom_function(fun = function(x) {sqrt(x)}, colour="black") +
+  geom_line(aes(dim_n, sqrt(dim_n) / dim_k, shape=param_alpha), linetype = 2) +
+  geom_point(size=1.5) +
+  #geom_function(fun = function(x) {sqrt(x)}, colour="black") +
   #scale_y_log10() +
   theme_bw() +
   labs(x=latex2exp::TeX("$n$"),
@@ -327,7 +335,9 @@ results05 |>
   ggplot(aes(dim_n, ssv, col=param_a, shape=param_alpha, interaction(param_a, param_alpha))) +
   geom_point(size=1.5) +
   geom_line() +
-  geom_function(fun = function(x) {sqrt(x)}, colour="black") +
+  geom_line(aes(dim_n, sqrt(dim_n) / dim_k, shape=param_alpha), linetype = 2) +
+  geom_point(size=1.5) +
+  #geom_function(fun = function(x) {sqrt(x)}, colour="black") +
   #scale_y_log10() +
   theme_bw() +
   labs(x=latex2exp::TeX("$n$"),
@@ -387,7 +397,9 @@ results06 |>
   ggplot(aes(dim_n, ssv * dim_k, col=param_a, shape=param_alpha, interaction(param_a, param_alpha))) +
   geom_point(size=1.5) +
   geom_line() +
-  geom_function(fun = function(x) {x^(0.5)}, colour="black") +
+  geom_line(aes(dim_n, sqrt(dim_n) / dim_k, shape=param_alpha), linetype = 2) +
+  geom_point(size=1.5) +
+  #geom_function(fun = function(x) {x^(0.5)}, colour="black") +
   #scale_y_log10() +
   theme_bw() +
   labs(x=latex2exp::TeX("$n$"),
@@ -422,3 +434,68 @@ results06 |>
 results <- list(results01, results02, results03, results04, results05, results06)
 save(results, file="results.RData")
 ```
+
+## Two dimensional example
+
+```{r k = n^alpha data generation with two dimensions, rate = 1}
+#| cache: true
+#| echo: false
+#| collapse: true
+ns <- seq(100, 1000, 100)
+a <- c(1, 1) / sqrt(2)
+a_norms <- seq(0, 20, 2)
+alphas <-  seq(0.1, 0.5, 0.1)
+
+set.seed(100)
+results07 <- data.frame(dim_n = integer(),
+                      dim_k = integer(),
+                      param_a = double(),
+                      param_alpha = double(),
+                      ssv = double())
+for (a_norm in a_norms) {
+  for (i in 1:length(ns)) {
+    for (j in 1:length(alphas)) {
+      n <- ns[i]
+      K <- floor(n^alphas[j])
+      if (!K > 0) next # skip if K is equal to zero
+      # use the default seed 1L
+      Q <- compute_matrix(seed=1L,
+                     a= a_norm * a,
+                     n = n,
+                     K = K,
+                     sample_X_fn = function(n) {matrix(rexp(2 * n), ncol = 2L)},
+                     fv = function(x) {dnorm(x, mean=0, sd=1)},
+                     Fv = function(x) {pnorm(x, mean=0, sd=1)},
+                     guard = 1e-12)
+      
+      ssv <- compute_minmax_sv(Q)[["smallest_singular_value"]]
+      
+      current_res <- data.frame(dim_n = n, dim_k = K, param_a = a_norm, param_alpha=alphas[j], ssv =ssv)
+      results07 <- rbind(results07, current_res)
+    }
+  }
+}
+```
+
+```{r k = n^alpha plotting, rate = 1}
+# plot the results
+results07 |>
+  filter(param_a %in% c(10, 20) & param_alpha < 0.12) |>
+  mutate(param_a = as.factor(param_a),
+         param_alpha = as.factor(param_alpha)) |>
+  group_by(param_a, param_alpha) |>
+  ggplot(aes(dim_n, ssv, col=param_a, shape=param_alpha, interaction(param_a, param_alpha))) +
+  geom_point(size=1.5) +
+  geom_line() +
+  geom_line(aes(dim_n, sqrt(dim_n) / dim_k, shape=param_alpha), linetype = 2) +
+  geom_point(size=1.5) +
+  #geom_function(fun = function(x) {sqrt(x)}, colour="black") +
+  #scale_y_log10() +
+  theme_bw() +
+  labs(x=latex2exp::TeX("$n$"),
+       y=latex2exp::TeX("Smallest singular value of $Q$"),
+       title=latex2exp::TeX("Smallest singular value of $Q$ with respect to $a$."),
+       subtitle = latex2exp::TeX(("Hyperparameter $k = n^{\\alpha}$. Black line is $\\sqrt{n}$, and $X \\sim Exp(1)$")),
+       colour=latex2exp::TeX("$a$"),
+       shape=latex2exp::TeX("$\\alpha$"))
+```