Inverse regression for ridge recovery: A data-driven approach for parameter reduction in computer experiments

TL;DR

This paper interprets inverse regression methods SIR and SAVE as estimating ridge functions for parameter reduction in computer experiments, analyzing their convergence and limitations in identifying true important directions.

Contribution

It provides a new interpretation of SIR and SAVE as ridge function estimators and analyzes their convergence and potential pitfalls in computer experiment contexts.

Findings

SIR and SAVE estimate matrices whose column spaces contain ridge directions.

Convergence of these column spaces improves with more data queries.

SIR and SAVE may misidentify important directions when functions are not true ridge functions.

Abstract

Parameter reduction can enable otherwise infeasible design and uncertainty studies with modern computational science models that contain several input parameters. In statistical regression, techniques for sufficient dimension reduction (SDR) use data to reduce the predictor dimension of a regression problem. A computational scientist hoping to use SDR for parameter reduction encounters a problem: a computer prediction is best represented by a deterministic function of the inputs, so data comprised of computer simulation queries fail to satisfy the SDR assumptions. To address this problem, we interpret SDR methods sliced inverse regression (SIR) and sliced average variance estimation (SAVE) as estimating the directions of a ridge function, which is a composition of a low-dimensional linear transformation with a nonlinear function. Within this interpretation, SIR and SAVE estimate matrices of integrals whose column spaces are contained in the ridge directions' span; we analyze and numerically verify convergence of these column spaces as the number of computer model queries increases. Moreover, we show example functions that are not ridge functions but whose inverse conditional moment matrices are low-rank. Consequently, the computational scientist should beware when using SIR and SAVE for parameter reduction, since SIR and SAVE may mistakenly suggest that truly important directions are unimportant.

Equations

$$y\,\rotatebox[origin={c}]{90.0}{$\models$}\,\mathbf{x}|\bm{A}^{\top}\mathbf{x},$$

$${\mathcal{S}}_{\mathrm{DRS}}={\mathcal{S}}_{\mathrm{DRS}}(\mathbf{A})=\mathrm{colspan}(\mathbf{A}).$$

$${\mathcal{S}}_{y|\mathbf{x}}=\bigcap {\mathcal{S}}_{\text{DRS}}$$

$$\mathbb{E}\left[\mathbf{x}\right]=0,\mathbb{C}\mathrm{ov}\left[\mathbf{x}\right]=\mathbf{I}.$$

$${\mathcal{S}}_{h(y)|\mathbf{x}}\subseteq {\mathcal{S}}_{y|\mathbf{x}},$$

$${\mathbf{C}}_{\mathrm{IR}}=\mathbb{C}\mathrm{ov}\left[\mathbb{E}\left[\mathbf{x}|y\right]\right]$$

$$\mathrm{colspan}({\mathbf{C}}_{\mathrm{IR}})\subseteq {\mathcal{S}}_{y|\mathbf{x}},$$

$$\underset{1\le i\le N}{\min }{y}_i={\tilde{y}}_0<{\tilde{y}}_1<\cdots <{\tilde{y}}_{R-1}<{\tilde{y}}_R=\underset{1\le i\le N}{\max }{y}_i,$$

$$h(y)=r\text{for}y\in J_r.$$

$${\mathbf{C}}_{\mathrm{SIR}}=\mathbb{C}\mathrm{ov}\left[\mathbb{E}\left[\mathbf{x}|h(y)\right]\right].$$

$$\mathrm{colspan}({\mathbf{C}}_{\mathrm{SIR}})\subseteq {\mathcal{S}}_{h(y)|\mathbf{x}}\subseteq {\mathcal{S}}_{y|\mathbf{x}}.$$

$${\hat{\mathbf{\mu }}}_h(r)=\frac{1}{N_r}\sum _{i\in {\mathcal{I}}_r}{\mathbf{x}}_i.$$

$${\hat{\mathbf{C}}}_{\mathrm{SIR}}=\frac{1}{N}\sum _{r=1}^R{N}_r{\hat{\mathbf{\mu }}}_h(r){\hat{\mathbf{\mu }}}_h(r{)}^{\top }.$$

$${\hat{\mathbf{C}}}_{\mathrm{SIR}}=\hat{\mathbf{W}}\hat{\mathbf{\Lambda }}{\hat{\mathbf{W}}}^{\top },$$

$$\mathbb{E}\left[\mathbb{C}\mathrm{ov}\left[\mathbf{x}|y\right]\right]$$

$$=\mathbf{I}-\mathbb{C}\mathrm{ov}\left[\mathbb{E}\left[\mathbf{x}|y\right]\right],$$

$$\mathbb{C}\mathrm{ov}\left[\mathbb{E}\left[\mathbf{x}|y\right]\right]$$

$$=\mathbb{E}\left[\mathbf{I}-\mathbb{C}\mathrm{ov}\left[\mathbf{x}|y\right]\right].$$

$${\mathbf{C}}_{\mathrm{AVE}}=\mathbb{E}\left[{\left(\mathbf{I}-\mathbb{C}\mathrm{ov}\left[\mathbf{x}|y\right]\right)}^2\right].$$

$$\mathrm{colspan}({\mathbf{C}}_{\mathrm{AVE}})\subseteq {\mathcal{S}}_{y|\mathbf{x}}.$$

$${\mathbf{C}}_{\mathrm{SAVE}}=\mathbb{E}\left[{\left(\mathbf{I}-\mathbb{C}\mathrm{ov}\left[\mathbf{x}|h(y)\right]\right)}^2\right].$$

$$\mathrm{colspan}({\mathbf{C}}_{\mathrm{SAVE}})\subseteq {\mathcal{S}}_{h(y)|\mathbf{x}}\subseteq {\mathcal{S}}_{y|\mathbf{x}}.$$

$${\hat{\mathbf{\mu }}}_h(r)=\frac{1}{N_r}\sum _{i\in {\mathcal{I}}_r}{\mathbf{x}}_i.$$

$${\hat{\mathbf{\Sigma }}}_h(r)=\frac{1}{N_r-1}\sum _{i\in {\mathcal{I}}_r}\left({\mathbf{x}}_i-{\hat{\mathbf{\mu }}}_h(r)\right){\left({\mathbf{x}}_i-{\hat{\mathbf{\mu }}}_h(r)\right)}^{\top }$$

$${\hat{\mathbf{C}}}_{\mathrm{SAVE}}=\frac{1}{N}\sum _{r=1}^R{N}_r{\left(\mathbf{I}-{\hat{\mathbf{\Sigma }}}_h(r)\right)}^2.$$

$${\hat{\mathbf{C}}}_{\mathrm{SAVE}}=\hat{\mathbf{W}}\hat{\mathbf{\Lambda }}{\hat{\mathbf{W}}}^{\top },$$

$$y=f(\mathbf{x}),\mathbf{x}\in {\mathbb{R}}^m,y\in \mathbb{R}.$$

$$y=f(\mathbf{x})=g({\mathbf{A}}^{\top }\mathbf{x}).$$

$$f(\mathbf{x}+\mathbf{w})$$

$$=g({\mathbf{A}}^{\top }\mathbf{x}+0)$$

$$=g({\mathbf{A}}^{\top }\mathbf{x})$$

$$=f(\mathbf{x}).$$

$$f(\mathbf{x})=g({\mathbf{A}}^{\top }\mathbf{x}),$$

Full text

1 Introduction and related literature

Advances in computing have enabled complex computer simulations of physical systems that complement traditional theory and experiments. We generically model this type of physics-based input/output system as a function $f(\mathbf{x})$, where $\mathbf{x}$ represents a vector of continuous physical inputs and $f$ represents a continuous physical output of interest. It is common to treat $f$ as a deterministic function, since a bug-free computational model produces the same outputs given the same inputs; in other words, there is no random noise in simulation outputs. This perspective is the de facto setup in computer experiments (Sacks et al., 1989; Koehler and Owen, 1996; Santner et al., 2003). The setup is similar in the field of uncertainty quantification (Smith, 2013; Sullivan, 2015; Ghanem et al., 2016); although the inputs may be modeled as random variables to represent aleatory uncertainty, the input/output map is deterministic.

Scientific studies based on deterministic computer simulations may require multiple evaluations of the model at different input values, which may not be feasible when a single model run is computationally expensive. In this case, one may construct a cheap approximation to the model—sometimes referred to as a response surface (Myers and Montgomery, 1995; Jones, 2001), surrogate model (Razavi et al., 2012; Allaire and Willcox, 2010), metamodel (Wang and Shan, 2006), or emulator (Challenor, 2012)—that can then be sampled thoroughly. However, such constructions suffer from the curse of dimensionality (Donoho, 2000; Traub and Werschulz, 1998); loosely, the number of function queries needed to construct an accurate response surface grows exponentially with the input parameter dimension. In computational science, this curse is exacerbated by the computational cost of each query, which often involves the numerical approximation of a system of partial differential equations. In high dimensions (i.e., several physical input parameters), the number of queries needed to reach error estimate asymptopia (or even construct the response surface) is generally considered infeasible. The computational scientist may attempt to simplify the model by fixing unimportant input parameters at nominal values, thus reducing the dimension of the input space to the point that response surface modeling becomes practical. This concept is generalized by active subspaces (Constantine, 2015), which seek out off-axis anisotropy in the function.

In regression modeling (Weisberg, 2005), the given data are predictor/response pairs $\{[\,\mathbf{x}_{i}^{\top}\,,\,y_{i}\,]\}$ that are assumed to be independent realizations of a random vector with a joint predictor/response probability density function. In this context, subspace-based dimension reduction goes by the name sufficient dimension reduction (SDR) (Cook, 1998; Li, 2018; Adragni and Cook, 2009). Techniques for SDR include sliced inverse regression (SIR) (Li, 1991), sliced average variance estimation (SAVE) (Cook and Weisberg, 1991), ordinary least squares (OLS) (Li and Duan, 1989), and principal Hessian directions (pHd) (Li, 1992)—among several others. These techniques seek a low-dimensional subspace in the predictor space that is sufficient to statistically characterize the relationship between predictors and response. SDR methods are gaining interest for parameter reduction in computational science (Li et al., 2016; Zhang et al., 2017; Pan and Dias, 2017). The first use we know applies OLS, SIR, and pHd to a contaminant transport model from Los Alamos National Labs, where the results revealed simple, exploitable relationships between transformed inputs and the computational model’s output (Cook, 1994b, Section 4). However, these papers do not carefully address the important mathematical differences that arise when applying a technique developed for statistical regression to study a deterministic function.

A data set comprised of point queries from a deterministic function (e.g., an experimental design in the physical inputs and associated computer model predictions) differs from the regression data set (independent draws of a joint predictor/response distribution). For example, the former does not admit a joint density on the input/output space. Thus, the data from computational science simulations fail to satisfy the regression assumptions. This failure has important consequences for interpreting error and convergence in the response surface. We do not explore these consequences in this paper, but the differences between statistical regression and function approximation are essential to our thesis. Instead, we focus on (i) the interpretation of dimension reduction subspaces from SIR or SAVE in the absence of random noise in the functional relationship between inputs and outputs and (ii) the proper meaning of sufficient dimension reduction for deterministic functions.

We evaluate and analyze the inverse regression methods SIR and SAVE as candidates for gradient-free subspace-based parameter reduction in deterministic functions. SIR- and SAVE-style dimension reduction follows from low-rank structure in the regression’s inverse conditional moment matrices. The appropriate way to interpret SIR and SAVE in this context is as methods for ridge recovery. A function of several variables is a ridge function when it is a composition of a low-dimensional linear transformation with a nonlinear function of the transformed variables (Pinkus, 2015)—i.e., $f(\mathbf{x})=g(\bm{A}^{\top}\mathbf{x})$ where $\bm{A}$ is a tall matrix. The goal of ridge recovery is to estimate the matrix $\bm{A}$ (more precisely, $\bm{A}$’s column space known as the ridge subspace) using only point queries of $f$. Recent work by Fornasier et al. (2012) and related work by Tyagi and Cevher (2014) develop ridge recovery algorithms that exploit a connection between a linear measurement of a gradient vector and the function’s directional derivative to estimate $\bm{A}$ with finite differences. Constantine et al. (2015) exploit a similar connection to estimate active subspaces with directional finite differences. In contrast, the inverse regression approach for ridge recovery is based on conditional moments of the predictors given the response; we detail how this difference (gradients versus inverse conditional moments) affects the methods’ ability to recover the ridge directions.

1.1 Main results and discussion

We develop and interpret sufficient dimension reduction theory in the context of parameter reduction for a deterministic scalar-valued function of several variables. By this interpretation, SDR methods translate naturally as ridge recovery methods; loosely, the noiseless analog of SDR is ridge recovery. This interpretation illuminates the theoretical limitations of inverse regression for parameter reduction. We show that, if the data-generating function is a ridge function, then the population (i.e., no error from finite sampling) SIR and SAVE subspaces are (possibly strict) subspaces of the function’s ridge subspace. Consequently, the rank of the inverse conditional moment matrices is bounded above by the dimension of the ridge subspace. We also show how the inverse conditional moments become integrals, so that the finite-sample SIR and SAVE algorithms can be viewed as numerical integration methods. By choosing input samples independently at random, SIR and SAVE become Monte Carlo methods; we analyze the convergence of SIR and SAVE subspaces as the number $N$ of samples increases, and we derive the expected $\mathcal{O}_{p}(N^{-1/2})$ rate, where $\mathcal{O}_{p}$ denotes convergence in probability and the constant depends inversely on associated eigenvalue gap. Moreover, this view enables more efficient numerical integration methods than Monte Carlo for SIR and SAVE (Glaws and Constantine, 2018).

1.2 Practical considerations and scope

Some readers want to know how these results impact practice before (or in lieu of) digesting the mathematical analysis. In practice, the computational scientist does not know whether her model’s input/output map is a ridge function, and she would like to use SIR and SAVE as a tool to test for possible dimension reduction; our theory gives some insight into this use case. First, in the idealized case of no finite-sampling (or numerical integration) error, SIR and SAVE can still mislead the practitioner, since there are non-ridge functions whose inverse conditional moment matrices are low-rank. Arguably, these may be pathological cases that do not arise (frequently) in practice, but quantifying this argument requires mathematically characterizing a typical problem found in practice. We believe these failure modes should be well known and guarded against. Second, finite-sample estimates (or numerical approximations) of the inverse conditional moment matrices will never have trailing eigenvalues that are exactly zero due to the finite sampling—even when the data-generating function is an exact ridge functions. Therefore, it may be difficult to distinguish between finite-sample noise in the eigenvalues and directions of small-but-nonzero variability in the function. The latter structure indicates anisotropic parameter dependence that is arguably more common in real world functions than exact ridge functions, and it (the anisotropy) can be exploited for ridge approximation (as opposed to ridge recovery). Because of such indistinguishability, we recommend exercising caution when using the estimated SIR and SAVE subspaces to build ridge approximations. Principally, heuristics such as bootstrapping or cross validation might help distinguish between finite sample noise and directions of relatively small variability. Third, we have found that breaking the elliptic symmetry assumption on the input density (see Section 4.1)—e.g., when the input density is uniform on a hyperrectangle—may produce nonzero eigenvalues in the inverse conditional moment matrices even when the data-generating function is an exact ridge function; distinguishing these effects from finite-sample effects or true variability in the data-generating function may be difficult.

1.3 Paper outline

The remainder is structured as follows. Section 2 reviews (i) the essential theory of sufficient dimension reduction in statistical regression and (ii) the SIR and SAVE algorithms. In Section 3, we review ridge functions and the ridge recovery problem. Section 4 translates SDR to deterministic ridge recovery and analyzes the convergence of the SIR and SAVE subspaces. Section 5 shows numerical examples that support the analyses using (i) two quadratic functions of ten variables and (ii) a simplified model of magnetohydrodynamics with five input parameters. Section 6 contains concluding remarks. To improve readability, all proofs are provided as supplemental material.

2 Sufficient dimension reduction

We review the essential theory of sufficient dimension reduction (SDR); our notation and development closely follow Cook’s Regression Graphics: Ideas for Studying Regressions through Graphics (Cook, 1998). The theory of SDR provides a framework for subspace-based dimension reduction in statistical regression. A regression problem begins with predictor/response pairs $\{[\,\mathbf{x}_{i}^{\top}\,,\,y_{i}\,]\}$, $i=1,\dots,N$, where $\mathbf{x}_{i}\in\mathbb{R}^{m}$ and $y_{i}\in\mathbb{R}$ denote the vector-valued predictors and scalar-valued response, respectively. These pairs are assumed to be independent realizations from the random vector $[\,\mathbf{x}^{\top}\,,\,y\,]$ with unknown joint probability density function $p_{\mathbf{x},y}$. The object of interest in regression is the conditional random variable $y|\mathbf{x}$; the statistician uses the given predictor/response pairs to estimate statistics of $y|\mathbf{x}$—e.g., moments or quantiles. SDR searches for a subspace of the predictors that is sufficient to describe $y|\mathbf{x}$ with statistics derived from the given predictor/response pairs.

The basic tool of SDR is the dimension reduction subspace (DRS).

Definition 1** **(Dimension reduction subspace

(Cook, 1996)).

Consider a random vector $[\,\mathbf{x}^{\top}\,,\,y\,]$, and let $\bm{A}\in\mathbb{R}^{m\times n}$ with $n\leq m$ be such that

[TABLE]

*where

$\models$

denotes independence of $y$ and $\mathbf{x}$. A dimension reduction subspace $\mathcal{S}_{\operatorname{DRS}}$ for $y|\mathbf{x}$ is*

[TABLE]

Equation (1) denotes the conditional independence of the response and predictors given $\bm{A}^{\top}\mathbf{x}$. We can define a new random variable $y|\bm{A}^{\top}\mathbf{x}$ that is the response conditioned on the $n$-dimensional vector $\bm{A}^{\top}\mathbf{x}$. If $\mathcal{S}_{\operatorname{DRS}}(\bm{A})$ is a DRS for $y|\mathbf{x}$, then the conditional CDF for $y|\mathbf{x}$ is the conditional CDF for $y|\bm{A}^{\top}\mathbf{x}$ (Cook, 1998, Chapter 6). If a regression admits a low-dimensional ($n<m$) DRS, then the predictor dimension can be reduced by considering $y|\bm{A}^{\top}\mathbf{x}$. Note that the matrix $\bm{A}$ in (2) does not uniquely define a DRS (Cook, 1994a). Any matrix with the same column space is a basis for $\mathcal{S}_{\operatorname{DRS}}$

A given regression problem may admit multiple DRSs. We require a uniquely-defined DRS to ensure a well-posed SDR problem. To this end, we consider the central dimension reduction subspace or simply the central subspace.

Definition 2** (Central subspace (Cook, 1996)).**

A subspace $\mathcal{S}_{y|\mathbf{x}}$ is a central subspace for $y|\mathbf{x}$ if it is a DRS and $\mathcal{S}_{y|\mathbf{x}}\subseteq\mathcal{S}_{\operatorname{DRS}}$, where $\mathcal{S}_{\operatorname{DRS}}$ is any DRS for $y|\mathbf{x}$.

When the central subspace exists it is the intersection of all other DRSs,

[TABLE]

for all $\mathcal{S}_{\operatorname{DRS}}$ of $y|\mathbf{x}$. The intersection in (3) always defines a subspace, but this subspace need not satisfy the conditional independence from Definition 1. Therefore, a regression need not admit a central subspace. However, when a central subspace exists, it is the unique DRS of minimum dimension for $y|\mathbf{x}$ (Cook, 1998, Chapter 6). This minimum dimension is referred to as the structural dimension of $y|\mathbf{x}$.

There are a variety of conditions that ensure the existence of $\mathcal{S}_{y|\mathbf{x}}$ for a given regression problem. We consider the following condition on the marginal density $p_{\mathbf{x}}$ of the predictors.

Theorem 1** (Cook (1998)).**

Suppose that $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ are two DRSs for $y|\mathbf{x}$ where the marginal density of the predictors $p_{\mathbf{x}}(\mathbf{x})$ has support over a convex set. Then, $\mathcal{S}_{1}\cap\mathcal{S}_{2}$ is also a DRS.

According to Theorem 1, if $p_{\mathbf{x}}$ has support over a convex set, then any intersection of DRSs will be a DRS—i.e., $\mathcal{S}_{y|\mathbf{x}}$ from (3) is a DRS and (hence) the central subspace. In regression practice, this existence condition may be difficult to verify from the given predictor/response pairs. Existence of the central subspace can be proven under other sets of assumptions, including a generalization of Theorem 1 to M-sets (Yin et al., 2008) and another based on location regressions (Cook, 1998, Ch. 6). The existence criteria in Theorem 1 is the most pertinent when we employ inverse regression for ridge recovery in Section 4.

There are two useful properties of the central subspace that enable convenient transformations. The first involves the effect of affine transformations in the predictor space. Let $\mathbf{z}=\bm{B}\mathbf{x}+\mathbf{b}$ for full rank $\bm{B}\in\mathbb{R}^{m\times m}$ and $\mathbf{x}\in\mathbb{R}^{m}$. If $\operatorname{colspan}(\bm{A})=\mathcal{S}_{y|\mathbf{x}}$ for some $\bm{A}\in\mathbb{R}^{m\times n}$, then $\operatorname{colspan}(\bm{B}^{-\top}\bm{A})=\mathcal{S}_{y|\mathbf{z}}$ (Cook, 1996). This allows us to make the following assumption about the predictor space without loss of generality.

Assumption 1** (Standardized predictors).**

Assume that $\mathbf{x}$ is standardized such that

[TABLE]

The second property involves mappings of the response. Let $h:\mathbb{R}\rightarrow\mathbb{R}$ be a function applied to the responses that produces a new regression problem, $\{[\,\mathbf{x}_{i}^{\top}\,,\,h(y_{i})\,]\}$, $i=1,\dots,N$. The central subspace associated with the new regression problem is contained within the original central subspace,

[TABLE]

with equality holding when $h$ is strictly monotonic (Cook, 2000). Equation (5) is essential for studying the slicing-based algorithms, SIR and SAVE, for estimating the central subspace, where the mapping $h$ partitions the response space; see Sections 2.1 and 2.2.

The goal of SDR is to estimate the central subspace for the regression from the given response/predictor pairs.

Problem 1** (SDR problem).**

Given response/predictor pairs $\{[\,\mathbf{x}_{i}^{\top}\,,\,y_{i}\,]\}$, with $i=1,\dots,N$, assumed to be independent draws from a random vector $[\,\mathbf{x}^{\top}\,,\,y\,]$ with joint density $p_{\mathbf{x},y}$, compute a basis $\bm{A}\in\mathbb{R}^{m\times n}$ for the central subspace $\mathcal{S}_{y|\mathbf{x}}$ of the random variable $y|\mathbf{x}$.

Next, we review two algorithms for Problem 1—sliced inverse regression (Li, 1991) and sliced average variance estimation (Cook and Weisberg, 1991). These algorithms use the given regression data to approximate population moment matrices of the inverse regression $\mathbf{x}|y$.

2.1 Sliced inverse regression

Sliced inverse regression (SIR) (Li, 1991) is an algorithm for approximating the matrix

[TABLE]

using the given predictor/response pairs. $\bm{C}_{\operatorname{IR}}$ is defined by the inverse regression function $\mathbb{E}\left[\mathbf{x}|y\right]$, which draws a curve through the $m$-dimensional predictor space parameterized by the scalar-valued response. If the given regression problem admits an elliptically symmetric marginal density $p_{\mathbf{x}}$, then it satisfies the linearity condition (Eaton, 1986). Under the linearity condition,

[TABLE]

which means the column space of $\bm{C}_{\operatorname{IR}}$ will uncover at least part of the central subspace for the regression (Li, 1991).

SIR approximates $\bm{C}_{\operatorname{IR}}$ using given predictor/response pairs by partitioning (i.e., slicing) the response space. Consider a partition of the observed response space,

[TABLE]

and let $J_{r}=[\tilde{y}_{r-1},\tilde{y}_{r}]$ denote the $r$th partition for $r=1,\dots,R$. Define the function

[TABLE]

Applying $h$ to the given responses creates a new regression problem $\{[\,\mathbf{x}_{i}^{\top}\,,\,h(y_{i})\,]\}$, where (6) becomes

[TABLE]

The notation $\bm{C}_{\operatorname{SIR}}$ emphasizes that this matrix is with respect to the sliced version of the original regression problem. Combining (7) and (5),

[TABLE]

The sliced partition of the response space and the sliced mapping $h$ bin the predictor/response pairs to enable sample estimates of $\mathbb{E}\left[\mathbf{x}|r\right]$ for $r=1,\dots,R$. This is the basic idea behind the SIR algorithm; see Algorithm 1. Note that if the response is discrete, then Algorithm 1 produces the maximum likelihood estimator of the central subspace (Cook and Forzani, 2009).

Eigenvectors of $\bm{C}_{\operatorname{SIR}}$ associated with nonzero eigenvalues provide a basis for the SIR subspace, $\operatorname{colspan}(\bm{C}_{\operatorname{SIR}})$. If the approximated eigenvalues $\hat{\lambda}_{n+1},\dots,\hat{\lambda}_{m}$ from Algorithm 1 are small, then $m\times n$ matrix $\hat{\bm{A}}$ approximates a basis for this subspace. However, determining the appropriate value of $n$ requires care. Li (1991) and Cook and Weisberg (1991) propose significance tests based on the distribution of the average of the $m-n$ trailing estimated eigenvalues. These testing methods also apply to the SAVE algorithm in Section 2.2.

For a fixed number of slices, SIR has been shown to be $N^{-1/2}$-consistent for approximating $\operatorname{colspan}(\bm{C}_{\operatorname{SIR}})$ (Li, 1991). In principle, increasing the number of slices may provide improved estimation of the central DRS. However, in practice, SIR’s success is relatively insensitive to the number of slices. The number of slices should be chosen such that there are enough samples in each slice to estimate the conditional expectations accurately. For this reason, Li (1991) suggests constructing slices such that the response samples are distributed nearly equally.

2.2 Sliced average variance estimation

SAVE uses the variance of the inverse regression. Li (1991) recognized the potential for $\mathbb{C}\operatorname{ov}\left[\mathbf{x}|y\right]$ to provide insights into the central subspace by noting that, under Assumption 1,

[TABLE]

which can be rewritten as

[TABLE]

The left side of (16) is precisely $\bm{C}_{\operatorname{IR}}$. This suggests that $\mathbb{E}\left[\bm{I}-\mathbb{C}\operatorname{ov}\left[\mathbf{x}|y\right]\right]$ may be useful in addressing Problem 1. Cook and Weisberg (1991) suggest using $\mathbb{E}\left[(\bm{I}-\mathbb{C}\operatorname{ov}\left[\mathbf{x}|y\right])^{2}\right]$, which has nonnegative eigenvalues. Define

[TABLE]

Under the linearity condition (Li, 1991) and the constant covariance condition (Cook, 2000), the column space of $\bm{C}_{\operatorname{AVE}}$ is contained within the central subspace,

[TABLE]

Both of these conditions are satisfied when the predictors have an elliptically symmetric marginal density (Eaton, 1986).

Using the partition (8) and the map (9),

[TABLE]

The notation $\bm{C}_{\operatorname{SAVE}}$ indicates application to the sliced version of the original regression problem. Combining (18) and (5),

[TABLE]

Algorithm 2 shows the SAVE algorithm, which computes a basis for the column span of the sample estimate $\hat{\bm{C}}_{\operatorname{SAVE}}$. This basis is a $N^{-1/2}$-consistent estimate of $\operatorname{colspan}(\bm{C}_{\operatorname{SAVE}})$ (Cook, 2000). Increasing the number of slices improves the estimate but suffers from the same drawbacks as SIR. In practice, SAVE performs poorly compared to SIR when few predictor/response pairs are available. This is due to difficulties approximating the covariances within the slices using too few samples. For this reason, Cook and Forzani (2009) suggest trying both methods to approximate the central DRS.

3 Ridge functions

The input/output map of a computer simulation is best modeled by a deterministic function,

[TABLE]

We assume the domain of $f$ is equipped with a known input probability measure $\pi_{\mathbf{x}}$. In an uncertainty quantification context, $\pi_{\mathbf{x}}$ encodes uncertainty in the physical inputs due to, e.g., experimental measurement error. We assume $\pi_{\mathbf{x}}$ admits a density function $p_{\mathbf{x}}:\mathbb{R}^{m}\rightarrow\mathbb{R}^{+}$ and that the $m$ inputs are independent. Common choices for $p_{\mathbf{x}}$ include a multivariate Gaussian and a uniform density over a hyperrectangle. The pair $f$ and $\pi_{\mathbf{x}}$ induces an unknown push-forward probability measure on the output space, which we denote by $\pi_{y}$.

Translating sufficient dimension reduction to deterministic functions naturally leads to ridge functions (Pinkus, 2015).

Definition 3** (Ridge function (Pinkus, 2015)).**

A function $f:\mathbb{R}^{m}\rightarrow\mathbb{R}$ is a ridge function if there exists a matrix $\bm{A}\in\mathbb{R}^{m\times n}$ with $n<m$ and a function $g:\mathbb{R}^{n}\rightarrow\mathbb{R}$ such that

[TABLE]

The columns of $\bm{A}$ are the directions of the ridge function and $g:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is the ridge profile.

A ridge function is nominally a function of $m$ inputs, but it is intrinsically a function of $n<m$ derived inputs. The term ridge function sometimes refers to the $n=1$ case, and $n>1$ is called a generalized ridge functions (Pinkus, 2015). We do not distinguish between the $n=1$ and $n>1$ cases and refer to all such functions as ridge functions for convenience.

Ridge functions appear in a wide range of computational techniques. For example, projection pursuit regression uses a regression model that is a sum of $n$ one-dimensional ridge functions, $\sum_{i=1}^{n}g_{i}(\mathbf{a}_{i}^{\top}\mathbf{x})$, where each $g_{i}$ is a spline or other nonparametric model (Friedman and Stuetzle, 1980). Neural network nodes use functions of the form $\sigma(\bm{W}^{\top}\mathbf{x}+\mathbf{b})$, where $\bm{W}$ is a matrix of weights from the previous layer of neurons, $\mathbf{b}$ is a bias term for the model, and $\sigma(\cdot)$ is an activation function (Goodfellow et al., 2016).

Ridge functions are good low-dimensional models for functions that exhibit off-axis anisotropic dependence on the inputs since they are constant along the $m-n$ directions of the input space that are orthogonal to the columns of $\bm{A}$. Consider a vector $\mathbf{w}\in\operatorname{null}(\bm{A}^{\top})$. Then, for any $\mathbf{x}\in\mathbb{R}^{m}$,

[TABLE]

We refer to the problem of finding the directions of a ridge function as ridge recovery.

Problem 2** (Ridge recovery).**

Given a query-able deterministic function $f:\mathbb{R}^{m}\rightarrow\mathbb{R}$ that is assumed to be a ridge function and an input probability measure $\pi_{\mathbf{x}}$, find $\bm{A}\in\mathbb{R}^{m\times n}$ with $n<m$ such that

[TABLE]

for some $g:\mathbb{R}^{n}\rightarrow\mathbb{R}$.

Recent papers in signal processing and approximation theory propose and analyze algorithms for ridge recovery. Cohen et al. (2012) consider the case of $n=1$ with conditions on the ridge direction. Fornasier et al. (2012) propose an algorithm for general $n<m$ under conditions on $\bm{A}$. Tyagi and Cevher (2014) extend the approximation results from Fornasier et al. to more general $\bm{A}$’s.

Note that the ridge recovery problem is distinct from ridge approximation (Constantine et al., 2017b; Hokanson and Constantine, 2017), where the goal is to find $\bm{A}$ and construct $g$ that minimize the approximation error for a given $f$.111Regarding nomenclature, neither ridge recovery nor ridge approximation is related to ridge regression, where Tikhonov regularization is applied to the regression model coefficients (Hastie et al., 2009, Chapter 3.4). Inverse regression may be useful for ridge approximation or identifying near-ridge structure in a given function, but pursuing these ideas is beyond the scope of this manuscript. Li et al. (2016) distinguish between an “ideal scenario”—which loosely corresponds to ridge recovery and where they can prove certain convergence results for inverse regression—and the “general setting”—where they apply inverse regression as a heuristic in numerical examples.

4 Inverse regression as ridge recovery

This section develops SIR (Algoroithm 1) and SAVE (Algorithm 2) as tools for ridge recovery (Problem 2). The next theorem connects the dimension reduction subspace (Definition 1) to the ridge directions (Definition 3).

Theorem 2**.**

Let $(\Omega,\Sigma,P)$ be a probability triple. Suppose that $\mathbf{x}:\Omega\rightarrow\mathbb{R}^{m}$ and $y:\Omega\rightarrow\mathbb{R}$ are random variables related by a measurable function $f:\mathbb{R}^{m}\rightarrow\mathbb{R}$ so that $y=f(\mathbf{x})$. Let $\bm{A}\in\mathbb{R}^{m\times n}$ be a constant matrix. Then y\,\rotatebox[origin={c}]{90.0}{\models}\,\mathbf{x}|\bm{A}^{\top}\mathbf{x} if and only if $y=g(\bm{A}^{\top}\mathbf{x})$ where $g:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is a measurable function.

Theorem 2 states that conditional independence of the inputs and output in a deterministic function is equivalent to the function being a ridge function. This provides a subspace-based perspective on ridge functions that uses the DRS as the foundation. That is, the directions of the ridge function are relatively unimportant compared to the subspace they span when capturing ridge structure of $f$ with sufficient dimension reduction.

The central subspace (Definition 2) corresponds to the unique subspace of smallest dimension that completely describes the ridge structure of $f$. Our focus on the subspace instead of the precise basis implies that we can assume standardized inputs $\mathbf{x}$ (as in Assumption 1) without loss of generality, which simplifies discussion of the inverse conditional moment matrices that define the SIR and SAVE subspaces for ridge recovery.

Theorem 1 guarantees existence of the central subspace in regression when the marginal density of the predictors has convex support. However, this condition is typically difficult or impossible to verify for the regression problem. In contrast, the deterministic function is accompanied by a known input measure $\pi_{\mathbf{x}}$ with density function $p_{\mathbf{x}}$. In practice, a relatively non-informative measure is used such as a multivariate Gaussian or uniform density on a hyper-rectangle defined by the ranges of each physical input parameter. Such choices for modeling input uncertainty satisfy Theorem 1 and guarantee the existence of the central subspace.

The input probability measure can influence the structure of the central subspace. For example, let $\mathbf{a},\mathbf{b}\in\mathbb{R}^{m}$ be constant vectors pointing in different directions and consider the function

[TABLE]

If $\pi_{\mathbf{x}}$ has a density function with support only for positive values of $\mathbf{x}$ (e.g., uniform over $[0,1]^{m}$), then the central subspace is $\operatorname{span}\{\mathbf{a}\}$. Alternatively, if the input density has support over all of $\mathbb{R}^{m}$ (e.g., multivariate Gaussian), then the central subspace is $\operatorname{span}\{\mathbf{a},\mathbf{b}\}$. For this reason, we write the central subspace for a deterministic function as $\mathcal{S}_{f,\pi_{\mathbf{x}}}$ to emphasize that this subspace is a property of the given function and input probability measure.

For deterministic functions, the translation of the inverse regression $\mathbf{x}|y$ (see Section 2) is the inverse image of $f$ for the output value $y$,

[TABLE]

Unlike the inverse regression $\mathbf{x}|y$ (which is a random vector), $f^{-1}$ is a fixed set determined by $f$’s contours. Furthermore, the inverse image begets the conditional probability measure $\pi_{\mathbf{x}|y}$, which is the restriction of $\pi_{\mathbf{x}}$ to the set $f^{-1}(y)$ (Chang and Pollard, 1997).

4.1 Sliced inverse regression for ridge recovery

For deterministic functions, we can write $\bm{C}_{\operatorname{IR}}$ from (6) as an integral:

[TABLE]

where the conditional expectation over the inverse image $f^{-1}(y)$ is

[TABLE]

This term represents the average of all input values that map to a fixed value of the output.

To understand how $\bm{C}_{\operatorname{IR}}$ can be used for dimension reduction in deterministic functions, consider the following. For $\mathbf{w}\in\mathbb{R}^{m}$ with unit norm,

[TABLE]

If $\mathbf{w}\in\operatorname{null}(\bm{C}_{\operatorname{IR}})$, then one possibility is that $\bm{\mu}(y)$ is orthogonal to $\mathbf{w}$ for all $y$. The following theorem relates this case to possible ridge structure in $f$.

Theorem 3**.**

Let $f:\mathbb{R}^{m}\rightarrow\mathbb{R}$ with input probability measure $\pi_{\mathbf{x}}$ admit a central subspace $\mathcal{S}_{f,\pi_{\mathbf{x}}}$, and assume $\pi_{\mathbf{x}}$ admits an elliptically symmetric and standardized density function. Then, $\operatorname{colspan}(\bm{C}_{\operatorname{IR}})\subseteq\mathcal{S}_{f,\pi_{\mathbf{x}}}$.

Theorem 3 shows that a basis for the range of $\bm{C}_{\operatorname{IR}}$ can be used to estimate $f$’s central subspace. However, this idea has two important limitations. First, by (33), we can write the inner product in the rightmost integrand in terms of the cosine of the angle between the vectors $\bm{\mu}(y)$ and $\mathbf{w}$,

[TABLE]

where $\theta(y)$ is the angle between $\bm{\mu}(y)$ and $\mathbf{w}$. Theorem 3 uses orthogonality of $\bm{\mu}(y)$ and $\mathbf{w}$ (i.e., $\cos(\theta(y))=0$ for all $y$) to show containment of the column space of $\bm{C}_{\operatorname{IR}}$ within the central subspace; however, the integrand in (34) also contains the squared 2-norm of $\bm{\mu}(y)$, which does not depend on $\mathbf{w}$. If $\bm{\mu}(y)=\bm{0}$ for all $y$, then $\mathbf{w}^{\top}\bm{C}_{\operatorname{IR}}\mathbf{w}=0$ for all $\mathbf{w}$. Consider the following example.

Example 1**.**

Assume $\mathbf{x}\in\mathbb{R}^{2}$ is weighted with a bivariate standard Gaussian. Let $y=f(\mathbf{x})=x_{1}x_{2}$. For any value of $y$, $\mathbf{x}\in f^{-1}(y)$ implies $-\mathbf{x}\in f^{-1}(y)$. Therefore, $\bm{\mu}(y)=\bm{0}$ for all $y$, and $\bm{C}_{\operatorname{IR}}=\bm{0}$. But $y$ is not constant over $\mathbb{R}^{2}$. Thus, $\{\bm{0}\}=\operatorname{colspan}(\bm{C}_{\operatorname{IR}})\subset\mathcal{S}_{f,\pi_{\mathbf{x}}}=\mathbb{R}^{2}$.

This example shows how $\bm{C}_{\operatorname{IR}}$, as a tool for ridge recovery, can mislead the practitioner by suggesting ridge structure in a function that is not a ridge function. This could lead one to ignore input space directions that should not be ignored. Note that if we shift the function such that $y=f(\mathbf{x})=(x_{1}+c_{1})(x_{2}+c_{2})$ for some constants $c_{1},c_{2}\neq 0$, then the symmetry is broken and $\bm{\mu}(y)\neq 0$ for all $y$. In this case, $\bm{C}_{\operatorname{IR}}$ will recover the central subspace (i.e., all of $\mathbb{R}^{2}$).

The second limitation of $\bm{C}_{\operatorname{IR}}$ for ridge recovery follows from the required elliptic symmetry of the input density $p_{\mathbf{x}}$. This assumption is satisfied if $p_{\mathbf{x}}$ is a multivariate Gaussian, but it is violated when the density is uniform over the $m$-dimension hypercube. If $f$ is a ridge function and $\mathbf{w}\in\operatorname{null}(\bm{A}^{\top})$, then $\mathbf{x}\in f^{-1}(y)$ implies $\mathbf{x}+\mathbf{w}\in f^{-1}(y)$ so that $f^{-1}(y)$ can be expressed as the union of lines parallel to $\mathbf{w}$. If the inputs are weighted by an elliptically symmetric density, then the expectation over $f^{-1}(y)$ will be centered such that $\bm{\mu}(y)$ is orthogonal to $\mathbf{w}$. If the inputs do not have an elliptically symmetric density, then the weighting can cause the conditional expectation to deviate in the direction of $\mathbf{w}$. The magnitude of this deviation also depends on the magnitude of the conditional expectation $||\bm{\mu}(y)||_{2}$.

Next, we examine the sliced approximation of $\bm{C}_{\operatorname{IR}}$. Recall the output partition from (8) and the slicing map $h(y)$ (9). Applying $h$ to the deterministic function $y=f(\mathbf{x})$ produces the discretized function $r=h(y)=h(f(\mathbf{x}))$, where $r\in\{1,\dots,R\}$. The output space of $h$ is weighted by the probability mass function

[TABLE]

Without loss of generality, we assume that the slices are constructed such that $\omega(r)>0$ for all $r$. If $\omega(r)=0$ for some $r$, then we can combine this slice with an adjacent slice. The conditional expectation for the sliced output is

[TABLE]

where $\pi_{\mathbf{x}|r}$ is the conditional measure defined over the set $f^{-1}(h^{-1}(r))=\{\,\mathbf{x}\in\mathbb{R}^{m}\,:\,h(f(\mathbf{x}))=h(y)=r\,\}$. Using (35) and (36), the sliced version of $\bm{C}_{\operatorname{IR}}$ is

[TABLE]

By Theorem 2, properties of the central subspace extend to the ridge recovery problem. This includes containment of the central subspace under any mapping of the output,

[TABLE]

By approximating $\bm{C}_{\operatorname{SIR}}$, we obtain an approximation of part of $f$’s central subspace. An important corollary of (38) is that the rank of $\bm{C}_{\operatorname{SIR}}$ is bounded above by the dimension of $\mathcal{S}_{f,\pi_{\mathbf{x}}}$.

Note that $\bm{C}_{\operatorname{SIR}}$ from (37) is a finite sum approximation of the integral in $\bm{C}_{\operatorname{IR}}$ from (31). Since $f^{-1}(h^{-1}(r))=\cup_{y\in J_{r}}f^{-1}(y)$, then $\bm{\mu}_{h}(r)$ is the average of the conditional expectations with $y\in J_{r}$. That is,

[TABLE]

Therefore, $\bm{C}_{\operatorname{SIR}}$ approximates $\bm{C}_{\operatorname{IR}}$ by a weighted sum of the average values of $\bm{\mu}(y)$ within each slice. If $\bm{\mu}(y)$ is continuous almost everywhere with respect to $\pi_{y}$, then $\bm{C}_{\operatorname{IR}}$ is Riemann-integrable (Folland, 1999, Ch. 2). This ensures that sum approximations using the supremum and infimums of $\bm{\mu}(y)$ over each slice converge to the same value. By the sandwich theorem, the average value will converge to this value as well (Abbott, 2001). Therefore, $\bm{C}_{\operatorname{SIR}}$ is a Riemann sum approximation of $\bm{C}_{\operatorname{IR}}$; as the number of slices $R$ increases, $\bm{C}_{\operatorname{SIR}}$ converges to $\bm{C}_{\operatorname{IR}}$.

We turn attention to asymptotic convergence of Algorithm 1 for ridge recovery. To generate the data for Algorithm 1, we choose $N$ points $\{\mathbf{x}_{i}\}$ in the input space consistent with $\pi_{\mathbf{x}}$. For each $\mathbf{x}_{i}$, we query the function to produce the corresponding output $y_{i}=f(\mathbf{x}_{i})$. In the computational science context, this corresponds to running the simulation model at particular sets of inputs. If we choose each $\mathbf{x}_{i}$ independently according to $\pi_{\mathbf{x}}$, then we can analyze SIR as a Monte Carlo method for estimating $\bm{C}_{\operatorname{SIR}}$ from (37). Given the input/output pairs, Algorithm 1 constructs the random matrix $\hat{\bm{C}}_{\operatorname{SIR}}$. To be clear, $\hat{\bm{C}}_{\operatorname{SIR}}$ is a random estimate of $\bm{C}_{\operatorname{SIR}}$ because of how we chose the points $\{\mathbf{x}_{i}\}$—not because of any randomness in the map $f$. Eigenpairs derived from $\hat{\bm{C}}_{\operatorname{SIR}}$ are also random, and the convergence analysis for Algorithm 1 is probabilistic.

The convergence depends on the smallest number of samples per slice over all the slices:

[TABLE]

where $N_{r}$ is from Algorithm 1. Recall that the slices are assumed to be constructed such that $\omega(r)>0$. Thus, $N_{r_{\min}}>0$ with probability 1 as $N\rightarrow\infty$. The following theorem shows that the eigenvalues of $\hat{\bm{C}}_{\operatorname{SIR}}$ converge to those of $\bm{C}_{\operatorname{SIR}}$ in a mean-squared sense.

Theorem 4**.**

Assume that Algorithm 1 has been applied to the data set $\{[\,\mathbf{x}_{i}^{\top}\,,\,y_{i}\,]\}$, with $i=1,\dots,N$, where the $\mathbf{x}_{i}$ are drawn independently according to $\pi_{\mathbf{x}}$ and $y_{i}=f(\mathbf{x}_{i})$ are point evaluations of $f$. Then, for $k=1,\dots,m$,

[TABLE]

where $\lambda_{k}(\cdot)$ denotes the $k$th eigenvalue of the given matrix.

In words, the mean-squared error in the eigenvalues of $\hat{\bm{C}}_{\operatorname{SIR}}$ decays at a $N_{r_{\min}}^{-1}$ rate. Since $\omega(r)>0$ for all $r$, $N_{r_{\min}}\rightarrow\infty$ as $N\rightarrow\infty$. Moreover, the convergence rate suggests that one should choose the slices in Algorithm 1 such that the same number of samples appears in each slice. This maximizes $N_{r_{\min}}$ and reduces the error in the eigenvalues.

An important consequence of Theorem 4 is that the column space of the finite-sample $\hat{\bm{C}}_{\operatorname{SIR}}$ is not contained in $f$’s central subspace. In fact, due to finite sampling, $\hat{\bm{C}}_{\operatorname{SIR}}$ is not precisely low-rank. With a fixed number of samples, it is difficult to distinguish the effects of finite sampling from actual variability in $f$. However, the eigenvalue convergence implies that one can devise practical convergence tests for low-rank-ness based on sets of samples with increasing size.

The next theorem shows the value of understanding the approximation errors in the eigenvalues for quantifying the approximation errors in the subspaces. We measure convergence of the subspace estimates using the subspace distance (Golub and Van Loan, 2013, Chapter 2.5),

[TABLE]

where $\bm{A},\hat{\bm{A}}$ are the first $n$ eigenvectors of $\bm{C}_{\operatorname{SIR}}$ and $\hat{\bm{C}}_{\operatorname{SIR}}$, respectively. The distance metric (42) is the principal angle between the subspaces $\operatorname{colspan}(\bm{A})$ and $\operatorname{colspan}(\hat{\bm{A}})$.

Theorem 5**.**

Assume the same conditions from Theorem 4, and let $\Delta_{n}=\lambda_{n}(\bm{C}_{\operatorname{SIR}})-\lambda_{n+1}(\bm{C}_{\operatorname{SIR}})$ denote the gap between the $n$th and $(n+1)$th eigenvalues of $\bm{C}_{\operatorname{SIR}}$. Then, for sufficiently large $N$,

[TABLE]

where $\bm{A},\hat{\bm{A}}$ are the first $n$ eigenvectors of $\bm{C}_{\operatorname{SIR}}$ and $\hat{\bm{C}}_{\operatorname{SIR}}$, respectively.

The subspace error decays with asymptotic rate $N_{r_{\min}}^{-1/2}$. The more interesting result from Theorem 5 is the inverse relationship between the subspace error and the magnitude of the gap between the $n$th and $(n+1)$th eigenvalues. That is, a large gap between eigenvalues suggests a better estimate of the subspace for a fixed number of samples. We do not hide this factor in the $\mathcal{O}$ notation to emphasize the importance of Theorem 4, which provides insights into the accuracy of the estimated eigenvalues of $\bm{C}_{\operatorname{SIR}}$.

4.2 SAVE for ridge recovery

Similar to (31), we express $\bm{C}_{\operatorname{AVE}}$ from (17) as an integral,

[TABLE]

The conditional covariance $\bm{\Sigma}(y)$ in (44) is an integral against the conditional probability measure $\pi_{\mathbf{x}|y}$,

[TABLE]

To see the relationship between the $\bm{C}_{\operatorname{AVE}}$ matrix and ridge functions, let $\mathbf{w}\in\mathbb{R}^{m}$ with unit norm:

[TABLE]

Equation (46) relates the column space of $\bm{C}_{\operatorname{AVE}}$ to the ridge structure in $f$.

Theorem 6**.**

Let $f:\mathbb{R}^{m}\rightarrow\mathbb{R}$ with input probability measure $\pi_{\mathbf{x}}$ admit a central subspace $\mathcal{S}_{f,\pi_{\mathbf{x}}}$, and assume $\pi_{\mathbf{x}}$ admits an elliptically symmetric and standardized density function. Then, $\operatorname{colspan}(\bm{C}_{\operatorname{AVE}})\subseteq\mathcal{S}_{f,\pi_{\mathbf{x}}}$.

This result shows the usefulness of $\bm{C}_{\operatorname{AVE}}$ for revealing ridge structure in deterministic functions: by estimating the column space of $\bm{C}_{\operatorname{AVE}}$, we obtain an estimate of a subspace of $f$’s central subspace. However, $\bm{C}_{\operatorname{AVE}}$ suffers two similar pitfalls as $\bm{C}_{\operatorname{IR}}$. First, $\bm{C}_{\operatorname{AVE}}$ can mislead the practitioner by suggesting ridge structure that does not exist—i.e., when $\operatorname{colspan}(\bm{C}_{\operatorname{AVE}})\subset\mathcal{S}_{f,\pi_{\mathbf{x}}}$—as the following example illustrates.

Example 2**.**

Assume $\mathbf{x}\in\mathbb{R}^{2}$ is weighted by a bivariate standard Gaussian. Let

[TABLE]

for some $0<r_{1}<r_{2}$. This functions looks like a bullseye with the central circle and outer ring mapping to $y_{1}$ and the middle ring mapping to $y_{2}$. If we choose $r_{1}$ and $r_{2}$ appropriately, then we can obtain

[TABLE]

Note that $\bm{\mu}(y_{1})=\bm{\mu}(y_{2})=0$ for all choices of $r_{1}$ and $r_{2}$. Then,

[TABLE]

and

[TABLE]

Thus, (48) holds when $r_{1}^{2}e^{-r_{1}^{2}/2}=r_{2}^{2}e^{-r_{2}^{2}/2}$, provided that $0<r_{1}<r_{2}$. Choosing $r_{1}$ and $r_{2}$ that satisfy these requirements results in

[TABLE]

However, we can see from inspection of (47) that $\mathcal{S}_{f,\pi_{\mathbf{x}}}=\mathbb{R}^{2}$. Thus, $\operatorname{colspan}(\bm{C}_{\operatorname{AVE}})\subset\mathcal{S}_{f,\pi_{\mathbf{x}}}$.

Example 2 shows one way that $\bm{C}_{\operatorname{AVE}}$ can be fooled into falsely suggesting low-dimensional structure in a function. Note that the rotational symmetry in $f(\mathbf{x})$ is the key feature of this function that results in a degenerate $\bm{C}_{\operatorname{AVE}}$ matrix. This symmetry also tricks the $\bm{C}_{\operatorname{IR}}$ matrix. In fact, it can be shown that, in general,

[TABLE]

which suggests that this sort of false positive is less likely to occur with $\bm{C}_{\operatorname{AVE}}$ than with $\bm{C}_{\operatorname{IR}}$ (Li, 2018, Chap. 5). The exhaustiveness of $\bm{C}_{\operatorname{AVE}}$ in capturing central subspace can be proven provided that at least one of $\mathbb{E}\left[\mathbf{w}^{\top}\mathbf{x}|y\right]$ and $\operatorname{Var}\left[\mathbf{w}^{\top}\mathbf{x}|y\right]$ are nondegenerate—i.e., explicitly depends on $y$—for all $\mathbf{w}\in\mathcal{S}_{f,\pi_{\mathbf{x}}}$. Notice that the function described in Example 2 violates this exhaustiveness condition.

The second limitation of $\bm{C}_{\operatorname{AVE}}$ arises from the elliptic symmetry requirement on the input density $p_{\mathbf{x}}$. When $p_{\mathbf{x}}$ is not elliptically symmetric, we cannot guarantee that the column space of $\bm{C}_{\operatorname{AVE}}$ is contained within the central subspace. Thus, a basis for the column space of $\bm{C}_{\operatorname{AVE}}$ could be contaminated by effects of $\pi_{\mathbf{x}}$—independent of whether or not $f$ is a ridge function.

Next, we consider the sliced version of $\bm{C}_{\operatorname{AVE}}$. We use the same slicing function $h$ from (9) to approximate $\bm{C}_{\operatorname{AVE}}$ from (45). The sliced approximation of $\bm{C}_{\operatorname{AVE}}$ is

[TABLE]

where $\omega(r)$ is the probability mass function from (35) and

[TABLE]

By containment of the central subspace,

[TABLE]

We can interpret $\bm{C}_{\operatorname{SAVE}}$ as a Riemann sum approximation of $\bm{C}_{\operatorname{AVE}}$ using a similar argument as in Section 4.1. An important corollary of (55) is that the rank of $\bm{C}_{\operatorname{SAVE}}$ is bounded above by the dimension of $f$’s central subspace.

Algorithm 2 computes a sample approximation of $\bm{C}_{\operatorname{SAVE}}$, denoted $\hat{\bm{C}}_{\operatorname{SAVE}}$, using given input/output pairs $\{[\,\mathbf{x}_{i}^{\top}\,,\,y_{i}\,]\}$. When the $\mathbf{x}_{i}$ are sampled independently according to $\pi_{\mathbf{x}}$ and each $y_{i}=f(\mathbf{x}_{i})$ is a deterministic function query, we can interpret $\hat{\bm{C}}_{\operatorname{SAVE}}$ as a Monte Carlo approximation to $\bm{C}_{\operatorname{SAVE}}$. Thus, $\hat{\bm{C}}_{\operatorname{SAVE}}$ and its eigenpairs are random—not because of any randomness in the map $f$ but because of the random choices of $\mathbf{x}_{i}$. The following theorem shows the rate of mean-squared convergence of the eigenvalues of $\hat{\bm{C}}_{\operatorname{SAVE}}$.

Theorem 7**.**

Assume that Algorithm 2 has been applied to the data set $\{[\,\mathbf{x}_{i}^{\top}\,,\,y_{i}\,]\}$, with $i=1,\dots,N$, where the $\mathbf{x}_{i}$ are drawn independently according to $\pi_{\mathbf{x}}$ and $y_{i}=f(\mathbf{x}_{i})$ are point evaluations of $f$. Then, for $k=1,\dots,m$,

[TABLE]

where $\lambda_{k}(\cdot)$ denotes the $k$th eigenvalue of the given matrix.

We note that the column space of $\hat{\bm{C}}_{\operatorname{SAVE}}$ is not contained in $f$’s central subspace because of finite sampling. However, using a sequence of estimates with increasing $N$, one may be able to distinguish effects of finite sampling from true directions of variability in $f$.

Next, we examine the convergence of the subspaces Algorithm 2 produces, where the subspace distance is from (42).

Theorem 8**.**

Assume the same conditions from Theorem 7, and let $\Delta_{n}=\lambda_{n}(\bm{C}_{\operatorname{SAVE}})-\lambda_{n+1}(\bm{C}_{\operatorname{SAVE}})$ denote the gap between the $n$th and $(n+1)$th eigenvalues of $\bm{C}_{\operatorname{SAVE}}$. Then, for sufficiently large $N$,

[TABLE]

where $\bm{A},\hat{\bm{A}}$ are the first $n$ eigenvectors of $\bm{C}_{\operatorname{SAVE}}$ and $\hat{\bm{C}}_{\operatorname{SAVE}}$, respectively.

The subspace error for Algorithm 2 decays asymptotically like $N_{r_{\min}}^{-1/2}$ with high probability. Similar to the estimated SIR subspace from Algorithm 1, the error depends inversely on the eigenvalue gap. If the gap between the $n$th and $(n+1)$th eigenvalues is large, then the error in the estimated $n$-dimensional subspace is relatively small for a fixed number of samples.

5 Numerical results

We apply SIR and SAVE to three ridge functions to study the methods’ applicability for ridge recovery and verify our convergence analysis. Since we are not concerned with using SIR and SAVE for ridge approximation, we do not study the methods’ behavior for functions that are not ridge functions. For each ridge function, we provide several graphics including convergence plots of estimated eigenvalues, eigenvalue errors, and subspace errors. One graphic is especially useful for visualizing the structure of the function relative to its central subspace: the sufficient summary plot (Cook, 1998). Sufficient summary plots show $y$ versus $\bm{A}^{\top}\mathbf{x}$, where $\bm{A}$ has only one or two columns that comprise a basis for the central subspace. Algorithms 1 and 2 produce eigenvectors that span an approximation of the SIR and SAVE subspaces, respectively. We use these eigenvectors to construct the low-dimensional inputs for sufficient summary plots. The label sufficient is tied to the precise definition of statistical sufficiency in sufficient dimension reduction for regression.

The convergence analysis from Sections 4.1 and 4.2 assume a fixed slicing of the observed $y$ range. However, Algorithms 1 and 2 are implemented using an adaptive slicing approach that attempts to maximize $N_{r_{\min}}$ for a given set of data. This is done as a heuristic technique for reducing eigenvalue and subspace errors.

The Python code used to generate the figures is available at https://bitbucket.org/aglaws/inverse-regression-for-ridge-recovery. The scripts require the dev branch of the Python Active-subspaces Utility Library (Constantine et al., 2016).

5.1 One-dimensional quadratic ridge function

We study a simple one-dimensional quadratic ridge function to contrast the recovery properties of SIR versus SAVE. Let $\pi_{\mathbf{x}}$ have a standard multivariate Gaussian density function on $\mathbb{R}^{10}$. Define

[TABLE]

where $\mathbf{b}\in\mathbb{R}^{10}$ is a constant vector. The span of $\mathbf{b}$ is the central subspace. First, we attempt to estimate the central subspace using SIR (Algorithm 1), which is known to fail for functions symmetric about $\mathbf{x}=\bm{0}$ (Cook and Weisberg, 1991); Figure 1 confirms this failure. In fact, $\bm{C}_{\operatorname{IR}}$ is zero since the conditional expectation of $\mathbf{x}$ for any value of $y$ is zero. Figure 1a shows that all estimated eigenvalues of the SIR matrix are nearly zero as expected. Figure 1b is a one-dimensional sufficient summary plot of $y_{i}$ against $\hat{\mathbf{w}}_{1}^{\top}\mathbf{x}_{i}$, where $\hat{\mathbf{w}}_{1}$ denotes the normalized eigenvector associated with the largest eigenvalue of $\hat{\bm{C}}_{\operatorname{SIR}}$ from Algorithm 1. If (i) the central subspace is one-dimensional (as in this case) and (ii) the chosen SDR algorithm correctly identifies the one basis vector, then the sufficient summary plot will show a univariate relationship between the linear combination of input evaluations and the associated outputs. Due to the symmetry in the quadratic function, SIR fails to recover the basis vector; the sufficient summary plot’s lack of univariate relationship confirms the failure.

Figure 2 shows results from applying SAVE (Algorithm 2) to the quadratic function (58). Figure 2a shows the eigenvalues of $\hat{\bm{C}}_{\operatorname{SAVE}}$ from Algorithm 2. Note the large gap between the first and second eigenvalues, which suggests that the SAVE subspace is one-dimensional. Figure 2b shows the sufficient summary plot using the first eigenvector $\hat{\mathbf{w}}_{1}$ from Algorithm 2, which reveals the univariate quadratic relationship between $\hat{\mathbf{w}}_{1}^{\top}\mathbf{x}$ and $y$.

5.2 Three-dimensional quadratic function

Next, we numerically study the convergence properties of the SIR and SAVE algorithms using a more complex quadratic function. Let $\pi_{\mathbf{x}}$ have a standard multivariate Gaussian density function on $\mathbb{R}^{10}$. Define

[TABLE]

where $\bm{B}\in\mathbb{R}^{10\times 2}$ and $\mathbf{b}\in\mathbb{R}^{10}$ with $\mathbf{b}\not\in\operatorname{colspan}(\bm{B})$. Figure 3a shows the eigenvalues of $\hat{\bm{C}}_{\operatorname{SIR}}$; note the gap between the third and fourth eigenvalues. Figure 3b shows the maximum squared eigenvalue error normalized by the largest eigenvalue,

[TABLE]

for increasing numbers of samples in 10 independent trials. We estimate the true eigenvalues using SIR with $10^{7}$ samples. The average error decays at a rate slightly faster than the $\mathcal{O}(N^{-1})$ from Theorem 4. The improvement can likely be attributed to the adaptive slicing procedure discussed at the beginning of this section. Figure 3c shows the error in the estimated three-dimensional SIR subspace (see (42)) for increasing numbers of samples. We use $10^{7}$ samples to estimate the true SIR subspace. The subspace errors decrease asymptotically at a rate of approximately $\mathcal{O}(N^{-1/2})$, which agrees with Theorem 5.

Figure 4 shows the results of a similar convergence study using SAVE (Algorithm 2). The eigenvalues of $\hat{\bm{C}}_{\operatorname{SAVE}}$ from (23) are shown in Figure 4a. Note the large gap between the third and fourth eigenvalues, which is consistent with the three-dimensional central subspace in $f$ from (59). Figures 4b and 4c show the maximum squared eigenvalue error and the subspace error, respectively, for $n=3$. The eigenvalue error again decays at a faster rate than expected in Theorem 7—likely due to the adaptive slicing implemented in the code. The subspace error decays consistently according to Theorem 8.

5.3 Hartmann problem

The following study moves toward the direction of using SIR and SAVE for parameter reduction in a physics-based model. The Hartmann problem is a standard problem in magnetohydrodynamics (MHD) that models the flow of an electrically-charged plasma in the presence of a uniform magnetic field (Cowling and Lindsay, 1957). The flow occurs along an infinite channel between two parallel plates separated by distance $2\ell$. The applied magnetic field is perpendicular to the flow direction and acts as a resistive force on the flow velocity. At the same time, the movement of the fluid induces a magnetic field along the direction of the flow. Figure 5 contains a diagram of this problem. We have recently used the following model as a test case for parameter reduction methods (Glaws et al., 2017).

The inputs to the Hartmann model are fluid viscosity $\mu$, fluid density $\rho$, applied pressure gradient $\partial p_{0}/\partial x$ (where the derivative is with respect to the flow field’s spatial coordinate), resistivity $\eta$, and applied magnetic field $B_{0}$. We collect these inputs into a vector,

[TABLE]

The output of interest is the total induced magnetic field,

[TABLE]

This function is not a ridge function of $\mathbf{x}$. However, it has been shown that many physical laws can be expressed as ridge functions by considering a log transform of the inputs, which relates ridge structure in the function to dimension reduction via the Buckingham Pi Theorem (Constantine et al., 2017a). For this reason, we apply SIR and SAVE as ridge recovery methods for $B_{\text{ind}}$ as a function of the logarithms of the inputs from (61). The log-transformed inputs are equipped with a multivariate Gaussian with mean $\bm{\mu}$ and covariance $\bm{\Sigma}$,

[TABLE]

Figure 6 shows the results of applying SIR (Algorithm 1) to the Hartmann model for the induced magnetic field $B_{\text{ind}}$. The eigenvalues of $\hat{\bm{C}}_{\operatorname{SIR}}$ from (13) with bootstrap ranges are shown in Figure 6a. Large gaps appear after the first and second eigenvalues, which indicates possible two-dimensional ridge structure. In fact, the induced magnetic field admits a two-dimensional central subspace relative to the log-inputs (Glaws et al., 2017). Figures 6c and 6d contain one- and two-dimensional sufficient summary plots of $B_{\text{ind}}$ against $\hat{\mathbf{w}}_{1}^{\top}\mathbf{x}$ and $\hat{\mathbf{w}}_{2}^{\top}\mathbf{x}$, where $\hat{\mathbf{w}}_{1}$ and $\hat{\mathbf{w}}_{2}$ are the first two eigenvectors of $\hat{\bm{C}}_{\operatorname{SIR}}$. We see a strong one-dimensional relationship. However, the two-dimensional sufficient summary plot shows slight curvature with changes in $\hat{\mathbf{w}}_{2}^{\top}\mathbf{x}$. These results suggest that ridge-like structure may be discovered using the SIR algorithm in some cases. Figure 6b shows the subspace errors as a function of the subspace dimension. Recall from Theorem 5 that the subspace error depends inversely on the eigenvalue gap. The largest eigenvalue gap occurs between the first and second eigenvalues, which is consistent with the smallest subspace error for $n=1$.

We perform the same numerical studies using SAVE. Figure 7a shows the eigenvalues of the $\hat{\bm{C}}_{\operatorname{SAVE}}$ from (23) for the induced magnetic field $B_{\text{ind}}$ from (62). Note the large gaps after the first and second eigenvalues. These gaps are consistent with the subspace errors in Figure 7b, where the one- and two-dimensional subspace estimates have the smallest errors. Figures 7c and 7d contain sufficient summary plots for $\hat{\mathbf{w}}_{1}^{\top}\mathbf{x}$ and $\hat{\mathbf{w}}_{2}^{\top}\mathbf{x}$, where $\hat{\mathbf{w}}_{1}$ and $\hat{\mathbf{w}}_{2}$ are the first two eigenvectors from $\hat{\bm{C}}_{\operatorname{SAVE}}$ in (23).

6 Summary and conclusion

We investigate sufficient dimension reduction from statistical regression as a tool for subspace-based parameter reduction in deterministic functions where applications of interest include computer experiments. We show that SDR is theoretically justified as a tool for ridge recovery by proving equivalence of the dimension reduction subspace and the ridge subspace for some deterministic $y=f(\mathbf{x})$. We interpret two SDR algorithms for the ridge recovery problem: sliced inverse regression and sliced average variance estimation. In regression, these methods use moments of the inverse regression $\mathbf{x}|y$ to estimate subspaces relating to the central subspace. In ridge recovery, we reinterpret SIR and SAVE as numerical integration methods for estimating inverse conditional moment matrices, where the integrals are over contour sets of $f$. We show that the column spaces of the conditional moment matrices are contained in the ridge subspace, which justifies their eigenspaces as tools ridge recovery.