# Measuring cluster masses with Planck and SPT-SZ CMB lensing

Alexandre Huchet  Université Paris-Saclay, CEA, IRFU, 91191, Gif-sur-Yvette,
France

###### Abstract

Galaxy cluster masses help to constrain cosmological parameters through the
halo mass function. To get rid of major biases in the mass measurement, we
directly probe the cluster gravitational potentials by observing their
gravitational lensing on the Cosmic Microwave Background (CMB). We measured
the average mass of a 468-cluster sample using SPT-SZ and Planck/HFI PR2 sky
maps and compared it with the masses provided in the SPT-SZ data set. We found
an average mass ratio of $\langle M_{\rm CMBlens}/M_{\rm SZ}\rangle=0.98\pm
0.19\rm{\leavevmode\nobreak\ (stat.)}\pm 0.03\rm{\leavevmode\nobreak\
(syst.)}$ [1], in agreement with SPT-SZ masses. We also showed that the CMB
large scales are important for the reconstruction of the lensing potential of
a cluster with a quadratic estimator.

##  1 Using CMB lensing to measure cluster masses

Galaxy clusters are the largest gravitationally bound structures. Their
distribution in mass and redshift is very sensitive to cosmological
parameters, particularly $\Omega_{\mathrm{m}}$ and $\sigma_{8}$, and can thus
be used to constrain them. Redshifts are obtained with optical or
spectroscopic surveys. Masses measurements can be obtained by various methods
and are subject to biases. Galaxy clusters actually contain only $\sim 3\%$ of
galaxies, $\sim 12\%$ of hot intracluster gas and $\sim 85\%$ of dark matter.
An optical or X-ray survey can map with good accuracy the galaxies or the hot
gas, but it will lack most of the mass, which is dark matter that does not
emit any light. To compensate this issue, one can use scaling relations that
link the mass measured through optical or X-rays to the total mass based on
hydro-dynamical simulations. One problem is that, to obtain a mass, some
assumptions need to be made, such as having a system in hydro-static
equilibrium – gravitation and pressure perfectly compensating each other.
These assumptions are idealistic and create biases in the final measurement.

To avoid those biases, one can choose to probe directly the gravitational
potentials of the clusters. This way, $100\%$ of the cluster mass can be
obtained by observing how the cluster distorts background images through
gravitational lensing. Two different types of background sources can be used:
the most used ones are galaxies, but one needs many galaxies to be able to
reconstruct the gravitational potential from statistical distortions of galaxy
images; the other one is the CMB. We chose the latter and focused on the
gravitational lensing of the CMB by galaxy clusters. Indeed, it gets harder to
find background galaxies at a redshift $z>1$, whereas the CMB was emitted at a
redshift $z\sim 1100$ and will always be behind the observed clusters.
Moreover, the signal-to-noise ratio (SNR) is expected to be constant in the
redshift range $0.5<z<2$ [2]. The gravitational lensing correlates Fourier
modes that were independent in the primary CMB. This is how we can manage to
retrieve the gravitational potential from this observable.

We used two complementary data sets in this analysis: SPT-SZ, a ground-based
telescope and Planck, a space-based telescope. On the one hand, the atmosphere
of the Earth does not allow the observation of the large scales of the CMB
because it fluctuates much faster than the telescope scanning speed. Planck
can provide these large scales, whereas the SPT-SZ data suffers from an
anisotropic transfer function that describes which CMB scales are lost. On the
other hand, it is hard to send large mirrors in space that would allow for the
observation of the CMB small scales. SPT-SZ and its 10-meter-diameter mirror
will allow to observe these scales. We used the 3 SPT-SZ and 6 Planck/HFI
frequency maps (95, 150 and 220 GHz for SPT-SZ and 100, 143, 217, 353, 545 and
857 GHz for Planck/HFI) and studied a sample of 468 clusters from SPT-SZ [3]
in the redshift range $0<z<1.7$.

##  2 Analysis

We ran the analysis pipeline on simulations of the 9 frequency maps and on
real data. For real data, we cut $10\times 10\leavevmode\nobreak\ \deg^{2}$
tangential maps from the All-sky maps for each frequency, and we directly
build flat maps for the simulations.

###  2.1 Simulations

The simulations were kept simple, as a mean to detect possible biases in our
analysis in a case with simple foregrounds. The first step was to create a
$10\times 10\leavevmode\nobreak\ \deg^{2}$ CMB map in $\Delta T/T$ units, from
the Planck PR2 CMB TT spectrum. This map was first periodic, i.e. created in
Fourier space and then transferred in real space. To make it more realistic,
we also did the same with a $50\times 50\leavevmode\nobreak\ \deg^{2}$ map and
then cut a non periodic $10\times 10\leavevmode\nobreak\ \deg^{2}$ map out of
it. We then added the gravitational lensing effect from a cluster with a
Navarro-Frenk-White (NFW) profile and its associated Sunyaev-Zel’dovich (SZ)
effect for each frequency (for a total of 9 maps). The latter is how the CMB
spectrum is modified because of the inverse Compton scattering of CMB photons
by hot intracluster free electrons. It can have a thermal (tSZ) or kinetic
origin (kSZ). Since kSZ has the same signature as the CMB, one can only
separate the tSZ. The latter depends on the integrated electronic pressure,
that we modelled with a generalised NFW profile, whereas the former was added
only in systematic checks for cluster peculiar velocities of $v\sim 300$
km.s-1. We then added the instrumental characteristics for each frequency,
i.e. the beam, transfer function and noise.

###  2.2 Pipeline

For each cluster (real or simulated), we studied flat $10\times
10\leavevmode\nobreak\ \deg^{2}$ maps at 9 different frequencies. Having the
same patch of sky at different frequencies helps to clean it from foregrounds
and noise. We used a constrained Internal Linear Combinations method [4]
(cILC) to get rid of the tSZ effect as well as mitigate the other
contaminants. This method takes advantage of the fact we know the frequency
signature of the tSZ effect as well as our instrumental characteristics. This
way, we built the best lensed CMB map we could from our data. Then, we looked
for a lensing signal at the cluster position using a quadratic estimator [5]
that extracts the lensing potential from the correlations between Fourier
modes in the CMB. We then have information on a Fourier mode
$\bm{K}=\bm{k^{\prime}}-\bm{k}$ from the correlations between
$\bm{k^{\prime}}$ and $\bm{k}$. We finally compared the obtained lensing
potential with the fiducial potential of a cluster with a NFW profile and the
mass given in the SPT-SZ data set. To do so, we used a matched filter [2] that
compares mode by mode the two profiles in Fourier space. The ratio of the
profiles is equal to the ratio of the masses.

The real maps required a few additional steps: the sky map was projected on a
tangential map using Healpy, the non periodic maps were apodised (the borders
were smoothed out, also for the non periodic simulations) and the point
sources were masked using a constrained Gaussian field [6] with the CMB
properties and local continuity. Moreover, to subtract systematic biases due
to the foregrounds, correlated noise or non periodicity, we ran the analysis
pipeline on 10 empty field (_off_ map) for each field with a cluster (_on_
map). The final result is _on_ \- $\langle\emph{off}\rangle_{10}$.

##  3 Results

We obtained consistent results for Planck only, SPT-SZ only and the
combination of both on periodic and non periodic simulations, and on real data
[1]. The combined result on real data is shown in Fig.1. There is an overall
good agreement between our observed CMB lensing mass and the mass derived by
SPT-SZ based on the tSZ signal. We observe a small increasing trend in the
ratio between our mass measurement and the normalized SPT-SZ measurement with
respect to the latter. This trend may be real or could be due to some
remaining bias in our analysis. After correcting for the bias due to the kSZ
estimated from the simulations ($\Delta\langle M_{\rm CMBlens}/M_{\rm
SZ}\rangle=-0.06$), we find a final average mass ratio of $\langle M_{\rm
CMBlens}/M_{\rm SZ}\rangle=0.98\pm 0.19\rm{\leavevmode\nobreak\ (stat.)}\pm
0.03\rm{\leavevmode\nobreak\ (syst.)}$ [1].

Figure 1: Recovered CMB lensing mass $M_{\rm CMBlens}$ (divided by the SZ mass
$M_{\rm SZ}$) as a function of SZ mass for the Planck \+ SPT-SZ combined
analysis. Black dots are individual measurements. The red dots are the
weighted means of the five mass bins, spaced logarithmically. The dotted
(resp. dashed) blue line is the zero (resp. unity) level. The green stripe is
the weighted mean within the error $\langle M_{\rm CMBlens}/M_{\rm
SZ}\rangle=0.919\pm 0.190$ (without correcting for the kSZ effect).

We also observe that the statistical error bars gain a factor of $\sqrt{2}$
between one data set ($\sigma_{\textit{Planck}}=0.27$, $\sigma_{\rm SPT-
SZ}=0.29$) and the combination of both ($\sigma_{\rm combi}=0.19$). It is not
obvious that two data sets bring independent information when they observe the
same sky. We show in Fig 2 the ratio $\sigma_{\rm combi}^{2}/\sigma_{i}^{2}$
with respect to the Fourier mode $K$, with $i$ representing either Planck or
SPT-SZ. We also plotted the sum of both ratios for comparison with purely
independent data sets. The figure shows that the Planck data set brings most
of the information at large scales ($K<50.9$ or $L<1830$), while the SPT-SZ
dataset provides information on the smaller scales ($K>145.0$ or $L>5208$).
The sum is the product between the combined $\sigma_{\rm SPT+Planck}^{2}$ term
and the inverse variance sum $1/\sigma_{\rm SPT}^{2}+1/\sigma_{\rm
Planck}^{2}$. It is below unity for $K>18.9$ or $L>680$ demonstrating the
significant gain of the combination with respect to individual measurements.

We therefore showed that information on CMB large scales is useful in the
reconstruction of the lensing potential of a cluster using a quadratic
estimator on temperature maps. In the future, Planck and/or the Small Aperture
Telescopes (SATs) data might complement the Large Aperture Telescopes (LATs)
on this matter.

Figure 2: Comparison of the variance $\sigma^{2}(K)$ of the lensing profile
reconstructed for SPT-SZ and Planck binned in $K=|\bm{K}|$ (or equivalently
$L=|\bm{l^{\prime}}-\bm{l}|$). The figure displays the ratios between the
combined $\sigma_{\rm SPT+Planck}^{2}(K)$ and the single datasets
$\sigma_{i}^{2}(K)$, and the sum of both ratios. See the main text for more
information.

## References

## References

  * [1] Alexandre Huchet and Jean-Baptiste Melin.  A measurement of cluster masses using Planck and SPT-SZ CMB lensing.  2024\.  Submitted to A&A.
  * [2] Jean-Baptiste Melin and James G. Bartlett.  Measuring cluster masses with CMB lensing: a statistical approach.  Astronomy & Astrophysics, 578:A21, June 2015.  arXiv:1408.5633 [astro-ph].
  * [3] S. Bocquet et al.  Cluster Cosmology Constraints from the 2500 deg2 SPT-SZ Survey: Inclusion of Weak Gravitational Lensing Data from Magellan and the Hubble Space Telescope.  The Astrophysical Journal, 878(1):55, June 2019.  Publisher: The American Astronomical Society.
  * [4] Mathieu Remazeilles, Jacques Delabrouille, and Jean-Francois Cardoso.  CMB and SZ effect separation with Constrained Internal Linear Combinations.  Monthly Notices of the Royal Astronomical Society, 410(4):2481–2487, February 2011.  arXiv:1006.5599 [astro-ph].
  * [5] Wayne Hu and Takemi Okamoto.  Mass Reconstruction with CMB Polarization.  The Astrophysical Journal, 574(2):566–574, August 2002.  arXiv:astro-ph/0111606.
  * [6] Yehuda Hoffman and Erez Ribak.  Constrained Realizations of Gaussian Fields: A Simple Algorithm.  The Astrophysical Journal, 380:L5, October 1991.  ADS Bibcode: 1991ApJ…380L…5H.


Case (i): Suppose $\mathbf{A}$ allocates the first arriving user to $a_2$. Then consider the following sequence, where user one arrives at time $0$, user two arrives at time $10$, and user three arrives at time $10.001$. Note that an optimal allocation for this sequence is to allocate user one on $a_1$, user two on $a_1$, and user three on $a_2$. Since user one departs from the system before user two arrives, the optimal allocation does not allocate two users concurrently on $a_1$. However, since $\mathbf{A}$ allocates user one on $a_2$, it must hold that both users two and three are allocated on $a_1$. To see this, we note that at the times users two and three arrive, user one is still on $a_2$ since the travel time of $a_2$ is strictly greater than the arrival time of user three. Since the capacity of $a_2$ is one, it must hold that $\mathbf{A}$ allocates both users two and three on $a_1$. Since $\Bar{d}$ can be arbitrarily large, we have that the ratio of the cost of such a deterministic algorithm to that of the optimal offline solution is unbounded.

Case (ii): Suppose $\mathbf{A}$ allocates the first arriving user to $a_1$. Then consider the following input sequence, where user one arrives at time $0$, user two arrives at time $0.015$, user three arrives at time $5.02$, and user four arrives at time $10.01$. Note that an optimal allocation for this input sequence is to allocate user one on $a_2$, users two and three on $a_1$, and user four on $a_2$. Again, as with case (i), the optimal allocation does not allocate two users concurrently on $a_1$. However, since $\mathbf{A}$ allocates user one on $a_1$, it must be that at least two of the last three users are concurrently allocated on $a_1$. To see this, we note that if $\mathbf{A}$ allocates the second arriving user on $a_1$, then we have two users concurrently occupying $a_1$ thereby incurring an arbitrarily large cost. Thus, we consider the case when $\mathbf{A}$ allocates user two to $a_2$. Now, we observe that both users three and four must be allocated to $a_1$ since both users' arrival time is before the departure of user two from $a_2$. Since the capacity of $a_2$ is one, it must be that $\mathbf{A}$ allocates both of the last two users on $a_1$. Thus, again we have that the ratio between the total travel cost of the online solution of $\mathbf{A}$ to that of the optimal allocation is unbounded since $\Bar{d}$ can be arbitrarily large.

§ GENERALIZATIONS OF OTR PROBLEM

§.§ Extending otr Problem to General (Non-Parallel) Road Networks

In this section, we present a generalization of the otr problem for parallel networks to the setting of general road networks. To model a general road network, we consider a directed graph $G = (V, E)$, where $V$ is its set of vertices and $E$ is its set of edges. Travelling through arc $a \in E$ requires $t_a$ time units, and each arc $a$ has a capacity $c_a$. We collect these parameters in $\mathcal{S} := \{ \{ c_a \}_{a \in E}, \{ t_a \}_{a \in E} \}$. We denote an input sequence of $n$ users by $\mathcal{I} := \{ \{ \theta_i \}_{i \in [n]}, \{ \tau_i \}_{i \in [n]}, \{ w_i \}_{i \in [n]} \}$, stating for each user $i$ its time of appearance $\tau_i\in\mathbb{R}_{\ge0}$, its value-of-time $\theta_i \in \Theta$, where $\Theta$ is finite, and its origin destination pair $w_i \in W$. Each origin destination pair $w \in W$ has $R_w$ routes and each route $r \in R_w$ is specified by an ordered sequence of $l$ edges $\{e_1^r, e_2^r, ..., e_l^r \}$. Here edge $e_1^r$ is a directed edge starting at the origin node of $w$ and edge $e_l^r$ is a directed edge ending at the destination node of $w$. With this notation, the otr (otr) problem for general road networks is as follows.

(otr (otr) for General Road Networks) Let the parameter set $\mathcal{S}$ and an input sequence $\mathcal{I}$ define a problem instance. Users arrive sequentially in the system and must be irrevocably assigned to a route $r \in R_{w_i}$ at their time of appearance $\tau_i$. The number of users concurrently traversing an arc is limited to the capacity $c_a$ of the arc, and each user traverses an arc $a$ for $t_a$ time units. Then, the otr problem for general road networks is to construct an assignment of users to routes such that the total cost or time that users spend to traverse their routes is minimal.

If all users are identical, i.e., $\theta_i = \theta_{i'}$ for each $i, i' \in [n]$, we refer to this problem variant as otr-i for general road networks.

We evaluate the efficacy of an algorithm $\mathbf{A}$ that solves Problem <ref> via its competitive ratio

$$CR(\mathbf{A}) = \max_{\mathcal{I} \in \Omega} \frac{ALG(\mathcal{I}, \mathbf{A})}{OPT(\mathcal{I})},$$

and say that an algorithm is $\alpha$-competitive if $CR(\mathbf{A}) \leq \alpha$. Here $OPT(\mathcal{I})$ denotes the optimal offline solution to the otr problem for general road networks and $ALG(\mathcal{I}, \mathbf{A})$ denotes the online solution achieved by algorithm $\mathbf{A}$ on input sequence $\mathcal{I}$. The input instances $\mathcal{I}$ belong to a set $\Omega$, which corresponds to user arrivals such that $\theta_i \geq 0$ for all users $i$ and $\tau_i \geq 0$, with the first user arriving at time $0$. We further assume without loss of generality that $\tau_1 = 0 < \tau_2 < ... < \tau_n $, i.e., no two users arrive in the system at the same point in time.

We compute $OPT(\mathcal{I})$ by solving the following optimization problem.
\beginmini!|s|[2]<b>
∀i ∈[n], r ∈{R_w}_w ∈W∑_i = 1^n ∑_a ∈E x_ir θ_i t_a δ_a, r

∑_r = 1^|R_w| x_ir=1   ∀i ∈[n]
x_ir∈{0, 1},   ∀r ∈[|R_w_i|, ∀i ∈[n] ]
∑_i=1^n ∑_r:a ∈r 1_τ_k ∈[τ_i + ∑_a' ∈r:a' < a t_a', τ_i + ∑_a' ∈r:a' ≤a t_a'] x_ir≤c_a,   ∀a ∈E, ∀τ_k ∈{τ_1, ..., τ_n }
\endmini!

\begin{align}
\min & \sum_{i = 1}^{n} \sum_{a \in E} x_{ir} \theta_i t_a \delta_{a, r}\label{eq:OPT-integer-Obj-general-road}\\
&\sum_{r = 1}^{|R_w|} x_{ir} = 1, \quad \forall i \in [n] \label{eq:prob-constraint-general-road} \\
&x_{ir} \in \{0, 1 \}, \quad \forall i \in [n], \forall r \in [|R_w|] \label{eq:binary-constraint-general-road}  \\
&\sum_{i=1}^n \sum_{r:a \in r} \mathbbm{1}_{\tau_{k} \in [\tau_{i} + \sum_{a' \in r:a' < a} t_{a'}, \tau_{i} + \sum_{a' \in r:a' \leq a} t_{a'}]} x_{ir} \leq c_a, \\ \forall a \in E, \forall \tau_k \in \{\tau_1, ..., \tau_n \} \label{eq:capacity-constraints-general-road}
\end{align}

Here we introduce binary variables $x_{ir}$ to denote whether user $i$ is assigned to route $r \in R_{w_i}$ ($x_{ir} = 1$) or not ($x_{ir}=0$). Further, we assume the constraints of this problem are feasible where, (<ref>) are assignment constraints, (<ref>) are binary constraints and (<ref>) are capacity constraints, and $\delta_{a, r}$ is the indicator function stating whether arc $a$ belongs to route $r$ ($\delta_{a, r}=1$) or not ($\delta_{a, r}=0$). In the capacity constraints, $\sum_{a' \in r:a' < a} t_{a'}$ denotes the sum of the travel times on edges preceding edge $a$ on route $r$, and $\sum_{a' \in r:a' \leq a} t_{a'}$ is the sum of the travel times on edges preceding as well as including edge $a$ on route $r$.

§.§ Extending OTR Problem to General Cost Functions

In this section, we present a generalization of the otr problem to the setting when the travel time of the arcs is modelled by a cost function. As with the otr problem, we consider a graph with $M \geq 2$ parallel arcs and two vertices representing an origin and a destination. Travelling through arc $a$ requires $f_a(y_a)$ time units, when $y_a$ users use arc $a$, and each arc $a$ has a capacity $c_a$. Here $f_a:\mathbb{R}_{\geq 0} \mapsto \mathbb{R}_{\geq 0}$ is non-decreasing in its argument. We collect these parameters in $\mathcal{S} := \{ \{ c_a \}_{a \in [M]}, \{ f_a \}_{a \in [M]} \}$. We denote an input sequence of $n$ users by $\mathcal{I} := \{ \{ \theta_i \}_{i \in [n]}, \{ \tau_i \}_{i \in [n]} \}$, stating for each user $i$ its time of appearance $\tau_i\in\mathbb{R}_{\ge0}$ and its value-of-time $\theta_i \in \Theta$, where $\Theta$ is finite. With this notation, the otr problem with general cost functions is as follows.

(Online Traffic Routing (otr) with General Cost Functions) Let the parameter set $\mathcal{S}$ and an input sequence $\mathcal{I}$ define a problem instance. Users arrive sequentially in the system and must irrevocably be assigned to an arc at their time of appearance $\tau_i$. The number of users concurrently traversing an arc is limited to the capacity $c_a$ of the arc, and each user traverses an arc $a$ for $f_a(y_a+1)$ time units, where $y_a$ users are using arc $a$ at time $\tau_i$. Then, the otr problem for general cost functions is to construct an assignment of users to arcs such that the total cost or time that users spend to traverse an arc is minimal.

If all users are identical, i.e., $\theta_i = \theta_{i'}\ \forall\ i, i' \in [n]$, we refer to this problem variant as otr-i with general cost functions.

We evaluate the efficacy of an algorithm $\mathbf{A}$ that solves Problem <ref> via its competitive ratio

$$CR(\mathbf{A}) = \max_{\mathcal{I} \in \Omega} \frac{ALG(\mathcal{I}, \mathbf{A})}{OPT(\mathcal{I})},$$

and say that an algorithm is $\alpha$-competitive if $CR(\mathbf{A}) \leq \alpha$. Here $OPT(\mathcal{I})$ denotes the optimal offline solution to the otr problem for general cost functions and $ALG(\mathcal{I}, \mathbf{A})$ denotes the online solution achieved by algorithm $\mathbf{A}$ on input sequence $\mathcal{I}$. The input instances $\mathcal{I}$ belong to a set $\Omega$, which corresponds to user arrivals such that $\theta_i \geq 0$ for all users $i$, and $\tau_i \geq 0$, with the first user arriving at time $0$. We further assume without loss of generality that $\tau_1 = 0 < \tau_2 < ... < \tau_n $, i.e., no two users arrive in the system at the same point in time, and that there exists a feasible assignment.

We compute $OPT(\mathcal{I})$ by solving the following optimization problem.
\beginmini!|s|[2]<b>
∀i ∈[n], a ∈[M]∑_i = 1^n ∑_a = 1^M x_ij θ_i f_a (x_ia + y_a(τ_i) )

∑_a = 1^M x_ia=1   ∀i ∈[n]
x_ia∈{0, 1},   ∀i∈[n], a ∈[M]
∑_i=1^n 1_τ_k ∈[τ_i, τ_i + f_j (x_ia + y_a(τ_i) )] x_ia≤c_a,   ∀a ∈[M], ∀τ_k ∈{τ_1, ..., τ_n }
y_a(τ_i)= ∑_k=1^n 1_τ_k < τ_i, τ_k+f_a(x_ka+y_a(τ_k))>τ_i  x_ia,   ∀a ∈[M], ∀τ_k ∈{τ_1, ..., τ_n }
\endmini!

\begin{align}
\min & \sum_{i = 1}^n \sum_{a = 1}^M x_{ij} \theta_i f_a \left(x_{ia} + y_a(\tau_i) \right) \label{eq:OPT-integer-Obj-general-congestion}\\
&\sum_{a = 1}^{M} x_{ia} = 1, \quad \forall i \in [n] \label{eq:prob-constraint-general-congestion} \\
&x_{ia} \in \{0, 1\}, \quad \forall i\in [n], a \in [M] \label{eq:binary-constraint-general-congestion}  \\
&\sum_i \mathbbm{1}_{\tau_{k} \in [\tau_{i}, \tau_{i} + f_j \left(x_{ia} + y_a(\tau_i) \right)]} x_{ia} \leq c_a, \quad \forall a \in [M], \forall \tau_k \in \{\tau_1, ..., \tau_n \} \label{eq:capacity-constraints-general-congestion} \\
& y_a(\tau_i) = \sum_{k} \mathbbm{1}_{\tau_{k} < \tau_i, \tau_k+f_a(x_{ka}+y_a(\tau_{k}))>\tau_i } x_{ia}, \quad \forall a \in [M], \forall \tau_k \in \{\tau_1, ..., \tau_n \} \label{eq:recursive-eqn}
\end{align}

Here, we introduce binary variables $x_{ia}$ to denote whether user $i$ is assigned to arc $a$ ($x_{ia} = 1$) or not ($x_{ia}=0$). Further, (<ref>) are assignment constraints, (<ref>) are binary constraints and (<ref>) are capacity constraints. Equations (<ref>) represent the total number of users on each arc $a$ at each time $\tau_i$ a user arrives.


# Recommender System Expedited Quantum Control Optimization

Priya Batra<EMAIL_ADDRESS>M. Harshanth Ram
<EMAIL_ADDRESS>T. S. Mahesh
<EMAIL_ADDRESS>Department of Physics and NMR Research Center,
Indian Institute of Science Education and Research, Pune 411008, India

###### Abstract

Quantum control optimization algorithms are routinely used to generate optimal
quantum gates or efficient quantum state transfers. However, there are two
main challenges in designing efficient optimization algorithms, namely
overcoming the sensitivity to local optima and improving the computational
speed. The former challenge can be dealt with by designing hybrid algorithms,
such as a combination of gradient and simulated annealing methods. Here, we
propose and demonstrate the use of a machine learning method, specifically the
recommender system (RS), to deal with the latter challenge of enhancing
computational efficiency. We first describe ways to set up a rating matrix
involving gradients or gate fidelities. We then establish that RS can rapidly
and accurately predict elements of a sparse rating matrix. Using this
approach, we expedite a gradient ascent based quantum control optimization,
namely GRAPE and demonstrate the faster performance for up to 8 qubits.
Finally, we describe and implement the enhancement of the computational speed
of a hybrid algorithm, namely SAGRAPE.

##  I Introduction

Quantum control is crucial for the trending quantum technology tasks such as
quantum sensing [1], scalable quantum information devices [2], quantum
simulations [3], quantum thermodynamics [4], quantum metrology [5], etc. In
general, quantum control optimization deals with finding the best
implementation of desired quantum dynamics on a given physical hardware [6,
7]. Such optimization algorithms are routinely used for control tasks like
unitary synthesis, state transfer, etc [8]. There has been a tremendous
progress in this area and numerous optimization algorithms have been
developed, such as gradient based algorithms [9, 10], variational principle
based algorithms [11, 12], chopped basis optimization [13, 14], and
metaheuristic algorithms [15, 16]. Lately, machine learning algorithms such as
reinforcement learning (RL), have also been used for the tasks like unitary
transformation [17], state preparation [18], robust controls [19], as well as
to control non-integrable quantum systems [20]. More recently, machine
learning protocols have also been used to control quantum thermal machines
[21].

An optimization method with an analytical expression for the gradient is an
efficient way to find a local optimum in the parameter space. One such class
of optimization algorithms for quantum control applications is based on
gradient ascent pulse engineering (GRAPE) [9] and its variants (for example,
[22, 23, 24]). However, it becomes problematic if there are too many sub-
optimal local solutions restricting the algorithm from reaching an optimal
solution. On the other hand, metaheuristic algorithms can escape the local
optima and therefore perform a robust search in the parameter space. The
metaheuristic algorithms such as the Nelder-Mead simplex method [15],
evolutionary algorithms [25, 16, 26, 27], simulated annealing (SA) [28], etc.
have been successfully adapted for quantum control optimization. However,
unlike the gradient methods, the metaheuristic algorithms often suffer from
slow convergence. To overcome this issue, recently there has been a proposal
to combine SA, a metaheuristic method with GRAPE, a local search method, to
realize the hybrid SAGRAPE algorithm [28]. Despite this progress, there
remains the problem of the poor computational efficiency of quantum
optimization algorithms. This problem becomes severe as the parameter space
grows exponentially with the size of the quantum system. This is where the
capabilities of machine learning can be useful (see Fig. 1).

Figure 1: Expediting gradient method, metaheuristic method, as well as a
hybrid of both by using machine learning.

In the field of machine learning, the recommender system (RS) algorithm is
widely used to recommend products to consumers [29, 30]. Here we propose RS
assisted speed up of GRAPE and SAGRAPE. The most time-expensive task in GRAPE
is the calculation of gradients for every segment of a control sequence, which
involves exponentiation of matrices. We show that, given a set of exactly
calculated gradients, RS can accurately and rapidly predict the remaining
ones. Thus, we can significantly improve the time-efficiency of the
optimization algorithm. SA involves a large number of function evaluations in
the neighborhood of an iterative solution, part of which can also be delegated
to RS. Here we demonstrate a significant time advantage of RS expedited GRAPE
as well as SAGRAPE, without sacrificing their convergence.

This article is arranged as follows. In Sec. II, we introduce the quantum
control problem along with two optimization algorithms, namely GRAPE and SA.
In Sec. III we introduce RS. In Sec. IV, we present the RS expedited GRAPE and
SAGRAPE, and discuss their performances. The article ends with a summary and
outlook in Sec. V.

##  II Quantum Control Optimization

We consider a quantum system with a constant internal Hamiltonian $H_{0}$ and
$m$ control Hamiltonians $\\{H_{j}\\}$, such that the time-dependent
Hamiltonian is given by

$H(t)=H_{0}+\sum_{j=1}^{m}\omega_{j}(t)H_{j}.$ (1)

Here $\omega_{j}(t)$ is the time-dependent strength of the $j$th control
Hamiltonian $H_{j}$. For a closed system undergoing Schrödinger evolution for
the time duration $T$, the corresponding unitary is given by

$U(T)=De^{-i\int_{0}^{T}H(t)dt}.$ (2)

Here $D$ is the Dyson time ordering operator and we have set $\hbar=1$.
Considering the difficulty in evaluating the propagator of a general time-
dependent Hamiltonian, we discretize the control function $\omega_{j}(t)$ by
dividing it into $N$ segments, each with a constant amplitude $\omega_{j,k}$
and duration $\delta=T/N$. The corresponding unitary for the $k$th segment
would be

$U_{k}=\exp\left[-i\delta\left(H_{0}+\sum_{j=1}^{m}\omega_{j,k}H_{j}\right)\right].$
(3)

The total time propagator $U(T)$ for the entire control sequence can be
expressed as the product of segment unitaries

$U(T)=U_{N}U_{N-1}\cdots U_{2}U_{1}.$ (4)

In this article, we mainly focus on unitary synthesis, although same methods
can be adapted for state preparation. Thus, the optimization procedure is
aimed at achieving a specific target unitary $U_{t}$ by numerically generating
the control sequence $\\{\omega_{j,k}\\}$. Optimization function is given by
the gate fidelity

$\displaystyle
F=\left|\mathrm{Tr}\left[U^{\dagger}(T)U_{t}\right]\right|^{2},$ (5)

which is the overlap of target unitary $U_{t}$ with the evolved unitary
$U(T)$. An improved optimization function would maximize the fidelity while
minimizing resources. Typically, it amounts to minimizing the power
consumption of control fields. Therefore the modified optimization function
$J$ can be cast as follows

$\displaystyle J=F-\sum_{j,k}\lambda_{j}\omega_{j,k}^{2},$ (6)

where $\lambda_{j}$ are the scalar penalty parameters. In the following, we
discuss two optimization methods, first a gradient method and the second a
metaheuristic method.

###  II.1 Gradient Ascent Pulse Engineering (GRAPE)

The GRAPE algorithm has following steps:

  * •

Start with a random control sequence $\\{\omega_{j,k}^{(0)}\\}$.

  * •

Forward propagate the initial unitary operator $U_{0}=\mathbbm{1}$ till $k$th
segment, i.e., $X_{k}=U_{k}U_{k-1}\cdots U_{1}U_{0}$.

  * •

Backward propagate the target unitary $U_{t}$ till $k$th segment, i.e.,
$P_{k}=U_{k+1}^{\dagger}U_{k+2}^{\dagger}\cdots U_{N}^{\dagger}U_{t}$.

  * •

Calculate the gradient for each segment using the first-order (in the norm
$||\delta H(t)||$) expression [9]

$\displaystyle
g_{j,k}=2\mathrm{Re}\left(-i\delta\mathrm{Tr}\left[P_{k}^{\dagger}H_{j}X_{k}\right]\mathrm{Tr}\left[X_{k}^{\dagger}P_{k}\right]\right).$
(7)

  * •

Update control amplitudes in the direction of gradients, i.e.,
$\omega_{j,k}\rightarrow\omega_{j,k}+\epsilon g_{j,k}$, where $\epsilon$ is
the step size.

GRAPE is remarkable to have such a simple analytical form for the gradient
function.

In a practical scenario, no physical hardware is perfect. For example, a
control field is typically associated with a distribution of amplitudes around
the nominal value [9]. We need to have a robust quantum control even with such
imperfect hardware. To this end, one minimizes the overall cost function
$1-\overline{J}=1-\sum_{i}p_{i}J_{i}$ obtained by summing over the costs
$1-J_{i}$ of individual elements in the distribution with respective
probabilities $p_{i}$.

###  II.2 Simulated Annealing (SA)

In metallurgy, annealing involves heating a material to high temperatures
followed by slow cooling, to allow the material reach a stable crystalline
form by finding its ground state. The same idea is adopted in the numerical
procedure namely, simulated annealing (SA) [31, 32]. Given an optimization
problem, SA starts with a high-temperature configuration, wherein even
nonoptimal solutions are selected with a certain probability. As the
iterations pass, the temperature parameter is gradually reduced, and optimal
solutions are increasingly favored. This stochastic process allows the
algorithm to overcome the local minima and reach the global minimum.

The steps for SA are as follows [28]:

  * •

Start with a random control sequence $\\{\omega_{j,k}^{(0)}\\}$. Set the
temperature to a high value $T^{0}$.

  * •

In the $i$th iteration, near the current solution $\\{\omega_{j,k}^{(i)}\\}$,
determine the control sequence $\\{\omega_{j,k}^{\prime}\\}$ with the best
fidelity $\overline{F}(\\{\omega_{j,k}^{\prime}\\})$ among a random set of
neighbourhood points.

  * •

If
$\Delta\overline{F}^{i}=\overline{F}(\\{\omega_{j,k}^{\prime}\\})-\overline{F}(\\{\omega_{j,k}^{(i)}\\})\geq
Q^{i}$, where
$Q^{i}=-\mathrm{min}\left[1,T^{i}\exp\left(\Delta\overline{F}^{i}/T^{i}\right)\right]$
then
$\\{\omega_{j,k}^{(i+1)}\\}=\\{\omega_{j,k}^{\prime}\\}$; else,
$\\{\omega_{j,k}^{(i+1)}\\}=\\{\omega_{j,k}^{(i)}\\}$.

  * •

After $i$th iteration, set the temperature to a reduced value
$T^{(i+1)}=\gamma T^{(i)}$, where $\gamma<1$ controls the cooling rate.

The above steps are iterated until the desired optimization function is
reached or the maximum number of iterations are completed. Note that, for
higher temperatures the algorithm may take a new solution
$\\{\omega_{j,k}^{\prime}\\}$ even if its fidelity is lower than the current
solution $\\{\omega_{j,k}^{(i)}\\}$. This stage is known as exploration. As
temperature goes down, the algorithm gradually switches to the exploitation
mode and it looks for solutions that are either better or slightly inferior
than the current solution. The combination of exploration and exploitation
helps SA to escape local optima and travel towards the global optimum.

In our previous work, we had combined SA with GRAPE to form the hybrid SAGRAPE
algorithm and demonstrated its superior convergence [28]. In the following, we
first describe recommender system, a particular type of machine learning
technique, and then explain how it can be used to expedite GRAPE as well as
SAGRAPE.

##  III Recommender System (RS)

Collaborative filtering is one of the most popular types of RS that is based
on experience of any particular consumer along with relative preferences among
all consumers [33, 34]. Here we use the matrix factorization algorithm for
collaborative filtering [35, 36]. It involves setting up a database ${\cal R}$
in the form of a rating matrix, wherein each row represents a particular
consumer and each column represents a particular product that is being
recommended [30, 37]. The rating matrix can be decomposed in terms of latent
vectors of the same dimension $f$, known as the number of features. Let the
parameter vector $\Theta^{(i)}\in\mathbb{R}^{f}$ and the feature vector
$\Phi^{(j)}\in\mathbb{R}^{f}$ encode latent vectors in real space
$\mathbb{R}^{f}$ for $i$th consumer and $j$th product. The predicted rating is
then modeled by scalar products

$r_{i,j}=\Theta^{(i)}\cdot\Phi^{(j)}.$ (8)

Depending on the problem at hand, one can specify products, users, as well as
ratings, and accordingly set up the rating matrix. One such example, for a
hypothetical quantum control problem of executing certain specific tasks with
various available control fields, is illustrated in Fig. 2. In this example,
tasks are users, control fields are products, and the ratings are different
levels of feasibility of implementing tasks. The job of RS is to efficiently
predict the unknown ratings.

Figure 2: Illustrating an example application of RS. Here the rating matrix
corresponds to different levels of feasibility of various tasks with certain
available control fields. The goal of RS is to predict the unknown ratings.

Given a sparse rating matrix ${\cal R}$, where missing elements correspond to
unknown preferences, our goal is to fill those missing elements with the help
of the collective information embedded in the overall database. Let
$\kappa={(i,j)}$ be the elements in the rating matrix for which actual ratings
${\cal R}_{i,j}$ are known. The discrepancy between predicted rating $r_{i,j}$
and actual rating ${\cal R}_{i,j}$ is quantified by the function $K_{0}$ as

$K_{0}=\sum_{(i,j)\in\kappa}(r_{i,j}-{\cal R}_{i,j})^{2}.$ (9)

Here, the objective is to minimize the cost function $K_{0}$. Generally two
regularization parameters $(\Lambda_{\Theta},\Lambda_{\Phi})$ are used to
avoid over-fitting, so that the altered cost function can be written as

$K=\frac{K_{0}}{2}+\frac{\Lambda_{\Theta}}{2}\sum_{i=1}^{m}\|\Theta^{(i)}\|+\frac{\Lambda_{\Phi}}{2}\sum_{j=1}^{n}\|\Phi^{(j)}\|.$
(10)

The latent vectors $\Theta^{(i)}$ and $\Phi^{(j)}$ are then determined by
minimizing the cost function using any standard minimization algorithm. In our
case, we use the Polack-Ribiere flavour of conjugate gradients to compute
search directions, and a line-search using quadratic and cubic polynomial
approximations for this purpose [38, 39, 40].

##  IV RS enhanced GRAPE and SAGRAPE

We shall now explain how we can use RS to speedup GRAPE as well as SAGRAPE.
For the sake of clarity, we use the context of quantum control of spin-
dynamics via nuclear magnetic resonance (NMR). To study RS enabled speedup, we
consider a spin system of $M=2+n$ on-resonant heteronuclear qubits, and
construct a controlled-NOT (CNOT) gate on the first two qubits which are
coupled by an indirect spin-spin interaction of strength ($J_{12}$). For
simplicity, we treat the remaining $n$ qubits as noninteracting spectator
spins, which contribute to the dimension of the overall Hilbert space, but not
to the complexity of the CNOT gate. This helps us to demonstrate the
efficiency of RS expedited algorithms for a varying size of qubit register
with same gate complexity. The internal NMR Hamiltonian for the system in the
multiply-rotating frame is given by

$H_{0}=2\pi J_{12}I_{z}^{1}I_{z}^{2},$ (11)

where $I_{z}^{i}$ represents the $z$-component of the $i$th spin operator. For
the $k$th segment of the $n+2$ channel control sequence, the total Hamiltonian

$H_{k}=H_{0}+\sum_{j=1}^{2+n}\omega_{x,j,k}I_{x}^{j}+\omega_{y,j,k}I_{y}^{j}.$
(12)

Here $\omega_{x,j,k}$ and $\omega_{y,j,k}$ are the $x$ and $y$ components of
radio frequency (RF) field on the $j$th nuclear species in $k$th segment of
the control sequence. In practice however, there exists a spatial RF
inhomogeneity (RFI) of amplitudes over the volume of the sample. One usually
models RFI with an $L$-point distribution and associates a scaling factor
$\\{\xi_{j,l}\\}$ with respective probabilities $\\{p_{j,l}\\}$. Thus, the
Hamiltonian for the $k$th control segment with the $l$th RFI element is

$H_{k}^{(l)}=H_{0}+\sum_{j=1}^{2+n}\xi_{j,l}\omega_{x,j,k}I_{x}^{j}+\xi_{j,l}\omega_{y,j,k}I_{y}^{j}.$
(13)

As explained in the last part of Sec. II.1, the overall cost function
$1-\overline{J}=1-\sum_{l}p_{l}J_{l}$ is obtained by the weighted sum of cost
functions of all the individual elements in the RFI distribution. In the
following, we describe how we can incorporate RS to expedite GRAPE.

###  IV.1 RS expedited GRAPE (RSGRAPE)

Let us first consider the GRAPE algorithm for the NMR context described above.
As explained in Sec. II.1, the $x(y)$ gradients for the $j$th channel in $k$th
segment with $l$th RFI element are expressed by the explicit form of Eq. 7,
i.e.,

$\displaystyle g_{x(y),j,k,l}$
$\displaystyle=2\mathrm{Re}\left(-i\delta\mathrm{Tr}\left[P_{k,l}^{\dagger}I_{x(y)}^{j}X_{k,l}\right]\mathrm{Tr}\left[X_{k,l}^{\dagger}P_{k,l}\right]\right).$
(14)

Here $X_{k,l}$ and $P_{k,l}$ are respectively the forward and backward
propagators for the $k$ segment with $l$th RFI element. In the traditional
GRAPE algorithm, one evaluates all the gradients and then calculates the
update values for the control amplitudes as described in Sec. II.1. This is
the heart of the algorithm and involves evaluating a large number of
propagators via matrix exponentiation. Accordingly, this is the bottleneck for
numerical efficiency, especially for the larger number of qubits. Here comes
the application of machine learning. Instead of evaluating all the gradients,
we only need to evaluate a fraction of the gradients, form a rating matrix,
and then let RS predict the rest of the gradients. In the language of RS, we
treat each segment as a consumer and each RF amplitude (corresponding to
indicies ${x(y),j,l}$) as a product. The corresponding rating matrix is in the
form of the TABLE 1.

RFI |  $l=1$ $\cdots$
---|---
Qubits | j=1 $\cdots$ | j=2+n
segments (k) phases | x | y | x | y
1 | $g_{x,1,1,1}$ | $g_{y,1,1,1}$ | ? | $g_{y,2+n,1,1}$
2 | $g_{x,1,2,1}$ | ? | $g_{x,2+n,2,1}$ | ?
3 | $g_{x,1,3,1}$ | $g_{y,1,3,1}$ | ? | $g_{y,2+n,3,1}$
4 | ? | $g_{y,1,4,1}$ | $g_{x,2+n,4,1}$ | ?
$\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$
Table 1: The rating matrix for RSGRAPE. Here rows correspond to various
segments ($k=1,2,\cdots,N$) and columns correspond to the $x(y)$ RF amplitudes
on all the heteronuclear qubits ($j=1,2,\cdots,2+n$) with various RFI
distribution elements ($l=1,2,\cdots,L$). The elements $g_{x(y),j,k,l}$ are
the gradients, and the goal of RS is to predict the unknown elements
(indicated by ‘?’).

We first study the dependence of RSGRAPE on the sparsity of the rating matrix.
To this end, we generate a two-qubit CNOT gate while varying the sparsity of
the rating matrix from 0% to 90%. In each case, we monitor, at the end of 500
iterations, the final infidelity and the time advantage

$\displaystyle\Gamma(\mbox{RSGRAPE})$
$\displaystyle=\tau(\mbox{RSGRAPE})/\tau(\mbox{GRAPE}),$ (15)

the ratio of computational times of RSGRAPE and GRAPE.

We now demonstrate RSGRAPE for generating an optimal control sequence
implementing a CNOT gate on the first two qubits of the $2+n$ qubit system. In
our analysis we have varied $2+n$ from $2$ to $8$ and in each case generated
an optimal control sequence with $N=200$ segments. For RS prediction we used
latent vectors of dimension 10. RFI is modeled by $L=5$ point distribution
with $\xi_{j,l}\in[0.8,0.9,1.0,1.1,1.2]$ and uniform probability
$p_{j,l}=0.2$. Firstly, we vary the sparsity (s) from $10\%$ to $90\%$ and
compare the gradients for the standard GRAPE ($g_{a}$) vs gradients predicted
by RS ($g_{p}$) as shown by red dots in Fig. 3 (a). For reference, also shown
are the ideal expected curves $g_{a}=g_{p}$ (in blue lines in Fig. 3 (a)).
It’s clear that RS is able to predict the gradients quite efficiently,
especially for sparsity values below $60\%$. For larger sparsity values, the
RS predictions become increasingly inaccurate and therefore worsens the
convergence of the algorithm.

Now, we iterate (for $i=1,\cdots,500$) GRAPE (as explained in Sec. II.1) as
well as RSGRAPE (as explained above) algorithms. In Fig. 3 (b), we plot the
final infidelity $1-\overline{F}$ (blue curves) as well as time advantage
$\Gamma(\mbox{RSGRAPE})$ (red curves) vs sparsity ($s$) for two (solid curves)
and four (dashed curves) qubit systems. Here for reliable analysis, all the
data points are obtained by averaging $10$ independent trials each starting
from a random initial guess. We find that the final infidelity remains low
($<0.005$) till about $60\%$ sparsity, and increases afterwards. Also, we
achieve time advantages by over a factor of 1.5 for both two and four qubits
systems.

We now set the sparsity of the rating matrix to 50%, meaning only half the
number of gradients randomly selected out of the total $2N(2+n)L$ elements
need to be evaluated using Eq. 14. The rest of the gradients in each iteration
are predicted by the RS algorithm. Again, for reliable analysis, we average
the results of 15 independent trials of RSGRAPE algorithm each starting from a
random initial guess. For comparison, we also carryout the standard GRAPE
algorithm in each case and monitor the time advantage
$\Gamma(\mbox{RSGRAPE})$. Fig. 3 (c) shows the infidelity ($1-\overline{F}$)
versus iteration number $i$ for various sizes $2+n$ of the qubit register. It
is clear that the convergence of RSGRAPE is not compromised compared to the
standard GRAPE algorithm, despite only 50% of the gradients being exactly
evaluated. In all the cases, the infidelity was well below 0.01.

Fig. 3 (d) shows the time-advantage for various sizes of the qubit register.
For small qubit registers (up to 4 qubits), the advantage is above 20 %, while
for larger registers (for 8 qubits), we find almost 100% time advantage. This
is because the RS overhead is dependent on (i) the dimension of the rating
matrix, which increases only linearly with the size of qubit register and (ii)
the dimensions of latent vectors. Therefore, as the complexity of GRAPE
algorithm increases exponentially with the size of the qubit register, the RS
overhead becomes increasingly insignificant, and the time-advantage improves.
However, for 50 % sparsity, $\Gamma(\mbox{RSGRAPE})$ remains bounded by a
factor of 2.

Figure 3: (a) Exact $g_{a}$ and predicted $g_{p}$ gradients vs sparsity (s in
%). (b) Infidelity $1-\overline{F}$ (blue) and time advantange
$\Gamma(\mbox{RSGRAPE})$ (red) vs sparsity ($s$ in %) of a CNOT gate generated
using RSGRAPE in two (solid line) and four (dashed line) qubit systems. (c)
Infidelity $1-\overline{F}$ vs number $i$ of iterations varying from 1 to 500
and for size $2+n$ of the qubit register varying from 2 to 8. Here solid-red
and dashed-blue lines correspond to RSGRAPE and GRAPE algorithms. (d)
$\Gamma(\mbox{RSGRAPE})$ vs size $2+n$ of the qubit register. Here (c) and (d)
are obtained with a sparsity value of 50 %.

###  IV.2 RS expedited SAGRAPE (RSSAGRAPE)

We now explain the RS expedited hybrid SAGRAPE algorithm. As explained in Sec.
II.2, an important step in SA in every iteration is to scan the neighborhood
points of the current solution. We add random deviation functions to the
current control sequence (of each heteronuclear qubit) to obtain a set of
neighborhood points. Since one needs to scan a large number of neighbourhood
points, this step forms a bottleneck in the standard SAGRAPE algorithm. This
is where RS can bring about a significant speedup.

The entire control sequence of current iteration can be represented in the
matrix form $\Omega^{(i)}\equiv\\{\omega^{(i)}_{x(y),j,k}\\}$ of dimension
$N\times 2(2+n)$. In our method, we first select a set of $Q$ spline vectors
$\\{s_{q}\\}$, each of dimension $N\times 1$. We also choose a set of $M$
random scaling vectors $\\{c_{m}\\}$, each of size $1\times 2(2+n)$. We now
setup the neighborhood points by adding the deviation function
$w^{(i)}s_{q}\cdot c_{m}$, i.e.,

$\Omega_{q,m}=\Omega^{(i)}+w^{(i)}c_{m}\cdot s_{q},$ (16)

where $w^{(i)}$ is a scalar weight parameter which can be gradually reduced
with iteration number to shrink the neighbourhood region. Now, we determine
the average fidelity $\overline{F}_{q,m}$ and form the rating matrix (see
TABLE. 2). In the RS expedited algorithm, we don’t need to evaluate all the
elements of the rating matrix, but only a set of randomly selected elements.
Rest of the elements are efficiently predicted by RS. The sequence
$\Omega_{q,m}$ corresponding to the maximum fidelity $\overline{F}_{q,m}$ is
then passed to the SA algorithm for comparing with the threshold function as
explained in Sec. II.2.

spline ($s_{q}$) scaling ($c_{m}$) | $c_{1}$ | $c_{2}$ | $c_{3}$ | $c_{4}\cdots$
---|---|---|---|---
$s_{1}$ | $\overline{F}_{1,1}$ | $\overline{F}_{1,2}$ | ? | $\overline{F}_{1,4}$
$s_{2}$ | $\overline{F}_{2,1}$ | ? | $\overline{F}_{2,3}$ | ?
$s_{3}$ | $\overline{F}_{3,1}$ | $\overline{F}_{3,2}$ | ? | $\overline{F}_{3,4}$
$s_{4}$ | ? | $\overline{F}_{4,2}$ | $\overline{F}_{4,3}$ | ?
$\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$
Table 2: The rating matrix for the RSSAGRAPE algorithm. Here rows correspond
to various spline functions $s_{q}$ and columns correspond to the scaling
factors $c_{m}$. Each element of rating matrix is the fidelity
$\overline{F}_{q,m}$ of the neighbourhood point $(q,m)$ obtained from the
current solution by adding the deviation function $w^{(i)}s_{q}\cdot c_{m}$.

Based on the procedure explained above, we now employ RSSAGRAPE to generate a
CNOT gate on a two-qubit system, using five iterations of GRAPE after every
iteration of SA. We used $N=200$ segment control sequence and scanned 100
neighbourhood points (using $Q=10$ spline functions and $M=10$ scaling
vectors) in each SA iteration. We again set the dimension of the latent
vectors to 10. The red dots in Fig. 4(a) correspond to the exact fidelities
($F_{a}$) plotted against the predicted fidelities ($F_{p}$) for a set of
random neighborhood points at various values of sparsity $s$. The blue lines
corresponding to expected distribution $F_{p}=F_{a}$ are also shown for
reference. Evidently, we see a good correlation of the predicted fidelities
with the exact values, thus confirming the effectiveness of RS, especially for
lower sparsity values.

Figure 4: (a) Exact ($F_{a}$) and predicted ($F_{p}$) fidelities of random
neighbourhood points versus sparsity ($s$ in %) in a particular SA iteration.
(b) Infidelity $1-\overline{F}$ vs number of iterations $i$ for a 2 qubit
system with GRAPE (dashed blue line) or RSSAGRAPE (solid red line). In (b) RS
was carried out with sparsity $s=50$ %.

We now set the sparsity of the rating matrix to 50%, meaning only half the
neighborhood points are evaluated exactly, while the remaining ones are
predicted by RS. Fig. 4(b) displays the infidelity $1-\overline{F}$ (averaged
over 10 independent trials each starting with a random initial guess) plotted
versus iteration number $i$ for SAGRAPE (dashed blue line) as well as
RSSAGRAPE (solid red line). It is clear from the average infidelity trajectory
that the convergence is not sacrificed by the partial prediction of the rating
matrix by RS. Furthermore, when compared with the SAGRAPE algorithm, we found
a time advantage

$\displaystyle\Gamma(\mbox{RSSAGRAPE})=\frac{\tau(\mbox{SAGRAPE})}{\tau(\mbox{RSSAGRAPE})}=1.9,$
(17)

meaning the RS enhancement has almost doubled the speed of the SAGRAPE
algorithm.

##  V Summary and outlook

Machine learning techniques are increasingly being utilized in almost every
field of science. Here we use recommender system (RS), a type of machine
learning technique to enhance the efficiency of quantum control algorithm,
particularly a gradient method (GRAPE) and a meta-heuristic method (simulated
annealing (SA)). Because of the analytical form for gradients, GRAPE is a
powerful tool, but it suffers from two problems. Firstly, being a local search
method, GRAPE is sensitive to local minima. Secondly, evaluating gradients for
all the constant-Hamiltonian segments is a computationally intensive task. The
former can be overcome by a hybrid algorithm such as SAGRAPE, which had
previously been demonstrated [28]. Here we address the latter issue by
employing RS to efficiently predict gradients and thereby to remarkably speed
up the GRAPE algorithm. We demonstrated the RSGRAPE algorithm for up to eight
qubits. Going further, we incorporated RS even in SA, for rapid evaluation of
a large set of random neighborhood points. Finally, by expediting both SA and
GRAPE simultaneously, we demonstrated almost doubling the speed of SAGRAPE.

The entirely different approaches of using RS in GRAPE and SA exemplifies the
flexibility and freedom of incorporating RS in quantum control problems. Note
that the particular approaches we have used are not unique. One can think of
different ways of encoding consumers and products to set up a rating matrix
for implementing RS. The generality of RS approach should allow its
application in conjunction with other gradient methods such as BFGS [22],
function-space control [23], etc. RS can also be used to enhance other meta-
heuristic algorithms as well as global search methods such as genetic
algorithm [26]. We believe, the present work encourages further applications
of machine learning protocols in quantum information tasks.

## Acknowledgements

PB acknowledges support from the Prime Minister’s Research Fellowship (PMRF)
of the Government of India. TSM acknowledges funding from
DST/ICPS/QuST/2019/Q67.

## References

  * Rembold  _et al._ [2020] P. Rembold, N. Oshnik, M. M. Müller, S. Montangero, T. Calarco, and E. Neu, AVS Quantum Science 2, 024701 (2020).
  * Dolde  _et al._ [2014] F. Dolde, V. Bergholm, Y. Wang, I. Jakobi, B. Naydenov, S. Pezzagna, J. Meijer, F. Jelezko, P. Neumann, T. Schulte-Herbrüggen, _et al._ , Nature communications 5, 1 (2014).
  * Georgescu  _et al._ [2014] I. M. Georgescu, S. Ashhab, and F. Nori, Reviews of Modern Physics 86, 153 (2014).
  * Cavina  _et al._ [2018] V. Cavina, A. Mari, A. Carlini, and V. Giovannetti, Phys. Rev. A 98, 012139 (2018).
  * Sekatski  _et al._ [2017] P. Sekatski, M. Skotiniotis, J. Kołodyński, and W. Dür, Quantum 1, 27 (2017).
  * Brif  _et al._ [2010] C. Brif, R. Chakrabarti, and H. Rabitz, New Journal of Physics 12, 075008 (2010).
  * d’Alessandro [2021] D. d’Alessandro,  _Introduction to quantum control and dynamics_ (Chapman and hall/CRC, 2021).
  * Glaser  _et al._ [2015] S. J. Glaser, U. Boscain, T. Calarco, C. P. Koch, W. Köckenberger, R. Kosloff, I. Kuprov, B. Luy, S. Schirmer, T. Schulte-Herbrüggen, _et al._ , The European Physical Journal D 69, 1 (2015).
  * Khaneja  _et al._ [2005] N. Khaneja, T. Reiss, C. Kehlet, T. Schulte-Herbrüggen, and S. J. Glaser, Journal of magnetic resonance 172, 296 (2005).
  * Machnes  _et al._ [2018] S. Machnes, E. Assémat, D. Tannor, and F. K. Wilhelm, Physical review letters 120, 150401 (2018).
  * Palao and Kosloff [2003] J. P. Palao and R. Kosloff, Phys. Rev. A 68, 062308 (2003).
  * Reich  _et al._ [2012] D. M. Reich, M. Ndong, and C. P. Koch, The Journal of chemical physics 136, 104103 (2012).
  * Doria  _et al._ [2011] P. Doria, T. Calarco, and S. Montangero, Physical review letters 106, 190501 (2011).
  * Egger and Wilhelm [2014] D. J. Egger and F. K. Wilhelm, Physical review letters 112, 240503 (2014).
  * Fortunato  _et al._ [2002] E. M. Fortunato, M. A. Pravia, N. Boulant, G. Teklemariam, T. F. Havel, and D. G. Cory, The Journal of chemical physics 116, 7599 (2002).
  * Zahedinejad  _et al._ [2014] E. Zahedinejad, S. Schirmer, and B. C. Sanders, Physical Review A 90, 032310 (2014).
  * An and Zhou [2019] Z. An and D. Zhou, EPL (Europhysics Letters) 126, 60002 (2019).
  * Zhang  _et al._ [2019] X.-M. Zhang, Z. Wei, R. Asad, X.-C. Yang, and X. Wang, npj Quantum Information 5, 1 (2019).
  * Niu  _et al._ [2019] M. Y. Niu, S. Boixo, V. N. Smelyanskiy, and H. Neven, npj Quantum Information 5, 1 (2019).
  * Bukov  _et al._ [2018] M. Bukov, A. G. R. Day, D. Sels, P. Weinberg, A. Polkovnikov, and P. Mehta, Phys. Rev. X 8, 031086 (2018).
  * Khait  _et al._ [2021] I. Khait, J. Carrasquilla, and D. Segal, arXiv preprint arXiv:2108.12441 (2021).
  * De Fouquieres  _et al._ [2011] P. De Fouquieres, S. Schirmer, S. Glaser, and I. Kuprov, Journal of Magnetic Resonance 212, 412 (2011).
  * Lucarelli [2018] D. Lucarelli, Physical Review A 97, 062346 (2018).
  * Batra  _et al._ [2020] P. Batra, V. R. Krithika, and T. S. Mahesh, Phys. Rev. Research 2, 013314 (2020).
  * Krause  _et al._ [1998] J. L. Krause, D. H. Reitze, G. D. Sanders, A. V. Kuznetsov, and C. J. Stanton, Physical Review B 57, 9024 (1998).
  * Bhole  _et al._ [2016] G. Bhole, V. Anjusha, and T. Mahesh, Physical Review A 93, 042339 (2016).
  * Zahedinejad  _et al._ [2015] E. Zahedinejad, J. Ghosh, and B. C. Sanders, Phys. Rev. Lett. 114, 200502 (2015).
  * Ram  _et al._ [2021] M. H. Ram, V. Krithika, P. Batra, and T. Mahesh, arXiv preprint arXiv:2112.01279 (2021).
  * Lü  _et al._ [2012] L. Lü, M. Medo, C. H. Yeung, Y.-C. Zhang, Z.-K. Zhang, and T. Zhou, Physics reports 519, 1 (2012).
  * Bobadilla  _et al._ [2013] J. Bobadilla, F. Ortega, A. Hernando, and A. Gutiérrez, Knowledge-based systems 46, 109 (2013).
  * Aarts and van Laarhoven [1987] E. H. Aarts and P. J. van Laarhoven, in  _Pattern recognition theory and applications_ (Springer, 1987) pp. 179–192.
  * Bertsimas and Tsitsiklis [1993] D. Bertsimas and J. Tsitsiklis, Statistical science 8, 10 (1993).
  * Herlocker  _et al._ [2004] J. L. Herlocker, J. A. Konstan, L. G. Terveen, and J. T. Riedl, ACM Transactions on Information Systems (TOIS) 22, 5 (2004).
  * Su and Khoshgoftaar [2009] X. Su and T. M. Khoshgoftaar, Advances in artificial intelligence 2009 (2009).
  * Singh and Gordon [2008] A. P. Singh and G. J. Gordon, in  _Joint European Conference on Machine Learning and Knowledge Discovery in Databases_ (Springer, 2008) pp. 358–373.
  * Koren  _et al._ [2009] Y. Koren, R. Bell, and C. Volinsky, Computer 42, 30 (2009).
  * Batra  _et al._ [2021] P. Batra, A. Singh, and T. S. Mahesh, Quantum 5, 598 (2021).
  * Polak, E. and Ribiere, G. [1969] Polak, E. and Ribiere, G., R.I.R.O. 3, 35 (1969).
  * Press  _et al._ [1992] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,  _Numerical Recipes in FORTRAN 77 Macintosh Diskette Version 2.0: The Art of Scientific Computing_ , Vol. 1 (Cambridge university press, 1992).
  * [40] Andrew Ng, Stanford Machine Learning Course, .


# An Unified Intelligence-Communication Model for Multi-Agent System Part-I:
Overview

Bo Zhang, Bin Chen, Jinyu Yang, Wenjing Yang, Jiankang Zhang B. Zhang is with
the Artificial Intelligence Research Center, National Institute of Defense
Technology Innovation, Beijing, P. R. China. E-mail:
<EMAIL_ADDRESS>B. Chen is with the Academy of Military Sciences,
Beijing, P. R. China. J. Yang and W. Yang are with the State Key Laboratory of
High Performance Computing, National University of Defense Technology,
Changsha, P. R. China. J. Zhang is with the Southampton Wireless Group,
Southampton University, SO17 1BJ, the United Kingdom.

###### Abstract

Motivated by Shannon’s model and recent rehabilitation of self-supervised
artificial intelligence having a ”World Model”, this paper propose an unified
intelligence-communication (UIC) model for describing a single agent and any
multi-agent system.

Firstly, the environment is modelled as the generic communication channel
between agents. Secondly, the UIC model adopts a learning-agent model for
unifying several well-adopted agent architecture, e.g. rule-based agent model
in complex adaptive systems, layered model for describing human-level
intelligence, world-model based agent model. The model may also provide an
unified approach to investigate a multi-agent system (MAS) having multiple
action-perception modalities, e.g. explicitly information transfer and
implicit information transfer.

This treatise would be divided into three parts, and this first part provides
an overview of the UIC model without introducing cumbersome mathematical
analysis and optimizations. In the second part of this treatise, case studies
with quantitative analysis driven by the UIC model would be provided,
exemplifying the adoption of the UIC model in multi-agent system.
Specifically, two representative cases would be studied, namely the analysis
of a natural multi-agent system, as well as the co-design of communication,
perception and action in an artificial multi-agent system. In the third part
of this treatise, the paper provides further insights and future research
directions motivated by the UIC model, such as unification of single
intelligence and collective intelligence, a possible explanation of
intelligence emergence and a dual model for agent-environment intelligence
hypothesis.

Notes: This paper is a Previewed Version, the extended full-version would be
released after being accepted.

###### Index Terms:

Information Theory, Multi-Agent System, Learning (artificial intelligence),
Communication Channels.

##  I A BRIEF INTRODUCTION OF THE UIC MODEL - FROM SHANNON’S PERSPECTIVE

This Section introduces the motivation of the UIC model and the basic
graphical description of the UIC model, illustrating the relationship between
an intelligent agent and the environment. This section would cover the
following aspects.

  * •

The environment is model as the generic communication channel between agents.

  * •

The agent is model as an communication transceiver. Specifically, the internal
world model is model as a source coder, which encodes the environment, while
the multi-modal perceiving and acting capability are abstracted as multiple
transmitter-receiver pairs.

Before exposing the proposed UIC model, let us first review the Shannon Model.

##  II A BRIEF REVIEW OF THE SHANNON MODEL

This section applies the Shannon Model to describe complex forms of
interactions between agents and environments.

The Shannon model proposed in the pioneering paper “A mathematical theory of
communication” may be illustrated as follows, and please note that the
Weiner’s feedback is not included for now.

Figure 1: Shannon Model without Feedback [1]

##  III COMPLEX INTERACTIONS MODELING: IS SHANNON MODEL SUFFICIENT?

This paper is definitely NOT the first effort to apply Shannon Model to model
and understand the interactions between agents and the environment.

An intuitive question arose immediately (especially the critics from the
communication society and most of the control society) and therefore may
prohibit the further efforts. Is the Shannon channel model sufficient to
describe complex interactions between agents and the environment? As an agent
not only transmits and receives waveforms, but also changes the physical
environment, e.g. moving objects, absorbs solar energy, etc.

Before answering the above question directly, the author may redirect to a
more fundamental question: What is the relationship between Information,
Substance and Energy? The answer to this question may help determine the range
of applications of Shannon Channel Model. Hence, the UIC model may be firmly
built based on the following statements:

  * •

Information is carried by waves. The waves may be electromagnetic waves,
mechanical waves, matter waves and other forms of waves, which are generated
by the four fundamental interactions.

  * •

The environment may be modeled as a possibly infinite set of waves. Each wave
carries a certain amount of information and energy. The different forms of
waves may interact, exchanging energy and information.

  * •

Any form of interactions may be contaminated by noise or uncertainty.
Therefore, the maximum amount of information that may be transferred by a wave
is quantified by the Shannon Channel Capacity.

Based on the above statements, we may come to the following statements: The
environment may be model as a set of multi-modal channels for different forms
of waves.

##  IV THE AGENT MODEL IN THE UIC MODEL

After modelling the environments, and based on the widely adopted philosophy
that an intelligent agent is an open system in the environment, the UIC model
incorporates the most important assumptions of an intelligent agent.

The UIC Assumption: As with the Environment, an agent is also a set of waves,
and the interactions between an agent and the environment may be model as
interactions between waves. The interaction process exchanges energy and
information, meanwhile incorporates uncertainties.

Therefore, the interactions between an agent and the environment may be
categorized by the widely used terminologies in describing the behaviour of an
intelligent agent:

###  IV-A Agent Perception

Perception may be considered as the measurement of energy and information
carried by a wave in the environment.

  * •

Based on the UIC Assumption, perception is modeled as the energy and
information transferring from the environment to the agent.

  * •

The perception module of an agent is then model as a set of channel receivers,
where each receiver extracts information from the waves, and an optimal
receiver extracts the maximum amount of information, as bounded by the Shannon
Capacity.

###  IV-B Agent Action

Action taken by an agent may be in the macro and micro level, and in various
domains, e.g. spatial, frequency, energy, etc.

  * •

Based on the UIC Assumption, any action may be model as a form of wave
interactions, action is model as the energy and information transferring from
the agent to the environment.

  * •

The action module of an agent is then model as a set of channel transmitters,
where each transmitter modulates a certain set of waves with information, and
a feedback from the receiver, as shown by Shannon and Wiener, may greatly
improve the optimal transmitter-receiver design.

We may now graphically depict the first UIC model figure, focusing on using
the Shannon model to interpret the interactions between the agent and the
environment.

Figure 2: Interactive Model between the Agent and the Environment

In the up-coming sections, we will dive deeper into the information source and
destination in the agent.

##  V MODELLING THE INFORMATION SOURCE AND DESTINATION IN THE AGENT

This Section introduces the model for information source and destination in
the UIC model as proposed in the last section. This section would cover the
following aspects.

  * •

The information source and destination is modelled, following the generic
learning agent model, which includes a performance element, a learning
element, a critic and a problem generator.

  * •

This Section reviews the different agent models, with an emphasis on the
design of learning agent model, and propose an unified, complete and operable
design for the learning agent that may lead to strong AI.

  * •

Along with the perception and action model proposed in the last section, we
may form the complete UIC model for an intelligent agent.

If we stick to the Shannon’s model of communication, the information source
and destination may be disconnected in a single agent. However, in order to
allow an agent to receive information from the outside world, to process
information and to change the world by transmitting information, a connection
between the information source and destination in the agent is inevitable.
Therefore, a question arose: How does an intelligent agent process the
information within itself? This paper is definitely NOT the first to answer
the question, so let us review the agent internal models.

##  VI A BRIEF REVIEW OF THE AGENT INTERNAL MODELS

An overview of the major five categories of agent models is given in Section
2.4 of [2]. Except for the simplest form of agent model, the rest of the agent
models includes an internal world model.

###  VI-A Agent Models without Learning

The agent without learning may be modelled as follows.

  * •

Simple Reflex Agent model (SRA model): The agent function is based on the
condition-action rule: if condition then action. So it may work well in a
fully-observable world.

  * •

Model-based, Reflex Agent model (MRA model): A world model is introduced to
maintain some kind of structure, which describes the part of the world which
cannot be observed.

  * •

Model-based, Goal-Based Agent model (MGA model): The goal information is
introduced to describe desirable situations and allow an agent to choose among
multiple possibilities, selecting the one which reaches a goal state. In the
complex adaptive system, Holland proposed to incorporates co-existence of
multiple micro SRA agents in an agent, and allowing the agent to select action
sequences[3].

  * •

Model-based, Utility-Based Agent model (MUA model): The utility function is
introduced to measure the space spanning from the goal state and the non-goal
state.

###  VI-B Learning Agent Models

Afore-mentioned agent models are static, as they do not explain how an agent
may learn and evolve. Therefore, the LA model is quite different by
introducing four conceptual components:

  * •

A performance element to select external actions, which is equivalent to a
static agent model that percepts from the world and decides external actions.

  * •

A learning element to make improvements by changing the performance element
and gathering knowledge from the performance element.

  * •

A critic to enable feedback on how the agent is doing and how the performance
element may be modified to do better.

  * •

A problem generator to select actions leading to new and informative
experiences.

Figure 3: Learning Agent model (LA model) [2]

##  VII RETHINK THE LEARNING AGENT MODEL WITH AN INTERNAL WORLD MODEL

Although it is not pointed out explicitly in [2], the learning agent model may
include a world model in its performance element. Therefore, the learning
element may change the world model.

If the design principle of the learning agent model is accepted (why NOT if we
have not come up with a more reasonable model), then the problem of designing
an intelligent agent may be further decomposed as follows:

  * •

Performance Element Design: What is inside the performance element and what is
a WM looks like?

  * •

Learning Element and Critic Design: How to gather knowledge and improve
performance based on the feedback provided by the critic in a complex dynamic
world? Also, how to update the WM and rest of Performance Element?

  * •

Problem Generator Design: How to find a GOOD tradeoff between short-term and
long-term reward in a changing world?

The answers to the above four questions may lead to the design principles of
Strong Artificial Intelligence. Therefore, the answers to these questions are
so attractive and many pioneering thoughts (NOT solutions) have been proposed
to light our way.

###  VII-A Agent Model in Complex Adaptive Systems

Professor Holland H. John is the pioneer of Complex Adaptive System Theory. In
his 1995 treatize [3], Holland proposed the following works:

  * •

Performance Element Design: Holland proposed the design of the performance
system, which has an input from sensors, an output to the actuators, while the
performance element is a collection of interactive micro-agents based on IF-
THEN rules.

  * •

Learning Element and Critic Design: A credit assignment mechanism was proposed
to allow multiple micro-agents to compete with each other. The micro-agent
with higher reward from the critic is given a larger credit, hence improving
performance of the agent.

  * •

Problem Generator Design: A rule discovery mechanism was proposed based on the
genetic algorithm, which was also invented by Holland.

###  VII-B Agent Model with Commonsense Thinking

In [5], Professor Marvin Minsky proposed to model a (human) brain into a
collection of resources, which may be activated by different emotion states.
The thinking activities may be broken down into a series of six-layered
actions: instinctive reactions, learned reactions, deliberate thinking,
reflective thinking, self-reflective thinking and self-conscious reflection.
Higher-level thinking is built upon the ones below, and the lowest instinctive
reaction level is driven by If-Do rules.

Specifically, we may re-organize Minsky’s thoughts as follows:

  * •

Performance Element Design: The instinctive reactions may be modelled as IF-DO
rules, which are equivalent to Holland’s IF-THEN rules.

  * •

Learning Element and Critic Design: The higher-layer actions are built upon
instinctive reactions. Minsky emphasized the limitation of IF-THEN rules
because the scalability of directly using IF-THEN rules for describing the
complex world is infeasible, therefore, representations, evaluations and
manipulations of higher-level abstraction is required. It should be noted that
Holland also dived much deeper later on following the same rule of agent
model, introducing the concept of dynamic finitely generated systems for
modelling complex interactions between agents [3].

  * •

Problem Generator Design: The layers higher than deliberate thinking motivate
an agent to criticize the reactive actions and explore un-explored world.

###  VII-C Self-Supervised Learning Agent Model

Very recently, Professor Yann Lecun proposed an agent model based on the
concept of Self-Supervised Learning with both a deep structure and a world
Model[5]. The agent may do reasoning by combining predicting and planning with
the aid of model-based deep reinforcement learning. Specifically, Lecun
considered the following strategies:

  * •

Performance Element Design: A world model is needed, accumulating background
knowledge about how the world works and perhaps include common sense. Deep
(artificial neural) networks may model complex interactions with the aid of
non-linear transformations within a manageable scale.

  * •

Learning Element and Critic Design: Latent-variable forward models is proposed
for planning and learning policies, but it is also far from a complete
solution.

  * •

Problem Generator Design: It seems Lecun did not consider it for the moment.

Figure 4: LeCun’s Agent Model with an Internal World Model[6]

##  VIII UNIFYING THE LEARNING AGENT MODELS

The models proposed by Holland[2][3], Minsky[4] and Lecun[5] may approach the
design of a learning agent model for strong AI from different perspective and
different emphases. With some non-difficult manipulation of their models, the
performance element, learning element and critic design for an unified
intelligent agent model may be proposed as follows:

  * •

Performance Element Design: An aggregation of micro-agents described by simple
reactive rules.

  * •

Learning Element and Critic Design: A collection of action modules for
building, manipulating and evaluating higher-level (likely to be non-linear)
representation of micro-agents’ behaviours.

However, it is not apparent to come up with a concrete and operable design for
the Problem Generator Design yet. Here, we try to solve the problem based on
the following basic understandings:

  * •

Non-Complete Representation as Motives: Even though an agent have non-linear
basis/generators for representing very complex interactions in the world, we
may still face up to non-complete representations, with limited resources.
Therefore, the agent may adopt a sub-optimal (or an action that cannot be
proved optimal with limited resource) in the world.

  * •

Multiple-Trial in the World Model: Though the agent cannot find the optimal
actions, the agent may come up with a set of sub-optimal candidate actions,
where the size of the set is limited based on available resources. After
taking previous actions in the world, these candidate actions could be tested
within the updated world model.

Based on the above insights, we may propose an operable pathway to the problem
generators for strong AI.

  * •

Problem Generator Design: A set of modules to evaluate non-complete
representation of complex action-reward interactions, to generate candidate
action set, and to trigger reflective actions in the world model..

We may now graphically represent our UIC model with the receivers for
modelling sensors’ perceptions, transmitters for modelling actuators’ actions,
as well as the performance element, learning element, critic and problem
generator for information source/destinations. Please note we not only merge
the different components in the learning agent model, but also provide
guidelines for designing each of them.

The details of the UIC design would be provided in the following sections,
through more elaborations, case studies and insightful discussions.

Figure 5: The Proposed UIC model with an Internal Model for Learning Agent

##  IX CONLUSIONS

In this paper, we first review the Shannon model for communications, and
propose an unified intelligence-communication model to describe the
interactions between intelligent agents and the environment, where the
perception and action of an agent may be modelled with information reception
and transmission.

Then, we review and discuss the relationship between different models for
intelligent agents with an emphasis on the learning agent model, and propose
an unified, complete and operable design for intelligent agents with learning
capabilities. Along with the perception and action model proposed, we may form
the complete UIC model for an intelligent agent.

In the second part of this treatize, we would show that case studies would be
provided, exemplifying the adoption of the UIC model. Specifically, we would
dive deep into two representative cases, the analysis of a natural bee-colony,
as well as the co-design of communication, perception and action in artificial
multi-robot system.

## References

  * [1] Shannon, Claude Elwood. “A mathematical theory of communication.” Bell system technical journal 27.3 (1948): 379-423.
  * [2] Russell, Stuart J. and Peter Norvig. Artificial intelligence: a modern approach. Pearson Education Limited, 2011.
  * [3] Holland, John H . Hidden Order: How Adaptation Builds Complexity. Leonardo, 1995.
  * [4] Holland, John H . Signals and boundaries : building blocks for complex adaptive systems. MIT Press, 2014.
  * [5] Minsky M . The Emotion Machine: Commonsense Thinking, Artificial Intelligence, and the Future of the Human Mind[M]. SIMON & SCHUSTER, 2007.
  * [6] Yann LeCun. Learning World Models: the Next Step towards AI. IJCAI KeyNote Speech, 2018