diff --git a/R/compute_derivative.R b/R/compute_derivative.R
new file mode 100644
index 0000000..ae90020
--- /dev/null
+++ b/R/compute_derivative.R
@@ -0,0 +1,41 @@
+source(here::here("R", "create_Q.R"))
+
+calc_deriv_a <- function(k, i, matrix_X, a, K, qFa, fv){
+  # Quantile levels
+  q1 <- k / K
+  q0 <- (k - 1) / K
+  
+  # F^{-1}_a at the two quantiles
+  z1 <- qFa(q1)
+  z0 <- qFa(q0)
+  
+  # a' X_j for all j
+  aTX <- as.vector(matrix_X %*% a)
+  
+  # f_v(z - a'X_j)
+  f1 <- fv(z1 - aTX)
+  f0 <- fv(z0 - aTX)
+  
+  # X_j - X_i
+  Xi <- matrix_X[i, ]
+  XminusXi <- sweep(matrix_X, 2, Xi, "-")
+  
+  # compute fraction
+  term1 <- (fv(z1 - aTX[i]) / sum(f1)) * colSums(XminusXi * f1)
+  term0 <- (fv(z0 - aTX[i]) / sum(f0)) * colSums(XminusXi * f0)
+  
+  # Difference
+  return(term1 - term0)
+  
+}
+X <- matrix(seq(1,15, 1), ncol=3)
+a <- c(1, -2, 1)
+Fv <- pnorm
+K <- 3
+
+Fa <- make_distribution_func(a=a, X_matrix = X, Fv=Fv)
+Q <- create_matrix_Q(Fa, a, K, Fv, matrix_X = X)
+
+qFa <- make_quantile_function(Fa)
+
+deriv <- calc_deriv_a(1, 1, X, a, K, qFa, dnorm)
diff --git a/R/create_Q.R b/R/create_Q.R
index 7cad462..af91a03 100644
--- a/R/create_Q.R
+++ b/R/create_Q.R
@@ -1,3 +1,82 @@
+#' Construct the matrix Q from graphon quantiles and covariates
+#'
+#' Constructs the matrix \(Q\) whose entries are given by differences of
+#' \(F_v\) evaluated at graphon quantiles shifted by the inner products
+#' \(X_j^\top a\). Specifically, for \(k = 1,\ldots,K\) and
+#' \(j = 1,\ldots,n\),
+#'
+#' \deqn{
+#' Q_{kj} =
+#' F_v\left(\hat F_a^{-1}\left(\frac{k}{K}\right)
+#' - X_j^\top a\right)
+#' -
+#' F_v\left(\hat F_a^{-1}\left(\frac{k-1}{K}\right)
+#' - X_j^\top a\right).
+#' }
+#'
+#' Here, \code{qgraphon} is a function that returns the graphon quantiles
+#' \(\hat F_a^{-1}(u)\), \code{Fv} is the distribution function \(F_v\),
+#' \code{a} is the parameter vector, and \code{matrix_X} contains the
+#' covariate vectors \(X_j\) as rows.
+#'
+#' If \code{scaled = TRUE}, the resulting matrix is multiplied by
+#' \(1/\sqrt{n}\).
+#'
+#' @param qgraphon A function that computes the graphon quantile function.
+#'   It must accept a numeric vector of probabilities in \([0,1]\) and return
+#'   the corresponding quantiles.
+#' @param a A numeric parameter vector. Its length must equal the number of
+#'   columns of \code{matrix_X}.
+#' @param K A positive integer specifying the number of intervals used to
+#'   construct the matrix \(Q\).
+#' @param Fv A function representing the distribution function \(F_v\).
+#'   It must accept numeric input and return values of the same length.
+#' @param matrix_X A numeric matrix whose rows contain the covariate vectors
+#'   \(X_j\). The number of columns must equal \code{length(a)}.
+#' @param scaled Logical indicating whether the resulting matrix should be
+#'   scaled by \(1/\sqrt{n}\), where \(n\) is the number of rows of
+#'   \code{matrix_X}. Defaults to \code{FALSE}.
+#'
+#' @return A numeric \(K \times n\) matrix. The \((k,j)\)-th entry is
+#'   \deqn{
+#'   F_v\left(\hat F_a^{-1}(k/K) - X_j^\top a\right)
+#'   -
+#'   F_v\left(\hat F_a^{-1}((k-1)/K) - X_j^\top a\right).
+#'   }
+#'   If \code{scaled = TRUE}, the matrix is multiplied by \(1/\sqrt{n}\).
+#'
+#' @examples
+#' n <- 100
+#' K <- 3
+#' a <- c(2.0, -0.5)
+#'
+#' X <- matrix(
+#'   rnorm(2 * n),
+#'   nrow = n,
+#'   ncol = 2
+#' )
+#'
+#' Fv <- function(x) {
+#'   pnorm(x, mean = 0, sd = 1)
+#' }
+#'
+#' qgraphon <- make_distribution_func(
+#'   a = a,
+#'   Fv = Fv,
+#'   X_matrix = X
+#' )
+#'
+#' Q <- create_matrix_Q(
+#'   qgraphon = qgraphon,
+#'   a = a,
+#'   K = K,
+#'   Fv = Fv,
+#'   matrix_X = X
+#' )
+#'
+#' dim(Q)
+#'
+#' @export
 create_matrix_Q <- function(
     qgraphon,
     a,
@@ -13,7 +92,7 @@ create_matrix_Q <- function(
   if (!is.function(Fv)) stop("'F_v' must be a function")
   if (!is.matrix(matrix_X)) stop("matrix_X must be a matrix")
   if (!is.logical(scaled)) stop("`scaled` must be a logical!")
-  if (ncol(matrix_X) != length(a)) {
+  if (!is.null(matrix_X) && ncol(matrix_X) != length(a)) {
     stop("Number of columns of `matrix_X` (", ncol(matrix_X), ") must equal length(a) (", length(a), ")")
   }