extfieldsofstudy:
Computer Science
Engineering
Mathematics

text:
Asynchronous Multi Agent Active Search

Active search refers to the problem of efficiently locating targets in an unknown environment by actively making data-collection decisions, and has many applications including detecting gas leaks, radiation sources or human survivors of disasters using aerial and/or ground robots (agents). Existing active search methods are in general only amenable to a single agent, or if they extend to multi agent they require a central control system to coordinate the actions of all agents. However, such control systems are often impractical in robotics applications. In this paper, we propose two distinct active search algorithms called SPATS (Sparse Parallel Asynchronous Thompson Sampling) and LATSI (LAplace Thompson Sampling with Information gain) that allow for multiple agents to independently make data-collection decisions without a central coordinator. Throughout we consider that targets are sparsely located around the environment in keeping with compressive sensing assumptions and its applicability in real world scenarios. Additionally, while most common search algorithms assume that agents can sense the entire environment (e.g. compressive sensing) or sense point-wise (e.g. Bayesian Optimization) at all times, we make a realistic assumption that each agent can only sense a contiguous region of space at a time. We provide simulation results as well as theoretical analysis to demonstrate the efficacy of our proposed algorithms.

Introduction
Active search (also referred to as active sensing, robotic sensing, seek and sample) defines the problem of efficiently locating targets in an unknown environment by interactively collecting data and finds use in applications such as detecting gas leaks, pollution sources or search and rescue missions [1,2,3]. In the field of active learning, most existing active search algorithms are developed for a single agent and are not extendable to multi agent scenarios (see related work in Section 1.2). As an example, [4] uses information greedy approaches to decide on best sensing actions for its agent. If we were to use multiple agents for this info-greedy method, all agents would make the same exact decision at each time step wasting resources of other agents. For other active learning algorithms that are extendable to multi agent scenarios, they usually need a central control system to coordinate the sensing actions of all agents [5,6]. In this paper, however, we are focusing on robotics applications where a central coordination of agents proves infeasible due to unreliable communication and limited onboard computing resources [7]. Such a constraint is very common in multi-robot applications of search and rescue and target detection or localization [8,9,10,11]. To clarify, one could still assume communication channels between agents to create coordination but a central coordinator that expects synchronicity is not feasible as any communication mishap could disrupt the entire process. While communication constraints have been highly discussed in robotics literature, there has not been many work on multi-robot data collection (task decomposition in robotics) without a central coordinator [8,7]. Another consideration of this paper is a realistic assumption on the sensing actions called region sensing initially introduced by [1]. Inspired by search robots, we assume that each agent senses an average value of a contiguous region (block) of the space at each time step. The size of the sensing block models the distance of the agent from the region. We furthermore model noise in our observations in accordance to this distance. Specifically, we assume sensing a larger contiguous region (modelling farther distance from region) inflicts a larger noise value on the resulting observation. Lastly, as an essential part of the real-world applications of active search, we assume targets are sparsely located around the environment.

Contributions
• We propose two novel algorithms, SPATS (Sparse Parallel Asynchronous Thompson Sampling) and LATSI (LAplace Thompson Sampling with Information gain), to actively locate targets in an unknown environment. SPATS is a completely nonparametric algorithm with a probabilistic exploration approach that does not need any prior information about the signal of interest. LATSI leverages the benefits of mutual information with probabilistic exploration in the search space.
• SPATS and LATSI have three main features that collectively distinguish them in robotics applications of active search. 1) They are multi agent methods where agents can asynchronously make independent data-collection decisions without a central coordinator. 2) They are developed given a practical region sensing assumption. 3) They consider sparse signal recovery.
• We demonstrate the efficacy of SPATS and LATSI with an extensive set of simulation results in an asynchronous multi agent setting. We provide theoretical analysis on the benefits of SPATS.
While there have been many sparse recovery algorithms proposed in the literature, to the best of our knowledge there is no algorithm proposed that develops sparse estimators for active learning methods along with multi agent structure and region sensing assumptions. In this paper, we show how sparsity in its nature limits the exploration factor in active learning methods and how a practical region sensing assumption exacerbates this situation. We propose SPATS and LATSI to strategically address such region sensing assumptions to successfully recover sparse signals.

Related Work
A prominent approach to estimating sparse signals is compressive sensing (CS) [12,13]. There has been a large number of work on adaptive CS that enables the ability to make online and adaptive measurements to estimate sparse signals and thus is applicable to active search problems [4,14,15,16,17]. Unfortunately, such adaptive CS methods are sequential and therefore not extendable to multi agent scenarios. Furthermore, CS algorithms in general assume that every measurement matrix can sense the entire environment with arbitrary coefficients which is not a practical assumption for active search problems with region sensing constraints.
Another area of work are multi-armed bandits. [18] and [19] are two multi-armed bandit algorithms that include a sparsity assumption on their hyperparameter. However, the focus of these algorithms are not on estimating the sparse parameter since they are solving for a different reward function. There have been other Bayesian Optimization (BO) and active learning methods proposed for active search. [20] uses BO to develop a spatial mapping of a region whereas we are interested in locating targeted signals. [21] uses BO for localization of single wireless devices but only focuses on point sensing actions. [1,22,23] aim at locating targets by optimizing some notion of Shannon information. Unfortunately, all of the aforementioned active learning algorithms are developed for single agent applications, and except for [1], they mostly lack any realistic assumptions on sensing actions.
Regarding multi-agent algorithms, one well-known approach is to use active learning performed in parallel computation settings. Such algorithms in general require a central control system to manage the actions of all the agents by optimizing a batch of actions at each time step and therefore are not applicable to our problem setting [5,6,24]. Another multi-agent area of work are mobile sensor networks (MSN) [25,26,27] where multiple mobile sensors/agents reconstruct a scalar map of sensory values in an entire area. MSNs typically consider some form of region sensing assumption on their actions, however, one substantial point of difference of these MSN algorithms is the presence of a very constricting sensor network where sensors only follow certain flocking patterns and can only communicate with their close neighbors. Such constrictions are unnecessary in our applications.
In robotics research, methods that deal with active search generally aim at autonomously building topological (identify obstacles and clearways) and/or spatial maps of a region. Our active search (a) Two-dimensional multi agent active search Here, the small vertical lines indicate the start of t'th task. In single agent, tasks start sequentially. In multi agent, task t can start before all previous t − 1 tasks are finished.
problem differs from topological mapping techniques such as SLAM [28,29] and can be most closely related to spatial mapping. For example, [2] identifies strong signals in environments with background information using trajectory planning with confidence intervals; but, unlike our problem setting, their algorithm is developed for a single agent performing point sensing observations. [30] formulates a multi agent game theoretic approach to coordinate unmanned aerial vehicles for cooperative search and surveillance. While they require the actions of neighboring agents for optimal action selection, we find such requirements to be non-essential in our probabilistic approach.
Notation Lowercase and uppercase boldface letters represent column vectors and matrices, respectively. For a matrix A, its transpose is A T . The 1 and 2 -norm of a are denoted by a 1 and a 2 .
The N × N identity matrix is denoted by I N . The Kronecker product is ⊗, and diag(a) is a square matrix with a on the main diagonal. For a set S, |S| denotes number of elements in that set.
2 Problem Formulation: Multi Agent Active Search with Region Sensing Figure 1a illustrates a multi agent active search problem for a two-dimensional environment. Our goal is to efficiently search for targets in an unknown environment by actively taking sensing actions given all the observations thus far. This can be thought of as an active learning problem (referred to as "Design of Experiment" in statistical literature) [31]. In particular, we are interested in recovering the sparse d-dimensional matrix B ∈ R n1×···×n d with minimum measurements. We have no knowledge on the true prior distribution of matrix B other than knowing it is sparse. Defining β ∈ R n as a flattened (vectorized) version of matrix B with n = n 1 ×...×n d , we can write each sensing operation at time step t as: Here, y t is the observation and vector x t ∈ R n is the sensing action at time step t. We call the set of (x t , y t ) the measurement at time step t. Our objective is to estimate the k-sparse vector β (k n) with as few number of measurements T as possible. Here, we are interested in rectangular sensing actions referred to as region sensing [1]. Precisely, in the original d-dimensional space, our sensing action will be a d-dimensional contiguous rectangle (region) with weights w t inside the rectangle and zeros outside. As an example, if d = 1, the sensing action becomes x t = [0, ..., 0, w t , ..., w t , 0, ..., 0] T . This constraint models a robot sensing a region of the search space as illustrated in Figure 1a. Furthermore, we dedicate a fixed amount of power to each sensing action by letting x t 2 = 1.This way sensing a larger region at a farther distance from the region would model larger distortion.
Multi Agent Setting To actively locate targets, at each time step t ≤ T , we choose a sensing action x t given all the available measurements thus far in the set D t−1 . For a single agent this procedure is sequential as in Figure 1b where at time step t the agent uses all previous sequential measurements ., t−1}} to make a decision. In this paper, however, we are interested in an asynchronous parallel approach with multiple agents independently making data-collection decisions as in Figure 1c. Here, asynchronicity means that agents don't wait on results from other agents; instead, an agent starts a new query immediately after its previous data acquisition is completed using all the measurements available thus far; e.g. in Figure 1c, second agent queries t = 6'th action before tasks 4 and 5 are completed using available measurements For easier computations, we can write a compact model of all the available measurements in D t−1 . For example for sequential D t−1 , by defining y = [y 1 , ..., y t−1 ] T , X = [x T 1 , ..., x T t−1 ] T we can write the model in (1) as:

Asynchronous Multi Agent Thompson Sampling
In order to develop an asynchronous multi agent algorithm without central coordination, we borrow ideas from [32] who develop asynchronous parallel Thompson Sampling (TS) for Bayesian Optimization (BO). We propose extending this idea to active search problems. TS is an exploration-exploitation algorithm originally introduced for clinical trials by [33] and later rediscovered for multi-armed bandits [34,35,36]. [32] shows that TS for BO can be developed for an asynchronous parallel setting where each agent makes independent and intelligent decision on the next action given the available measurements. This feature makes TS an ideal solution to our multi agent active search problem. As proven in this reference, asynchronous parallel TS outperforms existing parallel BO methods.
We now review TS for active learning. We start with the single agent setting as introduced in [37]. We are interested in recovering the n-dimensional vector β ∼ p 0 . We actively query actions x t and observe their outcome y t where the likelihood p(y t |x t , β) is known. To query the best action, we maximize a reward function λ(β , D t−1 ). As an example, the reward function can be where β is our belief of the true β, andβ(D t−1 ) is the estimated value of β given all the available measurements D t−1 (e.g. maximum likelihood estimate). We are interested in the myopic policy which selects action x that maximizes the expected reward of Here, the best reward would be the one that has access to the true value of β, i.e. λ + (β, D t−1 , x). Not knowing the true value of β, TS will sample it from the current posterior distribution of β conditioned on the measurements D t−1 , i.e. β ∼ p(β|D t−1 ). Hence, according to TS, the best action x t is one that maximizes the reward assuming that the sample β is the true value of β. The case of multi agent directly follows from [32]. Consider g agents planning on taking T measurements of an environment. Say an agent finishes making an observation and is ready to choose the t'th action. Using measurements available from all agents so far (|D t−1 | ≤ t−1), it will update and sample the posterior (posterior sampling), select its next sensing action that maximizes the reward (design), evaluate its action and share the observations with other agents. Algorithm 1 summarizes this process.

Thompson Sampling with Sparse Prior
In order to develop TS in Algorithm 1 for our active search problem in Section 2, we start by first establishing the prior p(β) and likelihood distribution p(y t |x t , β). As for the prior, our knowledge is limited to the presence of sparsity. Hence, we will assume β has a Laplace distribution with independent entries and a tunable parameter b, i.e. p(β) = 1 (2b) n exp(− β 1 b ). Laplace distribution translates to an 1 -norm regularization term in the cost function which has been shown to introduce sparsity into the estimator [38,39,40]. For the likelihood distribution, the sensing model in (2) gives p(y|X, β) = N (Xβ, σ 2 I t−1 ). Next step is to derive the posterior sampling and design stages of Algorithm 1 using this prior and likelihood. Appendix A provides a detailed derivation of these two stages. We call the resulting algorithm Laplace-TS.

Algorithm 1 Asynchronous Multi Agent Thompson Sampling
Assume: prior β ∼ p 0 and likelihood p(y t |x t , β) For t = 1, ..., T Wait for an agent to finish; for the freed agent: Design Observe y t given action x t ; update & share measurements D t = D t−1 ∪ (x t , y t ) Facing the Failure Mode of TS with Single Agent Unfortunately, Laplace-TS with single agent leads to poor performance that is on par with a point-wise algorithm that exhaustively searches all locations one at a time. We can associate this poor performance with one of the failure modes of TS discussed in Sec. 8.2 of the tutorial by [36]. According to this tutorial, TS faces a dilemma when solving certain kinds of active learning problems. One such scenario are problems that require a careful assessment of information gain. In general, by optimizing the expected reward, TS always restricts its actions to those that have a chance in being optimal which in our case are sparse sensing actions restricted further by the region sensing constraint. However, in active learning problems such as ours, suboptimal actions (i.e. nonsparse sensing actions) can carry additional information regarding the parameter of interest. Appendix B includes simulation results as well as an example to further illustrate the failure mode in active search problems. In the next section, we will modify Laplace-TS and propose two algorithms that can bypass this failure mode. Per our discussion in Section 3.1, introducing sparsity with Laplace prior into TS algorithm limited its ability to explore queries. With this in mind, one might conclude that choosing non-sparse samples in the posterior sampling stage of Algorithm 1 should solve this problem. However, this strategy will still face the failure mode of TS because it is the sparse estimator in the design stage that is limiting the feasible sensing actions. The next logical solution would then be to make both the estimator and posterior sampling procedures non-sparse. Even though with this strategy we will avoid the failure mode of TS, without taking advantage of the prior information about sparsity, the resulting non-sparse TS will be performing no better than exhaustively searching the entire space. To overcome this issue, we propose making an assumption on the prior distribution of both the sampling and estimation procedures that the neighbouring entries of the sparse vector β are temporally correlated, i.e. β is block sparse. Such temporal correlation creates the most compatible results to the region sensing constraint which only approves sensing actions with a single non-zero block of sensors. Furthermore, we expect block sparsity to introduce exploration ability while also keeping sparsity a useful information in the recovery process. In particular, by gradually reducing the length of the blocks from a large value, we gently trade exploration with exploitation capability over time.
In short, borrowing ideas from a block sparse Bayesian framework introduced in [41], we use a block sparse prior p(β) = N (0 n×1 , Σ 0 ), where: Σ 0 = diag ([γ 1 B 1 , ..., γ M B M ]) , with γ m and B m ∈ R L×L (m = 1, ..., M ) as hyperparameters. Here, γ m controls the sparsity of each block as is the case in sparse Bayesian learning methods [42,43], i.e. when γ m = 0, the corresponding block m is zero. Here, L is the length of the blocks that we will gradually reduce in the TS process. To avoid overfitting while estimating these hyperparameters, [44] suggests one matrix B to model all block covariances, namely Σ 0 = diag(γ) ⊗ B, where, γ is the vector containing all elements of γ m for m = 1, ..., M . Appendix C provides a detailed derivation of Algorithm 1 with this prior for the active search problem in Section 2. Algorithm 2 called SPATS summarizes our results in this section.

Algorithm 2 SPATS
Assume: Sensing model (1); sparse signal β; g agents; block length L Set: D 0 = ∅; L = n/g; γ m = 1; B : random highly correlated covariance matrix For t = 1, ..., T Wait for an agent to finish; for the freed agent: Sample β ∼ p(β|D t−1 , γ, B) from (14) in Appendix C Select x t = arg max x λ + (β , D t−1 , x) using (16) in Appendix C Observe y t given action x t ; update & share measurements D t = D t−1 ∪ (x t , y t ) Update hyper parameters γ and B using EM algorithm in (15) in Appendix C if t % g = 0 then L = L/2 SPATS has much lower computational cost than Laplace-TS algorithm since it does not require a Gibbs sampler. Furthermore, unlike Laplace-TS, SPATS is completely nonparametric and does not need to know the sparsity rate or any other prior information about the true signal β.

Theoretical Bounds for a Sparse Model
We now provide theoretical analysis testifying to the benefits of SPATS. SPATS has two aspects that distinguish it from a naïve TS developed for sparse signals. One is using a block sparse prior with varying block length and two is using multiple agents. In what follows, we introduce two theorems to investigate the benefits of each aspect separately. First in Theorem 1, for a sparse model with single agent setting we will compute and compare upper bounds on the expected regret of two TS algorithms with a 1-sparse and a 1-block sparse prior with one nonzero block. The 1-block sparse prior closely imitates SPATS's performance with a region sensing assumption. See proof in Appendix D.
Theorem 1. Consider an active search problem with a 1-sparse true parameter β ∈ R n and reward function R(x, β) = (x T β) 2 for action x ∈ R n chosen from set of actions X that satisfy region sensing in Section 2. Consider two single agent TS algorithms where one assumes a 1-sparse prior and another uses a 1-block sparse prior with varying block length as defined in Algorithm 2. Then, the A simple comparison of (4) and (5) in Theorem 1 shows that using TS with a block sparse prior and varying block length significantly reduces the regret bounds comparing to TS that is using the true 1-sparse prior. Next, we will compute and compare an upper bound on the expected regret of a single-agent and an asynchronous multi-agent TS algorithm. To the best of our knowledge, only theoretical analysis for asynchronous parallel TS has been provided by [32] which is limited to Gaussian Processes. In the following theorem, we provide theoretical guarantees for an asynchronous multi-agent active search problem with a sparse model with proof in Appendix E.
Theorem 2. Consider the active search problem in Theorem 1. Let us propose two TS algorithms with a 1-sparse prior where one is single agent and another uses g agents in an asynchronous parallel setting. Then, the expected regret as defined in Theorem 1 for each algorithms is: A simple analysis of (7) shows that for g n and g T (which is a reasonable assumption), the third term in the bound will be upper bounded by 2g + 1. As a result, the difference in expected regret between single agent and asynchronous multi agent is negligible in terms of number of measurements T . Hence, we can conclude that by dividing the same number of measurements T between g agents, multi agent algorithm achieves same regret g times faster than single agent setting.
Remark 1. Theorem 2 shows that our asynchronous multi agent algorithm performs on par with an optimal multi agent system with a central planner. This result is a consequence of the central planner's regret being bounded by the single agent in terms of number of measurements T [45].

LATSI: Laplace TS with Information gain
In Section 3.1, we discussed how single agent Laplace-TS fails due to a careless assessment of the information gain. To combat this issue, [46] proposes a new reward function which is a combination of expected regret and the mutual information between the optimal action and the next observation. They show that their algorithm called Information Directed Sampling (IDS) can considerably improve the performance of single agent Laplace-TS. Unfortunately, computing the mutual information as introduced in the IDS algorithm is not computationally feasible for our problem. Specifically, IDS proposes sampling techniques to approximate the mutual information between the optimal action and the next observation. However, sampling the optimal action would require computing it for each sample which considering the region sensing assumption is quite expensive.
Another information-theoretic active search algorithm is Region Sensing Index (RSI) by [1] that recognizes region sensing constraints. RSI searches for sparse signals by maximizing the mutual Algorithm 3 LATSI Assume: Sensing model (1); sparse signal β; g agents Set: D 0 =∅; τ : randomly initialized; α, η for tuning; p 0 : 1-sparse uniform distribution on β For t = 1, ..., T Wait for an agent to finish; for the freed agent: average λ + over all x using I from [1] and λ + in (12) Observe y t given action Update τ using EM algorithm in (11) and p(β|D t−1 , τ ) following RSI-A in [1]   information between the next observation and the true parameter β. RSI does not extend well to multi agent settings as without randomness in its reward function, all agents will solve for the same sensing action. However, inspired by the IDS algorithm, we propose combining the reward of Laplace-TS with the mutual information computed in RSI. We call this algorithm LATSI summarized in Algorithm 3 with additional details in Appendix F. We will provide simulation results in the next section, showing that while LATSI improves Laplace-TS, for k > 1 LATSI performs poorly comparing to SPATS due to RSI's poor approximation of mutual information for that case.

Numerical Results
We now compare the performance of our SPATS and LATSI against the information-theoretic approach called RSI proposed in [1]. In this section, we focus on 2-dimensional search spaces (d = 2), where we estimate a k-sparse signal β with length n = 8 × 16 and two sparsity rates of k = 1, 5. Here, β is generated with a randomly uniform sparse vector. We set the noise variance to σ 2 = 1. For LATSI, we set the tuning parameters α, η = 1. Note that neither SPATS nor LATSI are aware of the true uniform sparse prior or sparsity rate k. We then vary the number of measurements T and plot the mean and standard error of the full recovery rate over 50 random trials. The full recovery rate is defined as the rate at which an algorithm correctly recovers the entire vector β over random trials. To further demonstrate the efficacy of SPATS and LATSI, we provide additional experiments for d = 1, larger length n, k = 10 and a sensitivity analysis for LATSI in Appendix G.

Single Agent
In a single agent setting, Figure 2a shows that for k = 1, RSI and LATSI outperform SPATS. The reason is that RSI has a very accurate approximation of mutual information for k = 1 and consequently it is difficult for our SPATS to win over the information-optimal algorithms of RSI and LATSI. On the other hand, for higher sparsity rate of k = 5, SPATS outperforms RSI and LATSI. This is a result of poor approximation of mutual information for k > 1 by RSI. Specifically, for k > 1 RSI recovers the support of β by repeatedly applying RSI assuming k = 1. The authors use this strategy to avoid the large cost of computing mutual information for k > 1. Since our proposed LATSI is a combination of RSI and Laplace-TS, its performance is tied to that of both RSI and SPATS.  To better demonstrate the multi agent performance, we also plot full recovery rate as a function of time which is computed by dividing the number of measurements T by the number of agents g. In each sub-figure we vary the number of available agents between 1, 2, 4 and 8. Furthermore, in all subsequent plots, LATSI-g, RSI-g or SPATS-g indicate the corresponding algorithm with g agents available.
SPATS: We see in Figure 3 that by increasing the number of agents g, SPATS become g times faster as evident in the right two sub-figures. From the left two sub-figures, we can draw a similar conclusion. That is, increasing the number of agents from 1 to 2 to 4 and 8 hardly changes the total number of measurements required for a given recovery rate, i.e. the average number of sensing actions per agent is improved about g times. This result demonstrates that SPATS can efficiently perform active search in an asynchronous parallel fashion.

RSI:
We extend the RSI algorithm of [1] to multi agent setting by allowing each agent to independently choose its sensing action given RSI's acquisition function and utilizing the available measurements from other agents. Looking at Figure 4 for both k = 1 and k = 5, we see a significant deterioration in full recovery rate as a function of T as the number of agents increases. The reason is that without randomness in RSI's reward function, agents that are working at the same time are repeating the same sensing actions. For k = 5, this performance reduction is also obvious as a function of time. However, for k = 1 RSI performs slightly better in time by increasing agents. The reason for this contradicting behavior is that RSI's performance for k = 1 is so close to optimal (binary search) that it reaches recovery rate of 1 before the multi agent system can negatively affect it.
LATSI: Looking at Figure 5, we see that similar to SPATS, LATSI's multi-agent performance improves in time by increasing the number of agents. SPATS vs. LATSI vs. RSI: In Figure 2b, we plot all three algorithms against each other for 4 agents. Here, for k = 1, RSI and LATSI outperform SPATS due to their information-theoretic approach in computing the reward function. For k = 5, SPATS outperforms both RSI and LATSI. This is because SPATS is carefully designed to use randomness from TS in its reward function such that multiplying the number of agents would multiply its recovery rate. Furthermore, LATSI performs significantly better than RSI due to the probabilistic exploration aspect of TS in its reward function.

Conclusions
We have proposed two novel algorithms -SPATS and LATSI which are suitable for the recovery of sparse targets in a multi agent (parallel) asynchronous active search problem with a region sensing constraint. We have discussed the role of sparsity in our design principle and also compared the limitations of a purely information theoretic approach in this setting. An interesting direction for future work would involve constraining the total travel distance of the agents in addition to the number of measurements. Moreover, one could use continuous sensing along a trajectory as an objective rather than only sensing at stopping points.

Appendix: Asynchronous Multi Agent Active Search Abstract
This appendix includes derivations for the proposed Thompson Sampling algorithms Laplace-TS, SPATS and LATSI, discussion and example on failure mode of Laplace-TS, theoretical proofs for Theorem 1 and Theorem 2 and additional numerical results.
Notation Lowercase and uppercase boldface letters represent column vectors and matrices, respectively. For a matrix A, its transpose A T , and the kth row and th column entry is A k, . For a vector a, the kth entry is a k , and the sub-vector containing the ith to jth entries (excluding j) is a[i : j] (Python notation). The 1 and 2 -norm of a are denoted by a 1 and a 2 , and |a| represents its absolute value applied element-wise. The Kronecker product is ⊗, and the trace operator is tr(·). The N × N identity matrix is denoted by I N . diag(a) is a square matrix with a on the main diagonal. ∅ denotes an empty set and for a set S, |S| shows the number of elements in that set. A ∧ B depicts the logical AND between two sets A and B. We use symbols and 1(.) as the symbols for mathematical definition and indicator function, respectively.

A Laplace-TS: Multi Agent Thompson Sampling with Sparse Prior
We will now derive the posterior sampling and design stages of Laplace-TS algorithm introduced in Section 3.1. As detailed in this section, we use a Laplace sparse prior p(β) = 1 (2b) n exp(− β 1 b ) and likelihood distribution p(y|X, β) = N (Xβ, σ 2 I t−1 ) for this problem.
Posterior Sampling In this stage (also referred to as inference stage in active learning algorithms), we regularly update the posterior distribution given the available data sequence D t−1 . Using Bayes rule we have: Unfortunately, computing the normalization factor Z above is intractable for a Bayesian likelihood distribution with a Laplace prior. To compute this posterior, we borrow ideas from [47]. In particular, we substitute the Laplace prior distribution with a zero-mean Gaussian prior per entry p(β i |τ i ) = N (0, τ i ), where its variance τ i has an exponential hyper prior p(τ i ) = η 2 exp(− η 2 τ i ) with τ i ≥ 0. Integrating out τ i to compute p(β i ), we see that this new prior is equivalent to the original Laplace prior per entry β i : Here, √ η = 1/b is the scaling hyperparameter in Laplace distribution and for best performance needs to be tuned given the sparsity rate of the original signal β. If we were able to observe variance τ i , the posterior distribution of β given data sequence D t−1 would become: with, Here, τ is the vector containing all elements of τ i for i = 1, ..., n. Since we cannot observe τ , we use Gibbs sampling [48] to iteratively sample τ and β from their conditional posterior distributions. For the given prior, [49] has computed the conditional distributions to iteratively sample τ and β as follows: Design In this stage, we wish to compute the expected reward λ + (β , D t−1 , x) and optimize it for best sensing action x t . Let us use the reward function λ(β , In order to compute the expected reward in (3), we need to design an estimator β(D t−1 ∪ (x, y)) whose expectation we can compute. While there has been many sparse recovery algorithms proposed in the literature [50], for many of these well-known thresholding or iterative algorithms (e.g. [51,52,53,54,55]) computing the expectation in (3) is not tractable. Alternatively, we propose using maximum a posteriori (MAP) estimate of the posterior p(β|D t−1 , τ ) in (8) where we estimate τ using Expectation-Maximization (EM) proposed in [47] with the E-step and M-step for j = 1, ..., J iterations as follows: Since many elements ofβ (j) and consequently τ (j) are expected to be zero, to avoid inverting diag(τ ) we can rewrite the M-step as follows: Since the posterior is Gaussian, its MAP estimate is the mean of the distribution, i.e.β (J) in (11).
Adding (x, y) to the data sequence D t−1 , our MAP estimate becomes: Now, we can derive the expected reward in (3) as follows: With the posterior sampling and design stage developed, we can now run Algorithm 1 for the active search problem in Section 2. Let us call this algorithm Laplace-TS. In the next section, we will evaluate the performance of this algorithm.
B Failing Performance of Laplace-TS in Single Agent Setting We can associate this poor performance with one of the failure modes of TS discussed in Sec. 8.2 of the tutorial by [36]. According to this tutorial, TS faces a dilemma when solving certain kinds of active learning problems. One such scenario are problems that require a careful assessment of information gain. In general, by optimizing the expected reward, TS always restricts its actions to those that have a chance in being optimal. However, in certain active learning problems, suboptimal actions can carry additional information regarding the parameter of interest.
We provide the following example to better understand the situation in our problem formulation. Let us assume there is no noise in the sensing model (1) ( t = 0) and that there is only one non-zero element (k = 1). Under such conditions, the problem amounts to finding the location of the non-zero element in β which we callĩ. Thus, for every action x t , the observation y t = x T t β is non-zero ifĩ is in the support of x t and is zero, otherwise. Clearly, a binary search can locateĩ in log(n) steps. TS, however, in each time step estimates the sample β withβ and picks the sensing action x t =β. Or, in case of perfect estimation, we will have x t =β = β . So, unless β is the true β, observation y t would be zero and TS will only manage to eliminate support of β from the list of possibly correct supports. Therefore, TS ends up eliminating wrong supports one location at a time which explains its on par performance to an exhaustive point sensor. Similarly, for k > 1 TS is always limited to picking sensing actions x t that are in the support of the sparse estimates of the samples. Adding the region sensing constraint will only aggravate this problem by further shrinking this support set. In the next section, we will modify Laplace-TS and propose two algorithms that can bypass this failure mode.

C SPATS: Thompson Sampling with Block Sparse Prior
In what follows, we will derive the posterior sampling and design stages for SPATS algorithm proposed in Section 4.1. SPATS is an asynchronous multi agent TS algorithm using a block sparse prior for the problem formulation in Section 2.
As discussed in Section 4.1, we use the following prior and likelihood distributions. The likelihood function is p(y|X, β) = N (Xβ, σ 2 I t−1 ) from (2). For prior distribution, as discussed in Section 4.1 we use a block sparse prior p(β) = N (0 n×1 , Σ 0 ), where: with γ ∈ R M and B ∈ R L×L (M = n/L) as hyperparameters.
Posterior Sampling Similar to (8) and (9) in computing posterior distribution for Laplace-TS, the posterior with a Gaussian likelihood and Gaussian prior becomes: with, Using Expectation-Maximization proposed in [41], we can estimate the hyperparameters γ and B for j = 1, ..., J iterations as follows.

D Proof for Theorem 1
Since this theorem focuses on single agent settings, both proposed TS algorithms in the theorem have the theoretical guarantees as well as the limits of TS bounds as summarized in [36,Ch. 8]. As discussed in this chapter, there are different theoretical approaches proposed in literature for regret bounds on TS. For this proof, we use regret bounds computed via information theory proposed in [56]. According to this paper, the resulting bounds better describe the benefits of prior distribution which helps point to failure modes of TS and potential solutions. Defining π t as the distribution of actions for an online optimization policy at time t, [56] shows that for this policy, the expected regret defined in Theorem 1 is bounded by: where, H(x ) is the entropy for optimal action x andΨ t (π) is the average expected information ratio defined as follows.
where, ∆ t (π t ) and I t (π t ) are the single-period expected regret and information gain over policy π t at time step t defined below. For the single-period expected regret we have: where, R t (x, β) is the reward function for action x at time t. And, for information gain, we have: where, I t (x ; y|x) is the mutual information between optimal action x and observation y = (x T β) 2 at time t.
To compute the regret bound in (17), we need to first focus on computing single-period expected regret and information gain in (19) and (20). Computing these terms for a general problem formulation can be difficult. However, for an online optimization problem with a finite prior distribution and finite action set, [46, Algorithm 1] has provided the necessary formulation to compute these terms. Since our problem formulation falls under the category of this framework, we can take advantage of their formulations as follows. For single-period expected regret, we have: and, for the mutual information we have: where, with B(x) = {β : x = arg maxx E y|x,β [R(x, β)]} as the set containing all vectors β for which action x is optimal.
Next, we will compute this regret bound for our TS algorithms with 1-sparse and block sparse prior introduced in Theorem 1. Our emphasis here is on understanding the benefits of a block sparse prior distribution. Hence, to simplify the discussion we assume both algorithms are aware of the true 1-sparse prior for β.

1-sparse Prior
Let us here establish two parameters i and j that we will use throughout the proof. 1) Since the true parameter β is 1-sparse, we assume the location of the non-zero element of true β is i . Given the uniform 1-sparse prior on true β, i can be either of indices 1 through n with probability 1/n. 2) Since our TS algorithm here assumes a 1-sparse prior distribution, the only feasible sensing actions in X are 1-sparse with one nonzero element at location j.
With two parameters i and j, we can now compute the single-period expected regret in (21) for t = 1 as: Hence, the single-period expected regret ∆ 1 (π 1 ) over our policy becomes: Here, (a) follows from the fact that our TS algorithm at t = 1 will pick all sensing actions with the nonzero element j = 1 through j = n equally likely. Next, we need to compute single-period expected regret in (21) for t = 2. For t = 2, there are two policies for our TS algorithm depending on the action x 1 chosen in previous time step. If x 1 = x , then our TS algorithm will be exploiting the same optimal action at time t = 2 and set x 2 = x for which ∆ 2 (x) = ∆ 2 (π 2 ) = 0. The probability of finding x at t = 1 is 1/n. However, if x 1 = x , then our TS algorithm will be exploring the n − 1 feasible actions that are not observed. Let us assume j 1 shows the location of nonzero element in x 1 , then: Hence, for both policies the single-period expected regret ∆ 2 (π 2 ) becomes: Similarly, we can conclude that for any t > 2, we will have: We will next compute the information gain in (22) for time step t. Recall that our TS algorithm have two policies. One policy π t is to keep exploiting the optimal action if it has been discovered in the previous t−1 measurements D t−1 . Since ∆ t (π t ) = 0 for this policy, we have Ψ t (π t ) = ∆t(πt) 2 It(πt) = 0 and hence there is no need to compute the information term I t (π t ) for this scenario. The second policy of our TS method is used if optimal action x has not been found in the previous t − 1 measurements. For this scenario, at time t the TS algorithm picks the sensing action x t out of n − t + 1 measurements that have not been observed yet. For this policy, we can compute the probabilities in (23) as follows. First, for p(x ), we have: Second, for p(y|x) we have: where, {j 1 , ..., j t−1 } denotes the location of the nonzero elements in the previous t − 1 measurements D t−1 . Third, for p(x , y|x) we get: otherwise.
Using these three probability densities, the information gain in (22) becomes: resulting in: Now, using (24) and (25) , we can compute the average expected information ratio in (18) as: which together with an upperbound H(x ) < log (|X |) and (17) gives the following regret bound: Lastly, the regret bound we computed above is only applicable when T ≤ n. Since there are n 1-sparse actions available for a vector β with length n, after T = n actions the algorithm has surely found the optimal action. As a result, ∆ t (π t ) = 0 for t ≥ n which gives us the regret bound: We noticed a minor mistake in expected regret term (4) in Theorem 1 after paper submission. Equation (26) above describes the corrected term. This does not affect any of the conclusions in the main paper.

block sparse prior
We will again start with establishing two parameters i and m t that we will use throughout the proof. 1) We assume the location of the non-zero element of true β is i . Given the uniform 1-sparse prior on true β, i can be either of indices 1 through n with probability 1/n. 2) Since our TS algorithm here assumes a 1-block sparse prior distribution, the feasible sensing actions in With two parameters i and m t , we can now compute the single-period expected regret ∆ 1 (π 1 ) in (19) for t = 1 as: here, (b) follows from the fact that our TS algorithm at t = 1 will pick all sensing actions with the nonzero blocks m 1 = 1 through m 1 = M 1 equally likely, and (c) follows from (21). For t > 2, there are two policies depending on whether i has been discovered to belong to any of the elements of block m 1 or not. If i does not belong to block m 1 , then: Or if i does belong to block m 1 , similarly we ∆ 2 (π 2 ) = 1 − 1 L1 . For block length L 1 = n/2, we will have the regret ∆ 2 (π 2 ) = 1 − 1 n/2 for both policies which is the smallest regret value. Similarly, for any t > 2, we have: as long as L t = n 2 t . We will next compute the information gain in (22) for time step t = 1. For that, we will first compute the three probability densities in (23). For p(x ), we have: if x is 1-sparse 0, otherwise which follows the fact that the optimal action is 1-sparse. Second, for p(y|x) we have: Third, for p(x , y|x) we get: Using these three probability densities, the information gain in (22) becomes: The information term above is maximized when L 1 = n/2 which gives I 1 (x ; y|x) = log (2).
Next, let us discuss information gain for t = 2. Recall that for t = 2, there are two policies depending on whether i has been discovered to belong to elements of block m 1 or not. Similar computations to that of t = 1 shows that for the two policies we have: By setting L 1 = n/2 and L 2 = n/4, both policies will be maximized to give I 2 (x ; y|x) = log (2). Similarly, for any t > 2, the information gain I t (x ; y|x) is maximized with: as long as L t = n 2 t . The fact that the varying block length L t = n 2 t maximizes information gain I t (x ; y|x) and minimizes regret term ∆ t (π t ) shows that this varying block length is the optimal approach for our TS algorithm finding a 1-sparse signal of interest. This result in essence describes our strategy for varying block length L t = Lt−1 2 in SPATS algorithm. Finally, using (27) and (28) with the upperbound H(x ) < log (|X |), we can compute the expected regret in (17) as: Lastly, the regret bound we computed above is only applicable when T ≤ log 2 (n). Since the algorithm has a varying block length of L t = n 2 t , after T = log 2 (n) the block length is reduced to 1 which certainly includes nonzero element of β. As a result, after T = log 2 (n) actions the algorithm has surely found the optimal action. This result means that ∆ t (π t ) = 0 for t > log 2 (n) which gives us the regret bound: Remark 3. We noticed a minor mistake in expected regret term (5) in Theorem 1 after paper submission. Equation (29) above describes the corrected term. This does not affect any of the conclusions in the main paper.

E Proof for Theorem 2
Unlike Theorem 1, we use a direct approach to compute expected regret in this proof. Similar to Appendix D, to simplify the discussion, we assume both algorithms are aware of the prior distribution on true parameter β. We will first compute the expected regret for the single agent algorithm.

Single Agent
Since the true parameter β is 1-sparse, we assume the location of the non-zero element is i . Recall the equation for expected regret: where, (a) follows from assuming i is a fixed unknown location without loss of generality. To compute the first term R(x , β), note that the optimal action x is the 1-sparse vector with i as the nonzero element, i.e. R(x , β) = 1.
To compute the second term E xt [R(x t , β)], remember that x t for TS is selected by maximizing the reward function given the available measurements D t−1 . For a single agent setting, as depicted in Figure 1b, we have D t−1 = {(x j , y j )|j = 1, .., t−1}. Since our algorithm has a 1-sparse prior, only feasible sensing actions are 1-sparse. If x ∈ D t−1 , then TS will be exploiting the same optimal action at time t and set On the other hand, if x / ∈ D t−1 , then TS will select action x t out of n − t + 1 actions that are not observed as part of D t−1 . Hence, the reward R(x t , β) is 1 with probability 1 n−t+1 and 0 otherwise, Putting the information for the first and second term in (30), we have: where, (b) follows from p(x ∈ D t−1 ) = 1 − p(x / ∈ D t−1 ) = t−1 n given that there are t − 1 distinct measurements out of n distinct possible actions until x is found.
Lastly, the regret term we computed above is only applicable when T ≤ n. Since there are n 1-sparse actions available for a vector β with length n, after T = n actions the algorithm has surely found the optimal action. As a result, R(x , β) − R(x t , β) = 1 − 1 = 0 for t > n which gives us: T n (T n + 1) 2n , T n = min{T, n}.

Multi Agent
Now, we will compute the expected regret for TS with asynchronous multi agent setting. Similar to single agent, the equation for expected regret follows from (30) with i as the location of the nonzero element of β. Again, for the first term we have R(x , β) = 1 where the optimal action x is the 1-sparse vector with i as the nonzero element. For the second term E xt [R(x t , β)], we know that if x ∈ D t−1 , then TS will be exploiting the same optimal action at time t and set x t = x for which On the other hand, if x / ∈ D t−1 , then TS will select action x t out of actions that are not observed as part of D t−1 . Precisely: is not as straight forward as the single agent setting. In a single agent setting, |D t−1 | = t − 1 and all t − 1 actions are distinct until x is found. In an asynchronous multi agent setting, |D t−1 | < t − 1 and some of the actions in D t−1 might be the same. Specifically, as illustrated in Figure 1c an action t starts before all previous t − 1 actions are completed and as a result there is a chance that action t will equal the previous actions that were not included in D t−1 . While we cannot precisely compute p(x ∈ D t−1 ) in this asynchronous multi agent setting, in the following lemma we compute a lower bound for it with the proof provided in Appendix E.1.

Lemma 1.
Consider the active search problem in Theorem 2. For an asynchronous multi agent TS algorithm, p(x ∈ D t−1 ) is bounded by: Using the bound in Lemma 1 along with a naïve bound E xt [R(x t , β)|x / ∈ D t−1 ] ≥ 0, we can upper bound the expected regret in (31) as follows: Lastly, the regret term we computed above is only applicable when T ≤ n + g. If T > n + g, then |D t−1 | = n which includes all n 1-sparse actions available for a vector β with length n. Hence, after T = n + g actions, the algorithm has surely found the optimal action which results in R(x , β) − R(x t , β) = 1 − 1 = 0, and consequently: Proof. Since we need to upper bound the expected regret in (30), we will lower bound p(x ∈ D t−1 ).
To compute this bound, we need to better understand D t−1 . For that, we will look at the example illustrated in Figure 1c. For this example, D 0 through D 11 have the following measurement sets: Our first observation is that |D t−1 | = t − g which applies to any asynchronous example [32]. Next, let us compute p(x ∈ D t−1 ) as an example for D 6 : Here, (c) follows from the fact that on picking x 3 and x 2 , they depend on empty sets of measurements D 2 and D 1 . And, on choosing sensing action x 5 , the algorithm has access to measurements Consequently, in (d) we see that x 1 ,x 3 and x 2 have n possible actions to pick, while x 5 will pick out of n − 2 actions that had excluded the ones in D 4 .
As evident in (32), computing p(x ∈ D t−1 ) depends on the specifics of every example and is different for all D t−1 . However, this probability can be lower-bounded for any asynchronous setting as follows. We know that D t−1 has a cardinality of t − g. As a result, each asynchronous example translates to a set D t−1 with permutation of t − g measurements picked out of all t − 1 possible measurements {(x j , y j )|{j = 1, ..., t − 1}}. Out of all possible asynchronous examples, .., t − g}} is the worst case scenario in terms of finding the optimal action. In other words, for this measurement set, the algorithm has used the least amount of information to decide on each of the actions x 1 through x t . Hence, the asynchronous example with D W t−1 has the lowest probability to pick the optimal sensing action x and can be used to lower-bound p(x ∈ D t−1 ): Now, all we have to do is compute p(x ∈ D W t−1 ). For t < g, we have: For t ≥ g, we have:

F Additional Details on LATSI
In this section, we provide additional details on deriving LATSI algorithm we proposed in Section 4.2.
First, let us do a short review on IDS and RSI algorithms: Review of IDS (Information Directed Sampling) proposed by [46] There are certain online optimization problems such as our active search problem where traditional multi-armed bandit algorithms such as TS and Upper Confidence Bound (UCB) fail due to a careless assessment of the information gain (see Appendix B). IDS is an online optimization algorithm proposed by [46] that addresses this problem by introducing a novel reward function that balances between expected single-period regret and a measure of information gain. Specifically, defining π IDS t as the action sampling distribution of IDS at time t, at each time step t, IDS is computed by: where they call Ψ t (π) the information ratio of the action sampling distribution π, and ∆ t (π) and I t (π) are the single-period expected regret and information gain at time t defined below. For the single-period expected regret, we have: where, R t (x) is the reward function for action x at time t. And, for information gain, we have: where, I t (x ; y|x) is the mutual information between optimal action x and observations y at time t if sensing action x is chosen.
Once IDS computes the action sampling distribution at time t, it will choose an action by randomly sampling this distribution.

Review of RSI (Region Sensing Index) proposed by [1]
RSI is a single agent active search algorithm designed to locate sparse signals by actively making data-collection decisions. Similar to our problem formulation in Section 2, RSI makes a practical assumption that at each time step the agent senses a contiguous region of the space. To decide on their next action, RSI at each time step chooses the sensing action x t that maximizes the mutual information between the next observation y t and the signal of interest β, i.e.
Unfortunately, computing the mutual information I (β; y|x, D t−1 ) has high complexity for sparsity rates of k > 1. In order to reduce this complexity for k > 1 RSI recovers the support of β by repeatedly applying RSI assuming k = 1.

Deriving LATSI
Recall how IDS algorithm solves for the failure mode of TS by using a novel reward function that is the ratio of the expected single-period regret and information gain. Inspired by this algorithm, we propose revising the expected reward λ + (β , D t−1 , x) in (3) for Laplace-TS by adding a measure of information gain to it. For the information gain, we use the mutual information I (β; y|x, D t−1 ) computed in RSI algorithm. Thus, we can say that LATSI is the Laplace-TS algorithm where the design stage has been replaced by: where, Here, we normalize and add a tuning parameter α to best control the importance of each term on the overall reward R + .

G Additional Numerical Results
We now provide additional simulation results to support our analysis in Section 5.

1-dimensional Search Space
In this section we provide additional results for a 1-dimensional search space (d = 1), where we estimate a k-sparse signal β of length n = 128 with two sparsity rates of k = 1, 5. Here, we use the same set of parameters as those outlined in Section 5. Figure 7, 8 and 9 show the performance of SPATS, RSI and LATSI respectively. Our observations overall are similar to those in Section 5.2. Concretely, we observe that both SPATS and LATSI show an improvement in performance with time as more agents become available. The efficiency of this asynchronous search by SPATS and LATSI becomes more pronounced with higher sparsity rate k where there are more targets in the search space. On the other hand, multi agent RSI's performance worsens as we increase the number of agents (especially with larger sparsity rate k). This observation is similar to that in Section 5.2 and is explained by the lack of randomness in the information theoretic reward function of RSI.

2-dimensional Search Space
In this section we provide additional results for a 2-dimensional search space (d = 2), where we estimate a k-sparse signal β of length n = 16 × 16 and two sparsity rates of k = 1, 5. The same set of parameters outlined in Section 5 are followed. Figure 10, 11 and 12 show the performance of SPATS, RSI and LATSI respectively, in the multi agent setting. Our observations overall are similar to those in Section 5.2. The probabilistic decision-making nature of SPATS ensures that its performance multiplies by its number of agents. RSI relies on information theoretic decision making (no randomness). Hence, this lack of randomness in its reward function together with the poor approximation of mutual information for k > 1, results in a poor performance generally worsening with more agents. The multi agent performance of LATSI shows similar trends as SPATS given its partly probabilistic reward function. However, LATSI is not as efficient as SPATS given its poor approximation of mutual information on the info-greedy share of the reward function. These results further support the analysis in Section 5.

1-dimensional vs 2-dimensional Search Space
In Figure 13, we compare the performance of SPATS and LATSI for 1-dimensional (1d) and 2dimensional (2d) search spaces with the same length n = 128. We compare a 1d search space of length n = 128 with a 2d search space of length n = 8 × 16 (i.e. flattened length n = 128) for two sparsity rates k = 1, 5. The same set of parameters outlined in Section 5 are followed. As evident in this figure, while 1-dimensional SPATS and LATSI have a similar performance for k = 1, with larger k they perform better in a 2d space rather than in 1d. The reason for this behavior is our region sensing assumption. Specifically, region sensing in the 2-dimensional grid expands the available action space and allows for a larger feasible action set compared to the 1-dimensional grid. Mainly, the larger action set helps our algorithms to choose better sensing actions in the 2-dimensional space, particularly when multiple targets are present.

Sparsity Rate
In this section, we provide an additional set of results for sparsity rate of k = 10 in a 1-dimensional search space of length n = 128. Figure 14 illustrates the full recovery rate of SPATS, LATSI and RSI for one and four agents, respectively. Overall, we see that increasing the number of targets from k = 5 to k = 10 reduces the performance of all algorithms. Both figures confirm our previous observations that for sparsity rates of k > 1, SPATS algorithm outperforms LATSI and RSI. In the case of single agent, the poor performance of LATSI and RSI can be traced back to poor approximation of mutual information for k > 1 (see Section 5.1). In the four-agent scenario, RSI and LATSI are significantly worse than SPATS due to lack of randomness in information-greedy approaches (see Section 5.2 for details).

Sensitivity Analysis for LATSI
Recall tuning parameter α as the scaling factor in (33) which is used to combine Laplace-TS and RSI's reward functions in LATSI's reward function R + . As evident by (33), increasing the value of  Figure 15: Effect of α in full recovery performance of LATSI with multiple agents and multiple targets α in R + increases the weight of the reward computed by Laplace-TS. We experiment with different values of α ∈ {0, 0.1, 1, 10} to determine it's effect on the full recovery rate performance. Since it is our aim to exploit parallelization, we perform the simulations with 2 agents and k = 5 in a 1-dimensional search space of length n = 128. Plotting the results in Figure 15, we observe that the algorithm's performance is robust for a wide range of α. Based on this observation, we chose α = 1 for all experiments with LATSI.


extfieldsofstudy:
Materials Science
Computer Science

text:
GramGAN: Deep 3D Texture Synthesis From 2D Exemplars

We present a novel texture synthesis framework, enabling the generation of infinite, high-quality 3D textures given a 2D exemplar image. Inspired by recent advances in natural texture synthesis, we train deep neural models to generate textures by non-linearly combining learned noise frequencies. To achieve a highly realistic output conditioned on an exemplar patch, we propose a novel loss function that combines ideas from both style transfer and generative adversarial networks. In particular, we train the synthesis network to match the Gram matrices of deep features from a discriminator network. In addition, we propose two architectural concepts and an extrapolation strategy that significantly improve generalization performance. In particular, we inject both model input and condition into hidden network layers by learning to scale and bias hidden activations. Quantitative and qualitative evaluations on a diverse set of exemplars motivate our design decisions and show that our system performs superior to previous state of the art. Finally, we conduct a user study that confirms the benefits of our framework.

Introduction
Texture synthesis is an important field of research with a multitude of applications in both computer vision and graphics. While vision applications mainly focus on texture analysis and understanding, in computer graphics textures are most important in the field of image manipulation and photorealistic rendering. In this context, applications require synthesis techniques that can generate highly realistic textures. Due to the problem being highly complex, today's applications largely rely on techniques that generate textures based on real photos, which is a time-consuming process that requires expert knowledge. In this work, we present a novel texture synthesis framework that can generate high-quality 3D textures given an exemplar image as input (See Figure 1).
Techniques for texture synthesis can be classified into photogrammetric and procedural techniques. Photogrammetry relies on methods that capture the color appearance of physical objects to obtain (typically several) photos from the target texture. Next, texture mapping algorithms are employed to supply a digital 3D model, usually a polygon mesh, with texture coordinates that define a mapping from 3D surface points to 2D texture pixels. Such photogrammetric approaches yield high-quality textures but rely on time-consuming procedures that require expert knowledge. Moreover, the required projection from 3D to 2D introduces seam and distortion artifacts, for example when texturing spherical objects. Nevertheless, this is still a widely used technique, despite its shortcomings.
An entirely different class of techniques is referred to as procedural texturing, where textures are generated by an algorithm instead of relying on acquired data. Ideally the algorithm takes a 3D position as input and computes a texture color value as output, omitting the error-prone projection step entirely. Pioneering work in procedural texturing is the Perlin noise texture model [1], which dates back to 1985. The core idea is to linearly combine noise frequencies, or octaves, to generate different textures. An artist manually adjusts octave weights to obtain a specific look. While this early Figure 1: Results produce by our framework. Given an exemplar texture image (left), our model can faithfully resynthesize it infinitively (middle), suitable to texture 3D models (right). approach conveniently avoids both data acquisition and texture mapping, its applicability to natural textures is rather limited mainly due to the predefined, fixed noise frequencies. Moreover, a user has to manually tweak the parameters in order to obtain a desired look, which is a tedious process.
Recently, several techniques that leverage deep learning have been proposed to train generative models for procedural texture synthesis. Most of these leverage convolutional neural networks (CNNs) to generate texture pixel values as a function of spatially neighboring data. Recent examples in this direction either use generative adversarial networks (GANs) [2] or a style loss [3] to optimize their models. While these CNN-based approaches yield high-quality textures, in practice their application is limited to either 2D textures or very low resolution 3D textures due to the high memory footprint of 3D CNNs. Moreover, an unsolved problem with these techniques is the lack of diversity or rather stochasticity, i.e. such models struggle producing diverse results that represent the same texture. Despite attempts to mitigate this problem by explicitly minimizing a diversity term [4], true diversity has not been demonstrated by this approach.
Our work is inspired by to the model proposed by Henzler et al. [5], which features both memory efficiency and diversity. The former is achieved by formulating the texture synthesis problem as a point operation parametrized using a multilayer perceptron (MLP), while the latter is achieved by defining the model to be a function of noise frequencies, similarly to Perlin. However, this flexibility comes at the cost of lower image quality and a lack of similarity to the exemplar texture.
We present a point operation model that learns to generate textures as a function of noise frequencies, but we propose a novel loss function and architectural framework, which together improve image quality significantly. Moreover, our proposed network architecture inherently supports an extrapolation technique that improves generalization to unseen textures, that can be substantially different from the training set exemplars. In particular, we make the following contributions: (1) a deep MLP architecture that effectively combines learned noise frequencies; (2) a novel GAN loss inspired by style transfer to achieve both realistic image quality and high similarity to the input exemplar; and (3) a novel extrapolation strategy that significantly improves the results on unseen exemplars.

Related work
Texture synthesis techniques have a long history in computer vision and graphics. In 1985, Perlin [1] proposed a model to generate stochastic textures by linearly combining frequencies of noise in 2D or 3D. By using different combination weights, various visual appearances can be created. However, due to fixed, isotropic frequencies, the model expressiveness is limited and the space of generated textures is restricted to isotropic textural features, inhibiting anisotropic, e.g. elongated structures. Another major limitation of this model in practice is that the linear weights must be tediously tweaked to obtain an output that matches a desired texture. To overcome this limitation, a desired property of any texture synthesis method is the ability to generate textures that automatically match an exemplar texture at hand. Early works that tackle this problem are patch-based techniques [6][7][8][9][10][11], where texture exemplars are either padded, inpainted, or reproduced from scratch by considering local neighborhoods in the exemplar image. While these methods can faithfully reproduce a target texture, they are limited to 2D and are computationally expensive. With recent advances in deep learning, a multitude of learned texture synthesis models have been proposed in the vision community. Most such work either optimize a style loss [3], a GAN loss [2], or a hybrid of the two. As our work is closely related to these approaches, we give an overview of the respective methods in the following.
Texture synthesis using style loss In their groundbreaking work, Gatys et al. [3] demonstrated that statistical correlations of activations extracted from a pre-trained VGG-19 network [12] encode textural features. In particular, it has been shown that matching the Gram matrices of these deep activations successfully guides the synthesis network to generate samples similar to the target texture. Building upon these insights, several methods have been proposed to improve computation time [13,14] and diversity [15] of texture generation. All such models are parametrized using CNNs and therefore they do not scale well to 3D textures at acceptable texture resolutions. Moreover, these approaches all work on a single exemplar texture, which implies the need for training an individual model for each target texture. Finally, obtaining true diversity in generated samples is still an unsolved problem with these methods. A big leap forward has been achieved recently by Henzler et al. [5]. They propose a point operation model parametrized as an MLP instead of using CNNs, which enables efficient synthesis of high-resolution 3D textures. Moreover, true diversity is demonstrated thanks to combining noise frequencies directly, similarly to Perlin noise textures. Finally, their model is capable of synthesizing an entire space of textures, rather than a single exemplar, when trained on a set of exemplar textures. However, the resulting texture quality is not yet on par with state of the art CNN models and the synthesized textures are often relatively different from the respective target exemplars. In our work, we also model point operations that combine noise frequencies to efficiently support diverse 3D texture synthesis. However, we show that our novel GAN loss and network architecture significantly improve image quality and similarity to target textures. Moreover, we demonstrate superior generalization performance enabled by our proposed network architecture.
Texture synthesis using GANs GANs are another widely used class of models for texture synthesis, by training a CNN generator G to fool a CNN discriminator D in an adversarial minimax game [16][17][18][19]. Most such approaches learn purely generative models, often trained on individual exemplars. An effective technique to incorporate a user-provided exemplar is conditional GANs [20], as shown by Alanov et al. [21]. Here, G synthesizes a texture given an exemplar patch as input. To enforce G to respect its input, D discriminates the joint distribution of pairs of texture patches. Here, a sample of the generated distribution consists of a condition patch and a synthesized texture, and a sample of the ground truth distribution consists of two random patches coming from the same texture. A more widely used approach to incorporate a conditional input is to extend the generator loss by an additional style loss term [3] using an off-the-shelf feature extractor (VGG) [22][23][24]. While GAN-based models achieve astonishing texture quality, even surpassing style-loss based techniques [16], the proposed models all rely on CNNs, hindering their applicability to 3D textures due to computational limitations. While we also leverage GANs in our work to obtain high-quality textures, we propose a model based on point operations, suitable for 3D texture synthesis. Moreover, we incorporate ideas from style-loss based techniques to train a conditional GAN without the need of an additional off-the-shelf network, resulting in a unified loss function that is much more efficient than hybrid losses.

GramGAN
In this section we introduce our framework for 3D texture synthesis from 2D exemplars. At the core of our method is a texture sampler that is modeled as an MLP to efficiently learn to combine learned noise frequencies using a specialized layer architecture ( Figure 2). We train this model by optimizing a novel loss function that incorporates concepts from both GANs and style loss approaches to achieve high-quality outputs. In this work we consider two different target applications: (1) learning to synthesize highly realistic textures from a single exemplar (Section 3.1), and (2) learning a texture latent space given a set of many exemplars (Section 3.2). Later in Section 3.3 we propose an efficient strategy for improving the generalization performance.

Efficient, infinite, and highly realistic 3D textures from a single 2D exemplar
Given an exemplar image of a target texture, in order to learn a model that can generate new, highquality samples of the target texture in 3D, we train a transformation model and a texture sampler. Our transformation model learns to produce n noise frequencies, or octaves T i , i ∈ [0, n − 1] , given the 2D exemplar. These frequencies are linear 3D transformations T i : R 3 → R 3 defined as T i c = c i that map a spatial 3D coordinate c to a set of new locations {c i }, one location for each frequency.
The sampler S maps the spatial coordinate c to its RGB texture value. We first compute a noise value η i (c) for each octave by trilinearly sampling n 3D Gaussian noise tensor instancesη i ∼ N (0, 1) using the transformed locations c i , i.e., η i (c) =η i (c i ). This results in an input noise vector η(c) = (η 0 (c), η 1 (c), ..., η n−1 (c)) that is combined by S to form the final RGB output at location c: During training, we simultaneously optimize both the parameters of S and the n × 3 × 3 matrix entries in the transformation model.
We parametrize S using an MLP with 5 hidden layers, as shown in Figure 2. We find that directly providing the transformed noise η at the input layer yields subpar outputs. Therefore, we propose to inject η into the hidden layers of S and learn a constant input layer instead. Consider a dense layer defined as y = φ(W x + b), where W is the weight matrix, b is the bias vector, x is the layer input, φ is the activation function, and y stores the layer activations. We use dense layers of the form: where A is a learned m × n matrix for a layer with m hidden units. This concept resembles some similarities to StyleGAN [25], although the input noise has a different role here and network layers of our sampler do not relate to scale in any way. In practice we find it beneficial to force individual layers in S to specialize on a subset of noise frequencies. This reduces the number of parameters to optimize as well as minimizing the risk of S ignoring certain frequencies in computation of its output.
In particular, we partition the frequencies in equal bins of size n/4 and inject each bin into one hidden layer only. The last hidden layer thus does not receive any noise, which we find to reduce the risk of ignoring other frequencies further, and the dimension of the matrices A reduces to m × n/4.

Loss function
We optimize the parameters of our model using a GAN loss. However, while our sampler generates 3D textures, our critic D, parametrized as a standard CNN, discriminates 2D images. In particular, in each training iteration our model synthesizes 2D images by randomly slicing the 3D space. To optimize D, we use the WGAN-GP loss proposed by Gulrajani et al. [26], i.e., where f is a randomly oriented 2D slice produced by our model, r is a random patch from the exemplar texture image, and u is a data point uniformly sampled along the straight line connecting r and f . This is a crucial difference to [5], where only axis-aligned slices are sampled during training. To optimize our texture generator G, we introduce an additional style term L style to the WGAN loss: where L is the number of hidden layers in a feature extractor network F , M l (F (x)) is the Gram matrix of the activations in layer l when presenting x as input to F , N l is the number of filters in layer l, and M l is the number of pixels in the vectorized activations of layer l. Note that while previous work on style loss employ pretrained, off-the-shelf feature extractors (typically F = VGG-19 [3, 13-15, 5, 22-24]), our important difference is that we set F = D, i.e., we train F simultaneously with G in an adversarial manner! This concept has numerous advantages compared to the traditional hybrid approaches based on combinations of GAN and style losses. First, the activations in D are anyhow computed during GAN training, hence our approach significantly minimizes any computational overhead in loss function evaluation. Second and more importantly, such features explicitly trained for texture discrimination, as we show later, help generate textures of significantly higher quality compared to off-the-shelf features trained on image recognition as in works leveraging VGG. Note that we minimize L 1 distance between Gram matrices in L style as opposed to the L 2 distance commonly used in related work, since we find this to improve the stability of adversarial training.

Learning a latent texture space
While training a distinct texture synthesis model for each target texture is computationally feasible, in practice it is desirable to have a single conditional model that can generate an entire class, or space, of textures. Besides the reduced training computation time, such a model promises generalization capabilities such as interpolation and extrapolation, greatly simplifying the texture generation process from the user perspective. Our goal is to obtain a texture synthesis framework that allows a user to select a 2D image patch of a target texture, and the model can generate new samples of this target texture in 3D without retraining. For this purpose we extend our model proposed in Section 3.1 by an encoder E(x) = z that maps an input exemplar patch x to a compact latent representation z. Similar to Henzler et al. [5], we parametrize E using a CNN. The idea is now to condition both T i and S on the latent code z. We therefore train an additional 3-layer MLP Q(z) = {T i } to map the latent representation z to a set of noise frequencies {T i }. Moreover, we extend S to non-linearly combine noise frequencies given z, i.e., S(η(c)|z) = RGB(c)|z. Instead of presenting the condition z at the input layer of S, we find it beneficial to directly manipulate the hidden activations in S as a function of z (See Figure 2). Our approach is inspired by attention modules [27] and adaptive instance normalization techniques [28] for CNNs. First, we train a small 3-layer MLP f to learn a texture style w = f (z). Next, each hidden layer l learns an affine mapping g l (w) = (γ, δ) that computes a scale and bias for each unit in the respecitve layer. The definition from Equation 2 therefore extends to Note that in order to synthesize a texture sample given an input exemplar, a single forward pass through E, Q, f , and g l is required. Only S is evaluated for each RGB value to be synthesized.

Extrapolation to unseen data
While our proposed conditional model generalizes well to unseen samples that are sufficiently similar to the training data and interpolation in latent space produces plausible textures, we find that extrapolation, the generalization performance to entirely new textures, leaves a lot to be desired. This is not surprising since generalization in the sense of extrapolation is a major unsolved problem for deep learning models. Nevertheless, supported implicitly by the architectural design of our texture sampler S, we propose an extrapolation strategy that significantly improves texture quality on exemplars not belonging to the training data manifold. Our method is a mixture between zero-shot learning and domain-adaption.
Given a new exemplar texture patch x, our idea is to fine-tune only a fraction of the parameters of our model in order to help improve the quality of generated samples. For this purpose we use E and Q to predict an initial guess for the latent code z and the transformations {T i }. Next, we estimate a texture style w and from it a scale and bias vector γ, δ for each hidden layer in S. In addition, we also predict the noise weight matrix A for the first four hidden layers in S. Finally, we perform stochastic gradient descent by evaluating random slices of the generated texture in each iteration. During this optimization we solely alter the parameters θ = (T i , γ l , δ l , A l ), with l being the layer index in S, while keeping all other parameters fixed. While also keeping the parameters in D fixed, we minimize L style for a few hundred iterations.

Experiments and results
In this section we evaluate the effect of the components of our framework by ablating them and we compare our system with a state of the art 3D texture synthesis method [5]. Note that all synthesized 2D images shown in this work are in fact slices in 3D. For more results and analysis, we refer the reader to the appendix.  Figure 3. We report SIFID µ ±σ (×10 5 ) (lower is better). Implementation details We train all models using Adam optimizer [29] and we use (equalized [30]) learning rates of 2 × 10 −3 for D and 5 × 10 −4 for G (E, Q, and S). In each training iteration both G and D are updated once and sequentially. In L D we set λ = 10. An important factor is the choice of the hyperparameters α and β in L G . Although various settings produce plausible outputs, we achieved best results when setting α = 0.1 and β = 1. Note that training a conditional model is not feasible unless β > 0. Setting α = 0 leads to perfectly stable adversarial training in both the single exemplar and the conditional setting, which we find a rather surprising insight. In the single exemplar setting we set n = 16 and our conditional models use n = 32 noise frequencies. We train on texture patches of size 128 2 and use noise instances of resolution 64 3n . E maps texture patches to a latent representation of dimension 32. Our MLPs Q, S, and f consist of 1, 5, and 1 hidden layers with 128 units each. We do not use any normalization layers, neither in MLPs nor in CNNs.
During training, we sample completely random slices when training on isotropic materials (e.g. stone textures) by rotating a plane along all three axes. For more anisotropic structures like wood, we restrict random rotations along one axis, effectively minimizing the loss only along two axes and learning the third dimension implicitly. To efficiently compute texture values at many positions in a single forward pass, we implement S using 1 × 1 convolutions. Training our model on a single exemplar takes a few hours on a single NVIDIA 2080ti GPU and learning an entire texture space takes up to several days depending on the dataset.
Single exemplar For efficiency, we perform many ablation studies in the single-exemplar setting.
In order to draw robust conclusions we experimented with two very different exemplars: wood and granite (Figure 3 left). For comparisons we consider the following four variations of our framework: GAN+VGG: We replace our GramGAN style loss in L G with VGG style loss (similarly to [22][23][24]). NI, Noise-at-Input: Instead of injecting noise as in Section 3.1, we inject it at the input layer of S. α:1, β:0: We use a vanilla WGAN loss without any style loss. α:0, β:1: We solely optimize for GramGAN loss without traditional WGAN term.
In texture synthesis, the ground truth distribution is given by a single exemplar image and the quality of synthesized images is often evaluated by considering internal patch statistics. To quantitatively assess the differences between the approaches, we measure Single Image FID (SIFID) [16], a derivative of Fréchnet Inception Distance (FID) [31] tailored to quantify the distance between internal patch statistics of two images. For each method we create 50 random 3D textures (using different random seeds) and from each texture we synthesize a randomly oriented 2D slice at the resolution of the ground truth exemplar. We then compute SIFID for each of the 50 resulting sample images and report their mean and standard deviation in Table 1, with qualitative examples shown in Figure 3.
The examples show that all variants of our approach are capable of reasonably reproducing the target texture. However, we achieve significantly more plausible results with our proposed method (last column in Figure 3). In addition, we also train the model by Henzler et al. [5] on both exemplars and find that it produces significantly lower-quality results (second column in Figure 3). The numbers in Table 1 confirm these findings. While the results suggest that L style is more important for the wood exemplar, WGAN loss is more effective on granite. However, note that solely optimizing for L style , i.e. α:0, β:1, can cause artifacts that resemble Deep Dream [32] effects (tiny replicas of wood patterns when zoomed in Figure 3). For NI we use α:0.1, β:1 and for GAN+VGG we chose the hyperparameters to approximately result in the same weighting as ours with α:0.1, β:1.
Interestingly, we found that NI results in significantly less stable adversarial training and we had to lower the learning rate by a factor of 5 to 1 × 10 4 to achieve convergence. In conclusion, using our GramGAN loss (α:0.1, β:1) yields consistently competitive results without any artifacts for the given substantially different exemplars.
We show additional qualitative examples of our model in Figure 1 and 4, including a comparison to [5] on actual 3D renderings using the synthesized textures as diffuse albedo. The comparison in    Figure 4 shows that [5] produces implausible artifacts in 3D for anisotropic textures. This is a cause of sampling only axis-aligned slices during training, while we sample all possible slice orientations. The renderings are raytraced using mitsuba renderer [33] modified to perform inference on our model for texture look-ups. The cube model is rendered using perfectly diffuse albedo and light source, while the sphere model rendering features a more realistic scenario with a glossy coating material and an HDR environment map light source. The rendered images are high resolution, so please zoom in.
Texture space Next we report a quantitative evaluation of our conditional texture space model in comparison to the system proposed by Henzler et al. [5]. For this purpose we trained both models on a texture dataset consisting of 100 stone texture exemplars collected online. In this conditional setting, the input to both frameworks is an exemplar patch of resolution 128 2 and the systems learn to resynthesize the input texture. For evaluation we synthesized texture samples of size 128 2 from 22 reference input patches not seen during training, and measured their similarity to the references.
In Table 2 (first row) we report SIFID values averaged over the 22 samples generated by both frameworks. In addition, we also measured ground truth performance by sampling random crops from the respective training exemplars and compare them to the references. The numbers confirm that the outputs of our model are closer to the references in terms of SIFID. In addition, the comparison shows that SIFID for ground truth is not significantly lower. This can be explained by the fact that two random patches of the same texture exemplar image can look fairly different, especially for textures that feature very low-frequency variations. Note, however, that the computed SIFID numbers have rather high standard deviations due to the limited resolution of the texture images. In addition to SIFID we also measured average log-likelihood (ALL) reported in Table 2 (second row). An important property of generative models is to avoid producing artefacts and ALL has been proposed as a measure thereof [34], quantifying how probable generated samples are, given the ground truth density. ALL is typically estimated by using a Parzen window density estimate using a Gaussian kernel. Due to the stochastic nature of textures, we perform this evaluation over rather small kernel sizes to emphasize on the texture elements. In particular, we start with a large kernel size and reduce its size until there is a significant difference between the ALL of the methods. To set the second parameter, Gaussian kernel bandwidth, we performed a grid search using the ground truth test set. Our evaluation shows that our method is superior also in terms of ALL, although not as significant.
Finally, we conducted a formal user study to evaluate both the perceptual quality of the generated samples and their similarity to the reference. For this purpose participants were presented with a reference texture patch and three candidate patches: the ground truth patch, a patch generated by Table 2: Quantitative evaluation of our framework in comparison to [5] and ground truth. We report µ ±σ for SIFID (lower is better), ALL (higher is better), and user study ranking (lower is better).  our model, and a patch produced by [5]. The participants were asked to rank these according to their similarity to the reference patch. For sake of clarity, we elaborated similarity as how likely that the candidate is an image of the same material as in the reference image, e.g. perhaps a different region/section of the same sample. The presentation order of textures and candidates were randomized for each participant. We received rankings from 21 participants for 22 texture references. We report the average ranks of 462 rankings per method in Table 2 (third row). A paired Wilcoxon signed rank test confirms that the users prefer our results over [5] (with p < 2 × 10 −11 at 5% significance level). Moreover, the difference between our results and ground truth is not statistically significant (p = 0.56), whereas ground truth is significantly higher ranked than [5] (p < 1 × 10 −10 ). In Figure 4 we show a qualitative comparison of our texture space model on an exemplar texture patch with a texture significantly different than all training exemplars. This comparison shows that both [5] (first row) and our framework (second row) produce results fairly different from the input exemplar due to lacking extrapolation capabilities. Nevertheless, quantitative evaluation (SIFID over 50 random samples) indicate that our result (2.04 ±0.79) is significantly closer to the reference than [5] (11.31 ±0.79). Finally, we demonstrate that our extrapolation strategy proposed in Section 3.3 (forth row) further improves the output significantly (1.46 ±0.29). We also show an alternative of our extrapolation strategy (third row) that minimizes L style using VGG instead of D. This experiment demonstrates that our discriminator features are much better suited than VGG features for style-based texture synthesis (SIFID: 1.79 ±0.3).

Conclusions
In this work we presented GramGAN, an novel framework for 3D texture synthesis from 2D exemplars. Key to our method is a novel loss function that combines ideas from GANs as well as style transfer, to optimally guide a synthesis model. We demonstrate that our model significantly outperforms state of the art, with superior results compared to existing alternatives. An interesting future direction would be to extend our framework to non-diffuse albedos including specular and normal maps by leveraging differentiable renderers.

A Appendix
Here we show additional results and analysis of our framework for 3D texture synthesis. We first demonstrate the capability of our system to learn a continuous latent texture space when trained on a dataset consisting of diverse textures (Section A.1). Next, we present more qualitative results that demonstrate the benefits of our approach compared to the system by Henzler et al. [5] (Section A.2). In Section A.3, we tabulate the network architectures of the convolutional neural networks used in our experiments. Finally, we present detailed results from our user study and all the images used therein in Section A.4.  Here we demonstrate artifacts that occur due to non-axis-aligned slicing. If only axis-aligned slices are sampled during training ( [5]), distortion artifacts occur along non-axis-aligned slices. Left: Henzler et al. axis-aligned cube (first and third row) and 45 • rotated cube (second and fourth row). Right: ours axis-aligned cube (first and third row) and 45 • rotated cube (second and fourth row). Note how ours produces more consistent results for both axis-aligned and non-axis-aligned geometry. See Figure 10 (second row) and Figure 11 (third row) for exemplar images.