Experimental Design Networks:
A Paradigm for Serving Heterogeneous Learners under Networking Constraints ^††thanks: * Y. Liu and Y. Li contributed equally to the paper.

In the standard linear regression setting, a learner observes $n$ samples $(\boldsymbol{x}_{i},y_{i})$, $i=1,\ldots,n$, where $\boldsymbol{x}_{i}\in\mathbb{R}^{d}$ and $y_{i}\in\mathbb{R}$ are the feature vector and label of sample $i$, respectively. Labels are assumed to be linearly related to the features; in particular, there exists a model parameter vector $\boldsymbol{\beta}\in\mathbb{R}^{d}$ such that

and $\epsilon_{i}$ are i.i.d. zero mean normal noise variables with variance $\sigma^{2}$ (i.e., $\epsilon_{i}\sim N(0,\sigma^{2})$).

The learner’s goal is to estimate the model parameter $\boldsymbol{\beta}$ from samples $\{(\boldsymbol{x_{i}},y_{i})\}_{i=1}^{n}$. In Bayesian linear regression, it is additionaly assumed that $\boldsymbol{\beta}$ is sampled from a prior normal distribution with mean $\boldsymbol{\beta}_{0}\in\mathbb{R}^{d}$ and covariance $\boldsymbol{\Sigma}_{0}\in\mathbb{R}^{d\times d}$ (i.e., $\boldsymbol{\beta}\sim N(\boldsymbol{\beta}_{0},\boldsymbol{\Sigma}_{0})$). Under this prior, given dataset $\{(\boldsymbol{x_{i}},y_{i})\}_{i=1}^{n}$, maximum a posteriori (MAP) estimation of $\boldsymbol{\beta}$ amounts to [42]:

where $\boldsymbol{X}=[\boldsymbol{x}_{i}^{\top}]_{i=1}^{n}\in\mathbb{R}^{n\times d}$ is the matrix of features, $\boldsymbol{y}\in\mathbb{R}^{n}$ is the vector of labels, $\sigma^{2}$ is the noise variance, and $\boldsymbol{\beta}_{0},\boldsymbol{\Sigma}_{0}$ are the mean and covariance of the prior, respectively. We note that, in practice, the inherent noise variance $\sigma^{2}$ is often not known, and is typically treated as a regularization parameter and determined via cross-validation.

The quality of this estimator can be characterized by the covariance of the estimation error difference $\hat{\boldsymbol{\beta}}_{\mathtt{MAP}}-\boldsymbol{\beta}$ (see, e.g., Eq. (10.55) in [42]):

The covariance summarizes estimator quality in all directions in $\mathbb{R}^{d}$: given an unseen sample $(\boldsymbol{x},y)\in\mathbb{R}^{d}\times\mathbb{R}$, also obeying (1), the expected prediction error (EPE) is given by

Hence, the eigenvalues of Eq. (3) capture the overall variability of the expected prediction error in different directions.

In experimental design, a learner determines which experiments to conduct to learn the most accurate linear model. Formally, given $p$ possible experiment settings, each described by feature vectors $\boldsymbol{x}_{i}\in\mathbb{R}^{d}$, $i=1,\ldots,p$, the learner selects a total of $n$ experiments to conduct with these feature vectors, possibly with repetitions,¹¹1Note that, due to the presence of noise in labels, repeating the same experiment makes intuitive sense; formally, repetition of an experiment with features $\boldsymbol{x}_{i}$ reduces the EPE (4) in this direction. collects associated labels, and then performs linear regression on these sample pairs. In classic experimental design (see, e.g., [5]), the selection is formulated as an optimization problem minimizing a scalarization of the covariance (3). For example, in D-optimal design, the vector $\boldsymbol{n}=[n_{i}]_{i=1}^{p}\in\mathbb{N}^{p}$ of the number of times each experiment is to be performed is determined by minimizing

or, equivalently, by solving the maximization problem:

In other words, $\boldsymbol{n}\in\mathbb{N}^{p}$ is selected in such a way so that the $\log\det[\mathtt{cov}(\hat{\boldsymbol{\beta}}_{\mathtt{MAP}}-\boldsymbol{\beta})]$ is as small as possible. Intuitively, this amounts to selecting the experiments that minimize the product of the eigenvalues of the covariance;²²2Other commonly encountered scalarizations [5] behave similarly. E.g., E-optimality minimizes the maximum eigenvalue, while A-optimality minimizes the sum of the eigenvalues. alternatively, Prob. (5) also maximizes the mutual information between the labels $\boldsymbol{y}$ (to be collected) and $\hat{\boldsymbol{\beta}}_{\mathtt{MAP}}$; in both interpretations, the selection aims to pick experiments in a way that minimizes the variability of the resulting estimator $\hat{\boldsymbol{\beta}}_{\mathtt{MAP}}$.

We introduce here diminishing-returns submodularity:

The above definition generalizes the submodularity of set functions (whose domain is $\{0,1\}^{p}$) to functions over integer lattice (in the case of DR-submodularity), and continuous functions (in the case of continuous DR-submodularity). Particularly for continuous functions, if $f$ is differentiable, continuous DR-submodularity is equivalent to $\nabla f$ being antitone. Moreover, if $f$ is twice-differentialble, $f$ is continuous DR-submodular if all entries of its Hessian $\nabla^{2}f$ are non-positive. DR-submodularity is directly pertinent to D-optimal design:

For completeness, we provide a proof in Appendix -A. An immediate consequence of this lemma is that polynomial-time approximation algorithms exist to solve Prob. (5) with a $1-1/e$ guarantee (see, e.g., [31, 9]), although Prob. (5) is a classic NP-hard problem [9].

We model the above system as general multi-hop network with a topology represented by a directed acyclic graph (DAG) $\mathcal{G}(\mathcal{V},\mathcal{E})$, where $\mathcal{V}$ is the set of nodes and $\mathcal{E}\subset\mathcal{V}\times\mathcal{V}$ is the set of links. Each link $e=(u,v)\in\mathcal{E}$ can transmit data at a maximum rate (i.e., link capacity) $\mu^{e}\geq 0$. Sources $\mathcal{S}\subset\mathcal{V}$ of the DAG (i.e., nodes with no incoming edges) generate data streams, while learners $\mathcal{L}\subset\mathcal{V}$ reside at DAG sinks (nodes with no outgoing edges). We assume this for simplicity; we discuss how to remove this assumption, and how to generalize our analysis beyond DAGs, in Sec. VI.

We assume that every source generates labeled pairs $(\boldsymbol{x},y)\in\mathbb{R}^{d}\times\mathbb{R}$ of type $t$ according to a Poisson process of rate $\lambda_{\boldsymbol{x},t}^{s}\geq 0.$ Moreover, we assume that generated data follows a linear model (1); that is, for every type $t\in\mathcal{T}$, there exists a $\boldsymbol{\beta}_{t}\in\mathbb{R}^{d}$ s.t. $y=\boldsymbol{x}^{\top}\boldsymbol{\beta}_{t}+\epsilon_{t}$ where $\epsilon_{t}\in\mathbb{R}$ are i.i.d. zero mean normal noise variables with variance $\sigma^{2}_{t}>0$, independent across experiments and sources $s\in\mathcal{S}$.

This is satisfied, if, e.g., the network is a Kelly network [43] of $M/M/1$ queues, $M/M/c$ queues, etc., under FIFO, Last-In First-Out (LIFO), and processor sharing service disciplines, or other queues for which Burke’s theorem holds [42].

We consider a data acquisition time period $T$, at the end of which each learner $\ell\in\mathcal{L}$ estimates $\boldsymbol{\beta}_{t^{\ell}}$ based on the data it has received during this period via MAP estimation. Under Assumption 1, the arrivals of pertinent data pairs at learner $\ell$ form a Poisson process with rate $\lambda_{\boldsymbol{x}}^{\ell}$. Let $n_{\boldsymbol{x}}^{\ell}\in\mathbb{N}$ be the cumulative number of times that a pair $(\boldsymbol{x},y)$ of type $t^{\ell}$ was collected by learner $\ell$ during this period, and $\boldsymbol{n}^{\ell}=[n^{\ell}_{\boldsymbol{x}}]_{\boldsymbol{x}\in\mathcal{X}}$ the vector of arrivals across all experiments. Then,

for all $\boldsymbol{n}=[n_{\boldsymbol{x}}]_{\boldsymbol{x}\in\mathcal{X}}\in\mathbb{N}^{|\mathcal{X}|}$ and $\ell\in\mathcal{L}$. Motivated by standard experimental design (see Sec. III-B), we define the utility at learner $\ell\in\mathcal{L}$ as the following expectation:

where $G^{\ell}(\boldsymbol{n}^{\ell})\equiv G(\boldsymbol{n}^{\ell};\sigma_{t^{\ell}},\boldsymbol{\Sigma}_{\ell})$ and $G$ is given by (5a). We wish to solve the following problem:

Indexing flows by both type $t$ and features $\boldsymbol{x}$ implies that, to implement a solution $\boldsymbol{\lambda}\in\mathcal{D}$, routing decisions at intermediate nodes should be based on both quantities. Problem (16) is non-convex in general.⁵⁵5It is easy to construct instances of objective (16) that are non-concave. For example, when $|\mathcal{L}|=1$, $d=1$, $\mathcal{X}=\{0.1618,0.3116\}$, $\sigma=0.0422$, and $\Sigma_{\ell}=0.2962$, the Hessian matrix is not negative semi-definite. Nevertheless, we construct a polynomial time approximation algorithm in the next section.

Our first main result establishes the continuous DR-submodularity of the objective (16a):

The proof can be found in Section V-C; we establish the positivity of the gradient and non-positivity of the Hessian of $U$. We note that Theorem 1 identifies a new type of continuous relaxation to DR-submodular functions, via Poisson sampling; this is in contrast to the multilinear relaxation [31, 26, 28], which is ubiquitous in the literature and relies on Bernoulli sampling. Finally, though our objective is monotone and continuous DR-submodular, the constraint set $\mathcal{D}$ is not down-closed. Hence, the analysis by Bian et al. [14] does not directly apply, while using projected gradient ascent [28] would only yield a $1/2$ approximation guarantee.

Our algorithm is summarized in Algorithm 1. We follow the Frank-Wolfe variant for monotone DR-submodular function maximization by Bian et al. [14], deviating both in the nature of the constraint set $\mathcal{D}$, but also, most importantly, in the way we estimate the gradients of objective $U$.

We note that we face two challenges preventing us from computing the gradient of $\nabla U$ directly via. Eq. (17): (a) the gradient computation involves an infinite summation over $n\in\mathbb{N}$, and (b) conditional expectations in $\Delta^{\ell}_{\boldsymbol{x}}(\boldsymbol{\lambda}^{\ell},n)$ require further computing $|\mathcal{X}|-1$ infinite sums. Using (17) directly in Algorithm 1 would thus not yield a polynomial-time algorithm. To that end, we replace the gradient $\nabla U(\boldsymbol{\lambda}^{k})$ used in the standard Frank-Wolfe method by an estimator, which we describe next.

The proof can be found in Section V-D. Theorem 2 implies that, through an appropriate (but polynomial) selection of the total number of iterations $K$, the number of terms $n^{\prime}$ and samples $N$, we can obtain a solution $\boldsymbol{\lambda}$ that is within $1-e^{-1}\approx 0.63$ from the optimal. The proof crucially relies on (and exploits) the continuous DR-submodularity of objective $U$, in combination with an analysis of the quality of our gradient estimator, given by Eq. (19).

By the law of total expectation, we have:

Notably, $\frac{\partial U}{\partial\lambda^{\ell}_{\boldsymbol{x}}}=\frac{\partial U^{\ell}(\boldsymbol{\lambda}^{\ell})}{\partial\lambda^{\ell}_{\boldsymbol{x}}}$, for which the following is true:

where the last inequality is true because $G$ is monotone-increasing (Lemma 1).

Next, we compute the second partial derivatives $\frac{\partial^{2}U}{\partial\lambda^{\ell}_{\boldsymbol{x}}\partial\lambda^{\ell^{\prime}}_{\boldsymbol{x}^{\prime}}}$. It is easy to see that for $\textstyle\ell\neq\ell^{\prime}$, we have

For $\ell=\ell^{\prime}$ and $\boldsymbol{x}=\boldsymbol{x}^{\prime}$, it holds that $\frac{\partial^{2}U}{\partial(\lambda^{\ell}_{\boldsymbol{x}})^{2}}=\frac{\partial^{2}U^{\ell}(\boldsymbol{\lambda}^{\ell})}{\partial(\lambda^{\ell}_{\boldsymbol{x}})^{2}}$, where

and the last equality follows from the DR-submodularity of $G$ (Lemma 1).

For $\ell=\ell^{\prime}$ and $\boldsymbol{x}\neq\boldsymbol{x}^{\prime}$, it holds that $\frac{\partial^{2}U}{\partial\lambda^{\ell}_{\boldsymbol{x}}\partial\lambda^{\ell}_{\boldsymbol{x}^{\prime}}}=\frac{\partial^{2}U^{\ell}(\boldsymbol{\lambda}^{\ell})}{\partial\lambda^{\ell}_{\boldsymbol{x}}\partial\lambda^{\ell}_{\boldsymbol{x}^{\prime}}}$,

where the last inequality follows from the DR-submodularity of $G$ (Lemma 1). ∎

Our proof relies on a series of key lemmas; we state them below. Full proofs of all lemmas can be found in the appendix. We begin by associating the approximation guarantee of Algorithm 1 to the quality of gradient estimation $\widehat{\nabla U(\cdot)}$:

The proof, found in Appendix -B, relies on the continuous DR-submodularity of $U$, and follows [14]; we deviate from their proof to handle the additional issue that $\mathcal{D}$ is not downward closed (an assumption in [14]).

Next, we turn our attention to characterizing the quality of our gradient estimator. To that end, use the following subexponential tail bound:

The expression for $h(u)$ implies that the Poisson tail decays slightly faster than a standard exponential random variable (by a logarithmic factor). This lemma allows us to characterize the effect of truncating Eq. (17) in estimation quality. In particular, for $n^{\prime}\geq\lambda^{\ell}_{\boldsymbol{x}}T$, let:

Then, this is guaranteed to be within a constant factor from the true partial derivative:

Next, by estimating $\Delta^{\ell}_{\boldsymbol{x}}(\boldsymbol{\lambda}^{\ell},n)$ via sampling (see (20)), we construct our final estimator given by (19). Putting together Lemma 4 and along with a Chernoff bound [45], to attain a guarantee on sampling, we can bound the quality of our gradient estimator:

Theorem 2 follows by combining Lemmas 2 and 5. In particular, by Lemma 5 and a union bound, we have that (30) is satisfied for all iterations with probability greater than $1-2|\mathcal{X}||\mathcal{L}|\cdot e^{-\delta^{2}N/2T^{2}(n^{\prime}+1)}$. This, combined with Lemma 2, implies that

is satisfied with the same probability. This implies that for any $0<\epsilon_{0},\epsilon_{1}<1$ and $\epsilon_{2}>0$,

with probability $1-\epsilon_{0}$. From Eq. (27), the probability is an increasing function w.r.t. $\lambda^{\ell}_{\boldsymbol{x}}$, and $\lambda_{\mathrm{MAX}}$ is an upper bound for $\lambda^{\ell}_{\boldsymbol{x}}$. Letting

we have

where the last line holds because $u\ln u-u>u$ when $u$ is large enough, e.g., $u\geq e^{2}$. Thus, $n^{\prime}=O(\lambda_{\mathrm{MAX}}T+\ln\frac{1}{\epsilon_{1}})$.
We determine $K$ and $N$ by setting

and

Therefore, $K=O((\frac{\sqrt{2}}{2}|\mathcal{X}||\mathcal{L}|T\lambda_{\mathrm{MAX}}^{2}+2\mathcal{\lambda}_{\mathrm{MAX}})G_{\mathrm{MAX}}/\epsilon_{2})$, and $N=O(T^{2}n^{\prime}K^{2}\ln\frac{|\mathcal{X}||\mathcal{L}|K}{\epsilon_{0}})$. $\hfill\qed$

We consider a finite feature set $\mathcal{X}$ that includes randomly generated feature vectors with $d=100$, and a set $\mathcal{T}$ that of different Bayesian linear regression models with $\boldsymbol{\beta}_{t}$, $t\in\mathcal{T}$. Labels of each type are generated with Gaussian noise, whose variance $\sigma_{t}$ is uniformly at random (u.a.r.) chosen from 0.5 to 1. For each network, we u.a.r. select $|\mathcal{L}|$ learners and $|\mathcal{S}|$ data sources, and remove incoming edges of sources and outgoing edges of learners. Each learner has a target model $\boldsymbol{\beta}_{t^{\ell}}$, $t^{\ell}\in\mathcal{T}$ with a diagonal prior $\boldsymbol{\Sigma}_{t^{\ell}}$ generated as follows. First, we separate features into two classes: well-known and poorly-known. Then, we set the corresponding prior covariance (i.e., the diagonal elements in $\boldsymbol{\Sigma}_{t^{\ell}}$) to low (uniformly from 0 to 0.01) and high (uniformly from 100 to 200) values, for well-known and poorly-known features, respectively. The link capacity $\mu^{e},e=(u,v)\in\mathcal{E}$ is selected u.a.r. from 50 to 100, and source $s$ generates the data $(\boldsymbol{x},y)$ of type $t$ label with rate $\lambda_{\boldsymbol{x},t}^{s}$, uniformly distributed over [2,5].