Combinatorial Bandit Bayesian Optimization for Tensor Outputs
Authors: Jingru Huang
Haijie Xu
Jie Guo
Manrui Jiang
Chen Zhang
Venue: ICLR 2026 Poster (2026)
Keywords: Tensor data
Non-separable kernels
Gaussian process
Bayesian optimization
Combinatorial multi-arm bandit
Upper confidence bound
Primary area: probabilistic methods (Bayesian methods, variational inference, sampling, UQ, etc.)
TL;DR: We develop BO and combinatorial bandit BO frameworks for tensor-output systems, built upon the proposed tensor-output GP with a non-separable kernel as the surrogate model.

Abstract: Bayesian optimization (BO) has been widely used to optimize expensive and black-box functions across various domains. However, existing BO methods have not addressed tensor-output functions. To fill this gap, we propose a novel tensor-output BO framework. Specifically, we first introduce a tensor-output Gaussian process (TOGP) with two classes of tensor-output kernels as a surrogate model of the tensor-output function, which can effectively capture the structural dependencies within the tensor. Based on it, we develop an upper confidence bound (UCB) acquisition function to select query points. Furthermore, we introduce a more practical and challenging problem setting, termed combinatorial bandit Bayesian optimization (CBBO), where only a subset of the tensor outputs can be selected to contribute to the objective. To tackle this, we propose a tensor-output CBBO method, which extends TOGP to handle partially observed tensor outputs, and accordingly design a novel combinatorial multi-arm bandit-UCB2 (CMAB-UCB2) criterion to sequentially select both the query points and the output subset. We establish theoretical regret bounds for both methods, guaranteeing sublinear regret. Extensive experiments on synthetic and real-world datasets demonstrate the superiority of our methods.

Published as a conference paper at ICLR 2026
COMBINATORIAL BANDIT BAYESIAN OPTIMIZATION FOR TENSOR OUTPUTS
Jingru Huang", Haijie Xu", Jie GuoD, MMnrrui Jiang", Chen Zhang
" Department of Industrial Engineering, Tsinghua University, Beijing 100084, China b College of Economics and Management, Nanjing University of Aeronautics and Astronautics, Nanjing 211 106, China
Seggraahhang@tsiiggghu edu. cn, xu- han-hii22mmilsstsinnnnie tsiinguu....... edu. cn, guojie25enuaa
hwwp:/eeenaaaauuuuu.... edu. en, jiangmanrui @mail t tsinguua. edu. shangchenOl@tsinee edu. en
ABSTRACT
Bayesian optimization (BO) has been widely used to optimize expensive and black-box functions across various domains. However, existing BO methods have not addressed tensor-output functions. To fill this gap, we propose novel tensor- output BO framework. Specifically, process (TOGP) with two classes of censor-output kernels as a surrogate model of the tensor-output function, which can effectively capture the structural dependen- cies within the tensor. Based on it, we develop an upper confidence bound (UCB) acquisition function select query points. Furthermore, we introduce a more practical and challenging problem setting, termed combinatorial bandit Bayesian optimization (CBBO), where only a subset the tensor outputs can be selected to contribute to the objective. To tackle this, propose tensor-output CBBO method, which extends TOGP to handle partially observed tensor outputs, and ac- cordingly design novel combinatorial multi-arm bandit-UCBB (CMAB-UCB)) criterion to sequentially select both the query points and the output subset. We establish theoretical regret bounds for both methods, guaranteeing sublinear re- gret. Extensive experiments on synthetic and real-world datasets demonstrate the superiority of our methods.
first introduce a tensor-output Gaussian
INTRODUCTION
Bayesian optimization (BO) is a widely used strategy for optimizing expensive, black-box objective functions (Frazier, 2018, Wang et al. 2023). Its effectiveness has led to successful applications in various domains such as hyperparameter tuning, experimental design, and robotics (Snoek et al.) 2012 Shields et al. [022 Wang et al., 2022). Most existing BO methods focus on scalar outputs Bull 2011 Wu et al. 2017), while recent studies have extended BO to handle multi-output settings Chowdhury & Gopalan 2021 Tu et al. 2022. Maddox et al. 2021 Song et al. (022). However, to prior work has addressed BO with consorvallued outputs, where the system response is multi-mode tensor. In contrast, consor-allee data have been extensively ex-
the best of our knowledge,
plored in other areas, including tensor decomposition (Abed-Meraim et al. 2022 tensor regression [oock] 2018), and tensor completion (Song et al., 2019), among others.
In current outti- BO (MOBO) methods, surrogate model, typically multi-outpur Gaussian process (GP) or a collection of independent scalar-output GPs, is constructed from observed data, and acquisition function is then optimized to sequentially select query points by balancing explo- ration and exploitation. A straightforward way to handle tensor- tensor-valued outputs is to vectorize the tensor and apply existing MOBO methods. However, vectorization ignores the intrinsic structural correlations of tensor data, especially mode-wise dependencies that are critical in many applications. As a result, MOBO methods may become less effective when optimizing the acquisition function to identify the global optimum. As such, it is important to develop GP surrogates that model ten- sor structure directly. Existing tensor-output GP models typically adopt fully separable covariance
"oorrsponniig author.

---

Published as a conference paper at ICLR 2026
structures, in which the joint covariance factorizes Kronecker product of covariance matrices across the input and each output mode (Kia et al. 2018 Zhe et al. 2019 Belyaev et al. 2015). While computationally attractive, separability implies that correlations across tensor modes do not vary with the input, an assumption that is often violated in complex systems such as spatiotemporal processes (Hristopulos) 2023). This mismatch can lead to poor predictive accuracy, numerical in- stability in posterior inference, and ultimately degraded BO performance. Therefore, in this paper, first aim construct more flexible and scalable GP model that can capture input-dependent correlations within tensor outputs. Then, we aim to design an acquisition function and a sequential querying policy tailored for tensor-output BO.
Furthermore, we consider a more complex and practical setting in which only a subset of tensor entries contributes to the objective. This naturally transforms the problem into combinatorial multi-armed bandit (CMAB) setting. Specifically, each tensor element is treated as an base arm, and each round, select a subset of arms, referred to as a super-arm. The objective value associated with the chosen super-arm is then observed, which we interpret as the reward. We term this novel problem as combinatorial bandit Bayesian optimization (CBBO). The goal of tensor-output CBBO to jointly identify optimal input in the search space and the corresponding best super-arm over tensor outputs. Recent work has combined BO with multi-armed bandits (MAB), often termed bandit Bayesian optimization (BBO), handle mixed input spaces with both continuous and cafegorical variables (Nguyen et al.] 2020 Ru et al.) 2020)) Huang et al. 2022). In these settings, categorical variables can be viewed as indexing discrete modes, and choosing one category per variable amounts to selecting single arm along each mode, which is a special case of CBBO. However, existing BBO methods do not be directly extend our CBBO setting for two reasons. First, they typically model the effects of categorical choices using independent GPs, and thus cannot capture the rich structural correlations inherent in tensor outputs. Second, their selection strategies effectively decompose the decision into multiple independent bandit problems (one per categorical variable), whereas CBBO requires joint selection of multiple, potentially correlated arms within a super-arm. Thus, these methods are not well suited for our tensor-output CBBO framework.
To address the aforementioned challenges, we propose a tensor-output Bayesian optimization frame- work, termed TOBO, together with its extension to the CBBO setting, termed TOCBBO. Our main contributions are summarized as follows:
We propose tensor-output Gaussian process (TOGP) with two classes of tensor-output kernels. The proposed kernels explicitly incorporate tensor structure by extending the linear model of coregionalization from vector-valued censor-valued outputs, thereby capturing dependencies both across tensor modes and over the input domain.
Using the TOGP model as the surrogate, we develop a TOBO framework based on the up- per confidence bound (UCB) acquisition. We establish sublinear regret bound for TOBO, which provides the first Bayesian regret analysis tailored to valuue- outputs.
We formulate novel problem setting, referred to as CBBO. To address this problem, we design the TOCBBO framework by introducing a CMAB-UCB2 acquisition function, which couples UCB-based input selection with CMAB-UCB M-Ms--UBB-assed super-arm selection. We further establish a sublinear regret bound for TOCBBO.
We demonstrate the efficiency of our methods synthetic experiments and case studies.
Notably, compared with existing TOGP methods (Belyaev et al. 2019) our model provides a more general kernel-construutioon framework. Specifically, existing tensor-output kernels can be interpreted as special cases induced by particular low-rank tensor de- compositions, whereas our L.Me-based formulation allows a broad class of tensor constraints to be encoded through the coregionalization structure. Compared with standard BO methods (Srinivas) jet al. 2009 Belakaria et al.) 2019 Chowdhury & Gopalan 2021). our work is the first to establish a BO framework for tensor outputs and further extend it to the proposed CBBO setting. Moreover, our contributions lie in deriving regret bounds for both TOBO and TOCBBO based on concentration inequalities specialized the proposed TOGP under Bayesian framework.
2015 Kia et al. 2018 Zhe et al.
2

---

Published as a conference paper at ICLR 2026
2
RELATED WORKS
High-order Gaussian processes (HOGPs): Most existing studies for modeling HOGPs adopt sep- arable kernel structures. In particular, |elyyev et al. (2015) proposes tensor-variate GP with sepa- rable covariance across tensor modes, Kia et al. (2018) develops a scalable multi-task GP for tensor outputs by factorizing the cross-covariance kernel into mode-wise and input components, and ZZee] et al. 2019) introduces a scalable high-order GP framework based on Kronecker-structure kernels. Although these formulations simplify computation, they typically imply that correlations across tensor modes do not vary with the input, which limits their ability to capture input-dependent de- pendencies. Recent work has also developed non-separable kernels for multi-output GPs (MOGPs), including convolution based kernels (Fricker et al., 2013). linear model of coregionalization (LMC)
kernels (Li & Zhou 2016 Bruinsma et al. 2020), and linear damped harmonic oscillator-base
kernels (Hiisoouloss 2023). A closely related viewpoint is multi-task GPs, where each task cor- responds to one output component (Yu et al., 2018: Chowdhury & Gopalan [02T]] Maddox et al.) 2021). However, these methods are primarily designed for vector-valued outputs and therefore can- not fully exploit the intrinsic multi-mode structure of tensor data, see Section []
Multi-output Bayesian optimization (MOBO): MOBO typically refers to either multi-task BO or multi-objective BO (Frazierr 2018: Wang et al. 2023). A desirable property of BO algorithms is no-regret, namely achieving cumulative regret R(T = (T)) after T rounds (Srinivas et al. 2009 Chowdhury & Gopalan 2021). Recent work has developed various MOBO methods with theoreti- /0, Chowdhury & Gopalan (222)) : Dai et al.)
cally grounded acquisition functions. For multi-task BO,
(2020); Sessa al. (202)) employ UCB-type acquisition rules and establish O(VT) regret. For multn-bjecttive BO, Belakaria et al. (2019); Zhang et al. (2025) study max- max-value value entropy search
and obtain O (VT)) regret, while Daulton et al. (2022a) proposes a crust-region-basee criterion with
OV ) regret. In parallel, large body of MOBO methods emphasizes empirical perfor-
mance without regret guarantees, including improvement-basee criteria (Uhrenholt & Jensen 20199]
Daulton et al. 2020: 2022b), entropy-based search (Hernández-Lobate et al. 2014: 2016 Tu al.,
2022), and inf information gain- based methods (Chowdhury & Gopalan (0221) Although effective for correlations of tensor outputs and are therefore less suitable for our TOBO setting. Moreover, hypervolume- based and setrropy-eearc-bbaeeee multi-objective BO methods become computationally prohibitive as the number of objectives grows, which makes them impractical for tensor-output optimization.
multi-output problems, these methods typically do not exploit the intrinsic
Bandit Bayesian optimization (BBO): BBO combines BO with multi-armed bandit (MAB) al- gorithms to tackle optimization in mixed input spaces containing both continuous and categorical variables. Nguyen et al. 2020); Huang et al. (2022) study a setting with one continuous and one cat- egorical variable, coupling BO with Thompson sampling for the categorical choice and establishing O(VT*+I log T) regret. Ru et al. (2020) extends this idea to multiple categorical variables, named CoCaBO, by employing EXP3 (Auer et al. 2002) select categories and achieving O(VT log T') regret. Such settings can be viewed as special cases of CBBO, yet they do not directly extend to our problem due to the tensor-output structure and the need to jointly select multiple correlated arms.
Table D
summarizes a comprehensive comparison between our method and existing literature.
Table Comparison of related works and our proposed method
GP
Type
Liserature
Tensor structure
Seganbbee (Indepenee:: modes)
Noo-scurable (Crucs mode cocrelation)
DO (Regret
BBO (Regret)
CBBO (Regret)
Belvaw et a
2015
HOGP
m
A013
2m1
A034
MUGP
MTB0
/
v/ Wo
MTB0
C houcchrry & Chrolm
2011
x
yf (OIVTT)))
MOBU
Hglakan 21
6025
2019
x
x
v' VOO Zog
BBO
M0
x
/
x
v V10ty/ToTlo (ogT))
8B0
Xu ot al 5020
ViolVTlog Zi)
ViolyTing
TORO+TOCBRO
propoeed
V
V10lyTlog71)
3

---

Published as a conference paper at ICLR 2026
3 TENSOR-OUTPUT BAYESIAN OPTIMIZATION
In this section, we propose novel tensor-output Bayesian optimization (TOBO) framework for optimizing expensive black-box systems with tensor-valued responses. Let f X - y denote the objective function, where the input æ (II...., Td) lies in compact and convex domain X C Rd. and the output f(æ=) y C (g..xx...fm is m-mode tensor. Denote f..... fm (3) as the (i1...., )-th entry of the tensor, where i; = L..... & for I -,..., m, and let T I t =,......
be the total number of elements. To optimize tensor-output systems, an intuitive way to map the tensor-valued objective into a scalar function. To this end, we introduce bounded linear operator
I, E L( RR,, where L(), TR) denotes the set of bounded linear operators from y R. The optimization problem is thus
x = arg max (() ΦE
(1)
To solve this problem, the proposed TOBO aims to sequentially select inputs æ; and observe the corresponding tensor outputs,
= f (x;) + Ee.
Vi = 1,
where E X denotes i.i.d. measurement noise with vec(e,) ~ ~ and vec y R² is the vectorization operator.
Based on the collected data, we fit a tensor-output Gaussian process (TOGP) surrogate for f using two classes of tensor-output kernels, as detailed in Subsection 3.1 Then, we develop a UCB-based acquisition strategy efficiently identify the maximizer x", as presented in Subsection 3.2
3.1 TENSOR-OUTPUT GAUSSIAN PROCESS
Define the prior of f X + y as tensor-output Gaussian process (TOGP):
vec (f(x)) TOGP (z,,"
Va, x' E X.
(3)
where μ E IRT is the prior mean, o² > 0 is a variance hyperparameter, and K(x, x') E RTT is symmetric and positive semi-definite kernel function. The classes of K(xx, x') are discussed later.
Given n observations X,, = (2I,..., ;n) and Y, = =veecyy)) ..... vec(y,)
"e.,,,. T)² the posterior distribution of vec(f (x)) at new imput X is Gaussian in IR with mean and covariance
(x)
K (x) (K, nIuT)= (Y, 1, 1))
(4)
K (ax, 2)) o² K K( a)) K. (a:) (K,, + nITT) Kn (z))
(x))
(5)
where (x) RT is the block column matrix whose i-th block K(x,, x)), K,, E Ranx the block matrix with (i.j)-block K(x,, 1, is the n-vector of ones, and 17 r 12/022
To specify the tensor-output kernel in TOGP, we introduce two classes of kernels. The first class consists of non-separable tensor-output kernels constructed via tensor extension of the linear model of coregionalization (LMC) (Fricker et al. 2013, Li & Zhou 2016). Let vec(f(x))
m Σラ a₂ Mey (x), where aij RT are loading vectors and {ugj} are and covariance functions Cov(uej(a)) uey (2)))))
independent scalar GPs with zero krj( (xc, x²). We parameterize the loading vectors by tensor cores, namely arj = vee( Atj) with ArJ Rtix Rtix..xtm Then f(z) has zero mean and covariance Cov (vec(f(z)). vec(f(x")) EA Ei vec(Ae) vec(A.) kej( æ'). This yields a non-separable covariance structure since the cross-entry correlations are governed by mixture base kernels {kez} and therefore can vary with the input. The separable tensor-output kernel is defined as
Definition 1. Define the non-separable tensor-oulput kernel K(ær, x²) for any x. x' E X:
K(c, x') =
vcc(Ai)voc(A)) kiz (2, (x,
(6)
1=1 1=1
where Ar X.lxx..xt Xtm tensors and kij (xx x!) are base kernels X (e.g. Matérn).
4

---

Published as a conference paper at ICLR 2026
Furthermore, if all Wa (z) share the same covariance 22, for in [t.]. 1,.... ..... m, and the convolution degenerates to ri..... (x) = A:........
dm 5(2 - =), then the induced tutpor- kernel reduces to separable structure in which the correlations across tensor modes are independent of the input.
Definition 2. For any x. E X, define the separable lensor-oulput kernel
K(x, o²) vec(A) vec k(xx
where A E X2 is core tensor, and k(.,)) is base kernel on X.
Note that (xx, æ') is designed to capture both input correlations and dependencies among tensor entries. For the non-separable kernel in 66. each base kernel kı,(2; x') models correlations in the input space, while the matrix vec( (A vee vec(Ai) (At) encodes the dependency structure across tensor entries. For the separable kernel in (7). k(x, , x²) captures input correlations and vec(A) (e)eee())) specifies the dependency structure within the tensor output.
Proposition 1. The kernels in Definitions 7 Mandd
symmettrcc and positive semidefinite on X, i.e., (2) Vn E N, any It ;1....,Zay
(1) Va, x' E X, K(x., x') = K(x', x)) D. Eij=i vecly:) K( F. (vec(y)) 0.
En E X, and any y1...., E
Its proof is given in Appendix E Proposition I ensures that both classes of tensor-output kernels define valid covariance functions for TOGP, and hence the induced TOGP prior is well defined.
The core tensors {A in (6) and A in ( contain mil and free parameters, respectively. To reduce the parameter complexity low-rank decomposition based structures can be imposed on these cores, such as CANDECOMP//ARAFAC (CP) decomposition (Goulart et al. 2015 or lensor-train (TT) decomposition (Oseledets) 2011). Details are shown Appendix
Remark 1.
Denote the hyperparameters by e = a0, For the non-separable kernel in (6), 0
mad (0t are the parameters of the base kernel { and {ans {aug)i,s 1immm are the entries of (AI] For the kernel in 77. 0 (0I...., Oa} is the parameters of k(. .), and a {ai, Jaj mmth are the entries of A. We estimate @ by maximum likelihood. The detailed estimation, algorithm, and complexity analysis are provided in Appendix C
The ranks of the core tensors {Ae] are selected via cross-validation. Let R. (ri....,Tel Tel be the set of candidate ranks, we fit the TOGP model and evaluate its predictive per- formance held-out validation using the mean absolute error (MAE) criterion. We select
Remark 2.
f.- ||. The data-driven choice bal- 2i=1
r* = arg minFER, MAE(r), where MAE(r)
ances model flexibility and complexity, mitigating overfitting while preserving expressive capacity.
Remark 3. To improve scalability for large n., adopt Nyström low-rank approximation of the block covariance matrix Kn E 18T xn'," Specifically, it selects nh < v'l' landmark columns to construct a rank-nı approximation Kn Fs ÜÄÜT, which can be viewed as approxi- mating the leading spectral components of Kn. We choose the effective rank via cumulative ex- plained variance, i.e., l min {1o El....,, >/ d: 2. c). where {di} denote
the eigenvalues of Kn in non-increasing order. Then, evaluating (K,, + mTITT) can be reduced OnT n2 + ni)) WWilliams & Seeger 2000). Thus, the overall training complexity is reduced to O (InT ni n² n TT logn).
3.2 UPPER CONFIDENCE BOUND ACQUISITION STRATEGY
Building on TOGP as a surrogate for f, we develop a UCB-based acquisition strategy. At round n 1, given past observations X,, and Y,, we update the hyperparameters θ,, and the posterior mean (4) and covariance The UCB acquisition function for the scalarization- objective in is defined as
auca( D.) Ly (x) + 3n (x)1//
where B, > 0 balances exploration and exploitation. This criterion encourages exploration di- rections with greater predictive uncertainty under the tensor-output setting. The next query is then
5

---

Published as a conference paper at ICLR 2026
selected by maximizing (8
Pp+1 = arg max (u Dr). =EX
The complete TOBO algorithm and its complexity analysis given in Appendix D
We now study theoretical properties of TOBO under the following conditions. We first define two commonly used notions of regret.
Definition 3. At each round n,, the TOBO method selects queried input E X. The instantaneous regret is defined as Tn (x=) I, (m))), and the cumulative regret up to round N is defined as RN = En [Lf((=")= Lff(m,)).
The regret quantifies the gap from not knowing the objective in advance. A good strategy can achieve sub-linear cumulative regret so that the average regret per round converges zero N +
Definition 4. At round n, TOBO selects an input In E X. The instantaneous re- gret is In = (ff(x*) Lff(2,), and the cumulative regret up to round N Rx = E fff((=)==>>>)>)))())))))(()))[ L,
The regret quantifies the gap due to not knowing the objective in advance. A good strategy achieves sublinear cumulative regret that the average regret RN /N converges zero as N 30.
Assumption 1. The unknown function f is TOGP prior with kernel K defined in 66
Assumption 2. The scalarized objective a + [xx(( is L-Lipschitz on X, ie., for any xx E X, [f(({():) - (f([rr))) S L[(x)) f(æ;))] holds.
Assumption D specifies the Bayesian setting under which follows TOGP prior. Assumption 2 ensures that the scalarized objective is stable under small perturbations of the input.
Lemma 1. Let Qvec(f)/0zz E B be the gradient of vec(f) with respect to the j-th coordinate of E X. Then, Quee(f))/rz) is a GP with covariance K (Iz, z²y). Under Assumptions given data X,, and Y, with 2 1, there exist constants a, b >> such that for any L' > 0,
sup [[Qvec())/az T!' C+) <exx(__12118), j == ;;,,..,d ΦEx
Pr
(10)
where C = supzEx
tn RY (r),,))).
Its proof is given in Appendix Lemma shows that the derivative of the vectorized TOGP remains Gaussian and holds piobabbbbbylltyy confidence bounds.
Theorem 1. Under Assumptions 2 define Cn sup=Ex Cr(Kn(z, x))mx))mmmmxe( where (an)(())
(n) (x) is the largest eigenvalue of Kn( xz). Suppose X C 00, Then, for any 8 E (0, TOBO with 3n VOn + 2d logog(r (b.log((da/0) + Cv) holds that, Pr (RN L (VCINTN 1) BN ))) 1-8, where C1 0 a constant and iN (K, m 2 log det (INT +1 + KN) denotes the maximum information gain.
Its is given in Appendix G Theorem I implies that TOBO achieves sublinear cumulative regret with high probability.
Proposition If K(x,, x') is the separable kernel in Definition 2 then the maximum information gain satisfies 7n (K. 17) (T log(n)d+1) when k(x, x') is a Gaussian kernel, and "nn (K, n) log(n)) when k(x, ) is a Matérn kernel with v > 1. Details are provided
in Appendix H
TENSOR-OUTPOR COMBINATORIAL BANDIT BAYESIAN OPTIMIZATION
We consider challenging setting in which only k T entries of the tensor output contribute to the objective. Formally, define f(x A) e(A)vec(f(z)) D. where X

---

Published as a conference paper at ICLR 2026
(A1...., >T) is binary indicator vector in A f fo. [0.1)" = and e(X) E 10,1 is a binary selection matrix whose j-th row selects the 1 coordinate of vec(f(x)). The goal is to jointly identify x" E X and the optimal subset of k elements, encoded by X A, that maximize
(a' X+) = arg
max mEA,AEA
H f A),
(11)
where y is bounded linear operator H, C(Y, R). By interpreting each tensor element base arm. any X E corresponds to a super-arm S S [T] of size k. where j E S if > = 1 and j $ S otherwise. At each round i [N], the learner selects an input æ; together with a super-arm S: (i,,..., }- The observation is the partial output Y; E RR. indexed by Si, while entries with indices j $ S, remain unobserved.
In this section, we develop tensor-output combinatorial bandit Bayesian optimization (TOCBBO) framework solve (II). In Subsection 4.1 we extend TOGP for partially observed outputs. In Subsection 4.2 we develop an efficient AMAB-UCB2 acquisition strategy that couples UCB-based input selection with CMAB-UCB-basse super-arm selection.
4.1 PARTIALLY OBSERVED TENSOR-OUTPC" GAUSSIAN PROCESS
From the proposed TOGP in 3 the prior of f xxA + ý a partially observed TOGP (PTOGP):
f(z, ~ PTOgP(e()) re((A) x' )e Va, X. X. E A.
(12)
Let X,, (2,,...,,..)) TT,", A,, (AI, TT) and Y, (91,...,, Dn) be partial observa- tions, where X; encodes the selected super-arm S; and Dr = f( + with é; ~ ((H,+9Aa)
Then, for new pair (x, A), the posterior of f(x, A) is Gaussian in B with mean and covariance
(2, A) (()) + of (A)K, (z E E (Ye) E. (1, @ 11)
(13)
K2 (a. n X X = oe[e(A)K(a x' )e oee(A) (x) E, E,,,, (x,') (14)
where Σn o² o'E, K, + II and En E RokxnT n n block-diagonal matrix with the i-block given by e (X:). It is easy to verify that the posterior covariance is symmetric and pos- itive semidefinite. We estimate the hyperparameters of PTOGP via maximum likelihood. Detailed estimation, algorithm, and complexity analysis are presented in Appendix C
4.2 CMAB-UCB2 ACQUISITION STRATEGY
Building on the PTOGP, we develop TOCBBO to sequentially select query inputs f.i..... xx together with their associated super-arm indicators {'i...., Av). Directly optimizing over both æ and À is computationally intractable, since identifying the optimal super-arm of size k from T arms requires a combinatorial search over (1) candidates, which becomes prohibitive when cou- pled with continuous optimization over X. To address this challenge, we propose a CMAB-UCB2 criterion that decouples the joint optimization into two sequential steps.
At round n 1, let (a A denote the best pair observed far, that is, {
arg max (a. i= f(r,, A1). In the first step, we fix the super-arm to A and select the next input by maximizing a UCB acquisition function under this fixed super-arm:
En+1 arg max H p A1) Bm K (zz, x X X;)1 =EX
(15)
In the second step, given the chosen input In+1, selecting An+1 reduces to a combinatorial bandit problem. To this end, we adopt the CMAB-UCH criterion by constructing a UCB for each super-arm and selecting the one that maximizes the upper confidence value:
An+1 arg max H pr (Ix+1, A) + pn (T +1 X, A) XEA
(16)
Here, 3r and oa. are tuning parameters controlling the trade-off between exploration and exploita- tion. Detailed algorithm and its computational complexity analysis are provided in Appendix D
We further analyze the regret bound of TOCBBO under the following conditions. The regret for CBBO is defined follows.
7

---

Published as a conference paper at ICLR 2026
Definition 5. At round n. TOCBBO selects input an imput In X and super-arm S,, encoded by An A. The instantaneous regret is Tn Hff( X²) - Hy An), aand the cumulative regret up to round N is Rx =I #af(x* ( faf(
Assumption 3. The operator H in 11 is H H-Lisschitz with respect to f under the l2 norm, i.e., for Vo,,2; 2jj E X and As. X3 E A, I X.) y f >3)1 S H[l >:) f( X,)I holds.
Assumption 3 ensures that H varies smoothly to changes for partially observed tensor outputs.
Theorem 2. Under Assumption 1 and Assumption 2 denote Ön SuPxEX
(x, ..))))
A (2,
(x) Xmax Xn) is the largest eigenvalue of Kn w: Xn N,). For any 7 > O. TOCBBO with Bn = V.n + 2d rdn" /by/ log(da/ ()++)) 8) and Pr. 2log NTa. b
where
5 (0, 1) and
holds that, Pr Rx S H (VCoNTN (K, n)BN + of VCTTNNN (CTTNNN(K,M))))
+
2
where C3, Cs
0 are constants, Tn > O is a sequence such that Σn= 1/5tn = maxxycX 1 log det (IEN + E 2
and maXANCA 2 log det (I&N +1 EN is the maximum information gain for super-arms.
1, 7n (k,K))
7N (x(K))
The proof is provided in Appendix [ upper bound with high probability.
Theorem Soatremm@un immlies that TOCBBO achieves a sublinear regret
Proposition 3. Suppose the kernel K(x, x') follows the separable structure specified in Def- The maximum information gain 7nn (n,K,) n) and 5n(K, 17) are O (logg))++)) and
inition
T n 20 logn) for the Gaussian kernel and the Matérn kernel (v > 1), respectively, where d denotes the imput dimension. Details are provided in Appendix J
O
EXPERIMENTS
We evaluate TOBO and TOCBBO on both synthetic and real-case data, and compare them with baselines where the tensor output is vectorized and MOGPs used surrogate models. Specif- ically, we consider three GPs in the literature: (1) sMTGP: the scalable multi-task GP Kia et al. (2018) (2) MLGP: the multi-linear Yu et al. (2018;; and (3) MVGP: the multi-variate Chen ct al. (2020). For each GP, we examine two sequential BO sampling strategies: UCB-based BO and random sampling. In addition, we replace the UCB acquisition in TOBO and TOCBBO with random sampling to construct an ablation baseline, denoted as TOGP-RS. Details of all baseline settings are provided in Appendix K
5.1 SYNTHETIC EXPERIMENTS
We generate synthetic functions of the form f(x) B 81 U1 O2 ... 8m-1 Um-1 8m g(x), where each element of B E RR Pm is independently sampled from (0, 1). the ij-th element of
U, E xPi x1: defined as li cos(ijl/2) + sin(li), and g(x)) = (5x),, cos(x)) E R Pm.xTmmm Here
3, (Ti, T2,TTT Ts) (2,4,2), P2, P3) = 33,3,3); (2) m = 2, (T) = (3, 2). (PI. = (3,2), and (3) m = 3, (3, 3, 3). Observations are collected as y, f(x,) + Ei.
Pm d and x [0. 1]d. We consider three parameter settings: f(z): (1)
(2, (4, 2),, Ps,
where iid NO PI). For CBBO tasks, we set k = TT/6.
We generate Prrain 10d training samples and Niest 5d testing samples via a Latin hypereube design (Santner et al. 2003). The training set is used for hyperparameter estimation, and predictive performance is evaluated on the test set in terms of NLL, MAE, and ||Cov]]. with details provided in Appendix K To balance modeling flexibility and computational cost, use the separable tensor- output kernel in 7 for Settings (1) and (2), and the non-separable tensor-output kernel in (6) for Setting (3). The results are summarized in Table 2 As shown. our proposed method achieves the lowest NLL and MAE, indicating that the TOGP model provides the highest prediction accuracy.
Among the three baselines, ML.G performs the worst because its fitted covariance matrix becomes
8

---

Published as a conference paper at ICLR 2026
Table 2: The prediction performance of GPs in the three synthetic setrings.
Setting (1)
Setting
Setting (3)
NLL
MAE
||Cov]]
NLL
MAE
||Cov]|
NLL
MAE
[|Cov]
TOGP
503.0
0.1571
2.02
-18.1
0.1052
0.04
-3923.1
0.1372
2.82
MMGG
749.4
0.1684
1.44
-5.0
0.1566
0.06
-37430
0.1501
22.01
MLGP
707937.1
0,9428
67.00
7066.9
0.8789
5.12
-55800.7
1.1670
0.06
MVGP
11152.2
0.6746
22.20
46.54
0.1784
0.10
-2583.1
1.0000
142.72
singular under our experimental configuration. sMTGP outperforms MVGP because it considers modeling each mode of the tensor output, and MVGP ignores the tensor structure by vectorizing it.
Table summarizzes the optimization performance of different methods for the BO and CBBO tasks. We set N 10d and evaluate each method using MSE,,, MAE,y and Acc, defined in Appendix K For each surrogate, the UCB-based strategy consistently outperforms random sampling. This is
Table 3: The optimization performance of different methods in the three synthetic settings.
Setting (l)
Setting (2)
Setting (3)
MSE,
MAE,
Acc
MSE,
MAE,
Acc
MSE,
MAE,
Acc
TOBO
0,0000
0.0008
0.0003
0.0350
0.0001
0.0050
SMTGP- UCB
0.0001
0.0031
0.0003
0.0361
0.048
0.0590
MLGP-UCB
0,0433
0,3793
0.0512
0.9295
0.0342
0.6263
MVGP-UCB
0.0015
0.0523
0.0044
0,0550
0.0342
0.6263
BO
TOGP-RS
0.0893
0.3145
0.0026
0.0351
0.0106
0.1526
MMGG-RS
0.0251
0,3242
0.0206
0,3684
0.0084
0.1223
MLGP-RS
0,0433
0.3793
0.0435
0.7976
0.0075
0.0934
MVGP-RS
0.0148
0,2036
0.0157
0.2697
0.0084
0.1223
TOCBBO
0,0023
0.0172
1,00
0,0000
0.0000
1.00
0.0021
0.0145
1.00
<MTG--UCB
0,1832
0.5614
0.67
0,0000
0.0000
1.00
0.0075
0,2171
0.86
MLGP-UCB
0.0667
0.6527
0.33
0.1826
0.0779
1.00
0.2070
0.7105
0.43
CBBO
MVGP-UCB
0.0032
0.0285
1.00
0.0725
0.0312
1.00
0.0151
0.1988
0.71
TOGP-RS
0,0438
0,5319
0.67
0.0908
0.0395
1.00
0.0512
0.8793
0.43
SMGG-RRS
0,3151
0,5882
0.67
0.1826
0.0779
1.00
0.0117
0.9453
0.29
MLGP-RS 0.1053
0.6909
0.33
0.1826
0,0779
1.00
0.0117
0.9453
0.29
MVGP-RS
0.1313
0.5975
0.33
0.1489
0.0654
0.00
0.0512
0.8793
0.43
intuitive due to UCB's better theoretical guarantees. Across all GPs, TOBO and TOCBBO achieve the smallest MSE2 and MAE,, indicating that the selected pairs yield objective values closest the true optimum. Among the baseline methods, sMTGP-UCB delivers the second-best performance, followed by MVGP-UCB, while MI.GP-CBB performs the worst. This result is consistent with their modeling abilities shown in Table
Finally, we provide each round's logarithmic instantaneous regret for different methods for BO and CBBO in Figure I
We can observe that TOBO and TOCBBO consistently achieve the lowest
TOUPHO
(01g-42 MCPP-22
11
2
MTTE UCE
MESPICE
MEPREE
Figure I: Each round's logarithmic instantancous regret of different methods in the Setting (1) (L), (2) (M), and (3) (R) for BO (Top row) and CBBO (Bottom row).
instantaneous regret across all three settings, highlighting their superiority. Additional synthetic results are provided in Appendix M

---

Published as a conference paper at ICLR 2026
5.2 CASE STUDIES
We further apply apply TOBO and TOCBBO to four real-world datasets: (1) CHEM (Shields et al,
2021 input x: IR_ and output y E IR 14x3x3, (2) MAT (Wang al. 2020 input æ IT and
output y E Rx5xXx1, (3) PRINT (Zhai et al.
2023 y input œ: E R and output y E R (4)
REEN: input æ E R and output y R R Detailed description of these datasets is provided in Appendix K
Since the REEN dataset provides fully observed tensor outputs, we first evaluate surrogate modeling performance by randomly selecting 30 samples for training and 5 samples for testing. Table 4
reports the predictive performance of the four GP surrogates on the test set. TOGP achieves the best accuracy in terms of NLL and MAE. Figure 2 further compares optimization performance on both BO (left) and CBBO (right), where TOBO and TOCBBO consistently perform the best, demonstrating their effectiveness real-world black-box systems.
Table 4: The prediction performance of GPs in the REEN dataset.
TOGP
sMTGP
MI.GP
MVGP
NLL
15.6664
33.3198
88.2722
48.7167
MAE
0.0883
0.0993
0.0929
0,1054
[Cov]]
0.4918
0.3555
0.4318
0.5711
Objectovve vallue
Objectiovo value
3.7
9
3.4
3.1
10
15
20
25
30
5
10
15
20
25
30
Sequential samples
Sequential samples
Figure 2: Each round's optimal objective value REEN for BO (L) and CBBO (R).
The remaining three datasets contain only partially observed tensor outputs, we only evaluate op- timization performance under the CBBO setting. Figure 2 shows that TOCBBO can identify the optimal input consistently using fewer rounds than the baselines, further demonstrating its superior effectiveness.
TOCMO
MONAEN
MOUSE
Figure 3: Each round's optimal objective value in CHEM (L), MAT (M), and PRINT (R).
6 CONCLUSION AND DISCUSSION
In this work, we develop two BO methods for tensor-output systems. TOBO employs two classes of kernel-based TOGP surrogate model and selects query points via UCB acquisition function. TOCBBO extends TOGP to the partially observed tensor outputs and adopts a CMAB-UCB2 cri- terion to sequentially select both the query input and the super-arm. We establish theoretical regret bounds for both methods and demonstrate their effectiveness on extensive synthetic and real-world experiments. Future work could consider integrating the proposed tensor-output kernels with sparse techniques, such as sparse GPs (Snclson & Ghahramani] 2005) and scalable LMC (Bruinsma et al.) 2020), to improve the computational efficiency of TOGP. would be valuable to explore alterna- tive acquisition functions within this framework. For example, one may combine the TOGP model with improvement-baat based criteria such as El or PI (Frazier) 2018), and further extend them to the TOCBBO setting. Moreover, acquisition functions based on improvement (Uhrenholt & Jensen, 2019) or coormation-theoreee criteria (Tu et al. (022)) are also be considered in our framework, provided that the computational challenges for tensor outputs can be addressed. Finally, it is worth exploring more meaningful tensor structures for our proposed framework, such as a spatiotemporal system.
10

---

Published as a conference paper at ICLR 2026
ACKNOWLEDGMENTS
This work was supported by National Natural Science Foundation of China (Grant 200.72271138 and Ssinghua-Nationa University of Singapore Joint Funding (Grant No. No.20243080039, We grate- fully acknowledge this support.
REFERENCES
Karim
Nguyen Linh Trung, Adel Hafiane, et al. A contemporary and comprehensive survey on streaming tensor decomposition. IEEE Transactions Knowledge and Data Engi- neering,
2022.
Guglielmo Maria Accabi, Francesco Trovo, Alessandro Nuara, Nicola Gatti, and Marcello Restelli. When gaussian processes meet combinatorial bandits: Gcb. In 14th European Workshop Re- inforcement Learning, pp. 1-11, 2018.
Peter Auer, Nicolo Cesa-Bianchi, Yoav Freund, and Robert E Schapire. The nonstochastic multi- armed bandit problem. SIAM journal on computing, 32(1):48-77, 2002.
Syrine Belakaria, Aryan Deshwal, and Janardhan Rao Doppa. Max-value entropy search for multi- objective bayesian optimization. Advances in neural information processing systems, 32, 2019.
Mikhail Belyaev, Evgeny Burnaev, and Yermek Kapushev. Gaussian process regression for struc- tured data sets. In International Symposium on Statistical Learning and Data Sciences, pp. 106- 115. Springer, 2015.
Léon Bottou. Large-scale machine learning with stochastic gradient descent. Proc. of COMPSTAT, 01 2010.
Wessel Bruinsma, Eric Perim, William Tebbutt, Scott Hosking, Aro Solin, and Richard Turner. Scalable exact inference in multi-output gaussian processes. In International Conference Ma- chine Learning, pp. 1190-1201. PMLR, 2020.
Adam Bull. Convergence rates of efficient global optimization algorithms. Journal of Machine Learning Research, 12(10), 2011.
Zexun Chen, Bo Wang, and Alexander N Gorban. Multivariate gaussian and student-t process regression for multi multi-output prediction. Neural Computing and Applications, 32:3005-300 -3028, 2020.
Eric Chi and Tamara Kolda. On tensors, sparsity, and nonnegative factorizations. SIAM Journal on Matrix Analysis and Applications, 1772-1997, 2012.
Sayak Ray Chowdhury and Aditya Gopalan. No-regret algorithms for multi-task bayesian optimiza- tion. In International Conference on Artificial Intelligence and Statistics, pp. 1873-1881, PMLR,
2021.
Tinkle Chugh. Scalarizing functions in bayesian multiobjective optimization. In 2020 IEEE Congress on Evolutionary Computation (CEC), pp. 1-8. IEEE, 2020.
Sihui Dai, Jialin Song, and Yisong Yue. Multi-task bayesian optimization via gaussian process upper confidence bound. In ICML 2020 workshop on real world experiment design and active learning, volume 60, pp- 61, 2020.
Samuel Daulton, Maximilian Balandat, and Eytan Bakshy. Differentiable expected hypervolume improvement for parallel multi-objective bayesian optimization. Advances in Neural Information Processing Systems, 33:9851-9864, 2020.
Samuel Daulton, David Eriksson, Maximilian Balandar, and Eytan Bakshy. Multi-objective bayesian optimization over high-dimensional search spaces. In Uncertainty in Artificial Intelligence, pp- 507-517. PMLR, 2022a.
Samuel Daulton, Xingchen Wan, David Eriksson, Maximilian Balandat, Michael A Osborne, and Eytan Bakshy. Bayesian optimization over discrete and mixed spaces via probabilistic reparame- terization. Advances in Neural Information Processing Systems, 25:12760-12774 2022b.
11

---

Published as a conference paper at ICLR 2026
Vin De Silva and Lek-Heng Lim. Tensor rank and the ill ill-posedness of the best low-rank approxi- mation problem. SIAM Journal on Matrix Analysis and Applications, 3013::088-1117,, 2008.
Peter Frazier. A tutorial bayesian optimization. arXiv preprint
2018.
Thomas Fricker, Jeremy E Oakley, and Nathan M Urban. Multivariate gaussian process emulators with nonseparable covariance structures. Technometrics, 55(1):47-56, 2013.
José Henrique de M Goulart, Maxime Boizard, Rémy Boyer, Gérard Favier, and Pierre Comon. Tensor cp decomposition with structured factor matrices: Algorithms and performance. IEEE Journal of Selected Topics in Signal Processing, 10(4):757-769, 2015.
Daniel Hernández-Lobaad Jose Hernandez-Lobato Amar Shah, and Ryan Adams. Predictive en- tropy search for multi-objective bayesian optimization. In International conference on machine learning. pp. 1492-1501. PMLR, 2016.
José Miguel Hernández-Lobuue Matthew W Hoffman, and Zoubin Ghahramani. Predictive entropy search for efficient global optimization of black-box functions. Advances in neural information processing systems, 27, 2014.
Dionissios T Hristopulos. Non-separable covariance kernels for spatiotemporal gaussian processes based on hybrid spectral method and the harmonic oscillator. IEEE Transactions on Information Theory, 70(2):126) (122):226812221 1283, 2023.
Congfang Huang, Jaesung Lee, Yang Zhang, Shiyu Zhou, and Jiong Tang. Mixed-input bayesian optimization method for structural damage diagnosis. IEEE Transactions on Reliability, 72(2): 678-691, 2022.
Seyed Mostafa Kia, Christian F Beekmann, and Andre F Marquand. Scalable multi-task gaussian process tensor regression for normative modeling of structured variation in neuroimaging data. arXiv preprint 1808.00 000800080000,,, 2018.
Yongxiang Li and Qiang Zhou. Pairwise meta-modeling of multivariate output computer models using nonseparable covariance function. Technometrics, 58(4):483 494, 2016.
Eric F Lock. Tensor-on-tensoo regression. Journal of Computational and Graphical Statistics, 27 (3):638-647, 2018.
Wesley Maddox, Maximilian Balandat, Andrew Wilson, and Eytan Bakshy. Bayesian opti- mization with high-dimennsonal outputs. Advances in neural information processing systems, 34: 19274-19287, 2021.
Dang Nguyen, Sunil Gupta, Santu Rana, Alistair Shilton, and Svetha Venkatesh. Bayesian opti- mization for categorical and category-speeific continuous inputs. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 34, pp. 5256-5263, 2020.
Ivan V Oseledets. Tensor-train decomposition. SIAM Journal on Scientific Computing. 33(5):2295- 2011.
Omiros Papaspiliopoulos. High-dimensional probability: An introduction with applications in data science, 2020.
Binxin Ru, Absan Alvi, Nguyen, Michael A Osborne, and Stephen Roberts. Bayesian optimi- sation over multiple continuous and categorical inputs. In International Conference Machine Learning, pp. 8276-88 8285. PMLR, 2020.
Thomas J Santner, Brian J Williams, William I Notz, and Brain Williams. The design and analysis of computer experiments, volume 1. Springer, 2003.
Pier Giuseppe Sessa, Pierre Laforgue, Nicolò Cesa-Bianchi and Andreas Krause. Multitask leamn ing with no regret: from improved confidence bounds to active learning. Advances in Neural Information Processing Systems, 36:6770-6781, 2023.
12

---

Published as a conference paper at ICLR 2026
Benjamin J Shields, Jason Stevens, Jun Li, Marvin Parasram, Farhan Damani, Jesus I Martinez Alvarado, Jacob M Janey, Ryan Adams, and Abigail G Doyle. Bayesian reaction optimization as tool for chemical synthesis. Nature, 590(7844):89 2021.
Edward Snelson and 7oubin Ghahramani. Sparse gaussian processes using pseudo-inputs. Advances in neural information processing systems, 18, 2005.
Jasper Snoek, Hugo Larochelle, and Ryan P Adams. Practical bayesian optimization of machine learning algorithms. Advances in neural information processing systems, 25, 2012.
Lei Song, Ke Xue, Xiaobin Huang, and Chao Qian. Monte carlo tree search based variable selection for high frazier201 Stutorialdimensional bayesian optimization. Advances in Neural Information Processing Systems, 35:28488 28501, 2022.
Qingquan Song, Hancheng Ge, James Caverlee, and Xia Hu. Tensor completion algorithms in big data analytics. ACM Transuctions on Knowledge Discovery from Data (TKDD), 13(1):1-48, 2019.
Niranjan Srinivas, Andreas Krause, Sham M Kakade, and Matthias Seeger. cess optimization in the bandit setting: No regret and experimental design. arXiv:0912.3995, 2009.
Gaussian pro- arXiv preprint
Ben Tu, Axel Gandy, Nikolas Kantas, and Behrang Shafei. Joint entropy search for multi-objeccive bayesian optimization. Advances in Neural Information Processing Systems, 35:9922-9938, 2022.
Anders Kirk Uhrenholt and Bjoem Sand Jensen. Efficient bayesian optimization for target vector estimation. The 22nd International Conference on Artificial Intelligence and Statistics, pp- 2661-26700 PMLR, 2019.
Boqian Wang, Jiacheng Cai, Chuangui Liu, Jian Yang, and Xianting Ding. Harnessing novel machine-learing-asssstee evolutionary algorithm to co-optimize three characteristics of an elec- trospun oil sorbent. ACS Applied Materials & Interfaces, 12(38):42442-4284 2020.
Xilu Wang, Yaochu Jin, Sebastian Schmitt, and Markus Olhofer. Recent advances in bayesian optimization. ACM Computing Surveys, 55(13s):1 1-6, 2023.
Yan Wang, Meng Wang, Areej AlBahar, and Xiaowei Yue. Nested bayesian optimization for com- puter experiments. IEEE/ASME Transactions on Mechatronics, 28(1):440444,, 2022.
Christopher Williams and Matthias Seeger. Using the nyström method to speed up kernel machines.
Advances in neural information processing systems, 13, 2000.
Jian Wu, Matthias Poloczek, Andrew G Wilson, and Peter Frazier. Bayesian optimization with gradients. Advances in neural information processing systems, 30, 2017.
Rose Yu, Guangyu Li, and Yan Liu. Tensor regression meets gaussian processes. In International conference on artificial intelligence and statistics, pp. 482-490. PMLR, 2018.
Cuihong Zhai, Jianjun Wang, Yiliu Paul Tu, Gang Chang, Xiaolei Ren, and Chunfeng Ding. Ro- bust optimization of 3d printing process parameters considering process stability and production efficiency. Additive Manufacturing, 71:103588, 2023.
Yunchuan Zhang, Sangwoo Park, and Osvaldo Simeone. Multi-fidelity bayesian optimization with taso transferable max-value entropy search. IEEE Transactions on Signal Processing, 2025.
Shandian Zhe, Wei Xing, and Robert M Kirby. Scalable high-order gaussian process regression. In The 22nd International Conference on Artificial Intelligence and Statistics, pp. 261 1-2620. PMLR. 2019.
13

---

Published as a conference paper at ICLR 2026
A USE OF LLMS
Large Language Models (LLMs) were used to aid the writing and polishing of the manuscript. Specifically, we used an LLM to assist in refining the language, improving readability, and ensuring clarity in various sections of the paper. The model helped with tasks such as sentence rephrasing. grammar checking, and enhancing the overall flow of the text.
is important to note that the LLM was not involved in the ideation, research methodology. or experimental design. All research concepts, ideas, and analyses were developed and conducted by the authors. The contributions of the LLM were solely focused on improving the linguistic quality of the paper, with involvement in the scientific content data analysis.
The authors take full responsibility for the content of the manuscript, including any text generated or polished by the LLM. We have ensured that the LLM-generated text adheres to ethical guidelines and does not contribute to plagiarism or scientific misconduct.
B THE DETAILS OF THE TENSOR DECOMPOSITION IN REMARK 1
As discussed in Remark I directly estimating the entries of the core tensors {A or A re- quires mZ' parameters for the non-separable tensor-outpus kernel and T parameters for the separa- ble tensor-output kernel, which becomes intractable for large-scale or high-order tensors. To reduce the parameter complexity, we impose low-rank tensor structures and adopt different decompositions depending on the tensor order.
Low-order tensors (m S 3):
When the tensor order is small, we employ the CP decomposition:
Ri
A
Oirl airy ...... atrm for the non-separable kernel,
(17)
r=1
R
A
0,1 a,2 O
arm
for the separable kernel,
(18)
where airi, Oyi R for number of free parameters becomes Σm Σn= Rit, for the non-separable kernel and R > t. for the separable kernel, which scales linearly with each mode size t, Thus, CP provides a compact representation when m S 3. However, CP can be ill-posed for higher-order tensors because the set of tensors with bounded CP rank is not closed, and a best low-rank approximation may fail to exist (De Silva & Lim 2008). Moreover, the factor matrices can easily become ill-conditioned as m increases, leading to numerical instability (Chi Kolda 2012).
R and 1,..., for 1,.... R and i 1, 1,..., ............ The
m.
High-order tensors (m > 3): When the tensor order is large, we employ the TT decomposition:
A = Gn Gre(tt)) ... Gim(tm) for the non-separable kernel, A G2 ... G. (fm) for the separable kernel,
(19)
(20)
where each Gt,(t,) as X x XTI, three-mode tensor and the T'-ranks satisfy 11,0 = fl.m
for = 1,..., m. Similarly, (t is a-j x tú three-mode tensor that satisfies To I'm 1. The total number of free parameters is Ei Eg F1.j-165T1, for the non-separable tensor- output kernel in 6 and Σ ryyi for the separable tensor-output kernel in 7). It scales linearly with the tensor order m, instead of exponentially as in the full tensor. This makes TT decomposition highly suitable for high-order tensors (m >> 3), as it balances modeling flexibility with computational scalability and avoids the instability of CP.
In summary, we use CP for low-order cores and TT for high-order cores, ensuring both efficiency and numerical robustness across different tensor settings.
C THE ESTIMATION OF HYPERPARAMETER FOR THE TOGP AND PTOGP
In this appendix, we provide the details of hyperparameter estimation for both TOGP and PTOGP. Without loss of generality, we assume zero prior mean μ 0 in (3 and 12).
14

---

Published as a conference paper at ICLR 2026
C.I PARAMETER ESTIMATION FOR TRAINING TOGP
Given the training data X,, (#1,.... Fnd and Yn (vee((s)"....,.ee((yy) ..... the log marginal likelihood of TOGP is given by:
1 log En
log L(
2 Y E
(21)
where Σ, = K", I Then, 21 can be optimized by applying gradient- based optimization methods, such as L-BFGS algorithm.
The gradients of the log-likelihood function in 21 with respect to the hyperparameters 72, o², 0, and a is given by
i log L
lr E Kn.)
Y,
Dn
Kn YY.
(22)
i oog I 1
1
tr (m,11)
Y
Σ Σ YYn. ''n',,,
(23)
dg²
2
Olog L
72
72
OK n Σ, YY
Σm
Y Σ
(24)
00n
00is
2 n
Uitss
Olog tr
OK,
2 Y,
Σー UK
Σn
En YY.
(25)
dayz
days
2
Daus
where Big represents the scale parameters of ktj. the kernel function associated with the 1-th mode.
The matrices ok. and @K_ are the partial derivatives of the kernel matrix with respect to the
00.; Jany
corresponding kernel parameters.
The detailed algorithm for training TOGP is given follows:
Algorithm Parameter estimation for training TOGP
Input: Training data X,, and Y,. initial hyperparameters O, = {00) T0: 0 a
Initialize: o² <- oi 7² < Th. <- 0g 4 an;
1: while 7², o², 0, not converge do
2: Update 72 based on 22
3: Update o² based on 23
4: Update based on 24)
5: Update based 25
6: end while
Remark 4. For the non-separable tensor-output kernel in Definition 7 the total number of hy- perparameters to be estimated in TOGP is ma
2 + T + Eis Ens Rite when S 3 and 11.j-115715 when >> 3. For the separable tensor-outpul kernel in Def- = 2+ T +I= EØ=
inition 2 the total number of hypepprrameters be estimated in TOGP is my 2 R Ei 4 when < 3 and my. 2 T + Tー-tt,, when m > 3. The computational complexity of computing the gradient of log L(@) with respect to all my parameters O + 2 When using the L-BFGS algorithm to optimize the likelihood function results, the number of itera- tions typically scales as O (log(n) (Bottou 2010). Therefore, the overall computational complexity for training the TOGP takes O (og()) + n2 my log(n)) computational complexity.
C.2 PARAMETER ESTIMATION FOR TRAINING PTOGP
Given the partially observed training data X,, = (x1.... ... In), An (Ai...., An) and Y, (i1...., )n)) the log marginal likelihood of PTOGP is given by:
log L(0) log Σn ()ec(Yn)) "
(26)
where Σ, = E,E,KK,E, + I Then, 26) can also be optimized by applying gradient- based optimization methods.
15

---

Published as a conference paper at ICLR 2026
The gradients of the log- se-likelihood function in (26 with respect to the hyperparameters o². T², 0, and a given by
oi
Ir
K,
Y,
"
E,, E, nn
(27)
θr
OL
1 tr
_Yn En
1
T
E YYn
(28)
OG² 2
2
OL ?2
OKn
7² Y, E1-LE, En 00165
Okn E,, D,n 1,-Y
001j
tr
n.
En a0ts
En
2
(29)
aL 72
akn E
72 Y Σn En Dasi
Okn E, E
Σx
En
(30)
Janz
Dass
The detailed algorithm for training PTOGP is given as follows:
Algorithm 2 Parameter estimation for training PTOGP
Input: Training data X,,, A, and Y,, initial hyperparameters O, = (o0 T0 ao)
Initialize: o² <- 0ū. 7² +- Tō. <- 0o f- an:
1: while o², 72, 0, not converge do
2: Update 72 based on (7
3: Update o² based on 28
4: Update based on 229
5: Update based 30
6: end while
Remark 5. For the non-separable tensor-oulput kernel in Desinition 7 the total number of hy- perpurameters to estimated PTOGP ma 2 + Ei= E Rat; when m S and m. - - nition 2 the total number of hyperparameters to be estimated in PTOGP is m. = + R Eimh L when m < and m, = 2
when m For the separable tensor-output kernel in Defi-
-ーー when m 3. The computational complexity of computing the gradient of (()) with respect all my parameters is OKKTT + n n²=TTm)). When using the L-BFGS algorithm to optimize the likelihood function results, the number of iterations typically scales as (log(n)) Bottou (2010). Therefore, the overall computational complexity for training the PTOGP takes O (n' (oollol + n² log(n)) computational complexity.
D THE PROPOSED ALGORITHMS AND COMPUTATIONAL COMPLEXITY ANALYSIS
D.1 UCB-BASE TOBO ALGORITHM
The detailed procedure of the proposed TOBO method is given in Algorithm
Algorithm UCB JCB-based TOBO
Input: Total rounds N, initial dataset D, Ø, initial hyperparameters Ou; for round n 1..... N do
2: Update the posterior mean (4) and covariance 5 of f given On-
3: Select the next input æn < arg max our ( DD.-));
4: Evaluate the black-box box system and observe output yx;
5: Update dataset D,, + 1 U (In, Yn):
6: Update hyperparameters θ,, by maximizing
21 with L-BFGS;
7: end for
8: Identify i arg max]e[N] f(x:;:
9: Output: Optimal input 2;, and corresponding output y
16

---

Published as a conference paper at ICLR 2026
Remark 6. When using the TOBO method select x', the computational complexity of updating TOGP is ((n 1) at round n. Then, the computational complexity of querying the next point is O (m=1))122og(1))) by using the L-BFGS method. After updating the design dataset, the com- pulational complexity of updating the hyperparameters @ is O (n log(n) + mm log(n)) Thus, the computational complexity of TOOO at round is O (n log(n) T2 log(n) Suppose that there are N points needed to query, the total computational complexity of Algorithm [n log(n) + mm log(n))):
O
D.2 CMAB-UCB2-BASED TОСBBO ALGORITIM
The detailed algorithm of TOCBBO is provided in Algorithm
Algorithm 4 CMAB-UCB2-based TOCBBO
Input: Total rounds N, initial dataset Du 2, initial hyperparameters Ou:
1: for round n = ......N .N' do
2: Update the posterior mean (33) and covariance (14):
3: Select the next input on using 15)
4: Select the super-arm X2 using
5:
Evaluate under (In, An. and observe yn:
6: Update dataset D,, + Dn U {(In, An), 1:
7: Update hyperparameters On by maximizing log [(O) with L-BFGS;
8: Update incumbent solution {æ; A arg max;=1,..,, H (zx, A;);
9: end for
Identify = -rr max Ha, A:):
11: Output: Optimal input Vi., optimal super-arm As., and output y
Remark 7.
When using the proposed TOCBBO method to jointly select x² and X the compu- Then, based on the
tational complexity of updating the PTOGP at round n is O ( proposed CMAB-UCG UCB2 criterion, the computational complexity of querying the next input by us- ing L-BFGS method is O ((n-1)3k3 and selecting the next super-arm by using greedy Top-k method is OKTT)), respectively. After updating the current design dataset, the computational complexity of updating hyperparameters @ is O(n log(n) + nn mn log(n)). Therefore, at round n, the computational complexity of TOCBBO method is Onn*k3 log(n) + n KTTTm log(n) + KT Assuming a total of N rounds, the overall computational complexity of Algorithm
is
O
IN m=1 k log(n) nTTTm, log(n) + KTT)).
E THE PROOF OF PROPOSITION
I
Proof. We prove symmetry and positive semi-definite for the two kernel classes in Definitions 2
Non-separable kernel (Definition The tensor-output kernel is
Denote a: : voc(Ar) R (resp. a = voc(A) E R?).
K(x, x') = >>
a,a (z)(0,
x' E X,
(31)
8=13=1
where each kij is scalar positive semi-definite kernel.
For any x, x' E X, we have
K(æ,, T²)
(a k2)(2x, [ araa ke)(
Σ
aray kez (1)(2((x)) = K(x', x))
E
17

---

Published as a conference paper at ICLR 2026
This shows that the full tensor-output kernel K is symmetric.
Let {a>)i=1 C. X and {y:}1= C. IR be arbitrary. Then, we have
yi K( myyu
Σ
E
ye
(ara?
) ( m)
Ej
Σ
Σ
8ej.s (2 (82) 8e3.5
(32)
where Stj. = al yss E R. For each fixed (e,j), the matrix [Key +y)]" is positive semi-
definite, so Ec, Ség,i kez (z,, stj.j 0. Summing over (4,j) preserves nonnegativity, hence the Gram matrix induced by K is positive semi-definite.
Separable kernel (Definition
The kernel is
K (r. σ²) = a k(x,aa),
a, x' E X,
(33)
x)) where k is scalar positive semi-definite kernel. Since (aa (aaT)" aa and k(x, x²) k( have K(x,x') Eig (a T y:) k(a,, x;, (a y,) 2 0, because the Gram matrix [K( æ;)] is positive semi-definite.
K(x', z). For arbitrary {a, H and (y.)=- Ei.= y; K(2,; 2 )ı
Therefore, both kernel classes are symmetric and generate positive semi-definite Gram matrices for any finite set of inputs, i.e., they are valid tensor-output kernels on X.
F THE PROOF OF LEMMA
Proof. First, for tensor-output system vx)(f(x:) following a TOGP defined in
denote the
derivative field of v=)( to the j-coordinate element of a as gj (œ) == Övee(f(a
/dr where
(e)( Rd and {I.....d}. According to the derivative property of GP
Santner al.
2003 we have that (æ) E is multivariate-output GP with mean pa (xz) and covariance K.Y (xz 2 Cov(g;()), gj (x))) Øzy du o ("K(xx"")
day Pn pn(())
g²
At round n + 1. the observed data is denoted as X,, and Y,, then we have the posterior distribution of (x( is a T-dimensional Gaussian with mean and covariance
P (zz) = 0f)(()) Oz, (x,, OKK, (x.
(34)
dx, dz',
It is easy to verify that K. (Iy, 2ty) is positive semi-definite for every For any fixed x X. given X, and Yn. the random vector
( - p.y (z,) ~ N/O. K. (2, a;))).
(35)
Applying the Gaussian Lipschitz concentration (Proposition 2.5.2 and Theorem 5.2.2 in paspiliopoulos 22200) to the norm
Pa-
l|2 yields, for all 0,
Pr (1)(()) A. (2,)) 2 ( + +1) - exp
t²
(36)
2 X (xx,)
where (n) (æ; is the largest eigenvalue of K (æ, x)
Let DM (21,...,Ial Inr) C X' be any finite discretization. Using with t(x)
V (m) (xr)) log(M/5), and applying the union bound, we obtain with probability at least - 8,
(m) (z)()) n. (2)))) s Vur(K + Ambx(())) log(M/6),
Vo E Da. (37)
18

---

Published as a conference paper at ICLR 2026
Let CF == x / /tn( [ (z)) and Ao =- supzex Amaz (1z)) then (37) implies that, with
probability least 1 6,
[z) 1 || (z3)) + Co + V2AV log(M/6),
A x E DM-
(38)
Assume the kernel is sufficiently smooth so that g,(+) has almost surely Lipschitz sample paths.
Then there exist absolute constants a, b > 0 (depending on the Lipschitz modulus and covering- number bound of X) such that for any L' O,
exp (I.2
Pr sup MEX
B x))|
I' 1'+Cx) S
This follows by combining the net bound (38) with standard chaining argument pass from finite net to the full domain; the Gaussian tail is preserved with possible adjustment of absolute constants into a, b.
G THE PROOF OF THEOREM
Proof. The proof consists of two main parts. We first establish concentration inequality for (e)) for any æ E X, and then use this result to derive an upper bound for the regret.
Part Concentration inequality. We begin by proving concentration inequality for vec(f(x) evaluation on discrete set of points in the domain X. We then extend the result to neighborhoods of these discrete points, and ultimately to the entire space X.
We first prove basic concentration inequality for general T-dimensional Gaussian distribution, i.e. [u, K). If define U distribution V (O, IT), we can obtain
Z la and set another T-dimensional standard Gaussian
U L K2/2V
| |K//
Define the function f(U) |K / For any U,U', U', we have
((U))f(U))) - < | /2 (U - U'Il < K lop IU - U' = / Amax (K)
UU
where Amax (K) is the maximum eigenvalue of K. Therefore is Lipschitz function with constant Amax((K)
According to Proposition 2.5.2 and Theorem 5.2.2 in Papaspiliopoulos 2020). we have the follow-
ing concentration inequality for Lipschitz function of standard Gaussian distribution:
Pr (F(U) [[(() + s exp (-
212
Substituting f and L, obtain
t2 22max (K)
Pr ( 11l 2 E 1|ll + s exp
Since
E Z ||]] S / Z 11|22] vtr(K),
we get the final result for general Gaussian distribution:
Pr Z - اl| 2 (tr(K) + +)) < exp( (2)mms(())
(40)
Define the discrete set DM = { AM} C X. According to the above concentration inequality, have
Pr ([]voo(f(2)) (x)(1))
+ +z) S exp(
2² Va D
)
( (m))
19

---

Published as a conference paper at ICLR 2026
where [n] (a: is the maximum eigenvalue of (x)
By setting = = v (in (2)) log $84 we obtain
Pr [vec(f(x)) ( S Bn (x)
21-6,
Vx E DM,
(41)
where Bn V
Suppex
(x, X
V 2 22og M o/4
vCn /2log o/4
From Lemma we obtain that ( vec(f(ar))) S (I/ + (8) æ'1h holds with probability at least 1 - a e (12) for all x, x' EX.
Then, at round
the size of D Dan(n) as (Tn)", i.e., M(n) (T For (Min), have || - [ lh S nd where [ar],, the closest point in Dar to a.
Tu
Using the above aligns, we set L' -by log o. we obtain
[e)) - vec(f((==))))]] & (b log da + +cv)I | - 11 S (bylog de + cv) rd
with probability at least - 5 for all x E X.
Choosing Tn = rdn (6 log da 8 + CCV,, we have
Pr ([]vec((())) - vee vec(f((z)) S #/2 2 - 6,
Vo E X'
Combining the above results, obtain
Pr Pr((Ivo(f(az*) - p((==))))1 s 8n V (n) ((xx((2=))m) + n²
2 -
(b (lol///lg///5/ (6) -)-00))
where [x*], is the closest point in DM(()) to x", and Bn, vCn + 2o
Part 2. Regret bound. According to the Lipschitz property of h, we have Tn = h(xx)) - h(xx) LiIvec(f((x)) - vec(f(x..))I
(42)
where L > 0 is the Lipschitz constant.
Since
(xx) + Br V (n- (Ix)) (n=1((x01m) Bn 1 A ((x*)))) 2 vec(f(x*))
n
we have
In S I. 2,3 V X(m=x)) (1n-))((z()))) (Ixn) + nx
(43)
We first consider the first term:
4.2 Mix (x=1)) (xx,) (xx,) s 4:3 (1 (17 (xn=))(xx)) S 483 log (1 17 n1 (in= (xx))
< 4333//C log IT + +1 -'K,, 1 (xx, 2x,)1
where C2 is constant.
Define C = 2m/ Then
4,82((-))) S C 4,2 n max
Σ log Lr + 1) 1nn11(1,, an)) S (K,,))
where the last inequality holds by the definition of YN (K.
Since ΣN + S 6
we have the final result:
N
Σ On < L LLVVCNNNN(K.)) + +
(44)
7=
20

---

Published as a conference paper at ICLR 2026
H THE PROOF OF PROPOSITION 2
Proof. We start from the definition of the (maximum) information gain:
In (K, n) max 2 logdel del (InT + +1 1 K,), K, [K RRu
(45)
Under the separable kernel in Definition 2 we have
K (c, a²) = (vec( A Avee(A) l:x x²) =: B k(a,, 2');
(46)
where B
A) vee(A))ec(AA A) RTX.
Then, the Gram matrix factorizes Kronecker product
K K & B, K2 [k(a; my)]k E R
(47)
Let {a;) be the eigenvalues of K2 and {B,) be the eigenvalues B. By the spectral property of the Kronecker product, the eigenvalues of K B are {ay 1..... n, j 1........ D). Therefore, we have
T
log det +1 '0 2 B)) = Σ Σ
log(1 (1 + ca, Bg),
c= c==
(48)
For each fixed the function log(1 C.2 u) is non-decreasing for 2 0, hence
log (1 + ca;
Σ log(1 + ca. 95) s logo + c:2 Bmax),
(49)
3:0.30
where Bmax max; Bj- Summing over 1, ,....,, and multiplying by 1/2 gives
T
(, S
log 1 + c Ai) T Tm
(50)
where kt Punax & is simply rescaled version of k. Since rescaling by a positive constant does not change the asymptotic order of the information gain, we obtain the general comparison bound
7m (,n) n) S T 7 (n.kk, (k;
In our specific setting, the tightest bound is nn (K, 77) O (ks n)). Finally, we substitute known results for the scalar kernel k: (i) If & is the Gaussian (squared exponential) kernel in d dimensions, then (k, 7)) = (log n) dd+1, which gives 7n (K, 17) = (T d)d+1).
(ii) If l a Matérn kernel with smoothness parameter v > 1, then 7n (k, n) = on On log n
which gives 7n (K, 11) oT 2o at log
dtd 1)
I THE PROOF OF THEOREM
Denote h(æ, A) = H f(x A). then we have
h(x", (()= h(z An) = [1( X²) - h(xn, x-) + [[ X²) h(2 (m))
= Tin + T2n-
For the first item, according the Cauchy-Sehwarz inequality, we have
H H A) im (x, X) < H H1 A) (n A)-L-
(51)
Similar to the proof of Theorem 1 provide the following lemma.
21

---

Yn (K, n) max Published as a conference paper at ICLR 2026
Lemma 2. Under Assumption
3 suppose the noise vectors {vec(E:) are independently and
identically distributed in [(0,0 σ IR). Then, for any å E (0, with probability at least 1 - 6, [f A) (x A) - pnl (x, x; X, A)l holds uniformly over all æ E X and X A and
i L, where Bn V 2d log
rddn htt///oggd///4CCCC 6), Ön sup
tr(K, (.
and
Ann Ag)
O7 supzex V tr tr(Kk (z))
From Lemma 2 and 51 we have
H Hyf(x, H (z, X) S H H1lf(a X) (x, S ( A,
Then, have
Fin h(x', X²) h(x X*)
S Hy (x", x²)) + Bml Kn i (x""2")111 - x*)
$ 2 K (2 .n A X \11/1//2
Similar to Theorem
we obtain
RIN Σ Tin < 1
ani 2,))111 1/
< H (K) BN + =
where C3 > 0 constant.
For the second item, followed in Accabi et al. 2018) we have
R2N S H V C.TN5, n)px,
where C4 > 0 a constant.
Then we have the cumulative regret over N rounds is bounded by
I' S Fin -2n
311 11 (VCoTN(K, BN = + VCTNNN (K. m)xx)
J
THE PROOF OF PROPOSITION
Proof. Based on the definition of the (maximum) information gain, we have
X, 2 log det (Im + cKa),
c: 103
(52)
where the Gram matrix takes the form K2 = EK, E and K, = B K., with K, (,,,,, R=xn
=
(ee((A)ee(A) TXX
By Sylvester's identity, we have
dct (Ink + C,E,K,,,,, det (I,T + cK!/EE, E,, K1/2).
(53)
Since E,, selects rows, F.J Frn InT, hence
log det (Ink + C K,,KK, E..
log det [nT + è K!/"E, E, E. K2/2) log det [nr cK,). 1 we have
(54)
Then, followed in Proposition Nm (K,,, 11) = (T(log for the Gaussian kernel,
(55)
ا(باا) Ni [K,, 11) of I' नلक logn) for the Matérn kernl..
(56)
22

---

Published as a conference paper at ICLR 2026
From we also obtain that
in (Kn, n 11) 0 (log (logn) n) for the Gaussian kernel,
(57)
in (n,, 17) - of on''nnnnn (n log 12) for the Matérn kernel.
(58)
K DETAILED SETTINGS OF THE EXPERIMENTS
K.1 DETAILED SETTINGS OF SYNTHETIC EXPERIMENTS
The setting of simulations are as follows:
Compute resources: All experiments were run Windows 10 Pro (Build 19045) desk- top with an Intel Core 19-7900X CPU (3.10 GHz) and 32 GB RAM.
Kernel setting: We use the Matérn kernel function with the smoothing parameter 5/2 the input kernel for different GPs.
Data setting: For the GP prediction, generate Teraiin 10d training samples for timating hyperparameters and Nest = 5d testing samples for predicting. In the BO and CBBO framework, we generate no 5d initial design points based on Latin Hypercube Sampling (LHS), and N 10d sequential design points.
Criteria: To compare the GPs prediction performance of different methods, we use the following criteria:
(1) Negative LigeLikelihood (NLL): It measures how well a probabilistic model fits the observed data. For given data X,, and Y,, the NLL is defined as NLL log(>"K, + (r"K, + o I'T) -Yn
o² I'T Y
(2) Mean Absolute Error (MAE): MAE 1 Eis Yi=Yi
(3) Covariance Operator Norm (|CCoul):: |K,||| ii
K
For the BO framework, we use the the mean squared eтror of inputs (MSE.) and the absolute error of outputs MAE, compare the optimization performance of dif- ferent methods. Here MSE, = || 1y and MAE, - || over N rounds. For the CBBO framework, we add the Ace criterion compare the match between the optimal super arm and the super arm chosen over N rounds, that is, Acc = =Ix-=xx/ /k, where is indictor function.
K.2 DETAILED SETTINGS OF CASE STUDIES
The detailed datasets of studies are as follows:
Chemistry reaction (CHEM):: Reaction optimization is fundamental to synthetic chemistry, from improving yields in industrial processes to selecting reaction conditions for drug candidate synthesis. According to Shields et al. (2021), we aim to evaluate the reaction parameters (æą: concentration, Nz: temperature) improve the experimental yields (y 3 x 3) of palladium-catalysed direct arylation reaction under varying bases (Mode 1), ligands (Mode 2), and solvents (Mode 3).
(PS/PAN) materials are com- monly used as potential oil sorbents for marine oil spill remediation. From Wang et al. (2020), we aim optimize the fabrication parameters of PS/PAN materials, including spinneret speed (c) ), collector distance (z2), applied voltage (æ3), and fiber diameter (æ1), improve their oil absorption capacity (y : 5 x 4 x 4) under varying PS content (Mode 1), mass fraction (Mode 2), and SiO, content (Mode 3).
PS/PAN material (MAT): Electrospun polystreeneppoyyarylooiite
3D printing (PRINT): Material extrusion-bate based diree-dimeniional printed products have been widely used in aerospace, automotive, and other fields. Following Zhai et al. (202)), focus on selecting appropriate process parameters (æy: layer thickness, I2: platform temperature, æz: nozzle temperature, Z4: infill density, and Øs! printing speed) to reduce variations in part quality
23

---

Published as a conference paper at ICLR 2026
(y: x x 3) caused by different printer nozzles (Mode 1) and printing geometries (Mode 2). The quality (Mode 3) evaluated in terms of compression deformation, compressive strength, and the printing cost.
Renewable energy (REEN); Climate change affeets the availability and reliability of renewable energy sources such wind, solar, and hydropower. We employ the operational energy dataset from the Copernicus Climate Change Service (https://cds.c.) www.ate..mmmm.cccmmm...........e copernicus eu/) datasets/si en enezgy-derivec: httprrictorooarrassatttaaaeeeeee
explore the cli- mate conditions that are most beneficial to renewable energy generation in various European nations. The climate-related variables, used as imput features, include air temperature (æ1), precipitation (a2), surface incoming solar radiation (Iz), wind speed at 10 meters (zz)) and 100 meters (Iş), and mean sea level pressure (rg). The energy-related indicators (y 10 x 2) collected from ten European countries (Mode 1), used as outputs, include the capacity factor ratio of solar photovoltaic power generation and wind power generation onshore (Mode 2).
L THE COMPARISON OF COMPUTATIONAL COMPLEXITY FOR BASELINES
To provide clearer comparison, summarize the computational complexity of GP training, BO, and CBBO for the baseline methods sMTGP, MLGP, and MVGP, together with our proposed method, in the table below. Here mh, Mh1, Ma2, and maz for 1, m denote the numbers of
Table 5: The computational complexity of different methods for GP training, BO, and CBBO.
Task
Method
Computational complexity
TOGP
O + n logn)
sMTGP
((ni +
t
t + + + - am
mm log
GP Training
MLGP
((na
+
Σ
+ + + 2==1 ti fmm) log
mm
1=1
MVGP
((n)
7
+ +n² =n²ma + T2 log
TOGP
O + +n Tm2 ) log n
MMGG
(((n3
1=1 t + +n²ma + Σm tf log n
BO (Round n)
-
MI.G
((
21=1 ts + +n²m, + -i= 1m log
MVGP
0 (n 7 + +n²mAi + +T log
TOGP
0 (n mn=kTma log n KT
SMTGP
( + 13 + n2
Σm
1-1 tymms) log +
CBBO (Round n))
MLGP
(nn' + 13 + n² n²mn k
m
tum)s) log + k 1-1
MVGP O(n3 + T + n m + KT log KT
hyperparameters for the corresponding methods.
Furthermore, we report the empirical computational time for GP training across all methods under the three experimental settings. The results are summarized below.
Table 6: The runtime (s) of different methods for GP training in the three synthetic settings.
TOGP
sMTGP
MIGP
MVGP
Setting (1)
266.04
27.26
30.57
6,50
Setting (2)
69. 85
16.08
21.99
1.31
Setting (3) 900.13
851.34
891.73
769.50
We further report the empirical running time of BO and CBBO under all baseline methods, The results are summarized in the table below.
Finally, provide the prediction and optimization performance of TOGP based on the Nystr low-rank approximation (SVD-TOGP) in Remark 3 under the three synthetic settings. The results summarized below.
24

---

Published as a conference paper at ICLR 2026
Table 7: The runtime (s) of different methods for BO and CBBO in the three synthetic settings.
Setting (1)
Setting (2)
Setting (3)
Method/Task
BO
CBBO
BO
CBBO
BO
CBBO
TOBO/TOCBBO
6968.25
8017.08
1588.53
2348.16
6035.61
11699.09
sMTGP-UCB
1516.31
953.82
117.78
179.11
5291.26
9436.90
MLGP-UCB
1545.41
981.86
133.93
185.49
5452.90
9696.60
MVGP-UCB
340.44
437.29
14.86
20.76
3006.97
3199.47
TOGP-RS
6176.84
7990.44
1573.36
2338.97
5786.02
11501.01
sMTGP-RS
1503.24
939.27
106.93
164.95
5273.84
9421.08
MLGP-RS
1530.32
971.25
109.14
180.70
5252.32
9431.66
MVGP-RS
334.79
386.86
7.70
6.70
2997.75
3116.23
Table 8: Summary of TOGP and SVD-TOGP performance for GP training, BO and CBBO in the three synthetic settings.
Setting (1)
Setting (2)
Setting (3)
Task
Criterion
SVD-TOGP
TOGP
SVD-TOGP
TOGP
SVD-TOGP
TOGP
NLL
557.09
503.00
-18.39
-18.10
-3778.74
-3923.10
MAE
0.1756
0.1571
0.1099
0.1052
0.1436
0.1372
GP Training
CCol
3.9627
2.0200
0.6261
0,0400
0.0881
2.8200
Time
21.80
266.04
10.83
69.85
436.17
900.13
MSE,
0.0001
0.0000
0.0003
0,0003
0.0002
0.0001
MAE,
0.0041
0.0008
0.0356
0,0350
0.0051
0.0050
BO
Ins Regret
0.0040
0.0001
0.0002
0,0002
0.0402
0.0302
Time (8)
831.02
6968.25
162.78
1588.53
3604.66
6035.61
MSE,
0,0024
0.0023
0,0000
0.0000
0.0035
0,0021
MAE,
0,0180
0.0172
0,0000
0,0000
0.0246
0.0145
CBBO
Acc
1
1
1
Ins Regret Time (s)
1.0944
0.9807
0,0000
0.0000
1.8571
1.4406
655.83
8017.08
134.22
2348.16
4446.12
11699.09
M ADDITIONAL RESULTS OF SYNTHETIC EXPERIMENTS
Table 9: The additional optimization performance of different methods in three synthetic settings.
Setting (1)
Setting (2)
Setting (3)
k
MSE,
MAE,
Acc
MSE
MAE,
Acc
MSE,
MAE,
Acc
TОСBBO
0.0054
0.0509
1.00
0.0001
0.0036
1.00
0.0021
0.0312
1.00
SMTGP-UCB
0.0057
0.0942
1.00
0.0548
0.3336
1.00
0.0340
0.8616
0.71
MLGP-UCB
0.0677
0.9775
0.33
0.0434
0.6961
0.50
0.3549
0.8205
0.29
MVGP-UCB
0.0325
0.7530
0.33
0.0019
0.0191
1.00
0.0405
0.3647
0.71
T/3
TOGP-RS
0.0208
1.3279
0.83
0.0206
0.4347
0.50
0.0709
1.5521
0.57
SMTGP-RS
0.0203
1.2625
0.50
0.1201
0.4443
0.50
0.0254
1.0703
0.50
MIGGP-RS
0.0619
1.0296
0.67
0.1191
0.5081
1.00
0.1816
1.1735
0.64
MVGP-RS
0.0634
0.9733
0.67
0.0469
0.4444
0.50
0.0574
1.3837
0.43
TOCBBO
0.0006
0.0125
1.00
0.0010
0.0076
1.00
0.0103
0.1163
0.96
SMTGP- -UC
0.0066
0.0585
1.00
0.0012
0.0177
1.00
0.0122
1,9890
0.89
MLGP-UCB
0.0323
3.6747
0.82
0.2229
1.4532
0.75
0.0423
4.0294
0.81
MVGP-UCB
0.0047
0.4765
1.00
0.0022
0.0309
1.00
0.0324
1.9855
0.93
2T/3
TOGP-RS
0.0128
4.6512
0.82
0.0202
1.1872
1.00
0.0218
10.4919
0.70
SMTGP-RS
0.0128
4.6512
0.82
0.2229
1.4532
0.75
0.0110
6.0795
0.70
MLGP-RS
0.0323
3.6747
0.82
0.2229
1.4322
0.75
0.0423
4.0294
0.81
MVGP-RS
0.0128
4.6512
0.82
0.0439
0.8309
1.00
0.0423
4.0294
0.81
25

---

Published as a conference paper at ICLR 2026
TOMPRI
30
log Rugies
2.5
20
25
Seeq
Sequental
100930
MUCP-ACE
MYONUEs
Vig Rogiet
Sociee day
15
25
35
10
20
3
15
25
30
Sequartial sumppiss
Sequerbal parpphes
Sequarial
Figure 4: Each round's logarithmic instantaneous regret of different methods the Setting (1) (Left), (2) (Middle), and (3) (Right) when & = T/3 (Top row) and k 27/3 (Bottom row).
Table 10: Summary across CHEM, MAT, and PRINT datasets with 10 repetitions.
Method
CHEM
MAT
PRINT
TOBO
(0.00)
150.07 (0.00)
-25.33 (0.00)
sMTGP- UCB
99.91 (0.10)
147.97 (2.24)
-25.03 (0.00)
MLGP-UCB
99.29 (0.75)
147.90 (2.32)
(0.11)
MVGP-UCB
96.32 (4.89)
147.20 (3.44)
-25.27 (0.32)
TOGP-RS
93.84 (1.56)
147.83 (2.32)
-25.27 (0.23)
sMTGP RS
96.63 (3.74)
146.22 (4.30)
-25.03 (0.13)
MLGP-RS
99.41 (0.76)
148.28 (1.64)
-25.27 (0.23)
MVGP-RS
96.32 (3.72)
147.20 (2.80)
-25.27 (0.33)
In the BO setting, additionally include single-objective GP-UCB baseline, which directly sumes that f follows GP and uses the UCB acquisition function to select query points sequen- tially. The optimization results under the three settings in our numerical experiments are summarized below.
Table 11: The results of Single GP-UCB in three synthetic settings.
Criterion
Single GP-UCB
TOBO
MSE,
0.0004
0.0000
Setting (1)
Ins Regret
0.0459
0.0001
Time
457.34
6968.25
MSE,
0.4761
0.0003
Setting (2)
Ins Regret
0.2093
0.0002
Time
187.30
1588.53
MSE.
0.0084
0.0001
Setting (3)
Ins Regret
0.1964
0.0302
Time
415.61
6035.61
We include an independent GP-UCB baseline (Ind GP-UCB), which assumes that each tensor el- acquisition function to
ement follows its own GP model and applies the proposed sequentially select both query points and super-arms. The optimization results under the three syn- thetic settings are summarized below.
26

---

Published as a conference paper at ICLR 2026
Table 12: The results of the Ind GP-UCB in three synthetic settings.
Criterion
Ind GP-UCB
TOBO
MSE,
0.2121
0.0023
MAE,
0.6739
0.0172
Setting (1)
Acc
0.67
1
Ins Regret
36.301
0.9807
Time (8)
328.56
8017.08
MSE
0.0731
0.0000
MAFx
0.0486
0.0000
Setting (2)
Acc
1
Ins Regret
0.1561
0.0000
Time (s)
110.77
2348.16
MAE
0.0139
0.0021
Acc
0.86
1
Setting (3)
Regret
54.2714
1.4406
Time (s)
3296.72
11699.09
For non-separable multi-output GPs (Fricker et al. 2013). we additionally provide the results under the three synthetic settings below.
Table 13: The results of the non-separable MOGP in three synthetic settings.
Task
Criterion
Setting
Setting (2)
Setting (3)
NLL
507.19
-39.33
-3094.10
MAE
0.1923
0.1493
0.1591
|Cov]]
6.7398
0,5618
20.5800
Time (8)
453.11
151.31
1397.68
MSE,
0.0014
0.0005
0.0244
MAE,
0,0536
0.0503
0.6277
BO
Ins Regret
1.0284
0.0030
6,6283
Time (8)
9109.08
2101.45
7881.31
MSE,
0.0073
0.0683
0.0091
MAE, 0.0643
0.0390
0,1257
CBBO
Acc
0.86
Ins Regret 1.7389
0.1551
51.4646
Time (s)
10462.18
3107.82
14175.50
As shown, TOBO achieves consistently better optimization accuracy than Single GP-UCB, Ind GP-UCB, and non-separable MOGP-UCB across all settings. Although the runtime of TOBO is higher, we further introduce an efficient variant named SVD-TOBO, which employs a low-rank cigen-decomposition to approximate the TOGP covariance Remark 3 and significantly reduces runtime.
To evaluate the sensitivity our method with respect to different choicess Ly and Hy, we consider three types of operators defined as follows (Chugh) 2020).
Sum operator (used in the main manuscript): Lff(x) = Ef (x), (x,
H,f(x, Σ f (æ, A).
Weiggttd sum operator: Lyf(x)
E Wi fin..nmm Hff(x,
Ew fe....d (I, A), where Wi ~ U).
.
Exponential weighted operator: T.ff(x)
Σ Eeme -1
H,f(x, =) ePw, 11 (x.)) where p = and Wij...... U U(,
The results under Setting (1) in numerical experiments summarized below. It can be that our method robust across different choices of Ly and H,, and both TOBO and TOCBBO consistently achieve strong optimization performance.
27

---

Published as a conference paper at ICLR 2026
Table 14: The results of different operators in three synthetic settings.
Criterion
Sum operator
Weighted sum operator
Exponential weighted operator
MSE, 0.0000
0.0000
0.0000
BO
MAE, 0.0008
0.0083
0.0072
Ins Regret
0,0001
0.0044
0.0533
MSE, 0.0023
0.0037
0.0025
MAE, 0.0172
0.0329
0.0548
CBBO
Ace 0.1964
1
Ins Regret
0.0107
0.4972
N REAL-WORLD APPLICATION: SEMICONDUCTOR MANUFACTURING PROCESS
A motivating example arises from semiconductor manufacturing, where each wafer consists of nu- merous dies (chips) arranged in a two-dimensional grid. During the Chip Probing (CP) phase, a critical stage for functional quality control, each die is evaluated based multiple quality variables such as voltage, current, leakage, and power consumption. These variables are spatially correlated across neighboring dies due to physical effects such process variation and mechanical stress, forming a naturally structured tensor output. To ensure high yield and reliability, manufacturers aim adjust process control parameters (inputs) that all quality variables across the wafer remain within target specifications. This leads to an optimization problem where the output is a three-mode tensor: the first two modes index the die positions on the wafer, and the third mode captures multi- ple quality variables. Such scenario cannot be effectively modeled by scalar- or vector-output BO approaches, as they would lose essential structural information.
In practice each wafer may contain hundreds of dies. However, resource and time constraints often make it infeasible to measure all quality variables across all die positions on every wafer. Manufacturers instead selectively measure subset of output entries, such as centrally located dies, those more prone to failure, or historically most informative regions. Similarly, only a subset of quality variables may be measured if certain tests time-consuming costly. This results in partially observed tensor, where only part of the full output is available at each iteration.
We incorporate this example into the revised paper and conduct corresponding case study. The input is æ = (æ1, 22, 23) representing process parameters, and the output f(x) E TRx65x3 denotes die-wise quality variables the wafer. A black-box semiconductor simulator is employed as the true system. We generate observations as the initial design and then sequentially select 20 queried points based on different BO methods. Tables 15 and 16 compare the performance of TOBO and TOCBBO against several baselines. Our methods sigmificantly outperform the alternatives in terms of input accuracy (MAEg), regret, and final objective value, demonstrating their superiority.
TOBO
sMTGP-UCH
MLGP-UCB
MVGP-UCB
MAE,
0.0651 (0.0016)
0.1997 (0.0020)
0.1579 (0.0059)
0,2669 (0.0046)
Regret
0.1702 (0.0007)
0.2400 (0.0014)
0.1956 (0.0088)
0,3818 (0.0029)
Objective
0.8298 (0.0009)
0.7600 (0.0031)
0,8044 (0.0079)
0.6182 (0.0033)
Performance comparison of TOBO with baseline methods in terms of MAE,, regret, and objective (mean and standard deviation).
Table 15:
TOCBBO
sMTGP- UCB
MLGP-UCB
MVGP-UCB
MAE,
0.1453 (0.0466)
0.4313 (0.0592)
0.5702 (0.0747)
0.6758 (0.0844)
Regret
0.1431 (0.0388)
0.2483 (0.0522)
0.4338 (0.0918)
0.3461 (0.0709)
Objective
0.6085 (0.0436)
0.5290 (0.0500)
0.4282 (0.1147)
0.4794 (0.0912)
Performance comparison of TOCBBO with baseline methods in terms of MAE,, regret, Table 16:
and objective (mean and standard deviation).
28

---

Published as a conference paper at ICLR 2026
These results show that the proposed methods can effectively optimize under partially observed tensor data, highlighting their practical relevance for semiconductor manufacturing.
29


Reviews:
[{"review_id": "MOZg2rs7hA", "reviewer": "Reviewer_6Rc1", "summary": "This paper proposes a novel Bayesian optimization framework for systems with tensor-valued outputs, a setting previously unaddressed in the BO literature. The core idea involves a Tensor-output Gaussian Process (TOGP) with two new kernel classes that capture complex, input-dependent correlations within the tensor structure, avoiding the limitations of simple vectorization. The authors also introduce a more challenging Combinatorial Bandit BO (CBBO) setting where only a subset of tensor elements contributes to the objective, and propose TOCBBO to handle it via a CMAB-UCB2 acquisition function that decouples continuous input selection from combinatorial super-arm selection. Theoretical guarantees are provided for the sub-linear regret of both methods. Extensive synthetic and real-world experiments demonstrate superior performance over baselines that use standard multi-output GPs.", "questions": "See the weaknesses part above.", "limitations": "", "rating": 4, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "1. This paper introduces the first Bayesian optimization framework specifically designed for tensor-output systems, developing novel tensor-output Gaussian processes that capture complex structural dependencies through specialized kernel designs.\n2. This paper provides strong theoretical guarantees, establishing sublinear regret bounds for both the standard tensor-output setting and the more challenging combinatorial bandit scenario with partial observations.\n3. The proposed methods show robust performance and consistently outperform existing approaches across diverse real-world applications.", "weaknesses": "1. This paper does not provide computational complexity comparisons with existing methods. Additionally, the high O(n^3T^3) complexity of TOGP itself represents a significant limitation, as strategies for scaling to large tensor outputs are not discussed.\n2. The theoretical analysis relies on the assumption that the true function is a sample from the proposed TOGP. There is little discussion of whether this holds in practice or how to assess its validity in real-world applications.\n3. A more detailed explanation of Equation (3) in line 161 would be helpful. For instance, it would be useful to clarify whether the kernel function K(x,x′) is designed to capture both similarities in the input space X and correlations among the internal elements of the output tensor.\n4. One concern is that the modeling and algorithm design are claimed to target tensor-type data, yet the actual formulations (e.g., Equations (4), (5), and (8)) are carried out at the matrix level. Would this process compromise or distort the intrinsic tensor structure?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "qTSrXB7m5p", "reviewer": "Reviewer_9BYU", "summary": "- The paper introduces Tensor-Output Gaussian Processes (TOGPs) that model structured tensor-valued outputs with novel separable and non-separable kernels capturing input-dependent correlations.\n- Building on TOGP, the authors propose Tensor-Output Bayesian Optimization (TOBO) using a UCB acquisition function for efficient optimization of expensive black-box tensor-output functions.\n- They further extend the framework to a Combinatorial Bandit Bayesian Optimization (TOCBBO) setting where only a subset of tensor outputs can be observed.\n- Both TOBO and TOCBBO are theoretically analyzed, with sublinear regret bounds.\n- Extensive synthetic and real-world experiments demonstrate that the proposed methods outperform several baselines in predictive accuracy and optimization effectiveness.", "questions": "- The paper's primary contribution appears to lie in the Gaussian Process modeling rather than in the Bayesian Optimization framework itself. Could the authors clarify whether the main contribution is intended to advance GP modeling or BO methodology?\n- Given that ICLR allows for a revision, the authors might be able to revise or extend the work to more clearly strengthen its contribution to the BO literature?\n- While theoretical computational complexities are analyzed, no empirical runtime or efficiency results are provided. Could the authors report or estimate the practical computational overhead of TOGP and TOCBBO compared with existing multi-output GP or BO methods?\n- How does the improved predictive accuracy of TOGP quantitatively translate into better optimization results such as faster regret decay or fewer evaluations? Have the authors conducted any ablation or correlation analysis to demonstrate this link?\n- Could the authors consider extending their framework to more meaningful tensor structures, such as spatiotemporal tensors that capture both spatial and temporal dependencies?", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "- The manuscript is well written with polished language that enhances readability and effectively communicates complex technical ideas.\n- It introduces a novel Tensor-Output Gaussian Process (TOGP) framework with separable and non-separable kernels.\n- The framework extends Bayesian Optimization to handle tensor-valued functions.", "weaknesses": "- The paper's main novelty lies in the Gaussian Process modeling, while the Bayesian Optimization layer mainly adapts standard UCB and CMAB-UCB2 frameworks without introducing a fundamentally new acquisition strategy.\n- The experimental evaluation emphasizes predictive accuracy metrics (MSE, MAE, NLL) more than optimization-oriented ones, creating a mismatch with the paper’s stated BO focus.\n- The paper lacks runtime or computational efficiency results, reporting only asymptotic complexity analyses without empirical timing comparisons.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "neaCDkWUWj", "reviewer": "Reviewer_AnPz", "summary": "This paper presents two variants of GP-UCB for BO of tensor-valued objectives. Both model the objective using a tensor-output GP (TOGP) with the goal of capturing additional structure present in the tensor representation, with a kernel that is somewhat novel to capture correlation between modes. The first algorithm uses scalarization - basically a weighted sum of the ouputs - while the second does the same on a subset of outputs (partial observations) at each round. Convergence rates are analysed to show sub-linear regret, and some experimental results are presented.", "questions": "Relevant questions were answered last time I reviewed this paper.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "The major strength/novelty of this paper appears to be the TOGP model, in particular the kernel used. And the kernel used for this GP is novel, at least to my understanding.\n\nThat said, at least the first algorithm presented does not seem particularly novel to me: it is basically a scalarizing multi-objective BO algorithm. Each element of the tensor represents an objective, and a linearization operator takes the place of scalarization. The second algorithm is arguably more novel, as it assumes partial observations of the tensor components at each iteration; however, past this, it again boils down to a scalarized MOBO formulation.\n\nExperiments seem to have improved since last time I reviewed this paper, and experimental result seem sound, so I am inclined toward accepting.", "weaknesses": "See above", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "wmaBgZ0qXr", "reviewer": "Reviewer_f43k", "summary": "The paper investigates optimization problems involving tensor-output objective functions, where each function evaluation returns a multi-dimensional tensor instead of a vector, as in standard multi-output settings. The authors introduce a Tensor-Output Gaussian Process (TOGP) that captures both dependencies across inputs and structured correlations within and between tensor modes of outputs. The proposed TOGP generalizes and introduces two types of kernels: a separable kernel, which assumes constant correlations across tensor modes, and a non-separable kernel, which allows correlations within and across tensor modes to vary with the inputs. \n\nThe authors further apply the TOGP to Bayesian Optimization (BO) in two scenarios. In the fully observed setting (TOBO), all tensor entries are available at each iteration, and optimization proceeds via an Upper Confidence Bound acquisition function derived from the TOGP posterior. In the partially observed setting (TOCBBO), only a subset of tensor entries is observed per iteration, transforming the problem into a Combinatorial Bandit Bayesian Optimization (CBBO) framework. Here, the algorithm must select both the next input and the subset of tensor entries (the super-arm) to observe. The paper provides theoretical guarantees, including proofs of kernel validity and sublinear regret bounds for both TOBO and TOCBBO. Extensive evaluations include three synthetic and four real-world benchmark problems to compare against many multi-output GP baselines.", "questions": "1. What is the runtime trade-off of the proposed TOGP? Since TOGP models complex interactions across multiple tensor modes, it likely incurs higher computational costs than both single-task GPs and multi-output GPs based on the LMC, due to the larger number of hyperparameters involved. Additionally, given that the problem settings already require a scalarization operator (L_f or H_f) to map tensor outputs to scalar objectives, which allows single-task GP to work just fine, it would be useful to analyze the accuracy–efficiency trade-off - that is, how much performance gain the tensor-output model provides compared to the reduced computational cost relative to a simpler single-task GP baseline.\n2. Does the proposed TOBO algorithm account for the diversity of solutions? Given that the target problems have multiple outputs, some of which may be conflicting, can TOBO handle such cases, identify a set of solutions that satisfy different objectives or trade-offs? This is a common consideration in MOBO literature, where the goal is not only to find the best single solution but also to capture the Pareto diversity among competing objectives.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "-\tThe idea is sound, which generalizes the well-known LMC model in multi-output (vector-output) GP into the tensor-input GPs.\n-\tThe paper is clearly written, with theoretical proof for the kernel validity and regret bound under certain assumptions.\n-\tThe empirical performance shows that TOGP has potential in modelling tensor-output functions, when compared to vector-output GPs.", "weaknesses": "Some experiments are missing: \n- Lack of single-objective BO baselines: Since the target problem settings rely on a scalarization operator (L_f or H_f in the CBBO case) to convert tensor outputs into scalar objectives, standard single-objective BO methods (e.g., single-task GP combined with a UCB acquisition function) should, in principle, operate normally. The paper would be stronger if it compared TOGP-based optimization with such single-task GP baselines to demonstrate the benefit of modelling complex tensor dependencies rather than relying on a simpler approach.\n- Lack of single-objective BO baselines in the Bandit BO setting: Similarly, in the Bandit BO setting, the presence of the scalarization operator (L_f or H_f) allows single-objective BO algorithms to be applied directly. Prior works (e.g., Nguyen et al., 2020; Ru et al., 2020; Huang et al., 2022) already use this setup effectively, and including these as baselines would provide a fairer and more comprehensive comparison for TOCBBO.\n- Lack of multi-objective BO baselines: Because TOGP inherently models multiple correlated outputs, it would be valuable to assess its performance under a multi-objective optimization framework. Comparing against the SOTA MOBO methods using metrics such as the hypervolume indicator (e.g., Maddox et al., 2021; Daulton et al., 2022) would further complement the paper.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "4esbmGYxqG", "reviewer": "Reviewer_uCf6", "summary": "The paper introduces Tensor-Output Bayesian Optimization, a framework for optimizing black-box functions with tensor-valued outputs using a tensor-output Gaussian process. It proposes both separable and non-separable tensor-output kernels that capture input-dependent cross-mode correlations, and develops a UCB acquisition tailored to scalarized tensor objectives. The work further extends to a combinatorial setting where only a subset of tensor entries contributes to the objective, using a partially observed TOGP and a CMAB-UCB2 criterion for joint input and subset selection. Theoretical sublinear regret bounds are provided, and experiments on synthetic and real tasks show consistent improvements over vectorized multi-output GP baselines.", "questions": "- How are $L_f$ and $H_f$ instantiated across tasks, and how sensitive are results to these choices or to non-Lipschitz scalarizations encountered in practice?\n- How are tensor ranks chosen for the core tensors (line 208), and what is the performance/overfitting trade-off as ranks vary under separable vs non-separable kernels?", "limitations": "", "rating": 6, "confidence": 1, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "Modeling tensor outputs directly with non-separable, input-dependent kernels is a clear step beyond vectorized MOBO, enabling richer cross-mode structure to be exploited. The CBBO formulation and the PTOGP-based acquisition rule are principled. The regret analyses leverage information gain and cover both full and partial observation regimes, giving the method theoretical footing. Empirical results are broad, with synthetic and multiple real datasets, and show strong gains over standard multi-task and multivariate GP baselines. The appendices give thorough algorithmic and complexity details.", "weaknesses": "The computational complexity $O(n^3 T^3 log(n))$ of TOBO is cubic in both n **and** tensor size, which may severely limit practical scale. \n\nThe scalarization operators $L_f$ and $H_f$ are central in practice, but their choices, sensitivity, and domain-specific instantiations are only lightly discussed. \n\nBaselines omit state-of-the-art tensor GP surrogates with Kronecker structure and non-separable multi-output kernels beyond the selected MTGP/MVGP/MLGP such as proposed in \"Scalable High-Order Gaussian Process Regression\" (line 91).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)