gpmp.num module¶

The gpmp.num module is the numerical backend layer used by GPmp. It dispatches a NumPy-like API to the backend selected for the current Python process.

Backend object contract¶

GPmp functions operate on backend-native numerical objects:

With GPMP_BACKEND=numpy, backend arrays are NumPy ndarray objects.
With GPMP_BACKEND=torch, backend arrays are PyTorch Tensor objects.

Inputs passed to GPmp should be backend objects or objects convertible by gpmp.num.asarray. Outputs are backend objects unless a function explicitly documents conversion to NumPy. Use gpmp.num.to_np at plotting, reporting, or external-library boundaries when a NumPy array is required.

When writing custom mean functions, covariance functions, or selection criteria, use gpmp.num operations instead of direct numpy or torch calls when the code must be backend-independent.

Backend selection and dtype¶

The backend is selected at import time. Set GPMP_BACKEND to "numpy" or "torch" before importing GPmp, or call gpmp.config.set_backend(...) before importing gpmp.num. If no backend is requested, GPmp uses PyTorch when available and otherwise falls back to NumPy.

GPmp uses float64 only. The portable dtype is available through gpmp.num.get_dtype() and backend-specific constructors default to that dtype where applicable.

Core constructors and conversion¶

Numerical backend dispatcher for GPmp.

GPmp code uses backend-native numerical objects through this module. With the numpy backend, arrays are NumPy arrays. With the torch backend, arrays are PyTorch tensors. Public GPmp functions accept objects convertible by gpmp.num.asarray and may return backend-native objects unless their own docstring states that they convert outputs to NumPy.

Use this module, usually imported as gpmp.num as gnp, when writing custom mean functions, covariance functions, criteria, or examples that should work with both supported backends.

gpmp.num.asarray(x, dtype=None)[source]¶

gpmp.num.array(x, dtype=None)[source]¶

gpmp.num.zeros(shape, dtype=None)[source]¶

gpmp.num.ones(shape, dtype=None)[source]¶

gpmp.num.full(shape, fill_value, dtype=None)[source]¶

gpmp.num.eye(n, m=None, k=0, dtype=None)[source]¶

gpmp.num.to_np(x)[source]¶

gpmp.num.to_scalar(x)[source]¶

gpmp.num.get_dtype()[source]¶

Mathematical operations¶

These wrappers use NumPy-style argument names where possible. In particular, reductions use axis and ddof even when the torch backend is active.

gpmp.num.exp(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])¶

Calculate the exponential of all elements in the input array.

Parameters:

x (array_like) – Input values.
out (ndarray, None, or tuple of ndarray and None, optional) – A location into which the result is stored. If provided, it must have a shape that the inputs broadcast to. If not provided or None, a freshly-allocated array is returned. A tuple (possible only as a keyword argument) must have length equal to the number of outputs.
where (array_like, optional) – This condition is broadcast over the input. At locations where the condition is True, the out array will be set to the ufunc result. Elsewhere, the out array will retain its original value. Note that if an uninitialized out array is created via the default out=None, locations within it where the condition is False will remain uninitialized.
**kwargs – For other keyword-only arguments, see the ufunc docs.

Returns:

out – Output array, element-wise exponential of x. This is a scalar if x is a scalar.

Return type:

ndarray or scalar

See also

expm1: Calculate exp(x) - 1 for all elements in the array.
exp2: Calculate 2**x for all elements in the array.

Notes

The irrational number e is also known as Euler’s number. It is approximately 2.718281, and is the base of the natural logarithm, ln (this means that, if \(x = \ln y = \log_e y\), then \(e^x = y\). For real input, exp(x) is always positive.

For complex arguments, x = a + ib, we can write \(e^x = e^a e^{ib}\). The first term, \(e^a\), is already known (it is the real argument, described above). The second term, \(e^{ib}\), is \(\cos b + i \sin b\), a function with magnitude 1 and a periodic phase.

References

[1]

Wikipedia, “Exponential function”, https://en.wikipedia.org/wiki/Exponential_function

[2]

M. Abramovitz and I. A. Stegun, “Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables,” Dover, 1964, p. 69, https://personal.math.ubc.ca/~cbm/aands/page_69.htm

Examples

Plot the magnitude and phase of exp(x) in the complex plane:

>>> import numpy as np

>>> import matplotlib.pyplot as plt
>>> import numpy as np

>>> x = np.linspace(-2*np.pi, 2*np.pi, 100)
>>> xx = x + 1j * x[:, np.newaxis] # a + ib over complex plane
>>> out = np.exp(xx)

>>> plt.subplot(121)
>>> plt.imshow(np.abs(out),
...            extent=[-2*np.pi, 2*np.pi, -2*np.pi, 2*np.pi], cmap='gray')
>>> plt.title('Magnitude of exp(x)')

>>> plt.subplot(122)
>>> plt.imshow(np.angle(out),
...            extent=[-2*np.pi, 2*np.pi, -2*np.pi, 2*np.pi], cmap='hsv')
>>> plt.title('Phase (angle) of exp(x)')
>>> plt.show()

gpmp.num.log(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])¶

Natural logarithm, element-wise.

The natural logarithm log is the inverse of the exponential function, so that log(exp(x)) = x. The natural logarithm is logarithm in base e.

Parameters:

x (array_like) – Input value.
out (ndarray, None, or tuple of ndarray and None, optional) – A location into which the result is stored. If provided, it must have a shape that the inputs broadcast to. If not provided or None, a freshly-allocated array is returned. A tuple (possible only as a keyword argument) must have length equal to the number of outputs.
where (array_like, optional) – This condition is broadcast over the input. At locations where the condition is True, the out array will be set to the ufunc result. Elsewhere, the out array will retain its original value. Note that if an uninitialized out array is created via the default out=None, locations within it where the condition is False will remain uninitialized.
**kwargs – For other keyword-only arguments, see the ufunc docs.

Returns:

y – The natural logarithm of x, element-wise. This is a scalar if x is a scalar.

Return type:

ndarray

See also

log10, log2, log1p, emath.log

Notes

Logarithm is a multivalued function: for each x there is an infinite number of z such that exp(z) = x. The convention is to return the z whose imaginary part lies in (-pi, pi].

For real-valued input data types, log always returns real output. For each value that cannot be expressed as a real number or infinity, it yields nan and sets the invalid floating point error flag.

For complex-valued input, log is a complex analytical function that has a branch cut [-inf, 0] and is continuous from above on it. log handles the floating-point negative zero as an infinitesimal negative number, conforming to the C99 standard.

In the cases where the input has a negative real part and a very small negative complex part (approaching 0), the result is so close to -pi that it evaluates to exactly -pi.

References

[1]

M. Abramowitz and I.A. Stegun, “Handbook of Mathematical Functions”, 10th printing, 1964, pp. 67. https://personal.math.ubc.ca/~cbm/aands/page_67.htm

[2]

Wikipedia, “Logarithm”. https://en.wikipedia.org/wiki/Logarithm

Examples

>>> import numpy as np
>>> np.log([1, np.e, np.e**2, 0])
array([  0.,   1.,   2., -inf])

gpmp.num.sqrt(x, /, out=None, *, where=True, casting='same_kind', order='K', dtype=None, subok=True[, signature])¶

Return the non-negative square-root of an array, element-wise.

Parameters:

x (array_like) – The values whose square-roots are required.
out (ndarray, None, or tuple of ndarray and None, optional) – A location into which the result is stored. If provided, it must have a shape that the inputs broadcast to. If not provided or None, a freshly-allocated array is returned. A tuple (possible only as a keyword argument) must have length equal to the number of outputs.
where (array_like, optional) – This condition is broadcast over the input. At locations where the condition is True, the out array will be set to the ufunc result. Elsewhere, the out array will retain its original value. Note that if an uninitialized out array is created via the default out=None, locations within it where the condition is False will remain uninitialized.
**kwargs – For other keyword-only arguments, see the ufunc docs.

Returns:

y – An array of the same shape as x, containing the positive square-root of each element in x. If any element in x is complex, a complex array is returned (and the square-roots of negative reals are calculated). If all of the elements in x are real, so is y, with negative elements returning nan. If out was provided, y is a reference to it. This is a scalar if x is a scalar.

Return type:

ndarray

See also

emath.sqrt: A version which returns complex numbers when given negative reals. Note that 0.0 and -0.0 are handled differently for complex inputs.

Notes

sqrt has–consistent with common convention–as its branch cut the real “interval” [-inf, 0), and is continuous from above on it. A branch cut is a curve in the complex plane across which a given complex function fails to be continuous.

Examples

>>> import numpy as np
>>> np.sqrt([1,4,9])
array([ 1.,  2.,  3.])

>>> np.sqrt([4, -1, -3+4J])
array([ 2.+0.j,  0.+1.j,  1.+2.j])

>>> np.sqrt([4, -1, np.inf])
array([ 2., nan, inf])

gpmp.num.sum(a, axis=None, dtype=None, out=None, keepdims=<no value>, initial=<no value>, where=<no value>)¶

Sum of array elements over a given axis.

Parameters:

a (array_like) – Elements to sum.
axis (None or int or tuple of ints, optional) – Axis or axes along which a sum is performed. The default, axis=None, will sum all of the elements of the input array. If axis is negative it counts from the last to the first axis. If axis is a tuple of ints, a sum is performed on all of the axes specified in the tuple instead of a single axis or all the axes as before.
dtype (dtype, optional) – The type of the returned array and of the accumulator in which the elements are summed. The dtype of a is used by default unless a has an integer dtype of less precision than the default platform integer. In that case, if a is signed then the platform integer is used while if a is unsigned then an unsigned integer of the same precision as the platform integer is used.
out (ndarray, optional) – Alternative output array in which to place the result. It must have the same shape as the expected output, but the type of the output values will be cast if necessary.
keepdims (bool, optional) –
If this is set to True, the axes which are reduced are left in the result as dimensions with size one. With this option, the result will broadcast correctly against the input array.

If the default value is passed, then keepdims will not be passed through to the sum method of sub-classes of ndarray, however any non-default value will be. If the sub-class’ method does not implement keepdims any exceptions will be raised.
initial (scalar, optional) – Starting value for the sum. See ~numpy.ufunc.reduce for details.
where (array_like of bool, optional) – Elements to include in the sum. See ~numpy.ufunc.reduce for details.

Returns:

sum_along_axis – An array with the same shape as a, with the specified axis removed. If a is a 0-d array, or if axis is None, a scalar is returned. If an output array is specified, a reference to out is returned.

Return type:

ndarray

See also

ndarray.sum: Equivalent method.
add: numpy.add.reduce equivalent function.
cumsum: Cumulative sum of array elements.
trapezoid: Integration of array values using composite trapezoidal rule.

mean, average

Notes

Arithmetic is modular when using integer types, and no error is raised on overflow.

The sum of an empty array is the neutral element 0:

>>> np.sum([])
0.0

For floating point numbers the numerical precision of sum (and np.add.reduce) is in general limited by directly adding each number individually to the result causing rounding errors in every step. However, often numpy will use a numerically better approach (partial pairwise summation) leading to improved precision in many use-cases. This improved precision is always provided when no axis is given. When axis is given, it will depend on which axis is summed. Technically, to provide the best speed possible, the improved precision is only used when the summation is along the fast axis in memory. Note that the exact precision may vary depending on other parameters. In contrast to NumPy, Python’s math.fsum function uses a slower but more precise approach to summation. Especially when summing a large number of lower precision floating point numbers, such as float32, numerical errors can become significant. In such cases it can be advisable to use dtype=”float64” to use a higher precision for the output.

Examples

>>> import numpy as np
>>> np.sum([0.5, 1.5])
2.0
>>> np.sum([0.5, 0.7, 0.2, 1.5], dtype=np.int32)
np.int32(1)
>>> np.sum([[0, 1], [0, 5]])
6
>>> np.sum([[0, 1], [0, 5]], axis=0)
array([0, 6])
>>> np.sum([[0, 1], [0, 5]], axis=1)
array([1, 5])
>>> np.sum([[0, 1], [np.nan, 5]], where=[False, True], axis=1)
array([1., 5.])

If the accumulator is too small, overflow occurs:

>>> np.ones(128, dtype=np.int8).sum(dtype=np.int8)
np.int8(-128)

You can also start the sum with a value other than zero:

>>> np.sum([10], initial=5)
15

gpmp.num.mean(a, axis=None, dtype=None, out=None, keepdims=<no value>, *, where=<no value>)¶

Compute the arithmetic mean along the specified axis.

Returns the average of the array elements. The average is taken over the flattened array by default, otherwise over the specified axis. float64 intermediate and return values are used for integer inputs.

Parameters:

a (array_like) – Array containing numbers whose mean is desired. If a is not an array, a conversion is attempted.
axis (None or int or tuple of ints, optional) –
Axis or axes along which the means are computed. The default is to compute the mean of the flattened array.

If this is a tuple of ints, a mean is performed over multiple axes, instead of a single axis or all the axes as before.
dtype (data-type, optional) – Type to use in computing the mean. For integer inputs, the default is float64; for floating point inputs, it is the same as the input dtype.
out (ndarray, optional) – Alternate output array in which to place the result. The default is None; if provided, it must have the same shape as the expected output, but the type will be cast if necessary. See ufuncs-output-type for more details. See ufuncs-output-type for more details.
keepdims (bool, optional) –
If this is set to True, the axes which are reduced are left in the result as dimensions with size one. With this option, the result will broadcast correctly against the input array.

If the default value is passed, then keepdims will not be passed through to the mean method of sub-classes of ndarray, however any non-default value will be. If the sub-class’ method does not implement keepdims any exceptions will be raised.
where (array_like of bool, optional) –
Elements to include in the mean. See ~numpy.ufunc.reduce for details.

Added in version 1.20.0.

Returns:

m – If out=None, returns a new array containing the mean values, otherwise a reference to the output array is returned.

Return type:

ndarray, see dtype parameter above

See also

average: Weighted average

std, var, nanmean, nanstd, nanvar

Notes

The arithmetic mean is the sum of the elements along the axis divided by the number of elements.

Note that for floating-point input, the mean is computed using the same precision the input has. Depending on the input data, this can cause the results to be inaccurate, especially for float32 (see example below). Specifying a higher-precision accumulator using the dtype keyword can alleviate this issue.

By default, float16 results are computed using float32 intermediates for extra precision.

Examples

>>> import numpy as np
>>> a = np.array([[1, 2], [3, 4]])
>>> np.mean(a)
2.5
>>> np.mean(a, axis=0)
array([2., 3.])
>>> np.mean(a, axis=1)
array([1.5, 3.5])

In single precision, mean can be inaccurate:

>>> a = np.zeros((2, 512*512), dtype=np.float32)
>>> a[0, :] = 1.0
>>> a[1, :] = 0.1
>>> np.mean(a)
np.float32(0.54999924)

Computing the mean in float64 is more accurate:

>>> np.mean(a, dtype=np.float64)
0.55000000074505806 # may vary

Computing the mean in timedelta64 is available:

>>> b = np.array([1, 3], dtype="timedelta64[D]")
>>> np.mean(b)
np.timedelta64(2,'D')

Specifying a where argument:

>>> a = np.array([[5, 9, 13], [14, 10, 12], [11, 15, 19]])
>>> np.mean(a)
12.0
>>> np.mean(a, where=[[True], [False], [False]])
9.0

gpmp.num.var(a, axis=None, dtype=None, out=None, ddof=0, keepdims=<no value>, *, where=<no value>, mean=<no value>, correction=<no value>)¶

Compute the variance along the specified axis.

Returns the variance of the array elements, a measure of the spread of a distribution. The variance is computed for the flattened array by default, otherwise over the specified axis.

Parameters:

a (array_like) – Array containing numbers whose variance is desired. If a is not an array, a conversion is attempted.
axis (None or int or tuple of ints, optional) – Axis or axes along which the variance is computed. The default is to compute the variance of the flattened array. If this is a tuple of ints, a variance is performed over multiple axes, instead of a single axis or all the axes as before.
dtype (data-type, optional) – Type to use in computing the variance. For arrays of integer type the default is float64; for arrays of float types it is the same as the array type.
out (ndarray, optional) – Alternate output array in which to place the result. It must have the same shape as the expected output, but the type is cast if necessary.
ddof ({int, float}, optional) – “Delta Degrees of Freedom”: the divisor used in the calculation is N - ddof, where N represents the number of elements. By default ddof is zero. See notes for details about use of ddof.
keepdims (bool, optional) –
If this is set to True, the axes which are reduced are left in the result as dimensions with size one. With this option, the result will broadcast correctly against the input array.

If the default value is passed, then keepdims will not be passed through to the var method of sub-classes of ndarray, however any non-default value will be. If the sub-class’ method does not implement keepdims any exceptions will be raised.
where (array_like of bool, optional) –
Elements to include in the variance. See ~numpy.ufunc.reduce for details.

Added in version 1.20.0.
mean (array like, optional) –
Provide the mean to prevent its recalculation. The mean should have a shape as if it was calculated with keepdims=True. The axis for the calculation of the mean should be the same as used in the call to this var function.

Added in version 2.0.0.
correction ({int, float}, optional) –
Array API compatible name for the ddof parameter. Only one of them can be provided at the same time.

Added in version 2.0.0.

Returns:

variance – If out=None, returns a new array containing the variance; otherwise, a reference to the output array is returned.

Return type:

ndarray, see dtype parameter above

See also

std, mean, nanmean, nanstd, nanvar, ufuncs-output-type

Notes

There are several common variants of the array variance calculation. Assuming the input a is a one-dimensional NumPy array and mean is either provided as an argument or computed as a.mean(), NumPy computes the variance of an array as:

N = len(a)
d2 = abs(a - mean)**2  # abs is for complex `a`
var = d2.sum() / (N - ddof)  # note use of `ddof`

Different values of the argument ddof are useful in different contexts. NumPy’s default ddof=0 corresponds with the expression:

\[\frac{\sum_i{|a_i - \bar{a}|^2 }}{N}\]

which is sometimes called the “population variance” in the field of statistics because it applies the definition of variance to a as if a were a complete population of possible observations.

Many other libraries define the variance of an array differently, e.g.:

\[\frac{\sum_i{|a_i - \bar{a}|^2}}{N - 1}\]

In statistics, the resulting quantity is sometimes called the “sample variance” because if a is a random sample from a larger population, this calculation provides an unbiased estimate of the variance of the population. The use of \(N-1\) in the denominator is often called “Bessel’s correction” because it corrects for bias (toward lower values) in the variance estimate introduced when the sample mean of a is used in place of the true mean of the population. For this quantity, use ddof=1.

Note that for complex numbers, the absolute value is taken before squaring, so that the result is always real and nonnegative.

For floating-point input, the variance is computed using the same precision the input has. Depending on the input data, this can cause the results to be inaccurate, especially for float32 (see example below). Specifying a higher-accuracy accumulator using the dtype keyword can alleviate this issue.

Examples

>>> import numpy as np
>>> a = np.array([[1, 2], [3, 4]])
>>> np.var(a)
1.25
>>> np.var(a, axis=0)
array([1.,  1.])
>>> np.var(a, axis=1)
array([0.25,  0.25])

In single precision, var() can be inaccurate:

>>> a = np.zeros((2, 512*512), dtype=np.float32)
>>> a[0, :] = 1.0
>>> a[1, :] = 0.1
>>> np.var(a)
np.float32(0.20250003)

Computing the variance in float64 is more accurate:

>>> np.var(a, dtype=np.float64)
0.20249999932944759 # may vary
>>> ((1-0.55)**2 + (0.1-0.55)**2)/2
0.2025

Specifying a where argument:

>>> a = np.array([[14, 8, 11, 10], [7, 9, 10, 11], [10, 15, 5, 10]])
>>> np.var(a)
6.833333333333333 # may vary
>>> np.var(a, where=[[True], [True], [False]])
4.0

Using the mean keyword to save computation time:

>>> import numpy as np
>>> from timeit import timeit
>>>
>>> a = np.array([[14, 8, 11, 10], [7, 9, 10, 11], [10, 15, 5, 10]])
>>> mean = np.mean(a, axis=1, keepdims=True)
>>>
>>> g = globals()
>>> n = 10000
>>> t1 = timeit("var = np.var(a, axis=1, mean=mean)", globals=g, number=n)
>>> t2 = timeit("var = np.var(a, axis=1)", globals=g, number=n)
>>> print(f'Percentage execution time saved {100*(t2-t1)/t2:.0f}%')

Percentage execution time saved 32%

gpmp.num.std(a, axis=None, dtype=None, out=None, ddof=0, keepdims=<no value>, *, where=<no value>, mean=<no value>, correction=<no value>)¶

Compute the standard deviation along the specified axis.

Returns the standard deviation, a measure of the spread of a distribution, of the array elements. The standard deviation is computed for the flattened array by default, otherwise over the specified axis.

Parameters:

a (array_like) – Calculate the standard deviation of these values.
axis (None or int or tuple of ints, optional) – Axis or axes along which the standard deviation is computed. The default is to compute the standard deviation of the flattened array. If this is a tuple of ints, a standard deviation is performed over multiple axes, instead of a single axis or all the axes as before.
dtype (dtype, optional) – Type to use in computing the standard deviation. For arrays of integer type the default is float64, for arrays of float types it is the same as the array type.
out (ndarray, optional) – Alternative output array in which to place the result. It must have the same shape as the expected output but the type (of the calculated values) will be cast if necessary. See ufuncs-output-type for more details.
ddof ({int, float}, optional) – Means Delta Degrees of Freedom. The divisor used in calculations is N - ddof, where N represents the number of elements. By default ddof is zero. See Notes for details about use of ddof.
keepdims (bool, optional) –
If this is set to True, the axes which are reduced are left in the result as dimensions with size one. With this option, the result will broadcast correctly against the input array.

If the default value is passed, then keepdims will not be passed through to the std method of sub-classes of ndarray, however any non-default value will be. If the sub-class’ method does not implement keepdims any exceptions will be raised.
where (array_like of bool, optional) –
Elements to include in the standard deviation. See ~numpy.ufunc.reduce for details.

Added in version 1.20.0.
mean (array_like, optional) –
Provide the mean to prevent its recalculation. The mean should have a shape as if it was calculated with keepdims=True. The axis for the calculation of the mean should be the same as used in the call to this std function.

Added in version 2.0.0.
correction ({int, float}, optional) –
Array API compatible name for the ddof parameter. Only one of them can be provided at the same time.

Added in version 2.0.0.

Returns:

standard_deviation – If out is None, return a new array containing the standard deviation, otherwise return a reference to the output array.

Return type:

ndarray, see dtype parameter above.

See also

var, mean, nanmean, nanstd, nanvar, ufuncs-output-type

Notes

There are several common variants of the array standard deviation calculation. Assuming the input a is a one-dimensional NumPy array and mean is either provided as an argument or computed as a.mean(), NumPy computes the standard deviation of an array as:

N = len(a)
d2 = abs(a - mean)**2  # abs is for complex `a`
var = d2.sum() / (N - ddof)  # note use of `ddof`
std = var**0.5

Different values of the argument ddof are useful in different contexts. NumPy’s default ddof=0 corresponds with the expression:

\[\sqrt{\frac{\sum_i{|a_i - \bar{a}|^2 }}{N}}\]

which is sometimes called the “population standard deviation” in the field of statistics because it applies the definition of standard deviation to a as if a were a complete population of possible observations.

Many other libraries define the standard deviation of an array differently, e.g.:

\[\sqrt{\frac{\sum_i{|a_i - \bar{a}|^2 }}{N - 1}}\]

In statistics, the resulting quantity is sometimes called the “sample standard deviation” because if a is a random sample from a larger population, this calculation provides the square root of an unbiased estimate of the variance of the population. The use of \(N-1\) in the denominator is often called “Bessel’s correction” because it corrects for bias (toward lower values) in the variance estimate introduced when the sample mean of a is used in place of the true mean of the population. The resulting estimate of the standard deviation is still biased, but less than it would have been without the correction. For this quantity, use ddof=1.

Note that, for complex numbers, std takes the absolute value before squaring, so that the result is always real and nonnegative.

For floating-point input, the standard deviation is computed using the same precision the input has. Depending on the input data, this can cause the results to be inaccurate, especially for float32 (see example below). Specifying a higher-accuracy accumulator using the dtype keyword can alleviate this issue.

Examples

>>> import numpy as np
>>> a = np.array([[1, 2], [3, 4]])
>>> np.std(a)
1.1180339887498949 # may vary
>>> np.std(a, axis=0)
array([1.,  1.])
>>> np.std(a, axis=1)
array([0.5,  0.5])

In single precision, std() can be inaccurate:

>>> a = np.zeros((2, 512*512), dtype=np.float32)
>>> a[0, :] = 1.0
>>> a[1, :] = 0.1
>>> np.std(a)
np.float32(0.45000005)

Computing the standard deviation in float64 is more accurate:

>>> np.std(a, dtype=np.float64)
0.44999999925494177 # may vary

Specifying a where argument:

>>> a = np.array([[14, 8, 11, 10], [7, 9, 10, 11], [10, 15, 5, 10]])
>>> np.std(a)
2.614064523559687 # may vary
>>> np.std(a, where=[[True], [True], [False]])
2.0

Using the mean keyword to save computation time:

>>> import numpy as np
>>> from timeit import timeit
>>> a = np.array([[14, 8, 11, 10], [7, 9, 10, 11], [10, 15, 5, 10]])
>>> mean = np.mean(a, axis=1, keepdims=True)
>>>
>>> g = globals()
>>> n = 10000
>>> t1 = timeit("std = np.std(a, axis=1, mean=mean)", globals=g, number=n)
>>> t2 = timeit("std = np.std(a, axis=1)", globals=g, number=n)
>>> print(f'Percentage execution time saved {100*(t2-t1)/t2:.0f}%')

Percentage execution time saved 30%

gpmp.num.cov(m, y=None, rowvar=True, bias=False, ddof=None, fweights=None, aweights=None, *, dtype=None)¶

Estimate a covariance matrix, given data and weights.

Covariance indicates the level to which two variables vary together. If we examine N-dimensional samples, \(X = [x_1, x_2, ..., x_N]^T\), then the covariance matrix element \(C_{ij}\) is the covariance of \(x_i\) and \(x_j\). The element \(C_{ii}\) is the variance of \(x_i\).

See the notes for an outline of the algorithm.

Parameters:

m (array_like) – A 1-D or 2-D array containing multiple variables and observations. Each row of m represents a variable, and each column a single observation of all those variables. Also see rowvar below.
y (array_like, optional) – An additional set of variables and observations. y has the same form as that of m.
rowvar (bool, optional) – If rowvar is True (default), then each row represents a variable, with observations in the columns. Otherwise, the relationship is transposed: each column represents a variable, while the rows contain observations.
bias (bool, optional) – Default normalization (False) is by (N - 1), where N is the number of observations given (unbiased estimate). If bias is True, then normalization is by N. These values can be overridden by using the keyword ddof in numpy versions >= 1.5.
ddof (int, optional) – If not None the default value implied by bias is overridden. Note that ddof=1 will return the unbiased estimate, even if both fweights and aweights are specified, and ddof=0 will return the simple average. See the notes for the details. The default value is None.
fweights (array_like, int, optional) – 1-D array of integer frequency weights; the number of times each observation vector should be repeated.
aweights (array_like, optional) – 1-D array of observation vector weights. These relative weights are typically large for observations considered “important” and smaller for observations considered less “important”. If ddof=0 the array of weights can be used to assign probabilities to observation vectors.
dtype (data-type, optional) –
Data-type of the result. By default, the return data-type will have at least numpy.float64 precision.

Added in version 1.20.

Returns:

out – The covariance matrix of the variables.

Return type:

ndarray

See also

corrcoef: Normalized covariance matrix

Notes

Assume that the observations are in the columns of the observation array m and let f = fweights and a = aweights for brevity. The steps to compute the weighted covariance are as follows:

>>> m = np.arange(10, dtype=np.float64)
>>> f = np.arange(10) * 2
>>> a = np.arange(10) ** 2.
>>> ddof = 1
>>> w = f * a
>>> v1 = np.sum(w)
>>> v2 = np.sum(w * a)
>>> m -= np.sum(m * w, axis=None, keepdims=True) / v1
>>> cov = np.dot(m * w, m.T) * v1 / (v1**2 - ddof * v2)

Note that when a == 1, the normalization factor v1 / (v1**2 - ddof * v2) goes over to 1 / (np.sum(f) - ddof) as it should.

Examples

>>> import numpy as np

Consider two variables, \(x_0\) and \(x_1\), which correlate perfectly, but in opposite directions:

>>> x = np.array([[0, 2], [1, 1], [2, 0]]).T
>>> x
array([[0, 1, 2],
       [2, 1, 0]])

Note how \(x_0\) increases while \(x_1\) decreases. The covariance matrix shows this clearly:

>>> np.cov(x)
array([[ 1., -1.],
       [-1.,  1.]])

Note that element \(C_{0,1}\), which shows the correlation between \(x_0\) and \(x_1\), is negative.

Further, note how x and y are combined:

>>> x = [-2.1, -1,  4.3]
>>> y = [3,  1.1,  0.12]
>>> X = np.stack((x, y), axis=0)
>>> np.cov(X)
array([[11.71      , -4.286     ], # may vary
       [-4.286     ,  2.144133]])
>>> np.cov(x, y)
array([[11.71      , -4.286     ], # may vary
       [-4.286     ,  2.144133]])
>>> np.cov(x)
array(11.71)

Linear algebra and distances¶

gpmp.num.solve(a, b, lower=False, overwrite_a=False, overwrite_b=False, check_finite=True, assume_a=None, transposed=False)[source]¶

Solve the equation a @ x = b for x, where a is a square matrix.

If the data matrix is known to be a particular type then supplying the corresponding string to assume_a key chooses the dedicated solver. The available options are

diagonal	‘diagonal’
tridiagonal	‘tridiagonal’
banded	‘banded’
upper triangular	‘upper triangular’
lower triangular	‘lower triangular’
symmetric	‘symmetric’ (or ‘sym’)
hermitian	‘hermitian’ (or ‘her’)
symmetric positive definite	‘positive definite’ (or ‘pos’)
general	‘general’ (or ‘gen’)

Array argument(s) of this function may have additional “batch” dimensions prepended to the core shape. In this case, the array is treated as a batch of lower-dimensional slices; see linalg_batch for details.

Parameters:

a (array_like, shape (..., N, N)) – Square left-hand side matrix or a batch of matrices.
b ((..., N, NRHS) array_like) – Input data for the right hand side or a batch of right-hand sides.
lower (bool, default: False) – Ignored unless assume_a is one of 'sym', 'her', or 'pos'. If True, the calculation uses only the data in the lower triangle of a; entries above the diagonal are ignored. If False (default), the calculation uses only the data in the upper triangle of a; entries below the diagonal are ignored.
overwrite_a (bool, default: False) – Allow overwriting data in a (may enhance performance).
overwrite_b (bool, default: False) – Allow overwriting data in b (may enhance performance).
check_finite (bool, default: True) – Whether to check that the input matrices contain only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination) if the inputs do contain infinities or NaNs.
assume_a (str, optional) – Valid entries are described above. If omitted or None, checks are performed to identify structure so the appropriate solver can be called.
transposed (bool, default: False) – If True, solve a.T @ x == b. Raises NotImplementedError for complex a.

Returns:

x – The solution array.

Return type:

ndarray, shape (N, NRHS) or (…, N)

Raises:

ValueError – If size mismatches detected or input a is not square.
LinAlgError – If the computation fails because of matrix singularity.
LinAlgWarning – If an ill-conditioned input a is detected.
NotImplementedError – If transposed is True and input a is a complex matrix.

Notes

If the input b matrix is a 1-D array with N elements, when supplied together with an NxN input a, it is assumed as a valid column vector despite the apparent size mismatch. This is compatible with the numpy.dot() behavior and the returned result is still 1-D array.

The general, symmetric, Hermitian and positive definite solutions are obtained via calling ?GETRF/?GETRS, ?SYSV, ?HESV, and ?POTRF/?POTRS routines of LAPACK respectively.

The datatype of the arrays define which solver is called regardless of the values. In other words, even when the complex array entries have precisely zero imaginary parts, the complex solver will be called based on the data type of the array.

Examples

Given a and b, solve for x:

>>> import numpy as np
>>> a = np.array([[3, 2, 0], [1, -1, 0], [0, 5, 1]])
>>> b = np.array([2, 4, -1])
>>> from scipy.linalg import solve
>>> x = solve(a, b)
>>> x
array([ 2., -2.,  9.])
>>> a @ x == b
array([ True,  True,  True], dtype=bool)

Batches of matrices are supported, with and without structure detection:

>>> a = np.arange(12).reshape(3, 2, 2)   # a batch of 3 2x2 matrices
>>> A = a.transpose(0, 2, 1) @ a    # A is a batch of 3 positive definite matrices
>>> b = np.ones(2)
>>> solve(A, b)      # this automatically detects that A is pos.def.
array([[ 1. , -0.5],
       [ 3. , -2.5],
       [ 5. , -4.5]])
>>> solve(A, b, assume_a='pos')   # bypass structucture detection
array([[ 1. , -0.5],
       [ 3. , -2.5],
       [ 5. , -4.5]])

gpmp.num.solve_triangular(a, b, trans=0, lower=False, unit_diagonal=False, overwrite_b=False, check_finite=True)[source]¶

Solve the equation a @ x = b for x, where a is a triangular matrix.

The documentation is written assuming array arguments are of specified “core” shapes. However, array argument(s) of this function may have additional “batch” dimensions prepended to the core shape. In this case, the array is treated as a batch of lower-dimensional slices; see linalg_batch for details. Note that calls with zero-size batches are unsupported and will raise a ValueError.

Parameters:

a ((M, M) array_like) – A triangular matrix
b ((M,) or (M, N) array_like) – Right-hand side matrix in a x = b
lower (bool, optional) – Use only data contained in the lower triangle of a. Default is to use upper triangle.
trans ({0, 1, 2, 'N', 'T', 'C'}, optional) –
Type of system to solve:

trans

system

0 or ‘N’

a x = b

1 or ‘T’

a^T x = b

2 or ‘C’

a^H x = b
unit_diagonal (bool, optional) – If True, diagonal elements of a are assumed to be 1 and will not be referenced.
overwrite_b (bool, optional) – Allow overwriting data in b (may enhance performance)
check_finite (bool, optional) – Whether to check that the input matrices contain only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination) if the inputs do contain infinities or NaNs.

Returns:

trans	system
0 or ‘N’	a x = b
1 or ‘T’	a^T x = b
2 or ‘C’	a^H x = b

x – Solution to the system a x = b. Shape of return matches b.

Return type:

(M,) or (M, N) ndarray

Raises:

LinAlgError – If a is singular

Notes

Added in version 0.9.0.

Examples

Solve the lower triangular system a x = b, where:

     [3  0  0  0]       [4]
a =  [2  1  0  0]   b = [2]
     [1  0  1  0]       [4]
     [1  1  1  1]       [2]

>>> import numpy as np
>>> from scipy.linalg import solve_triangular
>>> a = np.array([[3, 0, 0, 0], [2, 1, 0, 0], [1, 0, 1, 0], [1, 1, 1, 1]])
>>> b = np.array([4, 2, 4, 2])
>>> x = solve_triangular(a, b, lower=True)
>>> x
array([ 1.33333333, -0.66666667,  2.66666667, -1.33333333])
>>> a.dot(x)  # Check the result
array([ 4.,  2.,  4.,  2.])

gpmp.num.cho_factor(a, lower=False, overwrite_a=False, check_finite=True)[source]¶

Compute the Cholesky decomposition of a matrix, to use in cho_solve

Returns a matrix containing the Cholesky decomposition, A = L L* or A = U* U of a Hermitian positive-definite matrix a. The return value can be directly used as the first parameter to cho_solve.

Warning

The returned matrix also contains random data in the entries not used by the Cholesky decomposition. If you need to zero these entries, use the function cholesky instead.

Parameters:

a ((M, M) array_like) – Matrix to be decomposed
lower (bool, optional) – Whether to compute the upper or lower triangular Cholesky factorization. During decomposition, only the selected half of the matrix is referenced. (Default: upper-triangular)
overwrite_a (bool, optional) – Whether to overwrite data in a (may improve performance)
check_finite (bool, optional) – Whether to check that the entire input matrix contains only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination) if the inputs do contain infinities or NaNs.

Returns:

c ((M, M) ndarray) – Matrix whose upper or lower triangle contains the Cholesky factor of a. Other parts of the matrix contain random data.
lower (bool) – Flag indicating whether the factor is in the lower or upper triangle

Raises:

LinAlgError – Raised if decomposition fails.

See also

cho_solve(): Solve a linear set equations using the Cholesky factorization of a matrix.

Notes

During the finiteness check (if selected), the entire matrix a is checked. During decomposition, a is assumed to be symmetric or Hermitian (as applicable), and only the half selected by option lower is referenced. Consequently, if a is asymmetric/non-Hermitian, cholesky may still succeed if the symmetric/Hermitian matrix represented by the selected half is positive definite, yet it may fail if an element in the other half is non-finite.

Examples

>>> import numpy as np
>>> from scipy.linalg import cho_factor
>>> A = np.array([[9, 3, 1, 5], [3, 7, 5, 1], [1, 5, 9, 2], [5, 1, 2, 6]])
>>> c, low = cho_factor(A)
>>> c
array([[3.        , 1.        , 0.33333333, 1.66666667],
       [3.        , 2.44948974, 1.90515869, -0.27216553],
       [1.        , 5.        , 2.29330749, 0.8559528 ],
       [5.        , 1.        , 2.        , 1.55418563]])
>>> np.allclose(np.triu(c).T @ np. triu(c) - A, np.zeros((4, 4)))
True

gpmp.num.cho_solve(c_and_lower, b, overwrite_b=False, check_finite=True)[source]¶

Solve the linear equations A x = b, given the Cholesky factorization of A.

Parameters:

(c (tuple, (array, bool)) – Cholesky factorization of a, as given by cho_factor
lower) (tuple, (array, bool)) – Cholesky factorization of a, as given by cho_factor
b (array) – Right-hand side
overwrite_b (bool, optional) – Whether to overwrite data in b (may improve performance)
check_finite (bool, optional) – Whether to check that the input matrices contain only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination) if the inputs do contain infinities or NaNs.

Returns:

x – The solution to the system A x = b

Return type:

array

See also

cho_factor: Cholesky factorization of a matrix

Examples

>>> import numpy as np
>>> from scipy.linalg import cho_factor, cho_solve
>>> A = np.array([[9, 3, 1, 5], [3, 7, 5, 1], [1, 5, 9, 2], [5, 1, 2, 6]])
>>> c, low = cho_factor(A)
>>> x = cho_solve((c, low), [1, 1, 1, 1])
>>> np.allclose(A @ x - [1, 1, 1, 1], np.zeros(4))
True

gpmp.num.cholesky_solve(A, b)[source]¶

gpmp.num.cholesky_inv(A)[source]¶

gpmp.num.scaled_distance(loginvrho: Any, x: Any, y: Any) → Any[source]¶

gpmp.num.scaled_distance_elementwise(loginvrho: Any, x: Any, y: Any | None) → Any[source]¶

Random numbers¶

Use gpmp.num.set_seed for backend-specific reproducibility. Random functions return backend-native objects.

gpmp.num.set_seed(seed: int) → None[source]¶: Set the global NumPy generator seed.

gpmp.num.rand(*shape: int) → Any[source]¶

gpmp.num.randn(*shape: int) → Any[source]¶

gpmp.num.choice(a: Any, size: int | None = None, replace: bool = True, p: Any | None = None) → Any[source]¶

gpmp.num.permutation(x: Any) → Any[source]¶

Differentiation helpers¶

The torch backend provides automatic differentiation. The NumPy backend uses finite-difference helpers where available. Higher-level parameter-selection code relies on DifferentiableSelectionCriterion to interface with SciPy optimizers.

gpmp.num.grad(f: Callable[[Any], Any]) → Callable[[Any], Any][source]¶

Return function that computes gradient of scalar f via finite differences.

Uses 5-point central difference formula for accuracy. Suitable for low to moderate dimensional problems.

Parameters:: f (callable) – Scalar-valued function taking an array and returning a scalar.
Returns:: Function grad_f(x) that computes nabla f(x) using finite differences.
Return type:: callable

gpmp.num.value_and_grad(f: Callable[[Any], Any], x: Any, *, h: float64 = 1e-05) → Tuple[Any, Any][source]¶: Returns (y, grad_y) where y = f(x) is scalar. Uses derivative_finite_diff on each coordinate (expects scalar input).

class gpmp.num.DifferentiableSelectionCriterion(crit: Callable[[Any, Any, Any], Any], x: Any, z: Any)[source]¶

class gpmp.num.BatchDifferentiableSelectionCriterion(crit: Callable[[Any, Any, Any], Any], loader: Iterable[Tuple[Any, Any]], reduction: str = 'mean', batches_per_eval: int = 0)[source]¶