Random Forest regression for manifold-valued responses

TL;DR

This paper introduces a non-parametric Random Forest-based method for predicting complex manifold-valued data, addressing challenges in modeling non-linear geometric spaces common in biomedical and computer vision applications.

Contribution

It develops a versatile distance-based Random Forest regression approach applicable to manifold-valued responses without requiring explicit manifold knowledge.

Findings

Superior predictive performance in simulations

Effective in image completion tasks

Applicable even when manifold structure is unknown

Abstract

An increasing array of biomedical and computer vision applications requires the predictive modeling of complex data, for example images and shapes. The main challenge when predicting such objects lies in the fact that they do not comply to the assumptions of Euclidean geometry. Rather, they occupy non-linear spaces, a.k.a. manifolds, where it is difficult to define concepts such as coordinates, vectors and expected values. In this work, we construct a non-parametric predictive methodology for manifold-valued objects, based on a distance modification of the Random Forest algorithm. Our method is versatile and can be applied both in cases where the response space is a well-defined manifold, but also when such knowledge is not available. Model fitting and prediction phases only require the definition of a suitable distance function for the observed responses. We validate our methodology using simulations and apply it on a series of illustrative image completion applications, showcasing superior predictive performance, compared to various established regression methods.

Tables

Table 1. Table 1: Test errors from the digit completion application. Classification Error was not applicable for RF and dRF. The number and percentage of badly reconstructed images was visually ascertained from Fig. 3 .

	MSE	Clas. Error	Bad Rec. No	Bad Rec. $\%$
1NN	3.3665	0.165	40	20
RF	2.7651	-	54	27
kRF	3.2675	0.21	48	24
dRF	3.3287	-	37	18.5

Equations

$$\begin{array}{r}{G}_{\mathcal{M}}(S_w)=\frac{1}{2N_w}\sum _{{\mathbf{y}}_i\in S_w}\sum _{{\mathbf{y}}_i\in S_w}{d}^2({\mathbf{y}}_i,{\mathbf{y}}_j)-\frac{1}{2N_{wl}}\sum _{{\mathbf{y}}_i\in S_{wl}}\sum _{{\mathbf{y}}_i\in S_{wl}}{d}^2({\mathbf{y}}_i,{\mathbf{y}}_j)-\frac{1}{2N_{wr}}\sum _{{\mathbf{y}}_i\in S_{wr}}\sum _{{\mathbf{y}}_i\in S_{wr}}{d}^2({\mathbf{y}}_i,{\mathbf{y}}_j),\end{array}$$

$${\hat{z}}_{new,k}=\frac{1}{\sqrt{{\lambda }_k}}\sum _{i=1}^N{u}_{ik}{\hat{k}}_{MDS}({\mathbf{y}}_{new},{\mathbf{y}}_i),k=1,\dots ,m,$$

$$\begin{array}{rl}{\hat{k}}_{MDS}({\mathbf{y}}_{new},& {\mathbf{y}}_i)=-\frac{1}{2}{d}^2({\mathbf{y}}_{new},{\mathbf{y}}_i)+\frac{1}{2(N+1)}\left(\sum _{j=1}^N{d}^2({\mathbf{y}}_j,{\mathbf{y}}_{new})\right)\\ & +\frac{1}{2(N+1)}\left(\left(\sum _{j=1}^N{d}^2({\mathbf{y}}_j,{\mathbf{y}}_i)\right)+d^2({\mathbf{y}}_{new},{\mathbf{y}}_i)\right)\\ & -\frac{1}{2(N+1{)}^2}\left(\left(\sum _{j,l=1}^N{d}^2({\mathbf{y}}_j,{\mathbf{y}}_l)\right)+2\sum _{j=1}^N{d}^2({\mathbf{y}}_{new},{\mathbf{y}}_j)\right).\end{array}$$

$$k_G(r,\mathbf{v},t,\mathbf{w}),r,t\in T^q=\{1,\dots ,q\},\mathbf{v},\mathbf{w}\in \mathcal{V}.$$

$$\sum _{i=1}^N\sum _{j=1}^N{a}_i{a}_j\sum _{r=1}^q\sum _{t=1}^q{k}_G(r,{\mathbf{v}}_i,t,{\mathbf{v}}_j)\ge 0,$$

$$g_{r\mathbf{v}}(\mathbf{w})={\left(k_G(r,\mathbf{v},1,\mathbf{w}),\dots ,k_G(r,\mathbf{v},q,\mathbf{w})\right)}^T.$$

$${\hat{\mathbf{y}}}_{new}=\hat{g}({\hat{\mathbf{z}}}_{new})=\sum _{i=1}^N\sum _{r=1}^q{C}_{ir}{g}_{r{\mathbf{z}}_i}({\hat{\mathbf{z}}}_{new}),$$

$$vec(\mathbf{C})=({\mathbf{K}}^G+\frac{N}{{\gamma }_g}\mathbf{I}{)}^{-1}vec(\mathbf{Y}),$$

$$k_G(r,\mathbf{v},t,\mathbf{w})=\{\begin{array}{l}\exp \left(-\frac{\Vert \mathbf{v}-\mathbf{w}{\Vert }^2}{{\sigma }_G}\right),r,t=1,\dots ,q,r=t\\ 0,r\ne t\end{array},$$

Full text

1 Introduction

Predictive modeling is an integral part of data analysis. It encompasses the process of developing models which can accurately predict yet-to-be-seen data. A multitude of regression models exist that can be used for prediction of univariate, or multivariate vectorial, responses. Nevertheless, less work has been done on the difficult problem of modeling and predicting more complex objects that possess additional structure, be it morphological, directional, or otherwise [13].

Images [12], shapes [17], graphs [19], deformation tensors [21] are examples of complex data types that appear naturally as responses in image analysis, computer vision, medical imaging and other application areas. While such objects are typically represented as points on very high-dimensional spaces, they can be meaningfully represented as points lying on smooth non-linear hyper-surfaces of lower dimensionality, a.k.a manifolds [3]. Manifolds can be understood as surfaces that locally resemble the Euclidean plane, but have different global structures [14]. Ideally, a suitable predictive modeling methodology should work under the additional constrains imposed by the data’s inherent manifold structure instead of simply treating these complex objects as points on Euclidean spaces. Major difficulties of course arise by forgoing the assumption of a Euclidean space, such as lack of coordinates, vectors and no analytical definitions of expected values [6]. In this work, we are additionally interested in experimental settings where the input observations may be very-high dimensional, which poses further requirements in the construction of a predictive algorithm.

Three main families of methodologies concerning with regression for manifold-valued responses can be found in the literature; intrinsic manifold regression, kernel-based and Manifold Learning (ML)-based methods. Intrinsic manifold regression models are generalizations of linear regression on manifolds [6, 21]. They require the analytical definition of the data manifold, since they fit a parametrized curve on the data.The assumption of an underlying manifold drives the choice of geometric elements employed during model definition and parameter estimation. Unfortunately, this requirement can be rarely satisfied, either due to the nature of data or the inability to define a suitable manifold. Furthermore, this methodology is not ideal for regression analysis with highly dimensional input observations.

Kernel methods first appeared in the literature as methodologies tailored for complex input objects, such as trees or graphs [15]. The basic idea behind these methods is that, if the input data lie on a non-linear space, then they can be implicitly mapped on a very high (or infinite) dimensional inner-product space, in which standard regression methods can be applied [15]. This implicit mapping is achieved through a kernel function that defines the objects’ inner-product in that space [15]. Kernel methods for complex responses have also been proposed [19, 7], but suffer from a number of issues. First, depending on the data at hand, a meaningful kernel function must be found or constructed, a process that is not always intuitive. Second, another problem is the formulation of a prediction methodology, which has to be described as a kernel minimization problem over the response space. This is most commonly solved by reducing the search space to the training dataset.

The last family of manifold regression methods is based on non-linear dimensionality reduction, a.k.a Manifold Learning. Given a set of observed points, ML methods aim to project these onto a space of lower dimensionality, whilst retaining as well as possible the original geometrical structure of the data. ML-based methods first apply ML on the response data, and subsequently use standard regression models trained on the learned response embedding. A suitable methodology must be formulated to map predicted points from the embedding to the manifold, a process that is referred to as backscoring [20]. Appropriate selection of the ML technique is largely based on intuition and previous experience. A decrease in model fitting accuracy, compared to an intrinsic manifold model, is to be expected, since the response embedding is not guaranteed to completely capture the structure of the underlying manifold. Furthermore, ML techniques generate discrete and one way maps from the manifold to the embedding, and the existence of an inverse and continuous map cannot always be guaranteed [20].

Here, we present a new methodology for predictive modeling of manifold-valued responses, which addresses a number of key issues listed above. Our objective was to propose a unified framework for regression and prediction of complex objects that is accurate, computationally efficient and can be readily deployed in a variety of applications. In the training phase, we employ our modified Random Forest (RF) regression algorithm that can be trained using only pairwise distances between response observations [16]. The non-parametric RF methodology is efficient and can handle high dimensional input spaces. In contrast to intrinsic manifold regression, no analytical definition of the underlying response manifold is required, apart from the construction of an appropriate distance function for the responses. Comparing our distance-based approach to kernel methods, we identify further advantages. First, a vast library of readily available distance metrics exists in the literature (see for example [5]). Second, when the manifold metric is not known, distances can be intuitively approximated, using for example the Isomap distance formulation [18].

In the prediction phase, a response point estimate for a new input observation is found as follows. Pairwise distances between the unseen response and all training observations are estimated. Using this set of distances, the response is predicted on a Euclidean embedding computed using the ML technique of metric multi-dimensional scaling (MDS) and is finally mapped to the response manifold through a standardized backscoring procedure. Our prediction methodology follows a two-step approach akin to ML-based methods, whilst offering two additional advantages. First, our regression model is trained on the original manifold, which enhances the quality of the fitted model. Second, due to the fact that manifold distances are known, MDS is employed for ML and predictions can be analytically computed and backscored to the response manifold.

The rest of this article is organized as follows. In section 2 we present the details of our manifold regression methodology. Our simulation and application experiments are included in section 3. We conclude this work in section 4.

2 Random Forest Predictive Methodology for Manifold-Valued Objects

2.1 Problem Definition

Let $S$ be a dataset of $N$ observed input-response pairs $\{(\mathbf{x}_{i},\mathbf{y}_{i})\}_{i=1}^{N}$, with inputs $\mathbf{x}_{i}=(x_{i1},\ldots,x_{ip})\in\mathbb{R}^{p}$ corresponding to responses $\mathbf{y}_{i}\in\mathcal{M}\subset\mathbb{R}^{q}$. $\mathcal{M}$ is a possibly unknown manifold, equipped with a distance metric $d(\cdot,\cdot)$. We refer to $\mathbb{R}^{q}$ as the response representation space. Let $\mathbf{D}$ be the $N\times N$ matrix of pairwise distances between the observed response points, with elements $D_{ij}=d(\mathbf{y}_{i},\mathbf{y}_{j})$. See Fig. 1(a) for a graphical illustration of the described data. We want to construct a predictive methodology that leverages these distances in order to ensure that, for any given input $\mathbf{x}_{new}$, the predicted $\mathbf{\hat{y}}_{new}$ lies on $\mathcal{M}$.

2.2 Distance Random Forest Regression

We first concentrate on the construction of a suitable manifold regression learning algorithm. Random Forest (RF) is a non-parametric, non-linear regression and classification algorithm [2]. In more detail, RF is a collection of Classification and Regression Trees (CARTs). A CART splits the input space recursively, according to a predefined split function, to small rectangular regions and then fits a simple model, commonly a constant value, in each one of them. See [2] for a detailed description of RF. In [16], we presented a modified distance Random Forest (dRF) algorithm, where the split criterion was formulated to depend only on pairwise distances between responses:

[TABLE]

with $\mathbf{y}_{i}\in\mathcal{M}$. Here $N_{w},N_{wl},N_{wr}$ are the cardinalities of the data subsets $S_{w}$, $S_{wl}$, $S_{wr}$ belonging to a node $w$ and its children nodes (left and right), respectively. The objective (1) is a generalization of the cost function used in standard regression RF, decoupled from the use of Euclidean norms and means, dependent only on pairwise response distances. As such, it enables the RF algorithm to be applied in general metric spaces. Previously, we used dRF for regression applications with graph and covariance-based response objects [16]. dRF lends itself naturally to manifold-valued data, whether the manifold in question is analytically defined or implied by the use of a specific distance metric.

2.3 Predictive Methodology for dRF

Traditionally, when working with responses lying in Euclidean spaces, an RF prediction in made simply by averaging response points; a new input observation follows the split decision rules learned during training and reaches a terminal node -leaf- in the tree. It is then assigned a suitable value on the response space, most often the average response of that leaf’s responses in the training set. This approach though is not valid under the assumption of non-vectorial manifold responses [11].

Our proposed methodology uses dRF for learning a family of trees from the data, which only requires an appropriate distance metric for the responses, as described above. The prediction phase is substantially different from standard RF. For a new input point $\mathbf{x}_{new}$, we use the trained dRF model to predict all pairwise distances between the unknown response, $\mathbf{y}_{new}$, and all observed responses in the training dataset, $\mathbf{y}_{i}$. Having estimated these distances, we then utilize them to to predict the response on an Euclidean embedding of the manifold, learned from the observed response set using metric MDS. Subsequently, a backscoring method is employed to project the point back to the original response manifold $\mathcal{M}$. Figure 1 provides an illustration of the proposed methodology, and the details are in order.

2.3.1 RF-based Distance Prediction

The first step is to estimate the set of distances $\{\hat{d}(\mathbf{y}_{new},\mathbf{y}_{i})\}_{i=1}^{N}$. We exploit the inherent ability of the dRF model to provide a measure of similarity between pairs of response observations (see Fig. 1(b)). When the new input point $\mathbf{x}_{new}$ is ‘dropped’ through each tree in the forest, it reaches a leaf associated with a subset of the training data. The $\mathbf{y}_{new}$ can be considered similar to that leaf’s training responses. To compute the dRF-based similarities, a vector $\mathbf{a}=(a_{1},\ldots,a_{N})$ is used, with each element $a_{i}$ corresponding to the similarity between $\mathbf{y}_{new}$ and $\mathbf{y}_{i}$. Initially $a_{i}$ is set to zero, for all $i=1,\ldots,N$. For each tree in RF, $a_{i}$ is incremented by one, each time $\mathbf{x}_{new}$ ends in the same node as $\mathbf{x}_{i}$. Normalization of similarities is performed by division with the number of trees.

Based on the similarity vector $\mathbf{a}$ and the training distance matrix $\mathbf{D}$, the response distances are predicted using the algorithmic procedure 1 ( see Fig. 1(c) ), which guarantees that the new point will lie in close proximity to at least its closest neighbor based on the dRF affinities, and that predicted distance values respect the triangular inequality property of metric $d$. In detail, the algorithm initially identifies the closest training point to the new response, according to $\mathbf{a}$, and assigns the minimum distance observed in the training set as the distance between these two responses. Subsequently, we iterate over the remaining responses, in decreasing order of distance from the first point, and assign a value for the distance to the new observation as follows: For each triplet including the considered point, the new point and a point for which the distance to $\mathbf{y}_{new}$ has been already estimated, we keep the maximum of the two known distances. Subsequently, we assign the minimum of all identified maximums as the predicted distance.

2.3.2 Prediction on the Response Embedding

Given $\{\hat{d}(\mathbf{y}_{new},\mathbf{y}_{i})\}_{i=1}^{N}$, we can predict a point estimate of the response on a Euclidean embedding of the manifold using MDS [4] (see Fig. 1(d) ).

MDS computes an approximation $\{\mathbf{z}_{i}|\mathbf{z}_{i}\in\mathbb{R}^{m}\}_{i=1}^{N}$ of the manifold-valued dataset that resides on a $m$-dimensional Euclidean space, by minimizing a stress function of the form $\sum_{i=1}^{N}\sum_{j=1}^{N}\left(d^{2}(\mathbf{y}_{i},\mathbf{y}_{j})-\|\mathbf{z}_{i}-\mathbf{z}_{j}\|^{2}\right)$. Let $\mathbf{D}^{(2)}$ be the $N\times N$ matrix of squared manifold distances and $\mathbf{K}=-\frac{1}{2}\mathbf{H}\mathbf{D}^{(2)}\mathbf{H}$, where $\mathbf{H}=\mathbf{I}_{N}-\frac{1}{N}\mathbf{e}\mathbf{e}^{T}$, with $\mathbf{I}_{N}$ the $N\times N$ identity matrix and $\mathbf{e}$ an $N\times 1$ column vector of all ones. Individual elements of $\mathbf{K}$ are given by $K_{ij}=-\frac{1}{2}(D^{(2)}_{ij}-\frac{1}{n}S_{i}-\frac{1}{n}S_{j}+\frac{1}{n^{2}}S),i,j=1,\ldots,N,$ where $S_{i}=\sum_{j}D^{(2)}_{ij},S_{j}=\sum_{i}D^{(2)}_{ij},S=\sum_{i,j}D^{(2)}_{ij}$ are the $i^{th}$ row, $j^{th}$ column and overall element-wise sums of $\mathbf{D}^{(2)}$, respectively. In [9], it was shown that, if $\mathbf{K}$ has rank $p\enspace(p\leq N-1)$, with $\lambda_{1},\ldots,\lambda_{p}$ the $p$ ordered non-zero eigenvalues of $\mathbf{K}$ with corresponding eigenvectors $\mathbf{u}_{1},\ldots,\mathbf{u}_{p}$, and $\mathbf{u}_{l}=(u_{l1},\ldots,u_{lN})$ for $l=1,\ldots,p$, then the solution to the MDS problem, assuming $m\leq p$, is given by $z_{il}=\sqrt{\lambda_{l}}u_{li}$, with $l=1,\ldots,m,\enspace i=1,\ldots,N$, and $\mathbf{z}_{i}=(z_{i1},\ldots,z_{im})$.

The result of the MDS decomposition is a discrete and one-way mapping $\mathbf{z}_{i}=f_{MDS}(\mathbf{y}_{i})$ defined on the training data set. An Out-Of-Sample (OOS) method which allows mapping of a new manifold observation on the learned space of the embedding was constructed in [1]. Let $Y\in\mathcal{M}$ be a random variable defined on the manifold surface and $\mathbf{y}_{a},\mathbf{y}_{b}\in\mathcal{M}$ two observations of $Y$. A continuous kernel function $k_{MDS}$, which gives rise to $\mathbf{K}$ under the training observations, is defined as $k_{MDS}(\mathbf{y}_{a},\mathbf{y}_{b})=-\frac{1}{2}(d^{2}(\mathbf{y}_{a},\mathbf{y}_{b})-E[d^{2}(Y,\mathbf{y}_{a})]-E[d^{2}(Y,\mathbf{y}_{b})]+E[d^{2}(Y,Y)]),$ where $E(\cdot)$ denotes expectation. Then, an OOS prediction $\hat{\mathbf{z}}_{new}=(\hat{z}_{new,1},\ldots,\hat{z}_{new,m})$ is given by [1]:

[TABLE]

with $\hat{k}_{MDS}$ denoting the mean estimator of $k_{MDS}$ under the augmented dataset $\{\mathbf{y}_{i}\}_{i=1}^{N}\cup\{\mathbf{y}_{new}\}$:

[TABLE]

The OOS formula does not depend on the actual value of $\mathbf{y}_{new}$, but rather on the distances between the new response and points on the training dataset. We can thus leverage the predicted $\{\hat{d}(\mathbf{y}_{new},\mathbf{y}_{i})\}_{i=1}^{N}$ on equation (2) in order to get a point estimate of our response $\hat{\mathbf{z}}_{new}$ on the embedding space $\mathbb{R}^{m}$. Now we are left with the task of mapping $\hat{\mathbf{z}}_{new}$ to the original manifold ( see figure 1(e) ).

2.3.3 Mapping from the Embedding to the Manifold

The mapping of the response from the embedding to the manifold - backscoring - can be formulated as an interpolation problem. Specifically, we are looking for a smooth continuous function $g:\mathbb{R}^{m}\rightarrow\mathbb{R}^{q}$, that minimizes the cost function $\enspace\gamma_{g}\sum^{N}_{i=1}(g(\mathbf{z}_{i})-\mathbf{y}_{i})^{2}+\|g\|^{2}_{G}$. $\gamma_{g}$ is a weight parameter balancing the smoothness of the interpolation and the adherence to the data and $G$ is a space of smooth functions equipped with the norm $\|\cdot\|_{G}$, which will be constructed in the following.

A solution to the interpolation problem was presented in [20]. Let $\mathcal{V}$ be a closed subset of $\mathbb{R}^{m}$ and $k_{G}$ a kernel function of the form

[TABLE]

Furthermore, assume that $k_{G}$ is semi-positive definite on $(T^{q}\times\mathcal{V})\times(T^{q}\times\mathcal{V})$, with

[TABLE]

for any finite set of points $\{\mathbf{v}_{i}|\mathbf{v}_{i}\in\mathcal{V}\}_{i=1}^{N}$ and real numbers $\{a_{i}|a_{i}\in\mathbb{R}\}_{i=1}^{N}$. For a fixed $(r,\mathbf{v})$, equation (4) defines a function from $\mathbb{R}^{m}$ to $\mathbb{R}^{q}$ by the formula

[TABLE]

Based on the above, let $G$ be the space of all finite linear combinations of functions of the form (6), as $(r,\mathbf{v})$ varies in $T^{q}\times\mathcal{V}$, and its closure w.r.t the scalar inner product $\langle g_{r\mathbf{v}},g_{t\mathbf{w}}\rangle=k_{G}(r,\mathbf{v},t,\mathbf{w})$. It follows that $\|g_{r\mathbf{v}}\|_{G}=\sqrt{\langle g_{r\mathbf{v}},g_{r\mathbf{v}}\rangle}$.

The interpolation problem can be solved over the space of functions $G$ as follows [20]. Let $\mathbf{Y}$ be the $N\times q$ matrix of training responses with rows $Y_{i\cdot}=\mathbf{y}_{i}$, and $\mathbf{K}^{G}$ the $Nq\times Nq$ matrix of kernel values $K_{N(r-1)+i,N(t-1)+j}^{G}=k_{G}(r,\mathbf{z}_{i},t,\mathbf{z}_{j})$. The minimizing function can be estimated as:

[TABLE]

where $C_{ir}$ are elements of the $N\times q$ coefficient matrix $\mathbf{C}$ given by

[TABLE]

with $vec(\cdot)$ denotes the vectorization of a matrix into a column vector.

It is clear, from equations (7) and (8), that the estimation of $\hat{g}(\hat{\mathbf{z}}_{new})$ requires just the knowledge of pairwise kernel values between the $N+1$ points $\mathbf{z}_{1}.\ldots,\mathbf{z}_{N},\hat{\mathbf{z}}_{new}$.

In our studies, we opted to simplify the minimizer (7), in favor of having just one tuning parameter, by choosing

[TABLE]

where $\sigma_{G}$ is a free parameter adjusting the bandwidth of the kernels.

We notice that the backscoring formulation does not take into explicit consideration the response manifold and $\hat{g}(\hat{\mathbf{z}}_{new})$ is not guaranteed to lie exactly on $\mathcal{M}$. Nevertheless, we justify our choice by pointing that since $\hat{g}$ is a smooth interpolating function from the embedding to the training responses, predictions should also adhere well to the manifold.

3 Experiments

In this section, we present a comparative simulation study to assess the ability of our methodology to cast predictions that adhere to the response manifold. Subsequently, we use our method in two illustrative image completion applications, where the objective is to predict one half of an image from the other half.

3.1 Simulation Study - Prediction on a Swiss-roll Manifold

We simulated $N=900$ paired input-response points $\{(\mathbf{x}_{i},\mathbf{y}_{i})\}_{i=1}^{N}$, with inputs $\mathbf{x}_{i}\in\mathbb{R}^{6}$ corresponding to responses $\mathbf{y}_{i}\in\mathcal{M}_{sr}\subset\mathbb{R}^{3}$, where $\mathcal{M}_{sr}$ denotes the 2-dimensional swiss roll manifold embedded in $\mathbb{R}^{3}$. Only the first two input dimensions were built to be predictive of the output. In detail, we first sampled $900$ points $\{t_{i}\}_{i=1}^{N}$ from the uniform distribution $\mathcal{U}(\pi,3\pi)$ and $\{u_{i}\}_{i=1}^{N}$ from $\mathcal{U}(0,21)$. The response points on the swiss-roll were computed as $y_{i1}=t_{i}\cos t_{i},\enspace y_{i2}=u_{i},\enspace y_{i3}=t_{i}\sin t_{i}$, while the input variables $x_{i1}$ and $x_{i2}$ were computed by mapping $u_{i}$ and $t_{i}$ in the unit circle: $x_{i1}=\frac{t_{i}-\bar{t_{i}}}{\max t_{i}}\sqrt{1-\frac{1}{2}\left(\frac{u_{i}-\bar{u_{i}}}{\max u_{i}}\right)^{2}},\enspace x_{i2}=\frac{u_{i}-\bar{u_{i}}}{\max u_{i}}\sqrt{1-\frac{1}{2}\left(\frac{t_{i}-\bar{t_{i}}}{\max t_{i}}\right)^{2}}$ . The values for the remaining input coordinate variables were drawn from the standard normal Gaussian distribution. Gaussian Noise, with $\Sigma=0.5\mathbf{I}_{3}$, was added on the response points.

We compared our dRF prediction methodology to $k$NN regression and standard RF, which do not take into consideration the special form of the response space, as well as kernel RF (kRF) [7], an RF method that employs a kernel function to capture the structure of the response space during training. The simulated dataset was split into $S_{train}$, comprising of $600$ input-response pairs and $S_{test}$, consisting of the remaining $300$. For $k$NN regression, the value of $k$ was selected to be $5$. All RFs were built with 150 trees, no pruning and the number of candidate split features in each node was set to $3$. For kRF, the training gram matrix was calculated using the Gaussian kernel $g(\mathbf{y}_{i},\mathbf{y})=\exp{\left(-\frac{\|\mathbf{y}_{i}-\mathbf{y}_{j}\|^{2}}{2\sigma^{2}}\right)}$ with $\sigma=2.5$. We followed the prediction methodology as described in [7], with $\mathbf{y}_{new}=\operatorname*{arg\,min}_{\mathbf{y}\in S_{train}}\left(g(\mathbf{y},\mathbf{y})-2\sum_{i=1}^{N}a(\mathbf{x}_{new},\mathbf{x}_{i})g(\mathbf{y},\mathbf{y}_{i})\right)$, where $a(\cdot,\cdot)$ is the RF-based affinity. The minimization problem was solved over the training set.

For dRF, the backscoring parameters were $\sigma_{\mathcal{G}}=100$ and $\gamma_{\mathcal{G}}=200$. Since there is no analytical form for computing distances on a swiss-roll manifold, we approximated manifold distances using the Isomap distance formulation [18]: A neighborhood graph $G$ was constructed, in which each point $\mathbf{y}_{i}$ was connected to its $k=7$ nearest neighbors in $\mathbb{R}^{q}$ . Graph edges were assigned weights equal to the Euclidean distances between the connected points in $\mathbb{R}^{q}$. For any two points $\mathbf{y}_{i}$ and $\mathbf{y}_{j}$ in $S$, $d(\mathbf{y}_{i},\mathbf{y}_{j})$ was then estimated by the shortest path connecting $\mathbf{y}_{i}$ and $\mathbf{y}_{j}$ in graph $G$.

Figure 2 shows test error vectors for the various methodologies used, projected on the $y_{\cdot 1}y_{\cdot 3}$ plane. We can see that $5$NN missed the goal of regressing on the manifold. For standard RF, it is clear that the model constantly underestimated the radius of the predicted points around the $y_{\cdot 1}$ axis. This can be justified by considering that predictions are taken as average points on the euclidean space $\mathbb{R}^{3}$, unaware of a possible structure in the response space. kRF preserved the manifold structure of the predicted points better than RF, but suffered significantly from the added noise in the responses. dRF outperformed the other methods both in terms of compliance to the underlying response manifold and regarding good robustness to the addition of noise.

3.2 Applications on Image Completion

In imaging analysis, it is common to assume that a collection of similar images lies on a manifold [12]. This assumption is guided by the complex nature of images as data objects, as well as the observation that the Euclidean metric and the corresponding geometric structure that it imposes on the space do not bode well with the human perception of similarity and difference between images [12]. Here, we used our manifold regression methodology to predict the bottom half of handwritten digits and human faces from their upper half.

3.2.1 Handwritten digits

For this application, we extracted $1000$ gray-scale images of handwritten digits, from the UCI Machine Learning Repository [8]. Each digit class, from [math] to $9$, was represented in the dataset by $100$ $8\times 8$ pixel images. Input data were constructed by vectorization of the $8\times 4$ upper half pixel intensities. The dataset was split into training and testing subsets with $N_{train}=800$ and $N_{test}=200$. The upper part of each images was taken as input for the predictive models. Responses comprised of the images’ bottom parts. The test set can be seen in Fig. 3(a).

We compared predictions from $1$NN, RF, kRF and dRF models. All RFs were built with $300$ trees and $5$ candidate split features in each node. For dRF, the distance matrix was constructed using the Isomap distance, with the number of neighbors set to $5$. A $25$-dimensional embedding space was used and the backscoring parameters were $\sigma_{\mathcal{G}}=3$ and $\gamma_{\mathcal{G}}=20$. For kRF, the training gram matrix was calculated using the Gaussian kernel with $\sigma=5$. The reconstructed test digits for the various models are shown in Fig. 3(b)-(e).

Table 1 summarizes the prediction results for the test data. In the first column we include the test Euclidean Mean Squared Errors (EMSE) for all models. In the case of $1$NN and kRF, which draw predictions from the training dataset, we are also able to report misclassification errors, a.k.a the percentage of predicted lower parts that mismatched the ground truth, which are shown on the second column of Table 1. Finally, in the last two columns, we report on the number and percentage of badly reconstructed test images from visual inspection, based on the following criteria: blurriness of the reconstructed image, smooth transition between the upper and bottom image parts, correct digit reconstruction. This qualitative performance measure is important due to the non-Euclidean nature of the data, which, as will be discussed below, makes the EMSE unsuitable for judging the predictive performance.

1NN and kRF cast predictions drawn directly from the training images. As such, no blurriness existed on the reconstructed digits of Fig. 3(b) and 3(d). Bad reconstructions were either misclassifications or reconstructions where the transition between the upper (input) and lower (predicted) image parts was not smooth. Surprisingly, $1$NN outperformed kRF in terms of classification error, as well as upon visual inspection of the images.

It is obvious from Fig. 3(c) that standard RF is ill-suited for the specific application. The RF prediction approach of averaging pixel intensities from various images resulted in a high number of blurry and nonsensical digit reconstructions. It is important to notice that while RF performed the worst, based on visual assessment of the reconstructed images, it had the lowest test EMSE, as shown in table 1. This observation highlights the unsuitability of EMSE as a measure of performance, in the case of manifold-valued responses.

Reconstructed digits from our dRF model are shown in Fig. 3(e). Our predictions were not drawn directly from the training set. Nevertheless, we notice that the large majority of reconstructions had no fuzziness and the transition from upper to lower parts was smooth. For cases where some blurriness could be noticed in the reconstructions, its effect was significantly less severe than in the RF predictions, resulting in the overall lowest number of badly reconstructed digits, based on visual inspection.

3.2.2 Faces

In the second application, we used $400$ images of faces from the Olivetti dataset, as included in the scikit-learn python package [10]. The dataset consisted of ten gray-scale $64\times 64$ pixel images for each of 40 distinct subjects, with same subject images taken at different times and with varying pose and facial expressions. Input data were constructed by vectorization of the $64\times 32$ upper half pixel intensities. The dataset was split into training and testing subsets with $N_{train}=300$ and $N_{test}=100$. Images of the same subject were only allocated either to the training or the testing set. Again, the upper parts of the images were taken as inputs, while bottom parts as responses.

We compared results from $1$NN, linear Regression (LR), RF and dRF. RFs were built with $300$ trees and $8$ candidate split features in each node. For dRF, the training distance matrix was constructed using the Isomap metric with the number of neighbors set to $5$, a $25$-dimensional embedding space was used and the backscoring parameters were $\sigma_{\mathcal{G}}=9$ and $\gamma_{\mathcal{G}}=50$.

Four test images and their reconstructions for the various models are included in Fig. 4. As we noted in the previous application, EMSE does not provide a suitable descriptor of performance. We rely again on visual inspection of the reconstructed images. 1NN exhibited a hits and miss behavior, with some predictions being similar to the original face, whilst others, such as the second and third depicted faces, being completely different. In addition, there was minimal smoothness in the transitions from the input to the predicted parts, accentuated specifically on the nasal and zygomatic edges. LR predictions were extremely blurry. RF also suffered from a large amount of blurriness, although transitions between the two parts of the faces looked more natural. Finally, dRF reconstructions exhibited the best transition smoothness of all methods, with the predicted half-images being well aligned to their inputs. Although the predictions were not completely free of blur artifacts, the effect was less severe and facial details, such as nasolabial folds and lips, were clearly portrayed.

4 Discussion and Conclusion

In this work we presented a predictive modeling approach for response objects occupying non-linear manifold spaces. The regression methodology is based on a distance modification of the RF algorithm that we previously published [16], which decouples the model’s training from the problem of response representation. For prediction purposes, we constructed a framework in which point estimates are first predicted on a Euclidean embedding of the response manifold, learned from the training dataset, and then projected back on the original space.

Our methodology, in contrast to intrinsic manifold methods, necessitates just the definition of a meaningful distance metric in the response space. This can be especially useful for real-life applications, such as image analysis, where the underlying manifolds are usually not known. Our distance-based regression algorithm draws similarities to the family of kernel-based methodologies. One benefit over kernel methods is the vast library of readily available distance metrics for a plethora of data objects. Furthermore, our methodology presents a unified framework, which deals with backscoring to the original response space, an issue that a lot of presented kernel methods do not tackle sufficiently well.

In the performed experiments, our method showed superior predictive performance in comparison to various regression methods, whilst being able to handle high-dimensional inputs and noise on the response observations. In the future, we aim to investigate the problem of automatic estimation of the Euclidean embedding’s dimensionality, as well as the use of more elaborate kernel functions in the backscoring formulation.

Acknowledgments

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors. The article is currently under consideration for publication at Patter Recognition Letters.