Evaluating the squared-exponential covariance function in Gaussian processes with integral observations

1 Introduction

The Gaussian process (GP) [18]

has over the past decade become an important and fairly standard tool for system identification. The flexibility of the GP and more specifically its inherent data-driven bias-variance trade off makes it a very successful model for linear systems

[16, 3, 17], where it in fact improves performance compared to more classical models. When it comes to nonlinear system the GP has also shown great promise, for example when combined with the state space model [6, 5, 20] and when used in a classic auto-regressive setting [13, 2].

A useful property of GP is that it is closed under linear operators, such as integrals [18, 15, 9]. This enables GP regression to be applied to data that is related to the underlying object of interest via a line integral. In this case, evaluating the resulting measurement covariance involves double line integrals of the covariance function. There are no analytical solutions available and we are forced to numerical approximations. To give a few examples of where this problem arises, we mention; optimization algorithms [11, 21], continuous-time system identification, quadrature rules [14], and tomographic reconstruction [12].

The contribution of this paper is a numerically efficient and robust method for solving these integrals when the squared-exponential covariance function is used. We also provide a thorough comparison of the accuracy and computational time to existing state of the art methods for this rather specific problem.

2 Problem Statement

Consider a set of $n$ line integral observations of the form

y_{i} = | | w_{i} | | \int_{0}^{1} f (w_{i} s + p_{i}) d s + e_{i},

(1)

that provide measurements of a scalar function $f$ over the input space $z \in R^{m}$ , where $e_{i} \sim N (0, σ_{n}^{2})$ ¹¹1The notation denotes that $e_{i}$

is normally distributed with mean

0

and variance

σ_{n}^{2}

.. By

p_{i}

we denote the start point of the line while the vector

w_{i}

specifies its direction and length; hence the end point is given by

w_{i} + p_{i}

, and any point on the line can be expressed as

w_{i} s + p_{i}

with

s \in [0, 1]

. Given this set we wish to make predictions

f_{*}

of the function for a new input

z_{*}

In GP regression this is accomplished through the specification of a prior covariance function $k (z, z^{'})$ which defines the correlation between any two function values and loosely speaking controls how smooth we believe the underlying function to be. Although there exists many possible covariance functions [18], here we will consider only the squared-exponential covariance function;

k (z, z^{'}) = exp (- 0.5 (z - z^{'})^{T} V (z - z^{'})),

(2)

where $V$ is a positive definite and symmetric scaling matrix.

The covariance between a measurement $y_{i}$ and the function at a test input $z_{*}$ is given by the single integral

K_{i *} = ∥ w_{i} ∥ 1 \int 0 exp (- 0.5 (w_{i} t + v_{i})^{T} V (w_{i} t + v_{i})) d t,

(3)

where $v_{i} = p_{i} - z_{*}$ . Two cases need considering.

Case 1: $∥ w_{i} ∥ = 0$ . Then the integral (3) reduces to

K_{i *} = ∥ w_{i} ∥ exp (- \frac{1}{2} v_{i}^{T} V v_{i}) = 0.

(4)

Case 2: $∥ w_{i} ∥ > 0$ . Then the solution is given by

	$\begin{matrix} K_{i *} & = ∥ w_{i} ∥ \sqrt{\frac{π}{2 c_{1}}} exp (\frac{c_{2}^{2}}{8 c_{1}} - \frac{1}{2} c_{3}) (erf (\frac{2 c_{1} + c_{2}}{2 \sqrt{2 c_{1}}}) - erf (\frac{c_{2}}{2 \sqrt{2 c_{1}}})), \end{matrix}$	(5a)
where $erf (\cdot)$ is the error function [1] and
	$\begin{matrix} c_{1} & = w_{i}^{T} V w_{i}, c_{2} = 2 w_{i}^{T} V v_{i}, c_{3} = v_{i}^{T} V v_{i} . \end{matrix}$	(5b)

Determining the covariance function between a pair of line integral measurements requires evaluating the following double integral

(6)

where $u_{i j} = p_{i} - p_{j}$ . This is analytically intractable and hence we are forced to numerical approximations.

3 Related Work

Here, we provide a brief overview of existing approaches and their benefits and disadvantages.

The need to perform a double integral can be avoided entirely by the use of reduced rank approximation schemes that approximate the squared-exponential by a finite set of basis functions. One such scheme is the Hilbert space approximation [19] that was used in [12] to perform inference of strain fields from line integral measurements. Although the problem of calculating the double line integral is removed, the number of basis functions required and the subsequent computational cost grows exponentially with the problem dimension; making this method infeasible for higher-dimensional problems. An additional drawback is that the expressiveness of the model is limited to the basis functions chosen which may make this approach undesirable in many applications. Although this approach may be appropriate for some applications, the rest of this paper focuses on solutions that do not approximate the covariance function.

A simple approach to evaluate the double line integral is the use of 2D numerical integration methods (e.g. using a 2D Simpson’s rule). Such an approach has a trade off between computation time and accuracy. Although the approach detailed later still requires the numerical evaluation of a single integral, we have found it to be more efficient for the accuracy provided.

Another approach is to transform the problem into the double integral of a bivariate normal distribution. The integral in (6) can be rewritten as


	$I_{2} = 1 \int 0 1 \int 0 e^{- \frac{1}{2} f (t, s)} d t d s,$		(7a)
where the function $f (t, s)$ is given by

	$f (t, s)$	$= (p_{i} - p_{j} + w_{i} t - w_{j} s)^{T} V$
		$(p_{i} - p_{j} + w_{i} t - w_{j} s)$	(7b)
		$= {(\begin{matrix} t s \end{matrix})}^{T} (\begin{matrix} a & b b & c \end{matrix}) (\begin{matrix} t s \end{matrix}) - 2 (\begin{matrix} d & e \end{matrix}) (\begin{matrix} t s \end{matrix}) + ¯ f,$	(7c)
with the following variables

	$Σ^{- 1} ≜ (\begin{matrix} a & b b & c \end{matrix})$		(7d)
	$d$	$= (p_{i} - w_{j})^{T} V w_{i},$	(7e)
	$e$	$= - (p_{i} - p_{j})^{T} V w_{j},$	(7f)
	$¯ f$	$= (p_{i} - p_{j})^{T} V (p_{i} - p_{j}) .$	(7g)

When the matrix $Σ^{- 1}$ is invertible, then it is possible to reformulate the integral and employ integration methods for bivariate normal distributions (see e.g. [7] and the references therein).

However, when $0 < | a c - b^{2} | < ϵ$ for some small value $ϵ > 0$ , then this approach fails to reliably deliver accurate solutions. This problem is most notable when $a c - b^{2} = 0$ , in which case $Σ$ is not invertible and the bivariate approach cannot be used. There is an alternate approach in this particular case, but the user is now faced with a choice to numerically determine the appropriate threshold to decide between the solution methods. This is most notably an issue for lines that are close to co-linear. Importantly, the case where $a c - b^{2} = 0$ occurs frequently, e.g. along the diagonal entries of the covariance matrix.

Transformation of the problem into the Bivariate Normal is presented in the Appendix of [11]. However the non-invertible case is not explicitly dealt with, therefore we have provided the details for applying this approach and the necessary cases in Appendix A.

4 Efficient Computation of the Double Line Integral

Here, a method to accurately evaluate the double line integral (6) in a computationally efficient manner is presented. The integral can be rewritten as

	$\begin{matrix} K_{i j} & = ∥ w_{i} ∥ ∥ ∥ w_{j} ∥ ∥ 1 \int 0 1 \int 0 exp (- \frac{1}{2} (a - b s + c t - d s t + e t^{2} + f s^{2})) d t d s, \end{matrix}$	(8a)
where
	$\begin{matrix} a = u_{i j}^{T} V u_{i j}, & b = 2 u_{i j}^{T} V w_{j}, & c = 2 u_{i j}^{T} V w_{i}, d = 2 w_{j}^{T} V w_{i}, & e = w_{i}^{T} V w_{i}, & f = w_{j}^{T} V w_{j} . \end{matrix}$	(8b)

The error function can be used to provide a solution to the integral over $t$ ;

\begin{matrix} K_{i j} & = ∥ w_{i} ∥ ∥ ∥ w_{j} ∥ ∥ \sqrt{\frac{π}{2 e}} \int_{0}^{1} Γ_{1} (s) Γ_{2} (s) d s, \end{matrix}

(9)

where

\begin{matrix} Γ_{1} (s) & = exp (\frac{1}{2} (b s - a - f s^{2} + \frac{(c - d s)^{2}}{4 e})) Γ_{2} (s) & = erf (\frac{c - d s + 2 e}{2 \sqrt{} 2 e}) - erf (\frac{c - d s}{2 \sqrt{2 e}}) . \end{matrix}

(10)

For numerical reasons the exponent should not be split into a product of $exp (- a)$ and $exp (\frac{1}{2} (b s - f s^{2} + \frac{(c - d s)^{2}}{4 e}))$ . The reason is that when $∥ ∥ V u_{i j} ∥ ∥$ is large, $exp (- a) \to 0$ while $exp (\frac{1}{2} (b s - f s^{2} + \frac{(c - d s)^{2}}{4 e})) \to \infty$ and rounding errors become a problem.

Numerical methods can then be used to evaluate the remaining integral. Care should be taken when considering the following two cases:

Case 1: Either $∥ w_{i} ∥ = 0$ or $∥ ∥ w_{j} ∥ ∥ = 0$ . Then the problem reduces to a single line integral and two sub cases should be considered in implementation.

Case 1a: $∥ w_{i} ∥ = 0$ and $∥ ∥ w_{j} ∥ ∥ > 0$ . Then Equation (9) is numerically unstable and two approaches can be taken. Either the solution to the single integral (5) can be used, which requires the user to determine a numerical threshold to decide between the two cases. Alternatively, the order of the integrals can be swapped, and Equation (9) can be applied.
Case 1b: $∥ w_{i} ∥ > 0$ and $∥ ∥ w_{j} ∥ ∥ = 0$ . Equation (9) is numerically stable, and hence the user can either ignore this case or define a numerical tolerance below which to use the solution to the single integral (5).

Case 2: $∥ w_{i} ∥ = 0$ and $∥ ∥ w_{j} ∥ ∥ = 0$ . Then Equation (9) is numerically unstable and the solution is instead given by

K_{i j} = ∥ w_{i} ∥ ∥ ∥ w_{j} ∥ ∥ exp (- \frac{1}{2} a) = 0.

(11)

Psudocode is provided in Algorithm 1 and a Matlab mex function implemented in C is provided at [10]. In the implementation and psuedocode, the order of the lines is swapped whenever $∥ ∥ {V w}_{j} ∥ ∥ > ∥ {V w}_{i} ∥$ as this means the $erf ()$ function is used to evaluate the larger interval and the numerical integration the shorter interval; which was found to require less function evaluations. The C implementation utilises the BLAS functions dgemm, ddot, and dnrm2 [4], and the GNU Scientific Library’s non-adaptive Gauss Konrod function gsl_integration_qng [8] with relative and absolute error tolerance set to square root of double precision epsilon.

1:procedure Double Integral(

u_{i j}

w_{i}

w_{j}

V

ϵ

)

2: if

| | V w_{i} | | < ϵ

| | V w_{2} | | < ϵ

then

K_{i j} := ∥ w_{i} ∥ ∥ ∥ w_{j} ∥ ∥ exp (- \frac{1}{2} u_{i j}^{T} V u_{i j})

4: else

5: if

∥ {V w}_{i} ∥ < ∥ ∥ {V w}_{j} ∥ ∥

then

a := u_{i j}^{T} V u_{i j}

b := - 2 u_{i j}^{T} V w_{i}

c := - 2 u_{i j}^{T} V w_{j}

d := 2 w_{i}^{T} V w_{j}

e := w_{j}^{T} V w_{j}

f := w_{i}^{T} V w_{i}

9: else

10:

a := u_{i j}^{T} V u_{i j}

b := 2 u_{i j}^{T} V w_{j}

11:

c := 2 u_{i j}^{T} V w_{i}

d := 2 w_{j}^{T} V w_{i}

12:

e := w_{i}^{T} V w_{i}

f := w_{j}^{T} V w_{j}

13: end if

14: Evaluate

K_{i j}

by applying a 1D numerical integration technique to equation (9)

15: end if

16: return

K_{i j}

17:end procedure

Algorithm 1 Psuedocode for evaluation of Equation 6

5 Results and Analysis

The performance of the methods was evaluated for several sets of double line integrals; containing 10,000 pairs each and with input dimension $m = 6$ . The sets were chosen to test a variety of cases and in particular to evaluate the performance on the transition between the corner cases. The following sets were used:

Set 1: A standard set; $V = I$ , $w_{i, k} \sim U (0, 1)$ , $w_{j, k} \sim U (0, 1)$ , and $u_{i j, k} \sim U (0, 1)$ for $k = 1, \dots, m$ ,
Set 2: Almost colinear lines; $V = I$ , $w_{i, k} \sim U (0, 1)$ , $w_{j, k} \sim w + U (0, 1 \times 10^{- 8})$ , and $u_{i j, k} \sim U (0, 1)$ for $k = 1, \dots, m$ ,
Set 3: Randomly selected diagonal scaling matrix; $V_{k k} \sim U (0, 1)$ , $w_{i, k} \sim U (0, 1)$ , $w_{j, k} \sim U (0, 1)$ , and $u_{i j, k} \sim U (0, 1)$ for $k = 1, \dots, m$ ,
Set 4: Reduced to single integral; $V = I$ , $w_{i} = 0$ , $w_{j, k} \sim U (0, 1)$ , and $u_{i j, k} \sim U (0, 1)$ for $k = 1, \dots, m$ ,
Set 5: Nicely scaled (large integral intervals and small random scaling); $V_{k k} \sim U (0, 0.01)$ , $w_{i, k} \sim U (0, 10)$ , $w_{j, k} \sim U (0, 10)$ , and $u_{i j, k} \sim U (0, 10)$ for $k = 1, \dots, m$ ,
Set 6: Poorly scaled (large integral intervals and large random scaling); $V_{k k} \sim U (0, 10)$ , $w_{i, k} \sim U (0, 10)$ , $w_{j, k} \sim U (0, 10)$ , and $u_{i j, k} \sim U (0, 10)$ for $k = 1, \dots, m$ ,
Set 7:: Almost reduced to single integral; $V = I$ , $w_{i, k} \sim U (0, 1 \times 10^{- 8})$ , $w_{j, k} \sim U (0, 1)$ , and $u_{i j, k} \sim U (0, 1)$ for $k = 1, \dots, m$ ,
Set 8: Almost reduced to no integral; $V = I$ , $w_{i, k} \sim U (0, 1 \times 10^{- 8})$ , $w_{j, k} \sim U (0, 1 \times 10^{- 8})$ , and $u_{i j, k} \sim U (0, 1)$ for $k = 1, \dots, m$ .

Here, the notation $v_{i, k} \sim U (a, b)$ denotes that the $k^{th}$ element of the vector $v_{i}$ is distributed uniformly between $a$ and $b$ .

Errors were calculated by comparison to Matlab’s Integral2 using absolute and relative error tolerance of double precision epsilon. The mean magnitude of the error for each set are recorded in Table 1. The computation times required by each method were consistent across the cases and so the average computation times are recorded in Table 2.

Method

Set 1

Set 2

Set 3

Set 4

Set 5

Set 6

Set 7

Set 8

Simpson’s Rule

p = 10

3.66 \times 10^{- 6}

3.49 \times 10^{- 6}

1.62 \times 10^{- 6}

0

1.63 \times 10^{- 4}

8.93 \times 10^{- 5}

3.14 \times 10^{- 14}

2.27 \times 10^{- 32}

p = 100

3.62 \times 10^{- 10}

3.44 \times 10^{- 10}

1.61 \times 10^{- 10}

0

1.61 \times 10^{- 8}

7.08 \times 10^{- 10}

3.12 \times 10^{- 18}

9.07 \times 10^{- 32}

p = 200

2.23 \times 10^{- 11}

2.15 \times 10^{- 11}

1.00 \times 12^{- 11}

0

1.01 \times 10^{- 9}

4.39 \times 10^{- 11}

1.95 \times 10^{- 19}

1.27 \times 10^{- 31}

Bivariate Normal

2.06 \times 10^{- 15}

1.84 \times 10^{- 1}

5.03 \times 10^{- 15}

0

4.81 \times 10^{- 13}

9.30 \times 10^{- 17}

3.26 \times 10^{- 17}

1.58 \times 10^{- 25}

Proposed

1.80 \times 10^{- 15}

6.39 \times 10^{- 16}

4.35 \times 10^{- 15}

0

4.15 \times 10^{- 13}

8.32 \times 10^{- 14}

9.56 \times 10^{- 24}

2.10 \times 10^{- 25}

Table 1: Mean magnitude of the error for each set. Errors are computed with respect to Matlab’s Integral2 function.

Method	Time ( $\times 10^{- 5}$ s)
Matlab Integral2	149.38*
Simpson’s Rule
$p = 10$	$1.24$
$p = 100$	$14.54$
$p = 200$	$54.86$
Bivariate Normal	$6.27$
Proposed	$1.21$

Table 2: Average computation time across all sets for each method. Note* the reported time is for Integral2 with absolute and relative tolerance of double precision epsilon, using the default tolerances reduces the computation time to

0.0011687 $ s $

A 2D Simpson’s rule was implemented with $p$ sub intervals along $t$ and $s$ ; results are shown for $p = 10$ , $p = 100$ and $p = 200$ . This 2D Simpson’s rule gives the lowest error for Set 8 where the integrand is close to a constant function and is quite fast a choice of small $p$ . Sets 1 through 6 it show relatively large errors with the errors decreasing as $p$ is increased. However, even with $p = 200$ the errors are still larger than for the Bivariate Normal method and the proposed method and the computation times are an order of magnitude greater.

The Bivariate Normal method has acceptable small errors across all sets excluding set 2. The Bivariate normal method required a different algorithm for the co-linear and non co-linear cases and Set 2 lies on the transition between these cases. Figure 1 provides a histogram of the error magnitudes for Set 2 as computed by each method. The histogram shows that the Bivariate normal method can result in large errors which would be prohibitive for many applications.

Our proposed method has acceptable errors across all sets. The worst errors, of $1 \times 10^{- 13}$ and $1 \times 10^{- 14}$ for Set 5 and Set 6 respectively, could potentially be improved by using an adaptive integration method that may be more suited to the scaling in these sets. The average computation time is the smallest reported.

Figure 1: Histogram of the error magnitudes for evaluating the double line integrals in set 2 (the almost co-linear case); (a) Simpson’s rule with $p = 10$ , (b) Simpson’s rule with $p = 200$ , (c) Bivariate Normal method, (d) our proposed method.

6 Conclusion

In this paper we have focussed on the specific problem of evaluating a double line integral over the squared-exponential covariance function. This problem arises in several applications of Gaussian process regression. Existing approaches to this problem and their advantages and disadvantages were outlined. An alternative method was proposed and its performance compared to the existing methods was evaluated for a number of sets of line integrals corresponding to both common cases and corner cases. The comparison shows the proposed method is computationally efficient while providing acceptable accuracy across all tested cases.

Acknowledgements

This research was financially supported by the Swedish Foundation for Strategic Research (SSF) via the project ASSEMBLE (contract number: RIT15-0012) and by the Swedish Research Council via the projects Learning flexible models for nonlinear dynamics (contract number: 2017-03807) and NewLEADS - New Directions in Learning Dynamical Systems (contract number: 621-2016-06079).

Appendix A Bivariate normal method

In this appendix we show how to compute the double integral in (6) by reformulating the problem as the double integral of a Bivariate normal distribution. The integral can be written as


	$I_{2} = \int_{0}^{1} \int_{0}^{1} e^{- \frac{1}{2} f (t, s)} d t d s,$		(12a)
where the function $f (s, t)$ is given by

	$f (t, s)$	$= (p_{i} - p_{j} + w_{i} t - w_{j} s)^{T} V (p_{i} - p_{j} + w_{i} t - w_{j} s)$
		$= {(\begin{matrix} t s \end{matrix})}^{T} (\begin{matrix} a & b b & c \end{matrix}) (\begin{matrix} t s \end{matrix}) - 2 (\begin{matrix} d & e \end{matrix}) (\begin{matrix} t s \end{matrix}) + ¯ f,$	(12b)
with the following variables

	$(\begin{matrix} a & b b & c \end{matrix})$	$= (\begin{matrix} w_{i}^{T} V w_{i} & - w_{i}^{T} V w_{j} - w_{j}^{T} V w_{i} & w_{j}^{T} V w_{j} \end{matrix}),$	(12c)
	$d$	$= - (p_{i} - p_{j})^{T} V w_{i},$	(12d)
	$e$	$= (p_{i} - p_{j})^{T} V w_{j},$	(12e)
	$¯ f$	$= (p_{i} - p_{j})^{T} V (p_{i} - p_{j}) .$	(12f)

There are two cases that need to be considered separately as the solutions to each case do not commute. The reason for needing two cases is that when $w_{i}$ is colinear with $w_{j}$ (as occurs along the block diagonal elements of the $K_{ℓ_{k}, ℓ_{k}}$ matrix), then the above matrix involving $a, b, c$ is not invertible.

a.1 Double integral with $a c - b^{2} = 0$

This case can occur in four distinct ways.

Case 1: $∥ w_{i} ∥ = 0$ and $∥ w_{j} ∥ > 0$ . Then $a = b = d = 0$ and $f (t, s) = ¯ f - 2 e s + c s^{2}$ and we can use (4) and (5) to solve this integral.

Case 2: $∥ w_{i} ∥ > 0$ and $∥ w_{j} ∥ = 0$ . Then $c = b = e = 0$ and $f (t, s) = ¯ f + 2 d t + a t^{2}$ and we can again use (4) and (5) to solve this integral.

Case 3: $∥ w_{j} ∥ = 0$ and $∥ w_{i} ∥ = 0$ . Then $a = b = c = d = e = 0$ and $f (t, s) = ¯ f$ and the solution is $I_{2} = e^{- \frac{¯ f}{2}} .$

Case 4: $∥ w_{i} ∥ > 0$ and $∥ w_{j} ∥ > 0$ . Since we are already in the case where $a c - b^{2} = 0$ , and $a > 0$ and $c > 0$ then it follows that $w_{i} = β w_{j}$ for some $β \neq 0$ since

	$a c - b^{2}$	$= (w_{i}^{T} V w_{i}) (w_{j}^{T} V w_{j}) - (w_{i}^{T} V w_{j})^{2},$
		$= (w_{i}^{T} V w_{i}) (β^{2} w_{i}^{T} V w_{i}) - (β w_{i}^{T} V w_{i})^{2},$
		$= β^{2} (w_{i}^{T} V w_{i})^{2} - β^{2} (w_{i}^{T} V w_{i})^{2}$
		$= 0$

as claimed. Therefore, by defining $β$ as

β ≜ \frac{w_{i}^{T} V w_{j}}{w_{i}^{T} V w_{i}}

(13)

then $f (t, s)$ can be expressed as

f (t, s)

= f + 2 t (p_{i} - p_{j})^{T} V w_{i} - 2 β s (p_{i} - p_{j})^{T} V w_{i} + t^{2} (w_{i}^{T} V w_{i}) - 2 β t s (w_{i}^{T} V w_{i}) + β^{2} s^{2} (w_{i}^{T} V w_{i}) .

If we make a change of variables for $s$ to

¯ s = β s ⟹ f (t, ¯ s) = f + a (t^{2} - 2 t ¯ s + s^{2})

and the integral becomes

	$I_{2}$	$= \frac{e^{- f / 2}}{β} \int_{0}^{β} \int_{0}^{1} e^{- \frac{a}{2} (t^{2} - 2 t ¯ s + {¯ s}^{2})} d t d ¯ s$
		$= \frac{e^{- f / 2}}{β} \int_{0}^{β} \int_{0}^{1} e^{- \frac{a}{2} (t - ¯ s)^{2}} d t d ¯ s$

By a change of variables $ϕ = \sqrt{\frac{a}{2}} (t - ¯ s)$ then

	$I_{2}$	$= \frac{\sqrt{2} e^{- f / 2}}{β \sqrt{a}} \int_{0}^{β} \int_{- ¯ s \sqrt{a / 2}}^{(1 - ¯ s) \sqrt{a / 2}} e^{- ϕ^{2}} d ϕ d% ¯ s$		(14)
		$= \frac{\sqrt{π} \sqrt{2} e^{- f / 2}}{2 β \sqrt{a}} \int_{0}^{β} (erf ((1 - ¯ s) \sqrt{a / 2}) - erf (- ¯ s \sqrt{a / 2})) d ¯ s$		(15)

where the second equality stems from the definition of the error function

erf (x) = \frac{2}{\sqrt{π}} \int_{0}^{x} e^{- u^{2}} d% u .

(16)

Next we can exploit the fact that the integral of the error function is given by (using integration by parts)

\int_{a}^{b} erf (x) d x = b erf (b) + \frac{e^{- b^{2}}}{\sqrt{π}} - a erf (a) - \frac{e^{- a^{2}}}{\sqrt{π}} .

(17)

To apply this to $I_{2}$ it is helpful to split the integral component via

$I_{2}$	$= \frac{\sqrt{π} \sqrt{2} e^{- f / 2}}{2 β \sqrt{a}} (I_{3} - I_{4}),$	(18)
$I_{3}$	$≜ \int_{0}^{β} erf ((1 - ¯ s) \sqrt{a / 2}) d ¯ s$	(19)
$I_{4}$	$≜ \int_{0}^{β} erf (- ¯ s \sqrt{a / 2}) d% ¯ s$	(20)

By change of variables $ζ = (1 - ¯ s) \sqrt{a / 2}$ then $I_{3}$ becomes

	$I_{3}$	$= \sqrt{\frac{2}{a}} \int_{(1 - β) \sqrt{a / 2}}^{\sqrt{a / 2}} % erf (ζ) d ζ$
		$= \sqrt{\frac{2}{a}} [\sqrt{\frac{a}{2}} erf (\sqrt{\frac{a}{2}}) + \frac{e^{- \frac{a^{2}}{2}}}{\sqrt{π}} - (1 - β) \sqrt{\frac{a}{2}} erf ((1 - β) \sqrt{\frac{a}{2}}) - \frac{e^{- (1 - β)^{2} \frac{a^{2}}{2}}}{\sqrt{π}}]$

And similarly for $I_{4}$ with change of variables $ψ = - ¯ s \sqrt{a / 2}$ to arrive at

	$I_{4}$	$= \sqrt{\frac{2}{a}} \int_{- β \sqrt{a / 2}}^{0} erf (ψ) d ψ$		(21)
		$= \sqrt{\frac{2}{a}} [\frac{1}{\sqrt{π}} + β \sqrt{\frac{a}{2}} erf (- β \sqrt{\frac{a}{2}}) - \frac{e^{- β^{2} \frac{}{a^{2}} 2}}{\sqrt{π}}]$		(22)

a.2 Double integral with $a c - b^{2} > 0$

By defining

	$z$	$= (\begin{matrix} t s \end{matrix}), Σ^{- 1} = (\begin{matrix} a & b b & c \end{matrix}), ν = (\begin{matrix} d e \end{matrix}), μ = Σ ν$
	$g (z)$	$≜ (z - μ)^{T} Σ^{- 1} (z - μ), h ≜ f - μ^{T} Σ^{- 1} μ$

we have that

f (t, s)

= g (z) + h

(23)

Therefore,

e^{- \frac{1}{2} f (t, s)}

= e^{- \frac{1}{2} g (z) - \frac{h}{2}} = 2 π e^{- \frac{h}{2}} \sqrt{det Σ} N (z | μ, Σ)

(24)

where $N (z | μ, Σ)$

denotes the probability density function for the normal distributed random variable

z

, with mean value

μ

and covariance

Σ

. With the change of variables

¯ z = z - μ

the integral can be expressed as

I_{2} = 2 π e^{- \frac{h}{2}} \sqrt{det Σ} \int_{- μ_{1}}^{1 - μ_{1}} \int_{- μ_{2}}^{1 - μ_{2}} N (¯ z | 0, Σ) d {¯ z}_{1} d {¯ z}_{2},

(25)

Note that the density in (25) is a bivariate normal distribution, which means that it can be written as

N (¯ z | 0, Σ)

= \frac{1}{2 π \sqrt{1 - ρ^{2}} \sqrt{c_{11} c_{22}}} exp (- \frac{α (¯ z)}{2 (1 - ρ^{2})})

where

		$[\begin{matrix} c_{11} & c_{12} c_{21} & c_{22} \end{matrix}] ≜ Σ = {[\begin{matrix} a & - b - b & c \end{matrix}]}^{- 1} = \frac{1}{a c - b^{2}} [\begin{matrix} c & b b & a \end{matrix}]$
	$α (¯ z)$	$= {(\frac{{¯ z}_{1}}{\sqrt{c_{11}}})}^{2} - 2 ρ (\frac{{¯ z}_{1}}{\sqrt{c_{11}}}) (\frac{{¯ z}_{2}}{\sqrt{c_{22}}}) + {(\frac{{¯ z}_{2}}{\sqrt{c_{22}}})}^{2}$
	$ρ$	$= \frac{c_{12}}{\sqrt{c_{11} c_{22}}} .$

By change of variables ${˜ z}_{1} = \frac{{¯ z}_{1}}{\sqrt{c_{11}}}$ and ${˜ z}_{2} = \frac{{¯ z}_{2}}{\sqrt{c_{22}}}$ the resulting integral as


	$I_{2} = 2 π e^{- \frac{h}{2}} \sqrt{det Σ} \int_{- \frac{μ_{1}}{\sqrt{c_{11}}}}^{\frac{1 - μ_{1}}{\sqrt{c_{11}}}} \int_{- \frac{μ_{2}}{\sqrt{} c_{22}}}^{\frac{1 - μ_{2}}{\sqrt{c_{22}}}} β (ρ) d {˜ z}_{1} d {˜ z}_{2} .$	(26a)
where

	$β (ρ) = \frac{1}{2 π \sqrt{1 - ρ^{2}}} exp (- \frac{{˜ z}_{1}^{2} - 2 ρ {˜ z}_{1} {˜ z}_{2} + {˜ s}_{2}^{2}}{2 (1 - ρ^{2})}) .$	(26b)

The integral (26a) has been widely studied in the literature, see e.g. [7] and the references therein. There is efficient code available²²2www.math.wsu.edu/faculty/genz/software/matlab/bvn.m to compute this integral.

References

[1] Larry C Andrews and Larry C Andrews. Special functions of mathematics for engineers. McGraw-Hill New York, 1992.
[2] Hildo Bijl, Thomas B. Schön, Jan-Willem van Wingerden, and Michel Verhaegen. System identification through online sparse Gaussian process regression with input noise. In IFAC Journal of Systems and Control, volume 2, pages 1–11, 2017.
[3] Tianshi Chen, Henrik Ohlsson, and Lennart Ljung.
On the estimation of transfer functions, regularizations and Gaussian processes–revisited.
Automatica, 48(8):1525–1535, 2012.
[4] Jack J Dongarra, Jeremy Du Croz, Sven Hammarling, and Iain S Duff. A set of level 3 basic linear algebra subprograms. ACM Transactions on Mathematical Software (TOMS), 16(1):1–17, 1990.
[5] Roger Frigola, Yutian Chen, and Carl E. Rasmussen. Variational Gaussian process state-space models. In Advances in Neural Information Processing Systems (NIPS), 2014.
[6] Roger Frigola, Fredrik Lindsten, Thomas B. Schön, and Carl E. Rasmussen. Bayesian inference and learning in Gaussian process state-space models with particle MCMC. In Advances in Neural Information Processing Systems (NIPS) 26, Lake Tahoe, NV, USA, December 2013.
[7] Alan Genz. Numerical computation of rectangular bivariate and trivariate normal and t probabilities. Statistics and Computing, 14(3):251–260, 2004.
[8] Brian Gough. GNU scientific library reference manual. Network Theory Ltd., 2009.
[9] Thore Graepel.
Solving noisy linear operator equations by Gaussian processes: Application to ordinary and partial differential equations.
In
International Conference on Machine Learning
, pages 234–241, 2003.
[10] Johannes Hendriks. Mex code for evaluating double line integrals of the scquared exponential covariance function. https://github.com/jnh277/lineIntSquaredExponential.git, 2018.
[11] Philipp Hennig and Martin Kiefel. Quasi-newton methods: A new direction. Journal of Machine Learning Research, 14, 06 2012.
[12] Carl Jidling, Johannes Hendriks, Niklas Wahlström, Alexander Gregg, Thomas B Schön, Chris Wensrich, and Adrian Wills. Probabilistic modelling and reconstruction of strain. Nuclear instruments and methods in physics research section B: Beam interactions with materials and atoms, 436:141–155, 2018.
[13] Juc Kocijan, Agathe Girard, Blaz Banko, and Roderick Murray-Smith. Dynamic systems identification with Gaussian processes. Mathematical and Computer Modelling of Dynamical Systems, 11(4):411—424, 2005.
[14] Thomas P Minka. Deriving quadrature rules from Gaussian processes. Technical report, Technical report, Statistics Department, Carnegie Mellon University, 2000.
[15] Athanasios Papoulis. Probability, random variables, and stochastic processes. 1965.
[16] Gianluigi Pillonetto and Giuseppe De Nicolao. A new kernel-based approach for linear system identification. Automatica, 46(1):81—93, 2010.
[17] Gianluigi Pillonetto, Francesco Dinuzzo, Tianshi Chen, Giuseppe De Nicolao, and Lennart Ljung. Kernel methods in system identification, machine learning and function estimation: a survey. Automatica, 50(3):657–682, 2014.
[18] Carl Edward Rasmussen and Christopher KI Williams. Gaussian process for machine learning. MIT press, 2006.
[19] Arno Solin and Simo Särkkä. Hilbert space methods for reduced-rank Gaussian process regression. arXiv preprint arXiv:1401.5508, 2014.
[20] Andreas Svensson and Thomas B. Schön. A flexible state space model for learning nonlinear dynamical systems. Automatica, 80:189–199, June 2017.
[21] Adrian G Wills and Thomas B Schön. On the construction of probabilistic newton-type algorithms. In Proceedings of the 56th IEEE Annual Conference on Decision and Control (CDC), pages 6499–6504, Melbourne, Australia, 2017.

Evaluating the squared-exponential covariance function in Gaussian processes with integral observations

Authors

Interpolation error of misspecified Gaussian process regression

Fast Kernel Approximations for Latent Force Models and Convolved Multiple-Output Gaussian processes

Numerical method of computing Hadamard finite-part integrals with a non-integral power singularity at the endpoint over a half infinite interval

On Lebesgue Integral Quadrature

A numerical method for Hadamard finite-part integrals with an integral power singularity at an endpoint

Computing the matrix fractional power based on the double exponential formula

Deep kernel learning for integral measurements

1 Introduction

2 Problem Statement

3 Related Work

4 Efficient Computation of the Double Line Integral

5 Results and Analysis

6 Conclusion

Acknowledgements

Appendix A Bivariate normal method

a.1 Double integral with $a c - b^{2} = 0$

a.2 Double integral with $a c - b^{2} > 0$

References

Evaluating the squared-exponential covariance function in Gaussian processes with integral observations

Authors

Interpolation error of misspecified Gaussian process regression

Fast Kernel Approximations for Latent Force Models and Convolved Multiple-Output Gaussian processes

Numerical method of computing Hadamard finite-part integrals with a non-integral power singularity at the endpoint over a half infinite interval

On Lebesgue Integral Quadrature

A numerical method for Hadamard finite-part integrals with an integral power singularity at an endpoint

Computing the matrix fractional power based on the double exponential formula

Deep kernel learning for integral measurements

This week in AI

1 Introduction

2 Problem Statement

3 Related Work

4 Efficient Computation of the Double Line Integral

5 Results and Analysis

6 Conclusion

Acknowledgements

Appendix A Bivariate normal method

a.1 Double integral with ac−b2=0

a.2 Double integral with ac−b2>0

References

a.1 Double integral with $a c - b^{2} = 0$

a.2 Double integral with $a c - b^{2} > 0$