Hierarchical Data Reduction and Learning

TL;DR

This paper introduces a hierarchical learning method for creating sparse data representations that effectively reduce data size while maintaining accuracy, with proven stability and convergence properties across various datasets.

Contribution

It presents a novel hierarchical approach for data reduction that guarantees stability, convergence, and efficient data reconstruction, applicable to diverse datasets.

Findings

Effective data reduction for synthetic and real datasets

Accurate data reconstruction with minimized prediction error

Stable and convergent approximation process

Abstract

This paper describes a hierarchical learning strategy for generating sparse representations of multivariate datasets. The hierarchy arises from approximation spaces considered at successively finer scales. A detailed analysis of stability, convergence and behavior of error functionals associated with the approximations are presented, along with a well chosen set of applications. Results show the performance of the approach as a data reduction mechanism for both synthetic (univariate and multivariate) and real datasets (geospatial and numerical model outcomes). The sparse representation generated is shown to efficiently reconstruct data and minimize error in prediction.

Equations

$$\underset{C_s\in {\mathbb{R}}^{l_s}}{\min }||f{|}_X-B^s{C}_s|{|}_2\le TOL$$

$$\underset{\tilde{f}\in \mathcal{Q}}{\min }||f-\tilde{f}|{|}_{\mathcal{V}}=||f-(Af)|{|}_{\mathcal{V}}$$

$$<f,K(x,\cdot ){>}_{\mathcal{H}}=f(x)x\in \Omega ,f\in \mathcal{H}$$

$$K(x,y)=<K(x,\cdot ),K(y,\cdot ){>}_{\mathcal{H}}=<{\delta }_x,{\delta }_y{>}_{{\mathcal{H}}^{\ast }}x,y\in \Omega$$

$${\Gamma }^s=span\{K^s(\cdot ,x_i):x_i\in X_s\}\approx span\{K^s(\cdot ,x_j):x_j\in X\},X_s\subseteq X$$

$$T=2(Diameter(X)/2{)}^2$$

$$G_s(x,y)=\exp \left(-\frac{||x-y|{|}^2}{{ϵ}_s}\right)\text{mbox};{ϵ}_s=\frac{T}{P^s}$$

$$R_{\delta}(G_{s})=\#\Bigg{(}j:\frac{\sigma_{j}(G_{s})}{\sigma_{0}(G_{s})}\geq\delta\Bigg{)}$$

$$R_{\delta}(G_{s})\leq\prod_{i=1}^{d}\Bigg{(}\frac{2|I_{i}|}{\pi}\sqrt{\epsilon_{s}^{-1}ln(\delta^{-1})}+1\Bigg{)}$$

$$R_{\delta }(G_s)\propto {ϵ}^{-d/2}\propto P^{sd}\propto s;\text{for }P=2$$

$$||f{|}_X-(A_{j}f){|}_X|{|}_2\le ||f{|}_X-(A_{i}f){|}_X|{|}_2$$

$$\underset{s\to S_c}{\lim }(A_{s}f){|}_X=f{|}_X$$

$$||f{|}_X-(A_{s}f){|}_X|{|}_2\le TOL$$

$$||E^s{|}_X||=||f{|}_X-(A_{s}f){|}_X||\le TOL$$

$$\begin{array}{rl}||E^s{|}_X|{|}_2& =||f{|}_X-R^{s}f{|}_X|{|}_2=||(I-R^s)f{|}_X|{|}_2\le ||I-R^s|{|}_2||f{|}_X|{|}_2\end{array}$$

$$\underset{s\to S_c}{\lim }||E^s{|}_X|{|}_2\to 0$$

$$(A_{s+1}f){|}_X=(A_{s}f){|}_X+{\alpha }_s\cdot E^s{|}_X$$

$${\Gamma }^s=span\{K^s(.,x_i):x_i\in X_s\}{X}_s\subseteq X$$

$$(A_{s+1}f){|}_X=(A_{s}f){|}_X+{\alpha }_s\cdot E^s{|}_X$$

$${\alpha }_s=\frac{<(A_{s+1}f-A_{s}f){|}_X,E^s{|}_X>}{<E^s{|}_X,E^s{|}_X>}=\frac{<(A_{s+1}f-A_{s}f){|}_X,E^s{|}_X>}{||E^s{|}_X|{|}^2}$$

$${\alpha }_s\le \frac{||(A_{s+1}f-A_{s}f){|}_X||\cdot ||E^s{|}_X||}{||E^s{|}_X|{|}^2}=\frac{||(A_{s+1}f-A_{s}f){|}_X||}{||E^s{|}_X||}\le 1+\frac{||E^{s+1}{|}_X||}{||E^s{|}_X||}$$

$$0\le {\alpha }_s\le 1+\frac{||E^{s+1}{|}_X||}{||E^s{|}_X||}$$

$$(A_{s+1}f-f){|}_X=(A_{s}f-f){|}_X+{\alpha }_s(f-A_{s}f){|}_X\text{ leads to}$$

$${\rho }_s=\frac{||(f-A_{s+1}f){|}_X||}{||(f-A_{s}f){|}_X||}=|1-{\alpha }_s|$$

$$\lim_{s\to S_{a}}\Big{|}\Big{<}(A_{s}f),E^{s}\Big{>}_{\mathcal{H}_{s}}\Big{|}\ \leq||C_{S_{a}}||_{\infty}\sqrt{n}(TOL)$$

$$\displaystyle\Big{|}\Big{<}(A_{s}f),E^{s}\Big{>}_{\mathcal{H}_{s}}\Big{|}$$

$$\displaystyle=\Big{|}\Big{<}\sum_{x_{j}\in X_{s}}K^{s}(\cdot,x_{j})c_{j},f\Big{>}_{\mathcal{H}_{s}}-\Big{<}\sum_{x_{j}\in X_{s}}K^{s}(\cdot,x_{j})c_{j},A_{s}f\Big{>}_{\mathcal{H}_{s}}\Big{|}$$

$$\displaystyle=\Big{|}\sum_{x_{j}\in X_{s}}f(x_{j})c_{j}-\sum_{x_{j}\in X_{s}}(A_{s}f)(x_{j})c_{j}\Big{|}\quad\text{(Reproducing Property)}$$

$$\displaystyle=\Big{|}\sum_{x_{j}\in X_{s}}c_{j}\Big{(}f(x_{j})-(A_{s}f)(x_{j})\Big{)}\Big{|}$$

$$\displaystyle\leq\sum_{x_{j}\in X_{s}}\Big{|}c_{j}\Big{(}f(x_{j})-(A_{s}f)(x_{j})\Big{)}\Big{|}\quad\text{(Triangle Inequality)}$$

$$\displaystyle\leq||C_{s}||_{\infty}\sum_{x_{j}\in X_{s}}\Big{|}\Big{(}f(x_{j})-(A_{s}f)(x_{j})\Big{)}\Big{|}$$

$$\displaystyle\leq||C_{s}||_{\infty}\sum_{x_{j}\in X}\Big{|}\Big{(}f(x_{j})-(A_{s}f)(x_{j})\Big{)}\Big{|}\quad(X_{s}\subseteq X)$$

$$=||C_s|{|}_{\infty }||E^s{|}_X|{|}_1$$

$$|<(A_{s}f),E^s{>}_{{\mathcal{H}}_s}|$$

Full text

{AMS}

68W25, 65D15, 33F05

1 Introduction

Hierarchical and multiscale modeling is widely used both in data driven and physics based models [24, 12, 15, 20]. In this paper, we focus on a strategy which analyzes data at different scales ($s$) by constructing corresponding basis representations $B^{s}\in\mathbb{R}^{n\times l_{s}}$ and inferring if these scale and data dependent bases are able to approximate the observed data $f|_{X}\in\mathbb{R}^{n}$ ($X\in\mathbb{R}^{n\times d}$) to satisfy:

[TABLE]

where, $l_{s}\in\mathbb{N}$ is a scale enumeration chosen by the algorithm, $TOL$ is a user defined tolerance level, and $n$ is the number of observations in a $d$ dimensional space ($x_{i}\in\mathbb{R}^{d}$, $1\leq i\leq n$). $C_{s}$ in (1) is computed as the optimal coordinate of projection on $B^{s}$. Besides constructing these bases at each scale, the algorithm identifies “representative” data points which form the corresponding sparse data $X_{s}$. $\{s,X_{s},C_{s}\}$ allow reconstruction of the original dataset and prediction (model) at any new point in the domain of interest resulting in both a sparse representation [17, 8, 22] and an efficient model [11, 5].

Data sparsification seeks to find a subset $X_{s}\subseteq X=\{x_{1},x_{2},..,x_{n}\}\subset\Omega$ ($x_{i}\in\mathbb{R}^{d}$ for $d$ dimensional inputs) for which we have data $f|_{X}={(f(x_{1}),f(x_{2}),...,f(x_{n}))}^{T}\in\mathbb{R}^{n}$, such that, we may obtain an acceptable reconstruction or approximation to $f$ as in equation (1). This subset with the corresponding projection coordinate $C_{s}$ with respect to the basis set $B^{s}$, is regarded as our obtained sparse representation ($D^{s}_{sparse}=[X_{s},C_{s}]$). Thus, this compressed representation of the dataset acts as a model for new prediction points $X_{*}\subset\Omega$ and reconstruction of the original dataset $f|_{X}$ at $X\subset\Omega$ as needed. The approximation problem [6, 19] is a more traditional one seeking to construct an approximation $(Af):\mathbb{R}^{d}\to\mathbb{R}$ for the underlying function $f$ which was discretely observed ($f|_{X}$) such that it minimizes the normed error measure $||(Af)-f||$. Following the standard notation, for a normed space $\mathcal{V}$ and any approximation $\tilde{f}$ to $f$ such that $f,\tilde{f}\in\mathcal{V}$, the quality of approximation $\tilde{f}$ is usually quantified by the error norm ${||f-\tilde{f}||}_{\mathcal{V}}$. $\tilde{f}$ here is usually constructed by first fixing a finite dimension approximation space $\mathcal{Q}\subseteq\mathcal{V}$ and then computing the optimal approximation $(Af)\in\mathcal{Q}$ by minimizing normed distance between $f$ and $\mathcal{Q}$:

[TABLE]

We now motivate our development of a data dependent bases ($B^{s}\equiv B(X)^{s}$) [9]. The Mairhuber-Curtis theorem 111Theorem 2 in [14] states that for $d\geq 2$ and $n\geq 2$, there is no Chebyshev system $B=\{s_{1},...s_{n}\}$ on $\mathbb{R}^{d}$, therefore for some data $X=\{x_{1},x_{2},..,x_{n}\}\subset\mathbb{R}^{d}$, the Vandermonde matrix approaches singular behavior affecting the quality of approximation. As a remedy, we choose a data dependent continuous kernel function $K:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}$ which is positive definite and symmetric on $\Omega$. Thus, for data points $X=\{x_{1},x_{2},....,x_{n}\}$, the kernel function becomes $K(x_{i},x_{j})_{1\leq i,j\leq n}$. From the literature on these functions [2, 4, 26] these kernels are reproducing in the native Hilbert Space $\mathcal{H}$ of functions on $\Omega$ in the sense

[TABLE]

The functionals $K(x,.)$ are the Reisz representers of the linear point functionals $\delta_{x}:f\mapsto f(x)$ in the dual space $\mathcal{H}^{*}$ of $\mathcal{H}$. Therefore we have the relationship

[TABLE]

For observed data at $X=\{x_{1},x_{2},...,x_{n}\}$, one can consider the space of trial function $\mathcal{D}_{X}=span\{K(\cdot,x_{j}):x_{j}\in X\}$ as a choice for the approximation space. However from work such as [9, 10] we know that the bases formed by translates of the kernel functions are highly dependent on the distribution of data points in X and hence numerically ill-conditioned for many problems. For dealing with this issue, at each scale s, our approach considers a subset of original trial functions as bases such that

[TABLE]

This is equivalent to identifying a lower dimension manifold in the space $\mathcal{D}_{X}$ (as compared to original dimension n for n data points, represented by the kernel matrix $K$). These linearly independent data derived basis functions are sampled using a Pivoted QR strategy as detailed in the following section (Algorithm 1). It should be noted here that at each scale we have a different kernel function $K^{s}$. Hence, effectively we are working in a different RKHS at every scale (from the uniqueness property [2] of kernel functions). Let this scale dependent RKHS be represented as $\mathcal{H}_{s}$. Therefore at any scale s, the function to be approximated ($f$) is assumed to be in the RKHS $\mathcal{H}_{s}$. However, our approach is aimed at finding an approximation $A_{s}f\in\Gamma^{s}$ to $f$ as in (5). Hence, if we look at the problem formulation (2), then $\mathcal{H}_{s}$ is the space $\mathcal{V}$ and $\Gamma^{s}$ is our approximation space $\mathcal{Q}$. While our overall approaches are somewhat related to statistical emulators [13], the notion of scale is rarely articulated in that domain.

The contributions of this paper can be summarized as follows. Firstly we introduce an approach which computes a scale ($s$) dependent sparse representation $D^{s}_{sparse}$ of a large dataset. Besides data reduction, the algorithm also provides a model for approximations at new data points in the observation space using just $D^{s}_{sparse}$, enabling any sort of learning from a compressed version of the dataset. We then provide detailed analysis on convergence, establishing the optimality of the solution obtained by the algorithm, computation of point-wise error functionals and their behavior which governs the performance of the approach. At each scale, besides providing confidence bounds which quantify our belief in the estimation, we also provide prediction intervals which estimates the bounds in which the algorithm believes any observations not included in the sparse representation as well as any new samples to be made in the future will lie. This also makes the sparse representation more useful. Towards the end, we have also provided stability estimates to further establish the dependability of the approach.

2 The Hierarchical Approach

In this approach, developing further the work of [3], we introduce a methodology for data reduction through efficient bases construction exploiting the multilevel nature of the correlation structure present in the data. The basic steps involved in the approach are presented as Algorithm 1, here.

2.1 Proposed Algorithm

The proposed approach (Algorithm 1) constructs a sequence of scale (s) dependent approximations $A_{1}f,A_{2}f,..,A_{s}f..$ to the unknown function $f:\mathbb{R}^{d}\mapsto\mathbb{R}$. Each of these approximations only use a subset of dataset $X_{1},X_{2},...,X_{s},..$ respectively. Algorithm begins by taking a dataset, where a data point $x_{i}\in\mathbb{R}^{d}$ is mapped to a functional value $f(x_{i})\in\mathbb{R}$. In matrix form $f|_{X}=\{f_{1},f_{2},....f_{n}\}$ values are obtained at data points $X=\{x_{1},x_{2},....,x_{n}\}$ ($f|_{X}\in\mathbb{R}^{n}$ and $X\in\mathbb{R}^{n\times d}$). The scalars $[T,P,TOL]\in\mathbb{R}^{3}$ are the algorithmic hyperparameters defined by the user. TOL is the simple $2$-$norm$ error tolerance, P is assumed to be 2 (Based on [3]). This choice of P reduces the length scale of the kernel (K) by a factor of 0.5 at each scale increment, providing an intuitive understanding of how the support of basis functions is adapted to scale variation. If we assume the diameter of the dataset to be distance between the most distant pair of datapoints, then T is given by

[TABLE]

Besides these hyperparameters, the algorithm also accepts the prediction points $X_{*}\subset\Omega\in\mathbb{R}^{n^{*}\times d}$, which represent the data points at which the user wants to approximate the underlying function. We also need to choose the positive definite function ($K:\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}$). For this paper we use the squared exponential (also referred to as Gaussian) kernel (7) for mapping the covariance structure and generating the space of trial functions $\Gamma^{s}$(equation 5) at each scale s.

[TABLE]

Here $\epsilon_{s}$ (also known as the length scale parameter) determines the width of the basis functions at a particular scale. This squared exponential kernel is very widely used in the Gaussian process literature [21] and is a popular choice for most spatial statistics problems. The learning phase (for generating $D^{s}_{sparse}$) of the proposed algorithm has been explained in $STEP-1$ and $STEP-2$ below. $STEP-3$ constitutes the prediction phase which uses the produced sparse representation for the dataset. It should be noted that, it is not required to wait till the convergence of the algorithm to go to the prediction phase. In fact, it is possible to predict from the sparse representation at each scale ($D^{s}_{sparse}$). We have shown the prediction phase separately for the purpose of clarity.

STEP-1 (Getting sparse data-$X_{s}$): Given the Dataset $D=[X,f|_{X}]$, the algorithm begins with the computation of the covariance operator $G_{s}$ (7). However, based on [9, 10], the distribution of the dataset might lead to ill-conditioning of this covariance kernel. Therefore we carry out a column pivoted QR decomposition to identify the space $\Gamma^{s}$ (at each scale) which represents the span of the trial functions $K^{s}(\cdot,x_{j}),(1\leq j\leq n)$ at scale s. The QR decomposition is carried out on W for obtaining the Permutation matrix $P_{R}$. $W$ is produced by the product of a random normal matrix $A$ with the $G_{s}$. Here we have $A\in\mathbb{R}^{k\times n}$ with $l_{s}=rank(G_{s})\leq k\leq n$. For our experiments we have assumed $k=l_{s}+8$ (as in [3]) which means we sample 8 additional rows to account for numerical round-offs during the QR decomposition. However, this is a conservative step and even without any additional sampled columns, the algorithm was found to perform well. The permutation matrix $P_{R}$ produced by the decomposition captures the information content of each column of W. $P_{R}$ is then used to extract independent columns with the biggest norm contributions along with the observation points these columns correspond to in the covariance kernel ($G_{s}$). This set of sampled observations from the original dataset is termed as the corresponding sparse data ($X_{s}$). The number of columns sampled from $G_{s}$ come from its numerical rank estimated by using strategies such as a $Rank\ Revealing-QR$ or a $SVD$ decomposition.

STEP-2 (Getting projection coordinate-$C_{s}$): Once, we have the relevant columns of $G^{s}$ (denoted as $B^{s}$) which also represent the approximation subspace $\Gamma^{s}$ (in the native RKHS) which spans $\mathcal{D}_{X}$, the algorithm proceeds to solve the over-constrained system ($B^{s}C_{s}=f|_{X}$). We can think of it as an orthogonal projection problem where $f|_{X}$ needs to be projected on the column space of $B^{s}$ and then it is required to compute the specific weighting of vectors in the basis matrix $B^{s}$ which produces this projection. The Algorithm computes the orthogonal projection ($(A_{s}f)|_{X}$) given the required coordinates and basis vectors in $B^{s}$. Once we have the scale $s=S_{a}$ at which the algorithm satisfies the $2$-$norm$ condition, it produces the sparse representation $D^{S_{a}}_{sparse}=\{X_{S_{a}},C_{{S_{a}}}\}$ and proceeds to the prediction stage if required by the user. Here $S_{a}$ is called the convergence scale as detailed in the following subsection.

STEP-3 (Prediction at $X_{*}$ from $D^{S_{a}}_{sparse}$): For computing functional values at unobserved location $X_{*}$, basis functions are constructed by computing $G^{*}_{S_{a}}$ (Gaussian kernel for the prediction points with the sparse data - $X_{S_{a}}$). The Prediction step weighs the constructed bases with coordinates of orthogonal projection $C_{s}$ ($C_{S_{a}}$, if the prediction is being made at the convergence scale) and linearly combines them to produce the required approximation.

Finally, before moving forward, it is worth mentioning here that the $2$-$norm$ criteria is just one of the many possible kinds of norms which can be used to measure the scale dependent fidelity of the model.

2.2 Critical Scale ($S_{c}$) and Convergence Scale ($S_{a}$)

In our work, the scale at which the kernel matrix becomes numerically full rank is referred to as Critical Scale. Working with proposition 3.7 in [3], if $\delta$ represents the precision of rank for the Gaussian kernel matrix, then we can define the numerical rank of the kernel as

[TABLE]

where $\sigma_{j}(G_{s})$ is the $j^{th}$ largest singular value of $G_{s}$. Also let $|I_{i}|$ represents the length of the bounding box of the data in $i^{th}$ ($i\in[1,d]$) dimension. Then given the length scale parameter $\epsilon_{s}$, the rank of the Gaussian kernel can be bounded above as

[TABLE]

Proposition 3.7 in [3] states that numerical rank of the Gaussian kernel matrix is proportional to the volume of the minimum bounding box $B=I_{1}\times I_{2}\times....\times I_{d}$ and to $\epsilon^{-d/2}$ . Therefore for a fixed data distribution, following relation holds

[TABLE]

Hence numerical rank of the Gaussian covariance kernel is directly proportional to the scale of study. Therefore, there exists a minimum scale $s=S_{c}$, at which the kernel becomes full and continue to stay full rank as the scale is further increased.This scale is referred to as the Critical Scale ($S_{c}$). Hence, based on the overall Algorithm 1, If $X_{1},X_{1}..X_{s}..$ and so on, are the sampled sparse representation at scale s, they satisfy $|X_{1}|\leq|X_{2}|\leq|X_{3}|\leq.....\leq|X_{S_{c}}|$ where $|\cdot|$ denotes the cardinality of a set. Therefore with increasing scales more data points are added leading to better approximation. Thus, for any scales $j$ and $i$ satisfying the relation $j\geq i$, the following relation holds

[TABLE]

From approximation theory [7], we can now state that for kernel at the Critical Scale and a finite point set $X\subset\Omega$: (a) $K^{S_{c}}$ is positive definite, (b) If $(A_{{S_{c}}}f)\in\Gamma^{S_{c}}$ vanishes on $X$, then $(A_{{S_{c}}}f)\equiv 0$, and, c) the approximation $(A_{{S_{c}}}f)$ satisfying the following relation is unique

[TABLE]

We can now define the Convergence scale ($S_{a}$) as the minimum scale $s$ which satisfies

[TABLE]

Therefore it is the scale at which Algorithm 1, stops. It is worth noting here that based on the definition of $S_{c}$ and $S_{a}$, we have $S_{a}\leq S_{c}$ and $|X_{S_{a}}|\leq|X_{S_{c}}|$

3 Convergence Properties

Following the notations introduced in the previous sections, again let $f$ be the function which needs to be approximated and $A_{s}f$ be the approximation produced by the Algorithm 1 at scale s. Correspondingly, let $E^{s}=f-A_{s}f$ be the error in approximation at scale $s$ with it’s $2$-$norm$ in $\mathbb{R}^{n}$ as $||E^{s}|_{X}||_{2}=||f|_{X}-(A_{s}f)|_{X}||_{2}$.

Lemma 3.1**.**

For any choice of $TOL\geq 0$ and any norm: $||\cdot||\in\mathbb{R}^{n}$ chosen by the user, the proposed algorithm (Algorithm 1) converges in the sense

[TABLE]

*for some finite $s\leq S_{c}$ and the corresponding produced approximation $A_{s}f$ *

Proof 3.2**.**

$(A_{s}f)|_{X}$* is the projection of $f|_{X}$ on the bases at scale $s$ with the projection operator $B^{s}{B^{s}}^{\dagger}\equiv R^{s}$, where ${B^{s}}^{\dagger}$ is the pseudo-inverse of basis matrix $B^{s}$. Therefore*

[TABLE]

From Equation (10), that numerical rank of bases $B^{s}$ increases monotonically with scale $s$. Hence $s\to S_{c}\Rightarrow R^{s}\to I_{n}$. Thus $||I-R^{s}||_{2}\to 0$ and hence $||E^{s}|_{X}||_{2}\to 0$.

[TABLE]

*This establishes the convergence for the proposed approach in $2$-$norm$ at some finite scale $S_{a}\leq S_{c}$. From equivalence of norms on finite dimensional vector space [16] this generalizes the convergence guarantee to all possible norm in $\mathbb{R}^{n}$. *

In the following theorem we analyze a particular hierarchical updating scheme.

Theorem 3.3**.**

If the projection update of Algorithm 1 is written in an iterative form

[TABLE]

*where $E^{s}=f-A_{s}f$. Then $\alpha_{s}$ follows the bounds $0<\alpha_{s}<2$ and convergence rate ($\rho_{s}$) satisfies the relation $\rho_{s}=|1-\alpha_{s}|$ *

Proof 3.4**.**

Since at each scale s, $A_{s}f$ is generated as projection of $f$ on space

[TABLE]

Therefore using the notation $E^{s}=f-A_{s}f$, we write this update in an iterative form

[TABLE]

From the definition of $E^{s}|_{X}$ it follows that any convergent algorithm requires $\alpha_{s}\geq 0$. Rearranging (16) and taking an inner product with respect to $E^{s}|_{X}$

[TABLE]

Using the Cauchy Schwarz and Triangle inequalities

[TABLE]

Therefore

[TABLE]

Now, for convergent approximations $||E^{s+1}|_{X}||\leq||E^{s}|_{X}||$. Therefore for all scales, $\alpha_{s}\in[0,2]$. The end members $\alpha_{s}=\{0,2\}$ imply the iterative scheme has converged and no improvement is needed (as reflected in the stopping criterion in Algorithm 1). Hence, we remove these end members obtaining our desired bounds. Rewriting equation (16)

[TABLE]

[TABLE]

Now we analyze the inner product of the error of the projection ($E^{s}=f-A_{s}f$) with respect to projection at scale $s$ in the native RKHS. The major power of the following theorem lies in the fact that as the proposed algorithm converges, we are able to upper bound this inner product as a function of the user defined tolerance ($TOL$). This provides the user direct control on the approximation properties on the algorithm even in the native RKHS. It should be noted again that at each scale $s$, it is assumed that the function $f$ to be approximated belongs to the same Hilbert space $\mathcal{H}_{s}$ and is approximated in $\Gamma^{s}$, justifying the use of the reproducing property.

Theorem 3.5**.**

For observed data $f|_{X}=(f(x_{1}),f(x_{2}),...f(x_{n}))^{T}\in\mathbb{R}^{n}$ on $X\subset\Omega$, the approximation produced by the $A_{s}f$ and its corresponding error $E^{s}$ satisfies

[TABLE]

*where C_{S_{a}}=\Big{(}|c_{1}|,|c_{2}|,...|c_{l_{S_{a}}}|\Big{)}^{T}\in\mathbb{R}^{l_{S_{a}}} contains the modulus of coordinates of projection at $S_{a}$, n is the number of observations and TOL is the $2$-$norm$ convergence error tolerance in the hierarchical algorithm. *

Proof 3.6**.**

Beginning with the Inner products in the RKHS \Big{|}\Big{<}(A_{s}f),E^{s}\Big{>}_{\mathcal{H}_{s}}\Big{|}

[TABLE]

Here $||C_{s}||_{\infty}$ is the $\infty$ norm of $C_{s}=[|c_{1}|,|c_{2}|,...,|c_{l_{s}}|]^{T}$.

[TABLE]

*However, from the convergence criteria in the proposed Algorithm, we know $||E^{s}|_{X}||_{2}\leq TOL$ as $s\to S_{a}$. Hence the theorem follows *

4 Approximation Properties and Confidence Intervals

This section provides estimates quantifying the quality of approximations generated by the proposed algorithm at each scale. The first result is a direct application of theorem (3.5)

4.1 Approximation Properties

Corollary 4.1**.**

For observed data $f|_{X}=(f(x_{1}),f(x_{2}),...f(x_{n}))^{T}\in\mathbb{R}^{n}$ on $X\in\Omega$, the approximation $(A_{S_{a}}f)$ produced by the hierarchical algorithm at the convergence scale ($S_{a}$), follows the Pythagoras Theorem in the limit $TOL\to 0$. i.e,

[TABLE]

Proof 4.2**.**

Let $f$ be the function to be approximated. As stated earlier, it was discretely observed at $X\in\Omega$ leading to the restricted function $f|_{X}$. Starting with the norm of the function in the Hilbert Space

[TABLE]

*Using theorem 3.5, at $s=S_{a}$ in the limit $TOL\to 0$, the result follows *

Now, we will provide results related to uniqueness and quality of solution.

Theorem 4.3**.**

For observed data $f|_{X}=(f(x_{1}),f(x_{2}),...f(x_{n}))^{T}\in\mathbb{R}^{n}$ on $X\in\Omega$, the approximation $(A_{s}f)\in\Gamma^{s}$ in the limit $TOL\to 0$ produced by the proposed algorithm at the convergence scale $s=S_{a}$ is

the unique orthogonal projection to $f$ 2. 2.

the unique best approximation to $f$ with respect to $||\cdot||_{\mathcal{H}_{s}}$

Proof 4.4**.**

(1): $A_{s}f$ is a unique orthogonal projection of $f$ on $\Gamma^{s}$ if $(f-A_{s}f)\perp\Gamma^{s}$. Again using the reproducing property of RKHS

[TABLE]

The rest of the proof is similar to steps in the proof of theorem 3.5 showing the considered native Hilbert space inner product vanishes as $TOL\to 0$

(2): Considering any $d_{s}\in\Gamma^{s}$. Then at $s=S_{a}$, $<d_{s}-A_{s}f,f-A_{s}f>_{\mathcal{H}_{s}}\to 0$ as $TOL\to 0$ (by part 1) Therefore,

[TABLE]

[TABLE]

*which establishes the optimality of the approximation *

Now, we will analyze the error functional and pointwise error bounds. At scale $s$, we are searching for a solution in the space $\Gamma^{s}$. Our approximation is of the form

[TABLE]

Any $A_{s}f\in\Gamma^{s}$ can also be expressed as

[TABLE]

$\delta_{x}$ is the evaluational functional for $f$, i.e. $\delta_{x}(f)=f(x)$ and $\lambda^{y}$ is application of linear functional $\lambda$ on $y$. Thus, we have the dual space \mathcal{H}^{*}_{s}=\Bigg{\{}\sum_{x_{j}\in X_{s}}c_{j}\delta_{x_{j}}\Bigg{\}}.

Next, we present a result which provides an inner product representation for approximation at a point $x\in\Omega$. The objective here is to show the optimality of the approximation $A_{s}f$ by requiring vanishing gradient for the norm of the Error Functional in $\mathcal{H}_{s}$.

Theorem 4.5**.**

If the error functional at scale $s$ is represented as $\varepsilon^{s}_{x}=\delta^{s}_{x}-M_{s}(x)^{T}\delta^{s}_{X}$ $\ni$

[TABLE]

Then the optimal value of $M_{s}(x)=[M^{1}_{s}(x),M^{2}_{s}(x),M^{3}_{s}(x)....,M^{n}_{s}(x)]^{T}\in\mathbb{R}^{n}$ which minimizes the error functional norm $||\varepsilon^{s}_{x}||^{2}_{{\mathcal{H}}_{s}}$ satisfies the inner product

[TABLE]

*establishing the optimality of approximation $(A_{s}f)$ at each scale s. *

Proof 4.6**.**

We begin by expressing the error functional norm

[TABLE]

Now, based on the property of dual space, we know at scale s,

[TABLE]

Also, let

[TABLE]

[TABLE]

Therefore, we have

[TABLE]

Differentiation yields the optimal $\hat{M}_{s}(x)$ that minimizes $||\varepsilon^{s}_{x}||^{2}_{\mathcal{H}_{s}}$, as the the solution of

[TABLE]

From Algorithm 1, we know $G_{s}$ need not be full rank. Using matrix A (as in Algorithm 1), constructing $W(=AG_{s})$ and a Column pivoted QR decomposition $WP_{s}=QR$, and, applying the permutation operator $P_{s}$ on equation (24) yields

[TABLE]

Since $G_{s}$ is a symmetric operator, sampling the first $l_{s}$ rows of $P_{s}G_{s}$, where $l_{s}$ represents the numerical rank of $G_{s}$, will produce ${B^{s}}^{T}$, the transpose of the bases considered at scale $s$. Correspondingly, sampling the respective values in $R_{s}(x)$ produces the restriction of $R_{s}(x)$ to set $X_{s}$ (represented as $R_{s}(x)|_{X_{s}}$). Therefore, restricting system (25) to equations only corresponding to $X_{s}$

[TABLE]

where we use the Moore-Penrose inverse for getting the optimal projection of $R_{s}(x)|_{X_{s}}$. Thus, optimal solution for $M_{s}(x)$ is given by

[TABLE]

Now computing the inner product $<f|_{X},\hat{M}_{s}(x)>$

[TABLE]

Now, we provide approximation results for optimal point evaluational functional. For the approximation at each scale, we can obtain the variance associated with that prediction by minimizing the squared error in evaluation at $x$ in the native reconstruction space. However, since there is no noise or uncertainty associated with the data, this variance will be zero.

Theorem 4.7**.**

For observed data $f|_{X}=(f(x_{1}),f(x_{2}),...f(x_{n}))^{T}\in\mathbb{R}^{n}$ on $X\subset\Omega$, there exists at each scale, a point approximation functional ($\psi^{s}_{x}$) for the optimal approximation $A_{s}f$ such that $\psi^{s}_{x}(f)=(A_{s}f)(x)$. This functional has a representation

[TABLE]

Also, it produces the minimized error functional in the sense

[TABLE]

Proof 4.8**.**

We begin with the expression for the norm of the error functional in equation (23)

[TABLE]

since, it is minimized when $M_{s}(x)$ is given by equation (27). Now, writing the optimal approximation functional at x as:

[TABLE]

which completes the proof for equation (28). Now, considering the error functional at x and putting the optimal $M_{s}(x)$ in equation (23), we get

[TABLE]

Recognizing $\psi^{s}_{x}$ from equation (28) and substituting in the above equation.

[TABLE]

*Which completes the proof for equation (29). *

Before moving forward, here we also provide a result for bounding the Error functional.

Corollary 4.9**.**

For the error functional defined in Theorem 4.5, the absolute error at x is bounded in the sense

[TABLE]

Proof 4.10**.**

Using the relation,

[TABLE]

again using equation (23)

[TABLE]

Now, since (23) is minimized when we satisfy (24). Also since, $||\delta^{s}_{x}||^{2}_{{\mathcal{H}}_{s}}=1$. Putting the values in (23) we get

[TABLE]

*Therefore result in equation (30) follows. *

4.2 Confidence and Prediction Intervals

Confidence intervals in general are a measure of our belief in the estimated approximation. Prediction intervals on the other hand refer to the bounds around the mean fit, where a future datapoint is expected to fall. This is crucial information in conjunction with the sparse representation, as even when we are not able to capture the function accurately at initial scales, we have an estimate of the expected behavior of the observations. Algorithm 1 makes predictions at each scale based on the corresponding sparse representation ($D^{s}_{sparse}$). Here if the error in approximation is greater than user defined tolerance, i.e. $||f|_{X}-(A_{s}f)|_{X}||_{2}\geq TOL$, then it signifies that Approximation $A_{s}f$ needs additional degrees of freedom to capture the underlying data generation process (as in principle it is the best possible approximation, given a basis set). For modeling this error value, we consider a model of form

[TABLE]

Therefore the sampling distribution for $f|_{X}$ would be given as

[TABLE]

Since, we know that the projection coordinate is given as $C_{s}=({B^{s}}^{T}B^{s})^{-1}{B^{s}}^{T}(f|_{X})$,

[TABLE]

[TABLE]

which shows an unbiased estimator. Hence $C_{s}$ is the best unbiased approximation at scale s given bases $B^{s}$. Now for computing the distribution of response $f$ at some $x_{*}\in\Omega$

[TABLE]

[TABLE]

For an estimated value of $\sigma$ ($\hat{\sigma}$), we can write the standard deviation of prediction at $x_{*}$

[TABLE]

Therefore if we use t-distribution for confidence bounds, we get the following $100(1-\alpha)\%$ confidence intervals for $E[A_{s}f(x_{*})]$

[TABLE]

Here $l_{s}$ is again the numerical rank of $G_{s}$ and $n$ is the original number of observations made. Now for prediction interval, we know, $var(f(x_{*})-(A_{s}f)(x_{*}))=var(\epsilon^{s})+var(A_{s}f(x_{*}))$, therefore,

[TABLE]

Hence the $100(1-\alpha)\%$ prediction intervals are given as

[TABLE]

For the estimated value of $\sigma^{2}$, we use its unbiased estimation at scale s given as

[TABLE]

5 Stability Properties

In this section we provide bounds related to stability of approximations obtained by the proposed algorithm. The first result bounds the approximation at scale $s$ with respect to the $L_{\infty}$ topology for some compact domain $\Omega\in\mathbb{R}^{d}$.

Theorem 5.1**.**

For observed data $f|_{X}=(f(x_{1}),f(x_{2}),...f(x_{n}))^{T}\in\mathbb{R}^{n}$ on $X\in\Omega$, the approximation $A_{s}f$ at any scale s is bounded in the $L_{\infty}(\Omega)$ norm as

[TABLE]

[TABLE]

where $\sigma_{smax}$ denotes the largest singular value and $D_{s}$ is obtained by applying the extension operator ${{B^{s}}^{\dagger}}^{T}$ on $G_{X_{s}}=R_{s}(x)|_{X_{s}}R_{s}(x)|^{T}_{X_{s}}$, i.e. $D_{s}={{B^{s}}^{\dagger}}^{T}G_{X_{s}}{B^{s}}^{\dagger}$

Proof 5.2**.**

Expressing $(A_{s}f)$ in terms of the inner product as in equation (22).

[TABLE]

Let $x^{*}\in\Omega$ be the data point at which $\sum_{j=1}^{n}|M^{j}_{s}(x)|$ is maximized. Therefore

[TABLE]

Now, we know

[TABLE]

Therefore

[TABLE]

For lower bound on $P^{s}_{\infty}$,

[TABLE]

For a tall thin matrix B, ${{B^{\dagger}}^{T}}^{\dagger}=B^{T}$. Therefore, $M_{s}(x)$ is the least squares solution of ${B^{s}}^{T}M_{s}(x)=R_{s}(x)|_{X_{s}}$, and,

[TABLE]

[TABLE]

where $\sigma_{smax}$ is the maximum singular value. Putting in equation (40).

[TABLE]

*Combining equation (39) and above result, the bounds on $P^{s}_{\infty}$ follow *

We will conclude this section by providing a bound on the approximation at any scale s at some point $x\in\Omega$

Theorem 5.3**.**

For observed data $f|_{X}=(f(x_{1}),f(x_{2}),...f(x_{n}))^{T}\in\mathbb{R}^{n}$ on $X\in\Omega$, the absolute value of the approximation produced by the hierarchical algorithm $(A_{s}f)$ at scale $s\leq S_{a}<S_{c}$ at any point $x\in\Omega$ is bounded as

[TABLE]

Where $D_{s}$ is defined as in Theorem (5.1). Also if the convergence happens at the critical scale ($S_{c}$), then at convergence, bound (42) can be simplified as

[TABLE]

*Where $\sigma_{max}$ is the maximum eigenvalue operator *

Proof 5.4**.**

Here again we begin with the absolute value of approximation produced by the proposed algorithm at scale s.

[TABLE]

The last step was carried out using equation (38). Now if the convergence happens at $s=S_{c}=S_{a}$, then $D_{S_{c}}=G_{S_{c}}^{-1}$ and therefore $D_{S_{c}}(j,j)=G_{S_{c}}^{-1}(j,j)$. Let $g(j,j)=G_{S_{c}}^{-1}(j,j)$

[TABLE]

Remark 5.5**.**

*We note that the bounds here are conservative and depend on the data set size $n$. The assumed global overlap of the basis functions leads to the loose upper bound. However, as Figure 4 (in the next section) shows the basis functions have a rapid decay and attained bounds in practice are much smaller. *

6 Results and Analysis

This section analyzes the behavior of the proposed approach on a variety of datasets under different conditions. The first subsection here studies the performance on synthetic datasets. This is important, as here we know the truth and hence quantification of performance becomes feasible. The following subsections deal with application on real datasets. Here we take 2 different applications which test the performance of the proposed hierarchical algorithm.

6.1 Analysis on Synthetic Datasets

Here we have chosen a set of 4 test functions (Figure 1) from literature [23] providing our proposed algorithm, the sampled data to learn the underlying function. These test functions have been shown in Figure 1, and are mathematically expressed as:

• •

Test Function 1 (TF1): Gramacy and Lee Test function

[TABLE]

• •

Test Function 2 (TF2): 1-D Schwefel Function

[TABLE]

• •

Test Function 3 (TF3): Dropwave Function

[TABLE]

• •

Test Function 4 (TF4): 2-D Schwefel Function

[TABLE]

After sampling the data from these test functions, for all axis (X and Y for TF1 and TF2, and X, Y and Z for TF3 and TF4 respectively) the values are normalized between -1 to 1 by dividing the measurements with the corresponding absolute maximum value along each axis. It should be noted here that for most of the analysis presented here, we have sampled 200 equidistant points for TF1 and TF2. And for TF3 and TF4, points are sampled on a $50\times 50$ grid. The motivation here is to sample data points from the test functions and reconstruct the functions from the sampled data using the generated sparse representation ($D^{s}_{sparse}$). These functions were specifically chosen as they have a lot of curvature changes and multiple local minima and maxima, which makes learning the function form difficult.

We begin with the analysis of the convergence of the algorithm on the test functions.

6.1.1 Convergence Behavior

Figure 2 shows the convergence behavior of the Algorithm 1 by studying the $2$-$norm$ error of the prediction with respect to the original observations. $TOL=10^{-2}$ was used in Algorithm 1 for generating these results. Following the notations used earlier, $S_{a}$ here represents the convergence scale with $S_{c}$ being the critical scale. Figure 2 also illustrates the proportion of dataset used at each scale $s$ for generating the approximation $A_{s}f$. It is essentially the proportion of the dataset used as the sparse representation. Therefore, from figure 2 we can make inferences like, for test function 1, at scale 6 with 23% of datapoints, the proposed algorithm was able to generate an approximation $A_{6}f$ which had a $2$-$norm$ error of less than $10^{-5}$. It should be noted here, that based on the structure, the error measure reduces in a unique manner for all 4 test functions which have different convergence scales.

Figure 3 shows the convergence bounds from Theorem 3.5. Here we have shown the results in $1$-$norm$, which state that at any scale s, in the Hilbert space

[TABLE]

Figure 3 plots the quantity on right in the above equation with increasing scales. The sharp drop in the this bound for all 4 test functions justifies the capacity of the algorithm to produce good approximations. Precisely, at higher scale in the native Hilbert space, the error in prediction approaches orthogonality with respect to the approximation $A_{s}f$.

6.1.2 Scale dependent $D^{s}_{sparse}$ and corresponding reconstructions

Figure 4 and Figure 5 show the behavior of the selected sparse representation and the approximation $A_{s}f$ produced with increasing scales. In this subsection we have only presented results for TF2 and TF3 for analyzing the performance of the algorithm in univariate and multivariate setting respectively. Starting with figure 4, the plots on the top row show that with increasing scale, the support of the basis functions becomes narrower. This corresponds to the increasing numerical rank of the kernel matrix $G_{s}$ with increasing $s$. The wider support of the basis functions in the initial scales also explains the corresponding over-smoothed approximations. This scenario is similar to behavior of approximation strategies with global bases (for example polynomial based approximation). In the bottom row, for every scale we have mentioned the number of points chosen as sparse representation (out of 200 points). It should be noted here that the blue points represent the $D^{s}_{sparse}$ sampled from the smaller green data points. Here, our motivation is not just to show that with few points, the algorithm is able to learn the underlying function. But also, that the algorithm has an inbuilt capacity to choose a small set of representative points which can appropriately capture the function structure.

Figure 5 shows the corresponding result for the 2-D wave functions. Since proper visualization of the basis functions for a surface is a little challenging, so here we have just shown the location of the points chosen in the sparse representation (in X-Y plane). One additional thing to be noted here is the higher density of sampling near the edges of the domain. This directly corresponds with the fact that at the edges, for matching the curvature appropriately, it needs more points as there is no scope of learning beyond the edges. The corresponding reconstruction also shows how the specific features are learned over the scales. Again, since points were sampled on a $50\times 50$ grid. So $D^{s}_{sparse}$ consists of data points sampled from 2500 design points.

6.1.3 Confidence and Prediction Intervals

The results for this section have been shown in figure 6. The analysis is shown for scale 0 to scale 5. The green points are the original data and the red points show the sampled ones for $D^{s}_{sparse}$. The thinner (bluish) bands show the 95% confidence interval on the estimated approximation at each scale. It should be noted here that, we have not carried out a full Bayesian analysis here. Instead, we have just used the fitting variance (equation 34) as a proxy for the variance and used it for scaling our interval. Specifically here we are using t-confidence bands which are suitable for smaller datasets (as compared to Gaussian bounds) and tend to be normal in the limit of larger datasets. The thicker band show the prediction interval. The main idea to be conveyed here is that if at a scale s, we have our sparse representation $D^{s}_{sparse}$, then along with these prediction intervals, we can make estimations of where all the deleted points and data points to be sampled in future would lie. One other way to say the same thing is, that if we are at a particular scale $s$, and if we fix the location on $x\in X$, then we can be 95% confident that the mean of $y$ values observed at that particular x will lie within the blue bounds. Similarly for a fixed $x$, the broader salmon color bands show the range in which any new observation to be made in the future will lie.

The results here also confirms the fact that with increasing scales both the bounds become very thin showing confidence in the approximation produced.

6.1.4 Importance metric for design points

This section aims at further exploring the application of the proposed algorithm. The results of the current study are presented in figure 7 and 8. The idea is based upon the requirement that besides just getting the sparse representation at each scale, sometimes we also need to define an ordering of data points in $D^{s}_{sparse}$ in decreasing order of importance with respect to efficient reconstruction of $f$. This is important because if a measurement at a design point is of very high priority, then more resources could be engaged to measure that particular observation accurately. Also if we need to further compress the data, then in what sequence the datapoints can be removed. A test case is shown in figure 7. It contains the importance ranking computed based on the order in which the points were sampled while getting the sparse representation. It shows that the most important point (rank 1) at scale 0 is near the center of the function (and hence should rank last in the ordering for deletion). This confirms our belief that an observation near the center is crucial to capture the behavior of the function to be approximated. The second and third most important points (rank 2 and 3 respectively) are found to be at the very end of the curve which again is logical based on the fact that algorithm needs precise information to capture the function at the edges. These results also make sense because of the inherent symmetry of the 1-D schwefel function under consideration. It should be noted that the location of these important points depends on the nature of the function under consideration, so for any other functions, the location of important points might be very different as compared to the ones obtained in figure 7.

In order to get a more detailed picture of importance metric, we ran 1000 simulations of Algorithm 1 and presented the distribution (for location) of the top 3 most important points. The results are presented in figure 8. Here the X-axis varies from 0 to 200 because we considered 200 sampled points from TF2 as our learning set. The analysis was run for scales 0, 1, 2, 3, 5, 10 for studying the behavior of the distribution. It could be very well seen here that all the weight of the distribution for the most important point is concentrated at the center of the domain for the initial scales (top row). For the same scales, the second (middle row) and third (bottom row) most important points have all their mass concentrated at the edges. However, if we move towards the critical scale (where all points are included in the sparse representation), all the points have approximately uniform importance density distribution. This is also expected as when all points are sampled, no single point is more important than the other.

6.1.5 For non-uniform sampling of data

We illustrate the performance of the proposed algorithm on non-uniformly sampled data (and fewer data points) in Figure 9. Here in plot (a), we can see the performance when the learning dataset composed of functional values at 40 randomly chosen locations. Here we have also shown the true function for visualizing the quality of the approximation generated. If we look closely at the reconstruction in figure 9(a) around $x=1.4$, then the ability of the proposed approach to capture the respective peak in the underlying function even though there was not enough data to reflect it is clearly visible. Besides 40, the reconstruction has also been shown (in figure 9) for datasets consisting of 50, 60 and 70 points (b, c and d respectively) with promising results.

Before moving forward with application on real datasets, the next subsection provides a comparison of the performance of Algorithm 1 with algorithm 4 in [3].

6.1.6 Comparison with Algorithm 4 in Bermanis et. al

As mentioned earlier, [3] was one of the major motivators for the current work. In that paper algorithm 4 for multiscale data sampling and function extension comes very close in behavior to our hierarchical approach. Briefly the idea there can be summarized as follows. Suppose the approximation at scale $s$ is represented as $H_{s}f$. Therefore starting with scale 0, $H_{0}f$ is the approximation to $f$ produced at scale [math]. Thus we can write

[TABLE]

However, when we move further, the error orthogonal to the search space at the previous scale becomes the target function for the next scale. Therefore,

[TABLE]

[TABLE]

[TABLE]

With a user defined error tolerance (err), the authors define the convergence scale ($s^{*}$) as the scale satisfying

[TABLE]

Therefore, final approximation to f is of the form

[TABLE]

Hence algorithm 4 in [3] has bases from all scales for the final approximation. If we take a closer look at our algorithm (Algorithm 1), we see that it also samples bases at each scale. However, it just uses the bases at convergence scale as the final bases. The behavior of the algorithm 4 from [3] is compared with Algorithm 1 in figure 10. Here the decay of error shows similar behavior for both the algorithms. For TF1, Algorithm 1 shows some faster convergence. However, for TF3 and TF4, the multiscale extension algorithm reduces the error to below machine precision faster than Algorithm 1. Although quite impressive, this doesn’t make Algorithm 4 from [3] any more useful because we are already at error levels of $10^{-13}$ at such higher scales.

With comparative learning behavior, we now move to the comparison of prediction capability. If we think of prediction at a new design point for the case of Algorithm 4 in [3], then we will have to keep track of points sampled at each scale from $s=0$ to $s=s^{*}$. Once we have that, we can combine the formulated bases of the prediction points points with respect to these points linearly using the projection coordinate at each scale as in equation (49). This is where our algorithm outperforms Algorithm 4 by only just requiring the bases formulated with respect to the sparse representation at the convergence scale. This characteristic of Algorithm 2.1 allows us to talk about sparse representation of the dataset which is not possible with the definition of Algorithm 4 in [3]. Figure 11 shows this more clearly. Here we have measured the time which each of the algorithms take for just prediction, with all the learning assumed to be performed beforehand. Here $|X|$ denotes the size of the data used for learning and $|X_{pred}|$ shows the number design points at which the prediction is to be made. Figure 11 visually shows the non-trivial improvement in the prediction latency by our hierarchical approach with repsect to the multiscale algorithm from [3].

6.2 Application on Real Data

In this section we have analyzed the performance of Algorithm 1 on datasets from some practical scenarios. Specifically we are dealing with two different datasets here. The first dataset is spatial in nature (figure 12 (a)) where the objective is to learn the sparse representation with the demonstration of capability to reconstruct the data from $D^{s}_{sparse}$. Here we discuss one additional application, that relates to improving the associated inherent spatial resolution.

The second category of dataset considered here comes from numerical modeling of gravity measurement changes (figure 12 (b)) observed over the Greenland Ice sheet. Here again we analyze the capacity of Algorithm 1 to construct sparse representation of the dataset and reconstruct the dataset from $D^{s}_{sparse}$.

6.2.1 Application on Spatial Dataset - generating $D^{s}_{sparse}$

For this paper we are considering a particular type of spatial dataset known as Digital Elevation Model (DEM). It is a topographical map of a particular region. The data is arranged on grids with each grid node $(x,y)$ associated with a height measurement. Here X and Y are projection coordinates on a horizontal plane (usually transformed from latitudes and longitudes). Here, we are considering the DEM dataset shown in figure 12 (a). The idea in this study is to generate a sparse representation of the DEM data and study the reconstructions produced by these representations as we move up the scale. Figure 13 show these results for scale 0,2 and 6. These specific scales were chosen so as to provide an idea of how the surface is evolving towards the starting scale and towards the end. The important thing to note here is that at scale 6, with $D^{s}_{sparse}$ only consisting of 763 points out of 4350 points, the algorithm was able to generate a reconstruction to the DEM where the prediction error in $\infty-norm$ is just 6.82 (compared to the range of variation observed as $\sim[250,500]$ in the colorbar in figure 12(a) ). The compression metric of ($763/4350\sim 17.5\%$) clearly shows the success of the algorithm in generating a sparse representation for the dataset. Here $\infty-norm$ was chosen as an error measure as it upper bounds the error at any individual point and gives an intuitive understanding of the performance of the algorithm.

6.2.2 Application on Spatial Dataset - Improving resolution

For many practical engineering problem like flow simulations, if the modeling region is relatively small, then the resolution of the DEM dataset for the topography plays a very crucial role in the accuracy and stability of the results. The resolution here is defined as the length of side of the square which determines one pixel value. This can also be stated as the length of a cell boundary in the gridded data. This problem of low resolution is also difficult to solve because even if we have DEM data of the same small region from a different source, it usually is from a different time epoch and so there is no guarantee that the topography would not have changed in this time duration.

The DEM dataset shown in figure 12(a) has a resolution of 60m. However for some particular analysis it might be required to approximate the topography at an even higher resolution. For this reason in this section we have presented the result obtained after interpolating the DEM to a 20m resolution. This result has been presented in figure 14(a). Although here we do not have a proper measure for the performance of the proposed algorithm for this particular task (as the ground truth is unknown), the improved clarity and smoothness of features from top panel to bottom panel in figure 14(a) shows promising nature of results.

6.3 Application on data from Numerical models

:

In this section we consider the output from a numerical model which has been widely used in the literature [25, 1] for determining the ice mass evolution of Antarctic and Greenland Ice sheets. [18] introduced a iterative strategy for generating the Gravity Recovery and Climate Experiment (GRACE) global solution of equal-area surface mass concentration parcels (also referred to as mascons) in equivalent height of water. Figure 12(b) shows the distribution of these mascons over Greenland ice sheet. The idea is to derive spatially and temporally distributed changes in the mass of land ice at 1 arc degree (approximately 100 km). These mascons are estimated directly from k-band range and range rate (KBRR) data for two co-orbitting satellites roughly 220 km apart. Here each of these approximately $100\times 100\ Km^{2}$ square regions have a time series for solution of each mascon in cm. equivalent of water height. For this study we are using the data product $V02.4$ (not corrected for glacial isostatic adjustments). In total, for the entire planet there are 41168 mascons divided broadly among land, ice and water. The time series associated with each of these mascons has 148 entries. We have assumed these observations to be associated with the middle of the mascon solution time window. However, for this study, we just use the spatial aspect of this dataset (just consider observations only at t=0 for all mascons) and showcase the capability of our approach to generate a sparse representation for this dataset. Figure 14(b) shows the sparse representation of the mascons at scale 0, 2, 4 and 5. In the bottom row, the corresponding reconstruction from the sparse representation has been shown as well. The main thing to notice here is that even at scale 0, with only a portion of the original mascons, our approach was able to learn the behavior of the dataset (this is evident from the comparison of reconstructions from the different scales)

7 Conclusion

In this paper, we have introduced a hierarchical method for learning a sequence of sparse representations for a large dataset. The hierarchy comes from approximation spaces defined based on a scale parameter. Principally, the proposed approach has been shown to be useful for data reduction applications coupled with learning a model for representing the data. The paper provides analysis that explains and studies the theoretical properties of the proposed approach. We derive bounds for stability, convergence and behavior of error functionals. In the results section, we have shown the performance of the approach as data reduction mechanism on both synthetic and real datasets (geo-spatial and numerical modeling datasets). The sparse model generated by the presented approach is also shown to efficiently reconstruct the data, while minimizing error in prediction.

Though the results shown in this paper depict the efficiency of the approach on a variety of datasets and settings, there are several areas in which the presented algorithm can be improved. Firstly, the implementation of the algorithm can be made more efficient by either optimizing the operations in the algorithm or by handling chunks of data at a time. Secondly, generation of sparse representation for noisy datasets poses another set of unique challenges which will be addressed in a companion paper under preparation. For this case, properly capturing the uncertainty in the generated sparse approximations becomes very crucial. Finally the hyperparameters like $P$ in (Algorithm 1) can be explored and studied further for making the approach behave better. This can be done by even utilizing additional information which is not directly observed but is inherently known (like the physics of a system). Hence, in essence the work presented addresses the need for efficient learning methods for large datasets and opens up new interesting approaches.