


\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.










