Comment on "Probabilistic Integration: A Role in Statistical Computation?"

TL;DR

This paper discusses the assumptions and implementation challenges of automatic Bayesian cubature, a probabilistic integration method that determines when a simulation has achieved a desired error tolerance with high confidence.

Contribution

It provides critical commentary on the modeling assumptions and practical issues in designing automatic Bayesian cubature methods.

Findings

Highlights key modeling assumptions in probabilistic integration.

Identifies implementation challenges in automatic Bayesian cubature.

Discusses the criteria for stopping simulations based on confidence levels.

Abstract

Probabilistic integration provides a criterion for stopping a simulation when a specified error tolerance is satisfied with high confidence. We comment on some of the modeling assumptions and implementation issues involved in designing an automatic Bayesian cubature.

Figures

Equations

$$\bigl{\lvert}\hat{\Pi}[f]-\Pi[f]\rvert\leq\lVert f\rVert_{\mathcal{H}}\lVert\mu(\hat{\pi})-\mu(\pi)\rVert_{\mathcal{H}},$$

$$\bigl{\lvert}\hat{\Pi}[f]-\Pi[f]\bigr{\rvert}\leq\varepsilon,$$

$$\displaystyle\operatorname{\mathbb{P}}\Bigl{[}\bigl{\lvert}\hat{\Pi}[f]-\Pi[g]\bigr{\rvert}\leq 2.58\sqrt{\operatorname{\mathbb{V}}_{n}[\Pi[g]]}\Bigr{]}=99\%,$$

$$V_n[\Pi [g]]=\Pi \Pi [c(\cdot ,\cdot )]-\Pi [\mathbf{c}(\cdot ,X)]{\mathbf{C}}^{-1}\Pi [\mathbf{c}(X,\cdot )].$$

$$c_{0}(\boldsymbol{x},\boldsymbol{x}^{\prime};\boldsymbol{\theta})=\prod_{l=1}^{d}\biggl{[}1-(-1)^{r}\gamma B_{2r}(|{x_{l}-x^{\prime}_{l}}|)\biggr{]},\\ \forall\boldsymbol{x},\boldsymbol{x}^{\prime}\in[0,1]^{d},\ \boldsymbol{\theta}=(r,\gamma),\ r\in\mathbb{N},\ \gamma>0,$$

$$\hat{\Pi }[f]$$

$$=\left(\frac{(1-1^T{\mathbf{C}}_0^{-1}\Pi [{\mathbf{c}}_0(X,\cdot )]1^T}{1^T{\mathbf{C}}_0^{-1}1}+\Pi [{\mathbf{c}}_0(\cdot ,X)]\right){\mathbf{C}}_0^{-1}\mathbf{f},$$

$$V_n[\Pi [g]]$$

$$\times (\frac{(1-\Pi [{\mathbf{c}}_0(\cdot ,X)]{\mathbf{C}}_0^{-1}1{)}^2}{1^T{\mathbf{C}}_0^{-1}1}$$

$$\displaystyle\qquad\qquad+\Pi\Pi[c_{0}(\cdot,\cdot)]-\Pi[\boldsymbol{c}_{0}(\cdot,X)]\boldsymbol{C}_{0}^{-1}\Pi[\boldsymbol{c}_{0}(X,\cdot)]\biggr{)},$$

$$\displaystyle\operatorname{\mathbb{P}}\Bigl{[}\bigl{\lvert}\hat{\Pi}[f]-\Pi[g]\bigr{\rvert}\leq t_{n-1,0.995}\sqrt{\operatorname{\mathbb{V}}_{n}[\Pi[g]]}\Bigr{]}=99\%.$$

$$\boldsymbol{\theta}_{\textup{EB}}=\operatorname*{argmin}_{\boldsymbol{\theta}}\biggl{\{}\log\left(\boldsymbol{f}^{T}\left[\boldsymbol{C}_{0}^{-1}-\frac{\boldsymbol{C}_{0}^{-1}\boldsymbol{1}\boldsymbol{1}^{T}\boldsymbol{C}_{0}^{-1}}{\boldsymbol{1}^{T}\boldsymbol{C}_{0}^{-1}\boldsymbol{1}}\right]\boldsymbol{f}\right)+\frac{1}{n}\log(\det(\boldsymbol{C}_{0}))\biggr{\}}.$$

$${\mathbf{C}}_0$$

$$\mathbf{V}$$

$$\mathbf{V}\text{ may be identified analytically},$$

$${\mathbf{v}}_1={\mathbf{V}}_1=1,$$

$$\tilde{\mathbf{b}}:={\mathbf{V}}^H\mathbf{b}\text{ requires only O(nlog(n)) operations }\forall \mathbf{b},$$

$$\Pi [c_0(\cdot ,\mathbf{x})]=1\forall \mathbf{x}.$$

$$\begin{gathered} \bullet \text{ Shift-invariant covariance functions, c0, which satisfy} \\ {c}_0(\mathbf{x},{\mathbf{x}}^{\prime })={\mathring{c}}_o(\mathbf{x}-{\mathbf{x}}^{\prime }\mathrm{mod}1)\forall \mathbf{x},{\mathbf{x}}^{\prime }\in [0,1{)}^d, \end{gathered}$$

$$\begin{gathered} \bullet \text{ Shifted rank-1 integration lattice node sets, {xi}i=1n, which satisfy} \\ \mathbf{x},{\mathbf{x}}^{\prime },{\mathbf{x}}^{\prime \prime }\in \{{\mathbf{x}}_i{\}}_{i=1}^n⟹\mathbf{x}+{\mathbf{x}}^{\prime }-{\mathbf{x}}^{\prime \prime }\mathrm{mod}1\in \{{\mathbf{x}}_i{\}}_{i=1}^n. \end{gathered}$$

$$\tilde{\mathbf{f}}$$

$$\hat{\Pi }[f]$$

$$V_n[\Pi [g]]$$

$${\mathbf{\theta }}_{\text{EB}}$$

$$k(\mathbf{x},{\mathbf{x}}^{\prime })=\sum _{i=1}^{\infty }{\varphi }_i(\mathbf{x}){\varphi }_i({\mathbf{x}}^{\prime }).$$

$$E({\Vert g\Vert }_{\mathcal{H}}^2)=E({\Vert \mathbf{A}\Vert }_2^2)={\Vert (1,1,\dots )\Vert }_2^2=\infty .$$

$$\hat{\Pi }[f]=\Pi [\hat{f}]=\Pi [\mathbf{k}(\cdot ,X)]{\mathbf{K}}^{-1}\mathbf{f},$$

$$\displaystyle\bigl{\lvert}\hat{\Pi}[f]-\Pi[f]\rvert^{2}\leq\lVert f-\hat{f}\rVert_{\mathcal{H}}^{2}\ \bigl{\{}\Pi\Pi[k(\cdot,\cdot)]-\Pi[\boldsymbol{k}(\cdot,X)]\boldsymbol{K}^{-1}\Pi[\boldsymbol{k}(X,\cdot)]\bigr{\}},$$

$$\Vert \hat{f}{\Vert }_{\mathcal{H}}^2={\mathbf{f}}^T{\mathbf{K}}^{-1}\mathbf{f}.$$

$$\hat{\Pi }[f]=E[\Pi [g]]=\Pi [{\mathbf{c}}_0(\cdot ,X)]{\mathbf{C}}_0^{-1}\mathbf{f},$$

$$\displaystyle\operatorname{\mathbb{V}}[\Pi[g]]=\frac{\boldsymbol{f}^{T}\boldsymbol{C}_{0}^{-1}\boldsymbol{f}}{n}\ \bigl{\{}\Pi\Pi[c_{0}(\cdot,\cdot)]-\Pi[\boldsymbol{c}_{0}(\cdot,X)]\boldsymbol{C}_{0}^{-1}\Pi[\boldsymbol{c}_{0}(X,\cdot)]\bigr{\}},$$

$$\displaystyle\operatorname{\mathbb{P}}\bigl{[}\bigl{\lvert}\hat{\Pi}[f]-\Pi[f]\rvert^{2}\leq t_{0.995,n}^{2}\operatorname{\mathbb{V}}[\Pi[g]]\bigr{]}=99\%.$$

$$\Vert f-\hat{f}{\Vert }_{\mathcal{H}}^2\le \frac{t_{0.995,n}^2\Vert \hat{f}{\Vert }_{\mathcal{H}}^2}{n}=\frac{t_{0.995,n}^2{\mathbf{f}}^T{\mathbf{K}}^{-1}\mathbf{f}}{n}.$$

Full text

Bayesian,

fast algorithms,

quasi-Monte Carlo,

T1This work is supported in part by NSF-DMS-1522687.

and

1 When to Stop?

In highlighting the possibilities of probabilistic integration, the authors of Briol et al. (2018+), henceforth abbreviated as PI, have suggested a useful stopping criterion for cubature. Numerical analysis provides an upper bound on the cubature error expressed as a product of the roughness of the integrand and the quality of our sampling scheme. For example, PI (5) quotes the error bound

[TABLE]

where

• •

the integrand, $f$, lies in a Hilbert space, $\mathcal{H}$,

• •

$\Pi[f]$ denotes the desired integral of $f$ defined in terms of the probability measure $\pi$, and

• •

$\hat{\Pi}[f]$ denotes a cubature defined in terms of the discrete measure $\hat{\pi}$.

The discrepancy between $\pi$ and $\hat{\pi}$ is defined as $\lVert\mu(\hat{\pi})-\mu(\pi)\rVert_{\mathcal{H}}$. As the sample size, $n$, increases, a well chosen sequence of discrete measures causes the discrepancy to tend to zero.

But, even if $\lVert\mu(\hat{\pi})-\mu(\pi)\rVert_{\mathcal{H}}$ can be computed efficiently, one typically does not have a good estimate or bound on $\lVert f\rVert_{\mathcal{H}}$. Therefore, it is impractical to use (1.1) to determine an $n$ satisfying the error criterion

[TABLE]

where $\varepsilon$ is the user-specified absolute error tolerance.

We believe that the practitioner would like an automatic cubature, i.e., an algorithm with a stopping criterion that guarantees (1.2) (with high probability). Probabilistic integration, and in particular Bayesian cubature, as espoused in PI, fulfills that wish.

Bayesian cubature, as explained in PI, assumes that the integrand, $f$, may be modeled by a Gaussian stochastic process, $g\sim\mathcal{G}\!\mathcal{P}(0,c)$, conditioned on $g$ having the same values as $f$ at the cubature nodes or states, $\{\boldsymbol{x}_{i}\}_{i=1}^{n}$. Thus, $\hat{\Pi}[g]=\hat{\Pi}[f]$. Furthermore, Bayesian cubature is designed to satisfy $\hat{\Pi}[g]=\operatorname{\mathbb{E}}_{n}[\Pi[g]]$. Here, $c$ is the covariance function (or kernel) for $g$. The definition of $g$ allows us to construct credible intervals for the cubature error via Proposition 1 in PI, namely,

[TABLE]

If the observed integrand, $f$, lies in the $99\%$ middle of the sample space for $g$, and not in the $1\%$ extreme, then increasing $n$ until $2.58\sqrt{\operatorname{\mathbb{V}}_{n}[\Pi[g]]}$ is no greater than $\varepsilon$ ensures that (1.2) holds with $99\%$ probability.

There are some practical obstacles to implementing this elegant recipe.

• •

How does one choose the covariance function $c$? While one may always choose the sample space large enough to include $f$, our use of the credible interval as a stopping criterion assumes that $f$ is not in the tails of the distribution $\mathcal{G}\!\mathcal{P}(0,c)$. We discuss this question in the next section.

• •

The computational cost of computing $\operatorname{\mathbb{V}}_{n}[\Pi[g]]$ involves matrix inversion, which requires $\mathcal{O}(n^{3})$ operations in general. This typically takes much more time than the $\mathcal{O}(n)$ operations required to compute the cubature, $\hat{\Pi}[f]$, unless obtaining an integrand value is quite time-consuming. We discuss how to circumvent this problem by matching covariance functions and cubature nodes in Section 3.

There are commonalities and differences in the deterministic and Bayesian approaches to numerical integration. We discuss some of these in Section 4.

2 Which Gaussian Process?

As mentioned above, using a credible interval as a stopping criterion requires a careful choice of the covariance function, $c$. The width of the credible interval in (1.3) depends on $\operatorname{\mathbb{V}}_{n}[\Pi[g]]$ given by (1.4). At first glance, nothing in (1.4) depends on the integrand data, \boldsymbol{f}=\bigl{(}f(\boldsymbol{x}_{i})\bigr{)}_{i=1}^{n}, although our intuition tells us that it should. The credible interval for the integral of $47f$ should be $47$ times as wide as the credible interval for the integral of $f$.

When constructing the confidence interval for the mean of a scalar random variable, $Y$, from independent and identically distributed (IID) data, one must estimate the variance of $Y$ by the sample variance. Analogously, when constructing the credible interval in (1.3) for the integral (mean) of a function, one must estimate the vertical scale factor inherent in the covariance function $c$.

We have recently explored Bayesian cubature as the basis for automatically selecting $n$ to satisfy the error criterion (1.2) in Jagadeeswaran and Hickernell (2018+), henceforth abbreviated as JH. As in Proposition 2 of PI, JH chooses the covariance function to take the form $c(\boldsymbol{x},\boldsymbol{x}^{\prime})=\lambda c_{0}(\boldsymbol{x},\boldsymbol{x}^{\prime};\boldsymbol{\theta})$, where $\lambda$ is the vertical scale factor, and the parameter $\boldsymbol{\theta}$ determines the smoothness and other properties of the covariance function. An example of $c_{0}$ is the following (JH (36)):

[TABLE]

where $B_{2r}$ is the Bernoulli polynomial of degree $2r$. The smoothness of the covariance function increases with $r$. Covariance functions of this form appear in Hickernell (1996); Dick, Kuo and Sloan (2013). Bernoulli polynomials are described in Chapter 24 of Olver et al. (2018).

To increase the possibility that our integrand $f$ lies in the middle of the sample space, we also allow the Gaussian process $g$ to have an arbitrary mean, $m$, so $g\sim\mathcal{G}\!\mathcal{P}(m,\lambda c_{0})$. One may imagine the situation where $f$ represents an option payoff. Then, $f$ is non-negative and its mean is non-negative. Assuming an improper prior on $(m,\lambda)$, the posterior marginal for $\Pi[g]$ is a Student-t distribution with $n-1$ degrees of freedom and with mean and variance both depending on the integrand data, $\boldsymbol{f}$ (JH (15–16)):

[TABLE]

Here, $\boldsymbol{1}$ is a vector of ones, and $t_{n-1,0.995}$ denotes the $99.5\%$ quantile of the Student-t distribution with $n-1$ degrees of freedom. For large $n$, $t_{n-1,0.995}\approx 2.58$. The expressions in (2.2) are similar to the conclusion of PI, Proposition 2. The differences are due to the mean of the Gaussian process being left unspecified in JH, which reduces the degrees of freedom by one, and adds additional terms to the expressions for $\operatorname{\mathbb{E}}_{n}[\Pi[g]]$ and $\operatorname{\mathbb{V}}_{n}[\Pi[g]]$.

Hidden in the definition of $c_{0}$ is the parameter $\boldsymbol{\theta}$. One may place a discrete prior on $\boldsymbol{\theta}$, but this strikes us as rather arbitrary. Thus, in JH we advocate estimating $\boldsymbol{\theta}$ by empirical Bayes, namely,

[TABLE]

JH also presents empirical Bayes as an alternative to assuming the improper prior on $(m,\lambda)$, as discussed in PI Section 4.1.3. Under empirical Bayes, the posterior marginal for $\Pi[g]$ has the same mean as for the full Bayes approach, but a somewhat smaller variance.

JH also discusses the alternative of generalized cross-validation for estimating the correct covariance function from the integrand data, as alluded to in Section 4.1.2 of PI. The formulas for the Bayesian cubature and the credible interval width are significantly different than those for full Bayes.

3 Speeding Up the Computation

Computing the estimate of $\boldsymbol{\theta}$ in (2.3) and then the credible interval according to (2.2) involve matrix factorization and computing a matrix determinant, which requires as many as $\operatorname{{\mathcal{O}}}(n^{3})$ operations. On the other hand, the computational cost of obtaining the integrand data, $\boldsymbol{f}$, is \operatorname{{\mathcal{O}}}(\(f)n)$, where$$(f)$ is the computational cost of a single integrand value.

If \(f)$is extraordinarily large compared to the expected sample size$n$, then the cost of obtaining integrand data dominates, and the$\operatorname{{\mathcal{O}}}(n^{3})$cost of matrix operations is unimportant. However, if$$(f)$is close to$\operatorname{{\mathcal{O}}}(1)$, then the cost of matrix operations may make Bayesian cubature prohibitively costly.

JH presents a scenario where the cost of matrix operations may be reduced to $\operatorname{{\mathcal{O}}}(n\log n)$ via fast transforms. The key is choosing covariance functions and cubature nodes that match. Let the matrix $\boldsymbol{C}_{0}$ be decomposed in terms of its eigenvectors, which comprise the columns of $\boldsymbol{V}$, and its eigenvalues, which comprise the diagonal elements of the diagonal matrix $\boldsymbol{\Lambda}$:

[TABLE]

Four assumptions are made regarding the covariance function, $c_{0}$, and the cubature nodes, $\{\boldsymbol{x}\}_{i=1}^{n}$ (JH (25, 27)):

[TABLE]

Here, $\boldsymbol{V}^{H}\boldsymbol{b}$ is called the fast transform of $\boldsymbol{b}$ because it takes fewer than the typical $\operatorname{{\mathcal{O}}}(n^{2})$ operations required for matrix-vector multiplication.

An example of matching covariance functions and cubature nodes is

[TABLE]

for some $\mathring{c}_{0}$ with period $1$ in each coordinate direction, and

[TABLE]

The covariance function in (2.1) is an example of a shift-invariant covariance function (Hickernell, 1998). Figure 1 (left) depicts a rank-1 integration node set Sloan and Joe (1994); Dick, Kuo and Sloan (2013). The reason that this family of covariance functions matches this family of cubature nodes and satisfies assumptions (3.1) is that the matrix $\boldsymbol{C}_{0}$ is circulant and $\boldsymbol{V}$ may be written in terms of complex exponentials.

Under assumptions (3.1) one may express (2.2) and (2.3) in terms of the fast transforms of the integrand data and the first column of the matrix $\boldsymbol{C}_{0}$ (JH Sections 3.2, 3.3):

[TABLE]

Apart from the computations in (3.2a), which require $\operatorname{{\mathcal{O}}}(n\log n)$ operations, all other calculations in (3.2) require only $\operatorname{{\mathcal{O}}}(n)$ operations. The expression for $\operatorname{\mathbb{V}}_{n}[\Pi[g]]$ excludes $i=1$ in the sum because we allow $g$ to have an arbitrary constant mean.

Section 5 of JH presents several numerical experiments for Bayesian cubature with shift-invariant covariance functions and lattice nodes sets. We reproduce one such experiment in Figure 1 (right). The integrand is the payoff of an arithmetic mean Asian option. The covariance function is the one given by (2.1) with $r=1$. A low degree of smoothness is chosen in view of the discontinuities in the partial derivatives of the integrand. The value of $\gamma$ is determined by empirical Bayes as in (3.2d). The sample size, $n$, is increased in a sequence of powers of $2$ until the stopping criterion implied by (2.2c), whose terms are computed quickly via (3.2), is satisfied. Four different values of the tolerance were tried, $\varepsilon=0.1,0.01,0.001,$ and $0.0001$. The aim is for the cubature error, $|\Pi[f]-\hat{\Pi}[f]|$, to be no greater than, but not too much less than, the prescribed tolerance nearly all the time. In this experiment, the error tolerance is always met. As expected, the computation time increases as the tolerance decreases.

In this example, and others provided in Section 5 of JH, the cost of evaluating the integrand is modest, and so the cost of obtaining the needed integrand data, $\boldsymbol{f}$, is on the same order as the matrix-vector operations required to compute the credible interval. JH also provides examples of the empirical Bayes and generalized cross-validation approaches to determining the parameters inherent in the covariance function and to using credible intervals as stopping criteria for Bayesian cubature. All of these approaches are successful, which suggests that they should be explored over a larger range of examples.

4 Bayesian Versus Deterministic Analysis

We return to the situation where the Gaussian process, $g$, has zero mean. Section 3.2 in PI sets the covariance function, $c$, identical to the reproducing kernel, $k$, of the Hilbert space containing the integrand, $f$, for “aesthetic” reasons. While this makes the application of several results from numerical analysis of deterministic cubature more readily transferable to Bayesian cubature, we think that such a correspondence muddies the waters.

Suppose that $\{\phi_{i}\}_{i=1}^{\infty}$ is an orthonormal basis for the Hilbert space $\mathcal{H}$ with reproducing kernel $k$. Then

[TABLE]

Moreover, the norm of any $f=\sum_{i=1}^{n}a_{i}\phi_{i}\in\mathcal{H}$ is $\left\lVert f\right\rVert_{\mathcal{H}}=\left\lVert\boldsymbol{a}\right\rVert_{2}$. If $g=\sum_{i=1}^{\infty}A_{i}\phi_{i}$, with $A_{i}\overset{\textup{IID}}{\sim}\mathcal{N}(0,1)$, then $g\in\mathcal{G}\!\mathcal{P}(0,c)$, where $c$ is identical to $k$ as defined in (4.1). While this may seem well and good, note that

[TABLE]

So setting $c$ identical to $k$ means that we are modeling an integrand with finite norm by a Gaussian process with an infinite expected squared norm. This seems counter-intuitive. Nevertheless, there are tantalizing mathematical similarities between the Bayesian and deterministic approaches to cubature.

When the optimal cubature weights are used, the deterministic error bound in (1.1) may be expressed as

[TABLE]

where $\hat{f}$ is the minimum Hilbert space norm interpolant of the integrand, $f$. The reason that $\lVert f\rVert_{\mathcal{H}}$ in (1.1) can be replaced by $\lVert f-\hat{f}\rVert_{\mathcal{H}}$ is that the cubature in (4.2a) integrates $\hat{f}$ exactly. Moreover, $f-\hat{f}$ is orthogonal to $\hat{f}$ under the Hilbert space inner product. Also note that

[TABLE]

Compare the error bound in (4.2) to Proposition 2 of PI, which implies that

[TABLE]

The formulas for the cubature $\hat{\Pi}[f]$ in both the deterministic and Bayesian senses are identical if $k$ is identical to $c_{0}$. They are also independent of the vertical scale factor multiplying inherent in the definition of $k$ or $c_{0}$, i.e., the reproducing kernels $k$ and $47k$ yield the same cubature rule. Likewise, the bound on the squared error in (4.2b) and the variance of the integral of the Gaussian process in (4.4c) both contain the common factor $\Pi\Pi[k(\cdot,\cdot)]-\Pi[\boldsymbol{k}(\cdot,X)]\boldsymbol{K}^{-1}\Pi[\boldsymbol{k}(X,\cdot)]$ if $k$ is identical to $c_{0}$.

Matching the deterministic error bound in (4.2) to the Bayesian credible interval in (4.4) when $k$ is identical to $c_{0}$ becomes possible if one applies (4.3) and asserts that

[TABLE]

Although this inequality is violated for some $f\in\mathcal{H}$, it holds for those $f\in\mathcal{H}$ that are well-modeled by their minimum norm interpolants, $\hat{f}$. Thus, one can mimic Bayesian cubature via a deterministic cubature which assumes that the integrand satisfies inequality (4.5).

5 Further Matters

A couple of matters deserve further investigation. How large a family of covariance functions must be considered for effective Bayesian cubature? A larger family increases the probability that the integrand in question lies in the middle of the space of Gaussian processes used to determine the stopping criterion. On the other hand, a larger family may require a more tedious choice of the underlying parameters $\boldsymbol{\theta}$.

The Bayesian approach to numerical integration assumes a Gaussian process. Do goodness-of-fit statistics confirm or discredit this assumption? How does the validity of this assumption affect the reliability of the proposed Bayesian automatic cubature? The alternative of Student-$t$ processes has been suggested by Shah, Wilson and Ghahramani (2014).

Finally, probabilistic numerics—including Bayesian cubature—deserves further participation from statisticians, numerical analysts, and software developers alike. Statisticians and numerical analysts should become more conversant in each other’s language and culture. Computational problems are better understood when one can look from multiple perspectives. Moreover, the algorithms that arise from probabilistic numerics should find their way into commonly used software libraries. Such libraries should be built using software engineering principles that are familiar to software developers, but perhaps not obvious to statisticians. Our Guaranteed Automatic Integration Library (GAIL) Choi et al. (2013–2017) is an example of such a library.

Acknowledgements

We thank the authors of PI for a provocative article and for urging the computational mathematics community to investigate more carefully probabilistic approaches. We also thank Sou-Cheng Choi for her comments. This work was supported in part by NSF-DMS-1522687 and NSF-DMS-1638521 (SAMSI).