\section{Introduction}

At a first-order phase transition, two distinct thermodynamic phases have the same free energy density
and thus coexist macroscopically. Is such phase coexistence still possible if the system contains
random disorder that locally favors one phase over the other (so-called random-$T_c$ or random-mass
disorder)? Building on earlier work for random fields \cite{ImryMa75}, Imry and Wortis
\cite{ImryWortis79} compared the possible free-energy gain from forming a domain that takes advantage
of the disorder with the free-energy cost
of its domain wall. They found that the formation of finite-size domains is favored in dimensions
$d\le 2$, even for arbitrarily weak disorder. (If the randomness breaks a continuous symmetry,
the marginal dimension is $d = 4$.)
 This destroys the macroscopic phase coexistence
and with it the first-order phase transition. After further work \cite{HuiBerker89},
this result was rigorously proven by Aizenman and Wehr \cite{AizenmanWehr89}.

What happens to a system whose first-order phase transition is unstable against disorder?
Is there an intermediate phase? Is the transition completely destroyed by smearing;
or does it become continuous? If the latter, in what universality class is the emerging
critical point?

The effects of disorder on first-order phase transitions have received less attention than the
corresponding effects on continuous transitions (see, e.g., Refs.\ \cite{Vojta06,Vojta10} for reviews).
As a result, the above questions are still being debated even for simple model systems such as
the two-dimensional Ashkin-Teller model. The classical $N$-color Ashkin-Teller model
\cite{AshkinTeller43,GrestWidom81,Fradkin84,Shankar85} is made up of $N$ Ising models
that are coupled via four-spin interactions. For more than two colors ($N>2$), the clean
Ashkin-Teller model undergoes a first-order phase transition from a magnetically ordered
phase at low temperatures to a paramagnetic phase at high temperatures.

In two dimensions, the
first-order transition cannot survive in the presence of disorder; the two-dimensional
$N$-color Ashkin-Teller model is thus a prototypical system for studying the questions posed above.
Perturbative renormalization group calculations \cite{Murthy87,Cardy96,Cardy99} predict
that disorder rounds the first-order transition of the clean-Ashkin-Teller model to a continuous
one which, somewhat surprisingly, is in the \emph{clean} Ising universality class with additional
logarithmic corrections to scaling. Early numerical simulations \cite{BellafardKatzgraberTroyerChakravarty12,BellafardChakravartyTroyerKatzgraber15}
of small systems reported nonuniversal critical behavior. However, a recent high-accuracy study
\cite{ZWNHV15} of the disordered three-color Ashkin-Teller model provided strong evidence in favor
of the scenario predicted by the perturbative renormalization group \cite{Murthy87,Cardy96,Cardy99}.

In the present paper, we test the universality of the emerging critical behavior in the
disordered Ashkin-Teller model by studying the four-color case and comparing it to the existing
simulations for three colors as well as the renormalization group findings. Specifically, we report results of large-scale
Monte Carlo simulations of the site-diluted two-dimensional four-color Ashkin-Teller model with up
to $2240^2$ sites. The paper is organized as follows. In Sec.\ \ref{sec:model}, we introduce
the model and summarize the predictions of the perturbative renormalization group.
Section \ref{sec:MC} is devoted to the Monte-Carlo simulation results. We summarize and
conclude in Sec.\ \ref{sec:conclusions}.

\section{Diluted Ashkin-Teller model}
\label{sec:model}

The two-dimensional $N$-color Ashkin-Teller model \cite{GrestWidom81,Fradkin84,Shankar85}
is defined on a square lattice of $L^2$ sites. Each lattice site $i$ contains $N$ Ising spins $S_i^\alpha=\pm 1$,
distinguished by the ``color''-index $\alpha=1\ldots N$. In the absence of disorder, the
Hamiltonian reads
\begin{equation}
H = -J \sum_{\alpha=1}^N \sum_{\langle ij \rangle} S_i^{\alpha} S_j^{\alpha}
    -\epsilon J \sum_{\alpha<\beta} \sum_{\langle ij \rangle} S_i^{\alpha} S_j^{\alpha} S_i^{\beta} S_j^{\beta}~.
\label{eq:HAT}
\end{equation}
It can be understood as $N$ identical Ising models that are coupled via their energy
densities. $\langle ij \rangle$ denotes the sum over pairs of nearest neighbor sites on the lattice;
$J>0$ is a ferromagnetic interaction; and $\epsilon\ge 0$ parameterizes the strength
of the inter-color coupling.
The clean two-color model (the original model proposed by Ashkin and Teller \cite{AshkinTeller43})
undergoes a continuous phase transition from a magnetically ordered (Baxter) phase at
low temperatures to a paramagnetic phase at high temperatures. In the Baxter phase, the spins of each color
order ferromagnetically  w.r.t. each other but the relative orientation of the colors is arbitrary (see, e.g., Ref.\ \cite{Baxter_book82}).
The critical behavior of the two-color model is nonuniversal, i.e., the critical exponents change continuously with $\epsilon$.
We are interested in the case of three or more colors for which the phase transition between the paramagnetic and Baxter phases
is of first order \cite{GrestWidom81,Fradkin84,Shankar85}.

We introduce quenched disorder into the Hamiltonian (\ref{eq:HAT}) by means of site dilution. This means,
a fraction $p$ of the lattice sites are randomly replaced by vacancies on which the spins $S_i^\alpha$ for
all colors are removed. Site dilution is a microscopic realization of random-$T_c$ or random mass disorder,
i.e., disorder that locally favors one phase over the other but does not break any of the spin symmetries.

Murthy \cite{Murthy87} and Cardy \cite{Cardy96,Cardy99} studied the phase transition of the
$N$-color Ashkin-Teller model with random-$T_c$ disorder by means of a perturbative renormalization group.
They found that the renormalization group flow on the critical $(\epsilon,\Delta)$-surface asymptotically
approaches the clean Ising fixed point $\epsilon=0,\Delta=0$. (Here, $\Delta$ is a measure of the disorder
strength.) This implies a continuous transition that is, surprisingly, in the clean 2D Ising universality
class, apart from logarithmic corrections analogous to those occurring in the disordered Ising model
\cite{DotsenkoDotsenko83,Shalaev84,Shankar87,Ludwig88}. The following finite-size scaling behavior
has been derived \cite{MazzeoKuhn99,HTPV08,KennaRuizLorenzo08}.
At criticality ($T=T_c$), the specific heat shows a characteristic double-logarithmic dependence
on the system size,
\begin{equation}
C \sim \ln \ln L~.
\label{eq:C_lnln}
\end{equation}
Magnetization and susceptibility at $T_c$ (averaged over all $N$ colors) behave as
\begin{eqnarray}
M  &\sim& L^{-\beta/\nu}\, (1+b_M/\ln L)~,
\label{eq:M_ln}\\
\chi &\sim& L^{\gamma/\nu}\,(1+b_\chi/\ln L)~,
\label{eq:chi_ln}
\end{eqnarray}
where $\gamma/\nu=7/4$ and $\beta/\nu=1/8$ as in the clean Ising model, and $b_M$ and $b_\chi$ are constant.
Any quantity $R$ of scale dimension zero behaves as
\begin{eqnarray}
R &=& R^\ast \,(1+b_R/\ln L)~,
\label{eq:R_ln}\\
dR/dT  &\sim& L^{1/\nu}(\ln L)^{-1/2}\, [1+O(1/(\ln L))]
\label{eq:dRdT_ln}
\end{eqnarray}
with the clean Ising exponent $\nu=1$. The finite-size scaling forms (\ref{eq:R_ln}) and (\ref{eq:dRdT_ln})
hold for the Binder cumulants
\begin{equation}
g_\textrm{av} = \left[1 - \frac {\langle m^4\rangle}{3\langle m^2\rangle^2} \right]_\textrm{dis}~, \quad
g_\textrm{gl} = 1 - \frac {[\langle m^4\rangle]_\textrm{dis}}{3[\langle m^2\rangle]^2_\textrm{dis}} ~.
\label{eq:Binder}
\end{equation}
Here, $\langle \ldots \rangle$ denotes the thermodynamic (Monte-Carlo) average while $[\ldots ]_\textrm{dis}$
stands for the disorder average. We need to distinguish average and ``global'' versions of this quantity, depending
on when the disorder average is performed.
The scaling forms (\ref{eq:R_ln}) and (\ref{eq:dRdT_ln})
also hold for the correlation length ratios $\xi_\textrm{av}/L$ and $\xi_\textrm{gl}/L$.
In our simulations, the correlation lengths are computed from the second moment of the spin correlation function
$G(\mathbf{r}) = (1/L^2) \sum_{i,j,\alpha} \langle S_i^\alpha S_j^\alpha \rangle \delta(\mathbf{r}-\mathbf{r}_{ij})$ \cite{CooperFreedmanPreston82,Kim93,CGGP01}.
They can be obtained efficiently from the Fourier transform $\tilde G(q)$ of the correlation function:
\begin{eqnarray}
\xi_\textrm{av} &=& \left[ \left(\frac{\tilde G(0) -\tilde G (q_\textrm{min})}{q_\textrm{min}^2 \tilde G(q_\textrm{min})} \right)^{1/2} \right]_\textrm{dis}~,
\label{eq:xi_av}\\
\xi_\textrm{gl} &=&  \left(\frac{[\tilde G(0) -\tilde G (q_\textrm{min})]_\textrm{dis}}{q_\textrm{min}^2 [\tilde G(q_\textrm{min})]_\textrm{dis}} \right)^{1/2}~.
\end{eqnarray}
Here, $q_\textrm{min}=2\pi/L$ is the minimum wave number that fits into a system of linear size $L$.

\section{Monte Carlo simulations}
\label{sec:MC}
\subsection{Method and overview}
\label{subsec:Method}

To simulate the thermodynamics of the four-color Ashkin-Teller model (\ref{eq:HAT}),
we employ an embedding algorithm analogous to that used in Refs.\ \cite{WisemanDomany95,ZWNHV15}.
It is based on a simple observation. If the spins of colors $\alpha=2,3$ and 4 are fixed,
the Hamiltonian (\ref{eq:HAT}) is equivalent to an Ising model for the spins $S_i^{(1)}$
with effective interactions $J_{ij}^\textrm{eff} = J+ \epsilon J ( S_i^{(2)}S_j^{(2)} + S_i^{(3)}S_j^{(3)} + S_i^{(4)}S_j^{(4)} )$.
This (embedded) Ising model can be simulated using any valid Monte-Carlo method.
Analogous embedded Ising models can be constructed for the spins $S_i^{(2)}$, $S_i^{(3)}$, and $S_i^{(4)}$.
By combining Monte-Carlo updates for all four embedded Ising models we obtain a valid
Monte-Carlo method for the entire Ashkin-Teller Hamiltonian.

Using this algorithm, we simulate systems with sizes from $35^2$ to $2240^2$ sites
with periodic boundary conditions and dilution $p=0.3$. All results are averaged over a large
number of disorder configurations (10,000 to 500,000), details will be given below.
We combine Wolff single-cluster updates \cite{Wolff89} with Swendsen-Wang multi-cluster
updates \cite{SwendsenWang87}. The latter help equilibrating small isolated clusters of sites
that occur for larger dilutions. Specifically, each full Monte Carlo sweep consists of
a Swendsen-Wang sweep for each color and a Wolff sweep (a number of single-cluster flips
such that the total number of flipped spins for each color equals the number of lattice sites).
We use 100 full Monte Carlo sweeps for equilibrating each sample (disorder configuration)
and 200 sweeps for measuring observables (one measurement per sweep).
The actual equilibration times are much shorter \cite{ZWNHV15}.
Biases in the observables due to the short measurement periods are overcome by using
improved estimators \cite{ZWNHV15}.

The Wolff and Swendsen-Wang algorithms are only valid as long as all effective interactions
of the embedded Ising models,
$J_{ij}^\textrm{eff} = J+ \epsilon J ( S_i^{(2)}S_j^{(2)} + S_i^{(3)}S_j^{(3)} + S_i^{(4)}S_j^{(4)} )$,
 are not negative. In the worst case, the term in parenthesis can take the value $-3$.
The coupling constant $\epsilon$ therefore must not exceed 1/3. In the production runs, we use the largest possible value,
$\epsilon=1/3$, because this leads to a strong first-order transition in the clean case.

\subsection{Results}
\label{subsec:Results}

To find the phase transition of the four-color Ashkin-Teller model with coupling $\epsilon=1/3$
and dilution $p=0.3$ we perform a series of simulation runs using linear system sizes $L=35$ to 2240.
The number of disorder realizations ranges from 500,000 for the smallest systems to 10,000 for the largest
ones. As usual, the critical temperature $T_c$ can be determined from the crossings of the Binder cumulant vs.
temperature curves for different system sizes (and analogously from the crossings of the reduced
correlation length curves $\xi/L$ vs. $T$).

The Binder cumulant $g_\textrm{gl}$ as a function of temperature $T$ is shown in Fig.\ \ref{fig:bindercrossing},
\begin{figure}
\includegraphics[width=\columnwidth]{Graph1.eps}
\caption{Binder cumulant $g_\textrm{gl}$ vs.\ temperature $T$ of the site-diluted two-dimensional four-color Ashkin-Teller model with $p=0.3$ and $\epsilon=1/3$
        for different linear system sizes $L$.  The shift towards higher temperatures of the crossing point with increasing $L$
        is caused by corrections to scaling.}
\label{fig:bindercrossing}
\end{figure}
and Fig.\ \ref{fig:xicrossing} presents the reduced correlation length $\xi_\textrm{gl}/L$ as a function of temperature $T$.
\begin{figure}
\includegraphics[width=\columnwidth]{Graph2.eps}
\caption{Reduced correlation length $\xi_\textrm{gl}/L$ vs.\ temperature $T$
        for different linear system sizes $L$.  With increasing $L$, the crossings move to lower $T$, again
        indicating corrections to scaling.}
\label{fig:xicrossing}
\end{figure}
Similar graphs can be created for the quantities $g_\textrm{av}$ and $\xi_\textrm{av}/L$. In all cases, the crossing points
of their curves for different $L$ move with increasing $L$. This indicates that corrections to the leading scaling behavior
are important for the studied system sizes. To find the true, asymptotic value of the critical temperature,
we therefore need to extrapolate the crossing points to infinite system size. To do so, we find the crossing temperature $T_x(L/2,L)$
of the $g_\textrm{gl}$ vs.\ $T$ curves for system sizes $L/2$ and $L$ as well as the analogous crossing temperatures
for   $g_\textrm{av}$, $\xi_\textrm{gl}/L$, and $\xi_\textrm{av}/L$. The resulting dependence of the crossing temperatures on the system
size is presented in Fig.\ \ref{fig:TxvsL}.
\begin{figure}
\includegraphics[width=\columnwidth]{Graph3.eps}
\caption{Crossing temperatures $T_x(L/2,L)$ vs.\ inverse system size $1/L$. The solid lines are fits to
$T_x(L/2,L) =T_c + a L^{-b}$ yielding $T_c=1.9639(3)$. The error bars of $T_x$ are about a symbol size for the smallest $L$, they become
much smaller with increasing $L$.}
\label{fig:TxvsL}
\end{figure}
The figure shows that all crossings approach the same temperature as $L$ increases. Fits to the heuristic relation
$T_x(L/2,L) =T_c + a L^{-b}$ yield a critical temperature of $T_c=1.9639(3)$ (the number in parentheses is an estimate of the error
of the last digits).

To determine the critical behavior, we now analyze the finite-size scaling properties of various observables right at $T_c$.
The data analysis follows Ref.\ \cite{ZWNHV15} and is based on the finite-size scaling forms (\ref{eq:C_lnln}) to (\ref{eq:dRdT_ln}).
Figure \ref{fig:CV} presents a double-logarithmic plot of specific heat $C$ vs.\ system size $L$.
\begin{figure}
\includegraphics[width=\columnwidth]{Graph4.eps}
\caption{Double-logarithmic plots of the specific heat $C$ as well as the slopes $L^{-1} d\ln(\xi_\textrm{gl}/L)/dT$
     and $L^{-1} d\ln(\xi_\textrm{av}/L)/dT$ vs.\
     system size $L$ at the critical temperature of $T_c=1.9639$. The error bars are smaller than the symbol sizes.
     The solid lines are fits to
     $a\ln[b\ln(cL)]$ for $C$ and $a[\ln(bL)]^{-1/2}$ for the slopes.}
\label{fig:CV}
\end{figure}
The specific heat increases more slowly than a power law with $L$, as indicated by the downward curvature of the graph. In agreement with the
prediction (\ref{eq:C_lnln}), the data can be fitted well to the form $a \ln[b \ln(cL)]$ over the entire
size range. The fit is of high quality, giving a reduced error sum $\bar\chi^2 \approx 0.8$. (The reduced error sum of a fit of $n$ data points
($x_i,y_i$) to a function $f(x)$ having $q$ fit parameters is defined as $\bar\chi^2 = 1/(n-q)\sum_i [y_i-f(x_i)]^2/\sigma_i^2$
where $\sigma_i$ is the standard deviation of $y_i$.) For comparison, we have also attempted to fit the specific heat to the simple logarithmic
form $a\ln(bL)$ and to $C_\infty - a L^{-b}$. The latter function corresponds to power-law scaling with a negative finite-size scaling
exponent $\alpha/\nu=-b$. The power-law fit is of poor quality with a reduced error sum $\bar\chi^2 \approx 4.3$, and the simple logarithmic fit
is completely off, giving $\bar\chi^2 \approx 420$.

Figure\ \ref{fig:CV} also presents the slopes $d\ln(\xi_\textrm{gl}/L)/dT$ and
$d\ln(\xi_\textrm{av}/L)/dT$ of the normalized correlation lengths vs.\ temperature curves at $T_c$. We have divided out the power law $d\ln(\xi/L)/dT \sim L$
of the clean Ising universality class to make the corrections predicted in eq.\ (\ref{eq:dRdT_ln})
more easily visible. The figure demonstrates that these corrections are not of power-law type, instead the data can be fitted well by the predicted logarithmic form
$a [\ln(bL)]^{-1/2}$,  yielding reduced error sums $\bar \chi^2$ of about 0.3 for $\xi_\textrm{av}$ and 0.9 for $\xi_\textrm{gl}$.

In addition, we study the system size dependence of the magnetization and the magnetic susceptibility at the critical temperature.
According to eqs.\ (\ref{eq:M_ln}) and (\ref{eq:chi_ln}), these quantities are predicted to follow the clean Ising power
laws with additive logarithmic corrections. We again divide out the clean Ising power laws and plot the resulting quantities,
viz., $M\, L^{1/8}$ and $\chi\, L^{-7/4}$ in Fig.\ \ref{fig:Mchi}.
\begin{figure}
\includegraphics[width=\columnwidth]{Graph5.eps}
\caption{Double-logarithmic plots of $M\, L^{1/8}$ and $\chi\,L^{-7/4}$ vs.\ $L$
     at the critical temperature of $T_c=1.9639$. The solid lines are fits to
     $a[1+b/\ln(cL)]$, as suggested by eqs.\ (\ref{eq:M_ln}) and (\ref{eq:chi_ln}).}
\label{fig:Mchi}
\end{figure}
The figure demonstrates that the corrections to the Ising universality class are not of power-law form.
Moreover, they are very weak, in particular for the magnetization where
they change the value by less than 1\% over the entire system size range. To capture this small correction,
we had to increase the number of disorder configurations significantly (500,000 for all system sizes)
to reduce the error bars. The larger numerical effort limits these simulations to system sizes between
$L=35$ and 800. Both  $M\, L^{1/8}$ and $\chi\, L^{-7/4}$ can be fitted well to the function
$a[1+b/\ln(cL)]$ predicted by the renormalization group results (\ref{eq:M_ln}) and (\ref{eq:chi_ln}).
The reduced error sums $\bar \chi^2$ are approximately 0.4 for the susceptibility and 0.3 for the
magnetization.

\section{Summary and conclusions}
\label{sec:conclusions}

In summary, we have carried out large-scale Monte Carlo simulations of the
site-diluted two-dimensional four-color Ashkin-Teller model. Our results confirm that the first-order
phase transition occurring in the undiluted (clean) model is rounded by the disorder,
as required by the Aizenman-Wehr theorem \cite{AizenmanWehr89}.

We have used finite-size scaling of the magnetization, magnetic susceptibility, specific heat
and correlation length to determine the universality class of the emerging continuous phase transition.
All our data agree very well with the results of the perturbative renormalization group
\cite{Murthy87,Cardy96,Cardy99} which predicts critical behavior in the clean two-dimensional
Ising universality class, but with logarithmic corrections similar to those occurring in
the two-dimensional disordered Ising model. These findings agree with those of extensive
high-accuracy simulations of the three-color case \cite{ZWNHV15}. Consequently, they provide strong evidence
for the universality of the critical behavior of the disordered Ashkin-Teller model.

Possible reasons for the discrepancies between our results (Ref.\ \cite{ZWNHV15} and the present paper)
and those of the earlier simulations \cite{BellafardKatzgraberTroyerChakravarty12,BellafardChakravartyTroyerKatzgraber15}
were discussed in detail in Ref.\ \cite{ZWNHV15}. Here, we just reiterate the main points:
Our systems are much larger, up to $2240^2$ sites (compared to just $32^2$ sites
in Ref.\ \cite{BellafardKatzgraberTroyerChakravarty12} and $128^2$ sites in Ref.\
\cite{BellafardChakravartyTroyerKatzgraber15}). This suggests that the earlier simulations
were not in the asymptotic regime, especially for weak disorder where the crossover from
the clean first-order transition to the disordered continuous one is slow. Note, however, that
a discrepancy exists already for the \emph{clean} three-color model. The clean phase diagram
of Ref.\ \cite{ZWNHV15}, which was determined and verified using three independent Monte Carlo algorithms,
coincides with older results by Grest and Widom  \cite{GrestWidom81} but disagrees with
Ref.\ \cite{BellafardKatzgraberTroyerChakravarty12}.

From a more general point of view, the renormalization groups results for the random-$T_c$
Ising \cite{DotsenkoDotsenko83,Shalaev84,Shankar87,Ludwig88} and Ashkin-Teller \cite{Cardy96} models
suggested a kind of ``super-universality'' of critical points in two-dimensional disordered systems.
This idea was initially supported by computer simulations of disordered Ising\cite{ADSW90,WSDA90a}, Ashkin-Teller\cite{WisemanDomany95},
and Potts\cite{WisemanDomany95,ChenFerrenbergLandau95} models as well as interface arguments
\cite{KSSD95}.
Numerical sudies of the disordered $q$-state Potts model\cite{JacobsenCardy98,ChatelainBerche98}
showed, however, that the finite-size scaling exponent $\beta/\nu$ differs from the Ising value and varies with $q$.
Moreover, the phase transition in the random-bond Blume-Capel model was also found to display complex non-Ising behavior, at least for
strong disorder \cite{MBHF09,MBHFP10,TheodorakisPanagiotisFytas12}.

Recently, the disordered quantum Ashkin-Teller spin chain has attracted lots of interest
because it serves as a paradigmatic model for studying disorder effects at first-order
\emph{quantum} phase transitions. As a quantum version of the Aizenman-Wehr theorem
has been proven \cite{GreenblattAizenmanLebowitz09}, the first-order character of the transition
in the clean problem must change upon the introduction of disorder.
Recent strong-disorder renormalization group approaches predict continuous transitions
governed by infinite-randomness critical points in different universality classes, depending on
the coupling strength $\epsilon$ \cite{GoswamiSchwabChakravarty08,HrahshehHoyosVojta12,Barghathietal14}.
Furthermore, in the case of two colors an exotic strong-disorder infinite-coupling phase \cite{HHNV14}
is predicted to appear for large $\epsilon$. These predictions can be tested by generalizations
of our Monte Carlo method to the  $(1+1)$-dimensional quantum case.
Some work along these lines is already in progress.

It may also be interesting to investigate the Ashkin-Teller model defined on a topologically
disordered lattice such as the random Voronoi-Delaunay lattice (see, e.g., \cite{OBSC_book00}).
Recent work \cite{BarghathiVojta14} has shown that the Imry-Ma argument does not hold for these
lattices because a topological constraint suppresses the disorder fluctuations. This leaves
open the possibility that the first-order phase transition survives in the presence of such
topological disorder.

This work was supported in part by the NSF under Grant Nos.\ DMR-1205803 and DMR-1506152.


\subsection{Need for modeling interaction}
We first show that if we ignore the interaction effect, we will get an
increasingly worse solution as the dimension of the problem
increases. We consider the following simple optimization problem for
our purposes.
\begin{equation}
\label{simple_QCQP}
\begin{aligned}
& \underset{\mathbf{x}}{\text{Minimize}} & &  -\mathbf{x}^T\mathbf{p} + \frac{\gamma}{2} \mathbf{x}^T\mathbf{x}\\
& \text{subject to}&  &\mathbf{x}'\mathbf{r}  \leq P\\
& & & 0 \leq \mathbf{x} \leq 1
\end{aligned}
\end{equation}
where the truth is $\mathbf{p} = -\mathbf{Q}_p \mathbf{x}$ and $\mathbf{r} = \mathbf{Q}_r \mathbf{x}$ for
positive definite matrices $\mathbf{Q}_p, \mathbf{Q}_r$ as defined in
Section \ref{sec:interact}. We solve the true optimization problem,
\begin{equation}
\label{truth}
\begin{aligned}
& \underset{\mathbf{x}}{\text{Minimize}} & &  \mathbf{x}^T\left( \mathbf{Q}_p + \frac{\gamma}{2}\mathbf{I}\right) \mathbf{x}\\
& \text{subject to}&  &\mathbf{x}'\mathbf{Q}_r\mathbf{x}  \leq P\\
& & & 0 \leq \mathbf{x} \leq 1
\end{aligned}
\end{equation}
to get ${\boldsymbol{x}}^*$. Now, ignoring the dependency structure we use $\mathbf{p} =
diag(\mathbf{Q}_p)$ and $\mathbf{r} = diag(\mathbf{Q}_r)$ to solve the
problem \eqref{simple_QCQP} to get $\hat{\mathbf{x}}$. Note that the
diagonals are just the single slot estimates.

We simulate two random instances of the single slot estimate and
create $diag(\mathbf{Q}_p), diag(\mathbf{Q}_r)$ as a decreasing function of the slot
position. We then simulate the rest of the matrix $\mathbf{Q}_p$ and $\mathbf{Q}_r$
according to the structure in Section \ref{sec:interact}. For
different values of sample size $n$, we calculate the relative error
as
\begin{align*}
err(n) = \frac{\hat{\mathbf{x}}^T\left( \mathbf{Q}_p
+ \frac{\gamma}{2}\mathbf{I}\right) \hat{\mathbf{x}} - (\mathbf{x}^*)^T\left( \mathbf{Q}_p
+ \frac{\gamma}{2}\mathbf{I}\right) \mathbf{x}^* }{(\mathbf{x}^*)^T\left( \mathbf{Q}_p
+ \frac{\gamma}{2}\mathbf{I}\right) \mathbf{x}^*}
\end{align*}
Figure \ref{fig:obj_val} shows the log of two different objective values. It
can be seen that we start getting increasing worse solutions as the
dimension increases. The fluctuations are because of the randomness in
the simulations.

\begin{figure}[!h]
\centering
\includegraphics[width = \linewidth]{obj_val2}
\caption{\label{fig:obj_val}
The black line shows the true objective value obtained from solving
problem \eqref{truth}. The blue line shows the objective value using
the solution from the ignored dependency problem.}
\end{figure}

Figure \ref{fig:err} show the log of the relative error as a function
of the dimension. It can be seen from our simulation, that we get an
average relative error of about $6 \times 10^8$ which is extremely
high of a cost to pay for ignoring the interaction effect.

\begin{figure}[!h]
\centering
\includegraphics[scale = 0.55]{rel_err}
\caption{\label{fig:err}
Plotting the log(relative error) against the dimension of the problem}
\end{figure}

\begin{table*}[!t]
\caption{\label{tab:obj_time}The optimal objective value and convergence time}
\begin{center}
 \begin{tabular}{|| c | c | c | c | c |c |c ||} \hline $n$ & Our
 method & Sampling on $[0,1]^n$ & Sampling on $\mathbb{S}^n$ & SDP
 relaxation & RLT relaxation & Exact \\[0.5ex]
\hline\hline
\multirow{2}{*} {5} & 3.00 & 2.99 & 2.95  & 3.07  & 3.08 & 3.07 \\
& (4.61s) & (4.74s) & (6. 11s) & (0.52s) & (0.51s) & (0.49) \\
 \hline
\multirow{2}{*} {10} & 206.85 & 205.21 & 206.5  & 252.88  & 252.88 & 252.88 \\
& (5.04s) & (5.65s) & (5.26s) & (0.53s) & (0.51s) & (0.51) \\
\hline
\multirow{2}{*} {20} & 6291.4 & 4507.8 & 5052.2  & 6841.6  & 6841.6 & 6841.6 \\
& (6.56s) & (6.28s) & (6.69s) & (2.05s) & (1.86s) & (0.54) \\
\hline
\multirow{2}{*} {50} & 99668 & 15122 & 26239 & $1.11\times10^5$  & $1.08\times10^5$ & $1.11\times10^5$ \\
& (15.55s) & (18.98s) & (17.32s) & (4.31s) & (2.96s) & (0.64) \\
\hline
\multirow{2}{*} {100} & $1.40\times10^6$ & 69746 & $1.24\times10^6$  & $1.62\times10^6$  & $1.52\times10^6$& $1.62\times10^6$ \\
& (58.41s) & (1.03m) & (54.69s) & (30.41s) & (15.36s) & (2.30s) \\
 \hline
\multirow{2}{*} {1000} & $2.24\times10^7$ & $8.34 \times 10^6$ & $9.02\times10^6$  & \multirow{2}{*} {NA}  & \multirow{2}{*} {NA} & \multirow{2}{*} {NA} \\
& (14.87m) & (15.63m) & (15.32m) &  &  & \\
\hline
\multirow{2}{*} {$10^5$} & $3.10\times10^8$ & $7.12 \times 10^7$ & $ 8.39\times10^7$  & \multirow{2}{*} {NA}  & \multirow{2}{*} {NA} & \multirow{2}{*} {NA} \\
& (25.82m) & (24.59m) & (27.23m) &  &  & \\
 \hline
\multirow{2}{*} {$10^6$} & $3.91\times10^9$ & $2.69 \times 10^8$ & $7.53\times10^8$  & \multirow{2}{*} {NA}  & \multirow{2}{*} {NA} & \multirow{2}{*} {NA} \\
& (38.30m) & (39.15m) & (37.21m) &  &  & \\
 \hline
 \end{tabular}
\end{center}
\end{table*}

\subsection{Comparative Study of QCQP Solutions}
We compare our proposed technique to the current state-of-the-art
solvers of QCQP. Specifically we compare it to the SDP and RLT
relaxation procedures. For small enough problems, we also compare our
method to the exact solution by interior point methods. Furthermore,
we provide empirical evidence to show that our sampling technique is
better than other simpler sampling procedures such as uniform sampling
on the unit square or on the unit sphere and then mapping it
subsequently to our domain as in Algorithm \ref{algo:simulate}. We
begin by considering a very simple QCQP for the form
\begin{equation}
\label{test_qcqp}
\begin{aligned}
& \underset{\mathbf{x}}{\text{Minimize}} & &  \mathbf{x}^T \mathbf{A} \mathbf{x}\\
& \text{subject to}&  &(\mathbf{x} - \mathbf{x}_0)^T \mathbf{B} (\mathbf{x} - \mathbf{x}_0)  \leq \tilde{b},\\
& & & \mathbf{l} \leq \mathbf{x} \leq \mathbf{u}
\end{aligned}
\end{equation}
We randomly sample $\mathbf{A}, \mathbf{B}, \mathbf{x}_0$ and $\tilde{b}$. The lower bound, $\mathbf{l}$ and upper bounds $\mathbf{u}$ are chosen in a way such that they intersect the ellipsoid. We vary the dimension $n$ of the problem and tabulate the final objective value as
well as the time taken for the different procedures to converge in
Table \ref{tab:obj_time}. Throughout our simulations we have chosen $\eta = 2$ and the number of optimal points as $N = \max (1024, 2^m) $, where $m$ is the smallest integer such that $2^m \geq 10n$.

Note that even though the SDP and the interior point methods converge very efficiently for small values of $n$, they cannot scale to values of $n \geq 1000$, which is where the strength of our method becomes evident. From Table \ref{tab:obj_time} we observe that the relaxation procedures SDP and RLT fail to converge within an hour of computation time for $n \geq 1000$, whereas all the approximation procedures can easily scale up to $n = 10^6$ variables.  We can further notice that SDP performs slightly better than RLT for higher dimensions.

Furthermore, our procedure gives the best approximation result when compared to the remaining two sampling
schemes. Lemma \ref{lem:step1} shows that if the both the objective
values are the same then we indeed get the exact solution. To see how much the approximation deviates from the truth, we also tabulate the relative error, i.e. $\norm{\mathbf{x}^*(N) - \mathbf{x}^*}/\norm{\mathbf{x}^*}$ for each
of the sampling procedures in Table \ref{tab:err}. We omit SDP and RLT results in Table \ref{tab:err} since both of them produce a solution very close to the exact minimum for small $n$. From the results in Table \ref{tab:err} it is clear that our procedure gets the smallest relative error compared to the other sampling schemes that we tried.

\begin{table}[!t]
\caption{\label{tab:err}The Relative Error : $\norm{\mathbf{x}^*(N) - \mathbf{x}^*}/\norm{\mathbf{x}^*}$}
\begin{center}
 \begin{tabular}{|| c | c | c | c ||}
 \hline
 $n$ & Our method & Sampling on $[0,1]^n$ & Sampling on $\mathbb{S}^n$  \\[0.5ex]
\hline\hline
5 &0.0615 	& 0.0828 & 0.0897\\
\hline
10&0.0714 	&	 0.1530&0.1229\\
\hline
20 &0.0895&	0.2455&	0.2368\\
\hline
50&0.3352&	3.8189&	1.0472\\
\hline
100&0.8768	&13.3709	&2.0849\\
 \hline
 \end{tabular}
\end{center}
\end{table}

\subsection{Choice of $\mathbf{Q}_p$ and $\mathbf{Q}_r$}
\label{sec:interact}
The choice of these matrices is a hard problem in practice since we do
not know the exact dependency structure of $\mathbf{p}$ and $\mathbf{r}$ on
$\mathbf{x}$. However, due to the presence of enormous empirical data for
clicks as well as complaints and dismisses, estimation of
$\mathbf{Q}_p, \mathbf{Q}_r$ is not a hard problem, if there exists a parametric
form. Since $\mathbf{Q}_p$ and $\mathbf{Q}_r$ can have similar parametric form,
without loss of generality we explain the form through $\mathbf{Q}_p$.

Throughout we assume $\mathbf{Q}_p = Diag( \mathbf{Q}_i, i = 1, \ldots, n)$,
i.e. the probabilities across the users are independent and each
$\mathbf{Q}_i$ is a $JK
\times JK$ dimensional positive definite matrix which we parametrize
below. We assume that all $K$ slots are viewable by a user, and the
chance of an event at any slot depends on all the rest. Moreover we
hide the $i$ in the subscript for notational simplicity.

Let us begin by introducing some more notation. For any $j \in \{1,
\ldots, J\}$, let $\mathbf{p}_j, \mathbf{x}_j$ denote $K$-dimensional vectors
comprising of elements $\{p_{jk}\}_{k=1}^K$ and $\{x_{jk}\}_{k=1}^K$
respectively. Let $\mathbf{Q}_{jj'}$ denote the $K \times K$ dependence
matrix between event $j$ and observation $j'$. This leads to the
following expression
\begin{align*}
\mathbf{p} = \left(
  \begin{tabular}{c}
  $\mathbf{p}_1$\\
  $\mathbf{p}_2$ \\
  $\vdots$\\
  $\mathbf{p}_J$
  \end{tabular}
\right) =
\begin{bmatrix}
    \mathbf{Q}_{11} & \mathbf{Q}_{12} & \mathbf{Q}_{13} & \dots  & \mathbf{Q}_{1J} \\
    \mathbf{Q}_{21} & \mathbf{Q}_{22} & \mathbf{Q}_{23} & \dots  & \mathbf{Q}_{2n} \\
    \vdots & \vdots & \vdots & \ddots & \vdots \\
    \mathbf{Q}_{J1} & \mathbf{Q}_{J2} & \mathbf{Q}_{J3} & \dots  & \mathbf{Q}_{JJ}
\end{bmatrix}
 \left(
  \begin{tabular}{c}
  $\mathbf{x}_1$\\
  $\mathbf{x}_2$ \\
  $\vdots$\\
  $\mathbf{x}_J$
  \end{tabular}
  \right)
\end{align*}

Now we characterize each $\mathbf{Q}_{jj'}$. Note that for any $j,k$, $p_{jk}$
is the probability of an event conditional of the fact that $x_{jk} =
1$. Thus, it depends on $\mathbf{x}_j$ only through $x_{jk}$. Hence it is easy
to see that for any $j \in \{1, \ldots, J\}$,
\begin{align}
\label{eq:q_diag}
\mathbf{Q}_{jj} = \tilde{p}_j \mathbf{I}_K
\end{align}
where $\tilde{p}_j$ is the prior event probability conditional on
seeing item $j$ irrespective of its position, and $\mathbf{I}_K$ is the $K
\times K$ identity matrix. The values of $\tilde{p}_j$ for $j = 1,
\ldots, J$ can be estimated from empirical data.

Now for any $j \neq j'$, we consider the dependency between the event
corresponding to item $j$ and the observation of $j'$. When we are
considering $p_{jk}$, since we know that $x_{jk} = 1$, we must have
that the contribution from $x_{j'k} = 0$ for $j' \neq j$. Moreover,
that is the only coefficient which is zero, since we let $p_{jk}$
depend on $x_{j'k'}$ for $k' \neq k$. Thus, we have,

\begin{align}
\label{eq:q_off_diag}
\mathbf{Q}_{jj'} = \begin{bmatrix}
    0 & a_{12} & a_{13} & \dots  & a_{1K} \\
    a_{21} & 0 & a_{23} & \dots  & a_{2K} \\
    \vdots & \vdots & \vdots & \ddots & \vdots \\
    a_{K1} & a_{K2} & a_{K3} & \dots  & 0
    \end{bmatrix}.
\end{align}
We estimate each of the $a_{\ell\ell'}$ using empirical data. Finally
to bring in symmetry in the problem, we assume that
\begin{align}
\label{eq:q_sym}
\mathbf{Q}_{j'j} = \mathbf{Q}_{jj'}^T
\end{align}
for all $j$ and $j'$. This condition will often be true for similar items, or for dissimilar items which have similar effects on one another (whether mutually synergistic or antagonistic). However, we also acknowledge that there will be certain cases where this condition will not hold, and that will make the problem non-convex.

Combining the structure in equations \eqref{eq:q_diag},
\eqref{eq:q_off_diag} and \eqref{eq:q_sym}, we get the complete
parametric form of the matrix $\mathbf{Q}$.

\subsection{Practical considerations}

If we assume the parametric form of $\mathbf{Q}$ as given in Section
\ref{sec:interact}, the postive-definiteness of $\mathbf{Q}$ is a concern since
we estimate each of the $a_{\ell\ell'}$ from the data. An easy fix
to this problem, is to add $\eta \mathbf{I}_{JK}$ to $\mathbf{Q}$ for some appropriate
value of $\eta$. This enforces that $\mathbf{Q}$ to be positive-definite by
only slightly increasing the elements corresponding to
$\tilde{p}_{j}$. Thus instead of using the estimated $\mathbf{Q}$ directly we
use,

\begin{align}
\label{eq:tilde_Q}
\tilde{\mathbf{Q}} := \begin{cases}
\mathbf{Q} + (-\lambda_1(\mathbf{Q}) + \epsilon)I_{JK} & \text{ if } \lambda_1(\mathbf{Q}) < 0\\
\mathbf{Q} & \text{ otherwise}
\end{cases}
\end{align}
where $\lambda_1(\mathbf{Q})$ denotes the minimum eigenvalue of $\mathbf{Q}$ and
$\epsilon > 0$ is arbitrary.

An advantage of using this technique from a computational point of
view is that since each $\mathbf{Q}$ is of small dimension ($JK$ being small),
the minimum eigenvalue calculation can be done extremely fast and can
be made highly parallel for different users $i$. Secondly, this
guarantees that $\mathbf{Q}_p = Diag(\tilde{\mathbf{Q}}_i, i = 1,\ldots,n)$ is also
positive-definite.

\subsection{Dependence of Parameters}
\label{sec:dep_p_r}
Whenever we are modeling the interaction effect, we have assumed that
vectors $\mathbf{p}$ and $\mathbf{r}$ are independent. However, there may arise problems
where that is not the case. For example, if we consider our original
problem \eqref{eq:orig_qp}, then $\mathbf{r} = {\boldsymbol{\$}} = \mathbf{p}\cdot\mathbf{c}$ and hence they are
clearly not independent. Here, it is clear that if we construct
$\mathbf{Q}_r = \mathbf{Q}_p \cdot \mathbf{c}$, where $\mathbf{Q}_p \cdot \mathbf{c}$ is a matrix whose $\ell$-th row is
the $\ell$-th row of $\mathbf{Q}_p$ multiplied by $c_\ell$, then $\mathbf{Q}_r$ loses
its positive-definiteness since it is not even symmetric. For a general dependency structure, we have the following result.

\begin{theorem}
Suppose we have an optimization problem of the form
\eqref{eq:dual_sol_easy}. Moreover, assume that $\mathbf{r} = f(\mathbf{p})$ for some
function $f$, that maps $\mathbb{R}^n \rightarrow \mathbb{R}^n$ and
let $\mathbf{Q}_p = Diag(\tilde{\mathbf{Q}}_i, i = 1,\ldots,n)$ where $\tilde{\mathbf{Q}}_i$'s
are generated by \eqref{eq:tilde_Q}. If there exists a
positive-definite matrix $\mathbf{Q}_r$ and a vector $\mathbf{c}_r$ such that
\begin{align}
\label{eq:cond1}
\mathbf{Q}_r \mathbf{x}  - 2 \mathbf{Q}_r \mathbf{c}_r = f(-\mathbf{Q}_p \mathbf{x})
\end{align}
for all $\mathbf{x}$, then there exists a convex QCQP formulation of
\eqref{eq:dual_sol_easy}.
\end{theorem}

\begin{proof}
Assume there exists a positive-definite matrix $\mathbf{Q}_r$ and $\mathbf{c}_r$ such
that such that \eqref{eq:cond1} holds. Thus, using the dependence
structure and modeling the interaction effect, we can write,
\begin{align*}
\mathbf{x}^T\mathbf{r} &= \mathbf{x}^Tf(\mathbf{p}) = \mathbf{x}^Tf(-\mathbf{Q}_p \mathbf{x})\\
&= \mathbf{x}^T\mathbf{Q}_r\mathbf{x} - 2\mathbf{x}^T\mathbf{Q}_r \mathbf{c}_r\\
&= (\mathbf{x} - \mathbf{c}_r)^T \mathbf{Q}_r (\mathbf{x} - \mathbf{c}_r) - \mathbf{c}_r^T\mathbf{Q}_r \mathbf{c}_r.
\end{align*}
Thus we have,
\begin{align*}
\mathbf{x}^T\mathbf{r} \leq P \Longleftrightarrow (\mathbf{x} - \mathbf{c}_r)^T \mathbf{Q}_r (\mathbf{x} - \mathbf{c}_r) \leq P +
\mathbf{c}_r'\mathbf{Q}_r \mathbf{c}_r.
\end{align*}
The result follows from observing the fact that since $\mathbf{Q}_r$ is
positive definite, $\mathbf{c}_r'\mathbf{Q}_r \mathbf{c}_r > 0$.
\end{proof}

\begin{remark}
Few simple functions $f$ for which the convex reformulation Theorem
holds include the linear shift operators and the positive scaling
operators. The exact characterization of the function class may be
hard and is beyond the scope this paper.
\end{remark}

\section{Introduction}
\label{sec:intro}
\input{Introduction}

\section{Feed Ranking as a MOO Problem}
\label{sec:problem_def}
\input{Problem_Defn}

\section{An Efficient Solution}
\label{sec:soln_prob}
\input{Problem_Soln}

\section{Modeling Interaction}
\label{sec:modeling_interaction}
\input{Modeling_Interaction}

\section{Solution to the QCQP}
\label{sec:sol_qcqp}
\input{Soln_QCQP}

\section{Experimental Results}
\label{sec:expt_results}
\input{Expt_results}

\section{Discussion}
\label{sec:disc}
\input{Discussion}

\section*{Acknowledgement}
We were sincerely like to thank Prof. Art Owen for the useful discussions and his help in providing some $(t,m,s)$-nets for our experiments.

\bibliographystyle{abbrv}

\subsection{Dual Solution}
The local constraints ${\boldsymbol{x}}_i \in \mathcal{K}_i$ for all $i = 1, \ldots, n$ can be
combined together to get a constraint like $\mathbf{K}\mathbf{x} \leq \mathbf{b}$, where $\mathbf{K} =
diag(\mathbf{K}_i, i = 1,\ldots, n)$. Thus with this notation we can transform
the problem \eqref{eq:orig_qp} as
\begin{equation}
\label{dual_sol}
\begin{aligned}
& \underset{x}{\text{Minimize}} & & -\mathbf{x}^T\mathbf{p} + \frac{\gamma}{2}\mathbf{x}^T\mathbf{x}\\
& \text{subject to} & & \mathbf{x}^T{\boldsymbol{\$}} \geq R \\
&  &  &\mathbf{x}^T\mathbf{d} \geq I\\
& & & \mathbf{K}\mathbf{x} \leq \mathbf{b}
\end{aligned}
\end{equation}
Writing the Lagrangian of the above problem, we have
\begin{align*}
L(\mathbf{x}, \mu_0, \mu_1, \boldsymbol{\eta} ) &= -\mathbf{x}^T\mathbf{p} + \frac{\gamma}{2}\mathbf{x}^T\mathbf{x} + \mu_0(R -
\mathbf{x}^T{\boldsymbol{\$}})\nonumber \\ & \;\;\; + \mu_1(I - \mathbf{x}^T\mathbf{d}) + \boldsymbol{\eta}^T(\mathbf{K}\mathbf{x} - \mathbf{b}).
\end{align*}
Finding the minimum with respect to $x$ we see get,
\begin{align*}
\hat{\mathbf{x}} = \frac{1}{\gamma}\left(\mu_0{\boldsymbol{\$}} + \mu_1\mathbf{d} + \mathbf{p} - \mathbf{K}^T\boldsymbol{\eta} \right.).
\end{align*}
Plugging this back into the Lagrangian we get,
\begin{align*}
L(\hat{\mathbf{x}}, \mu_0, \mu_1, \boldsymbol{\eta} ) &= \frac{\gamma}{2}\hat{\mathbf{x}}^T\hat{\mathbf{x}}
-\hat{\mathbf{x}}^T(\mu_0{\boldsymbol{\$}} + \mu_1\mathbf{d} + \mathbf{p} - \mathbf{K}^T\boldsymbol{\eta} )\\ &\;\;\; + \mu_0R +
\mu_1I - \boldsymbol{\eta}^T\mathbf{b} \\ & = -\frac{\gamma}{2}\hat{\mathbf{x}}^T\hat{\mathbf{x}} + \mu_0R +
\mu_1I - \boldsymbol{\eta}^T\mathbf{b}.
\end{align*}
Writing $\boldsymbol{\xi} = (R,I,-\mathbf{b}^T)^T, \;\; \mathbf{A} = ({\boldsymbol{\$}} : \mathbf{d} : -\mathbf{K}^T)$ and $\mathbf{y} =
(\mu_0, \mu_1, \boldsymbol{\eta})$ the above problem simplifies to
\begin{equation}
\label{dual_sol_simplified}
\begin{aligned}
& \underset{\mathbf{y}}{\text{Maximize}} & & \boldsymbol{\xi}^T\mathbf{y} -\frac{1}{2\gamma}(\mathbf{p} +
  \mathbf{A}\mathbf{y})^T(\mathbf{p} + \mathbf{A}\mathbf{y})\\ & \text{subject to} & & \mathbf{y} \geq 0
\end{aligned}
\end{equation}
which can be further simplified to
\begin{equation}
\label{dual_sol_final}
\begin{aligned}
& \underset{y}{\text{Minimize}} & & \frac{1}{2}\mathbf{y}^T\mathbf{M}\mathbf{y} - \mathbf{y}^T(\boldsymbol{\xi} -
  \mathbf{A}^T\mathbf{p}/\gamma)\\ & \text{subject to} & & \mathbf{y} \geq 0
\end{aligned}
\end{equation}
where $\mathbf{M} = \mathbf{A}^T\mathbf{A}/\gamma$. We solve this by the operator splitting
algorithm \cite{parikh2014block} to obtain $\mu_0$ and $\mu_1$.

\begin{remark}
  We cannot apply block splitting algorithm \cite{parikh2014block} to
  the dual problem since the objective function is no longer
  separable. We may apply it to the primal but obtaining the dual
  solution may be not be possible in all situations. For details, we
  refer to \cite{basu_user}.
\end{remark}

The ability to solve this problem on a large scale completely depends
on the sparsity of matrix $\mathbf{M}$.

\begin{definition}
Let $\mathbf{A}$ be any $m\times n$ matrix. We define the sparsity ratio,
$\psi(\mathbf{A})$ as $k$ if the number of non-zero entries in the matrix is
$kmn$.
\end{definition}

From previous experimentation we were able to solve the dual problem
having $n$ variables through operator splitting if $\psi(\mathbf{M})$ is of
order $O(1/n)$. We shall show in this setup too, $\psi(\mathbf{M})$ is of the
same order.

\begin{lemma}
\label{lem:sparse}
Consider the dual problem for the multi-slot optimization problem as
given in \eqref{dual_sol_final}. Let $\mathcal{J}_s$ and $\mathcal{J}_I$ denote
the sets of sponsored and impression important items. If there are $K$
slots, $J$ total items and $n$ users, then,
\[
\psi(\mathbf{M}) = \frac{1 + n(J + \beta + K(3 + \beta) + 7JK)}{(1+ nJ + nK + nJK)^2}.
\]
where $\beta = |\mathcal{J}_s| + |\mathcal{J}_I|$.
\end{lemma}
\begin{proof}
Note that, using the previous notation we can write,
\begin{align*}
\mathbf{M} = \frac{1}{\gamma} \mathbf{A}^T\mathbf{A} = \frac{1}{\gamma} \begin{bmatrix}
    {\boldsymbol{\$}}^T{\boldsymbol{\$}} & {\boldsymbol{\$}}^T\mathbf{d} & -{\boldsymbol{\$}}^T\mathbf{K}^T\\
    \mathbf{d}^T{\boldsymbol{\$}} & \mathbf{d}^T\mathbf{d} & -\mathbf{d}^T\mathbf{K}^T  \\
    -\mathbf{K}{\boldsymbol{\$}} & -\mathbf{K}\mathbf{d} & \mathbf{K}\Kb^T
\end{bmatrix}.
\end{align*}
Without loss of generality we assume that
$\mathcal{J}_s \cap \mathcal{J}_I \neq \emptyset$. Then the top-left $2 \times 2$
principal sub matrix of $M$ contains 4 non-zero terms. Otherwise only
the diagonals are positive and we only lose a count of 2 to the total
number of non-zero terms. Note that for each $i = 1, \ldots, n$ we can
write $\mathbf{K}_i^T$ as
\begin{align*}
\mathbf{K}^T_i = \begin{bmatrix}
    \mathbf{I}_{JK} & -\mathbf{I}_{JK} & \mathbf{B}^T & -\mathbf{B}^T & \mathbf{C}^T & -\mathbf{C}^T
\end{bmatrix},
\end{align*}
where,
\begin{align*}
\mathbf{B} = \begin{bmatrix}
    \mathbf{I}_{K} & \mathbf{I}_{K} & \ldots & \mathbf{I}_{K}
\end{bmatrix}_{K \times JK},\qquad
\mathbf{C} = \mathbf{I}_{J} \otimes 1_{K}^T,
\end{align*}
and $\mathbf{I}_n$ denotes the identity matrix of dimension $n$ and $\otimes$
denote the Kronecker product. Note that the matrices in $\mathbf{K}_i$
denotes the constraints for the $i$-user of optimization
problem \eqref{orig_problem}. Let $\mathrm{card}(\mathbf{A})$ denotes the
number of non-zero entries in a matrix $\mathbf{A}$. Now it is easy to see
that,
\begin{align}
\label{eq:s1}
\mathrm{card}(\mathbf{K}_i {\boldsymbol{\$}}_i) &= 2( K + |\mathcal{J}_s| + K|\mathcal{J}_s|)\\
\label{eq:s2}
\mathrm{card}(\mathbf{K}_i \mathbf{d}_i)  &= 2( K + |\mathcal{J}_I| + K|\mathcal{J}_I|).
\end{align}
To complete the proof, we note that,
\begin{align}
\label{eq:s3}
\mathrm{card}(\mathbf{I}_{JK}) &= \mathrm{card}(\mathbf{B}) = \mathrm{card}(\mathbf{C}) = JK,\\
\label{eq:s4}
\mathrm{card}(\mathbf{B}\Bb^T)  &= K, \;\;\mathrm{card}(\mathbf{C}\Cb^T)  = J, \;\;\mathrm{card}(\mathbf{B}\mathbf{C}^T)  = JK.
\end{align}
Now, using the structure is $\mathbf{K}\Kb^T$, and equations \eqref{eq:s3} and \eqref{eq:s4}, we get
\begin{align}
\label{eq:s5}
\mathrm{card}(\mathbf{K}\Kb^T) = 4n(J + K + 7JK).
\end{align}
The result follows from using \eqref{eq:s1}, \eqref{eq:s2}
and \eqref{eq:s5} and observing that the total number of entries in
the matrix $\mathbf{M}$ is $(2 + 2nK + 2nJ + 2nJK)^2$.
\end{proof}

From Lemma \ref{lem:sparse}, it is easy to see that $\psi(\mathbf{M}) =
O(1/(nJK))$. Thus, for example if we are able to solve the single slot
for 2 million users, having $K = 5$ and $J = 10$ we can solve this
multi-slot problem at the same time and same accuracy for 40000 users.

\subsection{Dual to Primal Trick}
The main idea behind this technique is to exploit the KKT conditions
of the problem, so that we can write the primal solution as a function
of the dual solution and the input parameters. To explain this
technique, let us start by rewriting the problem given in
\eqref{eq:orig_qp} as
\begin{equation}
\label{eq:dual_trick}
\begin{aligned}
& \underset{\mathbf{x}}{\text{Minimize}} & & -\mathbf{x}^T\mathbf{p} + \frac{\gamma}{2}\mathbf{x}^T\mathbf{x} +
  \mathds{1}_{\mathcal{K}_1 \times \mathcal{K}_2 \times \ldots \times \mathcal{K}_n}(\mathbf{x}) \\ & \text{subject
    to} & & \mathbf{x}^T{\boldsymbol{\$}} \geq R, \\ & & &\mathbf{x}^T\mathbf{d} \geq I,
\end{aligned}
\end{equation}
where $\mathcal{K}_1 \times \mathcal{K}_2 \times \ldots \times \mathcal{K}_n$ denotes the domain of
$\mathbf{x}$ and $\mathds{1}_{\mathcal{C}}(\mathbf{x}) = 0$ if $\mathbf{x} \in \mathcal{C}$ and $\infty$ otherwise. Introducing the dual variables $\mu_0$ and $\mu_1$ for the two
constraints we can write the Lagrangian as follows.
\begin{align*}
L(\mathbf{x}, \mu_0, \mu_1) = -\mathbf{x}^T\mathbf{p} + \frac{\gamma}{2}\mathbf{x}^T\mathbf{x} + \mu_0(R - \mathbf{x}^T{\boldsymbol{\$}}) +
\mu_1(I - \mathbf{x}^T\mathbf{d}),
\end{align*}
for $\mathbf{x} \in \mathcal{K}_1 \times \mathcal{K}_2 \times \ldots \times \mathcal{K}_n$. From this it is
easy to see that the optimal solution satisfies
\begin{align}
\label{primal_proj}
\mathbf{x}^{opt} = \Pi_{\mathcal{K}_1 \times \mathcal{K}_2 \times \ldots \times \mathcal{K}_n}
\left(\frac{1}{\gamma}(\mu_0{\boldsymbol{\$}} + \mu_1\mathbf{d} + \mathbf{p})\right),
\end{align}
where $\Pi_{\mathcal{C}}(\cdot)$ denotes the projection function onto
$\mathcal{C}$. Note that, since we have a product domain, we can write
this as
\begin{align*}
{\boldsymbol{x}}^{opt}_i = \Pi_{\mathcal{K}_i} \left(\frac{1}{\gamma}(\mu_0{\boldsymbol{\$}} + \mu_1\mathbf{d} +
\mathbf{p})_i\right)\;\; \forall \;\; i = 1,\ldots,n.
\end{align*}
Thus once we know $\mu_0$ and $\mu_1$ by
solving \eqref{dual_sol_final}, $x_i$ can be obtained simply by the
projection onto $\mathcal{K}_i$ which is a small quadratic problem and can
be solved fast and efficiently.

\subsection{Deterministic Serving Plan}
Here we present the algorithm to generate the serving plan for user
$i$ using the optimal probabilistic solution $x_{ijk}$. We know for
each slot $k$ the items follow a multinomial distribution with
probability $x_{ijk}$ for $j = 1, \ldots, J$. While serving we have
the added criteria that no two items can be repeated. Thus for each
slot $k$ we sample from $Mult( \{x_{ijk}\}_{j=1}^J)$ and we resample
if we get the same item for two slots. This guarantees that that
serving plan obeys the optimal serving distribution. The detailed
steps are written out in Algorithm \ref{algo:plan1}.

\begin{algorithm}
\caption{Deterministic Serving Plan}\label{algo:plan1}
\begin{algorithmic}[1]
\State \text{Input : Optimal primal solution $\mathbf{x}$}
\State \text{Output : Deterministic serving plan for each user}
\State \text{Set $S$ empty $n \times K$ matrix}
\For{$i = 1:n$}
\State $Sample = \emptyset$
\For{$k = 1:K$}
\State \label{sample-step} Pick $j \sim Mult( \{x_{ijk}\}_{j=1}^J)$
\If{$j \in Sample$}
\State Go to step \ref{sample-step}.
\Else
\State Add $j$ to Sample.
\EndIf
\EndFor
\State $S[i, ] = Sample$
\EndFor
\State \Return $S$
\end{algorithmic}
\end{algorithm}
\vspace{5em}

\subsection{QCQP to QP Approximation}
Here we describe the linearization technique to convert the quadratic
constraint into a set of $N$ linear constraints. The main idea behind
this approximation, is the fact that given any convex set in the
Euclidean plane, there exists a convex polytope that covers the
set.

Let us begin with a few notations. Let ${\cal P}$ denote the optimization
problem \eqref{eq:QCQP}. Let $\mathbf{x} \in {\mathbb R}^{s}$. Define,
\begin{align}
\label{eq:qc}
{\cal S} := \{ \mathbf{x} \in \mathbb{R}^s : (\mathbf{x} - \mathbf{b})^T\mathbf{B}(\mathbf{x} - \mathbf{b})  \leq \tilde{b} \}.
\end{align}
Let $\partial{\cal S}$ denote the boundary of the ellipsoid ${\cal S}$.

To generate the $N$ linear constraints for this one quadratic
constraint, we generate a set of $N$ points, $\mathbf{x}_1, \ldots, \mathbf{x}_N \in
\partial{\cal S}$. The sampling technique to select these points are given
in Section \ref{sec:sampling}. Corresponding to these $N$ points we
get the following set of $N$ linear constraints,
\begin{align}
\label{eq:lin_const}
(\mathbf{x} - \mathbf{b})^T \mathbf{B}(\mathbf{x}_j - \mathbf{b}) \leq \tilde{\mathbf{b}} \qquad \text{ for } j = 1, \ldots, N.
\end{align}
If we look at it geometrically, it is not hard to see that each of
these linear constraints are just tangent planes to ${\cal S}$ at $\mathbf{x}_j$
for $j = 1, \ldots, N$. Figure \ref{fig:qcqp_qp_lin} shows a set of
six linear constraints for a ellipsoidal feasible set in two
dimensions. Thus, using these $N$ linear constraints we can write the
approximate optimization problem, ${\cal P}(N)$, as follows.

\begin{equation}
\label{eq:QP_via_lin}
\begin{aligned}
& \underset{x}{\text{Minimize}} & &  (\mathbf{x} - \mathbf{a})^T\mathbf{A} (\mathbf{x} - \mathbf{a}) \\
& \text{subject to}&  &(\mathbf{x} - \mathbf{b})^T\mathbf{B}(\mathbf{x}_j - \mathbf{b})  \leq \tilde{b}\qquad \text{ for } j = 1, \ldots, N \\
& & & \mathbf{C}\mathbf{x} \leq \mathbf{c}.
\end{aligned}
\end{equation}
Thus, instead of solving ${\cal P}$, we solve ${\cal P}(N)$ for a large enough
value of $N$. Note that as we sample more points, our approximation
gets more and more accurate to the original solution.

\begin{figure}[!h]
\centering
\includegraphics[scale = 0.6]{qcqp_qp_lin.pdf}
\caption{\label{fig:qcqp_qp_lin} Converting a quadratic constraint
  into linear constraints. The tangent planes through the 6 points
  $\mathbf{x}_1, \ldots, \mathbf{x}_6$ create the approximation to ${\cal S}$.}
\end{figure}

\subsection{Sampling Scheme}
\label{sec:sampling}
Note that the accuracy of the solution of ${\cal P}(N)$ completely depends
on the choice of the $N$ points. The tangent planes to ${\cal S}$ at those
$N$ points create a cover of ${\cal S}$. Throughout this section, we are
going to use the notion of a bounded cover, which we define as
follows.

\begin{definition}
Let ${\cal T}$ be the convex polytope generated by the tangent planes to
${\cal S}$ at the points $\mathbf{x}_1, \ldots, \mathbf{x}_N \in \partial{\cal S}$. ${\cal T}$ is
said to be a bounded cover of ${\cal S}$ if,
$$ \sup_{\mathbf{t} \in {\cal T}} d(\mathbf{t}, {\cal S}) < \infty,$$ where $d(\mathbf{t}, {\cal S}) =
\inf_{\mathbf{x} \in {\cal S}} \norm{ \mathbf{t} - \mathbf{x} }$ and $\norm{\cdot}$ denotes the
Euclidean distance.
\end{definition}

The first result shows that we need a certain minimum number of points
to get a bounded cover.

\begin{lemma}
\label{lem:min_points}
Let ${\cal S}$ be an $s$ dimensional ellipsoid as defined in
\eqref{eq:qc}. Then we need at least $s+1$ points on $\partial{\cal S}$ to
be get a bounded cover.
\end{lemma}

\begin{proof}
Note that since ${\cal S}$ is a compact convex body in ${\mathbb R}^s$, there exists a $s$-dimensional simplex $T$ such that $ {\cal S} \subseteq T$. We can always shrink $T$ such that each edge touches ${\cal S}$ tangentially. Since there are $s+1$ faces, we will get $s+1$ points whose tangent surface creates a bounded cover.

To complete the proof we need to show that we cannot create a bounded cover using only $s$ or fewer points. Consider any set of $s$ points, $\mathbf{x}_1, \mathbf{x}_2, \ldots, \mathbf{x}_s$ on $\partial{\cal S}$. The equation of their tangent planes are,
\begin{align*}
(\mathbf{x} - \mathbf{b})^T \mathbf{B}(\mathbf{x}_j - \mathbf{b}) = \tilde{b}
\end{align*}
for $j = 1, \ldots, s$. Note that we can rewrite this as $\mathbf{A}\mathbf{x} = \boldsymbol{\gamma}$, where each row $\mathbf{a}_j = (\mathbf{x}_j -  \mathbf{b})^T\mathbf{B}$ and $\gamma_j = \tilde{b} + \mathbf{b}^T\mathbf{B}(\mathbf{x}_j - \mathbf{b})$ for $j = 1, \ldots, s$.

Without loss of generality we can assume that $\boldsymbol{\gamma} \in {\cal C}(\mathbf{A})$, otherwise the system of linear equations is inconsistent and this case, it is easy to see that the tangent planes do not create a compact polytope.
Now if $\mathbf{A}$ is not of full rank, there exists a continuum of solutions to $\mathbf{A}\mathbf{x} = \boldsymbol{\gamma}$. Hence, the polytope is not bounded. Finally if $\mathbf{A}$ is invertible, then there exists a unique solution to the system of equations, call it $\mathbf{x}_0$. Since all the planes intersect at a single point, the polytope is divergent. Thus, it is not possible to construct a bounded cover with only $s$ points. Similar proof holds for fewer than $s$ points. This completes the proof of the lemma.
\end{proof}

Lemma \ref{lem:min_points} states that we need a minimum of $s + 1$
points to create the bounding cover ${\cal T}$, but it does not guarantee
the its existence. Moreover, we wish to pick our points in a way such
that the constructed ${\cal T}$ is close as possible to ${\cal S}$, thus having
a better approximation. Formally, we introduce the notion of optimal
bounded cover.

\begin{definition}
${\cal T}^* = {\cal T}(\mathbf{x}_1^*, \ldots, \mathbf{x}_n^*)$ is said to be an optimal bounded cover, if
\begin{align*}
\sup_{\mathbf{t} \in {\cal T}^*} d(\mathbf{t}, {\cal S})  \leq \sup_{\mathbf{t} \in {\cal T}} d(\mathbf{t}, {\cal S})
\end{align*}
for any bounded cover ${\cal T}$ generated by any other $n$-point
sets. Moreover, $\{\mathbf{x}_1^*, \ldots, \mathbf{x}_n^*\}$ are defined to be the
optimal $n$-point set.
\end{definition}

Note that we can think of the optimal $n$-point set as that set of $n$
points which minimizes the maximum distance between ${\cal T}^*$ and
${\cal S}$. It is not hard to see that the optimal $n$-point set on the
unit circle are the $n$th roots of unity, unique upto rotation.

\subsubsection{Reisz Energy and Equidistribution}
\label{sec:reisz}
It has been shown that the $n$th roots of unity minimize the discrete Riesz energy for the unit circle \cite{gotz}. Riesz energy of a point set $A_n = \{\mathbf{x}_1, \ldots, \mathbf{x}_n\}$ is defined as
\begin{align}
\label{eq:riesz}
E_s(A_n) := \sum_{i \neq j = 1}^n \norm{\mathbf{x}_i - \mathbf{x}_j}^{-s}
\end{align}
for positive real parameter $s$. There is a vast literature on Riesz energy and its association with ``good" configuration of points. In fact, we can associate the optimal $n$-point set to the set of $n$ points that minimize the Riesz energy on ${\cal S}$. It is well known that the measures associated to the optimal point set that minimizes the Riesz energy on ${\cal S}$ converges to the normalized surface measure of ${\cal S}$. For sake of compactness we do not go deeper into this. For more details see \cite{grabner, hardin2005minimal} and the references therein.

Thus, to pick the optimal $n$-point set, we try to choose a point set, which has a very good equidistribution property that is lacking in random uniform sampling. One such point set in $[0,1]^s$ is called the $(t, m, s)$-net. To describe this point set we begin with a few definitions. Throughout these definitions $b\ge2$ is an integer
base, $s\ge1$ is an integer dimension and
$\mathbb{Z}_b=\{0,1,\dots,b-1\}$.

\begin{definition}
For $k_j\in\mathbb{N}_0$ and $c_j\in\mathbb{Z}_{b^{k_j}}$ for $j=1,\dots,s$, the set
$$
\prod_{j=1}^s\Bigl[ \frac{c_j}{b^{k_j}},\frac{c_j+1}{b^{k_j}}\Bigr)
$$
is a $b$-adic box of dimension $s$.
\end{definition}

\begin{definition}
For integers $m\ge t\ge0$,
the points $\mathbf{x}_1,\dots, $ $\mathbf{x}_{b^m}  \in[0,1]^s$
are a $(t,m,s)$-net in base $b$
if every $b$-adic box of dimension $s$ with volume $b^{t-m}$
contains precisely $b^t$ of the $\mathbf{x}_i$.
\end{definition}

The nets have good equidistribution (low discrepancy) because
boxes $[0,a]$ can be efficiently approximated by unions of
$b$-adic boxes. Digital nets can attain a discrepancy of
$O((\log(n))^{s-1}/n)$. There is vast literature on easy construction of these point sets. For more details on digital nets we refer to \cite{dick:pill:2010, nied92}.

\subsubsection{Area preserving map to $\partial{\cal S}$}
Now once we have a point set on $[0,1]^s$ we try to map it to $\partial{\cal S}$ using a measure preserving transformation so that the equidistribution property remains. We describe the mapping in two steps. First we map the point set from $[0,1]^s$ to the hyper-sphere $\mathbb{S}^s = \{ \mathbf{x} \in {\mathbb R}^{s+1} : \mathbf{x}^T\mathbf{x} = 1\}$. Then we map it to $\partial{\cal S}$. The mapping from $[0,1]^s$ to $\mathbb{S}^s$ is based on \cite{brauchart2012quasi}.

The cylindrical coordinates of the $d$-sphere, can be written as
\begin{align*}
\mathbf{x} = \mathbf{x}_s &= ( \sqrt{1 - t_s^s} \mathbf{x}_{s-1}, t_s)\\
& \ldots \\
 \mathbf{x}_2 &= (\sqrt{1 - t_2^2} \mathbf{x}_1, t_2)\\
\mathbf{x}_1 &= (\cos \phi, \sin \phi)
\end{align*}
where $0 \leq \phi \leq 2\pi, -1 \leq t_d \leq 1, \mathbf{x}_d \in \mathbb{S}^d$ and $d = (1,\ldots, s)$. Thus, an arbitrary point $\mathbf{x} \in \mathbb{S}^s$ can be represented through angle $\phi$ and heights $t_2, \ldots, t_s$ as,
\begin{align*}
\mathbf{x} = \mathbf{x}(\phi, t_2, \ldots, t_s), \qquad 0 \leq \phi \leq 2\pi, -1 \leq t_2, \ldots, t_s \leq 1.
\end{align*}
We map a point $\mathbf{y} = (y_1, \ldots, y_s) \in [0,1)^s$ to $\mathbf{x} \in \mathbb{S}^s$ using
\begin{align*}
\varphi_1(y_1) = 2\pi y_1, \qquad \varphi_d(y_d) = 1 - 2y_d \;\;\;(d = 2, \ldots, s)
\end{align*}
and cylindrical coordinates
\begin{align*}
\mathbf{x} = \Phi_s(\mathbf{y}) = \mathbf{x}( \varphi_1(y_1), \varphi_2(y_2), \ldots, \varphi_s(y_s)).
\end{align*}

\begin{lemma}
\label{lem:pres_1}
Under the above notation, $\Phi_s : [0,1)^s \rightarrow \mathbb{S}^s$ is an area preserving map.
\end{lemma}

\begin{proof}
See \cite{brauchart2012quasi}.
\end{proof}

\begin{remark}
Instead of using $(t,m,s)$-nets and mapping to $\mathbb{S}^s$, we could have also used spherical $t$-designs, the existence of which was proved in \cite{bondarenko2013optimal}. However, construction of such sets is still a hard problem in large dimensions. We refer to \cite{brauchart2015distributing} for more details.
\end{remark}

We consider the map from $\psi : \mathbb{S}^{s-1} \rightarrow\partial{\cal S}$ defined as follows.
\begin{align}
\label{eq:psi}
\psi(\mathbf{x}) =\sqrt{\tilde{b}}B^{-1/2}\mathbf{x} + b.
\end{align}
The next result shows that this also an area-preserving map, in the sense of normalized surface measures.

\begin{lemma}
\label{lem:pres_2}
Let $\psi$ be a mapping from $\mathbb{S}^{s-1} \rightarrow\partial{\cal S}$ as defined in \eqref{eq:psi}. Then for any set $A \subseteq \partial{\cal S}$,
\begin{align*}
\sigma_s(A) = \lambda_s(\psi^{-1}(A))
\end{align*}
where, $\sigma_s, \lambda_s$ are the normalized surface measure of $\partial{\cal S}$ and $\mathbb{S}^{s-1}$ respectively
\end{lemma}

\begin{proof}
Pick any $A \subseteq \partial{\cal S}$. Then we can write,
\begin{align*}
\psi^{-1}(A) = \left\{ \frac{1}{\sqrt{\tilde{b}}}B^{1/2}(\mathbf{x} - b) : \mathbf{x} \in A \right\}.
\end{align*}
Now since the linear shift does not change the surface area, we  have,
\begin{align*}
\lambda_s(\psi^{-1}(A)) &= \lambda_s\left(\left\{ \frac{1}{\sqrt{\tilde{b}}}B^{1/2}(\mathbf{x} - b) : \mathbf{x} \in A \right\}\right)  \\
&= \lambda_s\left(\left\{ \frac{1}{\sqrt{\tilde{b}}}B^{1/2}\mathbf{x} : \mathbf{x} \in A \right\}\right) = \sigma_s(A),
\end{align*}
where the last equality follows from the definition of normalized surface measures. This completes the proof.
\end{proof}

Following Lemmas \ref{lem:pres_1} and \ref{lem:pres_2} we see that the map $$\psi \circ \Phi_{s -1} : [0,1)^s \rightarrow \partial{\cal S},$$ is a measure preserving map. Using this map  and the $(t,m,s-1)$ net in base $b$, we derive the optimal $b^m$-point set on $\partial{\cal S}$. The
procedure is detailed as Algorithm \ref{algo:simulate}.

\begin{algorithm}
\caption{Point Simulation on $\partial{\cal S}$}\label{algo:simulate}
\begin{algorithmic}[1]
\State Input : $B, b, \tilde{b}$ to specify ${\cal S}$ and $N = k^m$ points
\State Output : $\mathbf{x}_1, \ldots, \mathbf{x}_N \in \partial{\cal S}$
\State Generate $\mathbf{y}_1, \ldots, \mathbf{y}_N$ as a $(t,m,s-1)$-net in base $k$.
\For {$i \in 1, \ldots, N$}
\State $\mathbf{x}_i = \psi \circ \Phi_{s-1}(\mathbf{y}_i)$
\EndFor
\State \Return {$\mathbf{x}_1, \ldots, \mathbf{x}_N$}
\end{algorithmic}
\end{algorithm}

Figure \ref{fig:net2sph2ell} shows how we transform a $(0,7,2)$-net in
base 2 to a sphere and then to an ellipsoid. For more general
geometric constructions we refer to \cite{basu:owen:2015, basu2015scrambled}.

\begin{figure*}
\centering
\includegraphics[width=\linewidth]{net2sph2ell_s.png}
\caption{\label{fig:net2sph2ell} The left panel shows a $(0, 7,
  2)$-net in base $2$ which is mapped to a sphere in 3 dimensions
  (middle panel) and then mapped to the ellipsoid as seen in the right
  panel.}
\end{figure*}

\subsection{Convergence guarantees}
Consider the problem ${\cal P}(N)$ as stated in \eqref{eq:QP_via_lin}. We
shall show that asymptotically as $N \rightarrow \infty$, we get back
the original problem ${\cal P}$ as stated in \eqref{eq:QCQP}. We shall also
prove some finite sample results to give some error bounds on the
solution to ${\cal P}(N)$. We begin with a few notation. Let $\mathbf{x}^*,
\mathbf{x}^*(N)$ denote the solution to ${\cal P}$ and ${\cal P}(N)$ respectively. Let
$f(\cdot)$ denote the objective function.

\begin{theorem}
\label{thm:converge}
Let ${\cal P}$ and ${\cal P}(N)$ be the optimization problems defined in
\eqref{eq:QCQP} and \eqref{eq:QP_via_lin} respectively. Then, ${\cal P}(N)
\rightarrow {\cal P}$ as $N \rightarrow \infty$ in the sense that
$$\lim_{N \rightarrow \infty} \norm{\mathbf{x}^*(N) - \mathbf{x}^*} = 0.$$
\end{theorem}

\begin{proof}
Fix any $N$. Let ${\cal T}_{N}$ denote the optimal bounded cover constructed with $N$ points on $\partial{\cal S}$. Note that to prove the result, it is enough to show that ${\cal T}_{N} \rightarrow {\cal S}$ as $N \rightarrow \infty$. This guarantees that linear constraints of ${\cal P}(N)$ converge to the quadratic constraint of ${\cal P}$, and hence the two problems match. Now since ${\cal S} \subseteq {\cal T}_N$ for all $N$, it is easy to see that ${\cal S} \subseteq \lim_{N \rightarrow \infty} {\cal T}_N.$

To prove the converse, let $t_0 \in \lim_{N \rightarrow \infty} {\cal T}_N$ but $t_0 \not \in {\cal S}$. Thus, $d(t_0, {\cal S}) > 0$. Let $t_1$ denote the projection of $t_0$ onto ${\cal S}$. Thus, $t_0 \neq t_1 \in \partial{\cal S}$. Choose $\epsilon$ to be arbitrarily small and consider any region $A_\epsilon$ on $\partial{\cal S}$ with diameter less than $\epsilon$. Since, we are working with a limiting case, there exists infinitely many points in $A_\epsilon$. Thus there exists a point $t^* \in A_\epsilon$, the tangent plane through which cuts the line joining $t_0$ and $t_1$. Thus, $t_0 \not\in \lim_{N \rightarrow \infty} {\cal T}_N$. Hence, we get a contradiction and the result is proved.
\end{proof}

Note that Theorem \ref{thm:converge} shows that $\lim_{N \rightarrow
  \infty} \norm{\mathbf{x}^*(N) - \mathbf{x}^*} = 0$ and hence, $\lim_{N \rightarrow
  \infty} \left| f(\mathbf{x}^*(N)) - f(\mathbf{x}^*)\right| = 0.$ Now we prove
state some finite sample results.

\begin{theorem}
\label{thm:upp_bound}
Let $g : \mathbb{N} \rightarrow {\mathbb R}$ such that $\lim_{n \rightarrow
  \infty} g(n) = 0.$ Further assume that $\norm{\mathbf{x}^*(N) - \mathbf{x}^*} \leq
C_1 g(N)$ for some constant $C_1 > 0$. Then,
\begin{align*}
\left| f(\mathbf{x}^*(N)) - f(\mathbf{x}^*) \right| \leq C_2 g(N)
\end{align*}
where $C_2 > 0$ is a constant.
\end{theorem}

\begin{proof}
Note that $f(\mathbf{x}) = (\mathbf{x} - {\boldsymbol{a}})^TA(\mathbf{x} - {\boldsymbol{a}})$ and
$\nabla f(\mathbf{x}) = 2A(\mathbf{x} - {\boldsymbol{a}})$. Now, note that we can write
\begin{align*}
f(\mathbf{x}) &= f(\mathbf{x}^*) + \int_{0}^1 \langle \nabla f(\mathbf{x}^* + t(\mathbf{x} - \mathbf{x}^*)), \mathbf{x} -\mathbf{x}^* \rangle dt \\
&= f(\mathbf{x}^*) +  \langle \nabla f(\mathbf{x}^*), \mathbf{x} - \mathbf{x}^* \rangle \\
&\;\;\; + \int_{0}^1 \langle \nabla f(\mathbf{x}^* + t(\mathbf{x} - \mathbf{x}^*)) - \nabla f(\mathbf{x}^*), \mathbf{x} -\mathbf{x}^* \rangle dt\\
&= I_1 + I_2 + I_3 \text{ (say) }.
\end{align*}
Now, we can bound the last term as follows. Observe that using Cauchy-Schwarz inequality,
\begin{align*}
\left|I_3\right| &\leq \int_{0}^1 \left|\langle \nabla f(\mathbf{x}^* + t(\mathbf{x} - \mathbf{x}^*)) - \nabla f(\mathbf{x}^*), \mathbf{x} -\mathbf{x}^* \rangle\right| dt\\
&\leq \int_{0}^1 \norm{ \nabla f(\mathbf{x}^* + t(\mathbf{x} - \mathbf{x}^*)) - \nabla f(\mathbf{x}^*)}\norm{ \mathbf{x} -\mathbf{x}^* } dt\\
&\leq 2  \sigma_{\max}(A) \int_{0}^1 \norm{t(\mathbf{x} -\mathbf{x}^*)}\norm{\mathbf{x} -\mathbf{x}^*} dt \\
&=  \sigma_{\max}(A)\norm{\mathbf{x} -\mathbf{x}^*}^2,
\end{align*}
where $\sigma_{\max}(A)$ denotes the highest singular value of $A$.
Thus, we have
\begin{align}
\label{eq:step1}
f(\mathbf{x}) &= f(\mathbf{x}^*) + \langle \nabla f(\mathbf{x}^*), \mathbf{x} - \mathbf{x}^* \rangle + \tilde{C} \norm{\mathbf{x} -\mathbf{x}^*}^2
\end{align}
where $|\tilde{C}| \leq \sigma_{\max}(A)$. Furthermore,
\begin{align}
\label{eq:step2}
|\langle \nabla f(\mathbf{x}^*), &\mathbf{x}^*(N) - \mathbf{x}^* \rangle| \nonumber\\
&= \left|\langle 2A(\mathbf{x}^* - a), \mathbf{x}^*(N) - \mathbf{x}^* \rangle\right| \nonumber\\
&\leq 2 \sigma_{\max}(A) (\norm{\mathbf{x}^*} + \norm{{\boldsymbol{a}}}) \norm{\mathbf{x}^*(N) - \mathbf{x}^*} \nonumber \\
&\leq 2 C_1 \sigma_{\max}(A) (s +  \norm{{\boldsymbol{a}}}) g(N)
\end{align}
Combining \eqref{eq:step1} and \eqref{eq:step2} we have,
\begin{align*}
\left| f(\mathbf{x}^*(N)) - f(\mathbf{x}^*) \right| \leq C_2 g(N)
\end{align*}
for some constant $C_2 > 0$. Thus, the result follows.
\end{proof}

Note that the function $g$ gives us an idea about how fast $\mathbf{x}^*(N)$
converges $\mathbf{x}^*$. To help, identify the function $g$ we state the
following results.

\begin{lemma}
\label{lem:step1}
If $f(\mathbf{x}^*) = f(\mathbf{x}^*(N))$, then $\mathbf{x}^* = \mathbf{x}^*(N)$. Furthermore, if
$f(\mathbf{x}^*) \geq f(\mathbf{x}^*(N))$, then $\mathbf{x}^* \in \partial{\cal U}$ and
$\mathbf{x}^*(N) \not\in {\cal U}$, where ${\cal U} = {\cal S} \cap \{\mathbf{x} : C\mathbf{x} \leq c \}$
is the feasible set for \eqref{eq:QCQP}
\end{lemma}

\begin{proof}
Let ${\cal V} = {\cal T}_N \cap \{\mathbf{x} : C\mathbf{x} \leq c \}$. It is easy to see that ${\cal U} \subseteq {\cal V}$. Assume $f(\mathbf{x}^*) = f(\mathbf{x}^*(N))$, but $\mathbf{x}^* \neq \mathbf{x}^*(N)$. Note that $\mathbf{x}^* , \mathbf{x}^*(N) \in {\cal V}$. Since ${\cal V}$ is convex, consider a line joining $\mathbf{x}^*$ and $\mathbf{x}^*(N)$. For any point $\lambda_t = t \mathbf{x}^* + (1 - t) \mathbf{x}^*(N)$,
\begin{align*}
f(\lambda_t) \leq t f(\mathbf{x}^*) + (1-t) f(\mathbf{x}^*(N)) = f(\mathbf{x}^*(N)).
\end{align*}
Thus, $f$ is constant on the line joining $\mathbf{x}^*$ and $\mathbf{x}^*(N)$. But, it is known that $f$ is strongly convex since $A$ is positive definite. Thus, there exists only one unique minimum. Thus, we have a contradiction, which proves $\mathbf{x}^* = \mathbf{x}^*(N)$

Now let us assume that $f(\mathbf{x}^*) \geq f(\mathbf{x}^*(N))$. Clearly, $\mathbf{x}^*(N) \not\in {\cal U}$. Suppose $\mathbf{x}^* \in \accentset{\circ}{{\cal U}}$. Let $\tilde{\mathbf{x}} \in \partial{\cal U}$ denote the point on the line joining $\mathbf{x}^*$ and $\mathbf{x}^*(N)$.  Clearly, $\tilde{\mathbf{x}} = t\mathbf{x}^* + (1-t)\mathbf{x}^*(N)$ for some $t > 0$. Thus,
\begin{align*}
f(\tilde{\mathbf{x}}) < t f(\mathbf{x}^*) + (1-t) f(\mathbf{x}^*(N)) \leq f(\mathbf{x}^*)
\end{align*}
But $\mathbf{x}^*$ is the minimizer over ${\cal U}$. Thus, we have a contradiction, which gives $\mathbf{x}^* \in \partial{\cal U}$. This completes the proof.
\end{proof}

\begin{lemma}
\label{lem:step2}
Following the notation of Lemma \ref{lem:step1}, if $\mathbf{x}^*(N) \not\in
{\cal U}$, then $\mathbf{x}^*$ lies on $\partial{\cal U}$ within the conic cap of ${\cal U}$
generated from $\mathbf{x}^*(N)$.
\end{lemma}

\begin{proof}
Since the gradient of $f$ is linear, the result is easy to see from the proof of the second assertion in Lemma \ref{lem:step1}.
\end{proof}

Now we can identify the function $g$ by considering the maximum
distance of the points lying on the intersection of the cone and
${\cal S}$. This is highly dependent on the shape of ${\cal S}$ and on the cover
${\cal T}_N$. Explicit calculation can give us explicit rates of
convergence, which we leave as future work.