On orthogonal projections for dimension reduction and applications in augmented target loss functions for learning problems

TL;DR

This paper explores the use of orthogonal projections in high-dimensional data for dimension reduction and introduces augmented target loss functions in deep learning, improving accuracy in clinical imaging and music classification.

Contribution

It analyzes the trade-offs in dimension reduction objectives and proposes a novel framework of augmented target loss functions for enhanced deep learning performance.

Findings

Projections often do not satisfy variance and distance preservation simultaneously.

Balanced projections improve classification results.

Augmented target loss functions increase accuracy in real-world tasks.

Abstract

The use of orthogonal projections on high-dimensional input and target data in learning frameworks is studied. First, we investigate the relations between two standard objectives in dimension reduction, preservation of variance and of pairwise relative distances. Investigations of their asymptotic correlation as well as numerical experiments show that a projection does usually not satisfy both objectives at once. In a standard classification problem we determine projections on the input data that balance the objectives and compare subsequent results. Next, we extend our application of orthogonal projections to deep learning tasks and introduce a general framework of augmented target loss functions. These loss functions integrate additional information via transformations and projections of the target data. In two supervised learning problems, clinical image segmentation and music information classification, the application of our proposed augmented target loss functions increase the accuracy.

Tables

Table 1. Table 4.1. Classification results of iris data, when using projected input data in support vector machine (SVM) and shallow neural networks (NN). Mean and Variance ( × 10 − 4 absent superscript 10 4 \times 10^{-4} ) of 1000 1000 1000 independent NN runs, 100 100 100 independent runs with 10 10 10 -fold cross-validation in SVM.

Input/Method	NN (10 HU)	NN (5 HU)	SVM
x	(97.6, 1.25)	(97.5, 2.29)	(96.7, 0.15)
p◇x	(98.4, 0.42)	(98.3, 2.15)	(97.3, 0.06)
p×x	(88, 1.73)	(87.9, 1.92)	(87.6, 0.63)
p${}_{\square }$x	(87.3, 9.56)	(86.9, 10.81)	(87.7, 0.42)
p${}_{○}$x	(96.8, 2.74)	(96.7, 1.77)	(96, 0.17)
p∗x (PCA)	(96.9, 1.36)	(96.5, 4.07)	(96, 0.37)
p⋆x	(62.1, 44.30)	(58.9, 70.78)	(56, 0.61)

Table 2. Table 6.1. Comparison of AUC values for photoreceptors segmentation and disruption detection.

Loss function	Photoreceptors	Disruptions
$L_{\mathrm{MSE}}$	0.9720	0.4399
$L_p$
$p=I_4$	0.9736	0.4686
$p_{{\lambda }_{2,4}}$	0.9746	0.4720
$p_{PCA}$	0.9716	0.5331
$p_{12}$	0.9755	0.5558

Table 3. Table 7.1. Classification results with different parameter choices. The standard inbuilt Tikhonov regularization ( ℓ 2 subscript ℓ 2 \ell_{2} -norm of θ 𝜃 \theta ) is weighted by β 𝛽 \beta . For α > 0 𝛼 0 \alpha>0 the feature transformations { T j } j = 1 16 superscript subscript subscript 𝑇 𝑗 𝑗 1 16 \{T_{j}\}_{j=1}^{16} are used in the loss function, either directly or weighted by a random projection p λ 6 , 16 subscript 𝑝 subscript 𝜆 6 16 p_{\lambda_{6,16}} . The accuracy of the model is measured by the number of correctly classified samples divided by the number of all samples.

$\alpha$	$\beta$	$p$	training	test data
0	0	-	0.5541	0.5716
0.01	0	$I_{16}$	0.5650	0.5683
0.01	0	$p_{{\lambda }_{6,16}}$	0.7722	0.7657
0	0.05	-	0.9771	0.9729
0.01	0.05	$I_{16}$	0.9849	0.9802
0.01	0.05	$p_{{\lambda }_{6,16}}$	0.9857	0.9833

Equations

$$G_{k,d}:=\{p\in {\mathbb{R}}^{d\times d}:p^2=p,p^{\top }=p,\mathrm{rank}(p)=k\}$$

$$\{\bar{x}+p(x_i-\bar{x}){\}}_{i=1}^m\subset \bar{x}+V,$$

$$\mathrm{tvar}(x):=\frac{1}{m-1}\sum _{i=1}^m\Vert x_i-\bar{x}{\Vert }^2.$$

$$\mathrm{tvar}(px)\le \mathrm{tvar}(x)$$

$$\mathrm{tvar}(x)$$

$$\mathrm{tvar}(px)$$

$$\frac{4\log (m)}{{ϵ}^2/2-{ϵ}^3/3}\le k,$$

$$(1-ϵ){\Vert x_i-x_j\Vert }^2\le \frac{d}{k}{\Vert p(x_i)-p(x_j)\Vert }^2\le (1+ϵ){\Vert x_i-x_j\Vert }^2$$

$$\Big{\{}\frac{d}{k}\frac{\|p(x_{i})-p(x_{j})\|^{2}}{\|x_{i}-x_{j}\|^{2}}:i<j\Big{\}}$$

$$\frac{(2+\tau )2\log (m)}{{ϵ}^2/2-{ϵ}^3/3}\le k,$$

$$\Big{\{}\frac{d}{k}\frac{\|P(x_{i})-P(x_{j})\|^{2}}{\|x_{i}-x_{j}\|^{2}}:i<j\Big{\}}\in[1-\epsilon,1+\epsilon]$$

$$(1-ϵ)\frac{k}{d}\mathrm{tvar}(x)\le \mathrm{tvar}(px)\le (1+ϵ)\frac{k}{d}\mathrm{tvar}(x),$$

$$\mathcal{M}(p,x)$$

$$\mathcal{V}(p,x)$$

$${\int }_{G_{k,d}}f(p)\mathrm{d}{\lambda }_{k,d}(p)={\int }_{G_{k,d}}f(p)\mathrm{d}\lambda (p),\text{for all }f\in {\mathrm{Pol}}_t({\mathbb{R}}^{d^2}),$$

$$\mathbb{E}\mathrm{tvar}(Px)=\frac{k}{d}\mathrm{tvar}(x),$$

$$\mathbb{E}\mathcal{M}(P,x)=1,$$

$$\mathbb{E}\mathcal{V}(P,x)=a_{k,d}(1-\frac{4}{m^2(m-1{)}^2}\sum _{\begin{array}{c}i<j\\ l<r\end{array}}⟨\frac{x_i-x_j}{\Vert x_i-x_j\Vert },\frac{x_l-x_r}{\Vert x_l-x_r\Vert }{⟩}^2),$$

$$\mathbb{E}\mathcal{V}(P,x)\le a_{k,d}⟶\frac{2}{k},\text{for}d\to \infty .$$

$$\mathrm{Corr}(\mathcal{M}(P,x),\mathrm{tvar}(Px))=\frac{\mathrm{Cov}(\mathcal{M}(P,x),\mathrm{tvar}(Px))}{\sqrt{\mathrm{Var}(\mathcal{M}(P,x))}\sqrt{\mathrm{Var}(\mathrm{tvar}(Px))}}.$$

$$\frac{{\min }_{i\ne j}{\Vert x_i-x_j\Vert }^2}{{\max }_{i\ne j}{\Vert x_i-x_j\Vert }^2}-\frac{m(m-1)}{2d}\cdot \frac{{\max }_{i\ne j}{\Vert x_i-x_j\Vert }^2}{{\min }_{i\ne j}{\Vert x_i-x_j\Vert }^2}.$$

$$\varrho (\{p_l{\}}_{l=1}^n):=\underset{p\in G_{k,d}}{\sup }\underset{1\le l\le n}{\min }\Vert p-p_l{\Vert }_F,$$

$$\frac{4\log (m)}{{ϵ}^2/2-{ϵ}^3/3}\le k\le d,$$

$$(1-\delta ){\Vert x_i-x_j\Vert }^2\le \frac{d}{k}{\Vert p_{l_0}(x_i)-p_{l_0}(x_j)\Vert }^2\le (1+\delta ){\Vert x_i-x_j\Vert }^2,i<j,$$

$$\Vert p_{l_0}x-px\Vert \le \Vert p_{l_0}-p{\Vert }_F\Vert x\Vert \le \varrho \Vert x\Vert ,x\in {\mathbb{R}}^d.$$

$${\varrho }_i\asymp t_i^{-1},$$

$$n_i\asymp t_i^{k(d-k)},i=1,2,\dots$$

$${\varrho }_i\asymp n_i^{-\frac{1}{k(d-k)}}.$$

$$\{\mathrm{tvar}(p_{1}x),\dots ,\mathrm{tvar}(p_{n}x)\},\{\mathcal{M}(p_1,x),\dots ,\mathcal{M}(p_n,x)\},$$

$$\mathrm{tvar}(p_{l}x)\approx s\cdot \mathcal{M}(p_l,x)+\gamma ,l=1,\dots ,n,$$

$$s$$

$$\gamma$$

$${\mathcal{V}}_{k,d}:=\{q\in {\mathbb{R}}^{k\times d}:q{q}^{\top }=I_k\}$$

Full text

1. Introduction

Linear dimension reduction is commonly used for preprocessing of high-dimensional data in complicated learning frameworks to compress and weight important data features. In contrast to nonlinear approaches, the use of orthogonal projections is computationally cheap, since it corresponds to a simple matrix multiplication. Conventional approaches apply specific projections that preserve essential information and complexity within a more compact representation. The projector is usually selected by optimizing distinct objectives, such as information preservation of the sample variance or of pairwise relative distances. Widely used orthogonal projections for dimension reduction are variants of the principal component analysis (PCA) that maximize the variance of the projected data, [39]. Preservation of relative pairwise distances asks for a near-isometric embedding, and random projections guarantee this embeddings with high probability, cf. [14, 6] and see also [1, 38, 5, 12, 32, 29]. The use of random projections is especially favorable for large, high-dimensional data ([50]), since the computational complexity is just $O(dkm)$, e.g. using the construction in [1], with $d,k\in\mathbb{N}$ being the original and lower dimensions and $m\in\mathbb{N}$ the number of samples. In contrast, PCA needs $O(d^{2}m)+O(d^{3})$ operations ([26]). Moreover, tasks that do not have all data available at once, e.g. data streaming, ask for dimension reduction methods that are independent of the data.

In the present manuscript, we study orthogonal projections regarding the interplay between

• O1)

preservation of variance,

• O2)

preservation of pairwise relative distances,

aiming for a sufficient lower-dimensional data representation. We shall consider the Euclidean distance exclusively since it is most widely used in applications, especially for error estimation. On manifolds, the geodesic distance is locally equivalent to the Euclidean distance. The two objectives O1) and O2) are directly addressed by PCA (O1) and random projections (O2). We achieve the following goals: first we clarify mathematically and numerically that the two objectives are competing, i.e. PCA and random projections preserve different kinds of information. Depending on the objectives we discuss beneficial choices of orthogonal projections and numerically find a balancing projector for a given data set. Finally, we define a general framework of augmented target (AT) loss functions for deep neural networks, that integrate information about target characteristics via features and projections. We observe that our proposed methodology can increase the accuracy in two deep learning problems.

In contrast to conventional approaches we study the joint behavior of the two objectives with respect to the entire set of orthogonal projectors. By analyzing the correlation between the variance and pairwise relative distances of projected data, we observe that O1) and O2) are competing and usually cannot be reached at the same time. In numerical learning experiments we investigate heuristic choices of projections applied to input features, for subsequent classification with support vector machine and shallow neural networks.

In view of learning frameworks, we utilize features and projections on target data. The class of augmented target loss functions incorporates suitable transformations and projections that provide beneficial representations of the target space. It is applied in two supervised deep learning problems dealing with real world data.

The first experiment is a clinical image segmentation problem in optical coherence tomography (OCT) data of the human retina. Related principles of dimension reduction for other clinical classification problems in OCT have already been successfully applied in [9]. In the second experiment we aim to categorize musical instruments based on their spectrogram, see [18] for related results. Our utilized augmented target loss functions can increase the accuracy in both experiments.

The outline is as follows. In Section 2 we address the analysis of the competing objectives and Theorem 2.5 yields the asymptotic correlation between variance and pairwise relative distances of projected data. Section 3 prepares for the numerical investigations by recalling $t$-designs as considered in [10], enabling subsequent numerics. Heuristic investigations on projected input used in a straightforward classification task, are presented in Section 4. Our framework of augmented target loss functions as modified standard loss functions for deep learning, is introduced in Section 5. Finally, in Sections 6 and 7 we present classification experiments on OCT images and musical instruments using aligned augmented target loss functions.

2. Dimension reduction with orthogonal projections

To reduce the dimension of a high-dimensional data set $x=\{x_{i}\}_{i=1}^{m}\subset\mathbb{R}^{d}$, we map $x$ into a lower-dimensional affine linear subspace $\bar{x}+V$, where $\bar{x}:=\frac{1}{m}\sum_{i=1}^{m}x_{i}$ is the sample mean and $V$ is a $k$-dimensional linear subspace of $\mathbb{R}^{d}$ with $k<d$. This mapping is performed by an orthogonal projector $p\in\operatorname{\mathcal{G}}_{k,d}$, where

[TABLE]

denotes the Grassmannian, so that the lower-dimensional data representation is

[TABLE]

with $\operatorname*{range}(p)=V$. A suitable choice of $p$ within $\operatorname{\mathcal{G}}_{k,d}$ depends on further objectives, i.e. which kind of information preservation shall be favored for subsequent analysis tasks. In the following, we consider two objectives associated to popular choices of orthogonal projectors for dimension reduction, in particular, random projectors from $\mathcal{G}_{k,d}$ and PCA. We will first observe that the two objectives are competing, especially in high dimensions, and then discuss consequences.

2.1. Objective O1)

The total sample variance111We use lower case letters for samples and upper case letters for random vectors/matrices. $\operatorname{tvar}(x)$ of $x=\{x_{i}\}_{i=1}^{m}\subset\mathbb{R}^{d}$ is the sum of the corrected variances along each dimension,

[TABLE]

PCA aims to construct $p\in\operatorname{\mathcal{G}}_{k,d}$, such that the total sample variance of (2.1) is maximized among all projectors in $\operatorname{\mathcal{G}}_{k,d}$. For other equivalent optimality criteria, we refer to [51].

The total sample variance of $px=\{px_{i}\}_{i=1}^{m}\subset V$ coincides with the one of (2.1) and satisfies

[TABLE]

for all $p\in\operatorname{\mathcal{G}}_{k,d}$. Thus, PCA achieves optimal variance preservation. The total variance (2.2) can also be expressed via pairwise absolute distances

[TABLE]

which reveals that PCA maximizes the sample mean of the projected pairwise absolute distances.

2.2. Objective O2)

In contrast to pairwise absolute distances, the Johnson-Lindenstrauss Lemma targets the global property of preservation of pairwise relative distances:

Lemma 2.1** (Johnson-Lindenstrauss, cf. [14, 38]).**

For any $0<\epsilon<1$, any $k\leq d,m\in\mathbb{N}$, with

[TABLE]

and any set $\{x_{i}\}_{i=1}^{m}\subset\mathbb{R}^{d}$, there is a projector $p\in\operatorname{\mathcal{G}}_{k,d}$ such that

[TABLE]

holds for all $i<j$.

For small $\epsilon>0$, the projector $p$ in Lemma 2.1 yields that all of the $\frac{m(m-1)}{2}$ pairwise relative distances

[TABLE]

are close to $1$, i.e. the projection $p$ preserves all scaled pairwise relative distances well. A good choice of $p$ in Lemma 2.1 is based on random projectors $P\sim\lambda_{k,d}$, where $\lambda_{k,d}$ denotes the unique orthogonally invariant probability measure on $\operatorname{\mathcal{G}}_{k,d}$. The following Theorem is essentially proved by following the lines of the proof of Lemma 2.1 in [14] after replacing the constant $4$ with $(2+\tau)2$ in the respective bound on $k$.

Theorem 2.2**.**

For any $0<\epsilon<1$, any $k\leq d,m\in\mathbb{N}$ and any $0<\tau$ with

[TABLE]

and any set $\{x_{i}\}_{i=1}^{m}\subset\mathbb{R}^{d}$, the random projector $P\sim\lambda_{k,d}$ satisfies

[TABLE]

with probability at least $1-\tfrac{1}{m^{\tau}}+\tfrac{1}{m^{\tau+1}}$.

2.3. Competing objectives

A projector $p$ satisfying the near-isometry property (2.5) implies

[TABLE]

so that the total variance of the projected data $px$ may not be maximized for $k<d$. In particular, with high probability a random projector $P\sim\lambda_{k,d}$ does not suit the objective of maximizing the total variance, and we even observe $\mathbb{E}\operatorname{tvar}(Px)=\frac{k}{d}\operatorname{tvar}(x)$, see (A.2) in the appendix. PCA does not guarantee any local geometric property and distances between pairs of points can be arbitrarily distorted [1], see [41] for more robust PCA. The preservation of larger distances is favored since PCA maximizes (2.4) among all $p\in\operatorname{\mathcal{G}}_{k,d}$ and $\|p(x_{i})-p(x_{j})\|\leq\|x_{i}-x_{j}\|$ holds for all $i<j$. Close but distinct points, could even be projected onto a single point, which violates the preservation of pairwise relative distances, see Figure 2.1.

To more quantitatively understand the relation between the two competing objectives, we consider the sample mean and the uncorrected sample variance of the pairwise relative distances (2.6),

[TABLE]

Recall that good preservation of the relative pairwise distances in (2.6) asks for $\mathcal{M}(p,x)$ being close to $1$ and the variance $\mathcal{V}(p,x)$ being small. In the following, we analyze $\operatorname{tvar}(px)$, $\mathcal{M}(p,x)$, and $\mathcal{V}(p,x)$ and their expectations for random $P\in\operatorname{\mathcal{G}}_{k,d}$.

In Figure 2.2 we see a simple numerical experiment, where we first create an independent, normally distributed fixed data set $\{x_{i}\}_{i=1}^{m}$ with $x_{i}\in\mathbb{R}^{d}$ for $i=1,\dots,m$ and $m=100$, $d=50$. We then compute PCA, for $k=10,20,30,40$, as well as $n=10000$ random projections $p$ distributed according to $\lambda_{k,50}$. In Figure 2.2 (a) - (d) we can see that the more $k$ differs from $d$, the more PCA and random projections differ concerning $\operatorname{tvar}(px)$ and $\mathcal{M}(p,x)$. Those differences may lead diverse behavior in subsequent data analysis. Moreover, we compare $\mathcal{M}(p,x)$ and $\mathcal{V}(p,x)$ in Figure 2.2 (e) - (h) for the different $k$. We can see that again when $k$ is much smaller than $d$, random projections and PCA differ more concerning the variance of pairwise distances $\mathcal{V}(p,x)$. For $k=10$ the variance for PCA is higher in comparison to random projections, see Figure 2.2 (e), for $k=40$ vice versa, see Figure 2.2 (h). Note that the theoretical bounds stated in Theorem 2.1 are much higher than the dimensions $k$ used in the experiments, but the projections still preserve relative pairwise distances very well. In [7] similar observations were made on empiric experiments with image and text data.

The amount of variance kept in the principal components comparing real world and random data has been experimentally studied, e.g. in [31] and [48]. Both studies determine that the difference occurs mainly in the first principal component.

Remark 2.3**.**

In the numerical example we compare random projections and PCA directly, serving as the corresponding projections to the objectives O1) and O2). We observe that even for not so high dimensional ($d=50$) data $x$ and $k\ll d/2$, PCA severely looses information in terms of total variance, i.e. more than $50\%$ for $k=10$, and more importantly, looses much more information on pairwise relative distances than random projections. If both types of information are of interest, pairwise relative distances and high total variance, one should therefore favor random projections over PCA for $k\ll d/2$ to balance the two objectives O1) and O2) and vice versa. Note that with a large amount of data one might still want to favor random projectors since their construction is computationally much cheaper and independent from the data. On the other hand, if objective O2) is negligible, e.g. tasks with very noisy data, then PCA would be the favorable choice for all $k$.

Information of data can be quantified and expressed in different ways. One crucial part in dimension reduction is the decision of what kind of information shall be kept, which depends on several parameters including the quality of the data and the analysis task. Variants of PCA, focusing on the preservation of variance, have been widely used in real world problems with big success, especially in denoising, when the preservation of all pairwise relative distances may be counterproductive, e.g. in dMRI imaging [53] and color filter array images [57]. Drawbacks are the necessity for all data being available from the start and the high computational costs. For very high dimensional and large data sets the computation of PCA is often not feasible. Besides the huge benefit of data independence and low computational cost when using random projections, the near-isometry property often allows to establish that the solution found in the low-dimensional space is a good approximation to the solution in the original space ([1], [37]).

Algorithms in machine learning often need or benefit from sufficient estimates of pairwise distances, e.g. approximate nearest-neighbor problems, supervised classification [29] and subspace clustering [28]. In [35] algorithmic applications of near-isometry embeddings have been introduced. In [7] random projections have been successfully applied to noisy and noiseless text and image data. The experimental studies include the comparison of preservation of pairwise distances between random projections and PCA. The results coincide with our observations, that for $k>d/2$ PCA is able to preserve the pairwise distances sufficiently, whereas for $k<d/2$ PCA distorts them. The smaller $k$ the worse the distortion, whereas random projections preserve similarities still well for very small $k$, while being computationally much cheaper than PCA. One should point out again that favoring preservation of pairwise distances relies on the accuracy of the original distances.

PCA and random projections are orthogonal projections favoring two different aims. We want to study in the context of the whole set of orthogonal projections if the two objectives O1) and O2) could be reached at the same time. We will see that the objectives act competing and therefore we suggest a balancing projector for tasks that benefit from both objectives.

2.4. Covariances and correlation between competing objectives

For further mathematical analysis we first introduce a more general class of probability measures on $\operatorname{\mathcal{G}}_{k,d}$ that resemble $\lambda_{k,d}$ sufficiently well:

Definition 2.4**.**

A Borel probability measure $\lambda$ on $\operatorname{\mathcal{G}}_{k,d}$ is called a cubature measure of strength $t$ if

[TABLE]

where $\operatorname{Pol}_{t}(\mathbb{R}^{d^{2}})$ denotes the set of multivariate polynomials of total degree $t$ in $d^{2}$ variables.

Existence of cubature measures is studied, for instance, in [16]. For random $P$, we now determine the expectation values for our $3$ quantities of interest: $\operatorname{tvar}(Px)$, $\mathcal{M}(P,x)$, and $\mathcal{V}(P,x)$. If $P\sim\lambda$ and $\lambda$ is a cubature measure of strength at least $2$, the identities (A.2) and (A.3) in the appendix and a short calculation yield

[TABLE]

where $a_{k,d}=\tfrac{2d(d-k)}{k(d-1)(d+2)}$. The expected sample variance in (2.12) satisfies

[TABLE]

This asymptotic bound relates to Theorem 2.2 and alludes to a near-isometry property of the type (2.7) for $k$ sufficiently large.

The following Theorem 2.5 provides a lower bound for random $P$ on the population correlation

[TABLE]

It holds for arbitrary dimensions $d$ and subsequently specifies the asymptotic behavior for $d\rightarrow\infty$:

Theorem 2.5**.**

Let $x=\{x_{i}\}_{i=1}^{m}\subset\mathbb{R}^{d}$ be pairwise different and let $P\sim\lambda$, with $\lambda$ being a cubature measure of strength at least $2$. For $d\geq\tfrac{m(m-1)}{2}$, the correlation (2.13) is bounded from below by

[TABLE]

Let $\{x_{i}\}_{i=1}^{m}\subset\mathbb{R}^{d}$ be random points, whose entries are independent, identically distributed with finite $4$-th moments, that are uniformly bounded in $d$. Then (2.14) converges towards $1$ in probability for $d\rightarrow\infty$.

The strong correlation for large dimensions $d$ in the second part of Theorem 2.5 suggests that increasing $\operatorname{tvar}(Px)$ may also lead to increasing $\mathcal{M}(P,x)$, see Figure 2.3 for illustration. Thus, large projected total variance $\operatorname{tvar}(Px)$ and the preservation of scaled pairwise distances, i.e. $\mathcal{M}(P,x)$ being close to $1$, are competing properties. As discussed in Section 2.3, the choice of which kind of information is favorable to preserve, depends on the data and the task; e.g. denoising (O1) and nearest neighbor classification (O2). PCA and random projections are extreme in preserving either O1) or O2). We will heuristically study the behavior of orthogonal projections balancing both objectives in the next section and will state a numerical experiment where a balancing projector yields highest classification accuracy.

Remark 2.6**.**

The second part of Theorem 2.5 relates to the well-known fact that random vectors in high dimensions are almost orthogonal, [4], and standard concentration of measure arguments may lead to more quantitative statements, cf. [54].

3. Preparations for numerical experiments

For the numerical experiments we need finite sets of projectors that represent the overall space well, i.e. cover $\operatorname{\mathcal{G}}_{k,d}$ properly.

3.1. Optimal covering sequences

Let the covering radius of a set $\{p_{l}\}_{l=1}^{n}\subset\operatorname{\mathcal{G}}_{k,d}$ be denoted by

[TABLE]

where $\|\cdot\|_{\operatorname*{F}}$ is the Frobenius norm. The smaller the covering radius, the better the set $\{p_{l}\}_{l=1}^{n}$ represents the entire space $\operatorname{\mathcal{G}}_{k,d}$. I.e., there are smaller holes and the points $\{p_{l}\}_{l=1}^{n}$ are better distributed within $\operatorname{\mathcal{G}}_{k,d}$. Following Lemma 2.1 we can connect finite sets of projections and their covering radius to the near-isometry property:

Lemma 3.1**.**

Let $\{p_{l}\}_{l=1}^{n}\subset\operatorname{\mathcal{G}}_{k,d}$ and denote $\varrho:=\varrho(\{p_{l}\}_{l=1}^{n})$. For any $0<\epsilon<1$, any $m,k,d\in\mathbb{N}$ with

[TABLE]

and any $\{x_{i}\}_{i=1}^{m}\subset\mathbb{R}^{d}$, there is $l_{0}\in\{1,\ldots,n\}$ such that

[TABLE]

where $\delta=\epsilon+2\varrho\sqrt{\frac{(1+\epsilon)d}{k}}+\frac{d}{k}\varrho^{2}$.

Proof.

Given an arbitrary projector $p\in\operatorname{\mathcal{G}}_{k,d}$, there is an index $l_{0}\in\{1,\dots,n\}$ such that

[TABLE]

From here, standard computations imply Lemma 3.1. We omit the details. ∎

The accuracy of the near-isometry property in (3.2) depends on the covering radius. Therefore, a set $\{p_{l}\}_{l=1}^{n}\in\operatorname{\mathcal{G}}_{k,d}$ with a small covering radius $\varrho$ is more likely to contain a projector with better preservation of pairwise relative distances. According to [11], it holds that 111We use the symbols $\lesssim$ and $\gtrsim$ to indicate that the corresponding inequalities hold up to a positive constant factor on the respective right-hand side. The notation $\asymp$ means that both relations $\lesssim$ and $\gtrsim$ hold. $\varrho\gtrsim n^{-\frac{1}{k(d-k)}}$, and we shall see next, how to achieve this lower bound.

A set of projectors $\{p_{l}\}_{l=1}^{n}\subset\operatorname{\mathcal{G}}_{k,d}$ is called a $t$-design if the associated normalized atomic measure $\frac{1}{n}\sum_{l=1}^{n}\delta_{p_{l}}$ is a cubature measure of strength $t$ (see Definition 2.4), see [46] for general existence results. Any sequence of $t_{i}$-designs $\{p^{i}_{l}\}_{l=1}^{n_{i}}\subset\operatorname{\mathcal{G}}_{k,d}$ with $t_{i}\rightarrow\infty$ satisfies

[TABLE]

and moreover, the bound $n_{i}\gtrsim t_{i}^{k(d-k)}$ holds, cf. [16, 11]. To relate $n_{i}$ with $\varrho_{i}$ via $t_{i}$, a sequence of $t_{i}$-designs $\{p^{i}_{l}\}_{l=1}^{n_{i}}\subset\operatorname{\mathcal{G}}_{k,d}$ is called a low-cardinality design sequence if $t_{i}\rightarrow\infty$ and

[TABLE]

For their existence and numerical constructions, we refer to [22] and [10, 11]. According to [11], see also (3.3) and (3.4), any low-cardinality design sequence $\{p^{i}_{l}\}_{l=1}^{n_{i}}$ covers asymptotically optimal, i.e.,

[TABLE]

Benefiting from the covering property, we will use low-cardinality design sequences as a representation of the overall space of orthogonal projectors $\operatorname{\mathcal{G}}_{k,d}$.

3.2. Linear least squares fit

With the linear least squares fit we can directly gain information about the relation between $\mathcal{M}(p,x)$ and $\operatorname{tvar}(px)$ for a given data set $x=\{x_{i}\}_{i=1}^{m}\subset\mathbb{R}^{d}$ when $p$ varies. Given the two samples

[TABLE]

the linear least squares fitting provides the best fitting straight line,

[TABLE]

where $s$ and $\gamma$ are determined by the sample variances and the sample covariance. If $\{p_{l}\}_{l=1}^{n}$ is a $2$-design, then the sample (co)variances coincide with the respective population (co)variances for $P\sim\lambda_{k,d}$, see Appendix A.3 for further details. It follows that

[TABLE]

The quantities $s$ and $\gamma$ can be directly computed, where $\operatorname{tvar}(x)$ is given by (2.2) and the covariances are stated in Corollary A.1. Note that (3.6) and (3.7) are now independent of the particular choice of $\{p_{l}\}_{l=1}^{n}$.

The correlation between the two samples (3.5) yields additional information about their relation. As before, if $\{p_{l}\}_{l=1}^{n}$ is a $2$-design, then the sample correlation coincides with the population correlation (2.13) for $P\sim\lambda_{k,d}$, cf. Appendix A.3. High correlation for a specific data set $x$ suggests that random projections and PCA preserve competing properties, whose benefits need to be assessed for the specific subsequent task.

4. Numerical experiments in pattern recognition

We investigate the impact on classification accuracy when applying specific orthogonal projections to input data. The real world data chosen yields a straightforward classification task, serving as a toy example for comparing the accuracy of several projected input data in simple learning frameworks. Projectors are chosen from a $t$-design in view of $\operatorname{tvar}(px)$ and $\mathcal{M}(p,x)$. For all computations made in this Section the ‘Neural Network’ and ‘Statistics and Machine Learning’ toolboxes in MatlabR2017a are used.

We use the publicly available iris data set from the UCI Repository of Machine Learning Database suitable for supervised classification learning. It consists of $3$ classes with $50$ instances each, where each class refers to a type of iris plant. The instances are described by $4$ features resulting in the input samples $\{x_{i}\}_{i=1}^{150}\subset\mathbb{R}^{4}$ and target samples $\{y_{i}\}_{i=1}^{150}\subset\{0,1\}^{3}$. For comparison we classify the diverse input data with support vector machine (SVM) and $3$-layer neural networks (NN) with $5$ and $10$ hidden units (HU).

4.1. Choice of orthogonal projection

In the experiment we use projections $p\in\operatorname{\mathcal{G}}_{2,4}$ reducing the original dimension from $d=4$ to $k=2$. As a finite representation of the overall space, we use a t-design of strength $14$ from a low-cardinality sequence (see Section 3.1) consisting of $8475$ orthogonal projectors. Note that the dimension reduction in practice takes place by applying $q\in\mathcal{V}_{k,d}$ with $q^{\top}q=p\in\operatorname{\mathcal{G}}_{k,d}$, where

[TABLE]

denotes the Stiefel manifold. When taking norms, $p$ and $q$ are interchangeable, i.e., $\|q(x)\|^{2}=\|p(x)\|^{2}$, for all $x\in\mathbb{R}^{d}$. Therefore we can use w.l.o.g. the theory developed for $p$.

The projections are chosen in a deterministic manner viewing the previously described competing properties. In Figure 4.1 the three quantities $\operatorname{tvar}(px)$, $\mathcal{M}(p,x)$ and $\mathcal{V}(p,x)$ are pairwise plotted for all projectors in $\{p_{l}\}_{l=1}^{8475}$. For comparison we choose the following projections $p\in\{p_{l}\}_{l=1}^{8475}\subset\operatorname{\mathcal{G}}_{2,4}$, see Figure 4.1(a) for a visualization.

• $p_{\times}$

closest to the expected values $1$ and $\tfrac{k}{d}\operatorname{tvar}(x)$ (see (2.10) and (2.11)),

• $p_{\Diamond}$

preserving $\mathcal{M}(p,x)\approx 1$ and maximizing $\operatorname{tvar}(px)$,

• p_{\text{\scalebox{0.8}{\square}}}

preserving $\mathcal{M}(p,x)\approx 1$ and minimizing $\operatorname{tvar}(px)$,

• p_{\text{\scalebox{0.6}{\bigcirc}}}

$\operatorname{tvar}(px)\approx\operatorname{tvar}(p_{\Diamond}x)$ and maximizing $\mathcal{M}(p,x)$,

• $p_{\Huge\star}$

minimal $\operatorname{tvar}(px)$,

• $p_{*}$

maximal $\operatorname{tvar}(px)$ (PCA).

4.2. Results

In Figure 4.1(b) we see the linear least squares fitting line, computed directly and via the slope and intercept as stated in (3.6) and $\eqref{gammapop}$. The correlation coefficient (2.13) is $0.98$, which suggests that preserving the two properties is highly competing and needs to be balanced.

In Table 4.1 the classification results of the iris data are presented. We can see that in this comparison the projector $p_{\Diamond}$, which corresponds to preserving $\mathcal{M}(p,x)\approx 1$ and maximizing $\operatorname{tvar}(px)$, yields the highest and most robust results. It even yields better results than working with the original input data. The projections that preserve $\mathcal{M}(p,x)\approx 1$ but do not take care of the magnitude of the total variance yield much worse results. On the other hand, the projections that just focus on high total variance still do not yield as high results as the projection $p_{\Diamond}$ that balances both properties.

Remark 4.1**.**

Given a data set $x$, the projector $p_{\Diamond}$ is a good choice to balance both objectives O1) and O2). It can be computed by directly analyzing $\{\operatorname{tvar}(p_{1}x),\ldots,\operatorname{tvar}(p_{n}x)\}$ and $\{\mathcal{M}(p_{1},x),\ldots,\mathcal{M}(p_{n},x)\}$ of a finite covering $\{p_{l}\}_{l=1}^{n}$ of $\operatorname{\mathcal{G}}_{k,d}$. For higher dimensions an accurate representation of $\operatorname{\mathcal{G}}_{k,d}$, in order to heuristically select $p_{\Diamond}$, requires large computational costs. The least squares regression line for a $2$-design, as stated in 3.7, can be directly computed with low computational cost. This offers helpful information about the interplay between O1) and O2).

5. Augmented target loss functions

In the previous section projectors were applied to input features of shallow neural networks. In more complex architectures, such as deep neural networks, the adaption of weights can be viewed as optimization of input features, e.g. arising features can be used for transfer learning [56]. Whereas the input data is processed and optimized in each iteration, the target data stays usually unchanged during the whole learning process, serving as measure of accuracy. The representation of the target data is one key property for successful approximation with neural networks. Here, we will introduce a general class of loss functions, i.e. augmented target (AT) loss functions, that use projections and features to yield beneficial representations of the target space, emphasizing important characteristics.

In optimization problems additional penalty terms are used for regularization or to enforce other beneficial constraints. In deep learning, weight decay (i.e. Tikhonov regularization) is a standard adaption of the loss function to that effect. Incorporating additional underlying information via features of the output/target data has been studied in diverse settings tailored to particular imaging applications. Perceptual loss functions have been used in [34] for image super-resolution, incorporating the comparison of high-level image features that arise from pretrained convolutional neural networks, i.e. the VGG-network [47]. Deep perceptual similarity metrics have been proposed in [19] for generating images, comparing image features instead of the original images. In [30] a similar approach was successfully used for style transfer and super-resolution, adding a network that defines loss functions. Anatomically constrained neural networks (ACNN) have been introduced in [42] and applied to cardiac image enhancement and segmentation. Their loss functions incorporate structural information by using autoencoders to gain features about lower dimensional parametrization of the segmentation. Brain segmentation was studied in [24], where information about the desired structure has been added in the loss function via an adjacency matrix. It was used for fine-tuning the supervised learned network with unlabeled data, reducing the number of abnormalities in the segmentation.

The information of certain target characteristics can be very powerful and even replace the need of annotations in some tasks. In [49] label-free learning is approached by using just structural information of the desired output in the loss function instead of annotated target values.

In the following, we will define a general framework of loss functions that add information of target characteristics via features and projections in supervised learning tasks.

5.1. General framework

Let the training data be input vectors $\{x_{i}\}_{i=1}^{m}\subset\mathbb{R}^{r}$ with associated target values $\{y_{i}\}_{i=1}^{m}\subset\mathbb{R}^{s}$. We consider training a neural network

[TABLE]

where $\theta\in\mathbb{R}^{N}$ corresponds to the vector of all free parameters of a fixed architecture. In each optimization step for $\theta$, the network’s output $\{\hat{y}_{i}=f_{\theta}(x_{i})\}_{i=1}^{m}\subset\mathbb{R}^{s}$ is compared with the targets $\{y_{i}\}_{i=1}^{m}$ via an underlying loss function $L$.

In contrast to ordinary learning problems with highly accurate target data, complicated learning tasks arising in many real world problems do not yield sufficient results when optimizing neural networks with standard loss functions $L$, such as the widely used mean least squares error

[TABLE]

The training data may include important information that is obvious for humans, but poorly represented within the original target data and therefore lacks consideration in the learning process. To overcome this issue, we propose to add information tailored to the particular learning problem represented by additional features of the outputs and targets.

First, we select transformations

[TABLE]

to enable error estimation in transformed output/target spaces. Note that the transformations $T_{j}$ are not required to be linear. However, they should be piecewise differentiable to enable subsequent optimization of the loss function with gradient methods. We shall allow for additional weighting of the transformations $T_{1},\dots,T_{d}$ to facilitate the selection of features for a specific learning problem. The previous sections suggest that orthogonal projections can provide favorable feature combinations, which essentially turns into a weighting procedure.

To enable suitable projections, we stack the $d$ output/target features

[TABLE]

so that applying a projector $p\in\operatorname{\mathcal{G}}_{k,d}$ to each column of $T(y_{i})$ yields $p(T(y_{i}))\in\mathbb{R}^{d\times t}$. We now define the augmented target loss function with projections by

[TABLE]

where $\alpha>0$ and $L$, $\tilde{L}$ correspond to conventional loss functions. Apparently, $L_{p}$ depends on the choice of $p\in\operatorname{\mathcal{G}}_{k,d}$. The projection $p(T(y_{i}))$ weighs the previously chosen feature transformations $T(y_{i})$. Standard choices of $L$ and $\tilde{L}$ are $L_{\operatorname{MSE}}$, in which case $L_{p}$ becomes

[TABLE]

Remark 5.1**.**

For $k=d$ the projector $p$ is the identity. In this case the transformations can map into different spaces, i.e.

[TABLE]

and we can now write the standard augmented target loss function by

[TABLE]

where $T_{1}$ corresponds to the identity function, $L^{1},\ldots,L^{d}$ are common loss functions and $\alpha_{1},\ldots,\alpha_{d}>0$ are weighting parameters.

It should be mentioned that $\alpha$ resembles a regularization parameter. The actual minimization of (5.1) among $\theta$ is usually performed through Tikhonov type regularization in many standard deep neural network implementations. The formulation (5.2) adds one further variational step for beneficial output data representation.

Remark 5.2**.**

Our proposed structure with target feature maps $T_{1},\ldots,T_{d}$ as in (5.4) relates to multi-task learning, which has been successfully used in deep neural networks [13]. It handles multiple learning problems with different outputs at the same time. In contrast to multi-task learning, we aim to solve a single problem but also penalize the error in transformed spaces enhancing certain target characteristics.

For the projected feature transformations in the augmented target loss function it is not possible to identify a balancing projection $p$ heuristically (such as $p_{\Diamond}$ in Section 4), because the output $y$ changes in each iteration when the loss function is called. In the following clinical numerical experiment we overcome this issue by using random projections and PCA in each optimization step and compare it to prior deterministic choices of projections.

6. Application to clinical image data

The first experiment is a clinical problem in retinal image analysis of the human eye, where the disruptions of the so-called photoreceptor layers need to be quantified in optical coherence tomography images (OCT). The photoreceptors have been identified as the most important retinal biomarker for prediction of vision from OCT in various clinical publications, see e.g. [25]. As OCT technology advances, clinicians are not able to look at each slice of OCT themselves (in mean they get 250 slices per patient and have 3-5 minutes/patients including their clinical examination). Therefore, automated classification of e.g. photoreceptor status is necessary for clinical guidance.

6.1. Data and objective

In this application, OCT images of different retinal diseases (diabetic macular edema and retinal vein occlusion) were provided by the Vienna Reading Center recorded with the Spectralis OCT device (Heidelberg Engineering, Heidelberg, Germany). Each patient’s OCT volume consists of $49$ cross-sections/slices ($496\times 512$ pixels) recorded in an area of $6\times 6$ mm in the center of the human retina, which is the part of the retina responsible for vision. Each of the slices was manually annotated by a trained grader of the reading center. This is a challenging and time-consuming procedure that is not feasible in clinical routine but only in a research setting. The binary pixelwise annotations serve as target values, enabling a supervised learning framework.

The objective is to accurately detect the photoreceptor layers and their disruptions pixelwise in each OCT slice by training a deep convolutional neural network with a suitable loss function. The learning problem is complicated by potentially inaccurate target annotations, as studies have shown that inconsistencies between trained graders are common, cf. [52]. Moreover, the learning task is unbalanced in the sense that there are many more slices showing none or very little disruptions. We shall observe that optimization with respect to standard loss functions performs poorly in regards to detecting disruptions. The augmented target loss function proposed in the previous section can enhance the detection.

6.2. Convolutional neural network learning

We implemented our experiments using Python 3.6 with Pytorch 1.0.0. A deep convolutional neural network $f_{\theta}$ is trained by applying the U-Net architecture reported in [45] with a sigmoid activation function and Tikhonov regularization. A set of $20$ OCT volumes ($980$ slices) from different patients with corresponding annotations are used for training, where $4$ volumes were used for calibration (validation set). Another $2$ independent test volumes were identified for evaluating the results, one without any disruptions in the photoreceptor layers, whereas the other one includes a high number of disruptions.

Each OCT slice is represented by a vector $x_{i}\in\mathbb{R}^{r}$ with $r=496\cdot 512$. The collection $\{x_{i}\}_{i=1}^{m}$ corresponds to all slices from the training volumes, i.e. $m=20\cdot 49$. Further matching the notation of the previous section, we have $r=s$ and $f_{\theta}:\mathbb{R}^{r}\rightarrow\mathbb{R}^{r}$ with binary target vectors $y_{i}\in\{0,1\}^{r}$. We observe that disruptions are not identified reliably when using the least squared loss function (5.1). To overcome this issues, we use the proposed augmented target loss function with least squared losses as stated in (5.3).

To enhance disruptions within the output/target space, we heuristically choose $d=4$ local features of the original representation. They are derived from convolutions with $2$ edge filters, $T_{1}$ (Prewitt) and $T_{2}$ (Laplacian of Gaussian), and from $2$ Gaussian highpass filters, yielding $T_{3}$ and $T_{4}$. Note that these feature transformations keep the same size, i.e. $T_{j}:\mathbb{R}^{r}\rightarrow\mathbb{R}^{r}$ for $j=1,\ldots,d$. See Figure 6.2 for example images.

We can derive several augmented target loss functions $L_{p}$ by choosing different $p\in\operatorname{\mathcal{G}}_{k,d}$ for (5.2). In this experiment we use the following projections:

• •

$p=I_{4}$,

• •

$\{p_{l}\}_{l=1}^{15}$, all projections from a t-design of strength $2$ $\subset\operatorname{\mathcal{G}}_{2,4}$ (see [10]),

• •

$p_{\operatorname*{PCA}}\in\operatorname{\mathcal{G}}_{2,4}$, projection determined by PCA on the training data,

• •

$p_{\lambda_{2,4}}$, random projection chosen according to $\lambda_{2,4}$ in each mini-batch.

6.3. Results

Since the detection problem is highly unbalanced we use precision/recall curves [15] for evaluating the overall performance of each loss function model. The area under the curve (AUC) was used as a numerical indicator of the success rate, [43]. The higher the AUC the better the classification.

The results of the different loss functions on the independent test set are stated in Table 6.1. Due to the imbalance within the data, the photoreceptor region is identified well, but disruptions are not identified reliably when using the least squared loss function (5.1). For $\alpha=0.1$ all proposed augmented target loss functions $L_{p}$ clearly increase the success rate of the disruption quantification. Note that all projections are independent from the actual data set, except PCA that was computed beforehand on the training data.

The features itself (i.e. $p=I_{4}$) improve the quantification and weighting them by projections increases the results even more: using the fixed projection $p_{12}$ from the t-design sequence $\{p_{l}\}_{l=1}^{15}$ on the output/target features yields the highest accuracy for photoreceptors and disruptions. This corresponds to the results of the previous sections, stating that depending on the particular data there are projections in the overall space acting beneficially.

Since this projection generally cannot be found beforehand, using random projections in each loss function evaluation step is easier possible in practice and independent from the data. They can be computed very efficiently and randomization can generalize and robusten the information, cf. [37]. In the following we will view a second classification problem based on spectrograms, where augmented target loss functions with random projections can improve the accuracy.

7. Application to musical data

Here, the learning task is a prototypical problem in Music Information Retrieval, namely multi-class classification of musical instruments. In analogy to the MNIST problem in image recognition, this classification problem is commonly used as a basis of comparison for innovative methods, since the ground truth is unambiguous and sufficient annotated data are available. The input to the neural network are spectrograms of audio signals, which is the standard choice in audio machine learning. Spectrograms are calculated from the time signal using a short-time Fourier transform and taking the absolute value squared of the resulting spectra, thus yielding a vector for each time-step and a two-dimensional array, like an image, cf. [17].

Reproducible code and more detailed information of our computational experiments can be found in the online repository [27].

7.1. Data and objective

The publicly available GoodSounds dataset [44] contains recordings of single notes and scales played by several single instruments. To gain equally balanced input classes we restrict the classification problem to $6$ instruments: clarinet, flute, trumpet, violin, alto saxophone and cello. Note that the recordings are monophonic, so that each recording yields one spectrogram that we aim to correctly assign to one of the $6$ instruments.

After removing the silence [3, 40], segments from the raw audio files are transformed into log mel spectrograms [20], so that we obtain images of time-frequency representations with size $100\times 100$. One example spectrogram for each class of instruments is depicted in Figure 7.1.

7.2. Convolutional neural network learning

We implemented a fully convolutional neural network $f_{\theta}:\mathbb{R}^{r}\to[0,1]^{s}$, cf. [36], where $r=100\times 100$ and $s=6$, in Python 3.6 using Keras 2.2.4 framework [21] and trained it on the Nvidia GTX 1080 Ti GPU. The data is split into $140\,722$ training, $36\,000$ validation and $36\,000$ independent test samples. We heuristically choose $d=16$ output features arising directly from the particular output class. The transformations $T_{1},\dots,T_{16}$, with $T_{j}:\mathbb{R}^{6}\to\mathbb{R}$ for $j=1,\dots,16$, are then given by the inner product of the output/target and the feature vectors. Amongst others the features are chosen from the enhanced scheme of taxonomy [55] and from the table of frequencies, harmonics and under tones [59]. We use the proposed augmented target loss function $L_{p}$ (5.2), where $L_{1}$ corresponds to the categorical-cross-entropy loss [58] and $L_{2}$ to the mean squared error as in (5.3). We consider here two choices of $p$: the identity $I_{16}$ and random projectors $p\sim\lambda_{6,16}$ in $\operatorname{\mathcal{G}}_{6,16}$.

The deep learning model is sensitive to various hyper-parameters, including $\alpha$ and $p$, in addition to conventional parameters, such as the number of convolutional kernels, learning rate and the parameter $\beta$ for Tikhonov regularization. To find the best choices in a fair trial we utilize a random hyper-parameter search approach, where we train $60$ models and select the $3$ best ones for a more precise search over different $\alpha$ in the augmented target loss function and $\beta$ for Tikhonov regularization. This results in $212$ models that are evaluated on the training and validation set. Finally, we select the best model based on the accuracy of the validation set and evaluate it on the independent test set. For comparison we also evaluate this model with no Tikhonov regularization, i.e. $\beta=0$, see Table 7.1.

7.3. Results

Table 7.1 shows that no regularization and no features provide the poorest results. It seems that adding features with random projections have a regularizing effect and improve the results significantly. As expected, it is important to include Tikhonov regularization on $\theta$. Further enhancement happens by adding features via the modified augmented target loss function with or without additional weighting from projections. All results are very stable and are generalizing very well from training to the independent test set, see [27] for further details.

8. Acknowledgment

This work was partially funded by the Vienna Science and Technology Fund (WWTF) through project VRG12-009, by WWTF AugUniWien/FA746A0249, by International Mobility of Researchers (CZ.02.2.69/0.0/0.0/16 027/0008371), and by project LO1401. For the research, infrastructure of the SIX Center was used.

Appendix A Proof of Theorem 2.5

A.1. Proof of (2.14) in Theorem 2.5

For $\{y_{i}\}_{i=1}^{M}\subset\mathbb{R}^{d}$ and $p\in\mathcal{G}_{k,d}$, we define

[TABLE]

Given two sets, $\{y_{i}\}_{i=1}^{M_{1}},\{z_{j}\}_{j=1}^{M_{2}}\subset\mathbb{R}^{d}$, suppose that $P\in\mathcal{G}_{k,d}$ is a random matrix, distributed according to a cubature measure of strength at least $2$. The covariance is given by

[TABLE]

Using the identity, cf.  [2],

[TABLE]

directly yields

[TABLE]

Following [8, Theorem 2.4, Section 3.1] we use that

[TABLE]

holds, where $q=(d-1)d(d+2)$, $\alpha_{1}=(d+1)k^{2}-2k$ and $\alpha_{2}=2k(d-k)$. This leads to the explicit formula of the population covariance

[TABLE]

with $a_{k,d}=\tfrac{2d(d-k)}{k(d-1)(d+2)}$.

For $y:=\{y_{i}\}_{i=1}^{M}\subset\mathbb{R}^{d}\setminus\{0\}$ we set $\hat{y}_{i}:=\frac{y_{i}}{\|y_{i}\|}$, for $i=1,\ldots,M$. The identity (A.4) enables us to compute the population correlation

[TABLE]

by the explicit formulas

[TABLE]

Since the variance is always nonnegative and $\tfrac{a_{k,d}}{d}>0$, the denumerator of $\operatorname{Corr}(f(P,y),f(P,\hat{y}))$ in (A.5) satisfies

[TABLE]

The enumerator of $\operatorname{Corr}(f(P,y),f(P,\hat{y}))$ in (A.5) is estimated by

[TABLE]

For $d\geq M$, a short calculation yields $\operatorname{Cov}(f(P,y),f(P,\hat{y}))\geq 0$, so that we obtain

[TABLE]

Since the correlation is scaling invariant the choice $y=\{x_{i}-x_{j}:1\leq i<j\leq m\}$ with $M=\frac{m(m-1)}{2}$ implies (2.14) in Theorem 2.5. Incorporating the correct scaling yields the following corollary:

Corollary A.1**.**

*For a given data set $x=\{x_{i}\}_{i=1}^{m}$ and for random $P\in\operatorname{\mathcal{G}}_{k,d}$ the (co)variances of $\operatorname{tvar}(Px)$ (2.4) and $\mathcal{M}(P,x)$ (2.8) are given by *

[TABLE]

where $M=\frac{m(m-1)}{2}$ and $a_{k,d}=\frac{2d(d-k)}{k(d-1)(d+2)}$.

A.2. Proof of the second part of Theorem 2.5

For fixed parameters $\mu>0,\sigma^{2}>0$, that do not depend on $d$, let $Y_{1}\in\mathbb{R}^{d}$ be a random vector, whose squared entries are independent, identically distributed with mean $\mathbb{E}Y_{1,l}^{2}=\mu$ and variance $\operatorname*{Var}(Y_{1,l}^{2})=\sigma^{2}$, for $l=1,\ldots,d$. We immediately observe

[TABLE]

For any $c>0$, Chebychev’s inequality yields

[TABLE]

Suppose that $Y_{2},\ldots,Y_{M}$ are copies of $Y_{1}$, not necessarily independent. Then the union bound

[TABLE]

implies that

[TABLE]

holds with probability at least $1-\frac{M}{c^{2}}$. Provided that $\sqrt{d}\mu\neq c\sigma$ and $0<\sqrt{d}\mu-c\sigma$, we deduce

[TABLE]

We can choose $c=\tfrac{\mu}{\sigma}\sqrt[4]{d}$, since $0<c\leq\tfrac{\sqrt[4]{d}\mu}{\sigma}\leq\tfrac{\sqrt{d}\mu}{\sigma}$. That directly yields

[TABLE]

holds with probability at least $1-\frac{\mu^{2}M}{\sigma^{2}\sqrt{d}}$.

It follows directly that $\frac{\min_{i}(\|Y_{i}\|)^{2}}{\max_{i}(\|Y_{i}\|)^{2}}$ converges towards $1$ in probability for $d\rightarrow\infty$,

The choice $\{Y_{1},\ldots,Y_{M}\}=\{X_{i}-X_{j}:1\leq i<j\leq m\}$ implies the second part of Theorem 2.5.

A.3. Calculations for population covariances

We notice that $\|p(x_{i}-x_{j})\|^{2}=\operatorname*{trace}(px_{i}x_{i}^{\top}-px_{j}x_{j}^{\top})$ is a polynomial of degree $1$ in $p$. Hence, $\operatorname{tvar}(px)$ in (2.4) is also a polynomial of degree $1$ in $p$. If $\{p_{l}\}_{l=1}^{n}$ is a $1$-design, then the sample mean of $\{\operatorname{tvar}(p_{1}x),\ldots,\operatorname{tvar}(p_{n}x)\}$ satisfies

[TABLE]

which is the population mean of $\operatorname{tvar}(Px)$, with $P\sim\lambda_{k,d}$. Similarly, the term $\|p(x_{i}-x_{j})\|^{4}$ is a polynomial of degree $2$ in $p$, so that $(\mathcal{M}(p,x))^{2}$ in (2.8) is a polynomial of degree $2$ in $p$. If $\{p_{l}\}_{l=1}^{n}$ is a $2$-design, then we derive

[TABLE]

with $P\sim\lambda_{k,d}$. In other words, the sample variance of $\{\mathcal{M}(p_{1},x),\ldots,\mathcal{M}(p_{n},x)\}$ coincides with the population variance $\operatorname*{Var}(\mathcal{M}(P,x))$. Analogously, we deduce that the sample covariance of (3.5) coincides with the population covariance $\operatorname{Cov}(\mathcal{M}(P,x),\operatorname{tvar}(Px))$ with $P\sim\lambda_{k,d}$.