def METHOD_NAME(
    self, num_samples: Optional[int] = None, dtype=np.float32
) -> Tensor:
    r"""
    Draw samples from the multivariate Gaussian distributions.
    Internally, Cholesky factorization of the covariance matrix is used:
        sample = L v + mu,
    where L is the Cholesky factor, v is a standard normal sample.
    Parameters
    ----------
    num_samples
        Number of samples to be drawn.
    dtype
        Data-type of the samples.
    Returns
    -------
    Tensor
        Tensor with shape (num_samples, ..., d).
    """
    def s(mu: Tensor, L: Tensor) -> Tensor:
        F = self.F
        samples_std_normal = F.sample_normal(
            mu=F.zeros_like(mu),
            sigma=F.ones_like(mu),
            dtype=dtype,
        ).expand_dims(axis=-1)
        samples = (
            F.linalg_gemm2(L, samples_std_normal).squeeze(axis=-1) + mu
        )
        return samples
    return _sample_multiple(
        s, mu=self.mu, L=self.L, num_samples=num_samples
    )