Many animals possess a remarkable omnidirectional sound localization ability enabled by subconsciously processing subtle features in the sounds received at the two ears from a common source location. For humans, these features arise due to the incoming acoustic wave scattering off the listener’s anatomical features (head, torso, pinnae) before reaching the eardrum. The spectral ratio between the sounds recorded at the eardrum and that would have been obtained at the center of the head in absence of the listener is known as the head-related transfer function (HRTF) ; HRTFs are thus specific to the individual’s anthropometry, wave direction, and contain other important cues such as the interaural time delay (ITD) and the interaural level difference (ILD) . Moreover, knowledge of individualized HRTFs allow for perceptually accurate D spatial audio synthesis [3, 4, 5].
We investigate the pre-image problem, namely how pairs of left and right ear HRTFs and functions of HRTFs (features based on them) map back to their measurement directions. This is related to the problem of sound-source localization (SSL) where under simple (anechoic) conditions, the direction of an acoustic event can be inferred from multi-receiver recordings of the sound spectrum by expressing the spectral cues solely in terms of the receiver’s transfer functions (independent of their actual content). This is of interest in robot perception (e.g. for event detection and localization [6, 7]), where the receiver’s transfer functions can be measured beforehand. For humans, this problem is restricted to two receivers (human ears) where functions of left and right pairs of HRTFs are mapped to their measurement directions in place of SSL directions. Thus, it possible to model this relation as either a classification or a regression problem between the two domains. Many works in literature have attempted similar tasks.
I-a Prior Works
Cue-mapping  uses ITD, ILD, and interaural envelope difference features paired with azimuth directions in a weighted kernel nearest-neighbor (NN) setting. A linear mapping between ITD, ILD, and HRTF notch frequency features to spherical coordinates can be learned 
. A self-organizing map between input ITD, spectral notches features and output horizontal and median plane coordinates can be trained10]. A probabilistic affine regression model between interaural transfer functions and the direction is possible .
Most closely related to our work are the source-cancellation and match-filtering algorithms [12, 13, 14, 15], where the binaural recordings ( left, right ears) are represented as convolutions of a common sound-source signal and the appropriate filters; for recording done in an anechoic space, these filters are the same-direction HRTFs ( left, right ears). The per-frequency domain representation is given by
where is element-wise product. The source-signal is removed by computing the ratio between left and right channel recordings (
). These binaural features, which are reduced to ratios of HRTFs, can be compared to those pre-computed from the subject’s collection of measured HRTFs; the measurement direction belonging to the maximally cross-correlated pair is reported as the sound-source direction. Such an approach can be interpreted as a nearest neighbor (NN) classifier where the binaural features and measurement directions are single class instances and labels respectively.
I-B Present Work
We propose a generalization of the match-filtering algorithm that addresses several deficiencies: While an NN classifier is accurate for a large number of training samples, it does not report out-of-sample spatial directions unless specified in a regression context. Linear regression methods via ordinary least squares (OLS) regressors111, for parameters often perform poorly due to inaccurate assumptions on the model complexity (number of parameters) and the linearity between predictors and outputs. Common issues include over-fitting the model to noise that arise from parametric OLS methods and under-fitting the training data from assumptions of linearity. Instead, we adopt a non-linear and non-parametric222Number of parameters is proportional to the number of data samples conditioned upon for inference. Gaussian process (GP) regression (GPR)  framework to address these issues.
GPR is a kernel method333Predictor variables are implicitly mapped to a reproducing kernel Hilbert space whose inner products are taken to be evaluations of a valid Mercer kernel or covariance function.
that places weak assumptions on the joint probability distribution444Normal distribution defined by prior mean and covariance functions of predictor variables (binaural features). of latent function realizations
that would model the output observations (spatial directions) in a Bayesian setting. Observations are drawn (realized) from a high-dimensional normal distribution that represents the joint probability density function of a collection of random variables indexed by their predictor variables. GPs have several attractive properties that are well-suited for SSL.
Based on the observation that HRTFs corresponding to different spatial directions covary smoothly with the considered binaural features (see sections III), we show they can be modeled via simple stationary GP covariance functions (see section IV
). The GP Bayesian formulation allows for the choice of the covariance function, which governs the smoothness between realizations at nearby predictors, to be automatically selected by evaluating a data marginal-likelihood criterion (goodness-of-fit); covariance functions belong to a function class and are specified by their “hyperparameters” (parameters that describe distributions). This allows the covariance function hyperparameters to be learned without the need for cross-validation and provides insights as to the intrinsic dimensionality of the high-dimensional feature space that the binaural features are mapped to. Most importantly, uncertainties in GP prediction are well-defined in terms of both prior and posterior distributions; the predicted variances at different inputs are tractable. Thus, GPR generalizes NN classifiers as it makes non-linear inferences to observations outside the training set. By the representer theorem, kernel methods such as support vector regression (SVR) and GPR make predictions expressible as linear combinations of non-linear covariance evaluations between the training features/observations and the test features.
In general, GPs perform better (make accurate inferences) with more observations (data) than other non-linear regression methods that do not encode and select for prior data-assumptions. The trade-off is its high computational costs ( operations for
number of observations) for both model-selection and inference; scaling GPs for for large datasets is an active field of research. Fortunately, the availability of high quality datasets, computational resources, and faster algorithmic formulations have allowed us to overcome these problems. In previous works, we have used several properties of HRTF datasets to to perform fast GP based HRTF interpolation and data-fusion . The current work is a major extension of our recent work on binaural SSL . For future references, we refer to GPs that predict SSL directions as GP-SSL models (see section IV for a complete derivation).
Ii Formulation of Problems
This work investigates two problems related to GP-SSL models (see Fig. 1 ). For notation, we refer to a binaural feature as a -dimensional vector ( is number of frequency bins), the measurement direction as the unit vector ( for the standard Cartesian basis), and collections of the aforementioned quantities ( number of samples) as concatenated into matrices and . The binaural features are independent of the sound-source content and thus strictly functions of the subject’s HRTFs (see section III). GP-SSL models are thereby specified and trained over known HRTFs and measurement directions belonging to CIPIC  database subjects.
Ii-a Feature subset-selection
Subset-selection for non-parametric methods such as NN and GPR is an important technique for reducing the model-order complexity and run-time costs for inference. SSL models that are trained with randomized subsets of samples trade measurement and prediction costs for localization accuracy. Increasing the density of measurement samples over the spherical grid results in a linear increase to both NN classification computational cost and accuracy, a quadratic and cubic increase to respective GP inference and training computational costs, and a non-linear increase to GP localization accuracy. We show how GP-SSL models using small and non-uniform subset-selected samples (which are most informative) make more accurate predictions over the full spherical grid than models evidenced on a randomized subset.
A simple greedy forward-selection (GFS) algorithm  that sequentially incorporates training samples into a subset without considerations in future iterations is implemented. It ranks all training samples outside the subset via a user-defined objective function (risk function) and adds the minimizer into the subset. We propose a class of risk functions that generalizes the GP prediction errors and show that the subset-selected GP-SSL models localize directions more accurately than models evidenced on randomized inputs (see section V); only a small fraction of training samples are required for reasonable accuracy ().
Ii-B Active-learning for individualizing HRTFs
, a number of works have sought indirect means for learning the subject’s HRTFs: regression models between the individual’s physically measured anthropometry and his/her HRTFs can be learned via neural-network and multiple non-linear regression models  but do not generalize well to test subjects. HRTFs can also be learned through listening tests [26, 27] by having an individual listen to a query HRTF x convolved with white Gaussian noise (WGN) (heard over a pair of headphones), localize the test signal (report a direction ), and then hand-tune the spectra of x or choose a new x out of a large candidate pool over a graphical user interface (GUI) as to move subsequent localizations towards a target direction
. The hand-tuning/selection step can be replaced by developing a recommendation system that selects for the query HRTF between rounds (steps) of localization. The listener can rank candidate HRTFs chosen from a genetic algorithm555Evaluates a fitness function w.r.t. localization accuracy of known u 
. HRTFs can also be tuned along a low-dimensional autoencoder space where u is unknown to the listener.
We propose to formulate the recommendation problem in an active-learning  context described as follows: given a finite set of candidate HRTFs sampled from a prior distribution (database or generative model), determine the HRTF from the that the listener would localize nearest to u within rounds of localizations. During round , the recommender selects a query x that the listener labels as without knowledge of u. The choice of x is referred to as the query-selection problem of minimizing the SSL error (SSLE) (modified cosine distance) given by
Unfortunately, the minimizer in Eq. 2 is unlikely to be found within rounds as can be large and must also be small as the cost of evaluating SSLE by the listener is high. It is more reasonable to model the SSLE function using an online regression model (adapting HRTFs predictors of SSLEs after each round) and select for x based on two competing strategies: query-selection exploits the online model by choosing x that the model predicts will have low SSLE and explores x that has high model uncertainty in its prediction; both concepts are trade-offs that require probabilistic treatments of model predictions. Fortunately, GPs are well-suited to this task as all predictions are expressed as probabilistic realizations sampled from normal distributions. Thus, we propose to solve the modeling problem via GP-SSLEs666GPs that predict the SSLE from HRTFs, and the query-selection problem using a method of GPs for the global optimization of smooth functions [31, 32] (see section VI). The relation between these methods and the GP-SSL models is also shown.
Iii Binaural Sound-Source Invariant Features
We consider several sound-source invariant features that can be extracted from short-time Fourier transforms of the left and right ear input channel streams in Eq.1 (see Table I and Fig. 2); it is useful to express the discrete Fourier transformed signals by their magnitude and phase representations where . The features are expressed as ratios between left and right ear channel recordings that remove the effects of the acoustic content in ; the remainder is strictly a per-frequency function of same-direction left and right ear HRTFs derived as follows:
|Avg. magnitude ratio|
|Magnitude pairs for flat|
Log-magnitude ratio (LMR) : While the source-cancellation method removes the dependence on signal , the resulting features are complex, noisy, and difficult to interpret. This can be avoided by considering the magnitude representation which gives the relative per-frequency energy between the channel signals. Adding a constant to the ratio prior to the log-transform penalizes the magnitude of the perturbation; adding a constant constrains the log-transform to be non-negative.
Phase difference (PD): Similarly, the per-frequency phase of the complex channel signal ratio can be expressed by the phase-difference between left and right HRTFs. For identical that differ by onset time-delays , the phase-difference is simply the constant delay across all frequencies; this ITD can be related to azimuth angles via Woodworth’s model . For arbitrary , the per-frequency phase-differences differ and are to be treated as independent variables in regression models.
Average magnitude ratio (AMR): The magnitude source-signal can also be removed by taking the ratio of left or right magnitude signals and the binaural average . Without the constant factor, the feature can be interpreted as the per-frequency contribution of the left or right magnitude HRTFs to the additive binaural magnitude response. Unlike log-magnitude ratio features that approaches a singularity as , these features are bounded in the interval and finite everywhere unless the binaural average is zero.
Magnitude pairs (MP): The magnitude pairs are the concatenation of the original left and right magnitude HRTFs that could be derived from convolution with a WGN with zero mean and unit variance. The power spectrum of is constant across all frequencies and so would be constant factors of magnitude HRTFs. Such conditions arise during listening tests where the source-signal can be specified; the test features can then be derived from per-frequency division given by and .
Binaural features extracted from CIPIC subjectHRTFs are shown for horizontal and median plane directions.
Iv Gaussian Process Regression for SSL
In a general regression problem, one predicts a scalar target variable from an input vector x of independent variables based on a collection of available observations. A common Bayesian approach for inference assumes that the observation is generated (realized) from a latent function given by
which is corrupted by additive Gaussian white noise with zero mean and constant variance. This latent function is given the form of a kernel regression where the function maps the inputs x into a high-dimensional space before computing the inner product with a vector of parameters realized from a collection of random variables with a prior multivariate normal distribution . Unlike linear regression, the parameters
are not explicitly found in order to perform inference but are marginalized in order to compute the first two moments (mean and covariance) of functiongiven by
The latent function realizations are thus drawn from a multivariate normal distribution with mean and variance . For , the inner product can be replaced with the covariance function which GPs generalize as follows:
A GP is a collection of random variables where any finite subset indexed at inputs has the joint multivariate normal distribution given by
and thus fully defined by the prior mean function and the prior covariance function . The prior mean function and vector are set to zero without loss of generality following Eq. 4. The covariance (Gram) matrix is characterized by the pairwise covariance function evaluations
; the covariance function is a positive semi-definite kernel (Mercer’s condition) that establishes the existence of the eigenfunction. This allows kernel methods such as SVR and GPR to omit computing the exact mapping as the inner products in the high-dimensional space, representing the similarity measure between input features , are well-defined.
GP inference at test inputs evidenced on training inputs and the observations in derives from the multivariate normal distribution of random variables conditioned on , . This is given by
where adjusts for the observation noise and are pair-wise covariance evaluations between training and test inputs. We refer to the distribution in Eq. 6 as the posterior GP defined by the posterior mean and posterior covariance functions and respectively. The former represents the vector of expected outputs (prediction) at
and the latter is gives the confidence intervals (diagonal of the matrix) of the predictions.
For the GP-SSL model, and are the respective binaural features in Table. I and their measurement directions (unit vectors where are values along the coordinate); test inputs refer to the binaural features extracted from test signals. While it is possible to model all output coordinates as a collection of independent GPs , a computationally cheaper alternative is to specify a common prior mean and covariance function shared by all GPs. Specifying a shared covariance model between GPs is reasonable as the original HRTFs are originally measured over the same physical topology of a human subject from a near-uniform spherical grid of directions. Thus for inference, we use three independent GPs, with shared priors, to model left-right, front-back, and top-down coordinate directions by either sampling from their posterior distribution or reporting their posterior means.
Iv-a Choice of Covariance Functions
The “smoothness/correlatedness” of realizations of for similar depends on the number of times that the covariance function is differentiable w.r.t. the input arguments. Consider the Matérn class of covariance functions where each function has varying orders of differentiation. For -dimensional inputs, we can specify the GP covariance function as the product of -independent Matérn covariance functions of identical class. Three common classes and the product covariance function are given as
for distance and hyperparameters . Covariance functions are times differentiable and stationary due to their dependence on . Each function contains a length-scale or bandwidth hyperparameter that represents a distance in the domain where outputs remain correlated; larger length-scales result in smoother .
A general hyperparameter is optimized by maximizing the data log-marginal likelihood (LMH) of the observations given the GP prior distributions; the derivation follows from integrating over the realizations by the product of data likelihoods (sampling from and sampling from the GP prior distribution). The LMH term and its partial derivative are both analytic and given by
where is the matrix of partial derivatives. A larger LMH represents a better goodness-of-fit of the data to the GP prior mean and covariances assumptions. Moreover, different covariance functions with optimized hyperparameters can be compared in this respect without resorting to domain-specific metrics.
Iv-B Model-Order and Cost Analysis
The GP model-order is proportional to the size of the GP prior distribution defined by the -dimensional multivariate normal distribution in Eq. 5 ( is the number of training samples). The associated costs of both conditioning on the GP prior distribution for inference and performing hyperparameter training is dominated by the inversion of the Gram matrix ( operations to compute and space to store). For large , exact GP becomes intractable and most practitioners rely on randomized sampling techniques  to reduce the costs at the expense of accuracy. Two types of analyses for evaluating this trade-off are given: first, empirical cross-validation experiments can demonstrate how data sampling (randomized and subset-selection) increases localization error. Second, the theoretical dimensionality of the feature space in Eq. 3
, despite not having been explicitly computed, can be estimated from an eigenanalysis of the GP Gram matrix. The distribution of eigenvalues (number of dominant ones) gives a minimum bound as to the number of input features whose mapping will contain most of the variances in the feature space.
To evaluate the dimensionality of
, we refer to the method of kernel principal component analysis of Gram matrix
. Its derivation expresses the eigenvectors(principal directions) and eigenvalues (measure of variance captured by ) of the sample covariance matrix of features in the high-dimensional space in the form of
where are the component scores between the feature mapping and the eigenvector. Applying the “kernel” trick allows to be reformulated in terms of the Gram matrix as a tractable eigendecomposition problem given by
which finds the eigenvalues and scores . Evaluating the contributions of the leading to the total energy estimates the number of eigenvectors that are relevant to .
GP-SSL models (input binaural features LMR, PD, AMR, and MP from Table. I belonging to CIPIC subject ) are trained (batch gradient descent of all covariance function hyperparameters via Eq. 8) for iterations. For a domain-metric, we use the angular separation distance between two directions u, u’ (predicted and reference directions) given by
Goodness-of-fit: GP-SSL models are specified/trained on the full set of inputs . The data LMHs in Table. II are computed for several covariance functions and feature types. The infinitely differentiable squared exponential gives the best-fit (highest LMH) across all features (latent functions modeling the SSL directions are smooth w.r.t. changes in the feature space). This confirms the fact that a finite collection of HRTFs approximates a sound-pressure field that is continuous in space. The best-fitting binaural features are the MPs (WGN sound-source) and AMRs (arbitrary sound-source); the LMH gap between the two suggest that GP-SSL models will perform more accurately when the recorded magnitude spectra match that of the HRTFs. The LMH gap between AMR and LMR suggests that relative contribution may be a better indicator of SSL than relative intensities. The low LMH of PD models suggests that phase may not be useful for SSL over the entire spherical coordinate system.
Eigenanaylsis of : The eigenvalues of the are computed for GP-SSL models trained/specified on the full dataset (). Fig. 3 shows the contribution of the leading eigenvalues to the total energy; specified by the four earlier features (LMR, PD, AMR, and MP) require respectively , , , and leading eigenvectors to capture of the total variance. The results suggest that feature mappings for MPs and PDs can be approximated with only a few samples while LMR and AMR feature mappings are more complex.
Cross-validation: GP-SSL models are trained on a randomized third of the available feature-direction pairs ( out of ); inference follows Eq. 6 at all available inputs () where only the posterior mean directions are reported. Table III shows the mean angular separation (Eq. 11) between predicted and reference directions for GP-SSL, NN classifier, OLS methods trained on the same data. Non-parametric methods (NN and GPR) outperform parametric methods (OLS) across all feature types. The MP and AMR features give the lowest errors across all methods (for a visual, see the first column of Fig. 4). OLS log-ratios perform the worse and suggest that the features are oversensitive linear predictors of change in localization direction. PD features, while useful for predictions on the horizontal plane, are insufficient for localizations over the full sphere.
V Feature Subset-Selection
Greedy feature selection is an efficient method for finding a subset of inputsthat best approximates a functional according to a user-specified risk function (measure of distance between and ). Determining the optimal subset via an combinatorial exhaustive search is prohibitive w.r.t. the number of evaluations of
. A greedy heuristic (rankingaccording to a point-inclusion in the risk evaluation and adding the minimizer into the subset without consideration in future iterations) reduces the search to a quadratic number of evaluations (see Algorithm 1). For GP-SSL, GFS approximates the GP posterior distribution (Eq. 6) evaluated on the full dataset () conditioned on a growing subset of inputs. We propose an efficient method for updating both GP prior and posterior distributions between point-inclusions in section V-A.
Specifying the risk function is more difficult as its evaluation costs must be low. Most risk functions that use second-order moments (e.g. GP posterior covariance in Eq. 6) are expensive and require approximations to remain tractable . Evaluating the GP posterior covariance requires space; its inverse and determinants are expensive to compute in sub-cubic time. Instead, we propose a cheaper class of risk functions that generalizes only the first-order moments (i.e. GP posterior mean in Eq. 6) in section V-B.
V-a Incremental GP Models
A point-update to a GP model can be defined in terms of changes to the first/second moments of the GP prior and posterior distributions (Eqs. 5, 6) and both the Gram matrix and its inverse generated from inputs in . While a point-update to simply contains an appended row and column of covariance function evaluations , its direct inverse would be expensive to compute. Instead, we define a recurrence relation with its previous inverse as follows.
Given a sample input-output pair for data index , let indices be the union with the subset indices . At iteration , append a row and column vector along the standard basis to the Gram matrix . The differences between and the appended are two rank- updates given by
where vectors , , , and is the
column of the identity matrix. The update in Eq.12 allows to follow from the modified Woodbury formulation  given by
which requires only two rank- updates. For a fixed set of test inputs , the updated posterior mean vector remains a matrix-vector product and the posterior variances are sums of diagonals given by
where matrix . The updated log-determinant is given by . The total computational costs of updating the GP prior and posterior distributions at iteration are and operations respectively.
V-B Gp Risk Function Criterions
We show how several risk functions can be derived from the distance between any two GP posterior mean functions evaluated at a possibly infinite sized set of test inputs . Given two GPs , defined over the subsets of inputs for indices and , the distance between their two GP posterior mean functions ( and ) is analytic under certain GP prior assumptions. For prior mean and the product of identical Matérn class covariance functions in Eq. 7, the errors evaluated at are given by
where vectors , are computed over training data. Updating the risk function evaluations between successive iterations is efficient as updating , need only rank- updates via Eq. 13. The associated matrices , , in Eq. 15 are sub-matrices of and can be pre-computed in operations. Computing depends on the following cases.
Finite Case: If is finite, then matrices , , and are the summation of outer-products whose entries are products of Matérn class covariance functions in Eq. 7.
Infinite Case: If is the full (unbounded) input domain, then matrices , , and contain improper integral entries. For a valid distance measure, the posterior mean functions converge to identical zero-mean priors at the limits and the improper integrals of the form given by
are shown to be finite (see Appendix Eq. 23). Several combinations of the distance are summarized as follows.
Prediction Error : The prediction error is taken between the GP posterior means at test inputs and the known sample pairs .
Generalized Error : The generalized error is taken between two GP posterior mean functions and evaluated at any finite (may be out-of-sample from ). For GFS, the two GPs are specified by subset-selected and the full set of inputs .
Normalized Error : The normalized error or ”frequentist“ risk is taken between two normalized GP posterior mean functions ( and ) evaluated at given uniform probability distribution over . The norm term is shown to be finite by setting either of the functions in Eq. 15 to zero. The two GPs are specified on subset-selected and the full set of inputs .
GFS selects for increasing subset sizes until it contains the full dataset. At each iteration , the incremental GP-SSL model infers directions (posterior means) along test inputs . The mean angular separation error (Eq. 11) between the predicted and the reference measurement directions are computed and shown in Fig. 5; intercepts with horizontal lines indicate subset sizes at and errors. The crossover points at the error line (localization accuracy) are achieved for MP and AMR features at a small fraction of the total input set (approximately and feature-direction pairs); decreases in localization error after randomized samples becomes logarithmic with diminishing returns. Moreover, GFS selected models generalize better than that of randomized selection in all but the PD features; a visual (second column plots in Fig. 4) shows that the former more accurately localizes directions further from the median plane.
Vi Active-Learner System
The active-learning process for inferring HRTFs is as follows. The collection of number of target directions is specified as . For rounds , a query HRTF (MP) is chosen from the candidate set and appended to form input matrix . The listener localizes , registers the direction over a GUI (see Fig. 6), and appends the directions to form matrix . The SSLEs w.r.t. are computed in s.t. . Last, the updated feature-direction pairs are added into the GP-SSLE models via incremental GPs (section V-A). The system components are organized below.
Vi-a Conditional Mixture of Gaussians Models
While it is possible to specify an entire HRTF database as the candidate set, it is reasonable to assume that most samples would not be localized near a target direction u; overt features arising from the reflections off the anthropometry may be a physical impossibility along all measurement directions. Conversely, choosing only HRTFs with measurement directions equivalent to u restricts the sample size to the number of subjects in the database. To address both issues, we model both the HRTFs and their corresponding measurement directions using a conditional mixture of Gaussians model (MoG) trained from the CIPIC database (see section VI-A). This allows for to be drawn from a distribution of HRTFs conditioned at any direction u.
The MoG models the joint distribution between input variables as if the samples are drawn from a latent set of normal distributions. The input variables consist of measurement directionsu and leading principal components (PCs)777PCs are computed from same-subject, mean-centered, log-magnitude pairs (concatenated left and right ear HRTFs). w associated with HRTFs along u. The joint distribution is modeled by a weighted sum of normal distributions with mean and covariances given by
are trained via the well-known expectation-maximization algorithm. The PCsw conditioned on u is also a MoG given by
where the conditional mean and covariance for the mixture are and respectively. The candidate set is given by PCs randomly sampled from the conditional MoG888Leading PCs are sampled (via Gibbs sampling) from one of multivariate normal distribution (randomly selected by weight). in Eq. 18 and decoded into HRTFs to form the candidate set. The non-individualized (directional-averaged) HRTFs are approximated by the sum of the weighted conditional mixture means.
Vi-B GPs for Modeling SSLE
GP-SSLE models () are specified by a common set of input MP features and output SSLEs for each of the number target directions in . Accurate modeling of the SSLE depends on the choice of GP prior mean and covariance functions. A zero mean prior is reasonable as reported directions v in the absence of localization should average to the zero vector. Choosing the GP covariance function is more difficult as the hyperparameters cannot be optimized in the absence of observations; inaccurate priors would result in poor generalizations error.
Fortunately, GP-SSLE models can be related to GP-SSL models when is the infinite set of target directions uniformly sampled over a unit sphere. Substituting the SSLE labels into Eq. 8, the GP-SSLE LMH is now given by
where matrix . As , the sample covariance of approaches a constant variance due to symmetry. The LMH in Eq. 19 reduces to
which is equivalent to that of GP-SSL models for MP features and directions .
The equivalence allows for the choice of the GP-SSLE model’s covariance function to approximated by that of GP-SSL models trained over known feature-direction pairs (e.g. CIPIC subject data). While these subjects are not identical to the listener, the trained GP-SSL models all share similar covariance functions as their hyperparameters are well-distributed (see Fig. 7); high frequency bands above kHz tend to be negligible while lower frequency sub-bands between and kHz are relevant.
We present GP based query-selection as a modification of a known algorithm  which is derived as follows. Consider the observed minimum SSLE for any u at round given by
Realizations of SSLEs () by the GP-SSLE posterior distribution (Eq. 6) at a candidate input will be normally distributed whose mean and variances represent the expected SSLE and uncertainty respectively. Thus, improvements (lowering) upon the global minimum
is given by the loss-functionwhose expectation can be computed via marginalizing over the .
The expected loss-function is analytic for any single u and so the weighted expected loss function (specified over each with independent GP-SSLE models) is given by
where weights can be set to a constant, GP-SSLE posterior mean and covariance functions at evidenced with are denoted by and , and the cumulative normal distribution of is denoted by . The query HRTF is chosen as the lowest scoring candidate or minimizer of the criterion Eq. 22 which balances local improvement through the posterior mean term with exploring uncertain predictions through the posterior variance term . The property is useful for proving the rate of convergence  to the true solution in Eq. 2.
GP-SSL active-learning trials: One method for fast and repeatable empirical validation substitutes the human listener for GP-SSL models trained on CIPIC subject data. Localizations at can be reported as either the GP posterior mean directions, or by sampling from the GP posterior distribution. This allows for large subsets of to be efficiently evaluated with little time costs. For coherence, we limit the query-selection criterion in Eq. 22 to single target directions u belonging to the CIPIC HRTF measurement directions (queries made for past u are discarded). GP-SSLE’s covariance hyperparameters are set to that of the GP-SSL mean hyperparameters (averaged across subject models); hyperparameters can be retrained after each round but is not necessary for improving the localization error. The variance term is set to .
In tests, the active-learner submits an initial non-individualized query HRTF for u and then proceeds through rounds of query-selection from a candidate HRTF set of samples drawn from a conditional MoG (Eq. 18). The nearest localized directions are shown to closer to their target directions than the non-individualized guesses (see Fig. 8). Non-individualized HRTFs are localized closer to the horizontal plane and towards the back of the head. Nearest localized directions accord with empirical studies of difficulties in front-back and up-down confusion with human subjects . The experiment is repeated across all GP-SSL CIPIC subject models (see Fig. 9). The improvement can be expressed as the mean ratio between the angular separation errors of the initial and nearest localized directions. The mean improvement is across all CIPIC measurement directions, for median plane directions, and for horizontal plane directions.
Human active-learning trials: For a human listener, we develop a simple GUI in Matlab that consists of an azimuth-elevation plot that the subject clicks to report . To introduce contrast in hearing, two test signals are alternatively played over headphones until the listener reports a direction. The first is a short burst of WGN independently generated for left and right ear channels. The second is the WGN convolved with the left and right min-phase HRTFs derived from the binaural MP features. The trials proceed as the listener localizes queries for rounds in each of the target directions ( on the horizontal and median planes each).
For sample human listeners, the initial and nearest (minimum) localization errors for each of the target direction are shown in Table IV and are compared to synthetic trials conducted with the GP-SSL CIPIC subject models. In both cases, the largest errors occur along the median plane direction . The mean percentage improvements of the nearest localizations over that of the non-individualized HRTFs are and for human and GP-SSL listeners respectively. GP-SSL localization errors are generally lower and more consistent across all direction than the human listeners; GP-SSL models can report a posterior mean direction whereas human listener exhibit variances in his/her localizations, even for identical test signals. It may be of interest in future work to both measure and model human localization variances via the GP-SSL’s variance term and by sampling localizations from the GP posterior distribution.
We developed a robust method for the SSL using sound-source invariant features derived from left and right ear HRTF measurements. Our GP-SSL models generalized NN based approaches and were shown to more accurate in both cases of randomized and subset-selected features; good spatialization accuracy () over the full sphere was possible using a fraction of the available features. For learning HRTFs in listening tests, we developed an active-learning method for query-selection using GP models. Both simulations with offline GP-SSL models and HRTFs recommended to real human listeners have shown large improvement in localization accuracy over non-individualized HRTFs.
Appendix A Matérn Product Integrals
Improper integrals in Eq. 16 have closed-formulations:
-  J. Blauert, Spatial hearing: the psychophysics of human sound localization. Cambridge, Massachusettes: MIT Press, 1997.
-  C. Cheng and G. Wakefield, “Introduction to head-related transfer functions (HRTFs): Representations of HRTFs in time, frequency, and space,” in Audio Engineering Society Convention 107, 1999.
-  A. Kulkarni and H. Colburn, “Role of spectral detail in sound-source localization,” Nature, vol. 396, no. 6713, pp. 747–749, 1998.
-  E. Wenzel, M. Arruda, D. Kistler, and F. Wightman, “Localization using nonindividualized head-related transfer functions,” JASA, vol. 94, p. 111, 1993.
-  G. Romigh, D. Brungart, R. Stern, and B. Simpson, “The role of spatial detail in sound-source localization: Impact on HRTF modeling and personalization.” in Proceedings of Meetings on Acoustics, vol. 19, 2013.
-  J. Hornstein, M. Lopes, J. Santos-Victor, and F. Lacerda, “Sound localization for humanoid robots-building audio-motor maps based on the HRTF,” in Intelligent Robots and Systems, 2006 IEEE/RSJ International Conference on. IEEE, 2006, pp. 1170–1176.
-  M. Rothbucher, D. Kronmüller, M. Durkovic, T. Habigt, and K. Diepold, “HRTF sound localization,” 2011.
-  T. Rodemann, M. Heckmann, F. Joublin, C. Goerick, and B. Scholling, “Real-time sound localization with a binaural head-system using a biologically-inspired cue-triple mapping,” in International Conference on Intelligent Robots and Systems. IEEE, 2006, pp. 860–865.
-  H. Nakashima and T. Mukai, “3D sound source localization system based on learning of binaural hearing,” in Systems, Man and Cybernetics, 2005 IEEE International Conference on, vol. 4. IEEE, 2005.
-  V. Willert, J. Eggert, J. Adamy, R. Stahl, and E. Körner, “A probabilistic model for binaural sound localization,” IEEE Transactions on Systems, Man, and Cybernetics–Part B: Cybernetics, vol. 36, no. 5, p. 1, 2006.
-  A. Deleforge and R. Horaud, “2D sound-source localization on the binaural manifold,” in Machine Learning for Signal Processing (MLSP), 2012 IEEE International Workshop on. IEEE, 2012, pp. 1–6.
-  F. Keyrouz, K. Diepold, and S. Keyrouz, “High performance 3D sound localization for surveillance applications,” in Advanced Video and Signal Based Surveillance, 2007. AVSS 2007. IEEE Conference on. IEEE, 2007, pp. 563–566.
-  F. Keyrouz, “Humanoid hearing: A novel three-dimensional approach,” in Robotic and Sensors Environments (ROSE), 2011 IEEE International Symposium on. IEEE, 2011, pp. 214–219.
-  F. Keyrouz and K. Diepold, “An enhanced binaural 3D sound localization algorithm,” in Signal Processing and Information Technology, 2006 IEEE International Symposium on. IEEE, 2006, pp. 662–665.
-  A. Pourmohammad and S. Ahadi, “TDE-ILD-HRTF-Based 3D entire-space sound source localization using only three microphones and source counting,” in Electrical Engineering and Informatics (ICEEI), 2011 International Conference on. IEEE, 2011, pp. 1–6.
-  C. E. Rasmussen and C. Williams, Gaussian Processes for Machine Learning. Cambridge, Massachusettes: MIT Press, 2006.
-  A. J. Smola and B. Schölkopf, “A tutorial on support vector regression,” Statistics and computing, vol. 14, no. 3, pp. 199–222, 2004.
-  Y. Luo, D. N. Zotkin, and R. Duraiswami, “Statistical analysis of head related transfer function (HRTF) data,” in International Congress on Acoustics, 2013.
-  ——, “Gaussian process data fusion for heterogeneous HRTF datasets,” in WASPAA, 2013.
-  ——, “Gaussian process models for HRTF based 3D sound localization,” in ICASSP, 2014.
-  V. R. Algazi, R. O. Duda, and C. Avendano, “The CIPIC HRTF Database,” in IEEE Workshop on Applications of Signal Processing to Audio and Acoustics, New Paltz, NY, 2001, pp. 99–102.
-  I. Guyon and A. Elisseeff, “An introduction to variable and feature selection,” Journal of Machine Learning Research, vol. 3, pp. 1157–1182, 2003.
-  D. Zotkin, J. Hwang, R. Duraiswaini, and L. S. Davis, “HRTF personalization using anthropometric measurements,” in Applications of Signal Processing to Audio and Acoustics, 2003 IEEE Workshop on. Ieee, 2003, pp. 157–160.
-  H. Hu, L. Zhou, H. Ma, and Z. Wu, “HRTF personalization based on artificial neural network in individual virtual auditory space,” Applied Acoustics, vol. 69, no. 2, pp. 163–172, 2008.
-  Q. Huang and Y. Fang, “Modeling personalized head-related impulse response using support vector regression,” J Shanghai Univ (Engl Ed), vol. 13, no. 6, pp. 428–432, 2009.
-  K. Fink and L. Ray, “Tuning principal component weights to individualize HRTFs,” in ICASSP, 2012.
-  A. Silzle, “Selection and tuning of HRTFs,” in Audio Engineering Society Convention 112. Audio Engineering Society, 2002.
-  P. Runkle, A. Yendiki, and G. Wakefield, “Active sensory tuning for immersive spatialized audio,” in Proc. ICAD, 2000.
-  Y. Luo, D. N. Zotkin, and R. Duraiswami, “Virtual autoencoder based recommendation system for individualizing head-related transfer functions,” in WASPAA, 2013.
-  B. Settles, “Active learning literature survey,” University of Wisconsin, Madison, 2010.
-  M. Osborne, R. Garnett, and S. Roberts, “Gaussian processes for global optimization,” in 3rd International Conference on Learning and Intelligent Optimization (LION3), 2009, pp. 1–15.
-  N. Srinivas, A. Krause, S. M. Kakade, and M. Seeger, “Gaussian process optimization in the bandit setting: No regret and experimental design,” arXiv preprint arXiv:0912.3995, 2009.
-  R. Woodworth and G. Schlosberg, Experimental psychology. Holt, Rinehard and Winston, 1962.
-  C. Williams and M. Seeger, “Using the Nyström method to speed up kernel machines,” in Advances in Neural Information Processing Systems, 2000.
-  S. Mika, B. Schölkopf, A. J. Smola, K.-R. Müller, M. Scholz, and G. Rätsch, “Kernel PCA and de-noising in feature spaces.” in NIPS, vol. 11, 1998, pp. 536–542.
-  M. Seeger, C. Williams, and N. Lawrence, “Fast forward selection to speed up sparse gaussian process regression,” in Artificial Intelligence and Statistics 9, no. EPFL-CONF-161318, 2003.
-  R. Saigal, “On the inverse of a matrix with several rank one updates,” University of Michigan Ann Arbor, Tech. Rep., 1993.