# An Improved Metarounding Algorithm via Frank-Wolfe

Ryotaro Mitsuboshi
Kyushu University/RIKEN AIP
<EMAIL_ADDRESS>
& Kohei Hatano
Kyushu University/RIKEN AIP
<EMAIL_ADDRESS>
& Eiji Takimoto
Kyushu University
<EMAIL_ADDRESS>
###### Abstract

Metarounding is an approach to convert an approximation algorithm for linear
optimization over some combinatorial classes to an online linear optimization
algorithm for the same class. We propose a new metarounding algorithm under a
natural assumption that a relax-based approximation algorithm exists for the
combinatorial class. Our algorithm is much more efficient in both theoretical
and practical aspects.

_Keywords_ Metarounding $\cdot$ Boosting $\cdot$ Frank-Wolfe

##  1 Introduction

Online decision-making problems using combinatorial objects arise in many
applications, e.g., routing, advertisements, and resource allocation tasks.
Examples of combinatorial objects include permutations [14], $k$-sets [26],
paths [12], matching [4, 14], spanning trees [20] and so on. Let
$\mathcal{C}\subset\mathbb{N}^{n}_{+}$ be a combinatorial class, i.e., a set
of combinatorial objects, where each combinatorial object is represented as a
positive vector. More formally, online combinatorial linear optimization over
a combinatorial class $\mathcal{C}$ is defined as a repeated game between a
player and an environment. The protocol is the following: For each round
$t=1,2,\dots$, (i) the player chooses $\bm{c}_{t}\in\mathcal{C}$. Then, (ii)
the environment reveals a loss vector $\bm{\ell}_{t}\in[0,1]^{n}$ and (iii)
the player incurs the loss $\bm{c}_{t}\cdot\bm{\ell}_{t}$.

The player’s goal is to minimize the $\alpha$-regret defined as

$R_{T}(\alpha)=\sum_{t=1}^{T}\bm{c}_{t}\cdot\bm{\ell}_{t}-\alpha\min_{\bm{c}\in\mathcal{C}}\sum_{t=1}^{T}\bm{c}\cdot\bm{\ell}_{t},$

for as small $\alpha\geq 1$ as possible.

There are two approaches for online combinatorial linear optimization over
$\mathcal{C}$. The first approach is to design algorithms for individual
classes [4, 14, 26, 27, 28]. The second approach is so called offline-to-
online conversion, which uses an offline approximation algorithm over
$\mathcal{C}$ to construct an online algorithm over $\mathcal{C}$ [16, 11, 9].
An offline approximation algorithm takes $\bm{\ell}$ as input and then outputs
a vector $\bm{c}\in\mathcal{C}$ satisfying
$\bm{c}\cdot\bm{\ell}\leq\alpha\min_{\bm{c}^{\prime}\in\mathcal{C}}\bm{c}^{\prime}\cdot\bm{\ell}$.
The offline-to-online conversion approach is generic in that once we design an
efficient conversion algorithm, we can just use known approximation algorithms
for different combinatorial classes. For example, Garber proposed an algorithm
that calls the sub-routine $O(\ln T)$ times for each round $t$ and achieves
$R_{T}(\alpha)=O(T^{2/3})$ [11]. The paper also proposed an algorithm that
calls the sub-routine $O(\sqrt{T}\ln T)$ per round and achieves
$R_{T}(\alpha)=O(\sqrt{T})$. The algorithm needs to know the number of total
rounds $T$ and the approximation ratio $\alpha$ a priori. In addition, the
number of sub-routine calls depends on $T$, so the conversion via the sub-
routine may slow for a huge $T$.

Figure 1: A conceptual diagram for online combinatorial linear optimization
via metarounding. In each round $t$, metarounding converts a vector
$\bm{x}_{t}$ into a combinatorial vector $\bm{c}_{t}\in\mathcal{C}$.

Another offline-to-online conversion approach is to use metarounding[2], which
uses a relaxation-based approximation algorithm. In fact, for many
combinatorial classes $\mathcal{C}$, there are relaxation-based approximation
algorithms. Thus, the requirement is natural. Metarounding converts a vector
$\bm{x}\in\mathop{\rm relax}(\mathcal{C})$ into a combinatorial vector
$\bm{c}\in\mathcal{C}$, where $\mathop{\rm relax}(\mathcal{C})$ is a convex
set satisfying $\mathcal{C}\subset\mathop{\rm relax}(\mathcal{C})$. The relax-
based approximation algorithms takes $\bm{\ell}$ as input and then outputs a
vector $\bm{c}\in\mathcal{C}$ satisfying
$\bm{c}\cdot\bm{\ell}\leq\alpha\min_{\bm{p}\in\mathop{\rm
relax}(\mathcal{C})}\bm{p}\cdot\bm{\ell}$. The only difference between
relaxation-based approximation and the approximation algorithm is the right-
hand-side of the condition. Using the relaxation-based approximation algorithm
as sub-routine, metarounding algorithms aim to output a vector
$\bm{c}\in\mathcal{C}$ randomly to satisfy
$\mathbb{E}_{\bm{c}}[\bm{c}]\cdot\bm{\ell}\leq\alpha\bm{x}\cdot\bm{\ell}$ for
all $\bm{\ell}$. Thus, combining the metarounding algorithm with an arbitrary
online linear optimization algorithm over $\mathop{\rm relax}(\mathcal{C})$,
one can achieve $R_{T}(\alpha)=\alpha(1\text{-regret of the OLO algorithm})$
[10]. Thus, once an efficient metarounding algorithm is obtained, we can fuse
it to an OLO algorithm to obtain a better regret bound. Figure 1 shows the
concept of metarounding for online combinatorial linear optimization over
$\mathcal{C}$. Assuming the existence of such relax-based approximation
algorithm, Fujita et al. [10] proposes a metarounding algorithm that finds a
distribution $\bm{\lambda}\in\mathcal{P}^{\mathcal{C}}$ satisfying
$\mathbb{E}_{\bm{c}}[\bm{c}]\cdot\bm{\ell}\leq(\alpha+\epsilon)\bm{x}\cdot\bm{\ell}$
for all $\bm{\ell}$ in $O(L^{2}D_{\infty}^{2}\ln(Ln)/\epsilon^{2})$
iterations, where $D_{\infty}^{2}$ is the maximal value of a combinatorial
vector $\bm{c}\in\mathcal{C}$ and $L$ is a constant not known in advance.
Thus, they need to run their algorithm for multiple $L$ to find a solution;
therefore, the actual complexity bound is unknown.

This paper proposes a new metarounding algorithm that finds a distribution
$\bm{\lambda}\in\mathcal{P}^{\mathcal{C}}$ over $\mathcal{C}$ that satisfies
$\mathbb{E}_{\bm{c}\sim\bm{\lambda}}[\bm{c}]\cdot\bm{\ell}\leq(\alpha+\epsilon)\bm{x}\cdot\bm{\ell}$
for all $\bm{\ell}\in[0,1]^{n}$ in $O(D_{\infty}^{2}\ln(n)/\epsilon^{2})$
rounds. This guarantee is quite better compared to the one by Fujita et al.
[10]. Our algorithm is designed based on Frank-Wolfe algorithms [7, 15]. Table
1 summarizes the previous works and our contribution.

Our technical contribution has two folds. The first contribution is a re-
formulation of the optimization problem for metarounding, based on a new
observation (Proposition 1). This formulation enables us to derive algorithms
with better iteration bounds. The second contribution is a new boosting-based
metarounding algorithm. Our algorithm is similar to ERLPBoost[25], but ours
employs a modified regularizer which is better suited for metarounding. In
fact, it can be viewed that our algorithm is a generalization of ERLPBoost for
the approximation factor $\alpha>1$. Furthermore, our analyses are based on
those for Frank-Wolfe, which are completely different from those for
ERLPBoost. As a result, our result significantly improves Fujita et al’s on
the number of approximation oracle calls per trial from
$O(L^{2}\mathcal{C}^{2}\ln(Ln)/\epsilon^{2})$ to
$O(\mathcal{C}^{2}\ln(n)/\epsilon^{2})$. We also observe its significant
improvements in the preliminary experiments.

Table 1: Comparison of the previous works. Our assumption is the same as the one in Fujita et al. [10]. Algorithm |  | Approx. alg.
---
oracle
| # of oracle calls
---
per trial
Regret
Kalai et al. [16] |  Any |  $O(T)$ | $O(T^{1/2})$
Garber [11] |  Any |  $O(n^{2}\ln T)$ | $O(T^{2/3})$
Garber [11] |  Any |  $O(n^{2}\sqrt{T}\ln T)$ | $O(T^{1/2})$
Fujita et al. [10] |  relax-based |  $O(L^{2}D_{\infty}^{2}\ln(Ln)/\epsilon^{2})$ | $O(T^{1/2})$
This work |  relax-based |  $O(D_{\infty}^{2}\ln(n)/\epsilon^{2})$ | $O(T^{1/2})$

##  2 Preliminary

Throughout this paper, we use $\mathcal{C}\subset\mathbb{N}^{n}$ to represent
the combinatorial set of our interest and let
$D_{\infty}=\max_{\bm{c}\in\mathcal{C}}\|\bm{c}\|_{\infty}$ as the maximal
value of the entry of $\bm{c}$. We use the notion $\bm{e}_{\bm{c}}$ to denote
the canonical basis vector over $\mathbb{R}^{\mathcal{C}}$. We denote
$\mathcal{P}^{n}=\\{\bm{p}\in[0,1]^{n}\mid\|\bm{p}\|_{1}=1\\}$ as the
probability simplex over $\mathbb{R}^{n}$. For abbreviation, we often writes
$\mathcal{P}^{\mathcal{C}}$ to represent the probability simplex over a set
$\mathcal{C}$. Given a set $\mathcal{C}\subset\\{0,1\\}^{n}$, To introduce
metarounding, We also define the barrier function
$I_{\mathcal{C}}:\mathbb{R}^{n}\to\\{0,\infty\\}$ over a set
$\mathcal{C}\subset\mathbb{R}^{n}$ such that $\bm{c}\in\mathcal{C}$ iff
$I_{\mathcal{C}}(\bm{c})=0$.

######  Definition 1 (Relax-based $\alpha$-approximation algorithm).

Algorithm $\mathcal{B}$ is said to be a relax-based $\alpha$-approximation
algorithm for $\mathcal{C}$ if $\mathcal{B}$ takes
$\bm{\ell}\in\mathbb{R}_{+}^{n}$ as input and then outputs a vector
$\bm{c}\in\mathcal{C}$ satisfying

$\displaystyle\bm{c}\cdot\bm{\ell}\leq\alpha\min_{\bm{p}\in\mathop{\rm
relax}(\mathcal{C})}\bm{p}\cdot\bm{\ell},$

where $\mathop{\rm relax}(\mathcal{C})$ is a convex set sastisfying
$\mathcal{C}\subset\mathop{\rm relax}(\mathcal{C})$.

######  Assumption 1.

We assume the existence of a polynomial time relax-based
$\alpha$-approximation algorithm $\mathcal{B}$ for $\mathcal{C}$ with a
relaxation set $\mathop{\rm relax}(\mathcal{C})$.

######  Definition 2 ($\alpha$-metarounding).

Algorithm $\mathcal{A}$ is said to be an $\alpha$-metarounding algorithm for
$(\mathcal{C},\mathop{\rm relax}(\mathcal{C}))$ if $\mathcal{A}$, when given
$\bm{x}\in\mathop{\rm relax}(\mathcal{C})$ as input, outputs a vector
$\bm{c}\in\mathcal{C}$ randomly such that for all
$\bm{\ell}\in\mathbb{R}_{+}^{n}$,
$\mathbb{E}_{\bm{c}}[\bm{c}\cdot\bm{\ell}]\leq\alpha\bm{x}\cdot\bm{\ell}$.

The existence of a relax-based $\alpha$-approximation algorithm implies the
existence of an $\alpha$-metarounding algorithm [3].

We say that a function $f$ is said to be $\eta$-strongly convex w.r.t. the
norm $\|\cdot\|$ if

$\displaystyle\forall\bm{x},\bm{y}\in\mathop{\rm dom}f,\quad f(\bm{y})\geq
f(\bm{x})+(\bm{y}-\bm{x})\cdot\nabla
f(\bm{x})+\frac{\eta}{2}\|\bm{y}-\bm{x}\|^{2}.$

Similarly, we say that a function $f$ is said to be $\zeta$-smooth w.r.t. the
norm $\|\cdot\|$ if

$\displaystyle\forall\bm{x},\bm{y}\in\mathop{\rm dom}f,\quad f(\bm{y})\leq
f(\bm{x})+(\bm{y}-\bm{x})\cdot\nabla
f(\bm{x})+\frac{\zeta}{2}\|\bm{y}-\bm{x}\|^{2}.$

Our analysis uses the Fenchel duality, so we introduce the Fenchel conjugate
function.

######  Definition 3 (Fenchel conjugate).

The Fenchel conjugate $f^{\star}$ of a function
$f:\mathbb{R}^{k}\to[-\infty,+\infty]$ is defined as

$\displaystyle
f^{\star}(\bm{\theta})=\sup_{\bm{\mu}\in\mathbb{R}^{k}}\left(\bm{\mu}\cdot\bm{\theta}-f(\bm{\mu})\right).$

The Fenchel conjugate of the max function $\bm{x}\mapsto\max_{i}x_{i}$ is a
barrier function of the probability simplex.

$\displaystyle f(\bm{x})=\max_{i\in[k]}x_{i}\quad\iff\quad
f^{\star}(\bm{y})=I_{\mathcal{P}^{k}}(\bm{y})=\begin{cases}0&\bm{y}\in\mathcal{P}^{k},\\\
+\infty&\bm{y}\notin\mathcal{P}^{k}.\end{cases}$

It is well known that if $f$ is a $1/\eta$-strongly convex function w.r.t. the
norm $\|\cdot\|$ for some $\eta>0$, $f^{\star}$ is an $\eta$-smooth function
w.r.t. the dual norm $\|\cdot\|_{\star}$. Further, if $f$ is a strongly convex
function, the gradient vector of $f^{\star}$ is written as $\nabla
f^{\star}(\bm{\theta})=\arg\sup_{\bm{\mu}\in\mathbb{R}^{k}}\left(\bm{\mu}\cdot\bm{\theta}-f(\bm{\mu})\right)$.
One can find the proof of these properties here [1, 22]. Now, we show the
duality theorem [1].

######  Theorem 1 (Fenchel duality theorem [1]).

Let $f:\mathbb{R}^{k}\to(-\infty,+\infty]$ and
$g:\mathbb{R}^{n}\to(-\infty,+\infty]$ be convex functions, and a linear map
$A:\mathbb{R}^{n}\to\mathbb{R}^{k}$. Define the Fenchel problems

$\displaystyle\gamma=\inf_{\bm{d}}f(\bm{d})+g(A^{\top}\bm{d}),\qquad\rho=\sup_{\bm{w}}-f^{\star}(-A\bm{w})-g^{\star}(\bm{w}).$
(1)

Then, $\gamma\geq\rho$ holds. Further, $\gamma=\rho$ holds if111For a convex
set $\mathcal{S}\subset\mathbb{R}^{n}$, $\mathop{\rm
int}(\mathcal{S})=\\{\bm{w}\in\mathcal{S}\mid\forall\bm{v}\in\mathbb{R}^{n},\exists
t>0,\forall\tau\in[0,t],\bm{w}+\tau\bm{v}\in\mathcal{S}\\}$ and $\mathop{\rm
dom}g-A^{\top}\mathop{\rm
dom}f=\\{\bm{w}-A^{\top}\bm{d}\mid\bm{w}\in\mathop{\rm
dom}g,\bm{d}\in\mathop{\rm dom}f\\}$.  $\bm{0}\in\mathop{\rm
int}\left(\mathop{\rm dom}g-A^{\top}\mathop{\rm dom}f\right)$. Furthermore,
points $\bar{\bm{d}}\in\mathcal{P}^{m}_{\nu}$ and
$\bar{\bm{w}}\in\mathcal{P}^{n}$ are optimal solutions for problems in (1)
respectively if and only if $-A\bar{\bm{w}}\in\partial f(\bar{\bm{d}})$ and
$\bar{\bm{w}}\in\partial g(A^{\top}\bar{\bm{d}})$.

##  3 Related Work

Previous work for this setting [10] tries to solve the problem with a
boosting-based approach designed based on SoftBoost [24]. The resulting
algorithm named MbB (the shorthand for “Metarounding-by-Boosting”) is
guaranteed to terminate in
$O\left(\frac{L^{2}D_{\infty}}{\epsilon^{2}}\ln(Ln)\right)$ iterations, where
$L$ is the maximal $L_{1}$-norm of the $\epsilon$-approximate solution, which
is not known in advance. The parameter $L$ is used for the MbB algorithm, so
they search for a proper $L$ by doubling; they need to run MbB for multiple
$L$ until a proper $L$ is found. This drawback is crucial, but they showed the
efficiency of MbB by numerical experiments against the ellipsoid method [3].
There are two main advantages for MbB. The first is that MbB does not need to
know the approximation ratio $\alpha$ a priori. The other is that the
iteration bound for MbB does not depend on the iteration $T$ of online
algorithms. Previous works take $O(\ln T)$ [11] to $O(T)$ [3] iterations for
every round of online combinatorial optimization, where $T$ is the number of
total rounds for the online algorithm.

###  3.1 Margin optimization boosting algorithms

Boosting is a common machine learning technique that combines multiple low-
accuracy predictors to yield a highly accurate one. Since the predictors to be
combined are chosen from a set of exponentially many predictors, one cannot
compute the combined predictor directory. The main idea of boosting is to
collect the most promising predictors one by one until some condition is
satisfied.

Gradient Boosting Machines [8], such as XGBoost [5] and LightGBM [17], are
popular algorithms that minimize some loss functions. For binary
classification, it is well known that a predictor that maximizes the quantity
so-called “margin” instead of minimizing loss functions is known to have
smaller losses for unseen data [21]. Thus, some algorithms, such as LPBoost
[6], SoftBoost [24], and ERLPBoost [25], are proposed to optimize the margin.

###  3.2 Frank-Wolfe algorithms

Frank-Wolfe algorithms (FWs for shorthand), a.k.a. conditional gradient
algorithms, are the class of first-order optimization algorithms for convex
objectives over convex sets. FW was originally invented by Frank & Wolfe [7]
and then studied extensively [15, 19, 23]. FW algorithms are known to converge
in $O(\eta/\epsilon)$ iterations if the objective function is $\eta$-smooth
w.r.t. some norm.

Algorithm 1 Metarounding by ERLPBoost

0: $\bm{x}\in\mathop{\rm relax}(\mathcal{C})$, accuracy $\epsilon>0$, and a
relax-based $\alpha$-approx. algorithm $\mathcal{B}$.

1: Initialize $\bm{\ell}_{0}=\frac{1}{n}\bm{1}$ and set
$H^{\star}(C^{(0)}\bm{\lambda}_{0})=\infty$, $\eta=2\ln(n)/\epsilon$.

2: for $k=0,1,2,\dots,K$ do

3: Call the approximation algorithm to obtain
$\bm{c}_{k+1}\in\mathcal{B}(\bm{\ell}_{k})$.

4: if $\epsilon_{k}:=H^{\star}(C^{(k)}\bm{\lambda}_{k})-\max_{j\leq
k}\bm{c}_{j+1}\cdot\bm{\ell}_{j}\leq\epsilon/2$ then

5: Set $K=k$. break.

6: end if

7: Solve Problem (7) over
$C^{(k+1)}=\begin{bmatrix}\bm{c}_{1}&\bm{c}_{2}&\dots&\bm{c}_{k+1}\end{bmatrix}$.
$\bm{\lambda}_{k+1}\leftarrow\arg\min_{\bm{\lambda}\in\mathcal{P}^{k+1}}H^{\star}(C^{(k+1)}\bm{\lambda})$

8: Set $\bm{\ell}_{k+1}=\nabla
H^{\star}(C^{(k+1)}\bm{\lambda}_{k+1})=\arg\max_{\bm{\ell}\in\mathcal{L}}\bm{\ell}^{\top}C^{(k+1)}\bm{\lambda}_{k+1}-\frac{1}{\eta}\Delta(\bm{\ell})$.

9: end for

10: Solve Problem (2) over $\bm{c}_{1},\bm{c}_{2},\dots,\bm{c}_{K+1}$ and let
the optimal solution be $\bm{\lambda}^{\star}$.

10: $\bm{\lambda}^{\star}\in\mathcal{P}^{\mathcal{C}}$

##  4 Main idea

First, we revisit the approach by Fujita et al. [10]. In MbB, they formulate
an $\alpha$-metarounding as an algorithm that solves the following
optimization problem.

$\displaystyle\min_{\beta,\bm{\lambda}}\;\beta\qquad\text{s.t.}\quad$
$\displaystyle\sum_{\bm{c}\in\mathcal{C}}\lambda_{\bm{c}}c_{i}\leq\beta
x_{i},\quad\forall i\in[n],$ (2)
$\displaystyle\bm{\lambda}\in\mathcal{P}^{\mathcal{C}}.$

Since Problem (2), has exponentially many variables, they consider the dual
problem.

$\displaystyle\max_{\gamma,\bm{\ell}}\;\gamma\qquad\text{s.t.}\quad$
$\displaystyle\bm{\ell}\cdot\bm{c}\geq\gamma,\quad\forall\bm{c}\in\mathcal{C},$
(3) $\displaystyle\bm{\ell}\cdot\bm{x}=1,\quad\bm{\ell}\geq\bm{0}.$

The resulting dual problem has $O(|C|)$ constraints, so MbB uses the column-
generation approach like soft margin boosting algorithms; That is, starting
from $\tilde{\mathcal{C}}=\emptyset$, they solve a problem under the
constraints in (3) over $\tilde{\mathcal{C}}$ and then append a new
combinatorial vector $\bm{c}\in\mathcal{C}$ to $\tilde{\mathcal{C}}$ one by
one until some conditions are met. One can see that Problem (3) can be seen as
so-called “Edge minimization” in boosting (See, e.g., [6, 24]). In fact,
setting $\bm{x}=\bm{1}$ in Problem (3) corresponds to margin optimization.
With this insight, MbB is proposed based on a boosting algorithm named
SoftBoost [24], but the resulting guarantee is not as good as the one in
SoftBoost.

We want to use the state-of-the-art soft margin optimization named ERLPBoost
[25] to solve Problem (3), but there is an issue. ERLPBoost assumes
$\bm{\ell}$ belongs to probability simplex, but Problem (3) does not.
Constraint $\bm{\ell}\cdot\bm{x}=1$ makes the direct application hard. With
this issue, we use a different approach, a Frank-Wolfe-based approach similar
to the one in Mitsuboshi et al. [18]. The resulting algorithm is a
generalization of the soft margin boosting in the following sense: Our
algorithm becomes ERLPBoost if we set $\bm{x}=\bm{1}$.

###  4.1 Bounding the dual variable.

Before getting into our main result, we prove the following proposition.

######  Proposition 1.

Let $M:=\max\\{\frac{1}{x_{i}}\mid x_{i}\neq 0,i\in[n]\\}$. There is an
optimal solution $(\gamma^{\star},\bm{\ell}^{\star})$ to Problem (3)
satisfying $\ell^{\star}_{i}\leq M$ for all $i\in[n]$.

###### Proof.

Let $(\beta^{\star},\bm{\lambda}^{\star})$,
$(\gamma^{\star},\bm{\ell}^{\star})$ be optimal solutions to Problem (2), (3),
respectively. Without loss of generality, we can restrict the combinatorial
set to
$\mathcal{C}^{\star}:=\\{\bm{c}\in\mathcal{C}\mid\lambda_{\bm{c}}>0\\}$. We
consider the following two cases.

  1. 1.

$\forall i\in[n],x_{i}\neq 0$. This case,

$\bm{\ell}^{\star}\cdot\bm{x}>\sum_{j\neq
i}\ell^{\star}_{j}x_{j}+Mx_{i}>\sum_{j\neq
i}\ell^{\star}_{j}x_{j}+M\min_{k\in[n]}x_{k}\geq 1,$

which contradicts to the constraint $\bm{\ell}^{\star}\cdot\bm{x}=1$, so that
$\ell^{\star}_{i}\leq M$ for all $i\in[n]$.

  2. 2.

$\exists i\in[n],x_{i}=0$. Since $(\beta^{\star},\bm{\lambda}^{\star})$ is a
feasible solution,
$\sum_{\bm{c}\in\mathcal{C}^{\star}}\lambda^{\star}_{\bm{c}}c_{i}\leq\beta^{\star}x_{i}=0$,
which implies $\lambda^{\star}_{\bm{c}}=0$ or $\bm{c}_{i}=0$ for all
$\bm{c}\in\mathcal{C}^{\star}$. By definition of $\mathcal{C}^{\star}$,
$\lambda^{\star}_{\bm{c}}>0$ for all $\bm{c}\in\mathcal{C}^{\star}$, so we
have $c_{i}=0$ for all $\bm{c}\in\mathcal{C}^{\star}$. Thus,
$\ell^{\star}_{i}c_{i}=0$ does not affect to the objective value
$\gamma^{\star}$. So we can set $\ell^{\star}_{i}=M$ without losing the
optimality.

∎

Thus, restricting $\bm{\ell}\in[0,M]^{n}$ does not affect the optimal solution
to Problem (3), so we consider the following problem

$\displaystyle\max_{\gamma,\bm{\ell}}\;\gamma\qquad\text{s.t.}\quad$
$\displaystyle\bm{\ell}\cdot\bm{c}\geq\gamma,\quad\forall\bm{c}\in\mathcal{C},$
(4) $\displaystyle\bm{\ell}\in\mathcal{L},$

where $M$ is the value defined in Proposition 1 and
$\mathcal{L}:=\\{\bm{\ell}\in[0,M]^{n}\mid\bm{\ell}\cdot\bm{x}=1\\}$.

###  4.2 The regularization

Regularizing a smooth function to a convex objective function does not make
the objective smooth while regularizing a strongly convex function makes it
strongly convex. So, we add a strongly convex regularizer to the objective
function of (4) then consider the primal problem. We use
$\Delta(\bm{\ell})=\sum_{i=1}^{n}\ell_{i}x_{i}\ln\frac{\ell_{i}x_{i}}{1/n}$ as
the regularizer222We define $\ell_{i}x_{i}\ln\ell_{i}x_{i}=0$ if
$\ell_{i}x_{i}=0$. . We solve Problem (4) by ERLPBoost-based algorithm, shown
in Algorithm 1. Our algorithm aims to find an $\epsilon/2$-approximate
solution to the following problem.

$\displaystyle\max_{\gamma,\bm{\ell}}\;\gamma-\frac{1}{\eta}\Delta(\bm{\ell})\qquad\text{s.t.}\quad$
$\displaystyle\bm{\ell}\cdot\bm{c}\geq\gamma,\quad\forall\bm{c}\in\mathcal{C},$
(5) $\displaystyle\bm{\ell}\in\mathcal{L}.$

Note that, by definition, $\Delta(\bm{\ell})\leq\ln n$ holds for all
$\bm{\ell}\in\mathcal{L}$.

Let $I_{\mathcal{L}}$ be the barrier function for $\mathcal{L}$,
$J(\bm{\theta})=\max_{\bm{c}\in\mathcal{C}}\theta_{\bm{c}}$ be the max
function over $\mathbb{R}^{\mathcal{C}}$, and let
$H=I_{\mathcal{L}}+\frac{1}{\eta}\Delta$. With these notations, one can
rewrite Problem (5) as

$\displaystyle\max_{\bm{\ell}}\;-J(-\bm{\ell}^{\top}C)-I_{\mathcal{L}}(\bm{\ell})-\frac{1}{\eta}\Delta(\bm{\ell})=:\max_{\bm{\ell}}\;-J(-\bm{\ell}^{\top}C)-H(\bm{\ell}).$
(6)

where
$C=\begin{bmatrix}\bm{c}_{1}&\bm{c}_{2}&\dots\end{bmatrix}\in\mathbb{R}^{n\times\mathcal{C}}$
is the matrix whose column vectors are elements of $\mathcal{C}$.

By Theorem 1, Problem (6) is equivalent to the following dual problem333One
can verify the equality by setting $f=J^{\star}$, $g=H^{\star}$, and
$A=C^{\top}$. .

$\displaystyle\min_{\bm{\lambda}}\;J^{\star}(\bm{\lambda})+H^{\star}(C\bm{\lambda})=\min_{\bm{\lambda}\in\mathcal{P}^{\mathcal{C}}}\;H^{\star}(C\bm{\lambda})$
(7)

Since $H$ is $1/\eta$-strongly convex w.r.t. $L_{1}$-norm, its conjugate
$H^{\star}$ is $\eta$-smooth w.r.t. $L_{\infty}$-norm. The objective function
$H^{\star}$ is a smoothified version of the one in Problem (2).

The smoothness of the objective function is a key property for our analysis,
so we aim to solve Problem (7) instead of Problem (2). Algorithm 1 summarizes
our method, which maintains the quantity $\epsilon_{k}$ defined in Line 4.

######  Lemma 1.

Let $\bm{x}\in\mathop{\rm relax}(\mathcal{C})$ be an input to Algorithm 1 and
let
$\epsilon_{k}=H^{\star}(C^{(k)}\bm{\lambda})-\max_{j\in[k]}\bm{c}_{j}\cdot\bm{\ell}_{j-1}$
be the optimality gap defined in Line 4 of Algorithm 1. If $\eta\geq
2\ln(n)/\epsilon$, $\epsilon_{k}\leq\epsilon/2$ implies
$\mathbb{E}_{\bm{c}\sim\bm{\lambda}_{k}}[\bm{c}]\cdot\bm{\ell}\leq(\alpha+\epsilon)\bm{x}\cdot\bm{\ell}$
for all $\bm{\ell}\in\mathcal{L}$.

###### Proof.

Recall that

$\bm{\ell}_{k}=\nabla
H^{\star}(C^{(k)}\bm{\lambda}_{k})=\arg\max_{\bm{\ell}\in\mathcal{L}}\bm{\ell}^{\top}C^{(k)}\bm{\lambda}_{k}-\frac{1}{\eta}\Delta(\bm{\ell}).$

By definition of $H^{\star}$, the following holds for all
$\bm{\ell}\in\mathcal{L}$.

$\displaystyle H^{\star}(C^{(k)}\bm{\lambda}_{k})=\;$
$\displaystyle\bm{\ell}_{k}^{\top}C^{(k)}\bm{\lambda}_{k}-\frac{1}{\eta}\Delta(\bm{\ell}_{k})$
$\displaystyle\geq\;$
$\displaystyle\bm{\ell}^{\top}C^{(k)}\bm{\lambda}_{k}-\frac{1}{\eta}\Delta(\bm{\ell})\geq\bm{\ell}^{\top}C^{(k)}\bm{\lambda}_{k}-\frac{1}{\eta}\ln
n\geq\mathbb{E}_{\bm{c}\sim\bm{\lambda}_{k}}\left[\bm{\ell}\cdot\bm{c}\right]-\frac{\epsilon}{2}.$

On the other hand, by the assumption of the $\mathcal{B}$, we have

$\displaystyle\max_{j\in[k]}\bm{\ell}_{j}\cdot\bm{c}_{j+1}\leq$
$\displaystyle\max_{j\in[k]}\alpha\min_{\bm{p}\in\mathop{\rm
relax}(\mathcal{C})}\bm{\ell}_{j}\cdot\bm{p}\leq\alpha\max_{j\in[k]}\bm{\ell}_{j}\cdot\bm{x},$

where the second inequality comes from the fact that $\bm{x}\in\mathop{\rm
relax}(\mathcal{C})$. Since Algorithm 1 chooses $\bm{\ell}_{j}$ to satisfy
$\bm{\ell}_{j}\cdot\bm{x}=1$, the right-hand-side is equal to $\alpha\cdot 1$.
Combining the above relations, we have

$\epsilon_{k}=H^{\star}(C^{(k)}\bm{\lambda}_{k})-\max_{j\in[k]}\bm{c}_{j+1}\cdot\bm{\ell}_{j}\geq\mathbb{E}_{\bm{c}\sim\bm{\lambda}_{k}}\left[\bm{\ell}\cdot\bm{c}\right]-\frac{\epsilon}{2}-\alpha,\qquad\forall\bm{\ell}\in\mathcal{L}.$

Thus, $\epsilon_{k}\leq\epsilon/2$ implies

$\forall\bm{\ell}\in\mathcal{L},\qquad\mathbb{E}_{\bm{c}\sim\bm{\lambda}_{k}}\left[\bm{\ell}\cdot\bm{c}\right]\leq\alpha+\epsilon=(\alpha+\epsilon)\bm{x}\cdot\bm{\ell},$

which is the relation we desired. ∎

Lemma 1 assures that the output $\bm{\lambda}^{\star}$ of Algorithm 1 is an
$\epsilon$-approximate solution for Problem (2). Thus, we can focus on proving
how quickly the optimality gap converges. The following lemma shows the
convergence rate.

######  Lemma 2.

For all $\\{\mu_{k}\\}_{k\geq 1}\subset[0,1]$, the sequence
$\\{\epsilon_{k}\\}_{k\geq 1}$ generated by Algorithm 1 satisfies the
following relation.

$\displaystyle\forall k\geq
1,\quad\epsilon_{k+1}\leq\left(1-\mu_{k}\right)\epsilon_{k}+2\eta\mu_{k}^{2}D_{\infty}^{2}.$

###### Proof.

First of all, $\bm{\lambda}_{k+1}$ is chosen to optimize

$\displaystyle\min_{\bm{\lambda}\in\mathcal{P}^{k+1}}H^{\star}(C^{(k+1)}\bm{\lambda}),\qquad
C^{(k+1)}=\begin{bmatrix}\bm{c}_{1}&\bm{c}_{2}&\dots&\bm{c}_{k+1}\end{bmatrix}\in\\{0,1\\}^{n\times(k+1)},$

so that $H^{\star}(C^{(k+1)}\bm{\lambda}_{k+1})\leq
H^{\star}(C^{(k+1)}\bar{\bm{\lambda}}_{k+1})$ holds for all interior point
$\bar{\bm{\lambda}}_{k+1}=\bm{\lambda}_{k}+\mu_{k}\left(\bm{e}_{\bm{c}_{k+1}}-\bm{\lambda}_{k}\right)$
with $\mu_{k}\in[0,1]$. Thus, using the $\eta$-smoothness of $H^{\star}$,

$\displaystyle\epsilon_{k}-\epsilon_{k+1}\geq\;$ $\displaystyle
H^{\star}(C^{(k)}\bm{\lambda}_{k})-H^{\star}(C^{(k+1)}\bm{\lambda}_{k+1})$
$\displaystyle\geq\;$ $\displaystyle
H^{\star}(C^{(k)}\bm{\lambda}_{k})-H^{\star}(C^{(k+1)}\bar{\bm{\lambda}}_{k+1})$
$\displaystyle\geq\;$ $\displaystyle-\mu_{k}\nabla
H^{\star}(C^{(k+1)}\bm{\lambda}_{k})^{\top}C^{(k+1)}\left(\bm{e}_{\bm{c}_{k+1}}-\bm{\lambda}_{k}\right)$
$\displaystyle-\frac{\eta}{2}\mu_{k}^{2}\|C^{(k+1)}(\bm{e}_{\bm{c}_{k+1}}-\bm{\lambda}_{k})\|_{\infty}^{2}$
$\displaystyle\geq\;$
$\displaystyle-\mu_{k}\bm{\ell}_{k}C^{(k+1)}\left(\bm{e}_{\bm{c}_{k+1}}-\bm{\lambda}_{k}\right)-2\eta\mu_{k}^{2}D_{\infty}^{2}$
$\displaystyle=\;$
$\displaystyle-\mu_{k}\bm{\ell}_{k}\cdot\bm{c}_{k+1}+\mu_{k}\bm{\ell}_{k}^{\top}C^{(k)}\bm{\lambda}_{k}-2\eta\mu_{k}^{2}D_{\infty}^{2}$
$\displaystyle\geq\;$
$\displaystyle-\mu_{k}\bm{\ell}_{k}\cdot\bm{c}_{k+1}+\mu_{k}\left[\bm{\ell}_{k}^{\top}C^{(k)}\bm{\lambda}_{k}-\frac{1}{\eta}\sum_{i=1}^{n}(\ell_{k,i}x_{i})\ln\frac{\ell_{k,i}x_{i}}{1/n}\right]$
$\displaystyle-2\eta\mu_{k}^{2}D_{\infty}^{2}$ $\displaystyle=\;$
$\displaystyle-\mu_{k}\bm{\ell}_{k}\cdot\bm{c}_{k+1}+\mu_{k}H^{\star}(C^{(k)}\bm{\lambda}_{k})-2\eta\mu_{k}^{2}D_{\infty}^{2}$
$\displaystyle\geq\;$
$\displaystyle-\mu_{k}\max_{j\in[k]}\bm{\ell}_{j}\cdot\bm{c}_{j+1}+\mu_{k}H^{\star}(C^{(k)}\bm{\lambda}_{k})-2\eta\mu_{k}^{2}D_{\infty}^{2}$
$\displaystyle=\;$
$\displaystyle\mu_{k}\epsilon_{k}-2\eta\mu_{k}^{2}D_{\infty}^{2}.$

Rearraging terms, we get the recurrence relation. ∎

With Lemma 2, we can derive the convergence rate for Algorithm 1.

######  Theorem 2.

Algorithm 1 terminates in $K=O\left(\frac{D_{\infty}^{2}}{\epsilon^{2}}\ln
n\right)$ iterations.

###### Proof.

By Lemma 2, we have the following recurrence relation.

$\displaystyle\epsilon_{k+1}\leq\left(1-\mu_{k}\right)\epsilon_{k}+2\eta\mu_{k}^{2}D_{\infty}^{2}$
(8)

Since $\\{\lambda_{k}\\}_{k\geq 1}\subset[0,1]$ is arbitrary, we choose
$\lambda_{k}=\frac{2}{k+2}$. Now we will show the following by induction on
$k$.

$\displaystyle\epsilon_{k}\leq\frac{8\eta D_{\infty}^{2}}{k+2}$ (9)

For the base case $k=1$,
$\epsilon_{1}\leq\left(1-\frac{2}{0+2}\right)\epsilon_{0}+2\eta\left(\frac{2}{0+2}\right)^{2}D_{\infty}^{2}\leq\frac{8\eta
D_{\infty}^{2}}{1+2}$ holds. Thus Eq. (9) holds for $k=1$. Now, assume that
Eq. (9) holds for $k\geq 1$. By the recurrence formula (8) and the induction
hypothesis, we have

$\displaystyle\epsilon_{k+1}\leq$
$\displaystyle\left(1-\frac{2}{k+2}\right)\frac{8\eta
D_{\infty}^{2}}{k+2}+2\eta\left(\frac{2}{k+2}\right)^{2}D_{\infty}^{2}$
$\displaystyle=$ $\displaystyle\frac{k}{k+2}\frac{8\eta
D_{\infty}^{2}}{k+2}+2\eta\left(\frac{2}{k+2}\right)^{2}D_{\infty}^{2}$
$\displaystyle=$ $\displaystyle\frac{8\eta
D_{\infty}^{2}}{k+2}\left[\frac{k}{k+2}+\frac{1}{k+2}\right]\leq\frac{8\eta
D_{\infty}^{2}}{k+2}\frac{k+2}{k+3}=\frac{8\eta D_{\infty}^{2}}{k+3}$

Thus, Ineq. (9) holds for all $k$.

Setting the right-hand-side of (9) to be at most $\epsilon/2$, we
$\epsilon_{K}\leq\epsilon/2$ is achieved after $K\geq\frac{16\eta
D_{\infty}^{2}}{\epsilon}\geq\frac{32D_{\infty}^{2}}{\epsilon^{2}}\ln
n=O\left(\frac{D_{\infty}^{2}}{\epsilon^{2}}\ln n\right)$ iterations. ∎

By Theorem 2, $\epsilon_{K}\leq\frac{\epsilon}{2}$ holds after
$K=O\left(D_{\infty}^{2}\ln(n)/\epsilon^{2}\right)$ iterations. Further, Lemma
1 guarantees that the output $\bm{\lambda}_{K}\in\mathcal{P}^{\mathcal{C}}$
satisfies
$\mathbb{E}_{\bm{c}\sim\bm{\lambda}}[\bm{c}]\cdot\bm{\ell}\leq(\alpha+\epsilon)\bm{x}\cdot\bm{\ell}$
for all $\bm{\ell}\in\mathcal{L}$. Therefore, Algorithm 1 is a matarounding
algorithm.

##  5 Experiment

Figure 2: Comparison of previous works and our work. Both figure uses
combinatorial vector of dimension $n=20$. Left: comparison of running time.
Right: comparison of cumulative losses. KKL ’09 is the algorithm proposed by
Kakade et al., 2009 [16]. In this experiment, we set $|\mathcal{C}|=100$.

We demonstrate numerical experiment to show the effectiveness of Algorithm 1
compared to the previous work 444The program is available at
https://github.com/rmitsuboshi/metarounding_mitsuboshi.. We performed our
experiment on Intel Xeon Gold 6124 CPU 2.60GHz processors.

Our experiment was demonstrated on a set cover instance, generated with the
same procedure in [10]. We choose a fixed item size $m=10$ and measure the
running time by varying the number of sets $n\in\\{10,50,100,200,500,1000\\}$.
The input for metarounding is a solution of relaxed set cover problem with a
random cost vector $\bm{x}\in[0,1]^{n}$. Here, the relaxation means that the
variable takes value in $[0,1]$ instead of $\\{0,1\\}$.

_Running time._ We compared the algorithm proposed by Fujita et al. [10] and
our algorithm. We use the Gurobi optimizer 9.0.1555https://www.gurobi.com/ to
solve the sub-problems. To solve the entropy minimization, we used the
sequential quadratic minimization approach; it first minimizes the quadratic
approximation of the objective, then the objective function is updated to the
approximation around the optimal solution. This procedure repeats until
convergent. As shown in the left side of Figure 2, our algorithm is extremely
fast compared to the others.

_Cumulative loss._ We compared the algorithm proposed by Kakade et al. [16]
and our algorithm. In each round, the loss vector is generated randomly. Our
algorithm uses the classic online gradient descent algorithm [13] over
$[0,1]^{n}$ as the online linear optimization. The algorithm achieves
$R_{T}(1)=O(\sqrt{T})$ over $\mathop{\rm relax}(\mathcal{C})$, so combining it
with our algorithm achieves $R_{T}(\alpha)=O(\sqrt{T})$. As the right side of
Figure 2, our algorithm is empirically superior to the other.

##  6 Conclusion

In this paper, we proposed a metarounding algorithm that finds an
$\epsilon$-approximate solution in $O(D_{\infty}^{2}\ln(n)/\epsilon^{2})$
rounds. Even in the numerical experiments, the proposed algorithm is faster
than previous works.

## Acknowledgements

The work was supported by JSPS KAKENHI Grant Numbers JP19H014174 and
JP19H04067, respectively.

## References

  * [1] Borwein, J.M., Lewis, A.S.: Convex Analysis, pp. 65–96. Springer New York (2006)
  * [2] Carr, D.R., Vempala, S.: Randomized metarounding. Random Structure and Algorithms 20(1), 343–352 (2002)
  * [3] Carr, R.D., Vempala, S.S.: Randomized metarounding. Random Struct. Algorithms 20(3), 343–352 (2002)
  * [4] Cesa-Bianchi, N., Lugosi, G.: Combinatorial Bandits. Journal of Computer and System Sciences 78(5), 1404–1422 (2012)
  * [5] Chen, T., Guestrin, C.: Xgboost: A scalable tree boosting system. In: Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. p. 785–794. KDD ’16, Association for Computing Machinery (2016)
  * [6] Demiriz, A., Bennett, K.P., Shawe-Taylor, J.: Linear Programming Boosting via Column Generation. Machine Learning 46(1-3), 225–254 (2002)
  * [7] Frank, M., Wolfe, P.: An algorithm for quadratic programming. Naval Research Logistics Quarterly 3(1-2), 95–110 (1956)
  * [8] Friedman, J.H.: Greedy function approximation: A gradient boosting machine. The Annals of Statistics 29(5) (2001)
  * [9] Fujita, T., Hatano, K., Takimoto, E.: Combinatorial Online Prediction via Metarounding. In: Proceedings of 24th Annual Conference on Algorithmic Learning Theory (ALT 2013). LNCS, vol. 8139, pp. 68–82 (2013)
  * [10] Fujita, T., Hatano, K., Takimoto, E.: Combinatorial online prediction via metarounding. In: Algorithmic Learning Theory - 24th International Conference, ALT 2013, Singapore, October 6-9, 2013. Proceedings. Lecture Notes in Computer Science, vol. 8139, pp. 68–82. Springer (2013)
  * [11] Garber, D.: Efficient online linear optimization with approximation algorithms. Math. Oper. Res. 46(1), 204–220 (2021)
  * [12] György, A., Linder, T., Lugosi, G., Ottucsák, G.: The On-Line Shortest Path Problem Under Partial Monitoring. Journal of Machine Learning Research 8, 2369–2403 (2007)
  * [13] Hazan, E.: Introduction to online convex optimization. CoRR abs/1909.05207 (2019), http://arxiv.org/abs/1909.05207
  * [14] Helmbold, D.P., Warmuth, M.K.: Learning Permutations with Exponential Weights. Journal of Machine Learning Research 10, 1705–1736 (2009)
  * [15] Jaggi, M.: Revisiting frank-wolfe: Projection-free sparse convex optimization. In: Proceedings of the 30th International Conference on Machine Learning, ICML 2013. JMLR Workshop and Conference Proceedings, vol. 28, pp. 427–435. JMLR.org (2013)
  * [16] Kakade, S., Kalai, A.T., Ligett, L.: Playing games with approximation algorithms. SIAM Journal on Computing 39(3), 1018–1106 (2009)
  * [17] Ke, G., Meng, Q., Finley, T., Wang, T., Chen, W., Ma, W., Ye, Q., Liu, T.: Lightgbm: A highly efficient gradient boosting decision tree. In: Guyon, I., von Luxburg, U., Bengio, S., Wallach, H.M., Fergus, R., Vishwanathan, S.V.N., Garnett, R. (eds.) Advances in Neural Information Processing Systems 30: Annual Conference on Neural Information Processing Systems 2017, December 4-9, 2017, Long Beach, CA, USA. pp. 3146–3154 (2017)
  * [18] Mitsuboshi, R., Hatano, K., Takimoto, E.: Boosting as frank-wolfe (2022)
  * [19] Pedregosa, F., Négiar, G., Askari, A., Jaggi, M.: Linearly convergent frank-wolfe without line-search. In: The 23rd International Conference on Artificial Intelligence and Statistics, AISTATS 2020, 26-28 August 2020, Online [Palermo, Sicily, Italy]. Proceedings of Machine Learning Research, PMLR (2020)
  * [20] Propp, J.G., Wilson, D.B.: How to Get a Perfectly Random Sample from a Generic Markov Chain and Generate a Random Spanning Tree of a Directed Graph. Journal of Algorithms 27(2), 170–217 (may 1998)
  * [21] Schapire, R.E., Freund, Y., Bartlett, P., Lee, W.S.: Boosting the margin: a new explanation for the effectiveness of voting methods. The Annals of Statistics 26(5), 1651–1686 (1998)
  * [22] Shalev-Shwartz, S., Singer, Y.: On the equivalence of weak learnability and linear separability: new relaxations and efficient boosting algorithms. Mach. Learn. 80(2-3) (2010)
  * [23] Tsuji, K., Tanaka, K., Pokutta, S.: Pairwise conditional gradients without swap steps and sparser kernel herding. In: Chaudhuri, K., Jegelka, S., Song, L., Szepesvári, C., Niu, G., Sabato, S. (eds.) International Conference on Machine Learning, ICML 2022, 17-23 July 2022, Baltimore, Maryland, USA. Proceedings of Machine Learning Research, vol. 162 (2022)
  * [24] Warmuth, M., Glocer, K., Rätsch, G.: Boosting Algorithms for Maximizing the Soft Margin. In: Advances in Neural Information Processing Systems 20 (NIPS 2007). pp. 1585–1592 (2007)
  * [25] Warmuth, M.K., Glocer, K.A., Vishwanathan, S.V.N.: Entropy regularized lpboost. In: Algorithmic Learning Theory, 19th International Conference, ALT 2008, Budapest, Hungary, October 13-16, 2008. Proceedings. Lecture Notes in Computer Science, vol. 5254, pp. 256–271. Springer (2008)
  * [26] Warmuth, M.K., Kuzmin, D.: Randomized Online PCA Algorithms with Regret Bounds that are Logarithmic in the Dimension. Journal of Machine Learning Research 9, 2287–2320 (2008)
  * [27] Yasutake, S., Hatano, K., Kijima, S., Takimoto, E., Takeda, M.: Online Linear Optimization over Permutations. In: Proceedings of the 22nd International Symposium on Algorithms and Computation (ISAAC 2011). LNCS, vol. 7074, pp. 534–543 (2011)
  * [28] Yasutake, S., Hatano, K., Takimoto, E., Takeda, M.: Online Rank Aggregation. In: NIPS 2011 Workshop on Computational Trade-offs in Statistical Learning (COST) (2011)


# Physics-informed neural networks for solving thermo-mechanics problems of
functionally graded material

Mayank Raj
Department of Mechanical Engineering
Indian Institute of Technology Madras
Chennai-600036
<EMAIL_ADDRESS>
& Pramod Kumbhar
Department of Mechanical Engineering
Indian Institute of Technology Madras
Chennai-600036
<EMAIL_ADDRESS>
& Ratna Kumar Annabattula
Department of Mechanical Engineering
Indian Institute of Technology Madras
Chennai-600036
<EMAIL_ADDRESS>
Corresponding author

###### Abstract

Differential equations are indispensable to engineering and hence to
innovation. In recent years, physics-informed neural networks (PINN) have
emerged as a novel method for solving differential equations. PINN method has
the advantage of being meshless, scalable, and can potentially be intelligent
in terms of transferring the knowledge learned from solving one differential
equation to the other. The exploration in this field has majorly been limited
to solving linear-elasticity problems, crack propagation problems. This study
uses PINNs to solve coupled thermo-mechanics problems of materials with
functionally graded properties. An in-depth analysis of the PINN framework has
been carried out by understanding the training datasets, model architecture,
and loss functions. The efficacy of the PINN models in solving thermo-
mechanics differential equations has been measured by comparing the obtained
solutions either with analytical solutions or finite element method-based
solutions. While R2 score of more than 99% has been achieved in predicting
primary variables such as displacement and temperature fields, achieving the
same for secondary variables such as stress turns out to be more challenging.
This study is the first to implement the PINN framework for solving coupled
thermo-mechanics problems on composite materials. This study is expected to
enhance the understanding of the novel PINN framework and will be seminal for
further research on PINNs.

_Keywords_ Physics-informed neural networks $\cdot$ Deep learning $\cdot$
Differential equations $\cdot$ Functionally graded material

##  1 Introduction

The term “neural network" encompasses a wide range of machine learning models
with one similarity, i.e., so-called “neurons" as the basic building block [1,
2]. Neural networks have seen many applications in fields not limited to
computer vision [3], banking and finance [4], bio-informatics [5, 6],
renewable energy systems [7, 8] and material science [9, 10, 11]. In other
words, if there is a function approximation task, neural networks are
indisputable models to use. The use cases are not limited to just function
approximation tasks. Rather, neural networks dominate in image captioning and
language translation [12], which are one-to-many or many-to-many mapping
tasks. So what is behind the success of neural networks? Certainly, the
flexibility they offer in designing composite parametric functions and the
ease of optimizing those parameters [13] lends to their success. Often called
universal function approximators [14], these neural networks can, in theory,
fit most of the relational mappings. However, the challenge lies in designing
a neural network architecture so that at least one set of optimum weights fits
the desired function. The so-called training of a neural network, i.e.,
estimation of optimal parameters using an optimization algorithm, can also be
challenging. However, after going through those steps, one gets a trained
neural network model that can predict/estimate the desired target variable on
a given set of inputs in seconds. Alongside all the success stories, neural
networks also have their own set of drawbacks. Training often requires
considerable data, careful data sampling to cover variance in the domain and
target space, difficult parallelization of the training process, and dedicated
hardware such as GPUs [15] for fast training.

Differential equations are indispensable to computational engineering.
Irrespective of engineering subdomains, differential equations appear in one
form or another. Often one needs to solve such equations on the domain of
interest numerically. The state of the art of technique for numerically
solving such differential equations is the finite element method. While the
finite element method has its pros and cons, the physics-informed neural
network (PINN) is a promising alternative to solving differential equations.

The family of a neural network is vast, and variety can be in various forms,
such as their architecture (arrangement of neurons and layers), learning
process (supervised or unsupervised), the purpose of the model (regression or
classification). Physics-informed neural network (PINN) is yet another variant
in this family. Presently in its infancy, PINN can potentially disrupt how
differential equations are solved. PINN has an apparent advantage of being a
meshless method over FEM. Unlike conventional supervised learning Deep Neural
Networks (DNNs), PINN models are trained without the target samples. While
PINN in its current form can take more time than commercial finite element
solvers such as Abaqus® to solve simple elasticity problems, it can do better
when simulations get complex and the number of simulations increases, such as
in an optimization scenario. PINN can use the neural network transfer learning
concept [16] to train on new samples in marginal time.

In an abstract sense, a physics-informed neural network tries to optimize the
parameters of a parametric function that approximates the solution of the
differential equation on a given domain. The network does the optimization by
minimizing an objective called the loss function in the machine learning
community. The prefix “Physics-informed" comes from the fact that the loss
function is obtained using the underlying physics of the problem represented
by the differential equations. Often, the loss is the functional corresponding
to the differential equation.

The physics-informed neural network was first introduced by Raissi et al. [17]
as an alternate method to solve partial differential equations. PINNs can be
classified into two categories based on the implementation of the loss
function. One of the approaches uses the collocation technique, where simply
the residual of the differential equation is minimized by training the neural
network [18, 19]. While on the other hand, there are energy-based loss
functions where the functional of the differential equation is minimized [20].
While the collocation approach is straightforward, it involves computing
higher-order derivatives, resulting in difficult training. On the other hand,
the energy-based method [20] requires lower order derivatives due to the
functional’s weak form of the differential equation. Hence, energy-based
methods have become more popular with numerous applications in the last few
years. It is used in the prediction of crack propagation using the phase-field
method [21, 22]. Further, the ease with which adaptive meshing can be done in
PINN was shown by Goswami et al. [23]. Li et al. [24] compared the
effectiveness of the PINN approach with other ML-based techniques in terms of
various metrics such as coefficient of determination. Recently, Dwivedi et al.
[25, 26] introduced Physics Informed Extreme Learning Machine, which improves
the efficiency of Physics Informed Neural Network by utilizing ideas of
Extreme Learning Machine. Dwivedi et al. also introduced Distributed Physics
Informed Neural Network [27, 28, 29] where the domain can be split into
multiple sub-domains, and separate PINNs can be trained over them with
ensuring continuity of variables on the interface of sub-domains.

Functionally graded materials (FGMs) are a unique class of composites
incorporating the material properties of two or more different homogeneous
materials [30]. The material property variation is usually continuous and
smooth such that it can be represented as explicit smooth mathematical
functions. Usually, these materials trump the conventional homogeneous
components in terms of superior thermal and mechanical capabilities. Hence, it
is not surprising that the FGMs find ubiquitous applications ranging from
mechanical to bio-medical applications [31].

While there are many papers published in PINNs in the last year, none of them
presents a comprehensive discussion on the pros and cons of the methodology.
Moreover, the application has been limited to isotropic material, and
discussions have mostly been centered around the primary variables such as
displacement fields and the damage parameter in phase-field modeling of
fracture. There has not been any discussion of secondary variables such as
strain and stress fields that depend on primary variables’ derivatives.
Moreover, there has been a knowledge gap in understanding the convergence of
energy-based PINN neural networks.

This paper starts with a discussion on using PINN models to solve one-
dimensional differential equations. We specifically solve coupled thermo-
mechanics problems in solid mechanics for the first time using PINN. We
consider materials with spatially varying properties and show the robustness
of PINN to arrive at the desired solution. We further discuss how and how not
to implement different kinds of boundary conditions. Further on, we point out
that the PINN model converges to a physically meaningful loss function value.
We then demonstrate the efficacy of PINN in solving the two-dimensional
thermo-mechanics problem with functionally graded properties. We achieve
remarkable PINN prediction accuracy of displacement and temperature fields.
However, we show that excellent accuracy in primary variables such as
displacement fields does not imply the same for secondary variables such as
stress fields. We conclude the paper by discussing further avenues of research
on PINN models.

##  2 Methods

###  2.1 Governing differential equations for thermo-elasticity

Figure 1: Schematic of a solid object with spatially varying properties lying
in a two-dimensional space.

Let us consider a solid body occupying spatial domain
$\Omega\subset\mathcal{R}^{d}$. The solid body is composed of material with
thermo-elastic properties (see Figure 1). The $d$ dimensional domain is
bounded by a surface $\Gamma$. The bounding surface is defined locally in
terms of the outward unit normal vector $\mathbf{n}$. Let the elastic modulus
and thermal conductivity of body vary as $E(\mathbf{x})$ and $k(\mathbf{x})$
where $\mathbf{x}$ is a point in the domain $\Omega$. The static thermo-
elastic response of the body, displacement field $\mathbf{u}:\Omega\to R^{d}$
and temperature field $T:\Omega\to R$ is governed by the following
differential equations:

$\nabla\cdot\mathbf{\sigma}(\mathbf{x})+\mathbf{f}(\mathbf{x})=0\quad\forall\mathbf{x}\in\Omega$
(1)

$\nabla\cdot\left(k(\mathbf{x})\nabla
T(\mathbf{x})\right)=0\quad\forall\mathbf{x}\in\Omega$ (2)

under the boundary conditions:

$\begin{split}\mathbf{u}(\mathbf{x})&=\mathbf{u_{bc}(x)}\quad\forall\mathbf{x}\in\Gamma_{u}\\\
\mathbf{\sigma(x)}\cdot\mathbf{n}&=\mathbf{t_{bc}(x)}\quad\forall\mathbf{x}\in\Gamma_{t}\\\
T(\mathbf{x})&=T_{bc}(\mathbf{x})\quad\forall\mathbf{x}\in\Gamma_{T}\end{split}$
(3)

$\mathbf{\sigma(x)}$ and $\mathbf{f(x)}$ represents stress field and body
force respectively. Moreover, the prescribed boundary conditions
$\mathbf{u_{bc}}$, $\mathbf{t_{bc}}$ and $T_{bc}$ are defined such that
$\left(\Gamma_{u}\cup\Gamma_{t}\cup\Gamma_{T}\right)\subset\Gamma$ and
$\Gamma_{u}\cap\Gamma_{t}=\Phi$. The stress field is related to strain field
which in turn is related to displacement and temperature fields as follows:

$\begin{split}\mathbf{\sigma(x)}&=\mathbf{C(x)}:\mathbf{\varepsilon(x)}\\\
\mathbf{\varepsilon(x)}&=\mathbf{\varepsilon}^{el}(\mathbf{x})-\alpha\Delta
T(\mathbf{x})\mathbf{I}\\\
\varepsilon^{el}_{ij}&=\dfrac{u_{i,j}+u_{j,i}}{2}\end{split}$ (4)

where $\mathbf{\varepsilon}^{el}$ is the elastic strain while $\alpha\Delta T$
is the thermal strain. $\alpha$ is the coefficient of thermal expansion ,
$\mathbf{I}$ is the identity matrix and $\Delta
T(\mathbf{x})=T(\mathbf{x})-T_{0}$ where $T_{0}$ is the temperature of the
body in the reference state. The elastic material constitutive matrix
$\mathbf{C(x)}$ is written in terms of elastic modulus and Poisson’s ratio
$\nu$ as:

$\mathbf{C(x)}=E(\mathbf{x})\begin{bmatrix}\dfrac{1}{1-\nu^{2}}&\dfrac{\nu}{1-\nu^{2}}&0\\\
\dfrac{\nu}{1-\nu^{2}}&\dfrac{1}{1-\nu^{2}}&0\\\
0&0&\dfrac{1}{(1+\nu)}\end{bmatrix}$ (5)

for two-dimensional object with plane-stress assumption.

The functional corresponding to differential equations governing the thermo-
elastic response of solid body is given as:

$\begin{split}\mathcal{F}&=W_{\text{int}}-W_{\text{ext}}\\\
&=\int_{\Omega}\frac{1}{2}\mathbf{\varepsilon}^{el}(\mathbf{x}):\mathbf{C(x)}:\mathbf{\varepsilon}^{el}(\mathbf{x})d\Omega+\int_{\Omega}\frac{1}{2}k(\mathbf{x})\nabla
T(\mathbf{x})\cdot\nabla
T(\mathbf{x})d\Omega-\left(\int_{\Omega}\mathbf{f(x)}\cdot\mathbf{n}\mkern
3.0mud\Omega+\int_{d\Gamma_{t}}\mathbf{t(x)}\cdot\mathbf{n}\mkern
3.0mud\Gamma\right)\end{split}$ (6)

The functional $\mathcal{F}$ can be viewed as internal energy
$(W_{\text{int}})$ less work done by the external forces $(W_{\text{ext}})$.

In this study, there are a few examples where we have used PINN to predict
only the elastic response of the object under consideration. In that scenario,
only Equation 1 governs the response of object. Moreover, the coefficient of
thermal expansion $\left(\alpha\right)$ is set to be zero, and the internal
energy corresponding to temperature gradient is removed from Equation 2. After
setting up the differential equations and boundary conditions governing the
problem at hand, PINN is used to minimize the functional, thereby obtaining
the desired response fields $\mathbf{u(x)}$ and $T(\mathbf{x})$. The details
of the PINN framework are discussed in the following subsection.

###  2.2 Physics informed neural network

Neural networks are composite parametric functions designed to learn any
relational mapping. While the traditional neural networks used to have only
one hidden layer, the abundance of computational resources has made it
possible to design and train models with numerous hidden layers. Hence the
term deep neural network (DNN) has become synonymous with neural networks. The
process of training a DNN is called deep learning. A DNN is built up using the
following building blocks: (1) neurons, (2) input layer, (3) output layer, (4)
hidden layers, (5) activation function/s, and (6) weights and biases.
Mathematically, a fully-connected deep neural network with straight-forward
architecture can be described as:

$\displaystyle\begin{split}\mathbf{h}^{1}&=\psi(\mathbf{W}^{1}\mathbf{x}+\mathbf{b}^{1})\\\
\mathbf{h}^{2}&=\psi(\mathbf{W}^{1}\mathbf{h}^{1}+\mathbf{b}^{1})\\\ \vdots\\\
\mathbf{h}^{i}&=\psi(\mathbf{W}^{i}\mathbf{h}^{i-1}+\mathbf{b}^{i})\\\
\vdots\\\
\mathbf{\tilde{v}}&=\mathbf{W}^{o}\mathbf{h}^{o-1}+\mathbf{b}^{o}\end{split}$
(7)

where $\mathbf{x}$ and $\mathbf{\tilde{v}}$ are respectively the input and
output vectors. $\mathbf{W}^{i}$ are the weight matrices, $\mathbf{b}^{i}$ are
the bias parameters, $\mathbf{h}^{i}$ are hidden layer vectors and $\psi(.)$
is the activation function used in different layers. Moreover, a further
transformation, $\zeta:\mathcal{R}^{m}\to\mathcal{R}^{m}$ is applied on
$\tilde{\mathbf{v}}$ in order to satisfy Dirichlet boundary conditions:

$\mathbf{\tilde{z}}(\mathbf{x})=\zeta\left(\mathbf{x},\mathbf{\tilde{v}(x)}\right)$
(8)

More on the transformation $\zeta(.)$ is discussed in Section 3. Hence, the
sequence of operations in Equations 7 and 8 can be represented as a DNN
$\mathbf{\tilde{z}}(\mathbf{x}):\mathcal{R}^{m}\to\mathcal{R}^{n}$, which maps
an $m$-dimensional real vector $\mathbf{x}$ to $n$-dimensional real vector
$\mathbf{\tilde{z}}$. In the Results and Discussion Section 3
$\mathbf{\tilde{u}}$ and ${\tilde{T}}$ will be often used in place of
$\mathbf{\tilde{z}}$. The vector $\mathbf{\tilde{z}}$ will In order to
optimize parameters of a conventional DNN, a loss function is defined as
$\mathcal{L}(\mathbf{\tilde{z}},\mathbf{z}):\mathcal{R}^{m\times
m}\to\mathcal{R}$, where $\mathbf{z}$ is true value of the target variable and
$\mathcal{L}(.)$ can represent a loss function such as mean squared error.
However, in case of PINN, the loss function is defined only in terms of
$\mathbf{\tilde{u}}$ and its derivatives with respect to $\mathbf{x}$, for
example:
$\mathcal{L}(\mathbf{\tilde{z}},\mathbf{\tilde{z}_{x}}):\mathcal{R}^{m\times
m\times n}\to\mathcal{R}$.

In order to solve the ODEs given in Equations 1 and 2 under the boundary
conditions in Equation 3, the functional given is Equation 6 must be
discretized to obtain the loss function $\mathcal{L}$ for training the PINN.
Hence, the discretized functional can be written as:

$\begin{split}\mathcal{L}&=\sum_{\mathbf{x}_{i}\in\Omega}\frac{1}{2}w(\mathbf{x}_{i})\mathbf{\tilde{\varepsilon}}^{el}(\mathbf{x}_{i}):\mathbf{C}(\mathbf{x}_{i}):\mathbf{\tilde{\varepsilon}}^{el}(\mathbf{x}_{i})\\\
&+\sum_{\mathbf{x}_{i}\in\Omega}\frac{1}{2}w(\mathbf{x}_{i})k(\mathbf{x}_{i})\nabla\tilde{T}(\mathbf{x}_{i})\cdot\nabla\tilde{T}(\mathbf{x}_{i})\\\
&-\left(\sum_{\mathbf{x}_{i}\in\Omega}w(\mathbf{x}_{i})\mathbf{f(\mathbf{x}_{i})}\cdot\mathbf{n}+\sum_{\mathbf{x}_{i}\in
d\Gamma_{t}}w^{\prime}(\mathbf{x}_{i})\mathbf{t}(\mathbf{x}_{i})\cdot\mathbf{n}\right)\end{split}$
(9)

where $w(\mathbf{x})$ and $w^{\prime}(\mathbf{x})$ represent the weights for
numerical integration on domain and boundary nodes respectively. The PINN
predicted solutions $\mathbf{\hat{u}(x)}$ and $\hat{T}(\mathbf{x})$ is
obtained by minimizing $\mathcal{L}$ such that:

$\begin{bmatrix}\mathbf{\hat{u}(x)}\\\ {\hat{T}\mathbf{(x)}}\
\end{bmatrix}=\mathbf{\tilde{z}}=\arg\inf_{\tilde{\mathbf{z}}}\mathcal{L}(\tilde{\mathbf{z}},\tilde{\mathbf{z}}_{x})$
(10)

###  2.3 Training physics informed neural network

The parameters of DNN: weights and biases are optimized using the
backpropagation [13] algorithm. The algorithm computes the derivative of the
loss function with respect to all parameters and subsequently updates them
using the gradient descent rule. The weight update process is repeated for a
number of epochs until the convergence criteria are satisfied. The parameters
of the PINN model are updated using the Adam optimizer [32], which is an
enhanced version of gradient descent. Computing derivatives in the
backpropagation algorithm is at the core of the training process. The
technique called automatic differentiation [33] is used for computing
derivatives and hence is at the core of any deep learning library such as
Tensorflow [34]. To understand automatic differentiation, let’s take a step
back to reflect on how neural networks (or any such parametric function) are
built in Tensorflow. To build a neural network, one uses the so-called
building block parametric operators in Tensorflow. Examples of such operators
are nothing but matrix multiplication, addition, and an activation function.
It is to be noted that the gradients to all such operators are already defined
in Tensorflow. For example, let’s look at the matrix multiplication operation.
Let us say there is a matrix $\mathbf{W}$, which when multiplied with vector
$\mathbf{x}$ of compatible size, returns the matrix multiplication result from
$\mathbf{y}=\mathbf{W}\mathbf{x}$. Hence, whenever one needs to compute the
gradient of this operation i.e, ${dy_{i}}/{dx_{j}}$, the result is the matrix
$\mathbf{W}$. In a nutshell, the parametric PINN model is built of such
fundamental operations whose derivative is well defined in terms of parameters
themselves. Hence, one can compute the derivative of the whole PINN function
by applying chain rule and using the derivatives of the basic operations.

###  2.4 Measure of PINN accuracy - coefficient of determination

The accuracy of PINN prediction is measured in terms of coefficient of
determination, popularly known as R2 score. It is a measure of the amount of
variance in the true target variable captured by the predicted values.
Mathematically, it can be defined as:

$R^{2}\text{
score}=1-\dfrac{\sum_{i\in\Omega}(u_{i}-\hat{u_{i}})^{2}}{\sum_{i\in\Omega}(u_{i}-\bar{u})^{2}}$
(11)

where $\Omega$ represents the domain. $u_{i}$, $\bar{u}$ and $\hat{u_{i}}$ are
scalar values and represent true displacement values, mean of true
displacement values, and predicted displacement values. The coefficient of
determination can similarly be obtained for any scalar variable of interest
such as a particular strain component. While R2 score is in general captures
the performance of PINN, it may not be a suitable measure in some cases. More
disucssion on R2 score is provided in Section 3.

##  3 Results and Discussion

###  3.1 1D-FGM-ELAS-DIRCH

Consider a bar with uniform cross-section and unit length lying along the
$x$-axis from $x=0$ to $x=1$. The bar is assumed to be functionally graded,
i.e. the elastic modulus varies spatially as $E(x)=1/(1+x)$. Moreover, the bar
is assumed to be fixed at one end, and a load is applied on the other to
observe the elastic response. Hence, for a displacement controlled analysis,
the boundary conditions are given as $u(0)=0,u(1)=1$, where $u(x)$ represents
the displacement field of the bar as a function of $x$ under the given load.
The differential equation governing the elastic response of object for given
properties and boundary conditions is given in Table 1. The solution $u({x})$
to the differential equation can be analytically derived to be:

$u(x)=\dfrac{x^{2}+2x}{3}$ (12)

Figure 2: Schematic of one-dimensional bar with spatially varying stiffness
and displacement node. Table 1: Differential equations and loss functions used
in training PINNs for solving different problems.

Problem Code | Differential Equation | Loss function
---|---|---
1D-FGM-ELAS-DIRCH | $\dfrac{d}{dx}\left(\dfrac{1}{1+x}\dfrac{du}{dx}\right)=0$ | $\dfrac{1}{N}\sum_{i}^{N}\dfrac{1}{2(x_{i}+1)}\left(\left.\dfrac{d{u}}{dx}\right\rvert_{x_{i}}\right)^{2}$
1D-FGM-ELAS-NEU | $\dfrac{d}{dx}\left(\dfrac{1}{1+x}\dfrac{du}{dx}\right)=0$ | $\dfrac{1}{N}\sum_{i}^{N}\dfrac{1}{2(x_{i}+1)}\left(\left.\dfrac{d{u}}{dx}\right\rvert_{x_{i}}\right)^{2}-t(1){u}(1)$
1D-ELAS-BF | $\dfrac{d^{2}u}{d^{2}x}+1=0$ | $\dfrac{1}{N}\sum_{i}^{N}\left(\dfrac{1}{2}\left(\left.\dfrac{d{u}}{dx}\right\rvert_{x_{i}}\right)^{2}-u_{i}\right)-t(1)u(1)$
1D-FGM-ELAS-THERMO | $\dfrac{d}{dx}\left(\dfrac{1}{1+x}\left(\dfrac{du}{dx}-\alpha(T-T_{0})\right)\right)=0$ | $\dfrac{1}{N}\sum_{i=1}^{N}\dfrac{1}{2(1+x_{i})}\left(\left.\dfrac{du}{dx}\right|_{x_{i}}-T(x_{i})\right)^{2}+\dfrac{10}{(1+x_{i})}\left(\left.\dfrac{dT}{dx}\right|_{x_{i}}\right)^{2}$
| $\dfrac{d}{dx}\left(k(x)\dfrac{\partial T}{\partial x}\right)=0$
Kirsch’s Problem | $\nabla\cdot\mathbf{\sigma(x)}=0$ | $\sum_{\mathbf{x}_{i}\in\Omega}{\dfrac{1}{2}w(\mathbf{x}_{i})\varepsilon(\mathbf{x}_{i}):\mathrm{C}:\varepsilon(\mathbf{x}_{i})}-\sum_{\mathbf{x}_{i}\in\Omega_{u}}w^{\prime}(\mathbf{x}_{i})\mathbf{t}(\mathbf{x}_{i})\cdot\mathbf{n}(\mathbf{x}_{i})$
2D-FGM-ELAS-NEU | $\nabla\cdot\mathbf{\sigma(x)}=0$ | $\sum_{\mathbf{x}_{i}\in\Omega}{\dfrac{1}{2}w(\mathbf{x}_{i})\varepsilon(\mathbf{x}_{i}):\mathrm{C}(\mathbf{x}_{i}):\varepsilon(\mathbf{x}_{i})}-\sum_{\mathbf{x}_{i}\in\Omega_{\Gamma_{t}}}\mathbf{t(x_{i})}\cdot\mathbf{n(x_{i})}$
2D-FGM-ELAS-DIRCH | $\nabla\cdot\mathbf{\sigma(x)}=0$ | $\sum_{\mathbf{x}_{i}\in\Omega}{\dfrac{1}{2}w(\mathbf{x}_{i})\varepsilon(\mathbf{x}_{i}):\mathrm{C}(\mathbf{x}_{i}):\varepsilon(\mathbf{x}_{i})}$
2D-FGM-THERMO-ELAS | $\nabla\cdot\mathbf{\sigma(x)}=0$ | $\sum_{\mathbf{x}_{i}\in\Omega}\left({\dfrac{1}{2}w(\mathbf{x}_{i})\varepsilon^{\text{el}}(\mathbf{x}_{i}):\mathrm{C}(\mathbf{x}_{i}):\varepsilon^{\text{el}}(\mathbf{x}_{i})}+\dfrac{1}{2}w(\mathbf{x}_{i})\dfrac{10}{1+x}\nabla T(\mathbf{x}_{i})\cdot\nabla T(\mathbf{x}_{i})\right)$
| $\nabla\cdot\left(k(\mathbf{x})\nabla T\right)=0$

It is to be noted that throughout this study, consistent units of various
physical quantities can be assumed. Hence, we do not explicitly mention units
in various figures and tables. Now, to solve the above differential equation
using a PINN, one needs to understand the three fundamental elements of the
PINN framework: (1) training dataset, (2) choice of neural network
architecture, (3) loss function. The training dataset or the input to the PINN
will be the spatial coordinate of nodes in the domain. Hence, the domain has
to be discretized into a set of nodes; however, unlike FEM, there is no notion
of elements in the PINN framework. While there is the freedom to non-uniformly
discretize the domain into nodes, for simplicity, the nodes are chosen to be
uniformly distributed in this example. It is to be noted that there is no
target dataset in training PINN models.

Figure 3: Representative PINN architecture for purely displacement controlled
analysis of 1D material.

The next step is to decide the neural network architecture. While the PINN
framework offers a lot of flexibility in determining the model architecture, a
few of such parameters get determined by the task at hand. For example, the
number of nodes in the input and output layers of the model must respectively
be equal to the dimensionality of the problem. In the case of the one-
dimensional solid mechanics problem, both the input and output layer has one
neuron each. The input neuron corresponds to the spatial position of nodes,
and the output neuron corresponds to the displacement. Moreover, in the case
of purely Dirichlet boundary conditions, the PINN model architecture can be
modified to satisfy displacement boundary conditions implicitly (see Equation
8).

Let us assume that there is a blackbox neural network function
$\tilde{v}(x):\mathrm{\Omega}\to\mathcal{R}$, as defined in Equation 7, which
maps points in the spatial domain $\mathrm{\Omega}$ to a real number. Now, the
Dirichlet boundary conditions as assumed in this example can be implemented by
the transformation $\tilde{u}(x):=x+x(1-x)\tilde{v}(x)$ thereby ensuring
$\tilde{u}(0)=0$ and $\tilde{u}(1)=1$. In Figure 3, connections in black
colour correspond to $\tilde{u}(x)$, while the connection in blue is the
transformation applied over $\tilde{v}(x)$ to get $\tilde{u}(x)$. The details
of neural network $\tilde{v}(x)$ is given in Table 2. It is to be noted that
one cannot do any transformation to implement Dirichlet BC such as
$\tilde{u}(x):=x^{2}+x^{2}(1-x^{2})\tilde{v}(x)$. This transformation will not
only ensure $\tilde{u}(0)=0$ and $\tilde{u}(1)=1$, but it will also
simultaneously put constraints on the displacement gradient which is
$\tilde{u}^{\prime}(0)=0$ for the given case. Hence, a transformation which
puts undesired constraints on gradients of primary variables must be avoided.

Table 2: Details of the neural network architecture
$(\mathbf{\tilde{v}}(\mathbf{x}))$ used in various problems.

Problem Code |  | No. of
---
input nodes
| No. of
---
output nodes
| No. of
---
hidden layers
| No. of nodes
---
in hidden layers
| Activation
---
function
1D-FGM-ELAS-DIRCH | 1 | 1 | 2 | [5, 5] | tanh(.)
1D-FGM-ELAS-NEU | 1 | 1 | 2 | [5, 5] | tanh(.)
1D-ELAS-BF | 1 | 1 | 2 | [5, 5] | tanh(.)
1D-FGM-THERMO-ELAS | 1 | (1,1) | (2,1) | [(5, 5), (10, 10), (5,5)] | tanh(.)
Kirch’s Problem | 2 | 2 | 3 | [100, 100, 100] | tanh2(.)
2D-FGM-ELAS-NEU | 2 | 2 | 3 | [100, 100, 100] | tanh2(.)
2D-FGM-ELAS-DIRCH | 2 | 2 | 3 | [100, 100, 100] | tanh2(.)
2D-FGM-THERMO-ELAS | 2 | 3 | 2 | [300, 300] | elu2(.)

The last and the most crucial element of the PINN framework is the loss
function used in training such models. As highlighted earlier, the loss used
sets the PINN model apart from the rest of the neural network family. The PINN
loss function is discretized form of the functional corresponding to the
underlying differential equation. Hence, the loss function will also change
from one differential equation to the other. It is to be noted that
discretization is done to numerically evaluate the functional in Equation 6
for training the PINN. The neural network model trains itself (fits the
differential equation) by minimizing the loss function and thereby finding a
function that minimizes the functional corresponding to the differential
equation. However, the tricky part is that the functional and hence the loss
function often depends not only on the primary target variables but also on
their derivatives. Differential equation and corresponding loss function can
be found in table Table 1.

In the context of solid mechanics, the loss function is the internal strain
energy minus work done by external forces on the object under consideration.
It is evident from Table 1, the loss involves computing strain values which
depends on gradient of $u(x)$ with respect to $x$. Since the neural network is
a parametric composite function, it is possible to compute the derivative of
the target variable $\tilde{u}(x)$ with respect to the input variable $x$ in
terms of parameters of the neural network. While computing such derivatives
can be a tedious task given the complexity of the neural network, the state-
of-the-art deep learning libraries turn out to be the savior. The derivative
can be easily computed using the automatic differentiation technique
(discussed in Section 2.3) implemented using the library Tensorflow. The
complete pipeline of the PINN framework is summarized in Figure 4.

Figure 4: Flowchart of PINN framework

After training the neural network with physics-informed loss, the given one-
dimensional differential equation is solved to obtain $\hat{u}(x)$ which is
the same as $\tilde{u}(x)$ at the end of training PINN. The comparison between
$\hat{u}(x)$, and analytical solution is shown in the Figure 5. It is
worthwhile to reiterate that the physics-informed loss must converge to the
minimum strain energy. The training of neural networks is an optimization task
where the global minima should ideally correspond to the unique minimum energy
solution of the system. R2 for different variables in given in Table 3.

Figure 5: Decline in physics-informed loss (left) and comparison between
analytical and predicted solution obtained at the end of training process
(right) for one-dimensional elasticity problem with displacement controlled
load (Section 3.1). Table 3: Coefficient of determination (R2 score) for
various variables of interest in different solid mechanics problems.

Problem Code | $u_{1}$ | $u_{2}$ | $\sigma_{11}$ | $\sigma_{12}$ | $\sigma_{22}$ | $\varepsilon^{el}_{11}$ | $\varepsilon^{el}_{12}$ | $\varepsilon^{el}_{22}$ | $T$
---|---|---|---|---|---|---|---|---|---
1D-FGM-ELAS-DIRCH | 99.94% | - | $<0$ | - | - | 95.14% | - | - | -
1D-FGM-ELAS-NEU | 99.99% | - | $<0$ | - | - | 99.99% | - | - | -
1D-ELAS-BF | 99.98% | - | 99.65% | - | - | 99.65% | - | - | -
1D-ELAS-THERMO-ELAS | 99.99% | - | $<0$ | - | - | 97.61% | - | - | 99.94%
Kirsch’s Problem | 99.89% | 99.02% | 97.52% | 99.29% | 89.75% | 97.32% | 99.29% | 88.05% | -
2D-FGM-ELAS-NEU | 99.83% | 99.99% | 83.01% | 80% | $<0$ | 97.11% | 64.11% | 96.85% | -
2D-FGM-ELAS-DIRCH | 99.92% | 99.99% | 81.78% | 90.81% | $<0$ | 98.16% | 87.34% | 96.75% | -
2D-FGM-THERMO-ELAS | 99.99% | 99.99% | 88.84% | 70.91% | 73.92% | 99.93% | 99.95% | 99.95% | 99.97%

###  3.2 1D-FGM-ELAS-NEU

The domain is again assumed to be spanning from $x=0$ to $x=1$, while the
cross-section is uniform. With the same functionally graded elastic modulus,
the differential equation governing the elastic response remains the same as
the last problem (see Table 1). While the left end of the bar is fixed, a unit
traction $t$ is applied on the right end such that $u(0)=0,\quad t(1)=1$. The
analytical solution for this case is given as:

$u(x)=\dfrac{2x+x^{2}}{2}$ (13)

The training dataset remains the same, i.e., the set of nodes on the domain.
While the neural network architecture implicitly ensured the Dirichlet
boundary conditions, implementation of Neumann boundary condition is ensured
in the definition of loss function given in Table 1. Since the loss function
is internal energy minus the work done by external forces, unlike the other
use cases of a neural network, the loss function does have a physical
significance in the case of PINN. The loss can be expected to converge to zero
only in situations where the internal energy stored is equal to the work done.
Otherwise, it can very well either saturate to a negative value if work done
by the external force is more than the energy stored or can saturate to a
positive value if the loss is only composed of internal energy. It is also to
be noted that loss function is not always internal energy minus the work done
by external force even when the later can be a non-zero quantity. For example,
in previous example the loss function was only composed of internal energy
even though work done by external force would be finite. So whenever, the
Dirichlet boundary condition is satisfied by a transformation abstracted in
Equation 8, the work done by the corresponding reaction force must not be
subtracted in the loss function.

The comparison between analytical and PINN solutions for this case is shown in
Figure 6. It can be observed that the loss converges to a negative value as
the work done by an external force is more than that of strain energy stored
in the bar. Moreover, a good match is observed between the analytical and PINN
predicted solution. The R2 score for different variables is given in Table 3.

Figure 6: Decline in physics-informed loss (left) and comparison between
analytical and predicted solution obtained at the end of training process
(right) for one-dimensional elasticity problem with force controlled load
(Section 3.2).

###  3.3 1D-ELAS-BF

In the last two examples, the use of PINN to implement Dirichlet and Neumann
boundary conditions was discussed. We consider one more example with body
force to demonstrate the robustness of the PINN framework. In this case, we
will assume the elastic modulus to be unity over the domain. Other properties
such as area of cross-section, mass density are assumed to be unity. The
governing equations are given in Table 1. Boundary conditions are applied such
that $u(0)=0,\quad t(1)=1$. The analytical solution for this case is given by:

$u(x)=\dfrac{4x-x^{2}}{2}$ (14)

Unsurprisingly, the PINN model can match the analytical solution in this case
as well. The comparison is shown in Figure 7 and R2 score was given in Table
3.

Figure 7: Decline in physics-informed loss (left) and comparison between
analytical and predicted solution obtained at the training process (right) for
1D mechanics problem with force controlled analysis and body force (Section
3.3).

It is quite remarkable that the PINN model can perfectly fit solutions to
differential equations under different boundary conditions. However, one
cannot just expect the black-box PINN model to fit the solution of all
possible differential equations. For example, a complex PINN model that is a
non-linear function can never fit a differential equation where the solution
$u$ is a constant value or linear in $x$. A practical example would be
estimating the displacement field in a coarsely discretized multi-phase
composite microstructure. Due to the coarse discretization, the true
displacement would be piece-wise linear functions. Moreover, the gradient at
the interface of two adjacent phases can be very high depending on the
difference in their stiffness. Hence, it will be difficult for PINN to fit
that displacement field with a high displacement gradient at the interface of
phases. However, this inability of PINN need not necessarily be a disadvantage
because coarse discretization is only a problem pertaining to the modeling of
materials. In reality, displacement fields in materials will be reasonably
smooth unless there is a discontinuity in the material due to cracks or
manufacturing defects.

It is to be noted that all the discussed examples were selected such that the
solution is quadratic in nature. While it is often challenging to comment on
the nature of the black-box PINN function, it can undoubtedly be said from the
excellent match that the models were just aptly parameterized to fit solutions
of quadratic nature.

###  3.4 1D-FGM-THERMO-ELAS

In this section, we demonstrate the ability of PINN to estimate thermo-
mechanical response of the functionally graded material. Let’s again consider
a 1D bar $\mathbf{\Omega}$ from $x=0$ to $x=1$. The thermo-elastic response of
the object is governed by the differetial equations given in Table 1. Let the
elastic modulus of object vary as $E(x)={1}/\left(1+x\right)$ and conductivity
to vary as $k(x)={10}/\left(1+x\right)$. The mechanical and thermal boundary
conditions are imposed such that $u(0)=0;\quad u(1)=1;\quad T(0)=0;\quad
T(1)=1$ where $u(x)$ and $T(x)$ are respectively the displacement and
temperature fields. The coefficient of thermal expansion $\alpha$ and
reference temperature $T_{0}$ are assumed to be unity and zero respectively.
The analytical solution for the purpose of comparison can be derived to be:

$u(x)=\dfrac{x^{3}}{9}+\dfrac{14}{27}x^{2}+\dfrac{10}{27}x$ (15)

and

$T(x)=\dfrac{2x+x^{2}}{3}$ (16)

The usual steps described in last section are followed in order to solve the
pair of differential equations using PINN. The domain is discretized into a
set of uniformly spaced $N=50$ nodes. The neural network
$\tilde{\mathbf{v}}(x)=\left(\tilde{v_{1}}(x)\quad\tilde{v_{2}}(x)\right)^{T}$,
is transformed as follows to ensure satisfaction of Dirichlet boundary
conditions:

$\begin{bmatrix}\tilde{u}(x)\\\
\tilde{T}(x)\end{bmatrix}=\begin{bmatrix}x+(1-x)\tilde{v}_{1}(x)\\\
x+(1-x)\tilde{v}_{2}(x)\\\ \end{bmatrix}$ (17)

The neural network architecture of $\mathbf{v}(x)$ is given in Table 2.
Further, the loss for training the PINN is given in Table 1. The model is
sufficiently trained and comparison between analytical and PINN predicted
solutions can be observed in Figure 8. The R2 score for primary and secondary
variables are given in Table 3.

Figure 8: Reduction is loss with epoch (left), comparison of analytical
solutions with predicted displacement field (middle) and temperature field
(right) obtained at the end of training process (Section 3.4).

###  3.5 Kirsch’s Problem

Let us consider an infinite rectangular plate with a circular hole. Let us
assume that the circular hole has a radius of $0.1$ and is centered at
$(0,0)$. Moreover, the plate’s thickness and the material’s elastic modulus
are assumed to be unity with appropriate units. Uniaxial traction
$\mathbf{t}=[1\quad 0]^{T}$ is applied in x-direction on both ends of the
plate. The problem can be solved by considering a quarter of the infinite
plate and applying Dirichlet boundary conditions as per symmetry. Moreover,
for the purpose of simulation, the domain is assumed to be finite, as shown in
Figure 9.

Figure 9: Domain, boundary conditions and applied load for Kirsch’s problem
(Section 3.5).

Let us say that
$\mathbf{\tilde{v}}(\mathbf{x}):\mathrm{\Omega}\to\mathcal{R}^{2}$ is the
neural network model which maps spatial coordinates vector $(\mathbf{x})$ to a
real-valued vector in $\mathcal{R}^{2}$. $\Omega$ represents the set of points
that constitute the plate while $\Gamma$ is the set of points on the boundary
of the plate. The details of neural network
$\mathbf{\tilde{v}}(\mathbf{x}):\mathrm{\Omega}\to\mathcal{R}^{2}$ is given in
Table 2. To ensure the satisfaction of Dirichlet boundary conditions, the
following transformation is applied
$\mathbf{\tilde{u}}(\mathbf{x})=\mathbf{x}\odot\tilde{v}(\mathbf{x})$, where
$\mathbf{\tilde{u}({x})}$ gives displacement vector and $\odot$ represents
Hadamard product. The governing equations and corresponding loss function is
given in Table 1. Under the plane stress conditions, the material constitutive
matrix is given in Equation 5. Distribution of nodes and corresponding weights
($w(\mathbf{x}))$ are shown in the Figure 10. It is to be noted that weights
are required for numerical integration in case nodes are non-uniformly
distributed over the domain. In the case of uniform distribution of nodes, as
in the case of other examples, the weights would be a constant value.

Figure 10: Distribution of nodes and weights in Kirsch’s Problem (Section
3.5). Non-uniform distribution of nodes is necessary to model curvature of the
hole.

The PINN model is trained for 2500 epochs in this case. The loss, which is the
internal strain energy - work done decreased rapidly at the beginning of the
training process and then saturated to a value of approximately $-0.512$. The
reduction in the loss vs. epochs is shown in Figure 11. In order to better
understand the training process, it is worthwhile to look at the distribution
of activation functions, weights and biases, and gradients for weight updates
throughout the training process. They are given in Appendix: Figures 26, 27
and 28. The comparison between predicted and analytical solutions is shown in
Figure 12, while the same for stress fields can be found in Figure 13. It is
observed that the predicted fields are in good agreement with the analytical
solution. The R2 score for various variables can found in Table 3.

Figure 11: Reduction of loss vs epoch as the PINN trained in Kirsch’s Problem
(Section 3.5).

(a) $u_{1}(\mathbf{x})$

(b) $u_{2}(\mathbf{x})$

Figure 12: Comparison between PINN predicted displacement field and
corresponding analytical solution (Section 3.5).

(a) $\sigma_{11}(\mathbf{x})$

(b) $\sigma_{12}(\mathbf{x})$

(c) $\sigma_{22}(\mathbf{x})$

Figure 13: Comparison between PINN predicted stress fields and corresponding
analytical solutions (Section 3.5).

###  3.6 2D-FGM-ELAS-NEU

Let us consider a rectangular domain $\Omega$ as shown in the Figure 14. As
shown in the figure, both the degrees of freedom of nodes at the bottom edge
are constrained to move. The elastic modulus of the plate is varying along
$x_{2}$ direction and is given as $E(x_{1},x_{2})={1}/\left(1+x_{2}\right)$.

Figure 14: Domain and boundary conditions for functionally graded material
object under force load defined in Section 3.6.

A uniform traction $\mathbf{t}$ is applied on the top boundary. Let us say
that $\tilde{v}(\mathbf{x}):\mathrm{\Omega}\to\mathcal{R}^{2}$ is the neural
network model which maps spatial coordinates vector $\mathbf{x}$ in domain
$\Omega$ to a real-valued vector in $\mathcal{R}^{2}$. $\Omega$ represents the
set of points that constitute the plate, while $\Gamma$ is the set of points
on the boundary of the plate. In order to satisfy Dirichlet boundary
conditions, the following transformation is applied:

$\mathbf{\tilde{u}(x)}=\begin{bmatrix}\tilde{u}_{1}\\\
\tilde{u}_{2}\end{bmatrix}=\begin{bmatrix}x_{2}\\\
x_{2}\end{bmatrix}\odot\begin{bmatrix}\tilde{v}_{1}\\\
\tilde{v}_{2}\end{bmatrix}$ (18)

where $u(\mathbf{x})$ gives displacement vector and $\odot$ represents
Hadamard product. Decrease in loss (internal strain energy - work done by an
external force) with epoch count as the PINN is trained shown in Figure 15.
Comparisons between displacement and stress fields predicted by PINN model and
the Abaqus FEM simulation are shown in Figure 15 and Figure 16, and 17. A good
match between the Abaqus and the PINN predicted solution is observed. The R2
score, for this case is given in Table 3. It is to be noted that for
$\sigma_{22}$ the R2 scores are negative, which usually means quite a bad
prediction except when the variance in the target variable is very close to
zero. The true value of $\sigma_{22}$ is more or less constant over the domain
where R2 score is negative. This results in small variance in true target
variable and hence a small value in denominator in calculation of R2 score
(see Equation 11). Therefore, even a small absolute error would result in a
negative R2 score.

Figure 15: Reduction in loss as the PINN is trained to predict elastic
response of two-dimensional functionally graded material object defined in
Section 3.6.

(a) $u_{1}\left(\mathbf{x}\right)$

(b) $u_{2}\left(\mathbf{x}\right)$

Figure 16: Comparison of PINN predicted displacement fields and corresponding
Abaqus solution for problem defined in Section 3.6.

(a) $\sigma_{11}\left(\mathbf{x}\right)$

(b) $\sigma_{12}\left(\mathbf{x}\right)$

(c) $\sigma_{22}\left(\mathbf{x}\right)$

Figure 17: Comparison of PINN predicted stress fields and corresponding Abaqus
solution for problem defined in section Section 3.6.

###  3.7 2D-FGM-ELAS-DIRCH

While we have seen how to train the PINN model on functionally graded material
with Neumann boundary conditions, let’s consider the same domain and material
properties with a different boundary condition. A uniform displacement
boundary condition is applied on the top edge as shown in the Figure 18. Since
there is no traction, in this case, the loss function will only have the
internal strain energy component (see Table 1).

Figure 18: Domain and boundary conditions for functionally graded material
object under displacement controlled load defined in Section 3.7.

$\mathbf{\tilde{u}(x)}=\begin{bmatrix}\tilde{u}_{1}\\\
\tilde{u}_{2}\end{bmatrix}=\begin{bmatrix}x_{2}\tilde{v}_{1}\\\
x_{2}/3L+(3L-x_{2})x_{2}\tilde{v}_{2}\end{bmatrix}$ (19)

Decrease in loss (internal strain energy) with epoch count as the PINN is
trained is shown in Figure 19. Comparisons between stress and displacement
fields obtained from PINN and using Abaqus simulation are done in Figures 20
and 21. R2 score for this case also was given in Table 3.

Figure 19: Reduction in loss vs epoch as the PINN is trained to predict
elastic response of two-dimensional functionally graded material object with
displacement load defined in Section 3.7.

(a) $u_{1}\left(\mathbf{x}\right)$

(b) $u_{2}\left(\mathbf{x}\right)$

Figure 20: Comparison between PINN prediction displacement fields and
corresponding Abaqus solution for problem defined in Section 3.7.

(a) $\sigma_{11}\left(\mathbf{x}\right)$

(b) $\sigma_{12}\left(\mathbf{x}\right)$

(c) $\sigma_{22}\left(\mathbf{x}\right)$

Figure 21: Comparison between PINN predicted stress fields and corresponding
Abaqus solution for problem defined in Section 3.7.

###  3.8 2D-FGM-THERMO-ELAS

Finally, let us consider a thermo-mechanics problem on 2D domain. The domain
and boundary conditions are shown in the Figure 22. Both elastic modulus
$E(x)$ and conductivity $k(x)$ varies along y-axis such that
$E(x_{1},x_{2})={1}/\left(1+x_{2}\right)$ and
$k(x_{1},x_{2})={10}/\left(1+x_{2}\right)$. The coefficient of thermal
expansion $\alpha$ is assumed to be unity and the reference temperature
$T_{0}$ is assumed to be zero.

Figure 22: Domain and boundary conditions for two-dimensional functionally
graded material object under displacement and temperature load defined in
Section 3.8.

The differential equations governing the thermo-elastic response of object is
given in Table 1. All primary variables $u_{1},u_{2}$ and $T$ are set to zero
on the bottom edge of the plate. The top edge is subjected to vertical
displacement of one unit i.e., $u_{2}(x_{2}=3L)=1$, while temperature of top
edge is fixed at one unit $T(x_{2}=3L)=1$. In order to implement boundary
conditions in PINN framework, the blackbox function $\tilde{v}(\mathbf{x})$ is
modified as:

$\mathbf{u}\left(\mathbf{x}\right)=\begin{bmatrix}\tilde{u}_{1}(\mathbf{x})\\\
\tilde{u}_{2}(\mathbf{x})\\\
\tilde{T}(\mathbf{x})\end{bmatrix}=\begin{bmatrix}x_{2}\tilde{v}_{1}(\mathbf{x})\\\
x_{2}/(3L)+x_{2}(3L-x_{2})\tilde{v}_{2}(\mathbf{x})\\\
x_{2}/(3L)+x_{2}(3L-x_{2})\tilde{v}_{3}(\mathbf{x})\end{bmatrix}$ (20)

The architecture of neural network $\mathbf{v}\left(x\right)$ is given in
Table 2. The loss function used in training the model is given in Table 1. The
decline in loss as the PINN model is trained can be seen in Figure 23.

Figure 23: Reduction in internal strain and thermal energy as the model is
trained to predict thermo-mechanical response of two-dimensional functionally
graded material object defined in Section 3.8.

After training the model, the PINN predictions are compared with solutions
obtained using Abaqus simulation. The comparison between primary fields
$u_{1}(\mathbf{x}),u_{2}(\mathbf{x}),T(\mathbf{x})$ and stress fields
$\sigma_{11}(x),\sigma_{22}(x)$ and $\sigma_{12}(x)$ can be observed in Figure
24 and Figure 25, respectively. The R2 score, for different variables is given
in Table 3.

(a) $u_{1}\left(\mathbf{x}\right)$

(b) $u_{2}\left(\mathbf{x}\right)$

(c) $T\left(\mathbf{x}\right)$

Figure 24: Comparison between PINN predicted displacement and temperature
fields and corresponding Abaqus simulation for problem defined in Section 3.8.

(a) $\sigma_{11}\left(\mathbf{x}\right)$

(b) $\sigma_{12}\left(\mathbf{x}\right)$

(c) $\sigma_{22}\left(\mathbf{x}\right)$

Figure 25: Comparision between PINN predicted stress fields and corresponding
Aabqus solution for problem defined in Section 3.8.

##  4 Conclusions

This study has used PINN to solve solid mechanics and coupled thermo-mechanics
problems of objects with functionally graded properties. Through various one
and two-dimensional examples, the efficacy of PINNs in accurately predicting
primary target variables has been established. On the other hand, it has also
been found that PINNs currently do not perform very well in predicting
secondary variables such as stress fields. When compared to the finite element
simulation method, PINN is advantageous in terms of being meshless and
flexible enough that various changes in the form of material properties,
domain, and boundary conditions can be incorporated with much ease.

At the same time, there are a few drawbacks of the PINN-based method. There is
no knowledge on how to design the neural network given a set of differential
equations. The choice of activation function used in PINN influences the kind
of functions it can fit. Hence, activation functions play a much significant
role in PINN as compared to other neural networks. The PINN framework is also
very sensitive to convergence criteria of losses. It is because a small change
in loss can sometimes result in a drastic change in the obtained solutions.

The current challenges in using PINN only open the door to further avenues of
research. Studying how to design PINN architectures, improving the accuracy of
secondary variables, and consistent convergence criteria will be beneficial.
Research along the lines of transfer learning using PINN and parallelized
training of distributed PINNs will be of great significance.

##  5 Appendix

A neural network is often debugged by looking at various variables such as
gradients, activation values, weights, and biases. The evolution of those
variables while training the model gives us about model architecture’s
efficacy. While training PINN, a slow training speed was observed. We consider
the Kirsch’s problem discussed earlier in the Section 3 and plot the
distribution of various variables to get insight into training. It is
worthwhile to mention that there are a few open-source libraries such as
Tensorboard (https://www.tensorflow.org/tensorboard), which can help in
generating such plots with much ease. Though, in this example, all the plots
are generated manually in python.

In Figure 26, we first look at the distribution of activation values in
different hidden layers as the model is trained. The distribution of
activation values at the start of the training process would be random due to
the random initialization of weights. However, as the training happens, the
distribution should be saturated. Moreover, mean activation close to zero with
a low standard deviation is not desired as it makes the gradients zero and
hampers the training process. In Figure 26, it is observed that mean of
activation values are not zero, which does not impede training. Moreover, the
activation values for Layer 2 and Layer 3 have not saturated till 2000 epochs.

Figure 26: Distribution of activation function values for different layers for
PINN trained on problem defined in Section 3.5. The thick line represents the
mean while the faint bars represents the spread about the mean of activation
values at various epochs.

In Figure 27, we also look at the means and standard deviations of weights and
biases in different layers. We observe that distributions of weights in
different layers have not saturated in 2000 epochs as the standard deviation
keeps on increasing.

Figure 27: Mean and standard deviation of weights and biases as the PINN is
trained for problem defined in Section 3.5

.

We finally look at Figure 28, which shows the distribution of gradients in
different layers. We see that the gradients are non-zero at the very start of
the training process. Hence, the PINN model trains rapidly for the initial few
epochs. However, after a few epochs, the gradients become very close to zero
results in slow training of the model.

Figure 28: Mean and standard deviation of gradients as the PINN trained for
problem defined in Section 3.5. It is to be noted that zero gradients would
result in no further weight updates.

## References

  * Goodfellow et al. [2016] I. Goodfellow, Y. Bengio, A. Courville, Deep Learning, The MIT Press, ISBN 0262035618, 2016\.
  * LeCun et al. [2015] Y. LeCun, Y. Bengio, G. Hinton, Deep learning, Nature 521 (7553) (2015) 436–444, ISSN 1476-4687.
  * Huang and Chen [2020] Y. Huang, Y. Chen, Survey of State-of-Art Autonomous Driving Technologies with Deep Learning, in: 2020 IEEE 20th International Conference on Software Quality, Reliability and Security Companion (QRS-C), 221–228, 2020.
  * Huang et al. [2020] J. Huang, J. Chai, S. Cho, Deep learning in finance and banking: A literature review and classification, Front. Bus. Res. China 14 (1) (2020) 13, ISSN 1673-7431.
  * Tang et al. [2019] B. Tang, Z. Pan, K. Yin, A. Khateeb, Recent Advances of Deep Learning in Bioinformatics and Computational Biology, Front. Genet. 10 (2019) 214, ISSN 1664-8021.
  * Ardabili et al. [2020] S. F. Ardabili, A. Mosavi, P. Ghamisi, F. Ferdinand, A. R. Varkonyi-Koczy, U. Reuter, T. Rabczuk, P. M. Atkinson, Covid-19 outbreak prediction with machine learning, Algorithms 13 (10) (2020) 249.
  * Peeketi et al. [2019] A. R. Peeketi, R. K. Desu, P. Kumbhar, R. K. Annabattula, Thermal analysis of large granular assemblies using a hierarchical approach coupling the macro-scale finite element method and micro-scale discrete element method through artificial neural networks, Comput. Part. Mech. 6 (4) (2019) 811–822.
  * Mosavi et al. [2019] A. Mosavi, M. Salimi, S. Faizollahzadeh Ardabili, T. Rabczuk, S. Shamshirband, A. R. Varkonyi-Koczy, State of the art of machine learning models in energy systems, a systematic review, Energies 12 (7) (2019) 1301\.
  * Masubuchi et al. [2020] S. Masubuchi, E. Watanabe, Y. Seo, S. Okazaki, T. Sasagawa, K. Watanabe, T. Taniguchi, T. Machida, Deep-learning-based image segmentation integrated with optical microscopy for automatically searching for two-dimensional materials, NPJ 2D Mater. Appl. 4 (1) (2020) 3, ISSN 2397-7132.
  * Raj et al. [2021] M. Raj, S. Thakre, R. K. Annabattula, A. K. Kanjarla, Estimation of Local Strain Fields in Two-Phase Elastic Composite Materials Using UNet-Based Deep Learning, Integr. Mater. Manuf. Innov. 10 (3) (2021) 444–460, ISSN 2193-9772.
  * Liu et al. [2021] B. Liu, N. Vu-Bac, T. Rabczuk, A stochastic multiscale method for the prediction of the thermal conductivity of Polymer nanocomposites through hybrid machine learning algorithms, Compos. Struct. 273 (2021) 114269\.
  * Merity et al. [2017] S. Merity, N. S. Keskar, R. Socher, Regularizing and Optimizing LSTM Language Models, CoRR abs/1708.02182.
  * Rumelhart et al. [1986] D. E. Rumelhart, G. E. Hinton, R. J. Williams, Learning representations by back-propagating errors, Nature 323 (6088) (1986) 533–536, ISSN 1476-4687.
  * Hornik et al. [1989] K. Hornik, M. B. Stinchcombe, H. L. White, Multilayer feedforward networks are universal approximators, Neural Networks 2 (1989) 359–366.
  * Oh and Jung [2004] K.-S. Oh, K. Jung, GPU implementation of neural networks, Pattern Recognition 37 (6) (2004) 1311–1314, ISSN 0031-3203.
  * Tan et al. [2018] C. Tan, F. Sun, T. Kong, W. Zhang, C. Yang, C. Liu, A Survey on Deep Transfer Learning, in: V. Kůrková, Y. Manolopoulos, B. Hammer, L. Iliadis, I. Maglogiannis (Eds.), Artificial Neural Networks and Machine Learning – ICANN 2018, Springer International Publishing, Cham, ISBN 978-3-030-01424-7, 270–279, 2018.
  * Raissi et al. [2019] M. Raissi, P. Perdikaris, G. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, J. Comput. Phys. 378 (2019) 686–707, ISSN 0021-9991.
  * Weinan and Bing [2018] E. Weinan, Y. Bing, The Deep Ritz Method: A Deep Learning-Based Numerical Algorithm for Solving Variational Problems, Commun. Math. Stat. 6 (1) (2018) 1–12, ISSN 2194-671X.
  * Anitescu et al. [2019] C. Anitescu, E. Atroshchenko, N. Alajlan, T. Rabczuk, Artificial neural network methods for the solution of second order boundary value problems, Comput. Mater. Contin. 59 (1) (2019) 345–359.
  * Samaniego et al. [2020] E. Samaniego, C. Anitescu, S. Goswami, V. M. Nguyen-Thanh, H. Guo, K. Hamdia, X. Zhuang, T. Rabczuk, An energy approach to the solution of partial differential equations in computational mechanics via machine learning: Concepts, implementation and applications, Comput. Methods Appl. Mech. Eng. 362 (2020) 112790.
  * Goswami et al. [2020a] S. Goswami, C. Anitescu, S. Chakraborty, T. Rabczuk, Transfer learning enhanced physics informed neural network for phase-field modeling of fracture, Theor. Appl. Fract. Mech. 106 (2020a) 102447\.
  * Goswami et al. [2020b] S. Goswami, C. Anitescu, T. Rabczuk, Adaptive fourth-order phase field analysis for brittle fracture, Comput. Methods Appl. Mech. Eng. 361 (2020b) 112808, ISSN 0045-7825.
  * Goswami et al. [2020c] S. Goswami, C. Anitescu, T. Rabczuk, Adaptive fourth-order phase field analysis using deep energy minimization, Theor. Appl. Fract. Mech. 107 (2020c) 102527.
  * Li et al. [2021] W. Li, M. Z. Bazant, J. Zhu, A physics-guided neural network framework for elastic plates: Comparison of governing equations-based and energy-based approaches, Comput. Methods Appl. Mech. Eng. 383 (2021) 113933.
  * Dwivedi and Srinivasan [2020a] V. Dwivedi, B. Srinivasan, Solution of Biharmonic Equation in Complicated Geometries With Physics Informed Extreme Learning Machine, Journal of Computing and Information Science in Engineering 20 (6), ISSN 1530-9827.
  * Dwivedi and Srinivasan [2020b] V. Dwivedi, B. Srinivasan, Physics Informed Extreme Learning Machine (PIELM)–A rapid method for the numerical solution of partial differential equations, Neurocomputing 391 (2020b) 96–118, ISSN 0925-2312.
  * Dwivedi et al. [2019] V. Dwivedi, N. Parashar, B. Srinivasan, Distributed physics informed neural network for data-efficient solution to partial differential equations, 2019.
  * Dwivedi et al. [2021] V. Dwivedi, N. Parashar, B. Srinivasan, Distributed learning machines for solving forward and inverse problems in partial differential equations, Neurocomputing 420 (2021) 299–316, ISSN 0925-2312.
  * Yadav et al. [0] G. K. Yadav, S. Natarajan, B. Srinivasan, Distributed PINN for linear elasticity - a unified approach for smooth, singular, compressible and incompressible media, Int. J. Comput. Methods 2021 (ja) (0) null.
  * Koizumi [1993] M. Koizumi, The concept of FGM, Ceram. Trans. 34 (1993) 3–10.
  * Pindera et al. [1997] M.-J. Pindera, J. Aboudi, A. M. Glaeser, S. M. Arnold, Use of composites in multi-phased and functionally graded materials, Compos. B. Eng. 28 (1-2).
  * Kingma and Ba [2017] D. P. Kingma, J. Ba, Adam: A Method for Stochastic Optimization, 2017.
  * Baydin et al. [2018] A. G. Baydin, B. A. Pearlmutter, A. A. Radul, J. M. Siskind, Automatic Differentiation in Machine Learning: a Survey, J Mach Learn Res. 18 (153) (2018) 1–43.
  * Abadi et al. [2015] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, X. Zheng, TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems, software available from tensorflow.org, 2015.