# Direct Constructions of Multiple Shift Complementary Sets of Flexible
Lengths

Abhishek Roy  Department of Mathematics
Indian Institute of Technology Patna
Bihta, Patna, Bihar-801106, India
<EMAIL_ADDRESS>Sudhan Majhi The work of Abhishek Roy was supported by
CSIR-SRF, Government of India, File no. 09/1023(0025)/2018-EMR-I. Department
of Electrical Communication Engineering
Indian Institute of Science
Bangalore, Karnataka-560012, India
<EMAIL_ADDRESS>
###### Abstract

Golay complementary set (GCS) plays a vital role in reducing peak-to-mean
envelope power ratio (PMEPR) in orthogonal frequency division multiplexing
(OFDM). A more general version of GCS is a multiple shift complementary set
(MSCS), where by relaxing the condition of zero auto-correlation sum
throughout all the non-zero time shifts to the integer multiples of some fixed
time shift, more sequence sets can be made available. In this paper, we
propose direct constructions of MSCSs with flexible and arbitrary lengths and
flexible set sizes, by using multivariable functions, which have not been
reported before.

###### Index Terms:

Golay complementary set (GCS), orthogonal frequency division multiplexing
(OFDM), multiple shift complementary set (MSCS), multivariable function, peak-
to-mean envelope power ratio (PMEPR).

##  I Introduction

Golay complementary pair (GCP) was first conceptualized in 1951 by Marcel J.
E. Golay [1]. It is a pair of sequences which has zero aperiodic auto-
correlation function (AACF) sum for all non-zero time shifts. Because of its
ideal AACF property, it has been widely used in orthogonal frequency division
multiplexing (OFDM) [2, 3], radar [4], channel estimation [5] etc. In OFDM
system, GCP carries out the role of reducing peak-to-mean envelope power ratio
(PMEPR) [6]. The first direct construction of GCP appears in [6], where
$2^{h}$-ary ($h\geq 1$ is an integer) GCPs were constructed by generalized
Boolean functions (GBFs). The idea of GCP was extended to Golay complementary
set (GCS) which is a set of more than two sequences having zero AACF sum for
all non-zero time shifts [7]. Although there are many constructions of GCS in
the literature [8, 9, 10, 11, 12, 13], Paterson et al. first proposed a direct
construction of GCS using GBFs [14]. Like GCPs, GCSs are also used in OFDM
system to reduce PMEPR, where it is upper bounded by the number of sequences
in the GCS.

The multiple shift complementary set (MSCS) is a more general version of GCS,
where the AACF sum is zero for multiples of some fixed time shift. It was
first introduced by Xin and Fair as an alternative to GCS [15]. Later Chen et
al. provided a direct construction of MSCS using GBFs [16] with sequence
length of the form of power-of-two. In [17], the authors studied even-shift
complementary pairs, which is a special case of MSCS, where AACF sum equals
zero when the time shift is even. Also, in this case the number of constituent
sequences in the set or the set size is $2$. Recently, Chen et al. proposed a
direct construction of MSCS of non-power-of-two length by using GBFs [18]. But
the length is of the form $2^{m-1}+2^{t}$, where $m\geq 2$ and $1\leq t\leq
m-1$. To the best of authors’ knowledge, there is no direct construction of
MSCS of arbitrary lengths and set sizes.

Motivated by this, in this paper, we propose constructions of MSCSs with
flexible and arbitrary lengths by using multivariable functions. The lengths
of the proposed MSCSs are of the form $p_{1}^{m_{1}}p_{2}^{m_{2}}\dots
p_{k}^{m_{k}}$ and $p_{1}^{m_{1}}p_{2}^{m_{2}}\dots p_{k}^{m_{k}}p_{k+1}$,
where $p_{i}$’s are prime numbers and $m_{i}\geq 1$ are integers $\forall
i=1,2,\dots,k$. The set sizes of the MSCSs are of the form $p_{1}p_{2}\dots
p_{k}$ and $p_{k+1}$.

The rest of the paper is organized as follows. In Section II, some definitions
are provided. Later in Section III, the main construction of the MSCSs are
proposed. The PMEPRs of the proposed constructions are also investigated in
this section. Finally, in Section IV, concluding remarks have been given.

##  II Preliminaries

###### Definition 1

Let $\mathbf{a}=(a_{0},a_{1},\dots,a_{L-1})$ and
$\mathbf{b}=(b_{0},b_{1},\dots,b_{L-1})$ be two complex-valued sequence of
length $L$. Then the aperiodic cross-correlation function (ACCF) at time shift
$\tau$ is defined by

$\rho(\mathbf{a},\mathbf{b})(\tau)=\begin{cases}\sum_{i=0}^{L-1-\tau}a_{i}b^{*}_{i+\tau},&0\leq\tau<L;\\\
\sum_{i=0}^{L-1+\tau}a_{i-\tau}b^{*}_{i},&-L<\tau<0,\\\ \end{cases}$ (1)

where $(\cdot)^{*}$ denotes the complex conjugate. When
$\mathbf{a}=\mathbf{b}$, then it is called AACF and denoted by
$\rho(\mathbf{a})(\tau)$.

######  Definition 2 (GCS)

A set of sequences
$\\{\mathbf{a}_{0},\mathbf{a}_{1},\dots,\mathbf{a}_{M-1}\\}$ with length $L$
is called a GCS if they satisfy

$\sum_{i=0}^{M-1}\rho(\mathbf{a}_{i})(\tau)=0,\forall\tau\neq 0.$ (2)

If $M=2$, then it is called a GCP.

######  Definition 3 (MSCS)

A set of sequences
$\\{\mathbf{a}_{0},\mathbf{a}_{1},\dots,\mathbf{a}_{M-1}\\}$ with length $L$
is called a MSCS if for some positive number $S$ they satisfy

$\sum_{i=0}^{M-1}\rho(\mathbf{a}_{i})(\tau)=0,\tau\neq 0\leavevmode\nobreak\
\&\leavevmode\nobreak\ \tau\leavevmode\nobreak\ \text{mod}\leavevmode\nobreak\
S=0.$ (3)

It is denoted by $(M,L,S)$-MSCS, where $M$ is called the set size. It should
be noted that when $S=1$, then it becomes a GCS.

######  Definition 4 (Type-II ZCS)

A set of sequences
$\\{\mathbf{a}_{0},\mathbf{a}_{1},\dots,\mathbf{a}_{M-1}\\}$ with length $L$
is called a type-II Z-complementary set (ZCS) if for some positive number $Z$
they satisfy

$\sum_{i=0}^{M-1}\rho(\mathbf{a}_{i})(\tau)=0,L-Z<\left|\tau\right|<L.$ (4)

Here $Z$ is called the zero correlation zone (ZCZ) width. The set may be
denoted as type-II $(M,L,Z)$-ZCS.

###  II-A Multivariable Function

Let $\mathbb{Z}_{p}=\\{0,1,\dots,p-1\\}$ be the set of integers modulo $p$. A
multivariable function can be defined as

$f:\mathbb{Z}_{p_{1}}^{m_{1}}\times\mathbb{Z}_{p_{2}}^{m_{2}}\times\dots\times\mathbb{Z}_{p_{k}}^{m_{k}}\rightarrow\mathbb{Z}_{\lambda}$

where $p_{1},p_{2},\dots,p_{k}$ are prime numbers, $m_{i}\geq 1$, $\forall
i=1,2,\dots,k$ and $\lambda$ is a positive integer.
$v_{p_{\alpha},1},v_{p_{\alpha},2},\dots,v_{p_{\alpha},m_{\alpha}}$ are the
$m_{\alpha}$ variables which takes values from $\mathbb{Z}_{p_{\alpha}}$ for
$\alpha=1,2,\dots,k$. The set of monomials of degree at most $r$ is given by
$A^{r}=\left\\{\prod_{\alpha=1}^{k}\prod_{\beta=1}^{m_{\alpha}}v_{p_{\alpha},\beta}^{j_{\beta}^{\alpha}}:0\leq\sum_{\alpha=1}^{k}\sum_{\beta=1}^{m_{\alpha}}j_{\beta}^{\alpha}\leq
r\right\\}$. A multivariable function of order $r$ is a
$\mathbb{Z}_{\lambda}$-linear combination of the monomials from $A^{r}$. Let
$\mathbf{u}_{p_{\alpha},i_{\alpha}}$ be the vector representation of
$i_{\alpha}$ with base $p_{\alpha}$, i.e.,
$\mathbf{u}_{p_{\alpha},i_{\alpha}}=(i_{\alpha,1},i_{\alpha,2},\dots,i_{\alpha,m_{\alpha}})$
where
$i_{\alpha}=\sum_{\gamma=1}^{m_{\alpha}}p_{\alpha}^{\gamma-1}i_{\alpha,\gamma}$.
The function rule for a multivariable function is defined as
$(\mathbf{u}_{p_{1},i_{1}},\mathbf{u}_{p_{2},i_{2}},\dots,\mathbf{u}_{p_{k},i_{k}})\longmapsto
f(\mathbf{u}_{p_{1},i_{1}},\mathbf{u}_{p_{2},i_{2}},\dots,\mathbf{u}_{p_{k},i_{k}})\mod\lambda$.
One can associate a $\mathbb{Z}_{\lambda}$-valued sequence of length
$p_{1}^{m_{1}}p_{2}^{m_{2}}\dots p_{k}^{m_{k}}$ corresponding to a
multivariable function $f$ as
$\mathbf{f}=\bigg{(}f(\mathbf{u}_{p_{1},i_{1}},\mathbf{u}_{p_{2},i_{2}},\dots,\mathbf{u}_{p_{k},i_{k}}):i_{1}=0,1,\dots,p_{1}^{m_{1}}-1;\dots;i_{k}=0,1,\dots,p_{k}^{m_{k}}-1\bigg{)}$.
Also, one can associate a complex-valued sequence of length
$p_{1}^{m_{1}}p_{2}^{m_{2}}\dots p_{k}^{m_{k}}$ in a similar manner as
$\psi(\mathbf{f})=\bigg{(}\omega_{\lambda}^{f(\mathbf{u}_{p_{1},i_{1}},\mathbf{u}_{p_{2},i_{2}},\dots,\mathbf{u}_{p_{k},i_{k}})}:i_{1}=0,1,\dots,p_{1}^{m_{1}}-1;\dots;i_{k}=0,1,\dots,p_{k}^{m_{k}}-1\bigg{)}$,
where $\omega_{\lambda}=\exp(2\pi\sqrt{-1}/\lambda)$.

.

##  III Construction

In this section, we propose constructions of MSCSs by using multivariable
functions.

###### Theorem 1

Let $\pi$ be a permutation on the set $\\{s,s+1,\dots,m\\}$ for integers
$m\geq 1$ and $s\geq 1$. Let $p$ be a prime and $\lambda$ be a positive
integer such that $p\mid\lambda$. Let
$f:\mathbb{Z}_{p}^{m}\rightarrow\mathbb{Z}_{\lambda}$ be defined such that

$\begin{split}f(v_{1},v_{2},\dots,v_{m})=&\frac{\lambda}{p}\sum_{i=s}^{m-1}v_{\pi(i)}v_{\pi(i+1)}+\sum_{i=1}^{m}g_{i}v_{i}+g\\\
&+h(v_{1},v_{2},\dots,v_{s-1}),\end{split}$ (5)

where $g_{i},g\in\mathbb{Z}_{\lambda}$ and $h(v_{1},v_{2},\dots,v_{s-1})$ is
any function $h:\mathbb{Z}_{p}^{s-1}\rightarrow\mathbb{Z}_{\lambda}$. For
$s=1$, we define $h=0$. Define
$a^{\gamma}:\mathbb{Z}_{p}^{m}\rightarrow\mathbb{Z}_{\lambda}$ such that

$a^{\gamma}=f(v_{1},v_{2},\dots,v_{m})+\frac{\lambda}{p}v_{\pi(s)}\gamma.$

Then $\\{\omega_{\lambda}^{a^{\gamma}}:\gamma\in\mathbb{Z}_{p}\\}$ is a
$(p,p^{m},p^{s-1})$-MSCS.

###### Proof:

As, $\rho(\mathbf{a})(-\tau)=\rho^{*}(\mathbf{a})(\tau)$, we shall only prove
for $\tau\geq 0$. We have to show that

$\sum_{\gamma=0}^{p-1}\sum_{i=0}^{p^{m}-1-\tau}\omega_{\lambda}^{(a^{\gamma})_{i}-(a^{\gamma})_{i+\tau}}=0,$
(6)

whenever $\tau\mod p^{s-1}=0$, where
$(a^{\gamma})_{i}=a^{\gamma}(i_{1},i_{2},\dots,i_{m})$, and
$(i_{1},i_{2},\dots,i_{m})$ is the $p$-ary vector representation of $i$. Let
$j=i+\tau$, where $\tau\mod p^{s-1}=0$, i.e., $\tau$ is a multiple of
$p^{s-1}$. We have

$(a^{\gamma})_{i}-(a^{\gamma})_{j}=\left(f_{i}-f_{j}\right)+\frac{\lambda}{p}\left(i_{\pi(s)}-j_{\pi(s)}\right)\gamma.$
(7)

Now, we have two cases.

  1. Case I:

Let, $i_{\pi(s)}\neq j_{\pi(s)}$. In this case, we have

$\sum_{\gamma=0}^{p-1}\omega_{\lambda}^{\frac{\lambda}{p}(i_{\pi(s)}-j_{\pi(s)})\gamma}=\sum_{\gamma=0}^{p-1}\omega_{p}^{(i_{\pi(s)}-j_{\pi(s)})\gamma}=0,$
(8)

as $\omega_{p}^{(i_{\pi(s)}-j_{\pi(s)})\gamma}$ are the $p$-th roots of $1$
for $\gamma=0,1,\dots,p-1$. So, we have

$\sum_{\gamma=0}^{p-1}\omega_{\lambda}^{(a^{\gamma})_{i}-(a^{\gamma})_{j}}=0.$
(9)

  2. Case II:

Let, $i_{\pi(s)}=j_{\pi(s)}$. But $i\neq j$ and $j=i+\tau$, where $\tau$ is a
multiple of $p^{s-1}$, i.e., $\tau=kp^{s-1}$ for some $1\leq k\leq
p^{m-s+1}-1$. Any $k$ in this range can be written as a $(m-s+1)$-tuple vector
representation form $(k_{1},k_{2},\dots,k_{m-s+1})$ with base $p$, where
$k=\sum_{i=1}^{m-s+1}k_{i}p^{i-1}$. So,
$\tau=\sum_{i=1}^{m-s+1}k_{i}p^{i+s-2}$, which implies $i_{l}=j_{l}$ for
$l=1,2,\dots,s-1$. Now, $i\neq j$ implies that $\exists$ some
$l\in\\{s,s+1,\dots,m\\}$ such that $i_{l}\neq j_{l}$. Let $\phi$ be the
smallest number such that $i_{\pi(\phi)}\neq j_{\pi(\phi)}$. Let $i^{\delta}$
be the integer whose $p$-ary vector representation is

$(i_{1},i_{2},\dots,i_{\pi(\phi-1)}-\delta,\dots,i_{m}),\leavevmode\nobreak\
\text{if}\leavevmode\nobreak\ i_{\pi(\phi-1)}-\delta\geq 0$

and

$(i_{1},i_{2},\dots,p+i_{\pi(\phi-1)}-\delta,\dots,i_{m}),\leavevmode\nobreak\
\text{if}\leavevmode\nobreak\ i_{\pi(\phi-1)}-\delta<0,$

where $\delta\in\\{1,2,\dots,p-1\\}$ and it differs form the $p$-ary vector
representation of $i$ only in the $\pi(\phi-1)$-th position. Similarly, we
take $j^{\delta}$. Now, for $i_{\pi(\phi-1)}-\delta\geq 0$ and
$j_{\pi(\phi-1)}-\delta\geq 0$, we have

$(a^{\gamma})_{i^{\delta}}-(a^{\gamma})_{i}=-\delta\left(\frac{\lambda}{p}i_{\pi(\phi-2)}+\frac{\lambda}{p}i_{\pi(\phi)}+g_{\pi(\phi-1)}\right)$
(10)

and

$(a^{\gamma})_{j^{\delta}}-(a^{\gamma})_{j}=-\delta\left(\frac{\lambda}{p}j_{\pi(\phi-2)}+\frac{\lambda}{p}j_{\pi(\phi)}+g_{\pi(\phi-1)}\right).$
(11)

But, we have

$\begin{split}&\left((a^{\gamma})_{i^{\delta}}-(a^{\gamma})_{j^{\delta}}\right)-\left((a^{\gamma})_{i}-(a^{\gamma})_{j}\right)\\\
&=\left((a^{\gamma})_{i^{\delta}}-(a^{\gamma})_{i}\right)-\left((a^{\gamma})_{j^{\delta}}-(a^{\gamma})_{j}\right)\\\
&=-\delta\frac{\lambda}{p}\left(i_{\pi(\phi)}-j_{\pi(\phi)}\right).\end{split}$
(12)

Also, for $i_{\pi(\phi-1)}-\delta<0$ and $j_{\pi(\phi-1)}-\delta<0$, we have

$\begin{split}&\left((a^{\gamma})_{i^{\delta}}-(a^{\gamma})_{j^{\delta}}\right)-\left((a^{\gamma})_{i}-(a^{\gamma})_{j}\right)\\\
&=(p-\delta)\frac{\lambda}{p}\left(i_{\pi(\phi)}-j_{\pi(\phi)}\right).\end{split}$
(13)

But
$\omega_{\lambda}^{(p-\delta)\frac{\lambda}{p}\left(i_{\pi(\phi)}-j_{\pi(\phi)}\right)}=\omega_{\lambda}^{-\delta\frac{\lambda}{p}\left(i_{\pi(\phi)}-j_{\pi(\phi)}\right)}$.
So, considering all the possibilities, we have

$\begin{split}&\sum_{\delta=1}^{p-1}\omega_{\lambda}^{\left((a^{\gamma})_{i^{\delta}}-(a^{\gamma})_{j^{\delta}}\right)-\left((a^{\gamma})_{i}-(a^{\gamma})_{j}\right)}=\sum_{\delta=1}^{p-1}\omega_{p}^{\delta(j_{\pi(\phi)}-i_{\pi(\phi)})}\\\
&\implies\sum_{\delta=1}^{p-1}\omega_{\lambda}^{\left((a^{\gamma})_{i^{\delta}}-(a^{\gamma})_{j^{\delta}}\right)-\left((a^{\gamma})_{i}-(a^{\gamma})_{j}\right)}=-1\\\
&\implies\sum_{\delta=1}^{p-1}\omega_{\lambda}^{\left((a^{\gamma})_{i^{\delta}}-(a^{\gamma})_{j^{\delta}}\right)}+\omega_{\lambda}^{\left((a^{\gamma})_{i}-(a^{\gamma})_{j}\right)}=0.\end{split}$
(14)

Hence, we have the result.

∎

In the following example, we illustrate the Theorem 1.

###### Example 1

Let, $p=3$, $m=3$, $s=2$, $q=6$ and $\pi$ be a permutation on $\\{2,3\\}$ such
that $\pi(2)=2$, $\pi(3)=3$. From Theorem 1, we can construct the function
$f:\mathbb{Z}_{3}^{3}\rightarrow\mathbb{Z}_{6}$, where
$f(v_{1},v_{2},v_{3})=2v_{2}v_{3}+5$. Then
$\\{\omega_{3}^{a^{\gamma}}:\gamma\in\mathbb{Z}_{3}\\}$ is a $(3,27,3)$-MSCS,
where $a^{\gamma}=f+2v_{2}\gamma$.

We can show that the set
$\\{\omega_{\lambda}^{a^{\gamma}}:\gamma=0,1,\dots,p-1\\}$ constructed in
Theorem 1 has some interesting property which is pretty straightforward and we
show it in the next corollary.

###### Corollary 1

Let $a^{\gamma}$ be the function defined in Theorem 1. Then we have

$\sum_{\gamma=0}^{p-1}\rho(\omega_{\lambda}^{a^{\gamma}})(\tau)=0,$ (15)

when $\left|\tau\right|>p^{s-1}$.

###### Proof:

We shall prove for $\tau>0$ as
$\rho(\mathbf{a})(-\tau)=\rho^{*}(\mathbf{a})(\tau)$. For $\tau>p^{s-1}$, let
$j=i+\tau$. As $\tau>p^{s-1}$, it is not possible that $i_{k}=j_{k},\forall
k\in\\{s,s+1,\dots,m\\}$. Now,
$(a^{\gamma})_{i}-(a^{\gamma})_{j}=(f_{i}-f_{j})+\frac{\lambda}{p}\left(i_{\pi(s)}-j_{\pi(s)}\right)\gamma$.
So, if $i_{\pi(s)}\neq j_{\pi(s)}$, then we have

$\sum_{\gamma=0}^{p-1}\omega_{\lambda}^{(a^{\gamma})_{i}-(a^{\gamma})_{j}}=\omega_{\lambda}^{f_{i}-f_{j}}\sum_{\gamma=0}^{p-1}\omega_{\lambda}^{\frac{\lambda}{p}\left(i_{\pi(s)}-j_{\pi(s)}\right)\gamma}=0.$
(16)

If $i_{\pi(s)}=j_{\pi(s)}$, then in a similar manner described in Theorem 1,
we find $i^{\delta}$ and $j^{\delta}$ for $\delta=1,2,\dots,p-1$. Now, arguing
similar to Theorem 1, we have

$\sum_{\gamma=0}^{p-1}\left(\sum_{\delta=1}^{p-1}\omega_{\lambda}^{\left((a^{\gamma})_{i^{\delta}}-(a^{\gamma})_{j^{\delta}}\right)}+\omega_{\lambda}^{\left((a^{\gamma})_{i}-(a^{\gamma})_{j}\right)}\right)=0.$
(17)

Hence, the result. ∎

###### Remark 1

Corollary 1 shows that the set $\\{a^{\gamma}:\gamma=0,1,\dots,p-1\\}$ is in
fact a type-II $(p,p^{m},p^{m}-p^{s-1})$-ZCS. Type-II ZCS has application in
wireless communication. For example, it can be applied in wideband wireless
communication system having large minimum interfering signal delay (ISD) for
removing asynchronous interference [19, 20]. Type-II ZCSs are preferable than
type-I ZCSs in these scenarios to mitigate inter-symbol interference [21].
Although, there are some indirect [22] and direct [20] constructions of type-
II ZCS in the literature, the proposed construction is direct, as well as
provides flexible ZCZ width when the length is power-of-prime.

Next, we generalize the Theorem 1 for MSCS.

###### Theorem 2

Let $p_{1},p_{2},\dots,p_{k}$ be $k$ primes and $\lambda$ be a positive
integer such that $p_{\alpha}\mid\lambda,\forall\alpha=1,2,\dots,k;$ and
$\pi_{\alpha}$ be permutations on the sets
$I_{\alpha}=\\{s_{\alpha},s_{\alpha}+1,\dots,m_{\alpha}\\}$ for
$s_{\alpha}\geq 1$ and $m_{\alpha}\geq 1$, where $\alpha=1,2,\dots,k$. Let
$f_{\alpha}:\mathbb{Z}_{p_{\alpha}}^{m_{\alpha}}\rightarrow\mathbb{Z}_{\lambda}$
be defined by

$\begin{split}&f_{\alpha}(v_{p_{\alpha},1},v_{p_{\alpha},2},\dots,v_{p_{\alpha},m_{\alpha}})\\\
&=\frac{\lambda}{p_{\alpha}}\sum_{i=s_{\alpha}}^{m_{\alpha}-1}v_{p_{\alpha},\pi_{\alpha}(i)}v_{p_{\alpha},\pi_{\alpha}(i+1)}+\sum_{i=s_{\alpha}}^{m_{\alpha}}g_{p_{\alpha},i}v_{p_{\alpha},i}+g_{p_{\alpha}}\\\
&\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
+h_{\alpha}(v_{p_{\alpha},1}v_{p_{\alpha},2},\dots,v_{p_{\alpha},s_{\alpha}-1}),\end{split}$
(18)

where $g_{p_{\alpha},i},g_{p_{\alpha}}\in\mathbb{Z}_{\lambda}$,
$h_{\alpha}(v_{p_{\alpha},1}v_{p_{\alpha},2},\dots,v_{p_{\alpha},s_{\alpha}-1})$
is any function
$h_{\alpha}:\mathbb{Z}_{p_{\alpha}}^{s_{\alpha}-1}\rightarrow\mathbb{Z}_{\lambda}$
and $h_{\alpha}=0$ when $s_{\alpha}=1$. We define

$a^{\bm{\gamma}}:\mathbb{Z}_{p_{1}}^{m_{1}}\times\mathbb{Z}_{p_{2}}^{m_{2}}\times\dots\times\mathbb{Z}_{p_{k}}^{m_{k}}\rightarrow\mathbb{Z}_{\lambda},$

for
$\bm{\gamma}=(\gamma_{1},\gamma_{2},\dots,\gamma_{k})\in\mathbb{Z}_{p_{1}}\times\mathbb{Z}_{p_{2}}\times\dots\times\mathbb{Z}_{p_{k}}$
such that

$\begin{split}a^{\bm{\gamma}}=\sum_{\alpha=1}^{k}f_{\alpha}+\sum_{\alpha=1}^{k}\frac{\lambda}{p_{\alpha}}v_{p_{\alpha},\pi_{\alpha}(s_{\alpha})}\gamma_{\alpha}.\end{split}$
(19)

Then the set
$\\{\omega_{\lambda}^{a^{\bm{\gamma}}}:\bm{\gamma}\in\mathbb{Z}_{p_{1}}\times\mathbb{Z}_{p_{2}}\times\dots\times\mathbb{Z}_{p_{k}}\\}$
is a
$\left(\prod_{\alpha=1}^{k}p_{\alpha},\prod_{\alpha=1}^{k}p_{\alpha}^{m_{\alpha}},\prod_{\alpha=1}^{k}p_{\alpha}^{s_{\alpha}-1}\right)$-MSCS.

Table I: Comparison Table for Parameters Construction | Method | Length |  Set Size |  Value of $S$ | Constraints
---|---|---|---|---|---
[16] | GBF | $2^{m}$ | $2$ | $2^{d}$ |  $m\geq 1$, $0\leq d<m$
[17] | Matrix operation |  $2N+1$, $2N+2$, $2N_{1}$, $N_{1}N$, $N_{1}+N_{2}$ | $2$ | $2$ |  $N$ is length of existing GCP, $N_{1}$ and $N_{2}$ are lengths of existing MSCS
[18] | GBF | $2^{m-1}+2^{t}$ | $2^{k-s+1}$ | $2^{s}$ |  $m\geq 2$, $1\leq t\leq m-1$, $k\leq m-1$, $1\leq s\leq k^{\prime}\leq k$
Proposed, Theorem 2 |  Multivariable function | $\prod_{\alpha=1}^{k}p_{\alpha}^{m_{\alpha}}$ | $\prod_{\alpha=1}^{k}p_{\alpha}$ | $\prod_{\alpha=1}^{k}p_{\alpha}^{s_{\alpha}-1}$ |  $p_{\alpha}$’s are prime $\forall\alpha$, $m_{\alpha}\geq 1,\forall\alpha$, $1\leq s_{\alpha}\leq m_{\alpha}$, $\forall\alpha$.
Proposed, Theorem 3 |  Multivariable function | $p_{k+1}\prod_{\alpha=1}^{k}p_{\alpha}^{m_{\alpha}}$ | $\prod_{\alpha=1}^{k}p_{\alpha}$ | $p_{k+1}$ |  $p_{\alpha}$’s are prime $\forall\alpha$, $m_{\alpha}\geq 1,\forall\alpha=1,2,\dots,k$, $s_{\alpha}=1$, $\forall\alpha=1,2,\dots,k$.

###### Proof:

We prove this by induction. From Theorem 1, it is evident that the statement
is true for $k=1$. We assume that the statement is true for $k=n$. Then we
have to prove that it is true for $k=n+1$. For $k=n+1$, we have

$\begin{split}a^{\bm{\gamma}}=&\sum_{\alpha=1}^{n+1}f_{\alpha}+\sum_{\alpha=1}^{n+1}\frac{\lambda}{p_{\alpha}}v_{p_{\alpha},\pi_{\alpha}(s_{\alpha})}\gamma_{\alpha}\\\
=&\left(\sum_{\alpha=1}^{n}f_{\alpha}+\sum_{\alpha=1}^{n}\frac{\lambda}{p_{\alpha}}v_{p_{\alpha},\pi_{\alpha}(s_{\alpha})}\gamma_{\alpha}\right)\\\
&+\left(f_{n+1}+\frac{\lambda}{p_{n+1}}v_{p_{n+1},\pi_{n+1}(s_{n+1})\gamma_{n+1}}\right)\\\
=&R+S\end{split}$ (20)

where

$R=\sum_{\alpha=1}^{n}f_{\alpha}+\sum_{\alpha=1}^{n}\frac{\lambda}{p_{\alpha}}v_{p_{\alpha},\pi_{\alpha}(s_{\alpha})}\gamma_{\alpha},$
(21)

and

$S=f_{n+1}+\frac{\lambda}{p_{n+1}}v_{p_{n+1},\pi_{n+1}(s_{n+1})\gamma_{n+1}}.$
(22)

Now,
$\omega_{\lambda}^{a^{\bm{\gamma}}}=\omega_{\lambda}^{S}\otimes\omega_{\lambda}^{R}$,
where $\otimes$ denotes the Kronecker product. It is easy to observe that the
length of each sequence in
$\\{\omega_{\lambda}^{a^{\bm{\gamma}}}:\bm{\gamma}\in\mathbb{Z}_{p_{1}}\times\mathbb{Z}_{p_{2}}\times\dots\times\mathbb{Z}_{p_{k}}\\}$
is $L=\prod_{\alpha=1}^{n+1}p_{\alpha}^{m_{\alpha}}$ and number of sequence is
$M=\prod_{\alpha=1}^{n+1}p_{\alpha}$. We need to show that

$\sum_{\bm{\gamma}\in\mathbb{Z}_{p_{1}}\times\mathbb{Z}_{p_{2}}\times\dots\times\mathbb{Z}_{p_{n+1}}}\rho(\omega_{\lambda}^{a^{\bm{\gamma}}})(\tau)=0$
(23)

when $\tau\mod\prod_{\alpha=1}^{n+1}p^{s_{\alpha}-1}=0$. We note that for two
sequences $\mathbf{a}$ and $\mathbf{b}$ having length $L_{1}$ and $L_{2}$,
respectively, we have

$\begin{split}\rho(\mathbf{a}\otimes\mathbf{b})(\tau)=&\rho(\mathbf{a})\left(\left\lfloor\frac{\tau}{L_{2}}\right\rfloor\right)\rho(\mathbf{b})\left(\tau\mod
L_{2}\right)\\\
+&\Delta_{L_{2}}\rho(\mathbf{a})\left(\left\lfloor\frac{\tau}{L_{2}}\right\rfloor+1\right)\rho(\mathbf{b})\left(\tau\\!\\!\\!\\!\\!\mod
L_{2}-L_{2}\right),\\\ \end{split}$ (24)

where

$\Delta_{L_{2}}=\begin{cases}0,&\tau\mod L_{2}=0;\\\
1,&\text{otherwise},\end{cases}$ (25)

and $\left\lfloor\cdot\right\rfloor$ is the floor function. So, we can write

$\begin{split}&\sum_{\bm{\gamma}\in\mathbb{Z}_{p_{1}}\times\mathbb{Z}_{p_{2}}\times\dots\times\mathbb{Z}_{p_{n+1}}}\rho(\omega_{\lambda}^{a^{\bm{\gamma}}})(\tau)\\\
&=\bigg{[}\sum_{\gamma_{n+1}\in\mathbb{Z}_{p_{n+1}}}\rho(\omega_{\lambda}^{S})\left(\left\lfloor\frac{\tau}{L_{2}}\right\rfloor\right)\\\
&\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
\times\sum_{\bm{\gamma}^{\prime}\in\mathbb{Z}_{p_{1}}\times\dots\times\mathbb{Z}_{p_{n}}}\rho(\omega_{\lambda}^{R})\left(\tau\mod
L_{2}\right)\bigg{]}\\\ &\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
+\bigg{[}\Delta_{L_{2}}\times\sum_{\gamma_{n+1}\in\mathbb{Z}_{p_{n+1}}}\rho(\omega_{\lambda}^{S})\left(\left\lfloor\frac{\tau}{L_{2}}\right\rfloor+1\right)\\\
&\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
\times\sum_{\bm{\gamma}^{\prime}\in\mathbb{Z}_{p_{1}}\times\dots\times\mathbb{Z}_{p_{n}}}\rho(\omega_{\lambda}^{R})\left(\tau\\!\\!\\!\\!\\!\mod
L_{2}-L_{2}\right)\bigg{]}\\\ \end{split}$ (26)

where
$\bm{\gamma}=(\bm{\gamma}^{\prime},\gamma_{n+1})$,$\bm{\gamma}^{\prime}=(\gamma_{1},\gamma_{2},\dots,\gamma_{n})$
and $L_{2}=\prod_{\alpha=1}^{n}p_{\alpha}^{m_{\alpha}}$. Now, let
$\tau\mod\prod_{\alpha=1}^{n+1}p_{\alpha}^{s_{\alpha}-1}=0$, i.e.,
$\tau=\beta\prod_{\alpha=1}^{n+1}p_{\alpha}^{s_{\alpha}-1}$ for some
$\beta\in\\{1,2,\dots,\prod_{\alpha=1}^{n+1}p_{\alpha}^{m_{\alpha}-s_{\alpha}+1}-1\\}$.
But either
$\beta=\beta_{0}\prod_{\alpha=1}^{n}p_{\alpha}^{m_{\alpha}-s_{\alpha}+1}$ for
some integer $\beta_{0}$, or
$\prod_{\alpha=1}^{n}p_{\alpha}^{m_{\alpha}-s_{\alpha}+1}\nmid\beta$. Now, we
have two cases.

  1. Case I:

Let, $\beta=\beta_{0}\prod_{\alpha=1}^{n}p_{\alpha}^{m_{\alpha}-s_{\alpha}+1}$
for some integer $\beta_{0}$. Then $\tau\mod L_{2}=0$ and $\Delta_{L_{2}}=0$.
Also, in this case,
$\left\lfloor\frac{\tau}{L_{2}}\right\rfloor=\beta_{0}p_{n+1}^{s_{n+1}-1}$,
i.e., $\left\lfloor\frac{\tau}{L_{2}}\right\rfloor\mod p_{n+1}^{s_{n+1}-1}=0$.
So, from our assumption, we have

$\sum_{\gamma_{n+1}\in\mathbb{Z}_{p_{n+1}}}\rho(\omega_{\lambda}^{S})\left(\left\lfloor\frac{\tau}{L_{2}}\right\rfloor\right)=0.$
(27)

  2. Case II:

Let, $\prod_{\alpha=1}^{n}p_{\alpha}^{m_{\alpha}-s_{\alpha}+1}$ does not
divide $\beta$. Then $\tau\mod L_{2}\neq 0$ and $\Delta_{L_{2}}=1$. Let
$\tau\mod L_{2}=\tau_{1}$, where $\tau_{1}=\tau-\beta_{1}L_{2}$ for some
integer $\beta_{1}$ and $0<\tau_{1}<L_{2}$. But
$\prod_{\alpha=1}^{n}p_{\alpha}^{s_{\alpha}-1}$ divides both $\tau$ and
$L_{2}$. Hence $\tau_{1}\mod\prod_{\alpha=1}^{n}p_{\alpha}^{s_{\alpha}-1}=0$,
i.e., $\left(\tau\mod
L_{2}\right)\mod\prod_{\alpha=1}^{n}p_{\alpha}^{s_{\alpha}-1}=0$. But from our
assumption, we have

$\sum_{\bm{\gamma}^{\prime}\in\mathbb{Z}_{p_{1}}\times\dots\times\mathbb{Z}_{p_{n}}}\rho(\omega_{\lambda}^{R})\left(\tau\mod
L_{2}\right)=0.$ (28)

Arguing in a similar manner, we can say

$\sum_{\bm{\gamma}^{\prime}\in\mathbb{Z}_{p_{1}}\times\dots\times\mathbb{Z}_{p_{n}}}\rho(\omega_{\lambda}^{R})\left(\tau\\!\\!\\!\\!\\!\mod
L_{2}-L_{2}\right)=0.$ (29)

Hence, the case for $k=n+1$ is proved and we have the result. ∎

###### Remark 2

If $s_{\alpha}=1,\forall\alpha=1,2,\dots,k$, then the set
$\\{\omega_{\lambda}^{a^{\bm{\gamma}}}:\bm{\gamma}\in\mathbb{Z}_{p_{1}}\times\mathbb{Z}_{p_{2}}\times\dots\times\mathbb{Z}_{p_{k}}\\}$
becomes a GCS of length $L=p_{1}^{m_{1}}p_{2}^{m_{2}}\dots p_{k}^{m_{k}}$ and
set size $M=p_{1}p_{2}\dots p_{k}$, which is given in [13]. So, the
construction of GCS in [13] is a special case of the proposed construction.

Now, we propose another construction of MSCS.

###### Theorem 3

Let $p_{1},p_{2},\dots,p_{k},p_{k+1}$ be $(k+1)$ distinct primes and $\lambda$
be a positive integer such that $p_{\alpha}|\lambda,\forall\alpha$. Let
$a^{\bm{\gamma}}:\mathbb{Z}_{p_{1}}^{m_{1}}\times\mathbb{Z}_{p_{2}}^{m_{2}}\times\dots\times\mathbb{Z}_{p_{k}}^{m_{k}}\rightarrow\mathbb{Z}_{\lambda}$
be the same function given in Theorem 2 for
$s_{\alpha}=1,\forall\alpha=\\{1,2,\dots,k\\}$. Then the set
$\\{b^{\bm{\gamma}}:\bm{\gamma}\in\mathbb{Z}_{p_{1}}\times\mathbb{Z}_{p_{2}}\times\dots\times\mathbb{Z}_{p_{k}}\\}$
is an
$(\prod_{\alpha=1}^{k}p_{\alpha},p_{k+1}\prod_{\alpha=1}^{k}p_{\alpha}^{m_{\alpha}},p_{k+1})$-MSCS,
where
$b^{\bm{\gamma}}:\mathbb{Z}_{p_{1}}^{m_{1}}\times\mathbb{Z}_{p_{2}}^{m_{2}}\times\dots\times\mathbb{Z}_{p_{k+1}}^{m_{k+1}}\rightarrow\mathbb{Z}_{\lambda}$
is defined by
$b^{\bm{\gamma}}=a^{\bm{\gamma}}+g_{p_{k+1},1}v_{p_{k+1},1}+g_{p_{k+1}}$,
where $g_{p_{k+1},1},g_{p_{k+1}}\in\mathbb{Z}_{\lambda}$.

###### Proof:

The proof follows form the fact that $b^{\bm{\gamma}}=f_{p_{k+1}}\otimes
a^{\bm{\gamma}}$, where $f_{p_{k+1}}=g_{p_{k+1},1}v_{p_{k+1},1}+g_{p_{k+1}}$,
and using (24) and Remark 2, we can get

$\begin{split}\sum_{\bm{\gamma}\in\mathbb{Z}_{p_{1}}\times\mathbb{Z}_{p_{2}}\times\dots\times\mathbb{Z}_{p_{k}}}\rho(\omega_{\lambda}^{a^{\bm{\gamma}}})(\tau\mod
L_{2})=&0,\\\
\sum_{\bm{\gamma}\in\mathbb{Z}_{p_{1}}\times\mathbb{Z}_{p_{2}}\times\dots\times\mathbb{Z}_{p_{k}}}\rho(\omega_{\lambda}^{a^{\bm{\gamma}}})(\tau\mod
L_{2}-L_{2})=&0,\\\ \end{split}$ (30)

when $\tau\mod p_{k+1}=0$ and
$L_{2}=\prod_{\alpha=1}^{k}p_{\alpha}^{m_{\alpha}}$. ∎

###  III-A Peak-to-Mean Envelope Power Ratio (PMEPR)

Let $\mathbf{x}=(x_{0},x_{1},\dots,x_{L-1})$ be a
$\mathbb{Z}_{\lambda}$-valued sequence. Then the OFDM signal is the real part
of the complex envelope

$P_{\mathbf{x}}(t)=\sum_{i=0}^{L-1}\omega_{\lambda}^{x_{i}+\lambda f_{i}t}$
(31)

where $f_{i}=f+i\Delta f$, $f$ is a constant frequency, $\Delta f$ is a
integer multiple of OFDM symbol rate. The term
$\frac{\left|P_{\mathbf{x}}(t)\right|^{2}}{L}$ is called instantaneous-to-
average power ratio (IAPR). PMEPR of the sequence $\mathbf{x}$ is defined as

$\text{PMEPR}(\mathbf{x})=\sup_{0\leq\Delta ft\leq
1}\frac{\left|P_{\mathbf{x}}(t)\right|^{2}}{L}.$ (32)

We shall use $\text{PMEPR}(\mathbf{x})$ and $\text{PMEPR}(\psi(\mathbf{x}))$
invariably whenever the context is clear. Similarly, for a set
$A=\\{\mathbf{x}_{0},\mathbf{x}_{1},\dots,\mathbf{x}_{M-1}\\}$ of sequences
PMEPR can be defined as

$\text{PMEPR}(A)=\max\\{\text{PMEPR}(\mathbf{x}_{i}):i=0,1,\dots,M-1\\}.$ (33)

###  III-B Bound for PMEPR

In this subsection, we calculate the PMEPR of the constructed MSCS. First we
state a result regarding the PMEPR which is available in the literature.

######  Lemma 1 ([15])

If $\mathbf{x}$ is a sequence of a $(2,L,S)$-MSCS, then
$\text{PMEPR}(\mathbf{x})$ is upper bounded by $2S$.

A similar statement can be made for a $(M,L,S)$-MSCS using methods similar to
[16], which we provide in the following lemma.

###### Lemma 2

If $\mathbf{x}$ is a sequence of a $(M,L,S)$-MSCS, then
$\text{PMEPR}(\mathbf{x})$ is upper bounded by $MS$.

###### Proof:

We briefly sketch the proof here. Let
$\\{\mathbf{a}_{0},\mathbf{a}_{1},\dots,\mathbf{a}_{M-1}\\}$ be an
$(M,L,S)$-MSCS, where
$\mathbf{a}_{i}=\left(a_{i,0}a_{i,1},\dots,a_{i,L-1}\right)$ for all $i$. We
let $\zeta=\exp\left(\frac{2\pi\sqrt{-1}}{S}\right)$ be the $S$-th primitive
root of $1$. We define the set of sequences
$\mathbf{a}_{i}^{u}=\left(a_{i,0}\zeta^{0u},a_{i,1}\zeta^{1u},\dots,a_{i,(L-1)}\zeta^{(L-1)u}\right)$
for $u=0,1,\dots,S-1$ and $i=0,1,\dots,M-1$. Then we have

$\begin{split}\sum_{u=0}^{S-1}\left|P_{\mathbf{a}_{i}^{u}}(t)\right|^{2}=&\sum_{u=0}^{S-1}\left|\sum_{k=0}^{L-1}a_{i,k}\zeta^{ku}z^{f_{k}}\right|^{2}\\\
=&LS+2\mathcal{R}\left(\sum_{k=0}^{L-1}\rho(\mathbf{a}_{i})(k)z^{f_{k}}\sum_{u=0}^{S-1}\zeta^{ku}\right),\end{split}$
(34)

where $z=\exp\left(2\pi t\sqrt{-1}\right)$ and $\mathcal{R}(\cdot)$ is the
real part of a complex number. As $\sum_{u=0}^{S-1}\zeta^{ku}$ equals $0$,
when $k\\!\\!\\!\mod S\neq 0$ and it equals $S$, otherwise, we have

$\begin{split}\sum_{i=0}^{M-1}\sum_{u=0}^{S-1}\left|P_{\mathbf{a}_{i}^{u}}(t)\right|^{2}=&MLS.\end{split}$
(35)

Hence, $\text{PMEPR}(\mathbf{a}_{i}^{u})\leq MS$. For $u=0$, we have
$\mathbf{a}_{i}^{u}=\mathbf{a}_{i}$ and it implies
$\text{PMEPR}(\mathbf{a}_{i})$ is upper bounded by $MS$. ∎

###### Example 2

Let $k=2$, $p_{1}=3$, $p_{2}=2$, $m_{1}=3$, $m_{2}=1$, $s_{1}=1$. Then, using
Theorem 3, we can directly construct a $(3,54,2)$-MSCS which has not been
reported before. We take
$f:\mathbb{Z}_{3}^{3}\times\mathbb{Z}_{2}\rightarrow\mathbb{Z}_{6}$, where
$f_{1}=2(v_{3,2}v_{3,3}+v_{3,3}v_{3,1})+2v_{3,1}+5v_{3,2}+v_{3,3}$ and
$f_{2}=3v_{2,1}$. In Fig. 1, we have shown plot of IAPR with respect to
$\Delta ft$ for all sequences of the MSCS. The PMEPR value is 5.9465 which
agree to the theoretical upper bound $6$ from Lemma 2.

Figure 1: IAPRs of constituent sequences of MSCS having length $54$

###  III-C Comparison with existing works

In the TABLE I, we have compared the proposed construction with the existing
ones [16, 17, 18] with respect to their parameters. It can be seen from TABLE
I that the proposed constructions provide flexible lengths and set sizes
compared to the existing constructions.

##  IV Conclusion

In this paper, we have proposed new constructions of MSCSs, which can be used
as alternatives to the conventional GCSs in OFDM. Although there are several
direct constructions of MSCSs in the literature, the sequence lengths and set
sizes are limited in those. The proposed constructions provide flexible as
well as arbitrary sequence lengths and flexible set sizes. For $k=1$, one of
the constructed MSCS reduces to type-II ZCS, which has applications in
wideband wireless communication systems.

## References

  * [1] M. J. E. Golay, “Static multislit spectrometry and its application to the panoramic display of infrared spectra.” _J. Opt. Soc. Am._ , vol. 41, no. 7, pp. 468–472, July 1951.
  * [2] B. Popovic, “Synthesis of power efficient multitone signals with flat amplitude spectrum,” _IEEE Trans. Commun._ , vol. 39, no. 7, pp. 1031–1033, 1991.
  * [3] R. Van Nee, “OFDM codes for peak-to-average power reduction and error correction,” in _1996 IEEE Global Telecommunications Conference (GLOBECOM)_ , vol. 1, 1996, pp. 740–744 vol.1.
  * [4] G. Welti, “Quaternary codes for pulsed radar,” _IRE Trans. Inf. Theory_ , vol. 6, no. 3, pp. 400–408, 1960.
  * [5] P. Spasojevic and C. Georghiades, “Complementary sequences for ISI channel estimation,” _IEEE Trans. Inf. Theory_ , vol. 47, no. 3, pp. 1145–1152, 2001\.
  * [6] J. Davis and J. Jedwab, “Peak-to-mean power control in OFDM, Golay complementary sequences, and Reed-Muller codes,” _IEEE Trans. Inf. Theory_ , vol. 45, no. 7, pp. 2397–2417, 1999.
  * [7] C.-C. Tseng and C. Liu, “Complementary sets of sequences,” _IEEE Trans. Inf. Theory_ , vol. 18, no. 5, pp. 644–652, 1972.
  * [8] K.-U. Schmidt, “Complementary sets, generalized Reed–Muller codes, and power control for OFDM,” _IEEE Trans. Inf. Theory_ , vol. 53, no. 2, pp. 808–814, 2007.
  * [9] C.-Y. Chen, “Complementary sets of non-power-of-two length for peak-to-average power ratio reduction in OFDM,” _IEEE Trans. Inf. Theory_ , vol. 62, no. 12, pp. 7538–7545, 2016.
  * [10] ——, “A novel construction of complementary sets with flexible lengths based on Boolean functions,” _IEEE Commun. Lett._ , vol. 22, no. 2, pp. 260–263, 2018.
  * [11] S. Das, S. Majhi, S. Budišin, and Z. Liu, “A new construction framework for polyphase complete complementary codes with various lengths,” _IEEE Trans. Signal Process._ , vol. 67, no. 10, pp. 2639–2648, 2019.
  * [12] S. Das, S. Majhi, and Z. Liu, “A novel class of complete complementary codes and their applications for APU matrices,” _IEEE Signal Process. Lett._ , vol. 25, no. 9, pp. 1300–1304, 2018.
  * [13] P. Sarkar, Z. Liu, and S. Majhi, “Multivariable function for new complete complementary codes with arbitrary lengths,” 2021. [Online]. Available: https://arxiv.org/abs/2102.10517
  * [14] K. Paterson, “Generalized Reed-Muller codes and power control in OFDM modulation,” _IEEE Trans. Inf. Theory_ , vol. 46, no. 1, pp. 104–120, 2000\.
  * [15] Y. Xin and I. Fair, “Multiple-shift complementary sequences and their peak-to-average power ratio values,” in _International Symposium on Information Theory, 2004. ISIT 2004. Proceedings._ , 2004, pp. 121–.
  * [16] W. Chen and C. Tellambura, “Identifying a class of multiple shift complementary sequences in the second order cosets of the first order Reed-Muller codes,” in _IEEE International Conference on Communications, 2005. ICC 2005. 2005_ , vol. 1, 2005, pp. 618–621 Vol. 1.
  * [17] B. Shen, Y. Yang, and Z. Zhou, “A construction of binary Golay complementary sets based on even-shift complementary pairs,” _IEEE Access_ , vol. 8, pp. 29 882–29 890, 2020.
  * [18] Y.-J. Lin, Z.-M. Huang, and C.-Y. Chen, “Golay complementary sets and multiple-shift complementary sets with non-power-of-two length and bounded PAPRs,” _IEEE Commun. Lett._ , vol. 25, no. 9, pp. 2805–2809, 2021.
  * [19] Z. Liu, U. Parampalli, and Y. L. Guan, “Optimal odd-length binary Z-complementary pairs,” _IEEE Trans. Inf. Theory_ , vol. 60, no. 9, pp. 5768–5781, 2014.
  * [20] R. Kumar, P. Sarkar, P. K. Srivastava, and S. Majhi, “A direct construction of asymptotically optimal type-II ZCP for every possible even length,” _IEEE Signal Process. Lett._ , vol. 28, pp. 1799–1802, 2021.
  * [21] W. C. Y. Lee, _Mobile Communications Design Fundamentals_ , 2nd ed. USA: John Wiley & Sons, Inc., 1992.
  * [22] R. Kumar, P. K. Srivastava, and S. Majhi, “A construction of type-II ZCCS of arbitrary sequence length with low PMEPR,” 2022. [Online]. Available: https://arxiv.org/abs/2207.02552


# Posterior Sampling of the Initial Conditions of the Universe from Non-linear
Large Scale Structures using Score-Based Generative Models

Ronan Legin,1,2,3 Matthew Ho,4 Pablo Lemos,1,2,3,5 Laurence Perreault-
Levasseur,1,2,3,5,6 Shirley Ho,5Yashar Hezaveh,1,2,3,5,6 and Benjamin
Wandelt4,5,7
1Department of Physics, Université de Montréal, Montréal, Canada
2Mila - Quebec Artificial Intelligence Institute, Montréal, Canada
3Ciela - Montreal Institute for Astrophysical Data Analysis and Machine
Learning, Montréal, Canada
4Sorbonne Université, CNRS, UMR 7095, Institut d’Astrophysique de Paris, 98
bis bd Arago, 75014 Paris, France
5Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue,
10010, New York, NY, USA
6Perimeter Institute for Theoretical Physics, Waterloo, Ontario, Canada, N2L
2Y5
7Sorbonne Université, Institut Lagrange de Paris, 98 bis boulevard Arago,
75014 Paris, France E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Reconstructing the initial conditions of the universe is a key problem in
cosmology. Methods based on simulating the forward evolution of the universe
have provided a way to infer initial conditions consistent with present-day
observations. However, due to the high complexity of the inference problem,
these methods either fail to sample a distribution of possible initial density
fields or require significant approximations in the simulation model to be
tractable, potentially leading to biased results. In this work, we propose the
use of score-based generative models to sample realizations of the early
universe given present-day observations. We infer the initial density field of
full high-resolution dark matter N-body simulations from the present-day
density field and verify the quality of produced samples compared to the
ground truth based on summary statistics. The proposed method is capable of
providing plausible realizations of the early universe density field from the
initial conditions posterior distribution marginalized over cosmological
parameters and can sample orders of magnitude faster than current state-of-
the-art methods.

###### keywords:

methods: statistical – early Universe – large-scale structure of Universe

††pubyear: 2023††pagerange: Posterior Sampling of the Initial Conditions of
the Universe from Non-linear Large Scale Structures using Score-Based
Generative Models–A

##  1 Introduction

In the standard model of cosmology, structure originates from quantum
fluctuations of a primordial density field, which are scaled to macroscopic
distances by a physical process called inflation (Guth, 1981; Albrecht &
Steinhardt, 1982; Linde, 1982, 1983). This initial density field represents
the seed to all structure seen in the Universe today. Furthermore, different
models of inflation predict various levels of non-Gaussianity in the
primordial density field (Acquaviva et al., 2002; Maldacena, 2003; Bartolo et
al., 2004). Therefore, accurate methods to reconstruct this primordial field
could shed light on the unknown mechanism behind inflation and guide our
search for new physics.

Beyond the study of early universe physics, knowledge of the initial
conditions of the Universe can be combined with a forward model to compute
predictions for any observable on our past lightcone. Such predictions can act
as discovery templates for new physical effects in cross-correlation with
external data, e.g., for the discovery of secondary anisotropies in the
microwave sky; be used to perform posterior predictive tests of the underlying
cosmological physics model; or give observational access to quantities that
were hitherto only accessible in simulations, such as the dynamical assembly
history of elements of the cosmic web, such as clusters, filaments, or voids
(Lavaux & Wandelt, 2010; Leclercq et al., 2015; Jasche & Lavaux, 2019;
Feldbrugge & van de Weygaert, 2022). For these reasons, substantial effort has
been put in the inference of initial conditions as one of the key problems in
cosmology. Currently, the primary constraints on the initial conditions come
from linear reconstruction (Komatsu et al., 2005; Yadav & Wandelt, 2005)
applied to observations of the Cosmic Microwave Background (CMB; Planck
Collaboration, 2020a, b).

Upcoming galaxy surveys will provide vast amounts of information on small
scales, begging the question: can we infer the small-scale initial conditions
from non-linearly evolved structure? This inverse problem compounds the
challenges of high dimensionality (representing the initial density field on a
3-D grid of $10^{7}$ voxels corresponds to a $10^{7}$ dimensional inference
problem) and the complexity of computing the non-linear forward mapping
between the primordial initial conditions and the present-day density field,
requiring, at the very least, high-resolution N-body simulations (Springel,
2005; Villaescusa-Navarro et al., 2020a) and further modeling of the
distribution of observable tracers of the underlying dark matter field
(Somerville & Davé, 2015).

Current state-of-the-art methods such as Bayesian Origin Reconstruction of
Galaxies (BORG; Jasche & Wandelt, 2013; Jasche & Lavaux, 2019) use Hamiltonian
Monte Carlo, a Markov Chain Monte Carlo (MCMC) algorithm for Bayesian
parameter inference of the initial conditions, in which one can exploit the
gradient of the likelihood to efficiently generate posterior samples. However,
due to the high computational cost of full N-body simulations, BORG relies on
approximate simulation methods such as second-order Lagrange Perturbation
Theory (2LPT) and Particle-Mesh (PM) simulations, which are inaccurate at
small scales. Moreover, they require fully-differentiable simulators, which
has so far precluded including non-differentiable operations—such as halo
finding—that are part of standard modeling pipelines that map dark-matter
density fields to galaxy surveys.

In the hopes of improving upon these limitations, machine learning has been
used to reconstruct the initial conditions from simulations (Modi et al.,
2021; Shallue & Eisenstein, 2023). Unfortunately, due to the high-dimensional
parameter space, these works have been limited to predicting a single point-
estimate, as modeling the full multi-million-dimensional posterior
distribution is intractable. Since these models do not produce samples of the
early universe, they do not provide any measures of uncertainty on the
reconstructions.

In this work, we propose the use of score-based generative diffusion models
(Song & Ermon, 2019; Ho et al., 2020; Song et al., 2020) to learn the
distribution of early universe density fields conditioned on the present-day
matter density field and to produce samples from it. We train a neural network
to predict the score of the posterior distribution using simulations from the
Quijote Latin Hypercube set (Villaescusa-Navarro et al., 2020b). We then use
the estimate of the score network to solve a reverse-diffusion stochastic
differential equation to sample the posterior distribution of the initial
conditions.

##  2 Methods

###  2.1 Problem Overview

Our goal is to infer the 3-D density field of the early universe $\bm{x}$
given observations $\bm{y}$ of the dark matter distribution at low redshift.
We can define this problem within a Bayesian framework, where we are
interested in sampling from the posterior distribution $p(\bm{x}|\bm{y})$.
Using Bayes’s theorem, the posterior distribution $p(\bm{x}|\bm{y})$ can be
written as,

$p(\bm{x}|\bm{y})=\frac{p(\bm{y}|\bm{x})p(\bm{x})}{p(\bm{y})},$ (1)

where $p(\bm{y}|\bm{x})$, $p(\bm{x})$, $p(\bm{y})$ represent the likelihood,
prior and evidence respectively. The prior distribution $p(\bm{x})$ reflects
our knowledge on the possible realizations of the early universe, while the
evidence $p(\bm{y})$ gives the probability for a realization of the data. The
likelihood distribution $p(\bm{y}|\bm{x})$ represents the distribution of
possible observations $\bm{y}$ given a fixed realization of the early universe
$\bm{x}$. It includes our cosmological forward simulator and additional
effects that resemble true observations (e.g., selection functions, galaxy
shot noise). For our problem, we let $\bm{y}$ be the simulated present-day
comoving dark matter density field with added noise sampled from a normal
distribution with standard deviation of 0.1. As detailed in section 2.5, the
early and late density fields, $\bm{x}$ and $\bm{y}$, are represented on 3-D
mesh grids with $128^{3}$ voxels. Therefore, the sampling space of $\bm{x}$ is
multi-million in dimension. Because of the high dimensionality of the
inference problem, modeling the posterior using density estimation techniques
such as normalizing flows is not feasible.

Instead of directly modelling the posterior $p(\bm{x}|\bm{y})$, we can model
the gradient of the log posterior distribution $\nabla_{\bm{x}}\log
p(\bm{x}|\bm{y})$. In comparison, $\nabla_{\bm{x}}\log p(\bm{x}|\bm{y})$ is
computationally tractable, as it does not depend on the normalization of the
posterior distribution. This means that it can be approximated by a simple
neural network that learns a function $s(\bm{x},\bm{y})$ mapping a set of
inputs ($\bm{x}$, $\bm{y}$) to an output prediction of the score
$\nabla_{\bm{x}}\log p(\bm{x}|\bm{y})$. In the following section, we describe
how we can train a conditional neural network to estimate the score
$\nabla_{\bm{x}}\log p(\bm{x}|\bm{y})$ and use it in the context of score-
based generative models (Song & Ermon, 2019; Song et al., 2020) to produce
samples from the posterior distribution $p(\bm{x}|\bm{y})$.

Figure 1: Top: Early universe dark matter density field at redshift $z=127$.
Bottom: Present-day dark matter density field at redshift $z=0$. The goal is
to infer the early universe (top cube) given the present-day dark matter
density field (bottom cube). Each face of the cubes shows the outermost pixels
for a $(1\ \text{Gpc}/h)^{3}$ simulation volume. The density cubes are
discretized on a $128^{3}$ grid.

###  2.2 Score-Based Generative Models

Score-based generative modeling is a framework designed to learn the
distribution of variables from a dataset, by approximating the score
$\nabla_{\bm{x}}\log p(\bm{x})$, which is typically modeled by a neural
network (Song et al., 2020). The procedure also provides the possibility to
generate new samples from the learned distributions. They have been used to
great success across a wide variety of domains (e.g., Song et al., 2021; Popov
et al., 2021; Gnaneshwar et al., 2022; Anand & Achim, 2022; Adam et al., 2022;
Mudur & Finkbeiner, 2022; Legin et al., 2023), and have surpassed previous
state-of-the-art methods such as GANs (e.g., Dhariwal & Nichol, 2021; Müller-
Franzes et al., 2022). These achievements have been made possible due to a
number of technical improvements in regard to the sampling strategy used to
generate samples. In part, one important breakthrough consisted of framing the
sampling method as a reverse-diffusion process, where the data $\bm{x}$ is
perturbed at various noise levels, and sampling new data points $\bm{x}$
consists of iteratively reversing this process starting from pure noise. In
Song et al. (2020), this process is defined as the stochastic differential
equation (SDE):

$d\bm{x}=f(\bm{x},t)dt+g(t)d\bm{w},$ (2)

where $f(\bm{x},t)$ is called the drift term, $d\bm{w}$ is a Wiener process
characterizing the random noise, and $g(t)$ is a scalar function that
determines the level of added noise. The key to generating samples $\bm{x}$
from $p(\bm{x})$ lies in reversing the diffusion process by solving the
reverse-SDE:

$d\bm{x}=\left(f(\bm{x},t)-g(t)^{2}\nabla_{\bm{x}}\log
p_{t}(\bm{x})\right)dt+g(t)d\bm{w},$ (3)

where $p_{t}(\bm{x})$ is the probability distribution of $\bm{x}$ at time $t$.
Alternatively, we can extend the previous equation to reverse a conditional
diffusion process by solving backward in time

$d\bm{x}=\left(f(\bm{x},t)-g(t)^{2}\nabla_{\bm{x}}\log
p_{t}(\bm{x}|\bm{y})\right)dt+g(t)d\bm{w},$ (4)

which requires the conditional score $\nabla_{\bm{x}}\log
p_{t}(\bm{x}|\bm{y})$ at time $t$.

In this work, we train a neural network conditioned on the observation
$\bm{y}$ and time $t$, denoted $s(\bm{x},\bm{y},t)$ to learn the score of the
posterior, $\nabla_{\bm{x}}\log p_{t}(\bm{x}|\bm{y})$, via score matching
(Hyvärinen, 2005; Vincent, 2011; Song et al., 2020). We then solve Equation 4
to sample from the posterior distribution of initial conditions
$p(\bm{x}|\bm{y})$ by replacing the score $\nabla_{\bm{x}}\log
p_{t}(\bm{x}|\bm{y})$ by its approximation from the score network,
$s(\bm{x},\bm{y},t)$.

###  2.3 Solving the Reverse-SDE

There exist many numerical methods to solve the reverse-SDE from Equation 4.
We choose the Euler-Maruyama method which discretizes the SDE as

$\bm{x}_{t+\Delta
t}=\bm{x}_{t}+\left(f(\bm{x}_{t},t)-g(t)^{2}\nabla_{\bm{x}_{t}}\log
p_{t}(\bm{x}_{t}|y)\right)\Delta t+g(t)\bm{z}_{t}\sqrt{-\Delta t},$ (5)

where $\Delta t=-1/N$ is the step size, $N$ is the number of steps and
$\bm{z}_{t}$ is sampled from a standard normal distribution with the same
dimensions as $\bm{x}_{t}$. We perturb the initial conditions $\bm{x}$ with
noise following the Variance Exploding Stochastic Differential Equation
(VESDE) proposed in Song et al. (2020). We then solve the discretized reverse-
SDE from Equation 5 to sample from the posterior distribution
$p(\bm{x}|\bm{y})$. Based on the geometric interpretation from Song & Ermon
(2020), we choose a maximum and minimum diffusion noise level of
$\sigma_{\text{max}}=100$ and $\sigma_{\text{min}}=0.01$ respectively, and
choose $N=1000$ steps.

###  2.4 Score Network Architecture and Training

The score network architecture is based on the Pytorch implementation by Song
et al. (2020) available on GitHub111https://github.com/yang-
song/score_sde_pytorch. The network follows the RefineNet architecture (Lin et
al., 2016) with five downsampling and upsampling levels each with 2 ResNet
blocks from the BigGAN model (Brock et al., 2018) and with outputs of 32
features maps for the first, fourth and fifth level and 64 feature maps for
the second and third level. Each ResNet block contains a Dropout layer
(Srivastava et al., 2014) with a dropout rate of $10\%$ during training. The
main difference with the implementation by Song et al. (2020) is that $\bm{y}$
is concatenated to $\bm{x}_{t}$ along the channel dimension and fed as input
to the network. A proof showing the validity of this input scheme to learn the
score of conditional probability distributions can be found in Batzolis et al.
(2021). We train the network for approximately 400 epochs with a batch size of
8 split across 4 NVIDIA H100 80GB Graphical Processing Units (GPU). We use the
Adam optimizer with a learning rate of $2\times 10^{-4}$ and clip the gradient
of the weights to a maximum gradient norm of 1. The duration of training under
these settings is approximately 24 hours.

###  2.5 Simulations

We use density fields from the Quijote Latin Hypercube set of N-body
simulations (Villaescusa-Navarro et al., 2020b) to train and test the score
network. Specifically, we use the set of 2,000 $(1\ \text{Gpc}/h)^{3}$ dark
matter N-body simulations without massive neutrinos. The simulations
initialized with different random seeds run with different values for the
cosmological parameters in the range $\Omega_{m}\in[0.1,0.5]$,
$\Omega_{b}\in[0.03,0.07]$, $h\in[0.5,0.9]$, $n_{s}\in[0.8,1.2]$ and
$\sigma_{8}\in[0.6,1.0]$. We use 1,900 N-body simulations for training and 100
for testing with the $\bm{x}$ and $\bm{y}$ density fields computed on a
$128^{3}$ grid at redshift $z=127$ and redshift $z=0$ respectively. An example
of initial conditions and present-day density field is shown in Figure 1. In
this work, the redshift $z=127$ density fields we train the score network with
are individually normalized by their own variance, resulting in generated
posterior samples of the early universe density field with variance of one. We
found that this improves training convergence and in practice, the generated
samples can be rescaled with the true variance predicted from the redshift
$z=0$ density field power spectrum from e.g. CAMB (Lewis et al., 2000;
Villaescusa-Navarro et al., 2020b). Note that the samples predicted at
redshift $z=127$ are within the linear regime on the scales we consider and
therefore can be related directly to an earlier density field (e.g., $z\sim
1000$) using linear theory.

##  3 Results

Figure 2: Left: The density field at redshift $z=0$ for the fiducial Planck
cosmology. Center: Initial conditions sampled from the posterior
$p(\bm{x}|\bm{y})$. Right: The true initial conditions. All three density
fields span a $1000\times 1000\times 125(h^{-1}\text{Mpc})^{3}$ region
averaged over the third axis. This example demonstrates the capability of
score-based generative models to sample highly detailed initial conditions
consistent with the ground truth. See Figure 3 for quantification of
uncertainty. Figure 3: Posterior variance of initial condition samples from
the fiducial Planck cosmology simulation per pixel averaged over the depth of
a volume of $1000\times 1000\times 125(h^{-1}\text{Mpc})^{3}$. Reconstruction
variance is high in patches that will collapse into large halos at the present
day (see Figure 2). As expected, posterior variance also increases near the
boundaries of the data volume.

We show results using our score network on density fields from the fiducial
Quijote set of simulations using the fiducial Planck cosmology with
$\Omega_{m}=0.3175$, $\Omega_{b}=0.049$, $h=0.6711$, $n_{s}=0.9624$ and
$\sigma_{8}=0.834$, and from six of the Quijote Latin Hypercube set of
simulations from our test set with different cosmological parameter values.
Note that our score network was not trained on density fields from fiducial
cosmology simulations. For each simulation, we solve Equation 4 to produce 100
samples from the posterior distribution $p(\bm{x}|\bm{y})$ using our trained
score network $s(\bm{x},\bm{y},t)$. In Figure 2 and in appendix A, we show
examples of sampled initial conditions for a simulation with the fiducial
cosmology. In Figure 3, we show a map of the posterior sample variance based
on the same present-day density field shown in Figure 2. We then compare the
power spectrum of the predicted samples and ground truth, compute the transfer
function and cross-correlation, and show the results in Figure 4 for the
fiducial cosmology case and in appendix A for the latin-hypercube simulations.

Moreover, we perform a posterior coverage test to verify if the distribution
of samples accurately cover the true posterior distribution. This consists of
computing the difference between the true initial conditions and the mean of
the initial conditions posterior samples, and then dividing by the standard
deviation of the posterior samples. If the distribution of the resulting
values is consistent with the statistics of a normal distribution (mean of
zero and standard deviation of one), then the generated posterior samples have
the correct first two marginal moments of the true posterior distribution
$p(\bm{x}|\bm{y})$. We show results of this coverage test in Figure 5 for the
fiducial simulation.

Figure 4: Test statistics to verify the accuracy of samples from the posterior
distribution of initial conditions against the ground truth from the fiducial
Planck cosmology simulation. Note that this cosmology was not used for
training the score network. Top: The power spectrum of the redshift $z=127$
density field vs the true density field. Middle: The cross-correlation $C(k)$
between every posterior sample and the true field. For comparison we show same
for the data and the true field, corresponding to the best achievable $C(k)$
for an optimal linear filter (green dashed). The posterior samples contain
significantly more information. Bottom: The transfer function between the
posterior samples and the true field with the Nyquist frequency shown as the
vertical black line. The shaded regions represent $1\sigma$ and $2\sigma$
errors. The results illustrate the high level of accuracy of the inferred
initial conditions. Figure 5: Posterior coverage test for the fiducial Planck
cosmology: we show the histogram of the difference between the true initial
condition $\bm{x}_{\text{true}}$ and the posterior mean $\mu_{\text{samples}}$
normalized by the posterior standard deviation $\sigma_{\text{samples}}$. If
the inferred posterior samples of initial conditions have the correct mean and
variance of the voxel-wise marginal posterior distributions, we expect this
histogram to have mean 0 and standard deviation 1.

##  4 Discussion

In this work, sampling possible initial conditions from the multi-million-
dimensional posterior distribution is done by solving a diffusion process
backwards in time. This requires the score $\nabla_{\bm{x}}\log
p_{t}(\bm{x}|\bm{y})$, which guides the diffusion towards the posterior
distribution while random noise diffuses the samples in order to explore the
sampling space. Compared to the posterior probability $p(\bm{x}|\bm{y})$,
modelling $\nabla_{\bm{x}}\log p_{t}(\bm{x}|\bm{y})$ is computationally
tractable as it is independent of the normalization of $p(\bm{x}|\bm{y})$,
which is intractable given the high dimensionality of the inference problem.
Therefore, methods that directly learn the posterior density using models such
as masked autoregressive flows (Papamakarios et al., 2017) are not suitable
for this task.

Efficient MCMC sampling methods in high dimensions such as Hamiltonian Monte
Carlo (HMC), as implemented in BORG (Jasche & Wandelt, 2013; Jasche & Lavaux,
2019), also do not need a normalized posterior distribution, but they do
require a differentiable forward model that must be run for each step in the
integration of the Hamiltonian trajectory. For our problem, each forward step
amounts to running an N-body simulation from initial conditions until the
present day. Each such run must be combined with an adjoint run to compute the
gradient. Therefore, the generation of each new sample from the posterior
requires tens of simulations. As in all MCMC techniques, accepted samples are
correlated leading to a burn-in stage that limits the parallelism of the
approach since each running chain must first converge to the posterior
distribution before generating useful samples. Moreover, the effective number
of samples in the chain is smaller than the number of accepted samples by a
factor proportional to the correlation length of the chain. As a result BORG
must run on the order of 100 N-body simulations for each approximately
independent sample drawn from the posterior distribution, although this could
be improved by fitting a variational approximation for the MCMC proposal
distribution (Modi et al., 2023). Because of this, BORG uses approximate
N-body simulators such as PM, which sacrifice accuracy for speed, as they are
typically orders of magnitude faster than full N-body simulations. Even so, it
can require upwards of 1000 CPU hours to generate a single independent
posterior sample of initial conditions. Given this, using a full N-body
simulator as the forward model for BORG is not computationally feasible.

In contrast, the method used in this work performs inference using full N-body
simulations and can generate independent samples of initial conditions from
the posterior distribution within minutes on a single GPU. The sampling
procedure can also be trivially parallelized by independently solving Equation
4 across multiple devices. Furthermore, our approach samples from the
posterior distribution of initial conditions marginalized over cosmological
parameters. There are straightforward generalizations to score-based sampling
that will allow sampling from the joint posterior distribution of initial
conditions and cosmological parameters. This will be explored in future work
once we move to more realistic models of the observables.

We now make a qualitative assessment of whether our posterior samples exploit
all the available information. As our work is the first to sample from the
posterior distribution of initial conditions using full N-body simulations, it
opens the door to using the full information content of future sky surveys for
inference. Figure 4 demonstrates that the posterior samples reproduce the true
power spectrum to an accuracy of better than 1% for all $k$ except within <10%
of the Nyquist frequency of the grid. The reconstructed fields account for 92%
of the variance of the true initial field; this level is expected given that
gravitational collapse destroys information within Lagrangian patches that
fall into halos. Using the Press-Schechter formalism and cutting off the power
spectrum at the grid frequency gives an estimated collapse fraction of 12% in
resolved halos. This is close to the measured fraction in the simulation of
16%. These calculations suggest that our samples saturate the information that
is physically available given non-linear gravitational collapse and coarse-
graining on the grid scale.

It bears mentioning that our analysis of DM density fields represents a more
stringent test of the network capability than what will be required for the
analysis of realistic observations, since these are far more sparsely sampled
and hence noisier than the DM density field.

A possible application of our methodology would be to run N-body simulations
from the samples of the inferred initial conditions; for existing samplers,
such as BORG, this is a natural byproduct (Leclercq et al., 2017). These
constrained realisations would open the possibility towards unlocking
important cosmological questions; we would obtain samples from the posterior
distribution of possible N-body simulations that result in the final
conditions. This would allow us to sample the posterior over properties of
galaxy clusters, such as the positions, masses, and velocities of all dark
matter halos in surveys. Furthermore, the possibility to sample constrained
high resolution N-body simulations would provide us with the means to uncover
the set of possible halo assembly histories. For example, this could be used
to sample possible histories of the Local Group and the Milky Way, as in
McAlpine et al. (2022). Initial condition inference with full N-body
simulations would help ensure accuracy since galaxy assembly histories are
sensitive to scales that are deeply in the non-linear regime today.

##  5 Conclusions

This work proposes score-based generative models for efficient sampling of the
posterior distribution of initial conditions, a problem that has up until now
been intractable using high-resolution N-body simulations. The key is to
sample by solving a reverse-diffusion process requiring only the gradient of
the noise-perturbed log posterior $\nabla_{\bm{x}}\log p_{t}(\bm{x}|\bm{y})$,
which can be learned using neural networks. The results show that we can
perform accurate inference of the initial conditions marginalized over
cosmological parameters at a fraction of the cost of state-of-the-art methods,
generating samples from the posterior within minutes on a single GPU. Our
tests indicate that the samples have the correct statistics at a level of
better than 1% across the relevant range of scales. In future work, we aim to
expand the proposed approach to infer initial conditions from simulated halo
and galaxy catalogs and ultimately apply it to real data.

## Acknowledgements

This work is supported by the Simons Collaboration on “Learning the Universe".
The Flatiron Institute is supported by the Simons Foundation. This work is
also supported by the NASA Research Opportunities in Space and Earth Science
(ROSES) program through grant number 12-EUCLID12-0004. The work is in part
supported by computational resources provided by Calcul Quebec and the Digital
Research Alliance of Canada. Y.H. and L.P. acknowledge support from the Canada
Research Chairs Program, the National Sciences and Engineering Council of
Canada through grants RGPIN-2020-05073 and 05102, and the Fonds de recherche
du Québec through grants 2022-NC-301305 and 300397. B.D.W. acknowledges
support by the ANR BIG4 project, grant ANR-16-CE23-0002 of the French Agence
Nationale de la Recherche; and the Labex ILP (reference ANR-10-LABX-63) part
of the Idex SUPER, and received financial state aid managed by the Agence
Nationale de la Recherche, as part of the programme Investissements d’avenir
under the reference ANR-11-IDEX-0004-02. R.L. is supported by a generous
donation by Eric and Wendy Schmidt with the recommendation of the Schmidt
Futures Foundation. R.L. thanks the Flatiron Institute for their hospitality
and the Centre for Research in Astrophysics of Quebec for their support. P.L
acknowledges support from the Simons Foundation. We thank Alexandre Adam for
reading the manuscript and providing important feedback.

## References

  * Acquaviva et al. (2002) Acquaviva V., Bartolo N., Matarrese S., Riotto A., 2002, arXiv preprint astro-ph/0209156
  * Adam et al. (2022) Adam A., Coogan A., Malkin N., Legin R., Perreault-Levasseur L., Hezaveh Y., Bengio Y., 2022, arXiv e-prints, p. arXiv:2211.03812
  * Albrecht & Steinhardt (1982) Albrecht A., Steinhardt P. J., 1982, Physical Review Letters, 48, 1220
  * Anand & Achim (2022) Anand N., Achim T., 2022, arXiv e-prints, p. arXiv:2205.15019
  * Bartolo et al. (2004) Bartolo N., Komatsu E., Matarrese S., Riotto A., 2004, Physics Reports, 402, 103
  * Batzolis et al. (2021) Batzolis G., Stanczuk J., Schönlieb C.-B., Etmann C., 2021, arXiv e-prints, p. arXiv:2111.13606
  * Brock et al. (2018) Brock A., Donahue J., Simonyan K., 2018, arXiv e-prints, p. arXiv:1809.11096
  * Dhariwal & Nichol (2021) Dhariwal P., Nichol A., 2021, arXiv e-prints, p. arXiv:2105.05233
  * Feldbrugge & van de Weygaert (2022) Feldbrugge J., van de Weygaert R., 2022, Cosmic Web ; Caustic Skeleton: non-linear Constrained Realizations – 2D case studies, doi:10.48550/ARXIV.2212.07840, https://arxiv.org/abs/2212.07840
  * Gnaneshwar et al. (2022) Gnaneshwar D., Ramsundar B., Gandhi D., Kurchin R., Viswanathan V., 2022, arXiv e-prints, p. arXiv:2203.04698
  * Guth (1981) Guth A. H., 1981, Physical Review D, 23, 347
  * Ho et al. (2020) Ho J., Jain A., Abbeel P., 2020, arXiv e-prints, p. arXiv:2006.11239
  * Hyvärinen (2005) Hyvärinen A., 2005, Journal of Machine Learning Research, 6, 695
  * Jasche & Lavaux (2019) Jasche J., Lavaux G., 2019, A&A, 625, A64
  * Jasche & Wandelt (2013) Jasche J., Wandelt B. D., 2013, Monthly Notices of the Royal Astronomical Society, 432, 894
  * Komatsu et al. (2005) Komatsu E., Spergel D. N., Wandelt B. D., 2005, The Astrophysical Journal, 634, 14
  * Lavaux & Wandelt (2010) Lavaux G., Wandelt B. D., 2010, MNRAS, 403, 1392
  * Leclercq et al. (2015) Leclercq F., Jasche J., Wandelt B., 2015, J. Cosmology Astropart. Phys., 2015, 015
  * Leclercq et al. (2017) Leclercq F., Jasche J., Lavaux G., Wandelt B., Percival W., 2017, J. Cosmology Astropart. Phys., 2017, 049
  * Legin et al. (2023) Legin R., Adam A., Hezaveh Y., Perreault Levasseur L., 2023, arXiv e-prints, p. arXiv:2302.03046
  * Lewis et al. (2000) Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473
  * Lin et al. (2016) Lin G., Milan A., Shen C., Reid I., 2016, arXiv e-prints, p. arXiv:1611.06612
  * Linde (1982) Linde A. D., 1982, Physics Letters B, 108, 389
  * Linde (1983) Linde A. D., 1983, Physics Letters B, 129, 177
  * Maldacena (2003) Maldacena J., 2003, Journal of High Energy Physics, 2003, 013
  * McAlpine et al. (2022) McAlpine S., et al., 2022, Monthly Notices of the Royal Astronomical Society, 512, 5823
  * Modi et al. (2021) Modi C., Lanusse F., Seljak U., Spergel D. N., Perreault-Levasseur L., 2021, arXiv e-prints, p. arXiv:2104.12864
  * Modi et al. (2023) Modi C., Li Y., Blei D., 2023, Journal of Cosmology and Astroparticle Physics, 2023, 059
  * Mudur & Finkbeiner (2022) Mudur N., Finkbeiner D. P., 2022, arXiv e-prints, p. arXiv:2211.12444
  * Müller-Franzes et al. (2022) Müller-Franzes G., et al., 2022, arXiv e-prints, p. arXiv:2212.07501
  * Papamakarios et al. (2017) Papamakarios G., Pavlakou T., Murray I., 2017, arXiv e-prints, p. arXiv:1705.07057
  * Planck Collaboration (2020a) Planck Collaboration 2020a, A&A, 641, A1
  * Planck Collaboration (2020b) Planck Collaboration 2020b, A&A, 641, A9
  * Popov et al. (2021) Popov V., Vovk I., Gogoryan V., Sadekova T., Kudinov M., 2021, arXiv e-prints, p. arXiv:2105.06337
  * Shallue & Eisenstein (2023) Shallue C. J., Eisenstein D. J., 2023, MNRAS, 520, 6256
  * Somerville & Davé (2015) Somerville R. S., Davé R., 2015, ARA&A, 53, 51
  * Song & Ermon (2019) Song Y., Ermon S., 2019, arXiv e-prints, p. arXiv:1907.05600
  * Song & Ermon (2020) Song Y., Ermon S., 2020, Advances in neural information processing systems, 33, 12438
  * Song et al. (2020) Song Y., Sohl-Dickstein J., Kingma D. P., Kumar A., Ermon S., Poole B., 2020, arXiv e-prints, p. arXiv:2011.13456
  * Song et al. (2021) Song Y., Shen L., Xing L., Ermon S., 2021, arXiv e-prints, p. arXiv:2111.08005
  * Springel (2005) Springel V., 2005, Monthly notices of the royal astronomical society, 364, 1105
  * Srivastava et al. (2014) Srivastava N., Hinton G., Krizhevsky A., Sutskever I., Salakhutdinov R., 2014, Journal of Machine Learning Research, 15, 1929
  * Villaescusa-Navarro et al. (2020a) Villaescusa-Navarro F., et al., 2020a, The Astrophysical Journal Supplement Series, 250, 2
  * Villaescusa-Navarro et al. (2020b) Villaescusa-Navarro F., et al., 2020b, ApJS, 250, 2
  * Vincent (2011) Vincent P., 2011, Neural Computation, 23, 1661
  * Yadav & Wandelt (2005) Yadav A. P., Wandelt B. D., 2005, Phys. Rev. D, 71, 123004

##  Appendix A Additional Figures

Figure 6: Summary statistics between generated posterior samples and ground
truth initial conditions for N-body simulations runned with different values
for the cosmological parameters. Note that the range of values for the
cosmological parameters is much larger than the constraints from current
observations. The results demonstrate that the score network is robust to
different cosmologies and provides accurate reconstructions across different
simulations.

Figure 7: Independent samples of initial conditions sampled from the posterior
$p(\bm{x}|\bm{y})$ based on the present-day density field $\bm{y}$ with
fiducial cosmology from Figure 2. Each sample shown here accurately represents
the true initial conditions, albeit with differences on the level of the
posterior variance shown in Figure 3.