# Deterministic Decoupling of Global Features
and its Application to Data Analysis

Eduardo Martínez-Enríquez, María del Mar González, and Javier Portilla E.
Martínez-Enríquez and J. Portilla are with the Instituto de Óptica, CSIC,
Spain.
E-mails<EMAIL_ADDRESS><EMAIL_ADDRESS>M.
González is with Dept. of Mathematics, Universidad Autónoma de Madrid and
ICMAT, Spain.
E-mail<EMAIL_ADDRESS>This work has been funded by the Spanish
Government grants FIS2016-75891-P, PID2020-118071GB- I00 and
PID2020-113596GB-I00. Additionally, Grant RED2018-102650-T funded by MCIN/AEI/
10.13039/501100011033, and the “Severo Ochoa Programme for Centers of
Excellence in R&D” (CEX2019-000904-S).

###### Abstract

We introduce a method for deterministic decoupling of global features and show
its applicability to improve data analysis performance, as well as to open new
venues for feature transfer. We propose a new formalism that is based on
defining transformations on submanifolds, by following trajectories along the
features’ gradients. Through these transformations we define a normalization
that, we demonstrate, allows for decoupling differentiable features. By
applying this to sampling moments, we obtain a quasi-analytic solution for the
orthokurtosis, a normalized version of the kurtosis that is not just decoupled
from mean and variance, but also from skewness. We apply this method in the
original data domain and at the output of a filter bank to regression and
classification problems based on global descriptors, obtaining a consistent
and significant improvement in performance as compared to using classical
(non-decoupled) descriptors.

###### Index Terms:

feature-based data analysis, feature redundancy, feature decoupling, nested
normalization, feature transfer, local de-correlation, orthokurtosis,
regression, classification.

##  1 Introduction

Data analysis relies on the statistical distribution of the observed samples.
Usually, some features (i.e., real functions) are applied to extract relevant
information from the observed data. In the Machine Learning field there are
two basic scenarios for data analysis: (i) the classical one, where the set of
features is chosen ad-hoc; (ii) the Deep Learning scenario, which involves
automatically learning the features from the input data. In all cases, when
features are used, the observed dependencies among them come both from the
statistical behavior of the data and from the features’ joint algebraic
structure.

To illustrate this problem imagine that we analyze vectors representing 1-D
signals, extracting some marginal sample moments and the sample auto-
correlation. We will find a strong dependency between skewness and kurtosis,
and also between consecutive correlation factors, e.g., $\rho(1)$ and
$\rho(2)$, even when the input data are i.i.d. samples. The reason is that the
two mentioned pairs of features, like many others, are algebraically coupled.
As a consequence, their joint range is not just the outer product of their
marginal ranges: some combinations of (independently) valid feature values are
incompatible for the same input data. E.g., a skewness of 10 and a kurtosis of
3, or $\rho(1)=0.9$ and $\rho(2)=0$.

Feature coupling, thus, produces spurious redundancy, a sort of feature
entanglement, regardless of (and in addition to) the redundancy derived from
input data statistics. It causes difficulties for analysis, processing, and
simulation. Having a joint feature range with a very intricate topology (full
of “holes” and complex boundaries’ structure) is an obstacle to interpreting
the role of each feature separately from the others. It also complicates the
manipulation of the samples, in case we wanted to average feature values [1],
study the effect of modifying the value of a particular feature independently
of the others, or simulating data by imposing some values onto their features
[2]. An interesting example of these problems appeared in [3], which addressed
the problem of simplifying the set of features used in [1] to visually
describe texture samples.

Despite its large negative impact on data analysis and processing, much less
effort has been devoted in the literature to study and reverse deterministic
feature coupling compared to statistical data modeling. Note that, when
algebraic coupling between features exists (as in the examples mentioned
above), conventional techniques such as PCA, ICA, or even more advanced non-
linear ICA (see, e.g., [4]), do not provide the right tools for disentangling
the involved joint feature vector structures. Even in the ideal scenario of
perfectly modeling all dependencies, a purely statistical approach applied to
the observed features would mix up the two sources of redundancy, namely,
statistical and algebraic. It is advantageous to address separately these two
redundancy sources, as it is usual to apply the same type of features for
diverse statistical distributions (even in ANNs, when doing transfer learning
[5]). Therefore, many different real problems on different data distributions
using the same features will benefit from their decoupling.

Here we propose a mathematical and algorithmic framework for decoupling a set
of given features, in the sense of finding another set of similar functions
with their gradients mutually orthogonal everywhere in the domain - and, as a
consequence, with their ranges decoupled. We study the mathematical conditions
under which that is possible. We also study less favorable scenarios where
only a limited and/or approximated decoupling is possible. We demonstrate the
practical application of the proposed method to several examples of data
analysis, namely, statistical regression and textured image classification.
Some of the seminal ideas and applied results presented has been published in
three conference proceedings [6, 7, 2]. In this work we provide a solid
framework, unifying, extending, and giving the necessary mathematical rigor to
our previous results.

As concrete study cases, here we have focused on marginal moments and on the
second-order moments at the output of a set of filters. Marginal moments are
widely used in the signal processing and statistics literature, for analysis
tasks such as estimation (e.g., the method of moments), detection, regression,
classification, etc. [8, 9, 10, 11, 12, 13], and also for synthesis-by-
analysis [14, 15, 1]. Typically they are used either implicitly, as empirical
marginal histograms, or in their standardized form, and up to fourth order:
sample mean, variance, skewness, and kurtosis. Whereas, as shown here, the
first three standardized moments are already mutually decoupled, that is not
the case for the skewness and kurtosis. The problem of the skewness-kurtosis
coupling has been pointed out by several authors [16, 17, 18], but it had not
been fully solved. In this respect, one of the main contributions of this
paper is presenting a normalized version of the fourth-order sample moment,
the orthokurtosis, which is not just decoupled from the sample mean and
variance, but also from the skewness. This new statistical function is fully
consistent with the previous standardized sample moments (mean, variance and
skewness), that result from applying our decoupling technique to the first
three raw moments. In addition, the orthokurtosis calculation has a modest
computational cost. By using this new fourth-order feature, instead of the
classical kurtosis, we obtain a dramatic accuracy gain in several regression
problems (see Subsection 6.3). Furthermore, higher-than-four order moments
have been very rarely used (see exceptions in, e.g., [19, 20]) because of
their instability and mutual redundancy. By decoupling the marginal moments
(exactly or approximately) here we are able to exploit very high order
decoupled moments (up to 10th order) and demonstrate their positive impact for
texture classification (see Subsection 6.4.1).

Banks of convolutional filters, on the other hand, are a classical tool in
signal processing, with a huge field of application, including early human
vision modelling and image/audio analysis, processing and synthesis. In
addition, they have also been incorporated [21, 22] into artificial neural
networks (ANNs) for signals having spatial dimensions (image, video, 3-D,
etc.) with tremendous impact. In neural science, they have long been used to
model the responses at early stages of animal and human visual and auditory
systems [23, 24, 25, 26]. The latter image/audio representations share the
feature of being redundant, thus avoiding the artifacts that plague
critically-sampled linear transformations (e.g., orthogonal or bi-orthogonal
wavelets [27]). However, redundancy in non-orthogonal linear representations
demands paying a high price, namely, the algebraic coupling of undecimated
sub-bands (outputs of the filters). Here we address the problem of
deterministically decoupling the second-order moments at the output of a
filter bank, with direct application, besides analysis, to transfer [28] and
synthesis [2]. In Section 6 we show how the gradients of the resulting
decoupled features are virtually orthogonal for white noise samples and very
close to orthogonal for photographic textured image patches. Furthermore, we
demonstrate how approximately decoupling not just variance, but also higher-
order moments, at the output of a bank of filters, results in an important
performance boost in texture classification (Subsection 6.4.2).

This paper is organized as follows. Section 2 sets the mathematical
foundations of the method, that allow, in favorable cases, to transform a
given feature by decoupling it from a set of other given features. Section 3
proposes algorithms (based on the Nested Normalization concept, NeN) to obtain
a hierarchically ordered set of mutually decoupled features and to transfer
features from one observed sample to another. Section 4 addresses in detail
two study cases of features for being decoupled, namely marginal moments, and
the second-order moments at the output of a filter bank. Section 5 addresses
analytically the local de-correlation effect of feature decoupling, and why
this improves parameter discrimination, in regression problems. Section 6 is
devoted to showing how the proposed method actually decouples the studied
features, and its practical impact when it is applied to data analysis
(regression and classification). Section 7 concludes the paper. In addition,
some technical and/or very detailed contents have been encapsulated in
appendices, for readability and reproducibility sake.

##  2 Deterministic Decoupling of Global Features

In this section we propose a method for, given a set of features ${\cal S}$
and another unrelated feature $g$, finding a transformed feature $\hat{g}$,
identical to $g$ on a high dimensional submanifold, that is decoupled to every
feature in ${\cal S}$.

###  2.1 Preliminary Concepts

####  2.1.1 Global Shift-Invariant Features

In this paper we associate a finite discrete signal made of $N$ samples with a
vector ${\mathbf{x}}\in{\mathbb{R}}^{N}$, possibly lexico-graphically
reordered, if the signal support is a multi-dimensional array. We will term a
feature $f$ of that vector ${\mathbf{x}}$ a differentiable real function
$f:\overline{\Omega}\subset{\mathbb{R}}^{N}\rightarrow{\mathbb{R}}$ for a
domain $\overline{\Omega}$. We define here global feature a feature that
depends on all vector’s coefficients.

In this paper we will focus on _shift-invariant_ features, a special case of
global features111The only non-global shift-invariant features are the trivial
functions $f({\mathbf{x}})=c$, where $c$ is a real constant.. Within them, we
will exemplify the application of our method to features of the form:

$f({\mathbf{x}})=\frac{1}{N}\sum_{n=1}^{N}{[{\mathbf{m}}({\mathbf{x}})]_{n}},$
(1)

with ${\mathbf{m}}:\overline{\Omega}\rightarrow{\mathbb{R}}^{N}$ being a
differentiable shift-equivariant (i.e., commuting with shift operations), or
shift-invariant, mapping, assuming a shift operation with boundary conditions
(e.g., circular) has been defined. This kind of functions, being averages,
play the role of sample statistics, like marginal moments, correlation
coefficients, moments at the output of filters, etc.

####  2.1.2 Decoupled Features

We say that two features $f_{i}$ and $f_{j}$ are algebraically decoupled (from
now on just decoupled) on a subset of $\overline{\Omega}$ iff

$\nabla f_{i}({\mathbf{x}})\cdot\nabla f_{j}({\mathbf{x}})=0,\quad\text{for
all }{\mathbf{x}}\text{ in that subset}.$ (2)

We extend this concept to a set of features ${\cal S}=\\{f_{j},j=1\dots M\\}$
by terming that the features of a set are decoupled, iff they are mutually
decoupled, i.e., iff all possible pairs of features within that set
$\\{(i,j)\,:\,i,j\in\\{1..M\\},i\neq j\\}$ are decoupled. Similarly, we say
that a feature is decoupled to a (decoupled or not) set of features iff it is
decoupled to each of the features in that set.

It is worth pointing out two special cases, namely, when features are
trivially decoupled and when they are trivially coupled. We term _trivially
decoupled_ features those for which there exists at least one orthogonal basis
where they have disjoint supports222Note that they can not have disjoint
supports in the original domain if they are both global. (e.g. Fourier,
orthogonal wavelets, etc.). Here we refer to support, in a given domain, as
the subset of vector indices the feature depends on. On the other extreme, a
feature map ${\mathbf{f}}$ is _trivially coupled_ iff it exists at least one
non-degenerate function $F:{\mathbb{R}}^{M}\rightarrow{\mathbb{R}}$ such that
$F({\mathbf{f}}({\mathbf{x}}))=0,\forall{\mathbf{x}}\in\overline{\Omega}$. In
this paper, we present methods for decoupling features assuming none of those
situations happens (for which decoupling is either unnecessary or impossible,
respectively).

####  2.1.3 Normalization map

The construction of a normalization will be key in our decoupling process.

Let ${\cal S}=\\{f_{j}:\overline{\Omega}\rightarrow{\mathbb{R}},j=1\dots M\\}$
be a set of features, $\Omega$ a subset of $\overline{\Omega}$,
$\hat{\mathbf{x}}_{\mathcal{S}}({\mathbf{x}}):\Omega\rightarrow\Omega$ be a
continuous non-constant mapping, and ${\mathbf{v}}^{ref}$ a vector made of
$\\{v_{j}^{ref},j=1\dots M\\}$ (the reference values), some jointly compatible
given reference values for the features in ${\mathcal{S}}$. We say that
$\hat{\mathbf{x}}_{\mathcal{S}}({\mathbf{x}};{\mathbf{v}}^{ref})$ is a
normalization w.r.t. ${\mathcal{S}}$ and ${\mathbf{v}}^{ref}$ in $\Omega$ iff
it holds that

  * (i)

$\\{f_{j}(\hat{\mathbf{x}}_{\mathcal{S}}({\mathbf{x}};{\mathbf{v}}^{ref}))=v_{j}^{ref}\\}_{j=1}^{M}$;

  * (ii)

if $\\{f_{j}({\mathbf{x}})=v_{j}^{ref}\\}_{j=1}^{M}$ then
$\hat{\mathbf{x}}_{\mathcal{S}}({\mathbf{x}};{\mathbf{v}}^{ref})={\mathbf{x}}$.

Note that previous conditions imply that every normalization is idempotent:

$\hat{\mathbf{x}}_{\mathcal{S}}(\hat{\mathbf{x}}_{\mathcal{S}}({\mathbf{x}};{\mathbf{v}}^{ref});{\mathbf{v}}^{ref})=\hat{\mathbf{x}}_{\mathcal{S}}({\mathbf{x}};{\mathbf{v}}^{ref}).$
(3)

We now set up some notation for this set of reference values which will be
useful in our exposition. Let
${\mathbf{f}}_{\mathcal{S}}:\overline{\Omega}\rightarrow{\mathbb{R}}^{M}$ be
the vector transformation made of the ordered features in ${\mathcal{S}}$,
$[f_{1}({\mathbf{x}})\dots f_{M}({\mathbf{x}})]$, and set ${\mathbf{v}}^{ref}$
to be an $M$-dimensional vector. We define a reference manifold as
${\mathcal{R}}_{\mathcal{S}}({\mathbf{v}}^{ref})={\mathbf{f}}_{\mathcal{S}}^{-1}({\mathbf{v}}^{ref})$,
i.e., the set of vectors ${\mathbf{x}}$ such that
${\mathbf{f}}_{\mathcal{S}}({\mathbf{x}})={\mathbf{v}}^{ref}$ .

For being a valid set of reference values for the features,
${\mathbf{v}}^{ref}$ must be made of jointly compatible values of the
functions in an algebraic sense, i.e.,
$\\{{\mathbf{x}}:{\mathbf{f}}_{\mathcal{S}}({\mathbf{x}})={\mathbf{v}}^{ref}\\}\neq\emptyset$.
In addition, we will assume a _non-degeneracy_ hypothesis, under which all
feature gradients are linearly independent at every point (this will be
condition C1 in Subsection 2.3.1). This is a stronger condition than the set
of features not being trivially coupled, and it implies that the dimension of
the reference manifold is everywhere $N-M$.

###  2.2 Motivating example: decoupling two features

In order to motivate the general algorithm, let us explain the method in the
case of two features.

####  2.2.1 Gradient systems

Let us fix one feature $f$, defined in a connected open set
$\overline{\Omega}$. We study the trajectories ${\mathbf{x}}(t)$ of the
initial value problem

$\left\\{\begin{aligned} &\frac{d{\mathbf{x}}}{dt}=-\nabla f({\mathbf{x}}),\\\
&{\mathbf{x}}(0)={\mathbf{x}}_{0}.\end{aligned}\right.$ (4)

This ODE is known as a gradient system. It is clear that moving along a (non-
constant) trajectory in the $+t$ (resp. $-t$) direction will strictly decrease
(resp. increase) the value of the function $f$ until it reaches the minimum
(resp. maximum) value or stabilize at a critical point of $f$ (see the
reference [29, Section 9.3] for a discussion of this type of systems).

Thus, in order to study the set of values that $f$ may take along a
trajectory, one needs to impose some constraints on its equilibrium points.
The precise set of conditions are given in Appendix A, and can be summarized
in:

  * B1.

Maxima and minima are all global extrema, not just local.

  * B2.

The set made of the basins of attraction of all saddle points, denoted by
$\Lambda$, is of lower dimension.

The first condition ensures that the trajectory will not stop at a local non-
global extreme. The second condition guarantees that saddle points are non-
degenerate and thus, essentially unstable, so it allows to circumvent them by
adding small perturbations on the initial condition ${\mathbf{x}}_{0}$ (see
Section 3.3.2).

We denote the trajectory that passes through the point ${\mathbf{x}}_{0}$, or
equivalently, the integral manifold of the system (4), by
${\mathcal{I}}({\mathbf{x}}_{0},f)$. The main property we will need is that
all possible $f$ values are reachable from any initial point
${\mathbf{x}}_{0}$ by moving along the gradient. As a consequence, fixed a
reference value $v_{f}^{ref}$ for $f$, conditions B1-B2 guarantee that, for
${\mathbf{x}}\not\in\Lambda$, each trajectory ${\mathcal{I}}({\mathbf{x}},f)$
reaches (only once) this value. See Appendix A for a technical discussion.

####  2.2.2 Decoupling via normalization

Assume that we are given two features $\\{f,g\\}$. We would like to replace
$g$ by a “similar” feature $\hat{g}$ that is decoupled from $f$ in the sense
given by (2), via normalization.

More precisely, fixed one feature $f$, we would like to construct a
normalization map $\hat{\mathbf{x}}_{f}({\mathbf{x}})$ as defined in Section
2.1.3. A possibility is to choose a point in the trajectory
${\mathcal{I}}({\mathbf{x}},f)$ that attains some reference value $v^{ref}$ of
the feature $f$. Thus it is natural to define the normalization by
$\hat{{\mathbf{x}}}_{f}({\mathbf{x}};v^{ref})$, this is, as the map that sends
a point ${\mathbf{x}}$ to the point where the trajectory
${\mathcal{I}}({\mathbf{x}},f)$ crosses the reference manifold
${\mathcal{R}}_{f}(v^{ref})$, which is unique by our previous discussion on
gradient systems.

Now, given another feature $g$, we define $\hat{g}$ by

$\hat{g}({\mathbf{x}})=g(\hat{\mathbf{x}}_{f}({\mathbf{x}})).$ (5)

Then $f$ and $\hat{g}$ are decoupled, this is, their gradients are orthogonal.
To show this, apply the chain rule to Eq. (5), to obtain
$\nabla\hat{g}_{f}({\mathbf{x}})={\mathbf{J}}_{\hat{{\mathbf{x}}}_{f}}^{T}({\mathbf{x}})\nabla
g(\hat{\mathbf{x}}_{f}({\mathbf{x}}))$, where
${\mathbf{J}}_{\hat{\mathbf{x}}_{f}}$ represents the Jacobian matrix of the
map $\hat{\mathbf{x}}_{f}$. Now, by pre-multiplying both terms by $\nabla
f({\mathbf{x}})^{T}$, we obtain:

$\nabla
f({\mathbf{x}})\cdot\nabla\hat{g}({\mathbf{x}})=({\mathbf{J}}_{\hat{\mathbf{x}}_{f}}({\mathbf{x}})\nabla
f({\mathbf{x}}))^{T}\nabla g(\hat{\mathbf{x}}_{f}({\mathbf{x}})).$ (6)

Finally, note that

${\mathbf{J}}_{\hat{\mathbf{x}}_{f}}({\mathbf{x}})\nabla
f({\mathbf{x}})=\left.\frac{d}{dt}\right|_{t=0}\hat{{\mathbf{x}}}_{f}({\alpha}(t)),$

where $\alpha(t)$ is the integral curve of (4) starting at the point
${\mathbf{x}}$. By construction, the map $\hat{\mathbf{x}}_{f}$ is constant
along this curve and, thus, the above expression vanishes. This shows that the
two gradients in expression (6) are orthogonal, as desired. Figure 1
graphically illustrates these concepts for two given features $f_{1},f_{2}$.

Figure 1: Decoupling two features through normalization. The proposed
normalization consists in finding the intersection of the invariance
submanifolds passing by ${\mathbf{x}}$ (parallel curved thick lines) with the
reference manifold
${\mathcal{R}}_{1}=\\{{\mathbf{x}}:f_{1}({\mathbf{x}})=v_{1}^{ref}\\}$ (black
dashed thick line). All vectors belonging to the same invariance submanifold
(e.g., the mustard dashed curve), like ${\mathbf{x}}_{1}$ and
${\mathbf{x}}_{2}$, have the same normalized vector $\hat{x}_{1}$ and, thus,
the same $\hat{f}_{2}$ value. Therefore, $\nabla\hat{f}_{2}$ must be locally
orthogonal everywhere to the invariance submanifolds (iso-level sets of
$\hat{f}_{2}$), and, as a consequence, also to $\nabla f_{1}$.

###  2.3 Decoupling features from a given set: multi-feature normalization

Now we consider the problem of decoupling a feature $g$ from a given set
${\mathcal{S}}=\\{f_{i}:\overline{\Omega}\subset{\mathbb{R}}^{N}\rightarrow{\mathbb{R}},i=1\dots
M\\}$ of $M$ features. The method follows the normalization scheme explained
in the two-feature case, by following trajectories given by the gradients of
all the $f_{i}$. Unlike the simple case of gradient systems, in order to build
integral manifolds from multiple feature gradients, additional conditions must
be fulfilled. Indeed, we will give a necessary and sufficient condition for
this method to apply (Proposition 2.2 below).

It is also clear that, in order to construct a normalization, we need to
restrict to a subdomain $\Omega$, obtained from $\overline{\Omega}$ by
removing its critical points corresponding to the set of given features
${\mathcal{S}}=\\{f_{j},j=1...M\\}$.

####  2.3.1 Invariant mapping with respect to a set of features

We first introduce the notion of a mapping being invariant with respect to a
set of features
${\mathcal{S}}=\\{f_{i}:\overline{\Omega}\subset{\mathbb{R}}^{N}\rightarrow{\mathbb{R}},i=1\dots
M\\}$, a concept that will greatly facilitate the decoupling of an arbitrary
feature $g$ from this set.

We say that a non-constant, differentiable mapping
${\mathbf{y}}_{\mathcal{S}}:\Omega\rightarrow\Omega$ is _invariant_ w.r.t.
${\mathcal{S}}$ iff

${\mathbf{J}}_{{\mathbf{y}}_{\mathcal{S}}}({\mathbf{x}})\nabla
f_{i}({\mathbf{x}})={\bf 0},\forall
f_{i}\in{\mathcal{S}},\forall{\mathbf{x}}\in\Omega,$ (7)

where ${\mathbf{J}}_{{\mathbf{y}}_{\mathcal{S}}}$ represents the Jacobian
matrix of ${\mathbf{y}}_{\mathcal{S}}$.

######  Proposition 2.1.

Obtaining decoupled features from invariant mappings. Let $g$ be an arbitrary
feature $g:\overline{\Omega}\rightarrow{\mathbb{R}}$, and
${\mathbf{y}}_{\mathcal{S}}$ an invariant mapping w.r.t. a set of features
${\mathcal{S}}$. From them we construct a new feature:

$\hat{g}_{\mathcal{S}}({\mathbf{x}})=g({\mathbf{y}}_{\mathcal{S}}({\mathbf{x}})).$
(8)

Then it holds that $\nabla\hat{g}_{\mathcal{S}}({\mathbf{x}})\cdot\nabla
f_{i}({\mathbf{x}})=0,\,\forall
f_{i}\in{\mathcal{S}},\forall{\mathbf{x}}\in\Omega$, i.e., the new feature
$\hat{g}_{\mathcal{S}}$ is decoupled from all features in ${\mathcal{S}}$.

###### Proof.

Mimicking the calculation in Eq. (6), we have

$\displaystyle\nabla
f_{i}({\mathbf{x}})\cdot\nabla\hat{g}_{\mathcal{S}}({\mathbf{x}})$
$\displaystyle=$
$\displaystyle({\mathbf{J}}_{{\mathbf{y}}_{\mathcal{S}}}({\mathbf{x}})\nabla
f_{i}({\mathbf{x}}))^{T}\nabla g({\mathbf{y}}_{\mathcal{S}}({\mathbf{x}}))$
(9) $\displaystyle=$ $\displaystyle 0,\,\forall
f_{i}\in{\mathcal{S}},\forall{\mathbf{x}}\in\Omega$

where the last equality holds because of Eq. (7). ∎

Now that the significance of having an invariant mapping has been established,
let us consider the problem of existence. We will need to assume that:

  * C1.

The gradients $\\{\nabla f_{i}({\mathbf{x}})\\}$ are linearly independent at
every point ${\mathbf{x}}$; and

  * C2.

they satisfy the Frobenius condition, which is an integrability condition for
several gradients $\nabla f_{i}$. The related technicalities are addressed in
Appendix B.

We will always assume condition C1 in order to avoid redundancy in the set of
features ${\mathcal{S}}$, even if not explicitly stated.

Now we show that the Frobenius condition C2 is necessary and sufficient for an
invariant mapping to apply. This gives a criterion for the possibility of
exact decoupling.

######  Proposition 2.2.

Given ${\mathcal{S}}$ as above, there exists an invariant mapping
${\mathbf{y}}_{\mathcal{S}}$ w.r.t. ${\mathcal{S}}$ iff the gradients
$\\{\nabla f_{i},\,i=1\dots M\\}$ satisfy the Frobenius condition C2 at each
point.

This proof will be given in two steps. The only if part will be postponed to
the Appendix (Section B.1) because it is rather technical and not relevant to
our study. Here we will concentrate in the if statement, which will be treated
in Subsection 2.3.3 below. Our proof is explicit, giving a precise
construction of the invariant map via normalization, which is the crucial
ingredient.

####  2.3.2 Invariance submanifolds

Given a set of features
${\mathcal{S}}=\\{f_{i}:\overline{\Omega}\rightarrow{\mathbb{R}},i=1\dots
M\\}$, the _invariance submanifold_
${\mathcal{I}}({\mathbf{x}}_{0},{\mathcal{S}})$ is an $M$-dimensional
submanifold passing through ${\mathbf{x}}_{0}$ whose tangent planes at each
point are spanned by the gradients $\\{\nabla f_{i},\,i=1\dots M\\}$.

To ensure that the invariance submanifold exists we need to assume conditions
C1 and C2 above, as it is explained in Appendix B. Indeed, Frobenius condition
C2 is a compatibility condition on the second derivatives of different $f_{i}$
which is needed for multi-feature integrability. Moreover, it is is vacuous if
$M=1$, as we only integrate along the gradient of a single feature, see Eq.
(4).

In addition, Frobenius theorem states that $\Omega$ is foliated by invariance
submanifolds.

####  2.3.3 Normalization of several features

Here we give the construction of a normalization of multiple features
generalizing the approach in Section 2.2 for the decoupling of two features,
and show that this normalization indeed yields an invariant mapping.

For this, we associate a single vector
${\mathbf{y}}_{\mathcal{S}}({\mathbf{x}})$ to each invariance submanifold
${\mathcal{I}}({\mathbf{x}},{\mathcal{S}})$, this is,

${\mathbf{y}}_{\mathcal{S}}({\mathbf{z}})={\mathbf{y}}_{\mathcal{S}}({\mathbf{x}}_{0}),\,\,\forall{\mathbf{z}}\in{\mathcal{I}}({\mathbf{x}}_{0},{\mathcal{S}}),$
(10)

and ensure that such mapping ${\mathbf{y}}_{\mathcal{S}}({\mathbf{x}})$ is
continuous and differentiable. The remaining question is, then, how to choose
a representative ${\mathbf{y}}_{\mathcal{S}}({\mathbf{x}}_{0})$ of each
invariance submanifold ${\mathcal{I}}({\mathbf{x}}_{0},{\mathcal{S}})$. In the
case of a normalization, we choose the vector belonging to
${\mathcal{I}}({\mathbf{x}}_{0},{\mathcal{S}})$ that attains some reference
values ${\mathbf{v}}_{\mathcal{S}}^{ref}$ in its features, values that we know
${\mathcal{I}}({\mathbf{x}}_{0},{\mathcal{S}})$ may take, for all
${\mathbf{x}}_{0}\in\Omega$.

The following proposition states that, after removing a lower-dimensional
subset $\Lambda$ from $\Omega$, we can attain a valid set of reference values
by moving along the invariance submanifold:

######  Proposition 2.3.

Let ${\mathbf{v}}^{ref}$ be any jointly compatible set of values of
${\mathbf{f}}_{\mathcal{S}}$. Then for all
${\mathbf{x}}_{0}\in\Omega\setminus\Lambda,$ there exists
${\mathbf{z}}\in{\mathcal{I}}({\mathbf{x}}_{0},{\mathcal{S}})$ that satisfies
${\mathbf{f}}_{\mathcal{S}}({\mathbf{z}})={\mathbf{v}}^{ref}$, and it is
unique in the connected component of $\Omega\setminus\Lambda$ where
${\mathbf{x}}$ belongs to. Thus, the solution set
${\mathcal{I}}({\mathbf{x}}_{0},{\mathcal{S}})\bigcap{\mathcal{R}}_{\mathcal{S}}({\mathbf{v}}^{ref})$
contains exactly one point in this connected component.

###### Proof.

By our assumptions on the critical points, the basin of attraction of critical
points that are not global maxima or minima is lower dimensional. Note that
the gradient flow of each $f_{i}$ starting at ${\mathbf{x}}_{0}$ is fully
contained in $\mathcal{I}({\mathbf{x}}_{0},{\mathcal{S}})$, and thus, the
$f_{i}$’s take all possible values along the flow unless ${\mathbf{x}}_{0}$
belongs to the basin of attraction of a saddle.

Next, there cannot be more than one point in
${\mathcal{I}}({\mathbf{x}}_{0},{\mathcal{S}})$ with exactly the same
reference values since, under our assumptions on critical points, the flow of
each gradient always strictly decreases (resp. increases) the value of the
corresponding $f_{i}$. ∎

Thanks to the previous Proposition, for
${\mathbf{x}}\in\Omega\setminus\Lambda$ we can define the normalization

$\hat{{\mathbf{x}}}_{\mathcal{S}}({\mathbf{x}};{\mathbf{v}}_{\mathcal{S}}^{ref})={\mathcal{I}}({\mathbf{x}},{\mathcal{S}})\bigcap{\mathcal{R}}_{\mathcal{S}}({\mathbf{v}}^{ref}).$
(11)

######  Proposition 2.4.

The normalization $\hat{{\mathbf{x}}}_{\mathcal{S}}$ constructed in (11) has
the following properties:

  * i.

$\hat{{\mathbf{x}}}_{\mathcal{S}}$ is an invariant mapping.

  * ii.

The Jacobian ${\mathbf{J}}_{\hat{{\mathbf{x}}}_{\mathcal{S}}}$, when evaluated
on the reference manifold ${\mathcal{R}}_{\mathcal{S}}$, is an orthogonal
projection map on ${\mathcal{R}}_{\mathcal{S}}$. Moreover, it is non-
degenerate, i.e.,
$rank({\mathbf{J}}_{\hat{{\mathbf{x}}}_{\mathcal{S}}}({\mathbf{x}}))=N-M$
$\forall{\mathbf{x}}\in\Omega$.

  * iii.

The pair
$\left({\mathbf{f}}_{\mathcal{S}}({\mathbf{x}}),\hat{{\mathbf{x}}}_{\mathcal{S}}({\mathbf{x}};{\mathbf{v}}^{ref})\right)$
carries the same information as ${\mathbf{x}}$. In particular, ${\mathbf{x}}$
can always be recovered from it by reversing the normalization, i.e.,
${\mathbf{x}}=\hat{{\mathbf{x}}}_{\mathcal{S}}(\hat{{\mathbf{x}}}_{\mathcal{S}}({\mathbf{x}};{\mathbf{v}}^{ref});{\mathbf{f}}_{\mathcal{S}}({\mathbf{x}}))$,
${\mathbf{x}}\in\Omega\setminus\Lambda$.

###### Proof.

The fact that $\hat{{\mathbf{x}}}_{\mathcal{S}}$ is an invariant mapping is
obvious from the construction. Indeed, by definition of Jacobian matrix,

${\mathbf{J}}_{\hat{{\mathbf{x}}}_{\mathcal{S}}}({\mathbf{x}})\nabla
f_{i}({\mathbf{x}})=\left.\frac{d}{dt}\right|_{t=0}\hat{{\mathbf{x}}}_{\mathcal{S}}({\alpha}_{i}(t))$

for a curve $\alpha_{i}(t)$ satisfying $\alpha_{i}(0)={\mathbf{x}}$ and
$\alpha_{i}^{\prime}(0)=\nabla f_{i}({\mathbf{x}})$. Since this curve can be
taken fully contained inside the invariant submanifold
${\mathcal{I}}({\mathbf{x}},{\mathcal{S}})$ (following, for instance, the
gradient flow of $f_{i}$), then
$\hat{{\mathbf{x}}}_{\mathcal{S}}(\alpha_{i}(t))$ is a constant function in
$t$ and thus, its derivative vanishes.

For the second statement note first that, by applying the chain rule in the
condition of Eq. (3) it immediately yields that

${\mathbf{J}}_{\hat{\mathbf{x}}_{\mathcal{S}}}(\hat{\mathbf{x}}_{\mathcal{S}}({\mathbf{x}})){\mathbf{J}}_{\hat{\mathbf{x}}_{\mathcal{S}}}({\mathbf{x}})={\mathbf{J}}_{\hat{\mathbf{x}}_{\mathcal{S}}}({\mathbf{x}}),\,\,\forall{\mathbf{x}}\in\Omega.$

Now, recall that for ${\mathbf{x}}$ in the reference manifold
${\mathcal{R}}_{\mathcal{S}}$ we have
$\hat{\mathbf{x}}_{\mathcal{S}}({\mathbf{x}})={\mathbf{x}}$, so the previous
equation reduces to

$({\mathbf{J}}_{\hat{\mathbf{x}}_{\mathcal{S}}})^{2}({\mathbf{x}})={\mathbf{J}}_{\hat{\mathbf{x}}_{\mathcal{S}}}({\mathbf{x}})\quad\text{for
all }{\mathbf{x}}\in{\mathcal{R}}_{\mathcal{S}}.$

Moreover, since $\hat{\mathbf{x}}_{\mathcal{S}}$ is an invariant mapping, the
rows of ${\mathbf{J}}_{\hat{\mathbf{x}}_{\mathcal{S}}}$ are orthogonal to
gradients $\\{\nabla f_{j},j=1\dots M\\}$, and then the projection is on the
space orthogonal to the linear span of the previous gradients. That is, on the
local tangent space in ${\mathcal{R}}_{\mathcal{S}}$.

In addition, the non-degeneracy of the Jacobian follows from a classical
result in linear algebra that states that the dimension $N$ is the sum of the
dimension of the kernel ($M$ in our case) plus the rank of the matrix.

The last statement is a consequence of $\Omega\setminus\Lambda$ being foliated
by invariance submanifolds. ∎

##  3 The Nested Normalization Method

So far we have proposed a method for decoupling a new feature with respect to
a given set of features. Here we apply the results of the previous analysis in
a particular hierarchical fashion, and study how and under which conditions we
can obtain a set of mutually decoupled features.

###  3.1 Analysis

Let us consider a set ${\mathcal{S}}$ of $M$ ordered non-trivially coupled
global features $\\{f_{i}({\mathbf{x}}),i=1..M\\}$. We propose here a
sequential algorithm that, starting by taking the first original feature
$\hat{f}_{1}=f_{1}$ unchanged, aggregates at each step $k$ a new feature
$\hat{f}_{k+1}$, as shown in Algorithm 1:

0: Coupled features $\\{f_{j},\,\,j=1,\dots,M\\}$

1: Initialization: $\hat{f}_{1}=f_{1}$

2: for $k=1$ to $M-1$ do

3: $\hat{f}_{k+1}\leftarrow\text{decouple}(f_{k+1},\\{\hat{f}_{i},i=1\dots
k\\})$

4: end for

5: return Decoupled features $\\{\hat{f}_{j},j=1,\dots,M\\}$

Algorithm 1 NeN: A hierarchical decoupling approach.

Our particular strategy involves constructing suitable normalization maps
following this sequential aggregation scheme. For this, we use hierarchically
nested reference manifolds:

${\mathcal{R}}_{M-1}\subset\dots\subset{\mathcal{R}}_{1}\subset{\mathcal{R}}_{0}=\Omega\subset{\mathbb{R}}^{N},$

where, in the notation of Section 2.1.3,
${\mathcal{R}}_{k}=\hat{\mathbf{f}}_{k}^{-1}({\mathbf{v}}^{ref})$, being
$\hat{\mathbf{f}}_{k}$ a map made from the ordered set of features
$\hat{\mathcal{S}}_{k}=\\{\hat{f}_{1},\dots,\hat{f}_{k}\\}$, and a
corresponding set of reference values ${\mathbf{v}}_{k}^{ref}$. At each step
$k$, we obtain a normalization map $\hat{\mathbf{x}}_{k}$ with respect to
$\hat{\mathcal{S}}_{k}$ and this is, precisely, what defines and gives its
name to the _Nested Normalization_ (NeN) method.

The proposed nested structure has some consequences. First, although each
normalization onto ${\mathcal{R}}_{k}$ imposes a new reference value to the
$k$ feature, it respects the previously normalized values for the features
$j=1\dots k-1$. Second, it implies that

${\mathcal{R}}_{k}=\hat{\mathbf{f}}_{k}^{-1}({\mathbf{v}}_{k}^{ref})={\mathbf{f}}_{k}^{-1}({\mathbf{v}}_{k}^{ref}),$
(12)

because $\hat{\mathbf{f}}_{k}({\mathbf{x}})={\mathbf{f}}_{k}({\mathbf{x}})$
when ${\mathbf{x}}\in{\mathcal{R}}_{k-1}$. This property shows that, under
these constraints, one can define reference values with respect to
${\mathbf{f}}_{k}$ or $\hat{\mathbf{f}}_{k}$, interchangeably.

We present two variants of the NeN algorithm: _Broad_ and _Narrow_ paths,
presented in Subsections 3.1.1 and 3.1.2, respectively, depending on the
choice of the integration path.

####  3.1.1 A broad path to normalization

This scheme is precisely explained in Algorithm 2. At the $k$-th step in the
Algorithm, we start with a set of features
$\hat{\mathcal{S}}_{k}=\\{\hat{f}_{1},\ldots,\hat{f}_{k}\\}$, constructed
inductively. Assume that these satisfy the Frobenius condition C2. Then the
broad path scheme yields a normalization $\hat{\mathbf{x}}_{k}$ with respect
to the $\hat{\mathcal{S}}_{k}$ by integrating along trajectories inside the
invariance submanifold of $\hat{S}_{k}$. Note that all trajectories made of
linear combinations of the features’ gradients imposing the desired
normalization values provide the same normalization result, as they belong to
the same integral manifold (which tells us that the order of integration does
not change the output normalization). Therefore, this method provides us with
valuable degrees of freedom for choosing convenient integration paths.

More formally, we look for suitable combinations of coefficients
$\alpha_{j,k}$, such that the initial value problem:

$\frac{d{\mathbf{y}}_{k}(t)}{dt}=\sum_{j=1}^{k}{\alpha_{j,k}(t)\nabla\hat{f}_{j}({\mathbf{y}}_{k}(t))},$
(13)

with ${\mathbf{y}}_{k}(0)={\mathbf{x}}$, can be integrated in
${\mathbf{y}}_{k}(t,\vec{\alpha}_{k}(t))$.

We can write, taking into account Eq. (12) for the choice of the reference
values,

$\begin{split}t_{s}&=\arg_{t}\left\\{{\mathbf{f}}({\mathbf{y}}_{k}(t,\vec{\alpha}_{k}(t)))={\mathbf{v}}_{k}^{ref}\right\\},\\\
\hat{\mathbf{x}}_{k}({\mathbf{x}};{\mathbf{v}}_{k}^{ref})&={\mathbf{y}}_{k}(t_{s},\vec{\alpha}_{k}(t_{s})),\end{split}$
(14)

with certainty that such a solution exists and is unique in a connected
domain, as it only depends on the reference values of the adjusted features,
and not on the choice of the $\alpha$ coefficients. In any case, $\alpha_{j}$
coefficients need to respect two constraints: (i) having all the sign of
$(v_{j}^{ref}-f_{j}({\mathbf{x}}))$, in order to go coordinately in the
direction of imposing the reference values to the features; and (ii) they
should not introduce any additional stationary solutions apart from the
already discussed admissible critical points of the features. Under these
constraints, each feature can be adjusted in its full range.

0: $\\{f_{j},v_{j}^{ref},\,\,j=1,\dots,M\\}$

1: Initialization: $\hat{f}_{1}({\mathbf{x}})=f_{1}({\mathbf{x}})$

2: for $k=1$ to $M-1$ do

3:
${\mathbf{f}}_{k}({\mathbf{x}})=[f_{j}({\mathbf{x}})],\,{\mathbf{v}}_{k}^{ref}=[v_{j}^{ref}],\,j=1,\dots,k$

4: ${\mathcal{R}}_{k}={\mathbf{f}}_{k}^{-1}({\mathbf{v}}_{k}^{ref})$

5: $\hat{\mathcal{S}}_{k}\,=\,\\{\hat{f}_{j},j=1\dots k\\}\,\,(*)\,\,$

6: Compute (ODEs)
$\hat{{\mathbf{x}}}_{k}({\mathbf{x}};{\mathbf{v}}_{k}^{ref})={\mathcal{I}}({\mathbf{x}},\hat{\mathcal{S}}_{k})\cap{\mathcal{R}}_{k}\,\,(*)\,\,$

7: $\hat{f}_{k+1}({\mathbf{x}})=f_{k+1}(\hat{{\mathbf{x}}}_{k}({\mathbf{x}}))$

8: end for

9: return $\\{\hat{f}_{j},j=1,\dots,M\\}$ (*) Modifications for the broad path
relaxation in Section 3.1.3: Substitute ${\mathcal{S}}_{k}$ by
$\hat{\mathcal{S}}_{k}$ and $f_{j}$ by $\hat{f}_{j}$ in Steps 5 and 6, for the
relaxed version of the NeN broad path.

Algorithm 2 Nested Normalization, Analysis - broad path.

The proposed Algorithm 1 (and its particular realization Algorithm 2) produces
a new set of features $\hat{f}_{1},\dots,\hat{f}_{M}$ such that each
$\hat{f}_{k+1}$ is decoupled from the previous ones
$\hat{f}_{1},\ldots,\hat{f}_{k}$. Unfortunately, it has several drawbacks that
often make its implementation difficult in practice.

A first obstacle in the method is the fact that the ODEs involved in computing
Step 6 of Algorithm 2 are typically difficult to solve analytically. In fact,
lacking an analytical solution for the normalization at the $k$-th iteration
translates into not being able to obtain the expression of the decoupled
feature for the $k+1$ iteration (Step 7), and beyond.

However, a more significant drawback of this scheme is that a new decoupled
feature is defined upon the previous ones. As a consequence, the loop stops
after a single decoupled feature no longer fulfills the requirements, meaning
that the next “decoupled features” in the loop simply do not exist. In
particular, Frobenius condition C2 is rather stringent, besides being usually
hard to verify since the gradients of the new features
$\hat{f}_{1},\ldots,\hat{f}_{M}$ will tend to have very convoluted
mathematical expressions.

In next subsection we develop a second version of the algorithm (_narrow
path_) that provides: i) an alternative method for feature decoupling that
does not require analytical calculations; and, ii) a way to demonstrate that
it is possible to relax the normalization at each step $k$ (also in the broad
path algorithm), by making it w.r.t. to the original feature set
${\mathcal{S}}_{k}$, instead of w.r.t. the decouple features set,
$\hat{\mathcal{S}}_{k}$. In Subsection 3.1.3 below we will discuss when both
approaches are equivalent.

####  3.1.2 The narrow path algorithm

Here we propose Algorithm 3 to construct a normalization
$\hat{\mathbf{x}}_{k}$ with respect to the features in
${\mathcal{S}}_{k}=\\{f_{1},\ldots,f_{k}\\}$. Such normalization is obtained,
at each step $k$, by moving along the gradient of each feature, projected over
the previous reference manifold. Thus our normalization $\hat{\mathbf{x}}_{k}$
is made of a sequence of 1-D integral manifolds, whose only requirement is
being each free from critical points (conditions B1-B2, see Subsection 2.2.1).

Assuming that Frobenius condition C2 holds for the gradients of
$\\{f_{1},\dots,f_{k}\\}$, then the normalization map
$\hat{\mathbf{x}}_{k}({\mathbf{x}})$ we obtain from the Algorithm 3 is an
invariant mapping with respect to the features in ${\mathcal{S}}_{k}$. This
fact follows from Proposition 2.4, since we have just provided an admissible
integration path inside the invariance submanifold
${\mathcal{I}}({\mathbf{x}},{\mathcal{S}}_{k})$ to reach
$\hat{\mathbf{x}}_{k}({\mathbf{x}})$ as defined in (11). Another consequence
of this Proposition (statement ii.) is that, for
${\mathbf{x}}\in\mathcal{R}_{k}$,

$\nabla\hat{f}_{k+1}({\mathbf{x}})=P_{\mathcal{R}_{k}}(\nabla
f_{k+1}({\mathbf{x}})),$ (15)

where we have denoted by $P_{\mathcal{R}_{k}}$ the (orthogonal) projection map
on the reference manifold. This in particular implies that
$\nabla\hat{f}_{k+1}$ is orthogonal to the linear span of the gradients
$\nabla f_{1},\dots,\nabla f_{k}$.

If, on the contrary, Frobenius condition does not hold for the original
gradients, we can still follow the direction of those projected gradients
(termed ${\mathbf{g}}_{k}$ in Algorithm 3) until we reach the desired feature
values. This will not give a new orthogonal gradient $\nabla\hat{f}_{k+1}$.
However, the algorithm still yields interesting results since it produces
approximate decoupling and improves performance in different applications (an
example is shown in Section 6).

0: ${\mathbf{x}}\in\Omega\setminus\Lambda$,
$\\{f_{j},v_{j}^{ref},\,\,j=1,\dots,M\\}$

1: Initialization: $\hat{f}_{1}=f_{1}$; ${\cal R}_{0}=\Omega\setminus\Lambda$;
$\hat{\mathbf{x}}_{0}({\mathbf{x}})={\mathbf{x}}$

2: for $k=1$ to $k=M-1$ do

3: Compute ${\mathbf{g}}_{k}=P_{{\mathcal{R}}_{k-1}}(\nabla f_{k})$

4: ${\mathbf{y}}_{k}(0)=\hat{\mathbf{x}}_{k-1}({\mathbf{x}})$

5: Follow ${\mathbf{g}}_{k}$ until $f_{k}({\mathbf{y}}_{k}(t))=v_{k}^{ref}$

6:
$\hat{\mathbf{x}}_{k}({\mathbf{x}};{\mathbf{v}}_{k}^{ref})={\mathbf{y}}_{k}(t)$

7: $\hat{f}_{k+1}({\mathbf{x}})=f_{k+1}({\hat{\mathbf{x}}_{k}})$

8: end for

9: return $\\{\hat{f}_{j}({\mathbf{x}}),j=1,\dots,M\\}$,
$\hat{x}_{M-1}({\mathbf{x}})$

Algorithm 3 Nested Normalization - narrow path.

Figure 2 illustrates the NeN algorithm in its narrow-path version, for three
dimensional vectors, defining two nested normalization levels.

Figure 2: Illustration of the NeN algorithm, in its narrow path version. The
original vector ${\mathbf{x}}$ is first normalized w.r.t. $f_{1}$, that is, it
is modified along the gradient of $f_{1}$ until reaching ${\mathcal{R}}_{1}$,
the reference manifold with all its vectors having
$f_{1}({\mathbf{x}})=v_{REF}^{(1)}$. That vector is
$\hat{\mathbf{x}}_{1}({\mathbf{x}})$, and there we can evaluate
$\hat{f}_{2}({\mathbf{x}})=f_{2}(\hat{\mathbf{x}}_{1}({\mathbf{x}}))$. From
there we follow the projection of the gradient of $f_{2}$ onto the local
hyperplane tangent to ${\mathcal{R}}_{1}$ until reaching ${\mathcal{R}}_{2}$,
the set of vectors having
$(f_{1}({\mathbf{x}})=v_{REF}^{(1)},f_{2}({\mathbf{x}})=v_{REF}^{(2)})$. That
vector is $\hat{\mathbf{x}}_{2}({\mathbf{x}})$, the normalization of
${\mathbf{x}}$ w.r.t. both $f_{1}$ and $f_{2}$. There we can evaluate
$\hat{f}_{3}({\mathbf{x}})=f_{3}(\hat{\mathbf{x}}_{2}({\mathbf{x}}))$ (and so
on).

####  3.1.3 From a narrow path to a broad relaxation

As we have seen in the previous subsection, the narrow path yields a
normalization by moving along the gradients of the original features, whereas
in the broad path we use the gradients of the modified features. These
approaches are equivalent if Frobenius condition C2 holds on the gradients of
the modified features. In fact that, under this condition, given ${\cal
S}_{k}=\\{f_{j},j=1\dots k\\}$ and $\hat{\cal S}_{k}=\\{\hat{f}_{j},j=1\dots
k\\}$ (the latter obtained with the Narrow Path Algorithm 3), then

$\hat{\mathbf{x}}_{{\cal{S}}_{k}}({\mathbf{x}})=\hat{\mathbf{x}}_{\hat{\cal{S}}_{k}}({\mathbf{x}}).$
(16)

We see that, in this favorable setting, the obtained features in
${\hat{\mathcal{S}}}_{k}$ are mutually decoupled. Moreover, since our
inductive scheme produces a $\hat{f}_{k+1}$ that is decoupled from the
previous ones, then the features in $\hat{S}_{k+1}$ will also be mutually
decoupled.

The proof of (16) is a consequence of our construction, since in Algorithms 2
and 3 the reference manifolds are the same thanks to Eq. (12). In addition, we
recall Eq. (15) that compares the gradients of the original and the modified
features when being on the reference manifolds. Thus, the solution constructed
by the narrow path algorithm is both a valid concatenation of 1-D integration
paths for the decoupled gradients (as we are following them), and for the
original gradients, as projected gradients are linear combinations of original
gradients.

As a corollary, given that the normalization result (and, thus, also the
resulting set of decoupled features) is unique, if it exists, for a given
ordered set ${\mathcal{S}}$ and their corresponding reference values
${\mathbf{v}}^{ref}$, then in the broad path Algorithm 2 we can simply
substitute $\hat{\mathcal{S}}_{k}$ by ${\mathcal{S}}_{k}$ in its Steps 5 and
6. By doing that we make it totally equivalent to the narrow path algorithm.
We call this change a relaxation of the broad path algorithm, and refer to
this modified version as its relaxed version. Besides being much easier to
implement, the broad path in its relaxed version (like the narrow path and
unlike the broad path in its original version), still provides useful results
(but not strictly decoupled) when the Frobenius condition does not hold on the
output features, as discussed below.

More generally, our previous discussion yields decoupling in the case when a
subset of the output features obtained using Algorithm 3 (or equivalent) have
gradients fulfilling the Frobenius condition:

######  Proposition 3.1.

If the set ${\cal S}_{k}=\\{f_{j},j=1\dots k\\}$ and the subset $\hat{\cal
S}_{k}=\\{\hat{f}_{j},j=1\dots k\\}$, $k\leq M$ (the latter obtained from
$\hat{\cal S}$ with Algorithm 3 or equivalent, for a given set of reference
values ${\mathbf{v}}^{ref}$) have gradients fulfilling Frobenius in ${\cal
S}_{k}$ and $\hat{\cal S}_{k}$, respectively, then all pairs
$(\hat{f}_{i},\hat{f}_{j})$, $i=1\dots k$, $j=1\dots M$, $i\neq j$, are
mutually decoupled in the whole domain $\Omega$.

Note that Frobenius condition is vacuous if $k=1$ and consequently,
$(\hat{f}_{1},\hat{f}_{j})$, for $j=2\dots M$, are decoupled unconditionally.

If Frobenius condition holds for the gradients of the original features, but
not for those of the transformed features, their corresponding gradients will
still be orthogonal on their corresponding reference manifolds, i.e.,
$\nabla\hat{f}_{i}({\mathbf{x}})\cdot\nabla\hat{f}_{j}({\mathbf{x}})=0,i\neq
j,$ for all ${\mathbf{x}}\in{\mathcal{R}}_{k-1}$, where $k=\max\\{i,j\\}$. The
proof of this fact follows from the nested structure. Indeed, assume without
loss of generality that $i<j$. Then, in ${\mathcal{R}}_{j-1}$,
$\nabla\hat{f}_{j}$ is the orthogonal projection over
${\mathcal{R}}_{j-1}\subset{\mathcal{R}}_{i-1}$ of the original gradient,
while $\nabla\hat{f}_{i}$ is orthogonal to ${\mathcal{R}}_{i-1}$ by Eq. (12).
The interest of this comes from the fact that many times, even if fully
decoupling is not possible, one can still have mutual decoupling over high
dimensional manifolds in $\Omega$. In such situations, and as a practical
consequence, when dealing with probability distributions it is convenient to
choose reference values that are the expected values of the density function.
This favors obtaining gradients close to mutually orthogonal (see Figs. 8 and
9, panels (a) and (d), in Subsection 6.2), as samples will be close to the
reference manifolds, where exact orthogonality holds.

Comparing narrow and broad path versions of the NeN algorithm, the principal
advantage of the broad path (especially in its much more convenient relaxed
form) is that it provides closed-form solutions in some favorable cases. This
usually translates on its solutions being easier to analyze and faster to
compute. On the other hand, narrow path is simpler and more systematic at the
implementation level, since it only requires explicit functions for the
original features’ gradients, and it just relies on numerical integration
along 1-D trajectories.

####  3.1.4 Normalization of homogeneous features

A key characteristic of the NeN algorithm as we have presented it so far, is
to establish a sequential order among the features to be decoupled. While
there are cases, such that of marginal moments, for which it is natural to
establish a hierarchical order for normalizing the features, there exist other
situations where this does not apply. An example is, given a bank of scaled
and rotated filters, to obtain features by applying some functions to the
filters’ outputs. In this case there is no reason to establish a hierarchy
among the features coming from the filters at the same scale, just rotated in
different angles. For these situations, a combined graph for the extracted
features, having both sequential and parallel nodes seems much more
appropriate than a purely sequential scheme.

Fortunately, the tools we have presented so far can be readily applied for (i)
“simultaneously” (the sequential order chosen for the integration does not
affect the result) changing an input sample by forcing a set of $P_{k}$
features to have their reference values (i.e., normalizing the vector w.r.t.
that set of features), and (ii) obtaining new features that are functions of
the normalized vector (applying Equation (5), and Section 2.3.3), which we
know that will be decoupled to those $P_{k}$ parallel features. Therefore, we
can easily adapt our algorithms just by changing some of the sequentially
adjusted features $f_{k}$ by $P_{k}$-dimensional features
$\overline{f}_{p}^{(k)}$ (not to be confused with the aggregated maps
${\mathbf{f}}_{k}$ from $j=1$ to $k$), without changing the underlying logic.
It must be noted, though, that this procedure does not mutually decouple these
parallel features 333Mutual decoupling can be obtained by imposing a
sequential order to these homogeneous features. Fully parallel (symmetrical)
decoupling is a hard problem requiring different techniques from the ones
presented here..

Algorithm 4 shows our strategy to “simultaneously” impose a set of reference
values to a set of homogeneous features, i.e., features having all the same
functional expression, but different parameters $\\{\vec{\lambda}_{j},j=1\dots
P_{k}\\}$. The underlying idea is simple: to express analytically the solution
of a single sequential ODE integration, one ODE excursion for each of the
homogeneous features, and then obtain an analytical expression concatenating
all these excursions. The time values $t_{j}$ are left as variables, that are
computed numerically in order to fulfill the normalization (or de-
normalization). Note the difference with Eqs. (14), where the normalization
was obtained by solving a single non-linear equation at a time (scalar $t$,
instead of a vector ${\mathbf{t}}$, like now). However, note as well that the
solution of Algorithm 4 can also be expressed in terms of those equations, by
choosing $\\{\alpha_{j}(t)\\}$’s that activate sequentially single-gradient
combinations, at times $\\{t_{j}\\}$. Finally, it must be also noted that a
variant of the previous algorithm can be used as well for the case of having
homogeneous features not in parallel, but hierarchically ordered (e.g.,
second-order moment at the output of a bank of filters). In that case we can
apply Algorithm 4 to hierarchically normalize nested subsets of features, as a
particular procedure for computing Step 6 in Algorithm 2.

0: ${\mathbf{x}}_{0}\in\Omega\setminus\Lambda,\,$
${\cal\tilde{S}}_{k}=\\{f_{j,k}({\mathbf{x}})=\tilde{f}_{k}({\mathbf{x}};\vec{\lambda}_{j})\\}$,
$\\{v_{j,k}^{ref}\\},\,\,j=1,\dots,P_{k}$

1: Solve for generic ODE ${\mathbf{y}}_{k}(t;\vec{\lambda},{\mathbf{x}})$
(analytically)

2: Initialization: ${\mathbf{x}}_{0,1}={\mathbf{x}}$

3: for $j=1$ to $P_{k}-1$ do

4:
${\mathbf{x}}_{0,j+1}={\mathbf{y}}_{k}(t_{j};\vec{\lambda}_{j},{\mathbf{x}}_{0,j})$

5: end for

6: ${\mathbf{t}}=[t_{1},\dots,t_{P_{k}}]$, ${\cal
L}=\\{\vec{\lambda}_{1},\dots,\vec{\lambda}_{P_{k}}\\}$

7: Solve for $\tilde{\mathbf{y}}_{k}({\mathbf{t}};{\cal
L},{\mathbf{x}})={\mathbf{x}}_{0,P_{k}}$ (analytically)

8: Solve for
$\hat{\mathbf{t}}=\arg_{{\mathbf{t}}}\\{f_{j,k}(\tilde{\mathbf{y}}_{k}({\mathbf{t}};{\cal
L},{\mathbf{x}}_{0}))=v_{j,k}^{ref}\\}_{j=1}^{P_{k}}$

9: return
$\hat{\mathbf{x}}_{{\cal\tilde{S}}_{k}}({\mathbf{x}}_{0})=\tilde{\mathbf{y}}_{k}(\hat{\mathbf{t}};{\cal
L},{\mathbf{x}}_{0})$

Algorithm 4 Normalization of homogeneous features

In Subsection 4.2 we study the different alternatives to approximately
decouple the second-order moments at the output of a set of hierarchically
ordered filters, among which Algorithm 4 provides the most systematic approach
still providing (partially) analytical solutions.

###  3.2 Feature Transfer

An essential characteristic of the integration along the direction of one or
several gradients is its reversibility. First, the adjustment of a set of
feature values, under the given constraints and assumptions, is always
possible for every ${\mathbf{x}}\in\Omega\setminus\Lambda$, whenever the set
of desired values are algebraically compatible (see Proposition 2.3). It is
also true, in particular, that we can de-normalize a normalized vector
(whenever the reference manifold itself is also contained in
$\Omega\setminus\Lambda$), not just for recovering the original vector (as
pointed out in the property (iii) from Proposition 2.4), but also for imposing
whatever new feature values we may aim for. Thus, the good properties of the
NeN analysis methodology presented so far allow us to change the role of
vector transformation (by integrating the gradient flows) from instrumental to
the main goal, and, as such, changing the focus from analysis to feature
transfer or even synthesis.

However, before going into how to do that, it is important to realize that a
set of mutually decoupled features have their joint range decoupled, as shown
next. Given two features $f_{i}$ and $f_{j}$ it is immediate to assess that
the decoupling condition of Eq. (2) implies that a local change in
${\mathbf{x}}$ along the gradient of one of the features does not affect the
value of the other feature. More precisely, let us assume
$\hat{{\mathcal{S}}}$ is a set of decoupled features defined in a set
$\hat{\Omega}\setminus\hat{\Lambda}$. Since we can modify each feature within
its whole range by navigating along its one-dimensional flow, independently of
the values of the other features, we obtain as a corollary that

$\operatorname{Rg}(\hat{\mathbf{f}})=\operatorname{Rg}(\hat{f}_{1})\times\ldots\times\operatorname{Rg}(\hat{f}_{M})$
(17)

in the set $\hat{\Omega}\setminus\hat{\Lambda}$.

The problems of finding the largest admissible domain $\hat{\Omega}$ for the
decoupled set of features (which implies knowing the location and subsequently
excluding all its critical points, C1 condition), and the set $\hat{\Lambda}$
of basins of all their saddles, as well as the range of each decoupled
feature, are not trivial, and depend on the particular set of features. Thus,
we leave that analysis for Section 4 and the Appendix C (C.1.1 and C.2.1),
where the decoupling of two particular sets of features is studied in detail.

####  3.2.1 Peeling the onion and covering it back with new layers

Now, the decoupling range property (17) allows to modify a signal by enforcing
arbitrary desired values (each within its valid range) for the decoupled
features without iterative corrections, opening up an unprecedented scenario.

We show below how to achieve this by means of our NeN algorithm. All we need
is applying the following 3-step procedure: (i) obtaining a set of desired
decoupled features ${\mathbf{v}}^{des}$ we want to transfer, either by
applying any of the NeN analysis algorithms we have presented so far to some
target vector data ${\mathbf{y}}$, or simply by choosing any combination of
decoupled features’ values within their valid range; (ii) normalizing the
source vector data ${\mathbf{x}}$ (the one we want to transform) up to the
$M-1$-level; and (iii) de-normalizing the previously normalized data, in
reverse order, until achieving for each decoupled feature $\hat{f}_{k}$ the
same value $v_{k}^{des}$ previously measured/chosen in step (i). More
formally:

$\begin{split}\text{(i)}&\quad{\mathbf{v}}^{des}=\hat{\mathbf{f}}({\mathbf{y}})\\\
\text{(ii)}&\quad\hat{\mathbf{w}}_{M-1}({\mathbf{x}})=\hat{\mathbf{x}}_{\mathcal{S}}({\mathbf{x}};{\mathbf{v}}_{M-1}^{ref})\\\
\text{(iii)}&\quad{\mathbf{z}}_{f}({\mathbf{x}};{\mathbf{v}}^{des})=\check{{\mathbf{x}}}_{\mathcal{S}}(\hat{\mathbf{w}}_{M-1}({\mathbf{x}});{\mathbf{v}}^{des}),\end{split}$

where we have used the symbol ${}_{\check{}}$ to indicate de-normalization.
Thanks to the reversibility of the whole process (based on the reversibility
of each ODE integration), the resulting transformed data vector
${\mathbf{z}}_{f}$ shares the same exact decoupled features
${\mathbf{v}}^{des}$ with the target vector, and the same normalized kernel
$\hat{\mathbf{w}}_{M-1}$ with the source vector ${\mathbf{x}}$:

$\displaystyle\hat{\mathbf{w}}_{M-1}({\mathbf{z}}_{f};{\mathbf{v}}_{M-1}^{ref})$
$\displaystyle=$
$\displaystyle\hat{\mathbf{w}}_{M-1}({\mathbf{x}};{\mathbf{v}}_{M-1}^{ref})$
$\displaystyle\hat{\mathbf{f}}({\mathbf{z}}_{f})$ $\displaystyle=$
$\displaystyle\hat{\mathbf{f}}({\mathbf{y}}),$

Figures 3 and 4 illustrate this process.

|
---|---
Analysis via Normalization | Transfer via De-normalization
Figure 3: Feature transfer can be done as a sequence of (i) extracting
decoupled features in a target vector data, (ii) normalization of the source
data (both (i) and (ii) represented on the left), and (iii) de-normalization
of the data normalized in (ii), by imposing the features ${\mathbf{v}}^{des}$
measured in (i), in reverse order (on the right). See Algorithm 5. Figure 4:
Three-step feature transfer using de-normalization, within the Nested
Normalization framework. The left half of the figure represent the
normalization+analysis process (first and second steps, in parallel), and the
right half of the figure (converging arrows) is the de-normalization+transfer
process (third step).

In Algorithm 5 we describe the de-normalization steps using the (reversed)
narrow-path algorithm.

This approach to feature transfer, termed Controlled Feature Adjustment in
[2], and applied to photo-realistic style transfer in [28], opens up a new
catalogue of transfer and synthesis possibilities, as explained in those
references.

0: ${\hat{\mathbf{w}}}_{M-1}\in\Omega\setminus\Lambda$,
${\mathbf{f}}^{ref}\in\operatorname{Rg}({\mathbf{f}}),{\mathbf{f}}^{des}\in\operatorname{Rg}(\hat{\mathbf{f}})$

1: Initialization: $\hat{z}_{M}=\hat{w}_{M-1}$

2: for $k=M$ to $1$ do

3: ${\mathbf{y}}(0)=\hat{z}_{k}$

4: Compute ${\mathbf{g}}_{k}=P_{{\mathcal{R}}_{k-1}}(\nabla f_{k})$

5: Follow ${\mathbf{g}}_{k}$ until $f_{k}({\mathbf{y}}_{k}(t))=v_{k}^{des}$

6: $\hat{z}_{k-1}={\mathbf{y}}(t)$

7: end for

8: return ${\mathbf{z}}_{f}=\hat{z}_{0}$

Algorithm 5 NeN, De-normalization (narrow path).

###  3.3 Critical points and perturbations

####  3.3.1 Critical points

In this subsection we look at the conditions (assuming B1, B2, C1, C2 hold)
for the critical points of the decoupled features obtained using the NeN
algorithm. This will be needed in order to establish the range of the newly
constructed features.

The critical points of a feature $\hat{f}_{k}$ are those points
${\mathbf{x}}_{k}^{*}$ satisfying $\nabla\hat{f}_{k}({\mathbf{x}}_{k}^{*})=0$.
As, when using the Nested Normalization algorithm, we have that
$\nabla\hat{f}_{k}=P_{{\mathcal{R}}_{k}}(\nabla f_{k})$ within
${\mathcal{R}}_{k-1}$, the critical point condition corresponds to the
orthogonality of $\nabla f_{k}$ to the reference manifold
${\mathcal{R}}_{k-1}$. Because, by the definition of reference manifold, the
local tangent space at ${\mathbf{x}}\in{\mathcal{R}}_{k-1}$ is the orthogonal
complement of the linear span ${\cal L}(\\{\nabla f_{j}\\}_{j=1}^{k-1})$,
having a null projection on that local plane implies in this case that $\nabla
f_{k}\in{\cal L}(\\{\nabla f_{j}\\}_{j=1}^{k-1})$, i.e., that there exist
$\\{\lambda_{j,k}\in{\mathbb{R}}\\}_{j=1}^{k-1}$ not all zero such that:

$\nabla f_{k}({\mathbf{x}}_{k}^{*})=\sum_{j=1}^{k-1}\lambda_{j,k}\nabla
f_{j}({\mathbf{x}}_{k}^{*}).$ (18)

(Note also that a trivial solution of Eq. (18) comes from having common
critical points of the original features for orders $j<k$, for which both
sides of the equation vanish.) This is precisely the same condition as C1.

By introducing the (known) structure of the gradients, we obtain a general
expression for the critical points (which, by the structure of the NeN method,
are not isolated, but entire submanifolds). In addition, by using the previous
equation plus the constraints derived from
${\mathbf{x}}_{k}^{*}\in{\mathcal{R}}_{k}$ we obtain the expression of the
critical points of $\hat{f}_{k}$ on ${\mathcal{R}}_{k}$ (see the corresponding
calculations in the two study cases, in Section 4).

Finally, it is important to note that decoupled features at level $j$ are not
defined at critical points of decoupled features at level $i$, for $i<j$. The
reason is that iso-level sets of the feature $j$ all cross orthogonally the
iso-level sets of the feature $i$, until converging all to a critical point
(maximum or minimum, a source or a drain of the gradient field of feature
$i$), and therefore the feature $j$ is not defined there. This is a direct
consequence of the normalization $\hat{\mathbf{x}}_{i}({\mathbf{x}})$
requiring the existence of a non-null gradient $\nabla f_{i}({\mathbf{x}})$
for initiating the ODE trajectory that adjusts the $i$-th feature to its
reference value, a previous step for computing
$\hat{f}_{i+1}({\mathbf{x}})=f_{i+1}(\hat{{\mathbf{x}}}_{i}({\mathbf{x}}))$.
Figure 5 (comparing the orthokurtosis to skewness and kurtosis) may help to
understand the involved concepts.

####  3.3.2 Introducing perturbations to avoid spurious basins

We have imposed that the domain in which the features in ${\mathcal{S}}$ are
defined, $\Omega$, is free from critical points. A practical way to choose the
common domain $\Omega$ for a set of features ${\mathcal{S}}$ is to find the
largest admissible set, consisting of the intersection of the domains on which
each feature is defined, and then to remove from it all critical points.

We have also assumed B1-B2 conditions, i.e., that the basins of attraction of
critical points other than the absolute maxima and minima are a set $\Lambda$
of submanifolds with a joint dimension lower than that of $\Omega$. In that
case, it is well known that a small perturbation of a point inside these
basins will take it back into $\Omega\setminus\Lambda$ with probability
one.444Although it is very easy to fabricate ad-hoc counter-examples, these
counter-examples will never happen by chance using continuous perturbation
densities in ${\mathbb{R}}^{N}$. We should also verify that for almost all
${\mathbf{x}}\in\Omega\setminus\Lambda$, the perturbed
$({\mathbf{x}}+\epsilon)$ will remain in $\Omega\setminus\Lambda$555Idem
previous footnote.. At which stage should we add a perturbation, then? The
simplest solution, because it does not even require to check if
${\mathbf{x}}\in\Lambda$, is to always add the perturbation just before
performing gradient integration.

In addition, in the context of a real application, there are two more concerns
on the practical impact of the above constraints and assumptions: (i) usually,
real-world signals are quantized, and, thus, they are not in
${\mathbb{R}}^{N}$, but in a finite, numerable subset. In particular, within a
quantized representation, it is no longer true that the probability of falling
on a lower dimension basin of attraction of a spurious critical point is zero.
Actually, quite the opposite, signal quantization will reduce the entropy and
favor the symmetry. (ii) Whereas being directly on a spurious basin of
attraction (of a saddle) will make our method to eventually get stuck, that is
not the only problem. Close-distance neighbor points, although outside those
basins of attraction, may be affected negatively in terms of the speed of
convergence of the differential equations, and should be avoided. Same can be
said about the adjustment of features to values too close to their absolute
extrema, especially when using numerical integration because of lacking
analytical solutions.

####  3.3.3 Choosing a suitable perturbation

Critical points, as we will see in the study cases (Section 4), are typically
points presenting low entropy configurations (high symmetry, a few repeated
values/patterns, etc.). The role of the added perturbation is to increase the
entropy of the signal in a suitable way, without affecting its relevant (e.g.,
perceptual) features. Here is a list of desirable characteristics of a
perturbation $\epsilon$ on a digital signal:

  * •

it should not cause a direct loss of information, i.e.,
${\mathbf{q}}({\mathbf{x}}+\epsilon)={\mathbf{q}}({\mathbf{x}})={\mathbf{x}}$,
where here ${\mathbf{q}}(\cdot)$ represents the quantization already present
in ${\mathbf{x}}$,

  * •

it should be reproducible,

  * •

it should increase the entropy (break the symmetry),

  * •

it should not affect the operational conditions (e.g., not noticeable under
human observation).

Among previous characteristics, the third and fourth ones discourage us from
naively using noise-like perturbations (i.i.d. pseudo-random coefficients),
because they may produce (i) some samples having very close values (just by
chance); and (ii) noticeable (e.g., visual, auditory, etc.) artifacts, because
of introducing unnecessarily large differences among neighbors. In appendix
subsections C.1.2 and C.2.1 we propose ad-hoc perturbations for two different
decoupling problems. In both studied cases the perturbation aims at increasing
the entropy of the signal, by increasing the number either of the distinct
signal values or of the active frequencies in the Fourier domain. As such,
these perturbations expand the decoupled feature ranges to their theoretically
broadest possible intervals.

##  4 Two study cases

###  4.1 Marginal Moments

We study here the problem of hierarchically decoupling the first $M$ sample
moments:

${\mathcal{S}}=\left\\{f_{j}({\mathbf{x}})=\frac{1}{N}\sum_{n=1}^{N}{x_{n}^{j}}\right\\}_{j=1}^{M}.$

First of all, we check that in this case the Frobenius condition C2 holds (see
Appendix B for a proof). Therefore, the corresponding set of gradients defines
a proper invariance submanifold and we can use the broad path from Subsection
3.1.1. This fact will be crucially exploited for computing in a quasi-explicit
way the fourth-order decoupled moment, that we termed the orthokurtosis. In
Section 6 we will show how the pair skewness-kurtosis is clearly inferior for
several analysis tasks than the couple skewness-orthokurtosis.

To fix notations, we will use for the normalization the moments of a zero-mean
uni-variate Gaussian distribution:

$v_{j}^{ref}=\begin{cases}(2j)!/(2^{j}j!)&\text{ for}\,j\,\text{even},\\\
0,&\text{ for}\,j\,\text{odd.}\end{cases}$

####  4.1.1 Analytic solutions: a path to the orthokurtosis

First, the gradient of the sample mean (the sample mean is both $f_{1}$ and
$\hat{f}_{1}$) is:

$\nabla\hat{f}_{1}({\mathbf{x}})=\frac{1}{N}{\mathbf{1}},$

and the normalization comes from solving the ODE

$\frac{d{\mathbf{y}}_{1}(t)}{dt}=\frac{1}{N}{\mathbf{1}}$

starting from ${\mathbf{y}}_{1}(0)={\mathbf{x}}$ and finding
$t_{1}=\arg_{t}\\{f_{1}({\mathbf{y}}_{1}(t))=0\\}$ (recall that
$\mathcal{R}_{1}={\mathbf{f}}_{1}^{-1}({\mathbf{v}}_{1}^{ref})=\\{{\mathbf{z}}\in\Omega\,:\,f_{1}({\mathbf{z}})=0\\}$).
The normalization will then be
$\hat{\mathbf{x}}_{1}({\mathbf{x}})={\mathbf{y}}_{1}(t_{1})$. In this case the
solution is straightforward:

$\displaystyle{\mathbf{y}}_{1}(t)$ $\displaystyle=$
$\displaystyle{\mathbf{x}}+\frac{1}{N}t$ $\displaystyle
f_{1}({\mathbf{y}}_{1}(t))$ $\displaystyle=$ $\displaystyle
f_{1}({\mathbf{x}})+\frac{1}{N}t$ $\displaystyle t_{1}$ $\displaystyle=$
$\displaystyle-Nf_{1}({\mathbf{x}})$
$\displaystyle\hat{\mathbf{x}}_{1}({\mathbf{x}})$ $\displaystyle=$
$\displaystyle{\mathbf{x}}-f_{1}({\mathbf{x}}){\mathbf{1}},$

which corresponds to the original vector with the sample mean subtracted to
every sample. Now we obtain the next decoupled feature,
$\hat{f}_{2}({\mathbf{x}})=f_{2}(\hat{\mathbf{x}}_{1}({\mathbf{x}}))$:

$\hat{f}_{2}({\mathbf{x}})=f_{2}({\mathbf{x}}-f_{1}({\mathbf{x}}){\mathbf{1}}),$

which is a (biased) version of the classical sample variance. Now we compute
its gradient,

$\displaystyle\nabla\hat{f}_{2}({\mathbf{x}})$ $\displaystyle=$
$\displaystyle\frac{2}{N}\hat{\mathbf{x}}_{1}({\mathbf{x}})$ (19)
$\displaystyle=$
$\displaystyle\frac{2}{N}({\mathbf{x}}-f_{1}({\mathbf{x}}){\mathbf{1}}).$

Because of the irrelevance of the factor $2/N$ for the subsequent
calculations, we drop it. Now we can modify $\hat{f}_{2}$ by moving along its
gradient without leaving ${\mathcal{R}}_{1}$, until reaching
${\mathcal{R}}_{2}={\mathbf{f}}_{2}^{-1}({\mathbf{v}}_{2}^{ref})=\\{{\mathbf{z}}\in\Omega:f_{1}({\mathbf{z}})=0,\,f_{2}({\mathbf{z}})=1\\}$.
Now ${\mathbf{y}}_{2}(0)=\hat{\mathbf{x}}_{1}({\mathbf{x}})$,

$\frac{d{\mathbf{y}}_{2}(t)}{dt}=\left({\mathbf{y}}_{2}(t)-f_{1}({\mathbf{y}}_{2}(t)){\mathbf{1}}\right).$

This ODE simplifies by noting that, when the gradient has zero sample mean,
the sample mean can not change when integrating it (the resulting curve
belongs to ${\mathcal{R}}_{1}$). As the initial value has zero mean,
$f_{1}({\mathbf{y}}_{2}(t))=0\,\,\forall t$, the ODE simplifies to:

$\frac{d{\mathbf{y}}_{2}(t)}{dt}={\mathbf{y}}_{2}(t),$

whose solution $\log({\mathbf{y}}_{2}(t))+C=t{\mathbf{1}}$ results in

${\mathbf{y}}_{2}(t)=\exp(t)\hat{\mathbf{x}}_{1}({\mathbf{x}}),$

by enforcing ${\mathbf{y}}_{2}(0)=\hat{\mathbf{x}}_{1}({\mathbf{x}})$. We see
that
$f_{2}({\mathbf{y}}_{2}(t))=f_{2}(\hat{\mathbf{x}}_{1}({\mathbf{x}}))\exp(2t)$
and the $t$ value intersecting with ${\mathcal{R}}_{2}$ is
$t_{2}=-1/2\log(f_{2}({\mathbf{x}}_{1}({\mathbf{x}})))$. Then

$\hat{\mathbf{x}}_{2}({\mathbf{x}})={\mathbf{y}}_{2}(t_{2})=\hat{\mathbf{x}}_{1}({\mathbf{x}})/\sqrt{\hat{f}_{2}({\mathbf{x}})}.$

We see that this second normalization is the standardization of
${\mathbf{x}}$. Now we can compute the next decoupled moment
$\hat{f}_{3}({\mathbf{x}})=f_{3}(\hat{\mathbf{x}}_{2}({\mathbf{x}}))$:

$\hat{f}_{3}({\mathbf{x}})=f_{3}\left(\frac{{\mathbf{x}}-f_{1}({\mathbf{x}}){\mathbf{1}}}{\sqrt{f_{2}({\mathbf{x}}-f_{1}({\mathbf{x}}){\mathbf{1}})}}\right),$
(20)

which is the sample skewness. Here again, we compute
$\nabla\hat{f}_{3}({\mathbf{x}})$ on ${\mathcal{R}}_{2}$ by differentiating in
Eq. (20). We obtain:

$\nabla\hat{f}_{3}({\mathbf{x}})=\frac{3}{N}\left(\hat{\mathbf{x}}_{2}({\mathbf{x}})^{\odot
2}-\hat{f}_{3}({\mathbf{x}})\hat{\mathbf{x}}_{2}({\mathbf{x}})-{\mathbf{1}}\right),$
(21)

where we use the symbol “$\odot$” for representing a pointwise scalar
operation for the vector coefficients (here a power).

Unlike in previous cases, now the resulting ODE equation
$\frac{d{\mathbf{y}}_{3}(t)}{dt}=\nabla\hat{f}_{3}({\mathbf{y}}_{3}(t))$ has
no known closed-form solution. However, we note two facts. First one, as
mentioned above, the original gradients fulfill the C1-C2 conditions and,
thus, define a proper invariance submanifold. Second, in this case we observe
that the linear span of $\\{\nabla\hat{f}_{j}({\mathbf{x}})\\}_{j=1}^{3}$ is
the same as the linear span of $\\{\nabla f_{j}({\mathbf{x}})\\}_{j=1}^{3}$,
$\forall{\mathbf{x}}\in\Omega$. As a consequence, both sets of features
produce the same invariance submanifolds. Thus to compute the normalization,
it is equivalent to use both broad path (Algorithm 2) or its relaxed version.
In particular, using the gradients of the original features instead of the
decoupled ones allows to find a closed-form solution for the normalization.
More precisely, we proceed as follows: first find a solution by following
$\nabla f_{3}$ until achieving zero skew, and then standardizing that result,
i.e., imposing also zero-mean and standard deviation one. The latter
adjustments are a shift and re-scaling that correspond to moving along $\nabla
f_{1}$ and $\nabla f_{2}$, admissible operations within the submanifold, that
do not affect the zero-skew condition. Thus, we can pose the much easier
Ricatti ODE equation obtained moving along the $\nabla f_{3}$,
$\frac{d{\mathbf{z}}_{3}(t)}{dt}=\nabla f_{3}({\mathbf{z}}_{3}(t))$666This is
not the only possibility for obtaining analytic solutions for the
normalization, although it seems the best, in this case. Note that other
integrable gradients, such as ${\mathbf{z}}_{3}(t)^{2}-1$ (giving rise to an
hiperbolic tangent solution) have a reduced range and converge to spurious
stationary solutions.:

$\frac{d{\mathbf{z}}_{3}(t)}{dt}=\frac{3}{N}({\mathbf{z}}_{3}(t))^{2},$

with ${\mathbf{z}}_{3}(0)=\hat{{\mathbf{x}}}_{1}({\mathbf{x}})$, whose
solution is (because of its irrelevance for the next calculations, we ignore
the $\frac{3}{N}$ factor):

${\mathbf{z}}_{3}(t)=\hat{{\mathbf{x}}}_{1}({\mathbf{x}})\odot/({\mathbf{1}}-t\hat{{\mathbf{x}}}_{1}({\mathbf{x}})).$
(22)

Then we find a numerical solution for

$t_{0}({\mathbf{x}})=\arg_{t}\\{\hat{f}_{3}(\hat{{\mathbf{x}}}_{1}({\mathbf{x}})\odot/({\mathbf{1}}-t\hat{{\mathbf{x}}}_{1}({\mathbf{x}})))=0\\},$
(23)

and we obtain the third-order normalization by standardizing
${\mathbf{z}}_{3}(t_{0}({\mathbf{x}}))$777It is easy to check that Eq. (22)
always has solution within the open interval
$(1/\min(\hat{{\mathbf{x}}}_{1}({\mathbf{x}}))),1/\max(\hat{{\mathbf{x}}}_{1}({\mathbf{x}})))$
(note that $sign(\min(\hat{{\mathbf{x}}}_{1}({\mathbf{x}})))\neq
sign(\max(\hat{{\mathbf{x}}}_{1}({\mathbf{x}})))$ for all
${\mathbf{x}}\neq{\mathbf{0}}$).:

$\hat{{\mathbf{x}}}_{3}({\mathbf{x}})=\hat{{\mathbf{x}}}_{2}\left(\hat{{\mathbf{x}}}_{1}({\mathbf{x}})\odot/({\mathbf{1}}-t_{0}({\mathbf{x}})\hat{{\mathbf{x}}}_{1}({\mathbf{x}}))\right).$

Finally, this allows us to define the fourth-order decoupled moment as
$\hat{f}_{4}({\mathbf{x}})=f_{4}(\hat{{\mathbf{x}}}_{3}({\mathbf{x}})).$ We
have termed this new decoupled sample moment $\hat{f}_{4}$ the orthokurtosis,
a function that, unlike the classical standardized fourth-order sample moment
(the kurtosis) is not just decoupled from the mean and variance, but also from
the skewness.

Note that the computation of the orthokurtosis includes a non-explicit
function, namely $t_{0}({\mathbf{x}})$. Although we could apply a similar
strategy for obtaining closed-form solutions (up to the integration parameter
value) for decoupling higher-order moments by integrating separately along
integer power gradients, that would not provide us with efficient solutions,
because we would still need to numerically find the parameter $t_{0}$ for
which the reference value for the decoupled moment is reached. For instance,
for computing the fifth order decoupled moment we need to normalize the
orthokurtosis. This implies evaluating every time in a loop this function,
which, in turn, requires the normalization in loop of the skewness. In
summary, once there are no fully analytic expressions, computationally
expensive nested search loops appear, and a piece-wise 1-dimensional ODE
integration strategy (as the one described in Algorithm 3) is preferable.

Finally, it is worth emphasizing how the first three decoupled moments
obtained using the NeN algorithm are precisely the classical standardized
moments: mean, variance and skewness. This clearly reflects how our method has
captured the pre-existing natural intuitions about decoupling features through
normalization. However, the next standardized moment, the kurtosis, breaks the
pattern of being decoupled from all previous standardized moments, as it is
algebraically coupled to skewness. This algebraic coupling has been previously
noted, and some solutions have been proposed (see, e.g., [16]). Previous
efforts have not aimed at orthogonalizing the involved gradients, and the few
proposed modifications of the kurtosis lack a theoretical foundation and have
proven inferior in their practical application to the solutions presented here
(see Fig.12 in Section 6).

Figure 5 illustrates the orthogonalization of the kurtosis with respect to the
skewness, giving rise to the orthokurtosis. It shows the actual iso-level
curves, for the case of $N=4$ (for visualization purposes, a dimension has
been removed, namely forcing that the solutions belong to the hyperplane
$\mu({\mathbf{x}})=0$). Each trajectory shown in the orthokurtosis
representation has been actually computed by integrating the projected
gradients, starting from a randomly perturbed maximum of the skewness (a
perturbation is necessary, because the new function is not defined at the
skewness’ critical points) and finishing in one minimum.

---
Figure 5: Visual represention of sample skewness (up), kurtosis (middle) and
orthokurtosis, the new fourth-order normalized moment (bottom). Dark curves
(at all three panels), and coloured curves (at the bottom), are actual iso-
level curves. Yellow/bright represents high values, and blue/dark low values.
At the bottom panel, the iso-orthokurtosis colored curves are drawn over the
iso-skewness curves (black) to show mutual orthogonality. Note how the
orthokurtosis is not defined at the extrema of the skewness.

####  4.1.2 Beyond orthokurtosis: higher-order approximately decoupled moments

As explained in Section 3, universal and exact feature decoupling is only
possible when the gradients obtained by any of the proposed algorithms fulfill
the Frobenius condition C2. The gradient of the orthokurtosis, together with
the gradients of its preceding decoupled moments (sample mean, variance and
skewness), no longer fulfill C2 condition. Therefore, an exact unconstrained
hierarchical decoupling solution is not possible beyond four-order moments.
Nevertheless, as shown in Section 3.1.3, gradients orthogonality can still
hold exactly within the reference manifolds and, as shown in Section 6.2,
approximately outside of them.

Thus, in a looser sense of “decoupling”, the lack of analytic solutions
fulfilling Frobenius beyond the fourth order is not an insurmountable obstacle
for computing higher-order approximately decoupled moments. In fact, we have
seen (Proposition 2.4 and posterior discussion, in Subsection 2.3.3) how easy
is to compute the gradient of a decoupled feature
$\nabla\hat{f}_{k}({\mathbf{x}})$ if we constrain ${\mathbf{x}}$ to belong to
the reference manifold ${\mathcal{R}}_{k-1}$. In that case it holds
$\nabla\hat{f}_{k}({\mathbf{x}})={\mathbf{g}}_{k}({\mathbf{x}})=P_{{\mathcal{R}}_{k-1}}(\nabla
f_{k}({\mathbf{x}}))$, where $P_{{\mathcal{R}}_{k-1}}$ is the projection
operator on the local tangent plane to ${\mathcal{R}}_{k-1}$ in
${\mathbf{x}}$, ${\mathbf{x}}\in{\mathcal{R}}_{k-1}$. As the orthogonal space
to that tangent plane is the linear span of the previous gradients $\nabla
f_{j},j=1\dots k-1$, the projection can be computed by finding a linear
combination of all gradients (including $\nabla f_{k}$) that is orthogonal to
the previous gradients. Because in our case the original gradients are made of
“monomial vector” of increasing orders, it is easy to solve the triangular
system of equations resulting from equating their inner products to zero [6],
yielding the projected gradients (a set of orthogonal vectors):

$\displaystyle{\mathbf{g}}_{1}({\mathbf{x}})$ $\displaystyle=$
$\displaystyle\frac{1}{N}({{\mathbf{1}}})$
$\displaystyle{\mathbf{g}}_{2}({\mathbf{x}})$ $\displaystyle=$
$\displaystyle\frac{2}{N}({\mathbf{x}}-f_{1}{{\mathbf{1}}})$
$\displaystyle{\mathbf{g}}_{3}({\mathbf{x}})$ $\displaystyle=$
$\displaystyle\frac{3}{N}\left({\mathbf{x}}^{\odot
2}-f_{2}{{\mathbf{1}}}-a_{2,3}{\mathbf{g}}_{2}({\mathbf{x}})\right)$
$\displaystyle\vdots$ $\displaystyle{\mathbf{g}}_{k}({\mathbf{x}})$
$\displaystyle=$ $\displaystyle\frac{k}{N}\left({\mathbf{x}}^{\odot
k-1}-\sum_{j=1}^{k-1}{a_{j,k}{\mathbf{g}}_{j}({\mathbf{x}})}\right),$ (24)

with

$\displaystyle c_{i,j}^{(\ell)}$ $\displaystyle=$
$\displaystyle\begin{cases}f_{i+j}-f_{i}f_{j}&\text{ if }\ell=1,\\\
c_{i,j}^{(\ell-1)}c_{\ell-1,\ell-1}^{(\ell-1)}-c_{\ell-1,i}^{(\ell-1)}c_{\ell-1,j}^{(\ell-1)}&\text{
if }\ell>1,\end{cases}$ $\displaystyle a_{j,k}$ $\displaystyle=$
$\displaystyle\begin{cases}f_{k-1}&\text{ if }j=1,\\\
c_{j-1,k-1}^{(j-1)}/c_{j-1,j-1}^{(j-1)}&\text{ if }j>1.\end{cases}$

We remind the reader that, in general, ${\mathbf{g}}_{k}=\nabla\hat{f}_{k}$
only for points belonging to their corresponding reference manifolds (the
$k$-th gradient is computed and applied in the ${\mathcal{R}}_{k-1}$
manifold). Some of them may get simpler expressions when imposing the
corresponding reference values (particularly, lower than order $k$ odd moments
vanish if we take as reference the moments of even symmetric pdf’s, e.g.,
Gaussian). They have the cross-invariance property, i.e., by integrating a
curve along the $k$-th decoupled gradient we do not change the previous
featuress $f_{j},\hat{f}_{j}\,j=1\dots k-1$. Although in this case we obtained
a closed-form (recursive) solution for the gradient projection, in case of
lacking close-form expressions we can always apply a purely numerical
orthogonalization method to the gradient vectors of the original features
(like Gram-Schmidt).

In our practical examples of applying decoupled moments to signal analysis in
Section 6, we demonstrate the usefulness of higher than four order decoupled
moments. This is, we believe, a relevant result, as the original moments of so
high order are very rarely used in the literature due to their instability and
high redundancy.

###  4.2 second-order moments at the output of a set of filters

We study here the problem of hierarchically decoupling second-order moments
measured at the output of a set of $M$ linearly independent band-pass (zero
DC-response) filters:

${\mathcal{S}}=\left\\{f_{j}({\mathbf{x}})=\frac{1}{N}\sum_{n=1}^{N}{[{\mathbf{x}}*{\mathbf{h}}_{j}]_{n}^{2}}\right\\}_{j=1}^{M}.$

Such a set of features provides an economical description of the auto-
correlation of ${\mathbf{x}}$. We will use for the normalization
$\\{v_{j}^{ref}=f_{j}({\mathbf{w}})\\}_{j=1}^{M}$, being
${\mathbf{w}}=N{\mathbf{\delta}}$, a scaled Kronecker delta. This is
equivalent to taking for reference values the expected value of
$f_{j}({\mathbf{y}})$ for ${\mathbf{y}}$ a vector made of i.i.d. zero-mean and
unit variance coefficients. This choice for reference values, being the values
of the functions applied to a given vector ${\mathbf{w}}$, guarantees the
algebraic compatibility of ${\mathbf{v}}^{ref}$, regardless of the chosen set
of filters $\\{{\mathbf{h}}_{j}\\}$. In addition, it makes the normalization
to whiten the input.

Same as in previous Subsection 4.1, we first check that in this case the
Frobenius condition C2 holds for the original gradients (see Appendix B for a
proof), and, therefore, the corresponding set of gradients defines a proper
invariance submanifold. This allows us to apply Algorithms 2 (in its relaxed
form), 3 or 4 for trying to decouple these features, and see a posteriori if
Frobenius condition also holds on the gradients of the obtained features.

####  4.2.1 General analytical approach

We first obtain the feature gradients and study their integrability. We have:

$\nabla
f_{j}({\mathbf{x}})=\frac{2}{N}{\mathbf{x}}*{\mathbf{h}}_{j}*\tilde{\mathbf{h}}_{j},$

where
$\tilde{{\mathbf{h}}}_{j}({\mathbf{n}})={\mathbf{h}}_{j}(-{\mathbf{n}})$. To
simplify this expression and its subsequent integration, we express it in the
Fourier domain, by doing the DFT of ${\mathbf{x}}$ and ${\mathbf{h}}$:

$G_{j}(X({\mathbf{\xi}}))={\cal F}\\{\nabla
f_{j}({\mathbf{x}})\\}=\frac{2}{N}|H_{j}({\mathbf{\xi}})|^{2}X({\mathbf{\xi}}),$
(25)

where ${\mathbf{\xi}}$ represents the (possibly vectorial, for $n$-D signals,
$n>1$) discrete frequencies, and, as usual, upper case letters represent the
Fourier transforms of their original lower-case counterparts (except for $G$,
that corresponds to the Fourier transform of the gradient of the features).

In order to normalize the first $k$ features, following the broad path
algorithm, we can write the relaxed version of Eq. (13) in the Fourier domain
as:

$\frac{dY_{k}({\mathbf{\xi}},t)}{dt}=\frac{2}{N}\left(\sum_{j=1}^{k}{\alpha_{j,k}|H_{j}({\mathbf{\xi}})|^{2}}\right)Y_{k}({\mathbf{\xi}},t),$

for some convenient choice of the integration path, encoded by the
coefficients $\vec{\alpha}_{k}$. Setting
$Y_{k}({\mathbf{\xi}},0)=X({\mathbf{\xi}})$ and calling
$L({\mathbf{\xi}};\vec{\alpha}_{k})$ the filter of the sum in brackets in
previous equation, the integration of the initial value problem for computing
the signal normalization is straightforward:

$Y_{k}({\mathbf{\xi}},t;\vec{\alpha}_{k})=X({\mathbf{\xi}})\exp\Big{(}\frac{2}{N}L({\mathbf{\xi}};\vec{\alpha}_{k})t\Big{)}.$
(26)

To obtain $\hat{\mathbf{x}}_{k}({\mathbf{x}},{\mathbf{v}}^{ref})$, we first
normalize $\alpha_{1,k}=1$ and then solve for

$(\vec{\alpha}_{k}^{ref},t_{k}^{ref})=\displaystyle\underset{\vec{\alpha},t}{\arg}\Big{\\{}\sum_{{\mathbf{\xi}}}{|H_{j}({\mathbf{\xi}})|^{2}\,|Y_{k}({\mathbf{\xi}},t;\vec{\alpha}_{k})|^{2}}=v_{j}^{ref}\Big{\\}}_{j=1}^{k}$
(27)

Finally,
$\hat{X}_{k}({\mathbf{\xi}})=Y_{k}({\mathbf{\xi}},t_{k}^{ref};\vec{\alpha}_{k}^{ref})\\}$,
and $\hat{\mathbf{x}}_{k}({\mathbf{x}},{\mathbf{v}}^{ref})={\cal
F}^{-1}\\{\hat{X}_{k}({\mathbf{\xi}})\\}$. Solving the non-linear system of
$k$ equations and $k$ unknowns of Eq. (27) will require a numerical
computation in a general case. However, it may also have simplified solutions
in some special cases, as we will see in next subsection.

An equivalent possibility for computing
$\hat{\mathbf{x}}_{k}({\mathbf{x}},{\mathbf{v}}^{ref})$ consists in following
the Algorithm 4: calculating analytically the result of a generic ODE
excursion using a single gradient $\nabla f_{j}$ from a given point, with
$\alpha_{j,k}=1$, as a function of $t_{j,k}$ (initially left as unknown), Step
1. Then we concatenate these 1-D trajectories (Steps 2-5) into a final
analytical solution depending on ${\mathbf{t}}_{k}=[t_{j,k}],j=1\dots k$ (Step
7). It is easy to see that such analytical solution has the same form of Eq.
(26) (but depending on a set of consecutive times $t_{j,k}$, instead of
$\alpha_{j,k}$), and that the Step 8 of that algorithm, which enforces the
solution achieving all $k$ reference values in
${\mathbf{f}}_{k}({\mathbf{x}})$, is equivalent to Eq. (27).

####  4.2.2 A case including two complementary filters

Let us consider a Parseval frame representation [30] using two filters,
$H_{1}({\mathbf{\xi}}),H_{2}({\mathbf{\xi}})$ fulfilling
$|H_{1}({\mathbf{\xi}})|^{2}+|H_{2}({\mathbf{\xi}})|^{2}=1$, e.g., a low and
high-pass kernels of a redundant wavelet transform. Here the features are
simply $f_{j}({\mathbf{x}})=\sigma_{j}^{2}({\mathbf{x}})$, $j=1,2$. Since the
Euclidean metric is preserved, we have
$\sigma_{1}^{2}+\sigma_{2}^{2}=\sigma^{2}$, the total signal variance (for
simplicity sake we assume here zero mean).

For normalizing the first feature, we integrate its corresponding gradient,
obtaining:

$Y_{1}({\mathbf{\xi}},t)=X({\mathbf{\xi}})\exp\Big{(}\frac{2}{N}|H_{1}({\mathbf{\xi}})|^{2}t\Big{)}.$

As usual, we then find the $t$ parameter providing us the desired reference
value:

$t_{1}^{ref}({\mathbf{x}})=\underset{t}{\arg}\Big{\\{}\sum_{{\mathbf{\xi}}}|H_{1}({\mathbf{\xi}})|^{2}\,|Y_{1}({\mathbf{\xi}},t)|^{2}=v_{1}^{ref}\Big{\\}}$

and we obtain $\hat{\mathbf{x}}_{1}({\mathbf{x}})={\cal
F}^{-1}\\{Y_{1}({\mathbf{\xi}},t_{1}^{ref}({\mathbf{x}}))\\}$, from which we
obtain the decoupled feature
$\hat{f}_{2}({\mathbf{x}})=f_{2}(\hat{\mathbf{x}}_{1}({\mathbf{x}}))$.

If we wanted to normalize also the second feature (e.g., in order to add more
decoupled features), we could move along the gradient of the second feature,
but, at the same time, control that the first feature does not change its
value (this is equivalent to project $\nabla f_{2}$ onto $\mathcal{R}_{1}$).
We can achieve that by dividing by the square root of the first feature
evaluated for each $t$. Such an adjustment is valid in this case because it
corresponds to moving along the sum of the two gradients (which in this case
corresponds to simply applying a scale factor to the vector):

$Y_{2}({\mathbf{\xi}},t)=\frac{\hat{X}_{1}({\mathbf{\xi}})\exp(\frac{2}{N}|H_{2}({\mathbf{\xi}})|^{2}t)}{\sqrt{(v_{1}^{ref})^{-1}\sum_{{\mathbf{\xi}}}|H_{1}({\mathbf{\xi}})|^{2}\,|\hat{X}_{1}({\mathbf{\xi}})\exp(\frac{2}{N}|H_{2}({\mathbf{\xi}})|^{2}t)|^{2}}},$
(28)

thus enforcing that
$\sum_{{\mathbf{\xi}}}|H_{1}({\mathbf{\xi}})|^{2}|Y_{2}({\mathbf{\xi}},t)|^{2}=v_{1}^{ref}$.
Same as before, we solve for the $t$ value that achieves the desired
normalization:

$t_{2}^{ref}({\mathbf{x}})=\underset{t}{\arg}\Big{\\{}\sum_{{\mathbf{\xi}}}|H_{2}({\mathbf{\xi}})|^{2}|Y_{2}({\mathbf{\xi}},t)|^{2}=v_{2}^{ref}\Big{\\}}$

and we obtain $\hat{\mathbf{x}}_{2}({\mathbf{x}})={\cal
F}^{-1}\\{Y_{2}({\mathbf{\xi}},t_{2}^{ref}({\mathbf{x}}))\\}$, from which we
could obtain another decoupled feature from any arbitrary (non-trivially
redundant) feature $g({\mathbf{x}})$; indeed,
$\hat{g}({\mathbf{x}})=g(\hat{\mathbf{x}}_{2}({\mathbf{x}}))$.

Note also that Eq. (28), although looking quite different from Eq. (26), is a
particular case of the latter with $\alpha_{1,2}=t_{1}^{ref}-\nu$,
$\alpha_{2,2}=t_{2}^{ref}-\nu$, $\nu$ being the logarithm of the square root
in the denominator of Eq. (28). No matter the adjustment method applied here
may seem arbitrary, the fulfillment of Frobenius conditions on the gradients
of the original features guarantees, jointly with the additional constraints
explained in Section 2, that the result of such normalization exists and is
unique, as it only depends on the reference values of the adjusted features,
and not on the choice for the coefficients $\alpha$’s in the linear
combinations of the feature gradients, in the ODEs.

Finally, although not mathematically proven here888Explicit expressions of the
second decoupled feature’s gradient can be obtained through implicit
derivation., it turns out that the resulting gradients of the new features
obtained in this case do not fulfill the Frobenius condition. This implies
that strict gradients’ orthogonality only holds for all pairs within their
reference manifolds, and in the whole domain for the pairs
$(\hat{f}_{1},\hat{f}_{j}),j=1\dots M$, as explained in Subsection 3.1.3.
Nevertheless, we have obtained gradients that are very close to being mutually
orthogonal also in the other cases when applying a set of bar and edge
detectors both to white noise and to textured patches of photographic images
(see Fig. 9 in Section 6.2).

###  4.3 A summary of feature-decoupling scenarios

Table I summarizes the different situations one may encounter when trying to
apply the decoupling framework proposed here to decouple a given set of
ordered features. From less favorable to more favorable, the first scenario is
when we do not have an explicit expression of the original features gradients.
The most extreme case would be that each of our features is a black box
function. In that case, a very computationally costly numerical procedure
(such as described in Section 6.2, see Eq. (33)) is the only option for
approximately computing the gradients. A better situation is when using
artificial neural networks (ANNs) with automatic differentiation for computing
the gradients of the cost functions. In that case we can evaluate the
gradients with a reasonable cost (and apply numerical integration, using the
narrow-path algorithm), but we may not be able to assess the fulfillment of
the Frobenius condition for the original features beyond the first layer. As
such, we should not expect a strict and universal decoupling beyond the second
feature (the second layer, if we decouple in a layer-wise fashion). For those
features, it becomes an empirical matter to test how far new gradients
typically are from being mutually orthogonal.

Second scenario corresponds to knowing the explicit expressions of our
features and their gradients, and knowing that they do not jointly fulfill the
Frobenius condition. In that case we can still apply the narrow-path version
of NeN and, again, test empirically how far the resulting gradients of the
(approximately) decoupled features are from being orthogonal. This
corresponds, for instance, to the case of higher-than-two order moments at the
output of a set of filters.

Third scenario is when we know the analytic expressions of the gradients of
the original features and they fulfill Frobenius. In that case the definitions
of Section 2 for the decoupled features apply for finding decoupled features
to the original ones, and we may end up finding explicit expressions for their
gradients. However, in this scenario these gradients turn out not fulfilling
Frobenius condition. First, we recall that the first decoupled feature is just
the first feature $\hat{f}_{1}({\mathbf{x}})=f_{1}({\mathbf{x}})$, as always,
and, for that reason ${\cal S}_{1}=\hat{\cal S}_{1}$, so the broad-path
algorithm applied to obtaining the second decoupled feature is the same as its
relaxed version. Furthermore, the second feature is exactly and universally
decoupled, as Frobenius condition becomes vacuous in 1-D. However, if we are
able to compute both gradients and they do not fulfill the Frobenius
condition, we then know that no added third or following features will be
universally and exactly decoupled from the two previous ones. Still, we may
find exact constrained decoupling solutions: the output features will be
exactly decoupled within the corresponding reference manifolds, but only
approximate outside them. Again, finding how close to be orthogonal are the
gradients outside those manifolds requires an empirical measurement. An
example of this situation is when decoupling the second-order moments at the
output of a set of filters overlapping in the Fourier domain.

Finally, the most favorable scenario corresponds to being able to obtain
explicit expressions for the input (coupled) and output features (by using the
broad-path NeN algorithm, Algorithm 2, either in its original form or in its
relaxed version) and their corresponding gradients, and that the two sets of
gradients fulfill Frobenius condition. In this case we obtain the full
decoupling solution, which we know is unique, universal and exact. In
addition, in this case there is a joint equivalence relationship between the
set of original features and the obtained decoupled set. E.g., when decoupling
the first three marginal moments, the decoupled result (sample mean, variance
and skewness) jointly carries exactly the same information as the original set
(first, second and third-order moments). When adding the fourth-order moment
we obtain another exact and universally decoupled feature (the orthokurtosis).
However, because of the gradient of the orthokurtosis no longer fulfills
Frobenius condition jointly with the lower order gradients, then: (1) the
joint equivalence relationship between coupled and decoupled moments up to
order four does not hold anymore, and (2) we can not further obtain exact
universal decoupling for any added feature to this set (and, in particular,
for any higher-order moments).

TABLE I:  Four feature decoupling scenarios.

| # S C E N A R I O
---|---
| 1 | 2 | 3 | 4
Original Feat. Gradients: | Analytic Expression | Non-Explicit | Yes | Yes | Yes
Frobenius | ? | No | Yes | Yes
Output Feat. Gradients: | Analytic Expression | ? | No | Yes (2) | Yes
Frobenius | ? | No | No | Yes
Joint Equivalence Coupled/Decoupled Set | No | No | No | Yes
Does accept one additional decoupled feature? | No | No | No | Yes
NeN Applicability: | Computation | Heavy/? | Medium | Medium | Light
Algorithm | 3 | 3 | 3,4 | 2,3,4
Decoupling | 1 layer | Approx/Constr | Approx/Constr | Exact&Univ.
Example (See Table II for acronyms’ description) | ANN | MF ($p>2$) | MF ($p=2$) | MM ($p\leq 3$)
See Secs.& Refs. | Future Work |  6.4, 6.4.2 &[6, 7] |  4.2, 6.2.2 &[2] |  4.1, 6.4.1, 6.2.1, 6.3

##  5 Deterministic Decoupling and Local De-Correlation

Here we show how features’ decoupling removes local covariance in the feature
space, and how this improves discrimination.

###  5.1 Covariance-free “balls” in the feature space

Let ${\mathbf{x}}\in\mathbb{R}^{N}$ be a random vector made of $N$ i.i.d.
samples obeying a probability distribution $p(x)$. Let us assume a feature set
${\cal S}$ of $M$ global features $\left\\{f_{j}\right\\}$. Define a vector
${\mathbf{c}}\in{\mathbb{R}}^{M}$ containing the expected value of the
features $f_{j}({\mathbf{x}})$ for different realizations of ${\mathbf{x}}$,
i.e., $c_{j}=\mathbb{E}\\{f_{j}({\mathbf{x}})\\},j=1\dots M$. Define a
manifold ${\cal A}({\mathbf{c}};{\cal S})={\mathbf{f}}^{-1}({\mathbf{c}})$,
i.e., the manifold containing the set of all vectors ${\mathbf{x}}$ having the
same ${\mathbf{c}}$ (note that it is non-degenerate by our initial
assumption). Describe vector samples as
${\mathbf{x}}_{i}={\mathbf{x}}_{0i}+{\mathbf{d}}_{i}$, where
${\mathbf{x}}_{0i}\in{\cal A}$, and ${\mathbf{d}}_{i}$ is a (relatively small)
sampling fluctuation, with $\mathbb{E}\\{d_{i,k}d_{i,l}\\}=0$,
$k,l\in[1,\dots,N]$, for all $k$-th and $l$-th components of the vector
${\mathbf{d}}_{i}$.

######  Proposition 5.1.

_(Decoupled features are locally uncorrelated)._ Under previous assumptions,
for $N$ large, decoupled features will have uncorrelated deviations from their
expected values, i.e.,
$\mathbb{E}\\{(\hat{f}_{n}({\mathbf{x}})-c_{n})(\hat{f}_{m}({\mathbf{x}})-c_{m})\\}\approx
0$, $n,m\in[1,\dots,M]$.

###### Proof.

For $N$ large, features’ values, as they behave like sample statistics (see
Eq. (1)), will not deviate much from their expected values and thus vector
samples $\\{{\mathbf{x}}_{i}\\}$ will be located in the vicinity of $\cal A$.
Thus, a first-order local approximation can be applied, which gives

$f_{j}({\mathbf{x}}_{i})\approx f_{j}({\mathbf{x}}_{0i})+\nabla
f_{j}({\mathbf{x}}_{0i})\cdot{\mathbf{d}}_{i}=c_{j}+\nabla
f_{j}({\mathbf{x}}_{0i})\cdot{\mathbf{d}}_{i}.$ (29)

Therefore,
$\mathbb{E}\\{(f_{n}({\mathbf{x}})-c_{n})(f_{m}({\mathbf{x}})-c_{m})\\}\approx\mathbb{E}\\{(\nabla
f_{n}({\mathbf{x}}_{0})\cdot{\mathbf{d}}({\mathbf{x}}_{0}))(\nabla
f_{m}({\mathbf{x}}_{0})\cdot{\mathbf{d}}({\mathbf{x}}_{0}))\\}$, yielding the
covariance:

$Cov(f_{n},f_{m})({\mathbf{x}})\approx\sigma_{d}^{2}\,\nabla
f_{n}({\mathbf{x}}_{0})\cdot\nabla f_{m}({\mathbf{x}}_{0}),$ (30)

where $\sigma_{d}^{2}$ is the expected quadratic dispersion of the features
fluctuations. ∎

In the decoupled features case gradients are mutually orthogonal, and thus
vector differences for the different features will be uncorrelated. In
contrast, when using coupled features, ${\mathbf{d}}_{i}$ is projected onto
non-orthogonal directions, leading to correlated sampling fluctuations in the
feature space, as illustrated in Figure 6.

Figure 6: A representation of the expected feature vector manifold in the
input space domain, showing some input samples, and the influence of feature
gradients’ correlation on the correlation in the feature domain.

###  5.2 Features’ covariance and discriminability

Let us assume now that our pdf depends on a parameter $\theta$, $p(x;\theta)$.
Consider also a global feature transformation
${\mathbf{f}}({\mathbf{x}}):\mathbb{R}^{N}\rightarrow\mathbb{R}^{M}$ meant to
be applied to vectors ${\mathbf{x}}(\theta)\in{\mathbb{R}}^{N}$ made of
samples from $p(x;\theta)$. How well can we discriminate samples coming from
similar values of $\theta$, based on ${\mathbf{f}}({\mathbf{x}}(\theta))$? For
studying this problem it is convenient to represent the samples $x(\theta)$
using an intermediate stochastic sample $x_{0}\sim p(x;\theta_{0})$ that does
not depend on $\theta$; then we obtain the final sample by applying a
deterministic invertible mapping $s_{\theta}:\mathbb{R}\rightarrow\mathbb{R}$
of $x_{0}$ depending on $\theta$: $x(\theta)=s_{\theta}(x_{0})$, such that
$x(\theta)\sim p(x;\theta)$ (re-parametrization trick [31, 32].

This allows us to study the dependency of the expected feature vector
$\bar{\mathbf{f}}$ on $\theta$, by expressing:

$\displaystyle\frac{d\bar{\mathbf{f}}(\theta)}{d\theta}$ $\displaystyle=$
$\displaystyle\frac{d\mathbb{E}\\{{\mathbf{f}}({\mathbf{x}}(\theta))\\}}{d\theta}$
(31) $\displaystyle=$ $\displaystyle\mathbb{E}\left\\{{\bf
J}_{{\mathbf{f}}}\frac{d{\mathbf{x}}(\theta)}{d\theta}\right\\}$
$\displaystyle=$ $\displaystyle\mathbb{E}\left\\{{\bf U}_{{\mathbf{f}}}{\bf
S}_{{\mathbf{f}}}{\bf
V}^{T}_{{\mathbf{f}}}\frac{d{\mathbf{x}}(\theta)}{d\theta}\right\\},$

where ${\bf J}_{{\mathbf{f}}}$ is the Jacobian matrix of ${\mathbf{f}}$ and
${\bf U}_{{\mathbf{f}}}{\bf S}_{{\mathbf{f}}}{\bf V}_{{\mathbf{f}}}^{T}$ is
its singular value decomposition, SVD (we have omitted here their dependency
on ${\mathbf{x}}(\theta)$ for brevity). On the other hand, from Eq. (30) we
can write the expected local covariance matrix ${\bf C}(\theta)$ of the
features fluctuations, as:

$\displaystyle{\bf C}(\theta)$ $\displaystyle\approx$
$\displaystyle\sigma_{d}^{2}\,\mathbb{E}\\{{\bf J}_{{\mathbf{f}}}{\bf
J}_{{\mathbf{f}}}^{T}\\}$ (32) $\displaystyle=$
$\displaystyle\sigma_{d}^{2}\mathbb{E}\\{{\bf U}_{{\mathbf{f}}}{\bf
S}_{{\mathbf{f}}}^{2}{\bf U}_{{\mathbf{f}}}^{T}\\}.$

Here it is crucial to notice that, under the assumptions made in previous and
current subsections, whereas ${\bf J}_{{\mathbf{f}}}({\mathbf{x}}(\theta))$
will heavily depend on ${\mathbf{x}}$, ${\bf
J}_{{\mathbf{f}}}({\mathbf{x}}(\theta)){\bf
J}_{{\mathbf{f}}}^{T}({\mathbf{x}}(\theta))$ will be much less sensitive to
${\mathbf{x}}$, as it only depends on the inner products of the different
features’ gradients (see Eq. (30)). Furthermore, in Subsection 6.2 we show how
these inner products (at least their correlation factor, which depends only on
their relative angle) are fairly stable, especially when inputs are samples
from pdf’s. Therefore, the ${\bf U}_{{\mathbf{f}}}({\mathbf{x}}(\theta))$ and
${\bf S}_{{\mathbf{f}}}({\mathbf{x}}(\theta))$ matrices, on their average
behavior, will determine both the direction of change of $\bar{\mathbf{f}}$
when changing $\theta$ (Eq. (31)) and the dominant direction of ${\bf
C}(\theta)$ (Eq. (32)), which, thus, will tend to coincide. Our observations
indicate that the eigenvalues of ${\bf C}(\theta)$ (the diagonal terms of
${\bf S}_{{\mathbf{f}}}^{2}$) are fairly concentrated in the studied cases.
This implies that the features’ coupling actually causes a worst case scenario
for discriminating between similar $\theta$’s: the features’ pdfs become (i)
elongated (due to eigenvalues’ concentration), and (ii) locally aligned with
the $\bar{\mathbf{f}}(\theta)$ curve. This causes strong overlapping of the
pdf’s having close $\theta$ values, and, as a result, poor discrimination.
Fig. 7(left) illustrates this phenomenon in a real experiment with real data.
Fig. 7(right) shows the effect of decoupling the kurtosis from the skewness
(orthokurtosis). We used 128 random vectors of 1024 i.i.d. samples each, from
$x(\theta)=x_{0}^{\theta}$, being $x_{0}\sim U(0,1)$. Ellipses correspond to a
Mahalanobis radius of 2, and $\theta=5$ (black), 6 (blue), and 7 (red).
Expected error probabilities, using a bi-variate Gaussian model, are 12.4%
(coupled) vs. 4.5% (decoupled).

|
---|---
Figure 7: Comparing Gaussian classes in the original feature space (left) and
in the decoupled feature space (right). See text for details.

To conclude this section, the techniques presented here attack the core of the
poor discrimination problem due to using coupled features, by orthogonalizing
the feature gradients, which has the effect of approximately diagonalizing the
local covariance matrix ${\bf C}(\theta)$. It is crucial to note that this
diagonalization is effective because it is local. A global diagonalization
(such as the classical Principal Component Analysis, PCA) would be useless for
reducing the pdfs’ overlapping corresponding to close $\theta$ values in the
feature space, as such overlapping is insensitive to global affine
transformations of that space. In contrast, a global linear correlation (as
the one shown in Figure 7(b) ) can be trivially removed, if needed, by
applying PCA after feature decoupling.

##  6 Experiments and Applications

###  6.1 Using two families of global features

In this section we define the two features’ families (sets) that will be
studied in the experiments, namely: (i) marginal moments of arbitrary order
$p$ (MM); (ii) $p$-th order moment at the output of a 2-D filter bank defined
from a filter $h$ (MF). Specifically, to extract the MF features, we first
applied the Translation Invariant Laplacian separable (TILs) representation
[33], a tight frame acting as bar and edge detector that provides nine
subbands, each with the same number of coefficients as pixels in the image.
Then, the $p$-th order moment was obtained for every subband. To obtain the
corresponding decoupled sets of the MM and MF families
($\text{DF}_{\text{MM}}$ and $\text{DF}_{\text{MF}}$ respectively), we used
the Nested Normalization-narrow path (Algorithm 3). For the MM set of
features, we used as reference values the $p$-th moment of a standardized
Gaussian distribution (i.e., $(p-1)!!$ for even $p$, 0 for odd [34]). For the
MF set of features, reference values corresponded to moments obtained by
convolving zero-mean univariate white Gaussian noise with a kernel $h$, which,
in our case, using use the set of kernels $\\{h_{j},j=1\dots 9\\}$ from the
TILs representation, are the same function of $p$ as for the MM family in all
subbands. Table II shows the original features for MM and MF. It also shows
their corresponding gradients expressions (ignoring scaling factors which do
not influence the result) and indicates in which cases the set of original
gradients fulfills the Frobenius condition.

TABLE II: Families of features (marginal moments, MM, and moments at the output of filters, MF), and their gradients. Family | $f({\mathbf{x}})$ | $\partial f({\mathbf{x}})/\partial x_{i}$ | Frobenius
---|---|---|---
MM | $\displaystyle 1/N\sum_{n=1}^{N}x_{n}^{p}$ | $x_{i}^{(p-1)}$ | For all $p$
MF | $\displaystyle 1/N\sum_{n=1}^{N}(x\ast h)_{n}^{p}$ | $(x\ast h)^{(p-1)}\ast\tilde{h}$ | Only for $p=2$

###  6.2 Measuring the amount of mutual coupling

In this section we evaluate how close to being mutually orthogonal are the
feature’s gradients, for two sets of coupled features and their corresponding
decoupled sets, namely: (i) a set composed of the first six orders of the
classical MM, in its standardized version: mean, variance, and the rest the
moments of the standardized sample to zero mean and unit variance (that is,
skewness, kurtosis, etc.). We will refer to this classical set of statistical
features by “MSM” (from Marginal Standardized Moments), and
$\text{DF}_{\text{MSM}}$ its corresponding decoupled set; (ii) a set composed
of the second-order moments at the output of a filter bank (“VF”, a particular
case of MF with $p=2$ and assuming zero mean) and its corresponding decoupled
set ($\text{DF}_{\text{VF}}$).

Let $\\{f_{j}({\mathbf{x}})\\}$ represent the original features and
$\\{\hat{f}_{j}({\mathbf{x}})\\}$ its corresponding decoupled set. Note that
$j=1,\dots,6$ for MSM and $j=1,\dots,9$ for VF (for the 9 subbands of the TILs
representation). Let $\mathbf{x_{0}}$ represent an $N$-D vector of i.i.d.
samples drawn from a Gaussian or an uniform distribution; or an $N$-D vector
that represents the pixel values of a texture patch extracted from an image of
the Broadtz database [35]. To obtain the gradient of a feature at
$\mathbf{x_{0}}$, we numerically calculated the partial derivatives with
respect to the $i$-th variable $x_{i}\in{\mathbf{x}}$ by adding a differential
perturbation $\epsilon$ to the $i$-th element of vector ${\mathbf{x}}$:

$\frac{\partial f(\mathbf{x_{0}})}{\partial x_{i}}=\lim_{\epsilon\to
0}\frac{f(\mathbf{x_{0}}+\epsilon\mathbf{e_{i}})-f(\mathbf{x_{0}}-\epsilon\mathbf{e_{i}})}{2\epsilon}.$
(33)

This expression yields the gradients for each feature of the MSM, the VF, and
their corresponding decoupled sets ($\text{DF}_{\text{MSM}}$ and
$\text{DF}_{\text{VF}}$ respectively). To evaluate the function
$\hat{f}_{j}({\mathbf{x}})$ in the DFs cases we used the Nested Normalization-
narrow path (Algorithm 3) using equation (4.1.2) for a fast analytical
computation of the gradient’s orthogonal projections.

Then we measured the angle $\alpha$ between pairs of gradient vectors of the
different features that belong to the MSM and VF sets
$\\{f_{j}({\mathbf{x}})\\}$ and that belong to the $\text{DF}_{\text{MSM}}$
and $\text{DF}_{\text{VF}}$ sets $\\{\hat{f}_{j}({\mathbf{x}})\\}$. The
deviation from orthogonality (DO) was obtained as the difference between 90
degrees (perfect orthogonality) and the actual calculated angles
($\text{DO}=90-\alpha$). As such, $\text{DO}>0$ indicates acute angles and
$\text{DO}<0$ obtuse angles. The number of samples were $N$=512 for the
Gaussian and uniform distributions and $N$=529 (23$\times$23 pixels) for
textures. The experiment was repeated $M$=256 times for the Gaussian and
uniform distributions. In the case where $\mathbf{x_{0}}$ came from textures,
we used a single patch for each of the $M$=112 different textures in the
Brodatz database. Table III shows the average of the absolute value of the DO
across the different pairs of feature’s gradients for the different
distributions tested, for the MSMs, VFs, $\text{DF}_{\text{MSM}}$ and
$\text{DF}_{\text{VF}}$. We excluded from the average calculation the mean and
the variance in the MSM and $\text{DF}_{\text{MSM}}$ cases, as they are
orthogonal by definition. Figures 8 and 9 show the DO results. Panels (a) and
(b) show the DO between different pairs of gradient’s features for the
original set and its decoupled counterpart (Gaussian and Textures cases
respectively). Panels (c-f) show the DO, in absolute value, between different
pairs of gradients. Panels (c) and (d) show results for the Gaussian case;
panels (e) and (f) show the results obtained for the Textures case. Blue color
indicates $|\text{DO}|$ close to 0 degrees (orthogonality), while yellow color
indicates a deviation from orthogonality close to 90 degrees (angle of 0
degrees). Note that the main diagonal only acts as a reference ($0$ degrees)
here.

####  6.2.1 Marginal moments

Let us analyze first the case of the marginal moments. In agreement with the
theory (the first three decoupled moments have gradients fulfilling Frobenius
condition, see Proposition 3.1) we see that the DO is exactly zero (perfect
orthogonality) for all standardized moments combined with orders 1 and 2 (note
that the MSM set is already a partially decoupled version of the original MM
family set), and for all decoupled moments combined with orders 1, 2 and 3. In
addition, it is close to zero on average in the rest of odd-even combinations,
in which case it presents a wide variance for MSM, and a much narrower one for
$\text{DF}_{\text{MSM}}$. For the odd-odd and even-even rest of the cases, MSM
presents highly acute angles (strong average coupling) with an extremely low
variance, whereas DF provides either perfect (3-5 case) or approximate (4-6
case) orthogonality, also with low variance. In summary, the decoupled moments
$\text{DF}_{\text{MSM}}$ provide close to orthogonal gradients also for orders
greater than 3, especially for the random sampling experiments (samples are
close to their expected values, in the reference manifolds, where exact
orthogonality holds), but also, to a lesser extent and with higher
variability, for photographic images (textures). In contrast, MSM presented in
those cases either very high deviations (with low variability) or low-to-
moderate average deviations with high variability (especially in real images).
In Table III we can see that, for real photographic images, the average
absolute DO has been reduced in a factor 6, approx., whereas for Gaussian
samples it has been reduced 16 times.

####  6.2.2 Second-order moments at the output of a filter bank

Figure 8 shows the empirical results obtained using the VF features. First, we
note that original features are fairly coupled, although not as much as in the
marginal moments’ case. Now the new $\text{DF}_{\text{VF}}$ features are
exactly and universally decoupled only for the eight pairs
$(\hat{f}_{1},\hat{f}_{j}),j=2\dots 9$, for which the theory tells us that the
gradients of the new features achieve perfect orthogonality when Frobenius
condition holds for the gradients of the original features (as it happens in
this case; see Proposition 3.1 and the note about the special case $\hat{\cal
S}_{1}={\cal S}_{1}$). Although this particular condition on the first eight
pairs of gradients is difficult to appreciate in Fig. 8(a) and (d), where
deviation from orthogonality seems approximately zero for all pairs of
features using white Gaussian samples, we measured a RMS value of the
deviation from orthogonality for the first eight couples of $1.6\times
10^{-5}$ degrees. This is a negligible numerical error exclusively due to the
numerical computation of the gradient (see Eq. (33)). For the rest of feature
couples we obtained an RMS of $0.42$ degrees, still small, but four orders of
magnitude larger. Aside from the first eight pairs of features, the excellent
practical decoupling of the other ones has been favored in this case by using
a random distribution (zero-mean uni-variate white Gaussian noise) for the
samples, with features whose expected values are precisely the reference
values used in the NeN decoupling algorithm. As explained in the theory
(Subsection 3.1.3), exact decoupling is also obtained for samples on the
reference manifolds, even if Frobenius condition does not hold for the output
features. For large enough samples, sample feature values do not deviate much
from their expected values, and, as a consequence of the features’ smoothness,
the gradients of these samples will also be close to orthogonal. To test how
the decoupling quality degrades when using real samples instead of pseudo-
random ones, we have tested the method, again, with the referred collection of
112 textured $23\times 23$ pixel patches. We first note that the amount of
mutual coupling between features is almost exactly the same as for the white
noise case, a relevant fact that adds support to the assumptions made in
Section 5: the angle between features is fairly constant for each couple of
features, both within a given distribution/collection (see the relatively
small amplitude of error bars) and also across different
distributions/collections (compare panel (a) with (b) and (c) with (e)). We
can say that, in this case, it is especially accurate to say that the local
covariance matrix is fairly independent of ${\mathbf{x}}$, behaving the
Jacobian pretty close to a moving frame. As for the decoupling quality for
this collection of real photographic patches, we now observe than in about
half of the pairs the decoupling is either perfect (first eight couples) or
almost perfect. For the other half, DO is still very moderate and much smaller
than for the original pairs, in the great majority of cases. In Table III we
can see that, for real photographic images, the average absolute DO has been
reduced in a factor 6, approx., whereas for Gaussian samples it has been
reduced 65 times.

TABLE III:  Average ($\pm\sigma$) absolute deviation from orthogonality ($|\text{DO}|$), in degrees, between feature’s gradients for the MSM, VF, $\text{DF}_{\text{MSM}}$ and $\text{DF}_{\text{VF}}$ sets of features. | Gaussian | Textures
---|---|---
MSM | 32$\pm$26 | 59$\pm$21
$\text{DF}_{\text{MSM}}$ | 2$\pm$3 | 9$\pm$17
VF | 13$\pm$10 | 12$\pm$11
$\text{DF}_{\text{VF}}$ | 0.2$\pm$0.3 | 2$\pm$3

(a) DO between different pairs of gradients. Gaussian case.

(b) DO between different pairs of gradients. Textures case.

(c) MSM set. Gaussian case.

(d) $\text{DF}_{\text{MSM}}$ set. Gaussian case.

(e) MSM set. Textures case.

(f) $\text{DF}_{\text{MSM}}$. Textures case.

Figure 8: Deviation from orthogonality (DO), in degrees, for the different
pairs of gradient’s features for the standardized marginal moments (MSM) and
the corresponding decoupled sets ($\text{DF}_{\text{MSM}}$). (a), (c) and (d)
show results for the Gaussian case, $N=512$. (b), (e) and (f) show results for
the Textures case, $N=529$ (23 $\times$ 23 pixels). In (a)-(b) positive values
indicate acute angle and negative obtuse. In (c-f) the color level represents
the average of the absolute value of DO.

(a) DO between different pairs of gradients. White Gaussian noise case.

(b) DO between different pairs of gradients. Textures case.

(c) VF, White Gaussian case.

(d) $\text{DF}_{\text{VF}}$, White Gaussian case.

(e) VF, Textures case.

(f) $\text{DF}_{\text{VF}}$, Textures case.

Figure 9: Deviation from orthogonality (DO), in degrees, for the different
pairs of gradient’s features for the second-order moment at the output of a
filter bank (VF), and the corresponding decoupled set
($\text{DF}_{\text{VF}}$). We used $N=529$ (23 $\times$ 23 pixels). (a), (c)
and (d) show results for the Gaussian case. (b), (e) and (f) show results for
the Textures case. In (a)-(b) positive values indicate acute angle and
negative obtuse. In (c-f) the color level represents the average of the
absolute value of DO.

###  6.3 Statistical regression

In this section we focus on the estimation of the parameters of a distribution
that best describe a dataset999Some results of this subsection have been
presented in [7].. Different approaches have been proposed in the literature
for different families of parametric distributions, such as the classical
maximum likelihood estimation or the method of moments. We approached the
estimation task as a regression problem, where the parameters of the
distribution are estimated from a set of global features obtained from the
observed data. Specifically, we tested, as global features: (i) a set of
classical MSM, and (ii) the corresponding decoupled set
($\text{DF}_{\text{MSM}}$). In order to thoroughly compare the descriptive
capabilities of both sets of features, we used several regression methods,
namely: linear regression models (LRM), regression trees (RT), support vector
regression (SVR), Gaussian process regression (GPR), ensembles of trees (ET)
and neural networks (NNR). All of these methods are implemented in the
Regression Learner App, ®Matlab. For reproducibility purposes, we used the
hyper-parameters set by default in the referred app. Specific information
about implementation, hyper-parameters selection and methodological details
can be found in [36].

In our experiments we used different statistical distributions, specifically:
generalized Gaussian distribution (GGD), Gamma distribution (GMD), and
absolute value of a Normal distribution raised to a positive number (GND). The
shape of these distributions depends on a shape parameter ($\beta$), and the
regression problem consists in estimating $\beta$ from an observed data set.
See Appendix D for the expressions of the probability density functions of
these distributions and their dependence with $\beta$.

Let ${\mathbf{x}}$ represent an $N$-D vector of i.i.d. samples drawn from a
GGD, GMD or GND distributions (we generated the samples following [37, 38] for
the GGD, and [39] for the GMD, using the ®Matlab function gamrnd.m),
normalized to have zero mean and unitary variance. For the GGD and GMD cases,
given that the kurtosis of these distributions changes very fast for small
values of the $\beta$ parameter, we defined $\beta=2^{A}$ and sampled
uniformly the exponent $A$ in $[-3,3]$, resulting in $\beta$ ranging in
$[1/8,8]$. In this way, we obtained a quite uniform distribution of kurtosis
values. In the GND case, $\beta$ was sampled uniformly in the range $[1,6]$.

MSM features and their corresponding decoupled features,
$\text{DF}_{\text{MSM}}$, were obtained from the third (skewness) to the sixth
order. This led to two sets of 4 dimensional predictors
$\\{f_{j}({\mathbf{x}}),j=3,\dots,6\\}$ and
$\\{\hat{f}_{j}({\mathbf{x}}),j=3,\dots,6\\}$ for each parameter $\beta$. We
compared the $\beta$ prediction accuracy of these two sets. We generated
$d=2048$ vectors ${\mathbf{x}}$ of different lengths $N$ of i.i.d. samples,
thus having $2048$ pairs of predictors ($\\{f_{j}\\}$ or
$\\{\hat{f}_{j}({\mathbf{x}})\\}$) and targets (known $\beta$ values). We
averaged 100 5-fold cross-validation runs to measure the accuracy of the
methods in terms of the RMSE in the estimation of the exponent
$A=\log_{2}\text{$\beta$}$ in the GGD and GMD cases, and the estimation of
$\beta$, in the GND case.

(a) GGD distribution

(b) GMD distribution

(c) GND distribution

Figure 10: RMSE of the estimated parameters as a function of the number of
samples $N$ for the NNR regression method applied to three density functions
(see text for details).

Figure 10 shows the RMSE of the results as a function of $N$ for the MSM and
$\text{DF}_{\text{MSM}}$ sets of descriptors, by using the NNR regression
method, the method providing the best results for almost all the tested cases
(see Appendix D for more detailed results). Figure 10(a) shows the results for
GGD, Figure 10(b) for GMD and Figure 10(c) for GND. The shadow area represents
the standard deviation across 100 repetitions of the experiment. The proposed
$\text{DF}_{\text{MSM}}$ clearly outperformed MSM for all the tested
distributions and sample sizes. For instance, in the GGD case, using the
proposed $\text{DF}_{\text{MSM}}$ the RMSE was reduced by factors of $0.50$,
$0.54$, $0.63$, $0.68$, $0.76$, and $0.86$, for $N=2048,1024,512,256,128$, and
$64$, respectively.

Tables VI, VII and VIII in Appendix D show the RMSE obtained using MSM and DF
descriptors for the different sample sizes $N$ and regression methods, for the
GGD, GMD and GND distributions, respectively. The best regression method for
each $N$ is highlighted in bold. Our decoupled descriptors
$\text{DF}_{\text{MSM}}$ outperformed the MSM across all the compared
regression methods and sample sizes in 105 out of 108 cases ($97.2\%$), which
shows the robustness and generalization ability of our approach.

###  6.4 Texture classification

In this section we apply the proposed method for texture classification in two
different settings: (i) comparing standardized moments (MSM) with their
corresponding decoupled features $\text{DF}_{\text{MSM}}$, both applied to a
set of subbands, the output of a filter bank (Section 6.4.1); (ii) using
features that directly are defined as marginal moments at the output of a
filter bank (MF), and their corresponding fully decoupled features,
$\text{DF}_{\text{MF}}$ (Section 6.4.2).

####  6.4.1 Texture classification based on marginally decoupled moments in a
filter bank

Here101010This subsection is a summary of the results published in [6]. we
compare the performance of two classifiers using features derived (i) from MSM
of order 2nd to 10th, (ii) from a generalized form of Ref. [16], and (iii)
from the corresponding decoupled features $\text{DF}_{\text{MSM}}$, all of
them at the output of subbands of the TILs representation [33]. We selected 54
textures of 640$\times$640 pixels from the Broadtz [35] database under the
criterion of looking homogeneous in 64$\times$64 pixel patches111111The list
of selected textures and the ®Matlab code of the experiments is available in
https://www.researchgate.net/project/Nested-Normalizations-for-Decoupling-
Global-Features.. Every texture was divided into $10\times 10$ disjoint
$64\times 64$ patches. The problem consisted in classifying patches in their
corresponding textures. To extract the features we first applied the TILs
representation [33] and discarded the low-pass band (having 8 subbands). Then,
for every subband of every patch, three sets of features were obtained: (i)
classical MSM features $\\{f_{j}({\mathbf{x}}),j=2,\dots,10\\}$ (something
commonly used to characterize textures, but to a lower order, see e.g.[40]);
(ii) modified moment set: same as MSM, except that now even order moments are
shift-minimized
(${\tilde{\mu}_{n}}=\min_{\alpha_{n}}{\mathbb{E}}\\{({\hat{x}}-\alpha_{n})^{n}\\}$,
${\hat{x}}$ being the standardized observation) [16]; and (iii) proposed
marginally decoupled moments $\text{DF}_{\text{MSM}}$
$\\{\hat{f}_{j}({\mathbf{x}}),j=2,\dots,10\\}$.

In order to quantify the redundancy between features, we estimated the mutual
information (MI) between pairs of features [41]. Some mean values across
subbands are shown in Table IV.

TABLE IV:  Averaged mutual information. Features | (3,4) | (3,5) | (4,5) | (3,6) | (4,6) | (5,6)
---|---|---|---|---|---|---
MSM | 0.26 | 0.61 | 1.39 | 0.55 | 3.10 | 3.11
Mod. [16] | 0.23 | 0.61 | 1.39 | 0.55 | 3.07 | 3.03
$\text{DF}_{\text{MSM}}$ | 0.22 | 0.04 | 0.08 | 0.09 | 0.54 | 0.09

We see that, except for the (3,4)-th case, the proposed
$\text{DF}_{\text{MSM}}$ features present a drastic reduction of redundancy
compared to the classical MSM and modified [16] ones. We recall the reader
that the aim of using decoupled features is not to compensate for the
statistics of the data (to which the decoupling method is totally
transparent), like non-linear ICA aimes for, but rather to avoid adding
spurious coupling in the processed data, leaving only the dependencies that
are effectively caused by the input data statistics. So, perfectly decoupled
features will generally reduce the mutual information among features, but will
not necessarily remove it, as that depends on the data statistics (see, e.g.,
what happened in the example of Fig. 7). Figure 11 shows a 3-D subset of the
MSM (left) and proposed $\text{DF}_{\text{MSM}}$ (right) features (i.e.,
before and after our decoupling method) for two texture classes, to illustrate
the decoupling effect on the data distribution in the feature space. Shown
data correspond to a single subband (number 5 in the representation), of