diff --git a/R/singular_values.R b/R/singular_values.R
index 41e1f71..8bf3ec3 100644
--- a/R/singular_values.R
+++ b/R/singular_values.R
@@ -44,6 +44,10 @@ source("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 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. 
+#'   The guard is used for the value k = 0, which can cause arithmetic errors.
 #'
 #' @return A numeric matrix **Q** of dimension `K × n`.  The \eqn{j}-th row
 #'         (for `j = 1,…,K`) contains the increments of the CDF evaluated at
@@ -100,7 +104,8 @@ compute_matrix <- function(
     K,
     sample_X_fn,
     fv,
-    Fv
+    Fv,
+    guard = sqrt(.Machine$double.eps)
 ) {
   ## 1.1 Check inputs ==========================================================
   if (!is.numeric(seed) || length(seed) != 1) stop("'seed' must be a single number")
@@ -127,6 +132,9 @@ compute_matrix <- function(
   
   ## 1.4 Compute the graphon quantiles =========================================
   k <- seq(0, K) / K
+  if (!is.null(guard)) {
+    k[1] <- guard
+  }
   graphon_quantiles <- qgraphon(k, a = a, Fv = Fv, X_matrix = X)
   
   ## 1.5 Build the matrix Q ====================================================