A Lipschitz Matrix for Parameter Reduction in Computational Science

TL;DR

The paper introduces the Lipschitz matrix, a generalization of the Lipschitz constant, to improve function approximation, integration, and optimization in high-dimensional computational science models by reducing uncertainty and enabling parameter reduction.

Contribution

It presents the Lipschitz matrix as a new tool for analyzing complex functions, enabling dimension-independent costs and parameter reduction through low-rank approximations.

Findings

Lipschitz matrix reduces worst-case approximation and integration costs.

Low-rank Lipschitz matrices enable dimension-independent analysis.

Parameter reduction via ridge approximation using dominant eigenvectors.

Abstract

We introduce the Lipschitz matrix: a generalization of the scalar Lipschitz constant for functions with many inputs. Among the Lipschitz matrices compatible a particular function, we choose the smallest such matrix in the Frobenius norm to encode the structure of this function. The Lipschitz matrix then provides a function-dependent metric on the input space. Altering this metric to reflect a particular function improves the performance of many tasks in computational science. Compared to the Lipschitz constant, the Lipschitz matrix reduces the worst-case cost of approximation, integration, and optimization; if the Lipschitz matrix is low-rank, this cost no longer depends on the dimension of the input, but instead on the rank of the Lipschitz matrix defeating the curse of dimensionality. Both the Lipschitz constant and matrix define uncertainty away from point queries of the function and by using the Lipschitz matrix we can reduce uncertainty. If we build a minimax space-filling design of experiments in the Lipschitz matrix metric, we can further reduce this uncertainty. When the Lipschitz matrix is approximately low-rank, we can perform parameter reduction by constructing a ridge approximation whose active subspace is the span of the dominant eigenvectors of the Lipschitz matrix. In summary, the Lipschitz matrix provides a new tool for analyzing and performing parameter reduction in complex models arising in computational science.

Tables

Table 1. Table 1: Comparison of Lipschitz matrix, Lipschitz constant, and corresponding volumes for a variety of test functions  [ 32 ] posed on the normalized domain 𝒟 = [ − 1 , 1 ] m 𝒟 superscript 1 1 𝑚 {\mathcal{D}}=[-1,1]^{m} with normalized outputs in [ 0 , 1 ] 0 1 [0,1] . Here we estimate 𝐋 𝐋 {\mathbf{L}} by minimizing the Frobenius norm using 1024 1024 1024 gradient samples including the 2 m superscript 2 𝑚 2^{m} corners and otherwise sampled randomly.

test problem	dim.	$L$		${\sigma }_{\text{max}}\left(𝐋\right)$		${\sigma }_{\text{min}}\left(𝐋\right)$		$\mathrm{vol}\left(L𝒟\right)$		$\mathrm{vol}\left(𝐋𝒟\right)$
Golinski volume [12]	$6$	$3$	$.0\cdot {10}^{-1}$	$3$	$.0\cdot {10}^{-1}$	$7$	$.3\cdot {10}^{-4}$	$4$	$.3\cdot {10}^{-2}$	$8$	$.4\cdot {10}^{-9}$
OTL circult [3]	$6$	$5$	$.4\cdot {10}^{-1}$	$5$	$.8\cdot {10}^{-1}$	$1$	$.4\cdot {10}^{-3}$	$1$	$.6\cdot {10}^0$	$1$	$.4\cdot {10}^{-5}$
piston [17]	$7$	$1$	$.0\cdot {10}^0$	$1$	$.0\cdot {10}^0$	$2$	$.6\cdot {10}^{-3}$	$1$	$.3\cdot {10}^2$	$2$	$.0\cdot {10}^{-6}$
borehole [13]	$8$	$9$	$.0\cdot {10}^{-1}$	$9$	$.2\cdot {10}^{-1}$	$4$	$.6\cdot {10}^{-6}$	$1$	$.1\cdot {10}^2$	$5$	$.0\cdot {10}^{-12}$
wing weight [10]	$10$	$3$	$.8\cdot {10}^{-1}$	$4$	$.1\cdot {10}^{-1}$	$4$	$.3\cdot {10}^{-3}$	$5$	$.9\cdot {10}^{-2}$	$2$	$.1\cdot {10}^{-10}$

Table 2. Table 2: Estimated maximum pointwise uncertainty over the domain using a 100-point minimax design with either isotropic or Lipschitz matrix metric. The maximum uncertainty is computed by extensively sampling the domain and measuring the diameter the Lipschitz constant or Lipschitz matrix uncertainty set. To aid comparison, each function’s range has been normalized to [ 0 , 1 ] 0 1 [0,1] .

test problem	dim.	isotropic Lip. const. uncertainty		isotropic Lip. matrix uncertainty		Lipschitz Lip. matrix uncertainty		ratio
Golinski volume [12]	$6$	$0$	$.69$	$0$	$.42$	$0$	$.15$	$4$	$.72$
OTL circult [3]	$6$	$1$	$.32$	$0$	$.61$	$0$	$.44$	$3$	$.01$
piston [17]	$7$	$3$	$.04$	$1$	$.26$	$0$	$.62$	$4$	$.90$
borehole [13]	$8$	$3$	$.12$	$1$	$.20$	$0$	$.49$	$6$	$.35$
wing weight [10]	$10$	$1$	$.72$	$0$	$.82$	$0$	$.39$	$4$	$.41$

Table 3. Table 3: Computational noise can lead to inaccurate estimates of the Lipschitz matrix. Using the ϵ italic-ϵ \epsilon -Lipschitz matrix instead, we can recover an accurate Lipschitz matrix. These Lipschitz matrices correspond to the partial trace example described in Section   8.1 .

Exact data	Noisy data	Noisy data with $ϵ$-Lipschitz
${\lambda }_1=3.50,{\lambda }_2=0.233$	${\lambda }_1=31.7,{\lambda }_2=31.4$	${\lambda }_1=3.48,{\lambda }_2=0.151$

Equations

$$\mathcal{L}(\mathcal{D},L):=\{f:\mathcal{D}\to \mathbb{R}:|f({\mathbf{x}}_1)-f({\mathbf{x}}_2)|\le L\Vert {\mathbf{x}}_1-{\mathbf{x}}_2{\Vert }_2,{\mathbf{x}}_1,{\mathbf{x}}_2\in \mathcal{D}\}.$$

$$\mathcal{L}(\mathcal{D},\mathbf{L}):=\{f:\mathcal{D}\to \mathbb{R}:|f({\mathbf{x}}_1)-f({\mathbf{x}}_2)|\le \Vert \mathbf{L}({\mathbf{x}}_1-{\mathbf{x}}_2){\Vert }_2,{\mathbf{x}}_1,{\mathbf{x}}_2\in \mathcal{D}\}.$$

$${\mathcal{L}}_{ϵ}(\mathcal{D},\mathbf{L}):=\{f:\mathcal{D}\to \mathbb{R}:|f({\mathbf{x}}_1)-f({\mathbf{x}}_2)|\le \Vert \mathbf{L}({\mathbf{x}}_1-{\mathbf{x}}_2){\Vert }_2+ϵ,{\mathbf{x}}_1,{\mathbf{x}}_2\in \mathcal{D}\}.$$

$$|f({\mathbf{x}}_1)-f({\mathbf{x}}_2)|\le L\Vert {\mathbf{x}}_1-{\mathbf{x}}_2{\Vert }_2=\Vert (L\mathbf{I})({\mathbf{x}}_1-{\mathbf{x}}_2){\Vert }_2.$$

$$|f({\mathbf{x}}_1)-f({\mathbf{x}}_2)|\le \Vert \mathbf{L}({\mathbf{x}}_1-{\mathbf{x}}_2){\Vert }_2\le \Vert \mathbf{L}{\Vert }_2\Vert {\mathbf{x}}_1-{\mathbf{x}}_2{\Vert }_2.$$

$$\underset{L\in {\mathbb{R}}_{+}}{\min }L\text{subject to}|f({\mathbf{x}}_1)-f({\mathbf{x}}_2)|\le L\Vert {\mathbf{x}}_1-{\mathbf{x}}_2{\Vert }_2{\mathbf{x}}_1,{\mathbf{x}}_2\in \mathcal{D}.$$

$$\Vert \mathbf{L}({\mathbf{x}}_1-{\mathbf{x}}_2){\Vert }_2=\Vert \mathbf{Q}\mathbf{P}({\mathbf{x}}_1-{\mathbf{x}}_2){\Vert }_2=\Vert \mathbf{P}({\mathbf{x}}_1-{\mathbf{x}}_2){\Vert }_2.$$

$$|f(\mathbf{x}+\delta \mathbf{p})-f(\mathbf{x}){|}^2\le \Vert \mathbf{L}(\mathbf{x}+\delta \mathbf{p}-\mathbf{x}){\Vert }_2^2={\delta }^2{\mathbf{p}}^{\top }{\mathbf{L}}^{\top }\mathbf{L}\mathbf{p}.$$

$$({\mathbf{p}}^{\top }\nabla f(\mathbf{x}){)}^2={\mathbf{p}}^{\top }\nabla f(\mathbf{x})\nabla f(\mathbf{x}{)}^{\top }\mathbf{p}\le {\mathbf{p}}^{\top }{\mathbf{L}}^{\top }\mathbf{L}\mathbf{p}$$

$$\nabla f(\mathbf{x})\nabla f(\mathbf{x}{)}^{\top }⪯{\mathbf{L}}^{\top }\mathbf{L}.$$

$$f(\mathbf{x})=g({\mathbf{U}}^{\top }\mathbf{x})\text{where}g:{\mathbb{R}}^n\to \mathbb{R},\mathbf{U}\in {\mathbb{R}}^{m\times n},m>n.$$

$$|f({\mathbf{x}}_1)-f({\mathbf{x}}_2)|\le \Vert \mathbf{L}({\mathbf{x}}_1-{\mathbf{x}}_2){\Vert }_2=0$$

$$\begin{array}{rl}|f({\mathbf{x}}_1)-f({\mathbf{x}}_2)|& =|g({\mathbf{U}}^{\top }{\mathbf{x}}_1)-g({\mathbf{U}}^{\top }{\mathbf{x}}_2)|\\ & \le \Vert {\mathbf{L}}_g[({\mathbf{U}}^{\top }{\mathbf{x}}_1)-({\mathbf{U}}^{\top }{\mathbf{x}}_2)]{\Vert }_2=\Vert \mathbf{U}{\mathbf{L}}_g{\mathbf{U}}^{\top }({\mathbf{x}}_1-{\mathbf{x}}_2){\Vert }_2.\end{array}$$

$$\mathbf{C}:=\mathbb{E}[\nabla f\nabla f^{\top }]={\int }_{\mathbf{x}\in \mathcal{D}}\nabla f(\mathbf{x})\nabla f(\mathbf{x}{)}^{\top }\mathrm{d}\mu (\mathbf{x})\in {\mathbb{R}}^{m\times m}.$$

$$\mathbf{C}={\int }_{\mathbf{x}\in \mathcal{D}}\nabla f(\mathbf{x})\nabla f(\mathbf{x}{)}^{\top }\mathrm{d}\mu (\mathbf{x})⪯{\mathbf{L}}^{\top }\mathbf{L}.$$

$$d_{\mathbf{L}}(\mathbf{x},\mathbf{y}):=\Vert \mathbf{L}(\mathbf{x}-\mathbf{y}){\Vert }_2\mathbf{x},\mathbf{y}\in \mathcal{D}.$$

$$N_{ϵ}(\mathcal{D};\mathbf{L}):=\underset{M,\{{\mathbf{x}}_j{\}}_{j=1}^M\subset \mathcal{D}}{\mathrm{argmin}}M\text{such that}\mathcal{D}\subseteq \bigcup _{j=1}^M{\mathcal{B}}_{ϵ}({\mathbf{x}}_j;\mathbf{L}),$$

$$\mathcal{X}(\mathcal{D},M,\mathbf{L}):=\underset{\{{\mathbf{x}}_j{\}}_{j=1}^M\subset \mathcal{D}}{\mathrm{argmin}}\underset{\mathbf{x}\in \mathcal{D}}{\max }\underset{j=1,\dots ,M}{\min }\Vert \mathbf{L}(\mathbf{x}-{\mathbf{x}}_j){\Vert }_2;$$

$${\left(\frac{1}{ϵ}\right)}^m\frac{\mathrm{vol}(\mathcal{D})}{\mathrm{vol}({\mathcal{B}}_1(\cdot ;\mathbf{L}))}\le N_{ϵ}(\mathcal{D};\mathbf{L})\le {\left(\frac{3}{ϵ}\right)}^m\frac{\mathrm{vol}(\mathcal{D})}{\mathrm{vol}({\mathcal{B}}_1(\cdot ;\mathbf{L}))},$$

$$\mathbf{L}\mathcal{D}:=\{\mathbf{L}\mathbf{x}:\mathbf{x}\in \mathcal{D}\}$$

$$N_{ϵ}(\mathcal{D};\mathbf{L})=N_{ϵ}(\mathbf{L}\mathcal{D};\mathbf{I}).$$

$$\mathcal{L}(\mathcal{D},\mathbf{L},{\mathcal{P}}_M):=\{f\in \mathcal{L}(\mathcal{D},\mathbf{L}):f({\mathbf{x}}_j)=y_j\{{\mathbf{x}}_j,y_j\}\in {\mathcal{P}}_M\}.$$

$$\mathcal{U}(\mathbf{x};\mathbf{L},{\mathcal{P}}_M):=\{f(\mathbf{x}):f\in \mathcal{L}(\mathcal{D},\mathbf{L},{\mathcal{P}}_M)\}.$$

$$y_j-\Vert \mathbf{L}(\mathbf{x}-{\mathbf{x}}_j){\Vert }_2\le f(\mathbf{x})\le y_j+\Vert \mathbf{L}(\mathbf{x}-{\mathbf{x}}_j){\Vert }_2.$$

$$f_{\min }(\mathbf{x};\mathbf{L},{\mathcal{P}}_M)$$

$$f_{\max }(\mathbf{x};\mathbf{L},{\mathcal{P}}_M)$$

$$\mathcal{U}(\mathbf{x};\mathbf{L},{\mathcal{P}}_M)=\left[f_{\min }(\mathbf{x};\mathbf{L},{\mathcal{P}}_M),f_{\max }(\mathbf{x};\mathbf{L},{\mathcal{P}}_M)\right].$$

$$f_{\mathrm{mid}}(\mathbf{x};\mathbf{L},{\mathcal{P}}_M):=\frac{1}{2}\left(f_{\min }(\mathbf{x})+f_{\max }(\mathbf{x})\right).$$

$$\underset{f\in \mathcal{L}(\mathcal{D},\mathbf{L},{\mathcal{P}}_M)}{\sup }|f(\mathbf{x})-f_{\mathrm{mid}}(\mathbf{x};\mathbf{L},{\mathcal{P}}_M)|\le \underset{f\in \mathcal{L}(\mathcal{D},\mathbf{L},{\mathcal{P}}_M)}{\sup }|f(\mathbf{x})-\tilde{f}(\mathbf{x})|.$$

$$\underset{f\in \mathcal{L}(\mathcal{D},\mathbf{L})}{\sup }\underset{\tilde{f}:\mathcal{D}\to \mathbb{R}}{\inf }\underset{\mathbf{x}\in \mathcal{D}}{\max }|f(\mathbf{x})-\tilde{f}(\mathbf{x})|\le ϵ\text{where}{\mathcal{P}}_M=\{{\mathbf{x}}_j,f({\mathbf{x}}_j){\}}_{j=1}^M$$

$$\underset{f\in \mathcal{L}(\mathcal{D},\mathbf{L})}{\sup }\underset{\tilde{f}:\mathcal{D}\to \mathbb{R}}{\inf }\underset{\mathbf{x}\in \mathcal{D}}{\max }|f(\mathbf{x})-\tilde{f}(\mathbf{x})|=\underset{f\in \mathcal{L}(\mathcal{D},\mathbf{L})}{\sup }\underset{\mathbf{x}\in \mathcal{D}}{\max }|f(\mathbf{x})-f_{\mathrm{mid}}(\mathbf{x};\mathbf{L},{\mathcal{P}}_M)|.$$

Full text

{AMS}

26B35, 62K05, 68Q25 {DOI}

1 Introduction

With the increasing sophistication of computer models, practitioners in science and engineering often confront the curse of dimensionality—the phenomena that for many relevant computational tasks, obtaining the desired solution requires a computational cost that grows exponentially with the number of input parameters [8, 34]. To mitigate this curse, practitioners may introduce a lower dimensional reparameterization of the model inputs yielding a similar output to the original. When applicable, this parameter reduction enables otherwise high-dimensional computations to exploit low-dimensional constructed parameters. Standard engineering practice uses a global sensitivity analysis [27] to identify a subset of input parameters that have relatively little affect on the output. Fixing these insignificant parameters at a nominal value provides a lower-dimensional parameterization using the remaining parameters [29]. Another approach reparameterizes the model using a few linear combinations of the input variables defining an active subspace along which the model varies [4]. Active subspaces include subset-based approaches since any subset can be encoded as the span of columns of the identity matrix. More generally, a nonlinear reparameterization of the input-output map can be used [15]. Here we introduce the Lipschitz matrix—a generalization of the scalar Lipschitz constant. The Lipschitz matrix can identify an active subspace, motivates a design of experiments, defines uncertainty away from function evaluations, and yields improved error bounds that can mitigate the curse of dimensionality.

1.1 Definition

We define the Lipschitz matrix analogously to the Lipschitz constant. Given a domain ${\mathcal{D}}\subset\mathbb{R}^{m}$, the scalar Lipschitz constant $L\in\mathbb{R}_{+}$ defines a set of scalar Lipschitz functions denoted ${\mathcal{L}}({\mathcal{D}},L)$ that satisfy

[TABLE]

The Lipschitz matrix changes this definition, moving the constant $L$ inside the norm and promoting it to a matrix ${\mathbf{L}}\in\mathbb{R}^{m\times m}$. This defines matrix Lipschitz functions

[TABLE]

In our notation, the type of the second argument of ${\mathcal{L}}(\cdot,\cdot)$ indicates whether this set refers to the scalar Eq. 1 or matrix Eq. 2 case. Additionally, we define $\epsilon$-Lipschitz functions that are near by Lipschitz functions:

[TABLE]

This function class is useful when analyzing functions with computational noise [20] or irrelevant oscillations; see examples in Section 8.

1.2 Equivalence

We refer to both scalar Eq. 1 and matrix Eq. 2 Lipschitz functions as simply Lipschitz functions as these two sets are nested. All scalar Lipschitz functions $f\in{\mathcal{L}}({\mathcal{D}},L)$ are also matrix Lipschitz functions with $f\in{\mathcal{L}}({\mathcal{D}},L{\mathbf{I}})$:

[TABLE]

Hence ${\mathcal{L}}({\mathcal{D}},L)={\mathcal{L}}({\mathcal{D}},L{\mathbf{I}})$. Similarly, all matrix Lipschitz functions $f\in{\mathcal{L}}({\mathcal{D}},{\mathbf{L}})$ are also scalar Lipschitz functions with $f\in{\mathcal{L}}({\mathcal{D}},\|{\mathbf{L}}\|_{2})$ as

[TABLE]

Hence ${\mathcal{L}}({\mathcal{D}},{\mathbf{L}})\subseteq{\mathcal{L}}({\mathcal{D}},\|{\mathbf{L}}\|_{2})$. This nesting also holds for $\epsilon$-Lipschitz functions.

1.3 Smallest Lipschitz Matrix

Given a particular function, we seek the smallest possible Lipschitz matrix to tighten our results. The challenge using the Lipschitz matrix is there is no natural ordering unlike the scalar Lipschitz constant. Hence we must impose an ordering.

For the scalar Lipschitz constant, the ordering of the real line provides a natural ordering for Lipschitz constants. The smallest Lipschitz constant for a function $f$ is

[TABLE]

Unlike scalars, matrices have no natural ordering. However by invoking the polar decomposition [16, Thm. 7.3.1] we can define a partial order for Lipschitz matrices. Any matrix ${\mathbf{L}}\in\mathbb{R}^{m\times m}$ has a polar decomposition into the product ${\mathbf{L}}={\mathbf{Q}}{\mathbf{P}}$ where ${\mathbf{Q}}\in\mathbb{R}^{m\times m}$ has orthonormal columns and ${\mathbf{P}}\in\mathbb{R}^{m\times m}$ is a symmetric positive semidefinite matrix, denoted ${\mathbf{P}}\in\mathbb{S}_{+}^{m\times m}$. As the 2-norm is unitarily invariant,

[TABLE]

Thus without loss of generality, we can assume a Lipschitz matrix ${\mathbf{L}}$ is symmetric positive semidefinite. Positive semidefinite matrices have a natural partial ordering: the Loewner partial order [16, Def. 7.7.1] where for ${\mathbf{A}},{\mathbf{B}}\in\mathbb{S}_{+}^{m\times m}$, we write ${\mathbf{A}}\preceq{\mathbf{B}}$ if ${\mathbf{B}}-{\mathbf{A}}$ is positive semidefinite. This is only a partial order as when ${\mathbf{B}}-{\mathbf{A}}$ is indefinite, we cannot order ${\mathbf{A}}$ and ${\mathbf{B}}$.

To define the smallest Lipschitz matrix, we choose a total order compatible with the partial order. Two convenient choices are the trace and determinant of ${\mathbf{L}}$ [16, Cor. 7.7.4d,e]. In the majority of our results, we use the squared Frobenius norm of ${\mathbf{L}}$, i.e., the squared sum of the eigenvalues of ${\mathbf{L}}$ which obeys the partial order by [16, Cor. 7.7.4c], as this yields an amenable optimization problem as discussed in Section 3.

1.4 Derivatives

The Lipschitz matrix bounds the derivatives of $f$. Suppose $f$ is differentiable at ${\mathbf{x}}$ in the interior of ${\mathcal{D}}$. Then for any ${\mathbf{p}}\in\mathbb{R}^{m}$ there is some $\delta>0$ such that ${\mathbf{x}}+\delta{\mathbf{p}}\in{\mathcal{D}}$. From the Lipschitz matrix constraint Eq. 2,

[TABLE]

Dividing by $\delta^{2}$ and taking the limit as $\delta\to 0$, we bound the gradient of $f$, $\nabla f$, by

[TABLE]

As this holds for all ${\mathbf{p}}$, we write this compactly using the Loewner partial order:

[TABLE]

If $f$ is not differentiable, a similar result holds for any Gateaux derivative. Note, due to the presence of $\epsilon>0$ in the definition of $\epsilon$-Lipschitz functions Eq. 3, the derivatives of $f$ do not provide a lower bound on the $\epsilon$-Lipschitz matrix.

1.5 Connection to ridge functions

A ridge function [23] depends only on a few linear combinations of its input variables; i.e., a function $f:\mathbb{R}^{m}\to\mathbb{R}$ is a ridge function if there is a function $g$ such that

[TABLE]

We call $g$ the ridge profile, $n$ the ridge dimension, and the range of ${\mathbf{U}}$ the active subspace. The Lipschitz matrix is intimately connected to ridge functions: informally, $f$ is a ridge function if and only if it has a low-rank Lipschitz matrix. The following theorem makes this precise.

Theorem 1.1**.**

*Suppose $f:{\mathcal{D}}\to\mathbb{R}$ is a Lipschitz function with a Lipschitz constant $L$. Then $f$ is a ridge function with ridge dimension $n$ if and only if there exists a ${\mathbf{L}}\in\mathbb{R}^{m\times m}$ with rank $n$ such that $f\in{\mathcal{L}}({\mathcal{D}},{\mathbf{L}})$. *

Proof 1.2**.**

Suppose $f\in{\mathcal{L}}({\mathcal{D}},{\mathbf{L}})$ where ${\mathbf{L}}$ is rank $n$. Let ${\mathbf{U}}\in\mathbb{R}^{m\times n}$ be an orthonormal basis for the range of ${\mathbf{L}}$. If ${\mathbf{x}}_{1},{\mathbf{x}}_{2}\in{\mathcal{D}}$ with ${\mathbf{U}}^{\top}({\mathbf{x}}_{1}-{\mathbf{x}}_{2})=0$, then

[TABLE]

as ${\mathbf{x}}_{1}-{\mathbf{x}}_{2}$ is in the nullspace of ${\mathbf{L}}$. Hence $f$ is constant in all directions in the nullspace of ${\mathbf{U}}$ and thus there is a $g:\mathbb{R}^{n}\to\mathbb{R}$ such that $f({\mathbf{x}})=g({\mathbf{U}}^{\top}{\mathbf{x}})$.

Suppose $f({\mathbf{x}})=g({\mathbf{U}}^{\top}{\mathbf{x}})$ where ${\mathbf{U}}\in\mathbb{R}^{m\times n}$ has orthonormal columns. As $f$ is Lipschitz, so too is $g$. Let ${\mathbf{L}}_{g}\in\mathbb{S}_{+}^{n\times n}$ be a full-rank Lipschitz matrix for $g$, then

[TABLE]

*Hence $f\in{\mathcal{L}}({\mathcal{D}},{\mathbf{U}}{\mathbf{L}}_{g}{\mathbf{U}}^{\top})$ and ${\mathbf{U}}{\mathbf{L}}_{g}{\mathbf{U}}^{\top}$ is rank $n$. *

Thus if we identify a low-rank Lipschitz matrix for a function, the range of the Lipschitz matrix defines an active subspace. However in our experience exact ridge functions are rare. More frequently, a function will be approximately a ridge function; we informally call these functions ridge-like. In the context of Lipschitz matrices, we say $f$ is ridge-like if it has an approximately low-rank Lipschitz matrix; then the dominant eigenspace defines an active subspace (see Theorem 4.1). Other approaches for identifying and approximating ridge-like functions include: polynomial ridge approximation [5, 14], sufficient dimension reduction techniques from statistical regression [11], and the mean gradient outer-product (MeGO) [4].

1.6 Comparison to MeGO

Given a probability measure $\mu$, the mean gradient outer-product (MeGO) [4, eq. (3.2)] is

[TABLE]

The dominant eigenspaces of ${\mathbf{C}}$ provide one way to identify the active subspace of $f$.

When $f$ is differentiable almost everywhere on its domain, we can bound the MeGO matrix by the Lipschitz matrix. From Eq. 10, the Lipschitz matrix bounds the gradient: $\nabla f({\mathbf{x}})\nabla f({\mathbf{x}})^{\top}\preceq{\mathbf{L}}^{\top}{\mathbf{L}}$. Then as $\mu$ is a probability measure

[TABLE]

This bound is tight when $f$ is a linear function $f({\mathbf{x}})={\mathbf{a}}^{\top}{\mathbf{x}}$ in which case ${\mathbf{C}}={\mathbf{a}}{\mathbf{a}}^{\top}$, ${\mathbf{L}}=\|{\mathbf{a}}\|_{2}^{-1}{\mathbf{a}}{\mathbf{a}}^{\top}$, and ${\mathbf{L}}^{\top}{\mathbf{L}}={\mathbf{a}}{\mathbf{a}}^{\top}$. Like the Lipschitz matrix, dominant eigenspace of ${\mathbf{C}}$ identifies an active subspace [4, Thm. 1]. As a corollary, if $f$ is a ridge function, then the nullspaces of ${\mathbf{C}}$ and ${\mathbf{L}}$ are the same.

1.7 Applications of the Lipschitz Matrix

The Lipschitz matrix provides both analytical tools and practical results. As we show in Section 2, replacing the Lipschitz constant with the Lipschitz matrix allows us to approximate, integrate, and optimize functions using fewer evaluations. If we can identify a low-rank Lipschitz matrix, then the cost of these tasks no longer scales with the number of parameters, but instead the rank of the Lipschitz matrix. In practice, there are few functions for which we can compute the Lipschitz matrix exactly. We show in Section 3 that we can estimate the Lipschitz matrix from arbitrary combinations of function evaluations $f({\mathbf{x}}_{i})$ and gradients $\nabla f({\mathbf{x}}_{k}^{\prime})$ by solving a semidefinite program. We can use the estimated Lipschitz matrix to then identify an active subspace (Section 4), guide the design of experiments (Section 5), construct ridge approximations for dimension reduction (Section 6), and quantify uncertainty (Section 7).

1.8 Reproducibility

Following the principles of reproducible research, we provide code implementing the algorithms described in this paper and scripts generating the data appearing in the figures and tables available at http://github.com/jeffrey-hokanson/PSDR/.

2 Algorithm Complexity

In this section we use results from information-based complexity [34, 35] to bound the worst-case optimal cost of approximation, integration, and optimization for Lipschitz matrix functions. These results parallel similar results for scalar Lipschitz functions [31] with one important distinction. For scalar Lipschitz functions on ${\mathcal{D}}\subset\mathbb{R}^{m}$, each of these three tasks requires $\mathcal{O}(\epsilon^{-m})$ function evaluations to obtain $\epsilon>0$ accuracy—this exponential growth in dimension is the curse of dimensionality [8]. For functions with a rank-$r$ Lipschitz matrix these tasks require only $\mathcal{O}(\epsilon^{-r})$ function evaluations. If $r<m$, then complexity is independent of dimension and we have mitigated the curse of dimensionality. The key ingredient is using the Lipschitz matrix to provide a (pseudo-)metric on the domain

[TABLE]

Using this metric, we show the complexity of approximation, integration, and optimization is proportional to $\epsilon$-internal covering number of ${\mathcal{D}}$:

[TABLE]

where ${\mathcal{B}}_{\epsilon}({\mathbf{x}}_{j};{\mathbf{L}})\coloneqq\{{\mathbf{x}}:\|{\mathbf{L}}({\mathbf{x}}-{\mathbf{x}}_{j})\|_{2}\leq\epsilon\}$ is the $\epsilon$-ball in $d_{{\mathbf{L}}}$. These points $\{{\mathbf{x}}_{j}\}_{j=1}^{M}$ have a special interpretation: they are an $M$-point minimax optimal design on ${\mathcal{D}}$:

[TABLE]

see proof in [6, Thm. 4.7]. By replacing the Lipschitz constant by the Lipschitz matrix fewer points are required cover the domain as illustrated in Fig. 1. The same can be seen in bounds for the covering number.

\DTLsettabseparator\DTLloaddb[noheader=false]coverscalardata/fig_cover_scalar.dat \DTLloaddb[noheader=false]covermatrixdata/fig_cover_matrix.dat \DTLmaketabspace

Theorem 2.1** ([36, Thm. 14.2]).**

If ${\mathcal{D}}\subset\mathbb{R}^{m}$ is convex and ${\mathbf{L}}$ is full rank

[TABLE]

*where $\operatorname{{\mathrm{v}ol}}$ denotes the Lebesgue measure in $\mathbb{R}^{m}$ and ${\mathcal{B}}_{1}(\cdot;{\mathbf{L}})$ the unit ball in $d_{{\mathbf{L}}}$. *

This theorem implies $N_{\epsilon}({\mathcal{D}};{\mathbf{L}})=\mathcal{O}(\epsilon^{-m})$ when ${\mathbf{L}}$ is full rank. If ${\mathbf{L}}$ is low-rank we can apply this theorem by interpreting the Lipschitz matrix as warping the domain

[TABLE]

and equipping this domain with the standard Euclidean metric $d({\mathbf{x}},{\mathbf{y}})=\|{\mathbf{x}}-{\mathbf{y}}\|_{2}$. Then covering number is

[TABLE]

When ${\mathbf{L}}$ is rank-$r$, then ${\mathbf{L}}{\mathcal{D}}$ is an $r$-dimensional subset of $\mathbb{R}^{m}$. Applying Theorem 2.1 on this $r$-dimensional problem, note $N_{\epsilon}({\mathbf{L}}{\mathcal{D}};{\mathbf{I}})=\mathcal{O}(\epsilon^{-r})$.

In the remainder of this section we first establish a definition of uncertainty for Lipschitz functions and then use this definition to derive complexity results for approximation, integration, and optimization. Then we discuss how the Lipschitz matrix yields both asymptotic and non-asymptotic reductions to the complexity of these tasks when compared to the Lipschitz constant.

2.1 Uncertainty

One use of the Lipschitz matrix is to provide constraints on what values a function can take away from points where know its values. We denote these point queries consisting of a point ${\mathbf{x}}_{j}\in{\mathcal{D}}$ and a response $y_{j}\in\mathbb{R}$ by ${\mathcal{P}}_{M}=\{{\mathbf{x}}_{j},y_{j}\}_{j=1}^{M}$. This allows us to define the space of all Lipschitz functions with the Lipschitz matrix ${\mathbf{L}}$ that interpolate these point queries:

[TABLE]

Using this notation we can then define the uncertainty set at a point ${\mathbf{x}}\in{\mathcal{D}}$

[TABLE]

This set is actually an interval for each ${\mathbf{x}}$. From the definition of Lipschitz matrix continuity Eq. 2, for any $f\in{\mathcal{L}}({\mathcal{D}},{\mathbf{L}},{\mathcal{P}}_{M})$ we have bounds

[TABLE]

As these bounds apply for each $\{{\mathbf{x}}_{j}\}_{j=1}^{M}$ these provide lower and upper bounds for each point ${\mathbf{x}}$:

[TABLE]

Since these are the only constraints on $f({\mathbf{x}})$, the uncertainty set is an interval

[TABLE]

An important tool in our results is the central approximation

[TABLE]

This function minimizes the worst-case pointwise error of all Lipschitz approximations that interpolate the point set.

Lemma 2.2**.**

Given a domain ${\mathcal{D}}\subset\mathbb{R}^{m}$, Lipschitz matrix ${\mathbf{L}}\in\mathbb{R}^{m\times m}$, and point queries ${\mathcal{P}}_{M}=\{{\mathbf{x}}_{j},y_{j}\}_{j=1}^{M}$ then for any fixed $\widetilde{f}:{\mathcal{D}}\to\mathbb{R}$

[TABLE]

Proof 2.3**.**

*This follows immediately from the uncertainty interval Eq. 27. Any other choice for $\widetilde{f}({\mathbf{x}})$ would have an error greater than $\frac{1}{2}[f_{\max}({\mathbf{x}})-f_{\min}({\mathbf{x}})]$ for some choice of $f$; cf., [34, pp. 12-13]. *

2.2 Complexity Results

The following results show the number of function queries are necessary to perform approximation, integration, and optimization to within a tolerance $\epsilon>0$ in the worst case depends on the $\epsilon$-covering of the domain.

2.2.1 Approximation

The following theorem shows that by querying a function at a minimax optimal design, the resulting central approximation yields the best approximation in the worst case.

Theorem 2.4**.**

Given a domain ${\mathcal{D}}\subset\mathbb{R}^{m}$ and a Lipschitz matrix ${\mathbf{L}}\in\mathbb{R}^{m\times m}$, the minimum number of point queries $M$ such that

[TABLE]

*is the $\epsilon$ internal covering number, $M=N_{\epsilon}({\mathcal{D}};{\mathbf{L}})$. The optimal point queries correspond to the $M$-point minimax optimal design, ${\mathcal{X}}({\mathcal{D}},M,{\mathbf{L}})$ and the optimal approximation $\widetilde{f}$ is the central approximation $f_{\mathrm{mid}}$. *

Proof 2.5**.**

Using Lemma 2.2 we can remove the middle optimization in Eq. 30 by replacing $\widetilde{f}$ with the central approximation $f_{\mathrm{mid}}$:

[TABLE]

The greatest uncertainty comes when $f({\mathbf{x}}_{j})=0$ for $j=1,\ldots,M$ and hence

[TABLE]

*Then the minimum $M$ to make the right hand side above less than $\epsilon$ is $N_{\epsilon}({\mathcal{D}},{\mathbf{L}})$ and the optimal queries correspond $N_{\epsilon}({\mathcal{D}},{\mathbf{L}})$-point minimax optimal design. *

2.2.2 Optimization

The argument for the complexity of optimization is essentially that of approximation: unless we can approximate $f$ within $\epsilon$, we cannot globally optimize $f$ within $\epsilon$.

Theorem 2.6**.**

*Suppose ${\mathcal{D}}\subset\mathbb{R}^{m}$ is compact and ${\mathbf{L}}\in\mathbb{R}^{m\times m}$. In the worst case, the minimum number of samples to find the maximum $f\in{\mathcal{L}}({\mathcal{D}},{\mathbf{L}})$ to within $\epsilon$ is the $\epsilon$-internal covering number $N_{\epsilon}({\mathbf{L}}{\mathcal{D}})$. *

Proof 2.7**.**

As ${\mathcal{D}}$ is compact and $f$ is Lipschitz, $f$ must have a finite maximizer $f^{\star}$ where $f^{\star}=f({\mathbf{x}}^{\star})$ for at least one ${\mathbf{x}}^{\star}\in{\mathcal{D}}$. We first establish an upper bound on the number of point queries required. Using $M=N_{\epsilon}({\mathcal{D}},{\mathbf{L}})$ point queries ${\mathcal{P}}_{M}$ in a minimax optimal design there is a ${\mathbf{x}}_{j}$ in ${\mathcal{P}}_{M}$ such that $\|{\mathbf{L}}({\mathbf{x}}^{\star}-{\mathbf{x}}_{j})\|_{2}\leq\epsilon$. By Eq. 24

[TABLE]

*and hence $f({\mathbf{x}}_{j})\geq f^{\star}-\epsilon$. To show this upper bound is obtained, consider the function $f\in{\mathcal{L}}({\mathcal{D}},{\mathbf{L}})$ where $f({\mathbf{x}})\geq 0$ for all ${\mathbf{x}}\in{\mathcal{D}}$, $f$ has optimizer $f({\mathbf{x}}^{\star})=2\epsilon=f^{\star}$, and $f$ has minimum integrand. If we choose $M<N_{\epsilon}({\mathcal{D}},{\mathbf{L}})$ point queries we can adversarially choose ${\mathbf{x}}^{\star}$ such that $\|{\mathbf{L}}({\mathbf{x}}-{\mathbf{x}}^{\star})\|_{2}>\epsilon$ and consequently $f({\mathbf{x}}_{j})<\epsilon$ for all $j=1,\ldots,M$. *

2.2.3 Integration

This result parallels the one-dimensional Lipschitz case presented by Traub and Werschulz [34, chap. 2].

Theorem 2.8**.**

Suppose ${\mathcal{D}}\subset\mathbb{R}^{m}$ be compact and ${\mathbf{L}}\in\mathbb{R}^{m\times m}$. Let $\phi^{\star}$ denote the quadrature rule from integrating the central approximation on an $M$-point minimax optimal design

[TABLE]

Let $\phi$ denote any other $M$ point quadrature rule integrating an approximation $\widetilde{f}$, then

[TABLE]

Proof 2.9**.**

From Lemma 2.2, the central approximation has the smallest pointwise worst case error of any approximation interpolating the point queries ${\mathcal{P}}_{M}$. Thus

[TABLE]

*Then as the points ${\mathcal{X}}({\mathcal{D}},M,{\mathbf{L}})$ minimize the worse case error in the central approximation, $\phi^{\star}$ is the worst-case optimal quadrature rule. *

Corollary 2.10**.**

In the setting of Theorem 2.8,

[TABLE]

Proof 2.11**.**

Let $M=N_{\epsilon}({\mathcal{D}},{\mathbf{L}})$. From the definition of $\phi^{\star}$ Eq. 33,

[TABLE]

*By Theorem 2.4, $|f_{\mathrm{mid}}({\mathbf{x}})-f({\mathbf{x}})|\leq\epsilon$; integrating this bound yields the result. *

2.3 Asymptotic Improvements

The Lipschitz matrix provides two asymptotic improvements over the Lipschitz constant. If the Lipschitz matrix is rank-$r$, then complexity of approximation, integration, and optimization grows like $\mathcal{O}(\epsilon^{-r})$—not $\mathcal{O}(\epsilon^{-m})$. This also applies to ridge functions with ridge dimension $r$ by Theorem 1.1. Even if the Lipschitz matrix is full rank, the Lipschitz matrix can still substantially reduce the covering number compared to the Lipschitz constant. Interpreting the Lipschitz matrix as transforming the domain, from Theorem 2.1 we have the bound

[TABLE]

For the Lipschitz matrix and Lipschitz constant

[TABLE]

Hence if $|\det{\mathbf{L}}|\ll L^{m}$, we have substantially reduced the covering number. This is not uncommon: for the test problems shown in Table 1, $|\det{\mathbf{L}}|$ is multiple orders of magnitude smaller than $L^{m}$.

2.4 Non-asymptotic Improvement

Before the asymptotic limit, the Lipschitz matrix can slow the growth of covering number. Denoting the singular values of ${\mathbf{L}}$ in non-increasing order as $\sigma_{j}({\mathbf{L}})$, when $\epsilon\gg\sigma_{j}({\mathbf{L}})$ then ${\mathbf{L}}{\mathcal{D}}$ is effectively $j$-dimensional as a single $\epsilon$-ball can cover the dimensions $j+1,\ldots,m$. This temporally slows the growth of the covering number as illustrated in Fig. 2.

3 Computing a Lipschitz Matrix

For a particular function $f$, our goal is to identify the smallest Lipschitz matrix ${\mathbf{L}}$ such that $f\in{\mathcal{L}}({\mathcal{D}},{\mathbf{L}})$. As discussed in Section 1.3, we must choose an ordering of positive semidefinite matrices. Here we minimize the Frobenius norm of the Lipschitz matrix yielding the program

[TABLE]

Ideally given $f$, we would identify this Frobenius-norm minimizing Lipschitz matrix in closed form. For most functions, this is infeasible. Instead, we approximate this program by discretizing the constraints using a finite number of point and/or gradient queries. This yields a semidefinite program to identify the Lipschitz matrix. Although tempting to enforce a low-rank constraint on the Lipschitz matrix in this program, we show this can yield in non-informative results. For the special case of quadratic functions, we provide a finite program exactly solving the full program Eq. 40. We finally illustrate how this discretized program converges with increasing queries.

3.1 Semidefinite Program for a Lipschitz Matrix

Both point queries ${\mathcal{P}}_{M}=\{{\mathbf{x}}_{j},f({\mathbf{x}}_{j})\}_{j=1}^{M}$ and gradient queries ${\mathcal{G}}_{N}=\{\nabla f({\mathbf{x}}_{k}^{\prime})\}_{k=1}^{N}$ provide constraints on the minimum Frobenius norm Lipschitz matrix. Recalling from Eq. 10 that the outer product of gradients is bounded above by ${\mathbf{L}}^{\top}{\mathbf{L}}$, the discretized minimum Frobenius-norm Lipschitz matrix solves

[TABLE]

This formulation has difficult non-convex quadratic constraints. Instead, we reformulate Eq. 41 in terms of the squared Lipschitz matrix ${\mathbf{H}}\coloneqq{\mathbf{L}}^{\top}{\mathbf{L}}$ and instead solve

[TABLE]

This, unlike Eq. 41, is convex semidefinite program. After parametrizing the space of symmetric matrices, this program has three sets of constraints: linear inequality constraints from the point queries, semidefinite constraints from the gradient queries, and the semidefinite constraint ${\mathbf{0}}\preceq{\mathbf{H}}$. Our numerical experiments use CVXOPT [1] to solve Eq. 42 and take ${\mathbf{L}}$ to be the symmetric positive definite square root of ${\mathbf{H}}$.

3.2 Semidefinite Program for an $\epsilon$-Lipschitz Matrix

We use a similar approach working in terms of the squared Lipschitz matrix ${\mathbf{H}}={\mathbf{L}}^{\top}{\mathbf{L}}$ to compute an $\epsilon$-Lipschitz matrix. However, because $\epsilon$-Lipschitz functions are not differentiable, gradient queries ${\mathcal{G}}_{N}$ do not constrain this matrix. This leaves only point queries which by Eq. 3 must satisfy

[TABLE]

Replacing the point query constraint in Eq. 42 with this expression yields

[TABLE]

3.3 Low-rank Solutions

It is tempting to impose a rank constraint on ${\mathbf{H}}$, and consequently ${\mathbf{L}}$, to reduce complexity. Unfortunately, a low-rank constraint can yield misleading results. For example, suppose we have point queries $\{{\mathbf{x}}_{j},y_{j}\}_{j=1}^{M}={\mathcal{P}}_{M}$ where $\{{\mathbf{x}}_{j}\}_{j=1}^{M}$ is in general position. Then for almost every vector ${\mathbf{a}}\in\mathbb{R}^{m}$, the projections onto $\operatorname{Span}{\mathbf{a}}$ are distinct: ${\mathbf{a}}^{\top}{\mathbf{x}}_{i}\neq{\mathbf{a}}^{\top}{\mathbf{x}}_{j}$ for $i\neq j$. Thus, for almost every ${\mathbf{a}}$ there is a rank-1 Lipschitz matrix whose range is $\operatorname{Span}{\mathbf{a}}$:

[TABLE]

Hence regardless of the actual structure of $f$ we have likely mistakenly identified it as a one-dimensional ridge function. Gradient constraints for the Lipschitz matrix do not fail in a similar way. If $\operatorname{Range}\{\nabla f({\mathbf{x}}_{k}^{\prime})\}_{j=1}^{N}$ is $n$-dimensional, then ${\mathbf{H}}$, and consequently ${\mathbf{L}}$, must be at least rank-$n$.

3.4 Using the Determinant

In Section 2.3 we saw that the complexity of approximation, integration, and optimization are proportional to $|\det{\mathbf{L}}|$. Why not use the determinant as the objective function instead of the Frobenius norm in Eq. 40? There are two important reasons. We loose convexity: $\|{\mathbf{L}}\|_{\mathrm{F}}^{2}=\operatorname{Tr}{\mathbf{H}}$ is convex whereas $|\det{\mathbf{L}}|^{2}=\det{\mathbf{H}}$ is concave. The other is that the determinant yields uninformative Lipschitz matrices. As illustrated in the previous subsection, given finite point queries we can always find a low-rank Lipschitz matrix. As the determinant of this matrix is zero it is an optimal, but uninformative, Lipschitz matrix.

3.5 Quadratic Functions

In the case of a quadratic function we can show that there are only finite number of active gradient constraints of the continuous problem Eq. 40 yielding a finite semidefinite program for its Lipschitz matrix. Suppose $f$ is a quadratic function:

[TABLE]

${\mathbf{A}}\in\mathbb{R}^{m\times m}$, and ${\mathbf{b}}\in\mathbb{R}^{m}$ whose gradient is $\nabla f({\mathbf{x}})=({\mathbf{A}}+{\mathbf{A}}^{\!{\top}}){\mathbf{x}}+{\mathbf{b}}$. The following theorem bounds the gradient outer-product above by points on the corners of the domain: those points where each coordinate takes on the value $-1$ or $1$.

Theorem 3.1**.**

Suppose $f:[-1,1]^{m}\to\mathbb{R}$ is a quadratic function as in Eq. 46. For any ${\mathbf{x}}\in{\mathcal{D}}$ there is a point ${\mathbf{c}}$ on a corner of ${\mathcal{D}}$ such that

[TABLE]

Proof 3.2**.**

Denote ${\mathbf{A}}_{S}={\mathbf{A}}+{\mathbf{A}}^{\top}$. Consider the difference of these two gradient outer-products for any ${\mathbf{c}}$ on the corner of ${\mathcal{D}}$ and ${\mathbf{x}}\in{\mathcal{D}}$:

[TABLE]

As ${\mathbf{c}}$ is on the corner of $[-1,1]^{m}$, then ${\mathbf{c}}{\mathbf{c}}^{\top}\succeq{\mathbf{x}}{\mathbf{x}}^{\top}$ and

[TABLE]

As ${\mathbf{c}}$ is on a corner, there is a ${\mathbf{c}}$ such that the entires of the entries of ${\mathbf{c}}-{\mathbf{x}}$ can have any combination of signs $\{+,0\}$ or $\{-,0\}$. Thus we can choose a corner such that ${\mathbf{A}}_{S}({\mathbf{c}}-{\mathbf{x}}){\mathbf{b}}^{\top}$ has nonnegative entries. Then the right hand side above is nonnegative and we conclude

[TABLE]

As a result of this theorem, the Frobenius-norm minimizing Lipschitz matrix for a quadratic function is the solution to a finite-dimensional semidefinite program

[TABLE]

where $\{{\mathbf{c}}_{k}\}_{k=1}^{2^{m}}$ are the $2^{m}$ corners of ${\mathcal{D}}=[-1,1]^{m}$.

3.6 Convergence in Queries

With increasing point and gradient queries, the solution the finite-constraint semidefinite program Eq. 41 converges to the continuous-constraint program Eq. 40 if the function $f$ is Lipschitz. Let ${\mathbf{H}}_{M,N}$ denote the solution to Eq. 42 with ${\mathcal{P}}_{M}$ and ${\mathcal{G}}_{N}$. If ${\mathcal{P}}_{M}\subset{\mathcal{P}}_{M+1}$ and ${\mathcal{G}}_{N}\subset{\mathcal{G}}_{N+1}$ then $\operatorname{Tr}{\mathbf{H}}_{M,N}\leq\operatorname{Tr}{\mathbf{H}}_{M+1,N}$ and $\operatorname{Tr}{\mathbf{H}}_{M,N}\leq\operatorname{Tr}{\mathbf{H}}_{M,N+1}$. As $f$ is Lipschitz there exists a solution ${\mathbf{H}}$ to the continuous-constraint problem and $\operatorname{Tr}{\mathbf{H}}_{M,N}\leq\operatorname{Tr}{\mathbf{H}}$. Then as $\operatorname{Tr}{\mathbf{H}}_{M,N}\to\operatorname{Tr}{\mathbf{H}}$, ${\mathbf{H}}_{M,N}\to{\mathbf{H}}$.

3.7 Convergence Rate

An important practical question is how fast do finite query approximations converge to the continuous-constraint Lipschitz matrix Eq. 40? Unfortunately if we use random sampling this rate can be very slow. Figure 3 shows a quadratic approximation to the OTL Circuit function. As this is a quadratic function, we know the Lipschitz matrix is determined by the gradient at the corners. The probability of querying within $\epsilon$ of these points when sampling randomly is $\mathcal{O}(\epsilon^{-m})$; hence convergence with gradient queries is $\mathcal{O}(\epsilon^{-m})$ and convergence with point queries is $\mathcal{O}(\epsilon^{-2m})$. This is very slow! However, this slow convergence is not surprising as the Lipschitz matrix tracks the maximum rate of change.

Although this result is distressing, hope is not lost. Many problems in science and engineering are approximately monotonic over their domains. Hence we can employ a similar argument to Theorem 3.1 to justify that querying gradients at the corners of the domain yields an accurate approximation of the Lipschitz matrix. We employ this strategy throughout the remainder of this manuscript to estimate the Lipschitz matrix, potentially in combination with other random point or gradient queries.

4 Ridge Approximation Error Bound

In the introduction, Theorem 1.1 showed that functions with low-rank Lipschitz matrices ridge functions. Here we provide a related result that bounds the error made when approximating Lipschitz function by a ridge function. This result is analogous to the error bound associated with the MeGO matrix [4, Thm. 4.3].

Theorem 4.1**.**

Suppose $f\in{\mathcal{L}}({\mathcal{D}},{\mathbf{L}})$, $\widetilde{f}\in{\mathcal{L}}({\mathcal{D}},{\mathbf{L}}{\mathbf{U}}{\mathbf{U}}^{\top})$ where ${\mathbf{U}}\in\mathbb{R}^{m\times n}$ and ${\mathbf{U}}^{\top}{\mathbf{U}}={\mathbf{I}}$, and $\{{\mathbf{x}}_{j}\}_{j=1}^{M}\subset{\mathcal{D}}$, then

[TABLE]

*where $\sigma_{\max}({\mathbf{A}})$ denotes the largest singular value of ${\mathbf{A}}$ and $\operatorname{diam}({\mathcal{D}})$ denotes the diameter of the set ${\mathcal{D}}$; $\operatorname{diam}({\mathcal{D}})\coloneqq\max_{{\mathbf{x}}_{1},{\mathbf{x}}_{2}\in{\mathcal{D}}}\|{\mathbf{x}}_{1}-{\mathbf{x}}_{2}\|_{2}$. *

Proof 4.2**.**

Suppose ${\mathbf{x}}\in{\mathcal{D}}$ is fixed. Let $j=\operatorname*{argmin}_{k}\|{\mathbf{L}}{\mathbf{U}}{\mathbf{U}}^{\top}({\mathbf{x}}_{k}-{\mathbf{x}})\|_{2}$ and define $\epsilon\coloneqq\max_{j=1,\ldots,M}|f({\mathbf{x}}_{j})-\widetilde{f}({\mathbf{x}}_{j})|$. Inserting an additive identity,

[TABLE]

Invoking each function’s Lipschitz matrix and again inserting the identify,

[TABLE]

Finally, using the triangle inequality in the last term,

[TABLE]

Then defining $\delta\coloneqq\max_{{\mathbf{x}}\in{\mathcal{D}}}\min_{j=1,\ldots,M}\|{\mathbf{L}}{\mathbf{U}}{\mathbf{U}}^{\top}({\mathbf{x}}-{\mathbf{x}}_{j})\|_{2}$ and bounding the second term by the diameter, we obtain a result independent of ${\mathbf{x}}$

[TABLE]

We can then use this theorem to motivate a particular choice for ${\mathbf{U}}$ and function evaluations ${\mathbf{x}}_{j}$. First note that each term in this theorem has an important interpretation:

[TABLE]

To minimize the Lipschitz matrix approximation error we should choose the leading eigenvectors of ${\mathbf{L}}$ as the columns of ${\mathbf{U}}$ (or right singular vectors of ${\mathbf{L}}$ is not symmetric positive definite). With ${\mathbf{U}}$ selected, we can minimize dispersion by constructing a minimax design of experiments ${\mathcal{X}}({\mathcal{D}},M,{\mathbf{L}}{\mathbf{U}}{\mathbf{U}}^{\top})$ as in Eq. 18. As discussed in Section 2.1, the number of function evaluations required to obtain a particular dispersion (the covering number) no longer grows with the dimension of the domain, but instead the rank of the Lipschitz matrix ${\mathbf{L}}{\mathbf{U}}{\mathbf{U}}^{\top}$. If $f$ is a ridge function then we can make the function approximation error zero; otherwise, this is not the case. If two coordinates have the same projection, ${\mathbf{U}}^{\top}{\mathbf{x}}_{i}={\mathbf{U}}^{\top}{\mathbf{x}}_{j}$, and $f({\mathbf{x}}_{i})\neq f({\mathbf{x}}_{j})$, then the function approximation error is at least $|f({\mathbf{x}}_{i})-f({\mathbf{x}}_{j})|/2$.

Before concluding, we note two corollaries when the function approximation error is zero.

Corollary 4.3**.**

In the setting of Theorem 4.1, let ${\mathbf{U}}$ be the identity matrix and choose $\widetilde{f}$ such that $f({\mathbf{x}}_{j})=\widetilde{f}({\mathbf{x}}_{j})$, then

[TABLE]

Corollary 4.4**.**

In the setting of Theorem 4.1, let $f$ be a ridge function and ${\mathbf{U}}$ be a basis for the range of ${\mathbf{L}}$. Taking $\widetilde{f}$ such that $f({\mathbf{x}}_{j})=\widetilde{f}({\mathbf{x}}_{j})$, then

[TABLE]

5 Design of Experiments

Given a particular function, where should we evaluate it to provide the most information? This is the subject of the design of computer experiments [28], a subfield of experimental design (see, e.g., [9]) distinguished by the assumption that observations are deterministic; e.g., $f({\mathbf{x}})$ returns only one value. From the ridge approximation error bound in Theorem 4.1 and the earlier complexity results in Section 2.2, a good experimental design of points $\{{\mathbf{x}}_{j}\}_{j=1}^{M}\subset{\mathcal{D}}$ should minimize the dispersion in the Lipschitz matrix metric:

[TABLE]

Minimizing the dispersion yields in a minimax optimal design; cf. Eq. 18:

[TABLE]

There is a substantial body of literature on solving this problem; see, e.g., [24] for a recent review. Here we provide a motivating example illustrating why an minimax design in the Lipschitz metric is important. Then we provide a brief description of how we construct locally optimal low-dispersion designs and illustrate the performance of this algorithm.

5.1 A Motivating Example

Consider a ridge function

[TABLE]

with a one-dimensional active subspace $\operatorname{Span}\{{\mathbf{1}}\}$ and a rank-one Lipschitz matrix ${\mathbf{L}}={\mathbf{1}}{\mathbf{1}}^{\top}$. If we randomly select points $\{{\mathbf{x}}_{j}\}_{j=1}^{M}$ with uniform probability over the domain, their sum tends to concentrate around the mean (zero) when projected onto the active subspace as seen in Fig. 4. The same is true with a Latin hypercube design where points are randomly selected so that projection onto the coordinate axes results in evenly spaced points. However, this is not true when we construct a minimax design under the Lipschitz matrix metric. This design has much lower dispersion and consequently we can construct much more accurate approximations according to Theorem 4.1.

5.2 Constructing Minimax Designs

Finding even an approximately optimal minimax design is challenging: it requires solving a deeply nested optimization problem Eq. 62 in a high dimensional space with many local minimizers with large dispersion. Here we briefly describe a combination of three algorithms that we use to construct the minimax designs appearing in the remainder of this paper. Each algorithm primarily uses bounded Voronoi vertices under the Lipschitz matrix metric. All such vertices can be computed using Qhull [2] or a random subset by sequential projection [19]; using sequential projection allows us to employ this combination algorithms spaces of moderate dimension, i.e., $m>4$.

To avoid finding a local minimizer with a large dispersion, the first two algorithm construct a good initial design that is then refined to local optimality. First, we construct a maximin coffeehouse design [21] where we greedily add new points to the design maximizing their distance from existing points in the design

[TABLE]

This optimizer ${\mathbf{x}}_{k}$ is a bounded Voronoi vertex that we can (approximately) identify by a finite minimization over (a subset of) these vertices. Second, we perform a few iterations of block coordinate descent [30] to maximize the minimum pairwise distance

[TABLE]

Again, we can (approximately) solve this optimization problem by a finite minimization over (a subset of) the bounded Voronoi vertices. Finally this design initializes a variant of Lloyd’s algorithm to find a locally optimal minimax design [33]. At each iteration we identify the bounded Voronoi vertices ${\mathbf{v}}_{i}^{\ell}$ associated with the design $\{{\mathbf{x}}_{j}^{\ell}\}_{j=1}^{M}$ and then move ${\mathbf{x}}_{j}^{\ell}$ to be the circumcenters of its Voronoi cell

[TABLE]

If this iteration converges, the design satisfies the local optimality conditions for a minimax design [6, Thm. 4.7]. Note that if the Lipschitz matrix is low-rank (i.e., rank one, two, or three) it feasible to compute all the bounded Voronoi vertices even if the ambient space is high dimensional.

It is tempting to think that a minimax design is only necessary when function evaluations are costly. However, as Fig. 5 illustrates increasing the number of random and Latin hypercube samples decreases dispersion only very slowly at $\mathcal{O}(M^{-1/6})$ whereas we see the minimax design temporally converges faster and has a substantially reduced dispersion (see discussion in subsections 2.3 and 2.4). Hence if we seek a low-dispersion design for accurate approximation, integration, or optimization it is critical to construct an approximate minimax design.

6 Parameter Reduction

Using the Lipschitz matrix, we can perform dimension reduction by constructing a ridge approximation

[TABLE]

Since $g$ is a ridge function, $g$ has at most a rank-$n$ Lipschitz matrix (Theorem 1.1) and complexity of tasks now scales with the ridge dimension $n$ and no longer the dimension of the ambient space $m$. There are many ways to construct ridge approximations; e.g., picking ${\mathbf{U}}$ from the dominant eigenvectors of the MeGO matrix and then fitting a polynomial [4] or directly picking ${\mathbf{U}}$ and polynomial $g$ via optimization [5, 14]. Here we construct a ridge approximation using the Lipschitz matrix in light of Theorem 4.1.

6.1 Building a Ridge Approximation

Our choice of the ridge subspace is clear from Theorem 4.1: the leading $n$ eigenvectors of the Lipschitz matrix form ${\mathbf{U}}\in\mathbb{R}^{m\times n}$, same as when using the MeGO matrix. To construct a ridge approximation satisfying Theorem 4.1 we need to choose a ridge function $\widetilde{f}\in{\mathcal{L}}({\mathcal{D}},{\mathbf{L}}{\mathbf{U}}{\mathbf{U}}^{\top})$. Due to the constraints its Lipschitz matrix, $\widetilde{f}$ may not interpolate point queries; i.e., $\widetilde{f}({\mathbf{x}}_{j})\neq f({\mathbf{x}}_{j})$. Hence to a compatible ridge approximation $\widetilde{f}$ we first solve an optimization problem for function values $\widetilde{y}_{j}=\widetilde{f}({\mathbf{x}}_{j})$ minimizing the $\ell_{1}$ error (or equivalently any other convex norm)

[TABLE]

Then these function values define point queries $\widetilde{\mathcal{P}}_{M}=\{{\mathbf{x}}_{j},\widetilde{y}_{j}\}_{j=1}^{M}$ and define $\widetilde{f}$ to be the central approximation of this data, cf. Eq. 28:

[TABLE]

As $\widetilde{f}$ approximates $f$ with a maximum error $\epsilon=\max_{j}|y_{j}-f({\mathbf{x}}_{j})|$, the uncertainty associated with $\widetilde{f}$ is, cf. Eq. 27,

[TABLE]

6.2 Typical Workflow

A typical workflow for using the Lipschitz matrix for dimension reduction starts with querying the gradient on the corners of the domain and using this data to estimate the Lipschitz matrix. Then, picking the active subspace to be the span of the leading eigenvector ${\mathbf{U}}\in\mathbb{R}^{m\times 1}$ of ${\mathbf{L}}$, we construct a 20-point minimax design under this metric ${\mathbf{L}}{\mathbf{U}}{\mathbf{U}}^{\top}$, $\{{\mathbf{x}}_{j}\}_{j=1}^{20}$. To better estimate the variability at each point, we take 5 maximin coffeehouse samples on each domain ${\mathcal{D}}_{j}=\{{\mathbf{x}}\in{\mathcal{D}}:{\mathbf{U}}^{\top}{\mathbf{x}}={\mathbf{U}}^{\top}{\mathbf{x}}_{j}\}$ in the Lipschitz metric ${\mathbf{L}}$. Finally, we construct the ridge approximation as described in the previous section. Figure 6 shows three examples of this workflow and the associated uncertainty.

7 Uncertainty Quantification

When working with expensive deterministic computer simulations it is often necessary to employ an approximation of certain quantities of interest, called a response surface or a surrogate. Supposing we have constructed this approximation using samples of $f$, it is natural ask: what are the range of possible values our approximation could take away from these samples? This is often called uncertainty in this setting. Gaussian processes provide one approach to define an uncertainty [25, sec. 2.2] and the Lipschitz constant provides another [26]. These two techniques are based on different assumptions and yield different results as illustrated in Fig. 7. The Gaussian process approach assumes that $f$ is a Gaussian process conditioned on the measurements $y_{j}=f({\mathbf{x}}_{j})$, maximizes the likelihood of a parameterized covariance kernel based on observations, and then computes the probability of observing $f({\mathbf{x}})$; the uncertainty is visualized as all function values above some probability threshold. In contrast, the Lipschitz approach assumes measurements $y_{j}=f({\mathbf{x}}_{j})$ comes from a Lipschitz function, chooses the smallest Lipschitz constant that consistent with this data, and defines uncertainty as the range of possible function values consistent with this Lipschitz constant and data; namely, the uncertainty set ${\mathcal{U}}$ Eq. 23. Replacing the Lipschitz constant by the Lipschitz matrix allows us to reduce uncertainty.

7.1 Decreasing Uncertainty

As the Lipschitz matrix more accurately encodes variability in the function, it can reduces uncertainty compared to the scalar Lipschitz constant. Figure 8 illustrates the improvement of the Lipschitz matrix over the Lipschitz constant on a two dimensional example. Table 2 demonstrates the same reduction in uncertainty occurs in higher dimensional examples. This example shows that uncertainty is reduced the most when Lipschitz matrix is used both when constructing the minimax design and evaluating the uncertainty set.

7.2 Visualizing Uncertainty on Shadow Plots

Shadow plots are an important tool for visualizing functions with a high-dimensional domain. To include the uncertainty in shadow plots, we generalize the uncertainty set Eq. 23 to take set-valued inputs

[TABLE]

In our implementation, we identify these lower and upper bounds by solving a sequential linear program [22]. Then to include the uncertainty in a one-dimensional shadow plot we evaluate the uncertainty set for ${\mathcal{S}}_{\alpha}=\{{\mathbf{x}}\in{\mathcal{D}}:{\mathbf{u}}^{\top}{\mathbf{x}}=\alpha\}$ for multiple $\alpha$ in the projection of the domain onto ${\mathbf{u}}$. Figure 9 provides an example of this projected uncertainty using both the Lipschitz constant and Lipschitz matrix as well as with an isotropic minimax design or a Lipschitz matrix minimax design. Note the uncertainty in Fig. 9 is the projection of a high-dimensional uncertainty Eq. 71 rather than an intrinsically low-dimensional uncertainty Eq. 70 in Fig. 6.

8 $\epsilon$-Lipschitz

In this section we explore three applications of the $\epsilon$-Lipschitz matrix exploiting its ability to ignore a small changes in the function.

8.1 Computational Noise

Many functions appearing in computational science and engineering often have computational noise [20]—a phenomenon emerging from many factors including convergence tolerances and mesh discretizations. With the addition of computational noise, functions that are ideally Lipschitz continuous can become discontinuous. The $\epsilon$-Lipschitz matrix enables us to perform dimension reduction for these discontinuous functions and remove the influence of noise. As an example, consider the partial trace function of Moré and Wild [20, Sec. 1]: the sum of the first five smallest eigenvalues of the parameterized matrix

[TABLE]

where ${\mathbf{A}}\in\mathbb{R}^{700\times 700}$ is the Trefethen_700 sparse matrix [7], and $\lambda_{i}$ denotes the $i$th eigenvalue in decreasing order. If we accurately evaluate this function it is approximately linear on the domain and yields an approximately low-rank Lipschitz matrix. To introduce computational noise, we use a very loose relative accuracy termination criteria of $0.1$ when computing these eigenvalues using ARPACK [18]. As evidenced in Table 3, the introduction of noise pollutes the estimate of the Lipschitz matrix, yielding a large, full-rank Lipschitz matrix. If we estimate an $\epsilon$-Lipschitz matrix instead with $\epsilon=6.47$ (the largest mismatch between accurate and noisy evaluations of this function) we are able to recover an accurate, approximately low-rank Lipschitz matrix. In this example we use samples on a $10\times 10$ grid and initialize the Lanczos iteration using the ones vector.

8.2 Distracting Oscillations

The active subspace identified using a Lipschitz matrix can yield counterintuitive results when a function has high frequency, low amplitude oscillations. Consider the “corrugated roof” function [5, eq. (26)]:

[TABLE]

Most of the variation in this function is due to $x_{1}$, but there is a high frequency oscillation in $x_{2}$. As $f$ is additive in the two coordinates, we can identify a Lipschitz matrix ${\mathbf{L}}=\begin{bmatrix}5&0\\ 0&10\pi\end{bmatrix}$ by considering the largest derivative in each coordinate. Then applying Theorem 4.1, we would choose the active subspace $\operatorname{Span}\{{\mathbf{e}}_{2}\}$. However, as illustrated in Fig. 10, this yields much higher uncertainty than using $\operatorname{Span}\{{\mathbf{e}}_{1}\}$ because the oscillations in $x_{2}$ have less influence on the output than the linear term in $x_{1}$. We can remove the influence of these oscillations by using the $\epsilon$-Lipschitz matrix. Taking $\epsilon=2$, we can ignore the sine term and identify an $\epsilon$-Lipschitz matrix ${\mathbf{L}}_{2}=\begin{bmatrix}5&0\\ 0&0\end{bmatrix}$ With this choice, we choose the active subspace $\operatorname{Span}\{{\mathbf{e}}_{1}\}$. Although this has increased the uncertainty at each point, it yields better active subspace which reduces overall uncertainty.

8.3 Dimension Reduction

Just as we can use an $\epsilon$-Lipschitz matrix to remove undesirable oscillations, we can also use it for dimension reduction. If there is a subspace ${\mathcal{U}}$ on which a function varies by less than $\epsilon$, then there is an $\epsilon$-Lipschitz matrix that has a nullspace of $\dim({\mathcal{U}})$. We have no guarantee that by minimizing the Frobenius norm with the given function evaluations we will identify this low-rank $\epsilon$-Lipschitz matrix. However Figure 11 illustrates that we can sometimes find a low-rank matrix. In this example we note that by ignoring about $20\%$ of the variation, we can identify a rank-1 $\epsilon$-Lipschitz matrix.

Acknowledgements

The authors would like to thank Akil Narayan and Drew Kouri for their comments that helped improve this manuscript.