# A Carleman-based numerical method for quasilinear elliptic equations with
over-determined boundary data and applications

Thuy T. Le Department of Mathematics and Statistics, University of North
Carolina at Charlotte, Charlotte, NC, 28223, USA<EMAIL_ADDRESS>Loc H.
Nguyen Department of Mathematics and Statistics, University of North Carolina
at Charlotte, Charlotte, NC, 28223, USA<EMAIL_ADDRESS>Hung V. Tran
Department of Mathematics, University of Wisconsin Madison, Madison, WI,
53706, USA<EMAIL_ADDRESS>(corresponding author).

###### Abstract

We propose a new iterative scheme to compute the numerical solution to an
over-determined boundary value problem for a general quasilinear elliptic PDE.
The main idea is to repeatedly solve its linearization by using the quasi-
reversibility method with a suitable Carleman weight function. The presence of
the Carleman weight function allows us to employ a Carleman estimate to prove
the convergence of the sequence generated by the iterative scheme above to the
desired solution. The convergence of the iteration is fast at an exponential
rate without the need of an initial good guess. We apply this method to
compute solutions to some general quasilinear elliptic equations and a large
class of first-order Hamilton-Jacobi equations. Numerical results are
presented.

Key words: numerical methods; Carleman estimate; linearization; boundary value
problems; quasilinear elliptic equations; Hamilton-Jacobi equations; viscosity
solutions; vanishing viscosity process.

AMS subject classification: 35D40, 35F30, 35J62, 35N25, 65N12, 78A46.

##  1 Introduction

The main aim of this paper is to develop a numerical method based on Carleman
estimates to solve quasilinear elliptic PDEs with over-determined boundary
data. We consider this new method the second generation of Carleman-based
numerical methods while the first generation is called the convexification,
which will be mentioned in detail later. Let $\Omega$ be an open and bounded
domain in $\mathbb{R}^{d}$, $d\geq 2$, with smooth boundary $\partial\Omega$.
Let $f$ and $g$ be two smooth functions defined on $\partial\Omega$. Let
$F:\overline{\Omega}\times\mathbb{R}\times\mathbb{R}^{d}\to\mathbb{R}$ be a
function in the class $C^{2}.$ Let
$A=(a_{ij})_{i,j=1}^{d}:\overline{\Omega}\to\mathbb{R}^{d\times d}$ be
$C^{2}$, symmetric, and positive definite, that is,

$\gamma|\xi|^{2}\leq a_{ij}({\bf
x})\xi_{i}\xi_{j}\leq\gamma^{-1}|\xi|^{2}\quad\text{ for all }{\bf
x}\in\overline{\Omega},\ \xi\in\mathbb{R}^{d},$

for some fixed $\gamma\in(0,1)$. The following problem is of our interests.

######  Problem 1.1.

Assume that the over-determined boundary value problem

$\left\\{\begin{array}[]{ll}-{\rm div}(A({\bf x})\nabla u({\bf x}))+F({\bf
x},u({\bf x}),\nabla u({\bf x}))=0&{\bf x}\in\Omega,\\\ u({\bf x})=f({\bf
x})&{\bf x}\in\partial\Omega,\\\ \partial_{\nu}u({\bf x})=g({\bf x})&{\bf
x}\in\partial\Omega\end{array}\right.$ (1.1)

has a solution $u^{*}$ in $C^{2}(\overline{\Omega})$. Compute the function
$u^{*}$.

Problem 1.1 is motivated by a class of nonlinear inverse problems in PDEs. One
important goal of inverse problems is to reconstruct the internal structure of
a domain from boundary measurements, which allow us to impose both Dirichlet
and Neumann data of the unknown in (1.1). Recently, a unified framework to
solve such inverse problems was developed by the research group of the first
and second authors, which has two main steps. In the first step, by
introducing a change of variables, one derives a PDE of the form (1.1) from
the given inverse problem, in which $f$ and $g$ can be computed directly by
the given boundary data. In the second step, one numerically solves (1.1) to
find $u^{*}$. The knowledge of $u^{*}$ directly yields that of the solution to
the corresponding inverse problem under consideration. See [14, 26] and the
references therein for some works in this framework. Moreover, this unified
framework was successfully tested with experimental data in [12, 13, 21].
Another motivation to study Problem 1.1 is to seek solutions to Hamilton-
Jacobi equations under the circumstance that the Neumann data of the unknown
can be computed by its Dirichlet data and the given form of the Hamiltonian,
see e.g., [23, Assumption 1.1 and Remark 1.1]. Since inverse problems are out
of the scope of this paper, we only focus on the applications in solving
quasilinear elliptic PDEs and first-order Hamilton-Jacobi equations.

A natural approach to solve (1.1) is based on optimization. That means one
sets the computed solution to (1.1) as a minimizer of a mismatch functional,
e.g.,

$v\mapsto J(v):=\int_{\Omega}\big{|}-{\rm div}(A({\bf x})\nabla v({\bf
x}))+F({\bf x},v({\bf x}),\nabla v({\bf x}))\big{|}^{2}\,d{\bf x}+\mbox{a
regularization term}$

subject to the Cauchy boundary conditions $v|_{\partial\Omega}=f$ and
$\partial_{\nu}v|_{\partial\Omega}=g.$ The methods based on optimization are
widely used in the scientific community, especially in computational
mathematics, physics and engineering. Although effective and popular, the
optimization-based approaches have some drawbacks:

  1. 1.

In general, it is not clear that the obtained minimizer approximates the true
solution to (1.1).

  2. 2.

The mismatch functional $J$ is not convex, and it might have multiple minima
and ravines (see an example in [36] for illustration). To deliver reliable
numerical solutions, one must know some good initial guesses of the true
solutions.

  3. 3.

The computation is expensive and time consuming.

Drawbacks # 1 and #2 can be treated by the convexification method, which is
designed to globalize the optimization methods. The main idea of the
convexification method is to employ some suitable Carleman weight functions to
convexify the mismatch functionals. The convexity of weighted mismatch
functionals is rigorously proved by Carleman estimates. Several versions of
the convexification method have been developed in [1, 14, 15, 16, 18, 21]
since it was first introduced in [20]. Moreover, we recently discovered that
the convexification method can be used to solve a large class of first-order
Hamilton-Jacobi equations [23].

In this paper, we introduce a new method to solve (1.1) based on linearization
and Carleman estimates. Like the convexification method, our method delivers a
reliable solution to (1.1) without requiring a good initial guess. This fact
is rigorously proved. Unlike the convexification method which is time
consuming, our new method quickly provides the desired solutions. Its converge
rate is $O(\theta^{n})$ as $n\to\infty$ for some $\theta\in(0,1).$

We find the numerical solution to (1.1) by repeatedly solving the
linearization of (1.1) by a new “Carleman weighted” quasi-reversibility
method. The classical quasi-reversibility method was first proposed in [24],
and it has been studied intensively since then (see [17] for a survey). By a
Carleman weighted quasi-reversibility method, we mean that we let a suitable
Carleman weight function involve in the cost functional suggested by the
classical quasi-reversibility method. The presence of the Carleman weight
function is the key for us to prove our convergence theorem. Our process to
solve Problem 1.1 is as follows. We first choose any initial solution that
might be far away from the true one. Denote this initial solution by the
function $u_{0}$. Linearizing (1.1) about $u_{0}$, we obtain a linear PDE. We
then solve this linear PDE by the Carleman weighted quasi-reversibility method
to obtain an updated solution $u_{1}$. Using the Carleman weighted quasi-
reversibility method rather than the classical one in this step is the key to
our success. By iteration, we repeat this step to construct a sequence
$\\{u_{n}\\}_{n\geq 0}$. The convergence of this sequence to the true solution
to (1.1) is proved by using Carleman estimates. We then apply this method to
numerically solve some quasilinear elliptic equations. It is important to note
that our approach works well for systems of quasilinear elliptic PDEs too.

Next, we use our method to solve some first-order Hamilton-Jacobi equations.
More precisely, to find viscosity solutions to the first-order equation

$F({\bf x},u({\bf x}),\nabla u({\bf x}))=0\quad{\bf x}\in\Omega,$

we use the vanishing viscosity procedure and consider, for $0<\epsilon_{0}\ll
1$,

$-\epsilon_{0}\Delta u({\bf x})+F({\bf x},u({\bf x}),\nabla u({\bf
x}))=0\quad{\bf x}\in\Omega,\\\ $

with given Cauchy boundary data. Our new method is robust in the sense that it
works for general nonlinearity $F({\bf x},u,\nabla u)$ that might not be
convex in $\nabla u.$ We refer the readers to [8, 7, 39] and the references
therein for the theory of viscosity solutions. A weakness of our new approach
in computing the viscosity solutions to Hamilton-Jacobi equations is that we
need to require both Dirichlet and Neumann data of the unknown $u$. See [23,
Remark 1.1] for some circumstances that this requirement is fulfilled. There
have been many important methods to solve Hamilton-Jacobi equations in the
literature. For finite difference monotone and consistent schemes of first-
order equations and applications, see [2, 9, 34, 37, 38] for details and
recent developments. If $F=F({\bf x},u,\nabla u)$ is convex in $\nabla u$ and
satisfies appropriate conditions, it is possible to construct some semi-
Lagrangian approximations by the discretization of the Dynamical Programming
Principle associated to the problem, see [10, 11] and the references therein.

The paper is organized as follows. In Section 2, we recall a Carleman estimate
and an example in which Problem 1.1 appears. In Section 3, we introduce the
new iterative method based on linearization and Carleman estimates. In Section
4, we prove the convergence of our method. Some numerical results for
quasilinear elliptic equations and first-order Hamilton-Jacobi equations are
presented in Section 5. Concluding remarks are given in Section 6.

##  2 Preliminaries

In this section, we recall a Carleman estimate, which plays a key role for the
proof of the convergence theorem in this paper. We then present an inverse
scattering problem in which Problem 1.1 appears.

###  2.1 A Carleman estimate

In this section, we present a simple form of Carleman estimates. Carleman
estimates were first employed to prove the unique continuation principle, see
e.g., [5, 35], and they quickly became a powerful tool in many areas of PDEs
afterwards. Let ${\bf x}_{0}$ be a point in
$\mathbb{R}^{d}\setminus\overline{\Omega}$ such that $r({\bf x})=|{\bf x}-{\bf
x}_{0}|>1$ for all ${\bf x}\in\Omega.$ For each $\beta>0$, define

$\mu_{\beta}({\bf x})=r^{-\beta}({\bf x})=|{\bf x}-{\bf
x}_{0}|^{-\beta}\quad\mbox{for all }{\bf x}\in\overline{\Omega}.$ (2.1)

We have the following lemma.

######  Lemma 2.1 (Carleman estimate).

There exist positive constants $\beta_{0},\lambda_{0}$ depending only on ${\bf
x}_{0}$, $\Omega$, $\gamma$, and $d$ such that for all function $v\in
C^{2}(\overline{\Omega})$ satisfying

$v({\bf x})=\partial_{\nu}v({\bf x})=0\quad\mbox{for all }{\bf
x}\in\partial\Omega,$ (2.2)

the following estimate holds true

$\int_{\Omega}e^{2\lambda\mu_{\beta}({\bf x})}|{\rm div}(A\nabla
v)|^{2}\,d{\bf x}\geq C\lambda\int_{\Omega}e^{2\lambda\mu_{\beta}({\bf
x})}|\nabla v({\bf x})|^{2}\,d{\bf
x}+C\lambda^{3}\int_{\Omega}e^{2\lambda\mu_{\beta}({\bf x})}|v({\bf
x})|^{2}\,d{\bf x}$ (2.3)

for all $\beta\geq\beta_{0}$ and $\lambda\geq\lambda_{0}$. Here, $C=C({\bf
x}_{0},\Omega,\gamma,d,\beta)>0$ depends only on the listed parameters.

###### Proof.

Lemma 2.1 is a direct consequence of [30, Lemma 5]. Let $0<R_{1}<1$ and
$R_{3}\gg 1$ such that $\Omega\Subset B({\bf
x}_{0},R_{3})\setminus\overline{B({\bf x}_{0},R_{1})}$. Here, $B({\bf
x}_{0},s)=\\{{\bf y}\in\mathbb{R}^{d}:|{\bf y}-{\bf x}_{0}|<s\\}$ for $s>0.$
Extend $v$ to the whole $\mathbb{R}^{d}$ such that $v({\bf x})=0$ for all
${\bf x}\in\mathbb{R}^{d}\setminus\Omega$. Using a change of variable ${\bf
x}\mapsto{\bf x}-{\bf x}_{0}$ and [30, Lemma 5], there exists a number
$\beta_{0}\geq 1$ depending on $\gamma$, $R_{1}$ and $R_{3}$ such that for all
$\beta\geq\beta_{0}$ and $|\lambda|\geq 2R_{3}^{-\beta}$, we have

$\int_{B({\bf x}_{0},R_{3})\setminus\overline{B({\bf x}_{0},R_{1})}}r({\bf
x})^{\beta+2}e^{2\lambda r({\bf x})^{-\beta}}|{\rm div}(A\nabla v)|^{2}\,d{\bf
x}\\\ \geq C\int_{B({\bf x}_{0},R_{3})\setminus\overline{B({\bf
x}_{0},R_{1})}}e^{2\lambda r^{-\beta}({\bf
x})}|\lambda|\beta\big{(}\beta^{3}\lambda^{2}r^{-2\beta-2}({\bf
x})|v|^{2}+|\nabla v|^{2}\big{)}\,d{\bf x}\\\ +C\int_{\partial(B({\bf
x}_{0},R_{3})\setminus\overline{B({\bf x}_{0},R_{1})})}|\lambda|\beta
e^{2\lambda r^{-\beta}({\bf x})}(r({\bf x})|\nabla v({\bf
x})|^{2}+\beta^{2}\lambda^{2}r^{-2\beta-1}({\bf x})|v({\bf
x})|^{2})\,d\sigma({\bf x})$ (2.4)

for some constant $C$ depending only on $d$ and $\gamma$. Since $v=0$ on
$\big{(}B({\bf x}_{0},R_{3})\setminus\overline{B({\bf
x}_{0},R_{1})}\big{)}\setminus\Omega$ and since $\Omega\Subset B({\bf
x}_{0},R_{3})\setminus\overline{B({\bf x}_{0},R_{1})}$, allowing $C$ to depend
on $\beta$ and $\Omega$, we deduce the Carleman estimate (2.3) from (2.4). ∎

An alternative way to obtain (2.3) is to apply the Carleman estimate in [25,
Chapter 4, §1, Lemma 3] for general parabolic operators. The arguments to
obtain (2.3) using [25, Chapter 4, §1, Lemma 3] are similar to that in [28,
Section 3] with the Laplacian replaced by the operator ${\rm
div}(A\nabla\cdot)$.

######  Remark 2.1.

We specially draw the reader’s attention to different forms of Carleman
estimates for all three main kinds of differential operators (elliptic,
parabolic and hyperbolic) and their applications in inverse problems and
computational mathematics [3, 4, 22]. It is worth mentioning that some
Carleman estimates hold true for all functions $v$ satisfying
$v|_{\partial\Omega}=0$ and $\partial_{\nu}v|_{\Gamma}=0$ where $\Gamma$ is a
part of $\partial\Omega$, see e.g., [23, 32], which can be used to solve
quasilinear elliptic PDEs with the boundary data partly given.

###  2.2 An inverse scattering problem

As mentioned in Section 1, Problem 1.1 arises from nonlinear inverse problems.
We present here an important example in the context of inverse scattering
problems in the frequency domain. Let $c:\mathbb{R}^{d}\to[1,\infty)$ be the
spatially distributed dielectric constant of the medium and
$[\underline{k},\overline{k}]$ be an interval of wavenumbers with
$\underline{k}>0$. Since the dielectric constant of the air is 1, we set
$c({\bf x})=1$ for all ${\bf x}\in\mathbb{R}^{d}\setminus\Omega.$ For each
$k\in[\underline{k},\overline{k}],$ let $w({\bf x},k)$, ${\bf
x}\in\mathbb{R}^{d}$, be the wave field generated by a point source at ${\bf
x}_{0}\in\mathbb{R}^{d}\setminus\Omega$ with wavenumber $k$. The function $w$
is governed by the Helmholtz equation and the Sommerfeld radiation condition

$\left\\{\begin{array}[]{ll}\Delta w({\bf x},k)+k^{2}c({\bf x})w({\bf
x},k)=-\delta({\bf x}-{\bf x}_{0})&{\bf x}\in\mathbb{R}^{d},\\\
(\frac{\partial}{\partial{|{\bf x}|}}-{\rm i}k)w({\bf x},k)=o(|{\bf
x}|^{(1-d)/2})&|{\bf x}|\to\infty.\end{array}\right.$ (2.5)

The inverse scattering problem is formulated as the problem of computing the
function $c({\bf x})$, ${\bf x}\in\Omega$, from the measurements of

$f_{1}({\bf x},k)=w({\bf x},k)\quad\mbox{and}\quad f_{2}({\bf
x},k)=\partial_{\nu}w({\bf x},k)$ (2.6)

for all ${\bf x}\in\partial\Omega$, $k\in[\underline{k},\overline{k}].$ The
knowledge of the function $c$ partly provides the internal structure of the
domain $\Omega$. In other words, solving the inverse scattering problem allows
us to examine a domain from external measurements, which has applications in
security, sonar imaging, geographical exploration, medical imaging, near-field
optical microscopy, nano-optics, see, e.g., [6] and references therein for
more details.

######  Remark 2.2.

In practice, measuring the wave field is expensive, and one does not measure
both $f_{1}$ and $f_{2}$. Rather, one measures the wave $w({\bf x},k)$ for all
${\bf x}$ on a surface $S$ far away from $\Omega.$ We call the information of
$w({\bf x},k)$ on $S\times[\underline{k},\overline{k}]$ the far field. Then,
we employ a technique, called “data propagation” to compute the near field on
$\partial\Omega\times[\underline{k},\overline{k}]$. This data propagation is
based on the angular spectrum representation, see, e.g., [33, Chapter 2] and
[29, Section 4.2.1] for details in implementation. This technique allows us to
compute the wave field in the region between the measurement surface $S$ and
$\partial\Omega$ where $c=1$. The knowledges of $f_{1}$ and $f_{2}$ follow.

There have been many important methods to solve inverse scattering problems in
the literature. Each method has its own advantages and disadvantages. A common
drawback of the widely-used method based on optimization to solve inverse
scattering problems is the need of a good initial guess of the true solution
$c$. We recall from [27] a method to solve the above inverse scattering
problem in which such a need is relaxed. Denote by

$z({\bf x},k)=\frac{1}{k^{2}}\log\Big{(}\frac{w({\bf x},k)}{w_{0}({\bf
x},k)}\Big{)}\quad\text{ for all }{\bf
x}\in\Omega,k\in[\underline{k},\overline{k}],$

where $w_{0}({\bf x},k)=\frac{e^{{\rm i}k|{\bf x}-{\bf x}_{0}|}}{4\pi|{\bf
x}-{\bf x}_{0}|}.$ Then, $z$ satisfies

$\Delta z({\bf x},k)+k^{2}|\nabla z({\bf x},k)|^{2}+\frac{2\nabla z({\bf
x},k)\cdot\nabla w_{0}({\bf x},k)}{w_{0}({\bf x},k)}=-c({\bf x})+1\quad\text{
for all }{\bf x}\in\Omega,k\in[\underline{k},\overline{k}].$

Let $\\{\Psi_{n}\\}_{n\geq 1}$ be the orthonormal basis of
$L^{2}(\underline{k},\overline{k})$ introduced in [19] and define

$z_{n}({\bf x})=\int_{\underline{k}}^{\overline{k}}z({\bf
x},k)\Psi_{n}(k)\,dk\quad\text{ for }n\geq 1,{\bf x}\in\Omega.$ (2.7)

We approximate

$v({\bf x},k)=\sum_{i=1}^{\infty}z_{i}({\bf
x})\Psi_{i}(k)\approx\sum_{i=1}^{N}z_{i}({\bf x})\Psi_{i}(k),$

for a suitable cut-off number $N\in\mathbb{N}$. Then, the vector
$Z_{N}=(z_{1},z_{2},\dots,z_{N})$ “approximately” satisfies the system

$\sum_{i=1}^{N}s_{li}\Delta z_{i}(\mathbf{x})+\sum_{i,j=1}^{N}a_{lij}\nabla
z_{i}(\mathbf{x})\cdot\nabla
z_{j}(\mathbf{x})+\sum_{i=1}^{N}B_{li}(\mathbf{x})\cdot\nabla
z_{i}(\mathbf{x})=0$ (2.8)

where

$\left\\{\begin{array}[]{l}s_{li}=\displaystyle\int_{\underline{k}}^{\overline{k}}\Psi_{i}^{\prime}(k)\Psi_{l}(k)\,dk,\\\
a_{lij}=\displaystyle
2\int_{\underline{k}}^{\overline{k}}\Big{(}k^{2}\Psi_{i}(k)\Psi_{j}^{\prime}(k)+k\Psi_{i}(k)\Psi_{j}(k)\Big{)}\Psi_{l}(k)\,dk,\\\
B_{li}(\mathbf{x})=\displaystyle
2\int_{\underline{k}}^{\overline{k}}\Big{(}\Psi_{i}^{\prime}(k)\frac{\nabla
w_{0}(\mathbf{x},k)}{w_{0}(\mathbf{x},k)}+\Psi_{i}(k)\partial_{k}\frac{\nabla
w_{0}(\mathbf{x},k)}{w_{0}(\mathbf{x},k)}\Big{)}\Psi_{l}(k)\,dk\end{array}\right.$

for all $i,j,l\in\\{1,\dots,N\\}$ and $\mathbf{x}\in\Omega$, see [27, Section
6] for details. The Dirichlet and Neumann boundary conditions for $Z_{N}$ can
be computed by the knowledges of $f_{1}$, $f_{2}$, and (2.7). Solving the
system of quasilinear elliptic equations (2.8) with the provided Dirichlet and
Neumann data is basically a goal of Problem 1.1. Doing so is the key step to
compute $c$. See [27] for convexification method to compute $Z_{N}$ and the
procedure to obtain $c$ from the knowledge of $Z_{N}.$

We refer the reader to [14, 13, 12] for another setup of the inverse
scattering problem and numerical results from experimental data. Since inverse
problems are out of the scope of this paper, we only mention it here to
explain the significance of Problem 1.1. In future works, we will use our
solver for Problem 1.1 to solve inverse problems.

##  3 The iteration and linearization approach for Problem 1.1

Our approach to solve (1.1) is based on linearization and iteration. Assume
that the solution $u^{*}$ to (1.1) is in the space $H^{p}(\Omega)$ for some
$p>\lceil d/2\rceil+2$ where $\lceil d/2\rceil$ is the smallest integer that
is greater than $d/2.$ Define the set of admissible solutions

$W=\big{\\{}\varphi\in
H^{p}(\Omega)\;:\,\varphi|_{\partial\Omega}=f,\partial_{\nu}\varphi|_{\partial\Omega}=g\big{\\}}.$
(3.1)

Then, the assumption in Problem 1.1 implies that $W\neq\emptyset$, and
$u^{*}\in W$. We now construct a sequence $\\{u_{n}\\}_{n\geq 0}$ that
converges to the solution $u^{*}$ to (1.1). Take a function $u_{0}\in W.$
Assume by induction that we have the knowledge of $u_{n}$ for some $n\geq 0$.
We find $u_{n+1}$ as follows. Assume that $u^{*}=u_{n}+h$ for some $h\in
H_{0}$ where

$H_{0}=\\{\varphi\in
H^{p}(\Omega)\,:\,\varphi|_{\partial\Omega}=0,\partial_{\nu}\varphi|_{\partial\Omega}=0\\}.$

Plugging $u^{*}=u_{n}+h$ into (1.1), we have

$-{\rm div}(A({\bf x})\nabla(u_{n}+h)({\bf x}))+F({\bf x},u_{n}({\bf
x})+h({\bf x}),\nabla u_{n}({\bf x})+\nabla h({\bf x}))=0\quad\mbox{for all
}{\bf x}\in\Omega.$ (3.2)

Heuristically, we assume at this moment that $h$ is small, that is,
$\|h\|_{C^{1}(\Omega)}\ll 1$. This assumption plays the role in establishing a
numerical scheme to find $u_{n+1}$ while it does not play a key role in the
proof of the convergence theorem. By Taylor’s expansion, we approximate (3.2)
as

$-{\rm div}(A\nabla u_{n})-{\rm div}(A\nabla h)+F({\bf x},u_{n}({\bf
x}),\nabla u_{n}({\bf x}))+D\mathcal{F}(u_{n})h=0\quad\mbox{for all }{\bf
x}\in\Omega$ (3.3)

where

$D\mathcal{F}(v)h=F_{s}({\bf x},v({\bf x}),\nabla v({\bf x}))h+\nabla_{{\bf
p}}F({\bf x},v({\bf x}),\nabla v({\bf x}))\cdot\nabla h({\bf x})$

for all $v\in H^{p}(\Omega)$. Here, $F_{s}$ and $\nabla_{{\bf p}}F$ are the
partial derivative of $F$ with respect to its second variable and its gradient
vector with respect to the third variable, respectively.

The next step is to compute a function $h\in H_{0}$ satisfying (3.3). Since
there is no guarantee for the existence of such a function $h$, we only
compute a “best fit” function $h$ by the Carleman-based quasi-reversibility
method described below. For each $v\in H^{p}(\Omega)$, $\lambda>1$, $\beta>0$
and $\eta>0$, define the functional
$J_{v}^{\lambda,\beta,\eta}:H_{0}\to\mathbb{R}$ as

$J_{v}^{\lambda,\beta,\eta}(\varphi)=\int_{\Omega}e^{2\lambda\mu_{\beta}({\bf
x})}\Big{|}-{\rm div}(A\nabla\varphi)-{\rm div}(A\nabla v)+F({\bf x},v({\bf
x}),\nabla v({\bf x}))+D\mathcal{F}(v)\varphi\Big{|}^{2}\,d{\bf
x}+\eta\|v+\varphi\|^{2}_{H^{p}(\Omega)}.$

Here, the function $\mu_{\beta}$ is defined in (2.1) and
$\eta\|v+\varphi\|^{2}_{H^{p}(\Omega)}$ is a regularization term.

For each $n\geq 0$, we minimize $J_{u_{n}}^{\lambda,\beta,\eta}(\varphi)$ on
$H_{0}$. The unique minimizer is the desired function $h.$ We then set

$u_{n+1}=u_{n}+h.$ (3.4)

The construction of the sequence $\\{u_{n}\\}_{n\geq 0}$ above is summarized
in Algorithm 1. We will prove that the sequence $\\{u_{n}\\}_{n\geq 0}$
converges to $u^{*}$ in Section 4 as $n\to\infty$ and $\eta\to 0$. The
presence of the Carleman weight $e^{2\lambda\mu_{\beta}({\bf x})}$ is a key
point for us to prove this convergence result.

Algorithm 1 The procedure to compute the numerical solution to (1.1)

1: Choose a threshold number $0<\kappa_{0}\ll 1.$ Choose an arbitrary initial
solution $u_{0}\in W.$

2: Set $n=0$.

3: Solve the linear equation (3.3) for a function $h\in H_{0}$ by minimizing
$J_{u_{n}}^{\lambda,\beta,\eta}(\varphi)$ on $H_{0}$. Set $u_{n+1}=u_{n}+h$.

4: if $\|u_{n+1}-u_{n}\|_{L^{2}(\Omega)}>\kappa_{0}$ then

5: Replace $n$ by $n+1.$

6: Go back to Step 3.

7: else

8: Set the computed solution $u_{\rm comp}=u_{n+1}.$

9: end if

##  4 The convergence analysis

In this section, we prove that the sequence $\\{u_{n}\\}_{n\geq 0}$ generated
by Algorithm 1 converges to the solution $u^{*}$ to (1.1). The following
result is the main theorem in this paper.

######  Theorem 4.1.

Assume that $\|F\|_{C^{2}(\Omega,\mathbb{R},\mathbb{R}^{d})}$ is a finite
number. Assume further that (1.1) has a unique solution $u^{*}\in W$. Let
$\\{u_{n}\\}_{n\geq 0}$ be the sequence generated by Algorithm 1, where
$u_{0}\in W$ is chosen arbitrarily. Then, there exist $\lambda_{0}>0$ and
$\theta\in(0,1)$ such that, for all $\lambda>\lambda_{0}$,

$\|u_{n+1}-u^{*}\|_{\lambda,\mu,\beta}^{2}\leq\theta^{n+1}\|u_{0}-u^{*}\|_{\lambda,\mu,\beta}^{2}+\eta\theta\frac{1-\theta^{n+1}}{1-\theta}\|u^{*}\|_{H^{p}(\Omega)}^{2}.$
(4.1)

Here,

$\|v\|_{\lambda,\mu,\beta}=\Big{[}\int_{\Omega}e^{2\lambda\mu_{\beta}({\bf
x})}\big{(}|v|^{2}+|\nabla v|^{2}\big{)}\,d{\bf
x}\Big{]}^{\frac{1}{2}}\quad\mbox{for all }v\in H^{1}(\Omega).$

###### Proof.

Fix $n\geq 0$. Let $h=u_{n+1}-u_{n}$. Due to Step 3 in Algorithm 1, $h$ is the
minimizer of $J_{u_{n}}^{\lambda,\beta,\eta}$. By the variational principle,
we have

$\int_{\Omega}e^{2\lambda\mu_{\beta}({\bf x})}\big{(}-{\rm div}(A\nabla
h)-{\rm div}(A\nabla u_{n})+F({\bf x},u_{n},\nabla
u_{n})+D\mathcal{F}(u_{n})h\big{)}\big{(}-{\rm
div}(A\nabla\varphi)+D\mathcal{F}(u_{n})\varphi\big{)}\,d{\bf x}\\\
+\eta\langle u_{n}+h,\varphi\rangle_{H^{p}(\Omega)}=0$ (4.2)

for all $\varphi\in H_{0}.$ Since $u_{n+1}=u_{n}+h$, we can rewrite (4.2) as

$\int_{\Omega}e^{2\lambda\mu_{\beta}({\bf x})}\big{(}-{\rm div}(A\nabla
u_{n+1})+F({\bf x},u_{n},\nabla
u_{n})+D\mathcal{F}(u_{n})(u_{n+1}-u_{n})\big{)}\big{(}-{\rm
div}(A\nabla\varphi)+D\mathcal{F}(u_{n})\varphi\big{)}\,d{\bf x}\\\
+\eta\langle u_{n+1},\varphi\rangle_{H^{p}(\Omega)}=0$ (4.3)

for all $\varphi\in H_{0}.$ As $u^{*}$ is the solution to (1.1), we have

$\int_{\Omega}e^{2\lambda\mu_{\beta}({\bf x})}\big{(}-{\rm div}(A\nabla
u^{*})+F({\bf x},u^{*},\nabla u^{*})\big{)}\big{(}-{\rm
div}(A\nabla\varphi)+D\mathcal{F}(u_{n})\varphi\big{)}\,d{\bf x}=0$ (4.4)

for all $\varphi\in H_{0}.$ It follows from (4.3) and (4.4) that for all
$\varphi\in H_{0}$

$\int_{\Omega}e^{2\lambda\mu_{\beta}({\bf x})}\big{(}-{\rm
div}(A\nabla(u_{n+1}-u^{*}))+F({\bf x},u_{n},\nabla u_{n})-F({\bf
x},u^{*},\nabla u^{*})+D\mathcal{F}(u_{n})(u_{n+1}-u_{n})\big{)}\\\
\big{(}-{\rm div}(A\nabla\varphi)+D\mathcal{F}(u_{n})\varphi\big{)}\,d{\bf
x}+\eta\langle u_{n+1},\varphi\rangle_{H^{p}(\Omega)}=0.$ (4.5)

Take $\varphi=u_{n+1}-u^{*}\in H_{0}$. It follows from (4.5) that

$\int_{\Omega}e^{2\lambda\mu_{\beta}({\bf x})}\Big{[}|{\rm
div}(A\nabla\varphi)|^{2}-D\mathcal{F}(u_{n})\varphi\,{\rm
div}(A\nabla\varphi)-\big{(}F({\bf x},u_{n},\nabla u_{n})-F({\bf
x},u^{*},\nabla u^{*})\big{)}{\rm div}(A\nabla\varphi)\\\ +\big{(}F({\bf
x},u_{n},\nabla u_{n})-F({\bf x},u^{*},\nabla
u^{*})\big{)}D\mathcal{F}(u_{n})\varphi-{\rm
div}(A\nabla\varphi)D\mathcal{F}(u_{n})(u_{n+1}-u_{n})\\\
+D\mathcal{F}(u_{n})(u_{n+1}-u_{n})D\mathcal{F}(u_{n})\varphi\Big{]}\,d{\bf
x}=-\eta\langle u_{n+1},\varphi\rangle_{H^{p}(\Omega)}.$ (4.6)

As $\|F\|_{C^{2}(\Omega,\mathbb{R},\mathbb{R}^{d})}=C<+\infty$, we can
estimate

$\displaystyle|D\mathcal{F}(u_{n})\varphi|\leq C(|\varphi|+|\nabla\varphi|),$
$\displaystyle|F({\bf x},u_{n},\nabla u_{n})-F({\bf x},u^{*},\nabla
u^{*})|\leq C(|u_{n}-u^{*}|+|\nabla(u_{n}-u^{*})|),$
$\displaystyle|D\mathcal{F}(u_{n})(u_{n+1}-u_{n})|=|D\mathcal{F}(u_{n})(u_{n+1}-u^{*}+u^{*}-u_{n})|$
$\displaystyle\hskip 56.9055pt\leq
C(|u_{n+1}-u^{*}|+|\nabla(u_{n+1}-u^{*})|)+C(|u_{n}-u^{*}|+|\nabla(u_{n}-u^{*})|),$
and $\displaystyle-\eta\langle
u_{n+1},\varphi\rangle_{H^{p}(\Omega)}=-\eta\langle\varphi+u^{*},\varphi\rangle_{H^{p}(\Omega)}=-\eta\|\varphi\|_{H^{p}(\Omega)}^{2}-\eta\langle
u^{*},\varphi\rangle_{H^{p}(\Omega)}\leq\frac{\eta}{2}\|u^{*}\|_{H^{p}(\Omega)}^{2}.$

These estimates, together with the inequality $|ab|\leq 4a^{2}+b^{2}/8$, imply

$\int_{\Omega}e^{2\lambda\mu_{\beta}({\bf x})}|{\rm
div}(A\nabla\varphi)|^{2}\,d{\bf x}\leq
C\Big{[}\int_{\Omega}e^{2\lambda\mu_{\beta}({\bf
x})}(|\varphi|^{2}+|\nabla\varphi|^{2})\,d{\bf x}\\\
+\int_{\Omega}e^{2\lambda\mu_{\beta}({\bf
x})}(|u_{n}-u^{*}|^{2}+|\nabla(u_{n}-u^{*})|^{2})\,d{\bf
x}\Big{]}+\frac{\eta}{2}\|u^{*}\|_{H^{p}(\Omega)}^{2}.$ (4.7)

Applying the Carleman estimate (2.3) for the function $\varphi$, we have

$\int_{\Omega}e^{2\lambda\mu_{\beta}({\bf x})}|{\rm
div}(A\nabla\varphi)|^{2}\,d{\bf x}\geq
C\lambda\int_{\Omega}e^{2\lambda\mu_{\beta}({\bf x})}|\nabla\varphi({\bf
x})|^{2}\,d{\bf x}+C\lambda^{3}\int_{\Omega}e^{2\lambda\mu_{\beta}({\bf
x})}|\varphi({\bf x})|^{2}\,d{\bf x}.$ (4.8)

Combining (4.7) and (4.8), we have

$\lambda\int_{\Omega}e^{2\lambda\mu_{\beta}({\bf x})}|\nabla\varphi({\bf
x})|^{2}\,d{\bf x}+\lambda^{3}\int_{\Omega}e^{2\lambda\mu_{\beta}({\bf
x})}|\varphi({\bf x})|^{2}\,d{\bf x}\leq
C\Big{[}\int_{\Omega}e^{2\lambda\mu_{\beta}({\bf
x})}(|\varphi|^{2}+|\nabla\varphi|^{2})\,d{\bf x}\\\
+\int_{\Omega}e^{2\lambda\mu_{\beta}({\bf
x})}(|u_{n}-u^{*}|^{2}+|\nabla(u_{n}-u^{*})|^{2})\,d{\bf
x}\Big{]}+\frac{\eta}{2}\|u^{*}\|_{H^{p}(\Omega)}^{2}.$ (4.9)

Letting $\lambda$ be sufficiently large, we can simplify (4.9) as

$\lambda\int_{\Omega}e^{2\lambda\mu_{\beta}({\bf
x})}(|\varphi|^{2}+|\nabla\varphi|^{2})\,d{\bf x}\leq
C\int_{\Omega}e^{2\lambda\mu_{\beta}({\bf
x})}(|u_{n}-u^{*}|^{2}+|\nabla(u_{n}-u^{*})|^{2})\,d{\bf
x}+\frac{\eta}{2}\|u^{*}\|_{H^{p}(\Omega)}^{2}.$ (4.10)

Recall that $\varphi=u_{n+1}-u^{*}$. We get from (4.10) that

$\|u_{n+1}-u^{*}\|_{\lambda,\mu,\beta}^{2}\leq\frac{C}{\lambda}\|u_{n}-u^{*}\|_{\lambda,\mu,\beta}^{2}+\frac{\eta}{2\lambda}\|u^{*}\|_{H^{p}(\Omega)}^{2}.$
(4.11)

Applying (4.11) for $n-1$ and denoting $\theta=C/\lambda\in(0,1)$, we have

$\displaystyle\|u_{n+1}-u^{*}\|_{\lambda,\mu,\beta}^{2}$
$\displaystyle\leq\theta\Big{[}\theta\|u_{n-1}-u^{*}\|_{\lambda,\mu,\beta}^{2}+\eta\theta\|u^{*}\|_{H^{p}(\Omega)}^{2}\Big{]}+\eta\theta\|u^{*}\|_{H^{p}(\Omega)}^{2}$
$\displaystyle=\theta^{2}\|u_{n-1}-u^{*}\|_{\lambda,\mu,\beta}^{2}+\eta\theta(1+\theta)\|u^{*}\|_{H^{p}(\Omega)}^{2}.$

By induction, we have

$\|u_{n+1}-u^{*}\|_{\lambda,\mu,\beta}^{2}\leq\theta^{n+1}\|u_{0}-u^{*}\|_{\lambda,\mu,\beta}^{2}+\eta\theta\sum_{i=0}^{n}\theta^{n}\|u^{*}\|_{H^{p}(\Omega)}^{2},$
(4.12)

which implies (4.1). The proof is complete. ∎

######  Remark 4.1 (Removing the boundedness condition of $F$ in $C^{2}$ in
Theorem 4.1).

In the case when
$\|F\|_{C^{2}(\Omega\times\mathbb{R}\times\mathbb{R}^{d+1})}=\infty$, we need
to assume that we know in advance that the true solution $u^{*}$ to (1.1)
belongs to the ball $B_{M}$ of $C^{1}(\Omega)$ for some $M>0$. This assumption
does not weaken the result since $M$ can be arbitrarily large. Define the cut-
off function $\chi\in C^{\infty}(\Omega\times\mathbb{R}\times\mathbb{R}^{d})$
as

$\chi({\bf x},s,{\bf p})=\left\\{\begin{array}[]{ll}1&(s^{2}+|{\bf
p}|^{2})^{1/2}<M,\\\ 0&(s^{2}+|{\bf p}|^{2})^{1/2}>2M\end{array}\right.$

and set $\widetilde{F}=\chi F.$ Since $({u^{*}}^{2}+|\nabla
u^{*}|^{2})^{1/2}<M$, it is obvious that $u^{*}$ solves the problem

$\left\\{\begin{array}[]{ll}-{\rm div}(A({\bf x})\nabla u({\bf
x}))+\widetilde{F}({\bf x},u({\bf x}),\nabla u({\bf x}))=0&{\bf
x}\in\Omega,\\\ u({\bf x})=f({\bf x})&{\bf x}\in\partial\Omega,\\\
\partial_{\nu}u({\bf x})=g({\bf x})&{\bf
x}\in\partial\Omega.\end{array}\right.$ (4.13)

Now, we can apply Algorithm 1 for (4.13) to compute $u^{*}$.

######  Remark 4.2.

Theorem 4.1 and estimate (4.1) rigorously guarantee that each sequence
generated by Algorithm 1 converges to $u^{*}$ at the exponential rate. This
fact is numerically confirmed by our numerical results in Section 5.

######  Remark 4.3.

As seen in the proof of Theorem 4.1, the efficiency of Algorithm 1 is
guaranteed by Carleman estimate (2.3). Therefore, we call the proposed method
described in Algorithm 1 the second generation of Carleman-based numerical
methods. The first generation of Carleman-based numerical method was developed
in [20], which is called the convexification. See [1, 23, 26, 27] for some
following up results. Like the convexification method, Algorithm 1 can be used
to compute solutions to nonlinear PDEs without requesting an initial good
guess. The advantage of our new method is the fast convergence rate, see
Remark 4.2.

##  5 Numerical study

In this section, we present some numerical results obtained by Algorithm 1.
For simplicity, we set $d=2$ and $\Omega=(-1,1)^{2}$. On $\Omega,$ we arrange
an $N\times N$ grid points

$\mathcal{G}=\\{(x_{i}=-1+(i-1)\delta,y_{j}=-1+(j-1)\delta):1\leq i,j\leq
N\\}$

where $\delta=2/(N-1).$ In our numerical scripts, $N=80$.

In Step 1 of Algorithm 1, we choose a function $u_{0}\in W$. It is natural to
find a function $u_{0}$ satisfying the equation obtained by removing from
(1.1) the nonlinearity $F({\bf x},u,\nabla u)$. We apply the Carleman-based
quasi-reversibility method to do so. That means, $u_{0}$ is the minimizer of

$J_{0}^{\lambda,\beta,\eta}(\varphi)=\int_{\Omega}e^{\lambda\mu_{\beta}({\bf
x})}|{\rm div}(A\nabla\varphi)|^{2}d{\bf
x}+\eta\|\varphi\|^{2}_{H^{2}(\Omega)}$ (5.1)

subject to the boundary conditions in (5.7). To simplify the efforts in
implementation, the norm in the regularization term is the
$H^{2}(\Omega)$-norm rather than the $H^{p}(\Omega)$-norm. This change does
not affect the performance of Algorithm 1. Algorithm 1 still provides
satisfactory solutions to (1.1).

######  Remark 5.1.

We employ the Carleman-based quasi-reversibility method to find $u_{0}$ for
the consistency to Step 3 of Algorithm 1. Since $u_{0}\in W$ can be chosen
arbitrarily, one can use the quasi-reversibility method without the presence
of the Carleman weight function $e^{\lambda\mu_{\beta}({\bf x})}$. In our
computation for all numerical examples below, $\lambda=4$, $\beta=10$ and
${\bf x}_{0}=(-4,0).$ The regularized parameter is $\eta=10^{-4}$. The
threshold number $\kappa_{0}=10^{-6}.$

We minimize $J_{0}^{\lambda,\beta,\eta}$ on $H_{0}$ by the least square MATLAB
command “lsqlin”. The implementation for the quasi-reversibility method to
minimize more general functionals than $J_{0}^{\lambda,\beta,\eta}$ was
described in [26, §5.3] and in [31, §5]. We do not repeat this process here.
In Step 3 of Algorithm 1, given $u_{n}\in W,$ we minimize the functional
$J^{\lambda,\beta,\eta}_{u_{n}}$ on $H_{0}$. Again, we refer the reader to
[26, §5.3] and in [31, §5] for details in implementation. The scripts for
other steps of Algorithm 1 can be written easily.

###  5.1 Quasilinear elliptic equations

The convexification method, first introduced in [20], was used to numerically
solve quasilinear elliptic equations in [1, 23, 27]. Our approach here is of
course different based on iteration and linearization. In particular, Step 3
in Algorithm 1 is very efficient as we only need to solve the linear equation
(3.3) as opposed to solving directly a nonlinear quasilinear elliptic
equation. In this subsection, we present two (2) numerical tests. In both
tests, we choose the matrix $A$ to be

$A=\left[\begin{array}[]{ll}2&1\\\ 1&2\end{array}\right].$

That means, ${\rm div}(A\nabla u)=2u_{xx}+2u_{xy}+2u_{yy}.$

In test 1, we solve (1.1) when

$F({\bf x},s,{\bf p})=s+|{\bf
p}|-\big{(}-x^{2}+2y^{2}+\sqrt{4x^{2}+16y^{2}}-4\big{)}$ (5.2)

for all ${\bf x}=(x,y)\in\Omega$, $s\in\mathbb{R}$ and ${\bf
p}\in\mathbb{R}^{2}$. The boundary conditions are given by

$u({\bf x})=-x^{2}+2y^{2},\quad\partial_{\nu}u({\bf
x})=\langle-2x,4y\rangle\cdot\nu$ (5.3)

for all ${\bf x}=(x,y)\in\partial\Omega.$ The true solution of (1.1) is the
function $u_{\rm true}(x,y)=-x^{2}+2y^{2}$.

(a) The function $u^{*}$

(b)  The relative error $\frac{|u^{*}({\bf x})-u_{\rm comp}({\bf
x})|}{\|u_{\rm true}\|_{L^{\infty}(\Omega)}}$

(c)  $\|u_{n}-u_{n-1}\|_{L^{2}(\Omega)}.$ The horizontal axis is the number of
iteration $n.$

Figure 1:  Test 1. The numerical solution to (1.1) where $F$ is given in (5.2)
and the boundary data are given in (5.3).

It is evident from Figure 1 that Algorithm 1 provides out of expectation
solution for test 1. The relative error $\|u_{\rm true}-u_{\rm
comp}\|_{L^{\infty}(\Omega)}/\|u_{\rm true}\|_{L^{\infty}(\Omega)}=1.23\times
10^{-5}$. One can see from Figure 1c that Algorithm 1 converges at the third
iteration.

In test 2, we solve (1.1) when

$F({\bf x},s,{\bf p})=|{\bf
p}|-\Big{[}\sqrt{\Big{[}\frac{\pi}{2}\cos\big{(}\frac{\pi}{2}(x+y)\big{)}+e^{x}\Big{]}^{2}+\frac{\pi^{2}}{4}\cos^{2}\big{(}\frac{\pi}{2}(x+y)\big{)}}\\\
+\frac{3\pi^{2}}{2}\sin\big{(}\frac{\pi}{2}(x+y)\big{)}-2e^{x}\Big{]}$ (5.4)

for all ${\bf x}=(x,y)\in\Omega$, $s\in\mathbb{R}$ and ${\bf
p}\in\mathbb{R}^{2}$. The boundary conditions are given by

$u({\bf
x})=\sin\big{(}\frac{\pi}{2}(x+y)\big{)}+e^{x},\quad\partial_{\nu}u({\bf
x})=\big{\langle}\frac{\pi}{2}\cos\big{(}\frac{\pi}{2}(x+y)\big{)}+e^{x},\cos\big{(}\frac{\pi}{2}(x+y)\big{)}\big{\rangle}\cdot\nu$
(5.5)

for all ${\bf x}=(x,y)\in\partial\Omega.$ The true solution of (1.1) is the
function $u_{\rm true}(x,y)=\sin\big{(}\frac{\pi}{2}(x+y)\big{)}+e^{x}$.

(a) The function $u^{*}$

(b)  The relative error $\frac{|u^{*}({\bf x})-u_{\rm comp}({\bf
x})|}{\|u_{\rm true}\|_{L^{\infty}(\Omega)}}$

(c)  $\|u_{n}-u_{n-1}\|_{L^{2}(\Omega)}$. The horizontal axis is the number of
iteration $n$

Figure 2:  Test 2. The numerical solution to (1.1) where $F$ is given in (5.4)
and the boundary data are given in (5.5).

It is evident from Figure 2 that Algorithm 1 provides out of expectation
solution for test 2. The relative error $\|u_{\rm true}-u_{\rm
comp}\|_{L^{\infty}(\Omega)}/\|u_{\rm true}\|_{L^{\infty}(\Omega)}=4.19\times
10^{-5}$. One can see from Figure 2c that Algorithm 1 converges at the fifth
iteration.

###  5.2 Application to first-order Hamilton-Jacobi equations

Our aim in this subsection is to solve numerically

$F({\bf x},u,\nabla u)=0\quad\mbox{for all }{\bf x}\in\Omega.$ (5.6)

with the boundary conditions

$u|_{\partial\Omega}=f\quad\mbox{and}\quad\partial_{\nu}u|_{\partial\Omega}=g.$
(5.7)

Basically, we use the vanishing viscosity process to approximate solutions to
(5.6). For $\epsilon>0$, consider

$-\epsilon\Delta u+F({\bf x},u,\nabla u)=0\quad\mbox{for all }{\bf
x}\in\Omega$ (5.8)

with boundary conditions (5.7). Again, (5.8)–(5.7) is an over-determined
boundary value problem. For $\epsilon>0$ sufficiently small, assume that
(5.8)–(5.7) has a solution $u^{\epsilon}\in W$. Then, $u^{\epsilon}$
approximates $u$, solution to (5.6)–(5.7), quite well under some appropriate
assumptions on $F$. In our numerical tests, we choose
$\epsilon=\epsilon_{0}=10^{-3}$.

In this part, we provide six (6) numerical tests, in which we compute the
viscosity solution to some Hamilton-Jacobi equations of the form (5.6) given
Cauchy boundary data. That means, by applying Algorithm 1, we numerically find
a function $u_{\rm comp}$ satisfying (5.8)-(5.7) when $F$, $f$ and $g$ are
given. The verification that $u^{*}$ is the correct viscosity solution can be
done in a similar manner as that in [39, 23].

Test 1. In this test, we solve (5.8)-(5.7) when

$F({\bf x},s,{\bf p})=s+|{\bf p}|+|x|-1.$ (5.9)

for all ${\bf x}=(x,y)\in\Omega,s\in\mathbb{R},{\bf p}\in\mathbb{R}^{2}$. The
boundary conditions are given by

$u({\bf x})=-|x|,\quad\partial_{\nu}u({\bf
x})=\langle-\mbox{sign}(x),0\rangle\cdot\nu$ (5.10)

for all ${\bf x}=(x,y)\in\partial\Omega$. In this case, the true solution is
$u^{*}=-|x|.$ The numerical result is displayed in Figure 3.

(a) The true solution $u^{*}$ to the Hamilton-Jacobi equation where the
Hamiltonian is given in (5.9).

(b) The initial solution $u_{0}$ computed by minimizing
$J_{0}^{\lambda,\beta,\eta}$ defined in (5.1)

(c) The computed solution $u_{\rm comp}$.

(d) $\|u_{n}-u_{n-1}\|_{L^{2}(\Omega)}$. The horizontal axis is the number of
iteration $n$.

(e) The true and computed solutions on the line $\\{(x=0,y)\in\Omega\\}$

(f)  The relative error $\frac{|u^{*}({\bf x})-u_{\rm comp}({\bf
x})|}{\|u_{\rm true}\|_{L^{\infty}(\Omega)}}.$

Figure 3:  Test 1. The true viscosity solution to (5.8)-(5.7) and the computed
one. The Hamiltonian and the boundary conditions are given in (5.9)-(5.10).

It is evident from Figure 3 that we successfully compute the solution $u_{\rm
comp}$. The procedure described in Algorithm 1 converges after four
iterations. The relative error $\frac{\|u^{*}-u_{\rm
comp}\|_{L^{\infty}(\Omega)}}{\|u^{*}\|_{L^{\infty}(\Omega)}}=5.33\%.$

Test 2. We find the viscosity solution to the eikonal equation. In this test,
we solve (5.8)-(5.7) when the Hamiltonian is

$F({\bf x},s,{\bf p})=|{\bf p}|^{2}-(1+(1+\mbox{sign}(x+y))^{2})$ (5.11)

for all ${\bf x}=(x,y)\in\Omega,s\in\mathbb{R},{\bf p}\in\mathbb{R}^{2}$. The
boundary conditions are given by

$u({\bf x})=-|x+y|-y,\quad\partial_{\nu}u({\bf
x})=(-\mbox{sign}(x+y),-\mbox{sign}(x+y)-1)\cdot\nu$ (5.12)

for all ${\bf x}=(x,y)\in\partial\Omega$. The true solution is $u^{*}({\bf
x})=-|x+y|-y$. The graphs of $u^{*}$ and $u_{\rm comp}$ are displayed in
Figure 4.

(a) The true solution $u^{*}$ to the Hamilton-Jacobi equation where the
Hamiltonian is given in (5.11).

(b) The initial solution $u_{0}$ computed by minimizing
$J_{0}^{\lambda,\beta,\eta}$ defined in (5.1)

(c) The computed solution $u_{\rm comp}$.

(d) $\|u_{n}-u_{n-1}\|_{L^{2}(\Omega)}$. The horizontal axis is the number of
iteration $n$.

(e) The true and computed solutions on the line $\\{(x=0.5,y)\in\Omega\\}$

(f)  The relative error $\frac{|u^{*}({\bf x})-u_{\rm comp}({\bf
x})|}{\|u_{\rm true}\|_{L^{\infty}(\Omega)}}.$

Figure 4:  Test 2. The true viscosity solution to (5.8)-(5.7) and the computed
one. The Hamiltonian and the boundary conditions are given in (5.11)-(5.12).

The relative error $\frac{\|u^{*}-u_{\rm
comp}\|_{L^{\infty}(\Omega)}}{\|u^{*}\|_{L^{\infty}(\Omega)}}=3.6\%.$

Test 3. We test the case when the Hamiltonian is not convex with respect its
third variable. The Hamiltonian in this test is given by

$F({\bf x},s,{\bf
p})=20s+|p_{1}|-|p_{2}|-\Big{[}20(-|x+0.5|+e^{\cos(2\pi(x^{2}+y^{2}))})\\\
+|\mbox{sign}(x+0.5)+4\pi
x\sin(2\pi(x^{2}+y^{2}))e^{\cos(2\pi(x^{2}+y^{2}))}|\\\ -|4\pi
y\sin(2\pi(x^{2}+y^{2}))e^{\cos(2\pi(x^{2}+y^{2}))}|\Big{]}$ (5.13)

for all ${\bf x}=(x,y)\in\Omega,s\in\mathbb{R},{\bf
p}=(p_{1},p_{2})\in\mathbb{R}^{2}$. The boundary conditions are given by

$u({\bf x})=-|x+0.5|+e^{\cos(2\pi(x^{2}+y^{2}))}\quad\mbox{for all }{\bf
x}=(x,y)\in\partial\Omega$ (5.14)

and

$\partial_{\nu}u({\bf x})=\big{(}-\mbox{sign}(x+0.5)-4\pi
x\sin(2\pi(x^{2}+y^{2}))e^{\cos(2\pi(x^{2}+y^{2}))},\\\ -4\pi
y\sin(2\pi(x^{2}+y^{2}))e^{\cos(2\pi(x^{2}+y^{2}))}\big{)}\cdot\nu$ (5.15)

for all ${\bf x}=(x,y)\in\partial\Omega$. The true solution is $u^{*}({\bf
x})=-|x+0.5|+e^{\cos(2\pi(x^{2}+y^{2}))}$. The graphs of $u^{*}$ and $u_{\rm
comp}$ are displayed in Figure 5.

(a) The true solution $u^{*}$ to the Hamilton-Jacobi equation where the
Hamiltonian is given in (5.13).

(b) The initial solution $u_{0}$ computed by minimizing
$J_{0}^{\lambda,\beta,\eta}$ defined in (5.1)

(c) The computed solution $u_{\rm comp}$.

(d) $\|u_{n}-u_{n-1}\|_{L^{2}(\Omega)}$. The horizontal axis is the number of
iteration $n$.

(e) The true and computed solutions on the line $\\{(x=0,y)\in\Omega\\}$

(f)  The relative error $\frac{|u^{*}({\bf x})-u_{\rm comp}({\bf
x})|}{\|u_{\rm true}\|_{L^{\infty}(\Omega)}}.$

Figure 5:  Test 3. The true viscosity solution to (5.8)-(5.7) and the computed
one. The Hamiltonian and the boundary conditions are given in (5.13)-(5.15).

The relative error $\frac{\|u^{*}-u_{\rm
comp}\|_{L^{\infty}(\Omega)}}{\|u^{*}\|_{L^{\infty}(\Omega)}}=10.65\%.$
Although the relative error in this test is large, the numerical result is
acceptable. In fact, we observe from Figure 5f that $\frac{\|u^{*}-u_{\rm
comp}\|_{L^{\infty}(\Omega)}}{\|u^{*}\|_{L^{\infty}(\Omega)}}$ is small almost
everywhere while it is large at only two small places near the left edge of
the domain.

Test 4. We test the case when the Hamiltonian is not increasing in the second
variable and is not convex in the third variable. The Hamiltonian in this test
is given by

$F({\bf x},s,{\bf p})=-40s+||{\bf
p}|-10|+40\Big{(}|x+y-0.5|+\sin\Big{(}\frac{x^{2}}{2}+y^{2}\Big{)}\Big{)}\\\
-\Big{|}\sqrt{\Big{[}\mbox{sign}(x+y-0.5)+x\cos\Big{(}\frac{x^{2}}{2}+y^{2}\Big{)}\Big{]}^{2}+\Big{[}\mbox{sign}(x+y-0.5)+2y\cos\Big{(}\frac{x^{2}}{2}+y^{2}\Big{)}\Big{]}^{2}}-10\Big{|}$
(5.16)

for all ${\bf x}=(x,y)\in\Omega,s\in\mathbb{R},{\bf
p}=(p_{1},p_{2})\in\mathbb{R}^{2}$. The boundary conditions are given by

$u({\bf x})=|x+y-0.5|+\sin\Big{(}\frac{x^{2}}{2}+y^{2}\Big{)}\quad\mbox{for
all }{\bf x}=(x,y)\in\partial\Omega$ (5.17)

and

$\partial_{\nu}u({\bf
x})=\Big{(}\mbox{sign}(x+y-0.5)+x\cos\Big{(}\frac{x^{2}}{2}+y^{2}\Big{)},\mbox{sign}(x+y-0.5)+2y\cos\Big{(}\frac{x^{2}}{2}+y^{2}\Big{)}\Big{)}\cdot\nu$
(5.18)

for all ${\bf x}=(x,y)\in\partial\Omega$. The true solution is $u^{*}({\bf
x})=|x+y-0.5|+\sin\Big{(}\frac{x^{2}}{2}+y^{2}\Big{)}$. The graphs of $u^{*}$
and $u_{\rm comp}$ are displayed in Figure 6.

(a) The true solution $u^{*}$ to the Hamilton-Jacobi equation where the
Hamiltonian is given in (5.16).

(b) The initial solution $u_{0}$ computed by minimizing
$J_{0}^{\lambda,\beta,\eta}$ defined in (5.1)

(c) The computed solution $u_{\rm comp}$.

(d) The relative error $\|u_{n}-u_{n-1}\|_{L^{2}(\Omega)}$. The horizontal
axis is the number of iteration $n$.

(e) The true and computed solutions on the line $\\{(x=0,y)\in\Omega\\}$

(f)  The relative error $\frac{|u^{*}({\bf x})-u_{\rm comp}({\bf
x})|}{\|u_{\rm true}\|_{L^{\infty}(\Omega)}}.$

Figure 6:  Test 4. The true viscosity solution to (5.8)-(5.7) and the computed
one. The Hamiltonian and the boundary conditions are given in (5.16)-(5.18).

The relative error $\frac{\|u^{*}-u_{\rm
comp}\|_{L^{\infty}(\Omega)}}{\|u^{*}\|_{L^{\infty}(\Omega)}}=0.95\%.$

Test 5. We solve a G-equation. The Hamiltonian in this test is given by

$F({\bf x},s,{\bf p})=5s+|{\bf
p}|-xp_{1}+\Big{[}5(|x-0.5|+|y|)-x\mbox{sign}(x-0.5)-\sqrt{2}\Big{]}$ (5.19)

for all ${\bf x}=(x,y)\in\Omega,s\in\mathbb{R},{\bf
p}=(p_{1},p_{2})\in\mathbb{R}^{2}$. The boundary conditions are given by

$u({\bf x})=-|x-0.5|-|y|\quad\mbox{for all }{\bf x}=(x,y)\in\partial\Omega$
(5.20)

and

$\partial_{\nu}u({\bf
x})=-\big{(}\mbox{sign}(x-0.5),\mbox{sign}(y)\big{)}\cdot\nu$ (5.21)

for all ${\bf x}=(x,y)\in\partial\Omega$. The true solution is $u^{*}({\bf
x})=-|x-0.5|-|y|$. The graphs of $u^{*}$ and $u_{\rm comp}$ are displayed in
Figure 5.

(a) The true solution $u^{*}$ to the Hamilton-Jacobi equation where the
Hamiltonian is given in (5.19).

(b) The initial solution $u_{0}$ computed by minimizing
$J_{0}^{\lambda,\beta,\eta}$ defined in (5.1)

(c) The computed solution $u_{\rm comp}$.

(d) The relative error $\|u_{n}-u_{n-1}\|_{L^{2}(\Omega)}$. The horizontal
axis is the number of iteration $n$.

(e) The true and computed solutions on the line $\\{(x=0,y)\in\Omega\\}$

(f)  The relative error $\frac{|u^{*}({\bf x})-u_{\rm comp}({\bf
x})|}{\|u_{\rm true}\|_{L^{\infty}(\Omega)}}.$

Figure 7:  Test 5. The true viscosity solution to (5.8)-(5.7) and the computed
one. The Hamiltonian and the boundary conditions are given in (5.19)-(5.21).

The relative error $\frac{\|u^{*}-u_{\rm
comp}\|_{L^{\infty}(\Omega)}}{\|u^{*}\|_{L^{\infty}(\Omega)}}=0.98\%.$

Test 6. The Hamiltonian in this test is given by

$F({\bf x},s,{\bf p})=20s+\min\\{|{\bf p}|,||{\bf
p}|-10|+6\\}-\Big{[}20(-|x|+\sin(\pi(x^{2}+y^{2})))\\\
+\min\Big{\\{}h(x,y),|h(x,y)-10|+6\Big{\\}}\Big{]}$ (5.22)

for all ${\bf x}=(x,y)\in\Omega,s\in\mathbb{R},{\bf
p}=(p_{1},p_{2})\in\mathbb{R}^{2}$ where

$h(x,y)=\sqrt{[-\mbox{sign}(x)+2\pi\cos(\pi(x^{2}+y^{2}))]^{2}+[2\pi\cos(\pi(x^{2}+y^{2})]^{2}}.$

The boundary conditions are given by

$u({\bf x})=-|x|+\sin(\pi(x^{2}+y^{2}))\quad\mbox{for all }{\bf
x}=(x,y)\in\partial\Omega$ (5.23)

and

$\partial_{\nu}u({\bf x})=\big{(}-\mbox{sign}(x)+2\pi
x\cos(\pi(x^{2}+y^{2})),2\pi y\cos(\pi(x^{2}+y^{2}))\big{)}\cdot\nu$ (5.24)

for all ${\bf x}=(x,y)\in\partial\Omega$. The true solution is $u^{*}({\bf
x})=-|x|+\sin(\pi(x^{2}+y^{2}))$. The graphs of $u^{*}$ and $u_{\rm comp}$ are
displayed in Figure 8.

(a) The true solution $u^{*}$ to the Hamilton-Jacobi equation where the
Hamiltonian is given in (5.22).

(b) The initial solution $u_{0}$ computed by minimizing
$J_{0}^{\lambda,\beta,\eta}$ defined in (5.1)

(c) The computed solution $u_{\rm comp}$.

(d) The relative error $\|u_{n}-u_{n-1}\|_{L^{2}(\Omega)}$. The horizontal
axis is the number of iteration $n$.

(e) The true and computed solutions on the line $\\{(x=0,y)\in\Omega\\}$

(f)  The relative error $\frac{|u^{*}({\bf x})-u_{\rm comp}({\bf
x})|}{\|u_{\rm true}\|_{L^{\infty}(\Omega)}}.$

Figure 8:  Test 6. The true viscosity solution to (5.8)-(5.7) and the computed
one. The Hamiltonian and the boundary conditions are given in (5.22)-(5.24).

The relative error $\frac{\|u^{*}-u_{\rm
comp}\|_{L^{\infty}(\Omega)}}{\|u^{*}\|_{L^{\infty}(\Omega)}}=4.8\%.$

######  Remark 5.2.

The $L^{\infty}$ relative errors in all tests above for first-order Hamilton-
Jacobi equations are compatible with
$\max\\{O(\sqrt{\eta}),O(\sqrt{\epsilon_{0}})\\}\simeq 3\%.$

##  6 Concluding remarks

We have developed a new globally convergent numerical method to solve over-
determined boundary value problems of quasilinear elliptic equations. The key
point of the method is to repeatedly solve the linearization of the given PDE
by the Carleman weighted quasi-reversibility method (Algorithm 1). As the
result, we obtain a sequence of functions converging to the solution thanks to
the Carleman estimate (Lemma 2.1). The strength of our new method includes (1)
the global convergence property and (2) the fast convergence rate, which is
described in Theorem 4.1. Some numerical results for quasilinear elliptic
equations and first-order Hamilton-Jacobi equations are presented to show the
effectiveness of our method.

## Acknowledgement

The works of TTL and LHN were supported by US Army Research Laboratory and US
Army Research Office grant W911NF-19-1-0044. The work of HT is supported in
part by NSF CAREER grant DMS-1843320 and a Simons fellowship.

## References

  * [1] A. B. Bakushinskii, M. V. Klibanov, and N. A. Koshev.  Carleman weight functions for a globally convergent numerical method for ill-posed Cauchy problems for some quasilinear PDEs.  Nonlinear Anal. Real World Appl., 34:201–224, 2017.
  * [2] G. Barles and P. E. Souganidis.  Convergence of approximation schemes for fully nonlinear second order equations.  Asymptotic Anal., 4(3):271–283, 1991.
  * [3] L. Beilina and M. V. Klibanov.  Approximate Global Convergence and Adaptivity for Coefficient Inverse Problems.  Springer, New York, 2012.
  * [4] A. L. Bukhgeim and Michael V. Klibanov.  Global uniqueness of a class of multidimensional inverse problems.  Soviet Math. Doklady, 17:244–247, 1981.
  * [5] T. Carleman.  Sur les systèmes linéaires aux derivées partielles du premier ordre a deux variables.  C. R. Acad. Sci. Paris, 197:471–474, 1933.
  * [6] David Colton and Rainer Kress.  Inverse acoustic and electromagnetic scattering theory. Applied Mathematical Sciences.  Springer, New York, 3rd edition, 2013.
  * [7] M. G. Crandall, L. C. Evans, and P.-L. Lions.  Some properties of viscosity solutions of Hamilton-Jacobi equations.  Trans. Amer. Math. Soc., 282(2):487–502, 1984.
  * [8] M. G. Crandall and P.-L. Lions.  Viscosity solutions of Hamilton-Jacobi equations.  Trans. Amer. Math. Soc., 277(1):1–42, 1983.
  * [9] M. G. Crandall and P.-L. Lions.  Two approximations of solutions of Hamilton-Jacobi equations.  Math. Comp., 43(167):1–19, 1984.
  * [10] M. Falcone and R. Ferretti.  Semi-Lagrangian schemes for Hamilton-Jacobi equations, discrete representation formulae and Godunov methods.  J. Comput. Phys., 175(2):559–575, 2002.
  * [11] M. Falcone and R. Ferretti.  Semi-Lagrangian approximation schemes for linear and Hamilton-Jacobi equations.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2014.
  * [12] V. A. Khoa, G. W. Bidney, M. V. Klibanov, L. H. Nguyen, L. Nguyen, A. Sullivan, and V. N. Astratov.  Convexification and experimental data for a 3D inverse scattering problem with the moving point source.  Inverse Problems, 36:085007, 2020.
  * [13] V. A. Khoa, G. W. Bidney, M. V. Klibanov, L. H. Nguyen, L. Nguyen, A. Sullivan, and V. N. Astratov.  An inverse problem of a simultaneous reconstruction of the dielectric constant and conductivity from experimental backscattering data.  Inverse Problems in Science and Engineering, 29(5):712–735, 2021\.
  * [14] V. A. Khoa, M. V. Klibanov, and L. H. Nguyen.  Convexification for a 3D inverse scattering problem with the moving point source.  SIAM J. Imaging Sci., 13(2):871–904, 2020.
  * [15] M. V. Klibanov.  Global convexity in a three-dimensional inverse acoustic problem.  SIAM J. Math. Anal., 28:1371–1388, 1997.
  * [16] M. V. Klibanov.  Global convexity in diffusion tomography.  Nonlinear World, 4:247–265, 1997.
  * [17] M. V. Klibanov.  Carleman estimates for the regularization of ill-posed Cauchy problems.  Applied Numerical Mathematics, 94:46–74, 2015.
  * [18] M. V. Klibanov.  Carleman weight functions for solving ill-posed Cauchy problems for quasilinear PDEs.  Inverse Problems, 31:125007, 2015.
  * [19] M. V. Klibanov.  Convexification of restricted Dirichlet to Neumann map.  J. Inverse and Ill-Posed Problems, 25(5):669–685, 2017.
  * [20] M. V. Klibanov and O. V. Ioussoupova.  Uniform strict convexity of a cost functional for three-dimensional inverse scattering problem.  SIAM J. Math. Anal., 26:147–179, 1995.
  * [21] M. V. Klibanov, T. T. Le, L. H. Nguyen, A. Sullivan, and L. Nguyen.  Convexification-based globally convergent numerical method for a 1D coefficient inverse problem with experimental data.  preprint, Arxiv:2104.11392, 2021.
  * [22] M. V. Klibanov and J. Li.  Inverse Problems and Carleman Estimates: Global Uniqueness, Global Convergence and Experimental Data.  to appear on De Gruyter, 2021.
  * [23] M. V. Klibanov, L. H. Nguyen, and H. V. Tran.  Numerical viscosity solutions to Hamilton-Jacobi equations via a Carleman estimate and the convexification method.  preprint Arxiv:2104.05870, 2021.
  * [24] R. Lattès and J. L. Lions.  The Method of Quasireversibility: Applications to Partial Differential Equations.  Elsevier, New York, 1969.
  * [25] M. M. Lavrent’ev, V. G. Romanov, and S. P. Shishat$\cdot$skiĭ.  Ill-Posed Problems of Mathematical Physics and Analysis.  Translations of Mathematical Monographs. AMS, Providence: RI, 1986.
  * [26] T. T. Le and L. H. Nguyen.  A convergent numerical method to recover the initial condition of nonlinear parabolic equations from lateral Cauchy data.  Journal of Inverse and Ill-posed Problems, DOI: https://doi.org/10.1515/jiip-2020-0028, 2020.
  * [27] T. T. Le and L. H. Nguyen.  The gradient descent method for the convexification to solve boundary value problems of quasi-linear PDEs and a coefficient inverse problem.  preprint, Arxiv:2103.04159, 2021.
  * [28] T. T. Le, L. H. Nguyen, T-P. Nguyen, and W. Powell.  The quasi-reversibility method to numerically solve an inverse source problem for hyperbolic equations.  Journal of Scientific Computing, 87:90, 2021.
  * [29] D-L. Nguyen, M. V. Klibanov, L. H. Nguyen, A. E Kolesov, M. A. Fiddy, and H. Liu.  Numerical solution of a coefficient inverse problem with multi-frequency experimental raw data by a globally convergent algorithm.  Journal of Computational Physics, 345:17–32, 2017.
  * [30] H. M. Nguyen and L. H. Nguyen.  Cloaking using complementary media for the Helmholtz equation and a three spheres inequality for second order elliptic equations.  Transaction of the American Mathematical Society, 2:93–112, 2015\.
  * [31] L. H. Nguyen.  A new algorithm to determine the creation or depletion term of parabolic equations from boundary measurements.  Computers and Mathematics with Applications, 80:2135–2149, 2020\.
  * [32] L. H. Nguyen, Q. Li, and M. V. Klibanov.  A convergent numerical method for a multi-frequency inverse source problem in inhomogenous media.  Inverse Problems and Imaging, 13:1067–1094, 2019.
  * [33] L. Novotny and B. Hecht.  Principles of Nano-Optics.  Cambridge University Press, Cambridge, UK, 2nd edition, 2012.
  * [34] S. Osher and R. Fedkiw.  Level set methods and dynamic implicit surfaces, volume 153 of Applied Mathematical Sciences.  Springer-Verlag, New York, 2003.
  * [35] M. H. Protter.  Unique continuation for elliptic equations.  Trans. Amer. Math. Soc., 95(1):81–91, 1960.
  * [36] J. A. Scales, M. L. Smith, and T. L. Fischer.  Global optimization methods for multimodal inverse problems.  J. Computational Physics, 103:258–268, 1992.
  * [37] J. A. Sethian.  Level set methods and fast marching methods, volume 3 of Cambridge Monographs on Applied and Computational Mathematics.  Cambridge University Press, Cambridge, second edition, 1999.  Evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science.
  * [38] P. E. Souganidis.  Approximation schemes for viscosity solutions of Hamilton-Jacobi equations.  J. Differential Equations, 59(1):1–43, 1985.
  * [39] H. V. Tran.  Hamilton–Jacobi equations: Theory and Applications.  Graduate Studies in Mathematics, Vol. 213. American Mathematical Society, Providence, RI, 2021.


utf8

# Speech Representation Analysis
based on Inter- and Intra-Model Similarities

###### Abstract

Self-supervised models have revolutionized speech processing, achieving new
levels of performance in a wide variety of tasks with limited resources.
However, the inner workings of these models are still opaque. In this paper,
we aim to analyze the encoded contextual representation of these foundation
models based on their inter- and intra-model similarity, independent of any
external annotation and task-specific constraint. We examine different SSL
models varying their training paradigm – Contrastive (Wav2Vec2.0) and
Predictive models (HuBERT); and model sizes (base and large). We explore these
models on different levels of localization/distributivity of information
including (i) individual neurons; (ii) layer representation; (iii) attention
weights and (iv) compare the representations with their finetuned
counterparts. Our results highlight that these models converge to similar
representation subspaces but not to similar neuron-localized concepts111A
concept represents a coherent fragment of knowledge, such as “a class
containing certain objects as elements, where the objects have certain
properties”[1]. We made the code publicly available for facilitating further
research, we publicly released our
code222.https://github.com/QCRIVoice/XSSL_speech .

Index Terms—  Self-Supervised Learning, Speech Models, Inter- and Intra-
Similarities

##  1 Introduction

Self-supervised Speech models like Wav2Vec2 [2] and HuBERT [3] have shown
remarkable advancements in a variety of speech processing tasks, including
speech recognition, emotion recognition, speaker verification, and language
identification [4, 5, 6] among others. This significant advancement over
supervised state-of-the-art methods and the opaqueness of these models has
sparked interest in understanding and exploring their internal mechanisms.

Several studies have aimed to understand the information these models capture
about different properties such as speaker characteristics [7, 8, 9, 10],
paralinguistic aspects [11, 12], articulatory features [13], acoustic-
linguistic elements [14], as well as accent features [15] among others.
Moreover, studies like [14, 16] have also shown how better model understanding
can lead to efficient fine-tuning strategies for downstream tasks.

A widely used interpretation technique includes training supervised
classifiers, aka probing classifiers [17, 8, 18], based on the learned
representations of given models, to predict various task properties. This
methodology has found application in various studies and showed the ability of
representations from different models to capture distinct properties.
Additionally, similarity-based methods are used to find associations at the
frame-, phoneme-, and word-levels. These methods utilizes metrics such as
projected-weighted canonical correlation analysis ($pwcca$) [19] and mutual
information without training classifiers. However, the effectiveness of this
approach is limited by the need for annotated data, requiring precise boundary
alignment, accurate phoneme transcription, coupled with consistent word
alignment, to ensure valid and reliable analysis and results.

In this study, we introduce inter- and intra-model similarity measures to
understand contextual representations within speech models. Instead of
focusing on a specific category or property of information, we focus on
exploring both inter- and intra-similarity across a spectrum of models. We
investigate localization/distributivity333Does every single neuron encode a
single concept or all concepts are spread across multiple neurons? [20]
properties in these models. We adopted a set of $5$ distinct similarity
measures, to explore the SSL models for localization/distributivity behavior
at individual neurons, layers, and attention mechanisms levels. This
comprehensive range of metrics allows us to capture the nuances in the
patterns embedded within the models, offering a granular view of their
structural and functional dynamics.

Our in-depth analysis reaffirms prior discoveries without the need for
external data or defined tasks. Moreover, our findings also reveal noteworthy
insights: (i) Speech SSL model neurons exhibit higher intra-model similarity
than inter-model similarity. (ii) Information encapsulated by neurons from one
layer can be represented as a linear combination of other layers. Models have
similar representation subspaces but different localized neuron concepts.
(iii) Lower and adjacent layers demonstrate a high degree of similarity across
diverse models. (iv) The training objective has a greater impact on
representation similarity than the size of the model architecture. (v)
Finally, we show how the similarity analysis can motivate efficient finetuning
for ASR, where freezing the bottom layers of models still maintains comparable
performance to finetuning the full network while reducing the finetuning time.

Fig. 1: Comparison of Heatmap Similarities Between HuBERT and Wav2Vec2.0
Models: $\text{neu-neu}_{\text{sm}}$, $\text{neu-lay}_{\text{sm}}$, and
$pwcca$ Similarities. Model boundaries are highlighted in yellow, and
noteworthy similarities are encircled in green.

##  2 Methodology

We analyzed $M$ pretrained speech SSL models for both localized and
distributed information using various widely accepted similarity metrics [21],
capturing different notations at individual neurons, layers, and attention
levels. We propose to remove any dependency on external labels or boundary
annotation by utilizing only frame-level representation for the study.

For each model $m$, ($m\in M$), we extracted frame level representations
$h^{(m)}_{l}$ $\in\mathbb{R}^{d_{m}}$, where $d_{m}=\\{768,1024\\}$, is
indicative of number of neurons, and attention weights $\alpha^{(m)}_{l}$ at
each layer $l$. We then exploit the extracted neuron/layer-level
representation and attention weight to find inter- and intra-similarities at
various levels of granularities.

###  2.1 Neuron-level Similarity

We adopted two different similarity measures: (i) neuron-neuron similarity
($\text{neu-neu}_{\text{sm}}$) – similarity between pairs of individual
neurons, and (ii) neuron-layer similarity ($\text{neu-lay}_{\text{sm}}$) –
similarity between a neuron in one model with a layer in another.

For a given neuron $h^{(m)}_{l}[k]$ of a model $m$, and a layer $l$,
$\text{neu-neu}_{\text{sm}}$ is defined as the maximum correlation between
$h^{(m)}_{l}[k]$ and another neuron $h^{(m^{\prime})}_{l}[k^{\prime}]$ of
layer $l^{\prime}$ of another model $m^{\prime}$:

$\displaystyle\widetilde{\text{neu-
neu}_{\text{sm}}}(h^{(m)}_{l}[k],h^{(m^{\prime})}_{l^{\prime}})$
$\displaystyle=\max_{k^{\prime}}\rho(h^{(m^{\prime})}_{l^{\prime}}[k^{\prime}],h^{(m)}_{l}[k])$
(1)

Then, we average over all neurons in layer $l$ of the model $m$,

$\displaystyle\text{neu-
neu}_{\text{sm}}(h^{(m)}_{l},h^{(m^{\prime})}_{l^{\prime}})$ $\displaystyle=$
(2) $\displaystyle\frac{1}{d_{m}}\times\sum_{k}\widetilde{\text{neu-
neu}_{\text{sm}}}(h^{(m)}_{l}[k],h^{(m^{\prime})}_{l^{\prime}}).$

where $\rho$ is the Pearson correlation.

$\text{neu-neu}_{\text{sm}}$ is designed to assess the localization of
information, reflecting higher values when two layers exhibit pairs of neurons
that demonstrate similar behavioral patterns.

In contrast, $\text{neu-neu}_{\text{sm}}$ assesses how can a neuron
$h^{(m)}_{l}[k]$ be expressed as linear regression of neurons of another layer
$l^{\prime}$ of another model $m^{\prime}$, and measures the quality of
regression fit which is defined as:

$\displaystyle\widetilde{\text{neu-
lay}_{\text{sm}}}(h^{(m)}_{l}[k],h^{(m^{\prime})}_{l^{\prime}})$
$\displaystyle=lstsq(h^{(m^{\prime})}_{l^{\prime}},h^{(m)}_{l}[k]).r$ (3)

$lstsq$ denotes linear least-squares, and $r$ represents the associated
r-value. As before, this is extended to the layer level:

$\displaystyle\text{neu-
lay}_{\text{sm}}(h^{(m)}_{l},h^{(m^{\prime})}_{l^{\prime}})$
$\displaystyle=\frac{1}{d_{m}}\sum_{k}\widetilde{\text{neu-
lay}_{\text{sm}}}(h^{(m)}_{l}[k],h^{(m^{\prime})}_{l^{\prime}})$ (4)

$\text{neu-lay}_{\text{sm}}$ reflects how the localized information in a
particular neurons in $m$ are distributed across the the layers of models
$m^{\prime}$.

###  2.2 Representation-level Similarity

For layer-level representation analysis, we focus on canonical correlation
analysis (CCA) similarities measures. Despite previous work that focuses on
using projected-weighted CCA similarity ($pwcca$)[14, 16] and singular vector
CCA ($svcca$)[22], we focus on examining similarities among frame
representations instead of frame representation with other information such as
phonemes, words, and boundaries, which typically require extensive annotation
and linguistic expertise. These similarity measures underscore the
distributive nature of information across layers which highlights scenarios
where two layers exhibit similar behaviors across all their neurons,
emphasizing the collective patterns rather than relying solely on individual
neuron matching.

###  2.3 Attention-level Similarity

Similar to $\text{neu-neu}_{\text{sm}}$, $\text{attention}_{\text{sm}}$
similarity identifies the most “correlated” other attention heads within the
model $m$ and across $m^{\prime}$. This measure captures the behavior
similarity indicating the focus alignment. Given two attention heads,
$\alpha_{l}^{m}[k]$ and $\alpha_{l^{\prime}}^{m^{\prime}}[k^{\prime}]$, we
calculate their similarity based on their Pearson correlation, then we average
over the heads in layer $l$ as in Section 2.1.

##  3 Experimental Setup

#### SSL Models

We adopt widely used self-supervised speech models, HuBERT (hub) and
Wav2Vec2.0 (w2v)444Available here: https://huggingface.co/collections/facebook
as reported in Table 1. Both models share similar architectures. The encoder
network consists of blocks of temporal convolution layers with $512$ channels,
and the convolutions in each block have strides and kernel sizes that compress
about $25$ms of $16$kHz audio every $20$ms. The context network consists of
${12\text{ (base)}}$ and $24$ (large) blocks with model dimension ${768\text{
(base)}}$ and $1024$ (large) and attention heads of ${12\text{ (base)}}$ and
$16$ (large). The underlying difference in the models lies in their training
objectives; w2v undergoes training through Contrastive Predictive Coding (CPC)
loss, employing masking techniques, thereby classifying it as a contrastive
model. On the other hand, hub, follows a different approach by attempting to
predict discrete targets of masked regions using Cross-Entropy (CE) loss,
classifying the model as a predictive model.

#### Dataset

We use the extensively employed TIMIT dataset [23] in research, studies on
phone recognition, phone segmentation, and speaker recognition. TIMIT
comprises $5.4$ hours of clean data manually transcribed. Despite its limited
size, the dataset features a diverse set of approximately $630$ speakers
delivering phonetically rich sentences, rendering it favorable for our task.
For our task focuses on studying similarities, we exclusively utilize the
official training set. Given that we employ frame-level embeddings in this
context, each of the selected models yields over $700K$ embeddings from each
layer (each $20$ ms corresponds to an embedding).

Models | (Abbreviation) | Training Data
---|---|---
HuBERT (BASE) | hub-base | Librispeech 960hrs
HuBERT (LARGE) | hub-large | Libri-Light
Wav2Vec 2.0 (BASE) | w2v-base | Librispeech 960hrs
Wav2Vec 2.0 (LARGE) | w2v-large | Libri-Light

Table 1: Examined Pretrained Speech SSL Models

##  4 Analysis & Discussion

#### A. Neuron Intra-Model Similarity.

Figure 1 illustrates heatmaps showing similarities between different neurons
and layers across various models using $\text{neu-neu}_{\text{sm}}$,
$\text{neu-lay}_{\text{sm}}$, and $pwcca$ similarities. The $\text{neu-
neu}_{\text{sm}}$ reveals a distinctive diagonal pattern within each model,
indicating that neurons within a specific model $m$ and layer $l$ tend to
exhibit similarity to their counterparts neurons in adjacent layers within the
same model $m$. However, individual neurons are very different when comparing
a given model $m$ individual neurons to other models neurons. This observation
suggests that neurons exhibit significantly higher intra-model similarities
than inter-model similarities. A similar pattern was found in contextualized
language models [21]. Furthermore, the identified similarity pattern in lower
layer neurons is consistent across all examined models, potentially linked to
their proximity to CNN feature extraction layers equivalent to spectrogram
features, as demonstrated in previous work [24].

Fig. 2: Heatmap of $svcca$ similarity Fig. 3: Comparing Heatmap Similarities
between HuBERT and Wav2Vec2.0 Models in Relation to their ASR Finetuning
Variations: $\text{neu-neu}_{\text{sm}}$ and $pwcca$ Similarities. Noteworthy
similarities are encircled in green.

#### B. Layer Inter-Model Similarity.

While individual pairs exhibit distinct characteristics across different
models, $\text{neu-lay}_{\text{sm}}$ in Figure 1 reveals strong inter-model
similarity, suggesting that the representations of different models converge
to similar subspaces. Furthermore, the results also suggest that the
individual neurons of a model can be represented as a linear combination of
neurons from other layers of the model. These cross-model similarities are
also observed using representation-level similarities $pwcca$, and $svcca$.

#### C. Models within the same family behaves similarly.

Notably, $\text{neu-lay}_{\text{sm}}$ similarity reveals that neurons concepts
in the top layers of hub based models (base and large) are less similar to
lower layers, and vice versa, which supports the fact that the lower layer
captures different fine-grained concepts, whereas the higher layers are
capturing more abstract information as seen in [14, 16, 8]. In contrast, w2v
models show a different trend. A notable similarity is observed in the higher
layers of both base (layer L8 - L10) and large (L20 - L23) models with respect
to all the layers within the model. These intra-model similarities are seen
using both $pwcca$ and $svcca$ similarity measures. Despite the high
similarities in inter-model layer representation (as shown in Section 4.B),
the final layer of the hub and w2v (both base and large) are very distinct.
Our observations indicate that the models within the same family (base and
large) exhibit behavioral similarities in representation. We hypothesize the
uniqueness in representation between the family – hub vs w2v is more likely
attributed to the training objective of self-supervised models rather than the
architecture’s number of layers which aligns with the findings reported in
[17].

#### D. Adjacent and Lower Layers Similarity.

All the heatmaps in Figure 1, including the $svcca$ similarity in Figure 2,
show a bright diagonal and its neighboring areas. This brightness suggests
that neighboring layers share similar representations, indicating that
adjacent layers in the models encapsulate similar information subspaces.
Similar patterns are observed in both language models and vision networks
[25]. Additionally, $pwcca$ discloses that lower layers subspaces demonstrate
similarity across studied models. This alignment with expectations, and with
previous findings in Section 4.A, where lower layers closely resemble CNN
layers functioning as feature extractors, and these features exhibit
equivalence representations across the considered models.

#### E. Attention Weights Similarity.

Examining Figure 4, we observe high similarities between the attention heads
in the upper layers of the models with respect to the lower layers. These high
similarities in the upper layer could indicate redundancy in design. However,
it is important to note that attention-based similarity measures are not
reliable and are harder to interpret, as fine-grained patterns are less
noticeable in the similarity-based analysis.

Fig. 4: Heatmap of $\text{attention}_{\text{sm}}$ similarity

#### F. ASR Finetuning Effect.

Figure 3 depicts the similarities in heatmaps between w2v-large, and hub-
large, along with their ASR fine-tuned counterparts w2v-large-ft and hub-
large-ft, on Librispeech dataset. The analysis utilizes $pwcca$ and
$\text{neu-neu}_{\text{sm}}$ similarities to explore potential changes in
information across different layers and the localized information within
neurons. Results indicate that hub-based models primarily undergo significant
changes after fine-tuning only at the last few layers in comparison to other
layers (similarity between the foundation and its finetuned counterpart is
less than $0.5$ in upper layers). For w2v model, we observe that a large
number of upper layers has changed significantly at both the neuron and layer
levels, in alignment with the findings in [14]. Such findings indicate that
finetuning exclusively the upper layers can be as effective as finetuning the
full model. Our observation is aligned with the findings in [16] where this
conclusion was drawn by examining phoneme-level $cca$ following the fine-
tuning of only the layer $16$ in w2v-large and layer $20$ in hub-large, which
yielded comparable results to finetuning all parameters. Note in [16] used
human annotation for phoneme boundaries, whereas our proposed method gave the
same conclusion without relying on any external annotation.

#### Key Points.

Our study highlights how different models trained with distinct objectives can
converge toward similar representations and concepts. We observed that neurons
in one layer can be expressed as linear combinations of neurons from other
layers in different models. Importantly, this convergence is driven more by
the distributivity nature of representations than by neuron concept
localization. In other words, individual neurons learn different localized
concepts, but overall, they contribute to similar subspaces across layers.

##  5 Conclusion

The paper introduces both annotation- and task-independent approaches for
analyzing various speech SSL models. Our in-depth analysis explores both
Wav2Vec2.0 and HuBERT model families, revealing intricate convergence patterns
in inter- and intra-model neurons, layers, and attention weights similarities.
Our finding suggests that models share similar distributional representations
but different localized concepts, and the training objective emerges as a
pivotal factor, outweighing the influence of model size. Hence, signaling how
understanding the inner workings of such large models can facilitate effective
and parameter-efficient design decisions for both foundation and downstream
models.

## References

  * [1] Stock Wolfgang G,  “Concepts and semantic relations in information science,”  Journal of the American Society for Information Science and Technology, vol. 61, 2010.
  * [2] Baevski Alexei et al.,  “wav2vec 2.0: A framework for self-supervised learning of speech representations,”  Advances in neural information processing systems, 2020.
  * [3] Hsu Wei-Ning et al.,  “Hubert: Self-supervised speech representation learning by masked prediction of hidden units,”  IEEE/ACM Transactions on Audio, Speech, and Language Processing, 2021.
  * [4] Mohamed Abdelrahman et al.,  “Self-supervised speech representation learning: A review,”  IEEE Journal of Selected Topics in Signal Processing, 2022.
  * [5] Shon Suwon al.,  “Slue: New benchmark tasks for spoken language understanding evaluation on natural speech,”  in International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2022.
  * [6] Borgholt Lasse et al.,  “A brief overview of unsupervised neural speech representation learning,”  arXiv preprint arXiv:2203.01829, 2022.
  * [7] Fan Zhiyun et al.,  “Exploring wav2vec 2.0 on speaker verification and language identification,”  arXiv preprint arXiv:2012.06185, 2020.
  * [8] Chowdhury Shammur Absar et al.,  “What do end-to-end speech models learn about speaker, language and channel information? a layer-wise and neuron-level analysis,”  Computer Speech & Language, 2023.
  * [9] Shammur A Chowdhury, Ahmed Ali, Suwon Shon, and James R Glass,  “What does an end-to-end dialect identification model learn about non-dialectal information?,”  in INTERSPEECH, 2020.
  * [10] Feng Chi-Luen et al.,  “Silence is sweeter than speech: Self-supervised model using silence to store speaker information,”  arXiv preprint arXiv:2205.03759, 2022.
  * [11] Shah Jui et al.,  “What all do audio transformer models hear? probing acoustic representations for language delivery and its structure,”  arXiv preprint arXiv:2101.00387, 2021.
  * [12] Li Yuanchao et al.,  “Exploration of a self-supervised speech model: A study on emotional corpora,”  in Spoken Language Technology Workshop (SLT). IEEE, 2023.
  * [13] Ji Hang et al.,  “Predicting within and across language phoneme recognition performance of self-supervised learning speech pre-trained models,”  arXiv preprint arXiv:2206.12489, 2022.
  * [14] Pasad Ankita et al.,  “Layer-wise analysis of a self-supervised speech representation model,”  in Automatic Speech Recognition and Understanding Workshop (ASRU). IEEE, 2021.
  * [15] Yang Mu et al.,  “What can an accent identifier learn? probing phonetic and prosodic information in a wav2vec2-based accent identification model,”  arXiv preprint arXiv:2306.06524, 2023.
  * [16] Pasad Ankita et al.,  “Comparative layer-wise analysis of self-supervised speech models,”  in International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2023.
  * [17] Chung Yu-An et al.,  “Similarity analysis of self-supervised speech representations,”  in International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2021.
  * [18] Belinkov Yonatan and James Glass,  “Analysis methods in neural language processing: A survey,”  Transactions of the Association for Computational Linguistics, 2019.
  * [19] Morcos Ari et al.,  “Insights on representational similarity in neural networks with canonical correlation,”  Advances in neural information processing systems, vol. 31, 2018.
  * [20] James L McClelland, David E Rumelhart, and Geoffrey E Hinton,  “The appeal of parallel distributed processing,”  MIT Press, Cambridge MA, 1986.
  * [21] Wu John et al.,  “Similarity analysis of contextual word representation models,”  in Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics, Online, July 2020, pp. 4638–4655, Association for Computational Linguistics.
  * [22] Raghu Maithra et al.,  “Svcca: Singular vector canonical correlation analysis for deep learning dynamics and interpretability,”  Advances in neural information processing systems, vol. 30, 2017.
  * [23] Garofolo John S,  “Timit acoustic phonetic continuous speech corpus,”  Linguistic Data Consortium, 1993, 1993.
  * [24] Wu Felix et al.,  “Performance-efficiency trade-offs in unsupervised pre-training for speech recognition,”  in International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2022.
  * [25] Kornblith Simon et al.,  “Similarity of neural network representations revisited,”  in International conference on machine learning. PMLR, 2019, pp. 3519–3529.