\section{Introduction}
Experimental quantum computing has made great strides in the past few decades. From the basic quantum logic gates realized by Monroe et al. \cite{monroe:1995} to the 32 quantum volume machines recently announced by IBM \cite{IBM:QV32}, quantum computing systems have grown in size, computational power and general interest. Publications on new or improved algorithms point to the potential of quantum computing. Yet, several key obstacles must be overcome if today’s envisioned applications become tomorrow’s reality. One such obstacle is the infidelity of two-qubit gates.

Two-qubit gate operations are fundamental to quantum computing. It has been shown that a quantum universal gate set can be made with one and two-qubit gates \cite{divincenzo:1995} \cite{barneco:1995} and as a result two-qubit gates are ubiquitous in quantum algorithms.
However, two-qubit gate fidelity must be sufficiently high to ensure quantum algorithms yield useful results.

Characterizing two-qubit gate errors is therefore paramount to obtaining high quality computational results
from NISQ era quantum processors. Gate error data is used to evaluate the general performance and hence selection of a quantum computing system. Furthermore, it is used in compiler optimization techniques such as noise-aware circuit compilation \cite{murali:2019} and variability-aware qubit allocation and movement \cite{tannu:variability:2019}. Inaccuracies in the characterization of gate operation can yield unexpected results and poor computational performance.

At IBM the Controlled-NOT (CNOT) gate is characterized by a single error rate. However, an exploration of the pulse-level quantum controls used by IBM to execute a CNOT gate uncovered two different microwave pulse schedules for
a CNOT realized with a given pair of qubits, depending on which qubit is the control.
If the error characteristics of the CNOT configurations differ, then the systems, controls and compilers that make use of this data may not be fully optimized.

In this paper, we explore the performance of the different CNOT configurations and attempt to quantify an observed asymmetry. Our results suggest that taking this asymmetry into account when mapping a circuit onto physical gates could reduce the error in the computation.

\section{Background}
\begin{comment}
Our findings and discussion require knowledge of the physical implementation of the CNOT gate on the IBM quantum processors. We therefore spend a brief amount of time reviewing the cross-resonance gate, the CNOT gate and the IBM quantum system.

\subsection{The Cross-Resonance Gate}

The cross-resonance (CR) gate is a microwave-only gate that stimulates a ZX interaction by driving the control qubit with microwave pulses at the frequency of the target qubit \cite{alexander:2020}. As described in \cite{alexander:2020} and \cite{magesan:2020}, the complete CR Hamiltonian has the structure
\begin{equation}\label{cr_hamiltonian}
\begin{split}
   & \mathcal{H} = \frac{Z \bigotimes A}{2} + \frac{I \bigotimes B}{2} \\
    & A = \omega_{ZI}I + \omega_{ZX}X + \omega_{ZY}Y + \omega_{ZZ}Z \\
    & B = \omega_{IX}X + \omega_{IY}Y + \omega_{IZ}Z
\end{split}
\end{equation}
Using echo and cancellation pulses the CR gate can be tuned to a ZX gate described by the unitary%
\begin{equation}\label{cr_unitary}
U_{ZX}(\theta_{ZX}) = e^{\left(-i\theta_{ZX}\frac{ZX}{2}\right)}
\end{equation}
When \(\theta_{ZX}\) is set to \(\frac{\pi}{2}\) it can be used as a generator for the CNOT gate when combined with additional single qubit gates \cite{sheldon:2016}.

\subsection{The Controlled-NOT Gate}

The controlled-NOT (CNOT) is a two-qubit logic gate with one qubit acting as the control and the other as the target. When the control qubit is in the \(\ket{1}\) state the target qubit is flipped. When the control qubit is in the \(\ket{0}\) state the target qubit is left alone. If \(q_{a}\) and \(q_{b}\) are representative of two qubits in a quantum processor, the action of the CNOT in the standard basis can be written as
\begin{equation}\label{cnot_ab_matrix}
\text{CNOT}_{a,b} =
\begin{bmatrix}
0 & 1 & 0 & 0 \\
0 & 0 & 0 & 1 \\
0 & 0 & 1 & 0
\end{bmatrix},
\end{equation}
\begin{equation}\label{cnot_ba_matrix}
\text{CNOT}_{b,a} =
\begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 \\
0 & 0 & 1 & 0 \\
0 & 1 & 0 & 0 \\
\end{bmatrix},
\end{equation}
where \(q_{a}\) is the control in the first CNOT and \(q_{b}\) is the control in the second CNOT.

\subsection{Two-Qubit Gates on IBM Q Processors}
\end{comment}
IBM uses superconducting transmon quantum processors with fixed frequency qubits, and uses the cross-resonance (CR) gate as a two-qubit entangling gate \cite{chow:2011}.
The CR gate stimulates a $ZX$ interaction by driving the control qubit with microwave pulses at the frequency of the target qubit \cite{alexander:2020}.
To generate a CNOT gate, the CR gate is combined with an additional frame change $S$ on the control qubit and a \(X_{\frac{\pi}{2}}\)on the target qubit \cite{IBM:faq}. IBM combines the CNOT gate with the single qubit $U$ gate, defined by \(U(\theta, \psi, \lambda) = R_z(\phi)R_y(\theta)R_z(\lambda)\), to form a universal quantum gate set \cite{zulehner:2019}.
\begin{figure}[b]
\begin{subfigure}[t]{0.24\textwidth}
\centering
\includegraphics[width=0.7\textwidth]{nonoptimal-cnot}
\caption{}
\label{fig:nonoptimal cnot}
\end{subfigure}
\begin{subfigure}[t]{0.24\textwidth}
\centering
\includegraphics[width=0.7\textwidth]{optimal-cnot}
\caption{}
\label{fig:optimal cnot}
\end{subfigure}
\begin{subfigure}[t]{0.48\textwidth}
\centering
\includegraphics[width=0.4\textwidth]{inverted-cnot}
\caption{}
\label{fig:inverted cnot}
\end{subfigure}
\caption{Three examples of a CNOT gate. The circuits in (a) and (c) are equivalent. If a circuit is programmed to execute the $\mathrm{CNOT_{1,0}}$ seen in
part (a)
but it is not physically allowed, a compiler may swap it out for the circuit in part (c).
This effectively takes the $\mathrm{CNOT_{0,1}}$ to the $\mathrm{CNOT_{1,0}}$; the gate originally intended by the circuit designer.}
\label{fig:cnot replacements}
\end{figure}

In general, the CR gate used in the IBM quantum processors is restricted to use the higher frequency qubit as the control for the lower frequency qubit, since this results in a shorter entangling time \cite{tripathi:2019} \cite{magesan:2020} \cite{IBM:backend}. However, this restriction is not enforced at the software level and circuit designers are free to configure the CNOT gate to use the lower frequency qubit as the control. To overcome this conflict, circuit or microwave pulse compilers replace the disallowed CNOT with an equivalent gate or pulse sequence that adheres to the underlying physical restriction. An example of one such gate sequence can be seen in FIG. \ref{fig:inverted cnot}. This means that the gate sequence for a CNOT may vary with the configuration of the control and target qubits. If these varied gate sequences have different error characteristics, it is plausible that CNOT performance may vary with configuration as well.

We endeavored to explore this idea within the IBM quantum system. However, before we could begin, we first needed to verify whether the microwave pulse sequences used in the IBM control systems varied with CNOT configuration. Using Qiskit \cite{Qiskit}, an open source software development kit used for working with the IBM quantum processors, we built two circuits. The first contained a single CNOT gate configured to use the higher frequency qubit as the control. The second circuit contained a single CNOT gate configured to use the lower frequency qubit as the control. We then used the \texttt{schedule()} \cite{IBM:schedule} function available in Qiskit to generate a pulse schedule for each circuit (using the \texttt{ibmq\char`_almaden} quantum processor as the backend). As expected, we found the schedules to differ, as detailed in TABLE \ref{table:cnot pulse comparison}. With this established, we then moved on to designing our experiments.

\begin{table}[t]
\centering
    \begin{tabular}{|m{2cm}|m{2cm}|m{2cm}|}
        \hline
        \textbf{CNOT Configuration} & \textbf{Number of Gates} & \textbf{Schedule Duration}  \\
        \hline
        $\mathrm{CNOT_{hfreq}}$ & 13 & 348 ns \\
        \hline
        $\mathrm{CNOT_{lfreq}}$ & 21 & 384 ns \\
        \hline
    \end{tabular}
    \caption{The number of gates and the duration of the pulse schedule was dependent on CNOT configuration. When the CNOT was configured to use the lower frequency qubit as the control ($\mathrm{CNOT_{lfreq}}$) the schedule included more gates and had a longer duration.}
    \label{table:cnot pulse comparison}
\end{table}

\section{Methods}

To evaluate the symmetry of a CNOT, we conducted an experiment which consisted of multiple circuits, with each circuit consisting of $n$ repetitions of an identity operation, where $n= 1,2, …, 6$. The identity operation, which can be seen in FIG. \ref{fig:identity operation}, is composed of a Bell State and its inverse. A barrier is inserted to ensure the circuit is not converted to an identity gate during compilation.
This experiment was run twice for each pair of qubits included in our evaluation. In the first run
the CNOT gates were configured to use $q_0$ as the control and $q_1$ as the target. In the second run
the CNOT gates were configured to use $q_1$ as the control and $q_0$ as the target.
Note, here $q_0$ and $q_1$ are being used generically to represent two different qubits and do not explicitly represent the qubits indexed with 0 or 1 on any particular quantum processor.

In this paper we refer to the circuit with $n$ repetitions of the identity operation as the $n$-stage circuit. The 3-stage circuits can be seen in FIG. \ref{fig:3-stage circuit}. Note that barriers are inserted between each stage to ensure the compiler does not combine identity operations.

\begin{figure}[b]
    \centering
    \includegraphics[width=0.40\textwidth]{identity-operation}
    \caption{The identity operation is composed of a Bell state and its inverse. The $n$-stage circuit consists of $n$ repetitions of this identity operation. }
    \label{fig:identity operation}
\end{figure}

\begin{figure*}[t]
    \centering
    \begin{subfigure}{0.9\textwidth}
        \centering
        \includegraphics[width=0.9\textwidth]{3-stage-circuit}
        \caption{}
        \label{fig:3-stage q0_1}
    \end{subfigure}
    \begin{subfigure}{0.9\textwidth}
        \centering
        \includegraphics[width=0.9\textwidth]{3-stage-circuit-reverse}
        \caption{}
        \label{fig:3-stage q1_0}
    \end{subfigure}
    \caption{In the 3-stage circuit, the identity operation is repeated 3 times. Barriers are placed between the identity operations and between the Bell state and its inverse to prevent compiler optimizations that may alter the gate sequence. The direction of the CNOT and placement of the Hadamard gates are reversed in \ref{fig:3-stage q0_1} and \ref{fig:3-stage q1_0} to evaluate the effect of CNOT direction on circuit performance.}
    \label{fig:3-stage circuit}
\end{figure*}

For each circuit in the experiment, the selected qubits were first prepared in the \(\ket{00}\) state after which the $n$ identity operations of the $n$-stage circuit were run. In the ideal scenario each identity operation would return the qubits to the \(\ket{00}\) state, and therefore we would also expect to measure the qubits in the \(\ket{00}\) state at the end of the $n$-stage circuit. However, given the presence of noise and systematic errors in the physical system, actual results varied and included all states in the computational space. We benchmarked the performance of each CNOT direction on a qubit pair by calculating the percentage of shots that were measured in the \(\ket{00}\) state. A CNOT that exhibited symmetry in implementation direction should produce ground state percentages invariant of implementation direction. Even in the presence of noise and systematic error we expect the percentages to be statistically similar if the CNOT is symmetric.

We define the asymmetry of an $n$-stage circuit as the absolute magnitude of the difference between the fraction of correct answers for each implementation:

\begin{equation}\label{asymmetry}
    f(n) = \left|g_{01}(n) - g_{10}(n)\right|,
\end{equation}
\begin{equation}\label{percent in ground}
    g_{c,t}(n) = \frac{G_{c,t}(n)}{T_{c,t}(n)},
\end{equation}
where \textit{\(G_{c,t}(n)\)} and \(T_{c,t}(n)\) are the number of shots measured in the ground state and the total number of shots, respectively, for the $n$-stage circuit with CNOTs configured with \textit{c} as the control and \textit{t} as the target

We classified the CNOT pair as asymmetric if $f(n)\geq 0.02$ for any value of $n$ considered (i.e., for $n$ between 1 and 6). Since measurement errors are typically of order 2-3\%, this is roughly equivalent to classifying the CNOT pair as asymmetric if the asymmetry is larger than the read error.

Specific qubit pairs on each processor were evaluated using the experiment described above. The two qubits used in each experiment were nearest neighbors and were connected by a bus resonator to ensure SWAP gates were not included during circuit compilation. Each $n$-stage circuit was executed 3 times with a shot count of 4,096 shots. Results were aggregated so the total shot count for an $n$-stage circuit was 12,288.

Experiments were conducted on \texttt{ibmq\char`_valencia}, \texttt{ibmq\char`_almaden} and \texttt{ibmq\char`_paris}, which are 5, 20 and 27 qubit processors, respectively. Our work focused primarily on \texttt{ibmq\char`_almaden} on which we tested every qubit pair in the coupling map. To determine if the asymmetry was present on other quantum processors we ran additional experiments on \texttt{ibmq\char`_valencia} and \texttt{ibmq\char`_paris}. We were able to confirm the presence of CNOT asymmetry on these machines through a limited number of experiments on a subset of the available qubit pairs.
The qubit pairs evaluated on each machine can be found in Appendix \ref{appendix:experiments}.

\section{Results}

\subsection{\texttt{ibmq\char`_almaden}}
 An asymmetry was observed in 71.8\% of the experiments conducted on \texttt{ibmq\char`_almaden}. Of the 23 qubit couplings available on \texttt{ibmq\char`_almaden}, 21 exhibited an asymmetry at some point in our evaluation. In 9.6\% of tests, the asymmetry was seen in the shallow depth 1-stage circuit.  The 6-stage circuit, the deepest circuit tested, exhibited the asymmetry in 67.12\% of tests. The divergence in performance was most pronounced for the 6-stage circuits which on average showed a 5.94\% difference in the fraction of correct answers for the different CNOT implementations.

 Table \ref{tab:almaden results} displays the percentage of experiments in which an asymmetry was observed for each $n$-stage circuit.

\begin{table*}[t]
    \centering
    \begin{tabular}{|c|c|}
    \hline
        \textbf{Circuit} & \textbf{Percentage of Experiments Exhibiting an Asymmetry} \\
         \hline
        1-stage circuit & 9.6\% \\
        \hline
        2-stage circuit & 30.14\% \\
        \hline
        3-stage circuit & 38.36\% \\
        \hline
        4-stage circuit & 52.05\% \\
        \hline
        5-stage circuit & 63.01\% \\
        \hline
        6-stage circuit & 67.12\% \\
        \hline
    \end{tabular}
    \caption{The percentage of experiments that exhibited an asymmetry for each of the $n$-stage circuits. The asymmetry increased with circuit depth.}
    \label{tab:almaden results}
\end{table*}

In Fig. \ref{fig:cnot asymmetry} we plot the fraction of correct answers for each CNOT implementation as a function of circuit stage. Fig. \ref{fig:almaden 5 and 6} illustrates the asymmetric performance of a CNOT on qubits 5 and 6 on \texttt{ibmq\char`_almaden}. The fraction of correct answers diverges substantially as $n$ increases. The asymmetric behavior can even be seen at the shallow-depth 2-stage circuit. This suggests the possibility of non-optimal performance for circuits which incorporated the \(\text{CNOT}_{5,6}\) gate on \texttt{ibmq\char`_almaden} on the day this experiment was conducted.

In contrast, Fig. \ref{fig:almaden 1 and 2} shows the performance of a CNOT gate configured to use qubits 1 and 2 on \texttt{ibmq\char`_almaden} that did not display an asymmetry. In this plot we do not see a divergence in the fraction of correct answers and performance of the \(\text{CNOT}_{1,2}\) on \texttt{ibmq\char`_almaden} on the day of experimentation appeared to be invariant to CNOT configuration.

\begin{figure}[h]
    \centering
    \begin{subfigure}{0.49\textwidth}
        \centering
        \includegraphics[width=0.8\textwidth]{figures/2020-04-23-Almaden-5_6.png}
        \caption{}
        \label{fig:almaden 5 and 6}
    \end{subfigure}
    \begin{subfigure}{0.49\textwidth}
        \centering
        \includegraphics[width=0.8\textwidth]{2020-04-04-Almaden-1_2}
        \caption{}
        \label{fig:almaden 1 and 2}
    \end{subfigure}
    \caption{An asymmetry of the CNOT gate was observed across multiple pairs of qubits on \texttt{ibmq\char`_almaden}. The asymmetry increased with circuit depth as seen in (a). The performance of CNOT gates which did not exhibit an asymmetry was invariant to CNOT configuration as seen in (b). }
    \label{fig:cnot asymmetry}
\end{figure}

Interestingly, the CNOT asymmetry was even observed on qubits with relatively low error rates. In one such example we tested qubits 3 and 4 on \texttt{ibmq\char`_almaden}. Daily calibration data published by IBM on the day of experimentation showed the error rates for qubits 3 and 4 to be some of the lowest on the \texttt{ibmq\char`_almaden} quantum processor (see Appendix \ref{appendix:calibration data} for specific values). Again, the percentage of shots measured in the ground state diverged as $n$ increased, an indicator of the performance asymmetry.
This suggests that CNOT asymmetry characterization would be a useful addition to the calibration data.

A performance difference was also observable when measurement error mitigation techniques were applied. On IBM quantum processors, read operations are typically more error prone than single-qubit and two-qubit gates \cite{tannu:2019}. To determine if read errors factored into the observed asymmetry, we applied a measurement error mitigation technique and compared the mitigated results to the raw results. For these tests we used the technique described in \cite{IBM:measurement}. In general, measurement error mitigation did not adequately address CNOT asymmetry, and in some cases even exacerbated the disparity in performance between CNOT configurations. Fig. \ref{fig:measurement error mitigation} illustrates one example. In this case, measurement error mitigation improved overall CNOT performance which can be seen in the vertical translation of the mitigated plot relative to the unmitigated plot. However, the asymmetry increased when the measurement error mitigation technique was applied. Specifically, the absolute magnitude defined in Eq. \ref{asymmetry} increased from 5.55\% to 6.18\%, constituting an 11.35\% increase in asymmetry.

\begin{figure}[h]
    \centering
    \includegraphics[width=0.45\textwidth]{2020-04-23-Almaden-2_3-unmitigated-and-mitigated}
    \caption{Measurement error mitigation improved overall CNOT performance, but did not reduce CNOT asymmetry.}
    \label{fig:measurement error mitigation}
\end{figure}

\subsection{\texttt{ibmq\char`_valencia} and \texttt{ibmq\char`_paris}}

The asymmetry was also observable on \texttt{ibmq\char`_valencia} and \texttt{ibmq\char`_paris}, despite a smaller sample size
when compared to \texttt{ibmq\char`_almaden}. We observed an asymmetry on several CNOT gates, one of which can be seen in  Fig. \ref{fig:paris 5 and 6}.  The CNOT gate configured to use qubits 5 and 6 on \texttt{ibmq\char`_paris} displayed an asymmetry for $n\geq2$.

The pattern of asymmetry observed on \texttt{ibmq\char`_valencia} was markedly different than the asymmetry observed on other machines. A periodicity can be seen in Fig. \ref{fig:valencia 3 and 4} which is not present on any of the other plots. This pattern was observed only on the Q3, Q4 pair
and suggests the asymmetry may be due to a systematic error.

\begin{figure}
    \centering
    \begin{subfigure}{0.49\textwidth}
        \centering
        \includegraphics[width=0.7\textwidth]{2020-05-07-Paris-5_6}
        \caption{}
        \label{fig:paris 5 and 6}
    \end{subfigure}
    \begin{subfigure}{0.49\textwidth}
        \centering
        \includegraphics[width=0.7\textwidth]{2020-04-23-Valencia-3_4}
        \caption{}
        \label{fig:valencia 3 and 4}
    \end{subfigure}
    \caption{An asymmetry of the CNOT gate was also observed across multiple pairs of qubits on \texttt{ibmq\char`_valencia} and \texttt{ibmq\char`_paris}. The asymmetry generally increased with circuit depth, although CNOTs using Q3 and Q4 on \texttt{ibmq\char`_valencia} exhibited a non-monotonic dependence on circuit depth, as shown in part (b).}
    \label{fig:my_label}
\end{figure}

\section{Conclusions and Future Work}

In this work we presented experimental evidence showing an asymmetry in CNOT operation. Specifically, the performance of the CNOT gate on the IBM Q processors may vary depending upon the configuration of the control and target qubits. Some qubit pairs showed what appears to be a periodic variation in the asymmetry, suggesting the presence of systematic rather than purely random errors.

Based on these findings, we suggest two separate CNOT error characterizations, one for each configuration of control and target qubits, be considered when evaluating CNOT gate operation.  This additional information would better characterize the true performance of a CNOT gate and provide insightful statistics for research into quantum controls and compiler optimizations.

Our results suggest that
noise-aware circuit compilation can be improved by incorporating configuration-specific CNOT error rates into optimization objective functions. Prior work in noise-aware quantum circuit compilation has incorporated daily calibration data into optimization techniques for qubit mapping \cite{murali:2019}. The efficacy of these techniques is dependent on the accuracy of the reported gate errors. Since current calibration data from IBM characterizes CNOT gate errors with a single error rate, these prior works assume the CNOT gate operation to be invariant to control and target configuration. Our work shows this is not always the case.
Therefore we believe if stated objective functions are updated to include multiple CNOT error rates, the performance of noise-aware compilers can be further improved.

Possibilities for future work include additional analysis of the CNOT gate and the underlying CR gate to better understand the underpinnings of the observed asymmetry. In particular, performing state tomography on the output of the multi-stage circuit may give insight into the source of both systematic and random errors.
Access to the IBM Q Network processors was obtained through the IBM Q Hub at NC State.

\bibliographystyle{./bibliography/IEEEtran}


\section{Supplemental Material}
\subsection{First principles theoretical approach}
We solve the linearized Boltzmann transport equation for phonons for the non-equilibrium distribution function, $n_\lambda = n^0_\lambda + n^1_\lambda$, resulting from a small applied temperature gradient, $\nabla T$, where $n^0_\lambda$ is the Bose distribution function, and $n^1_\lambda$ is the small deviation from $n^0_\lambda$ generated by $\nabla T$. Here, $\lambda\sim\left(j, \mathbf{q}\right)$ designates the phonon mode, where $j$ is the phonon branch and $\mathbf{q}$ is the phonon wave vector. The linearized Boltzmann equation is:
\begin{align}
\mathbf{v}_\lambda\cdot\nabla T\frac{\partial n^0_\lambda}{\partial T} = \frac{\partial n_\lambda}{\partial t}\Bigg|_{\mathrm{collisions}}	\label{SI:PBE}
\end{align}
where $\mathbf{v}_\lambda$ is the phonon velocity of the mode $\lambda$, and the right-hand side is the collision term, which includes three-phonon scattering, four-phonon scattering and phonon scattering by the isotopic mass disorder. Harmonic interatomic force constants (IFCs), and anharmonic third- and fourth-order IFCs are determined using the Density Functional Theory (DFT) as implemented in Quantum Espresso.  Expressing $n^1_\lambda = n^0_\lambda\left(n^0_\lambda + 1\right)\mathbf{F}_\lambda\cdot\left(-\nabla T\right)$ allows Eq.~\ref{SI:PBE} to be recast into an integral equation for the function, $\mathbf{F}_\lambda$, which is solved numerically. The thermal conductivity is then obtained as:
\begin{align}
\kappa_{\alpha\beta} = \frac{k_B}{V}\sum_\lambda\frac{\partial n^0_\lambda}{\partial T}v_{\lambda, \alpha}F_{\lambda, \alpha}
\end{align}
where $k_B$ is the Boltzmann constant, $V$ is the crystal volume, and $\kappa_{\alpha\beta}$ is the thermal conductivity tensor for heat current flow along the Cartesian direction, $\alpha$, resulting from a temperature gradient along the direction, $\beta$. For the cubic structures considered in the present work, the thermal conductivity tensor is diagonal: $\kappa_{\alpha\beta} = k\delta_{\alpha\beta}$. Further details of the computational approach including computation scheme for IFCs, expressions for the phonon scattering rates and implementation of the solution of the phonon Boltzmann equation are given in Ref.~\cite{ravichandran_unified_2018}. All the calculations in this study are performed using norm-conserving pseudopotentials with the local density approximation (LDA) for the exchange correlation. The converged parameters used in the first principles calculations in this work for BP and SiC, such as the energy cutoffs for the DFT calculations, $\mathbf{q}$-grids used in the solution of the Boltzmann equation and the cutoffs for the harmonic, cubic and quartic IFCs, can be found in Ref.~\cite{ravichandran_phonon-phonon_2020}. The pressure in our calculations is calculated by obtaining the derivative of the Helmholtz free energy, correct to fourth-order in anharmonicity ($F_{4^{\mathrm{th}}-\mathrm{order}}$), with respect to crystal volume ($V$) at each temperature ($T$) using the expression: $P(a) = -\frac{\partial F_{4^{\mathrm{th}}-\mathrm{order}}}{\partial V}\Bigg|_{T, a} \approx -\frac{F_{4^{\mathrm{th}}}(a + \Delta a) - F_{4^{\mathrm{th}}}(a - \Delta a)}{V(a + \Delta a) - V(a - \Delta a)}\Bigg|_T$, where $a$ is the crystal lattice constant and $\Delta a \sim 0.05\%$ of $a$. The complete expression for the Helmholtz free energy, correct to fourth-order in anharmonicity of the crystal potential, can be found in  section 10 of the supplemental materials in Ref.~\cite{ravichandran_non-monotonic_2019}.
\subsection{Supplemental figures to the main text for BP}
\begin{figure*}[!ht]
\begin{center}
\includegraphics*[scale=0.325]{SI_Fig1_Aug30_2020.png}
\end{center}
\caption{Process-wise three-phonon phase space of the acoustic phonons calculated for boron arsenide [BAs - \textbf{(a)}], cubic gallium arsenide [c-GaN - \textbf{(b)}], boron phosphide [BP - \textbf{(c)}] and silicon carbide [SiC - \textbf{(d)}], organized according to decreasing mass ratio of the heavy to light constituent atoms, which correlates with a decreasing size of the frequency gap between acoustic and optic phonons (A-O gap).  The four compounds, BAs, c-GaN, BP and SiC show similar sharp dips  in the phase space for scattering between three acoustic phonons (\emph{AAA} scattering, red stars).  The large A-O gap in BAs almost completely removes the \emph{AAO} scattering (black squares) and fully exposes the dip in the \emph{AAA} scattering phase space.  In contrast, \emph{AAO} scattering processes increasingly dominate in the frequency region of the \emph{AAA} dips in c-GaN, BP and SiC as the size of the A-O gap decreases.}	\label{SI:fig_1}
\end{figure*}
\begin{figure*}[!ht]
\begin{center}
\includegraphics*[scale=0.4]{SI_Fig2_Aug29_2020.pdf}
\end{center}
\caption{\textbf{(a)} Schematic of phonon dispersions with a frequency gap between acoustic and optic phonons, $\Delta\omega_{A-O}$, and an optic phonon bandwidth of $\Delta\omega_O$. \textbf{(b)} Schematic of the corresponding three-phonon phase space for \emph{AAO} and \emph{AOO} processes illustrating the upper and lower cut-offs at $\Delta\omega_O$ and $\Delta\omega_{A-O}$, respectively.  If $\Delta\omega_{A-O} > \Delta\omega_O$, a frequency window opens in which only \emph{AAA} scattering can occur.}	\label{SI:fig_2}
\end{figure*}
\begin{figure*}[!ht]
\begin{center}
\includegraphics*[scale=0.285]{SIFig_10BPvs11BP_PhonSpec.png}
\end{center}
\caption{\textbf{(a)} Phonon dispersions of $^{10}$BP and $^{11}$BP at ambient pressure. \textbf{(b)} Phase space for \emph{AAA} and \emph{AAO} three-phonon scattering processes for $^{10}$BP and $^{11}$BP. \textbf{(c)} Scattering rates for \emph{AAA} and \emph{AAO} three-phonon processes (top panel) and spectral contributions to $\kappa$ [$\kappa (\nu)$] (bottom panel) for $^{10}$BP and $^{11}$BP. The slightly stiffer optic phonons in $^{10}$BP compared to $^{11}$BP shifts \emph{AAO} processes to higher frequencies, thereby exposing more of the \emph{AAA} phase space dip, causing weaker \emph{AAO} scattering rates and, as a result, higher $\kappa$ contributions from a narrow frequency region (shown by the blue shaded region in \textbf{(c)}) in $^{10}$BP.}	\label{SI:fig_3}
\end{figure*}
\begin{figure*}[!ht]
\begin{center}
\includegraphics*[scale=0.5]{SIFig_B11PvsMgO_100-300K_kOverk0_Aug30_2020Data.pdf}
\end{center}
\caption{Calculated $\kappa (P)$ of $^{11}$BP and MgO, scaled by their zero pressure values, $\kappa_0$, at temperatures of 100 K, 200 K and 300 K.  With decreasing temperatures, $\kappa (P)/\kappa_0$ decreases in $^{11}$BP, in striking contrast to the behavior in MgO which shows an increase in $\kappa (P)/\kappa_0$.}	\label{SI:fig_4}
\end{figure*}
\subsection{Increasing strength of four-phonon scattering in BP with increasing pressure}
Figs. 1 and 3(a) in the main text show that, as $P$ increases, inclusion of four-phonon scattering increasingly suppresses the $\kappa (P)$ of BP, and that, this suppression becomes even larger at higher $T$.  This behavior is driven by strengthening of \emph{AAAA} four-phonon scattering with increasing $P$, as shown in Fig.~\ref{SI:fig_5}.  Hydrostatic pressure stiffens the chemical bonding, resulting in increased interatomic force constants (IFCs).  These increases in the fourth-order IFCs produce corresponding increases in the four-phonon scattering matrix elements acting to increase the scattering rates of all four-phonon channels.  At the same time, the increases in second-order IFCs shift optic phonons to higher frequencies, thereby reducing their populations and reducing the scattering rates of all four-phonon processes involving optic phonons, such as \emph{AAAO}, \emph{AAOO} and \emph{AOOO} processes.  As a result, \emph{AAAA} four-phonon scattering rates are increased preferentially over all other four-phonon scattering rates, which contain at least one optic phonon. \\
\begin{figure*}[!ht]
\begin{center}
\includegraphics*[scale=0.4]{BP_4phProcwiseScatRate_300K_0to90GPa_4plots_InkscapeEdited.png}
\end{center}
\caption{Process-wise four-phonon scattering rates for BP at $T$ = 300 K at different pressures.  \emph{AAAA} scattering rates increase with increasing $P$ and dominate at high $P$. }		\label{SI:fig_5}
\end{figure*}
\pagebreak
\newpage
\clearpage
\subsection{Hidden influence of the selection rules on the $\kappa$ of SiC}
The results for the pressure and temperature-dependent phonon dispersions, three-phonon and four-phonon scattering rates, and $\kappa$ for SiC are presented in this section. Figure~\ref{SI:fig_6} (a) shows the stiffening of the longitudinal acoustic and all optic phonon branches, and the weak softening of the transverse acoustic phonon branch of SiC with pressure. Figure~\ref{SI:fig_6} (b) shows the evolution of the process-wise separated three-phonon and total four-phonon scattering rates with pressure at room temperature. The two plots show the striking similarity of the trends found for SiC and for BP (in the main text). Similarly, application of hydrostatic pressure causes a peak and drop in $\kappa$ in SiC (Fig.~\ref{SI:fig_7}) as in BP (in the main text). Furthermore, the observed temperature- and pressure-dependence of the competition between \emph{AAA} and \emph{AAO} scattering rates in BP is also seen in SiC (Fig.~\ref{SI:fig_8}). Finally, the dominant all-acoustic \emph{AAAA} four-phonon scattering channel (see Fig.~\ref{SI:fig_9}) drives the reduction in $\kappa$ of SiC at high $P$, also similar to BP.\\
\begin{figure*}[!ht]
\begin{center}
\includegraphics*[scale=0.255]{SIFig_SiC_DispersionPrDep_ScatRates_300K.png}
\end{center}
\caption{\textbf{(a)} Phonon dispersions for SiC at 300 K and different pressures. \textbf{(b)} Process-wise separated three-phonon and total four-phonon scattering rates of SiC at 300 K and different pressures.}	\label{SI:fig_6}
\end{figure*}
\begin{figure*}[!ht]
\begin{center}
\includegraphics*[scale=0.275]{SIFig_SiC_kappa_34ph_IsoPureNat.pdf}
\end{center}
\caption{Pressure dependence of the calculated thermal conductivity of isotopically enriched and naturally occurring cubic SiC scaled by its zero pressure values for different temperatures.}	\label{SI:fig_7}
\end{figure*}
\begin{figure*}[!ht]
\begin{center}
\includegraphics*[scale=0.5]{SIFig_SIC_3phProcwise_4phTotal_ScatRates_100K-1000K_0to90GPa.png}
\end{center}
\caption{Process-wise three-phonon and total four-phonon scattering rates for SiC at $T$ = 100 K and $T$ = 1000 K, and different pressures.}	\label{SI:fig_8}
\end{figure*}
\begin{figure*}[!ht]
\begin{center}
\includegraphics*[scale=0.4]{SiC_4phProcwiseScatRate_300K_0to90GPa_4plots_InkscapeEdited.png}
\end{center}
\caption{Process-wise four-phonon scattering rates for SiC at $T$ = 300 K and different pressures.  \emph{AAAA} scattering rates increase with increasing $P$ and dominate at high $P$, qualitatively similar to the behavior in BP.}	\label{SI:fig_9}
\end{figure*}
\end{document}


\section{Introduction}

Since the early days of the online advertising, the industry has faced the challenge of selling and distributing relevant ads to users at scale. The predominant solution has been to facilitate these efforts through online auctions organized by ad exchanges. The most widely used online auction format has been the generalized second-price auction (SPA), a powerful medium for advertisers to bid for reaching their target audiences with relevant product/service ads. These second-price auctions have become the engine for the online advertising industry, used for over 20 years, and driving worldwide advertiser revenues to over \$250 Billion in 2018 \cite{emarketer2018}, with more than \$100 billion in the US alone\footnote{\url{https://www.iab.com/wp-content/uploads/2019/05/Full-Year-2018-IAB-Internet-Advertising-Revenue-Report.pdf}, accessed January 2020}.

In the second-price auction, when a user places a bid, they are charged the price of the next highest competing bidder plus (usually) a penny; so for example, if they bid \$5.00, and the next highest bidder is \$2.50, then they are charged \$2.51. SPAs exhibit the Vickrey property \cite{wilkens:cinderella}, which states that bidding ones true value is a dominant strategy - meaning that this will produce the best possible payoff regardless of actions taken by other bidders. In SPAs advertisers may simply focus
on calculating the value of the incoming impression and then bidding that true value.

However, midway through 2017, this situation suddenly changed. Many ad exchanges began switching to generalized first-price auctions (FPAs), a different bidding format that was very popular in the first years of online advertising~\cite{edelman2007internet}. In January 2017 there were no known FPAs used by major display ad exchanges \cite{kitts:fpa}, while between January and December 2019, the percent of auctions running on FPAs had risen from an astonishing 40\% \cite{benes:prices} to nearly 100\% \cite{google:rtb}.
Several factors conspired to drive the industry towards the adoption of the FPA,
with the most important ones being increased demand for transparency and accountability \cite{chari1992us,sluis:guardiansuesrubicon,getintent:rtb1}.
In the first-price auction, there would be no possibility of an ad exchange improperly manipulating the clearing prices, as the price charged would always be exactly equal to the price offered. FPAs also allowed ad exchanges to capture all of the revenue from buyers since there would no longer be any discounting; this was an extremely attractive option for ad exchanges.

However, FPAs also carry some significant problems that were very notable in the early days of online advertising.
They introduce considerable complexity to bidding systems, in particular because they do not have the Vickrey property and have no strategy equilibrium (the easiest way to win an auction is to be fast and frequent in \textit{revising} bids with respect to the competitors bids). They are, thus, susceptible to system gaming strategies by different parties, which easily create volatile prices that in turn cause allocative inefficiencies \cite{edelman2007internet}.

Most importantly, in first-price auctions there is the fear of overpaying.
If an advertiser bids \$5.00 for an impression, and the next highest competing bid is \$2.50, they will get charged \$5.00 - but a strategic bid of \$3.00 would save the advertiser money. This means that the bidder has to introduce a whole new system - after calculating the value of the impression, the bidder then has to adjust their submitted bid downwards, so that if they win, they are charged only the minimum necessary to win their bid. Bidding too low though, increases the chances of not winning the auction at all, hence there is an intrinsic trade-off between the likelihood of winning and the payment in case it is won that needs to be taken into account. This practice of making a final adjustment to the bid price is referred to in the auction literature as \textit{Bid Shading}, and it is a frequent practice in many auction systems such as cattle auctions~\cite{zulehner2009bidding,pownall2013bidding}, US treasury auctions~\cite{hortaccsu2018bid} and FCC spectrum auctions~~\cite{chakravorti1995auctioning}.
The difference between the advertiser's private value and shaded bid for won auctions is called the \textit{bid surplus}, and optimization of this quantity is a major objective for the advertiser.

Online advertising auctions are dynamic, affected by a range of external factors, not to mention the arrival and departure of different bidders, bid shading is thus a very difficult problem.

In this study, we discuss a machine learning based bid shading approach for \textit{open (non-censored) first-price auctions} -- a type of auctions where there is a feedback containing minimum bid to win sent to all participants regardless of an auction outcome.
Unlike closed FPAs where the price of the highest competing bidder is censored and SPAs where the price is known only for won auctions, in the \textit{open FPAs} it is always shared.
Moreover, as in the early days of online advertising, \textit{open FPAs} are becoming a dominant type of first-price online auctions with the largest ad exchanges adopting it early.
The main benefit of this type of auctions is that they shed light on the bidding landscape and competition in such a manner that demand platforms can fairly and transparently optimize their surplus and submit appropriate bids.

The approach this study discusses learns from historical non-censored auction data using features available at the ad opportunity, so as to estimate the \textit{optimal bid shading ratio} defined as ratio of the highest competing bid and calculated bid value. This approach is successfully deployed in a major Demand Side Platform (DSP), and we demonstrate its effectiveness in a real production environment, compared against more traditional approaches, and using a range of relevant metrics in offline and online settings.

Contribution of this study are enlisted below:
\begin{itemize}
    \item We propose and characterize the problem and challenges of estimating the optimal bid shading ratio for the open (non-censored) first-price online auctions.
    \item We propose efficient Factorization Machines based approach for optimizing the bidding surplus to model the historical auction data and estimate the shading ratio in real time. To the best of our knowledge this paper is the first to describe a working bid shading algorithm for bidding systems.
    \item Improvements of total surplus show 18.57\% and 20.54\% for offline results and online results as compared to the existing models in production.
\end{itemize}

\section{Bid Shading Problem Definition, Recent Work and Challenges}
\label{sec:background}

\subsection{Canonical bidding optimization objective}
\label{sec:bidding}

In order to show ads to consumers, advertisers rely on DSPs (such as Google DoubleClick or Verizon Media DSP) to deploy highly efficient auction bidding infrastructures and run campaigns and lines on their behalf. Campaigns and lines can target an activity (such as click or conversion) and can optimize a range of objectives such are cost-per-X (view, click, mille, acquisition, install, etc.). Furthermore, each line is evaluated with key performance indicators (KPIs) which in addition to line objectives can measure total spend, win rate, etc.
Demand platforms are responsible for finding the best ad opportunities to deliver advertisers' ads, bidding for those opportunities, while spending as much as possible advertisers' budget to maximize delivery.

Online ad auction participants can submit a single bid in dollar value trying to win the opportunity to display an ad to a user. To achieve this, the following optimization objective (simplified for brevity) is formalized across all ad opportunities $i$:
\begin{equation} \label{eq:bidding}
\begin{aligned}
    & \underset{b_i}{\arg\max} \sum_{i = 1}^{N} \mathbf{I} (b_i) * v_i \\
    & \text{subject to:} \sum_{i = 1}^{N} \mathbf{I} (b_i) * c_i \leq B,
\end{aligned}
\end{equation}
where $b_i$ is the bid value, $c_i$ is the payment to the exchange if the auction is won, indicated by $\mathbf{I} (b_i)$, and $v_i$ is the expected value to the advertiser of showing the intended ad for impression opportunity $i$.
Typically, that value is defined in terms of an action associated to the ad impression, such as a click, view, conversion, etc.
For example, for conversion lines $v_i = pCVR_i *event\_value$, where $pCVR_i$ is the estimated probability of a conversion, and $event\_value$ is the monetary value of such conversion event. Finally, the constraint is placed on the optimization function to ensure that the total sum of impression costs to the advertiser $c_i$ does not exceed the budget $B$ defined for the same time window.

To solve the problem statement in (\ref{eq:bidding}),
$b_i$ is often calculated as:
\begin{equation}
	b_i = \alpha * v_i,
\end{equation}
where $\alpha$ is summarizing multiplicative parameter that control the delivery of a line, i.e. with respect to the bid shading, delivery pacing and KPIs. Depending on the auction type, $b_i$ can be unshaded (annotated with $b_i^u$) or shaded.
Bid shading, thus, plays a key role in deciding on the bid value as a part of multiplicative term $\alpha$.

\subsection{Bid Shading Maximization}
\label{sec:bid_shading}

For first-price auctions, the optimal bidding strategy $f$ attempts to  decrease (shade) the bid price as much as possible to the minimum bid necessary to win. This process is called bid shading.

The optimization function of the bid shading can be defined through maximizing the surplus:
\begin{equation} \label{eq:surplus_per_response}
\text{bid surplus} = \sum_{i=1}^{N}(b_i^u - b_i)\mathbf{I}(b_i),
\end{equation}
where bid $b_i^u$ is the unshaded bid.
It should be noted that this formulation of the problem implies that bid shading is the process of decreasing the bid value only, as increasing the bid value would negatively coact with the process of estimating the value of an impression.

\subsection{Prior bid shading algorithms}

\subsubsection{Supply-side bid shading algorithms}

The shift from second to first-price auctions was extremely rapid, and many advertisers may have difficulty calculating shading adjustments. To ease the transition to first-price, a range of ad exchanges have begun to offer "default bid shading services" to help advertisers who might struggle to implement their own algorithm. Examples include Google's ADX ``Bid Translation Service''~\cite{google:rtb}, Rubicon's ``Estimated Market Rate''~\cite{rubicon:EMR}, and AppNexus's ``Bid Price Optimization'' system ~\cite{appnexus:BPO}. These services are designed to shade high, private-value-like bids down to a price at which the advertiser continues to win at an equivalent rate, but with a less egregious price. While this is a useful outcome for supply-side platforms, since it ensures impressions continue to sell, whilst avoiding what might be otherwise be considered ``price gouging'', it is most certainly not utility maximizing for advertisers. These services are only sometimes available, and discontinuation of these services is imminent once demand platforms' bid shading efforts have become sufficiently mature.

\subsubsection{Segmented, non-linear shading algorithm}
\label{sec:non_linear_bid_shading}

One approach that is simple but effective, is to estimate the optimal bid shading ratio using parametric functions of the unshaded bid. The function parameters are either tuned manually, or algorithmically, based on a feedback mechanism that tries to maximize surplus. Parameters are trained for each unique inventory segment, which is represented as a combination of salient inventory properties such as seller exchange, top level domain, and so on. For each of these inventory segments, historical surplus data is fed back to an algorithm that adjusts parameter values iteratively, attempting to maximize surplus for that segment.
As a benchmark for the new methodology proposed in this paper, the following parametric function used in production by the DSP system was considered for both offline and online experiments:
\begin{equation}
    b_i=f(b_i^u) = \begin{cases}
            \log \frac{1 + u_1 * u_2 * b_i^u}{u_2} & u_2 > 0\\
            b_1 * b_i^u & \text{otherwise}
        \end{cases}
\end{equation}
where $u_1, u_2, b_1$ are the parameters, and $b_i^u$ is the unshaded bid. The algorithm used to adjust the parameters was a recursive least square methodology and segments were defined as the combination of exchange, top level domain, device and layout.

One drawback of the above segment-based, nonlinear approach, is that segments have to be predetermined and finding a suitable segment definition requires substantial analysis without being able to provide optimality guarantees.  Another disadvantage is that information across segments is not shared, which is a problem for segments that do not have enough traffic.

\subsubsection{Related work on bidding landscape prediction}
Several approaches were developed for modeling censored data of SPAs to estimate the bidding landscape~\cite{wu2015predicting,wu2018deep}. These approaches have a goal to directly estimate minimum bid to win, using their bid price as lower and upper bound when auction is lost and won, respectively, through distributional assumptions (such as Normal or Gamma) on bid shading landscape.
The distributional assumption can be a very limiting factor in predicting bid landscape as we show in the following section, and the studies have not shown that any distribution tested is distinctly the best. Moreover, as these approaches are developed for the SPAs it would require extending them to the use case of FPA in order to conduct comparisons.

\subsection{Modeling challenges}
\label{sec:challenges}
For traditional machine learning problems such as CTR or CVR predictions a predictor is trained on established training data, and the assumption is that available features provide sufficient richness to describe the problem. Estimating the optimal bid shading factor for an ad opportunity, on the other hand, relies on not only drawing insights from historical auction data, but also foreseeing how the other bidders would \textit{play the game}. The following challenges are often barriers when we formulate and solve bid shading problems:
\begin{itemize}

    \item \textit{A large variability of the highest competing bid price}. Knowing the highest competing bid prices as provided by open bid auctions when bidding for their inventory is crucial piece of information to know the optimal bid ratio. However, the distribution of the highest competing bid prices can vary significantly across different sets as shown in Figure~\ref{fig:mb2w_examples}. Each graph in Figure~\ref{fig:mb2w_examples} shows the empirical distribution of the highest competing bid price (per mille) for a set of fixed page and user attributes, and the differences among different sets can provide insight into the variability of the problem bid shading algorithm needs to model.
        \begin{figure} [t!]
        \centering
        \subfloat{\includegraphics[width=0.25\textwidth]{mb2w_perc_2.png}}
        \subfloat{\includegraphics[width=0.25\textwidth]{mb2w_perc_5.png}}\\
        \subfloat{\includegraphics[width=0.25\textwidth]{mb2w_perc_6.png}}
        \subfloat{\includegraphics[width=0.25\textwidth]{mb2w_perc_17.png}} \\
        \caption{Diversity of distributions of the highest competing bid prices for four fixed sets of Web page and user attributes.}
        \label{fig:mb2w_examples}
    \end{figure}

    \item \textit{Data collection}. In order to learn a bid shading strategy, it is important to log not only won but also lost bids information which increases storage requirements exponentially.
    \item \textit{One can never step into the same auction twice}. Every auction is unique. Machine learning approaches often make assumption that the test data will be generated from the same distribution that generated the training data. However, in an online ad auction system, ad inventory, end user, budget, and participating bidders are constantly changing, which poses challenges for any machine-learning approach.
    \item \textit{Many other factors contribute to the final bid price}. Bid shading is only one step of many in a complex DSP system. For example, targeting efforts, CTR/CVR predictions, pacing control, etc., attempt to make adjustments to the final bid price. A bid shading algorithm is only useful if it's resilient to noises and changes in other parts of the DSP system.
\end{itemize}

\section{emand platform Bid shading system overview}
\label{sec:system_overview}
We provide a system overview of the Verizon Media DSP which is one of the major DSPs in the US and the World with more than 200 billion daily requests received from more than 40 different ad exchanges. The DSP's production ad serving system workflow (of which the \textit{bid shading module} is part of) is shown in Figure~\ref{fig:system_overview}. On the left side of the DSP bidder we can see the sequence of events that occur for each ad opportunity generated by a user. User generates visits a web page, Supply Side Platform (SSP) generates ad request and sends it to several DSP's, the DSP will provide a first-price bid and it will receive information of minimum bid to win regardless of the bidding outcome. If bid was won, ad impression is made and logged in the data store together with the ad exchange feedback. On the backend side of the demand platform, a feature generating process is run with training machine learning model to predict optimal shading factor for different impressions once a day. The model is stored into a model file loaded by host machines in the bid shading module of the DSP bidder to be used in the following bid requests made by different SSP's.
\begin{figure}[ht!]
	\centering
	\includegraphics[width=0.48\textwidth]{bidshading_1.png}
	\caption{Overview of the bid shading in the bid generation system.}
	\label{fig:system_overview}
	\vspace{-15pt}
\end{figure}

\paragraph{Requirements and limitations for predictive approaches}
This particular system setup allows for the models to be trained on very large datasets and stored in formats optimized for efficient serving in a very limited bid request response time frame.
During the total allowable 20ms response time, multiple modules needed to go through one by one including targeting, click/conversion prediction, internal auction, value capping, bid shading and others. The current programmatic bidding system in Verizon Media does not support parallel processing of these modules, which results in strict latency constraints of each module that we take into consideration during machine learning modeling phase.

\section{Predicting the Optimal Shading Ratio}
\label{sec:methodology}

\subsection{Methodologies}
In this section we describe the Factorization Machines (FM) model whose objective is to estimate the optimal bid shading ratio defined as the bid price divided by the minimum bid to win on historical bidding data.
We then provide details of the loss function used designed to ensure surplus maximization.

\subsubsection{Factorization Machines}
Factorization Machines are a popular and efficient model used in recommendation systems. They have been very successfully used in applications involving sparse data such as CTR~\cite{pan2018field,juan2016field} and CVR~\cite{pan2019predicting} estimation outperforming popular second-degree polynomial models, while they maintain a terrific property of having only quadratic floating-operations complexity $O(N^2)$ in the space of data fields ($N$). With their high efficiency and effectiveness they are the approach that performs while complying with strict deployment restrictions discussed in Section~\ref{sec:system_overview}.

FMs learn a weight $w_i$ and an embedding vector $\bm{v}_i \in \mathbb{R}^K$ for each feature to model both its main effect and its interaction effects with other features, respectively.
The interaction between two features $i$ and $j$ is modeled as the dot product $\langle \bm{v}_i, \bm{v}_j \rangle$ of their corresponding embedding vectors $\bm{v}_i$, $\bm{v}_j$. The formula is shown in Eq.~\ref{eq:fm}:

\begin{equation}
    \Phi_{FM}\left((\bm{w},\bm{v}),\bm{x} \right)=  w_0 + \sum_{i=1}^m x_i w_i +\sum_{i=1}^{m}\sum_{j=i+1}^m x_{i}x_{j} \langle \bm{v}_{i}, \bm{v}_{j}\rangle
    \label{eq:fm}
\end{equation}
Here, $K$ is the dimension of the embedding vectors, and it's usually a small integer, in our experiments $10$.

The reason behind FMs success is that it can always learn some meaningful embedding vector for each feature as long as the feature itself appears enough times in the data, which makes the dot product a good estimate of the interaction effect of two features even if they have never or seldom occurred together in the data.

\subsubsection{Loss function and optimization}
Provided that the model's objective is to estimate the optimal bid shading ratio as its target $y$, a natural choice for a loss function was mean squared error loss:
\begin{equation}
    \mathcal{L} = \frac{1}{N}\sum_{i=1}^{N} \left( y_i - \Phi_{FM}(\bm{w},\bm{v}),\bm{x}_i) \right)^2
\end{equation}

However, the goal of the algorithm like this would be only to match the optimal bid shading ratio, whereas our goal is also to \textit{maximize the bidding surplus}, \textit{spend more},
and \textit{improve win rate}.
To that end we designed an asymmetric loss function \cite{tian2009forecasting} that:
\begin{itemize}
    \item Treats win/loss differently by penalizing more when we overshade and loose the bid, thus optimizing win rate and spend
    \item Has different loss weights for different bid requests such as penalizing less the higher bid surplus accumulation is achieved, thus optimizing the bidding surplus

    \item Has a ready capping hyperparameter that allows control over the penalty -- the more penalty when losing the bids will result in overall higher shading factor, higher win rate, and more spend. And vice versa.
\end{itemize}

We extend the simple quadratic loss function by imposing asymmetry between over-prediction and under-prediction as:
\begin{equation}
    \mathcal{L} = \frac{1}{N} \sum_{i=1}^{N} \left[(y_i - \Phi_{FM}(\bm{w},\bm{v}),\bm{x}_i))^2 * |\mathbb{I}_{\Phi_i < y_i} + \alpha| \right],
\end{equation}
where
\begin{equation}
    \mathbb{I}_{x} = \left\{\begin{matrix} 1 & \text{if } x \geq 0\\   -1 & \text{if } x < 0 \end{matrix}\right. ,
\end{equation}
and the optimization function becomes:
\begin{equation}
    \mathcal{L} = \frac{1}{N}\sum_{i=1}^{N} \left[(y_i - \Phi_{FM}(\bm{w},\bm{v}),\bm{x}_i))^2 * \left|\left\{\begin{matrix} 1 + \alpha & \text{if loss of the bid} \\  1 - \alpha & \text{if win of the bid} \\  \end{matrix}\right. \right| \right].
\end{equation}

The parameter $\alpha \in (0,1)$ indicates the degree of asymmetry. For instance, if $\alpha > 0$ the model suffers more loss from under-prediction than from over-prediction. The further $\alpha$ is from 0, the more asymmetric the model's preference between under-prediction and over-prediction is. In our definition, $\alpha$ is obtained as:
\begin{equation}
    \alpha = \text{min} (1, \text{max} (\text{opt\_surplus}, \gamma)), \gamma \in (0, 1),
\end{equation}
where
\begin{equation}
    \text{opt\_surplus = price before shading - minimum price to win}
\end{equation}
is normalized to [0, 1] scale to match the range of $\alpha$.
The hyperparameter $\gamma$ used for flooring $\alpha$ by allowing a small manual control over the shading factor predictions in the production. We thus retain control over some performance metrics not included in the loss function above to preserve its generalization across different demand platform optimization goals.

Parameters $(\bm{w},\bm{v})$ can be easily trained end-to-end through standard gradient descent approaches as the proposed asymmetric loss function preserves the smoothness of the MSE loss.

\section{Experimental evaluation}
\label{sec:experiments}
A detailed analysis in both offline and online setting of the proposed methodology as compared against several baselines across different metrics is provided here.

All algorithms have been trained on a sample of 7 days of data and evaluated on the next day using won bids only. The datasets sample approximately has 9 million auction feedback per day. We are not using the feedback from auctions that were lost, because in the current bidding system definition (Section~\ref{sec:bidding}), the bid shading function is supposed to optimize only the bid surplus given the bid value determined through other modules, thus only lowering the bid price is allowed (i.e. we do not wish to pay for opportunity more that it is worth).

For input features algorithms use as inputs to predict the optimal shading ratio for each ad opportunity we are using 13 publisher and context fields available during the ad call:
\begin{itemize}
    \item Publisher fields: $Page\_TLD$, $Subdomain$, $Publisher\_ID$ and $Request\_publisher\_ID$,
    \item Context fields: $Country\_ID$, $Day\_of\_Week$, $Hour\_of\_Day$, $Device\_type\_ID$, $App\_name$, $Is\_new\_user$, $Target\_deal\_ID$, $Layout\_ID$ and $Ad\_position\_ID$.
\end{itemize}

These 13 fields amount to hundreds of thousands of total categorical features used by the algorithms.

\subsection{Baselines}
For the machine learning approaches, we considered the following: 1) Linear regression (LR), 2) Multi Layer Perceptron regression (MLP) and 3) Factorization Machines regression (FM). All algorithms are optimized using the aforedescribed asymmetric loss function.

\subsection{Metrics}
Below we describe a spectrum of metrics used to evaluate performance of the bid shading algorithms from perspectives of regression performance and demand platform performance metrics.

\subsubsection{Regression performance metrics}
Selected metrics for evaluating regression performance of proposed algorithms in this study are:
\begin{itemize}
\item Mean Squared Error (MSE):
\begin{equation}
 MSE = \frac{1}{N}\sum_{i=1}^{n}(y_i - \hat{y}_i)^2
\end{equation}
\item Coefficient of determination ($r^2$):
\begin{equation}
 r^2 = \frac{n\Sigma_{xy} - \Sigma_x \Sigma_y}{\sqrt{(n\Sigma_{x^2} - (\Sigma_x)^2)(n\Sigma_{y^2} - (\Sigma_y)^2)}}
\end{equation}
\end{itemize}

\subsubsection{Demand platform performance metrics}
The proposed regression task metrics do not encapsulate the most important aspects of the performance of the bid shading system. As discussed in Section~\ref{sec:bid_shading}, the bid shading is actually a process of optimizing the eq.~\ref{eq:surplus_per_response}. To that end we propose the following metrics of evaluation (using the consistent notation) that we will consider the main performance metrics:
\begin{itemize}
    \item Total surplus per bid response:
    \begin{equation}
        surplus = \sum_{i=1}^{N}(b_i^u - b_i)\mathbf{I}_{B^*}(b_i)
    \end{equation}
    \item Total spend:
    \begin{equation}
        total\_spend = \sum_{i=1}^{N}b_i\mathbf{I}_{B^*}(b_i)
    \end{equation}
    \item Win rate:
    \begin{equation}
        win\_rate = \frac{1}{N}\sum_{i=1}^{N}\mathbf{I}_{B^*}(b_i)
    \end{equation}
    \item Cost per mille (CPM):
    \begin{equation}
        CPM = \frac{1}{N}\sum_{i=1}^{N}b_i\mathbf{I}_{B^*}(b_i)
    \end{equation}
\end{itemize}
Results reported on these metrics are percentage improvements or increases as compared to production's non linear bid shading algorithm (Section~\ref{sec:non_linear_bid_shading}) as absolute numbers could not be disclosed.

\subsection{Offline experiments and analysis}
In this section, we will analyze the performance of the machine learning based approaches using both regression metrics and demand platform performance metrics in an offline setting using the retrospective analysis. This will allow us to very easily characterize algorithmic performance on the available data, as well as to do a deeper dive into goal type performance analysis and fine-tuning of hyperparameters.

\subsubsection{Machine learning based approaches and their performance}
On the described dataset we train the three mentioned machine learning based approaches, characterize their regression performance and finally compare their individual improvements against the production baseline.
The regression performance is aimed to evaluate how well the algorithms fit the data, which, in the setup proposed in this study, means predicting the competitive landscape, a very difficult task indeed given a large number of demand platforms and highly variable number of external lines and ads that could potentially target the current opportunity. This is clearly reflected in lower $r^2$ values on all algorithms, however its positivity shows that models can still capture some of the large variability of the bid shading signal.

\begin{wrapfigure}{r}{0.30\textwidth}
  \begin{center}
    \includegraphics[width=0.30\textwidth]{exp_main_regression.png}
  \end{center}
  \caption{MSE and $r^2$ performance metrics for the three machine learning based approaches.}
  \label{fig:exp_main_regression}
  \vspace{-15pt}
\end{wrapfigure}
For the regression metrics (Figure~\ref{fig:exp_main_regression}) we can see that there is a distinction among the algorithms (especially in terms of the $r^2$ metric) where the MLP approach had the worst performance, while the proposed FM approach had the best performance. We suspect that the combination of existing features does not provide as much as predictive information as does modeling heterogeneity of each feature (i.e. through the feature learning layer in the FM model).

More interesting results are shown in Table~\ref{tbl:exp_main_performance}, where the three machine learning based approaches are directly compared against the non-linear bid shading used in production on surplus, spend, win rate and CPM metrics. Across the board we can see that each machine learning based approach outperforms the nonlinear bid shader on surplus, total spend and win rate.
The CPM metric is naturally increased with the increase of spend and win rate and in production environment its value is controlled in several goal type lines such as eCPM which we will show in our online experiments.  Moreover we additionally control the balance between spend and CPM increase through hyperparameters of the algorithms (in particular the $\gamma$ parameter) which we discuss in Section~\ref{sec:gamma_tuning}.
\begin{table}[h!]
\resizebox{0.48\textwidth}{!}{%
\begin{tabular}{lcccc}
\toprule
& Surplus imp. & Spend imp. & CPM inc. & Win rate impr.       \\
\midrule
LR           & 17.20\%            & 128.63\%   & \textbf{67.56\%}   & \textbf{36.40\%} \\
MLP          & 16.18\%            & \textbf{136.67\%}   & 73.85\%   & 36.12\% \\
FM           & \textbf{18.57\%}            & \textit{133.02\%}   & \textit{70.82\%}   & \textbf{36.40\%} \\
\bottomrule
\end{tabular}}%
\caption{Improvements across surplus, total spend, win rate and increase of the CPM expressed in percentages over the nonlinear bid shading production baseline as obtained by the three machine learning based approaches.}
\label{tbl:exp_main_performance}
\vspace{-10pt}
\end{table}
The highest surplus and win rate improvements are obtained with the factorization machine approach yielding second highest spend at a bit higher CPM increase than of the linear regression.

\subsubsection{Performance analysis per demand platforms' goal type}
Further characterization of the performance of different models can be done per goal types lines are optimizing towards, primarily because of the two following reasons: 1) the bid price calculation is slightly different for each goal type and 2) each goal type can have different control signals.
The goal types of particular interest on the platform with their prevalence on the test set are \textit{None} goal type with 40.5\% and five CP-X goals -- CPC (click) with 12.12\%, CPA (action) with 38.95\%, CPCV (completed view) with 0.57\%, CPViewI (viewed impressions) with 0.78\% and eCPM (mille) with 3.60\%.

Again, we separate the evaluation of algorithms into regression and platform performance metrics and show the performance of the proposed FM model as the best performing model. Regression metrics per goal type are shown in Fig.~\ref{fig:exp_fm_goal_type_regression2}. We can see that regression performance does vary per line type, where the performance is most aligned with the overall results for two dominant goal types -- None and CPA. It is expected for the performance numbers to drop for goal types that have a very small yet non-negligible share of the serving opportunities due to the inherent noise of the problem.
\begin{figure}[ht!]
	\centering
	\includegraphics[width=0.5\textwidth]{exp_fm_goal_type_regression2.png}
	\caption{MSE and $r^2$ performance metrics for the Factorization Machine algorithm broken down by the goal type.}
	\label{fig:exp_fm_goal_type_regression2}
	\vspace{-10pt}
\end{figure}

Furthermore, once we analyze the performance of the FM-based bid shading approach against the nonlinear production baseline across surplus, total spend, win rate and CPM we can clearly see that there are improvements on the key metrics.
\begin{figure}[h!]
	\centering
	\includegraphics[width=0.5\textwidth]{exp_fm_goal_type_performance.png}
	\caption{Improvements across surplus, total spend, win rate and increase of the CPM expressed in percentages over Nonlinear bid shading production baselines obtained by the Factorization Machine model broken down by the goal type.}
	\label{fig:exp_fm_goal_type_performance2}
\end{figure}
Moreover, as the highest increase in surplus we observe is for the CPA goal which typically has the highest prices, we can concur that better shading can deliver gains in practice for advertisers targeting sophisticated, rare, and high-priced goals. But also for goal types with smaller values such as eCPM, or even CPC, where shading has fewer degrees of freedom to pick and chose how to adjust the price we can clearly see that there is added value.
As a result of this analysis, we can conclude that the FM model provides a clear advantage over the production model across different goal types.

\subsubsection{Using parameter $\gamma$ to control the performance metrics}
\label{sec:gamma_tuning}
For successful deployment in production environment we allowed a control over the algorithm through parameter $\gamma$ acting as a floor to asymmetric parameter $\alpha$. By changing $\gamma$ we can have a control over models performance, i.e. choosing one with the best bid surplus or lowest CPM. This parameter is meant to select the model parameters that will be aligned with the control signals of production system.
For the FM algorithm, we searched for the best $\gamma$ parameter by increments of $0.2$ (results shown in Table~\ref{tbl:gamma_impact}).
\begin{table}[h!]
\resizebox{0.48\textwidth}{!}{%
\begin{tabular}{lllllll}
\toprule
$\gamma$ value    & 0      & 0.2    & 0.4    & 0.6    & 0.8    & 1      \\
\midrule
Surplus impr. & 17.77\%  & \textbf{18.57\%}  & 17.51\%  & 17.31\%  & 17.20\%  & 17.49\%  \\
Spend impr.   & 144.02\% & 133.02\% & 146.57\% & 148.27\% & \textbf{149.30\%} & 147.28\% \\
CPM inc.              & 75.24\%  & \textbf{70.82\%}  & 75.65\%  & 76.01\%  & 76.25\%  & 75.83\%  \\
Win rate impr.      & 39.25\%  & 36.40\%  & 40.36\%  & 41.03\%  & \textbf{41.42\%}  & 40.61\% \\
\bottomrule
\end{tabular}}%
\caption{Impact of $\gamma$ parameter evaluated across 6 values on 0.2 increment for the Factorization Machine algorithm.}
\label{tbl:gamma_impact}
\vspace{-15pt}
\end{table}
In this experiment we obtained the highest surplus improvement and the lowest CPM increase at $\gamma = 0.2$ operating point, while the highest total spend and win rate are obtained at $\gamma = 0.8$ operating point, at the cost of the highest CPM increase and lowest surplus improvement. With respect to the criteria described above, we selected $\gamma = 0.2$ as a parameter for all algorithms throughout the experiments.

Using the asymmetric loss function with $\gamma = 0.2$ compared to symmetric mean square error yielded an increase of 0.1\% and decrease of 1.2\% for MSE and $r^2$, respectively, while it provided a increase of 6\% for total surplus, and a 6.6\% decrease of CPM (with <10\% decrease in win rate and spend suggesting that this approach focuses more on higher value bids).

\subsection{Online experiments and analysis}
Offline experiments were able to show retroactive performance of the proposed bid shading approach while lacking any interaction with the other modules of the system.
Finally, we discuss performance of the proposed bid shading strategy in the online setting. The two LR and two FM models were deployed on distinct servers (we removed MLP model due to its high computational cost and poor performance in offline experiments), serving independent traffic and performance was monitored over several days.

\subsubsection{Online A/B test}
In terms of online AB tests, we aggregated results for both LR and FM models and compared them to the randomly selected traffic volume of same scale. Results obtained by the two models as compared to non-linear bid shading approach are given in Figure~\ref{fig:exp_online_ab}. Online results are similar to the offline results in terms of relative ranking of models, with the proposed FM approach being the best one. However, thanks to the online serving control the online setting provided similar or better surplus improvements but at a much lower CPM increase ($< 10\%$ compared to $> 70\%$ in the offline counterparts). This behavior is very desirable as the machine learning based bid shading module is capable of playing well into the system providing desired improvements on key metrics.
\begin{figure}[t!]
	\centering
	\includegraphics[width=0.43\textwidth]{online_results.png}
	\caption{Performance metrics improvements and increases (expressed in percentages) for LR and FM model as compared to the nonlinear bid shared collected from separate serving buckets (A/B test).}
	\label{fig:exp_online_ab}
	\vspace{-15pt}
\end{figure}

\subsubsection{Online A/A/B/B tests}
To test the stability, we separately aggregated performance of the two algorithms across different servers and thus created an A/A/B/B test. In Fig.~\ref{fig:exp_online_aa} we can see that the performance increases compared to the non linear bid shader is very similar between servers running the same algorithm, be it LR or FM model. We can thus strengthen the claim that the performance measured of any of the algorithms was not a fluke.
\begin{figure}[h!]
	\centering
	\includegraphics[width=0.46\textwidth]{online_results_aa_test.png}
	\caption{Performance metrics improvements and increases (expressed in percentages) for LR and FM model as compared to the nonlinear bid shading collected from two different servers each showing stability of performance of ML-based approach (A/A/B/B test).}
	\label{fig:exp_online_aa}
	\vspace{-15pt}
\end{figure}

Finally we ran a similar setup of experiments only per goal type (similarly to the offline counterpart) with sufficient traffic available for the best performing FM model. Figure~\ref{fig:exp_online_goal_type_aab} shows that across different goal types and across two different servers, the FM model provides very stable and high increase in bid surplus, total spend and win rate, while maintaining relatively low CPM increase ($<10\%$).
\begin{figure}[h!]
	\centering
	\subfloat{\label{fig:assym_regression}\includegraphics[width=80mm]{online_results_goal_type_a.png}}\hspace{2mm
	\vspace{-10pt}
	\subfloat{\label{fig:assym_performance}\includegraphics[width=80mm]{online_results_goal_type_b.png}}
	\caption{Performance metrics improvements and increases (expressed in percentages) for the FM model across 6 available goal types as compared to the nonlinear bid shared collected from three separate serving buckets (A/A/B/B test).}
	\label{fig:exp_online_goal_type_aab}
	\vspace{-10pt}
\end{figure}
With the last online A/A/B/B test we conclude that the machine learning based bid shading approach using the FM model is superior to any other approach considered for the demand platform's bidding system, and that it provides a significant value of allowing the strategy equilibrium in the non-censored first-price auctions.

\section{Conclusions and future directions}
\label{sec:conclusions}
In this study, we focused on describing the problem of estimating optimal shading factor using machine learning models designed for the non-censored type of first-price auctions, which are becoming the widespread type for online advertising. The bid shading problem under such circumstances becomes a point estimate problem with many opportunities and challenges. We described an approach that fitted well within DSP's production environment and performed exceptionally. The proposed FM approach with asymmetric loss function provided valuable improvements over other approaches both non-linear shading and machine learning based in both offline and online settings. The proposed approach brought improvements in estimated surplus of $18.57\%$ and $20.54\%$ in offline and online setting respectively, while providing stability of performance per goal type and also on A/A tests in an online setting.
However many challenges still remain: Different goal types weight metrics differently which could be taken into account when modeling; furthermore, using information from bids lost
could be used to improve total spend and performance for lines by modeling the bidding competition.

\balance
\bibliographystyle{ACM-Reference-Format}


\section{Introduction}
Neural networking has achieved outstanding performance in a wide range of areas such as computer vision and natural language processing (see the review article \cite{lecunbh15} and recent book \cite{goodfellowbc16}  for more background), and the question of why they work so well have naturally attracted much attention.
\cite{barron94,cybenko89,funahashi89,horniksw89} have already answered the question about a network with a single hidden layer in detail. A network with a single hidden layer can approximate any continuous function with compact support to arbitrary accuracy when the width goes to infinity. This result is referred to as \textit{the universal approximation theorem}.

This classical universal approximation theorems focused on the width increasing to infinity, namely the fat network. However, the recent tremendous success of the neural network originates from the larger and deeper network structure. The literature has reported the theoretical support for the deep neural network  \cite{bauerk19,cohenss16,eldans16,lupwhw17,linj18,mhaskarp16,ohnk19,telgarsky16,hanins18,hanin19}.

Many literature error estimates for approximators exist \cite{bauerk19,hanin19,montanellid19,yarotsky17}. Yarotsky provides the $L^\infty$-error of approximation of functions belonging to a certain Sobolev space with a fully connected ReLU network \cite{yarotsky17}. Hanin gives quantitative depth estimates on only the length of hidden layers for the rate of approximation of any continuous function by fully connected ReLU networks with a bounded width \cite{hanin19}. Also, Bauer and Kohler show that least squares estimates based on multilayer feedforward neural networks allow to circumvent the curse of dimensionality in nonparametric regression in \cite{bauerk19}.

In recent practical applications of deep learning, special network structures involving several local connections, such as CNNs, were also widely adopted.
Many works have tried to answer the question for the case of feedforward neural networks with full connection; however, not many literatures  have been answered the question yet for neural networks with a local connection \cite{zhou20nn,zhou20uni}.
 The main contributions of this paper are summarized as follows:
\begin{itemize}
\item We analyze the approximation error of the specific neural network architecture with a local connection.
\item Notably, the analytical results derived in this paper show that the locally-connected networks (CNNs) are universal function approximators.
\item Through our analysis, we provide new, interesting, useful, and helpful engineering insights to understand the approximation performance of the network architecture with local connections, and why such network architecture is working well.
\end{itemize}

\subsection{Notations}
Let $\mathbb Z_k=\{-k,-k+1,\cdots ,k\}$, $g_1(x)=\ReLU(1-|x|)$ where $ \ReLU(x)=\max\{x,0\}$ and $g_n(\xx)=\prod_{j=1}^ng_1(x_j),\xx=(x_1,x_2\cdots,x_n)\in \RR^n$ (see Figure \ref{fig:h}).
 \begin{figure*}
\centering
\includegraphics[width=0.25\textwidth] {1d}
\includegraphics[width=0.25\textwidth] {2d}
\includegraphics[width=0.2\textwidth] {3d}
\includegraphics[width=0.03\textwidth] {range}
\caption{ Left: $g_1(x)$ Middle: $g_2(\xx)$ Right: $g_3(\xx)$}
\label{fig:h}
\end{figure*}
The network architecture $(L,\mathbf n)$ consists of several hidden layers $L\in\mathbb N$ and a width vector $\mathbf n=(n_0,n_1,\cdots,n_{L+1})\in\mathbb N^{L+2}.$ A neural network with network architecture $(L,\mathbf n)$ is then any function of the form
\begin{equation}\label{eq:netwrokfunction}
\begin{array}{ll}
f_{(L,\mathbf n)}:\RR^{n_0}\to \RR^{n_{L+1}},\\
f_{(L,\mathbf n)}(\xx)=W_{L+1} \circ \ReLU \circ W_L \circ \ReLU \circ \cdots \circ  W_2 \circ \ReLU \circ  W_1(\xx)
\end{array}
\end{equation}
where $\ReLU(\xx)=(\ReLU(x_1),\ReLU(x_2),\cdots,\ReLU(x_n))$, $W_i=(w^{(i)}_{lm})\in\RR^{n_i}\times\RR^{n_{i-1}+1}$ and $W(\xx):\RR^{n_{i-1}}\to \RR^{n_i}$ is defined by
$$
\begin{array}{ll}
\displaystyle W_i(\xx)=\left(\sum_{m=1}^{n_{i-1}}w_{1m}x_{m}+w_{1,n_{i-1}+1},\sum_{m=1}^{n_{i-1}}w_{2m}x_{m}+w_{2,n_{i-1}+1},\right.\left.\cdots,\sum_{m=1}^{n_{i-1}}w_{n_im}x_{m}+w_{n_i,n_{i-1}+1}\right).
\end{array}
$$

Here we set $n_0=n$ and $n_{L+1}=1$.
To this end, we define the space of network functions with the given network architecture
$$
\mathcal F(L,\mathbf n):=\{f \mbox{ of the form \eqref{eq:netwrokfunction}}\}.
$$
\section{Main Results}

\begin{thm}\label{thm:error}
Let $f\in C^1(\RR^n)$ with $\operatorname{supp}f\subset[-1,1]^n$ and
$$
(f*g_n)_k(\xx)=\sum^{}_{\mathbf{i}\in\mathbb Z_k^n}f(\mathbf{i}k^{-1})g_{n}(k\xx-\mathbf i).
$$
Then we have
$$
|f(\xx)-(f*g_n)_k(\xx)|\le \max_{j=1,2,\cdots,n}||\partial_{x_j}f(\cdot)||_\infty n (2k)^{-1}.
$$
\end{thm}
Appendix \ref{sec:proof} presents the proof.
Let us denote by $||f||_\infty=\max\{ |f(\xx)|:\xx\in[0,1]^n\}$ and  $\lceil x\rceil=[x]+1$, where $[x]$ is a Gauss symbol.

\begin{thm}
Let $ f\in C^1(\RR^n)$ with $\operatorname{supp}f\subset[0,1]^n$ and $m,k\in\mathbb N$.
For $n\ge2$, there is a network function $f_{(L,\mathbf n)}\in \mathcal F((m+5)\lceil \log_2 n\rceil +3,(n,3n(2k+1), n(2k+1),6n(2k+1)^n,\cdots,6n(2k+1)^n ,(2k+1)^n,1)) $   such that
\begin{equation}\label{eq:estimate}
\begin{array}{ll}
\displaystyle \sup_{\xx\in[0,1]^n}|f(\xx)-f_{(L,\mathbf n)}(\xx)|\\
\le \displaystyle  ||f||_\infty3^n2^{-m}(2k+1)+ \max_{j=1,2,\cdots,n}||\partial_{x_j}f(\xx)||_\infty n (2k)^{-1},
\end{array}
\end{equation}
where $L=(m+5)\lceil \log_2 n\rceil +3$ and $\mathbf n=(n,3n(2k+1), n(2k+1),6n(2k+1)^n,\cdots,6n(2k+1)^n ,(2k+1)^n,1)$.
For $n=1$, there is a network function $f_{(2,(1,3(2k+1),2k+1,1))}\in \mathcal F(2,(1,3(2k+1),2k+1,1)) $ such that
\begin{equation*}
\sup_{\xx\in[0,1]^n}|f(\xx)-f_{(2,(1, 3(2k+1),2k+1,1))}(\xx)|\le ||f'||_\infty (2k)^{-1}.
\end{equation*}
\end{thm}
This \eqref{eq:estimate} can be written as
\begin{equation*}\label{eq:estimateO}\sup_{\xx\in[0,1]^n}|f(\xx)-f_{(L,\mathbf n)}(\xx)|=O(2^{-L/\lceil \log_2 n\rceil}k)+O((k^n)^{-1/n}).\end{equation*}
Here, $L$ is the  number of hidden layers, and $6n(2k+1)^n$ is the number of nodes of hidden layers.
If $k^2\approx 2^{L/\lceil \log_2 n\rceil}$, then
\begin{equation*}\label{eq:estimate0O}\sup_{\xx\in[0,1]^n}|f(\xx)-f_{(L,\mathbf n)}(\xx)|=O(2^{-L/(2\lceil \log_2 n\rceil)}).\end{equation*}

Theorem 2 shows that the approximation error (or accuracy) of the constructed network architecture depends on two factors: the number $L$ of the hidden layers (i.e., depth) of the network, and the parameter $k$ that is related to the maximum number $6n(2k+1)^n$ among the nodes of the hidden layers. Specifically, as shown in \eqref{eq:estimate}, the approximation error is proportional to these two factors and, more precisely, the higher increase of the number $L$ of the hidden layers than the number of nodes of the hidden layers is more efficient to achieve sufficiently high accuracy for function approximation. Notably, Theorem 2 implies that the specific network architecture with local connections, which should provide higher use than the one with the full connection because the locally-connected one is 1D CNN, is a universal function approximator.

\begin{rmk}
We can have a sharper estimate using $\operatorname{supp}f\subset[0,1]^n$.
When computing $(f*g_n)_k$, we need only half, i.e.
$$
(f*g_n)_k(\xx)=\sum^{}_{\mathbf{i}\in\mathbb Z_k^n}f(\mathbf{i}k^{-1})g_{n}(k\xx-\mathbf i)=\sum^{}_{\mathbf{i}\in(\mathbb N_k^{*})^n}f(\mathbf{i}k^{-1})g_{n}(k\xx-\mathbf i),
$$
where $\mathbb N_k^*=\{0,1,\cdots ,k\}$.
Thus, there is a network function of $f_{(L,\mathbf n)}\in \mathcal F((m+5)\lceil \log_2 n\rceil +3,(n,3n(k+1), n(k+1),6n(k+1)^n,\cdots,6n(k+1)^n ,(k+1)^n,1)) $ with
$$
\begin{array}{ll}
\displaystyle\sup_{\xx\in[0,1]^n}|f(\xx)-f_{(L,\mathbf n)}(\xx)|\le \displaystyle ||f||_\infty3^n2^{-m}(k+1)+ \max_{j=1,2,\cdots,n}||\partial_{x_j}f(\xx)||_\infty n (2k)^{-1}.
\end{array}
$$
\end{rmk}
\section{Proof}

\subsection{1-dimensional input}
To convey our idea clearer, we first construct a ReLU network with $n_{}=1$ and $2$ hidden layers with a width of at most $3(2k+1)$.
We then construct a ReLU network with a general input dimension $n$.

Let $W_{1,temp}:\RR^2\to \RR^3$ and $W_{2,temp}:\RR^3\to\RR$ with
$$
W_{1,temp}(x)=\left(x-1,x,x+1\right
$$
and
$$
W_{2,temp}(x_1,x_2,x_3)=x_1-2x_2+x_3.
$$
Then the function represented by $W_{2,temp}\circ\ReLU\circ W_{1,temp}$ is
$$
\begin{array}{ll}
W_{2,temp}\circ\ReLU\circ W_{1,temp}(x)
=\ReLU(x-1)-2\ReLU(x)+\ReLU(x+1)
=\left\{\begin{array}{ll}0\quad&\mbox{if}\quad |x|>1,\\
1-|x|\quad&\mbox{if}\quad -1<x<1,
\end{array}\right.
\end{array}
$$
which is equal to $g_1(x)$

By Theorem \ref{thm:error}, we have
$$
\sup_{-{1}\le x\le{1}}\left|f(x)- \sum^k_{i=-k}f(ik^{-1})g_1(kx-i )\right|<||f'||_\infty (2k)^{-1}.
$$

Now we contruct a feedforward neural network $(L,\mathbf n)$ with ReLU activations, input and output dimensions 1, and 2 hidden layer width at most $3(2k+1)$ whose repesent $(f*g_n)_k$: Then
 $$
\sup_{-1\le x\le1} \left|f(x)- f_{(L,\mathbf n)}(x)\right|< ||f'||_\infty (2k)^{-1}.
 $$
 Let $W_{1}:\RR\to \RR^{(2k+1)\times3}$, $W_{2}:\RR^{(2k+1)\times3}\to \RR^{2k+1},$ and $W_{3}:\RR^{2k+1}\to \RR$ be defined by
 $$
 \begin{array}{ll}
W_{1}(x)=\left(\begin{array}{ll}\displaystyle kx+k-1,kx+k,kx+k+1\\
\displaystyle kx+k-2,kx+k-1,kx+k\\
\qquad\vdots\\
\displaystyle kx-i-1,kx-i,kx-i+1\\
\qquad\vdots\\
\displaystyle kx-k-1,kx-k,kx-k+1\\
\end{array}\right),\\
 W_{2}(\xx)=\left(\begin{array}{ll}\displaystyle x_{11}-2x_{12}+x_{13}\\
\displaystyle x_{21}-2x_{22}+x_{23}\\
\qquad\vdots\\
\displaystyle x_{i1}-2x_{i2}+x_{i3}\\
\qquad\vdots\\
\displaystyle x_{k1}-2x_{k2}+x_{k3}\\
\end{array}\right)
\mbox{and}\quad
W_{3}(\xx)=\left(\begin{array}{ll}
\displaystyle \sum_{i=-k}^kf(ik^{-1})x_i\\
\end{array}\right)\mbox{(see Figure \ref{fig:n=1})}.
\end{array}
$$
Here, we express the vector in the $\RR^{(2k+1)\times3}$ to the $(2k+1)\times 3$ matrix for easy conveyance of the idea.
Then we have
$$
W_2\circ \ReLU\circ  W_1(x)=\left(\begin{array}{c}
 \displaystyle g_1(kx+k)\\
\displaystyle g_1(kx+k-1)\\
\vdots\\
\displaystyle g_1(kx-i)\\
\vdots\\
\displaystyle g_1(kx-k)\\
\end{array}\right)
$$
and thus
$$
\begin{array}{ll}W_3\circ\ReLU\circ W_2\circ \ReLU\circ W_1(x)
\displaystyle =\sum^k_{i=-k}f(ik^{-1})g_1(kx-i),
\end{array}
$$
because for a non-negative function $h:\RR\to \RR$, $\ReLU(h)=h$.
\begin{figure}
\begin{center}
\begin{tikzpicture}[scale=0.55]
\tikzset{myptr/.style={decoration={markings,mark=at position 1 with
    {\arrow[scale=2,>=stealth]{>}}},postaction={decorate}}}
   \draw (5,15) circle (0.8);
  \draw (5,13) circle (0.8);
  \draw (5,11) circle (0.8);
  \draw (5,9) circle (0.8);
  \draw (5,5) circle (0.8);
  \draw (5,7) circle (0.8);
   \draw (5,3) circle (0.8);
   \draw (5,1) circle (0.8);
    \draw (5,-1) circle (0.8);
  \draw (10,13) circle (0.8);
  \draw (10,7) circle (0.8);
   \draw (10,1) circle (0.8);
    \draw (15,7) circle (0.8);
  \draw[myptr] (5.8,15)--(9.2,13);
  \draw[myptr] (5.8,13)--(9.2,13);
  \draw[myptr] (5.8,11)--(9.2,13);
  \draw[myptr] (5.8,9)--(9.2,7);
  \draw[myptr] (5.8,7)--(9.2,7);
  \draw[myptr] (5.8,5)--(9.2,7);
  \draw[myptr] (5.8,3)--(9.2,1);
    \draw[myptr] (5.8,1)--(9.2,1);
      \draw[myptr] (5.8,-1)--(9.2,1);
  \draw[myptr] (10.8,7)--(14.2,7);
  \draw[myptr] (10.8,13)--(14.2,7);
    \draw[myptr] (10.8,1)--(14.2,7);
   \draw[->] (11.5,16.3)--(13.8,16.3);
   \node at (0,7) { $x$};
  \node at (5,-1) {{\small}  $\sum \ReLU$};
    \node at (5,1) {{\small}  $\sum \ReLU$};
      \node at (5,3) {{\small}  $\sum \ReLU$};
 \node at (5,15) {{\small}  $\sum \ReLU$};
 \node at (5,13) {{\small}  $\sum \ReLU$};
 \node at (10,13) {{\small}  $\sum \ReLU$};
 \node at (10,7) {{\small}  $\sum \ReLU$};
    \node at (5,11) {{\small}  $\sum \ReLU$};
  \node at (5,9) {{\small}  $\sum \ReLU$};
  \node at (5,7) {{\small}  $\sum \ReLU$};
  \node at (5,5) {{\small}  $\sum \ReLU$};
  \node at (10,1) {{\small}  $\sum \ReLU$};
  \node at (15,7) {{\small}  $\sum \ReLU$};
  \draw (0,7) circle (0.8) ;
  \draw[myptr] (0.8,7) -- (4.2,15);
  \draw[myptr] (0.8,7) -- (4.2,13);
  \draw[myptr] (0.8,7) -- (4.2,11);
  \draw[myptr] (0.8,7) -- (4.2,9);
  \draw[myptr] (0.8,7) -- (4.2,7);
  \draw[myptr] (0.8,7) -- (4.2,5);
  \draw[myptr] (0.8,7) -- (4.2,-1);
  \draw[myptr] (0.8,7) -- (4.2,1);
  \draw[myptr] (0.8,7) -- (4.2,3);
  \draw[->] (0.8,16.3)--(3.3,16.3);
  \draw[->] (6.7,16.3)--(8.3,16.3);
  \node at (5,16.3) {\small hidden layer};
  \node at (10,16.3) {\small hidden layer};
   \node at (15,16.3) {\small output};
   \node at (0,16.3) {\small input};
   \node at (2,16.7) {$W_1$};
    \node at (7.5,16.7) {$W_2$};
   \node at (12.5,16.7) {$W_3$};
\end{tikzpicture}
\end{center}
\caption{neural network with $n=1$ and $k=1$}\label{fig:n=1}
\end{figure}
\subsection{General-dimensional input}
As in a 1-dimensional input case, we first construct the approximate identity $g_n(\xx)$ for $n$ dimensions.
Let $W_{1,temp}(\xx):\RR^n\to \RR^{n \times 3}$ and $W_{2,temp}(\xx): \RR^{n\times 3}\to\RR^{n}$ be defined by
$$
W_{1,temp}(\xx)=\left(\begin{array}{cc}x_1-1,x_1,x_1+1\\
x_2-1,x_2,x_2+1\\
\vdots\\
x_n-1,x_n,x_n+1\end{array}\right)
$$
and
$$
 W_{2,temp}(\xx)=\left(\begin{array}{cc}x_{1,1}-2x_{1,2}+x_{1,3}\\
x_{2,1}-2x_{2,2}+x_{2,3}\\
\vdots\\
x_{n,1}-2x_{n,2}+x_{n,3}
\end{array} \right).
$$
Then we have
$$
 W_{2,temp}\circ \ReLU\circ W_{1,temp}(\xx)=\left(\begin{array}c g_1(x_1)\\g_1(x_2)\\\cdots\\g_1(x_n)\end{array}\right)
$$

Now, we construct the function $g_n$, the multiplication of $g_1$ that is presented in \cite{schmidthieber17,telgarsky16,yarotsky17}.
\begin{lem}[Lemma A3 in \cite{schmidthieber17}]\label{lem:mul}
There exists a network $\operatorname{Mult}_m^r\in \mathcal F((m+5)\lceil \log_2r \rceil,(r,6r,6r,\cdots,6r,1))$ such that
$$
\left|\operatorname{Mult}_m^r(\xx)-\displaystyle\prod_{j=1}^rx_j\right|\le 3^r2^{-m}
$$
\end{lem}
for all $\xx=(x_1,x_2,\cdots,x_r)\in[0,1]^r.$

Applying Lemma \ref{lem:mul}, we have a network $$\operatorname{AppId}_{m,temp}^n\in \mathcal F((m+5)\lceil \log_2n \rceil+2,(n,3n, n,6n,6n,\cdots,6n,1))$$ such that
$$
\displaystyle\left|\operatorname{AppId}_{m,temp}^n(\xx)-\prod_{j=1}^ng_1(x_j)\right|\le 3^n2^{-m},
$$
for all $\xx=(x_1,x_2,\cdots,x_n)\in[0,1]^n.$
Now, we construct a feedforward neural network $(L,\mathbf n)$ with ReLU activations, input dimension $n$, and output dimension 1.
Let $W_1(\xx):\RR^n\to \RR^{(2k+1)\times n \times 3}$ and $W_2(\xx): \RR^{(2k+1)\times n\times 3}\to\RR^{n\times (2k+1)}$ be defined by
$$
\begin{array}{ll}
W_1(\xx)=\left(\begin{array}{cc}\left(\begin{array}{cc}kx_1+k-1,kx_1+k,kx_1+k+1\\
kx_1+k-2,kx_1+k-1,kx_1+k\\
\vdots\\
kx_1-k-1,kx_1-k,kx_1-k+1\end{array}\right)\\
\left(\begin{array}{cc}kx_2+k-1,kx_2+k,kx_2+k+1\\
\vdots\\
kx_2-k-1,kx_2-k,kx_2-k+1\end{array}\right)\\
\vdots\\
\left(\begin{array}{cc}kx_n+k-1,kx_n+k,kx_n+k+1\\
\vdots\\
kx_n-k-1,kx_n-k,kx_n-k+1\end{array}\right)
\end{array}\right)
\quad\mbox{and}\\
\displaystyle W_2(\xx)=\left(\begin{array}{cccc}x_{1,1}-2x_{1,2}+x_{1,3},&\cdots,&x_{n,1}-2x_{n,2}+x_{n,3}\\
x_{n+1,1}-2x_{n+1,2}+x_{n+1,3},&\cdots,&x_{2n,1}-2x_{2n,2}+x_{2n,3}\\
\vdots&\vdots&\vdots\\
x_{n^2-n+1,1}-2x_{n^2-1,2}+x_{n^2-1,3},&\cdots,&x_{n^2,1}-2x_{n^2,2}+x_{n^2,3}\end{array}
 \right)\\%_{j=1,\cdots,n,\mathbf i=(i_1,\cdots,i_n)\in \mathbb Z_k^n}\\
\mbox{(see Figure \ref{fig:n=2})}.
 \end{array}
$$
Then we have
$$
\begin{array}{ll}
 W_2\circ \ReLU\circ W_1(\xx)\\
 =\left(\begin{array}{llll}g_1(x_1+k)&g_1(x_1+k-1)&\cdots&g_1(x_1-k)\\
g_1(x_2+k)&g_1(x_2+k-1)&\cdots&g_1(x_2-k)\\
&\vdots&&\\
g_1(x_n+k)&g_1(x_n+k-1)&\cdots&g_1(x_n-k) \end{array}\right).
\end{array}
$$

\begin{figure}
\begin{center}
\begin{tikzpicture}[scale=0.55]
\tikzset{myptr/.style={decoration={markings,mark=at position 1 with
    {\arrow[scale=2,>=stealth]{>}}},postaction={decorate}}}
  \draw (0,-5) circle (0.8) ;
  \draw (0,13) circle (0.8) ;
  \node at (0,13) {$x_1$};

  \node at (0,-5) {$x_2$};
  \draw (5,21) circle (0.8) ;
  \draw (5,19) circle (0.8) ;
  \draw (5,17) circle (0.8) ;
  \draw (5,15) circle (0.8) ;
  \draw (5,13) circle (0.8) ;
  \draw (5,11) circle (0.8) ;
  \draw (5,9) circle (0.8) ;
  \draw (5,7) circle (0.8) ;
  \draw (5,5) circle (0.8) ;
  \draw (5,-3) circle (0.8) ;
  \draw (5,-5) circle (0.8) ;
  \draw (5,1) circle (0.8) ;
  \draw (5,-1) circle (0.8) ;
  \draw (5,3) circle (0.8) ;
  \draw (5,-7) circle (0.8) ;
  \draw (5,-9) circle (0.8) ;
  \draw (5,-11) circle (0.8) ;
  \draw (5,-13) circle (0.8) ;

  \node at (5,15) {{\small}  $ \sum\ReLU$};
  \node at (5,9)  {{\small}  $ \sum\ReLU$};
  \node at (5,13) {{\small}  $ \sum\ReLU$};
  \node at (5,3)  {{\small}  $ \sum\ReLU$};
  \node at (5,-5) {{\small}  $ \sum\ReLU$};
  \node at (5,-7) {{\small}  $ \sum\ReLU$};
  \node at (5,5)  {{\small}  $ \sum\ReLU$};
  \node at (5,-3) {{\small}  $ \sum\ReLU$};
  \node at (5,1)  {{\small}  $ \sum\ReLU$};
  \node at (5,-1) {{\small}  $ \sum\ReLU$};
  \node at (5,11) {{\small}  $ \sum\ReLU$};
  \node at (5,7)  {{\small}  $ \sum\ReLU$};
  \node at (5,17) {{\small}  $ \sum\ReLU$};
  \node at (5,19)  {{\small}  $ \sum\ReLU$};
  \node at (5,21) {{\small}  $ \sum\ReLU$};
  \node at (5,-9)  {{\small}  $ \sum\ReLU$};
  \node at (5,-11) {{\small}  $ \sum\ReLU$};
  \node at (5,-13) {{\small}  $ \sum\ReLU$};
  \draw[myptr] (0.8,13) -- (4.2,15);
  \draw[myptr] (0.8,13) -- (4.2,13);
  \draw[myptr] (0.8,13) -- (4.2,11);
  \draw[myptr] (0.8,13) -- (4.2,7);
  \draw[myptr] (0.8,13) -- (4.2,9);
  \draw[myptr] (0.8,13) -- (4.2,5);
  \draw[myptr] (0.8,13) -- (4.2,17);
  \draw[myptr] (0.8,13) -- (4.2,19);
  \draw[myptr] (0.8,13) -- (4.2,21);
  \draw[myptr] (0.8,-5) -- (4.2,1);
  \draw[myptr] (0.8,-5) -- (4.2,3);
  \draw[myptr] (0.8,-5) -- (4.2,-3);
  \draw[myptr] (0.8,-5) -- (4.2,-1);
  \draw[myptr] (0.8,-5) -- (4.2,-5);
  \draw[myptr] (0.8,-5) -- (4.2,-7);
  \draw[myptr] (0.8,-5) -- (4.2,-11);
  \draw[myptr] (0.8,-5) -- (4.2,-9);
  \draw[myptr] (0.8,-5) -- (4.2,-13);
  \draw (10,13) circle (0.8) ;
  \draw (10,7) circle (0.8) ;

  \draw (10,-5) circle (0.8) ;
  \draw (10,1) circle (0.8) ;
  \draw (10,19) circle (0.8) ;
  \draw (10,-11) circle (0.8) ;
  \draw (15,20) circle (0.9) ;
  \draw (15,16) circle (0.9) ;
  \draw (15,12) circle (0.9) ;
  \draw (15,8) circle (0.9) ;
  \draw (15,4) circle (0.9) ;
  \draw (15,0) circle (0.9) ;
  \draw (15,-4) circle (0.9) ;
  \draw (15,-8) circle (0.9) ;
  \draw (15,-12) circle (0.9) ;

  \node at (10,13) {{\small}  $ \sum\ReLU$};
  \node at (10,7)  {{\small}  $ \sum\ReLU$};
  \node at (10,1)  {{\small}  $ \sum\ReLU$};
  \node at (10,-5) {{\small}  $ \sum\ReLU$};
  \node at (10,19)  {{\small}  $ \sum\ReLU$};
  \node at (10,-11) {{\small}  $ \sum\ReLU$};
  \node at (15,20) {{\small}  $ \sum\ReLU\atop{(\operatorname{AppId}^2_{m,\mathbf i})}$};
  \node at (15,16) {{\small}  $ \sum\ReLU\atop{(\operatorname{AppId}^2_{m,\mathbf i})}$};
  \node at (15,12)  {{\small}  $ \sum\ReLU\atop{(\operatorname{AppId}^2_{m,\mathbf i})}$};
  \node at (15,8)  {{\small}  $ \sum\ReLU\atop{(\operatorname{AppId}^2_{m,\mathbf i})}$};
  \node at (15,4) {{\small}  $ \sum\ReLU\atop{(\operatorname{AppId}^2_{m,\mathbf i})}$};
  \node at (15,0) {{\small}  $ \sum\ReLU\atop{(\operatorname{AppId}^2_{m,\mathbf i})}$};
  \node at (15,-4)  {{\small}  $ \sum\ReLU \atop{(\operatorname{AppId}^2_{m,\mathbf i})}$};
  \node at (15,-8)  {{\small}  $ \sum\ReLU\atop{(\operatorname{AppId}^2_{m,\mathbf i})}$};
  \node at (15,-12) {{\small}  $ \sum\ReLU\atop{(\operatorname{AppId}^2_{m,\mathbf i})}$};
  \draw[myptr] (5.8,15)--(9.2,13);
  \draw[myptr] (5.8,13)--(9.2,13);
  \draw[myptr] (5.8,11)--(9.2,13);

  \draw[myptr] (5.8,21)--(9.2,19);
  \draw[myptr] (5.8,19)--(9.2,19);
  \draw[myptr] (5.8,17)--(9.2,19);

  \draw[myptr] (5.8,7)--(9.2,7);
  \draw[myptr] (5.8,9)--(9.2,7);
  \draw[myptr] (5.8,5)--(9.2,7);

  \draw[myptr] (5.8,-1)--(9.2,1);
  \draw[myptr] (5.8,3)--(9.2,1);
  \draw[myptr] (5.8,1)--(9.2,1);

  \draw[myptr] (5.8,-7)--(9.2,-5);
  \draw[myptr] (5.8,-5)--(9.2,-5);
  \draw[myptr] (5.8,-3)--(9.2,-5);

  \draw[myptr] (5.8,-9)--(9.2,-11);
  \draw[myptr] (5.8,-11)--(9.2,-11);
  \draw[myptr] (5.8,-13)--(9.2,-11);

  \node at (12.5,13) { $\cdots$};
  \node at (12.5,7)  { $\cdots$};
  \node at (12.5,1)  { $\cdots$};
  \node at (12.5,-5) { $\cdots$};
  \node at (12.5,-11) { $\cdots$};
  \node at (12.5,19) { $\cdots$};
  \draw (20,4) circle (0.8) ;
  \node at (20,4) {{\small} $ \sum\ReLU$};
  \draw[myptr] (15.8,20)--(19.2,4);
  \draw[myptr] (15.8,16)--(19.2,4);
  \draw[myptr] (15.8,12)--(19.2,4);
  \draw[myptr] (15.8,8)--(19.2,4);
  \draw[myptr] (15.8,4)--(19.2,4);
  \draw[myptr] (15.8,0)--(19.2,4);
  \draw[myptr] (15.8,-4)--(19.2,4);
  \draw[myptr] (15.8,-8)--(19.2,4);
  \draw[myptr] (15.8,-12)--(19.2,4);
  \node at (5,22.3) {\small hidden layer};
  \node at (10,22.3) {\small hidden layer};
  \node at (15,22.3) {\small hidden layer};
  \node at (12.5,22.3) { $\cdots$};
  \node at (20,22.3) {\small output};
  \node at (0,22.3) {\small input};
  \draw[->] (0.8,22.3)--(3.3,22.3);
  \draw[->] (6.7,22.3)--(8.3,22.3);
  \draw[->] (16.7,22.3)--(18.8,22.3);
  \node at (2,22.7) {$W_1$};
  \node at (7.5,22.7) {$W_2$};
  \node at (17.7,22.7) {$W_{Last}$};
\end{tikzpicture}
\end{center}
\caption{neural network with $n=2$ and $k=1$}
\label{fig:n=2}
\end{figure}

Now we apply Lemma \ref{lem:mul} to a set of the $(j,i_j)$ elements of $W_2\circ \ReLU\circ W_1$ for each $\mathbf i=(i_1,i_2,\cdots,i_n)\in\mathbb Z_k^n$:
Then we have a network
$\operatorname{AppId}_m^n\in \mathcal F(L ,(n,3n(2k+1), n(2k+1),6n(2k+1)^n,\cdots,6n(2k+1)^n,(2k+1)^n))$
$(\,L=(m+5)\lceil \log_2(n(2k+1))\rceil +2 \,)$
such that for $\xx\in [0,1]^n$ and $\mathbf i\in\mathbb Z_k^n$
$$
\begin{array}{ll}
\displaystyle\max_{\mathbf i\in \mathbb Z_k^n}\left|\operatorname{AppId}_{m,\mathbf i}^n(\xx)-\prod_{j=1}^ng_1(x_j-i_j)\right|\le 3^n2^{-m}\quad\mbox{and}\quad \\
\operatorname{AppId}_m^n(\xx)=(\operatorname{AppId}_{m,\mathbf i}^n(\xx))_{\mathbf i\in\mathbb Z_k^n}
\end{array}
$$

To make the network approximating $(f*g_n)_k$, let $W_{Last}:\RR^{(2k+1)^n}\to\RR$ be defined by
$$
W_{Last}(\xx)=\sum_{\mathbf i\in \mathbb Z_k^n}f(\mathbf i k^{-1})x_{\mathbf i}.
$$
Therefore, our network is
$
W_{Last}\circ\ReLU\circ \operatorname{AppId}_m^n  \in \mathcal F(L,(n,3n(2k+1), n(2k+1),6n(2k+1)^n,\cdots,6n(2k+1)^n,(2k+1)^n,1))
$
such that
$$
\begin{array}{ll}
&|W_{Last}\circ\ReLU\circ\operatorname{AppId}_m^n(\xx)-(f*g_n)_k|\\
&\displaystyle \le\left |\sum_{\mathbf i\in \mathbb Z_k^n}f(\mathbf i k^{-1})\operatorname{AppId}_{m,\mathbf i}^n(\xx)-\sum^{}_{\mathbf{i}\in\mathbb Z_k^n}f(\mathbf{i}k^{-1})g_{n}(k\xx-\mathbf i)\right|\\
&\displaystyle = \sum_{\mathbf i\in \mathbb Z_k^n}|f(\mathbf i k^{-1})||\operatorname{AppId}_{m,\mathbf i}^n(\xx)-g_{n}(k\xx-\mathbf i)|\\
&\displaystyle \le ||f||_\infty 3^n2^{-m}(2k+1).
\end{array}
$$
By Theorem \ref{thm:error}, we have
$$
\begin{array}{ll}
&|f(\xx)-W_{Last}\circ\ReLU\circ\operatorname{AppId}_m^n(\xx)|\\
&\le|f(\xx)-(f*g_n)_k(\xx)|
+|(f*g_n)_k(\xx)-W_{Last}\circ\ReLU\circ\operatorname{AppId}_m^n(\xx)|\\
&\le \displaystyle  \max_{j=1,2,\cdots,n}||\partial_{x_j}f(\xx)||_\infty n (2k)^{-1}+ ||f||_\infty3^n2^{-m}(2k+1).
\end{array}
$$