\documentclass[acmtog]{acmart}
\acmSubmissionID{126}


\usepackage{booktabs} % For formal tables

\usepackage{amsmath}
\usepackage{nicematrix}
\usepackage{tikz}
\usetikzlibrary {arrows.meta}
\usepackage[normalem]{ulem}
\usepackage{wrapfig}
\usepackage{algpseudocodex}
\usepackage{multirow}

\usetikzlibrary{decorations.pathreplacing}
\usetikzlibrary{spy} % LATEX and plain TEX

% TOG prefers author-name bib system with square brackets
\citestyle{acmauthoryear}
%\setcitestyle{nosort,square} % nosort to allow for manual chronological ordering


\usepackage{overpic}

\usepackage[ruled]{algorithm2e} % For algorithms
\renewcommand{\algorithmcfname}{ALGORITHM}

\SetAlFnt{\small}
\SetAlCapFnt{\small}
\SetAlCapNameFnt{\small}
\SetAlCapHSkip{0pt}

\DeclareMathSymbol{:}{\mathord}{operators}{"3A}
\DeclareMathOperator{\diag}{diag}
\DeclareMathOperator{\Id}{Id}


% Metadata Information
\setcopyright{cc}
\setcctype{by}
\acmJournal{TOG}
\acmYear{2025} \acmVolume{44} \acmNumber{4} \acmArticle{} \acmMonth{8} \acmPrice{}\acmDOI{10.1145/3730821}





\usepackage[many]{tcolorbox}
\tcbuselibrary{hooks, breakable}

\newtcolorbox{prop}[2][]{boxsep=0mm,bottom=1pt,top=0pt,toptitle=0mm, bottomtitle=0mm,detach title,after upper={\par\texttt{\tcbtitle}\hfill},arc=0mm,outer arc=0mm, boxrule=1pt, toprule=1pt,boxsep=0pt,colbacktitle=white,colback=white,coltitle=black,
title={#2},#1}
%% https://tex.stackexchange.com/questions/568782/allowing-paragraph-continuation-after-a-tcolorbox
\makeatletter
\tcbset{
  after app={%
    \ifx\tcb@drawcolorbox\tcb@drawcolorbox@breakable
    \else
      % add only when not breakabel
      \@endparenv
    \fi
  }
}
% for breakable
\appto\tcb@use@after@lastbox{\@endparenv\@doendpe}
\makeatother

% Document starts
\begin{document}
% Title portion
\title{Sobol' Sequences with Guaranteed-Quality 2D Projections}

% DO NOT ENTER AUTHOR INFORMATION FOR ANONYMOUS TECHNICAL PAPER SUBMISSIONS TO SIGGRAPH 2019!
\author{Nicolas Bonneel}
\orcid{0000-0001-5243-4810}
 \affiliation{
   \institution{CNRS, Université Claude Bernard Lyon 1, INSA Lyon}
     \country{France}
     }
 \email{nicolas.bonneel@liris.cnrs.fr} 
\author{David Coeurjolly}
\orcid{0000-0003-3164-8697}
\affiliation{
    \institution{CNRS, Université Claude Bernard Lyon 1, INSA Lyon}
    \country{France}
     }
 \email{david.coeurjolly@cnrs.fr}


\author{Jean-Claude Iehl}
\orcid{0000-0001-6877-2398}
 \affiliation{
     \institution{Université Claude Bernard Lyon 1, CNRS, INSA Lyon}
     \country{France}
     }
 \email{jean-claude.iehl@liris.cnrs.fr}

\author{Victor Ostromoukhov}
\orcid{0009-0004-3123-9388}
 \affiliation{
     \institution{Université Claude Bernard Lyon 1, CNRS, INSA Lyon}
     \country{France}
     }
 \email{victor.ostromoukhov@liris.cnrs.fr}



%\renewcommand\shortauthors{Bonnneel, N. et al}

\setlength{\fboxsep}{.2pt}
\newcommand{\be}{{{\small [e]}}}
\newcommand{\btwoe}{{{\small [2e]}}}
\newcommand{\bfoure}{{{\small [4e]}}}
\newcommand{\bet}{{{\small [\tilde{e}]}}}
\newcommand{\btwoet}{{{\small [2\tilde{e}]}}}
\newcommand{\bfouret}{{{\small [4\tilde{e}]}}}

\definecolor{OliveGreen}{RGB}{63, 126, 49}
\definecolor{BurntOrange}{RGB}{255, 142, 39}


\newcommand{\dummyfig}[3]{
\begin{tikzpicture}
\draw[color=blue,fill=blue!20] (0,0) rectangle (#2,#3) node[pos=.5] {#1};
\end{tikzpicture}}

\begin{teaserfigure}
    \centering
    \includegraphics[width=12cm]{figures_siggraph25/teaser_e5_10.pdf}
    \begin{overpic}[width=17cm]{figures_siggraph25/JKvsOneTwoSeq.png}
 	 \put(-2,8.5) {\rotatebox{90}{Sobol'}}
   	 \put(-2,2) {\rotatebox{90}{Ours}}
        \put(-2.2,-3){\rotatebox{90}{dims.}}
        \put(1.5,-2){$(0,1)$}
        \put(8.61538,-2){$(2,3)$}
        \put(15.7308,-2){$(4,5)$}
        \put(22.8462,-2){$(6,7)$}
       \put(29.9615,-2){$(8,9)$}        
       \put(37.0769,-2){$(10,11)$}        
       \put(44.1923,-2){$(12,13)$}        
       \put(51.3077,-2){$(14,15)$}        
      \put(58.4231,-2){$(16,17)$}        
      \put(65.5385,-2){$(18,19)$}        
      \put(72.6538,-2){$(20,21)$}        
      \put(79.7692,-2){$(22,23)$}        
      \put(86.8846,-2){$(24,25)$}        
      \put(94,-2){$(26,27)$}        
    \end{overpic}
  \vspace{0.5cm}
    \caption{We demonstrate that 2D Sobol' sequences constructed with polynomials $p$ and $p^2+p+1$ have a characteristic matrix $K = M_{p^2+p+1} M_{p}^{-1}$ that can be obtained with a simple recursive algorithm.
    This is illustrated using polynomials of degrees $e=5$ and $2e=10$, where each colored block has dimensions $e \times e$.
     They produce high-quality $(1, 2)$-sequences (with quality factor $t=1$) under mild conditions on $K$'s blocks and $p$ (bottom). The quality factor $t$ of each point set is indicated in the upper-right corner of each patch (only 256 points are shown).
    We use these $(1, 2)$-sequences to construct higher-dimensional low-discrepancy sequences with high-quality 2D and 4D projections. }
    \label{fig:teaser}
\end{teaserfigure}

\begin{abstract}
Low-discrepancy sequences, and more particularly Sobol' sequences are gold standard for drawing highly uniform samples for quasi-Monte Carlo applications. 
They produce so-called $(t,s)$-sequences, that is, sequences of $s$-dimensional samples whose uniformity is controlled by a non-negative integer quality factor $t$.
The Monte Carlo integral estimator has a convergence rate that improves as $t$ decreases. Sobol' construction in base 2 also provides extremely fast sampling point generation using efficient xor-based arithmetic.
Computer graphics applications, such as rendering, often require high uniformity in consecutive 2D projections and in higher-dimensional projections at the same time.
However, it can be shown that, in the classical Sobol' construction, only a single 2D sequence of points  (up to scrambling), constructed using irreducible polynomials $x$ and $x+1$, achieves the ideal $t=0$ property. 
Reusing this sequence in projections necessarily loses high dimensional uniformity.
We prove the existence and construct many 2D Sobol' sequences having $t=1$ using irreducible polynomials $p$ and $p^2+p+1$. 
They can be readily combined to produce higher-dimensional low discrepancy sequences with a high-quality $t=1$, guaranteed in consecutive pairs of dimensions. 
We provide the initialization table that can be directly used with any existing Sobol' implementation, along with the corresponding generator matrices, for an optimized 692-dimensional Sobol' construction. 
In addition to guaranteeing the $(1,2)$-sequence property for all consecutive pairs, we ensure that $t \leq 4$ for consecutive 4D projections up to $2^{15}$ points.


\end{abstract}



%
% The code below should be generated by the tool at
% http://dl.acm.org/ccs.cfm
% Please copy and paste the code instead of the example below.
%
\begin{CCSXML}
<ccs2012>
   <concept>
       <concept_id>10002950.10003648.10003670.10003687</concept_id>
       <concept_desc>Mathematics of computing~Random number generation</concept_desc>
       <concept_significance>500</concept_significance>
       </concept>
   <concept>
       <concept_id>10010147.10010371.10010372</concept_id>
       <concept_desc>Computing methodologies~Rendering</concept_desc>
       <concept_significance>300</concept_significance>
       </concept>
   <concept>
       <concept_id>10002950.10003714.10003715.10003718</concept_id>
       <concept_desc>Mathematics of computing~Computations in finite fields</concept_desc>
       <concept_significance>100</concept_significance>
       </concept>
   <concept>
       <concept_id>10002950.10003714.10003740</concept_id>
       <concept_desc>Mathematics of computing~Quadrature</concept_desc>
       <concept_significance>300</concept_significance>
       </concept>
 </ccs2012>
\end{CCSXML}

\ccsdesc[500]{Mathematics of computing~Random number generation}
\ccsdesc[300]{Computing methodologies~Rendering}
\ccsdesc[100]{Mathematics of computing~Computations in finite fields}
\ccsdesc[300]{Mathematics of computing~Quadrature}

%
% End generated code
%


\keywords{Quasi-Monte Carlo, Sobol', low discrepancy sequences, irreducible polynomials}




\maketitle

\section{Introduction}


Integrating a function using quasi-Monte Carlo consists in evaluating it at well-chosen uniformly distributed sample points, and averaging these values.
Sobol' sequences arise as the cornerstone of quasi-Monte Carlo, by producing extremely well-distributed sampling points whose uniformity drastically improves the convergence rate, compared to classical Monte Carlo methods.
 These points are constructed using a fast and compact recursive algorithm involving polynomials and matrices.
 Sobol' sequences have thus been widely adopted in computer graphics, notably for rendering where efficiency is paramount. Their mathematical beauty, their connections to Pascal matrices (or Sierpinsky triangle) and Galois theory also make them appealing to the mathematician, but their difficulty comprehending may discourage others. This paper gives new fundamental mathematical insights on Sobol' sequences, exploring particular pairs of polynomials and new Sobol' matrix constructions, with practical and provable benefits in terms of quasi-Monte Carlo integration error.

Sobol' sequences produce a sequence of samples in arbitrary dimension by multiplying Sobol' matrices with a base-$b$ representation of the sample index, where one Sobol' matrix is used per dimension. To compute a Sobol' matrix, a polynomial $p$ of degree $e$ and an initialization matrix of size $e \times e$ (often called \textit{direction vectors} in prior work) are required. The (triangular) Sobol' matrix is recursively constructed column by column, where the next column is computed as a linear combination of the previous $e$ columns 
weighted by the polynomial coefficients (plus a shifted column), with all operations performed modulo $b$ (or over Galois Field $GF(b)$).
The uniformity of the produced sample points is determined by a non-negative integer parameter $t$, where $t=0$ corresponds to the best achievable quality, and the quasi-Monte Carlo integral estimator converges with a rate roughly proportional to $b^t$~(see \cite{Lemieux2009MCQMCSampling} page 157, and \cite{niederreiter1988low}). 

An $s$-dimensional  set with $b^m$ samples with quality $t$ is called a $(t,m,s)$-net (see Fig.~\ref{fig:t1}). If such a point set is a $(t,m,s)$-net for all $m$, then it is a  $(t,s)$-sequence. 
We are particularly interested in the case $s=2$, not only for producing 2D points for 2-dimensional problems, but most importantly to control 2-dimensional projections of higher-dimensional problems~\cite{joe2008constructing} arising for example in computer graphics such as rendering.

It has been demonstrated that, in base $b=2$, and when considering 2D points, matrices that generate $(0,2)$-sequences are inherently related by a Pascal matrix~\cite{HoferSuzuki2019,ahmed2023analysis}. 
Consequently, the only pair of polynomials that produce $t=0$ using Sobol' construction are $x$ and $ x+1$. 
The space of $t=0$ sequences is thus extremely limited, when $b=2$.
One possible solution is to increase $b$, as suggested in~\cite{Quad2024}, which allows for additional polynomials that generate $t=0$ sequences.
However, the effect on the integration error for $t\neq 0$ becomes more significant due to the $b^t$ factor in the convergence rate. 
In addition, base $b=2$ allows for extremely fast implementations using vectorized xor-based arithmetic, which is not the case when $b>2$. For practical reasons, Sobol' polynomials and initialization matrices have thus been optimized mostly in base $b=2$~\cite{joe2008constructing}, though consecutive dimensions typically produce increasing $t$ values as dimension increases which limits their use for rendering~\cite{christensen2018renderman}.

Focusing on the case $b=2$, in this paper we show that there exist many $(1,2)$-sequences that can be constructed from pairs of irreducible polynomials $p$ and $p^2+p+1$. 
These 2D sequences can be combined to produce higher dimensional $(t,s)$-sequences of high-quality $t=1$ in consecutive 2D projections, which is the best quality achievable for $b=2$ aside from the single $t=0$ pair mentioned above. 
We also observed that, in practical integration problems, the quasi-Monte Carlo convergence rate did not differ significantly between $t=0$ and $t=1$ in base $2$ (see Fig.~\ref{fig:correlation}), making $t=1$ a compelling compromise.
Our use of a standard $b=2$ Sobol' framework makes our sequences readily usable in production renderers already using Sobol' sequences, by simply replacing existing polynomials and initializations with ours.
In our quest to prove the quality of our polynomials, we discovered a new recursive construction of Sobol' matrices, derived by iteratively squaring polynomials. 
This, in turn, led to the identification of interesting patterns that characterize these sequences.
This paper explores the depths of this new construction and demonstrates how our sequences can be applied to quasi-Monte Carlo rendering. Code and data are available at \href{https://github.com/liris-origami/OneTwoSobolSequences}{https://github.com/liris-origami/OneTwoSobolSequences}.
\begin{figure}[!h]
    \begin{overpic}[width=1.5cm]{figures_siggraph25/tms/t0/0_0_4.pdf}

        \put(-25,0){\rotatebox{90}{$(0,4,2)-$net}}
        \put(-25,-130){\rotatebox{90}{$(1,4,2)-$net}}
        \put(-25,-255){\rotatebox{90}{$(2,4,2)-$net}}

    \end{overpic}
    \includegraphics[width=1.5cm]{figures_siggraph25/tms/t0/0_1_3.pdf}
    \includegraphics[width=1.5cm]{figures_siggraph25/tms/t0/0_2_2.pdf}
    \includegraphics[width=1.5cm]{figures_siggraph25/tms/t0/0_3_1.pdf}
    \includegraphics[width=1.5cm]{figures_siggraph25/tms/t0/0_4_0.pdf}\\
    \vspace{.2cm}
    \hrule
    \vspace{.2cm}
    \includegraphics[width=1.5cm]{figures_siggraph25/tms/t1/1_0_3.pdf}
    \includegraphics[width=1.5cm]{figures_siggraph25/tms/t1/1_1_2.pdf}
    \includegraphics[width=1.5cm]{figures_siggraph25/tms/t1/1_2_1.pdf}
    \includegraphics[width=1.5cm]{figures_siggraph25/tms/t1/1_3_0.pdf}
    \vspace{.2cm}
    \hrule
    \vspace{.2cm}
    \includegraphics[width=1.5cm]{figures_siggraph25/tms/t2/2_0_2.pdf}
    \includegraphics[width=1.5cm]{figures_siggraph25/tms/t2/2_1_1.pdf}
    \includegraphics[width=1.5cm]{figures_siggraph25/tms/t2/2_0_2.pdf}
   
    \caption{The definition of a $(t,m,2)$-net of $n=2^m$ samples is that all dyadic intervals of area $2^t/2^m$ contain $2^t$ points. 
    As such, a $(0,4,2)$-net (top) has exactly one (=$2^0$) point in each interval of volume $1/{2^4}$, a $(1,4,2)$-net has exactly two points (=$2^1$) in each interval of volume $1/{2^3}$ (middle), and a $(2,4,2)$-net has four points on intervals of volume $1/{2^2}$ (bottom). 
    Increasing $t$ thus reduces uniformity constraints and produces larger gaps and clusters in the distribution.
    \label{fig:t1}}
\end{figure}

\begin{figure*}[!ht]
       \includegraphics[width=\textwidth]{scripts/t-correlation.pdf}
       \caption{Experimental validation in 2D of the impact of the $t$ value of a Sobol' sequence on various metrics (from left to right, the generalized $L_2$ discrepancy \cite{hickernell1998generalized} and Monte Carlo integration errors on random Gaussians and Heaviside functions). We randomly select a pair of Sobol' polynomials from the first thousand entries of Joe and Kuo  \shortcite{joe2008constructing}, we  evaluate the metrics for each sample count and plot error distribution (box plots) per values $t$ (box plot colors) observed for that sample count (64 realizations for each $t$). We observe a strong correlation between the observed $t$ value of the point set and its quality for Monte Carlo integration, while $t=1$ appears as a good compromise in quality compared to $t=0$. \label{fig:correlation}}
  \end{figure*}
  
\section{Related work}

We summarize fundamentals about Sobol' construction~\cite{Sobol1967distribution} and refer the reader to reference books~\cite{Niederreiter1992random,Lemieux2009MCQMCSampling,Dick_Pillichshammer_2010} for in-depth discussions.

\paragraph{Sobol' construction}
To construct $(t,s)$-sequences, Sobol' proposed a solution based on primitive polynomials. Given $s$ primitive polynomials $p_0, \dots, p_{s-1}$ in the Galois Field of prime base $b$ called $GF(b)$ (think of it as the set of integers modulo $b$), Sobol' constructs $s$ upper triangular matrices $M_{p_0}, \dots, M_{p_{s-1}}$, one per dimension, which are used to obtain each coordinate of the $i^{th}$ sample point $\mathbf{x}_i$ in the sequence. Specifically, the $k^{th}$ coordinate of the $i^{th}$ sample point, $\mathbf{x}_i^k$, is obtained using a matrix-vector product between matrix $M_{p_k}$ and the base $b$ decomposition of the sample point index $i$ interpreted as a column vector, denoted $\bar{i}$ hereafter, $\bar{q} = M_{p_k} \bar{i}$ with $i = \sum_j \bar{i}_{j} b^j$. The sample point coordinate is now 
$$\mathbf{x}_i^k = \sum_j \bar{q}_j b^{-j}\,.$$
The construction of the upper triangular invertible matrix $M_{p_k}$ using the primitive polynomial $p_k$ uses a recursive formula to obtain a column given its $e$ previous columns, where $e$ is the degree of the polynomial $p_k$. 
The initialization of this recursion, an $e \times e$ upper triangular matrix at the top-left corner of matrix $M_{p_k}$, provides additional degrees of freedom in addition to the chosen polynomial. 
We base our construction on matrix blocks instead of the more common column-wise recursion, as proposed by Faure and Lemieux~\shortcite{faure2016irreducible}, that we briefly describe in Sec.~\ref{sec:blocklemieux}.

Joe and Kuo numerically optimized top-left $e\times e$ blocks, resulting in improved Sobol' sequences on consecutive projections~\cite{joe2008constructing}. 
Matrices can be directly obtained without Sobol' recurrence using an integer linear program solver, but this limits their use to only moderately large problem~\cite{matbuilder22}.

Faure and Lemieux showed that the larger set of irreducible polynomials can be used instead of primitive polynomials~\cite{Sloane2001,faure2016irreducible}. 
Irreducible polynomials are similar to prime numbers, meaning they cannot be factored into products of other non-constant polynomials. 
Faure and Lemieux showed that the parameter $t$ of the resulting $(t,s)$-sequence is bounded by the sum of the polynomial degrees minus one.
A simple way to obtain $(0,b)$-sequences in base $b$ consists of using the first $b$ irreducible polynomials $p_0(x)=x$, $p_1(x)=x+1$, $\dots$, $p_{b-1}(x)=x+b-1$, each of degree $1$. 
The theorem by Faure and Lemieux then shows that $0 \leq t \leq \sum_i \left(deg(p_i)-1\right) = 0$. 
However, this produces a unique sequence (up to sample permutations), related to those produced by Faure~\shortcite{Faure1982}, which limits its use in more general settings that require sample diversity. 
Ostromoukhov et al.~\shortcite{Quad2024} used a construction with quadruplets of irreducible polynomials in base $b=3$ to achieve progressive point sets of excellent consecutive projections.

\paragraph{LDS and projective LDS in Computer Graphics.} In rendering applications, low-discrepancy sequences can have a significant impact on path-tracing performance~\cite{keller2004myths,keller2013quasi,christensen2018renderman,jarosz2019orthogonal}. 
When the sampling pattern defined on the canonical domain $[0,1)^s$ is mapped to a pixel (or a group of pixels), decorrelating the pattern across different pixels typically requires a scrambling procedure. 
Owen's scrambling is usually considered, as it preserves the $t$ value of the point set \cite{owen1995randomly}. 
Due to the nature of the rendering equation, several authors have explored projective strategies aimed at achieving highly uniform consecutive 2D projections. 
Achieving high-quality in 2D projections often comes at the cost of degrading uniformity in higher dimensions~\cite{kollig2002efficient,perrier2018sequences,paulin2021cascaded,ahmed2020screen}. Notably, the ZeroTwo sequence uses the first two Sobol' dimensions repeatedly with random permutations~\cite{pharr2023physically}, while padded 4D  Sobol' repeats and shuffles samples of the first four dimensions~\cite{burley2020practical}. These provide ideal behavior in consecutive 2 or 4D projections, but behave similarly to white noise in higher dimensions. 
Some methods are dedicated to generating point sets rather than sequences~\cite{matbuilder22,Quad2024}, or are not low discrepancy \cite{reinert2016projective,paulin2020sliced}.  
Our new construction enables the definition of complete $(t,s)-$sequences with guaranteed high-quality (\emph{i.e.} $(1,2)$-sequences) 2D projections.


\section{Overview of our construction}



We first focus on 2-dimensional Sobol' sequences. Our goal is thus to obtain Sobol' matrices $M_{p_0}$ and $M_{p_1}$ for two polynomials $p_0$ and $p_1$. In our framework, we may use standard Sobol' construction to generate  these matrices using respectively polynomials $p_0$ and $p_1$ and initialization matrices $D_{p_0}$ and $D_{p_1}$: our contribution is to provide simple conditions for these initialization matrices to yield high-quality parameter $t=1$ Sobol' sequences (see Fig.~\ref{fig:mainIteration}).
This section describes an overview of this process, while Sec.~\ref{sec:mathconstruction} details proofs. In the following, we assume modulo 2 arithmetic, and all values involved are binary. % (we then have $x+x=0$ and $-x=x$).

We consider a matrix $K = M_{p_1} M_{p_0}^{-1}$ that uniquely represents a 2D Sobol' sequence up to a permutation of points, called characteristic matrix (denoted $C$ by Ahmed et al.~\shortcite{ahmed2023analysis}). 
We show that when $p_1 = p_0^2+p_0+1$, where $p_0$ is a degree $e$ polynomial and $p_1$ thus has degree $2e$, $K$ has a peculiar form, and can be obtained recursively from a decomposition into square blocks $A$, $B$ and $C$:  % i starts at 1 here


$$
K^{(i)} = 
\left[\begin{array}{cc}
A & B\\
0 & C
\end{array}\right]
\rightarrow K^{(i+1)} =
\left[\begin{array}{cc|cc}
A & B & A+B & A\\
0& C & C & 0 \\ \hline
0&0& A & A+B\\
0& 0&0  & C
\end{array}\right]\,.
$$ 
where each block $A$, $B$ or $C$ has size $2^{i-1}e$. This produces a matrix $K$ of arbitrary size, doubling its size at each iteration. Our paper introduces conditions on $K^{(1)}$ and $K^{(2)}$, that lead to conditions on the Sobol' initialization matrices $D_{p_0}$ and $D_{p_1}$, for our Sobol' sequences to be $(1,2)$-sequences.
We note that the initialization $K^{(1)}$, of size $2e \times 2e$, is $ K^{(1)}= D_{p_1} \widetilde{D}_{p_0}^{-1}$, where we denote $\widetilde{D}_{p_0}$ the $2e \times 2e$ matrix obtained by extending the $e \times e$ matrix $D_{p_0}$ to $2e \times 2e$ using standard Sobol' iterations.

To obtain Sobol' initialization matrices, we thus first generate a random invertible triangular $e \times e$ matrix $D_{p_0}$ which we extend to $2e \times 2e$ using Sobol' iterations. We then use a matrix $K^{(1)}$ satisfying our conditions (see next), and easily obtain $2e \times 2e$ initialization matrix $D_{p_1} = K^{(1)} \widetilde{D}_{p_0}$. Matrices $M_{p_0}$ and $M_{p_1}$, and 2D sample coordinates are then obtained using standard Sobol' procedures, from $p_0$, $p_1=p_0^2+p_0+1$, $D_{p_0}$ and $D_{p_1}$.


\begin{figure}

\definecolor{dx}{RGB}{255, 0, 255}
\definecolor{ddx}{RGB}{255, 128, 128}
\definecolor{dy}{RGB}{128, 0, 128}

    \begin{tikzpicture}
\node (itera) at (0,0)  {\includegraphics[width=\linewidth]{figures_siggraph25/mainIteration.pdf}};
\node[text=dx] (dp) at (-4, -1) {$D_{p}$};
\node[text=ddx] (dpp) at (-2.7, -1) {$D_{p^2}$};
\node (dppp) at (-1.7, -2) {${\color{dy}D_{p^2+p+1}} = K^{(1)}D_{p^2}$};
\node (k1) at (-3.5,2.2) {$K^{(1)}$};
\node (k2) at (-1.6,2.2) {$K^{(2)}$};
\node (k3) at (1.8,2.2) {$K^{(3)}$};

\draw [draw=gray,decorate,decoration={brace,amplitude=5pt,raise=4ex}]
  (-0.2,2) -- (4.5,2) node[color = gray,midway,yshift=3em]{Only used for proofs (Theorem \ref{theorem:constructionRanks})};


  
\draw [draw=gray,very thick] (-4.4,2.5) rectangle (-.4,-.3);
\node[text=gray] at (-.85,2.3) {Alg.~\ref{algo:CharacteristicMatrices}};
  
\draw [draw=gray,very thick] (-4.4,-.6) rectangle (-.4,-2.5);
\node[text=gray] at (-.85,-.8) {Alg.~\ref{algo:ConstructingDxDy}};

  \draw [-{Stealth[length=2mm]}](-3.6,.8) to [bend right=15] (dp);
%\draw [-{Stealth[length=2mm]}] (-3.4,.8) to [bend left=15] (dppp);
\draw [-{Stealth[length=2mm]}] (dp) -- (dpp);
\draw [-{Stealth[length=2mm]}] (dpp) -- (dppp);

\draw [draw=gray,very thick] (-4.4,-2.7) rectangle (4,-7);
\node[text=gray] at (2,-2.9) {Classical Sobol' construction};
\node (mp) at (-3.8, -3.5) {$M_{p}$};
\node (m1) at (-2,-5.2) {\includegraphics[width=3cm]{figures_siggraph25/matrixM1.pdf}};
\node (mppp) at (0.2, -3.5) {$M_{p^2+p+1}$};
\node (m2) at (2,-5.2) {\includegraphics[width=3cm]{figures_siggraph25/matrixM2.pdf}};
\draw [-{Stealth[length=2mm]}] (dp) to [bend right=15] (mp);
\draw [-{Stealth[length=2mm]}] (dppp) -- (mppp);

\node (pts) at (0,-9) {\includegraphics[width=3cm]{figures_siggraph25/pts56.pdf}};

\draw [very thick] (1.55,-10.45) -- (1.55,-7.55);
\draw [-{Stealth[length=2mm]}] (m2) to [bend left=25] (1.65,-9);

\draw [very thick] (-1.5,-7.45) -- (1.5,-7.45);
\draw [-{Stealth[length=2mm]}] (m1) to [bend left=15] (0,-7.4);

\end{tikzpicture}
\caption{Overview. We introduce a recursive construction of the characteristic matrix associated with a pair of polynomials ($p,p^2+p+1$). We use it in proofs to obtain conditions for generating $(1,2)$-sequences based on the first iteration alone of this recursion.
From characteristic matrices meeting these conditions, we 
derive Sobol' initialization matrices $D_{p}$ and $D_{p^2+p+1}$,  which in turn allows to construct the corresponding Sobol' matrices $M_{p}$ and $M_{p^2+p+1}$ generating $(1,2)-$sequences in base 2.}\label{fig:mainIteration} 
\end{figure}


To generate higher-dimensional Sobol' sequences, we combine pairs of dimensions but further require that $p_0$ and $p_1 = p_0^2+p_0+1$ be irreducible polynomials so as to guarantee that the resulting $s$-dimensional sequence is a $(t,s)$-sequence~\cite{faure2016irreducible}.


We claim the following contributions.


\begin{theorem}\label{theorem:construction}
The sequence of iterations 
\begin{equation}
K^{(i)} = 
\left[\begin{array}{cc}
A & B\\
0 & C
\end{array}\right]
\rightarrow K^{(i+1)} =
\left[\begin{array}{cc|cc}
A & B & A+B & A\\
0& C & C & 0 \\ \hline
0&0& A & A+B\\
0& 0&0  & C
\end{array}\right], 
\label{eq:KiKiplusone}
\end{equation}
where each matrix block $A$, $B$, $C$, is of size $2^{i-1}e \times 2^{i-1}e$, produces the characteristic matrix of a 2D Sobol' sequence given by a degree $e$ polynomial $p$, and degree $2e$ polynomial $p^2+p+1$.
\end{theorem}

\begin{corollary}\label{corollary:coeffs}

A 2D Sobol' sequence given by polynomials $p$ and $p^2+p+1$ only depends on the degree of the polynomial and initialization matrices, and does not depend on the coefficients of $p$ themselves, up to a permutation of samples. 
\end{corollary}

This corollary is readily justified since the recursive construction of theorem~\ref{theorem:construction} does not involve polynomial coefficients, but merely polynomial degrees.
As such, a pair $p_0$ of degree $e$ and $p_1 = p_0^2+p_0+1$ with given initialization matrices would result in the same sequence as $p_0'$ of degree $e$ and $p_1' = p_0'^2+p_0'+1$ for \textit{another pair} of initialization matrices.

%$p_0(x)=x^e$ and $p_1(x)=x^{2e}+x^e+1$ produces the same Sobol' sequence of samples as any other pair of polynomials of the same degree satisfying $p_1=p_0^2+p_0+1$ if initialization matrices are the same. \nb{PAS BON! A REVOIR PROPREMENT}


\begin{theorem}\label{theorem:constructionRanks}
Iterations $K^{(i)} \rightarrow K^{(i+1)}$ characterize Sobol' $(1,2)$-sequences if and only if both conditions are met:
\begin{itemize}
\item $(\mathcal{P})$: $\emph{\text{corank}}(T) \leq 1$ for any rectangular $(w-1) \times w$ submatrix $T$ of $K^{(2)}$ anchored at its first row
\item $(\mathcal{Q})$: $\emph{\text{corank}}(C') \leq 1$ for any square submatrix $C'$ of block $C$ in $K^{(1)}$ obtained by removing any consecutive set of $k$ columns and the last $k$ rows of $C$.
\end{itemize}
\end{theorem}

These results allow us to pre-compute many matrices $K^{(1)}$ satisfying the conditions of theorem~\ref{theorem:construction} for a given polynomial degree $e$ and infer pairs of irreducible polynomials and initialization matrices that produce Sobol' $(1, 2)$-sequences.

In the process of proving these theorems, we discovered new insights on Sobol' constructions for more general polynomials:

\begin{itemize}
\item Polynomials $p$ and $p^2$ will produce the same Sobol' matrices, given proper initialization matrices (see lemma~\ref{Lemma1:doubling}).
\item This property allows building general Sobol' matrices in a recursive way by doubling their size at each iteration via polynomial squaring.
\item The characteristic matrix can also be constructed recursively.
\end{itemize}

The goal of Sec.~\ref{sec:mathconstruction} is to provide mathematical proofs for the claims we summarized in this overview section. The practitioner may thus skip Sec.~\ref{sec:mathconstruction} at first read and jump to the description of our algorithms in Sec.~\ref{sec:experimental}.
Specifically, sections ~\ref{sec:blocklemieux} and \ref{sec:squaring} prove lemma~\ref{Lemma1:doubling} related to polynomials $p$ and $p^2$ yielding the same matrices. This in turns helps proving in Sec.~\ref{sec:blockchar} that characteristic matrices can be constructed recursively. When applied to polynomials $p$ and $p^2+p+1$, Sec.~\ref{sec:pp2p1} proves that the characteristic matrix has the recursive construction of Eq.~\ref{eq:KiKiplusone} and thus proves Theorem~\ref{theorem:construction}. Finally, Sec.~\ref{sec:ranks} and our supplementary materials prove Theorem~\ref{theorem:constructionRanks} that explicits conditions under which the characteristic matrices generate $(1, 2)$-sequences.



\section{Construction of $(1,2)$-sequences in base $2$}
\label{sec:mathconstruction}

We first recall a construction of Sobol' matrices based on matrix blocks by Faure and Lemieux~\shortcite{faure2016irreducible} in Sec.~\ref{sec:blocklemieux}, which will serve as a basis for the next sections describing our proof. Our proof first consists of introducing a new recursive construction of Sobol' matrices by squaring polynomials in Sec.~\ref{sec:squaring}. We then show that a similar squaring procedure can be obtained for characteristic matrices in Sec.~\ref{sec:blockchar}. We then show, using this construction, that the characteristic matrix for polynomials $p$ and $p^2+p+1$ has a specific form exhibiting a self-similar pattern in Sec.~\ref{sec:pp2p1}. We finally show that ranks of characteristic matrices with this self-similar pattern are necessarily such that the produced 2D sequences are $(1,2)$-sequences in Sec.~\ref{sec:ranks}.


\subsection{Block-based recursive construction}
\label{sec:blocklemieux}
We first briefly describe a block-based 1-D Sobol' construction as described by Faure and Lemieux~\shortcite{faure2016irreducible}. 

For a given irreducible polynomial $p(x) = x^e + \sum_{i=0}^{e-1} a_i x^i$ and upper $e\times e$ invertible triangular initialization matrix $D_p$, Faure and Lemieux~\shortcite{faure2016irreducible} rewrite Sobol' iterations in terms of block matrices:

$$
M_p = 
\begin{pmatrix}
    B_{1,1} & B_{1,2} & B_{1,3} & \dots   \\
    0 & B_{2,2} & B_{2,3} & \dots   \\
    0 & 0 & B_{3,3} & \dots   \\
    \vdots & \vdots & \vdots & \ddots 
\end{pmatrix},
$$
where the blocks $B_{i,j}$, of size $e\times e$, are defined according to the following recursive procedure:
$$
B_{1,1} = D_p ;  \hspace{10pt} B_{i,i}  = D_p F_p^{i-1}  \\
$$
\begin{equation}
B_{i,j} = \left\{ 
		\begin{array}{lll} 
			B_{i,j-1} Q_p F_p && \text{when} \hspace{5pt} i=1 \\
			B_{i,j-1} Q_p F_p +  B_{i-1,j-1} F_p && \text{elsewhere}.
		\end{array}
	\right.
\label{eq:Blr}
\end{equation}

Here, $Q_p$ is an $e\times e$ lower triangular Toeplitz matrix involving the coefficients $(a_i)_{i}$ of polynomial $p$:
\begin{equation}
Q_p = 
\begin{pmatrix}
    a_{0} & 0 & 0 & \dots & 0 \\
    a_{1} & a_{0} & 0 & \dots  & 0 \\
    a_{2} & a_{1} &a_{0} & \dots  & 0 \\
    \vdots & \vdots & \vdots & \ddots & \vdots \\
    a_{e-1} & a_{e-2} & a_{e-3} & \dots  & a_{0} 
\end{pmatrix}.
\label{eq:Qp}
\end{equation}
Matrix $F_p$ of size $e\times e$ is defined as 
\begin{equation}
F_p = ({Id}_{\small e} + R_{p,2} ) ({Id}_{\small e} + R_{p,3} ) \dots ({Id}_{\small e} + R_{p,e} ),
\label{eq:F}
\end{equation}
where $Id_{\small e}$ is an identity matrix of size $e\times e$, and $R_{p, k}$ are matrices of size $e\times e$  
with zeros everywhere except in the first $k-1$ entries of the $k$-th column, given by the coefficients $(a_{e-(k-1)},\dots ,a_{e-1})$ of polynomial $p$.



We introduce a new recursion to build Sobol' matrices inspired by the construction of Faure and Lemieux~\shortcite{faure2016irreducible}.




\subsection{Sobol' construction by squaring polynomials}
\label{sec:squaring}

In the following, we introduce a squared superscript notation to clarify matrix sizes when appropriate, e.g., $M_{p}^\btwoe$ denoting the Sobol' matrix of polynomial $p$ restricted  to the first $2e$ rows and $2e$ columns.


We observe that the Sobol' matrix $M_{p^2}$ of a squared polynomial (although not irreducible) is identical to the Sobol' matrix $M_p$ of the original polynomial, provided that the initialization matrix $D_{p^2}$ coincides with the top-left corner of $M_p$. We formalize this:

\begin{lemma}\label{Lemma1:doubling}
The Sobol' matrix $M_p$ generated by a polynomial $p$ of degree $e$ and initialization matrix $D_p$ is identical to the Sobol' matrix $M_{p^2}$  of polynomial $p^2$ and initialization matrix 
\begin{equation}
    D_{p^2} = M^\btwoe_{p} = \begin{pmatrix}
    D_p & 0 \\
    0 & D_p 
\end{pmatrix} \begin{pmatrix}
    Id_{\small e} & Q_p F_p \\
    0 & F_p 
\end{pmatrix},\label{eq:lemma41}
\end{equation}
corresponding to the top-left $2e \times 2e$ submatrix of $M_p$.

Matrices $Q_{p^2}$ and $F_{p^2}$ of Faure and Lemieux~\shortcite{faure2016irreducible} can be obtained by applying a Kronecker product (or tensor product) with the matrix ${Id}_2 = {\begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}}$ to matrices $Q_p$ and $F_p$:
\begin{equation}
Q_{p^2} = Q_{p}  \otimes {Id}_2, \hspace{20pt} F_{p^2} = F_{p}  \otimes {Id}_2.
\label{eq:FQ-prime}
\end{equation}
where $Q_{p^2}$ and $F_{p^2}$ are of size $2e \times 2e$.
\end{lemma}

\begin{figure*}
    \begin{tikzpicture}
        \draw[very thick] (2,-1.5) rectangle (3.5,0) node[pos=.5] {0};
        \draw[very thick] (3.5,-1.5) rectangle (5,0) node[pos=.5] {${\color{teal}M_p^\be} Y^{(1)}$};
        \draw [very thick](3.5,0) rectangle (5,1.5) node[pos=.5] {${\color{teal}M_p^\be} X^{(1)}$};

        \draw[fill=teal,color=teal,very thick] (2,0) rectangle (3.5,1.5) node[pos=.5] {\color{white}\small$M_p^\be=D_p$};

        \draw[very thick] (5,-1.5) rectangle (8,1.5) node[pos=.5] {${\color{magenta}M_p^\btwoe} X^{(2)}$};
        \draw[very thick] (5,-4.5) rectangle (8,-1.5) node[pos=.5] {${\color{magenta}M_p^\btwoe} Y^{(2)}$};
        \draw[very thick] (2,-4.5) rectangle (5,-1.5) node[pos=.5] {0};
        \node at (3.5,2) {${\color{magenta}M_p^\btwoe}$};
        \draw[color=magenta,very thick] (2,-1.5) rectangle (5,1.5);
        \node at  (8.5,1.5) {\ldots\ldots};
        \node at  (2.0,-4.9) {\rotatebox{90}{\ldots\ldots}};
        \node at (3.5,-5) {$X^{(1)} = Q_pF_p$};
        \node at (6.8,-5) {$X^{(i)} = X^{(i-1)}\otimes Id_2$};
        \node at (3.3,-5.5) {$Y^{(1)} = F_p$};
        \node at (6.8,-5.5) {$Y^{(i)} = Y^{(i-1)}\otimes Id_2$};

        \draw [color=gray,decorate,decoration={brace,amplitude=5pt,raise=.5ex}]
        (2,0) -- (2,1.5) node[midway,xshift=-1em]{$e$};
        \draw [color=gray,decorate,decoration={brace,amplitude=5pt,raise=1ex}]
        (1.8,-1.5) -- (1.8,1.5) node[midway,xshift=-1.5em]{$2e$};
        \draw [color=gray,decorate,decoration={brace,amplitude=5pt,raise=1ex}]
        (1.4,-4.5) -- (1.4,1.5) node[midway,xshift=-1.5em]{$4e$};

    \end{tikzpicture}\hspace{.5cm}
\begin{tikzpicture}
    \draw[very thick] (2,-1.5) rectangle (3.5,0) node[pos=.5] {0};
    \draw[very thick] (3.5,-1.5) rectangle (5,0) node[pos=.5] {${\color{teal}K_{q,r}^\bet} Y'^{(1)}$};
    \draw[very thick] (3.5,0) rectangle (5,1.5) node[pos=.5] {${\color{teal}K_{q,r}^\bet} X'^{(1)}$};

    \draw[fill=teal,color=teal,very thick] (2,0) rectangle (3.5,1.5) node[pos=.5] {\tiny\color{white}$K_{q,r}^\bet=D_rD_q^{-1}$};

    \draw[very thick] (5,-1.5) rectangle (8,1.5) node[pos=.5] {${\color{magenta}K_{q,r}^\btwoet} X'^{(2)}$};
    \draw[very thick] (5,-4.5) rectangle (8,-1.5) node[pos=.5] {${\color{magenta}K_{q,r}^\btwoet} Y'^{(2)}$};
    \draw[very thick] (2,-4.5) rectangle (5,-1.5) node[pos=.5] {0};
    \node at (3.5,2) {${\color{magenta}K_{q,r}^\btwoet}$};
    \draw[very thick,color=magenta,very thick] (2,-1.5) rectangle (5,1.5);
    \node at  (8.5,1.5) {\ldots\ldots};
        \node at  (2.0,-4.9) {\rotatebox{90}{\ldots\ldots}};
    \node at (4.2,-5) {$X'^{(1)} = D_{q} (Q_{q} + Q_{r}F_{r}F_{q}^{-1} ) D_{q}^{-1}$};
    \node at (7.8,-5) {$X'^{(i)} = X'^{(i-1)}\otimes Id_2$};
    \node at (3.7,-5.5) {$Y'^{(1)} = D_{q} ( F_{r}F_{q}^{-1} ) D_{q}^{-1}$};
    \node at (7.8,-5.5)  {$Y'^{(i)} = Y'^{(i-1)}\otimes Id_2$};
    \draw [color=gray,decorate,decoration={brace,amplitude=5pt,raise=.5ex}]
    (2,0) -- (2,1.5) node[midway,xshift=-1em]{$e$};
    \draw [color=gray,decorate,decoration={brace,amplitude=5pt,raise=1ex}]
    (1.8,-1.5) -- (1.8,1.5) node[midway,xshift=-1.5em]{$2e$};
    \draw [color=gray,decorate,decoration={brace,amplitude=5pt,raise=1ex}]
    (1.4,-4.5) -- (1.4,1.5) node[midway,xshift=-1.5em]{$4e$};
\end{tikzpicture}
        \caption{In contrast to the column-by-column Sobol’ approach~\cite{Sobol1967distribution}, and the block formulation of Faure and Lemieux~\shortcite{faure2016irreducible} (Sec.~\ref{sec:blocklemieux}), we present a novel recursive construction for Sobol’ matrices (left) and characteristic matrices (right) by squaring polynomials.
         The construction of matrix $M_p$ follows from lemma~\ref{Lemma1:doubling}, Eq.~(\ref{eq:lemma41}), while the construction of $K^\btwoet_{q,r}$ is obtained from Eq.~(\ref{eq:C}). Matrices $X^{(i)}$, $Y^{(i)}$, $X'^{(i)}$ and $Y'^{(i)}$  provide compact representations of expressions involving $D_p$, $F_p$, and $Q_p$, using the distributivity of the Kronecker product over matrix multiplication.} 
        \label{fig:doublingMK}
\end{figure*}


This lemma brings a new recursive construction of Sobol' matrices, doubling their sizes at each iteration by squaring polynomials, illustrated in Fig.~\ref{fig:doublingMK}. 
%To double the size of an $m \times m$ Sobol' matrix, where $m=2^k e$, we consider this matrix as an initialization matrix $D^{\small [m]}_{p^{m/e}}$ for a polynomial of degree $m$, and use $Q^{\small [m]}_{p^{m/e}} = Q^\be_{p} \otimes {Id}_{m/e}$ and  $F^{\small [m]}_{p^{m/e}} = F^\be_{p} \otimes {Id}_{m/e}$ in a single iteration of the block-based formulation of Faure and Lemieux~\shortcite{faure2016irreducible} to obtain the corresponding three $m \times m$ blocks required to double its size.


In our development, we first note that, while $F_p$ can have a complicated form, its inverse can be expressed much more easily, as an upper triangular Toeplitz matrix:
\begin{equation}
F_p^{-1} = 
\begin{pmatrix}
    1 & a_{e-1} & a_{e-2} & \dots & a_{1} \\
    0 & 1 & a_{e-1} & \dots   & a_{2} \\
    0 & 0 &1  & \dots  & a_{3}  \\
    \vdots & \vdots & \vdots  & \ddots & \vdots \\
    0 & 0 & 0 & \dots  & 1 
\end{pmatrix},
\label{eq:Q}
\end{equation}
where we ignore signs in modulo 2 arithmetic. This is obtained by observing that matrices ${Id}_e + R_{p,k}$ are their own inverses, and that inverting $F_p$ then merely reverses the order of the multiplication: $F_p^{-1} = ({Id}_{\small e} + R_{p,e} ) \dots ({Id}_{\small e} + R_{p,2} )$.


\begin{proof}[Proof of Lemma \ref{Lemma1:doubling}] 
The relation between $D_{p^2}$ and $D_{p}$ is simply obtained by running Faure and Lemieux' iterations for one block of columns.
We may now assume that the top-left $2e \times 2e$ submatrices of $M_p$ and $M_{p^2}$ coincide.
We also note that for a polynomial $p(x) = x^e + \sum_{i=0}^{e-1} a_i x^i$, given modulo-2 arithmetic cancelling odd degrees, $p^2(x) = x^{2e} + \sum_{i=0}^{e-1} a_i x^{2i}$.
This, in turn, leads to $Q_{p^2} = Q_{p} \otimes Id_2$, and similarly, to $F^{-1}_{p^2} = F^{-1}_{p} \otimes {Id}_2$ and thus $F_{p^2} = F_{p} \otimes {Id}_2$ (the inverse of the Kronecker product being the Kronecker product of the inverse matrices), and finally, $Q_{p^2} F_{p^2} = (Q_{p} F_{p})\otimes {Id}_2$.
The recursive construction of Eq.~(\ref{eq:Blr}) thus produces the same matrices whether using $Q_{p^2}$ and $F_{p^2}$ or $Q_{p}$ and $F_{p}$.

\end{proof}



\subsection{Block-based characteristic matrices}
\label{sec:blockchar}


For a pair of dimensions  with Sobol' matrices $M_{q}$ and $M_{r}$, we base our analysis on the characteristic matrix $K = M_{r} M_{q}^{-1}$ as defined by Ahmed et al.~\shortcite{ahmed2023analysis}
, which uniquely characterizes Sobol' 2D sequences up to a permutation of samples.
When the polynomials $q$ and $r$ are of the same degree $\tilde{e}$, we may build a characteristic matrix of size $2\tilde{e}$ by applying one iteration of Faure and Lemieux' block construction:

\begin{align}
        K^{\btwoet}_{q, r} &=  M^\btwoe_{r}  \left(M^{\btwoe}_{q}\right)^{-1}\nonumber \\
         &= {\begin{pmatrix} D_{r} & 0 \\ 0 & D_{r}\end{pmatrix}} {\begin{pmatrix} Id_{\small \tilde{e}} & Q_{r} F_{r} \\ 0 & F_{r}\end{pmatrix}} \left(  {\begin{pmatrix} D_{q} & 0 \\ 0 & D_{q}\end{pmatrix}} {\begin{pmatrix} Id_{\small \tilde{e}} & Q_{q} F_{q} \\ 0 & F_{q}\end{pmatrix}} \right)^{-1} \nonumber \\	
        &=  {\begin{pmatrix} D_{r} & 0 \\ 0 & D_{r}\end{pmatrix}} {\begin{pmatrix} Id_{\small \tilde{e}} & Q_{q} + Q_{r} F_{r} F_{q}^{-1} \\ 0 & F_{r} F_{q}^{-1}\end{pmatrix}}    {\begin{pmatrix} D_{q} & 0 \\ 0 & D_{q}\end{pmatrix}}^{-1} \nonumber \\
         %= {\begin{pmatrix} D_{r} & 0 \\ 0 & D_{r}\end{pmatrix}}  {\begin{pmatrix} R_{11} & R_{12} \\ 0 & R_{22} \end{pmatrix}}    {\begin{pmatrix} D_{q}^{-1} & 0 \\ 0 & D_{q}^{-1}\end{pmatrix}}  \\
         &= {\begin{pmatrix} K^{\bet}_{q, r} & 0 \\ 0 & K^{\bet}_{q, r}\end{pmatrix}}   {\begin{pmatrix} Id_{\small \tilde{e}} & D_{q} (Q_{q} + Q_{r}F_{r}F_{q}^{-1} ) D_{q}^{-1} \\ 0 & D_{q} ( F_{r}F_{q}^{-1} ) D_{q}^{-1}\end{pmatrix}.} 	 
     \label{eq:C}
    \end{align}
where we used the identity $
\begin{array}{lll} 
	\begin{pmatrix} A & B \\ 0 & C\end{pmatrix}^{-1} = 
	\begin{pmatrix} A^{-1} & -A^{-1}  B C^{-1} \\ 
	0 & C^{-1} \end{pmatrix}
\end{array}
$ for invertible $A$ and $C$, and where $K^{\bet}_{q, r} = D_{r}  D_{q}^{-1}$ is the $\tilde{e} \times \tilde{e}$ initialization block for this characteristic matrix, obtained from the initializations of  $M_{q}$ and $M_{r}$.

This construction considers polynomials $q$ and $r$ of the same degree, but we can use lemma~\ref{Lemma1:doubling} to square our first polynomial and then consider $q = p^2$ and $r = p^2+p+1$ of the same degree. We thus consider this construction using $\tilde{e}=2e$.

Also, this construction merely produces a matrix of size $2\tilde{e}$ given polynomials of degree $\tilde{e}$. However, it can be used recursively by considering lemma~\ref{Lemma1:doubling}, doubling the size of the matrix at each iteration by squaring polynomials. By lemma~\ref{Lemma1:doubling}, the effect of squaring polynomials on all matrices involved is merely a tensor product with $Id_2$. We also illustrate this process in Fig.~\ref{fig:doublingMK}.



\subsection{The special case $(p, p^2+p+1)$}
\label{sec:pp2p1}

The construction for characteristic matrices in Sec.~\ref{sec:blockchar} is general, and applies to any pair of polynomials. In this section, we show the special case when $q = p^2$, and $r = p^2+p+1$.

First, we note that applying lemma~\ref{Lemma1:doubling} for using ${q} = p^2$ in place of ${q} = p$ leads to an extended initialization matrix $D_{p^2}$ which inverse can be expressed: 
\begin{equation}
D_{p^2}^{-1} = \left(M^{\btwoe}_p\right)^{-1} =  \begin{pmatrix}
    Id_{\small e} & Q_p  \\
    0 & F_p^{-1} 
\end{pmatrix}\begin{pmatrix}
    D_p^{-1} & 0 \\
    0 & D_p^{-1} 
\end{pmatrix}. \label{eq:invDp2}
\end{equation}
We thus now consider that our polynomials have degree $2e$ and we seek to apply our characteristic matrix construction to obtain a matrix 	$K^{\bfoure}_{q, r}$ of size $4e \times 4e$.

It can then be shown (see Appendix \ref{app:proofsR}) that $F_{p^2+p+1} F_{p^2}^{-1}$ has a particular form:
\begin{equation}
F_{p^2+p+1} F_{p^2}^{-1} = \begin{pmatrix}
    Id_{\small e} & F_p  \\
    0 & Id_{\small e} 
\end{pmatrix}.\label{eq:FrinvFp}
\end{equation}
The resulting matrix is also its own inverse: $F_{p^2+p+1} F_{p^2}^{-1} = F_{p^2} F_{p^2+p+1}^{-1}$.
Combining Eq.~(\ref{eq:invDp2}) and (\ref{eq:FrinvFp}) (see Appendix \ref{app:proofs1}), we have: 
\begin{align} 
F_{p^2+p+1} F_{p^2}^{-1} D_{p^2}^{-1} &= D_{p^2}^{-1} \begin{pmatrix}
    Id_{\small e} & Id_{\small e}  \\
    0 & Id_{\small e} 
\end{pmatrix}.\label{eq:FrinvFpDp2}
\end{align}
Similarly, we have 
\begin{equation}
(Q_{p^2} + Q_{p^2+p+1} F_{p^2+p+1} F_{p^2}^{-1}) D_{p^2}^{-1} = D_{p^2}^{-1}  \begin{pmatrix}
    Id_{\small e} & Id_{\small e}  \\
    Id_{\small e} & 0 
\end{pmatrix},\label{eq:Qp2plus}
\end{equation}
where we further use the properties of the product of our Toeplitz matrices, see Appendix \ref{app:proofs2} for details.

Putting all together, and considering the initially given matrix $K^{\btwoe}_{p^2, p^2+p+1}$ has size $2e$ since we considered $p^2$ and $p^2+p+1$ of degree $2e$, we see that the recursive construction of $K$ becomes 
\begin{equation}
\begin{array}{lll} 
	K^{\bfoure}_{p^2, p^2+p+1} 
	 = {\begin{pmatrix} K^{\btwoe}_{p^2, p^2+p+1} & 0 \\ 0 & K^{\btwoe}_{p^2, p^2+p+1}\end{pmatrix}}   {\begin{pmatrix} Id_{\small 2e} &\begin{pmatrix}
    Id_{\small e} & Id_{\small e}  \\
    Id_{\small e} & 0 
\end{pmatrix} \\ 0 & \begin{pmatrix}
    Id_{\small e} & Id_{\small e}  \\
    0 & Id_{\small e} 
\end{pmatrix} \end{pmatrix}}\,. 	 
\end{array}
 \label{eq:C2}
\end{equation}

Lemma~\ref{Lemma1:doubling} indicates that this procedure becomes recursive by squaring polynomials, allowing to double the size of $K$ at each iteration. We rewrite these iterations using the following notation involving blocks $A$, $B$, $C$ of size $2^ie \times 2^ie$, doubling their size at each iteration:

$$
K^{(i)} = 
\left[\begin{array}{cc}
A & B\\
0 & C
\end{array}\right]
\rightarrow K^{(i+1)} =
\left[\begin{array}{cc|cc}
A & B & A+B & A\\
 & C & C & 0 \\ \hline
 &  & A & A+B\\
 &  &  & C
\end{array}\right] 
$$ 
with the $2e \times 2e$ initial matrix 
\begin{align*}
K^{(1)} &= K^{\btwoe}_{p^2, p^2+p+1} =  D_{p^2+p+1} D_{p^2}^{-1} \\
&= 
    D_{p^2+p+1} 
\left(\begin{pmatrix}
    D_{p} & 0 \\
    0 & D_{p} 
\end{pmatrix} \begin{pmatrix}
    Id_{\small e} & Q_p F_p \\
    0 & F_p 
\end{pmatrix}\right)^{-1}, \end{align*}
for any given $e \times e$ upper triangular invertible matrix $D_p$ and $2e \times 2e$ upper triangular invertible matrix $D_{p^2+p+1}$ giving the degrees of freedom for the generated sequences. $D_{p^2}$ is here obtained by running standard Sobol' iterations to extend the $e \times e$ matrix $D_{p}$ to obtain the next $e$ rows and columns.

\subsection{Rank of submatrices}
\label{sec:ranks}

Let us denote by $\bar{T}^{j,w}$ a square $w \times w$ submatrix of $K$ 
starting at column $1\leq j<m-w$.
Ahmed et al.~\shortcite{ahmed2023analysis} showed that $(M_p,M_q)$ is a $(0,2)-$sequence if and only if all submatrices $\bar{T}^{j,w}$ have nonzero determinant (modulo 2). In our settings, we define $T^{j,w}$ 
\begin{wrapfigure}{r}{1.9cm} 
\vspace{-10pt}
\begin{tikzpicture}
\draw [draw=blue!50,very thick] (0,0) rectangle (2,2);
\draw [draw=green!50!black,very thick] (0.28,.95) rectangle (1.7,2);
\draw [draw=orange,very thick] (0.3,1) rectangle (1.3,2);
\draw [decorate,decoration={brace,amplitude=2pt,mirror,raise=4ex}] (0.3,1.4) -- (1.3,1.4) node[midway,yshift=-2.5em]{$w$};
\draw [decorate,decoration={brace,amplitude=3pt,mirror,raise=4ex}] (1.3,1.4) -- (1.7,1.4) node[midway,yshift=-2.5em]{$1$};
\node at (0.3,2) [draw=orange, fill=orange,circle,fill,inner sep=1.5pt]{};
\node at (0.25,2.3) [text=orange]{$j$};
\node at (.8,1.5) [text=orange]{$\bar{T}^{j,w}$};
\node at (1.65,2.3) [text=green!50!black]{$T^{j,w}$};
\node at (0.3,0.2) [text=blue!50]{$K$};
\end{tikzpicture}
\end{wrapfigure}
 as the rectangular submatrix of $K$ starting at column $j$ and of size $(w-1) \times w$. A first claim is that 
a $(1,2)$-sequence is characterized by a matrix $K$ having all submatrices $T^{j,w}$ rank-deficient by at most $1$ (see Appendix \ref{app:corank}). 
If we denote $\text{corank}(T^{j,w}) := w-1 - \text{rank}(T^{j,w})$ the rank deficiency of matrix $T^{j,w}$, a $(1,2)$-sequence is thus characterized by the property that all submatrices $T^{j,w}$ of $K$ have a corank of at most $1$.  We call this property $\mathcal{P}$. Ranks need to be computed in $GF(2)$, e.g., the matrix 
${\left[\begin{array}{cccc}
    1 & 1 & 0 & 0\\
    1 & 0 & 1 & 1\\
    0 & 1 & 1 & 1
    \end{array}\right]} $
 has rank 2 in $GF(2)$ (because the third column is the sum of the previous two, modulo 2) although it has rank 3 over the integers. This can be obtained numerically using Gaussian elimination.


Since matrix $K$ is an infinite-sized matrix, systematic numerical evaluation of ranks for all possible submatrices of $K$ quickly becomes intractable.

We instead benefit from our recursive construction of $K$ to propagate properties across iterations. We show by induction that if property $\mathcal{P}$ holds for matrix $K^{(2)}$, and an additional property $\mathcal{Q}$ holds for the block $C$ of $K^{(1)}$, then properties $\mathcal{P}$ and $\mathcal{Q}$ necessarily hold for all $K^{(i)}$, $i\geq 1$.

Property $\mathcal{Q}$ states that all submatrices ${C'}$ of $C$ obtained by removing $1 \leq t < m$ consecutive columns and the last $t$ rows have $\text{corank}(C')\leq 1$. It is easy to verify that property $\mathcal{Q}$ holds for $K^{(i)}$, $i\geq 1$, if it holds for $K^{(1)}$, since the recursive procedure transforms $C$ into a block triangular matrix $\begin{pmatrix}
    A & A+B \\
    0 & C 
\end{pmatrix}$, where $A$ is full rank.

Verifying by induction that property $\mathcal{P}$ holds for $K^{(i)}$, $i\geq 1$, provided that it holds for $K^{(2)}$ (and thus $K^{(1)}$) is more involved. Given a matrix of the form $K^{(i)}$, we iterate our construction to obtain $K^{(i+1)}$ and $K^{(i+2)}$ ; $K^{(i+2)}$ has $8 \times 8$ blocks, each of size $2^{i-1}e \times 2^{i-1}e$. From $K^{(i+2)}$, we extract matrices $T^{j,w}$ that overlap any number $1 \leq b \leq 8$ of consecutive blocks horizontally and either $b$ or $b-1$ blocks vertically. We symbolically perform Gaussian elimination on $T^{j,w}$ to exhibit block triangular structures for which ranks can be easily obtained~\cite{meyer1973generalized}. Specifically, we seek to have any number of blocks on the diagonal with full rank, and, at most, $1$ block fulfilling property $\mathcal{P}$ or $\mathcal{Q}$ by hypothesis to conclude that $\text{corank}(T^{j,w})\leq 1$ .
Given the sheer number of cases, we refer the reader to the supplementary materials for the exhaustive list of cases, and show one typical case in Fig.~\ref{fig:rank} on a submatrix $T = 
T^{j,w}$ overlapping $4$ blocks (out of 8) of $K^{(i+2)}$. The base case of the induction is tested numerically on matrix $K^{(2)}$ (if it holds for $K^{(2)}$ it also holds for $K^{(1)}$).

\begin{figure}
\includegraphics[width=\linewidth]{blocks/bloc44_3h.pdf}
\caption{
Example of an $(s-1) \times s$ submatrix $T$ (in red) of $K^{(i+2)}$ overlapping 4 consecutive blocks horizontally and vertically. Gaussian elimination is performed (here by addition and permutation of blocks of columns) to exhibit a block triangular structure, where one block (in orange) has non-zero determinant, and another block (in pink) whose corank is smaller than 1 since $K^{(i+1)}$ satisfies $\mathcal{P}$ by hypothesis. This proves that this particular submatrix $T$ also has $\text{corank}(T) \leq 1$. In supplementary materials, we exhaustively list all cases for rectangular submatrices $T$ included in $K^{(i+2)}$ to conclude that $K^{(i+2)}$ then satisfies $\mathcal{P}$.
}
\label{fig:rank}
\end{figure}


\section{Experimental results}
\label{sec:experimental}

In this section, we outline practical aspects of our method. First, we provide a brief overview of the process for generating polynomials and initialization matrices for Sobol' sequences, with a focus on ensuring high-quality 2D projections. Next, we present a concrete example demonstrating how the proposed method can be integrated into a typical physically-based rendering (PBR) framework, using PBRT~\cite{pharr2023physically} as a case study.



\subsection{Constructing projective $(1,2)$-sequences }

In our construction, we cannot use all irreducible polynomials as Faure and Lemieux do~\cite{faure2016irreducible}, because we focus on pairs of polynomials $(p_{i}, p_{i+1})$ such that $p_{i+1} = p_{i}^2 + p_{i} + 1$.
We found 346 such pairs of irreducible polynomials of degree up to $e=16$ ($2e = 32$). This allows for 692D sampling with guaranteed-quality 2D projections, and more precisely $(1,2)$-sequences for consecutive pairs of dimensions.
The property of $(1,2)$-sequences for each pair is guaranteed by Theorem~\ref{theorem:construction}, provided that the appropriate initialization matrices are provided.

First, we precompute a set of candidate characteristic matrices $\mathcal{K}_e$ for each degree $e$ of the polynomials we are considering. 
Note that by Corollary~\ref{corollary:coeffs}, at this stage we do not need the specific polynomials involved, but only their degrees. 
The construction of this collection consists in a random search for matrices $A$, $B$, $C$ of size $e \times e$ by verifying that $K^{(2)}$ satisfies property $\mathcal{P}$ and that $C$ satisfies property $\mathcal{Q}$ (see Sec.~\ref{sec:ranks} and  Algorithm~\ref{algo:CharacteristicMatrices}). 
Exploring the space of all matrices $A$, $B$, $C$ that satisfy $\mathcal{P}$ and $\mathcal{Q}$, respectively, becomes infeasible for large $e$, as the search space grows exponentially with $e$. 
For $e\in\{1,2,3,4,5\}$, we found 2, 6, 40, 1688, and 727 matrices, respectively. For higher degrees, we leverage the fact that doubling a matrix in $\mathcal{K}_e$ by squaring the polynomial as described in Eq.~(\ref{eq:KiKiplusone}), provides a candidate matrix for $\mathcal{K}_{2e}$. These matrices are available in the supplementary material.


\begin{figure*}[!h]
    \begin{tabular}{ccc}
    \begin{overpic}[width=4.5cm]{figures_siggraph25/fig-2DProjections_sobol.pdf}
    \put(-15,102) {dims} 
    \put(-7,3) {11}
    \put(-7,13) {10}
    \put(-7,21) {9}
    \put(-7,30) {8}
    \put(-7,39) {7}
    \put(-7,48) {6}
    \put(-7,57) {5}
    \put(-7,66) {4}
    \put(-7,75) {3}
    \put(-7,83) {2}
    \put(-7,93) {1}
    \put(3,102) {0}
    \put(12,102) {1}
    \put(20,102) {2}
    \put(29,102) {3}
    \put(37,102) {4}
    \put(46,102) {5}
    \put(54,102) {6}
    \put(63,102) {7}
    \put(72,102) {8}
    \put(81,102) {9}
    \put(90,102) {10}
\end{overpic}&
    \includegraphics[width=4.5cm]{figures_siggraph25/fig-2DProjections_Lemieux.pdf}&
    \includegraphics[width=4.5cm]{figures_siggraph25/fig-2DProjections_zt.pdf}\\
    (a) Sobol' & (b) Faure-Lemieux & (c) ZeroTwo\\
    \includegraphics[width=4.5cm]{figures_siggraph25/fig-2DProjections_Cascaded.pdf}&
    \includegraphics[width=4.5cm]{figures_siggraph25/fig-2DProjections_IrreducibleGF3.pdf}&
    \includegraphics[width=4.5cm]{figures_siggraph25/fig-2DProjections_Ours.pdf}\\
    (d) Cascaded Sobol'& (e) QuadOptimized GF(3) & (f) Ours
    \end{tabular}
\caption{
We compare consecutive 2D projections for the first 12 dimensions of several constructions:
(a) Sobol' with Joe and Kuo initializations~\cite{joe2008constructing},
(b) Faure and Lemieux~\shortcite{faure2016irreducible,faure2019implementation}, 
(c) the first two Sobol' dimensions, repeated with a random permutation of sample indices~\cite{pharr2023physically}, 
(d) the Cascaded Sobol' approach of Paulin et al.~\shortcite{paulin2021cascaded} (not sequence) 
(e) the Quad-optimized LDS in GF(3) by Ostromoukhov et al.~\shortcite{Quad2024},
and (f) our approach. Here, orange squares designate guaranteed $(0,2)$-progressive or $(0,2)/(1,2)$-sequence properties.
 Blue squares designate optimized 4-tuples of dimensions. 
 Green squares designate additional optimizations, supported by our optimization process (See details in Sec.~\ref{sec:FurtherImprovements}).
For low discrepancy projections, the factor $t$ of each point set is numerically computed and indicated in the upper-right corner of each patch.
}
\label{fig:2dpatches}
\end{figure*}

\begin{figure*}
\centering \setlength{\tabcolsep}{2pt}
\begin{tabular}{ccccc}
   & $m=6$ & $m=8$ & $m=12$ & max $\forall m\in \{2...12\}$\\
\raisebox{.6cm}{\rotatebox{90}{Sobol' / Joe\&Kuo}}&
    \begin{tikzpicture}[spy using outlines={circle, magnification=4, size=2cm, connect spies}]
    \node { \includegraphics[width=4cm]{figures_siggraph25/heatmaps/heatmap_100Dx100D_SobolJK_m6-crop.pdf} };
    \spy [red] on (-.2,-.2)  in node [left] at (1.5,1);
    \begin{overpic}{figures_siggraph25/empty.png}
        %\put(-6000,6000) {dims}
        \put(-5500,5000) {$\longrightarrow$}
        \put(-5700,5300) {$i$}
        \put(-6600,4450) {$j$}
        \put(-6300,4700) {\rotatebox{270}{{$\longrightarrow$}}}
        \end{overpic}
\end{tikzpicture}

&\begin{tikzpicture}[spy using outlines={circle, magnification=4, size=2cm, connect spies}]
    \node { \includegraphics[width=4cm]{figures_siggraph25/heatmaps/heatmap_100Dx100D_SobolJK_m8-crop.pdf} };
    \spy [red] on (-.2,-.2)  in node [left] at (1.5,1);
\end{tikzpicture}
&\begin{tikzpicture}[spy using outlines={circle, magnification=4, size=2cm, connect spies}]
    \node { \includegraphics[width=4cm]{figures_siggraph25/heatmaps/heatmap_100Dx100D_SobolJK_m12-crop.pdf} };
    \spy [red] on (-.2,-.2)  in node [left] at (1.5,1);
\end{tikzpicture}
&\begin{tikzpicture}[spy using outlines={circle, magnification=4, size=2cm, connect spies}]
    \node { \includegraphics[width=4cm]{figures_siggraph25/heatmaps/heatmap_Sobol_max-crop.pdf} };
    \spy [red] on (-.2,-.2)  in node [left] at (1.5,1);
\end{tikzpicture}
\\
\raisebox{1.5cm}{\rotatebox{90}{Ours}}&

\begin{tikzpicture}[spy using outlines={circle, magnification=4, size=2cm, connect spies}]
    \node { \includegraphics[width=4cm]{figures_siggraph25/heatmaps/heatmap_100Dx100D_OneTwo_m6-crop.pdf} };
    \spy [red] on (-.2,-.2)  in node [left] at (1.5,1);
\end{tikzpicture}
&\begin{tikzpicture}[spy using outlines={circle, magnification=4, size=2cm, connect spies}]
    \node { \includegraphics[width=4cm]{figures_siggraph25/heatmaps/heatmap_100Dx100D_OneTwo_m8-crop.pdf} };
    \spy [red] on (-.2,-.2)  in node [left] at (1.5,1);
\end{tikzpicture}
&\begin{tikzpicture}[spy using outlines={circle, magnification=4, size=2cm, connect spies}]
    \node { \includegraphics[width=4cm]{figures_siggraph25/heatmaps/heatmap_100Dx100D_OneTwo_m12-crop.pdf} };
    \spy [red] on (-.2,-.2)  in node [left] at (1.5,1);
\end{tikzpicture}
&
\begin{tikzpicture}[spy using outlines={circle, magnification=4, size=2cm, connect spies}]
    \node { \includegraphics[width=4cm]{figures_siggraph25/heatmaps/heatmap_OneTwo_max-crop.pdf} };
    \spy [red] on (-.2,-.2)  in node [left] at (1.5,1);
\end{tikzpicture}\\
 \raisebox{.3cm}{\rotatebox{90}{$|i-j|\leq 4$}}& \includegraphics[width=4cm]{figures_siggraph25/heatmaps/hist_OneTwo_m6.pdf}&\includegraphics[width=4cm]{figures_siggraph25/heatmaps/hist_OneTwo_m8.pdf} &\includegraphics[width=4cm]{figures_siggraph25/heatmaps/hist_OneTwo_m12.pdf} &  \includegraphics[width=4cm]{figures_siggraph25/heatmaps/hist_max.pdf}\\
 \raisebox{.5cm}{\rotatebox{90}{All pairs}} &\includegraphics[width=4cm]{figures_siggraph25/heatmaps/hist_OneTwo_m6-all.pdf}&\includegraphics[width=4cm]{figures_siggraph25/heatmaps/hist_OneTwo_m8-all.pdf}&\includegraphics[width=4cm]{figures_siggraph25/heatmaps/hist_OneTwo_m12-all.pdf}&\includegraphics[width=4cm]{figures_siggraph25/heatmaps/hist_max-all.pdf}
\end{tabular}
\caption{Up to 100 dimensions, we show the experimental $t$ values for each  2D projection pair $(i,j)$ of Sobol' sequences with Joe and Kuo (top) and our $(1,2)$-sequences (bottom) for $m=6$, $m=8$, and $m=12$. The last column corresponds to the maximum $t$ values over $m\in\{2,..,12\}$. The histograms  highlight the distributions of $t$ values for close pairs (\emph{i.e.}, $|i-j|\leq 4$) and all pairs. While most pairs have comparable $t$ values with only a small degradation far from the diagonal, our construction shows a significant improvement for close consecutive pairs (with $t\leq 1$ by construction for pairs $(2i,2i+1)$).\label{fig:heatmaps}}
\end{figure*}


\begin{algorithm}
    \caption{Constructing  $\mathcal{K}_e$.}
    \KwResult{A set of candidate characteristic matrices $\mathcal{K}_e$.}
    \While{true}{
        Draw a random upper triangular matrix  $A$ with $1$ on the diagonal and a random square matrix $B$, both of size $e \times e$ \;
        Draw a random triangular matrix $C$ of size $e \times e$ satisfying the property $\mathcal{Q}$\;
        Construct $K^{(2)}$ of size $4e \times 4e$ using Eq.~(\ref{eq:KiKiplusone})\;
        \If{ $K^{(2)}$ satisfies the property $\mathcal{P}$}{
            Append $K^{(1)}$ to $\mathcal{K}_e$\;
        }
 	\label{algo:CharacteristicMatrices}
   }
\end{algorithm}

Then, each characteristic matrix in $\mathcal{K}_e$ is used to define two Sobol' matrices for each pair of irreducible polynomials $(p,p^2+p+1)$ of degrees $e$ and $2e$ respectively (Theorem~\ref{theorem:construction}). We construct the initialization matrices $D^{\be}_p$ and $D^{\btwoe}_{p^2+p+1}$ as follows:  we draw a random non-singular upper triangular matrix $D^\be_p$ of size $e \times e$, and expand it to $D^\btwoe_{p^2}$ using standard Sobol' iterations for polynomial $p$, and construct $D^{\btwoe}_{p^2+p+1}$ using $D^{\btwoe}_{p^2+p+1} = K^\btwoe D^\btwoe_{p^2}$ where the characteristic matrix $K^\btwoe$ is drawn from $\mathcal{K}_e$ (see Algorithm~\ref{algo:ConstructingDxDy}).


Finally, we convert the initialization matrices $D^{\be}_p$ and $D^{\btwoe}_{p^2+p+1}$ into a set of direction vectors for Sobol' construction, which is compatible with the format of Joe and Kuo~\shortcite{joe2008constructing}.



\begin{algorithm}
    \caption{Constructing many $(1,2)-$sequence initialization matrices} \label{alg:diffowen}
    \KwData{a degree $e$, a set of candidate characteristic matrices $\mathcal{K}_e$.}
    \KwResult{A collection of tuples $\left\{\left(p, D^{\be}_p and D^{\btwoe}_{p^2+p+1}\right)\right\}$}
    \While{true}{
    \ForAll{pairs of irreducible polynomials $p$ and $p^2+p+1$ of degrees $e$ and $2e$ respectively}{
       Create random non-singular upper triangular matrix $D^\be_p$  of size $e \times e$ \;
       Expand $D^\be_p$ to $D^\btwoe_{p^2}$ using Sobol' construction with $p$\; 
       Draw a characteristic matrix $K^\btwoe$ from $\mathcal{K}_e$\;
       Compute $D^\btwoe_{p^2+p+1}  = K^\btwoe  D^\btwoe_{p^2}$\;
        Append $\left(p, D^\be_p, D^\btwoe_{p^2+p+1}\right)$ to the result\;
    }
    \label{algo:ConstructingDxDy}
   }
\end{algorithm}


\subsection{Further improvements}
\label{sec:FurtherImprovements}

For each pair of polynomials $(p_{i}, p_{i+1})$ we can generate a large number of possible initializations, as outlined in Algorithm~\ref{algo:ConstructingDxDy}, which all satisfy our conditions for generating $(1, 2)$-sequences. 
Consequently, we enforce additional criteria to enhance our multi-dimensional construction. 
In the context of computer graphics, we aim to achieve higher quality not only for consecutive pairs of dimensions but also for 4-tuples of dimensions, which group consecutive pairs. 
We select only solutions with guaranteed reasonably-good $t \leq 4$ for 4D projections up to $2^{15}$ points.
We further seek to achieve low $t$ for dimensions that are close to $(i, i+1)$.
Specifically, pairs of dimensions are progressively added by proposing
pairs of matrices generated from characteristic matrices. Accepting a new pair of matrices requires that, within the 6D block of dimensions involving the last 4D block and the new pair, all-pairs 2D values of $t\leq 3$ for any $m\leq 8$. Further, for the 4D block involving the last pair of dimensions and the new pair, the 4D value of $t \leq 4$ for $m \leq 15$ and $t \leq 3$ for $m \leq 10$.
Pairs of polynomials of degree lower than $e=12$ (involving the first 36 dimensions) were further inspected manually to ensure high quality. 
This optimization process is inspired by the pioneering works of Joe and Kuo~\shortcite{joe2008constructing} and Faure and Lemieux~\shortcite{faure2019implementation}.
It is also close to the optimization described by Ostromoukhov et al.~\shortcite{Quad2024}.
Visualization of 2D projections for our resulting sequence can be seen in Fig.~\ref{fig:2dpatches} while discrepancy and integration errors for 2D and 4D projections can be seen in Fig.~\ref{fig:discrepanceAll2D4D}. In Figure \ref{fig:heatmaps}, we further analyze the experimental $t$ values any 2D projections, for various sample counts, up to 100D. While our construction provides better $t$ values for nearly consecutive pairs (see histograms), the experimental $t$ values for distant polynomials are only slightly worse than Sobol's. 

It is important to note that, aside from the optimization criteria, our construction behaves like any other Sobol' construction. 
Specifically, some remote pairs of dimensions or n-tuples beyond the optimized 4-tuples mentioned earlier may exhibit  ``good'' or ``bad'' values of $t$, which fall outside the control of our optimization process. 
This limitation is also present in the aforementioned optimizations~\cite{joe2008constructing,faure2019implementation,Quad2024}.

For dimensions greater than 692, the standard Joe and Kuo initializations can be used, provided they do not reuse our optimized polynomials. To assist with this, we provide a complementary initialization table for reference, along with the corresponding initialization matrices, integrating Joe and Kuo's dimensions for dimensions greater than 692 that excludes our polynomials.


\begin{figure*}[!h]
    \begin{tabular}{c|c}
      &  \quad\quad \quad\quad Discrepancy \quad\quad\quad\quad\quad\quad\quad\quad\quad MSE Integration (Gaussians) \quad\quad\quad\quad\quad\quad MSE Integration (Heavisides)\\
  &\includegraphics[width=.95\linewidth]{figures_siggraph25/GraphsDiscrepancyAndIntegration/2D_Discrepancy+MSEGauss+MSEHeaviside_dims_0_1.pdf} \\
  \raisebox{3cm}{\rotatebox{90}{2D}}&\includegraphics[width=.95\linewidth]{figures_siggraph25/GraphsDiscrepancyAndIntegration/2D_Discrepancy+MSEGauss+MSEHeaviside_dims_14_15.pdf}\\  
  \hline
  &\includegraphics[width=.95\linewidth]{figures_siggraph25/GraphsDiscrepancyAndIntegration/4D_Discrepancy+MSEGauss+MSEHeaviside_dims_0_1_2_3.pdf} \\
  \raisebox{3cm}{\rotatebox{90}{4D}}&\includegraphics[width=.95\linewidth]{figures_siggraph25/GraphsDiscrepancyAndIntegration/4D_Discrepancy+MSEGauss+MSEHeaviside_dims_12_13_14_15.pdf}  
    \end{tabular}
      \caption{In 2D and 4D, we evaluate the samplers quality with respect to the generalized $L_2$ discrepancy measure \cite{hickernell1998generalized} and integration errors (MSE) for random Gaussians and random Heavisides integrands (results averaged over 64 Owen-scrambled point sets). 
      Although Sobol'~\shortcite{Sobol1967distribution}/Joe and Kuo~\shortcite{joe2008constructing} and Faure and Lemieux \shortcite{faure2016irreducible} sequences are of high quality for the pair $(0,1)$ and the quadruple $(0,1,2,3)$, higher discrepancies and integration errors can be observed for the pair $(14,15)$ and the quadruple $(12,13,14,15)$. In contrast, quad-optimized LDS in $GF(3)$~\cite{Quad2024} and our sequences show comparable results, with our sequences more easily computed in $GF(2)$.}
  \label{fig:2ddiscr}
  \end{figure*}

% VO: typical nTrials per npts {npts,nTrials}: {{16,1024},{32,1024},{64,1024},{128,512},{256,256},{512,128},{1024,64},{2048,64},{4096,64},{8192,64},{16384,64},{32768,64},{65536,64}} 

\begin{figure*} \setlength{\tabcolsep}{0pt}
    \begin{tabular}{cc}
    \includegraphics[width=8cm]{figures_siggraph25/superpos/discrepance-all-2D.pdf}&
    \includegraphics[width=8cm]{figures_siggraph25/superpos/discrepance-all-4D.pdf}\\
    2D & 4D
    \end{tabular}
    \caption{Generalized $L_2$ discrepancy \cite{hickernell1998generalized} of consecutive 2D pairs (left) and quadruples of dimensions (right) of the first 36 dimensions of Sobol' using tables of Joe and Kuo~\shortcite{joe2008constructing} (red), 
    quad-optimized projection in $GF(3)$~\cite{Quad2024} (blue), Faure and Lemieux \shortcite{faure2016irreducible} (magenta), and our sequences (green). We observe comparable results to the quad-optimized projection in $GF(3)$ while staying in $GF(2)$, both improving over Joe and Kuo.}
\label{fig:discrepanceAll2D4D}
\end{figure*}

\begin{figure}
    \begin{overpic}[width=8cm]{figures_siggraph25/new-graphs_2D+4D+8D+16D+32D.pdf}
        \put(91,9) {2D}
        \put(91,15) {4D}
        \put(91,25) {8D}
        \put(89,35) {16D}
        \put(89,48) {32D}
    \end{overpic}
    \caption{For $(t,s)$-sequences only, we compare their generalized $L_2$ discrepancy in higher dimensions  (from 2D to 32D, by increasing order of the polynomials). We observe similar results for all LDS sequences, while our sequence has highly uniform 2D projections  (see Fig.~\ref{fig:discrepanceAll2D4D}-left). ZeroTwo~\cite{pharr2023physically} and Padded 4D~\cite{burley2020practical} are not LDS in higher dimensions and thus do not offer the same convergence rate. } \label{fig:alldisc}
\end{figure}

\subsection{Renderings}

We  evaluate our sequence with PBRT-v4~\cite{pharr2023physically} used as a per-pixel path tracer. 
PBRT constructs paths by combining 1D and 2D random variables. 
When sampling 1D variables, we sample 2 of our dimensions and keep one of them cached for the next 1D variable.
Constructing a path involves sampling a pixel (2D), time (1D) and the lens (2D). 
Evaluating direct lighting additionally requires selecting the light source (1D) and a point on that light source (2D). 
Evaluating one bounce of indirect lighting requires selecting the material (1D) and sampling a direction from it (2D). 
In this setting, rendering with direct lighting uses 11 dimensions involving 6 optimized pairs, % from 6 sequences, 
while rendering with one bounce of indirect lighting requires 17 dimensions (9 pairs), two bounces of indirect lighting require 23 dimensions (12 pairs), etc. We did not use Russian roulette nor spectrum sampling. 
We compare our results to those of other samplers in Fig.~\ref{fig:rendering}, focusing on rendering error. We used Owen scrambling for all methods.

\begin{figure*}[!tbh]
    \begin{tabular}{rl}
		\raisebox{.2cm}{\begin{overpic}[width=6.7cm]{results/rendering/cornell/cornell}
		\end{overpic}}
	&
		\begin{overpic}[width=7cm]{results/rendering/cornell/mse}
            \put(-1,26) {\small\rotatebox{90}{MAE}}
        \end{overpic}
	\\
    \raisebox{.2cm}{\begin{overpic}[width=6.7cm]{results/rendering/dining-room/dining-room}
		\end{overpic}}
	&
		\begin{overpic}[width=7cm]{results/rendering/dining-room/mse}
            \put(-1,26) {\small\rotatebox{90}{MAE}}
		\end{overpic}\\
		\raisebox{.2cm}{\begin{overpic}[width=6.7cm]{results/rendering/spaceship/spaceship}
		\end{overpic}}
	&
		\begin{overpic}[width=7cm]{results/rendering/spaceship/mse}
            \put(30,-3) {\small{Number of samples per pixel}}
            \put(-1,26) {\small\rotatebox{90}{MAE}}
        \end{overpic}
    \end{tabular}
    \caption{Rendering results in 17D with 2 bounces (top and middle), and 35D with 5 bounces (bottom). We provide convergence graphs as a function of the number of samples per pixel (mean absolute error -- MAE) showing that rendering can benefit from high-quality projections ($t=1$ in our case) while being sequence for progressive rendering, contrary to Cascaded Sobol'~\cite{paulin2021cascaded} (additional results in supplementary material). \textit{Breakfast Room by blendswap user Wig42 and Spaceship by thecali, compiled by Benedikt Bitterli}.
	\label{fig:rendering}}
\end{figure*}


Our sequences with guaranteed $t=1$ 2D projections perform similarly to the base-3 progressive point sets of quad-optimized $GF(3)$ \cite{Quad2024} and the base-2 point sets of Cascaded Sobol'~\cite{paulin2021cascaded}. 
This result is in agreement with other discrepancy and integration results in Fig.~\ref{fig:discrepanceAll2D4D} and Fig.~\ref{fig:2ddiscr}. Padding 4D Sobol' samples with random shuffling~\cite{burley2020practical}  yields better results than padding in 2D (ZeroTwo~\cite{pharr2023physically}). While our high-dimensional behavior is guaranteed low-discrepancy and padded 4D Sobol' has poor discrepancy convergence (see Fig.~\ref{fig:alldisc}), our renderings remain similar in most cases.

Working with $GF(2)$ arithmetic is also faster than $GF(3)$. 
Additions in $GF(2)$ can be computed with a binary \textit{xor} in parallel on 32 values whereas $GF(3)$ requires modulo arithmetic and tabulated operations on scalar values~\cite{Quad2024}. 
Generating 8D points is roughly four time slower with quad-optimized $GF(3)$ (798ms vs 201ms respectively, for 16M samples, on a Ryzen 3900X). In our tests, when rendering a Cornell box at 256spp at 1k resolution, sampling (in base-2) takes at least 75\% of the total render time in PBRT, while the more complex SanMiguel scene results in 15\% of the time spent in sampling.
Easy-to-use precompiled matrices and fast point generation functions are available in the supplementary materials, as well as the modified PBRT source code.

\section{Conclusions}

We designed a theoretical construction of 2D Sobol' sequences with $t=1$ using $p$ and $p^2+p+1$, while remaining low-discrepancy in higher dimensions.
In practice, we found many 
solutions of unique characteristic matrices, in contrast to the unique solution for $t=0$. We used 346 such pairs to produce a 692D sequence having at most $t=1$ in 2D consecutive projections. 
In the process of proving $t=1$, we discovered a new recursive construction for Sobol' matrices and for characteristic matrices.
However, the availability of pairs of irreducible polynomials in the form $p$ and $p^2+p+1$ is limited, and their degrees  quickly increase. 
In practice, we use polynomials of up to degree $2e = 32$ to produce 692 dimensions, while the construction of Faure and Lemieux~\shortcite{faure2016irreducible} uses at maximum polynomials of degree 13 to produce 1377 dimensions. 
While low-degree polynomials may appear desirable since they are guaranteed to reduce $t$ for high-dimensional integration problems, as $t$ is bounded by sums of polynomial degrees, 
this does not mean that $t$ is necessarily large when the degree is large (in fact, our solution could lead to $t=1$ in 2D for arbitrarily large polynomial degrees). 
We have found that the quality of our sequence remains competitive  for moderately high-dimensional integration problems arising in path tracing, despite our use of higher-degree polynomials. 
Our use of a base-2 construction remains an advantage in rendering where efficiency is critical, and base-2 allows for both efficient sampling and Owen scrambling~\cite{owen1995randomly,burley2020practical}. 
Our sampler produces a sequence, which is ideal for progressive rendering. Our use of standard Sobol' construction makes integration into existing renderers already supporting Sobol' extremely lightweight.
We nevertheless intend to explore $b>2$ within our framework to discover $t=0$ sequences in higher dimensions, which remains a gold standard for numerical integration.


\begin{acks}
    This work was partially funded by ANR-20-CE45-0025 (MoCaMed), by ANR-22-CE46-000 (StableProxies), and donations from Adobe Inc.
\end{acks}

% Bibliography
\bibliographystyle{ACM-Reference-Format}
\bibliography{references}    


% Appendix
\appendix
\section{Additional derivations}
\label{app:details}



\subsection{Proofs of Eq.~(\ref{eq:FrinvFp})}
\label{app:proofsR}

Starting from eq.~\ref{eq:F} and $F_p^{-1} = ({Id}_{\small e} + R_{p,e} ) \dots ({Id}_{\small e} + R_{p,2} )$, we have:
\begin{align*}
&F_{p^2+p+1} F_{p^2}^{-1}=\\ 
&({Id}_{\small 2e} + R_{p^2+p+1,2} ) \dots ({Id}_{\small 2e} + R_{p^2+p+1,2e} )({Id}_{\small 2e} + R_{p^2,2e} ) \dots ({Id}_{\small 2e} + R_{p^2,2} )
\end{align*}

To simplify the notations, we  denotes $R'_{k}=({Id}_{\small 2e} + R_{p^2+p+1,k})$  and $R''_{k}=({Id}_{\small 2e} + R_{p^2,k})$. Hence, we have:

\begin{align}
&F_{p^2+p+1}F_{p^2}^{-1}=  \nonumber\\
&\quad \underbrace{R'_{2} \dots  R'_{e-1}}_{(i)}\underbrace{R'_{e} \dots  R'_{2e-1}R'_{2e} R''_{2e} R''_{2e-1}\dots R''_{e}}_{(ii)}\underbrace{R''_{e-1}\dots R''_{2}}_{(iii)}\label{eq:appcases}
\end{align}
In the following, we will use this illustration for $R'_{k}$ (the column of index $k$ contains the $(k-1)$ highest degree coefficients of $p^2+p+1$):
\begin{center}
\begin{overpic}[width=2cm]{figures_siggraph25/blockmatbis.pdf}
 \put(30,60){$Id$}
 \put(80,10){$Id$}
 \put(82,60){0}
 \put(61,49){\rotatebox{90}{$p^2+p$}}
 \put(30,12){0}
 \put(65,12){0}
 \put(65,27.5){1}
\end{overpic}
\end{center}
Note that by definition of $R'_k$ and $R''_k$ matrices, we can only consider polynomials  $p^2+p$ and $p^2$ respectively, as the constant factor is dropped by construction. 

Let us first consider the first innermost product in part $(ii)$ of Eq.~(\ref{eq:appcases}) involving the $p^2+p$ and $p^2$ polynomials $R'_{2e} R''_{2e}$:
\begin{center}
\begin{overpic}[width=2cm]{figures_siggraph25/blockmat.pdf}
 \put(40,60){$Id$}
 \put(84,40){\rotatebox{90}{$p^2+p$}}
 \put(90,5){1}
 \put(40,5){0}
\end{overpic}
\begin{overpic}[width=2cm]{figures_siggraph25/blockmat.pdf}
 \put(40,60){$Id$}
 \put(84,50){\rotatebox{90}{$p^2$}}
 \put(90,5){1}
 \put(40,5){0}
\end{overpic}
\raisebox{1cm}{=} \begin{overpic}[width=2cm]{figures_siggraph25/blockmatpplusone.pdf}
 \put(40,60){$Id$}
 \put(87,68){\rotatebox{90}{$p$}}
 \put(89,30){0}
 \put(90,5){1}
 \put(40,5){0}
 \put(103,68){$e$}
\end{overpic}
\end{center}
as $p^2+p^2$ coefficients cancel out for rows greater or equal to $e$. We denote by $U_1$ the resulting matrix. Let us now consider the product  $U_2 = R'_{2e-1}\,U_1\,R''_{2e-1}$:
\begin{center}
\begin{overpic}[width=2cm]{figures_siggraph25/blockmater.pdf}
 \put(40,60){$Id$}
 \put(74,38){\rotatebox{90}{$p^2+p$}}
 \put(90,2){1}
 \put(79,13){1}
 \put(79,2){0}
 \put(90,13){0}
  \put(40,2){0}
 \put(40,2){0}
 \put(90,50){0}
\end{overpic}
\begin{overpic}[width=2cm]{figures_siggraph25/blockmatpplusone.pdf}
 \put(40,60){$Id$}
 \put(87,72){\rotatebox{90}{$p$}}
 \put(89,30){0}
 \put(90,5){1}
 \put(40,5){0}
\end{overpic}
\begin{overpic}[width=2cm]{figures_siggraph25/blockmater.pdf}
 \put(40,60){$Id$}
 \put(74,50){\rotatebox{90}{$p^2$}}
  \put(40,13){0}
    \put(40,2){0}
  \put(90,2){1}
 \put(79,13){1}
 \put(79,2){0}
 \put(90,13){0}
   \put(90,50){0}
\end{overpic}
\end{center}
which simplifies to
\begin{center}
\begin{overpic}[width=2cm]{figures_siggraph25/blockmater2.pdf}
 \put(40,60){$Id$}
  \put(87,64){\rotatebox{90}{$p$}}
 \put(72,72){\rotatebox{90}{$p$}}
 \put(90,2){1}
 \put(75,18){1}
 \put(90,18){0}
 \put(90,38){0}
 \put(75,38){0}
 \put(75,2){0}
  \put(40,2){0}
  \put(40,18){0}
\end{overpic}
\end{center}
If we repeat this process for all triplets of matrices $R'_{k}\, U_{2e-k}\,R''_k$ for the $k$ indices of $(ii)$, we end up with the matrix $U_e$:
\begin{center}
\begin{overpic}[width=2cm]{figures_siggraph25/blockmatU.pdf}
 \put(20,70){$Id$}
 \put(65,22){$Id$}
  \put(65,70){$F^{-1}_p$}
 \put(22,22){0}
\end{overpic}
\end{center}
Indeed, for each product $R'_{k}\, U_{2e-k}\,R''_k$, all $p^2$ coefficents vanish, leading to a triangular upper-right block  with shifted $p$ coefficients as in Eq.~(\ref{eq:Q}).

Let us now consider the product between $(ii)$ and the $(i)$ and $(iii)$ parts in Eq.~(\ref{eq:appcases}).
First, we observe that
\begin{center}
\raisebox{.9cm}{$R''_{e-1}R''_{e-2}=$}~~~~ \begin{overpic}[width=2cm]{figures_siggraph25/blockmati.pdf}
 \put(20,70){$Id$}
 \put(65,22){$Id$}
  \put(75,70){0}
 \put(22,22){0}
 \put(40,68){\rotatebox{90}{\small $p^2$}}
 \put(35,85){\rotatebox{90}{\small $p^2$}}
  \put(48,51){\tiny 1}
  \put(40,51){\tiny 0}
  \put(40,57){\tiny 1} 
\end{overpic}
\end{center}
By doing such products for all matrices of $(iii)$, we obtain an upper-left block which corresponds to the first $e\times e$ entries of $ F^{-1}_{p^2+p}$.
$$
R''_{e-1}\ldots R''_2 =  \begin{pmatrix}
    F^{-1}_{p^2} & 0  \\
    0 & Id_e
  \end{pmatrix}\,.
$$
For products in $(i)$, we use the fact that 
$$
\left(R'_{2}\ldots {R'}_{e-1}\right)^{-1} = {R'}_{e-1}^{-1}\ldots {R'}_{2}^{-1} =  R'_{e-1}\ldots R'_{2} =\begin{pmatrix}
    F^{-1}_{p^2+p} & 0  \\
    0 & Id_e
  \end{pmatrix}\,,
$$
as $R'_{k}$ is its own inverse and using a similar construction as for $(iii)$. Thus, using the inverse of a block matrix, we obtain
$$
R'_{2}\ldots {R'}_{e-1} = \begin{pmatrix}
    F_{p^2+p} & 0  \\
    0 & Id_e
  \end{pmatrix}\,,$$
  which equals to $\begin{pmatrix}
    F_{p^2} & 0  \\
    0 & Id_e
  \end{pmatrix}$ as no  coefficients for the $p$ term are present in the upper-left block.

  We finally have
\begin{align*}
  F_{p^2+p+1}F_{p^2}^{-1}&= \begin{pmatrix}
    F_{p^2} & 0  \\
    0 & Id_e
  \end{pmatrix}
  \begin{pmatrix}
    Id_e & F^{-1}_p  \\
    0 & Id_e
  \end{pmatrix}
  \begin{pmatrix}
    F^{-1}_{p^2} & 0  \\
    0 & Id_e
  \end{pmatrix}\\
  &=\begin{pmatrix}
    F_{p^2} & 0  \\
    0 & Id_e
  \end{pmatrix}
  \begin{pmatrix}
    F^{-1}_{p^2}  & F^{-1}_p  \\
    0 & Id_e
  \end{pmatrix}\\
  &= \begin{pmatrix}
    Id_{\small e} & F_p  \\
    0 & Id_{\small e} 
\end{pmatrix}.
\end{align*}
which concludes Eq.~(\ref{eq:FrinvFp}). The last step uses the observation that, for the upper-right block,  $F_{p^2}F_p^{-1} = \left(\left(F_{p^2} F^{-1}_{p}\right)^{-1}\right)^{-1} = \left(F_p F_{p^2}^{-1}\right)^{-1} = \left(F_p F_{p}^{-1}F_{p}^{-1}\right)^{-1} = F_p$.

\subsection{Proofs of Eq.~(\ref{eq:FrinvFpDp2})}
\label{app:proofs1}

First, combining Eq.~(\ref{eq:invDp2}) and (\ref{eq:FrinvFp}), we have: 
\begin{align}
F_{p^2+p+1} F_{p^2}^{-1} D_{p^2}^{-1} &=\begin{pmatrix}
    Id_{\small e} & F_p  \\
    0 & Id_{\small e} 
\end{pmatrix}\begin{pmatrix}
    Id_{\small e} & Q_p  \\
    0 & F_p^{-1} 
\end{pmatrix}\begin{pmatrix}
    D_p^{-1} & 0 \\
    0 & D_p^{-1} 
\end{pmatrix} \\
&= \begin{pmatrix}
    Id_{\small e} & Q_p + Id_{\small e}  \\
    0 & F_p^{-1} 
\end{pmatrix} \begin{pmatrix}
    D_p^{-1} & 0 \\
    0 & D_p^{-1} 
\end{pmatrix} \label{eq:Gx} \\
&= \begin{pmatrix}
    Id_{\small e} & Q_p   \\
    0 & F_p^{-1} 
\end{pmatrix} \begin{pmatrix}
    D_p^{-1} & D_p^{-1} \\
    0 & D_p^{-1} 
\end{pmatrix} \\
&=\begin{pmatrix}
    Id_{\small e} & Q_p  \\
    0 & Id_{\small e} 
\end{pmatrix}\begin{pmatrix}
    D_p^{-1}  & 0  \\
    0 & D_p^{-1}  
\end{pmatrix}\begin{pmatrix}
    Id_{\small e} & Id_{\small e} \\
    0 & Id_{\small e} 
\end{pmatrix} \\
 &= D_{p^2}^{-1} \begin{pmatrix}
    Id_{\small e} & Id_{\small e}  \\
    0 & Id_{\small e} 
\end{pmatrix},
\end{align}
leading to Eq.~(\ref{eq:FrinvFpDp2}). Starting from Eq.~(\ref{eq:Gx}), we also have
\begin{align}
F_{p^2+p+1} F_{p^2}^{-1} D_{p^2}^{-1} &= \begin{pmatrix}
    Id_{\small e} & Q_p + Id_{\small e}  \\
    0 & F_p^{-1} 
\end{pmatrix} \begin{pmatrix}
    D_p^{-1} & 0 \\
    0 & D_p^{-1} 
\end{pmatrix} \\
&= \begin{pmatrix}
    Id_{\small e} & Q_{p+1}  \\
    0 & F_p^{-1} 
\end{pmatrix} \begin{pmatrix}
    D_p^{-1} & 0 \\
    0 & D_p^{-1} 
\end{pmatrix} \label{eq:GxBis}\,,
\end{align}
that will be used later.



\subsection{Proofs of Eq. (\ref{eq:Qp2plus})}
\label{app:proofs2}

Let us now prove the following statement:
$$
(Q_{p^2} + Q_{p^2+p+1} F_{p^2+p+1} F_{p^2}^{-1}) D_{p^2}^{-1} = D_{p^2}^{-1}  \begin{pmatrix}
    Id_{\small e} & Id_{\small e}  \\
    Id_{\small e} & 0 
\end{pmatrix}\,.
$$

From now on, we make explicit the size of the matrices using $[e]$ or $[2e]$ superscripts. First, by definition, $Q^\btwoe_{p^2}$ is
\begin{equation}
    Q^\btwoe_{p^2} = 
    \begin{pmatrix}
        a_{0} & 0 & 0 & \dots & 0 \\
        0 & a_{0} & 0 & \dots  & 0 \\
        a_{1} & 0 &a_{0} & \dots  & 0 \\
        0 & a_1 & 0 & \dots  & 0 \\
        \vdots & \vdots & \vdots & \ddots & \vdots\\ 
        a_{e-1} & 0 & a_{e-2} & \dots & a_0 
    \end{pmatrix}\,.
\end{equation}
We decompose $Q^\btwoe_{p^2}$ into $e\times e$ blocks:
\begin{equation}
 Q^\btwoe_{p^2} =  \begin{pmatrix}
        Q^\be_{p^2} & 0 \\
        \widetilde{Q}^\be_{p^2} & Q^\be_{p^2}
    \end{pmatrix}.\label{eq:decQtwoe}
\end{equation}




Similar to $Q^\btwoe_{p^2}$, $\widetilde{Q}^\be_{p^2}$  is also Toeplitz.
Furthermore, we have 
$$
Q^\be_{p+q} = Q^\be_p + Q^\be_q\quad\text{and}\quad Q^\be_{pq} = Q^\be_p Q^\be_q\,,
$$
for any polynomial $p$ and $q$ of  degree $e$. The same holds for $\widetilde{Q}^\be_{p+q}$ and $\widetilde{Q}^\be_{pq}$ matrices.
Now,
%
%
\begin{align}
  (Q_{p^2} &+ Q_{p^2+p+1} F_{p^2+p+1} F^{ -1}_{p^2}) D_{p^2}^{-1} =\nonumber\\
 &=Q_{p^2}  D_{p^2}^{-1} + Q_{p^2+p+1} F_{p^2+p+1} F^{ -1}_{p^2}D_{p^2}^{-1}\nonumber\\
 \text{using Eq.~(\ref{eq:invDp2})}& \text{ and \ref{eq:GxBis} with $T=\begin{pmatrix}
    D_p^{-1} & 0 \\
    0 & D_p^{-1} 
 \end{pmatrix}$ }\nonumber\\
 &=\left(\begin{pmatrix}
 Q^\be_{p^2} & 0 \\
 \widetilde{Q}^\be_{p^2} & Q^\be_{p^2}
\end{pmatrix}\begin{pmatrix}
     Id_{\small e} & Q^{\be}_{p}  
     \\
     0 & \widetilde{Q}^{\be}_{p}
 \end{pmatrix} \right.\nonumber\\
&\left.\quad + \begin{pmatrix}
    Q^\be_{p^2+p+1} & 0 \\
    \widetilde{Q}^\be_{p^2+p+1} & Q^\be_{p^2+p+1}
   \end{pmatrix}\begin{pmatrix}
    Id_{\small e} & Q^{\be}_{p+1}  \\
    0 & \widetilde{Q}^{\be}_{p+1}
\end{pmatrix}\right)T\,. \label{eq:Tfactor}
\end{align}
First, we observe that  $\widetilde{Q}^{\be}_{p} = \widetilde{Q}^{\be}_{p+1}$. The first factor can be rewritten
\begin{align*}   
    &\begin{pmatrix}
        Q^\be_{p^2} &  Q^\be_{p^3} \\
        \widetilde{Q}^\be_{p^2} & \widetilde{Q}^\be_{p^2}Q^\be_{p}+Q^\be_{p^2}\widetilde{Q}^\be_{p}
       \end{pmatrix}\\
&\quad\quad+  
    \begin{pmatrix}
        Q^\be_{p^2+p+1} &  Q^\be_{p^3+1}  \\
        \widetilde{Q}^\be_{p^2+p+1} & \widetilde{Q}^\be_{p^2+p+1}Q^\be_{p+1}+{Q}^\be_{p^2+p+1}\widetilde{Q}^{\be}_{p+1}
       \end{pmatrix}\,,
\end{align*}
since $Q^\be_{p^2}Q^\be_{p}=Q^\be_{p^3}$ and $Q^\be_{p^2+p+1}Q^\be_{p}=Q^\be_{p^3+p^2+p}$. Furthermore, for any polynomial $p$ and $q$ of  degree $e$, we have
\begin{align*}
Q^\btwoe_{pq} &= Q^\btwoe_{p}Q^\btwoe_{q}\\
&= \begin{pmatrix}
    Q^\be_{p} & 0 \\
    \widetilde{Q}^\be_{p} & Q^\be_{p}
\end{pmatrix}\begin{pmatrix}
    Q^\be_{q} & 0 \\
    \widetilde{Q}^\be_{q} & Q^\be_{q}
\end{pmatrix}\\
&= \begin{pmatrix}
    {Q}^\be_{pq} &0 \\
    \widetilde{Q}^\be_{p}{Q}^\be_{q} + {Q}^\be_{p}\widetilde{Q}^\be_{q}&{Q}^\be_{pq}
\end{pmatrix}\,.
\end{align*}


Hence, the first factor of Eq.~(\ref{eq:Tfactor}) is 
\begin{align*}   
    \begin{pmatrix}
        Q^\be_{p^2} &  Q^\be_{p^3} \\
        \widetilde{Q}^\be_{p^2} & \widetilde{Q}^\be_{p^3}
       \end{pmatrix}
+  
    \begin{pmatrix}
        Q^\be_{p^2+p+1} &  Q^\be_{p^3+1}  \\
        \widetilde{Q}^\be_{p^2+p+1} & \widetilde{Q}^\be_{p^3+1}
       \end{pmatrix}
  &=  \begin{pmatrix}
        Q^\be_{p+1} &  Q^\be_{1} \\
        \widetilde{Q}^\be_{p+1} & \widetilde{Q}^\be_{p^3}+\widetilde{Q}^\be_{p^3+1}
       \end{pmatrix}\\
       &=  \begin{pmatrix}
        Q^\be_{p+1} & Id_e \\
        F_p^{-1}& 0
       \end{pmatrix}\,,
\end{align*}
using the fact that  $\widetilde{Q}^{\be}_{p}= F_p^{-1}$ from the construction of both matrices. 


Finally, 
\begin{align*}
    (Q_{p^2} + Q_{p^2+p+1}&F_{p^2+p+1} F^{ -1}_{p^2}) D_{p^2}^{-1} =\begin{pmatrix}
        Q^\be_{p+1} & Id_e \\
        F_p^{-1} & 0
       \end{pmatrix} T\\
       &=\begin{pmatrix}
        Q^\be_{p+1} & Id_e \\
        F_p^{-1} & 0
       \end{pmatrix} \begin{pmatrix}
    D_p^{-1} & 0 \\
    0 & D_p^{-1} 
 \end{pmatrix}\\
 &=\begin{pmatrix}
    Q^\be_{p} + Id_e & Id_e \\
    F_p^{-1} & 0
   \end{pmatrix} \begin{pmatrix}
D_p^{-1} & 0 \\
0 & D_p^{-1} 
\end{pmatrix}\\
&=\begin{pmatrix}
    Id_e&Q^\be_{p} + Id_e  \\
    0& F_p^{-1}
   \end{pmatrix}\begin{pmatrix}
    0 &1  \\
    1&0
   \end{pmatrix} \begin{pmatrix}
D_p^{-1} & 0 \\
0 & D_p^{-1} 
\end{pmatrix}\\
 &=D_{p^2}^{-1} \begin{pmatrix}
    Id_{\small e} & Id_{\small e}  \\
    Id_{\small e} &0 
\end{pmatrix} \,,
\end{align*}   
which concludes the proof of Eq.~(\ref{eq:Qp2plus}).

\subsection{$(1,2)-$sequences and corank 1 submatrices of $K$}
\label{app:corank}
Let us consider two Sobol' matrices $M_p$  and $M_q$ of size $m\times m$ forming a $(t,m,2)$-net. We denote $K=M_qM_p^{-1}$. First we remind that the pairs of matrices $(M_p,M_q)$ and $(Id_m, K)$ generate the same point set (up to indices permutation). Let  $\mathcal{K}_k^t$   the $(m-t)\times m$  matrix consisting of the first $k$ rows of $Id_m$ and the first $m-k-t$ rows of $K$:
\vspace{.5cm}
\begin{center}
   \raisebox{1cm}{$\mathcal{K}_k^t$ =}
\begin{overpic}[width=3cm]{figures_siggraph25/tmat.pdf}
\put(10,50) {$Id_k$}
\put(10,15) {$K'$}
\put(60,15) {$K''$}
\put(60,50) {0}
\put(10,70) {$k$}
\put(105,50) {$k$}
\put(50,70) {$m-k$}
\put(105,15) {$m-k-t$}
\end{overpic}
\end{center}

\begin{lemma}[Niederreiter \shortcite{Niederreiter1992random} (p. 73) and Paulin et al \shortcite{paulin2022generator}]
    $M_p$ and $M_q$ is a $(t,m,2)-$net if and only if for all $k\in \{1,\dots,m\}$, $\mathcal{K}_k^t$ has corank $t$.
\end{lemma}
From block-wise rank computation (the corank of a block triangular matrix with one full rank diagonal block is the corank of the other diagonal block \cite{meyer1973generalized}), we have
$$
\text{corank}(\mathcal{K}_k^t) = \text{corank}(K'')\,.
$$

Focusing on $(1,2)$-sequences, matrices $K''$ for all $m$ and all $k$ of size $(m-k-1)\times(m-k)$ are exactly the $T^{j,w}$ matrices involved in the property $\mathcal{P}$ (see Sect.~\ref{sec:ranks}). As a consequence, if each such matrices $T$ has corank 1, we can conclude that $K$ characterizes a $(1,2)-$sequence.


\end{document}
