# Poissonian Blurred Image Deconvolution by Framelet based Local Minimal Prior

Reza Parvaz111 Corresponding author<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Image production tools do not always create a clear image, noisy and blurry
images are sometimes created. Among these cases, Poissonian noise is one of
the most famous noises that appear in medical images and images taken in
astronomy. Blurred image with Poissonian noise obscures important details that
are of great importance in medicine or astronomy. Therefore, studying and
increasing the quality of images that are affected by this type of noise is
always considered by researchers. In this paper, in the first step, based on
framelet transform, a local minimal prior is introduced, and in the next step,
this tool together with fractional calculation is used for Poissonian blurred
image deconvolution. In the following, the model is generalized to the blind
case. To evaluate the performance of the presented model, several images such
as real images have been investigated.

Department of Mathematics, University of Mohaghegh Ardabili, 56199-11367
Ardabil, Iran.

Keywords: Deblurring; Framelet; Poissonian noise; Medical images; Fractional
calculation.

##  1 Introduction

Image is one of the most important tools in the development of various
sciences, such as medicine and astronomy. In these sciences, the details of
the image are very important in diagnosing the disease or analyzing an image.
Therefore, removing a detail from the image can cause unfavorable image
analysis. Also despite the development of imaging tools in astronomy and
medicine, it has not been possible to create a completely clear image. Noise
and blur are two elements that affect image quality. Deblurring problem is
influenced by the point spread function (PSF) and is divided into two types:
non-blind and blind problems. In the first type, we have information about the
PSF. And in the second case, only the blurred image is available. There are
different types of PSF in this discussion, such as average, disk, gaussian,
motion. These types of functions are defined in the image processing software
such as Matlab or OpenCV and can be easily used. Like the PSF, there are
different types of image noise that reduce image quality, such as the
Poissonian, salt and pepper and periodic noises. Poissonian noise also called
Photon noise has a Poissonian distribution and appears in medical images as
X-ray image [1] and astronomical images [2]. The reader can find more
information about the structure of this type of noise in [3]. The general
structure of Poissonian blurred image is written as follows

$\displaystyle Y=\mathbf{P}(K\circledast X),$ (1.1)

where $Y$ and $X\in\mathbb{R}^{n\times m}$ denote blurred and clear images,
respectively, $K\in\mathbb{R}^{l\times s}$ shows the point spread function and
$\circledast$ represents two dimensional convolution operator. $\mathbf{P}$
also represents the process of adding Poissonian noise to the blurred image.
This type of problem can be expressed as a system of linear equations as
follows

$\displaystyle y=\mathbf{P}(kx),$ (1.2)

where $y,x\in\mathbb{R}^{nm\times 1}$ and $k\in\mathbb{R}^{nm\times nm}$. This
system of equations can not be solved directly due to the large volume and
ill-conditioned. Also, the structure of the blurred kernel matrix depends on
the boundary conditions of the problem as zeros and periodic bounders. It is
well known that this matrix can be obtained by the fast Fourier transform(FFT)
under the periodic condition. Using the Fourier method to create matrices has
the advantage of reducing computational volume and program execution time. The
first attempt to solve this type of problem can be to use the singular value
decomposition, although this may not be an efficient method despite Poissonian
noise. One of the efficient methods to solve this type of problem is to use
the total variation (TV) method [4], which is presented by the first article
and is used in most articles. The results for the nonblind and blind types of
problems using this method show the efficiency of the TV method in obtaining
the appropriate solution for deblurring problem. With the expansion of the
discussion on fractional derivatives, the TV method is developed and various
articles are used based on fractional derivatives [5, 6]. Another tool used to
improve the TV method is the use of framelet transform [7]. This transform
provides a sparse approximation and has been used in various articles as [8,
9, 10]. In recent years, various methods have been used for Poissonian blurred
image deconvolution. The frame based algorithm by a nonstationary accelerating
alternating direction method is introduced in [11], an algorithm based on
regularized by a denoiser constraint and deep image prior is used in [12],
also the non-local total variation is introduced in [13]. One of the efficient
tools used to solve the blurred image deconvolution problem is the dark
channel metric [14, 15]. Although this method has a good solution, but the
volume of calculations of this method is high. To solve this problem, in [16],
an algorithm based on the local minimal pixels over non overlapping patches is
introduced. Due to the fact that in this method all pixels of the image are
not evaluated, the computational volume is considerably reduced. In the
proposed algorithm, this method is improved by using framelet theory.
According to the simulated results that are compared in the next section, the
difference between blurred and clear images is more visible compared to the
method in [15, 16]. And this sharp discrepancy caused by the use of framelet
provides a good tool for solving the blurred image deconvolution. In addition
to the framelet based local minimal prior, the fractional derivative based on
the discrete form of the G-L definition is used in the proposed algorithm.
Using this method becomes an efficient tool for approximating the edges of a
restored image. In the real world, we only know the information inside the
picture frame. Lack of data outside the frame causes problems in image
deblurring. One of these problems that appears in the restored image is the
ringing artifacts near the image borders. Different methods to solve this
problem are studied in different articles. In [17], no-boundary convolution by
choosing the valid option in Matlab’s convolution is proposed. Also, in [18] a
method based on Fourier domain restoration filter and extrapolated image is
proposed. In this article, the method that is introduced in [18] is used in
the proposed algorithm. One of the most important advantages of this method is
that the bounders of the restored image are not lost.

Outline: The organization of this paper is as follows: Motivation and
preliminaries of the paper are studied in Sect. 2. Also in this section,
framelet based local minimal prior which has an important role in the proposed
model is introduced. In Sect. 3, the proposed model for Poissonian blurred
image deconvolution is introduced and a numerical algorithm is proposed to
solve this model. The results of the numerical algorithm for the proposed
model is shown in Sect. 4. The conclusion of the paper is given in Sect. 5.
Notations: The notations used are summarized as: $F(\cdot)$ and
$F^{-1}(\cdot)$ denote the fast Fourier transform (FFT) and the inverse fast
Fourier transform (IFFT), $\circledast$ and $\langle\cdot,\cdot\rangle$ are
used for two dimensional convolution operator and inner product, respectively.
And separable Hilbert space is shown by $H$.

##  2 Motivation and preliminaries

In this section, we introduce the tools and ideas used and study them.

###  2.1 Framelet theory and fractional derivative

A countable set as $X$ in a separable Hilbert space $H$ is named a tight frame
if for all $f\in H$ the following relation holds

$\displaystyle f=\sum_{\psi\in X}\langle f,\psi\rangle\psi.$

The collection $X(\Psi)=\\{2^{j/2}\psi_{i}(2^{j}.-t);1\leq i\leq
n,~{}t,j\in\mathbb{Z}\\}$ is called wavelet system for
$\Psi=\\{\psi_{0},\cdots,\psi_{n}\\}$. If the wavelet system be a tight frame,
then this system is called a tight wavelet frame and each element of this
system is named a framelet. For a compactly supported scaling function like
$\phi$ a wavelet tight frame is obtained by

$\displaystyle
F\big{(}\phi_{i}(2\omega)\big{)}=F(h_{i})F\big{(}\phi(\omega)\big{)},~{}~{}i=0,\ldots,n.$

The piecewise linear B-spline framelets with low-pass and high-pass filters as
$h_{0}=\frac{1}{4}[1,2,1],~{}h_{1}=\frac{\sqrt{2}}{4}[1,0,-1]$ and
$h_{2}=\frac{1}{4}[-1,2,-1]$, respectively, is used in proposed model. It is
suggested that the reader find more information about this concept from [7,
19].
The concept of derivative can be considered as the most widely used tools in
the image processing. With the development of the concept of ordinary
derivative to fractional derivatives, the use of this tool has been used in
various fields of image processing such as image denoising and image
deblurring [20, 21]. Despite having a fixed definition for a ordinary
derivative, various definitions have been developed for the fractional
derivative as Riemann-Liouville, Caputo, Caputo-Fabrizio and Grünwald-Letnikov
(G-L) [22, 23]. Due to the discrete structure of G-L definition, this type of
definition is most commonly used in image processing. In the proposed
algorithm, this definition is used. For a real function
$u:\Omega\subseteq\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$, the fractional
gradient of order $\alpha\in\mathbb{R}$ is defined as

$\displaystyle\nabla^{\alpha}u:=\big{(}\nabla^{\alpha}_{h}u,\nabla^{\alpha}_{v}u\big{)},$

where

$\displaystyle\nabla^{\alpha}_{h}u(i,j)=\sum^{L-1}_{l=0}\frac{(-1)^{l}\Gamma(\alpha+1)}{\Gamma(l+1)\Gamma(\alpha-l+1)}u(i-l,j),$
$\displaystyle\nabla^{\alpha}_{v}u(i,j)=\sum^{L-1}_{l=0}\frac{(-1)^{l}\Gamma(\alpha+1)}{\Gamma(l+1)\Gamma(\alpha-l+1)}u(i,j-l).$

###  2.2 Framelet based local minimal prior

For an image $I\in\mathbb{R}^{m\times n\times c}$, the framelet patch-wise
minimal pixels (FPMP) is defined as a vector of local minimal pixels over
nonoverlapping patches on the framelet transform and obtained by using
following formula

$\displaystyle
P_{w}(I)(i)=\min_{(x,y)\in\Omega_{i}}\big{(}\min_{c\in{\\{r,g,b\\}}}WI(x,y,c)\big{)},$

where $\Omega_{i}$ means the $i$-th patch of image. In the above formula,
after the framelet transform, the collection of local minimal pixels are
calculated for each of the low-pass and high-pass filters. Therefore, if the
size of patch is considered as $r\times r$ then
$P_{w}(I)\in\mathbb{R}^{9\varpi}$ where
$\varpi=\lceil\frac{m}{r}\rceil.\lceil\frac{n}{r}\rceil$. Similarly, the
framelet dark channel (FDC) is defined as

$\displaystyle
D_{w}(I)(i,j)=\min_{(x,y)\in\Omega_{i,j}}\big{(}\min_{c\in{\\{r,g,b\\}}}WI(x,y,c)\big{)},$

where $\Omega_{i,j}$ is an image patch centered at $(i,j)$-th pixel. By using
this definition, we can see that $D_{w}(I)\in\mathbb{R}^{3m\times 3n}$. Using
the same process as [15], the following properties are obtained.

######  Property 2.1.

Let $P_{w}(B)$ and $P_{w}(I)$ show the FPMP of the blurred and clear images,
respectively. Then

$\displaystyle P_{w}(B)\geq P_{w}(I).$

######  Property 2.2.

Let $D_{w}(B)$ and $D_{w}(I)$ show the FDC of the blurred and clear images,
respectively. Then

$\displaystyle D_{w}(B)\geq D_{w}(I).$

In addition to the above properties, the histogram plots are drawn in Fig.1,
which shows that the number of zero elements in the FPMP and FDC for the
blurred image is greater than these values for the clear image. Also in this
figure, a comparison between the FPMP, FDC, PMP and DC is done. Based on the
results, it can be seen that the difference in the number of non-zero elements
between the clear image and the blurred image in the FPMP and FDC is larger
than the difference of these value in the PMP and DC methods. Therefore, using
the $l_{0}$-norm along with the FPMP and FDC can be more efficient method in
the deblurring problems; remember that the $l_{0}$-norm counts the nonzero
elements. Most deblurring algorithms are iterative algorithms and the FDC is
large in size then this increases the computational time of the algorithm. To
avoid this problem the FPMP is used in the proposed algorithm.

(a) clear image($I$)

Figure 1: Histogram plots for FPMP, FCD, PMP and DC of clear and blurred
images.

##  3 Poissonian blurred image deconvolution

One of the most common ways to find a clear image $x$ is to use the Poisson
assumption for problem (1.2). The Poisson distribution for blurred problems
can be written as

$\displaystyle\mathbf{P}(y|kx)=\prod^{N}_{i=1}\frac{(kx)_{i}^{y_{i}}e^{(kx)_{i}}}{y_{i}!}.$

Based on the negative log-likelihood, we have

$\displaystyle-\log\mathbf{P}(y|kx)=\langle
1,kx\rangle-\big{\langle}y,\log(kx)\big{\rangle}+\sum^{N}_{i=1}\log(y_{i}!)$

where $1$ represents a vector that all elements are equal to one. The above
problem is an ill-conditioned problem, therefore, special methods are used to
solve this type of problem. In a common method, two main methods,
regularization method and maximum a posteriori (MAP) are used to solve this
type of problem In many algorithms, the following variational problem is used
to estimate clear image $x$.

$\displaystyle\min_{x\geq
0}\mu\big{\langle}1,kx-y\log(kx)\big{\rangle}+\varphi(x),$

where $\varphi(x)$ and $\mu>0$ denote for regularizer and regularization
parameters, respectively.

###  3.1 Non-blind deconvolution problem

The following model is proposed to solve the Poissonian blurred image
deconvolution.

$\displaystyle\min_{x\geq
0}\mu\big{\langle}1,kx-y\log(kx)\big{\rangle}+\lambda\|P_{w}(x)\|_{0}+\|\nabla^{\alpha}x\|_{l}+\delta_{\mathbb{R}^{N}_{+}}(x),$
(3.1)

where $\delta_{\mathbb{R}^{N}_{+}}$ denotes the indicator function over
$\mathbb{R}^{N}_{+}$. In this model $l$ is considered as $0$ and $1$.
Therefore, two different norms for the proposed model are studied in the
proposed model. In the continuation of this paper the model based on $l_{0}$
norm is named $l_{0}$-FPMP and the model based on $l_{1}$ norm is called
$l_{1}$-FPMP. In order to approximate the answer, the alternating direction
method of multipliers (ADMM) is used. By using of $z$, $v$ and $m$ as
auxiliary variables, Eq. (3.1) is reformulated as

$\displaystyle\min_{x\geq 0}\mu\langle 1,v-y\log
v\rangle+\lambda\|P_{w}(x)\|_{0}+\|z\|_{l}+\delta_{\mathbb{R}^{N}_{+}}(m),$
$\displaystyle\text{s.t.}~{}~{}v=kx,~{}z=\nabla^{\alpha}x,~{}m=x.$

Then augmented Lagrangian is written as

$\displaystyle L(x,v,z,m;p_{1},p_{2},p_{3})$ $\displaystyle=\mu\langle
1,v-y\log
v\rangle+\lambda\|P_{w}(x)\|_{0}+\|z\|_{l}+\delta_{\mathbb{R}^{n}_{+}}(m)$
$\displaystyle+\frac{\gamma}{2}\|kx-v+\gamma^{-1}p_{1}\|^{2}_{2}+\frac{\eta}{2}\|x-m+\eta^{-1}p_{2}\|^{2}_{2}$
$\displaystyle+\frac{\beta}{2}\|\nabla x-z+\beta^{-1}p_{3}\|^{2}_{2},$ (3.2)

where $\gamma$, $\eta$, $\beta$ are positive penalty parameters and
$p_{1},p_{2}$, $p_{3}$ are the Lagrangian multipliers. The iterative scheme
based on ADMM is given as

$\displaystyle
v^{j+1}=\arg\min_{v>0}L(x^{j},v,z^{j},m^{j};p^{j}_{1},p^{j}_{2},p^{j}_{3})=\arg\min_{v}\mu\langle
1,v-y\log v\rangle+\frac{\gamma}{2}\|kx^{j}-v+\gamma^{-1}p^{j}_{1}\|^{2}_{2},$
(3.3) $\displaystyle x^{j+1}$
$\displaystyle=\arg\min_{x}L(x,v^{j+1},z^{j},m^{j};p^{j}_{1},p^{j}_{2},p^{j}_{3})=\arg\min_{x}\frac{\gamma}{2}\|kx-v^{j+1}+\gamma^{-1}p^{j}_{1}\|^{2}_{2}+\frac{\eta}{2}\|x-m^{j}$
$\displaystyle+\eta^{-1}p^{j}_{2}\|^{2}_{2}+\lambda\|P_{w}(x^{j})\|_{0}+\frac{\beta}{2}\|\nabla^{\alpha}x-z^{j}+\beta^{-1}p^{j}_{3}\|^{2}_{2},$
(3.4) $\displaystyle
z^{j+1}=\arg\min_{z}L(x^{j+1},v^{j+1},z,m^{j};p^{j}_{1},p^{j}_{2},p^{j}_{3})=\arg\min_{z}\|z\|_{l}+\frac{\beta}{2}\|\nabla^{\alpha}x^{j+1}-z+\beta^{-1}p^{j}_{3}\|^{2}_{2},$
(3.5) $\displaystyle
m^{j+1}=\arg\min_{m}L(x^{j+1},v^{j+1},z^{j+1},m;p^{j}_{1},p^{j}_{2},p^{j}_{3})=\arg\min_{m}\delta_{\mathbb{R}^{n}_{+}}(m)+\frac{\eta}{2}\|x^{j+1}-m+\eta^{-1}p^{j}_{2}\|^{2}_{2},$
(3.6) $\displaystyle p^{j+1}_{1}=p^{j}_{1}+\gamma(kx^{j+1}-v^{j+1}),$
$\displaystyle p^{j+1}_{2}=p^{j}_{2}+\eta(x^{j+1}-m^{j+1}),$ $\displaystyle
p^{j+1}_{3}=p^{j}_{3}+\beta(\nabla^{\alpha}x^{j+1}-z^{j+1}).$

We note that proposed algorithm is based on the ADMM and the convergence
analysis of this method is studied in [24, 25, 26]. Also, in the numerical
section, the convergence of the proposed algorithm is studied using numerical
results.

###  3.2 Solving subproblems

In the previous subsection, the subproblems that are required for image
deblurring are introduced. In this subsection, the methods of solving these
subproblems are explained. It should be noted that in the following, for
simplicity, we avoid using the power symbols $j$ and $j+1$. In a similar way
to the quadratic formula for the roots of the general quadratic equation for
subproblem (3.3), we get

$\displaystyle v=\frac{\gamma kx+p_{1}-\mu+\sqrt{(\mu-\gamma kx-
p_{1})^{2}+4\mu\gamma g}}{2\gamma}.$ (3.7)

Finding the answer for (3.1) is not as simple as (3.3) and we need to design a
special method. By introduced auxiliary variable $n$, (3.1) is written as

$\displaystyle(x,n)$
$\displaystyle=\arg\min_{x,n}\frac{\gamma}{2}\|kx-v+\gamma^{-1}p_{1}\|^{2}_{2}+\frac{\eta}{2}\|x-m+\eta^{-1}p_{2}\|^{2}_{2}+\lambda\|n\|_{0}$
$\displaystyle+\frac{\beta}{2}\|\nabla^{\alpha}x-z+\beta^{-1}p_{3}\|^{2}_{2}+\frac{\rho}{2}\|n-P_{w}(x)\|^{2}_{2},$
(3.8)

The new problem (3.2) is divided into two subproblems as

$\displaystyle
x=\arg\min_{x}\frac{\gamma}{2}\|kx-\xi_{1}\|^{2}_{2}+\frac{\eta}{2}\|x-\xi_{2}\|^{2}_{2}+\frac{\beta}{2}\|\nabla^{\alpha}x-\xi_{3}\|^{2}_{2}+\frac{\rho}{2}\|n-P_{w}(x)\|^{2}_{2},$
(3.9) $\displaystyle
n=\arg\min_{n}\lambda\|n\|_{0}+\frac{\rho}{2}\|n-P_{w}(x)\|^{2}_{2},$ (3.10)

where
$\xi_{1}:=v-\gamma^{-1}p_{1},~{}~{}\xi_{2}:=m-\eta^{-1}p_{2},~{}~{}\xi_{3}:=z-\beta^{-1}p_{3}$.
Now, we define $x^{w}_{P}:=P_{w}^{T}(P_{w}(x))=Wx\circ M$ and
$\widehat{x}^{w}_{P}:=(1-M)\circ Wx$. Where $M$ is the binary matrix and
defined as

$\displaystyle M_{i,j}=\begin{cases}1,&\text{if (i,j) be a minimal pixel in
patch,}\\\ 0,&\text{other wise}.\end{cases}$

It is easy to see that $x=W^{T}\big{(}x^{w}_{P}+\widehat{x}^{w}_{P}\big{)}$.
For simplicity in performing calculations, we define $x_{p}:=W^{T}x^{w}_{P}$
and $\widehat{x}_{p}:=W^{T}\widehat{x}^{w}_{P}$, therefore
$x=x_{p}+\widehat{x}_{p}$. Based on these definitions, subproblem (3.9) can be
written as the two following subproblems.

$\displaystyle
x_{p}=\arg\min_{x_{p}}\frac{\gamma}{2}\|kx_{p}-\xi_{4}\|^{2}_{2}+\frac{\eta}{2}\|x_{p}-\xi_{5}\|^{2}_{2}+\frac{\beta}{2}\|\nabla^{\alpha}x_{p}-\xi_{6}\|^{2}_{2}+\frac{\rho}{2}\|x_{p}-\xi_{7}\|^{2}_{2},$
(3.11)
$\displaystyle\widehat{x}_{p}=\arg\min_{\widehat{x}_{p}}\frac{\gamma}{2}\|k\hat{x}_{p}-\xi_{8}\|^{2}_{2}+\frac{\eta}{2}\|\widehat{x}_{p}-\xi_{9}\|^{2}_{2}+\frac{\beta}{2}\|\nabla^{\alpha}\widehat{x}_{p}-\xi_{10}\|^{2}_{2},$
(3.12)

where $\xi_{i},~{}i=4,\cdots,10$ are defined as

$\displaystyle\xi_{4}:=\xi_{1}-k\widehat{x}_{p},~{}~{}\xi_{5}:=\xi_{2}-\hat{x}_{p},~{}~{}\xi_{6}:=\xi_{3}-\nabla^{\alpha}\widehat{x}_{p},~{}~{}\xi_{7}:=W^{T}P^{T}_{w}(n),$
$\displaystyle\xi_{8}=\xi_{1}-kx_{p},~{}~{}\xi_{9}=\xi_{2}-x_{p},~{}~{}\xi_{10}=\xi_{3}-\nabla^{\alpha}x_{p}.$

By using the optimal condition for subproblems (3.11)-(3.12), we get

$\displaystyle\big{(}\gamma
k^{T}k+(\eta+\rho)I+\beta(\nabla^{\alpha})^{T}\nabla^{\alpha}\big{)}x_{p}=\gamma
k^{T}\xi_{4}+\eta\xi_{5}+\beta(\nabla^{\alpha})^{T}\xi_{6}+\rho\xi_{7},$
$\displaystyle\big{(}\gamma k^{T}k+\eta
I+\beta(\nabla^{\alpha})^{T}\nabla^{\alpha}\big{)}\widehat{x}_{p}=\gamma
k^{T}\xi_{8}+\eta\xi_{9}+\beta(\nabla^{\alpha})^{T}\xi_{10}.$

These subproblems can be easily solved by assuming the periodic boundary
conditions and using fast Fourier transform. Therefore, $x_{p}$ and
$\widehat{x}_{p}$ are obtained by

$\displaystyle x_{p}=F^{-1}\Big{(}\frac{\gamma F(k^{T}\xi_{4})+\eta
F(\xi_{5})+\beta F\big{(}(\nabla^{\alpha})^{T}\xi_{6}\big{)}+\rho
F(\xi_{7})}{\gamma F^{*}(k)F(k)+\eta+\rho+\beta
F^{*}(\nabla^{\alpha})F(\nabla^{\alpha})}\Big{)},$ (3.13)
$\displaystyle\widehat{x}_{p}=F^{-1}\Big{(}\frac{\gamma F(k^{T}\xi_{8})+\eta
F(\xi_{9})+\beta F\big{(}(\nabla^{\alpha})^{T}\xi_{10}\big{)}}{\gamma
F^{*}(k)F(k)+\eta+\beta F^{*}(\nabla^{\alpha})F(\nabla^{\alpha})}\Big{)}.$
(3.14)

As a final step in approximating $x$, by using the proximal mapping of
$l_{0}$-norm, the solution of the subproblem (3.10) is obtained as follows.

$\displaystyle
n=\begin{cases}P_{w}(x),&|P_{w}(x)|^{2}\geq\frac{2\lambda}{\rho},\\\
0,&\text{other wise}.\end{cases}$ (3.15)

In the following algorithm, the solution of the subproblem (3.5) for $l=1$ can
be find by using the shrinkage operator as

$\displaystyle
z=\text{shrink}\big{(}\nabla^{\alpha}x+\beta^{-1},\beta^{-1}\big{)},$ (3.16)

where $\text{shrink}(a,b)=\operatorname{sign}(a)\max\big{(}|a|-b,0\big{)}$.
Also for $l=0$, the solution of the subproblem (3.5) can written as

$\displaystyle
z=\begin{cases}\nabla^{\alpha}x,&|\nabla^{\alpha}x|^{2}\geq\frac{2}{\beta},\\\
0,&\text{other wise}.\end{cases}$ (3.17)

And finally, $m$ in the subproblem (3.6) is obtained by using the following
formula.

$\displaystyle m=\max\big{(}x+\eta^{-1}p_{2},0\big{)}.$ (3.18)

The proposed algorithm, which is discussed in details in the above, is
summarized as Algorithm 1.

Algorithm 1: Non-blind image restoration algorithm.

Input: Blurred image $y$, blur kernel $k$.

Choose: $\mu,\lambda,\gamma,\beta,\rho$.

Initialize: $x^{0}=y,n^{0}=W(y),z=D(y),$

$~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}m^{0}=z^{0}=p^{0}_{1}=p^{0}_{2}=p^{0}_{3}=b^{0}=0$.

$j=0$

repeat

obtain $v^{j+1}$ by (3.7).

obtain $x_{p}$ by (3.13).

obtain $\widehat{x}_{p}$ by (3.14).

consider $x^{j+1}=x_{p}+\widehat{x}_{p}$.

obtain $n$ by (3.15).

obtain $z^{j+1}$ by (3.16) for $l_{1}$-FPMP or (3.17) for $l_{0}$-FPMP.

obtain $m^{j+1}$ by (3.18).

update $p^{j+1}_{1},p^{j+1}_{2}$ and $p^{j+1}_{3}$ by

$~{}~{}~{}p^{j+1}_{1}=p^{j}_{1}+\gamma(kx^{j+1}-v^{j+1}),$

$~{}~{}~{}p^{j+1}_{2}=p^{j}_{2}+\eta(x^{j+1}-m^{j+1}),$

$~{}~{}~{}p^{j+1}_{3}=p^{j}_{3}+\beta(\nabla^{\alpha}x^{j+1}-z^{j+1})$.

$j=j+1$.

until Stop condition is met

Output: Deblurred image $x$.

###  3.3 Blind deconvolution problem

In the previous subsections, non-blind problem is studied. But in the real
world, sometimes we have no information about the PSF, and we are faced with
the non-blind type of problem. In this subsection, we deal with this type of
problem. The proposed model for this type is as follows:

$\displaystyle\min_{x\geq 0,k\in\omega}\mu\langle
1,kx-y\log(kx)\rangle+\lambda\|P_{w}(x)\|_{0}+\|\nabla^{\alpha}x\|_{1}+\varrho\|\nabla
k\|_{1}+\delta_{\mathbb{R}^{N}_{+}}(x),$

where $\omega:=\\{k\in\mathbb{R}^{l\times s}|k_{i,j}\geq
0,\sum_{i,j}k_{i,j}=1\\}$. This problem is divided into two sub-problems. In
the first problem, the original image is approximated, and in the second
problem, the blurred kernel is approximated. The first problem is studied in
the previous subsection and the second problem can be expressed as follows

$\displaystyle\min_{k\in\omega}\mu\langle
1,kx-y\log(kx)\rangle+\varrho\|\nabla k\|_{1},$

In $(j+1)$th iterative procedure, the solution of this problem is obtained as
(see [27])

$\displaystyle
k_{t}^{j+1}=\frac{k^{n}}{1-\frac{\varrho}{\mu}div\big{(}\frac{\nabla
k^{j}}{|\nabla
k^{j}|}\big{)}}\big{(}X^{j}\big{)}^{T}\Big{[}\frac{y}{X^{j}k^{j}}\Big{]},$
(3.19) $\displaystyle k^{j+1}=\frac{k_{t}^{j+1}}{\sum k_{t}^{j+1}}.$ (3.20)

It should be noted that as mentioned in the introduction section to solve
ringing artifacts near the image borders, the method based on Fourier domain
restoration filter and extrapolated image [18] is used in this algorithm. The
proposed blind deblurring algorithm is summarized in Algorithm 2.

Algorithm 2: Blind image restoration algorithm.

Input: Blurred image $y$.

Choose: $\mu,\lambda,\gamma,\beta,\rho,\varrho$.

Initialize: $x^{0}=y,n^{0}=W(y),z=D(y),$

$~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}m^{0}=z^{0}=p^{0}_{1}=p^{0}_{2}=p^{0}_{3}=b^{0}=0$,

Consider $k^{0}$ as average kernel.

prepare the blurred image using algorithm in [18].

$j=0$

repeat

update $x^{j+1}$ by Algorithm 1.

update $k^{j+1}$ by (3.19)-(3.20)

$j=j+1$.

until Stop condition is met

Output: Deblurred image and kernel.

##  4 Experiment results

In the previous section, a model for Poissonian blurred image deconvolution is
introduced and then a numerical algorithm is proposed to solve this model. In
this section, in order to show the efficiency of the algorithm, some examples
are studied. To simulate the algorithm, the system based on Windows 10-64bit
and Intel(R) Core(TM) i3-5005U<EMAIL_ADDRESS>and MATLAB 2014b are used. In
the simulation part, `fspecial` and `poissrnd` are used to generate PSF and
Poissonian noise, respectively. Also different values are used in the
following examples for the peak to check the noise intensity. In the numerical
results by using Frobenius norm, the stop condition is considered as
$Error=\|x^{j+1}-x_{j}\|_{F}/\|x^{j}\|_{F}\leq tol$, where $10^{-4}$ and
$10^{-3}$ are chosen as $tol$ for nonblind and blind problems, respectively.
In this section to check the algorithm, the peak signal to noise ratio (PSNR)
and the structural similarity (SSIM) values are calculated. The differences
and similarities between the original and restored image is measured by PSNR
an obtained as

$\displaystyle PSNR=10\log 10\frac{p^{2}}{MSE(u,U)},$

where $u$ and $U$ denote original and restored image, respectively, $p$ shows
the maximum pixel value and $MSE$ is mean squared error. Also $SSIM$ is
calculated by

$\displaystyle
SSIM=\frac{(\mu_{o}\mu_{r}+C_{1})(2\sigma_{or}+C_{2})}{(\mu^{2}_{o}+\mu^{2}_{r}+C_{1})(\sigma_{o}+\sigma_{r}+C_{2})},$

where $\mu_{i},\sigma_{i},~{}i=o,r$ are the means and variances of the clear
and deblurred images, respectively, and $\sigma_{or}$ shows covariance between
the clear and deblurred images. The following values are used in the
simulation section and the best output are selected after calculating. patch
size$\in\\{12,15,25,35\\}$, $\alpha\in\\{0.1i|i=1,\cdots,10\\}$,
$\lambda,\rho,\eta\in\\{0.25,0.5,1,2,4\\}$,
$\beta,\gamma\in\\{10e-i|i=1,2,3,4,5\\}$, $\mu\in\\{5,100,200,300\\}$ and
$\varrho/\mu\in\\{0.1,1,5,20\\}$.

###  4.1 Non-blind problem

For the first example, the image of Cameraman with peak $255$ is blurred by
Gaussian PSF with size $9$ and standard deviation $\sqrt{3}$. The results of
the proposed algorithm is shown and compared in Fig. 2. By using results in
this figure, it is observed that the use $l_{1}$ norm compared to $l_{0}$ norm
makes the restored image smoother. Also, for a more accurate comparison, the
value of the PSNR is given in Table 1 for two images, Cameraman and Galaxy. In
these cases, the previous PSF is considered. According to these results, as
seen in Fig. 2, the use of $l_{0}$ norm does not always create a suitable
answer, but the use of $l_{1}$ norm causes a more appropriate result with the
compared algorithms. In order to examine the astronomical image, the results
of the proposed algorithm for Galaxy in Fig. 3 are compared with other
algorithms. Based on the simulation results, we can see that the proposed
method compares favorably or even better against compared algorithms.

(a) clear image($I$)

Figure 2: Visual comparison of Cameraman with peak $255$.

(a) clear image($I$)

Figure 3: Visual comparison for Galaxy with peak $255$. Table 1: Compared results of PSNR for different images and peaks. Imgae | Peak |  RLTV [28] |  NPTool [29] |  BM3D [30] |  FOTV [17] |  $l_{0}$-FPMP |  $l_{1}$-FPMP
---|---|---|---|---|---|---|---
| 255 | 24.49 | 23.74 | 24.56 | 24.59 | 24.33 | 25.62
Cameraman | 127.5 | 23.87 | 23.90 | 24.20 | 24.14 | 23.78 | 25.72
| 51 | 22.55 | 23.21 | 23.58 | 23.40 | 21.24 | 24.32
| 25.5 | 20.96 | 22.49 | 23.09 | 22.75 | 19.56 | 23.67
| 255 | 25.57 | 24.64 | 25.04 | 25.34 | 25.73 | 26.65
Galaxy | 127.5 | 24.99 | 24.28 | 24.66 | 24.90 | 25.43 | 26.07
| 51 | 24.24 | 24.12 | 24.14 | 24.38 | 24.82 | 25.54
| 25.5 | 23.28 | 23.83 | 23.77 | 23.86 | 22.28 | 24.41

In the following, in order to investigate the algorithm, we consider four
images, including fluo cells, hill, moon and saturn images. These images are
shown in Fig. 4. For fluo cells and hill motion PSF with 20 pixels and an
angle of 45 degree in a counterclockwise direction is used. Also Gaussian PSF
with size $9$ and standard deviation $\sqrt{3}$ is used for moon and saturn
images. The results of PSNR and SSIM values are given in Table 2. Numerical
results are given in this table for different values of peak and compared with
other algorithms. In this table, as seen in the previous examples, the results
of using $l_{0}$ norm have less accuracy than the algorithm based on $l_{1}$
norm. The results for the restored images for different peaks for fluo cells
image are given in Fig 5. These results are compared with other algorithms.
The results show the efficiency of the proposed method. In order to
investigate the convergence of the proposed algorithm, the figures of Energy,
PSNR, Error and SSIM curves for fluo cells image with peak$=255$ are drawing
in Fig. 6. As it is clear from these results, after a number of repetitions,
the resulting value approaches a certain value.

(a) clear image($I$)

Figure 4: Test images.

(a) clear image($I$)

Figure 5: Visual comparison for fluo cells image with different peaks. Table 2: Compared results of PSNR and SSIM for different images and peaks. |  |  method in [31] |  FOTV [17] |  $l_{0}$-FPMP |  $l_{1}$-FPMP
---|---|---|---|---|---
Image | Peak | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM
| 255 | 24.92 | 0.69 | 26.84 | 0.74 | 26.16 | 0.61 | 27.28 | 0.75
fluo cells | 127.5 | 24.11 | 0.81 | 26.98 | 0.84 | 22.73 | 0.59 | 27.26 | 0.84
| 51 | 22.48 | 0.92 | 26.41 | 0.93 | 21.98 | 0.78 | 26.78 | 0.94
| 25.5 | 16.55 | 0.82 | 18.04 | 0.82 | 18.14 | 0.59 | 18.18 | 0.83
| 255 | 22.38 | 0.45 | 22.67 | 0.45 | 22.14 | 0.45 | 23.21 | 0.47
hill | 127.5 | 22.04 | 0.63 | 23.03 | 0.64 | 22.10 | 0.48 | 23.40 | 0.65
| 51 | 21.40 | 0.85 | 22.88 | 0.86 | 20.72 | 0.69 | 22.94 | 0.85
| 25.5 | 21.11 | 0.94 | 22.81 | 0.94 | 17.51 | 0.49 | 22.99 | 0.94
| 255 | 25.02 | 0.67 | 25.65 | 0.70 | 26.68 | 0.72 | 27.70 | 0.73
moon | 127.5 | 24.33 | 0.79 | 25.63 | 0.82 | 25.08 | 0.69 | 26.19 | 0.82
| 51 | 23.55 | 0.91 | 25.57 | 0.93 | 23.16 | 0.83 | 25.97 | 0.93
| 25.5 | 23.36 | 0.96 | 25.16 | 0.96 | 20.26 | 0.63 | 25.47 | 0.97
| 255 | 32.26 | 0.89 | 31.24 | 0.90 | 31.83 | 0.90 | 33.34 | 0.91
saturn | 127.5 | 31.12 | 0.93 | 30.94 | 0.93 | 31.37 | 0.93 | 31.80 | 0.94
| 51 | 29.64 | 0.96 | 30.44 | 0.97 | 30.11 | 0.96 | 31.14 | 0.97
| 25.5 | 28.42 | 0.97 | 29.53 | 0.98 | 27.89 | 0.96 | 30.27 | 0.98

(a) clear image($I$)

Figure 6: Energy, PSNR, Error and SSIM curves for fluo cells image with
peak$=255$, (a-d) $l_{0}$-FPMP, (e-h) $l_{1}$-FPMP.

###  4.2 Blind problem

Now the proposed algorithm for blind case is studied in this part. As first
example, the spin image is blurred by PSF $k_{1}$ and $k_{2}$ where $k_{1}$ is
consider as Gaussian PSF with size $9$ and standard deviation $20$ and $k_{2}$
is obtained as motion PSF with 11 pixels and an angle of 50 degree. The
results of proposed algorithm for PSNR value are compared in Table 3 for
different kernels and peaks. Also, the results related to restored images are
given in Fig.8 and compared These results indicate the efficiency of the
proposed method in image restoring. In the last example of this section, we
study the real image using the proposed algorithm. For this example, different
values is chosen for the kernel size as $9\times 9$ (Fig.7-(b,c)), $12\times
12$ (Fig.7-(d)) and $15\times 15$ (Fig.7-(e)). The results are shown in Fig.
7. Comparing the results with the original image, it can be seen that the
restored images obtained from the proposed algorithm reveal more details about
the surface of the moon that is hidden in the original image. In these
figures, Fig.7-(b,g) show the restored image without using algorithm in [18].
Comparing this result with the images in which the algorithm [18] is used, it
is observed that the ringing artifacts is created near the image borders.

Table 3: Compared results of PSNR for satellite image. Imgae | Peak |  |  EM [17] |  |  FOTV [17] |  | proposed
---|---|---|---|---|---|---|---
$k_{1}$ | 255 |  | 14.92 |  | 21.44 |  | 22.18
| 127.5 |  | 14.72 |  | 21.26 |  | 22.56
$k_{2}$ | 255 |  | 16.80 |  | 23.44 |  | 24.56
| 127.5 |  | 16.63 |  | 20.98 |  | 24.36

(a) clear image($I$)

Figure 7: Visual results for real image.

(a) clear image($I$)

Figure 8: Visual comparison for satellite images.

##  5 Conclusion

In this paper the nonblind and blind Poissonian blurred image deconvolution
algorithm based on the local minimal prior are introduced. In the local
minimal prior, Framelet transform is used to find the difference between the
clear and the blurred images. Also, the alternating direction method of
multipliers (ADMM) is used to solve the proposed model. In the blind
deconvolution algorithm, Fourier domain restoration filter is used to remove
the ringing artifacts in the restored image. The results of simulation
algorithms have been compared with different methods and these results confirm
the efficiency of the proposed algorithms.

## References

  * [1] de Bruijne, Marleen, Philippe C. Cattin, Stéphane Cotin, Nicolas Padoy, Stefanie Speidel, Yefeng Zheng, and Caroline Essert, eds. Medical Image Computing and Computer Assisted Intervention-MICCAI 2021: 24th International Conference, Strasbourg, France, September 27-October 1, 2021, Proceedings, Part III. Vol. 12903. Springer Nature, 2021.
  * [2] Gendler, Robert. "Lessons from the Masters." Lessons from the Masters: Current Concepts in Astronomical Image Processing, The Patrick Moore Practical Astronomy Series, Volume 179. ISBN 978-1-4614-7833-1. Springer Science+ Business Media New York, 2013 1 (2013).
  * [3] Sikula, Josef, Michael Levinshtein, and Mikhail Efimovich Levinshteĭn, eds. Advanced experimental methods for noise research in nanoscale electronic devices. Vol. 151. Springer Science & Business Media, 2004.
  * [4] Rudin, Leonid I., Stanley Osher, and Emad Fatemi. "Nonlinear total variation based noise removal algorithms." Physica D: nonlinear phenomena 60, no. 1-4 (1992): 259-268.
  * [5] Chen, Gao, Gang Li, Yu Liu, Xiao-Ping Zhang, and Li Zhang. "SAR image despeckling based on combination of fractional-order total variation and nonlocal low rank regularization." IEEE Transactions on Geoscience and Remote Sensing 58, no. 3 (2019): 2056-2070.
  * [6] Ren, Zemin, Chuanjiang He, and Qifeng Zhang. "Fractional order total variation regularization for image super-resolution." Signal Processing 93, no. 9 (2013): 2408-2421.
  * [7] Han, Bin. "Framelets and wavelets." Algorithms, Analysis, and Applications, Applied and Numerical Harmonic Analysis. Birkhäuser xxxiii Cham (2017).
  * [8] Parvaz, Reza. "Point spread function estimation for blind image deblurring problems based on framelet transform." The Visual Computer (2022): 1-17.
  * [9] Liu, Jingjing, Yifei Lou, Guoxi Ni, and Tieyong Zeng. "An image sharpening operator combined with framelet for image deblurring." Inverse Problems 36, no. 4 (2020): 045015.
  * [10] Chen, Yingpin, Yuming Huang, Lingzhi Wang, Huiying Huang, Jianhua Song, Chaoqun Yu, and Yanping Xu. "Salt and pepper noise removal method based on stationary Framelet transform with non-convex sparsity regularization." IET Image Processing 16, no. 7 (2022): 1846-1865.
  * [11] Huang, Jie, and Ting-Zhu Huang. "A nonstationary accelerating alternating direction method for frame-based Poissonian image deblurring." Journal of Computational and Applied Mathematics 352 (2019): 181-193.
  * [12] Feng, Yayuan, Yu Shi, and Dianjun Sun. "Blind poissonian image deblurring regularized by a denoiser constraint and deep image prior." Mathematical Problems in Engineering 2020 (2020).
  * [13] Shi, Yu, Jie Song, and Xia Hua. "Poissonian image deblurring method by non-local total variation and framelet regularization constraint." Computers & Electrical Engineering 62 (2017): 319-329.
  * [14] He, Kaiming, Jian Sun, and Xiaoou Tang. "Single image haze removal using dark channel prior." IEEE transactions on pattern analysis and machine intelligence 33, no. 12 (2010): 2341-2353.
  * [15] Pan, Jinshan, Deqing Sun, Hanspeter Pfister, and Ming-Hsuan Yang. "Blind image deblurring using dark channel prior." In Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 1628-1636. 2016.
  * [16] Wen, Fei, Rendong Ying, Yipeng Liu, Peilin Liu, and Trieu-Kien Truong. "A simple local minimal intensity prior and an improved algorithm for blind image deblurring." IEEE Transactions on Circuits and Systems for Video Technology 31, no. 8 (2020): 2923-2937.
  * [17] Chowdhury, Mujibur Rahman, Jing Qin, and Yifei Lou. "Non-blind and blind deconvolution under poisson noise using fractional-order total variation." Journal of Mathematical Imaging and Vision 62, no. 9 (2020): 1238-1255.
  * [18] Liu, Renting, and Jiaya Jia. "Reducing boundary artifacts in image deconvolution." In 2008 15th IEEE International Conference on Image Processing, pp. 505-508. IEEE, 2008.
  * [19] Averbuch, Amir Z., Pekka Neittaanmäki, and Valery A. Zheludev. Spline and spline wavelet methods with applications to signal and image processing. Vol. 1. New York: Springer, 2014.
  * [20] Shukla, Anil K., Rajesh K. Pandey, and P. K. Reddy. "Generalized fractional derivative based adaptive algorithm for image denoising." Multimedia Tools and Applications 79, no. 19 (2020): 14201-14224.
  * [21] Guo, Zhichang, Wenjuan Yao, Jiebao Sun, and Boying Wu. "Nonlinear fractional diffusion model for deblurring images with textures." Inverse Problems & Imaging 13, no. 6 (2019): 1161.
  * [22] Yang, Xiao-Jun. General fractional derivatives: theory, methods and applications. Chapman and Hall/CRC, 2019.
  * [23] Cafagna, Donato. "Fractional calculus: A mathematical tool from the past for present engineers [Past and present]." IEEE Industrial Electronics Magazine 1, no. 2 (2007): 35-40.
  * [24] Gabay, Daniel, and Bertrand Mercier. "A dual algorithm for the solution of nonlinear variational problems via finite element approximation." Computers & mathematics with applications 2, no. 1 (1976): 17-40.
  * [25] Glowinski, Roland, and Americo Marroco. "Sur l’approximation, par éléments finis d’ordre un, et la résolution, par pénalisation-dualité d’une classe de problèmes de Dirichlet non linéaires." ESAIM: Mathematical Modelling and Numerical Analysis-Modélisation Mathématique et Analyse Numérique 9, no. R2 (1975): 41-76.
  * [26] Eckstein, Jonathan, and Dimitri P. Bertsekas. "On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators." Mathematical Programming 55, no. 1 (1992): 293-318.
  * [27] Dey, Nicolas, Laure Blanc-Feraud, Christophe Zimmer, Pascal Roux, Zvi Kam, Jean-Christophe Olivo-Marin, and Josiane Zerubia. "Richardson-Lucy algorithm with total variation regularization for 3D confocal microscope deconvolution." Microscopy research and technique 69, no. 4 (2006): 260-266.
  * [28] Dupé, F-X., Mohamed-Jalal Fadili, and J-L. Starck. "Image deconvolution under Poisson noise using sparse representations and proximal thresholding iteration." In 2008 IEEE International Conference on Acoustics, Speech and Signal Processing, pp. 761-764. IEEE, 2008.
  * [29] Landi, Germana, and E. Loli Piccolomini. "An efficient method for nonnegatively constrained Total Variation-based denoising of medical images corrupted by Poisson noise." Computerized Medical Imaging and Graphics 36, no. 1 (2012): 38-46.
  * [30] Azzari, Lucio, and Alessandro Foi. "Variance stabilization in Poisson image deblurring." In 2017 IEEE 14th International Symposium on Biomedical Imaging (ISBI 2017), pp. 728-731. IEEE, 2017.
  * [31] Jon, Kyongson, Jun Liu, Xiaoguang Lv, and Wensheng Zhu. "Poisson noisy image restoration via overlapping group sparse and nonconvex second-order total variation priors." Plos one 16, no. 4 (2021): e0250260.


\frac{\sigma}{{\overline{e}}(t)}\partial_{x}w(t,1)+{\overline{v}}w(t,1)&=0,&\text{for}\leavevmode\nobreak\
t\in\mathbb{R}_{+}^{*},\\\
\frac{\sigma}{{\overline{e}}(t)}\partial_{x}w(t,0)&=0,&\text{for}\leavevmode\nobreak\
t\in\mathbb{R}_{+}^{*},\\\
w(0,x)&=w^{0}(x):=g^{0}(x{\overline{e}}_{0}),&\text{for}\leavevmode\nobreak\
x\in(0,1),\end{aligned}\right.$ (B.1)

and the associated notion of weak $L^{2}$ solution:

######  Definition 7.

A function $w\in\mathcal{C}^{0}([0,+\infty),L^{2}(0,1))$ is said to be a
$L^{2}$-weak solution of (B.1) if for any $T>0$, it satisfies:

$\int_{0}^{T}\int_{0}^{1}w(t,x)\left[\partial_{t}\tilde{v}(t,x)+\frac{\sigma}{{\overline{e}}(t)^{2}}\partial_{xx}^{2}\tilde{v}(t,x)-\frac{{\overline{v}}}{{\overline{e}}(t)}x\partial_{x}\tilde{v}(t,x)-\left(\lambda+\frac{{\overline{v}}}{{\overline{e}}(t)}\right)\tilde{v}(t,x)\right]dxdt+\int_{0}^{1}w(0,x)\tilde{v}(0,x)dx=0,$

for any test function $\tilde{v}$ that satisfies:

  * •

$\tilde{v}\in\left(L^{2}\left((0,T);H^{2}(0,1)\right)\right)\cap\mathcal{C}^{0}([0,T],L^{2}(0,1°))$,

  * •

$\partial_{t}\tilde{v}\in\left(L^{2}\left((0,T);L^{2}\right)\right)$,

  * •

$\tilde{v}(T,\cdot)=0$,

  * •

$\sigma\partial_{x}\tilde{v}(t,0)=0,\leavevmode\nobreak\ \forall t\in(0,T)$,

  * •

$\sigma\partial_{x}\tilde{v}(t,{\overline{e}}(t))=0,\forall t\in(0,T)$.

Note that the two definitions 4 and 7 are equivalent: from the latter to the
former, take a test function of the form
${\overline{e}}(t)v(t,{\overline{e}}(t)x)$ where $v$ satisfies the assumptions
in 4. The other way around, take a test function of the form
$\frac{1}{{\overline{e}}(t)}\tilde{v}(t,\frac{x}{{\overline{e}}(t)})$ where
$\tilde{v}$ satisfies the previous assumptions.

The problem is uniformly parabolic: for any $0\leq t\leq T$

$\frac{\sigma}{{\overline{e}}(t)^{2}}\geq\frac{\sigma}{{\overline{e}}(T)^{2}}>0,$

so that we expect the classical parabolic estimates and well-posedness in
$\mathcal{C}^{0}([0,T],H^{k}(0,1))$ for any $k\in\mathbb{N}$ as soon as
$w^{0}\in H^{k}(0,1)$. In fact, we have the following a priori estimates

######  Lemma 9.

Assume $w^{0}\in L^{2}((0,1))$. Any smooth solution $w$ to (B.1) must satisfy
the energy estimate:

$\begin{split}&\frac{1}{2}\|w\|_{L^{\infty}(0,T;L^{2})}^{2}+\frac{\sigma}{{\overline{e}}(T)^{2}}\|\partial_{x}w\|_{L^{2}(0,T;L^{2}(0,1))}^{2}\\\
+&\frac{{\overline{v}}}{2{\overline{e}}(T)}\int_{0}^{T}w(t,1)^{2}+\left(\frac{{\overline{v}}}{2{\overline{e}}(T)}+\lambda\right)\|w\|_{L^{2}(0,T;L^{2}(0,1))}^{2}\leq\|w^{0}\|_{L^{2}(0,1)}^{2}.\end{split}$
(B.2)

Furthermore, it must satisfy the stability estimate: for any $0\leq t\leq T$

$\|w(t)\|_{L^{2}(0,1)}\leq e^{-\lambda t}\|w^{0}\|_{L^{2}(0,1)}$ (B.3)

###### Proof.

Multiply the first equation in (B.1) by $w$ and integrate by parts in space at
time $t$. One first obtains:

$\frac{1}{2}\frac{d}{dt}\|w(t)\|_{L^{2}}^{2}+\frac{\sigma}{{\overline{e}}(t)^{2}}\int_{0}^{1}(\partial_{x}w)^{2}+\left(\frac{{\overline{v}}}{2{\overline{e}}(t)}+\lambda\right)\int_{0}^{1}w^{2}-\frac{\sigma}{{\overline{e}}(t)^{2}}\partial_{x}w(t,1)w(t,1)-\frac{{\overline{v}}}{2{\overline{e}}(t)}w(t,1)^{2}=0.$
(B.4)

Then using the boundary conditions in (B.1), one gets:

$\begin{split}\frac{1}{2}\frac{d}{dt}\|w(t)\|_{L^{2}}^{2}+\frac{\sigma}{{\overline{e}}(t)^{2}}\int_{0}^{1}(\partial_{x}w)^{2}+\left(\frac{{\overline{v}}}{2{\overline{e}}(t)}+\lambda\right)\int_{0}^{1}w^{2}+\frac{{\overline{v}}}{2{\overline{e}}(t)}w(t,1)^{2}=0.\end{split}$
(B.5)

It comes in particular:

$\frac{1}{2}\frac{d}{dt}\|w(t)\|_{L^{2}(0,1)}^{2}\leq-\lambda\|w(t)\|_{L^{2}(0,1)}^{2},$

from which we conclude to (B.3) with the Gronwall lemma. Integrating (B.5)
with respect to time in $[0,T]$, one finds (LABEL:eq:energy-estimate).

∎

From these estimates, we define a notion of _energy solution_ for the problem.

######  Definition 8.

A function $w\in\mathcal{C}^{0}([0,T],L^{2}(0,1))\cap L^{2}((0,T);H^{1}(0,1))$
such that $\partial_{t}w\in L^{2}((0,T);(H^{1}(0,1))^{\prime})$ is an energy
solution to (B.1) if, for almost any time $0\leq t\leq T$ and any function
$\tilde{v}\in H^{1}(0,1)$, it satisfies

$\langle\partial_{t}w(t),\tilde{v}(t)\rangle_{(H^{1})^{\prime},H^{1}}+a(t;w,\tilde{v})=0,$
(B.6)

where the bilinear form $a$ is given by:

$a(t;w,\tilde{v})=\int_{0}^{1}\left(\frac{\sigma}{{\overline{e}}(t)^{2}}\partial_{x}w\partial_{x}\tilde{v}-\frac{{\overline{v}}}{{\overline{e}}(t)}w\left(\tilde{v}+x\partial_{x}\tilde{v}\right)+\lambda
w\tilde{v}\right)dx.$ (B.7)

By construction of the weak formulation, such solutions still satisfy the
previous estimates. In particular and by linearity, such a solution is unique,
if it exists. The existence follows from the Galerkin method (see [10][Theorem
10.9, p.341] for a general result). This is summarized in the following
proposition

######  Proposition 4.

Let $w^{0}\in L^{2}(0,1)$. There exists a unique energy solution to (B.1).
This solution satisfies (LABEL:eq:energy-estimate) and (B.3).

Proposition 4 gives existence to a weak $L^{2}$ solution since the energy
solution is a particular one, and this solution satisfies in particular (B.3).
But one cannot directly conclude that _any_ weak $L^{2}$ satisfies (B.3) and
deduce uniqueness from the estimate. Indeed, such an estimate cannot be
deduced directly from the weak formulation in (4). Instead, we will first
prove uniqueness from the weak formulation, then deduce that the only weak
$L^{2}$ solution satisfies indeed the estimate.

######  Lemma 10.

There is at most one weak $L^{2}$ solution in the sense of Definition 7.

###### Proof.

Consider two such solutions $w_{1},w_{2}\in\mathcal{C}^{0}([0,T];L^{2}(0,1))$.
Then the difference satisfies, for any test function $\tilde{v}$ that
satisfies the assumptions of Definition 7:

$\int_{0}^{T}\int_{0}^{1}(w_{1}-w_{2})\left(\partial_{t}\tilde{v}+\frac{\sigma}{{\overline{e}}(t)^{2}}\tilde{v}-\frac{{\overline{v}}}{{\overline{e}}(t)}x\partial_{x}\tilde{v}-\left(\lambda+\frac{{\overline{v}}}{{\overline{e}}(t)}\right)\tilde{v}\right)dxdt=0$

Now fix $S\in L^{2}((0,T);L^{2}(0,1))$ and consider the inhomogeneous dual
problem with source term $S$:

$\left\\{\begin{aligned}
\partial_{t}\tilde{v}+\frac{\sigma}{{\overline{e}}(t)^{2}}\partial^{2}_{xx}\tilde{v}-\frac{{\overline{v}}}{{\overline{e}}(t)}x\partial_{x}\tilde{v}-\left(\lambda+\frac{{\overline{v}}}{{\overline{e}}(t)}\right)\tilde{v}&=S,&\text{for}\leavevmode\nobreak\
(t,x)\in[0,T]\times(0,1),\\\
\partial_{x}\tilde{v}(t,1)&=0,&\text{for}\leavevmode\nobreak\ t\in[0,T],\\\
\partial_{x}\tilde{v}(t,0)&=0,&\text{for}\leavevmode\nobreak\ t\in[0,T],\\\
\tilde{v}(T,x)&=0,&\text{for}\leavevmode\nobreak\
x\in(0,1).\end{aligned}\right.$ (B.8)

Up to time reversal $t\to T-t$, this is a classical parabolic problem with
smooth coefficients. Therefore, there exists a (unique) solution
$\tilde{v}\in\left(L^{2}\left((0,T);H^{2}(0,1)\right)\right)\cap\mathcal{C}^{0}([0,T],L^{2}(0,1))$
and such that
$\partial_{t}\tilde{v}\in\left(L^{2}\left((0,T);L^{2}(0,1)\right)\right)$ (the
regularity is limited by $S\in L^{2}$). Consequently, $\tilde{v}$ can be taken
as a test function against $(w_{1}-w_{2})$ and it holds:

$\int_{0}^{T}\int_{0}^{1}(w_{1}-w_{2})Sdxdt=0.$

Since this is true for any $S\in L^{2}((0,T);L^{2}(0,1))$, $w_{1}=w_{2}$ and
uniqueness is proved.

∎

Proposition 2 follows from Proposition 4 and 10.

##  Appendix C Formal derivation of the backstepping kernel problems

As explained in the remarks in Section 4.4, the derivation is done assuming
that the kernels explicitly depend on time $t$, then we show they do not have
to depend on $t$. In the spirit of Section 4.3, we assume that all the
functions are smooth and we differentiate (4.2) at $x=0$. It gives:

$\partial_{x}g(t,0)=-k(t,0,0)\zeta(t,0),$

which suggests, since $\zeta(t,0)$ is undetermined, that the kernel $k$ should
be supplied with the condition

$k(t,0,0)=0,$ (C.1)

for the boundary condition at $x=0$ in (4.10) to be satisfied.

At this stage, it is unclear how the kernel depends on time. We make it
explicit by applying a rescaling in the space variable into a fixed domain.
More precisely, we consider the rescaled versions of problems (4.1) and
(4.10). The latter one was defined in (B.1) while the former is given by:

$\left\\{\begin{aligned}
\partial_{t}z-\frac{1}{{\overline{e}}(t)^{2}}\sigma\partial^{2}_{xx}z-\frac{{\overline{v}}}{{\overline{e}}(t)}x\partial_{x}z&=0,&\text{for}\leavevmode\nobreak\
(t,x)\in\mathbb{R}_{+}^{*}\times(0,1),\\\
\frac{\sigma}{{\overline{e}}(t)}\partial_{x}z(t,1)+{\overline{v}}z(t,1)&=\delta\psi(t),&\text{for}\leavevmode\nobreak\
t\in\mathbb{R}_{+}^{*},\\\
\frac{\sigma}{{\overline{e}}(t)}\partial_{x}z(t,0)&=0,&\text{for}\leavevmode\nobreak\
t\in\mathbb{R}_{+}^{*},\\\
z(0,x)&=\zeta^{0}(x{\overline{e}}(t)),&\text{for}\leavevmode\nobreak\
x\in(0,1).\end{aligned}\right.$ (C.2)

We consider the backstepping transformation associated to these rescaled
problems: for any $t\geq 0$, $x\in(0,1)$:

$w(t,x):=z(t,x)-\int_{0}^{x}{\tilde{k}}(t,x,y)z(t,y)dy,$ (C.3)

where $z$ is a solution to (C.2), $w$ a solution to (B.1) and ${\tilde{k}}$ an
unknown function defined in $D_{1}:=\\{0<y\leq x<1\\}$. Assume everything is
smooth and compute derivatives:

$\displaystyle\partial_{x}w(t,x)$
$\displaystyle=\partial_{x}z(t,x)-{\tilde{k}}(t,x,x)z(t,x)-\int_{0}^{x}\partial_{x}{\tilde{k}}(t,x,y)z(t,y)dy,$
(C.4) $\displaystyle\partial^{2}_{xx}w(t,x)$
$\displaystyle=\partial_{xx}^{2}z(t,x)-\partial_{x}\left({\tilde{k}}(t,x,x)z(t,x)\right)-\partial_{x}{\tilde{k}}(t,x,x)z(t,x)-\int_{0}^{x}\partial^{2}_{xx}{\tilde{k}}(t,x,y)z(t,y)dy.$
(C.5)

and:

$\displaystyle\partial_{t}w(t,x)$
$\displaystyle=\partial_{t}z-\int_{0}^{x}{\tilde{k}}(t,x,y)\partial_{t}z(t,y)dy-\int_{0}^{x}\partial_{t}{\tilde{k}}(t,x,y)z(t,y)dy$
$\displaystyle=\frac{\sigma}{{\overline{e}}(t)^{2}}\partial_{xx}^{2}z+\frac{{\overline{v}}}{{\overline{e}}(t)}x\partial_{x}z-\int_{0}^{x}{\tilde{k}}(t,x,y)\left(\frac{\sigma}{{\overline{e}}(t)^{2}}\partial_{yy}^{2}z+\frac{{\overline{v}}}{{\overline{e}}(t)}y\partial_{y}z\right)dy-\int_{0}^{x}\partial_{t}{\tilde{k}}(t,x,y)z(t,y)dy$

Now use integration by parts for the integral terms in the middle. It holds:

$\int_{0}^{x}{\tilde{k}}(t,x,y)y\partial_{x}z(t,y)dy={\tilde{k}}(t,x,x)xz(t,x)-\int_{0}^{x}\partial_{y}(y{\tilde{k}}(t,x,y))z(t,y)dy,$

and, using $\partial_{x}z(0)=0$:

$\displaystyle-\int_{0}^{x}{\tilde{k}}(t,x,y)\partial_{yy}^{2}z$
$\displaystyle=-{\tilde{k}}(t,x,x)\partial_{x}z(t,x)+\int_{0}^{x}\partial_{y}{\tilde{k}}(t,x,y)\partial_{x}zdy$
$\displaystyle=-{\tilde{k}}(t,x,x)\partial_{x}z(t,x)+\left(\partial_{y}{\tilde{k}}(t,x,x)z(t,x)-\partial_{y}{\tilde{k}}(t,x,0)z(t,0)\right)-\int_{0}^{x}\partial_{yy}^{2}{\tilde{k}}(t,x,y)z(t,y)dy,$

so that it holds:

$\begin{split}\partial_{t}w(t,x)=&\frac{\sigma}{{\overline{e}}(t)^{2}}\left(\partial_{xx}^{2}z-{\tilde{k}}(t,x,x)\partial_{x}z(t,x)+\left(\partial_{y}{\tilde{k}}(t,x,x)z(t,x)-\partial_{y}{\tilde{k}}(t,x,0)z(t,0)\right)\right.\\\
&-\left.\int_{0}^{x}\partial_{yy}^{2}{\tilde{k}}(t,x,y)z(t,y)dy\right)\\\
&+\frac{{\overline{v}}}{{\overline{e}}(t)}\left(x\partial_{x}z+\int_{0}^{x}\partial_{y}(y{\tilde{k}}(t,x,y))z(t,y)dy-{\tilde{k}}(t,x,x)xz(x)\right)-\int_{0}^{x}\partial_{t}{\tilde{k}}(t,x,y)z(t,y)dy\end{split}$
(C.6)

Now we insert (C.4),(C.6) into the equation satisfied by $w$:

$\partial_{t}w-\frac{\sigma}{{\overline{e}}(t)^{2}}\partial_{xx}^{2}w-\frac{{\overline{v}}}{{\overline{e}}(t)}x\partial_{x}w+\lambda
w=0.$

After cancellations, it remains, for any $x\in(0,1)$:

$\begin{split}0=\frac{\sigma}{{\overline{e}}(t)^{2}}\left(2\frac{d}{dx}{\tilde{k}}(t,x,x)+\lambda\right)z(x)-\partial_{y}{\tilde{k}}(t,x,0)z(0)\\\
-\int_{0}^{x}z\left(\partial_{t}{\tilde{k}}-\frac{\sigma}{{\overline{e}}(t)^{2}}(\partial^{2}_{xx}{\tilde{k}}-\partial_{yy}^{2}{\tilde{k}})-\frac{{\overline{v}}}{{\overline{e}}(t)}(x\partial_{x}{\tilde{k}}+\partial_{y}(y{\tilde{k}}))+\lambda{\tilde{k}}\right)dy,\end{split}$

which leads to the following problem for $t>0$

$\left\\{\begin{aligned}
\partial_{t}{\tilde{k}}-\frac{\sigma}{{\overline{e}}(t)^{2}}\left(\partial_{yy}^{2}{\tilde{k}}(t,x,y)-\partial_{xx}^{2}{\tilde{k}}(t,x,y)\right)-\frac{{\overline{v}}}{{\overline{e}}(t)}\left(x\partial_{x}{\tilde{k}}+y\partial_{y}{\tilde{k}}+{\tilde{k}}\right)+\lambda{\tilde{k}}&=0&(x,y)\in\\{0<y\leq
x<1\\},\\\
\frac{\sigma}{{\overline{e}}(t)^{2}}\partial_{y}{\tilde{k}}(t,x,0)&=0&x\in(0,1),\\\
\frac{2\sigma}{{\overline{e}}(t)^{2}}\frac{d}{dx}{\tilde{k}}(t,x,x)&=-\lambda&x\in(0,1).\end{aligned}\right.$
(C.7)

Now we look for a solution with separate variables under the form ($k$ does
not depend explicitly on time):

${\tilde{k}}(t,x,y)={\overline{e}}(t)k(x{\overline{e}}(t),y{\overline{e}}(t)).$
(C.8)

Inserting (C.8) into (C.7), the terms in ${\overline{e}}(t)$ cancel each other
and coming back to the original domain we obtain

$\left\\{\begin{aligned}
\partial_{xx}^{2}k(x,y)-\partial_{yy}^{2}k(x,y)&=\frac{\lambda}{\sigma}k(x,y)&(x,y)\in\\{0<y\leq
x<{\overline{e}}(t)\\},\\\ \partial_{y}k(x,0)&=0&x\in(0,{\overline{e}}(t)),\\\
k(x,x)&=-\frac{\lambda}{2\sigma}x&x\in(0,{\overline{e}}(t)),\end{aligned}\right.$
(C.9)

Moreover, inserting (C.8) into (C.3), it is clear with a change of variables
that the $k$ defined from (C.8) enables to recover the original kernel in
(4.2)-(4.4) we were looking for.

##  Appendix D Failure of the basic quadratic Lyapunov approach

We show why the common approach of directly using a basic quadratic Lyapunov
function would fail to provide exponential stabilization. In order to have a
proper basic quadratic Lyapunov function, we work on the rescaled system
(C.2). A basic quadratic Lyapunov function for the $L^{2}$ norm has the form,
for some positive function $f\in C^{2}((0,1))\cap C^{1}([0,1])$,

$V(z(t,\cdot))=\int_{0}^{1}f(x)z(t,x)^{2}dx.$ (D.1)

Let us take $V$ as a Lyapunov function candidate. By differentiating along
$C^{2}$ solutions of (C.2), we have

$\frac{d}{dt}V(z(t,\cdot))=\int_{0}^{1}2f(x)z\left[\frac{\sigma}{{\overline{e}}(t)^{2}}\partial_{xx}^{2}z+\frac{{\overline{v}}}{{\overline{e}}(t)}x\partial_{x}z\right]dx,$
(D.2)

Integrating by parts the first term gives

$\left[2f\frac{\sigma}{{\overline{e}}(t)^{2}}z\partial_{x}z\right]_{0}^{1}-\frac{\sigma}{{\overline{e}}(t)^{2}}\int_{0}^{1}\left[2f(x)(\partial_{x}z)^{2}+2f^{\prime}(x)z\partial_{x}z\right]dx,$

while the second terms gives

$\left[fx\frac{{\overline{v}}}{{\overline{e}}(t)}z^{2}\right]_{0}^{1}-\frac{2{\overline{v}}}{{\overline{e}}(t)}\int_{0}^{1}\left(z^{2}\left(f+xf^{\prime}\right)+fxz\partial_{x}z\right)dx.$

Note that the last term in the previous equation is the same term we
integrated by parts. Therefore putting everything together we obtain:

$\begin{split}\frac{d}{dt}V(z(t,\cdot))=\left(\left[2f\frac{\sigma}{{\overline{e}}(t)^{2}}z\partial_{x}z\right]_{0}^{1}+\left[fx\frac{{\overline{v}}}{{\overline{e}}(t)}z^{2}\right]_{0}^{1}\right)-\frac{\sigma}{{\overline{e}}(t)^{2}}\int_{0}^{1}&\left[2f(x)(\partial_{x}z)^{2}+2f^{\prime}(x)z\partial_{x}z\right]dx\\\
&-\frac{{\overline{v}}}{{\overline{e}}(t)}\int_{0}^{1}z^{2}\left(f+xf^{\prime}\right)dx,\end{split}$
(D.3)

which gives, using again an integration by parts for the second term in the
first integral and the boundary conditions of (C.2)

$\begin{split}\frac{d}{dt}V(z(t,\cdot))=&\left(\left[2f\frac{\sigma}{{\overline{e}}(t)^{2}}z\partial_{x}z\right]_{0}^{1}+\left[fx\frac{{\overline{v}}}{{\overline{e}}(t)}z^{2}\right]_{0}^{1}-\left[z^{2}f^{\prime}\frac{\sigma}{{\overline{e}}(t)^{2}}\right]_{0}^{1}\right)\\\
&-\frac{\sigma}{{\overline{e}}(t)^{2}}\int_{0}^{1}\left[2f(x)(\partial_{x}z)^{2}-f^{\prime\prime}(x)z^{2}\right]dx-\frac{{\overline{v}}}{{\overline{e}}(t)}\int_{0}^{1}z^{2}\left(f+xf^{\prime}\right)dx.\\\
=&\left(\frac{2}{{\overline{e}}(t)}f(1)\delta\psi(t)z(t,1)-\frac{{\overline{v}}}{{\overline{e}}(t)}z^{2}(t,1))+\frac{\sigma}{{\overline{e}}(t)^{2}}\left(z(t,0)^{2}f^{\prime}(0)-z(t,1)^{2}f^{\prime}(1)\right)\right)\\\
&-\frac{\sigma}{{\overline{e}}(t)^{2}}\int_{0}^{1}\left[2f(x)(\partial_{x}z)^{2}-f^{\prime\prime}(x)z^{2}\right]dx-\frac{{\overline{v}}}{{\overline{e}}(t)}\int_{0}^{1}z^{2}\left(f+xf^{\prime}\right)dx.\end{split}$
(D.4)

To have a Lyapunov function ensuring an exponential stability estimate, there
has to exist $\gamma>0$ such that the right-hand side is lower or equal than
$-\gamma V$ for any $t\in[0,T]$ and any solution of (C.2). From that point one
would typically require in the Lyapunov approach that for all
$t\in[0,+\infty)$ and $\;Z\in C^{2}([0,1])$,

$\begin{split}&\left(\frac{2}{{\overline{e}}(t)}f(1)\delta\psi(t)Z(1)-\frac{{\overline{v}}}{{\overline{e}}(t)}Z^{2}(1))+\frac{\sigma}{{\overline{e}}(t)^{2}}\left(Z(0)^{2}f^{\prime}(0)-Z(1)^{2}f^{\prime}(1)\right)\right)\\\
&-\int_{0}^{1}\left[2f(x)\frac{\sigma}{{\overline{e}}(t)^{2}}(\partial_{x}Z)^{2}+\left(\frac{{\overline{v}}}{{\overline{e}}(t)}(f+xf^{\prime})-f^{\prime\prime}(x)\frac{\sigma}{{\overline{e}}(t)^{2}}-\gamma
f(x)\right)Z^{2}\right]dx\leq 0,\end{split}$ (D.5)

In particular this would be true for any $Z\in C^{2}([0,1])$ with compact
support which implies that

$\begin{split}\int_{0}^{1}\left[2f(x)\frac{\sigma}{{\overline{e}}(t)^{2}}(\partial_{x}Z)^{2}+\left(\frac{{\overline{v}}}{{\overline{e}}(t)}(f+xf^{\prime})-f^{\prime\prime}(x)\frac{\sigma}{{\overline{e}}(t)^{2}}-\gamma
f(x)\right)Z^{2}\right]dx\geq 0.\end{split}$ (D.6)

Since this has to be true for any time and any $Z\in C_{c}^{2}([0,1])$, and
since ${\overline{e}}(t)\rightarrow+\infty$ when $t\rightarrow+\infty$, this
implies that for any $x\in(0,1)$,

$\begin{split}xf^{\prime}(x)-\gamma f(x)\geq 0,\end{split}$ (D.7)

but as $f\in C^{1}([0,1];(0,+\infty))$ this is impossible: indeed, denoting
$M=\sup_{[0,1]}(f^{\prime})\in\mathbb{R}$ and $m=\inf_{[0,1]}(f)>0$ this would
imply in particular that

$xM\geq\gamma m>0,\;\forall x\in(0,1),$ (D.8)

which would lead to a contradiction. Note that, although Lyapunov functionals
of the form (D.1) fail here, some other Lyapunov functionals ([22, 47] or
quadratic functionals with time-dependent weights) may manage to provide rapid
stabilization results for this system.

## References

  * [1] Luca Alasio, Helene Ranetbauer, Markus Schmidtchen, and Marie-Therese Wolfram.  Trend to equilibrium for systems with small cross-diffusion.  ESAIM: Mathematical Modelling and Numerical Analysis, 54(5):1661–1688, September 2020.
  * [2] Jean Auriol, Federico Bribiesca-Argomedo, David Bou Saba, Michael Di Loreto, and Florent Di Meglio.  Delay-robust stabilization of a hyperbolic pde–ode system.  Automatica, 95:494–502, 2018.
  * [3] Jean Auriol and Florent Di Meglio.  Minimum time control of heterodirectional linear coupled hyperbolic PDEs.  Automatica J. IFAC, 71:300–307, 2016.
  * [4] Jean Auriol and Florent Di Meglio.  Robust output feedback stabilization for two heterodirectional linear coupled hyperbolic pdes.  Automatica, 115:108896, 2020.
  * [5] Antonello Baccoli, Alessandro Pisano, and Yury Orlov.  Boundary control of coupled reaction–diffusion processes with constant parameters.  Automatica, 54:80–90, April 2015.
  * [6] Athmane Bakhta and Virginie Ehrlacher.  Cross-diffusion systems with non-zero flux and moving boundary conditions.  ESAIM: Mathematical Modelling and Numerical Analysis, 52(4):1385–1415, July 2018.
  * [7] Andras Balogh and Miroslav Krstic.  Infinite dimensional backstepping-style feedback transformations for a heat equation with an arbitrary level of instability.  European journal of control, 8(2):165–175, 2002.
  * [8] Georges Bastin, Jean-Michel Coron, and Amaury Hayat.  Diffusion and robustness of boundary feedback stabilization of hyperbolic systems.  preprint, 2022.
  * [9] Dejan M. Bošković, Andras Balogh, and Miroslav Krstić.  Backstepping in infinite dimension for a class of parabolic distributed parameter systems.  Math. Control Signals Systems, 16(1):44–75, 2003.
  * [10] Haïm Brezis.  Analyse fonctionnelle: théorie et applications.  Mathématiques appliquées pour le master. Dunod, Paris, nouvelle présentation, [nachdr.] edition, 2010.
  * [11] Brezis, Haïm.  Functional Analysis, Sobolev Spaces and Partial Differential Equations.  Springer International Publishing.
  * [12] Martin Burger, Marco Di Francesco, Jan-Frederik Pietschmann, and Bärbel Schlake.  Nonlinear Cross-Diffusion with Size Exclusion.  SIAM Journal on Mathematical Analysis, 42(6):2842–2871, January 2010.
  * [13] Christopher I. Byrnes and Alberto Isidori.  New results and examples in nonlinear feedback stabilization.  Systems Control Lett., 12(5):437–442, 1989.
  * [14] Li Chen and Ansgar Jüngel.  Analysis of a parabolic cross-diffusion population model without self-diffusion.  Journal of Differential Equations, 224(1):39–59, May 2006.
  * [15] Xiuqing Chen and Ansgar Jüngel.  When do cross-diffusion systems have an entropy structure?  arXiv:1908.06873 [math], August 2019.
  * [16] Jean-Michel Coron and Brigitte d’Andréa Novel.  Stabilization of a rotating body beam without damping.  IEEE Trans. Automat. Control, 43(5):608–618, 1998.
  * [17] Jean-Michel Coron, Ludovick Gagnon, and Morgan Morancey.  Rapid stabilization of a linearized bilinear 1-D Schrödinger equation.  J. Math. Pures Appl. (9), 115:24–73, 2018.
  * [18] Jean-Michel Coron, Amaury Hayat, Shengquan Xiang, and Christophe Zhang.  Stabilization of the linearized water tank system.  arXiv preprint arXiv:2103.08293, 2021.
  * [19] Jean-Michel Coron, Long Hu, Guillaume Olive, and Peipei Shang.  Boundary stabilization in finite time of one-dimensional linear hyperbolic balance laws with coefficients depending on time and space.  J. Differential Equations, 271:1109–1170, 2021.
  * [20] Jean-Michel Coron and Qi Lü.  Fredholm transform and local rapid stabilization for a Kuramoto-Sivashinsky equation.  J. Differential Equations, 259(8):3683–3729, 2015.
  * [21] Jean-Michel Coron and Hoai-Minh Nguyen.  Null Controllability and Finite Time Stabilization for the Heat Equations with Variable Coefficients in Space in One Dimension via Backstepping Approach.  Archive for Rational Mechanics and Analysis, 225(3):993–1023, September 2017.
  * [22] Jean-Michel Coron and Emmanuel Trélat.  Global steady-state controllability of one-dimensional semilinear heat equations.  SIAM journal on control and optimization, 43(2):549–569, 2004.
  * [23] L. Desvillettes, T. Lepoutre, A. Moussa, and A. Trescases.  On the Entropic Structure of Reaction-Cross Diffusion Systems.  Communications in Partial Differential Equations, 40(9):1705–1747, September 2015.
  * [24] Nicolás Espitia, Andrey Polyakov, Denis Efimov, and Wilfrid Perruquetti.  Boundary time-varying feedbacks for fixed-time stabilization of constant-parameter reaction–diffusion systems.  Automatica, 103:398–407, 2019.
  * [25] Lawrence Evans.  Partial Differential Equations, volume 19 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, Rhode Island, second edition, March 2010.
  * [26] Ludovick Gagnon, Amaury Hayat, Shengquan Xiang, and Christophe Zhang.  Fredholm transformation on laplacian and rapid stabilization for the heat equations.  preprint, 2021.
  * [27] Ludovick Gagnon, Pierre Lissy, and Swann Marx.  A fredholm transformation for the rapid stabilization of a degenerate parabolic equation.  preprint, 2020.
  * [28] Katharina Hopf and Martin Burger.  On multi-species diffusion with size exclusion.  Nonlinear Analysis, 224:113092, 2022.
  * [29] Long Hu, Florent Di Meglio, Rafael Vazquez, and Miroslav Krstic.  Control of homodirectional and general heterodirectional linear coupled hyperbolic PDEs.  IEEE Trans. Automat. Control, 61(11):3301–3314, 2016.
  * [30] Long Hu, Rafael Vazquez, Florent Di Meglio, and Miroslav Krstic.  Boundary exponential stabilization of 1-dimensional inhomogeneous quasi-linear hyperbolic systems.  SIAM J. Control Optim., 57(2):963–998, 2019.
  * [31] Mojtaba Izadi, Javad Abdollahi, and Stevan S. Dubljevic.  PDE backstepping control of one-dimensional heat equation with time-varying domain.  Automatica, 54:41–48, April 2015.
  * [32] Mojtaba Izadi and Stevan Dubljevic.  Backstepping output-feedback control of moving boundary parabolic pdes.  European Journal of Control, 21:27–35, 2015.
  * [33] Ansgar Jüngel.  The boundedness-by-entropy principle for cross-diffusion systems.  Nonlinearity, 28(6):1963–2001, June 2015.
  * [34] Ansgar Jüngel.  Entropy Methods for Diffusive Partial Differential Equations.  SpringerBriefs in Mathematics. Springer International Publishing, Cham, 2016.
  * [35] Ansgar Jüngel and Ines Viktoria Stelzer.  Entropy structure of a cross-diffusion tumor-growth model.  Mathematical Models and Methods in Applied Sciences, 22(07):1250009, July 2012.
  * [36] Ansgar Jüngel and Ines Viktoria Stelzer.  Existence Analysis of Maxwell–Stefan Systems for Multicomponent Mixtures.  SIAM Journal on Mathematical Analysis, 45(4):2421–2440, January 2013.
  * [37] Simon Kerschbaum and Joachim Deutscher.  Backstepping control of coupled linear parabolic pdes with space and time dependent coefficients.  IEEE Transactions on Automatic Control, 65(7):3060–3067, 2019.
  * [38] Daniel E Koditschek.  Adaptive techniques for mechanical systems.  Proc. 5th. Yale University Conference, page pp. 259–265, 1987\.
  * [39] Miroslav Krstic and Andrey Smyshlyaev.  Backstepping boundary control for first-order hyperbolic PDEs and application to systems with actuator and sensor delays.  Systems Control Lett., 57(9):750–758, 2008.
  * [40] Miroslav Krstic and Andrey Smyshlyaev.  Boundary Control of PDEs: A Course on Backstepping Designs, volume 16 of Advances in Design and Control.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2008.
  * [41] Weijiu Liu.  Boundary feedback stabilization of an unstable heat equation.  SIAM journal on control and optimization, 42(3):1033–1043, 2003\.
  * [42] Jeremy Quastel.  Diffusion of color in the simple exclusion process.  Communications on Pure and Applied Mathematics, 45(6):623–679, 1992\.
  * [43] Andrey Smyshlyaev and Miroslav Krstic.  Closed-form boundary state feedbacks for a class of 1-D partial integro-differential equations.  IEEE Trans. Automat. Control, 49(12):2185–2202, 2004.
  * [44] Andrey Smyshlyaev and Miroslav Krstic.  On control design for pdes with space-dependent diffusivity or time-dependent reactivity.  Automatica, 41(9):1601–1608, 2005.
  * [45] John Tsinias.  Sufficient Lyapunov-like conditions for stabilization.  Math. Control Signals Systems, 2(4):343–357, 1989.
  * [46] Rafael Vazquez, Jean-Michel Coron, Miroslav Krstic, and Georges Bastin.  Local exponential $H^{2}$ stabilization of a $2\times 2$ quasilinear hyperbolic system using backstepping.  50th IEEE Conference on Decision and Control and European Control Conference, Orlando, pages 1329–1334, 2011.
  * [47] Shengquan Xiang.  Quantitative rapid and finite time stabilization of the heat equation.  arXiv preprint arXiv:2010.04696, 2020.
  * [48] Nicola Zamponi and Ansgar Jüngel.  Analysis of degenerate cross-diffusion population models with volume filling.  Annales de l’Institut Henri Poincaré C, Analyse non linéaire, 34(1):1–29, January 2017.
  * [49] Christophe Zhang.  Internal rapid stabilization of a 1-D linear transport equation with a scalar feedback.  Preprint, October 2018.  working paper or preprint.
  * [50] Christophe Zhang.  Finite-time internal stabilization of a linear 1-D transport equation.  Systems Control Lett., 133:104529, 8, 2019.


This paper introduces the Variational Determinant Estimator (VDE), a variational extension of the recently proposed determinant estimator discovered by Sohl-Dickstein, 2020. Our estimator significantly reduces the variance even for low sample sizes by combining (importance-weighted) variational inference and a family of normalizing flows which allow density estimation on hyperspheres. In the ideal case of a tight variational bound, the VDE becomes a zero variance estimator, and a single sample is sufficient for an exact (log) determinant estimate.

§ INTRODUCTION

§ THE VARIATIONAL DETERMINANT ESTIMATOR

§ RELATED WORK

§ RESULTS

In this section, we demonstrate the performance of our method in determinant estimation. We consider two cases: estimating the determinant of randomly sampled $10 \times 10$ dense matrices and estimating the determinant of a convolutional layer.

Nr. of samples     $10^2$     $10^3$     $10^4$     $10^5$

VDE det. (ours)     3.4 $\pm$ 2.1 %      1.7 $\pm$ 0.6 %     1.6 $\pm$ 1.3 %     0.3 $\pm$ 0.3 %

MC det.     533 $\pm$ 660  %     348 $\pm$ 262 %     104 $\pm$ 30 %      59 $\pm$  43 %

Comparison of determinant estimates. Results are in mean absolute relative difference for $5$ unit Gaussian sampled $10 \times 10$ matrices. Deviations are given in one standard deviation.

§.§ Dense Matrix

The determinant is estimated for five $10 \times 10$ matrices where the entries are sampled from unit Gaussians, see the Appendix <ref> for the specific matrices. The spherical flow utilizes a Moebius transformation for the circle part with $N_C$ = 12 number of centers and Neural Spline flows [Durkan et al., 2019] with $N_B = 16$ number of bins for the interval part, see [Rezende et al., 2020] for details of the architecture and the parameters. Furthermore, we stacked $N_F=8$ flows on top, used coupling layers, and trained the models for $10k$ iterations with a batch size of $1024$. Flows based on autoregressive masking are also possible.

The results can be seen in Table <ref>, where we present the mean of the relative absolute differences of the estimated determinant in comparison to the true absolute determinant. The variational determinant estimator achieves even for a low sample size of $10^2$ low relative differences, whereas in contrast, the naïve Monte Carlo estimate still has high variance throughout all sample sizes.

[Absolute determinant estimates.]

[Relative variation of determinant estimates]
Determinant estimates of a structured $16 \times 16$ matrix. Left: absolute determinant estimated values. Right: the absolute relative difference to the true determinant value. Both plots are on log-scale.

Nr. of samples      $10^2$     $10^3$     $10^4$     $10^5$      true (log) determinant

VDE det. (ours)     7.62     7.71     7.64     7.70      2*7.71

MC det.      481.09      74.37      39.08      13.21

VDE log det. (ours)      2.03     2.04     2.03     2.04      2*2.04

MC log det.     6.18      4.31      3.67     2.58

VDE Rel. diff of det. (ours)      1.1 %     0.05 %     0.9 %     0.1 %     2*0 %

MC Rel. diff of det.     6144 %      865 %      407 %      71 %

Comparison of variational and Monte Carlo determinant estimates. First four rows show results in absolute numbers and log-space. The last two rows illustrate the relative difference of the estimates to the true determinant.

§.§ Convolutional Layer

In this experiment, the (log) determinant of a convolutional layer is estimated. The reason for this experiment is that these types of sparse linear transformations often occur in deep learning, and they typically have cheap matrix-vector products. We manually reconstructed the $16 \times 16$ equivalent matrix $W$ of a convolution of a $3 \times 3 $ filter with an $ 4 \times 4$ image, see Appendix <ref> for the filter. The determinant of $W$ equals the determinant of the convolution operation. We chose the same architecture as in our previous experiment but ran the experiment this time for $40k$ iterations with the same batch size. Note that the parametrization of our spherical flow is fully connected, and better optimization behavior is expected when the parametrization would be convolutional.

Figure <ref> and Table <ref> show the results of the experiment in terms of absolute and log abs. determinant estimates. We observe the same behavior as in the previous experiment: The variational determinant estimator achieves low errors with already low sample sizes, whereas the MC determinant estimator is not able to capture the determinant correctly with even high sample sizes.

§ CONCLUSION

In this paper, we introduced the Variational Determinant Estimator, which achieves with low sample sizes high accuracy in estimating a determinant of a linear operator. Interestingly, the estimator in its original variant and the VDE allows estimation if only matrix-vector products are available.

In our experiments, we considered an offline setting where a large number of samples are required first to optimize the model. However, in future work, VDE could also be applied in an online, moving target settings. In this case, small updates to the matrix $A$ would only require small updates to the density model $q(s)$. A perpendicular direction for VDE would be to estimate the Jacobian determinant of a function $f$. The Jacobian would depend on the input $x$, and the density model $q(s|x)$ can be amortized. Additionally, not the entire Jacobian but only Jacobian-vector products would be required to estimate the determinant.

Code to reproduce the results and to enable further research concerning the Variational Determinant Estimator and spherical normalizing flows will be made public in the git repository[https://github.com/P4ppenheimer] of one of the authors at a later time point.

[Burda et al., 2015]
Yuri Burda, Roger Grosse, and Ruslan Salakhutdinov.
Importance weighted autoencoders.
arXiv preprint arXiv:1509.00519, 2015.

[Chen et al., 2019]
Tian Qi Chen, Jens Behrmann, David Duvenaud, and Jörn-Henrik Jacobsen.
Residual flows for invertible generative modeling.
In Neural Information Processing Systems, NeurIPS, pages
9913–9923, 2019.

[Dinh et al., 2016]
Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio.
Density estimation using real nvp.
arXiv preprint arXiv:1605.08803, 2016.

[Durkan et al., 2019]
Conor Durkan, Artur Bekasov, Iain Murray, and George Papamakarios.
Neural spline flows.
In Advances in Neural Information Processing Systems, pages
7511–7522, 2019.

[Gemici et al., 2016]
Mevlana C. Gemici, Danilo Jimenez Rezende, and Shakir Mohamed.
Normalizing flows on riemannian manifolds.
CoRR, abs/1611.02304, 2016.

[Goliński et al., 2019]
Adam Goliński, Frank Wood, and Tom Rainforth.
Amortized monte carlo integration.
arXiv preprint arXiv:1907.08082, 2019.

[Hesterberg, 1988]
Timothy Classen Hesterberg.
Advances in importance sampling.
PhD thesis, Citeseer, 1988.

[Kingma and Welling, 2014]
Diederik P. Kingma and Max Welling.
Auto-encoding variational bayes.
In 2nd International Conference on Learning Representations,
ICLR, 2014.

[Kingma et al., 2016]
Diederik P Kingma, Tim Salimans, Rafal Jozefowicz, Xi Chen, Ilya Sutskever, and
Max Welling.
Improved variational inference with inverse autoregressive flow.
In Advances in Neural Information Processing Systems, pages
4743–4751, 2016.

[Kingma and Dhariwal, 2018]
Durk P Kingma and Prafulla Dhariwal.
Glow: Generative flow with invertible 1x1 convolutions.
In Advances in neural information processing systems, pages
10215–10224, 2018.

[Müller et al., 2019]
Thomas Müller, Brian McWilliams, Fabrice Rousselle, Markus Gross, and Jan
Neural importance sampling.
ACM Transactions on Graphics (TOG), 380 (5):0
1–19, 2019.

[Owen, 2013]
Art B Owen.
Monte carlo theory.
Methods and Examples, 665, 2013.

[Perugachi-Diaz et al., 2020]
Yura Perugachi-Diaz, Jakub M. Tomczak, and Sandjai Bhulai.
CoRR, abs/2010.02125, 2020.

[Rezende and Mohamed, 2015]
Danilo Jimenez Rezende and Shakir Mohamed.
Variational inference with normalizing flows.
arXiv preprint arXiv:1505.05770, 2015.

[Rezende et al., 2014]
Danilo Jimenez Rezende, Shakir Mohamed, and Daan Wierstra.
Stochastic backpropagation and approximate inference in deep
generative models.
In Proceedings of the 31th International Conference on Machine
Learning, ICML, 2014.

[Rezende et al., 2020]
Danilo Jimenez Rezende, George Papamakarios, Sébastien Racanière,
Michael S Albergo, Gurtej Kanwar, Phiala E Shanahan, and Kyle Cranmer.
Normalizing flows on tori and spheres.
arXiv preprint arXiv:2002.02428, 2020.

[Sohl-Dickstein, 2020]
Jascha Sohl-Dickstein.
Two equalities expressing the determinant of a matrix in terms of
expectations over matrix-vector products.
arXiv preprint arXiv:2005.06553, 2020.

[Tabak and Turner, 2013]
Esteban G Tabak and Cristina V Turner.
A family of nonparametric density estimation algorithms.
Communications on Pure and Applied Mathematics, 660
(2):0 145–164, 2013.

§ EXPERIMENTAL DETAILS

§.§ Cover Density

The cover image of this paper illustrates the learned proposal distribution $q(s)$ corresponding to the matrix

\begin{align*}
A = \begin{bmatrix}
\begin{array}{rrrr}
-0.7056  & 0.6741 &-0.5454  \\
0.9107  & 1.0682 & 0.1424 \\
-1.2754  & -0.1769 & 1.0084
\end{array}
\end{bmatrix},
\end{align*}

which is created with  and torch manual seed $15$. The optimal proposal distribution $q^* \propto \lVert As\rVert^{-n}$ is illustrated in Figure <ref>. We trained the model for $10$k iterations and in contrast to the architecture in Section <ref>, we used $N_F = 6$ flows with autoregressive masking and Neural Spline flows for both the spherical part and the interval part of $\mathbb{S}^1 \times [-1,1]$ with $N_B = 32$ bins, see again [Rezende et al., 2020] for details.

Optimal proposal distribution corresponding to $A$.

§.§ Dense $10 \times 10$ Matrices

In this section we publish the absolute determinants in Table <ref> and the $10 \times 10$ matrices of our experiment in Section <ref>. Numbers are rounded to two decimals.

$A_1$     $A_2$     $A_3$     $A_4$     $A_5$

Absolute det.     520.36    748.68     945.02     3000.5     252.29

Absolute determinants of matrices of the experiment in Section <ref>.

\begin{align*}
A_1 &= \begin{bmatrix}
\begin{array}{rrrrrrrrrr}
-1.08 & -0.6 & 0.06 & 0.71 & -0.81 & 0.57 & 0.69 & 0.51 & -0.94 & 0.18\\
-0.55 & 1.5 & 1.39 & -0.18 & -0.56 & -0.05 & 0.98 & 1.82 & 1.48 & 0.01\\
-0.26 & -2.07 & -1.12 & -0.27 & -1.03 & 0.97 & -1.84 & -0.5 & -0.47 & -1.17\\
1.01 & -1.25 & 1.71 & 1.24 & -0.79 & -0.17 & -1.05 & 0.44 & 0.02 & 0.04\\
1.24 & -0.31 & -0.18 & -0.74 & -0.43 & 0.29 & -0.67 & 1.43 & -1.01 & -0.17\\
-0.49 & -1.17 & 0.43 & 1.4 & 1.28 & 1.8 & -0.45 & 1.67 & -0.93 & -1.72\\
0.78 & 1.19 & 0.02 & -0.06 & 0.72 & -1.24 & -1.19 & -0.71 & 1.73 & 0.81\\
0.53 & 1.56 & -1.09 & 0.33 & -0.29 & -0.47 & 1.02 & 1.67 & -0.17 & 0.26\\
1.16 & -0.18 & 0.86 & 0.94 & 0.26 & -1.64 & -0.38 & -0.31 & -0.79 & 1.31\\
0.54 & 1.39 & -0.21 & -0.12 & 0.14 & 0.8 & 0.78 & 0.85 & -1.3 & -0.41\\
\end{array}
\end{bmatrix} \\[3pt]
A_2 &= \begin{bmatrix}
\begin{array}{rrrrrrrrrr}
-1.92 & -0.19 & 0.34 & 0.41 & -0.58 & -2.08 & 0.29 & -0.46 & -1.37 & -0.45\\
-0.56 & 0.71 & 0.06 & 0.17 & 1.44 & -1.81 & -1.19 & 1.02 & -2.84 & 2.28\\
1.64 & 0.14 & -1.86 & 0.23 & 0.85 & 1.33 & -0.88 & -0.73 & -0.53 & 2.09\\
-0.11 & -0.43 & 0.68 & -1.45 & 0.08 & 0.81 & 0.53 & 0.41 & 0.41 & -0.27\\
-0.05 & 0.05 & 0.7 & -1.09 & 1.77 & -0.79 & -0.35 & 1.71 & 0.85 & 0.8\\
1.24 & -0.22 & 0.41 & -1.02 & -0.64 & -0.21 & -1.25 & 0.71 & 0.6 & -0.75\\
0.71 & -0.91 & -0.11 & 0.18 & 1.13 & -0.48 & 1.85 & -0.03 & 0.29 & -1.25\\
0.52 & -1.06 & 0.48 & -2.26 & 1.52 & -0.63 & 1.26 & -1.42 & -0.02 & -1.66\\
-1.01 & -1.23 & 0.42 & -0.37 & 1. & -0.04 & -0.32 & 0.52 & -1.91 & -1.78\\
0.89 & -0.1 & -0.39 & -0.52 & 0.21 & -0.99 & 0.48 & 0.22 & 0.77 & -0.19\\
\end{array}
\end{bmatrix} \\[3pt]
A_3 &= \begin{bmatrix}
\begin{array}{rrrrrrrrrr}
-0.15 & -1.65 & -0.95 & 0.26 & 1.35 & -0.1 & 0.37 & 0.45 & 0.23 & -1.12\\
0.61 & -1.81 & -0.68 & 0.58 & 0.94 & 2.36 & -0.49 & 0.04 & 0.86 & 0.52\\
1.91 & -1.44 & -0.51 & 0.96 & -2.56 & -0.01 & -1.13 & 0.19 & -2.5 & 0.68\\
0.93 & -1.3 & -0.65 & -1.9 & -0.09 & 0.24 & 0.69 & 1.28 & -0.57 & -0.39\\
0.55 & -1.34 & -1.37 & -1.29 & -0.28 & -0.67 & 0.77 & -0.25 & 0.85 & -2.89\\
-0.86 & 1.95 & -1.33 & 0.68 & 0.27 & 0.25 & 0.29 & -1.29 & 2.05 & 0.11\\
0.82 & 0.52 & -0.71 & -0.59 & -1.57 & -1.05 & 0.46 & -0.61 & 0.63 & 2.02\\
0.76 & 0.01 & -0.06 & -0.43 & 1.12 & 1.05 & -1.35 & -0.04 & -0.62 & -0.35\\
-2.13 & -0.8 & 1.12 & 1.77 & -0.79 & -0.1 & 1.17 & -1.06 & -0.37 & 0.01\\
-0.86 & 1.44 & -0.55 & 1.19 & 2.52 & 0.81 & -0.36 & -0.61 & 1.24 & -0.06\\
\end{array}
\end{bmatrix} \\[3pt]
A_4 &= \begin{bmatrix}
\begin{array}{rrrrrrrrrr}
-0.31 & -0.68 & -0.22 & -0.28 & 1.64 & -0.41 & -0.66 & -0.59 & 1.57 & 0.38\\
-0.15 & 0.6 & 1.08 & 1.29 & -0.12 & 1.89 & -1.85 & -0.11 & 1.5 & 0.72\\
0.88 & -1.71 & 0.69 & -1.75 & -0.06 & 0.9 & 0.08 & -0.11 & -0.21 & 1.75\\
0.31 & -1.4 & -1.79 & 0.17 & 0.57 & -0.86 & 1.64 & -1.55 & 0.91 & -2.06\\
1.1 & -1.19 & 0.47 & -0.84 & 0.37 & 0.25 & 0.03 & -0.23 & 1.32 & 0.36\\
0.43 & -0.02 & -0.04 & 1.19 & 0.2 & -1.13 & 1.36 & 1.23 & -0.01 & 2.08\\
-0.8 & 0.48 & -1.57 & 0.6 & -0.19 & -0.18 & -0.88 & -1.53 & -0.66 & -0.83\\
1.32 & -1.09 & 0.71 & 1.04 & 1.02 & -0.09 & 1.51 & -0.51 & -0.73 & -0.82\\
0.21 & -2.07 & 0.61 & 0.29 & 1.41 & -1.93 & -2.06 & 0.23 & -0.09 & 0.24\\
-1.36 & 0.3 & 0.15 & 1.33 & -1.1 & -0.72 & 0.37 & 0.09 & -0.56 & 2.81\\
\end{array}
\end{bmatrix} \\[3pt]
A_5 &= \begin{bmatrix}
\begin{array}{rrrrrrrrrr}
0.36 & -0.95 & 0.12 & 0.85 & -0.4 & -0.1 & 0.58 & -0.48 & 0.79 & 0.12\\
0.11 & -0.02 & -0.66 & -0.98 & -0.28 & -1.61 & -0.82 & 1.13 & 1.18 & 0.33\\
-0.7 & 0.65 & -1.5 & -0.33 & -0.18 & -0.6 & -0.84 & -0.43 & -0.42 & 1.12\\
-2.33 & -0.49 & 0.61 & 0.88 & -0.85 & -0.68 & 0.38 & 0.53 & 0.34 & 1.59\\
0.43 & 1.61 & -0.14 & 1.15 & -1.25 & 2.28 & -0.32 & -0.36 & -2.1 & 0.98\\
-0.68 & -0.54 & -0.88 & 1.55 & 0.7 & -1.34 & 0.15 & -0.27 & -0.86 & 1.35\\
-0.83 & -0.52 & -0.83 & -1.98 & 1.79 & -0.86 & 0.05 & 1.29 & 0.1 & 1.17\\
-1.34 & -0.66 & 0.12 & -0.95 & -0.46 & 2.15 & -0.67 & -0.77 & 1.87 & 1.4\\
0.54 & -0.51 & 0.16 & 1.38 & 1.49 & 0.61 & 0.22 & 0.64 & -0.27 & -0.47\\
0.62 & -0.24 & -0.11 & 0.27 & -0.48 & 0.75 & 0.59 & 0.41 & -0.81 & 0.07\\
\end{array}
\end{bmatrix}
\end{align*}

§.§ Filter of Convolutional Layer

Here we present the $3 \times 3$ filter $k$ which was used to create the $16 \times 16$ matrix $W$ of our experiment in section <ref>. The matrix $W$ is the equivalent matrix of the convolution $k \star x$ where $x$ is an arbitrary $4 \times 4$ image.

\begin{align*}
k = \begin{bmatrix}
\begin{array}{rrr}
-0.107 & -0.689 & -0.027\\
0.226 & 1.393 & -0.544\\
-0.28 & -0.467 & 0.024\\
\end{array}
\end{bmatrix}
\end{align*}