Release notes#

1.4.0 (2026-09-30)#

A point marked Upgrading says what a script written for 1.3.3, or the environment it runs in, may have to change.

The warning with which a single-point gpyreg.GP.update() falls back to a full recomputation of a posterior names the line that called update, as it did up to 1.3.2; in 1.3.3 it named a line of gpyreg.
gpyreg.GP takes the keyword raise_on_cholesky_failure, False by default. With True, a Cholesky factorization of the training covariance that fails raises numpy.linalg.LinAlgError at its first attempt, in the posteriors and in the objective of gpyreg.GP.fit() (gpyreg.GP.log_likelihood() and gpyreg.GP.log_posterior()), as in MATLAB BADS (CholAttempts = 0). By default such a factorization is tried again with the noise multiplied tenfold, up to ten times, and the posterior keeps the multiplier, which its predictions apply to the noise while gpyreg.GP.get_hyperparameters() returns the noise as fitted. With the switch on, the space-filling design of fit ranks a starting point whose factorization fails last, as one of infinite value, and starts no optimization from it unless every point failed; a failure in the optimization, in the sampling or at the fitted hyperparameters raises from fit, which leaves the GP as it was. A copy or a pickle of the GP keeps the switch. Without it nothing changes, to the last bit.
gpyreg.covariance_functions.SquaredExponential, gpyreg.covariance_functions.Matern and gpyreg.covariance_functions.RationalQuadraticARD take the keyword periods: one period per input dimension, which makes the kernel periodic along that dimension, or np.inf for a dimension that is not periodic. Along a dimension of period p, the squared difference delta**2 of two inputs is replaced by the squared chord (p / pi)**2 * sin(pi * delta / p)**2, the squared distance between the two points mapped onto a circle of circumference p: to the kernel, coordinates a whole number of periods apart are the same, and the length scale keeps the units of the input. The periods are fixed constants, not hyperparameters, so the hyperparameters and their recommended bounds are those of the kernel without periods. A kernel without periods, as by default, or with periods that are all infinite, computes what it computed in 1.3.3, to the last bit. The isotropic kernels (gpyreg.isotropic_covariance_functions.MaternIsotropic and gpyreg.isotropic_covariance_functions.SquaredExponentialIsotropic) refuse periods, and so does gpyreg.GP.quad(), whose Gaussian integrals assume the non-periodic squared exponential kernel; gpyreg.GP refuses a covariance function whose number of periods is not the dimension of the GP.
A Gaussian prior whose centre lies some 38 of its scales or more outside the bounds of its hyperparameter has a finite log prior, and the space-filling design of gpyreg.GP.fit() draws finite values of that hyperparameter inside its bounds. The prior is renormalized by its mass inside the bounds, the difference of its cumulative distribution function (or of its survival function) at the two bounds, and so far in one tail both values underflowed to zero, and the mass with them: the log prior was infinite at every hyperparameter, so that gpyreg.GP.log_posterior() was infinite and the objective of the fit minus infinity, and the design mapped its draws through the same zero values, to infinity. Where the mass underflows to zero its log is taken in log space, and a draw that would be infinite comes from the prior truncated to the bounds; every other computation is unchanged, to the last bit.
The gradients of gpyreg.covariance_functions.RationalQuadraticARD and gpyreg.covariance_functions.Matern with respect to the length scales compute the factor that does not depend on the dimension once, where they computed it again for each dimension. With its gradient, the rational-quadratic kernel takes about a third less time at 10 dimensions and 150 inputs and a quarter less at 6 dimensions and 110 inputs, and the Matern kernels 15 to 40% less at those sizes, which speeds up gpyreg.GP.fit() with them. Every value is unchanged, to the last bit.
The ARD kernels (gpyreg.covariance_functions.SquaredExponential, gpyreg.covariance_functions.Matern and gpyreg.covariance_functions.RationalQuadraticARD) compute their gradient without a new array for each intermediate result, with the squared differences of every dimension taken at once where no dimension is periodic, and the squared exponential and rational-quadratic kernels compute their values the same way; gpyreg.GP.predict() scales and squares its intermediate arrays in place; the objective of gpyreg.GP.fit() solves its triangular systems by direct calls of LAPACK and sums its gradient in one array; and the masses of Gaussian priors inside their bounds, and the starting points that fit draws from them, take the standard normal’s functions from scipy.special rather than through scipy.stats. With one BLAS thread, at 3 to 10 dimensions and 50 to 150 inputs, with SciPy 1.17 and 1.18, the kernels with their gradient take 35 to 55% less time, the squared exponential kernel between the inputs and 2048 test points 55 to 65% less and the rational-quadratic one 45 to 55% less (10 to 25% at the smallest size), predict at 2048 test points 25 to 35% less (5 to 15% at 50 inputs), fit 8 to 20% less, and gpyreg.GP.set_bounds(), which computes those masses, a fourteenth to a twentieth of the time. Every value is unchanged, to the last bit, with any SciPy up to 1.18, for inputs of any floating type, infinite coordinates included.
gpyreg imports the modules of SciPy that it uses, scipy.stats among them. With SciPy older than 1.9, which loads a module only where it is imported, the methods of gpyreg.GP that use scipy.stats, gpyreg.GP.fit() among them, and gpyreg.f_min_fill raised AttributeError unless the caller had imported it.
The recommended bounds of the kernels and of gpyreg.mean_functions.NegativeQuadratic (their get_bounds_info) are computed without the warnings of NumPy (RuntimeWarning) that they raised on training inputs with a column without spread, as every column of a single input is: on a logarithm of zero, and on the standard deviation of a single input. gpyreg.GP.fit() computes these bounds even where the caller gives every bound, so that a fit on such inputs printed the warnings whatever its bounds. Every value is unchanged, to the last bit: the bounds built from a column without spread are -inf, which gpyreg.GP.get_recommended_bounds() refuses unless the caller gives a finite lower bound.
Installing gpyreg no longer installs pytest, pytest-rerunfailures and numdifftools, which only its tests use. The extra test installs them (pip install "gpyreg[test]"), and the extra dev includes them. Upgrading: an environment that runs gpyreg’s tests, or that relied on gpyreg to install one of these packages, installs the extra test or the packages themselves.

1.3.3 (2026-09-24)#

A point marked Upgrading says what a script written for 1.3.2 may have to change.

update(), update(hyp=...), gpyreg.GP.set_hyperparameters(), fit(X, y) and fit() run on a GP whose s2 holds a number, as a caller may assign it, which the noise function adds at every input (an update given new data needs s2 as an array, as in every release). The checks that a GP holds as many noise variances as inputs read the number of rows of s2 and raised AttributeError or IndexError for such a value; a value without rows is not counted, and the calls run as they did in 1.3.1, to the last bit.
A call of gpyreg.GP.fit() or gpyreg.GP.update(), and of gpyreg.GP.set_hyperparameters() through update, that raises leaves the GP as it was before the call, with the training data, bounds, priors and posteriors that it held; fit also refuses a column of inputs without spread, as gpyreg.GP.get_recommended_bounds() does, before it changes anything. Both stored the data they were given first, and fit then the bounds filled from them, so that a call that raised after that left the GP holding the data of the call without posteriors that match them, from which gpyreg.GP.predict() raised, or mixed the new inputs with the old posteriors: a fit that raised on such a column, on a starting point that holds NaN (which a fit without a space-filling design or an optimization keeps, and whose posterior update refuses), on a failed factorization of the training covariance, or in the sampler, after the optimization; and an update whose factorization of the training covariance failed, which left the posteriors replaced by empty ones or, where it failed in the full recomputation of a posterior to which a single-point update falls back, the other posteriors already extended. A fit that raised after filling its bounds also left those bounds, which a later fit that takes the bounds of the GP, as it does by default, kept. Every fit and update that completes is unchanged, to the last bit. Upgrading: a caller that reads a GP after a failed fit or update finds it as it was before the call, where 1.3.2 left it holding the new data (and, after a failed fit, the bounds the fit had set) with no usable posterior; a script that goes on with what the failed call stored gives it to the GP itself: the data through update or the next fit, the bounds through gpyreg.GP.set_bounds().
gpyreg.GP.update(), and gpyreg.GP.set_hyperparameters() through it, refuse hyperparameters that are NaN (not set), where the update would compute the posteriors from them, before it changes anything. The update stored the data it was given and replaced the posteriors with empty ones first, so that a GP given a NaN hyperparameter kept no posterior, and gpyreg.GP.predict() then raised AttributeError, and a GP without hyperparameters given data held them after the refusal. Upgrading: a script that gives data to a GP without hyperparameters through update and catches this refusal passes compute_posterior=False, with which the update stores the data without computing a posterior.

1.3.2 (2026-09-24)#

A point marked Upgrading says what a script written for 1.3.1 may have to change.

gpyreg.GP.get_recommended_bounds(), through which gpyreg.GP.fit() fills the bounds that the caller leaves unset, refuses with ValueError training inputs without spread in a column, as a single training point has none in any column, where a hyperparameter whose recommended bounds take their scale from that column has no finite lower bound from the caller, and the message names the columns and the hyperparameters concerned. The recommended bounds of a length scale of the kernel, and of the scale of gpyreg.mean_functions.NegativeQuadratic, take their scale from the width of the inputs, and a column of width zero gave them the pair (-inf, -inf). With the lower bound of such a hyperparameter left unset or given as -inf, the optimizer of fit ended with KeyError where the upper bound was unset or -inf as well, as it did with a lower bound of +inf. Beside a finite upper bound or one of +inf, the fit returned -inf for the logarithm of the scale, a length scale (or scale of the mean) of zero, from which its predictions were NaN, except on some fits to a single training point, which ended at the upper bound, or at NaN where that was +inf and then raised ValueError; with hyperparameter samples, its slice sampler raised ValueError. A finite lower bound from the caller is taken, and an upper bound left unset collapses onto it, as before: the fit runs as it did, to the last bit. Upgrading: a script that catches the KeyError catches ValueError, and a script that fits on such inputs, or reads the recommended bounds for them, gives the hyperparameters that the message names a finite lower bound first.
gpyreg.GP.set_priors() leaves the GP as it was when it refuses its argument, and marks the GP as having priors only where a coordinate of some block has one. On a GP without priors, a refused call marked the GP as having priors, and so did a call whose blocks had no prior in any coordinate, such as ("student_t", (nan, nan, nan)), so that str reported them as present and gpyreg.GP.fit() added to its objective their log prior, which is zero except at a point off the value of a hyperparameter that equal bounds fix, where it is -inf. A fit without a space-filling design (init_N=0) ranks its starting points, those it is given or those of the GP, before it moves them into the bounds: it ranked such a point last, and ranks it by its likelihood, and where that point comes first, the fit starts from it, moved onto the fixed value, and its result changes. Upgrading: a script that reads from str whether a GP set with such blocks alone has priors finds none, and one that compares the result of such a fit with values stored from 1.3.1 stores them again.
gpyreg.GP.get_priors() returns every prior that gpyreg.GP.set_priors() takes in a form that set_priors writes back as it was, so set_priors(get_priors()) changes nothing, the mark of the GP as having priors or none included. A Student’s t block whose degrees of freedom mix zero with NaN or a number, such as [0, nan] or [0, 3], came back as None, which dropped the prior; degrees of freedom infinite throughout, or zero on the coordinates that have a prior and NaN on the others, came back as the Gaussian family, which set_priors writes with zero; and a family set on a block with no prior in any coordinate came back as None. Such a block, whose location and sigma are NaN throughout, comes back under the Gaussian or the Student’s t family that its degrees of freedom name: as "gaussian" where it was set with "gaussian" or "smoothbox", or with "student_t" or "smoothbox_student_t" and zero degrees of freedom throughout; as None where it was set with one of the latter two and NaN degrees of freedom throughout; and as "student_t" where it was set with one of them and other degrees of freedom. get_priors raises ValueError, with set_priors’ message, for priors that set_priors refuses, as priors written into hyper_priors directly can be, and those of a GP pickled by gpyreg 1.2.x or earlier, whose set_priors took any sigma, where it returned a block that set_priors then refused (for a negative sigma, among others), or None (for a NaN sigma beside a location). Upgrading: a script that compares what get_priors returns finds a Student’s t family with infinite degrees of freedom, or with degrees of freedom of zero where it has a prior and NaN elsewhere, under its own name, not as the Gaussian family, a Student’s t block whose degrees of freedom mix zero with NaN or a number as "student_t", not None, and a block without a prior in any coordinate, unless its degrees of freedom are NaN throughout, under the family just given, not None; a script that writes priors into hyper_priors directly sets them through set_priors instead; and one that calls get_priors on a GP pickled by gpyreg 1.2.x or earlier whose priors set_priors refuses sets the priors of that GP again with set_priors first.
gpyreg.GP.set_bounds() refuses with ValueError a lower bound above the upper bound of the same hyperparameter, naming it, as gpyreg.GP.get_recommended_bounds() and gpyreg.GP.fit() already did, and leaves the bounds as they were. It stored the inverted pair, which the next fit refused where it took both its lower and its upper bounds from the GP, as it does by default, and ignored where it was given its own lower_bounds and upper_bounds. Upgrading: a script that set an inverted pair, and never fitted or fitted with bounds of its own, gives the pair in order.
gpyreg.GP.fit() takes its option thin as a whole number of an integer or a float type, as gpyreg.slice_sample.SliceSampler.sample() takes it, where a whole float such as 2.0 raised TypeError from the fit’s own use of it. It refuses with ValueError, before it changes anything, a thin that is not a whole number greater than zero: a fraction, which raised TypeError, or the sampler’s ValueError on the burn-in where the default burn, thin times n_samples, is a fraction; zero or a negative number, which raised another ValueError after the optimization; a bool, which ran as the integer it stands for; and, in a fit without hyperparameter samples, which does not use it, any such value, which it took. Upgrading: a script that passes thin=True passes 1; one that catches the TypeError of a fractional thin catches ValueError; and one that gives a fit with n_samples=0 a thin that is not a whole number greater than zero leaves it out.
gpyreg.GP.fit() takes its other counts, n_samples, opts_N and init_N, as whole numbers of an integer or a float type, where a whole float such as 2.0 raised TypeError, and refuses with ValueError, before it changes anything, one that is not a whole number of at least zero: a fraction, which raised TypeError or another ValueError after the optimization; a negative number, with which opts_N and init_N ran as zero and n_samples raised after the optimization; a NaN init_N, which ran as zero; and a bool, which ran as the integer it stands for, except init_N=True, which raised TypeError. gpyreg.slice_sample.SliceSampler.sample() takes its number of samples N the same way, and refuses with ValueError one that is not a whole number greater than zero, where a whole float, a fraction or a bool raised TypeError and zero returned an empty chain; it refuses a bool as thin or burn as well, which ran as the integer it stands for. The counts of fit and of sample, thin and burn included, take a 0-d array that holds a whole number as that number, and refuse one that holds anything else as they refuse what it holds. An array holding an integer ran as that integer, except as thin or burn of sample, and so as burn of fit, which sample refused; one holding a whole float raised TypeError, or that refusal. Upgrading: a script that passes a negative opts_N or init_N, or a NaN init_N, passes 0; one that passes a bool as a count passes the integer; one that catches the TypeError of a count that is not a whole number, such as opts_N=1.5 or init_N=32.5 in fit or N=5.5 or N=True in sample, catches ValueError; and one that asks sample for zero samples does not call it.
gpyreg.GP.fit() checks its option burn with its other counts, before it changes anything and whether or not it draws hyperparameter samples, by the rule of gpyreg.slice_sample.SliceSampler.sample(): a whole number of at least zero, taken as the other counts are, or None, which leaves the burn-in to the sampler. It refuses any other value with ValueError. The fit passed burn to the sampler unchecked: with samples, the sampler refused such a value with ValueError after the optimization, raised TypeError for a string, and ran a bool as the integer it stands for; without samples, the fit took any value. Upgrading: a script that catches the TypeError of a string burn in a fit that draws samples catches ValueError, and one that gives a fit with n_samples=0 a burn that is neither None nor a whole number of at least zero leaves it out.
Where the smallest noise variance at the training inputs is below 1e-6, as it can be after a fit on noiseless targets, the posterior holds the negative inverse of the training covariance (the low-noise representation). gpyreg.GP.predict() and gpyreg.GP.predict_full() form the predictive covariance there from a Cholesky factor of the training covariance, which the posterior keeps beside the inverse (the attribute L_factor of gpyreg.gaussian_process.Posterior, one more N by N matrix per hyperparameter sample). Formed from the inverse, as gplite forms it, the covariance carried a rounding error that grows as the noise shrinks: with 30 training inputs on [-2, 2], a squared exponential kernel of unit length and output scales and a noise standard deviation of 1e-6, the variances at the training inputs, of order 1e-12, were off by about 1e-4, many of them clamped to zero, and the covariance of predict_full had eigenvalues of order -1e-3. A single-point gpyreg.GP.update() of such a posterior takes from the factor the predictive variance of the new point and the solve that it divides by that variance, and extends the factor. With the variance formed from the inverse it recomputed the posterior in full where rounding had clamped that variance, and otherwise extended the posterior with a wrong one, which moved the predictive mean away from that of the posterior computed in full: over the last ten of the 30 inputs of the GP above, added one at a time, with targets of range 2 and noise standard deviations from 1e-7 to 1e-4, by up to 1e-2, and by up to 0.3 with a length scale of 3, as the noise and the rounding of the BLAS build decided. gpyreg.GP.random_function() draws from the kept factor instead of factoring the training covariance at each call: its draws from a posterior that no single-point update has extended are unchanged, and after such updates they change with the posterior. A posterior pickled by an earlier version has the factor computed again where it is needed. The Cholesky representation is unchanged, to the last bit. Upgrading: the predictions of a GP in the low-noise representation, and its posterior after single-point updates, change; a script that compares them with values stored from 1.3.1 stores them again. PyBADS reaches this representation: at its default options the lower bound of the noise of its GP is a variance of about 1.4e-7, below the 1e-6 at which the representation starts. It makes no single-point updates and reads only the hyperparameters of the posteriors, so for any fit whose noise variance ends below 1e-6, as on a noiseless target, its predictive means are unchanged and its predictive variances change, which can take its search elsewhere.
gpyreg.GP.quad() with compute_var=True forms the variance of an integral from the same Cholesky factor in the low-noise representation, where it formed it from the inverse and carried its rounding: with two training points 1e-3 apart, a squared exponential kernel of unit length and output scales, a noise standard deviation of 1e-6 and a Gaussian measure of standard deviation 0.5 centred between the two points, a variance of 0.016 was off by 6e-11. The means of the integrals, and both in the Cholesky representation, are unchanged, to the last bit. Upgrading: the variances of integrals of a GP in the low-noise representation change.
gpyreg.GP.fit() on a GP without training data, given neither to it nor held from an earlier fit or update, raises ValueError that says so and names what is missing, before it changes anything. It raised AttributeError from a covariance function or from the check of the input shapes, or, given X alone, the ValueError of gpyreg.GP.get_recommended_bounds() after storing X. Upgrading: a script that catches the AttributeError catches ValueError.
gpyreg.GP.update() and gpyreg.GP.fit() refuse with ValueError, before they change anything and with a message that gives the numbers, a call that would make the number of targets, or of noise variances, that the GP holds differ from its number of inputs. update refuses targets y_new or noise variances s2_new given without the inputs X_new to a GP that holds no inputs, or whose inputs have their targets, or their variances; X_new without y_new to a GP that holds targets; and X_new with y_new to a GP that holds inputs without targets. fit refuses X given without y, or without s2, where the targets, or the variances, that the GP holds are one per input it holds and not one per row of X; it counts only the variances that the noise function reads (as gpyreg.noise_functions.GaussianNoise with user_provided_add does), and keeps the others as before. On a GP without inputs, update raised AttributeError. Elsewhere both stored the data, and the computation of the posterior, or the objective of the fit, then raised ValueError and left the GP without posteriors that match its data, so that the next gpyreg.GP.predict() raised. They raised nothing, and left the numbers different, where no posterior was computed (compute_posterior=False, or a GP without targets), where only noise variances that the noise function does not take differed, as after s2_new given alone to data that hold their variances, and where NumPy broadcast a single value against the others: on a GP of a single training point, and where one input and its target were given to a GP that holds inputs without targets, whose posterior took that target at every input. Targets given alone, one per input, to a GP that holds inputs without targets, and noise variances given alone to one that holds data without variances, are taken as before: no call that ran to its end and left the GP with as many targets, and as many noise variances, as inputs, or none of them, is refused. Upgrading: a script that catches the AttributeError catches ValueError, and one that relied on a call that left the numbers different (without a posterior, in variances that the noise function does not take, or through such a broadcast) gives the targets and the variances with their inputs, one per input.
Documentation: the Raises section of gpyreg.GP.fit() names the ValueError it passes on from gpyreg.GP.get_recommended_bounds(), from the check of the shapes of its training data and from gpyreg.slice_sample.SliceSampler, and the ValueError and LinAlgError it passes on from gpyreg.GP.update(). Those of gpyreg.GP.predict(), gpyreg.GP.predict_full(), gpyreg.GP.quad() and gpyreg.GP.random_function() name the LinAlgError of the factorization that a posterior pickled by an earlier version in the low-noise representation needs. Those of update, predict and predict_full name the ValueError of the check of the shapes of the data they are given. Those of fit, update, predict and predict_full name the TypeError of a noise variance that is neither an array, a number nor None, and that of predict the ValueError of return_lpd=True without y_star.

1.3.1 (2026-09-23)#

A point marked Upgrading says what a script written for 1.3.0 may have to change.

The gradient of the log prior is zero, not NaN, for a hyperparameter whose lower and upper bounds are equal and that has no prior, or a smooth-box prior whose box holds its value. The NaN reached gpyreg.GP.log_posterior() with compute_grad=True, and, where another hyperparameter has a prior, the optimizer of gpyreg.GP.fit(), which stopped within an iteration, short of the optimum.
The probability that a hyperprior puts inside the bounds of its hyperparameter, by which gpyreg.GP.log_posterior() renormalizes the prior, is computed from the survival function where both bounds lie above the centre of the prior. As a difference of two values of the cumulative distribution function it was zero with both bounds far in the upper tail (beyond about 8.3 scales of a Gaussian prior), which made the log posterior infinite everywhere and sent gpyreg.GP.fit() to a poor point. Where the lower bound is not above the centre, the value is the same as before, to the last bit. gpyreg.f_min_fill has the survival functions of the two smooth-box families, smoothbox_sf and smoothbox_student_t_sf.
The space-filling design of gpyreg.GP.fit(), drawn by gpyreg.f_min_fill, draws the starting values of a hyperparameter whose bounds both lie above the centre of its prior through the survival function of the prior and its inverse, where it drew them through the cumulative distribution function and the percent point function. With both bounds far in the upper tail, every draw of that hyperparameter lay at infinity (the starting points the fit was given stayed finite), and a fit with hyperparameter samples raised ValueError because the widths of the slice sampler were NaN. Where the lower bound is not above the centre, the design is the same as before, to the last bit. gpyreg.f_min_fill has the inverse survival functions of the two smooth-box families, smoothbox_isf and smoothbox_student_t_isf.
The space-filling design of gpyreg.GP.fit(), drawn by gpyreg.f_min_fill, gives a hyperparameter whose lower and upper bounds are equal their value at every point when the hyperparameter has a prior, as it did for one without. Mapped through the prior’s quantile function, the value came back an ulp or two off, where the log prior is -inf, so that the objective was infinite at every point of the design but the starting points the fit was given, and the ranking of the design that picks the starts of the optimization was lost. A GP that PyVBMC builds meets this when the noise is held at its lower bound (targets of a range below about 3e-3) and the option noise_size moves the centre of the noise prior away from that bound.
gpyreg.GP.set_priors() refuses, with a message that says what is wrong, a coordinate whose sigma is finite beside a location that is not: an infinite or NaN mu of a Gaussian or Student’s t prior, or an infinite or NaN end a or b of a smooth box, as it refuses a sigma that is not finite and positive. 1.3.0 took such a prior, and the log posterior was NaN with the bounds that gpyreg.GP.fit() fills. A coordinate of a block without a prior is written, as before, with NaN for both its location and its sigma. Upgrading: a script that wrote a non-finite mu for no prior, which gplite reads that way, sets the prior of a hyperparameter none of whose coordinates has a prior to None, and gives a coordinate of a block without a prior NaN for both its location and its sigma, as the message says. A smooth box with an infinite or NaN end has no such replacement: the script writes the prior it means.
gpyreg.GP.set_priors() refuses, with a message that says what is wrong, a smooth box of either family whose lower end a is above its upper end b. 1.3.0 took such a box, whose normalization constant is then below one or negative, and gave a log prior that was wrong or NaN. A smooth box with a == b, which has no plateau and is the Gaussian or the Student’s t centred at a, is taken as before. Upgrading: a script that wrote an inverted box writes the box it means, with a <= b.

1.3.0 (2026-09-23)#

A point marked Upgrading says what a script written for 1.2.1 may have to change.

The recommended bounds of the hyperparameters take the statistics of the training inputs per input dimension. The starting value of each length scale is the logarithm of the standard deviation of its own dimension, and the location and the scale of gpyreg.mean_functions.NegativeQuadratic have, per dimension, a hard box, a plausible box and a starting value built from that dimension’s minimum, maximum, median, width and standard deviation. They were taken over all the entries of X at once, one number for every dimension, which mixes the locations of the dimensions with their spreads: where the dimensions lie apart, every recommendation came out as wide as the widest gap between them. The per-dimension statistics are those of MATLAB’s gplite. With one input dimension the two agree, and nothing changes. In more than one, a gpyreg.GP.fit() with that mean function takes its plausible box from the recommendation, and its hard box wherever the caller leaves it unset, so its space-filling design, its default starting point and its slice sampler change, and its hyperparameters move; a fit that starts from given hyperparameters (hyp0, or those the GP already has), with no space-filling design (init_N = 0) and no hyperparameter samples (n_samples = 0), changes only where its start or its optimizer meets a hard bound that moved. fit does not start from the recommended length scales, which reach a caller that reads them from get_bounds_info; the bounds of the length scales were per dimension already. PyVBMC 1.0.4 is such a caller: it starts the length scales of its GP from the recommendations for the points of highest density of its training set, and where those points share one value in a column, the start of that length scale is -inf, or tens below zero where rounding leaves their standard deviation in it above zero, and from -inf the first GP fit of the run fails, where it completed with gpyreg 1.0.3 to 1.2.1; PyVBMC 1.5 takes the statistics of such a column from the whole training set.
The recommended bounds and starting length scale of the isotropic kernels take the means of the logarithms of the per-dimension widths and standard deviations, as MATLAB’s gplite does, where they took the logarithm of the mean width and the standard deviation pooled over all entries of X. With one input dimension the two agree; in more than one, a fit with gpyreg.isotropic_covariance_functions.SquaredExponentialIsotropic or gpyreg.isotropic_covariance_functions.MaternIsotropic takes a different plausible box for the length scale than in 1.2.1, and a different hard box wherever the caller leaves it unset.
gpyreg.GP.fit() takes its starting points as a copy of the space-filling design, so the default widths of the slice sampler are the standard deviation of the design as it was drawn; a fit that draws hyperparameter samples with a noise hyperparameter and 1 < opts_N < init_N, as at the defaults, sees its samples move.
The plausible upper bound of the shape parameter of gpyreg.covariance_functions.RationalQuadraticARD is its own (5) and the output scale keeps the range of the targets; the shape’s line wrote into the output scale’s slot. The plausible box of a fit with that kernel changes, and with it the space-filling design and the default starting point.
Slice sampling with a burn-in of two or three iterations no longer returns the same point once per requested sample: a coordinate whose burn-in variance estimate is not positive keeps the width it has instead of being floored to zero, and a window of fewer than two iterations adapts nothing.
An infinite entry of SliceSampler(widths=...), which the documentation of widths allows for an unbounded coordinate, no longer comes back after the burn-in and fills the chain with NaN.
gpyreg.covariance_functions.Matern and gpyreg.isotropic_covariance_functions.MaternIsotropic of degree 1 can be fitted: the gradient with respect to the length scale is zero, not NaN, where two inputs coincide, so the gradient of the marginal likelihood is finite.
A training set whose targets are all equal and finite is given a range of one, with a warning, instead of bounds of -inf that ended gpyreg.GP.fit() with KeyError: (-inf, -inf); such a fit now completes.
gpyreg.GP.quad() with compute_var=True normalizes by the noise scale the stored Cholesky factor carries, so the variance of an integral is right after a rank-one gpyreg.GP.update() that appended a point of lower total noise, where it used to be clamped to machine epsilon.
gpyreg.GP.random_function() draws a sample where the predictive covariance is numerically singular, as it is on a dense one-dimensional grid: eigenvalues that are negative but of rounding size count as the zeros they are, where the draw used to collapse onto the predictive mean without a word. The rounding is measured against the prior variance at the test points, from which the predictive covariance is formed by subtraction, so a grid inside the training data draws as well. This holds in both representations of the posterior. Where the smallest noise variance at the training inputs is below 1e-6, as it can be after a default fit on noiseless targets, the posterior holds the inverse of the training covariance instead of its Cholesky factor (the low-noise representation). The predictive covariance that 1.2.1 formed from that inverse carried a rounding error that grows as the noise shrinks, so that the draws of 1.2.1 from such a posterior often came out as the predictive mean, inside the training data and outside it; the draw now forms the covariance from a Cholesky factor of the training covariance. It raises LinAlgError where the covariance it computes has a negative eigenvalue beyond the rounding of the prior variance, and, in the low-noise representation, where the Cholesky decomposition of the training covariance fails even with the noise raised. The eigenvalue fallback also fixes the signs of whole eigenvectors and uses the symmetric eigensolver, so the factor it builds is a factor of the matrix it was given; a draw at closely spaced or duplicated test points came from the wrong covariance before. Draws through the eigenvalue fallback, and draws from a posterior in the low-noise representation, change for a given generator.
gpyreg.GP.predict_full() with add_noise=True adds the observation noise on the diagonal; with a noise that varies from point to point the returned matrix was neither symmetric nor a covariance matrix, while its diagonal was already right. A constant noise function is unaffected.
A hyperprior whose degrees of freedom are infinite or NaN is the Gaussian family it names, MATLAB’s HPRIOR.nu = Inf convention, where it used to contribute no prior at all; gpyreg.GP.fit() fills NaN degrees of freedom with its option df_base for its own duration (below). A smooth-box prior set on a block of several hyperparameters has one normalization constant per coordinate of the block, where it doubled the log density with two coordinates on different sides of the box and raised with three. gpyreg.GP.get_priors() returns, in the form gpyreg.GP.set_priors() reads back, such a smooth-box prior, where it raised, and a prior with NaN degrees of freedom or a block with a coordinate that has no prior, where it returned None; so set_priors(get_priors()) keeps them.
A single-observation gpyreg.GP.update() of a posterior in the low-noise representation falls through to a full recomputation, with a warning, where rounding drives the latent variance of the new point to zero or below, as it can for an observation at an existing training input, and its predictive variance is the noise alone; the Cholesky representation already did so for its own stability test.
gpyreg.GP.fit() collapses an inverted plausible pair onto its upper bound, inside the hard box, before the space-filling design, which needs the pair ordered. For training targets whose standard deviation is below 1e-3 (the noise’s plausible lower bound) and whose range is above 1e-6 (below that the hard pair collapses first), the noise’s recommended plausible pair is inverted, and in 1.2.1 a fit with a space-filling design whose noise has no prior raised AssertionError.
The convergence diagnostics of gpyreg.slice_sample.SliceSampler recognize a parameter whose recorded chain did not move by its range, whatever value it is frozen at: R and eff_N are NaN for it, and for a free parameter exit_flag is -3, with the “did not move” message. Before, this held only for a constant whose mean is exact, and a parameter fixed by LB == UB at, say, 0.3 reported finite R and eff_N that were rounding noise.
The normalization constant of a smooth-box Student’s t prior, and the cumulative distribution and quantile functions of that distribution, are computed through the logarithms of the gamma functions, so degrees of freedom above about 340 no longer make the log posterior NaN.
gpyreg.GP.log_likelihood() and gpyreg.GP.log_posterior() accept the dictionary of hyperparameters their docstrings document.
gpyreg.GP.fit() reads the sampler under sampler_name, the name it documents, as well as under sampler, and fills df_base into a copy of the priors for its own duration, so a fitted GP keeps the priors the caller set and a second fit with another value uses it. A prior left with NaN degrees of freedom therefore reads as Gaussian outside the fit, in gpyreg.GP.log_posterior() after it for one, where 1.2.1 wrote df_base into the GP’s priors and read a Student’s t from then on. Upgrading: a script that relies on that Student’s t after the fit gives the prior its degrees of freedom.
A single-observation gpyreg.GP.update() on a GP that carries no posterior factors, after gpyreg.GP.clean() or after an update with compute_posterior=False, recomputes the posterior in full instead of raising TypeError.
The gradient of the marginal likelihood no longer raises where the noise function carries a scale for a user-provided variance that is never given.
gpyreg.GP.quad() reads a one-dimensional mu or sigma of length D as one measure in D dimensions, as MATLAB’s gplite does, where 1.2.1 raised.
Inputs that 1.2.1 took in silence are refused with a message: a coordinate that has a prior needs a finite, positive sigma (a hyperparameter without a prior is written as None, and a coordinate of a block without one as NaN location and sigma, as 1.2.1 accepted it); gpyreg.GP.set_priors() and gpyreg.GP.set_bounds() reject a hyperparameter name the model has not; gpyreg.GP.get_recommended_bounds() takes any array_like and rejects a bound pair given inverted, and so does gpyreg.GP.fit(), which fills its bounds through it; gpyreg.GP.update() with hyp checks the width of the hyperparameter row, and with compute_posterior=True refuses hyperparameters that are NaN (never set), naming them; gpyreg.GP.quad() refuses a mean function it cannot place, and a mu with more columns than input dimensions, for which 1.2.1 ignored the extra columns, broadcast them into a wrong integral or raised a broadcasting error, depending on the mean function and the dimension; and a noise variance given as an array whose row count is not that of the inputs is refused, where 1.2.1 reshaped any array of N entries into a column and failed on any other. Inputs on which 1.2.1 failed with an error from inside the computation are refused with a message that names the problem: gpyreg.GP.quad() refuses a GP without training data or posterior factors, where 1.2.1 raised AttributeError or TypeError, and a mu with fewer columns than input dimensions, where it raised IndexError; and gpyreg.slice_sample.SliceSampler.sample() refuses a thin or burn that is not a whole number, an infinite one included, with ValueError instead of raising TypeError from range. A noise variance may be any number or 0-d array; quad takes a sigma of one column, one standard deviation per measure in every dimension, as 1.2.1 and gplite take it, and refuses a sigma of any other width than one or the number of input dimensions; and SliceSampler.sample takes a thin or burn that is a whole number of any type, a float such as 2.0 included, where 1.2.1 took integers only. A failed Cholesky decomposition reports LinAlgError in both noise representations, where the low-noise one raised TypeError, and the checks of the shape of the inputs raise ValueError where they raised AssertionError. Upgrading: each refusal of an input that 1.2.1 took in silence can stop a script that ran under 1.2.1, and so can the new type of an exception where a script catches the old one; such a script changes as follows:
- a name the model has not, which set nothing, is left out; PyVBMC 1.0.4 is such a caller: at its uncertainty level 1 (uncertainty_handling on, and a target that does not return its own noise estimate) it sets a prior on noise_provided_log_multiplier, a hyperparameter its GP does not have, so that every such run stops at its first GP fit; it needs gpyreg 1.2.1 or earlier, or PyVBMC 1.5;
- a prior whose sigma is infinite, zero or negative, or NaN beside a location, is given a positive sigma, or is replaced by None where no prior is meant (1.2.1 took the absolute value of a negative sigma in the log prior);
- a bound pair given inverted, which 1.2.1 collapsed onto its lower bound, is given in order (equal bounds fix a hyperparameter);
- update(hyp=...) is given one column per hyperparameter of the GP, where 1.2.1 stored a row of another width and read it by offset;
- gp.update(X_new=X, y_new=y) on a GP whose hyperparameters were never set, which 1.2.1 completed with NaN posterior factors, gives the data to gp.fit(X, y) instead. A following gp.fit() found hyperparameters in 1.2.1 only where it drew no hyperparameter samples (n_samples=0) and started from a space-filling design (init_N above zero, as at the default); for such a fit, compute_posterior=False in the update keeps the pattern working;
- quad with another mean function, whose integral 1.2.1 computed as if that mean were a constant equal to its first hyperparameter, has no replacement;
- quad with a mu of more columns than input dimensions, whose extra columns 1.2.1 ignored or broadcast into a wrong integral, is given one column per dimension;
- a noise variance given as a row of N is given as a column;
- a script that catches the exception of a check whose type changed catches the new one: ValueError for a shape check (an X that is not two-dimensional, or whose number of columns is not the GP’s D), which raised AssertionError, for quad on a GP without training data or posterior factors (AttributeError or TypeError) or with a mu of fewer columns than input dimensions (IndexError), and for a thin or burn that is not a whole number (TypeError); LinAlgError for a failed Cholesky decomposition in the low-noise representation (TypeError).
The covariance kernels refuse compute(compute_diag=True, compute_grad=True) with ValueError. Upgrading: that combination returned the diagonal beside the gradient of the full matrix, a pair that meant nothing; ask for the two separately.
The undocumented Metropolis step of gpyreg.slice_sample.SliceSampler is removed, with its options and the attributes metropolis_pdf, metropolis_rnd and metropolis_flag. Upgrading: its options never turned the step on (their key was misspelled), but setting the three attributes directly did run it; a script that sets them now samples without the step, with no error.
uuinv returns NaN for a p outside [0, 1] in every case, not only in the general one, and its documentation states the mixture it implements: the weight outside the plausible box is spread over the two tails in proportion to their lengths.
Documentation: the isotropic kernels have a page; AbstractKernel.compute documents the (N, 1) shape of the diagonal it returns; gpyreg.noise_functions.GaussianNoise documents the np.spacing(1) nugget it keeps with constant_add=False; the default of add_noise of gpyreg.GP.predict_full(), that its diagonal is not clamped and can be negative on a nearly singular posterior; which variance gpyreg.GP.predict() returns by default, and how it pools the variances of several hyperparameter samples; that the log predictive density always carries the observation noise and over several hyperparameter samples is the density of the Gaussian with the pooled mean and variance; that gpyreg.GP.log_posterior() renormalizes each prior over its bounds; what NaN degrees of freedom of a prior mean, in gpyreg.GP.set_priors() and in the df_base option of gpyreg.GP.fit(); and, in the docstring of the sampler’s effective sample size, that its autocorrelations are paired from lag 0.

1.2.1 (2026-09-14)#

gpyreg.GP.quad() returns the correct integral variance when the training noise is heteroskedastic, for example with user-provided observation variances. The variance previously collapsed to machine epsilon in that case. quad also handles SquaredExponentialIsotropic in more than one dimension, where it previously read the hyperparameters in the wrong order. Results with homoskedastic noise, and with the isotropic kernel in one dimension, are unchanged.
gpyreg.GP.update() applies replacement hyperparameters passed together with a single new observation. The rank-one shortcut is taken only when the hyperparameters stay in place; it previously discarded the new values. Calls without hyp are unchanged.
The rank-one update of a single appended observation extends the Cholesky factor with the noise scale the factor was built with, stored as the gpyreg.gaussian_process.Posterior attribute sl. It is therefore exact for heteroskedastic noise, where it previously assumed the new point’s noise equal to the smallest training noise, and it works with output-dependent noise, which previously raised an error. The shortcut now also applies when s2_new is supplied, which previously triggered a full recomputation; the result is the same up to floating-point rounding. The homoskedastic case is unchanged.
gpyreg.GP.update() keeps the stored user-provided noise aligned with the training inputs. Points added without s2_new to a GP that stores s2, or earlier points of a GP that receives s2_new for the first time, get a variance of zero. Previously s2 could end up with fewer rows than X, which made the next posterior computation fail, and in the second case a single new variance was silently applied to every training point.
gpyreg.slice_sample.SliceSampler accepts bounds and widths given as lists. A coordinate with equal list bounds is now recognized as fixed; previously it went undetected and list widths raised an error. The sampler also logs its no-violation message when the diagnostics report success, which previously never happened.
gpyreg.slice_sample.SliceSampler with step_out=True in more than one dimension evaluates the stepping-out brackets on the current coordinate line. Previously the bracket ends of a later coordinate carried an earlier coordinate’s shrunk bracket edge. Sampling with the default step_out=False is unchanged.
The convergence diagnostics of gpyreg.slice_sample.SliceSampler report a positive, bounded effective sample size, following Geyer’s initial positive sequence with the floor used by Stan and ArviZ; it could previously be negative. A coordinate fixed by LB == UB is excluded from the checks, and undefined diagnostics or fewer than eight recorded samples give exit_flag -3 instead of reporting success. R and the effective sample size are unchanged for chains that mix.
The documentation gains a section on BLAS threading, and the attributes of gpyreg.gaussian_process.Posterior are documented.

1.2.0 (2026-09-13)#

gpyreg.GP.predict() accepts the keyword-only option return_cross_covariance=True. It appends a tuple containing the latent training-to-prediction kernel matrix for each hyperparameter sample. Downstream calculations can reuse these matrices instead of computing them again. Ordinary calls retain their existing return values.
Kernel matrices remain separate when predictions are averaged across hyperparameter samples. They exclude observation noise and are unconditioned kernels, rather than posterior predictive covariances. Prior-only GPs return None for each sample. See Gaussian processes for examples, shapes and memory considerations.
The prediction reference corrects the documented default of add_noise to False and the averaged-output shape to (M, 1). These are documentation corrections; the corresponding behavior is unchanged.

1.1.0 (2026-09-05)#

gpyreg.GP.fit(), gpyreg.GP.random_function(), gpyreg.slice_sample.SliceSampler and f_min_fill accept rng=. Pass a NumPy generator to share a stream with the caller, or a seed to create a generator. The default rng=None continues to use NumPy’s global stream controlled by np.random.seed. See Random number generation.
Prediction and hyperparameter fitting have lower computational overhead, including batched evaluation of bundled mean functions and reuse of a factorization when only mean hyperparameters change during fitting. Custom mean and covariance implementations retain their supported paths.
gpyreg.GP.log_likelihood() and gpyreg.GP.log_posterior() correctly return (value, gradient) when compute_grad=True; these calls previously raised TypeError.
Slice samplers using the default stream remain copyable and picklable. Their saved state does not capture NumPy’s global random state. Samplers with an explicit generator preserve its state when copied or pickled; older sampler saves resume using the global stream.

Release notes

Contents