## 2 p-Refined Multilevel Quasi-Monte Carlo
### Algorithm
1\(L=0\), B \(\leftarrow\infty\), \(Var\leftarrow\infty\) ;
2 while\(\frac{\epsilon}{\sqrt{2}}<B\)do
3 if\(L>2\)then
4 while\(\frac{\epsilon}{\sqrt{2}}<Var\)do
5 Find the level \(\ell\) with largest \(\mathcal{V}_{\ell}/\text{C}_{\ell}\) among all \(\ell\in\{0,1,...,L\}\) ;
6\(N_{\ell}\gets N_{\ell}*1.2\) ;
7 Take \(N_{\ell}\) samples on level \(\ell\) for each random shift \(\Xi_{r},r=1,2,...,R_{L}\) ;
8 Calculate \(\mathcal{E}_{\ell}\) and \(\mathcal{V}_{\ell}\) ;
9\(Var\leftarrow\sum_{\ell=0}^{L}\mathcal{V}_{\ell}\) ;
10 Calculate \(B\) ;
11
12 else
13 Take \(N_{L}\) samples on level \(L\) for each random shift \(\Xi_{r},r=1,2,...,R_{L}\) ;
14 Compute \(\text{C}_{L}\), and calculate \(\mathcal{E}_{L}\) and \(\mathcal{V}_{L}\) according to Equation (17) and (18);
15\(Var\leftarrow\sum_{\ell=0}^{L}\mathcal{V}_{\ell}\) ;
16\(L\gets L+1\) ;
17
18 error \(\leftarrow\sqrt{Var+B^{2}}\) ;
19\(E\leftarrow\sum_{\ell=0}^{L}\mathcal{E}_{\ell}\) ;
20 return \(E\), error
```  
**Algorithm 1**p-MLQMC. We will now describe an algorithm which generates the location of the points in which the discrete values of the random field are to be computed. Figures 5 and 6, show the location of these points \(\bullet\). As input to this algorithm, we require a set of Gauss quadrature points in natural coordinates  
Figure 5: Quadrature points \(\square\) and points where the discrete values of the random field are to be generated on one element \(\bullet\) for model problem 1. Figure 6: Quadrature points \(\triangle\) and points where the discrete values of the random field are to be generated on one element \(\bullet\) for model problem 2.  
for one element of the finest considered mesh. An algorithm for generating the point sets on each level is given in Algorithm 2. The output consists of a vector of points in which the discrete values of the random field are to be computed. This vector contains all the locations of the points on the finest level. Because the vector is sorted, taking subsets hereof gives the locations of the points on the coarser levels. Figure 7 shows in which sequence the points are chosen on a unit square and on a unit triangle. For example, for model problem 2, the evaluation points at the coarsest level are points \(\{1,2,3,4,5,6,7\}\), the second coarsest level points \(\{1,...,7,8,9,10,11,12,13,14,15,16\}\), the third coarsest level points \(\{1,...,16,17,18,19,20,21,22,23,24,25,26,27,28\}\), and at the finest level points \(\{1,...,28,29,30,31,32,33,34,35,36,37\}\). See also Tables 1 and 2 which, among other, list the number of quadrature points used per level. Once the locations of the points for one element in natural coordinates for each level have been obtained, the points are mapped to each element of the mesh by means of a transformation matrix. This results in a point set, listing the locations in global coordinates where the discrete values of the random field for the model problem are to be computed. As such, Algorithm 2 only works for two dimensional problems. Algorithm 2 consists of two main loops, the first loop starts on line 3, and returns a set of points which are closest to points of a given quadrature set on level \(\ell\). The second loop, starting on line 12, reorders the obtained set of closest points as in Figure 7, following the principle illustrated in Figure 8. Extending Algorithm 2 for three dimensions can be done be rewriting lines 12 to 24 to accommodate a reordering in three dimensions or to not consider a reordering at all. This would mean that lines 12 to 24 have to be omitted when computing the locations for the points in three dimensions. In its current form, Algorithm 1 is level-adaptive, in the sense that levels are added to the estimator (\(L\gets L+1\)) only when needed to reduce the bias. This means that the maximum level \(L\) that satisfies the bias constraint is not known a priori.