A search for short-period Tausworthe generators over Fb with application to Markov chain quasi-Monte Carlo

Abstract

A one-dimensional sequence $u_0,u_1,u_2,\dots \in [0,1)$ is said to be completely uniformly distributed (CUD) if overlapping s-blocks $(u_i,u_{i+1},\dots ,u_{i+s-1})$, $i=0,1,2,\dots$, are uniformly distributed for every dimension $s\ge 1$. This concept naturally arises in Markov chain quasi-Monte Carlo (QMC). However, the definition of CUD sequences is not constructive, and thus there remains the problem of how to implement the Markov chain QMC algorithm in practice. Harase [A table of short-period Tausworthe generators for Markov chain quasi-Monte Carlo. J Comput Appl Math. 2021;384:Paper No. 113136, 12.] focussed on the t-value, which is a measure of uniformity widely used in the study of QMC, and implemented short-period Tausworthe generators (i.e. linear feedback shift register generators) over the two-element field ${\mathbb{F}}_2$ that approximate CUD sequences by running for the entire period. In this paper, we generalize a search algorithm over ${\mathbb{F}}_2$ to that over arbitrary finite fields ${\mathbb{F}}_b$ with b elements and conduct a search for Tausworthe generators over ${\mathbb{F}}_b$ with t-values zero (i.e. optimal) for dimension s = 3 and small for $s\ge 4$, especially in the case where b = 3, 4, and 5. We provide a parameter table of Tausworthe generators over ${\mathbb{F}}_4$, and report a comparison between our new generators over ${\mathbb{F}}_4$ and existing generators over ${\mathbb{F}}_2$ in numerical examples using Markov chain QMC.

Keywords:

• Pseudorandom number generation
• Quasi-Monte Carlo
• Markov chain Monte Carlo
• Bayesian inference
• Linear regression

1. Introduction

We study the problem of calculating the expectation $E_{\pi }\left[f\right(\mathbf{X}\left)\right]$ using Markov chain Monte Carlo (MCMC) methods for a target distribution π on a state space $\mathcal{X}$ and some function $f:\mathcal{X}\to \mathbb{R}$, where $\mathbf{X}$ is a π-distributed random variable on $\mathcal{X}$. We are interested in improving the accuracy by replacing IID uniform random points with quasi-Monte Carlo (QMC) points. However, traditional QMC points (e.g.Sobol', Niederreiter–Xing, Faure, and Halton) are not straightforwardly applicable. Motivated by a simulation study conducted by Liao [1], Owen and Tribble [2] and Chen et al. [3] theoretically showed that Markov chain QMC remains consistent if the driving sequences are completely uniformly distributed (CUD). A one-dimensional sequence $u_0,u_1,u_2,\dots$ $\in [0,1)$ is said to be CUD if overlapping s-blocks $(u_i,u_{i+1},\dots ,u_{i+s-1})$, $i=0,1,2,\dots$, are uniformly distributed for every dimension $s\ge 1$. Levin [4] proposed some constructions of CUD sequences, but they are not suitable to implement. Thus, there remains the problem of how we implement the Markov chain QMC algorithm, in particular, how we construct suitable driving sequences in practice.

Tribble and Owen [5] and Tribble [6] proposed an implementation method to obtain point sets that approximate CUD sequences by using short-period linear congruential and Tausworthe generators (i.e.linear feedback shift register generators over the two-element field ${\mathbb{F}}_2:=\{0,1\}$) that run for the entire period. Moreover, Chen et al. [7] implemented short-period Tausworthe generators in terms of the equidistribution property, which is a coarse measure of uniformity in the area of pseudorandom number generation [8].

In a previous study, Harase [9] implemented short-period Tausworthe generators over ${\mathbb{F}}_2$ that approximate CUD sequences in terms of the t-value, which is a central measure in the theory of $(t,m,s)$-nets and $(t,s)$-sequences. The key technique was to use a polynomial analogue of Fibonacci numbers and their continued fraction expansion, which was originally proposed by Tezuka and Fushimi [10]. More precisely, we can view Tausworthe generators as a polynomial analogue of Korobov lattice rules with a denominator polynomial $p\left(x\right)\in {\mathbb{F}}_2\left[x\right]$ and a numerator polynomial $q\left(x\right)\in {\mathbb{F}}_2\left[x\right]$ (cf. [11,12]), and hence, the t-value is zero (i.e. optimal) for dimension s = 2 if and only if the partial quotients in the continued fraction of $q\left(x\right)/p\left(x\right)$ are all of degree one [10,13]. By enumerating such pairs of polynomials $\left(p\right(x),q(x\left)\right)$ efficiently, Harase [9] conducted an exhaustive search of Tausworthe generators over ${\mathbb{F}}_2$ with t-values zero for s = 2 and small (but not zero) for $s\ge 3$, and demonstrated the effectiveness in numerical examples using Gibbs sampling.

From the theoretical and practical perspective, the most interesting case is the t-value zero. However, Kajiura et al. [14] proved that there exists no maximal-period Tausworthe generator over ${\mathbb{F}}_2$ with t-value zero for s = 3. In fact, in finite fields ${\mathbb{F}}_b$ of prime power order $b\ge 3$, we can find maximal-period Tausworthe generators with t-value zero for s = 3, for some combinations of b and m.

In this paper, our aim is to conduct an exhaustive search of maximal-period Tausworthe generators over ${\mathbb{F}}_b$ with t-values zero for dimension s = 3, in addition to s = 2, especially in the case where b = 3, 4, and 5. For this purpose, we generalize the search algorithms of Tezuka and Fushimi [10] and Harase [9] over ${\mathbb{F}}_2$ to those over ${\mathbb{F}}_b$. We provide a parameter table of Tausworthe generators over ${\mathbb{F}}_4$ with t-values zero for s = 3 and small for $s\ge 4$ to implement the Markov chain QMC algorithm. Accordingly, we report a comparison between our new Tausworthe generators over ${\mathbb{F}}_4$ and existing generators [7,9] over ${\mathbb{F}}_2$ in numerical examples using Markov chain QMC.

The rest of this paper is organized as follows: In Section 2, we recall the definitions of CUD sequences, Tausworthe generators, and the t-value, and recall a connection between the t-value and continued fraction expansion. In Section 3, we discuss our main results: In Section 3.1, we investigate the number of polynomials $q\left(x\right)$ for which the partial quotients of the continued fraction expansion of $q\left(x\right)/p\left(x\right)$ all have degree one for a given irreducible polynomial $p\left(x\right)$ over ${\mathbb{F}}_b$. In Section 3.2, we describe a search algorithm of Tausworthe generators over ${\mathbb{F}}_b$. In Section 3.3, we conduct an exhaustive search in the case where b = 3, 4, and 5, and provide tables. In Section 4, we present numerical examples, such as Gibbs sampling and a simulation of a queuing system, in which both Tausworthe generators over ${\mathbb{F}}_2$ and ${\mathbb{F}}_4$ optimized in terms of the t-value perform comparable to or better than Tausworthe generators [7] optimized in terms of the equidistribution property. In Section 5, we conclude this paper.

2. Preliminaries

We refer the reader to [8,9,11,13,15] for general information.

2.1. Discrepancy and completely uniformly distributed sequences

Let $P_s=\{{\mathbf{u}}_0,{\mathbf{u}}_1,\dots ,{\mathbf{u}}_{N-1}\}\subset [0,1{)}^s$ be an s-dimensional point set of N elements in the sense of a multiset. Let us recall the definition of discrepancy $D_N^{\ast s}\left(P_s\right)$ as a measure of uniformity.

Definition 2.1

Discrepancy

For a point set $P_s=\{{\mathbf{u}}_0,{\mathbf{u}}_1,\dots ,{\mathbf{u}}_{N-1}\}\subset [0,1{)}^s$, the (star) discrepancy is defined as $D_N^{\ast s}\left(P_s\right):=\underset{J}{sup}|\frac{\nu (J;P_s)}{N}-\mathrm{vol}(J\left)\right|,$where the supremum is taken over all intervals J of the form $\prod _{j=1}^s[0,t_j)$ for $0<t_j\le 1$, $\nu (J;P_s)$ denotes the number of i with $0\le i\le N-1$ for which ${\mathbf{u}}_i\in J$, and $\mathrm{vol}\left(J\right):=\prod _{j=1}^s{t}_j$ denotes the volume of J.

We define the CUD property for a one-dimensional sequence $\{u_i{\}}_{i=0}^{\mathrm{\infty }}$ in $[0,1)$, which is known as one of the definitions of random number sequences in [16, Chapter 3.5].

Definition 2.2

CUD sequences

A one-dimensional infinite sequence $u_0,u_1,$ $u_2,\dots$ $\in [0,1)$ is said to be completely uniformly distributed (CUD) if overlapping s-blocks satisfy $\underset{N\to \mathrm{\infty }}{lim}{D}_N^{\ast s}\left(\right(u_0,\dots ,u_{s-1}),(u_1,\dots ,u_s),\dots ,(u_{N-1},\dots ,u_{N+s-2}\left)\right)=0$for every dimension $s\ge 1$; in short, the sequence of s-blocks $(u_i,\dots ,u_{i+s-1}),$ $i=0,1,\dots$, is uniformly distributed in $[0,1{)}^s$ for every dimension $s\ge 1$.

In the study of Markov chain QMC, it is desirable that $D_N^{\ast s}$ converges to zero as fast as possible if $N\to \mathrm{\infty }$ (cf. [17,18]). As a necessary and sufficient condition of Definition 2.2, Chentsov [19] proved the following theorem:

Theorem 2.3

[19]

A one-dimensional infinite sequence $u_0,u_1,$ $u_2,\dots$ $\in [0,1)$ is CUD if and only if non-overlapping s-blocks satisfy (1) $\underset{N\to \mathrm{\infty }}{lim}{D}_N^{\ast s}\left(\right(u_0,\dots ,u_{s-1}),(u_s,\dots ,u_{2s-1}),\dots ,(u_{(N-1)s},\dots ,u_{\mathrm{Ns}-1}\left)\right)=0$(1) for every dimension $s\ge 1$.

We thus use a sequence $\{u_i{\}}_{i=0}^{\mathrm{\infty }}$ in $[0,1)$ for Markov chain QMC in this order.

2.2. Tausworthe generators over ${\mathbb{F}}_b$

Let ${\mathbb{F}}_b$ be a finite field with b elements, where b is a prime power, and perform addition and multiplication over ${\mathbb{F}}_b$. We define Tausworhe generators over ${\mathbb{F}}_b$, which are usually defined over ${\mathbb{F}}_2$ [8,20–22].

Definition 2.4

Tausworthe generators over ${\mathbb{F}}_b$

Let $p\left(x\right):=x^m-c_1{x}^{m-1}-\cdots -c_{m-1}x-c_m\in {\mathbb{F}}_b\left[x\right]$, where $c_m\ne 0$. We consider the linear recurrence over ${\mathbb{F}}_b$ given by (2) $a_i:=c_1{a}_{i-1}+\cdots +c_m{a}_{i-m}\in {\mathbb{F}}_b,i=0,1,2,\dots ,$(2) whose characteristic polynomial is $p\left(x\right)$. Let σ be a step size and w a digit number. We define the output $u_i\in [0,1)$ at step i as (3) $u_i:=\sum _{j=0}^{w-1}\eta \left(a_{\mathrm{i\sigma }+j}\right)b^{-j-1}\in [0,1),i=0,1,2,\dots ,$(3) where $\eta :{\mathbb{F}}_b\to {\mathbb{Z}}_b:=\{0,1,\dots ,b-1\}$ is a bijection with $\eta \left(0\right)=0$. If $p\left(x\right)$ is primitive, $(a_0,\dots ,a_{m-1})\ne (0,\dots ,0)$, $0<\sigma <b^m-1$, and $gcd(\sigma ,b^m-1)=1$, then the sequences (2(2) $a_i:=c_1{a}_{i-1}+\cdots +c_m{a}_{i-m}\in {\mathbb{F}}_b,i=0,1,2,\dots ,$(2) ) and (3(3) $u_i:=\sum _{j=0}^{w-1}\eta \left(a_{\mathrm{i\sigma }+j}\right)b^{-j-1}\in [0,1),i=0,1,2,\dots ,$(3) ) are both purely periodic with maximal period $b^m-1$. Throughout this paper, we assume these maximal-period conditions. We call a generator in such a class a Tausworthe generator over ${\mathbb{F}}_b$ (or a linear feedback shift register generator over ${\mathbb{F}}_b$).

Similar to the case of ${\mathbb{F}}_2$, Tausworthe generators over ${\mathbb{F}}_b$ can be viewed as a polynomial analogue of linear congruential generators (LCGs): $\begin{array}{rl}{X}_i\left(x\right)& :=q\left(x\right)X_{i-1}\left(x\right)\mathrm{mod}p\left(x\right),\\ X_i\left(x\right)/p\left(x\right)& =a_{\mathrm{i\sigma }}{x}^{-1}+a_{\mathrm{i\sigma }+1}{x}^{-2}+a_{\mathrm{i\sigma }+2}{x}^{-3}+\cdots \in {\mathbb{F}}_b\left(\right(x^{-1}\left)\right)\end{array}$where $X_i\left(x\right)\in {\mathbb{F}}_b\left[x\right],i=0,1,2,\dots$, is a sequence of polynomials, $p\left(x\right),q\left(x\right)\in {\mathbb{F}}_b\left[x\right]$ represent a modulus and multiplier, respectively, and the step size σ satisfies $q\left(x\right)=x^{\sigma }\mathrm{mod}p\left(x\right)$ and $0<\sigma <b^m-1$. Then, the output $u_i$ in (3(3) $u_i:=\sum _{j=0}^{w-1}\eta \left(a_{\mathrm{i\sigma }+j}\right)b^{-j-1}\in [0,1),i=0,1,2,\dots ,$(3) ) is expressed as $u_i={\nu }_w\left(X_i\right(x\left)/p\right(x\left)\right)$, where a map ${\nu }_w:{\mathbb{F}}_b\left(\right(x^{-1}\left)\right)\to [0,1)$ is given by $\sum _{j=j_0}^{\mathrm{\infty }}{k}_j{x}^{-j-1}\mapsto \sum _{j=max\{0,j_0\}}^{w-1}\eta \left(k_j\right)b^{-j-1}$, which transforms a formal power series in ${\mathbb{F}}_b\left(\right(x^{-1}\left)\right)$ into a b-adic expansion with w digits in $[0,1)$.

Moreover, similar to LCGs, Tausworthe generators have a lattice structure. Let $N:=b^m$. We consider a sequence (4) $u_0,u_1,\dots ,u_{N-2},u_{N-1}=u_0,u_1,\dots \in [0,1)$(4) generated by a Tausworthe generator (3(3) $u_i:=\sum _{j=0}^{w-1}\eta \left(a_{\mathrm{i\sigma }+j}\right)b^{-j-1}\in [0,1),i=0,1,2,\dots ,$(3) ) with period length N−1. We set s-dimensional overlapping points ${\mathbf{u}}_i=(u_i,\dots ,u_{i+s-1})$ for $i=0,1,\dots ,N-2$, that is, ${\mathbf{u}}_0=(u_0,\dots ,u_{s-1}),{\mathbf{u}}_1=(u_1,\dots ,u_s),\dots ,{\mathbf{u}}_{N-2}=(u_{N-2},u_0,\dots ,u_{s-2})$. We construct a QMC point set (5) $P_s=\left\{\mathbf{0}\right\}\cup \{{\mathbf{u}}_i{\}}_{i=0}^{N-2}\subset [0,1{)}^s,$(5) adding the origin $\left\{\mathbf{0}\right\}$. Note that the cardinality is $\mathrm{\#}{P}_s=b^m$. Then, a point set $P_s$ in (5(5) $P_s=\left\{\mathbf{0}\right\}\cup \{{\mathbf{u}}_i{\}}_{i=0}^{N-2}\subset [0,1{)}^s,$(5) ) can be viewed as a polynomial analogue of Korobov lattice rules: (6) $P_s=\{{\nu }_w\left(\frac{h\left(x\right)}{p\left(x\right)}\right(1,q\left(x\right),q(x{)}^2,\dots ,q(x{)}^{s-1}\left)\right)|\mathrm{deg}(h\left(x\right))<m\},$(6) where the map ${\nu }_w$ is applied component-wise and $m=\mathrm{deg}\left(p\right(x\left)\right)$; see [11,12] for details.

A pair of polynomials $\left(p\right(x),q(x\left)\right)$ is a parameter set of Tausworthe generators. Thus, in accordance with Definition 2.2, we would like to find a pair of polynomials $\left(p\right(x),q(x\left)\right)$ with small discrepancy $D_N^{\ast s}\left(P_s\right)$ for each $s\ge 1$.

2.3. t-value and continued fraction expansion

A point set $P_s$ in (5(5) $P_s=\left\{\mathbf{0}\right\}\cup \{{\mathbf{u}}_i{\}}_{i=0}^{N-2}\subset [0,1{)}^s,$(5) ) and (6(6) $P_s=\{{\nu }_w\left(\frac{h\left(x\right)}{p\left(x\right)}\right(1,q\left(x\right),q(x{)}^2,\dots ,q(x{)}^{s-1}\left)\right)|\mathrm{deg}(h\left(x\right))<m\},$(6) ) generated by a Tausworthe generator (3(3) $u_i:=\sum _{j=0}^{w-1}\eta \left(a_{\mathrm{i\sigma }+j}\right)b^{-j-1}\in [0,1),i=0,1,2,\dots ,$(3) ) is a digital net. Hence, we can compute the t-value, which is closely related to $D_N^{\ast s}\left(P_s\right)$ for $N=b^m$.

Definition 2.5

$(t,m,s)$-nets

Let $s\ge 1$ and $0\le t\le m$ denote integers. A point set $P_s$ of $b^m$ points in $[0,1{)}^s$ is said to be a $(t,m,s)$-net in base b if every interval of the form $E=\prod _{j=1}^s[r_j/b^{d_j},(r_j+1)/b^{d_j})$ in $[0,1{)}^s$ with integers $d_j\ge 0$ and $0\le r_j<b^{d_j}$ and of volume $b^{t-m}$ contains exactly $b^t$ points from $P_s$.

For a given dimension s, the smallest value t for which $P_s$ is a $(t,m,s)$-net is said to be the t-value. For a $(t,m,s)$-net $P_s$ in base b, we have an upper bound $D_N^{\ast s}\left(P_s\right)\le C_{b,s}{b}^t(\log N{)}^{s-1}/N$, where the constant $C_{b,s}>0$ only depends on s and b; hence a small t-value is desirable. Therefore, we adopt the t-value as a measure of uniformity instead of the direct calculation of $D_N^{\ast s}\left(P_s\right)$ to obtain low-discrepancy point sets.

Furthermore, in the case s = 2, there is a connection between the t-value of a polynomial Korobov lattice rule (6(6) $P_s=\{{\nu }_w\left(\frac{h\left(x\right)}{p\left(x\right)}\right(1,q\left(x\right),q(x{)}^2,\dots ,q(x{)}^{s-1}\left)\right)|\mathrm{deg}(h\left(x\right))<m\},$(6) ) and the continued fraction expansion of $q\left(x\right)/p\left(x\right)$. Let $q\left(x\right)/p\left(x\right)$ be a rational function over ${\mathbb{F}}_b$ with $gcd\left(p\right(x),q(x\left)\right)=1$ and $\mathrm{deg}\left(p\right(x\left)\right)\ge 1$. Then, $q\left(x\right)/p\left(x\right)$ has a unique regular continued fraction expansion $\begin{array}{rl}q\left(x\right)/p\left(x\right)& =A_{v+1}\left(x\right)+1/\left(A_v\right(x)+1/(A_{v-1}\left(x\right)+\cdots +1/A_1\left(x\right)\left)\right)\\ & =\left[A_{v+1}\right(x);A_v(x),A_{v-1}(x),\dots ,A_1(x\left)\right]\end{array}$with a polynomial part $A_{v+1}\left(x\right)\in {\mathbb{F}}_b\left[x\right]$ and partial quotients $A_k\left(x\right)\in {\mathbb{F}}_b\left[x\right]$ satisfying $\mathrm{deg}\left(A_k\right(x\left)\right)\ge 1$ for $1\le k\le v$. Under this condition, we put $K\left(q/p\right):=\underset{1\le k\le v}{max}\mathrm{deg}\left(A_k\right(x\left)\right)$. We have the following theorem:

Theorem 2.6

[10,13]

Let $p\left(x\right)\in {\mathbb{F}}_b\left[x\right]$ with $m=\mathrm{deg}\left(p\right(x\left)\right)$. Let $q\left(x\right)\in {\mathbb{F}}_b\left[x\right]$ with $\mathrm{deg}\left(q\right(x\left)\right)<m$. Suppose that $gcd\left(p\right(x),q(x\left)\right)=1$. Then, the two-dimensional point set (7) $P_2=\{{\nu }_w\left(\frac{h\left(x\right)}{p\left(x\right)}\right(1,q\left(x\right)\left)\right)|\mathrm{deg}(h\left(x\right))<m\}$(7) is a $(t,m,2)$-net in base b with $t=K\left(q/p\right)-1$, which is exactly the t-value. In particular, $P_2$ has the t-value zero if and only if $K\left(q/p\right)=1$, so $\mathrm{deg}\left(A_k\right(x\left)\right)=1$ for all $1\le k\le v$ and v = m.

Using the continued fraction expansion based on the above theorem, Tezuka and Fushimi [10] proposed an algorithm to search for Tausworthe generators over ${\mathbb{F}}_2$ having pairs of polynomials $\left(p\right(x),q(x\left)\right)$ with t-value zero for s = 2 and small for $s\ge 3$. Harase [9] recently indicated that their technique is applicable to QMC points that approximate CUD sequences and conducted an exhaustive search over ${\mathbb{F}}_2$ removing some conditions.

Remark 2.1

In previous studies, L'Ecuyer and Lemieux [12,23] constructed short-period Tausworthe generators for QMC numerical integration in general-purpose use. To assess the uniformity of QMC points, they took into account the quality of the projections and developed several figures of merit using the equidistribution property, which are often used for selecting pseudorandom number generators with very long period [20,21,24]. These figures of merit are implemented in LatNet Builder [25], a software tool to find good parameters, and are probably useful in our study, but they are not so closely related to the discrepancy $D_N^{\ast s}$ as the t-value because the condition of the equidistribution property is sometimes weaker than that of the t-value. The CUD sequences are defined via the discrepancy $D_N^{\ast s}$ in Definition 2.1, so we adopt the t-value as a primary criterion. We also note that our study is aimed at an application to Markov chain QMC, not usual pseudorandom number generation.

3. Main results

In the theory of $(t,m,s)$-nets and $(t,s)$-sequences, the most interesting case is the t-value zero. Kajiura et al. [14] proved that there exists no maximal-period Tausworthe generator over ${\mathbb{F}}_2$ with t-value zero for dimension s = 3. Thus, we conduct a search of Tausworthe generators over ${\mathbb{F}}_b$ with t-value zero for s = 3, especially in the case where b = 3, 4, and 5.

3.1. Orthogonal multiplicity

To obtain a pair of polynomials $\left(p\right(x),q(x\left)\right)$ with t-value zero for $s=3$, it is necessary to satisfy at least $K\left(q/p\right)=1$ in Theorem 2.6. Thus, we first investigate how many polynomials $q\left(x\right)$ satisfying $K\left(q/p\right)=1$ exist for each irreducible polynomial $p\left(x\right)\in {\mathbb{F}}_b\left[x\right]$. For a given irreducible $p\left(x\right)$, we define the number $M\left(p\right):=\mathrm{\#}\left\{q\right(x)\in {\mathbb{F}}_b[x\left]\right|\mathrm{deg}\left(q\right(x\left)\right)<\mathrm{deg}\left(p\right(x\left)\right)\text{ and }K\left(q/p\right)=1\}.$The number $M\left(p\right)$ is called the orthogonal multiplicity of $p\left(x\right)$ in [26]. Specializing the proof for the case ${\mathbb{F}}_2$, Mesirov and Sweet [27] proved that every irreducible polynomial $p\left(x\right)\in {\mathbb{F}}_2\left[x\right]$ has exactly $M\left(p\right)=2$ for $\mathrm{deg}\left(p\right(x\left)\right)\ge 2$, that is, there exist only two polynomials $q\left(x\right)$ for which the partial quotients of $q\left(x\right)/p\left(x\right)$ have all degree one. Moreover, such polynomials are $q\left(x\right)$ and its inverse element $q^{-1}\left(x\right)\mathrm{mod}p\left(x\right)$, and hence, they yield exactly the same lattice point set $P_s$. This result asserts the existence of $P_2$ with t-value zero for every irreducible polynomial $p\left(x\right)\in {\mathbb{F}}_2\left[x\right]$ in Theorem 2.6 but also asserts that there is no degree of freedom to select such $q\left(x\right)$ for each $p\left(x\right)$.

In fact, in the case ${\mathbb{F}}_b$ for $b\ge 3$, Blackburn [26] indicated that the situation is different far from the case ${\mathbb{F}}_2$. More precisely, the orthogonal multiplicities $M\left(p\right)$ are not always the same number but are often much greater than two. Figure 1 shows some histograms of orthogonal multiplicities $M\left(p\right)$ for all monic irreducible polynomials $p\left(x\right)\in {\mathbb{F}}_b\left[x\right]$ with $\mathrm{deg}\left(p\right(x\left)\right)=m$. No clear regularity as in ${\mathbb{F}}_2$ has been observed. Additionally, for arbitrary ${\mathbb{F}}_b$ with $b\ge 3$, it is not even known whether there exist irreducible polynomials $p\left(x\right)\in {\mathbb{F}}_b\left[x\right]$ with $M\left(p\right)=0$ in general (see Remark 3.1). Thus, using computer calculations, we checked the existence of $M\left(p\right)>0$ as follows:

Figure 1. Distribution of orthogonal multiplicities $M\left(p\right)$ for all monic irreducible polynomials $p\left(x\right)\in {\mathbb{F}}_b\left[x\right]$ with $\mathrm{deg}\left(p\right(x\left)\right)=m$.

Theorem 3.1

Let ${\mathbb{F}}_b$ be a finite field with b elements. Every monic irreducible polynomial $p\left(x\right)\in {\mathbb{F}}_b$ with $\mathrm{deg}\left(p\right(x\left)\right)=m$ has $M\left(p\right)>0$, at least under the following conditions:

• $1\le m\le 15$ for b = 3;

• $1\le m\le 12$ for b = 4;

• $1\le m\le 10$ for b = 5.

Remark 3.1

Let $p\left(x\right)\in {\mathbb{F}}_b\left[x\right]$ be an irreducible polynomial over ${\mathbb{F}}_b$. Assume that $0<\mathrm{deg}\left(p\right(x\left)\right)<b$, that is, $\mathrm{deg}\left(p\right(x\left)\right)$ is less than the order b of ${\mathbb{F}}_b$. Under this condition, Friesen [28, Theorem 2] proved that every irreducible $p\left(x\right)\in {\mathbb{F}}_b\left[x\right]$ has $M\left(p\right)>0$, that is, every irreducible $p\left(x\right)$ has $M\left(p\right)>0$ provided the order b of ${\mathbb{F}}_b$ is sufficiently large. This result is an improvement of that in the study by Blackburn [26, Theorem 2]. However, the assumption $\mathrm{deg}\left(p\right(x\left)\right)<b$ is significantly restrictive compared with the numerical results, and there has been no progress on the study of orthogonal multiplicities $M\left(p\right)$ since Friesen's study. Thus, we numerically checked the existence of $M\left(p\right)>0$ only in the range required for our study.

3.2. A search algorithm using Fibonacci polynomials over ${\mathbb{F}}_b$

Tausworthe generators associated with $\left(p\right(x),q(x\left)\right)$ attain the t-value zero for s = 3 only if $K\left(q/p\right)=1$ in Theorem 2.6. Our strategy is to choose $\left(p\right(x),q(x\left)\right)$ with t-value zero for s = 3 among pairs satisfying $K\left(q/p\right)=1$. Thus, we generalize the search algorithms of Tezuka and Fushimi [10] and Harase [9] over ${\mathbb{F}}_2$ to those over arbitrary finite fields ${\mathbb{F}}_b$.

Recall that Fibonacci numbers $F_k,k=1,2\dots$, are defined by the recurrence $F_k=F_{k-1}+F_{k-2}$, where we choose the starting values $F_{-1}=0,F_0=1$. Then, the continued fraction expansion of the ratio of two successive Fibonacci numbers $F_{k-1}/F_k$ is given by $F_{k-1}/F_k=[0;1,1,\dots ,1]$ with partial quotients that are all one. As a polynomial analogue, we consider a sequence of polynomials $F_k\left(x\right),k=1,2,\dots$, defined as (8) $\begin{array}{rl}& F_k\left(x\right)=A_k\left(x\right)F_{k-1}\left(x\right)+F_{k-2}\left(x\right),\end{array}$(8) (9) $\begin{array}{rl}& F_{-1}\left(x\right)=0,F_0\left(x\right)=1,\end{array}$(9) (10) $\begin{array}{rl}& A_k\left(x\right)=\text{β}x+\gamma ,\end{array}$(10) where $\text{β}\in {\mathbb{F}}_b^{\ast }:={\mathbb{F}}_b\mathrm{\setminus }\left\{0\right\}$ and $\gamma \in {\mathbb{F}}_b$ so that $\mathrm{deg}\left(A_k\right(x\left)\right)=1$. Similarly, we have the continued fraction expansion $F_{k-1}\left(x\right)/F_k\left(x\right)=[0;A_k(x),A_{k-1}(x),\dots ,A_1(x\left)\right]$, so $K\left(F_{k-1}/F_k\right)=1$ holds. The polynomials $F_k\left(x\right),k=0,1,2,\dots$, are called Fibonacci polynomials over ${\mathbb{F}}_b$ (cf. [29]). Figure 2 shows an example of the initial part of a tree of Fibonacci polynomials over ${\mathbb{F}}_3$. Note that there exist $\left\{\right(b-1)b{\}}^m$ different pairs $\left(F_m\right(x),F_{m-1}(x\left)\right)$ of Fibonacci polynomials over ${\mathbb{F}}_b$. Among pairs $\left(F_m\right(x),F_{m-1}(x\left)\right)$, we choose a suitable pair of $\left(p\right(x),q(x\left)\right)$ with t-values zero for s = 3 and small for $s\ge 4$ satisfying Definition 2.4.

Figure 2. Initial part of the tree of Fibonacci polynomials over ${\mathbb{F}}_3$.

We now generalize the algorithms [9,10] over ${\mathbb{F}}_2$ to those over ${\mathbb{F}}_b$. Let $\mathrm{lc}\left(F_m\right(x\left)\right)$ denote the leading coefficient of $F_m\left(x\right)$ and $s_{\max }$ denote a given maximum dimensionality. Our algorithm proceeds as follows:

Before we begin our algorithm, we create a lookup table of primitive polynomials over ${\mathbb{F}}_b$ in advance to avoid repeated computation in Step 3. In Step 2, $F_m\left(x\right)$ generated by (8(8) $\begin{array}{rl}& F_k\left(x\right)=A_k\left(x\right)F_{k-1}\left(x\right)+F_{k-2}\left(x\right),\end{array}$(8) ) is not always monic over arbitrary finite fields ${\mathbb{F}}_b$ except for ${\mathbb{F}}_2$, so it is necessary to divide $F_m\left(x\right)$ and $F_{m-1}\left(x\right)$ by the leading coefficient $\mathrm{lc}\left(F_m\right(x\left)\right)$. In Steps 5 and 6, we compute the t-values using Gaussian elimination [30]. For some combinations of b and m, $\left(F_m\right(x),F_{m-1}(x\left)\right)$ might not exist with t-value zero for s = 3 in Step 5. In this case, we skip Steps 6–9 and terminate the algorithm.

Remark 3.2

Tezuka and Fushimi [10] and Harase [9] dealt with the search algorithms that are similar to Algorithm 1 but restricted to the special case ${\mathbb{F}}_2$. We now note that there are several differences between the cases ${\mathbb{F}}_2$ and ${\mathbb{F}}_b$ for $b\ge 3$. With regard to Equation (10(10) $\begin{array}{rl}& A_k\left(x\right)=\text{β}x+\gamma ,\end{array}$(10) ), we have only two polynomials $A_k\left(x\right)$ with degree one in the case ${\mathbb{F}}_2$, that is, $A_k\left(x\right)=x$ or x + 1; but we have many polynomials with degree one in general cases ${\mathbb{F}}_b$ (e.g. see Figure 2). Thus, the patterns of continued fraction expansions $[0;A_k(x),A_{k-1}(x),\dots ,A_1(x\left)\right]$ drastically increase as opposite to the case ${\mathbb{F}}_2$. Moreover, every polynomial over ${\mathbb{F}}_2$ is always monic, and hence, Step 2 in Algorithm 1 is not appeared in the existing algorithms. Once again, as mentioned in Section 3.1, there are only two polynomials $q\left(x\right)$ for which the partial quotients of $q\left(x\right)/p\left(x\right)$ have all degree one over ${\mathbb{F}}_2$, but many polynomials $q\left(x\right)$ with such property exist in the case ${\mathbb{F}}_b$. Therefore, our generalization would not be straightforward and simple when we search for parameters in practice.

3.3. Specific parameters

We conduct an exhaustive search of short-period Tausworthe generators over ${\mathbb{F}}_3$, ${\mathbb{F}}_4$, and ${\mathbb{F}}_5$ using Algorithm 1. We set $s_{\max }=20$. If b is a prime number (i.e, b = 3 or 5), we identify ${\mathbb{F}}_b$ with ${\mathbb{Z}}_b$ and set a bijection $\eta :{\mathbb{F}}_b\to {\mathbb{Z}}_b$ as the identity map. If b = 4, we set ${\mathbb{F}}_4=\{0,1,\alpha ,{\alpha }^2\}$ with ${\alpha }^2=\alpha +1$ and ${\alpha }^3=1$ and set a bijection $\eta :{\mathbb{F}}_4\to {\mathbb{Z}}_4$ consisting of $0\mapsto 0,1\mapsto 1,\alpha \mapsto 2,{\alpha }^2=\alpha +1\mapsto 3.$Table 1 summarizes the number of maximal-period Tausworthe generators with t-value zero for dimension s = 3. We observe that a very few pairs of polynomials $\left(p\right(x),q(x\left)\right)$ exist over ${\mathbb{F}}_3$; however, many pairs exist over ${\mathbb{F}}_4$ and ${\mathbb{F}}_5$, at least within the range described in the table. From the viewpoint of applications, we tabulate specific parameters of pairs of polynomials $\left(p\right(x),q(x\left)\right)$ over ${\mathbb{F}}_4$ and step sizes σ for $2\le m\le 11$ in Table 2. In Table 2, each first and second row shows the coefficients of $p\left(x\right)$ and $q\left(x\right)$ respectively; for example, ${\alpha }^{2}11$ means ${\alpha }^2+x+x^2$. Table 3 shows the t-values in the range of $1\le s\le 20$. Throughout our search, we find several parameters with the same t-values, so we choose one from them.

Table 1. Number of pairs of polynomials $\left(p\right(x),q(x\left)\right)$ that attain maximal-period Tausworthe generators with t-value zero for dimension $s=3$.

Number of Tausworthe generators over ${\mathbb{F}}_3$ with t-value zero.
m	2	3	4	5	6	7	8	9	10	11	12	13
Num.	8	6	0	0	8	6	0	0	0	0	0	0
Number of Tausworthe generators over ${\mathbb{F}}_4$ with t-value zero.
m	2	3	4	5	6	7	8	9	10	11
Num.	32	72	128	1296	2016	7648	4640	5328	4176	4560
Number of Tausworthe generators over ${\mathbb{F}}_5$ with t-value zero.
m	2	3	4	5	6	7	8
Num.	32	480	1056	16,800	38,720	514,640	706,496

Table 2. Specific parameters of pairs of polynomials $\left(p\right(x),q(x\left)\right)$ over ${\mathbb{F}}_4$ and step sizes σ.

m = 2	${\alpha }^{2}11$
	$\alpha 1$  ($\sigma =8$)
m = 3	${\alpha }^2{\alpha }^2{\alpha }^{2}1$
	$1\alpha {\alpha }^2$  ($\sigma =47$)
m = 4	${\alpha }^2{\alpha }^2{\alpha }^{2}01$
	${\alpha }^{2}11{\alpha }^2$  ($\sigma =131$)
m = 5	${\alpha }^2{\alpha }^2\alpha 101$
	$\alpha {\alpha }^2{\alpha }^2{\alpha }^2{\alpha }^2$  ($\sigma =724$)
m = 6	${\alpha }^{2}101101$
	$11{\alpha }^2{\alpha }^{2}1\alpha$  ($\sigma =2267$)
m = 7	$\alpha {\alpha }^{2}0\alpha {\alpha }^2\alpha \alpha 1$
	$00{\alpha }^2{\alpha }^2\alpha {\alpha }^{2}1$  ($\sigma =1633$)
m = 8	$\alpha {\alpha }^{2}110\alpha 001$
	$111100\alpha {\alpha }^2$  ($\sigma =16423$)
m = 9	${\alpha }^2{\alpha }^2\alpha 01\alpha \alpha 101$
	$\alpha 11{\alpha }^2{\alpha }^2{\alpha }^2\alpha 01$  ($\sigma =36887$)
m = 10	$\alpha {\alpha }^2\alpha 01{\alpha }^{2}00{\alpha }^{2}01$
	${\alpha }^{2}00\alpha 101111$  ($\sigma =1030108$)
m = 11	${\alpha }^2$  α  1  ${\alpha }^2$  α  ${\alpha }^2$  1  ${\alpha }^2$  ${\alpha }^2$  1  α  1
	${\alpha }^2$  α  ${\alpha }^2$  α  α  ${\alpha }^2$  1  ${\alpha }^2$  1  1  α  ($\sigma =3144209$)

Table 3. The t-values for good Tausworthe generators over ${\mathbb{F}}_4$.

$m\mathrm{\setminus s}$	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16	17	18	19	20
2	0	0	0	0	0	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1
3	0	0	0	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1
4	0	0	0	1	1	1	1	2	2	2	2	2	2	2	2	2	2	2	2	2
5	0	0	0	1	1	2	2	2	2	2	2	2	2	2	2	2	2	2	2	2
6	0	0	0	1	2	2	2	2	3	3	3	3	3	3	3	3	3	3	3	3
7	0	0	0	1	2	2	2	3	3	3	3	3	3	4	4	4	4	4	4	4
8	0	0	0	1	2	4	4	4	4	4	4	4	4	4	4	4	4	4	4	4
9	0	0	0	1	3	3	3	3	3	4	4	4	4	4	4	5	5	5	5	5
10	0	0	0	2	2	3	3	3	4	4	4	5	5	6	6	6	6	6	6	6
11	0	0	0	2	3	3	3	4	5	5	5	5	5	5	5	5	5	6	6	6

For the implementation, we introduce a reasonably fast algorithm to generate the output values (3(3) $u_i:=\sum _{j=0}^{w-1}\eta \left(a_{\mathrm{i\sigma }+j}\right)b^{-j-1}\in [0,1),i=0,1,2,\dots ,$(3) ) from Tausworthe generators over ${\mathbb{F}}_4$. Assume that $m\le w$. Let ${\tilde{\mathbf{x}}}_i=(a_{\mathrm{i\sigma }},a_{\mathrm{i\sigma }+1},\dots ,a_{\mathrm{i\sigma }+m-1},a_{\mathrm{i\sigma }+m},\dots ,a_{\mathrm{i\sigma }+w-1}{)}^{\mathrm{\top }}\in {\mathbb{F}}_b^w$ denote a state vector at step i (${}^{\mathrm{\top }}$ means ‘transposed’). We can define a state-space representation ${\tilde{\mathbf{x}}}_{i+1}=\tilde{\mathbf{B}}{\tilde{\mathbf{x}}}_i$, where (11) $\tilde{\mathbf{B}}=\left(\begin{array}{ccccccc}{\tilde{\mathbf{b}}}_0& {\tilde{\mathbf{b}}}_1& \cdots & {\tilde{\mathbf{b}}}_{m-1}& \mathbf{0}& \cdots & \mathbf{0}\end{array}\right)$(11) is a $w\times w$ state transition matrix in ${\mathbb{F}}_b$ that consists of m column vectors ${\tilde{\mathbf{b}}}_0,{\tilde{\mathbf{b}}}_1,\dots ,$ ${\tilde{\mathbf{b}}}_{m-1}\in {\mathbb{F}}_b^w$ and w−m zero column vectors $\mathbf{0}\in {\mathbb{F}}_b^w$. We now set b = 4 and decompose $a_i\in {\mathbb{F}}_4$, $i=0,1,\dots$, into $a_i={\overline{a}}_i\alpha +{\underset{\_}{a}}_i$ with ${\overline{a}}_i,{\underset{\_}{a}}_i\in {\mathbb{F}}_2$ since a set $\{1,\alpha \}$ is a basis of ${\mathbb{F}}_4$ over ${\mathbb{F}}_2$. Then, we can write ${\tilde{\mathbf{x}}}_{i+1}={\overline{a}}_{\mathrm{i\sigma }}\left(\alpha {\tilde{\mathbf{b}}}_0\right)+{\underset{\_}{a}}_{\mathrm{i\sigma }}{\tilde{\mathbf{b}}}_0+{\overline{a}}_{\mathrm{i\sigma }+1}\left(\alpha {\tilde{\mathbf{b}}}_1\right)+{\underset{\_}{a}}_{\mathrm{i\sigma }+1}{\tilde{\mathbf{b}}}_1+\cdots +{\overline{a}}_{\mathrm{i\sigma }+m-1}\left(\alpha {\tilde{\mathbf{b}}}_{m-1}\right)+{\underset{\_}{a}}_{\mathrm{i\sigma }+m-1}{\tilde{\mathbf{b}}}_{m-1}$, that is, a linear combination of $2m$ column vectors $\alpha {\tilde{\mathbf{b}}}_j,{\tilde{\mathbf{b}}}_j\in {\mathbb{F}}_4^w$, $j=0,1,\dots ,m-1$, with coefficients in ${\mathbb{F}}_2$. From this, we can calculate ${\tilde{\mathbf{x}}}_{i+1}$ by only adding vectors $\alpha {\tilde{\mathbf{b}}}_j$ if ${\overline{a}}_{\mathrm{i\sigma }+j}=1$ and ${\tilde{\mathbf{b}}}_j$ if ${\underset{\_}{a}}_{\mathrm{i\sigma }+j}=1$ for each j. Moreover, the elements $0,1,\alpha ,\alpha +1\in {\mathbb{F}}_4$ can be represented as column vectors $(00{)}^{\mathrm{\top }},(01{)}^{\mathrm{\top }},(10{)}^{\mathrm{\top }},(11{)}^{\mathrm{\top }}\in {\mathbb{F}}_2^2$, respectively, and hence, $\alpha {\tilde{\mathbf{b}}}_j,{\tilde{\mathbf{b}}}_j,{\tilde{\mathbf{x}}}_{i+1}\in {\mathbb{F}}_4^w$ can be viewed as column vectors in ${\mathbb{F}}_2^{2w}$. Using this property, we can generate $\{u_i{\}}_{i=0}^{\mathrm{\infty }}$ in (3(3) $u_i:=\sum _{j=0}^{w-1}\eta \left(a_{\mathrm{i\sigma }+j}\right)b^{-j-1}\in [0,1),i=0,1,2,\dots ,$(3) ) with reasonable speed, as if we performed additions over ${\mathbb{F}}_2$. The sample code is available at https://github.com/sharase/cud-f4.

Remark 3.3

Kajiura et al. [14] proved that there exists no Tausworthe generator over ${\mathbb{F}}_2$ with both maximal periodicity and t-value zero for s = 3 if $m\ge 3$. More precisely, they proved that ${\mathbb{F}}_2$-linear generators, which are a general class of linear pseudorandom number generators over ${\mathbb{F}}_2$ including Tausworthe generators (cf. [8,11]), have the t-value zero for $s=3$ only if the period length is exactly three. Their proof was specialized for the case ${\mathbb{F}}_2$; for example, they used the property ${\mathcal{L}}_m\cap {\mathcal{U}}_m=\left\{{\mathbf{I}}_m\right\}$ in [14, Proof of Theorem 1], where ${\mathbf{I}}_m$ denotes the identity matrix of order m, and ${\mathcal{L}}_m$ and ${\mathcal{U}}_m$ denote a set of non-singular $m\times m$ lower-triangular and upper-triangular matrices, respectively. This is false in the fields ${\mathbb{F}}_b$ except for ${\mathbb{F}}_2$. Indeed, Harase [9] obtained Tausworthe generators over ${\mathbb{F}}_2$ with t-value two or three for s = 3, but they were not optimal with respect to the t-value. Thus, we conducted a search over ${\mathbb{F}}_b$, whose restrictions are looser than those over ${\mathbb{F}}_2$.

Remark 3.4

We consider a reason why there are a very few pairs of polynomials $\left(p\right(x),q(x\left)\right)$ with the t-value zero for s = 3 over ${\mathbb{F}}_3$ in Table 1. Assume that $m\le w$. Let ${\mathbf{A}}_0$ denote the $(m\times m)$-transpose companion matrix of $p\left(x\right)$ in ${\mathbb{F}}_b$ given by ${\mathbf{A}}_0=\left(\begin{array}{cccc}& 1& & \\ & & \ddots & \\ & & & 1\\ c_m& c_{m-1}& \cdots & c_1\end{array}\right),$where blank entries in this matrix mean zeros. We set an $(m\times m)$-matrix $\mathbf{A}={\mathbf{A}}_0^{\sigma }$. According to [8, Section 5.1], we can obtain the state transition matrix $\tilde{\mathbf{B}}$ in (11(11) $\tilde{\mathbf{B}}=\left(\begin{array}{ccccccc}{\tilde{\mathbf{b}}}_0& {\tilde{\mathbf{b}}}_1& \cdots & {\tilde{\mathbf{b}}}_{m-1}& \mathbf{0}& \cdots & \mathbf{0}\end{array}\right)$(11) ) by expanding $\mathbf{A}$, that is, if m = w, then we put $\tilde{\mathbf{B}}=\mathbf{A}$, and if m<w, for $j=m+1,\dots ,w$, we attach $(c_1^{\left(j\right)},\dots ,c_m^{\left(j\right)})$ as the jth row vector, where the coefficients $c_i^{\left(j\right)}$ are given by the relation $a_{\mathrm{i\sigma }+j}=c_1^{\left(j\right)}{a}_{i-1}+\cdots +c_m^{\left(j\right)}{a}_{i-m}$, and add w−m columns of the zero vector $\mathbf{0}$. Thus, the t-value for dimension s is determined by the the maximum number of linear independence of leading row vectors of s generating matrices (${\mathbf{I}}_m$, $\mathbf{A},\dots ,{\mathbf{A}}^{s-1})$; see [13, Theorem 4.28] or [15, Theorem 4.52] for details. In our construction scheme, one can only change the parameter values $c_1,\dots ,c_m$ and σ, so that the search space is restricted.

As an alternative, we have conducted a numerical experiment for which we discard the structure of Tausworthe generators and take general $(m\times m)$-full rank matrices $\mathbf{A}$, not given by ${\mathbf{A}}_0^{\sigma }$, as described in [14, Equ. (4)]. (Here, we may assume without loss of generality that the row vectors $(c_1^{\left(j\right)},\dots ,c_m^{\left(j\right)})$ are arbitrary.) Our goal here is to find a full rank matrix $\mathbf{A}$ such that a digital net generated by $(\mathbf{I},\mathbf{A},{\mathbf{A}}^2)$ has the t-value zero for s = 3 and the multiplicative order of $\mathbf{A}$ is $b^m-1$. For this, we generate full rank matrices $\mathbf{A}$ at random and check the above conditions. In computer search, we have confirmed the existence of such $\mathbf{A}$ in ${\mathbb{F}}_3$ for $2\le m\le 10$. It might be expected that the existence holds true for every m in arbitrary ${\mathbb{F}}_b$ except for ${\mathbb{F}}_2$. However, this approach seems to be significantly inefficient and time-consuming because it is not so easy to find matrices $\mathbf{A}$ that generate the digital nets with the t-value zero even for s = 2 if m is large for small b. Therefore, it would be desirable to design some mathematical structure of $\mathbf{A}$ in advance before we conduct a search. In contrast, our algorithm always ensures the t-value zero for s = 2. In this paper, we conduct an exhaustive search, but our algorithm has the advantage that we can easily switch from an exhaustive search to a random search by generating $A_k\left(x\right)$ in (10(10) $\begin{array}{rl}& A_k\left(x\right)=\text{β}x+\gamma ,\end{array}$(10) ) randomly.

4. Numerical examples

We provide numerical examples to confirm the performance of Markov chain QMC. In our examples, we estimate the expectation $E_{\pi }\left[f\right(\mathbf{X}\left)\right]$ and compare the following driving sequences:

1. New: our new Tausworthe generators over ${\mathbb{F}}_4$;

2. Harase: Tausworthe generators over ${\mathbb{F}}_2$ developed by Harase [9];

3. Chen: Tausworthe generators over ${\mathbb{F}}_2$ developed by Chen et al. [7]; and

4. IID: Mersenne Twister [31].

We briefly explain how to use Tausworthe generators over ${\mathbb{F}}_b$. Recall that $N=b^m$ and the period length is N−1. For the output values (3(3) $u_i:=\sum _{j=0}^{w-1}\eta \left(a_{\mathrm{i\sigma }+j}\right)b^{-j-1}\in [0,1),i=0,1,2,\dots ,$(3) ) generated by Tausworthe generators, if $\gcd (s,N-1)=1$, we simply define s-dimensional non-overlapping points starting from the origin: (12) $(0,\dots ,0),(u_0,\dots ,u_{s-1}),(u_s,\dots ,u_{2s-1}),\dots ,(u_{(N-2)s},\dots ,u_{(N-1)s-1}).$(12) If $\gcd (s,N-1)=d>1$, instead of (12(12) $(0,\dots ,0),(u_0,\dots ,u_{s-1}),(u_s,\dots ,u_{2s-1}),\dots ,(u_{(N-2)s},\dots ,u_{(N-1)s-1}).$(12) ), we generate d distinct short loops of s-dimensional points, that is, (13) $(u_j,\dots ,u_{j+s-1}),(u_{j+s},\dots ,u_{j+2s-1}),\dots ,(u_{j+\left(\right((N-1)/d)-1)s},\dots ,u_{j+\left(\right(N-1\left)/d\right)s-1}),$(13) for $j=0,\dots ,d-1$, and concatenate them starting from the origin $(0,\dots ,0)$ in this order. For these points, we apply b-adic digital shifts, that is, we add $(z_1,\dots ,z_s)$ to each s-dimensional point using the digit-wise addition ${\oplus }_b$ (see Remark 4.1), where $z_1,\dots ,z_s$ are IID samples from $\mathcal{U}(0,1)$, that is, the continuous uniform distribution over $(0,1)$. We use the resulting points as input for Markov chain QMC; see Remark 4.3 and [2,5,6,32] for more details. We set w = 32 over ${\mathbb{F}}_2$ and w = 16 over ${\mathbb{F}}_4$ as a digit number in Definition 2.4.

Remark 4.1

We recall the definition of digital shifts. For $x=\sum _{j=0}^{\mathrm{\infty }}{\xi }_j{b}^{-j-1}\in [0,1)$ and $z=\sum _{j=0}^{\mathrm{\infty }}{\zeta }_j{b}^{-j-1}\in (0,1)$ with ${\xi }_j,{\zeta }_j\in {\mathbb{Z}}_b$, we define the b-adic digitally shifted point $\tilde{x}\in (0,1)$ as $\tilde{x}=x{\oplus }_{b}z:=\sum _{j=0}^{\mathrm{\infty }}{\psi }_j{b}^{-j-1}$, where ${\psi }_j:=\eta \left({\eta }^{-1}\right({\xi }_j)+{\eta }^{-1}({\zeta }_j\left)\right)$ with ${\psi }_j\ne b-1$ for infinitely many j and ‘+’ represents the addition in ${\mathbb{F}}_b$. For higher dimensions s>1, let $\mathbf{z}=(z_1,\dots ,z_s)\in (0,1{)}^s$. For $\mathbf{x}=(x_1,\dots ,x_s)$, we similarly define the b-adic digitally shifted point $\tilde{\mathbf{x}}\in (0,1{)}^s$ as $\tilde{\mathbf{x}}=\mathbf{x}{\oplus }_b\mathbf{z}:=(x_1{\oplus }_b{z}_1,\dots ,x_s{\oplus }_b{z}_s)$.

4.1. Gaussian Gibbs sampling

Our first example is a systematic Gibbs sampling scheme to generate the s-dimensional multivariate Gaussian (normal) distribution $\mathbf{X}=(X_1\cdots X_s{)}^{\mathrm{\top }}\sim \mathcal{N}(\mathit{\mu },\mathbf{\Sigma })$ for a mean vector $\mathit{\mu }=({\mu }_1\cdots {\mu }_s{)}^{\mathrm{\top }}$ and covariance matrix $\mathbf{\Sigma }=\left({\sigma }_{\mathrm{ij}}\right)$. This can be implemented as (14) $X_k|{\mathbf{X}}_{-k}\sim \mathcal{N}({\mathit{\mu }}_k+{\mathbf{\Sigma }}_{k,-k}{\mathbf{\Sigma }}_{-k,-k}^{-1}({\mathbf{X}}_{-k}-{\mathit{\mu }}_{-k}),{\mathbf{\Sigma }}_{k,k}-{\mathbf{\Sigma }}_{k,-k}{\mathbf{\Sigma }}_{-k,-k}^{-1}{\mathbf{\Sigma }}_{-k,k}),$(14) for $k=1,\dots ,s$, which reduces to the iteration of the calculation of the one-dimensional normal distribution. (Here, for simplicity of notation, the indices ${}_k$ and ${}_{-k}$ represent the kth component and the components except for the kth component, respectively; e.g.${\mathbf{\Sigma }}_{k,-k}=({\sigma }_{k,1},\dots ,{\sigma }_{k,k-1},{\sigma }_{k,k+1},\dots ,{\sigma }_{k,s})$ and so on.) Thus, we apply the inverse transform method in (14(14) $X_k|{\mathbf{X}}_{-k}\sim \mathcal{N}({\mathit{\mu }}_k+{\mathbf{\Sigma }}_{k,-k}{\mathbf{\Sigma }}_{-k,-k}^{-1}({\mathbf{X}}_{-k}-{\mathit{\mu }}_{-k}),{\mathbf{\Sigma }}_{k,k}-{\mathbf{\Sigma }}_{k,-k}{\mathbf{\Sigma }}_{-k,-k}^{-1}{\mathbf{\Sigma }}_{-k,k}),$(14) ). We set the parameter values s = 3 and $\mathit{\mu }=\left(\begin{array}{c}0\\ 0\\ 0\end{array}\right),\mathbf{\Sigma }=\left(\begin{array}{ccc}1& 0.3& -0.2\\ 0.3& 1& 0.5\\ -0.2& 0.5& 1\end{array}\right),$which were used in [6, Chapter 6.1].

First, we estimate $E\left[X_1\right],E\left[X_2\right]$, and $E\left[X_3\right]$ with true value 0 by taking the sample mean. Figure 3 shows a summary of the root-mean-square errors (RMSEs) in ${\log }_2$ scale for sample sizes N from $2^{10}$ to $2^{20}$ using 300 digital shifts. In all cases, the Tausworthe generators (labeled ‘New’ and ‘Harase’) optimized in terms of the t-value have almost the same accuracy and outperform Chen's generators.

Figure 3. RMSEs for $E\left[X_1\right]$, $E\left[X_2\right]$, and $E\left[X_3\right]$ with true value 0.

Furthermore, we estimate the second-order moments $E\left[X_1{X}_2\right]$, $E\left[X_1{X}_3\right]$, $E\left[X_2{X}_3\right]$ and the third-order moment $E\left[X_1{X}_2{X}_3\right]$ using 300 digital shifts, respectively. Figures 4 and 5 show summaries of the RMSEs. In Figure 5, we observe that Chen's generators are unstable and have several bumps when we estimate $E\left[X_1{X}_2{X}_3\right]$.

Figure 4. RMSEs for $E\left[X_1{X}_2\right]$, $E\left[X_1{X}_3\right]$, and $E\left[X_2{X}_3\right]$ with true values 0.3, −0.2, and 0.5.

Figure 5. RMSEs for $E\left[X_1{X}_2{X}_3\right]$ with true value 0.

4.2. M/M/1 queuing system

Our second example is an M/M/1 queuing model, which has the same setting as that in [32, Chapter 8.3.2]. Consider a single-server queuing model, where the customers arrive as a Poisson process with intensity $\lambda >0$ and the service time is exponentially distributed with intensity $\mu >0$. Assume that $\mu >\lambda$ for system stability. Let $W_j$ denote the waiting time of the jth customer, $S_j$ denote the service time of the jth customer, and $T_j$ denote the time interval between the jth customer and the $(j-1)$th customer. Then, we have the Lindley recurrence: (15) $\begin{array}{rl}& W_0=0,W_j=max(W_{j-1}+S_{j-1}-T_j,0),\end{array}$(15) (16) $\begin{array}{rl}& S_{j-1}\sim \mathcal{E}\mathit{xp}\mathit{\left(}\mu \mathit{\right)}\mathit{,}{\mathit{T}}_{\mathit{j}}\sim \mathcal{E}\mathit{xp}\mathit{\left(}\lambda \mathit{\right)}\mathit{,}\end{array}$(16) for $j\ge 1$, where $\mathcal{E}\mathit{xp}\mathit{(}\cdot \mathit{)}$ denotes the exponential distribution. Under stationarity, the average waiting time is known as (17) $E\left[W_j\right]=\frac{\lambda }{\mu (\mu -\lambda )}$(17) (cf. [33]). We estimate the average waiting time (17(17) $E\left[W_j\right]=\frac{\lambda }{\mu (\mu -\lambda )}$(17) ) by taking the sample mean $(W_1+\cdots +W_N)/N$ via Equations (15(15) $\begin{array}{rl}& W_0=0,W_j=max(W_{j-1}+S_{j-1}-T_j,0),\end{array}$(15) )–(16(16) $\begin{array}{rl}& S_{j-1}\sim \mathcal{E}\mathit{xp}\mathit{\left(}\mu \mathit{\right)}\mathit{,}{\mathit{T}}_{\mathit{j}}\sim \mathcal{E}\mathit{xp}\mathit{\left(}\lambda \mathit{\right)}\mathit{,}\end{array}$(16) ). Note that we need $2N$ random points $(S_0,T_1),\dots ,(S_{N-1},T_N)$ for N customers. Note also that the function $max(\cdot ,0)$ is unsmooth at 0.

We set the parameters $\lambda =0.5$ and $\mu =1$. Figure 6 shows the RMSEs for the average waiting time using 300 digital shifts. The three types of QMC points have almost the same performance, except for Chen's generator at $N=2^{13}$. The bump of Chen's generator coincides with that in [32, Figure 8.2].

Figure 6. RMSEs for the average waiting time of the M/M/1 queuing model.

4.3. A linear regression model

In the third example, we consider a linear regression model $y_i={\mathbf{x}}_i^{\mathrm{\top }}\text{β}+{ϵ}_i,{ϵ}_i\stackrel{\mathrm{IID}}{\sim }\mathcal{N}(0,{\tau }^2)(i=1,\dots ,n),$where $y_i$ is the ith observation on the response variable, ${\mathbf{x}}_i=(1,x_{i,1},\dots ,x_{i,k}{)}^{\mathrm{\top }}$ is a $(k+1)\times 1$ vector of 1 and ith observations on the k explanatory variables, $\text{β}=({\text{β}}_0,{\text{β}}_1,\dots ,{\text{β}}_k{)}^{\mathrm{\top }}$ is a $(k+1)\times 1$ vector of regression coefficients, and the error term ${ϵ}_i$ is IID normal with mean zero and common variance ${\tau }^2$. Let $\mathbf{X}=({\mathbf{x}}_1,\dots ,{\mathbf{x}}_n{)}^{\mathrm{\top }}$ be an $n\times (k+1)$ design matrix (with rank k + 1<n) and $\mathbf{y}=(y_1,\dots ,y_n{)}^{\mathrm{\top }}$ an $n\times 1$ vector.

We now consider Bayesian inference as follows: We assume that the parameters $\text{β}$ and ${\tau }^2$ are independent and have the prior distributions $\begin{array}{rl}\text{β}& \sim \mathcal{N}({\mathbf{b}}_0,{\mathbf{B}}_0),\\ {\tau }^2& \sim \mathcal{IG}(\frac{n_0}{2},\frac{s_0}{2}),\end{array}$where $\mathcal{IG}({\alpha }_1,{\alpha }_2)$ denotes the inverse gamma distribution with shape parameter ${\alpha }_1>0$ and rate parameter ${\alpha }_2>0$, and $(k+1)$-dimensional mean vector ${\mathbf{b}}_0$, $\left(\right(k+1)\times (k+1\left)\right)$-covariance matrix ${\mathbf{B}}_0$, $n_0$, and $s_0$ are hyperparameters. Then, according to [34,35], sampling from the joint posterior distribution of $({\text{β}}^{\mathrm{\top }},{\tau }^2{)}^{\mathrm{\top }}$ can be generated through sampling from the full conditional distributions (18) $\begin{array}{rl}\text{β}|{\tau }^2,\mathbf{y}& \sim \mathcal{N}({\mathbf{b}}_1,{\mathbf{B}}_1),\end{array}$(18) (19) $\begin{array}{rl}{\tau }^2|\text{β},\mathbf{y}& \sim \mathcal{IG}(\frac{n_1}{2},\frac{s_1}{2}),\end{array}$(19) where $\begin{array}{rl}{\mathbf{b}}_1& ={\mathbf{B}}_1({\mathbf{B}}_0^{-1}{\mathbf{b}}_0+{\tau }^{-2}{\mathbf{X}}^{\mathrm{\top }}\mathbf{y}),{\mathbf{B}}_1^{-1}={\mathbf{B}}_0^{-1}+{\tau }^{-2}{\mathbf{X}}^{\mathrm{\top }}\mathbf{X},\\ n_1& =n_0+n,s_1=s_0+(\mathbf{y}-\mathbf{X}\text{β}{)}^{\mathrm{\top }}(\mathbf{y}-\mathbf{X}\text{β}).\end{array}$Thus, we calculate $E\left[\text{β}\right]$ and $E\left[{\tau }^2\right]$ by taking the sample mean using the Gibbs sampler based on (18(18) $\begin{array}{rl}\text{β}|{\tau }^2,\mathbf{y}& \sim \mathcal{N}({\mathbf{b}}_1,{\mathbf{B}}_1),\end{array}$(18) ) and (19(19) $\begin{array}{rl}{\tau }^2|\text{β},\mathbf{y}& \sim \mathcal{IG}(\frac{n_1}{2},\frac{s_1}{2}),\end{array}$(19) ). We generate $\text{β}$ in (18(18) $\begin{array}{rl}\text{β}|{\tau }^2,\mathbf{y}& \sim \mathcal{N}({\mathbf{b}}_1,{\mathbf{B}}_1),\end{array}$(18) ) via $\mathbf{L}\left({\mathrm{\Phi }}^{-1}\right(u_j),$ $\dots ,{\mathrm{\Phi }}^{-1}\left(u_{j+k}\right){)}^{\mathrm{\top }}$ for $u_j,\dots ,u_{j+k}\in (0,1)$, where ${\mathbf{B}}_1=\mathbf{L}{\mathbf{L}}^{\mathrm{\top }}$ is the Cholesky decomposition and $\mathrm{\Phi }(\cdot )$ is the cumulative distribution function of the standard normal distribution.

As a numerical example, we use the Boston housing data analysed in [36]. To investigate the demand for clean air, Harrison and Rubinfeld [36] built a linear regression model given by (20) $\begin{array}{rl}\log \left(\mathrm{MEDV}\right)& ={\text{β}}_0+{\text{β}}_1\mathrm{CRIM}+{\text{β}}_2\mathrm{ZN}+{\text{β}}_3\mathrm{INDUS}+{\text{β}}_4\mathrm{CHAS}+{\text{β}}_5{\mathrm{NOX}}^2\\ & +{\text{β}}_6{\mathrm{RM}}^2+{\text{β}}_7\mathrm{AGE}+{\text{β}}_8\log \left(\mathrm{DIS}\right)+{\text{β}}_9\log \left(\mathrm{RAD}\right)+{\text{β}}_{10}\mathrm{TAX}\\ & +{\text{β}}_{11}\mathrm{PTRATIO}+{\text{β}}_{12}\mathrm{B}+{\text{β}}_{13}\log \left(\mathrm{LSTAT}\right)+ϵ,ϵ\sim \mathcal{N}(0,{\tau }^2),\end{array}$(20) where the housing price MEDV is a response variable, and $\mathrm{CRIM}$, $\mathrm{ZN}$, $\mathrm{INDUS}$, $\mathrm{CHAS}$, $\mathrm{NOX}$, $\mathrm{RM}$, $\mathrm{AGE}$, $\mathrm{DIS}$, $\mathrm{RAD}$, $\mathrm{TAX}$, $\mathrm{PTRATIO}$, $\mathrm{B}$, and $\mathrm{LSTAT}$ are 13 explanatory variables (i.e. k = 13); see [36, Table IV] for more details about the variables. In our experiment, we estimate the same linear regression model as in (20(20) $\begin{array}{rl}\log \left(\mathrm{MEDV}\right)& ={\text{β}}_0+{\text{β}}_1\mathrm{CRIM}+{\text{β}}_2\mathrm{ZN}+{\text{β}}_3\mathrm{INDUS}+{\text{β}}_4\mathrm{CHAS}+{\text{β}}_5{\mathrm{NOX}}^2\\ & +{\text{β}}_6{\mathrm{RM}}^2+{\text{β}}_7\mathrm{AGE}+{\text{β}}_8\log \left(\mathrm{DIS}\right)+{\text{β}}_9\log \left(\mathrm{RAD}\right)+{\text{β}}_{10}\mathrm{TAX}\\ & +{\text{β}}_{11}\mathrm{PTRATIO}+{\text{β}}_{12}\mathrm{B}+{\text{β}}_{13}\log \left(\mathrm{LSTAT}\right)+ϵ,ϵ\sim \mathcal{N}(0,{\tau }^2),\end{array}$(20) ). Note that the state vector $({\text{β}}^{\mathrm{\top }},{\tau }^2{)}^{\mathrm{\top }}$ has 15 dimensions (i.e.s = 15).

We set the hyperparameter values ${\mathbf{b}}_0=\mathbf{0}$, ${\mathbf{B}}_0=100{\mathbf{I}}_s$, $n_0=5$, $s_0=0.01$, and run the Gibbs sampler for 5000 iterations using random numbers as a burn-in period. Then, we calculate $E\left[\text{β}\right]$ and $E\left[{\tau }^2\right]$ by running the Gibbs sampler for $N=2^{12}$, $2^{14}$, $2^{16}$, and $2^{18}$ iterations. Table 4 shows a summary of sample variances of posterior mean estimates using 300 digital shifts. In both cases ${\mathbb{F}}_2$ and ${\mathbb{F}}_4$, Tausworthe generators optimized in terms of the t-value provide comparable to or better results than Chen's Tausworthe generators optimized in terms of the equidistribution property, excluding some exceptions (e.g. ${\text{β}}_6,\dots ,{\text{β}}_{13}$ for $N=2^{14}$ estimated by Tausworthe generators over ${\mathbb{F}}_2$). Furthermore, we plot the histograms of ${\text{β}}_0$ and ${\text{β}}_8$ using $2^{14}$ IID uniform random points and QMC points generated by our new generator over ${\mathbb{F}}_4$ in Figure 7. In the case ${\text{β}}_0$, the sampling using our new QMC points tends to converge to the posterior distribution faster than the sampling using IID uniform random points, but in the case ${\text{β}}_8$, the difference seems to be unclear. From this, the estimation of ${\text{β}}_8$ might be more difficult than that of ${\text{β}}_0$ when we apply QMC. Overall, our experiment implies that the t-value is a good measure of uniformity in the study of Markov chain QMC.

Figure 7. Histograms of ${\text{β}}_0$ and ${\text{β}}_8$ using $2^{14}$ IID uniform random points and QMC points generated by our new generator over ${\mathbb{F}}_4$.

Table 4. Variances of posterior mean estimates for $\text{β}=({\text{β}}_0,\dots ,{\text{β}}_{13}{)}^{\mathrm{\top }}$ and ${\tau }^2$.

$N=2^{12}$
Parameter	${\text{β}}_0$	${\text{β}}_1$	${\text{β}}_2$	${\text{β}}_3$	${\text{β}}_4$	${\text{β}}_5$	${\text{β}}_6$
IID	6.51e−06	3.56e−10	6.31e−11	1.26e−09	2.53e−07	3.28e−06	3.79e−10
Chen	1.15e−09	6.29e−14	1.03e−14	2.52e−13	4.55e−11	5.38e−10	8.95e−14
Harase	2.98e−10	1.41e−14	2.77e−15	5.60e−14	1.18e−11	1.27e−10	1.72e−14
New	5.35e−10	3.59e−14	5.40e−15	1.40e−13	2.99e−11	3.85e−10	7.83e−14
Parameter	${\text{β}}_7$	${\text{β}}_8$	${\text{β}}_9$	${\text{β}}_{10}$	${\text{β}}_{11}$	${\text{β}}_{12}$	${\text{β}}_{13}$	${\tau }^2$
IID	7.17e−11	2.92e−07	9.45e−08	3.60e−12	5.69e−09	3.07e−12	1.48e−07	1.22e−09
Chen	1.37e−14	5.02e−11	1.60e−11	7.22e−16	1.18e−12	4.97e−16	3.19e−11	2.59e−13
Harase	2.03e−15	1.18e−11	3.33e−12	1.59e−16	2.05e−13	1.16e−16	6.78e−12	1.27e−13
New	7.26e−15	2.89e−11	1.06e−11	4.11e−16	6.26e−13	3.09e−16	2.09e−11	1.34e−13
$N=2^{14}$
Parameter	${\text{β}}_0$	${\text{β}}_1$	${\text{β}}_2$	${\text{β}}_3$	${\text{β}}_4$	${\text{β}}_5$	${\text{β}}_6$
IID	1.34e−06	1.17e−10	1.57e−11	3.13e−10	7.18e−08	7.61e−07	1.12e−10
Chen	1.24e−10	9.66e−15	1.21e−15	3.27e−14	2.14e−11	6.62e−11	8.80e−15
Harase	2.23e−11	1.75e−15	2.40e−16	5.19e−15	1.07e−12	1.26e−11	8.91e−15
New	1.31e−11	6.83e−16	1.12e−16	4.55e−15	9.99e−13	1.09e−11	1.09e−15
Parameter	${\text{β}}_7$	${\text{β}}_8$	${\text{β}}_9$	${\text{β}}_{10}$	${\text{β}}_{11}$	${\text{β}}_{12}$	${\text{β}}_{13}$	${\tau }^2$
IID	1.83e−11	6.36e−08	2.36e−08	9.15e−13	1.66e−09	6.75e−13	4.37e−08	3.00e−10
Chen	1.46e−15	5.42e−12	1.81e−12	7.96e−17	1.29e−13	5.71e−17	3.32e−12	2.88e−14
Harase	2.46e−15	6.10e−12	2.74e−12	9.44e−17	1.27e−13	7.52e−17	4.07e−12	2.34e−14
New	1.29e−16	5.94e−13	1.69e−13	8.36e−18	1.26e−14	6.15e−18	3.43e−13	5.58e−15
$N=2^{16}$
Parameter	${\text{β}}_0$	${\text{β}}_1$	${\text{β}}_2$	${\text{β}}_3$	${\text{β}}_4$	${\text{β}}_5$	${\text{β}}_6$
IID	3.07e−07	2.74e−11	4.13e−12	8.43e−11	1.64e−08	1.72e−07	2.22e−11
Chen	9.28e−12	6.59e−16	1.07e−16	2.37e−15	4.71e−13	4.73e−12	6.69e−16
Harase	1.40e−12	9.12e−17	1.82e−17	3.55e−16	1.71e−13	7.11e−13	7.65e−17
New	1.73e−12	9.86e−17	1.71e−17	4.20e−16	7.31e−14	9.01e−13	1.24e−16
Parameter	${\text{β}}_7$	${\text{β}}_8$	${\text{β}}_9$	${\text{β}}_{10}$	${\text{β}}_{11}$	${\text{β}}_{12}$	${\text{β}}_{13}$	${\tau }^2$
IID	4.26e−12	1.54e−08	5.87e−09	2.25e−13	3.48e−10	1.89e−13	1.09e−08	7.35e−11
Chen	1.14e−16	4.70e−13	1.67e−13	6.91e−18	1.08e−14	4.52e−18	2.34e−13	2.36e−15
Harase	8.86e−18	4.47e−14	1.42e−14	1.18e−18	1.36e−15	6.42e−19	3.02e−14	9.24e−16
New	1.63e−17	7.71e−14	3.20e−14	1.36e−18	2.28e−15	1.20e−18	4.68e−14	1.00e−15
$N=2^{18}$
Parameter	${\text{β}}_0$	${\text{β}}_1$	${\text{β}}_2$	${\text{β}}_3$	${\text{β}}_4$	${\text{β}}_5$	${\text{β}}_6$
IID	9.16e−08	5.65e−12	9.53e−13	1.85e−11	4.39e−09	5.28e−08	6.81e−12
Chen	5.65e−13	3.33e−17	5.70e−18	1.36e−16	2.55e−14	2.73e−13	4.09e−17
Harase	4.62e−14	2.48e−18	4.95e−19	1.07e−17	2.44e−15	2.40e−14	4.26e−18
New	7.02e−14	5.55e−18	9.27e−19	2.10e−17	3.93e−15	4.46e−14	5.85e−18
Parameter	${\text{β}}_7$	${\text{β}}_8$	${\text{β}}_9$	${\text{β}}_{10}$	${\text{β}}_{11}$	${\text{β}}_{12}$	${\text{β}}_{13}$	${\tau }^2$
IID	9.29e−13	3.75e−09	1.48e−09	5.92e−14	9.52e−11	3.47e−14	2.59e−09	1.52e−11
Chen	6.90e−18	2.74e−14	9.00e−15	3.51e−19	5.77e−16	2.69e−19	1.56e−14	2.31e−16
Harase	6.80e−19	2.49e−15	8.37e−16	3.47e−20	5.11e−17	2.51e−20	1.48e−15	3.08e−17
New	8.56e−19	4.08e−15	1.19e−15	5.45e−20	7.38e−17	3.27e−20	1.89e−15	2.56e−17

Remark 4.2

According to some heuristic arguments in [6, Chapter 7.1], it is expected that Markov chain QMC drastically improves the rate of convergence when the dependence of states on the past decays quickly. To investigate such phenomena, we plot the sample paths and autocorrelation functions (ACFs) of ${\text{β}}_0,{\text{β}}_1,{\text{β}}_2$, and ${\tau }^2$ in Figure 8 using $2^{14}$ IID uniform random points, after a burn-in period with 5000 iterations. The ACF plots imply that the dependence of states on the past decays very quickly (i.e. at one step) and has negligible effect on the current state. Moreover, it is believed that QMC methods in high-dimensional problems are successful especially in the case where the problems are dominated by the first leading variables or well approximated by a sum of functions of at most one or two variables (cf. [37]). Our linear regression example is probably included in such a class of problems, and hence, all the three Tausworthe generators drastically improve the rate of convergence. On the other hand, in more complicated applications in practice, the difference among these generators might not become clear, but we expect that Tausworthe generators optimized in terms of the t-value would be at worst superior to IID uniform random points.

Figure 8. Sample paths and ACFs of ${\text{β}}_0,{\text{β}}_1,{\text{β}}_2$, and ${\tau }^2$ using $2^{14}$ IID uniform random points.

Remark 4.3

The generation scheme in (13(13) $(u_j,\dots ,u_{j+s-1}),(u_{j+s},\dots ,u_{j+2s-1}),\dots ,(u_{j+\left(\right((N-1)/d)-1)s},\dots ,u_{j+\left(\right(N-1\left)/d\right)s-1}),$(13) ) was originally used in [2,5]. However, if $\gcd (s,N-1)=d>1$, then we have d−1 skips in (13(13) $(u_j,\dots ,u_{j+s-1}),(u_{j+s},\dots ,u_{j+2s-1}),\dots ,(u_{j+\left(\right((N-1)/d)-1)s},\dots ,u_{j+\left(\right(N-1\left)/d\right)s-1}),$(13) ) through the entire period. To avoid such irregular skips, Tribble [6] and Chen [32] suggested a strategy in which the skips are the same between every pair of non-overlapping s-blocks: Let r be the smallest integer $r\ge s$ such that $gcd(r,N-1)=1$. Then, instead of (13(13) $(u_j,\dots ,u_{j+s-1}),(u_{j+s},\dots ,u_{j+2s-1}),\dots ,(u_{j+\left(\right((N-1)/d)-1)s},\dots ,u_{j+\left(\right(N-1\left)/d\right)s-1}),$(13) ), we can consider s-dimensional non-overlapping points starting from the origin: (21) $\begin{array}{rl}(0,\dots ,0),& (u_0,\dots ,u_{s-1}),(u_r,\dots ,u_{r+s-1}),\\ & (u_{2r},\dots ,u_{2r+s-1}),\dots ,(u_{(N-2)r},\dots ,u_{(N-2)r+s-1}),\end{array}$(21) which maintain balance (i.e. every r steps) in each coordinate by discarding r−s points between each block. We also implemented the strategy (21(21) $\begin{array}{rl}(0,\dots ,0),& (u_0,\dots ,u_{s-1}),(u_r,\dots ,u_{r+s-1}),\\ & (u_{2r},\dots ,u_{2r+s-1}),\dots ,(u_{(N-2)r},\dots ,u_{(N-2)r+s-1}),\end{array}$(21) ) and conducted the same experiments as in Section 4. We obtained almost similar results with a slight fluctuation. In this paper, we optimized Tausworthe generators in terms of the t-value for consecutive output values, so we adopted the scheme (13(13) $(u_j,\dots ,u_{j+s-1}),(u_{j+s},\dots ,u_{j+2s-1}),\dots ,(u_{j+\left(\right((N-1)/d)-1)s},\dots ,u_{j+\left(\right(N-1\left)/d\right)s-1}),$(13) ) without discarding r−s points, which seems to be closer to the condition (1(1) $\underset{N\to \mathrm{\infty }}{lim}{D}_N^{\ast s}\left(\right(u_0,\dots ,u_{s-1}),(u_s,\dots ,u_{2s-1}),\dots ,(u_{(N-1)s},\dots ,u_{\mathrm{Ns}-1}\left)\right)=0$(1) ). We refer the reader to [6, Chapter 5] and [32, Chapter 8.2] for more details.

5. Conclusion

We attempted to search for short-period Tausworthe generators over arbitrary finite fields ${\mathbb{F}}_b$ for Markov chain QMC in terms of the t-value. To achieve this, we generalized the search algorithms [9,10] over ${\mathbb{F}}_2$ to those over ${\mathbb{F}}_b$. We conducted an exhaustive search, especially in the case where b = 3, 4, and 5, and implemented Tausworthe generators over ${\mathbb{F}}_4$ with t-values zero for dimension s = 3, in addition to s = 2, and small for $s\ge 4$. We also reported numerical examples in which both Tausworthe generators over ${\mathbb{F}}_2$ and ${\mathbb{F}}_4$ optimized in terms of the t-value perform comparable to or better than Tausworthe generators [7] optimized in terms of the equidistribution property.

The two-element field ${\mathbb{F}}_2$ is the most important finite field in applications, but has some restrictions that do not occur over other fields ${\mathbb{F}}_b$. Therefore, in future work, it would be interesting to study implementations of other types of QMC points over ${\mathbb{F}}_b$, such as polynomial lattice rules [13,15,30] and irreducible Sobol'–Niderreiter sequences [38], which are closely related to the t-value. Furthermore, we are also planning to apply our new generators, including [9], to a large variety of Bayesian computation using real-life data.

To conclude this paper, we mention some recent related works. In past a decade, the application of QMC methods to computational statistics has received a lot of attention for researchers and many novel studies have been proposed. For example, Gerber and Chopin [39] present a class of algorithms where a sequential Monte Carlo strategy is implemented with QMC. Buchholz and Chopin [40] derive approximate Bayesian computation (ABC) algorithms based on QMC for dealing with models with an intractable likelihood. As another direction, research on kernel density estimation using QMC has been actively conducted. We refer the reader to the survey paper [41] for recent progress in this topic.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by JSPS KAKENHI Grant Numbers JP22K11945, JP18K18016.