# MrSARP: A Hierarchical Deep Generative Prior for SAR Image Super-resolution

Tushar Agarwal, Nithin Sugavanam and Emre Ertin T. Agarwal, Nithin Sugavanam,
and E. Ertin are with the Department of Electrical and Computer Engineering,
The Ohio State University, Columbus, OH, 43210 USACorresponding Author: T.
Agarwal (agarwal.270@buckeyemail.osu.edu)This work has been submitted to the
IEEE for possible publication. Copyright may be transferred without notice,
after which this version may no longer be accessible.

###### Abstract

Generative models learned from training using deep learning methods can be
used as priors in inverse under-determined inverse problems, including imaging
from sparse set of measurements. In this paper, we present a novel
hierarchical deep-generative model MrSARP for SAR imagery that can synthesize
SAR images of a target at different resolutions jointly. MrSARP is trained in
conjuction with a critic that scores multi resolution images jointly to decide
if they are realistic images of a target at diffferent resolutions. We show
how this deep generative model can be used to retrieve the high spatial
resolution image from low resolution images of the same target. The cost
function of the generator is modified to improve its capability to retrieve
the input parameters for a given set of resolution images. We evaluate the
model’s performance using the three standard error metrics used for evaluating
super-resolution performance on simulated data and compare it to upsampling
and sparsity based image sharpening approaches.

###### Index Terms:

Deep Learning, Super-Resolution, Compressive sensing.

##  I Introduction

Synthetic aperture radar (SAR) imagery captures the physical aspects of the
target differently compared to the Optical imagery because of multi-path
reflections, specular nature of reflectors, and Imaging geometry effects
leading to overlay, shadowing. In this work, we present a generative model
that captures SAR phenomenology at different resolutions exhibited at SAR
magnitude imagery. Specifically, we propose a specialized hierarchical
architecture of the generative model, called MrSARP, that jointly models the
data manifold of multiple-resolutions of a SAR image. We show that such a
generative model acts as the projection operator to a lower dimensional
manifold and can be directly used for super-resolving magnitude SAR images
from low resolution magnitude images. The super-resolution performance of
MrSARP is evaluated by comparing it with LASSO (or $L_{1}$ recovery) and
Nearest Neighbor Upsampling using empirical data.

Next, we describe important notation used throughout the text. The lower-case
character $x$ represents a fixed sample value of a random vector (rv) $X$ (the
corresponding upper-case character), $M:N$ is the enumeration of natural
numbers from $M$ to $N$, $g(X_{1:N};w)$ denotes function $g$ of rv’s $X_{1:N}$
with deterministic parameters $w$ and an index-based $x[v]$ represents the
$v^{th}$ value from the ordered set $\\{x[i]:i\in\\{1:N\\}\\}$. Formal
definitions of terms Generative Modeling and Compressed Sensing, as used in
the context of this paper, are stated below.

######  Definition 1 (Generative Modeling).

Generative Modeling aims to estimate the distribution $p(X)$ using a
parameterized distribution family $q(X;w)$ and a set $\mathcal{D}$ of samples
from the distribution $p(X)$. The stochasticity in sampling is typically
achieved by sampling a low-dimensional latent rv $Z$, internal to $q(X;w)$,
from a simple distribution (e.g. standard gaussian). Therefore, the goal is to
generate $N$ realistic samples $\\{\hat{X}[v]\\}_{v=1}^{N}$ from $q(X;w)$ as
if they were from $p(X)$.

######  Definition 2 (Compressed Sensing).

Given measurements $y\in\mathbb{C}^{d}$ obtained using a known measurement (or
forward) operator $F$ of an underlying signal $x\in\mathbb{C}^{D}$, Compressed
Sensing seeks to recover this underlying signal under the model

$\displaystyle y=F(x)+\eta,\qquad x\in\Omega,$ (1)

where $\eta$ is the measurement noise and $\Omega$ represents the constraint
set that is non-convex. It is an under-determined system of equations, i.e.,
the no. of measurements are less than the signal dimension and so the signal
structure must be exploited through appropriate constraints to obtain a unique
solution for the problem:

$\displaystyle\min_{x\in\Omega}\qquad$ $\displaystyle L(F({x});{y})$ (2)

Where $L(F({x});{y})$ is an appropriate loss function to be minimized.

We want to use the range space of the generative model $q(X;w)$ as the
constraint set $\Omega$ that can super-resolve SAR data as well as potentially
be used to solve a compressed sensing problem for SAR data. Specifically, we
want to find a prior function as an ANN based generative model $G$ with
parameters $w_{G}$, that generates random samples of $X$ (SAR magnitude
images) using a low-dimensional latent rv $Z$, i.e. sampling from $G(Z;w_{G})$
is a good approximation of sampling from $p(X)$. This implicitly constrains
samples $X$ on a low-dimensional manifold while having the flexibility to
adapt the basis to any dataset, unlike commonly used sparsity priors.
Therefore, $G$ can be used as the learned projection function
$\mathcal{P}_{\Omega}$ in a Projected Gradient Descent optimization to find
the solution in the constraint set $\Omega$. We show that a specially designed
$G$ can itself be used for super-resolving SAR images. The dataset of all
samples is split into 3 subsets $\mathcal{D}_{train}$, $\mathcal{D}_{val}$ and
$\mathcal{D}_{test}$ for model training, cross-validation and evaluation of
$G$ respectively.

###  I-A Related Work

A survey of trained and untrained DL methods to solve the Inverse problems is
given in [1]. Data-driven models can be broadly classified as end-to-end
models that are agnostic to the measurement operator and deep generative
priors that model the distribution representing the constraint set estimated
from a corpus of data. We’ll discuss the approaches based on the latter
briefly.

Theoretical guarantees are established on the number of measurements and
conditions on the generator network for successful reconstruction in [2]. A
generative model is estimated on a large corpus of natural images in [3].
Given the generator, which is implemented as a generative adversarial network
or variational auto-encoder, and the measurements we can infer the image by
solving (2). The projection operator $\mathcal{P}_{\Omega}$ for the set of
natural images $\Omega$ is learned from a large corpus of data using an
adversarial network in [4]. This projection operator is used in solving any
inverse problem in a plug and play manner since the projection operator has
enough capacity to model the complicated non-convex set of natural images.
Theoretical guarantees for the generative prior are established in [5]. It is
shown that if the entries of the generator function is near-random and the
number of weights in each layer increases with the width then the loss
function landscape contains descent directions to the global optima. The deep-
geometric prior is extended to MR imaging in [6]. Invertible generators are
used. Since the dimensionality of the latent code is large and equivalent to
the signal dimension, a block-wise structure is imposed on the structure of
the latent variable. The idea of generative prior is extended to video-
sequences in a novel way to synthesize video from a sub-sample of image frames
in [7]. The image sequences in time are denoted by
$X=\mathbf{x}_{1},\mathbf{x}_{2},\cdots,\mathbf{x}_{T}$. It is shown that if
all the images are obtained from the same generator than the network
parameters $\theta$ can be fixed and for each image in the sequence a latent
code can be estimated to succinctly represent the video. The latent code
sequence $Z=[\mathbf{z}_{1},\mathbf{z}_{2},\cdots,\mathbf{z}_{T}]$ compactly
represents the image. The smoothness in the image sequences of a video can be
imposed by imposing a smoothness constraint or a low-rank constraint on the
latent code $Z$ while jointly estimating in the prediction step. Furthermore
if images are dissimilar then the network parameters $\theta$ can also be re-
trained to capture the variability. It is shown that missing frames can be
synthesized using interpolation in the latent space. Next, we consider the
regime when no or limited training data is available, where the generator is
estimated for each image.

Deep image prior [8] presents a non-trained version of the generator discussed
in previous section. The latent code is assumed fixed and can be chosen
arbitrarily but the network parameters are optimized to represent the image.
It is hypothesized that the structure of the network imposes a strong
regularizer or prior on the image. This strong prior is shown to have high
impedance to noise and uncorrelated samples. Therefore, the optimization
problem solved is

$\displaystyle\min_{\theta}L(g(\mathbf{z},\theta),\mathbf{y},F).$ (3)

The main problem is over-fitting it is shown that early-stopping is necessary
to capture the de-noised image. This number of iterations serves as a hyper-
parameter that needs to be carefully chosen. Different neural network
architectures are explored and it is empirically shown that as the network is
over-parameterized the performance improves because the capacity of network is
flexible to learn the image.

The works in [9], and [10] propose an alternative loss function termed as the
backprojection loss and demonstrate the efficacy theoretically as well as
empirically for multiple imaging linear inverse problems such as de-blurring
and super-resolution as a function of the condition-number of the measurement
operator. The loss function considered is

$\displaystyle L(\theta)=\lVert
F^{+}(\mathbf{y}-Fg(\mathbf{x},\theta))\rVert^{2},$ (4)

where $F^{+}$ denotes the pseudo-inverse of the linear measurement operator.
Theoretical guarantees for compressive sensing are established for deep priors
for [11], and [12]. It is shown that the number of measurements required are
similar to the traditional requirements established in compressive sensing.
Empirically, these methods are shown to perform better than $\ell_{1}$ and
$TV$ norm regularization on fastMRI dataset.

##  II Method

Building upon eq. 2, squared Euclidean $L_{2}$ norm is commonly used for the
loss function, i.e., $L(F({x});{y})=\lVert{y}-F({x})\rVert^{2}_{2}$ under the
assumption that $\eta$ is additive Gaussian noise. Sparsity in some known
basis $\Phi$, such as Fourier or Wavelet basis, is the most widely used
constraint on $x$, achieved by adding $L_{1}$ regularization to sparse vector
$\tilde{x}$ where $\tilde{x}$ is such that $x=\Phi\tilde{x}$. This derivation
involves a convex relaxation as well as $\Phi$ must follow the restricted
isometry property and [13] can be referred for details. Instead of this
sparsity assumption, we follow a different approach. The problem in eq. 2 can
also be solved by finding a convex relaxation of the constraint set $\Omega$
and using a projected gradient method as

$\displaystyle\hat{x}_{mp}$
$\displaystyle=\arg\,\min_{x}\lVert{y}-F({x})\rVert^{2}_{2}$
$\displaystyle=\arg\,\min_{m,p}\lVert{y}-F({m\odot\exp(jp)})\rVert^{2}_{2}$
(5) $\displaystyle\hat{x}_{\Omega}$
$\displaystyle=\mathcal{P}_{\Omega}(\hat{x}_{mp})=\arg\,\min_{x\in\Omega}\lVert\hat{x}_{mp}-x\rVert^{2}_{2}$
(6)

Where $m,p$ are magnitude and phase of complex-valued $x$ respectively
alongwith subscripts (if any), $\mathcal{P}_{\Omega}$ denotes the projection
operator onto set $\Omega$ Alternating between these two steps will find the
desired solution if $\Omega$ and $F$ were convex starting from an appropriate
initial condition $x_{0}$. However, like most $x$ of interest, we don’t know
$\Omega$ for SAR data. Therefore, similar to Projected Gradient Descent GAN by
Shah and Hegde [14], we propose to learn $\Omega$ as the range-space of an ANN
based Generative model $G(Z;w_{G})$. Here, $Z$ as much lower dimensionality
than $x$. Moreover, we only constrain the magnitude $m$ of $x$ using
$G(Z;w_{G})$ and allow phase $p$ to be unconstrained without projection.
Assuming that we have obtained such a $G(Z;w_{G})$ that well approximates the
probability distribution $P(m)$, projection $\hat{x}_{\Omega}$ becomes

$\displaystyle\hat{x}_{\Omega}$
$\displaystyle=\mathcal{P}_{\Omega}(\hat{m}_{mp}\odot\exp(j\hat{p}_{mp}))$
$\displaystyle=\mathcal{P}_{G}(\hat{m}_{mp})\odot\exp(j\hat{p}_{mp})$
$\displaystyle=\hat{m}_{G}\odot\exp(j\hat{p}_{mp})$ (7)

Where $\hat{m}_{G}$ is the projection of magnitude $\hat{m}_{mp}$ derived
using $G$. Here, the magnitude $m$ of the structured signal $x$ is assumed to
be sampled from a low-dimensional manifold with a latent variable $z$. To find
$\hat{m}_{G}$, we first optimize over $z$ starting from initial condition
$z_{0}=0$ as

$\displaystyle z^{*}_{z}$
$\displaystyle=\arg\,\min_{z}\lVert\hat{m}_{mp}-G(z;w_{G})\rVert^{2}_{2}$ (8)

However, since $G$ is a non-linear ANN, range-space of $G(Z;w_{G})$ is non-
convex. Therefore, solving optimization problem in eq. 8 only guarantees
local-minima. Hence, it will be sensitive to initialization and may yield
inconsistent projections. Motivated from works of Bojanowski et al. [15] and
Wu et al. [16], we attempt to alleviate this problem by making $G(Z;w_{G})$
aware of this inversion task during its training phase. This is described in
section III.

Additionally, we borrow important ideas from similar work of IAGAN [17]. Even
after significant progress in Deep Generative Modeling, such models are still
an approximation of the true distribution due to limited representation
capabilities of ANN. Therefore, magnitudes $m$ of many samples $x$ may not
even belong to the range-space of $G(Z;w_{G})$. To mitigate this, IAGAN
proposes to do image-adaptive projections, i.e. optimizing over both latent
vector $z$ and ANN weights $w$ when projecting. This adds the following step
to eq. 8 as

$\displaystyle(z^{*}_{zw},w^{*}_{G})$
$\displaystyle=\arg\,\min_{z,w}\lVert\hat{m}_{mp}-G(z;w)\rVert^{2}_{2}$ (9)

Where $z$, $w$ are initialized as $z_{0}=z^{*}_{z}$, $w_{0}=w_{G}$
respectively while optimizing eq. 9. Finally, the magnitude projection
required for eq. 7 is
$\hat{m}_{G}=\mathcal{P}_{G}(\hat{m}_{mp})=G(z^{*}_{zw};w^{*}_{G})$. We call
our complete algorithm as Image Adaptive Projected Gradient Descent WGAN and
it is mentioned in algorithm 1.

Algorithm 1 IAPGD WGAN

1:$n,n_{mp},n_{z},n_{zw}\geq 0,y,F,G(z;w_{G}),x_{0}$

2:while $n\neq 0$ do

3:
$\hat{x}_{mp}\leftarrow\left[{\arg\,\min_{m,p}\lVert{y}-F({m\odot\exp(jp)})\rVert^{2}_{2}}\right]_{n_{mp}}$

4: $\hat{m}_{mp}\leftarrow\absolutevalue{\hat{x}_{mp}}$

5: $z_{0}\leftarrow 0$

6:
$z^{*}_{z}\leftarrow\left[{\arg\,\min_{z}\lVert\hat{m}_{mp}-G(z;w_{G})\rVert^{2}_{2}}\right]_{n_{z}}$

7: $(z_{0},w_{0})\leftarrow(z^{*}_{z},w_{G})$

8:
$(z^{*}_{zw},w^{*}_{G})\leftarrow\left[{\arg\,\min_{z,w}\lVert\hat{m}_{mp}-G(z;w)\rVert^{2}_{2}}\right]_{n_{zw}}$

9: $\hat{m}_{G}\leftarrow G(z^{*}_{zw};w^{*}_{G})$

10: $\hat{x}_{\Omega}\leftarrow\hat{m}_{G}\odot\exp(j\hat{p}_{mp})$

11: $x_{0}\leftarrow\hat{x}_{\Omega}$

12: $n\leftarrow n-1$

13:end while

Where the notation $\left[{\arg\,\min_{x}\mathcal{L}(x)}\right]_{n_{x}}$
refers to minimizing the loss $\mathcal{L}(x)$ for $n_{x}$ steps iteratively
using a variant of Stochastic Gradient Descent (SGD).

We propose a special hierarchical architecture of a WGAN that jointly
generates magnitude of multiple resolutions of the same SAR image. This idea
of jointly modeling multiple resolutions was inspired from the Progressive GAN
by Karras et al. [18] though their motivations were different. We aim to
exploit this hierarchical structure for super-resolution i.e. to find a higher
resolution image given its lower resolutions. Suppose we are given a dataset
containing 4 exponentially increasing resolution images’ magnitude
$m^{r_{1}},m^{r_{2}},m^{r_{3}},m^{r_{4}}$ where resolution of $m^{r_{i+1}}$ is
twice of $m^{r_{i}}$. Then our WGAN $G(Z;w_{G})$ models the joint probability
distribution $P(m^{r_{1}},m^{r_{2}},m^{r_{3}},m^{r_{4}})$. Now if we are given
a new sample from $P(m^{r_{1}},m^{r_{2}},m^{r_{3}})$ i.e. of the 3 lower
resolutions, we can use $G(Z;w_{G})$ and steps 5 to 9 of algorithm 1 to
project it onto the joint data manifold of
$P(m^{r_{1}},m^{r_{2}},m^{r_{3}},m^{r_{4}})$ by finding a common
$(z^{*}_{zw},w^{*}_{G})$ pair. The highest resolution $r_{4}$ image is then
obtained by a simple forward pass
$\hat{m}^{r_{4}}_{G}=G(z^{*}_{zw},w^{*}_{G})$. The overall algorithm is
therefore Projecting from a Multi-Resolution SAR Prior (or MrSARP) and is
summarized in algorithm 2. Note that $G(z;w_{G})^{r_{i}}$ denotes ${r_{i}}$
resolution output from $G(z;w_{G})$ and $G(z;w_{G})^{r_{4}}$ is simply dropped
in step 3.

Algorithm 2 Projecting from MrSARP

1:$n_{z},n_{zw}\geq 0,m^{r_{1}},m^{r_{2}},m^{r_{3}},G(z;w_{G})$

2:$z_{0}\leftarrow 0$

3:$z^{*}_{z}\leftarrow\left[{\arg\,\min_{z}\Sigma_{i=1}^{3}\lVert
m^{r_{i}}-G(z;w_{G})^{r_{i}}\rVert^{2}_{2}}\right]_{n_{z}}$

4:$(z_{0},w_{0})\leftarrow(z^{*}_{z},w_{G})$

5:$(z^{*}_{zw},w^{*}_{G})\leftarrow\left[{\arg\,\min_{z,w}\Sigma_{i=1}^{3}\lVert
m^{r_{i}}-G(z;w)^{r_{i}}\rVert^{2}_{2}}\right]_{n_{zw}}$

6:$\hat{m}^{r_{4}}_{G}\leftarrow G(z^{*}_{zw};w^{*}_{G})^{r_{4}}$

##  III Model

(a) Generator

(b) Critic

Figure 1: Hierarchical Model Architecture of MrSARP

MrSARP consists of a WGAN-GP [19] with a Hierarchical architecture. Our
architecture is inspired from the CIFAR-10 ResNet architecture used by
Gulrajani et al. [19] and ProGAN [18]. We use the ResNet block from the former
and the idea of FromImage/ToImage layers at various ResNet input/output
features from latter. The schematic of the hierarchical model is shown in
figure 1. The output shapes of layers are specified. The ResNet blocks in the
generator use Batch-Normalization layers and have nearest-neighbour upsampling
after the input. The ResNet blocks in the critic uses Layer-Normalization
layers and have an average pooling based downsampling before the final output,
and channel-wise concatenation in case of multiple inputs. The FromImage layer
is a 2D convolution layer with 1 channel input (image) and 2 channel output
(features). The ToImage layer comprises slicing first 2 channels of input
features, Batch-Normalization, ReLU activation and 2D convolution layer with 1
channel output (image) after tanh activation. All convolutional layers in the
generator have kernel-size as 3 and ReLU activation function. All
convolutional layers in the critic have kernel-size as 3 and LeakyReLU
activation function with 0.2 slope. The dense (or fully-connected) layers in
both, generator and critic have linear activation and the GAP referes to
Global Average Pooling operation. The generator $G(Z;w_{G})$ of MrSARP
generates samples from the unconditional joint distribution of 4 resolutions
of SAR images $P(m^{r_{1}},m^{r_{2}},m^{r_{3}},m^{r_{4}})$ where resolution of
$m_{r_{i}}$ is $2^{i+3}\times 2^{i+3}$ and bandwidth is $125.2^{i}Mhz.$.

WGAN-GP was chosen because of its advantages over traditional GANs, especially
meaningful loss curves for cross-validation and reliable training as
demonstrated by Gulrajani et al. [19]. WGAN-GP, like all GANs, require an
extra critic ANN $D(\cdot,\cdot;w_{D})$ to be learned simultaneously to aid
the learning of the generator. Both, the critic and generator are trained
alternately using gradient descent steps. Latent variable
$Z\sim\mathcal{N}(0,I)$ as in a conventional GAN.

Since MrSARP would primarily be used for inverting and projecting on data-
manifold instead of sampling, we think it is essential to inform the generator
$G$ about this task during training. Inspired from Wu et al’s work [16], we
add Model Agnostic Meta Learning (MAML) based regularizing loss term to the
overall loss functions used for training $G$. Proposed by Finn et al. [20],
MAML is a general-purpose meta-learning method to adapt parameters $w$ of a
statistical model to a number of tasks as long as the loss function
$\mathcal{L}(\mathcal{T};w)$ for the task $\mathcal{T}$ is differentiable.
Since our primary task is inversion using step 3 of algorithm 2, we perform a
small number $n_{z}=5$ iterations of this optimization step to get $z^{*}$ and
use residual error on all 4 channels as our MAML loss $\mathcal{L}_{MAML}$.
Since the $Z\sim\mathcal{N}(0,I)$ with high dimensionality $d_{Z}$, the
samples lie near the hypersphere of $\sqrt{d_{Z}}$. To enforce such a
constraint on $z^{*}$, we additionally project $z^{*}$ on to the
$\sqrt{d_{Z}}$ hypersphere after every SGD step using projection operator
$\mathcal{P_{S}}$. The complete loss functions are as follows.

$\displaystyle\mathcal{L}$
$\displaystyle=\begin{cases}-\mathcal{L}_{W}+\lambda_{1}\mathcal{L}_{GP}&\text{for
}D\\\ \mathcal{L}_{W}+\lambda_{2}\mathcal{L}_{MAML}&\text{for }G\end{cases}$
(10)

$\displaystyle\mathcal{L}_{W}$
$\displaystyle=\operatorname{\mathbb{E}}_{m}\left[{D(m);w_{D}}\right]-\operatorname{\mathbb{E}}_{Z}\left[{D(G(Z;w_{G});w_{D})}\right]$
$\displaystyle\mathcal{L}_{GP}$
$\displaystyle=\operatorname{\mathbb{E}}_{\bar{m}}(\left\lVert{\nabla_{\bar{m}}D(\bar{m};w_{D})}\right\rVert_{2}-1)^{2}$
$\displaystyle\mathcal{P_{S}}(z^{*})$
$\displaystyle=\frac{z^{*}}{\min(\left\lVert{z^{*}}\right\rVert_{2},\sqrt{d_{Z}})}.\sqrt{d_{Z}}$
(11) $\displaystyle z^{*}_{S}$
$\displaystyle=\left[{\mathcal{P_{S}}\left({\arg\,\min_{z}\Sigma_{i=1}^{3}\lVert
m^{r_{i}}-G(z;w_{G})^{r_{i}}\rVert^{2}_{2}}\right)}\right]$ (12)
$\displaystyle\mathcal{L}_{MAML}$
$\displaystyle=\operatorname{\mathbb{E}}_{m}\left[{\Sigma_{i=1}^{4}\lVert
m^{r_{i}}-G(z^{*}_{S};w_{G})^{r_{i}}\rVert^{2}_{2}}\right]$ (13)

Where $m\sim P(m^{r_{1}},m^{r_{2}},m^{r_{3}},m^{r_{4}})$, $\bar{m}\sim
r(\bar{M})$ is uniform-sampling along straight lines between pairs of sampled
points $m$, $G(z;w_{G})$. All expectations are approximated using
corresponding empirical means. Motivated by [19], we employ negative critic
loss as our primary metric for model selection. Hence, the weights of the
final model are set to their values corresponding to the epoch where the
smallest negative critic loss on $\mathcal{D}_{val}$ was achieved. For the
hyperparameters, we use $\lambda_{1}=10$ (as in [19]) and set
$\lambda_{2}=1000$ to approximately balance the two regularizing loss term
$L_{MAML}$ with the WGAN loss term $\mathcal{L}_{W}$ in magnitude.

##  IV Experimental Setup

We use the Tensorflow (2.1) [21] deep-learning Python library for our ANN
implementations. Unless mentioned otherwise, the Adam optimizer from
Tensorflow with default parameters is used for most optimizations. We use
NVIDIA GeForce RTX 2080 Ti GPU alongwith Intel Xeon CPU as our primary
computation hardware.

###  IV-A CVDomes Dataset

We use the phase histories obtained from simulation of back-scattered energy
from civilian vehicles in [22]. We consider only the HH polarization
measurement in our experiment. The image at a bandwidth B and azimuth span
$\Delta\theta$ is obtained by backprojection method using a Hamming window to
suppress the side-lobes. We spotlight on a square patch of $9m\times 9m$. To
generate 4 different resolutions $r{1:4}$, we repeat this process for 4
bandwidths $125.2^{i}$ and corresponding pixel resolutions $2^{i+3}\times
2^{i+3}$ for $i\in{1,2,3,4}$.

###  IV-B Pre-Processing

We work with magnitude images only for MrSARP. Hence, we first find the
absolute values of complex-valued SAR images and then perform min-max
normalization of every image individually to restrict their values in range
$[-1,1]$. We then upsample all lower resolution images to the highest
resolution of $128\times 128$ with nearest-neighbor upsampling and combine all
4 resolutions into a single 4 channel image. The lower resolution images are
appropriately downsampled using average pooling before their input into the
critic. This is only done to simplify multi-resolution architecture
implementation and is based on the simple fact that nearest-neighbor
upsampling followed by average-pooling gives identity function.

##  V Results

We present some results of this study in this section. After training
generator $G$ as described in section III above, we use algorithm 2 to find
the super-resolved images $\hat{m}_{G}^{r_{4}}$ for the unseen lower
resolution samples in $\mathcal{D}_{test}$. We compare these images with the
available ground truth images $m^{r_{4}}$ qualitatively as well as using three
quantitative metrics viz. PSNR (Peak Signal-to-Noise Ratio), NMSE (Normalised
Mean Squared Error) as defined in [23] and SSIM (Structural Similarity Index)
proposed in [24]. The use of these metrics is motivated from the evaluation
schemes used in existing literature on super-resolution of SAR images e.g.
[25]. These metric used are defined as follows.

$\displaystyle MSE$
$\displaystyle=\frac{1}{N_{m}}\left\lVert{m^{r_{4}}-\hat{m}^{r_{4}}}\right\rVert_{F}^{2}$
$\displaystyle PSNR$
$\displaystyle=10\log_{10}\left({\frac{(\max(m^{r_{4}})-\min(m^{r_{4}}))^{2}}{MSE}}\right)$
(14) $\displaystyle NMSE$
$\displaystyle=\frac{\left\lVert{m^{r_{4}}-\hat{m}^{r_{4}}}\right\rVert_{F}^{2}}{\left\lVert{m^{r_{4}}}\right\rVert_{F}^{2}}$
(15) $\displaystyle SSIM$ $\displaystyle=l(x,y).c(x,y).s(x,y)$ (16)
$\displaystyle l(x,y)$
$\displaystyle=\left({\frac{2\mu_{x}\mu_{y}+C_{1}}{\mu_{x}^{2}+\mu_{y}^{2}+C_{1}}}\right)$
$\displaystyle c(x,y)$
$\displaystyle=\left({\frac{2\sigma_{x}\sigma_{y}+C_{2}}{\sigma_{x}^{2}+\sigma_{y}^{2}+C_{2}}}\right)$
$\displaystyle s(x,y)$
$\displaystyle=\left({\frac{\sigma_{xy}+C_{3}}{\sigma_{x}\sigma_{y}+C_{3}}}\right)$

where $N_{m}$ are total number of pixels in image $m^{r_{4}}$,
$\mu_{x},\mu_{y}$ are empirical means of patches $x,y$ respectively,
$\sigma_{x},\sigma_{y}$ are sample standard deviations of patches $x,y$
respectively, $\sigma_{xy}$ is sample cross-correlation of $x,y$ after mean
subtraction and $C_{1:3}$ are small constants added for numerical stability.
SSIM is calculated on smaller local patches $x,y$ of size $7\times 7$ and the
mean SSIM is calculated for every image comparison pair
$m^{r_{4}},\hat{m}^{r_{4}}$.

We also compare $r_{4}$ resolution estimate $\hat{m}_{N}^{r_{4}}$ obtained
from nearest neighbor upsampling of $m^{r_{3}}$ and estimates
$\hat{m}_{L}^{r_{4}}$ obtained from using the popular LASSO method. The LASSO
method is used to obtain the sparse scattering center representation of the
vehicle at the resolution ${r_{3}}$. The sparse representation is projected
back into the phase history measurement domain using the SAR forward operator.
These measurements are converted to SAR imagery using the backprojection
method with the hamming window at resolution ${r_{4}}$. All the evaluations
are done on 504 images from the unseen $\mathcal{D}_{test}$.

Figure 2 shows some samples from $\mathcal{D}_{test}$ that are super-resolved
using different methods. There is a qualitative similarity between the Nearest
Neighbor upsampled and $L_{1}$ recovered images. Both methods differ
significantly from resulting images of MrSARP.

(a) Image 1

(b) Image 2

(c) Image 3

Figure 2: Qualitative comparison of samples super-resolved using different
methods

These quantitative results are presented in table I. Addition of Image-
Adaptive steps to MrSARP improves performance quantitatively for all 3
metrics. In fact, among all the methods tested, Image-Adaptive MrSARP results
perform best in terms of both, $NMSE$ and $PSNR$. However, they are inferior
to both Nearest Neighbor Upsampling and $L_{1}$ recovery in terms of SSIM
indicating potential bias in amplitude for MrSARP.

TABLE I: Quantitative evaluation of Super-resolution performance. Method | $NMSE\downarrow$ | $PSNR\uparrow$ | $SSIM\uparrow$
---|---|---|---
Nearest-Neighbor Upsampling | 9.071 | 27.27 | 0.931
$L_{1}$ recovery | 8.698 | 27.638 | 0.937
MrSARP | 6.92 | 29.785 | 0.9
Image-Adaptive MrSARP | 6.003 | 30.963 | 0.918

##  VI Conclusion

In this paper, we showed how a GAN with hierarchical architecture can jointly
model the distribution of multiple resolutions of magnitude SAR images. We
further showed how this GAN, called MrSARP, be used to for super-resolving SAR
images. We saw some improvements over baselines of nearest neighbor upsampling
as well as $L_{1}$ recovery in terms of $NMSE$ and $PSNR$ values but there is
significant scope of further improvements in terms of perceptual quality. This
was indicated by MrSARP’s inferior performance in terms of SSIM. Furthermore,
we plan to utilize this generative model and the proposed algorithm 1 to
regularize inverse problems in SAR imaging with structured interrupted
measurements.

## Acknowledgements

This research was partially supported by NSF grants CNS-1823070 CBET-2037398
and NIH Grant P41EB028242

## References

  * [1] G. Ongie, A. Jalal, C. A. M. R. G. Baraniuk, A. G. Dimakis, and R. Willett, “Deep learning techniques for inverse problems in imaging,” _IEEE Journal on Selected Areas in Information Theory_ , 2020.
  * [2] Z. Liu and J. Scarlett, “Information-theoretic lower bounds for compressive sensing with generative models,” _IEEE Journal on Selected Areas in Information Theory_ , 2020.
  * [3] A. Bora, A. Jalal, E. Price, and A. G. Dimakis, “Compressed sensing using generative models,” in _International Conference on Machine Learning_ , 2017, pp. 537–546.
  * [4] J. H. R. Chang, C. Li, B. Póczos, B. V. K. Vijaya Kumar, and A. C. Sankaranarayanan, “One network to solve them all — Solving linear inverse problems using deep projection models,” in _2017 IEEE International Conference on Computer Vision (ICCV)_ , 2017, pp. 5889–5898.
  * [5] P. Hand and V. Voroninski, “Global guarantees for enforcing deep generative priors by empirical risk,” _IEEE Transactions on Information Theory_ , vol. 66, no. 1, pp. 401–418, 2019.
  * [6] V. A. Kelkar, S. Bhadra, and M. A. Anastasio, “Compressible latent-space invertible networks for generative model-constrained image reconstruction,” _arXiv preprint arXiv:2007.02462_ , 2020.
  * [7] R. Hyder and M. S. Asif, “Generative models for low-dimensional video representation and reconstruction,” _IEEE Transactions on Signal Processing_ , vol. 68, pp. 1688–1701, 2020.
  * [8] V. Lempitsky, A. Vedaldi, and D. Ulyanov, “Deep image prior,” in _2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2018, pp. 9446–9454.
  * [9] T. Tirer and R. Giryes, “Back-projection based fidelity term for ill-posed linear inverse problems,” _IEEE Transactions on Image Processing_ , vol. 29, pp. 6164–6179, 2020.
  * [10] J. Zukerman, T. Tirer, and R. Giryes, “BP-DIP: A backprojection based deep image prior,” _arXiv preprint arXiv:2003.05417_ , 2020.
  * [11] R. Heckel and M. Soltanolkotabi, “Compressive sensing with un-trained neural networks: Gradient descent finds the smoothest approximation,” _arXiv preprint arXiv:2005.03991_ , 2020.
  * [12] R. Heckel, “Regularizing linear inverse problems with convolutional neural networks,” _arXiv preprint arXiv:1907.03100_ , 2019.
  * [13] R. G. Baraniuk, V. Cevher, M. F. Duarte, and C. Hegde, “Model-Based Compressive Sensing,” _IEEE Transactions on Information Theory_ , vol. 56, no. 4, pp. 1982–2001, Apr. 2010. [Online]. Available: http://arxiv.org/abs/0808.3572
  * [14] V. Shah and C. Hegde, “Solving Linear Inverse Problems Using Gan Priors: An Algorithm with Provable Guarantees,” in _2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , Apr. 2018, pp. 4609–4613.
  * [15] P. Bojanowski, A. Joulin, D. Lopez-Paz, and A. Szlam, “Optimizing the Latent Space of Generative Networks,” May 2019. [Online]. Available: http://arxiv.org/abs/1707.05776
  * [16] Y. Wu, M. Rosca, and T. Lillicrap, “Deep Compressed Sensing,” May 2019. [Online]. Available: http://arxiv.org/abs/1905.06723
  * [17] S. A. Hussein, T. Tirer, and R. Giryes, “Image-Adaptive GAN based Reconstruction,” Nov. 2019. [Online]. Available: http://arxiv.org/abs/1906.05284
  * [18] T. Karras, T. Aila, S. Laine, and J. Lehtinen, “Progressive Growing of GANs for Improved Quality, Stability, and Variation,” Feb. 2018\. [Online]. Available: http://arxiv.org/abs/1710.10196
  * [19] I. Gulrajani, F. Ahmed, M. Arjovsky, V. Dumoulin, and A. Courville, “Improved Training of Wasserstein GANs,” _arXiv:1704.00028 [cs, stat]_ , Dec. 2017. [Online]. Available: http://arxiv.org/abs/1704.00028
  * [20] C. Finn, P. Abbeel, and S. Levine, “Model-Agnostic Meta-Learning for Fast Adaptation of Deep Networks,” _arXiv:1703.03400 [cs]_ , Jul. 2017. [Online]. Available: http://arxiv.org/abs/1703.03400
  * [21] M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard, M. Kudlur, J. Levenberg, R. Monga, S. Moore, D. G. Murray, B. Steiner, P. Tucker, V. Vasudevan, P. Warden, M. Wicke, Y. Yu, and X. Zheng, “{}TensorFlow{}: A System for {}Large-Scale{} Machine Learning,” in _12th USENIX Symposium on Operating Systems Design and Implementation (OSDI 16)_ , 2016, pp. 265–283. [Online]. Available: https://www.usenix.org/conference/osdi16/technical-sessions/presentation/abadi
  * [22] K. E. Dungan, C. Austin, J. Nehrbass, and L. C. Potter, “Civilian vehicle radar data domes,” in _Algorithms for Synthetic Aperture Radar Imagery XVII_ , vol. 7699. SPIE, Apr. 2010, pp. 242–253. [Online]. Available: http://www.spiedigitallibrary.org/conference-proceedings-of-spie/7699/76990P/Civilian-vehicle-radar-data-domes/10.1117/12.850151.full
  * [23] N. Karimi and M. R. Taban, “Nonparametric blind SAR image super resolution based on combination of the compressive sensing and sparse priors,” _Journal of Visual Communication and Image Representation_ , vol. 55, pp. 853–865, 2018.
  * [24] Z. Wang, A. Bovik, H. Sheikh, and E. Simoncelli, “Image quality assessment: From error visibility to structural similarity,” _IEEE Transactions on Image Processing_ , vol. 13, no. 4, pp. 600–612, Apr. 2004.
  * [25] C. He, L. Liu, L. Xu, M. Liu, and M. Liao, “Learning Based Compressed Sensing for SAR Image Super-Resolution,” _IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing_ , vol. 5, no. 4, pp. 1272–1281, Aug. 2012.


# Functional data learning using convolutional neural networks

J Galarza, T Oraby  School of Mathematical and Statistical Sciences, The
University of Texas Rio Grande Valley, Edinburg, TX 78539, USA

###### Abstract

In this paper, we show how convolutional neural networks (CNN) can be used in
regression and classification learning problems of noisy and non-noisy
functional data. The main idea is to transform the functional data into a 28
by 28 image. We use a specific but typical architecture of a convolutional
neural network to perform all the regression exercises of parameter estimation
and functional form classification. First, we use some functional case studies
of functional data with and without random noise to showcase the strength of
the new method. In particular, we use it to estimate exponential growth and
decay rates, the bandwidths of sine and cosine functions, and the magnitudes
and widths of curve peaks. We also use it to classify the monotonicity and
curvatures of functional data, algebraic versus exponential growth, and the
number of peaks of functional data. Second, we apply the same convolutional
neural networks to Lyapunov exponent estimation in noisy and non-noisy chaotic
data, in estimating rates of disease transmission from epidemic curves, and in
detecting the similarity of drug dissolution profiles. Finally, we apply the
method to real-life data to detect Parkinson’s disease patients in a
classification problem. The method, although simple, shows high accuracy and
is promising for future use in engineering and medical applications.

  * September 2023

Keywords: Functional Data Learning, Deep Learning, Convolutional Neural
Networks, Regression, Classification

##  1 Introduction

Functional data (FD) are functions observed for each unit over certain
intervals; see (?, ?). FD appears in various scientific fields, such as
engineering, geology, biology, medicine, pharmacology, and chemistry. It
involves analyzing data in the form of continuous vector functions or curves,
which could be treated as realizations of stochastic processes; see (?, ?).
Functional data analysis (FDA) provides methods for extracting intrinsic
information from infinite-dimensional and irregular observation data; see (?,
?). FDA combines statistics, spatial analysis, and multivariate modeling tools
to analyze and predict functional data. It provides advantages over
traditional pointwise estimation methods by using irregularly sampled data in
space, time, and depth to fit space-time functional models (?).

Earlier stages of FDA were developed by (?) to study the relationship between
the ability of an examinee and his or her probability of correctly selecting
an option in a standardized item response test. First, (?) used kernels to
correctly fit the observations to one dimension and later extended the idea in
(?) to fit the data to larger dimensions. Then, (?) introduced Principal
Differential Analysis to find an approximate solution that solves the
differential equation $Lu=0$ in which we can pick the order of the
differential operator $L$ and a basis for the weights, using splines or other
adequate functions. This technique was applied to approximate Chinese writing
in (?) using a second-order operator. Later, complexity was added to this type
of multilevel model (?). For more information and applications of FD, we
recommend reading the full work on Functional Analysis case studies of
Ramsay’s work in (?) and finding more mathematical foundation in (?).

Functional data learning (FDL) has also received a lot of attention recently
in the domain of machine learning and deep learning. Different approaches were
used to establish predictive models of functional data like in (?) where
gradient descent was derived for classical Neural Networks using Fréchet
derivates and the Reisz representation theorem to establish a deep learning
method for FDL. In (?), the support vector machine method was applied to
feature vectors obtained by voxels transformation. In (?), some FDA R
libraries were compared to machine learning techniques such as random forest
and were shown to be outperformed by the latter. In (?), a new methodology was
used to learn the bases of different functional subspaces to model functional
data before applying learning methods. In (?), orthogonal B-spline bases for
FDA were shown to be more efficient than other Fourier-based methods.

Other deep learning and machine learning approaches include the work in (?),
in which inputs are fed directly into a layer composed of nodes of nested
neural networks. The output of these neural networks is the basis function
themselves. Inspired by the square root velocity, the characteristics in (?)
were obtained from functional data using derivatives. It was found to work
well along with classification methods such as logistic regression and support
vector machine. Other methods of FDL such as manifold learning have recently
been introduced, see, e.g. (?, ?, ?).

Convolutional neural network (CNN) has been used extensively in image
recognition. The architecture of Neural Networks has been evolving since the
introduction of very deep CNNs (?), ResNets (?), and MobileNets (?, ?). The
efficiency of CNN is affected by its architecture and hardware (?), as well as
training and cross-validation sets and scaling (?). Several other advances
have been made by adding other algorithms to CNN, such as the ”You Only Look
Once” approach proposed by (?, ?). To our knowledge, convolutional neural
networks have not been used before for functional data learning.

In this paper, we explore how to use CNNs in functional data learning. In
Section 2, we introduce our method and the CNN architecture. In most of the
problems, we use the same CNN architecture except for a few that will be
described in place of their application. In Section 3, we start with some
regression and classification problems of curves that represent different
functions with and without random noise. In particular, we examine the ability
of CNN to estimate parameters of exponential and trigonometric functions. In
addition, we employed the same CNN to estimate the magnitude and width of the
curve peaks. We also use CNN for the classification problems of increasing
versus decreasing function, concave versus convex function, algebraic versus
exponential growth, and discerning curves of 0, 1, and 2 peaks. In Section 4,
we examine the ability of CNN to estimate parameters of characteristics of
dynamical systems. In particular, we use it to estimate Lyapunov exponents of
some chaotic curves with and without random noise. We apply the same CNN to
estimate transmission rates from epidemic curves. Then, we applied a Siamese
CNN to test the similarity between two drug dissolution curves (profiles). In
Section 5, we apply the same methods to classify the actual functional data of
cases of Parkinson’s disease and control.

##  2 Methods

Let $\\{(x_{i},f(x_{i})):i=1,2,\ldots,n\\}$ be the data points of a graph with
equidistance sampled points $x_{i}$, $i=1,2,\ldots,n$. The data pre-processing
step in functional data learning via CNNs creates an input image using
functional data points. We assume that $f$ is Min-Max normalized, so we can
assume that $f:\mathbb{R}\to[0,1]$. A signed distance matrix $D$ is defined,
so its elements $(i,j)$, $d_{ij}=f(x_{i})-f(x_{j})$. The matrix $D$ represents
a grayscale intensity value that is used to produce a 28 by 28 grayscale
image; see Figure 1.

Next, we describe a typical architecture as presented in Mathworks’ documents
(?), which is based on Lenet networks in (?). The first multilayer of CNN
consists of batch normalization, a RELU activation function, and an average
normalization layer. The first output is $13\times 13\times 8$ layers. The
second multilayer is the same as the first with an output of $5\times 5\times
16$ layers. The third multilayer does not contain an average pooling layer and
has an output of $3\times 3\times 32$ layers. The final multilayer has no
average pooling layer; however, it has a dropout layer of 20% with a fully
connected layer followed by a regression/classification layer. All codes and
simulations were performed in MATLAB and using its deep learning toolbox (?).

Figure 1: Diagram and procedure proposed for the regression and classification
problems.

We use the same CNN’s architecture for all of the regression and
classification problems. We tested the procedure on different functional types
and we start first by discussing our findings for seven case studies of
functional data.

##  3 Results

In this section, we discuss various regression and classification problems of
functional data with and without random noise. We use randomly generated
values for the functional parameters from specified ranges to produce
training, validation, and testing data sets of sizes 1000, 100, and 100,
respectively. We use $n=100$ data points for each curve. In the case of width
estimation, height estimation, and number of peaks classification, we use
$10000$ curves for training, and we use $n=1000$ functional data points for
each curve. This would ensure that there will be a wide variety of heights and
widths and enough peaks for each classification.

###  3.1 Regression Problems

####  3.1.1 Exponential Function

Let the exponential curve be $y=\exp(\omega x)$ with parameter $\omega$
representing the rate at which $y$ increases or decreases. The first task is
to use our method to estimate $\omega$. See Figure 2 for the curve of
$y=\exp(-0.27x)$ with and without noise as an example and the corresponding 28
by 28 image. For training, validation, and testing data, the parameter
$\omega$ is sampled from a uniform distribution over the interval $[-1,1]$ for
$1000$, $100$, and $100$ times, respectively.

Figure 2: (a) The curve of $y=\exp(\omega x)$ when $\omega=-.27$. (b) The 28
by 28 image that corresponds to the curve in (a). (c) The curve of
$y=\exp(\omega x)+\sigma z$ when $\omega=-.27$ and $\sigma=.1$ where $z$ is a
standard normal random variable. (d) The 28 by 28 image that corresponds to
the curve in (c). (f) The curve of $y=\exp(\omega x)+\sigma z$ when
$\omega=-.27$ and $\sigma=1$ where $z$ is a standard normal random variable.
(f) The 28 by 28 image that corresponds to the curve in (e).

Figure 3 shows the predicted values versus the estimated value of the
exponential function parameter that closely follows the diagonal line with no
intercept and slope of one. Table 1 shows a strong diagonal linear
relationship between the true and predicted values of the rate.

Figure 3: Results of the exponential function regression using the test dataset without noise (a), with noise of magnitude $\sigma=.1$ (b), and with noise of magnitude $\sigma=1$ (c). Table 1: Correlation coefficient, intercept, and slope of a simple regression line of the actual and predicted values, and their p-values for three cases. Case | Correlation Coefficient ($r$) | Intercept (p-value) | Slope (p-value)
---|---|---|---
Without Noise | 0.999 | -0.000 (0.867) | 0.994 (0.049)
With Noise ($\sigma=0.1$) | 0.995 | 0.002 (0.692) | 0.995 (0.657)
With Noise ($\sigma=1$) | 0.970 | -0.011 (0.415) | 1.012 (0.655)

####  3.1.2 Sine and Cosine Functions

Consider the functions $y=\sin(\omega x)$ and $y=\cos(\omega x)$ with the
frequency parameter $\omega$. See some example curves in Figures 4 and 5 for
$y=\sin(1.06x)$ and $y=\cos(0.96x)$, respectively. For training, validation,
and testing data, the parameter $\omega$ is randomly uniformly sampled from
the interval $[0,3]$ for $1000$, $100$, and $100$ times, respectively.

Figure 4: (a) The curve of $y=\sin(\omega x)$ when $\omega=1.06$. (b) The 28
by 28 image that corresponds to the curve in (a). (c) The curve of
$y=\sin(\omega x)+\sigma z$ when $\omega=1.06$ and $\sigma=.1$ where $z$ is a
standard normal random variable. (d) The 28 by 28 image corresponds to the
curve in (c). (e) The curve of $y=\sin(\omega x)+\sigma z$ when $\omega=1.06$
and $\sigma=1$ where $z$ is a standard normal random variable. (f) The 28 by
28 image that corresponds to the curve in (e).

Figure 5: (a) The curve of $y=\cos(\omega x)$ when $\omega=.96$. (b) The 28 by
28 image that corresponds to the curve in (a). (c) The curve of $y=\cos(\omega
x)+\sigma z$ when $\omega=.96$ and $\sigma=.1$ where $z$ is standard normal
random variable. (d) The 28 by 28 image corresponds to the curve in (c). (f)
The curve of $y=\cos(\omega x)+\sigma z$ when $\omega=.96$ and $\sigma=1$
where $z$ is a standard normal random variable. (e) The 28 by 28 image that
corresponds to the curve in (e).

Figure 6 shows the predicted values versus the estimated value of the
parameter of the sine and cosine functions that closely follow the diagonal
line without an intercept and slope of one. Tables 2 and 3 show a strong
diagonal linear relationship between the true and predicted bandwidth values.

Figure 6: Results of the sine function regression using the test dataset without noise (a), with noise of magnitude $\sigma=.1$ (c), and with noise of magnitude $\sigma=1$ (e). Results of the cosine function regression using the test dataset without noise (b), with noise of magnitude $\sigma=.1$ (d), and with noise of magnitude $\sigma=1$ (f). Table 2: Correlation Coefficient, Intercept, and Slope for Sine Data with P-values Case | Correlation Coefficient ($r$) | Intercept (p-value) | Slope (p-value)
---|---|---|---
Without Noise | $\approx 1$ | 0.0015 (0.8306) | 0.9954 (0.9121)
With Noise ($\sigma=0.1$) | $\approx 1$ | 0.0171 (0.0836) | 0.9911 (0.4942)
With Noise ($\sigma=1$) | $\approx 1$ | 0.0011 (0.9649) | 1.0108 (0.2765)
Table 3: Correlation Coefficient, Intercept, and Slope for Cosine Data with P-values Case | Correlation Coefficient ($r$) | Intercept (p-value) | Slope (p-value)
---|---|---|---
Without Noise | $\approx 1$ | -0.0059 (0.3546) | 0.9979 (0.5825)
With Noise ($\sigma=0.1$) | $\approx 1$ | 0.0143 (0.3497) | 1.0045 (0.6084)
With Noise ($\sigma=1$) | $\approx 1$ | 0.0639 (0.0305) | 0.9647 (0.0259)

####  3.1.3 Estimation of the Height and Width of Peaks of Functions

To produce a number of peaks with different heights and widths, we used a
mixture of Gaussian curves, see Figures 7. The maximum height of a curve is
the height of the peak at the global maximum of the curve. The width of a
curve is the horizontal distance of the contour located half of the prominence
of that peak. It is important to note that the curves need to be normalized
with the highest peak of all the generated curves to have grayscale values
representative of all of the heights. Otherwise, the highest grayscale value
of 1 will be assigned to the peak height of each curve and thus will not be
able to train the model properly. On the other hand, with the width estimation
and peak detection (see peak classification section), we use local
normalization, since they do not necessarily depend on the height.

To generate the curves, we use a mixture of Gaussian kernels given by
$y=G(x)=\sum^{2}_{k=1}H_{k}\exp({-(\frac{x-P_{k}}{W_{k}}})^{2})$ for
$x\in[0,50]$. Then the heights $H_{k}$ are randomly and uniformly selected
from $[0,2200]$. The shape parameters $W_{k}$ are randomly generated using
$W_{k}=\lfloor 50\cdot U_{1}+1\rfloor$ where $U_{1}$ is a uniformly
distributed random variable on [0,1]. Similarly, position parameters are
randomly generated using $P_{k}=\lfloor 50\cdot U_{2}+1\rfloor$ where $U_{2}$
is a uniformly distributed random variable over [0,1]. See Figure 7 for an
example of a mixture of Gaussian kernels with and without random noise.

Figure 7: (a) The curve is an example of a mixture of Gaussians $G(x)$. (b)
The 28 by 28 image that corresponds to the curve in (a). (c) The curve of
$G(x)+\sigma z$ when $\sigma=.1$ where $z$ is a standard normal random
variable. (d) The 28 by 28 image that corresponds to the curve in (c). (e) The
curve of $G(x)+\sigma z$ when $\sigma=1$ where $z$ is a standard normal random
variable. (f) The 28 by 28 image that corresponds to the curve in (e).

Figure 8 shows that the normalized width predictions closely follow the actual
values. In this part, Table 4 shows a strong indication of accurate prediction
of the width of the peaks of the curves. Table 5 supported by Figure 9,
however, shows that the CNN could not retrieve the actual values of height on
average unless it has more noise.

Figure 8: Results of the maximum width regression using the test dataset
without noise (a), with noise of magnitude $\sigma=.1$ (b), and with noise of
magnitude $\sigma=1$ (c).

Figure 9: Results of the maximum height regression using the test dataset without noise (a), with noise of magnitude $\sigma=.1$ (b), and with noise of magnitude $\sigma=1$ (c). Table 4: Correlation Coefficient, Intercept, and Slope for Width Estimation Data with P-values Case | Correlation Coefficient ($r$) | Intercept (p-value) | Slope (p-value)
---|---|---|---
Without Noise | 0.974 | -0.0017 (0.6247) | 1.0057 (0.4453)
With Noise ($\sigma=0.1$) | 0.980 | 0.0022 (0.4370) | 0.9964 (0.5736)
With Noise ($\sigma=1$) | 0.977 | -0.0071 (0.0178) | 1.0138 (0.0463)
Table 5: Correlation Coefficient, Intercept, and Slope for Height Estimation Data with P-values Case | Correlation Coefficient ($r$) | Intercept (p-value) | Slope (p-value)
---|---|---|---
Without Noise | 0.951 | -0.0169 (0.0013) | 1.1195 (0.0000)
With Noise ($\sigma=0.1$) | 0.941 | 0.0084 (0.1244) | 1.1044 (0.0000)
With Noise ($\sigma=1$) | 0.932 | -0.0057 (0.3473) | 1.0327 (0.0105)

###  3.2 Classification Problems

In this subsection, we examine the capabilities of CNN in functional data
classification. For this type of problem, the same CNN architecture is used
except that after the drop-out layer, the fully connected layer has hidden
nodes equal to the number of classes, followed by a softmax layer followed by
a classification Layer.

####  3.2.1 Increasing versus Decreasing Curves

We use CNN to classify the monotonicity of curves. Curves
$y=e^{w_{1}(x-w_{2})}$ are used to generate increasing or decreasing
exponential curves for training, validation, and training datasets, see Figure
10. We also use the following random variables $w_{1}=sign(U_{1}-.5)$ and
$w_{2}=2U_{2}+2.5$ where $U_{1}$ and $U_{2}$ are uniformly distributed random
variables on [0,1]. We found that the accuracy of the classification of the
functional test data is 100%.

Figure 10: (a)The curve of $y=e^{-(x-2)}$ . (b) The 28 by 28 image that
corresponds to the curve in (a). (c) The curve of $e^{-(x-2)}+\sigma z$ when
$\sigma=.1$ where $z$ is standard normal random variable. (d) The 28 by 28
image corresponds to the curve in (c). (e) The curve of $e^{-(x-2)}+\sigma z$
when $\sigma=1$ where $z$ is standard normal random variable. (f) The 28 by 28
image corresponds to the curve in (e).

.

####  3.2.2 Convex versus Concave Curves

We also examined the ability of CNN to classify the curvature of the curves as
convex or concave. Curves $y=w_{1}(x-w_{2})^{2}$ are used to generate convex
and concave curves for training, validation, and training datasets; see Figure
11. We also use the following random variables $w_{1}=sign(U_{1}-.5)$ and
$w_{2}=2U_{2}+2.5$ where $U_{1}$ and $U_{2}$ are uniformly distributed random
variables on [0,1]. We found that the accuracy of the classification of the
functional test data is 100%.

Figure 11: (a)The curve of $y=(x-2)^{2}$ . (b) The 28 by 28 image that
corresponds to the curve in (a). (c) The curve of $(x-2)^{2}+\sigma z$ when
$\sigma=.1$ where $z$ is a standard normal random variable. (d) The 28 by 28
image corresponds to the curve in (c). (e) The curve of $(x-2)^{2}+\sigma z$
when $\sigma=1$ where $z$ is a standard normal random variable. (f) The 28 by
28 image corresponds to the curve in (e).

####  3.2.3 Exponential versus Algebraic Growth Curves

We examined the capabilities of CNN in classifying curve growth as exponential
(in the form of $y=e^{cx}$ for $c>0$) or algebraic (in the form of $y=x^{c}$
for $c>0$), see Figure 12. Training, validation, and training curves are
generated using a random $c=3U+1$ where $U$ is a uniformly distributed random
variable over [0,1]. We found that the classification accuracy of the
functional test data is 100%.

Figure 12: (a)The curve of an algebraic growth represented by $y=x^{3}$. (b)
The 28 by 28 image that corresponds to the curve in (a). (c) The curve of
$y=x^{3}+\sigma z$ when $\sigma=.1$ where $z$ is a standard normal random
variable. (d) The 28 by 28 image corresponds to the curve in (c). (e) The
curve of $y=x^{3}+\sigma z$ when $\sigma=1$ where $z$ is a standard normal
random variable. (f) The 28 by 28 image corresponds to the curve in (e).

####  3.2.4 Finding the Number of Peaks of Curves

The same CNN architecture that is used to estimate the magnitude of the
maximum height and the width of the peaks is used for classification with the
number of peaks as mentioned in Subsection 3.1.3. The curves are generated
with the same scheme as mentioned in Subsection 3.1.3. In this case, we use a
local normalization; that is, the highest grayscale value of that particular
curve is 1, and the lowest is 0. This helps to identify the peaks; since we
are not interested in the height, the grayscale values do not need to be
proportional to the height. A zero number of peaks is possible if the randomly
generated position is outside the interval [0,50]. We found that the
classification accuracy of the functional test data is 98. 0% for noise-free
data, 97. 2% for noisy data with magnitude $\sigma=.1$ and 96. 2% for noisy
data with magnitude $\sigma=1$.

##  4 Applications to Dynamical Systems

In this section, we apply a CNN to dynamical systems. First, we show how CNN
can estimate the Lyapunov exponent from one curve of the three variables
solving the Lorenz system. Second, we show that CNN can be used to estimate
transmission rates from epidemic curves. Third, we use CNN to test the
similarity and dissimilarity of drug dissolution profiles.

###  4.1 Estimating Lyapunov Exponent

The study of human motion has been the focus of interest in the medical field
to determine which exercise and range of motion are stable. Different
perspectives have been performed to study such stability. Although (?) pointed
out the link between movement variability and stability. (?) emphasized that
variability should not be confused with instability, as it can be observed in
healthy and unhealthy subjects. (?) helped develop the standard procedure for
analyzing stability using Lyapunov exponents estimated by the Rosenstein
method (?). Lyapunov exponents are the average exponential divergence of the
orbits between nearest trajectories, which are called nearest neighbors. The
more unstable the system, the higher the value of the Lyapunov exponent. An
alternative parameter is the mean of the absolute value of the Floquet
multipliers. (?) is the first to extensively use Floquet multipliers for the
same purpose. Floquet multipliers are the eigenvalues of the Jacobian matrix
that measure the separation between orbits of the system (?) that could be
calculated using Poincare maps. If the mean of the magnitudes of the
eigenvalues is less than 1, the orbits are considered stable.

To examine our method in estimating the Lyapunov exponent we use the prototype
of chaotic systems, Lorenz systems with parameters $\alpha,\beta,\rho$ defined
as the system of ordinary differential equations:

$\displaystyle\frac{dx}{dt}$ $\displaystyle=$ $\displaystyle\alpha(y-x)$
$\displaystyle\frac{dy}{dt}$ $\displaystyle=$ $\displaystyle x(\rho-z)-y$
$\displaystyle\frac{dz}{dt}$ $\displaystyle=$ $\displaystyle xy-\beta z$

Initial values are $x(0)=1$, $y(0)=1$, and $z(0)=1$. Figure 13 shows an
example image of the attractor with $\alpha=2.8029$, $\beta=1.1114$, and
$\rho=11.9620$.

We randomly simulate the values $\alpha=10\,U_{1}$, $\beta=8/3\,U_{2}$, and
$\rho=20\,U_{3}$, where $U_{i}$ are independent uniformly distributed random
variables in $[0,1]$ for $i=1,2,3$. We ran the Lorenz system simulation over
the time interval $[0,1]$ using the Runge-Kutta hybrid order 4 and 5 numerical
method in MATLAB. We use component $x$ only for the training, validation, and
testing of the CNN, see Figure 14.

Figure 13:  An example of Lorenz attractor at $\alpha=2.8029$, $\beta=1.1114$,
and $\rho=11.9620$ without noise (a), with noise of magnitude $\sigma=.1$ (b),
and with noise of magnitude $\sigma=1$ (c).

Figure 14: (a) The curve of component $x$ of the Lorenz attractor. (b) The 28
by 28 image that corresponds to the curve in (a). (c) The curve of $x+\sigma
w$ when $\sigma=.1$ where $w$ is a standard normal random variable. (d) The 28
by 28 image corresponds to the curve in (c). (e) The curve of $x+\sigma w$
when $\sigma=1$ where $w$ is a standard normal random variable. (f) The 28 by
28 image corresponds to the curve in (e).

Figure 15 shows the predicted values versus the estimated value of the
Lyapunov exponents that closely followed the diagonal line without an
intercept and a slope of one. Table 6 shows a strong diagonal linear
relationship between the true and predicted values of the Lyapunov exponents.

Figure 15: Results of Lyapunov exponent estimation using the test dataset without noise (a), with noise of magnitude $\sigma=.1$ (b), and with noise of magnitude $\sigma=1$ (c). Table 6: Correlation Coefficient, Intercept, and Slope for Noise-Free Testing Data with P-values Case | Correlation Coefficient ($r$) | Intercept (p-value) | Slope (p-value)
---|---|---|---
Without Noise | 0.916 | -0.0146 (0.8550) | 0.9600 (0.3477)
With Noise ($\sigma=0.1$) | 0.949 | 0.0109 (0.8734) | 1.0908 (0.0152)
With Noise ($\sigma=1$) | 0.910 | -0.1140 (0.0900) | 0.9369 (0.1476)

Also, we tested the capability that a CNN model trained with noise-free data
in estimating the Lyapunov exponent in noisy data. Figure 16 and Table 7 show
relatively good results.

Figure 16: Results of Lyapunov exponent estimation when the CNN is trained only using the noise-free data for the test dataset without noise (a), with noise of magnitude $\sigma=.1$ (b), and with noise of magnitude $\sigma=1$ (c). Table 7: Correlation Coefficient, Intercept, and Slope for Lyapunov Testing Data with P-values when CNN is Trained with Noise-Free Data. Case | Correlation Coefficient ($r$) | Intercept (p-value) | Slope (p-value)
---|---|---|---
Without Noise | 0.944 | 0.0567 (0.3275) | 0.8856 (0.0004)
With Noise ($\sigma=0.1$) | 0.928 | -0.1326 (0.0733) | 0.9871 (0.7472)
With Noise ($\sigma=1$) | 0.867 | -0.3390 (0.0015) | 1.0815 (0.1970)

An endurance test was also performed using 10000 test curves to estimate the
Lyapunov exponent and found that the CNN is approximately 600 times faster
(2.7628 seconds) compared to the MATLAB Rosenstein method (1692.6 seconds).

###  4.2 Estimating the Transmission Rates from Epidemic Curves

Estimation of transmission rates and exponential growth curve rates (see
Subsection 3.1.1) are important for emerging epidemics (?) such as at the
beginning of COVID-19 (?). Some epidemics grow algebraically and not
exponentially (?, ?) and so it is also helpful to discern them through
classification; see Subsection 3.2.3. One of the main goals in epidemiology is
to estimate the basic reproduction number $R_{0}$ which is almost always
proportional to the transmission rate $\beta$. If $R_{0}<1$, then the disease
diminishes; otherwise, there is a chance that it will become endemic.

We use a susceptible-infected-recovered (SIR) to produce epidemic curves with
different transmission rates $\beta$ and estimate that parameter. The disease
dynamics of a susceptible-infected-recovered (SIR) compartmental model follows
the system of differential equations

$\displaystyle\frac{dS}{dt}$ $\displaystyle=$
$\displaystyle\mu-\beta\,S\,I-\mu\,S$ $\displaystyle\frac{dI}{dt}$
$\displaystyle=$ $\displaystyle\beta\,S\,I-(\mu+\gamma)\,I$
$\displaystyle\frac{dR}{dt}$ $\displaystyle=$ $\displaystyle\gamma\,I-\mu\,R$

where $S$, $I$, and $R$ are the proportion of susceptible, infected, and
recovered individuals in the population at time $t$, such that $S+I+R=1$.
Initial values are $S(0)=.99$, $I(0)=.01$ and $R(0)=0$. The parameter $\mu$ is
the per capita birth/death rate, $\beta$ is the transmission rate, and
$\gamma$ is the recovery rate. The basic reproduction number in the above SIR
model is given by $R_{0}=\beta/(\mu+\gamma)$.

We assume values $\mu=1/(365*50)$ days-1, $\gamma=1/28$ days-1 and $\beta$
randomly selected from a uniform distribution over $(.01,1)$ to reflect a
basic reproduction number in the range of $(.28,28)$. The SIR model
simulations are run for $50$ days using the Runge-Kutta hybrid order 4 and 5
numerical method in MATLAB. See Figure 17 for a simulated epidemic curve $I$.

Figure 17: (a) A simulation curve of $I$. (b) The 28 by 28 image that
corresponds to the curve in (a). (c) The curve of $I+\sigma z$ when
$\sigma=.1$ where $z$ is a standard normal random variable. (d) The 28 by 28
image corresponds to the curve in (c). (f) The curve of $I+\sigma z$ when
$\sigma=1$ where $z$ is a standard normal random variable. (f) The 28 by 28
image corresponds to the curve in (e).

Figure 18 and Table 8 show strong prediction of transmission rates.

Figure 18: Results of SIR regression using the test set. (a) No noise data.(b) noisy data with $\sigma=.1$. (c) noisy data with $\sigma=1$. Table 8: Correlation Coefficient, Intercept, and Slope for Transmission Rate Testing Data with P-values Case | Correlation Coefficient ($r$) | Intercept (p-value) | Slope (p-value)
---|---|---|---
Without Noise | 0.999 | -0.0023 (0.2573) | 1.0048 (0.1642)
With Noise ($\sigma=0.1$) | 0.997 | -0.0066 (0.1318) | 1.0177 (0.0229)
With Noise ($\sigma=1$) | 0.979 | 0.0005 (0.9686) | 0.9920 (0.7029)

###  4.3 Detecting Similarity of Drug Dissolution Profiles

The problem of drug release or dissolution profiles is important for the
pharmaceutical industry. Regulatory guidelines seek to advise coherent
characteristics of drug dissolution prior to their approval; see, e.g. (?, ?).
Many statistical approaches have been developed to test the similarity between
dissolution curves or profiles, including cluster analysis, decision trees,
and linear models; see, e.g. (?, ?, ?, ?, ?, ?). That also includes
nonparametric measures such as the two measures adopted by the US Food and
Drug Administration (FDA) and the European Medicines Agency (EMA),

$f_{1}=\frac{\sum_{i=1}^{n}|R_{i}-S_{i}|}{\sum_{i=1}^{n}R_{i}}\times 100\%$

and

$f_{2}=50\log_{10}\left(\left[1+\frac{1}{n}\sum_{i=1}^{n}(R_{i}-S_{i})^{2}\right]^{-0.5}\times
100\right)$

which are widely used to detect dissimilarity and similarity, respectively,
between two curves $\\{(t_{i},R_{i}):i=1,2,\ldots,n\\}$ and
$\\{(t_{i},S_{i}):i=1,2,\ldots,n\\}$. If $f_{1}$ is between 0 and 15 and
$f_{2}$ is between 50 and 100, then the two curves are considered similar;
see, for example, (?, ?) for a complete set of models and measures, as well as
FDA & EMA guidelines.

There are several mathematical models of drug dissolution; see, for example,
an important model of drug dissolution is the logistic curve
$f=\frac{100}{1+\exp(-c(t-6))}$ for some release rate $c>0$, see (?). We
follow (?) by using the logistic curve at the sampling time points, t = 5, 10,
15, 20, 30, 60, 90, and 120 minutes. To generate a number of similar and
dissimilar profiles, we use release rates $c_{1}=.01+.001*z$ and
$c_{2}=.03+.001*z$ where $z$ are generated randomly and independently from the
standard normal distribution. We can see the simulated curves in Figure 19 for
(a) dissimilar and (b) similar curves. In addition, we have the histogram for
$f_{1}$ and $f_{2}$ for the dissimilar (c) and similar (d) curves in the
training set. Figure also shows the histogram for $f_{1}$ and $f_{2}$ for the
dissimilar (e) and similar (f) curves in the test set.

Figure 19: (a) Example of dissimilar curves. (b) Examples of similar curves.
(c) Training data histogram of $f_{1}$ and $f_{2}$ for the dissimilar curves.
(d) Training data histogram of $f_{1}$ and $f_{2}$ for the similar curves. (e)
Test data histogram of $f_{1}$ and $f_{2}$ for the dissimilar curves. (f) Test
data histogram of $f_{1}$ and $f_{2}$ for the similar curves.

To detect the similarity of any two drug dissolution curves, we use a Siamese
CNN in which two parallel CNN’s final layers are inputs to a cross-entropy
measure of the two input images. Some minor changes are made to the overall
CNN architecture, in which we avoid using batch normalization and change the
average pooling layer to a max pooling layer. In the end, there is a dense
layer with $28^{2}$ hidden nodes. Also, the weights are initialized by
sampling from a normal distribution of mean zero and a standard deviation of
0.01. Following (?), we use the cross-entropy in the output layer to identify
the similarity between the images of the dissolution curves. The test resutls
can be seen in the confusion matrix in Figure 20.

Figure 20: Confusion matrix shows the true predicted positives and the true
predicted negatives without resulting in false positives or negatives when
comparing 50 similar and 50 dissimilar pairs of curves.

It is important to note that the number of hidden nodes in the last layer
helps in convergence. We notice that it is possible to use one hidden node;
however, convergence is not consistently guaranteed. The number of nodes of
the hidden layer therefore should be considered an important hyperparameter.
Furthermore, below $28^{2}$ nodes, there was no considerable change in
convergence time.

##  5 Real-Life Application: Detecting Parkinson’s Disease

Parkinson’s disease is a progressive neurodegenerative disorder that results
in motor and non-motor symptoms such as tremors, rigidity, and impaired
movement control. Detecting Parkinson’s disease involves a thorough physical
examination to assess motor skills, reflexes, muscle strength, and
coordination, and searching for characteristic signs of Parkinson’s disease.

Our technique could be successfully applied to detect Parkinson’s disease
using motor tests. We use a dataset introduced by (?, ?) in which 62
Parkinson’s patients and 15 healthy subjects draw a spiral curve on a tablet.
The original test was divided into three parts: a static test, a dynamic test,
and a circular motion test. In the static test, subjects draw a certain fixed
shape. In the dynamic test, the subjects draw a shape that disappears and
appears at certain times, so the subjects need to memorize their location and
continue drawing. In the circular motion test, subjects draw a circle at a red
reference point on the tablet. Here, we use the data set for the static
drawing test of Parkinson’s patients by (?), which contains at each time stamp
the positions $X$, $Y$, $Z$ of the drawing, the pressure that we call $P$, and
the grip angle that we call $A$, see Figure 21.

Figure 21: (a) A spiral drawing of a normal subject. (b) A spiral drawing of a
Parkinson’s patient.

(?) used three ML/DL techniques, random forest, logistic regression, and an
artificial neural network, to classify subjects into Parkinson’s and non-
Parkinson’s. While the original data set had only 77 subjects, (?) used tens
of thousands of instances. Among the three ML/DL techniques, the artificial
neural network showed 100% precision in all cases.

(?) used logistic regression, random forest, support vector machines, and
K-Nearest Neighbors using unbalanced data from 25 patients and 15 controls.
Logistic regression was the best model using the AUC. With the same idea, (?)
used logistic regression, a support vector machine, and a restricted Boltzmann
machine followed by a neural network for the same data set. The restricted
Boltzmann machine with the neural network model achieved an accuracy of 95%.
See also a review on the use of ML/DL techniques to detect Parkinson’s disease
(?).

Regarding the data set used, due to the small size of the data, we use
augmentation (see (?) for more details on data augmentation) to make a larger
training and testing data set. We made a tensor, in which the third dimension
represents the fixed signals of $X$ and $Y$, while we checked a combination of
$Z$, the pressure measurement $P$, and the grip angle $A$. So, we have 8
tensors with size 1024 by 1024 by 2 plus the combination of the last 3
variables to determine which has the best results. Using MATLAB data
augmentation (see (?) functions for more details), we made a random
translation of $X$ within $[-25,25]$, random reflections, and translation of
$Y$ within $[-50,50]$, adding a random translation for $Z$ within $[-.5,.5]$.
We also performed a random translation of the pressure within $[0,50]$ and a
random translation of the angle within $[0,25]$. Data augmentation was
performed four times more for the control class of the set since there are
only 15 control images compared to the 61 case images. By doing that
augmentation, we prepared training data of 300 control images and 305 disease
images. In combination, we split it into 80% for the training and validation
set, and the rest 20% is used in testing. See examples of Parkinson’s and
control patients in Figure 22.

Figure 22: Curves with their respective image transformation. (a) Example of a
control subject’s $x$ spiral component time series. (b) The 28 by 28 image
corresponds to the curve in (a). (c) Example of a Parkinson’s patient’s $x$
spiral component time series. (d) The 28 by 28 image corresponds to the curve
in (c).

Our method gave 100% validation and testing accuracy for all of the
combinations of the features $X$, $Y$, $Z$ coordinates, the pressure $P$, and
the griping angle of $A$. Using the simplest model with $X$ and $Y$ gives the
confusion matrix in Figure 23.

Figure 23: Confusion matrix shows the true predicted positives and the true
predicted negatives without resulting in false positives or negatives when
using $xy$ combination only.

##  6 Discussion

In this paper, we tested our new method of functional data learning using
convolutional neural networks (CNN) with various examples and applications.
CNN performance was very close to perfect, as evidenced by the test curves of
the functional cases in both regression and classification problems. We also
see that the variation starts to appear in practical cases such as the chaotic
system of Lorentz’s attractor and when estimating transmission rates from
epidemic curves of the SIR system. Also, we found that training the CNN with
noisy data improves CNN performance. These results show that in both cases the
new method is robust to noise and would handle different cases of functional
data. While some of the p-values of the results show slopes and intercepts
that are statistically significant from one to zero, respectively, their
estimates and the correlation coefficients show close values to one and zero,
indicating an excellent effect size. This conflict might be the result of
using large data sizes when testing the trained CNN.

On the practical side, a pre-trained CNN can be used in some of the
applications. A pre-trained CNN could estimate Lyapunov exponents and assess
the stability of some systems. For example, it could be used in the medical
field to determine the stability of human motion or walking gait. This
methodology provides a more practical approach in which, with moderate noise,
the CNN performs well and is approximately 2 orders of magnitude faster. As
such, the measured data can be used as input without the need to filter before
using the CNN since it is robust to noise. A CNN pre-trained on epidemic
curves can be used to estimate the transmission rate or directly estimate the
basic reproduction number $R_{0}$ and discern exponential growth from
algebraic growth.

When it comes to classification problems, the new method always gave an
accuracy of 100%. We tested the new method for the classification of curves
according to their monotonicity and curvature, as well as their type of
growth. Furthermore, we simulated drug dissolution profiles to train a Siamese
CNN which accurately managed to determine their similarity and dissimilarity.
Finally, using real-life data, CNN trained with functional motion data of a
few cases of Parkinson’s disease and even fewer controls discerned cases with
an accuracy of 100%.

##  7 Conclusion

In this paper, we show a simple method to convert any curve into an image.
Using convolution neural networks (CNN), we trained on a number of those
images, together with a validation set of images in various regression and
classification problems. The same technique could be used for regression and
classification problems in gene analysis and other medical sciences. Other
areas where one can explore are multiple output problems in which several
parameters are estimated, or classification and regression problems are
combined to find the type of curve and estimate its parameters. Extension of
the method to allow other types of kernel to produce those images might be a
viable extension to the main idea in this paper. However, the presented
technique might require large functional data that we could not find at the
time of writing this paper. In that case, we had to perform a data
augmentation. New methods for efficient functional synthetic data might be
needed to handle small functional data learning problems.

## Acknowledgment

This work was supported by the U.S. Department of Defense Manufacturing
Engineering Education Program (MEEP) under Award No. N00014–19-1–2728.

The authors also thank Salman Rahman and Harrinson Arrubla for their early
discussions and trials on this project. The authors did not receive any
funding for the research performed in this article. The authors declare that
they have no conflict of interest.

## Data Availability Statement

This study did not involve the creation, collection, or generation of new
data. The study mainly relied on publicly available data previously published.

## Conflict of interest

The authors declare that they have no conflict of interest.

## Code details

The codes used in this study can be accessed at
https://github.com/jesusgl86/FDAP01.

## References

## References

  * [1]
  * [2] [] Abend A M, Hoffelder T, Cohen M J, Van Alstine L, Diaz D A, Fredro-Kumbaradzi E, Reynolds J, Zheng Y, Witkowski K & Heimbach T 2023 AAPS Journal .
  * [3]
  * [4] [] Abraham A, Pedregosa F, Eickenberg M, Gervais P, Mueller A, Kossaifi J, Gramfort A, Thirion B & Varoquaux G 2014 Frontiers in neuroinformatics 8, 14.
  * [5]
  * [6] [] Akyol K 2017 Int J Inf Technol Comput Sci 9, 45–51.
  * [7]
  * [8] [] Basna R, Nassar H & Podgórski K 2021 arXiv preprint arXiv:2103.07453 .
  * [9]
  * [10] [] Boonpatcharanon S, Heffernan J M & Jankowski H 2022 PloS ONE 17(6), e0269306.
  * [11]
  * [12] [] Chowell G, Viboud C, Hyman J M & Simonsen L 2015 PLoS Currents .
  * [13]
  * [14] [] Costa P & Sousa Lobo J M 2001 European Journal of Pharmaceutical Sciences .
  * [15]
  * [16] [] Dingwell J B & Cusumano J P 2000 Chaos: An Interdisciplinary Journal of Nonlinear Science 10(4), 848–863.
  * [17]
  * [18] [] Dingwell J B & Kang H G 2006 Journal of Biomechanical Engineering 129(4), 586–593.
  * [19]
  * [20] [] Enǎchescu D 2010 Biocybernetics and Biomedical Engineering .
  * [21]
  * [22] [] Górecki T, Krzyśko M, Waszak Ł & Wołyński W 2018 Statistical Papers 59, 153–182.
  * [23]
  * [24] [] He K, Zhang X, Ren S & Sun J 2016 in ‘Proceedings of the IEEE conference on computer vision and pattern recognition’ pp. 770–778.
  * [25]
  * [26] [] Hernández-Roig H A, Aguilera-Morillo M C & Lillo R E 2020 in ‘Functional and High-Dimensional Statistics and Related Fields 5’ Springer pp. 115–122.
  * [27]
  * [28] [] Hernández-Roig H A, Aguilera-Morillo M C & Lillo R E 2021 Mathematics 9(4), 406.
  * [29]
  * [30] [] Howard A G, Zhu M, Chen B, Kalenichenko D, Wang W, Weyand T, Andreetto M & Adam H 2017 arXiv preprint arXiv:1704.04861 .
  * [31]
  * [32] [] Hurmuzlu Y & Basdogan C 1994 Journal of biomechanical engineering 116(1), 30–36.
  * [33]
  * [34] [] Isenkul M, Sakar B, Kursun O et al. 2014 in ‘The 2nd international conference on e-health and telemedicine (ICEHTM-2014)’ Vol. 5 pp. 171–175.
  * [35]
  * [36] [] Isenkul M, Sakar B & PhD. P 2017 ‘Parkinson Disease Spiral Drawings Using Digitized Graphics Tablet’ UCI Machine Learning Repository.  DOI: https://doi.org/10.24432/C5Q01S.
  * [37]
  * [38] [] James G, Witten D, Hastie T, Tibshirani R et al. 2013 An introduction to statistical learning Vol. 112 Springer.
  * [39]
  * [40] [] Kamble M, Shrivastava P & Jain M 2021 Measurement: Sensors 16, 100047.
  * [41]
  * [42] [] Koch G, Zemel R, Salakhutdinov R et al. 2015 in ‘ICML deep learning workshop’ Lille.
  * [43]
  * [44] [] Kolebaje O T, Vincent O R, Vincent U E & McClintock P V 2022 Communications in Nonlinear Science and Numerical Simulation .
  * [45]
  * [46] [] LeCun Y, Bottou L, Bengio Y & Haffner P 1998 Proceedings of the IEEE 86(11), 2278–2324.
  * [47]
  * [48] [] Maggio R M, Castellano P M & Kaufman T S 2008 European Journal of Pharmaceutical Sciences .
  * [49]
  * [50] [] MATLAB 2022 version 9.12.0 (R2022a) The MathWorks Inc. Natick, Massachusetts.
  * [51]
  * [52] [] Mei J, Desrosiers C & Frasnelli J 2021 Frontiers in aging neuroscience 13, 633752.
  * [53]
  * [54] [] Mughal M R, Koch H, Huang J, Chiaromonte F & DeGiorgio M 2020 PLoS genetics 16(8), e1008896.
  * [55]
  * [56] [] Paixão P, Gouveia L F, Silva N & Morais J A 2017 European Journal of Pharmaceutics and Biopharmaceutics .
  * [57]
  * [58] [] Pfisterer F, Beggel L, Sun X, Scheipl F & Bischl B 2019 arXiv preprint arXiv:1911.07511 .
  * [59]
  * [60] [] Polson N & Sokolov V 2020 Wiley Interdisciplinary Reviews: Computational Statistics 12(5), e1500.
  * [61]
  * [62] [] Pourmohamad T & Ng H K T 2023 Pharmaceutical Statistics .
  * [63]
  * [64] [] Rafajłowicz W & Rafajłowicz E 2021 in ‘Artificial Intelligence and Soft Computing: 20th International Conference, ICAISC 2021, Virtual Event, June 21–23, 2021, Proceedings, Part I 20’ Springer pp. 485–495.
  * [65]
  * [66] [] Ramsay J O 1991 Psychometrika 56(4), 611–630.
  * [67]
  * [68] [] Ramsay J O 1995 Psychometrika 60(3), 323–339.
  * [69]
  * [70] [] Ramsay J O 1996 Journal of the Royal Statistical Society Series B: Statistical Methodology 58(3), 495–508.
  * [71]
  * [72] [] Ramsay J O 2000 Journal of the American Statistical Association 95(449), 9–15.
  * [73]
  * [74] [] Ramsay J O 2002 Modeling intraindividual variability with repeated measures data: Methods and applications pp. 87–102.
  * [75]
  * [76] [] Ramsay J O & Silverman B W 2002 Applied functional data analysis: methods and case studies Springer.
  * [77]
  * [78] [] Ramsay J O & Silverman B W 2005a Functional Data Analysis Springer New York.
  * [79]
  * [80] [] Ramsay J O & Silverman B W 2005b Springer New York New York, NY pp. 1–18.
*#1
  * [81]
  * [82] [] Redmon J, Divvala S, Girshick R & Farhadi A 2016 in ‘Proceedings of the IEEE conference on computer vision and pattern recognition’ pp. 779–788.
  * [83]
  * [84] [] Redmon J & Farhadi A 2017 in ‘Proceedings of the IEEE conference on computer vision and pattern recognition’ pp. 7263–7271.
  * [85]
  * [86] [] Rosenstein M T, Collins J J & De Luca C J 1993 Physica D: Nonlinear Phenomena 65(1-2), 117–134.
  * [87]
  * [88] [] Sandler M, Howard A, Zhu M, Zhmoginov A & Chen L C 2018 in ‘Proceedings of the IEEE conference on computer vision and pattern recognition’ pp. 4510–4520.
  * [89]
  * [90] [] Simonyan K & Zisserman A 2014 arXiv preprint arXiv:1409.1556 .
  * [91]
  * [92] [] Stergiou N & Decker L M 2011 Human movement science 30(5), 869–888.
  * [93]
  * [94] [] Tan M & Le Q 2019 in ‘International conference on machine learning’ PMLR pp. 6105–6114.
  * [95]
  * [96] [] Thakur M, Dhanalakshmi S, Kuresan H, Senthil R, Narayanamoorthi R & Lai K W 2023 Journal of Ambient Intelligence and Humanized Computing 14(1), 175–189.
  * [97]
  * [98] [] Tuite A R & Fisman D N 2020 ‘Reporting, epidemic growth, and reproduction numbers for the 2019 novel coronavirus (2019-NCoV) epidemic’.
  * [99]
  * [100] [] Vranić E, Mehmedagić A & Hadžović S 2002 Bosnian Journal of Basic Medical Sciences .
  * [101]
  * [102] [] Wilson C, Simpson S E, Van Emmerik R E & Hamill J 2008 Sports biomechanics 7(1), 2–9.
  * [103]
  * [104] [] Xiaoying W, Qian S & Jialiang G 2021 in ‘2021 2nd International Conference on Education, Knowledge and Information Management (ICEKIM)’ IEEE pp. 562–569.
  * [105]
  * [106] [] Yao J, Mueller J & Wang J L 2021 in ‘International Conference on Machine Learning’ PMLR pp. 11898–11908.
  * [107]
  * [108] [] Yarger D, Stoev S & Hsing T 2022 The Annals of Applied Statistics 16(1), 216–246.
  * [109]
  * [110] [] Zhang C, Yan H, Lee S & Shi J 2021 Technometrics 63(3), 370–383.
  * [111]
  * [112] [] Zhao J 2012 Mathematical Methods in the Applied Sciences 35(17), 2111–2121.
  * [113]