# Total Positivity of Quasi-Riordan Arrays

Tian-Xiao He  Department of Mathematics, Illinois Wesleyan University,
Bloomington, IL 61702-2900, Indiana, USA  Roksana Słowik  Faculty of Applied
Mathematics, Silesian University of Technology, Kaszubska 23, 44-100 Gliwice,
Poland

###### Abstract

In this paper the total positivity of quasi-Riordan arrays is investigated
with use of the sequence characterization of quasi-Riordan arrays. Due to the
correlation between quasi-Riordan arrays and Riordan arrays, this study is an
in-depth discussion of the total positivity of Riordan arrays.

AMS Subject Classification: 05A15, 05A05, 15B36, 15A06, 05A19, 11B83.

Key Words and Phrases: Riordan array, Riordan group, Quasi-Riordan array,
Almost-Riordan array, Quasi-Riordan group.

##  1 Introduction

Following Karlin [13] and Pinkus [19], an infinite matrix is called totally
positive (abbreviately, TP), if its minors of all orders are nonnegative. An
infinite sequence of nonnegative real numbers $(a_{n})_{n\geq 0}$ is called a
Pólya frequency sequence (abbreviately, PF), if its Toeplitz matrix

$\left[a_{i-j}\right]_{i,j\geq 0}=\left[\begin{array}[]{lllll}a_{0}&&&&\\\
a_{1}&a_{0}&&&\\\ a_{2}&a_{1}&a_{0}&&\\\ a_{3}&a_{2}&a_{1}&a_{0}&\\\
\vdots&\vdots&\vdots&\vdots&\ddots\end{array}\right]$

is TP. We say that a finite sequence $(a_{0},a_{1},\ldots,a_{n})$ is PF if the
corresponding infinite sequence $(a_{0},a_{1},\ldots,a_{n},0,\ldots)$ is PF.
Denote by ${{\mathbb{N}}}$ the set of all nonnegative integers. A fundamental
characterization for PF sequences is given by Schoenberg et al.[1, 2] (see
also Karlin [13]), which states that a sequence $(a_{n})_{n\geq 0}$ is PF if
and only if its generating function can be written as

$\sum_{n\geq 0}a_{n}t^{n}=Ct^{k}e^{\gamma t}\frac{\Pi_{j\geq
0}(1+\alpha_{j}t)}{\Pi_{j\geq 0}(1-\beta_{j}t)},$ (1)

where $C>0$, $k\in{{\mathbb{N}}}$, $\alpha_{j}$, $\beta_{j}$, $\gamma\geq 0$,
and $\sum_{j\geq 0}(\alpha_{j}+\beta_{j})<\infty$. In this case, the above
generating function is called a Pólya frequency formal power series. For more
preliminary materials and more relevant results see, for example, Brenti [6],
Mao, Mu, and Wang [17] and Pinkus [19].

Let $A$ and $B$ be $m\times m$ and $n\times n$ matrices, respectively. Then we
define the direct sum of $A$ and $B$ by

$A\oplus B=\left[\begin{matrix}A&0\\\
0&B\end{matrix}\right]_{(m+n)\times(m+n)}.$ (2)

Let us recall some basic definitions and introduce the notation.

First, the Riordan arrays (or Riordan matrices) are infinite, lower triangular
matrices defined by the generating function of their columns. With the matrix
multiplication, the set of all Riordan arrays form a group, called the Riordan
group (see Shapiro, Getu, Woan and Woodson [22]).

More formally, let us consider the set of formal power series ring
${\mathcal{F}}={{\mathbb{K}}}[\\![$$t$$]\\!]$, where ${{\mathbb{K}}}$ is the
field ${{\mathbb{R}}}$ or ${{\mathbb{C}}}$. The _order_ of
$f(t)\in{\mathcal{F}}$, $f(t)=\sum_{k=0}^{\infty}f_{k}t^{k}$
($f_{k}\in{{\mathbb{K}}}$), is the minimum number $r\in{{\mathbb{N}}}$ such
that $f_{r}\neq 0$, where ${{\mathbb{N}}}=\\{0,1,2\ldots\\}$ is the set of all
natural numbers. We denote by ${\mathcal{F}}_{r}$ the set of formal power
series of order $r$. Let $g(t)\in{\mathcal{F}}_{0}$ and
$f(t)\in{\mathcal{F}}_{1}$; the pair $(g,\,f)$ defines the (proper) Riordan
array $R=(d_{n,k})_{n,k\in{{\mathbb{N}}}}=(g,f)$ having

$d_{n,k}=[t^{n}]g(t)f(t)^{k}$ (3)

or, in other words, having $gf^{k}$ as the generating function of the $k$th
column of $(g,f)$. Hence,

$(g,f):=(g,gf,gf^{2},gf^{3},\ldots),$ (4)

where $g$, $gf$, $gf^{2}$, $gf^{3}\cdots$ are the generating functions of the
$0$th, $1$st, $2$nd, $3$rd, $\cdots$, columns of the matrix $(g,f)$,
respectively. The first fundamental theorem of Riordan arrays can be
represented as

$(g(t),f(t))h(t)=g(t)(h\circ f)(t),$

which can also be simplified to $(g,f)h=gh(f)$. Thus we immediately see that
the usual row-by-column product of two Riordan arrays is also a Riordan array:

$(g_{1},\,f_{1})(g_{2},\,f_{2})=(g_{1}g_{2}(f_{1}),\,f_{2}(f_{1})).$ (5)

The Riordan array $I=(1,\,t)$ is the identity matrix because its entries are
$d_{n,k}=[t^{n}]t^{k}=\delta_{n,k}$.

Let $(g\left(t\right),\,f(t))$ be a Riordan array. Then its inverse is

$(g\left(t\right),\,f(t))^{-1}=\left(\frac{1}{g(\overline{{f}}(t))},\overline{f}(t),\right)$
(6)

where $\overline{f}(t)$ is the compositional inverse of $f(t)$, i.e.,
$(f\circ\overline{f})(t)=(\overline{f}\circ f)(t)=t$. In this way, the set
$\mathcal{R}$ of all proper Riordan arrays forms a group (see [22]) called the
Riordan group.

Here is a list of six important subgroups of the Riordan group (see [22, 21,
4]).

  * •

the Appell subgroup $\\{(g(t),\,t):g\in{\mathcal{F}}_{0}\\}$.

  * •

the Lagrange (associated) subgroup $\\{(1,\,f(t):f\in{\mathcal{F}}_{1})\\}$.

  * •

the $k$-Bell subgroup $\\{(g(t),\,t(g(t))^{k}):g\in{\mathcal{F}}_{0}\\}$,
where $k$ is a fixed positive integer.

  * •

the hitting-time subgroup
$\\{(tf^{\prime}(t)/f(t),\,f(t)):f\in{\mathcal{F}}_{1}\\}$.

  * •

the derivative subgroup $\\{(f^{\prime}(t),\,f(t)):f\in{\mathcal{F}}_{1}\\}$.

  * •

the checkerboard subgroup
$\\{(g(t),\,f(t)):g\in{\mathcal{F}}_{0},f\in{\mathcal{F}}_{1}\\},$ where $g$
is an even function and $f$ is an odd function.

The $1$-Bell subgroup is referred to as the Bell subgroup for short, and the
Appell subgroup can be considered as the $0$-Bell subgroup if we allow $k=0$
to be included in the definition of the $k$-Bell subgroup.

An infinite lower triangular matrix $[d_{n,k}]_{n,k\in{{{\mathbb{N}}}}}$ is a
Riordan array if and only if a unique sequence
$A=(a_{0}\not=0,a_{1},a_{2},\ldots)$ exists such that for every
$n,k\in{{{\mathbb{N}}}}$

$d_{n+1,k+1}=a_{0}d_{n,k}+a_{1}d_{n,k+1}+\cdots+a_{n}d_{n,n}.$ (7)

This is equivalent to

$f(t)=tA(f(t))\quad\text{or}\quad t=\bar{f}(t)A(t).$ (8)

Here, $A(t)$ is the generating function of the $A$-sequence. The first formula
of (8) is also called the second fundamental theorem of Riordan arrays.

Moreover, there exists a unique sequence $Z=(z_{0},z_{1},z_{2},\ldots)$ such
that every element in column $0$ can be expressed as the linear combination

$d_{n+1,0}=z_{0}d_{n,0}+z_{1}d_{n,1}+\cdots+z_{n}d_{n,n},$ (9)

or equivalently,

$g(t)=\frac{1}{1-tZ(f(t))},$ (10)

in which and thoroughly we always assume $g(0)=g_{0}=1$, a usual hypothesis
for proper Riordan arrays. From (10), we may obtain

$Z(t)=\frac{g(\bar{f}(t))-1}{\bar{f}(t)g(\bar{f}(t))}.$

$A$\- and $Z$-sequence characterizations of Riordan arrays were introduced,
developed, and/or studied in Merlini, Rogers, Sprugnoli, and Verri [18], Roger
[20], Sprugnoli and the author [11], Jean-Louis and Nkwanta [12], Luzón,
Morón, Prieto-Martinez [15, 16], and [10] as well as their references.

Secondly, we recall the definition of quasi-Riordan arrays given in [9].

######  Definition 1.1.

[9] Let $g\in\mbox{$\mathcal{F}$}_{0}$ with $g(0)=1$ and
$f\in\mbox{$\mathcal{F}$}_{1}$. We call the following matrix a quasi-Riordan
array and denote it by $[g,f]$.

$[g,f]:=(g,f,tf,t^{2}f,\ldots),$ (11)

where $g$, $f$, $tf$, $t^{2}f\cdots$ are the generating functions of the
$0$th, $1$st, $2$nd, $3$rd, $\cdots$, columns of the matrix $[g,f]$,
respectively. It is clear that $[g,f]$ can be written as

$[g,f]=\left(\begin{matrix}g(0)&0\\\ (g-g(0))/t&(f/t,t)\end{matrix}\right),$
(12)

where $(f,t)=(f,tf,t^{2}f,t^{3}f,\ldots)$ and $(f/t,t)$ is an Appell Riordan
array. Particularly, if $f=tg$, then the quasi-Riordan array
$[g,tg]=(g,t)$,the Appell-type Riordan array.

The relationship between the Riordan arrays and quasi-Riordan arrays can be
presented in the following theorem.

######  Theorem 1.2.

[9] Let $(g,f)=(d_{n,k})_{n,k\geq 0}$ be a Riordan array, and let
$([1]\oplus(g,f))$ and $[g,f]$ be defined by (2) and (11), respectively. Then
$(g,f)$ has the expression

$(g,f)=[g,f]([1]\oplus(g,f)).$ (13)

For an integer $r\geq 0$, denote the $r$-th truncation of a formal power
series $h=\sum_{n\geq 0}h_{n}t^{n}$ by $h|_{r}:=\sum^{r}_{n=0}h_{n}t^{n}$. For
an integer $n\geq 0$, denote by $(g,f)_{n}$ the $n$-th order leading principle
submatrix of the Riordan array $(g(t),f(t))_{n}$. This notation is introduced
by Luzón et al in [14] and to the one used later by Luzón and her
collaborators. Similarly, $[g,f]_{n}$ denotes the $n$-th order leading
principle submatrix of the quasi-Riordan array $[g,f]$ defined by (11), namely

$[g,f]_{n}=\left(\begin{matrix}g(0)&0\\\
((g-g(0))/t)|_{n-1}&(f,t)_{n-1}\end{matrix}\right),$ (14)

where the $n-1$-st truncation of $(g-g(0)/t)|_{n-1}$ is

$(g-g(0)/t)|_{n-1}=\sum^{n-1}_{k=1}g_{k}t^{k-1},$

and $(f,t)_{n-1}=(f,tf,t^{2}f,\ldots,t^{n-1}f).$ We call $[g,f]_{n}$ the
$n$-th truncation of the quasi-Riordan array $[g,f]$.

######  Corollary 1.3.

[9, 17] Let $(g,f)=(d_{n,k})_{n,k\geq 0}$ be a Riordan array with
$g=\sum_{n\geq 0}g_{n}t^{n}$ and $f=\sum_{n\geq 1}f_{n}t^{n}$, and let
$(g,f)_{n}$ be the $n$-th order leading principle submatrix of $(g,f)$. Then
we have the recurrence relation of $(g,f)_{n}$ in the following form:

$(g,f)_{n}=[g,f]_{n}([1]\oplus(g,f)_{n-1}).$ (15)

Let $R=(g(t),f(t))$ be a Riordan array, where $g(t)=\sum_{n\geq 0}g_{n}t^{n}$
and $f(t)=\sum_{n\geq 1}f_{n}t^{n}$. If the lower triangular matrix

$Q=\left[\begin{array}[]{llllll}g_{0}&0&0&0&0&\cdots\\\
g_{1}&f_{1}&0&0&0&\cdots\\\ g_{2}&f_{2}&f_{1}&0&0&\cdots\\\
g_{3}&f_{3}&f_{2}&f_{1}&0&\cdots\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\end{array}\right]=\left(g(t),f(t),tf(t),t^{2}f(t),\ldots\right)$
(16)

is totally positive (TP), then so is $R$ [17]. This can be proved easily by
using Corollary 1.3 and mathematical induction.

The question is when $Q$ is totally positive. It is obvious that the Pólya
frequency property for the formal power series $g$, $f$ is not sufficient to
ensure TP of $Q$.

######  Example 1.4.

Let $g(t)=(1+t)^{2}$ and $f(t)=t/(1-t)$. Then both $g(t)$ and $f(t)$ are Pólya
frequency, and $(g,f)$ is TP. However, the minor of $[g,f]$

$M^{1,2,3}_{0,1,2}=\det\left(\begin{array}[]{lll}2&1&0\\\ 1&1&1\\\
0&1&1\end{array}\right)=-1<0.$

Hence, $[g,f]$ is not TP although $(g,f)$ is TP.

In this paper, we search for the conditions of TP quasi-Riordan arrays and the
sequence characterization of TP quasi-Riordan arrays. From the last example,
we know that the formal power series $g(t)$ plays an important role in
identifying whether a quasi-Riordan array is TP. Hence, in the next section,
we discuss the conditions and their equivalents for TP quasi-Riordan arrays
$[g,f]$ with some special power series $g$. In Section 3, we study the
sequence characterization of TP quasi-Riordan arrays. In particular, we
present some generalizations of results given in [7]. More specifically,
Theorems 3.2 and 3.4 in this section are generalizations of Theorem 2.1 and
Theorem 2.3 of [7] for quasi-Riordan arrays. In Section 4 we will present some
more examples with various $g$’s for TP and non-TP quasi-Riordan arrays
$[g,f]$, respectively. From a list of examples showing in this section and the
results in the previous sections, one may see a kind of independence of a TP
quasi-Riordan array $[g,f]$ and the properties of that its formal power series
$g$ and $f$ are the Pólya frequency formal power series.

##  2 TP Quasi-Riordan Arrays

Let us consider two special types of Riordan arrays and discuss their TP
property and the TP property for the corresponding quasi-Riordan array.

If $g(t)=1$, then from (12) we have

$[g,f]=\left(\begin{matrix}1&0\\\ 0&(f/t,t)\end{matrix}\right).$ (17)

Hence,

######  Proposition 2.1.

Let $f\in\mathcal{F}_{1}$. The following conditions are equivalent.

  1. 1.

$f$ is a Pólya frequency sequence.

  2. 2.

The Riordan array $(f/t,t)$ is TP.

  3. 3.

The Riordan array $(1,f)$ is TP.

  4. 4.

The quasi-Riordan array $\left[1,\frac{f(t)}{t}\right]$ is TP.

###### Proof.

Since $(f/t,t)$ is a Riordan array, from the definition of the TP matrices, we
immediately obtain Proposition 2.1. ∎

Consider the case when $g(t)=1+g_{1}t$ and $f(t)=f_{1}t+f_{2}t^{2}$. Then

$[g,f]=\left[\begin{array}[]{llll}1&0&&\\\ g_{1}&f_{1}&0&\\\ &f_{2}&f_{1}&\\\
&&f_{2}&\\\ &\ddots&\ddots&\ddots\end{array}\right].$

Hence, from the proof of Proposition 2.6 of [7] we have

######  Corollary 2.2.

The quasi-Riordan array $[1+g_{1}t,f_{1}t+f_{2}t^{2}]$ is TP if and only if
$g_{1},f_{1},f_{2}\geqslant 0$.

Proposition 2.1 can be somewhat generalized to the following.

######  Proposition 2.3.

If $[g,f]$ is $TP$, then $f$ is a Pólya frequency sequence.

Clearly, the above necessary condition is absolutely not sufficient for total
positivity of $[g,f]$. Primarily, if the terms of $g$ are not nonnegative,
then $[g,f]$ is definitely not TP, even if $f$ is a Pólya frequency sequence.
Moreover, it should also be emphasized again that, although the fact that both
$g$ and $f$ are Pólya frequency sequences implies that the Riordan array
$(g,f)$ is TP [8], the analogous result for quasi-Riordan arrays does not hold
as we have seen in Example 1.4.

##  3 Sequence Characterization of the TP of Quasi-Riordan Arrays

Let us now bring the attention to the fact that the quasi-Riordan array
$[g,f]$ can be treated as almost-Riordan array
$\left(g|\frac{f(t)}{t},t\right)$ [5]. From [3], we obtain the generating
functions of the $A$-, $Z$-, and $W$-sequences of the latter one as follows:

$\displaystyle A(t)=1,$ (18) $\displaystyle
Z(t)=\frac{f(t)-z_{0}tg(t)}{f(t)}+z_{0},$ (19) $\displaystyle
W(t)=\frac{(1-w_{0}t)g(t)-1}{f(t)}+w_{0},$ (20)

where $z_{0}=f_{1}$ and $w_{0}=g_{1}$. For instance, for quasi-Riordan array
$[1/(1-t),t/(1-t)]$, $A(t)=Z(t)=W(t)=1$.

We call the following matrix the production matrix of the quasi-Riordan matrix
$[g,f]$:

$J=\left[\begin{array}[]{llllll}w_{0}&z_{0}&&&&\\\ w_{1}&z_{1}&a_{0}&&&\\\
w_{2}&z_{2}&a_{1}&a_{0}&&\\\ w_{3}&z_{3}&a_{2}&a_{1}&a_{0}&\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\end{array}\right]$ (21)

######  Proposition 3.1.

Let $[g,f]$ be a quasi-Riordan array, and let $J$ be the matrix defined by
(21). Then,

$[g,f]J=\overline{[g,f]},$ (22)

where $\overline{[g,f]}$ is the matrix by deleting the first row of $[g,f]$.

###### Proof.

Equation (22) can be checked straightforwardly by using matrix multiplication
and definition of quasi-Riordan array. ∎

######  Theorem 3.2.

Let $J$ be a matrix defined by (21) and let $R$ be the quasi-Riordan array
defined by $A$-, $Z$-, and $W$-sequences. If $J$ is TP, then so is $R$.

###### Proof.

It can be proved by using (22) and mathematical induction. ∎

######  Corollary 3.3.

Let $[g,f]$ and $J$ be two matrices defined by (11) and (21), respectively. If
$J$ is TP, then so are $[g,f]$ and $(g,f)$.

###### Proof.

The corollary can be proved by using Theorems 1.2 and 3.2. ∎

Clearly, by (18)-(20) the matrix $J$ shown above is of a special form. Namely

$J=\left[\begin{array}[]{lllllll}w_{0}&z_{0}&&&&&\\\ w_{1}&z_{1}&1&&&&\\\
w_{2}&z_{2}&0&1&&&\\\ w_{3}&z_{3}&0&0&1&&\\\ w_{4}&z_{4}&0&0&0&1&\\\
\vdots&\vdots&&&&&\ddots\\\ \end{array}\right].$ (23)

Then, straight from the definition of total positivity one derives:

######  Theorem 3.4.

The matrix $J$ given by $(\ref{eq:J_for_[g,f]})$ is TP if and only if the two
followings hold:

  1. 1.

$w_{k}=z_{k}=0$ for all $k\geqslant 2$ and $w_{0},w_{1},z_{0},z_{1}\geq 0$

  2. 2.

$w_{0}z_{1}-w_{1}z_{0}\geqslant 0$.

###### Proof.

A matrix is TP if and only if its all minors are nonnegative.

Since minors of order $1$ must be nonnegative, we must have
$w_{k},z_{k}\geqslant 0$.

If a minor does not consist neither the $0$-th nor the $1$-st row, then it is
a minor of an identity matrix with an extra $0$ row on the top which is TP.

Now notice that choosing a minor consisting of the first and $k$-th
($k\geqslant 2$) row and $0$-th and second column we get

$\left|\begin{array}[]{ll}w_{1}&1\\\ w_{k}&0\\\
\end{array}\right|=-w_{k}\geqslant 0,$

so comparing with the previous observation $w_{k}=0$ for all $k\geqslant 2$.

Analogously for $z_{k}$.

Thus $J$ takes now the form

$\left[\begin{array}[]{lllllll}w_{0}&z_{0}&&&&&\\\ w_{1}&z_{1}&1&&&&\\\
0&0&0&1&&&\\\ 0&0&0&0&1&&\\\ 0&0&0&0&0&1&\\\ \vdots&\vdots&&&&&\ddots\\\
\end{array}\right]$

Consider the minors containing the $0$-th column and some columns $j_{2}$,
$\ldots$, $j_{n}$, where $1<j_{2}<\cdots<j_{n}$.

  1. 1.

If a minor does not contain any of the rows $0$ or $1$, then its first column
(coming after the $0$th one) is equal to $0$ and so is our minor.

  2. 2.

If a minor contains the $0$-th row, then it is triangular and equal either to
$w_{0}$ or $0$.

  3. 3.

If a minor does not contain the $0$-th row, then exapanding it along this row
we get that it is equal either to $0$ or $w_{1}$.

For the minors containing the first column and some columns $j_{2}$, $\ldots$,
$j_{n}$, where $1<j_{2}<\cdots<j_{n}$, the conclusions are the same.

Consider now minors containing both columns $0$ and $1$.

  1. 1.

If a only the $0$-th row or the first row, then the first two columns are
linearly dependent, so the minor is equal to $0$.

  2. 2.

Suppose that a minor contains the $0$-th and the first row. If its order is
$2$, then it is simply equal to $w_{0}z_{1}-w_{1}z_{0}$. Otherwise, it is a
minor of a partitioned matrix of the form

$\left[\begin{array}[]{l|l}A&B\\\ \hline\cr 0&C\\\
\end{array}\right],\quad\textrm{with
}A=\left[\begin{array}[]{ll}w_{0}&z_{0}\\\ w_{1}&z_{1}\\\ \end{array}\right]$

and the first column of $C$ is zero, so the whole minor is equal to $0$ as
well.

∎

Joining the results of Theorems 3.2, 3.4 and
$(\ref{eq:A_Z_W-1})$-$(\ref{eq:A_Z_W-3})$ one can derive the following.

######  Corollary 3.5.

Suppose $w_{0},w_{1},z_{0},z_{1}\geqslant 0$ and
$w_{0}z_{1}-w_{1}z_{0}\geqslant 0$. Then the quasi-Riordan array

$\left[\frac{1-z_{1}t}{(w_{0}z_{1}-w_{1}z_{0})t^{2}-(w_{0}+z_{1})t+1},\frac{z_{0}t}{(w_{0}z_{1}-w_{1}z_{0})t^{2}-(w_{0}+z_{1})t+1}\right]$

is TP.

###### Proof.

Substituting $A(t)=1$, $Z(t)=z_{0}+z_{1}t$ and $W(t)=w_{0}+w_{1}t$ into
$(\ref{eq:A_Z_W-2})$ and (20) yields the system

$\begin{cases}z_{0}tg(t)+(z_{1}t-1)f(t)=0\\\ (w_{0}t-1)g(t)+w_{1}tf(t)=-1\\\
\end{cases}$

which has the solution

$\begin{array}[]{l}g(t)=\frac{1-z_{1}t}{(w_{0}z_{1}-w_{1}z_{0})t^{2}-(w_{0}+z_{1})t+1},\\\
f(t)=\frac{z_{0}t}{(w_{0}z_{1}-w_{1}z_{0})t^{2}-(w_{0}+z_{1})t+1}.\\\
\end{array}$

By Theorem 3.4 $J$ is TP. Thus, by Proposition 3.2 so is the obtained quasi-
Riordan array. ∎

######  Example 3.6.

By Corollary 3.5 the quasi-Riordan array (with $w_{0}=z_{0}=1$, $w_{1}=2$,
$z_{1}=3$):

$\left[\frac{1-3t}{t^{2}-4t+1},\frac{t}{t^{2}-4t+1}\right]=\left[\begin{array}[]{llllll}1&&&&&\\\
1&1&&&&\\\ 3&4&1&&&\\\ 11&15&4&1&&\\\ 41&56&15&4&1&\\\ \vdots&&&&&\ddots\\\
\end{array}\right]$

is TP.

Note that since the discriminant of the denominator $D(t)$ of the obtained $g$
and $f$ is equal to

$\Delta=w_{0}^{2}+z_{1}^{2}-2w_{0}z_{1}+4w_{1}z_{0}=(w_{0}-z_{1})^{2}+4w_{1}z_{0},$

all the coefficients are nonnegative, so $D(t)$ has two real roots. Moreover,
according to the assumption of $w_{0}z_{1}-w_{1}z_{0}\geqslant 0$, their sum
$\frac{w_{0}+z_{1}}{w_{0}z_{1}-w_{1}z_{0}}$ and their product
$\frac{1}{w_{0}z_{1}-w_{1}z_{0}}$ are nonnegative. Then, they must be
nonnegative as well. Thus, $f$ is of form (1), so by [1, 2, 13] this confirms
that $f$ must be Pólya frequency formal power series. If $z_{1}=0$, then by
$w_{0}z_{1}-w_{1}z_{0}\geq 0$ we have $w_{1}z_{0}=0$, which implies $w_{1}=0$
due to $z_{0}\not=0$. Hence, if $z_{1}=0$, then $z_{0}>0$, $w_{1}=0$, and
$w_{0}>0$, and all those make
$[g,f]=\left[\frac{1}{1-w_{0}t},\frac{z_{0}t}{1-w_{0}t}\right]$, where both
$f$ and $g$ are Pólya frequency. If $z_{1}>0$, then $g$ is not a Pólya
frequency formal power series, although both $[g,f]$ and, hence, $(g,f)$ are
TP. This gives another counterexample for Item 2 in [24]. As what we have
seen, if $z_{1}>0$, then $g$ is not a Pólya frequency formal power series– its
numerator has a positive root. However, as Proposition 3.2 holds,
$[g(t),f(t)]$ is TP, so the coefficients of the series $g(t)$ are nonnegative.
Moreover, we have already seen in Example 1.4 that, in general, fact that both
series $g(t)$ and $f(t)$ are Pólya frequency is not sufficient for $[g,f]$ to
be TP, so one could mistrust the necessity of this condition as well.

##  4 Counterexamples

Note that in Example 1.4 we have already presented a non-TP quasi-Riordan
array $[g,f]$, where $g$ and $f$ are Pólya frequency power series. Here, in
the first subsection, we give some more examples like that. In the second
subsection, we give some more examples for TP of quasi-Riordan arrays $[g,f]$
with non-Pólya frequency power series $g$.

###  4.1 Quasi-Riordan arrays that are not TP

######  Proposition 4.1.

Let $k_{1},k_{2}$ be natural numbers such that $k_{2}>k_{1}$, $\alpha>0$, and
let $f(t)=\sum_{k\geqslant 1}f_{k}t^{k}$ be a Pólya frequency power series. If
for some positive natural number $n$ it holds
$\alpha^{k_{2}-k_{1}}f_{k_{1}-n+1}>f_{k_{2}-n+1}$, then the quasi-Riordan
array $\left[\frac{1}{1-\alpha t},f(t)\right]$ is not TP, although both series
$g$ and $f$ are Pólya frequency power series.

###### Proof.

Consider a minor of $\left[\frac{1}{1-\alpha t},f(t)\right]$ that consists of
the entries from the rows $k_{1}$, $k_{2}$ and the columns $0$ and $n$. It is
equal to

$M_{0,n}^{k_{1},k_{2}}=\alpha^{k_{1}}f_{k_{2}-n+1}-\alpha^{k_{2}}f_{k_{1}-n+1},$

so by the assumption it is negative, and so $\left[\frac{1}{1-\alpha
t},f(t)\right]$ is not TP. ∎

######  Example 4.2.

Let

$f(t)=\frac{t}{(1-2t)^{2}}=t+4t^{2}+12t^{3}+32t^{4}+80t^{5}+192t^{6}+\cdots$

Choosing $k_{1}=3$, $k_{2}=4$, $n=1$, for $\alpha=3$ we have

$\alpha^{4-3}f_{3}>f_{4}.$

Hence,

$\left[\frac{1}{1-3t},\frac{t}{(1-2t)^{2}}\right]=\left[\begin{array}[]{llllll}1&&&&&\\\
3&1&&&&\\\ 9&4&1&&&\\\ 27&12&4&1&&\\\ 81&32&12&4&1&\\\ \vdots&&&&&\ddots\\\
\end{array}\right]$

is not TP, as

$M^{0,1}_{3,4}=\left|\begin{array}[]{ll}27&12\\\ 81&32\\\
\end{array}\right|=-108<0,$

although $g$ and $f$ are Pólya frequency power series.

In general, Proposition 4.1 implies the following.

######  Corollary 4.3.

For every Pólya frequency power series $f$, there exists a constant $\alpha>0$
such that the quasi-Riordan array $\left[\frac{1}{1-\alpha t},f(t)\right]$ is
not TP, although both $\frac{1}{1-\alpha t}$ and $f(t)$ are Pólya frequency
power series.

###### Proof.

Clearly, fixing arbitrarily indices $k_{1},k_{2},n$, $k_{1}<k_{2}$, $n>0$, it
suffices to take

$\alpha>\sqrt[k_{2}-k_{1}]{\frac{f_{k_{2}-n}}{f_{k_{1}-n}}}.$ (24)

∎

######  Example 4.4.

Let

$f(t)=\frac{t(t+1)}{1-2t}=t+\sum_{n=2}^{\infty}\left(2^{n-2}+2^{n-1}\right)t^{n}.$

To find $\alpha$ appearing in Proposition 4.1 we consider the inequality

$\alpha>\sqrt[k_{2}-k_{1}]{\frac{f_{k_{1}-n}}{f_{k_{2}-n}}}.$

Taking, for instance $k_{2}=n+2$, $k_{1}=n+1$, we get

$\alpha>\frac{f_{2}}{f_{1}}=3.$

Thus, for any $\alpha>3$ the quasi-Riordan $\left[\frac{1}{1-\alpha
t},f(t)\right]$ is not TP.

Series $g$, that in from Proposition 4.1 is equal to $\frac{1}{1-\alpha t}$,
can be generalized to various other forms. A formula for this general form
seems difficult to be found, yet let us present one more example.

######  Example 4.5.

Let $f(t)$ be a Pólya frequency power series and let

$g(t)=\frac{1}{(1-\alpha t)(1-\beta
t)}=\sum_{n=0}^{\infty}\frac{\beta^{n+1}-\alpha^{n+1}}{\beta-\alpha}t^{n}\quad\textrm{with
}\beta>\alpha,$

i.e. $g(t)$ is also a Pólya frequency power series.

Looking at the minor

$M^{1,2}_{0,1}=\left|\begin{array}[]{ll}\frac{\beta^{2}-\alpha^{2}}{\beta-\alpha}&f_{0}\\\
\frac{\beta^{3}-\alpha^{3}}{\beta-\alpha}&f_{1}\\\ \end{array}\right|$

it can be observed that it is negative if and only if

$\frac{\alpha^{2}+\alpha\beta+\beta^{2}}{\alpha+\beta}>\frac{f_{1}}{f_{0}}\quad\Leftrightarrow\quad\alpha^{2}+\beta^{2}+\alpha\beta-\frac{f_{1}}{f_{0}}\alpha-\frac{f_{1}}{f_{0}}\beta>0.$

The above inequality has infinitely many solutions. One of the examples can be
the pair $\alpha=\frac{f_{1}}{2f_{0}}$, $\beta=\frac{f_{1}}{f_{0}}$. More
precisely, the set of all solutions is the exterior of the ellipse

$\frac{\left(x-\frac{\sqrt{2}f_{1}}{3f_{0}}\right)^{2}}{\frac{2f_{1}^{2}}{9f_{0}^{2}}}+\frac{y^{2}}{\frac{2f_{1}^{2}}{3f_{0}^{2}}}=1$

rotated by the angle $-\frac{\pi}{4}$. The example of this area for
$\frac{f_{1}}{f_{0}}=2$ is shown on Fig.1.

Figure 1: The area depicting pairs $(\alpha,\beta)$ for which
$\left[\frac{1}{(1-\alpha t)(1-\beta t)},f(t)\right]$ is not TP,
$\frac{f_{1}}{f_{0}}=2$. Note that the fact that point $(\alpha,\beta)$ lies
inside the ellipse does not necessarily mean that $[g,f]$ is TP.

###  4.2 Quasi-Riordan arrays that are TP

######  Example 4.6.

Let $f(t)=\frac{t}{1-2t}$, $g(t)=1+t+t^{2}$. Then $f$ is a Pólya frequency
power series, $g$ is not. Consider

$R=[g,f]=\left[\begin{array}[]{llllll}1&&&&&\\\ 1&1&&&&\\\ 1&2&1&&&\\\
0&4&2&1&&\\\ 0&8&4&2&1&\\\ \vdots&&&&&\ddots\\\ \end{array}\right]$

and its minors.

Clearly, those minors that do not involve the $0$th column are the minors of
$\left[1,f(t)\right]$, so they are nonnegative. Moreover, only those that are
either of order $1$ or those that are triangular minors are positive, the
other ones are equal to $0$.

Consider now a minor of order $n$ ($n\geqslant 2$) containing the $0$th
column. If it contains the $0$th row, then expanding along the $0$th row, we
get that it is equal to some minor of $\left[1,f(t)\right]$, so it is
nonnegative. Suppose that it does not not contain the $0$th row. Then let us
expand it along the $0$th column. If it contains one of the $1$st row and the
$2$ row, expanding along the $0$th column, we again get a minor of
$\left[1,f(t)\right]$. Suppose that it contains both $1$st and $2$nd row.
Again, expanding along the $0$th column, we get then the difference of two
minors of $\left[1,f(t)\right]$. If $n=2$, then we deal with the minor

$\left|\begin{array}[]{ll}1&r_{k\>1}\\\ 1&r_{k\>2}\\\
\end{array}\right|\in\left\\{\left|\begin{array}[]{ll}1&1\\\ 1&2\\\
\end{array}\right|,\;\left|\begin{array}[]{ll}1&0\\\ 1&1\\\
\end{array}\right|,\;\left|\begin{array}[]{ll}1&0\\\ 1&0\\\
\end{array}\right|\right\\}$

that is nonnegative.

If $n\geqslant 3$, both minors are of order greater or equal to $3$, so they
are equal to $0$, and so is their difference. Hence, $R$ is TP.

Example 4.6 can be extended to the following.

######  Corollary 4.7.

Let $\alpha>0$, and let $g_{0},g_{1},g_{2}>0$ be such that
$g_{1}^{2}-4g_{0}g_{2}<0$. If $g_{1}\alpha-g_{2}\geqslant 0$, then
$R=\left[g_{0}+g_{1}t+g_{2}t^{2},\frac{t}{1-\alpha t}\right]$ is TP, although
$g_{0}+g_{1}t+g_{2}t^{2}$ is not a Pólya frequency power series.

###### Proof.

Analogously as in the example, the only minors that are nonzero, are the
triangular ones and the minor

$\left|\begin{array}[]{ll}r_{10}&r_{11}\\\ r_{20}&r_{21}\\\
\end{array}\right|=\left|\begin{array}[]{ll}g_{1}&1\\\ g_{2}&\alpha\\\
\end{array}\right|=g_{1}\alpha-g_{2}$

which, by assumption, is nonnegative. ∎

Note that replacing $g(t)=1+t+t^{2}$ from Example 4.6 by, for instance,
$g(t)=t^{2}+1$ (that does not satisfy the assumption $g_{1}>0$ of Corollary
4.7) yields the array of the form

$\left[\begin{array}[]{lllll}1&&&&\\\ 0&f_{1}&&&\\\ 1&f_{2}&f_{1}&&\\\
0&f_{3}&f_{2}&f_{1}&\\\ \vdots&&&&\ddots\\\ \end{array}\right]$

whose minor $M^{1,2}_{0,1}$ is equal to $-f_{1}<0$, so the above array is not
TP.

##  5 Closing comments

In the present work we have studied total positivity of quasi-Riordan arrays.
As we have mentioned, these are, in fact, some particular cases of almost-
Riordan arrays. Thus, it is natural question whether the presented results can
be generalized for almost-Riordan arrays. Considering the capacity of the
paper, the TP of almost-Riordan arrays will be discussed in a subsequent
paper. For the same reason, some other related topics, such as analogues of
$TPr$ Riordan arrays in [7], $TP_{r}$ quasi-Riordan arrays and $TPr$ almost-
Riordan arrays, will also be discussed in a subsequent paper.

##  6 Data availability statement

All data generated or analyzed during this study are included in this article.

## References

  * [1] M. Aissen, A. Edrei, I.J. Schoenberg and A.M. Whitney, On the generating functions of totally positive sequences, Proc. Nat. Acad. Sci. U.S.A. 37 (1951), 303–307.
  * [2] M. Aissen, I.J. Schoenberg and A.M. Whitney, On the generating functions of totally positive sequences I, J. Analyse Math. 2 (1952), 93–103.
  * [3] Y. Alp and E. G. Kocer, Sequence characterization of almost-Riordan arrays, Linear Algebra Appl. 664 (2023), 1–23.
  * [4] P. Barry, Riordan Arrays: A Primer, LOGIC Press, Kilcock, 2016.
  * [5] P. Barry, On the group of almost-Riordan arrays, arXiv:1606.05077v1.
  * [6] F. Brenti, Combinatorics and total positivity, J. Combin. Theory Ser. A, 71 (1995) 175–218.
  * [7] X. Chen, H. Liang and Y. Wang, Total positivity of Riordan arrays, European J. Combin. 46 (2015) 68–74.
  * [8] X. Chen and Y. Wang, Notes on the total positivity of Riordan arrays, Linear Algebra Appl. 569 (2019) 156–161.
  * [9] T.-X. He, The vertical recursive relation of Riordan arrays and its matrix representation, J. Integer Seq. 25 (2022), no. 9, Art. 22.9.5, 22 pp.
  * [10] T. -X. He, Matrix characterizations of Riordan matrices, Linear Algebra Appl. 465 (2015), 15-42.
  * [11] T. -X. He and R. Sprugnoli. Sequence Characterization of Riordan Arrays, Discrete Math., 309 (2009), 3962-3974.
  * [12] C. Jean-Louis and A. Nkwanta, Some algebraic structure of the Riordan group, Linear Algebra Appl. 438 (2013), 2018–2035.
  * [13] S. Karlin, Total Positivity, Vol.1, Stanford University Press, 1968.
  * [14] A. Luzó n, D. Merlini, M. A. Moró n, L. F. Prieto-Martinez, and R. Sprugnoli, Some inverse limit approaches to the Riordan group, Linear Algebra Appl. 491 (2016), 239–262.
  * [15] A. Luzó n, M. A. Moró n, and L. F. Prieto-Martinez, A formula to construct all involutions in Riordan matrix groups. Linear Algebra Appl. 533 (2017), 397–417.
  * [16] A. Luzó n, M. A. Moró n, and L. F. Prieto-Martinez, The group generated by Riordan involutions, Revista Matemática Complutense https://doi.org/10.1007/s13163-020-00382-8, 2020.
  * [17] J. Mao, L. Mu, and Y. Wang, Yet another criterion for the total positivity of Riordan arrays, Linear Algebra Appl. 634 (2022), 106–111.
  * [18] D. Merlini, D. G. Rogers, R. Sprugnoli, and M. C. Verri, On some alternative characterizations of Riordan matrices, Canadian J. Math., 49 (1997), 301–320.
  * [19] A. Pinkus, Totally Positive Matrices, Cambridge University Press, Cambridge, 2010.
  * [20] D. G. Rogers, Pascal triangles, Catalan numbers and renewal matrices, Discrete Math., 22 (1978), 301–310.
  * [21] L. W. Shapiro, Bijections and the Riordan group, Random generation of combinatorial objects and bijective combinatorics, Theoret. Comput. Sci. 307 (2003), no. 2, 403–413.
  * [22] L. W. Shapiro, S. Getu, W. J. Woan and L. Woodson, The Riordan group, Discrete Appl. Math. 34(1991) 229–239.
  * [23] L. W. Shapiro, R. Sprugnoli, P. Barry, G.-S. Cheon, T.-X. He, D. Merlini, and W. Wang, The Riordan Group and Applications, Springer, 2022.
  * [24] R. Slowik, Some (counter)examples on totally positive Riordan arrays, Linear Algebra Appl. 594 (2020), 117–123.


# Residual-based attention and connection to information bottleneck theory in
PINNs

Sokratis J. Anagnostopoulos  Juan Diego Toscano
Nikolaos Stergiopulos  George Em Karniadakis

###### Abstract

Driven by the need for more efficient and seamless integration of physical
models and data, physics-informed neural networks (PINNs) have seen a surge of
interest in recent years. However, ensuring the reliability of their
convergence and accuracy remains a challenge. In this work, we propose an
efficient, gradient-less weighting scheme for PINNs, that accelerates the
convergence of dynamic or static systems. This simple yet effective attention
mechanism is a function of the evolving cumulative residuals and aims to make
the optimizer aware of problematic regions at no extra computational cost or
adversarial learning. We illustrate that this general method consistently
achieves a relative $L^{2}$ error of the order of $10^{-5}$ using standard
optimizers on typical benchmark cases of the literature. Furthermore, by
investigating the evolution of weights during training, we identify two
distinct learning phases reminiscent of the fitting and diffusion phases
proposed by the information bottleneck (IB) theory. Subsequent gradient
analysis supports this hypothesis by aligning the transition from high to low
signal-to-noise ratio (SNR) with the transition from fitting to diffusion
regimes of the adopted weights. This novel correlation between PINNs and IB
theory could open future possibilities for understanding the underlying
mechanisms behind the training and stability of PINNs and, more broadly, of
neural operators.

###### keywords:

Residual-based attention , PINNs convergence accuracy , information bottleneck
theory, self-adaptive weights.

††journal: arXiv

[inst1]organization=Laboratory of Hemodynamics and Cardiovascular Technology,
EPFL,city=Lausanne, postcode=1015, state=VD, country=Switzerland

[inst2]organization=School of Engineering, Brown University,city=Providence,
postcode=02912, state=RI, country=USA [inst3]organization=Division of Applied
Mathematics, Brown University,city=Providence, postcode=02912, state=RI,
country=USA label1label1footnotetext: The first two authors contributed
equally to this work

##  1 Introduction

Physics-informed neural networks (PINNs) offer an alternative to traditional
numerical methods for solving partial differential equations (PDEs) in forward
and inverse problems[1]. In the PINN approach, a neural network is trained to
approximate the solution of a PDE by minimizing a composite loss function that
includes terms related to the physical laws and data from initial boundary
conditions, simulations, or experiments [2]. Keeping a balance between the
loss terms is crucial to avoid convergence and accuracy issues, especially
when dealing with highly nonlinear, multi-scale, or chaotic behavior
problems[3, 4, 5]. To address this imbalance, researchers have developed
several approaches that can be roughly classified into three categories:
neural network modifications, PDE transformations, and adaptive weighting
strategies.

A neural network modification addresses the multi-layer perceptron
optimization capabilities via model reparametrizations [6, 7], input dimension
expansions [8, 9], sequential learning [10, 11], or adaptive activation
functions [12]. On the other hand, PDE transformations attempt to simplify the
problem by enforcing the physical laws in their homologous forms. For
instance, by implementing the PDE on its constrained expression, the number of
losses can be reduced, and the optimizer can focus on a single condition [13,
14, 15]. Similarly, by using alternate or auxiliary physical laws, it is
possible to reduce the number of equations or decrease the PDE order [16, 17],
which eases the optimization process.

Adaptive weighting strategies address the loss imbalances by iteratively
altering the contribution of particular terms or regions of the domain during
the training process. This modification can be done indirectly by adaptively
resampling points in crucial areas [18, 19, 20] or directly by assigning
multipliers that adjust the contribution of each error function. For instance,
Wang et al.,[6] proposed a learning rate annealing algorithm that balances the
losses’ weights based on their back-propagated gradients. Wang et al. [3] also
defined a causal parameter for time-dependent problems that force the model to
learn sequentially based on the time steps. Global multipliers generally
modify the average contribution of each loss term [6, 3, 21, 22]. Similarly,
these multipliers can be applied locally (i.e., per training point); McCLenny
et al. [4] proposed a self-adaptive (SA) method where the loss weights are
trained via adversarial training. Zhang et al. [23] extended this concept and
proposed a differentiable adversarial self-adaptive (DASA) pointwise weighting
scheme that uses a subnetwork to find the optimal weights. These weighting
strategies enhance the capabilities of PINNs and enable them to address
challenging problems in diverse fields. Nevertheless, their implementation can
be expensive since they require coupling auxiliary networks [23], training
additional multipliers [4], or computing and processing gradients [6].
Moreover, some of these approaches grow unboundedly [4, 17] and do not relate
to the neural network training capabilities, which may result in over-refining
some regions while ignoring others, as pointed out by [23].

To address these problems, we propose a bounded residual-based attention (RBA)
scheme that efficiently computes adaptive weights for the collocation points
based on the residuals of the given PDEs. The suggested weighting approach
considers the neural network training dynamics and adapts itself during the
training process. Furthermore, as the optimizer significantly exceeds the
accuracy threshold of the vanilla PINN, we observe two distinct learning
phases, evident in the evolution of weight distributions. We perform
subsequent gradient analysis, and we associate this behavior with the
information bottleneck (IB) theory, which proposes an optimal method to devise
a condensed representation of an input while preserving the majority of
information related to an output [24, 25]. This approach employs the concept
of mutual information, which measures the knowledge gained about one random
variable through the observation of another. According to the IB theory, a
well-functioning model should retain essential output information while
discarding insignificant input details, thereby creating an “information
bottleneck” [26, 27, 28].

This paper is organized as follows. In Section 2, we provide an overview of
PINNs, introduce the RBA weighting scheme, and review additional modifications
that enhance the model performance. Then, in Section 3, we present the
implementation of our method through two benchmark problems and conduct an
ablation study to interpret the relative contribution of each modification. In
Section 4, we analyze the evolution of our multipliers during training and
link their behavior with the bottleneck theory. Finally, in Section 5, we
summarize the findings and discuss our future work. In the appendix, we
include details regarding the implementation.

##  2 Methods

###  2.1 Physics-Informed Neural networks

In the context of Physics-Informed Neural Networks (PINNs), a general partial
differential equation (PDE) can be expressed as:

$\mathcal{L}\\{u(x,t)\\}=f(x,t),$ (1)

where $u(x,t)$ is the unknown function we wish to approximate, $\mathcal{L}$
denotes the differential operator, and $f(x,t)$ is a given source or forcing
function that introduces external influences to the system. The differential
operator $\mathcal{L}$ depends on the specific PDE under consideration. For
instance, in the case of the heat equation,
$\mathcal{L}=\frac{\partial}{\partial t}-\kappa\nabla^{2}$, where $\kappa$ is
the thermal diffusivity and $\nabla^{2}$ is the Laplace operator.

Additionally, the problem may be constrained by several types of boundary
conditions. These conditions are typically categorized into Dirichlet,
Neumann, Robin, or Periodic, and are represented as follows:

3 &u(x_b,t) = g_d(x_b,t), (Dirichlet)
∂u∂n(x_b,t) = g_n(x_b,t), (Neumann)
αu(x_b,t) + β∂u∂n(x_b,t) = g_r(x_b,t), (Robin)
u(x+L,t) = u(x,t), (Periodic)

where $x_{b}$ denotes the boundary points. The Dirichlet condition defines the
function value at the boundary, the Neumann condition specifies the normal
derivative at the boundary, and the Robin condition is a weighted combination
of the function and its normal derivative. For the Periodic condition, the
function value repeats after a certain length $L$.

The loss function in PINNs is designed to encompass the deviation of the
neural network solution from the initial conditions, boundary conditions, and
the PDE itself. The influence of these components in the loss function is
balanced using Lagrange multipliers $\lambda_{r}$, $\lambda_{ic}$ and
$\lambda_{bc}$:

$\mathcal{L}=\lambda_{ic}\cdot\mathcal{L}_{ic}+\lambda_{bc}\cdot\mathcal{L}_{bc}+\lambda_{r}\cdot\mathcal{L}_{r},$
(2)

where $\mathcal{L}_{ic}$, $\mathcal{L}_{bc}$, and $\mathcal{L}_{r}$ represent
the losses associated with the initial conditions, boundary conditions, and
the residuals of the PDE, respectively. In general, each loss term
$\mathcal{L}_{i}$ takes the form:

$\mathcal{L}_{i}=\langle h_{i}(u_{NN})\rangle,\ \ for\ \ i\in\\{ic,bc,r\\}$
(3)

where $\langle\cdot\rangle$ is the mean operator and $h_{i}$ is a function of
the neural network $u_{NN}$, associated with the error between the prediction
and the ground truth.

Each of these terms encourages the neural network to conform to the respective
physical laws or initial/boundary conditions encoded in the problem. The
Lagrange multipliers can be global (e.g., scalars balancing each individual
loss term) or local (e.g., vectors balancing each collocation or boundary
point). In the most basic PINN definition, $\lambda=1$ for all loss terms.

The training of the neural network weights is performed using an optimization
method like stochastic gradient descent (SGD). In each iteration $k$, an
update is applied to a weight $w$ as follows:

$w^{k+1}=w^{k}-\eta\cdot\nabla_{w}\mathcal{L},$ (4)

where $\eta$ is the learning rate, and $\nabla_{w}\mathcal{L}$ is the gradient
of the loss function with respect to $w$. This process iteratively reduces the
loss, pushing the solution towards satisfying the given PDE and
boundary/initial conditions.

###  2.2 Residual-based attention (RBA) scheme

One of the inherent challenges in obtaining accurate results with PINNs is
that the residuals of key collocation points can get overlooked by the mean
calculation of the objective function (Eq. 3). Consequently, despite a
decrease in total loss during training, certain spatial or temporal
characteristics might not be fully captured. This issue becomes particularly
pronounced in multiscale problems, where such troublesome areas could not only
result in a lack of detail but could also impair the flow of important
information from the initial and boundary conditions into the domain of
interest, thereby hindering convergence. Although mini-batching techniques
present a potential solution to this issue, the difficult task of a priori
identifying the problematic regions remains.

A workaround for this problem is selecting a set of Lagrange multipliers,
either global or local, which can be adjustable during training. The aim of
the global multipliers is to balance different terms of the loss function,
while local multipliers aim at weighting the influence of specific collocation
points of the domain. Establishing a systematic rule for updating these
multipliers can significantly aid the optimizer, with several successful
strategies documented in the literature, demonstrating major improvements to
the vanilla PINNs performance [3, 4, 23].

Furthermore, as opposed to classical numerical methods, which can guarantee
stability, the training convergence of PINNs is affected by the evolution of
residuals which will have some degree of stochasticity, leading to a
corresponding degree of instability. Driven by the aforementioned challenges,
the initial objective of this work is the development of a simple, gradient-
less weighting scheme based on the rolling history of the cumulative residuals
to extend the attention span of the optimizer. The proposed scheme aims to
increase attention to the challenging regions of information propagation in
both space and time dimensions defining the problem.

The update rule for the proposed residual-based Lagrange multipliers for any
collocation point $i$ on iteration $k$ is given by:

$\lambda_{i}^{k+1}\leftarrow\gamma\lambda_{i}^{k}+\eta^{*}\frac{|r_{i}|}{\max_{i}({|r_{i}|})},\
\ i\in\\{0,1,...,N_{r}\\},$ (5)

where $N_{r}$ is the number of collocation points, $r_{i}$ is the PDE residual
for point $i$, $\gamma$ is the decay parameter and $\eta^{*}$ is the learning
rate, respectively. Note that the learning rate of the optimizer $\eta$ and
the learning rate of the weighting scheme $\eta^{*}$ are two different
hyperparameters.

This is a convergent linear homogeneous recurrence relation. Given that
$\frac{|r_{i}|}{\max_{i}({|r_{i}|})}\in[0,1]$ and $\lambda_{i}^{0}\neq
0,\forall i\in\\{0,1,...,N_{r}\\}$, its bounds are given by:

$\lambda_{i}^{k}\in(0,\frac{\eta^{*}}{1-\gamma}],\ \ as\ k\rightarrow\infty.$
(6)

These local multipliers scale with the normalized cumulative residuals, but
they also behave dynamically during training due to the decay $\gamma$ as
follows: since $\gamma$ lies between 0 and 1, $\lambda^{k}$ will gradually
decrease in its contribution to $\lambda^{k+1}$, as k increases.
Simultaneously, $\eta^{*}\frac{|r_{i}|}{\max{(|r_{j}|)}}$ is being added in
each iteration, leading to an eventual equilibrium. However, as different
points $i$ may be targeted depending on the stage of the training process, the
distribution of $\lambda_{i}$ will dynamically vary across the domain of
interest.

The main advantages of the proposed residual weighting scheme are summarized
as follows:

  1. 1.

Deterministic operation scheme, bounded by the $\gamma$ and $\eta^{*}$
parameters, which ensure the absence of exploding multipliers.

  2. 2.

No training or gradient calculation is involved, leading to negligible
additional computational cost.

  3. 3.

Scaling with the cumulative residuals guarantees increased attention on the
solution fronts where the PDEs are unsatisfied in both spatial and temporal
dimensions.

The modified training process for an indicative standard Gradient Descent
optimization method is outlined in Algorithm 1. We define the RBA weights as
$\lambda_{i}$ and initialize them with 0 values for all collocation points,
noting that $\lambda_{i}=\sqrt{\lambda_{r}}$ (from Eq. 2). For all training
examples, the RBA learning rate and decay hyperparameters are set as
$\eta^{*}=0.01$ and $\gamma=0.999$, respectively. A constant $c$ may be added
to the updated weights to raise the lower bound and adjust the ratio
$max(\lambda_{i})/min(\lambda_{i})$ (for this study $c=0$).

Data: Training set, learning rates $\eta$, $\eta^{*}$, decay $\gamma$, lower
bound $c$

Result: Optimized parameters $w_{i},\ \ i\in\\{0,1,...,N_{r}\\}$

Initialization: Xavier [29] $\forall$ $w_{i}$ and 0 $\forall$ $\lambda_{i}$;

while  _stopping condition not met_ do

for  _each batch of examples in the training set_ do

$\lambda_{i}^{k+1}\leftarrow\gamma\lambda_{i}^{k}+\eta^{*}\frac{|r_{i}|}{\max_{i}({|r_{i}|})}$;

$\lambda_{i}^{k+1}\leftarrow\lambda_{i}^{k+1}+c$;

$\mathcal{L}_{r}^{*}=\langle(\lambda_{i}\cdot r_{i})^{2}\rangle$;

$\mathcal{L}=\lambda_{ic}\mathcal{L}_{ic}+\lambda_{bc}\mathcal{L}_{bc}+\mathcal{L}_{r}^{*}$;

$w_{i}^{k+1}=w_{i}^{k}-\eta\cdot\nabla_{w}\mathcal{L}$;

end for

end while

Algorithm 1 Residual-based attention Gradient Descent

###  2.3 Additional PINN enhancements

####  2.3.1 Modified multi-layer perceptrons

For the design of the network architecture, we leverage the recently
introduced modified multi-layer perceptron (mMLP) [6], which draws its
inspiration from the attention mechanisms prominent in transformers [30]. The
mMLP aims at augmenting the efficacy of PINNs by embedding the input variables
$x$ into the hidden states of the network. Initially, these inputs are encoded
in a feature space by employing two distinct encoders, $U$ and $V$. They are
then assimilated within each hidden layer of a conventional MLP by point-wise
multiplication. The encoders are given by:

$U=\sigma(x^{0}w^{U}+b^{U}),\quad V=\sigma(x^{0}w^{V}+b^{V})$ (7)

We implement a variation of this concept by performing the encoding operation
right before each activation function, as it performed better for our
experiments. Thus, each forward pass is defined as:

$\alpha^{l}(x)=\alpha^{l-1}(x)w^{l}+b^{l},\quad for\ \
l\in\\{1,2,...,layers\\}$ (8)

$\alpha^{l}(x)=(1-\alpha^{l}(x))\odot U+\alpha^{l}(x)\odot V$ (9)

$\alpha^{l}(x)=\sigma(\alpha^{l}(x))$ (10)

where $x$ is the input, $\alpha^{l}$ and $w^{l}$ are the neurons and weights
of layer $l$, $\sigma$ is the activation function and $\odot$ is the element-
wise product.

####  2.3.2 Exact imposition of boundary conditions

The inexact imposition of boundary conditions can negatively affect the
accuracy and training stability of neural networks [31, 6, 32].

##### Dirichlet Boundary Conditions

A recent study by Sukumar et al. [13] showed how to impose Dirichlet, Neuman,
and Robin boundary conditions in PINNs by using approximate distance functions
(ADF). The constrained expression for Dirichlet boundary conditions is defined
in equation 11.

$u(\mathbf{x})=g(\mathbf{x})+\phi(\mathbf{x})u_{NN}(\mathbf{x}),$ (11)

where $g(\mathbf{x})$ is a function that satisfies $u$ along the boundaries,
$u_{NN}(\mathbf{x})$ is the output of our neural network, and
$\phi(\mathbf{x})$ is a composite distance function that equals zero when
evaluated on the boundaries. If the boundary is composed of $M$ partitions
$[S_{1},...S_{M}]$, the composite distance function for Dirichlet boundary
conditions can be described as:

ϕ(ϕ_1,ϕ_2,…,ϕ_M)=Π_i^M ϕ_i,

where $[\phi_{1},..,\phi_{M}]$ are the individual distance functions. Notice
that if $\mathbf{x}\in S_{i}$, then $\phi(\mathbf{x})=0$ and the neural
network approximation exactly satisfies the boundary conditions, i.e.,
$u(\mathbf{x})=g(\mathbf{x})$.

##### Periodic Boundary Conditions

Similarly, it is possible to enforce periodic boundary conditions as hard
constraints strictly by constructing Fourier feature embeddings of the input
data [3, 31]. In particular, the periodic constraint of a smooth function
$u(x)$ can be encoded in a neural network using a one-dimensional Fourier
feature embedding (See equation 12).

$v(x,y)=[1,cos(\omega_{x}x),sin(\omega_{x}x),...,cos(m\omega_{x}x),sin(m\omega_{x}x)]$
(12)

In the same way, a neural network can incorporate 2D periodic constraints by
utilizing a two-dimensional Fourier feature embedding [3].

$p(x,y)=\bmatrix
cos(\omega_{x}x)cos(\omega_{y}y),...,cos(n\omega_{x}x)cos(m\omega_{y}y)\\\
cos(\omega_{x}x)sin(\omega_{y}y),...,cos(n\omega_{x}x)sin(m\omega_{y}y)\\\
sin(\omega_{x}x)cos(\omega_{y}y),...,sin(n\omega_{x}x)cos(m\omega_{y}y)\\\
sin(\omega_{x}x)sin(\omega_{y}y),...,sin(n\omega_{x}x)sin(m\omega_{y}y)\\\ $
(13)

where, $w_{x}=\frac{2\pi}{P_{x}}$, $w_{y}=\frac{2\pi}{P_{y}}$, $m,n$ are a
non-negative integers, and $P_{x}$ and $P_{y}$ are the periods in the $x$ and
$y$ directions. Dong et al. [31] proved that any network representation
$u_{NN}(v(x))$ and $u_{NN}(p(x,y))$ are periodic in the $x$ and $x,y$
directions, respectively.

##  3 Results

###  3.1 Dynamic case: 1D Allen-Cahn equation

The Allen-Cahn equation is a compelling benchmark for PINNs due to its
challenging characteristics. It is regarded as a “stiff” PDE, a term that
describes equations that require careful numerical handling to avoid
instability in their solutions. This complexity arises from its ability to
produce solutions with pronounced, sharp transitions in both spatial and
temporal dimensions. The 1D Allen-Cahn PDE is given by:

$\frac{\partial u}{\partial t}=10^{-4}\frac{\partial^{2}u}{\partial
x^{2}}-5u^{3}+5u,$ (14)

with initial condition and periodic boundary conditions:

$u(0,x)=x^{2}\cos(\pi x),\quad\forall x\in[-1,1],$ (15)
$u(t,x+1)=u(t,x-1),\quad\forall t\geq 0\quad\text{and}\quad x\in[-1,1].$ (16)

To obtain the solution of Allen-Cahn PDE in Figure 1, we train a PINN for
$3\cdot 10^{5}$ iterations with the standard Adam optimizer [33] and an
exponential learning rate scheduler on 25600 collocation points. For the
periodic boundary conditions, we use the formulation of section 2.3.2 and set
$\lambda_{ic}=100$, following the benchmark setup of the recent literature
[3]. The rest of the parameters are mentioned in the appendix. The benchmark
comparison of the final relative $L^{2}$ error obtained by different methods
is presented in Table 1, where the error for the RBA weights represents the
average of five runs with random seeds. It should be noted that only the best
performing methods of Table 1 ([3], [23]) utilize the PINN enhancements
discussed in section 2.3. The relative $L^{2}$ is defined as:

$relative\ \ L^{2}=\frac{\lVert Exact-Predicted\rVert_{2}}{\lVert
Exact\rVert_{2}}$ (17)

Figure 1: Exact solution of the 1D Allen-Cahn with the corresponding network prediction and the absolute error difference. Method | Relative $L^{2}$ error
---|---
Vanilla PINN [1] | 4.98e - 1
Self-adaptive weights [4] | 2.1e - 2
Time marching [34] | 1.68e - 2
Causal training [3] | 1.39e - 4
DASA-PINNs [23] | 8.57e - 5
RBA weights | 5.7e - 5
Table 1: Relative $L^{2}$ comparison between different methods of the
literature for the 1D Allen-Cahn equation. The RBA run represents the average
of five randomly chosen seeds.

####  3.1.1 Ablation Study for Allen-Cahn

An ablation study involves measuring the performance of a system after
removing one or more of its components to measure the relative contribution of
the ablated components [35]. To investigate the effect of the adopted methods
on the relative $L^{2}$ convergence, we perform a single-component ablation
study and report the results in Table 2. The best-performing model of Table 1
is selected to represent the full model. For each run, a component is removed
to identify the effect of the omitted method on the performance. Figure 2
shows the convergence histories of each run. As evident from these results,
the coupling of the RBA scheme along with the Fourier feature embedding is the
most important component to achieving low relative $L^{2}$. Moreover, the RBA
weights initiate a steep convergence trajectory at 12000 iterations, achieving
a final $L^{2}$ of $3.16\cdot 10^{-3}$ without any additional component.
Combined with the modified MLP and the Fourier features, the full model leads
to the best $L^{2}$ at $4.5\cdot 10^{-5}$. The adopted mMLP accelerates the
convergence in the first 60000 iterations, while it does not provide
significant accuracy gains in the long run. The noise associated with the
best-performing methods indicates that the optimizer is “jumping” over local
minima of the loss landscape and effectively not getting stuck in a bad
solution.

Method | Relative $L^{2}$ error
---|---
RBA | 3.16 - 3
RBA+mMLP | 3.79e - 3
mMLP+Fourier | 3.2e - 4
RBA+Fourier | 5.19e - 5
RBA+mMLP+Fourier | 4.55e - 5
Table 2: Final relative $L^{2}$ error for each case of the single-component
ablation study on the 1D Allen-Cahn. Figure 2: Ablation study convergence for
the 1D Allen-Cahn: Progression of convergence for each experiment. The results
clearly demonstrate that the integration of the RBA approach and the Fourier
feature embedding is crucial for attaining a minimal relative $L^{2}$. When
combined with a modified MLP architecture, an optimal $L^{2}$ of $4.5\cdot
10^{-5}$ is reached. The implemented mMLP aids in speeding up the convergence
during the initial 60000 iterations. Additionally, the RBA weights address the
instability issues in the standard PINN, achieving a relative $L^{2}$ of
$3.16\cdot 10^{-3}$ without supplementary components. The noise linked with
the top-performing methods suggests that the optimization process successfully
avoids becoming trapped in sub-optimal solutions by effectively “leaping” past
local minima in the loss landscape.

###  3.2 Static case: 2D Helmholtz equation

The Helmholtz equation is commonly employed to model wave and diffusion
phenomena, capturing changes over time in either a spatial or combined
spatial-temporal domain. The 2D Helmholtz PDE is expressed as follows:

$\frac{\partial^{2}u}{\partial x^{2}}+\frac{\partial^{2}u}{\partial
y^{2}}+k^{2}u-q(x,y)=0$ (18)

$\split q(x,y)=&-(a_{1}\pi)^{2}sin(a_{1}\pi x)sin(a_{2}\pi y)\\\
-(a_{2}\pi)^{2}sin(a_{1}\pi x)sin(a_{2}\pi y)\\\ +ksin(a_{1}\pi x)sin(a_{2}\pi
y),$ (19)

where $q(x,y)$ is the forcing term which leads to the analytical solution
$u(x,y)=sin(a_{1}\pi x)sin(a_{2}\pi y)$ [4]. The boundary conditions are
defined as follows:

$u(-1,y)=u(1,y)=u(x,-1)=u(x,-1)=0$ (20)

To allow a direct comparison with previous studies, we set $a_{1}=1$,
$a_{2}=4$, $k=1$ and follow the experimental setup described in [23], setting
$\lambda_{bc}=100$. We initially apply a combination of Dirichlet and Periodic
boundary conditions, with the latter enforced using the Fourier feature
embedding. However, to avoid encoding the PDE analytical solution (i.e.,
$\sin(\pi x)\sin(4\pi y)$), we expand our inputs (i.e., $x,y$) using a
truncated embedding with $m=n=5$ (Eq. 21):

$p(x,y)=\bmatrix
cos(\omega_{x}x)cos(\omega_{y}y),...,cos(n\omega_{x}x)cos(n\omega_{y}y)\\\
cos(\omega_{x}x)sin(\omega_{y}y),...,cos(n\omega_{x}x)sin(n\omega_{y}y)\\\
sin(\omega_{x}x)cos(\omega_{y}y),...,sin(n\omega_{x}x)cos(n\omega_{y}y)\\\ $
(21)

Then, we approximate the solution of the PDE as $u_{F}=u_{NN}(p(x,y))$, where
$u_{NN}$ is the modified MLP. Additionally, to demonstrate the flexibility of
RBA weights, we examine an additional case applying only Dirichlet boundaries
with approximate distance functions (ADF). Towards this end, we split the BCs
into four partitions, i.e., $y=1$, $y=-1$, $x=1$, and $x=-1$, and use
$\phi=(x^{2}-1)(y^{2}-1)$ as the composite distance function. Based on
equation 20 we define $g(\mathbf{x})=0$ and approximate the solution of the
Helmholtz equation as,

$u_{ADF}=(x^{2}-1)(y^{2}-1)u_{NN}(x,y),$ (22)

which exactly satisfies the BCs. We train our models (i.e., $u_{F}$ and
$u_{ADF}$) using the Adam optimizer for $2\cdot 10^{4}$ iterations, followed
by L-BFGS for $3\cdot 10^{3}$ iterations, noting that the RBA weights can be
deployed for both optimizers. Figure 3 illustrates the exact solution of the
1x4 Helmholtz PDE alongside the corresponding prediction obtained from the
PINN. Furthermore, in Table 3, we compare the average $L^{2}$ errors for the
Fourier and ADF exact boundary imposition methods with the current state-of-
the-art methods.

Figure 3: Analytical solution for the 2D Helmholtz equation and the corresponding network prediction with the absolute error difference. Method | Relative $L^{2}$ error
---|---
Vanilla PINN [6] | 8.14e - 2
Self-adaptive weights [4] | 3.20e - 3
Learning-rate annealing [6] | 1.49e - 3
DASA-PINNs [23] | 5.35e - 5
RBA weights (ADF) | 8.04e - 5
RBA weights (Fourier) | 1.46e - 5
Table 3: Relative $L^{2}$ comparison between different methods of the
literature for the 2D Helmholtz equation. The RBA run represents the average
of five randomly chosen seeds.

####  3.2.1 Ablation Study for Helmholtz

In this ablation study for the 2D Helmholtz, we present the convergence
histories for each ablated system (Fig. 2) and summarize the final $L^{2}$
errors in Table 5 and 5. For the first case (i.e., $u_{F}$), the Fourier
feature embedding has the most significant influence in reaching a minimized
relative $L^{2}$. The accuracy further improves paired with either the RBA
scheme or the mMLP. The best-performing seed leads to an $L^{2}$ of $8.91\cdot
10^{-6}$; however, we note that this method is employed for comparison
purposes since the Fourier features incorporate the trigonometric traits of
the analytical solution. Therefore, we also perform the same study using
distance functions ($u_{ADF}$), which indicates that the essential component
is the ADF (Table 5). By exactly imposing boundary conditions, the number of
loss terms is reduced, allowing the optimizer to focus on minimizing a single
error function.

Figure 4: Ablation study convergence for the 2D Helmholtz: The convergence
trajectory for each experiment is illustrated here. For this case, the Fourier
feature embedding is critical in achieving a low relative $L^{2}$. Coupled
with the RBA weights and a modified MLP architecture, the model reaches an
optimal $L^{2}$ of $8.91\cdot 10^{-6}$ after $2\cdot 10^{4}$ Adam and $3\cdot
10^{4}$ L-BFGS training steps.

Method | Relative $L^{2}$ error
---|---
RBA | 1.25e - 3
RBA+mMLP | 7.10e - 4
mMLP+Fourier | 1.63e - 5
RBA+Fourier | 3.04e - 5
RBA+mMLP+Fourier | 8.91e - 6

Table 4: Final relative $L^{2}$ error for each case of the ablation study on
the 2D Helmholtz. The exact symmetry boundary condition is enforced using the
Fourier feature embedding.

Method | Relative $L^{2}$ error
---|---
RBA | 1.33e - 3
RBA+mMLP | 7.40e - 4
mMLP+ADF | 6.60e - 5
RBA+ADF | 3.36e - 4
RBA+mMLP+ADF | 5.33e - 5

Table 5: Final relative $L^{2}$ error for 2D Helmholtz where the exact
Dirichlet boundary condition is enforced using ADF.

###  3.3 RBA weight evolution

Figure 5 shows the evolution of RBA weights for the previously analyzed cases.
For each run, the maximum value is bounded by Eq. 6, while the upper bound for
this study equals $10$. The distribution of weights significantly varies
depending on the stage of training; hence the maximum values fluctuate while
approaching the upper bound. The fluctuation of the mean values is less
pronounced until they finally converge to $\approx 20\%$ of the upper bound.
This behavior indicates that, on average, the total magnitude of weights
remains constant while their distribution varies as the optimizer focuses on
different parts of the domain. Thus, by introducing the decay factor $\gamma$,
we can effectively prevent exploding weight values. At the same time, this
allows the weights to be flexibly distributed, adapting as needed throughout
different stages of the solution process.

Figure 5: Evolution of RBA weights. For each case, the peak value is limited
as per Eq. 6, with this study’s upper bound equaling $10$. Notice that for the
Helmholtz (Fourier) case, the RBA weights were only updated during Adam,
indicating that the model did not fully converge. However, for the ADF the
weights are updated during L-BFGS which pushes the maximum value to the upper
bound after 20000 iterations. In general, the maxima fluctuate until they
approach the upper limit during the final learning stages. On the other hand,
the mean values show a less noticeable fluctuation before eventually
stabilizing at about $\approx 20\%$ of the upper bound. This trend indicates
that, on average, the total magnitude of weights remains constant while their
distribution adjusts dynamically as the optimizer shifts its focus across
different domain sections.

##  4 Information bottleneck (IB) theory

The information bottleneck theory was developed as a framework for
understanding the trade-off between compression and prediction in supervised
learning [24, 25]. It formulates an idealized principle for creating a
compressed, or “bottlenecked,” representation of an input variable that
retains as much information as possible about an output variable. The theory
utilizes the concept of mutual information, which measures the amount of
information obtained about one random variable by observing another. It
proposes that a good model will retain all the relevant information about the
output while discarding irrelevant details about the input, creating an
“information bottleneck.”

###  4.1 Fitting and Diffusion Phases of Learning

Recently, it has been proposed that the process of deep learning can be
divided into two primary phases: the fitting phase and the diffusion phase
[26, 28].

  1. 1.

Fitting Phase: In this phase, the model learns to correctly classify the
training data, gradually reducing the empirical error. The model can memorize
the data, building complex representations that fit the prediction set. This
phase tends to decrease both the training and test errors and is characterized
by a high Signal-to-Noise Ratio (SNR), indicating that the model is capturing
more of the useful signal (important features of the data) compared to the
noise (irrelevant or less meaningful details). In this context, “noise” is the
random or unwanted information that does not contribute to the learning task.

  2. 2.

Diffusion Phase: After the model has learned to fit the data, it continues to
learn, but in a subtly different way. The network weights start to diffuse to
improve the model’s generalization capabilities. The model simplifies the
internal representations of the learned patterns by keeping the important
features and discarding the irrelevant ones, thus effectively reducing the
complexity of the representations. The training error remains the same in this
phase, but the test error can decrease. The SNR also decreases as the model
becomes more focused on preserving the meaningful signals and is less
influenced by the noise. These two distinct phases can be seen as a process
where the model first fits the data (captures relevant information) and then
compresses it (discards irrelevant information), further enhancing its
generalization ability.

###  4.2 Connection of IB theory to PINNs

The benchmark examples demonstrate that the proposed residual-based method can
extend the training accuracy beyond the original PINN formulation. Upon
inspection of the weight evolution during the training process, we observed a
behavior that largely met our expectations; the weights indeed form an ordered
pattern consistent with the flow of information, primarily dictated by the
initial conditions (Figs. 6,7). Yet, a surprising element emerged as we noted
an abrupt shift to chaotic patterns, occurring once the model descended below
a specific test loss threshold (defined by the relative $L^{2}$ error). Hence,
two definite phases, akin to the “fitting” and “diffusion” from the IB theory,
became apparent.

Following the formulation proposed by the IB theory, we performed a gradient
analysis to quantify the transitions from high to low SNR. The normalized
signal, $\lVert\mu\rVert$ and noise, $\lVert\sigma\rVert$ of the stochastic
gradient distributions for each hidden layer $l$ are given by:

$\lVert\mu_{l}\rVert=\frac{\lVert\langle\nabla_{w}\mathcal{L}_{i}\rangle\rVert_{F}}{\lVert
w_{i}\rVert_{F}}$ (23)

$\lVert\sigma_{l}\rVert=\frac{\lVert
std(\nabla_{w}\mathcal{L}_{i})\rVert_{F}}{\lVert w_{i}\rVert_{F}}$ (24)

$SNR=\frac{\lVert\mu_{l}\rVert}{\lVert\sigma_{l}\rVert},$ (25)

where $i$ is the mini-batch and $\lVert\cdot\rVert_{F}$ is the Forbenius norm.
To apply this formulation, we divide our collocation points into 100 mini-
batches and perform gradient calculations per batch for every 50 iterations of
the optimizer without performing a backward pass. Hence, we deduce the batch-
wise gradient information while the models are trained with the full batch.

We present the findings for the indicative 1D Allen-Cahn (Fig. 6) and 2D
Helmholtz with $a_{1}=a_{2}=6$ (Fig. 7) runs. We define the three regimes of
RBA weights evolution as Phase I (fitting), Phase II (transitional), and Phase
III (diffusion). During these phases, the RBA weights undergo a swift change
from order to disorder which is accompanied by a high to low $SNR$ transition
of the layer gradients, respectively. We also plot the normalized Loss and
normalized relative $L^{2}$ history to compare their convergence rates. In
both cases, the convergence history of the normalized Loss and $L^{2}$
indicates that the model continues to learn after the Loss has reached a
plateau by further decreasing the $L^{2}$ error. In addition, the magnitudes
of the network’s weights $\lVert w\rVert$ begin their final convergence
trajectory at the beginning of the diffusive phase. For both examples, a fully
connected MLP is used as it produces more pronounced transitions. It is
important to mention that the convergence histories show the monotonic
decrease of the Loss and $L^{2}$ errors. The animations of the presented
cases, along with the relevant code, will be available on
https://github.com/soanagno/rba-pinns.

Figure 6:  Three phases of learning for Allen-Cahn (top row): Distribution of
RBA weights for Fitting (I), Transitional (II), and Diffusion (III) phases for
the Allen-Cahn problem. The weights undergo a swift change from order to
disorder at $2\cdot 10^{4}$ iterations, where they stop following the flow of
information imposed by the initial condition. Gradient distributions of
network weights (middle row): Corresponding $\lVert\mu\rVert$,
$\lVert\sigma\rVert$ and $SNR$ curves for each layer, which indicate a high
$SNR$ during the fitting phase followed by a low $SNR$ during the diffusion
phase, as proposed by the IB theory. Convergence history (bottom row):
Convergence of normalized Loss/relative $L^{2}$ and network weights $\lVert
w\rVert$. Even though the Loss reaches a plateau after phase II, the relative
$L^{2}$ decreases quickly while the network weights begin their final
convergence trajectory. Notably, the learning continues far beyond the
convergence of the Loss. Figure 7: Three phases of learning for Helmholtz (top
row): Fitting (I), Transitional (II), and Diffusion (III) phases, during which
the RBA weights distribution is presented. This case is run for $10^{5}$
iterations to capture the final learning states fully. After 1400 iterations,
the weights experience a rapid transition from an ordered to a disordered
state, breaking away from the information flow dictated by the boundary
conditions. Gradient distributions of network weights (middle row): Gradient
distributions of the network weights, displaying the corresponding
$\lVert\mu\rVert$, $\lVert\sigma\rVert$ and $SNR$ curves for each layer. This
clearly indicates a high $SNR$ in the fitting phase and a subsequent drop
during the diffusion phase, aligning with the IB theory’s predictions.
Convergence history (bottom row): Convergence of the normalized Loss/relative
$L^{2}$ and network weights $\lVert w\rVert$. Although the normalized Loss and
relative $L^{2}$ reach a plateau relatively early (after 5000 iterations), the
Loss converges first around phase II. This is also concurrent with the final
convergence trajectory onset of the network weights.

##  5 Summary

In this work, we propose a fast, residual-based attention scheme which can
enhance the accuracy of PINNs in both static and dynamic cases, competing with
the state-of-the-art performing models of the literature. Moreover, as the
ablation study points out, a critical aspect of the solution accuracy is the
formulation of exact boundary conditions. Additionally, even though the RBA
multipliers provide considerable accuracy improvements, we believe that an
equally important element of this type of weighting strategy is in the realms
of convergence. For example, the 1D Allen-Cahn equation cannot be solved by
the vanilla PINN due to its stiff nature, but RBA can initiate an early
convergence trajectory. Furthermore, the versatility of the proposed method
allows for its application with a few additional lines of code.

In addition, given that the adopted methodology can overcome a certain
accuracy threshold, we provide evidence that the learning process is dominated
by two distinct learning regimes, as suggested by the IB theory. In the first
fitting phase, the model learns the most important parts of the solution. The
distribution of the RBA weight values follows the main shape of the solution
while it sets the stage for the next learning phase. In the diffusion phase,
the weights stop following the flow of information from the initial conditions
and start focusing on unregulated parts of the domain, resembling a diffusion
process. Based on the IB theory, it is during this final phase that the model
finds the best balance between detail and simplicity by tweaking what was
learned in the fitting phase, hence refining its understanding of the problem.

This novel connection between IB theory and PINNs could open new possibilities
for quantifying the performance of the many versions of PINNs and even neural
operators. Furthermore, in applications where accuracy is essential, utilizing
methods that collectively aid the optimizer in overcoming the threshold of the
fitting phase could be a critical convergence criterion. Planned future work
is focused on expanding the application of RBA weights to a broader spectrum
of problems and testing the sensitivity of the associated hyperparameters.
Additionally, we aim to undertake comprehensive parametric analyses of the
bottleneck phase transition to enhance our understanding of its occurrence
timing, underlying causes, and possible strategies to instigate it.

## Acknowledgements

SJA and NS thank the Swiss National Science Foundation grant “Hemodynamics of
Physiological Aging” (Grant nr $205321\\_197234$). JDT and GEK acknowledge
support by the DOE SEA-CROGS project (DE-SC0023191), the MURI-AFOSR
FA9550-20-1-0358 project, and the ONR Vannevar Bush Faculty Fellowship
(N00014-22-1-2795). Finally, we extend our gratitude to Prof. N. Sukumar from
the University of California, Davis, for engaging in insightful discussions on
the approximate distance functions.

##  Appendix A Implementation Details

The chosen architecture for all benchmarks is an MLP of 6 hidden layers and
128 neurons per layer, using the hyperbolic tangent ($tanh$) activation
function. The weights are initialized with Xavier normal distribution [29].
The rest of the case-specific information and their corresponding training
time is tabulated in the following sections. However, the reported training
times should be used only as an inner reference since the codes were not
optimized for speed.

###  A.1 Allen-Cahn

Implementation Details
---
Sampling method | Latin Hypercube
Number of collocation points | 25600
Number of initial condition points | 512
Initial learning rate | 1e-3
Decay rate | 0.9
Decay step | 5000
ML framework | PyTorch
GPU | GeForce RTX $3090$
Table 6: Allen Cahn Implementation Details Method |  Training time (ms/iteration)
---|---
Vanilla |  17
RBA |  17
RBA+mMLP |  58
RBA+Fourier |  19
mMLP+Fourier |  60
RBA+mMLP+Fourier |  61
Table 7: Comparison of computational speed for different combinations of the
ablation study for the 1D Allen-Cahn equation.

###  A.2 Helmholtz

Implementation Details
---
Sampling method | Uniform grid
Number of collocation points | 1201
Number of boundary points (Fourier) | 200
Initial learning rate | 5e -3
Decay rate | 0.7
Decay step | 1000
ML framework | Jax
GPU | GeForce RTX $3090$
Table 8: Helmholtz ($a_{1}=1$,$a_{2}=4$) Implementation Details Method |  Training time (ms/iteration)
---|---
Vanilla |  13.8
RBA |  13.8
RBA+mMLP |  18.6
RBA+Fourier |  22.6
mMLP+Fourier |  37.8
RBA+mMLP+Fourier |  34.4
RBA+ADF |  14.2
mMLP+ADF |  18.4
RBA+mMLP+ADF |  18.6
Table 9: Training time for 2D Helmholtz ($a_{1}=1$,$a_{2}=4$)

###  A.3 2D Helmholtz-Bottleneck

Implementation Details
---
Number of collocation points | 90400
Initial learning rate | 5e -3
Decay rate | 0.7
Decay step | 1000
ML framework | Jax
GPU | GeForce RTX $3090$
Training time per iteration | 88.4ms
Table 10: Helmholtz ($a_{1}=6$,$a_{2}=6$) Implementation Details

##  Appendix B Snapshots of weights and residuals

This section displays indicative snapshots of the RBA weights and the
corresponding residual evolution for the Allen-Cahn and Helmholtz models. As
stipulated by the information bottleneck theory, during the training process
we can observe the fitting, transition, and diffusion learning phases. The
interplay between weights and residuals is evident as the former keep track of
the short-term residual history. Therefore, the RBA weights may differ from
the network’s residuals during the early iterations but as the latter
stabilize in the final stages of convergence, the weights tend to follow suit,
reflecting the trends observed in the PDE residuals.

Figure 8: Snapshots of weights and residuals for the Allen-Cahn PDE. We show
results up to $3\cdot 10^{5}$ iterations to capture fitting, transition, and
diffusion learning phases. The fitting phase ends around iteration 20000 while
the weights experience a transition from order to disorder, breaking away from
the information flow dictated by the initial condition. After the transition
phase is over, the diffusion state is predominant around i=40000. During this
chaotic phase, the SNR decreases and the model randomly focuses on refining
scattered parts of the domain. Figure 9: Snapshots of weights and residuals
for Helmholtz ($a_{1}=6,a_{2}=6$). Results are shown for a total of $2\cdot
10^{4}$ iterations, capturing the fitting, transition, and diffusion learning
phases. During the first 800 iterations, the fitting phase is dominant and the
weights follow a specific order dictated by the trigonometric shape of the PDE
solution. After that, the transition phase starts and the RBA weights rapidly
switch to a disordered state. From 5000 iterations and onward, the system is
primarily ruled by diffusion, causing the RBA weights to become chaotic and
focus on optimizing unregulated regions of the domain.

## References

  * [1] M. Raissi, P. Perdikaris, G. E. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational physics 378 (2019) 686–707.
  * [2] S. Cai, Z. Mao, Z. Wang, M. Yin, G. E. Karniadakis, Physics-informed neural networks (pinns) for fluid mechanics: A review, Acta Mechanica Sinica 37 (12) (2021) 1727–1738.
  * [3] S. Wang, S. Sankaran, P. Perdikaris, Respecting causality is all you need for training physics-informed neural networks, arXiv preprint arXiv:2203.07404 (2022).
  * [4] L. McClenny, U. Braga-Neto, Self-adaptive physics-informed neural networks using a soft attention mechanism, arXiv preprint arXiv:2009.04544 (2020).
  * [5] V. A. Es’ kin, D. V. Davydov, E. D. Egorova, A. O. Malkhanov, M. A. Akhukov, M. E. Smorkalov, About optimal loss function for training physics-informed neural networks under respecting causality, arXiv preprint arXiv:2304.02282 (2023).
  * [6] S. Wang, Y. Teng, P. Perdikaris, Understanding and mitigating gradient flow pathologies in physics-informed neural networks, SIAM Journal on Scientific Computing 43 (5) (2021) A3055–A3081.
  * [7] T. Salimans, D. P. Kingma, Weight normalization: A simple reparameterization to accelerate training of deep neural networks, Advances in neural information processing systems 29 (2016).
  * [8] W. Guan, K. Yang, Y. Chen, Z. Guan, A dimension-augmented physics-informed neural network (dapinn) with high level accuracy and efficiency, arXiv preprint arXiv:2210.13212 (2022).
  * [9] S. Wang, H. Wang, P. Perdikaris, On the eigenvector bias of fourier feature networks: From regression to solving multi-scale pdes with physics-informed neural networks, Computer Methods in Applied Mechanics and Engineering 384 (2021) 113938.
  * [10] A. Krishnapriyan, A. Gholami, S. Zhe, R. Kirby, M. W. Mahoney, Characterizing possible failure modes in physics-informed neural networks, Advances in Neural Information Processing Systems 34 (2021) 26548–26560.
  * [11] C. L. Wight, J. Zhao, Solving allen-cahn and cahn-hilliard equations using the adaptive physics informed neural networks, arXiv preprint arXiv:2007.04542 (2020).
  * [12] A. D. Jagtap, K. Kawaguchi, G. E. Karniadakis, Adaptive activation functions accelerate convergence in deep and physics-informed neural networks, Journal of Computational Physics 404 (2020) 109136.
  * [13] N. Sukumar, A. Srivastava, Exact imposition of boundary conditions with distance functions in physics-informed deep neural networks, Computer Methods in Applied Mechanics and Engineering 389 (2022) 114333.
  * [14] C. Leake, D. Mortari, Deep theory of functional connections: A new method for estimating the solutions of partial differential equations, Machine learning and knowledge extraction 2 (1) (2020) 37–55.
  * [15] L. Lu, R. Pestourie, W. Yao, Z. Wang, F. Verdugo, S. G. Johnson, Physics-informed neural networks with hard constraints for inverse design, SIAM Journal on Scientific Computing 43 (6) (2021) B1105–B1132.
  * [16] X. Jin, S. Cai, H. Li, G. E. Karniadakis, Nsfnets (navier-stokes flow nets): Physics-informed neural networks for the incompressible navier-stokes equations, Journal of Computational Physics 426 (2021) 109951.
  * [17] S. Basir, Investigating and mitigating failure modes in physics-informed neural networks (pinns), arXiv preprint arXiv:2209.09988 (2022).
  * [18] L. Lu, X. Meng, Z. Mao, G. E. Karniadakis, Deepxde: A deep learning library for solving differential equations, SIAM review 63 (1) (2021) 208–228.
  * [19] C. Wu, M. Zhu, Q. Tan, Y. Kartha, L. Lu, A comprehensive study of non-adaptive and residual-based adaptive sampling for physics-informed neural networks, Computer Methods in Applied Mechanics and Engineering 403 (2023) 115671.
  * [20] B. Zapf, J. Haubner, M. Kuchta, G. Ringstad, P. K. Eide, K.-A. Mardal, Investigating molecular transport in the human brain from mri with physics-informed neural networks, Scientific Reports 12 (1) (2022) 15475.
  * [21] Z. Xiang, W. Peng, X. Liu, W. Yao, Self-adaptive loss balanced physics-informed neural networks, Neurocomputing 496 (2022) 11–34.
  * [22] D. Liu, Y. Wang, A dual-dimer method for training physics-constrained neural networks with minimax architecture, Neural Networks 136 (2021) 112–125.
  * [23] G. Zhang, H. Yang, F. Zhu, Y. Chen, et al., Dasa-pinns: Differentiable adversarial self-adaptive pointwise weighting scheme for physics-informed neural networks, SSRN (2023).
  * [24] N. Tishby, F. C. Pereira, W. Bialek, The information bottleneck method, arXiv preprint physics/0004057 (2000).
  * [25] N. Tishby, N. Zaslavsky, Deep learning and the information bottleneck principle, in: 2015 ieee information theory workshop (itw), IEEE, 2015, pp. 1–5.
  * [26] R. Shwartz-Ziv, N. Tishby, Opening the black box of deep neural networks via information, arXiv preprint arXiv:1703.00810 (2017).
  * [27] A. M. Saxe, Y. Bansal, J. Dapello, M. Advani, A. Kolchinsky, B. D. Tracey, D. D. Cox, On the information bottleneck theory of deep learning, Journal of Statistical Mechanics: Theory and Experiment 2019 (12) (2019) 124020.
  * [28] Z. Goldfeld, Y. Polyanskiy, The information bottleneck problem and its applications in machine learning, IEEE Journal on Selected Areas in Information Theory 1 (1) (2020) 19–38.
  * [29] X. Glorot, Y. Bengio, Understanding the difficulty of training deep feedforward neural networks, in: Proceedings of the thirteenth international conference on artificial intelligence and statistics, JMLR Workshop and Conference Proceedings, 2010, pp. 249–256.
  * [30] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, I. Polosukhin, Attention is all you need, Advances in neural information processing systems 30 (2017).
  * [31] S. Dong, N. Ni, A method for representing periodic functions and enforcing exactly periodic boundary conditions with deep neural networks, Journal of Computational Physics 435 (2021) 110242.
  * [32] J. Chen, R. Du, K. Wu, A comparison study of deep galerkin method and deep ritz method for elliptic problems with different boundary conditions, arXiv preprint arXiv:2005.04554 (2020).
  * [33] D. P. Kingma, J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980 (2014).
  * [34] R. Mattey, S. Ghosh, A novel sequential method to train physics informed neural networks for allen cahn and cahn hilliard equations, Computer Methods in Applied Mechanics and Engineering 390 (2022) 114474.
  * [35] J. D. Toscano, C. Zuniga-Navarrete, W. D. J. Siu, L. J. Segura, H. Sun, Teeth mold point cloud completion via data augmentation and hybrid rl-gan, Journal of Computing and Information Science in Engineering 23 (4) (2023) 041008.