Optimizing computationally expensive black-box functions is challenging because each function evaluation incurs a high computational cost. Bayesian optimization (BO) addresses this problem by using Gaussian process surrogate models to guide the selection of new evaluation points. Within the Bayesian optimization with Pareto-based selection (BOPS) framework, the predictive mean and standard deviation define a bi-objective optimization problem that balances exploitation and exploration. A batch of solutions is then selected from its approximate Pareto set. This paper presents a systematic empirical study of Pareto front approximation and Pareto-based solution selection within the BOPS framework. We investigated three space-filling designs (Sobol sequences, Latin hypercube sampling, and uniform random sampling) in combination with five Pareto-based solution selection strategies. Experiments were conducted on the 20-dimensional CEC2022 benchmark suite (12 functions) and a real-world robotic manipulation problem using a batch size of $q=4$, an evaluation budget of 800 function evaluations, and 30 independent runs. BOPS variants were compared with standard acquisition-based batch BO methods, including $q$-expected improvement, $q$-upper confidence bound, and $q$-knowledge gradient. The results showed that optimization performance depends primarily on the solution selection mechanism rather than the Pareto approximation strategy once a sufficiently accurate Pareto set is obtained. Moreover, several BOPS variants using simple space-filling-based Pareto approximation achieve optimization performance comparable to or better than standard acquisition-based BO. These findings provide empirical insights into Pareto-based batch Bayesian optimization, demonstrate the effectiveness of simple space-filling-based Pareto approximation, and provide a unified framework for developing and evaluating future Pareto-based solution selection strategies.

We introduce a generating-function framework for analyzing the no-feedback card-guessing game after $k$ Gilbert–Shannon–Reeds riffle shuffles. We show that the distribution of the card appearing in position $i$ can be expressed as a structured mixture of $2^k$ tractable components, each corresponding to a sum of independent Bernoulli trials. From this decomposition, we derive an explicit closed-form expression for the probability generating function, represented as a product of binomial-type polynomials with a clear and systematic structure, valid for any number of cards $n$ and any number of shuffles $k$. This formulation replaces recursive convolutions with a single analytic expression, enabling efficient computation and revealing the combinatorial–probabilistic structure underlying riffle shuffles. Beyond exact evaluation, the framework connects optimal no-feedback strategies with the generating functions and suggests asymptotic behavior in both the fixed- $k$, large- $n$ and fixed- $n$, large- $k$ regimes.