**A Self-Adaptive Shrinking Projection Method with an Inertial Technique for Split Common Null Point Problems in Banach Spaces**

### **Chibueze Christian Okeke 1, Lateef Olakunle Jolaoso 2,\* and Regina Nwokoye <sup>3</sup>**


Received: 27 October 2020; Accepted: 20 November 2020; Published: 2 December 2020

**Abstract:** In this paper, we present a new self-adaptive inertial projection method for solving split common null point problems in *p*-uniformly convex and uniformly smooth Banach spaces. The algorithm is designed such that its convergence does not require prior estimate of the norm of the bounded operator and a strong convergence result is proved for the sequence generated by our algorithm under mild conditions. Moreover, we give some applications of our result to split convex minimization and split equilibrium problems in real Banach spaces. This result improves and extends several other results in this direction in the literature.

**Keywords:** split common null point; strong convergence; resolvent; metric resolvent; split minimization problem; split equilibrium problem; Banach space

### **1. Introduction**

Let *H*<sup>1</sup> and *H*<sup>2</sup> be real Hilbert spaces and *C* and *Q* be nonempty, closed and convex subsets of *H*<sup>1</sup> and *H*2, respectively. We consider the Split Common Null Point Problem (SCNPP) which was introduced by Byrne et al. [1] as follows:

$$\text{Find } \quad z \in H\_1 \quad \text{such that} \quad z \in A^{-1}(0) \bigcap T^{-1}(B^{-1}(0)), \tag{1}$$

where *<sup>A</sup>* : *<sup>H</sup>*<sup>1</sup> <sup>→</sup> <sup>2</sup>*H*<sup>1</sup> and *<sup>B</sup>* : *<sup>H</sup>*<sup>2</sup> <sup>→</sup> <sup>2</sup>*H*<sup>2</sup> are maximal monotone operators and *<sup>T</sup>* : *<sup>H</sup>*<sup>1</sup> <sup>→</sup> *<sup>H</sup>*<sup>2</sup> is a linear bounded operator. The solution set of SCNPP (1) is denoted by Ω. The SCNPP contains several important optimization problems such as split feasibility problem, split equilibrium problem, split variational inequalities, split convex minimization problem, split common fixed point problems, etc., as special cases (see, e.g., [1–5]). Due to their importance, several researchers have studied and proposed various iterative methods for finding its solutions (see, e.g., [1,4–9]). In particular, Byrne et al. [1] introduced the following iterative scheme for solving SCNPP in real Hilbert spaces:

$$\begin{cases} \mathbf{x}\_0 \in H\_1, \lambda > 0, \\ \mathbf{x}\_{n+1} = f\_\lambda^A(\mathbf{x}\_n + \lambda T^\*(f\_\lambda B)T\mathbf{x}\_n), \quad n \ge 0, \end{cases} \tag{2}$$

where *J<sup>A</sup> <sup>λ</sup> <sup>x</sup>* = (*<sup>I</sup>* <sup>+</sup> *<sup>λ</sup>A*)−1*x*, for all *<sup>x</sup>* <sup>∈</sup> *<sup>H</sup>*1. They also proved that the sequence {*xn*} generated by (2) converges weakly to a solution of SCNPP provided the step size *λ* satisfies

$$
\lambda \in \left(0, \frac{2}{L}\right),
$$

where *L* is the spectral radius of *T*. Furthermore, Kazmi and Rizvi [10] proposed a viscosity method which converges strongly to a solution of (1) as follows:

$$\begin{cases} \mathbf{x}\_{\text{0}} \in H\_{1}, \lambda > 0, \\ u\_{n} = J\_{\lambda}^{A} (\mathbf{x}\_{n} + \lambda \mathbf{T}^{\*} (I\_{\lambda}^{B} - I) A \mathbf{x}\_{n}), \\ \mathbf{x}\_{n+1} = a\_{\text{H}} f(\mathbf{x}\_{\text{H}}) + (1 - a\_{\text{H}}) S u\_{\text{H}}, \quad n \ge 0, \end{cases} \tag{4}$$

where {*αn*} ⊂ (0, 1) satisfies some certain conditions and *S* : *H*<sup>1</sup> → *H*<sup>1</sup> is a nonexpansive mapping. It is important to emphasize that the convergence of (4) is achieved with the aid of condition (3). Other similar results can be found, for instance, in [11,12] (and references therein). However, it is well known that the norm of bounded linear operator is very difficult to find (or at least estimate) (see [13–15]). Hence, it becomes necessary to find iterative methods whose step size selection does not require prior estimate of the norm of the bounded linear operator. Recently, some authors have provided breakthrough results in the framework of real Hilbert spaces (see, e.g., [13–15]).

On the other hand, Takahashi [8,16] extends the study of SCNPP (1) to uniformly convex and smooth Banach spaces as follows: Let *E*<sup>1</sup> and *E*<sup>2</sup> be uniformly convex and uniformly smooth real Banach spaces with dual *E*∗ <sup>1</sup> and *E*<sup>∗</sup> <sup>2</sup> , respectively, and *T* : *E*<sup>1</sup> → *E*<sup>2</sup> be a bounded linear operator. Let *<sup>A</sup>* : *<sup>E</sup>*<sup>1</sup> <sup>→</sup> <sup>2</sup>*E*<sup>∗</sup> <sup>1</sup> and *<sup>B</sup>* : *<sup>E</sup>*<sup>2</sup> <sup>→</sup> <sup>2</sup>*E*<sup>∗</sup> <sup>2</sup> be maximal monotone operators such that *<sup>A</sup>*−1(0) <sup>=</sup> <sup>∅</sup>, *<sup>B</sup>*−1(0) <sup>=</sup> <sup>∅</sup> and *<sup>Q</sup><sup>μ</sup>* is a metric resolvent operator with respect to *<sup>B</sup>* and parameter *<sup>μ</sup>* <sup>&</sup>gt; 0. Takahashi and Takahashi [17] introduced the following shrinking projection method for solving SCNPP in uniformly convex and smooth Banach spaces:

$$\begin{cases} \mathbf{x}\_1 \in \mathbb{C}, \mu\_1 > 0, \\ \mathbf{z}\_n = \mathbf{x}\_n - f\_{\lambda\_n} f\_{E\_1}^{-1} T^\* f\_{E\_2} (T \mathbf{x}\_n - Q\_{\mu\_n} T \mathbf{x}\_n), \\ \mathbf{C}\_{n+1} = \{ \mathbf{z} \in \mathbb{C}\_n : \langle z\_n - z, f\_{E\_1} (\mathbf{x}\_n - z\_n) \rangle \ge 0 \}, \\ \mathbf{x}\_{n+1} = P\_{\mathbb{C}\_{n+1}} \mathbf{x}\_1, \quad \text{for all } n \in \mathbb{N}\_\prime \end{cases} \tag{5}$$

where *JEi* are the normalized duality mapping with respect to *Ei* for *i* = 1, 2 (defined in the next section). They proved a strong convergence result with the condition that the step size satisfies

$$0 < a \le \lambda\_n ||T||^2 < b < 1 \quad \text{and} \quad 0 < c \le \mu\_n \quad \text{for all } n \in \mathbb{N}.$$

Furthermore, Suantai et al. [18] introduced a new iterative scheme for solving SCNPP in a real Hilbert space *H* and a real Banach space *E* as follows:

$$\begin{cases} \mathbf{x}\_{1} \in H, \\ y\_{n} = f\_{\lambda\_{n}}^{A} (\mathbf{x}\_{n} + \lambda\_{n} T^{\*} f\_{E} (Q\_{\mu\_{n}} - I) \, \mathbf{x}\_{n}), \\ \mathbf{x}\_{n+1} = \mathbf{a}\_{n} f(\mathbf{x}\_{n}) + \beta\_{n} \mathbf{x}\_{n} + \gamma\_{n} y\_{n}, \quad n \ge 1, \end{cases} \tag{6}$$

where {*αn*}, {*βn*}, {*γn*} ⊂ (0, 1) such that *α<sup>n</sup>* + *β<sup>n</sup>* + *γ<sup>n</sup>* = 1 and *f* : *H* → *H* is a contraction mapping. They also proved a strong convergence result under the condition that the step size satisfies

$$0 < \lambda\_n \|\|T\|\|^2 < 2.$$

*Axioms* **2020**, *9*, 140

Recently, Takahashi [19] introduced a new hybrid method with generalized resolvent operators for solving the SCNPP in real Banach spaces as follows:

$$\begin{cases} \boldsymbol{z}\_{n} = \boldsymbol{J}^{-1}(f\_{E}\mathbf{x}\_{n} - \boldsymbol{r}\_{n}T^{\*}\left(f\_{F}T\mathbf{x}\_{n} - f\_{F}Q\_{\mu\_{n}}T\mathbf{x}\_{n}\right)),\\ \boldsymbol{y}\_{n} = f\_{\lambda\_{n}}\boldsymbol{z}\_{n},\\ \boldsymbol{\mathsf{C}}\_{n} = \{\boldsymbol{z} \in \boldsymbol{E} : 2\langle\mathbf{x}\_{n} - \boldsymbol{z}, f\_{E}\mathbf{x}\_{n} - f\_{E}\mathbf{z}\_{n}\rangle \geq r\_{n}\rho\_{F}(T\mathbf{x}\_{n}, Q\_{\mu\_{n}}T\mathbf{x}\_{n})\},\\ \boldsymbol{D}\_{n} = \{\boldsymbol{z} \in \boldsymbol{E} : \langle\boldsymbol{y}\_{n} - \boldsymbol{z}, f\_{E}\mathbf{z}\_{n} - f\_{E}\mathbf{y}\_{n}\rangle \geq 0\},\\ \boldsymbol{Q}\_{n} = \{\boldsymbol{z} \in \boldsymbol{E} : \langle\mathbf{x}\_{n} - \boldsymbol{z}, f\_{E}\mathbf{x}\_{1} - f\_{E}\mathbf{x}\_{n}\rangle \geq 0\},\\ \boldsymbol{x}\_{n+1} = \boldsymbol{II}\_{\mathbb{C}\_{n}\cap\mathbb{D}\_{N}}\mathbf{x}\_{1}, \text{ for all } n \in \mathbb{N}. \end{cases} \tag{7}$$

He also proved that the sequence generated by Algorithm (7) converges strongly to a solution of SCNPP provided the step sizes satisfy

$$0 < a \le r\_n \le \frac{1}{||T||^{2^\prime}} \quad \text{and} \ 0 < b \le \lambda\_n, \mu\_n \quad \text{for all } n \in \mathbb{N}.$$

It is evident that the above methods and other similar ones (see, e.g., [6,9,20]) require prior knowledge of the operator norm, which is very difficult to find. Thus, the following natural question arises.

**Problem 1.** *Can we provide a new iterative method for solving SCNPP in real Banach spaces such that the step size does not require prior estimate of the norm of the bounded linear operator?*

Let us also mention the inertial extrapolation process which is considered as a means of speeding up the rate of convergence of iterative methods. This technique was first introduced by Polyak [21] as a heavy-ball method of a two-order time dynamical system and has been employed by many authors recently (see, e.g., [22–27]). Moreover, Dong et al. [27] introduced a modified inertial hybrid algorithm for approximating the fixed points of non-expansive mappings in real Hilbert spaces as follows:

$$\begin{cases} \mathbf{x}\_{0}, \mathbf{x}\_{1} \in \mathbb{C}, \\ w\_{n} = \mathbf{x}\_{n} + \theta\_{n} (\mathbf{x}\_{n} - \mathbf{x}\_{n-1}), \\ z\_{n} = (1 - \beta\_{n}) w\_{n} + \beta\_{n} T w\_{n}, \\ \mathbb{C}\_{n} = \{\mathbf{x} \in \mathbb{C} : ||z\_{n} - \mathbf{x}||^{2} \le ||\mathbf{x}\_{n} - \mathbf{x}||^{2}\}, \\ Q\_{n} = \{\mathbf{x} \in \mathbb{C} : \langle \mathbf{x}\_{n} - \mathbf{x}, \mathbf{x}\_{0} - \mathbf{x}\_{n} \rangle \ge 0\}, \\ \mathbf{x}\_{n+1} = P\_{\mathbb{C}\_{n} \cap Q\_{n}} \mathbf{x}\_{0}. \end{cases} \tag{8}$$

where {*θn*} ⊂ [*a*1, *a*2], *a*<sup>1</sup> ∈ (−∞, 0], *a*<sup>2</sup> ∈ [0, +∞), {*βn*} ⊂ (0, 1) are suitable parameters.

More recently, Cholamjiak et al. [28] introduced an inertial forward-backward algorithm for finding the zeros of sum of two monotone operators in Hilbert spaces as follows:

$$\begin{cases} \mathbf{x}\_{0}, \mathbf{x}\_{1} \in H, \mathbf{r}\_{n} > 0, \\ \mathbf{y}\_{n} = \mathbf{x}\_{n} + \theta\_{n} (\mathbf{x}\_{n} - \mathbf{x}\_{n-1}), \\ z\_{n} = \alpha\_{n} y\_{n} + (1 - \alpha\_{n}) \mathbf{T} y\_{n}, \\ \mathbf{z}\_{n} = \beta\_{n} z\_{n} + (1 - \beta\_{n}) I\_{r\_{n}}^{B} (I - r\_{n} A) z\_{n\prime} \\ \mathbf{C}\_{n+1} = \{\mathbf{v} \in \mathbf{C}\_{n} : \|\mathbf{v}\_{n} - \mathbf{v}\|^{2} \le \|\mathbf{x}\_{n} - \mathbf{v}\|^{2} + \mathbf{K}\_{n} \}, \\ \mathbf{x}\_{n+1} = P\_{\mathbf{C}\_{n+1}} \mathbf{x}\_{1}, n \ge 1, \end{cases} \tag{9}$$

where *Kn* = 2*θ*<sup>2</sup> *<sup>n</sup>xn* <sup>−</sup> *xn*−1 − <sup>2</sup>*θnxn* <sup>−</sup> *<sup>z</sup>*, *xn*−<sup>1</sup> <sup>−</sup> *xn*, *<sup>J</sup><sup>B</sup> rn* = (*<sup>I</sup>* <sup>+</sup> *rnB*)−1, {*θn*} ⊂ [0, *<sup>θ</sup>*] for some *θ* ∈ [0, 1) and {*αn*}, {*βn*} are sequences in [0, 1]. The authors proved that the sequence {*xn*} generated by (9) converges strongly to a solution *<sup>x</sup>* <sup>∈</sup> (*<sup>A</sup>* <sup>+</sup> *<sup>B</sup>*)−1(0) under some mild conditions.

Motivated by the above results, in this paper, we aim to provide an affirmative answer to Problem 1. We introduce a new inertial shrinking projection method for solving SCNPP in *p*-uniformly convex and uniformly smooth real Banach spaces. The algorithm is designed such that its step size is determined by a self-adaptive technique and its convergence does not require prior knowledge of the norm of the bounded operator. We also prove a strong convergence result and provide some applications of our main theorem to solving other nonlinear optimization problems. This result improves and extends the results in [6,8,9,11,12,16,19,20] and many other recent results in the literature.

### **2. Preliminaries**

Let *E* be a real Banach space with dual *E*<sup>∗</sup> and norm ·. We denote the duality pairing between *f* ∈ *E* and *g*<sup>∗</sup> ∈ *E*<sup>∗</sup> as *f* , *g*∗. The weak and strong convergence of {*xn*} ⊂ *E* to *a* ∈ *E* are denoted by *xn a* and *xn* → *a*, respectively, ∀ by "for all" and ⇔ by "if and only if". The function *δ<sup>E</sup>* : [0, 2] → [0, 1] defined by

$$\delta\_E(\alpha) = \inf \left\{ 1 - \frac{\|f + \mathbf{g}\|}{2} : \|f\| = 1 = \|\mathbf{g}\|\_{\prime} \|f - \mathbf{g}\| \ge a \right\}$$

is called the modulus of convexity of *E*. The Banach space *E* is said to be uniformly convex if *δE*(*α*) > 0. If there exists a constant *Cp* <sup>&</sup>gt; 0 such that *<sup>δ</sup>E*(*α*) <sup>≥</sup> *Cpα<sup>p</sup>* for any *<sup>α</sup>* <sup>∈</sup> (0, 2], then we say *<sup>E</sup>* is *<sup>p</sup>*-uniformly convex. In addition, the function *ρE*(*β*) : [0, ∞) → [0, +∞) defined by

$$\rho\_E(\beta) = \left\{ \frac{||f + \beta \mathbf{g}|| + ||f - \beta \mathbf{g}||}{2} - 1 : ||f|| = ||\mathbf{g}|| = 1 \right\}$$

is called the modulus of smoothness of *E*. The Banach space *E* is said to be uniformly smooth if lim*β*→+<sup>∞</sup> *<sup>ρ</sup>E*(*β*) *<sup>β</sup>* <sup>=</sup> 0. If there exists a constant *Dq* <sup>&</sup>gt; 0 such that *<sup>ρ</sup>E*(*β*) <sup>≤</sup> *Dqβ<sup>q</sup>* for any *<sup>β</sup>* <sup>&</sup>gt; 0, then *<sup>E</sup>* is called *<sup>q</sup>*-uniformly smooth Banach space. Let 1 <sup>&</sup>lt; *<sup>q</sup>* <sup>≤</sup> <sup>2</sup> <sup>≤</sup> *<sup>p</sup>* with <sup>1</sup> *<sup>p</sup>* <sup>+</sup> <sup>1</sup> *<sup>q</sup>* = 1. We Remark that a Banach space *E* is *p*-uniformly convex if and only if its dual *E*∗ is *q*-uniformly smooth. Examples of *q*-uniformly smooth Banach spaces include Hilbert spaces, *Lq*(or *lp*) spaces, 1 < *p* < ∞ and the Sobolev spaces, *W<sup>p</sup> <sup>m</sup>*, 1 < *p* < ∞ (see [29]). Moreover, the Hilbert spaces are uniformly smooth while

$$L\_p(or \, l\_p) \text{ or } \mathcal{W}\_m^p \text{ is } \begin{cases} p-\text{uniformly smooth if } 1 < p \le 2 \\ 2-\text{uniformly smooth if } p \ge 2. \end{cases}$$

Let *<sup>ϕ</sup>* : <sup>R</sup><sup>+</sup> <sup>→</sup> <sup>R</sup><sup>+</sup> be a continuous strictly increasing function. *<sup>ϕ</sup>* is called a gauge function if

$$
\varphi(0) = 0, \quad \lim\_{t \to \infty} \varphi(t) = +\infty.
$$

The duality mapping with respect to *ϕ*, i.e., *J<sup>ϕ</sup>* : *E* → *E*<sup>∗</sup> is defined by

$$J\_{\mathfrak{g}}(\mathbf{x}) = \{ j \in E^\* : \langle \mathbf{x}, j \rangle = ||\mathbf{x}|| ||j||\_{\*\*} ||j||\_\* = \mathfrak{g}(||\mathbf{x}||) \}, \ \mathbf{x} \in E.$$

When *ϕ*(*t*) = *t*, then we call *J<sup>ϕ</sup>* = *J* a normalized duality mapping. In addition, if *ϕ*(*t*) = *tp*−<sup>1</sup> where *p* > 1, then, *J<sup>ϕ</sup>* = *Jp* is called a generalized duality mapping defined by

$$J\_p(\mathfrak{u}) = \{ f \in E^\* : \langle \mathfrak{u}, f \rangle = ||\mathfrak{u}|| \| f \|\_{\ast \ast} \| f \|\_{\ast} = ||\mathfrak{u}||^{p-1} \}, \ \mathfrak{x} \in E.$$

In the sequel, *C* is a nonempty closed convex subset of *E* and *F*(*T*) = {*x* ∈ *C* : *Tx* = *x*} is the set of fixed point of *T* : *C* → *C*.

**Definition 1.** *Ref. [30] Let E be a Banach space, J<sup>ϕ</sup>* : *E* → *E*<sup>∗</sup> *a duality mapping with gauge function ϕ*, *and C a nonempty subset of E*. *A mapping T* : *C* → *E is said to be*

*(i) ϕ-firmly non-expansive if*

$$\langle \!\langle Tu - Tv, f\_{\mathfrak{P}}(Tu) - f\_{\mathfrak{P}}(Tv) \rangle \rangle \leq \langle Tu - Tv, f\_{\mathfrak{P}}(u) - f\_{\mathfrak{P}}(v) \rangle \rangle$$

*for all u*, *v* ∈ *C*.

*(ii) ϕ-firmly quasi-non-expansive if F*(*T*) = ∅ *and*

$$\langle Tu - z\_\prime J\_\varphi(u) - J\_\varphi(Tu) \rangle \ge 0$$

*for all u in C and z in F*(*T*).

**Definition 2.** *Given a Gâteaux differentiable and convex function f* : *<sup>E</sup>* <sup>→</sup> <sup>R</sup>, *the function*

$$\Delta\_f(u,\upsilon) := f(\upsilon) - f(u) - \langle f'(u), \upsilon - u \rangle, \text{ for all } u, \upsilon \in E \tag{10}$$

*is called the Bregman distance of u to v with respect to the function f* .

Moreover, since *J p <sup>E</sup>* is the derivative of the function *fp*(*u*) = <sup>1</sup> *<sup>p</sup> up*, in that case, the Bregman distance with respect to *fp* becomes

$$\begin{aligned} \Delta\_p(\boldsymbol{\mu}, \boldsymbol{\upsilon}) &= \quad \frac{1}{q} ||\boldsymbol{\mu}||^p - \langle f\_E^p \boldsymbol{\mu}, \boldsymbol{\upsilon} \rangle + \frac{1}{p} ||\boldsymbol{\upsilon}||^p \\ &= \quad \frac{1}{p} (||\boldsymbol{\upsilon}||^p - ||\boldsymbol{\mu}||^p) + \langle f\_E^p \boldsymbol{\mu}, \boldsymbol{\iota} - \boldsymbol{\upsilon} \rangle \\ &= \quad \frac{1}{q} (||\boldsymbol{\mu}||^p - ||\boldsymbol{\upsilon}||^p) - \langle f\_E^p \boldsymbol{\mu} - f\_E^p \boldsymbol{\upsilon}, \boldsymbol{\upsilon} \rangle. \end{aligned}$$

**Remark 1.** *It follows from the Definition of* Δ*<sup>p</sup> that*

$$
\Delta\_{\mathcal{P}}(u,\upsilon) = \Delta\_{\mathcal{P}}(u,z) + \Delta\_{\mathcal{P}}(z,\upsilon) + \langle z - \upsilon, l\_{\mathcal{E}}^{p}u - l\_{\mathcal{E}}^{p}z \rangle, \text{ for all } u,\upsilon,z \in \mathcal{E},\tag{11}
$$

*and*

$$
\Delta\_p(u, \upsilon) + \Delta\_p(\upsilon, u) = \langle u - \upsilon, \mathbb{J}\_\mathbb{E}^p u - \mathbb{J}\_\mathbb{E}^p \upsilon \rangle, \text{ for all } u, \upsilon, z \in \mathbb{E}. \tag{12}
$$

Although the Bregman is not symmetrical, it however has the following relationship with · distance:

$$\|u(\|u-\upsilon\|)^p \le \Delta\_p(u,\upsilon) \le \langle u-\upsilon, l\_{\to}^p u - l\_{\to}^p \upsilon \rangle\_{\prime} \quad \text{for all } u,\upsilon \in E,\ u > 0. \tag{13}$$

This indicates that Bregman distance is non-negative.

**Definition 3.** *The Bregman projection mapping* Π*<sup>C</sup>* : *E* → *C is defined by*

$$\Pi\_{\mathbb{C}} u = \arg\min\_{v \in \mathbb{C}} \Delta\_p(u, v), \text{ for all } u \in E. \tag{14}$$

The Bregman projection can also be characterized by the following inequality

$$
\langle f\_E^p u - f\_E^p \Pi\_\mathbb{C} u\_\prime z - \Pi\_\mathbb{C} u \rangle \le 0,\text{ for all } z \in \mathbb{C},\tag{15}
$$

This is equivalent to

$$
\Delta\_{\mathbb{P}}(\Pi\_{\mathbb{C}}u, z) \le \Delta\_{\mathbb{P}}(u, z) - \Delta\_{\mathbb{P}}(u, \Pi\_{\mathbb{C}}u), \text{ for all } z \in \mathbb{C}. \tag{16}
$$

**Lemma 1.** *Ref. [31] Let E be a q-uniformly smooth Banach space with q-uniformly smoothness constant cq* > 0. *For any u*, *v* ∈ *E*, *the following inequality holds:*

$$||\mu - \upsilon||^q \le ||\mu||^q - q\langle \upsilon, f\_E^q \mu \rangle + c\_{\eta}||\upsilon||^q.$$

**Definition 4.** *A mapping T* : *C* → *C is said to be closed or has a closed graph if a sequence* {*xn*} ⊂ *C converges strongly to a point x* ∈ *C and Txn* → *y*, *then Tx* = *y*.

**Lemma 2.** *Ref. [29] It is known that the generalized duality has the following properties:*


**Lemma 3.** *Ref. [32] For any* {*xn*} ⊂ *<sup>E</sup>*, {*tn*} ⊂ (0, 1) *with* <sup>∑</sup>*<sup>N</sup> <sup>n</sup>*=<sup>1</sup> *tn* = 1, *the following inequality holds:*

$$\Delta\_{\mathcal{P}}(f\_{\mathbb{E}^{\mathbf{x},\prime}}^{\mathcal{G}}(\sum\_{n=1}^{N}t\_{\mathcal{U}}f\_{\mathbb{E}}^{\mathcal{P}}(\mathbf{x}\_{\mathbb{H}})),\mathbf{x}) \leq \sum\_{n=1}^{N}t\_{\mathcal{U}}\Delta\_{\mathcal{P}}(\mathbf{x}\_{\mathbb{H}},\mathbf{x}) \text{ for all } \mathbf{x} \in E.$$

We now define some important operators which play key role in our convergence analysis.

**Definition 5.** *Let <sup>A</sup>* : *<sup>E</sup>* <sup>→</sup> <sup>2</sup>*E*<sup>∗</sup> *be a multi-valued mapping. We define the effective domain of A by* D(*A*) = {*<sup>x</sup>* ∈ *<sup>E</sup>* : *Ax* = <sup>0</sup>} *and range of <sup>A</sup> by* (*A*) = *<sup>x</sup>*∈D(*A*) *Ax*. *The operator <sup>A</sup> is said to be monotone if <sup>x</sup>* <sup>−</sup> *<sup>y</sup>*, *<sup>u</sup>*<sup>∗</sup> <sup>−</sup> *<sup>v</sup>*∗ ≥ <sup>0</sup> *for all <sup>x</sup>*, *<sup>y</sup>* <sup>∈</sup> <sup>D</sup>(*A*), *<sup>u</sup>*<sup>∗</sup> <sup>∈</sup> *Ax and <sup>v</sup>*<sup>∗</sup> <sup>∈</sup> *Ay*. *When the graph of <sup>A</sup> is not properly contained in the graph of any other monotone operator, then we say that A is maximally monotone.*

*Let <sup>E</sup> be a smooth, strictly convex, and reflexive Banach space and <sup>A</sup>* : *<sup>E</sup>* <sup>→</sup> <sup>2</sup>*E*<sup>∗</sup> *be a maximal monotone operator. The metric resolvent operator with respect to A is defined by Q<sup>ϕ</sup> <sup>r</sup>* (*u*)=(*I* + *r J*−<sup>1</sup> *<sup>ϕ</sup> A*)−1(*u*). *It is easy to see that*

$$0 \in J\_{\boldsymbol{\theta}}(\boldsymbol{Q}\_{r}^{\boldsymbol{\theta}}(\boldsymbol{u}) - \boldsymbol{u}) + r A \boldsymbol{Q}\_{r}^{\boldsymbol{\theta}}(\boldsymbol{u}),\tag{17}$$

*and F*(*Q<sup>ϕ</sup> <sup>r</sup>* ) = *A*−10 *for all r* > 0 *(see, e.g., [20]). Moreover, by the monotonicity of A*, *we can show that*

$$<\langle Q\_r^{\mathfrak{g}}(\mathfrak{u}) - Q\_r^{\mathfrak{g}}(\mathfrak{v}), I\_{\mathfrak{g}}(\mathfrak{u} - Q\_r^{\mathfrak{g}}(\mathfrak{u})) - I\_{\mathfrak{g}}(\mathfrak{v} - Q\_r^{\mathfrak{g}}(\mathfrak{v})) \rangle \ge 0 \tag{18}$$

*for all u*, *<sup>v</sup>* <sup>∈</sup> *<sup>E</sup>*. *In addition, if A*−10 <sup>=</sup> <sup>∅</sup>, *then*

$$
\langle Q\_r^{\varphi}(u) - z, l\_{\varphi}(u - Q\_r^{\varphi}(u)) \rangle \ge 0 \tag{19}
$$

*for all <sup>u</sup>* <sup>∈</sup> *<sup>E</sup> and <sup>z</sup>* <sup>∈</sup> *<sup>A</sup>*−10. *In the case <sup>ϕ</sup>*(*t*) = *<sup>t</sup>p*−<sup>1</sup> *with <sup>p</sup>* <sup>∈</sup> (1, <sup>+</sup>∞), *we denote <sup>Q</sup><sup>ϕ</sup> <sup>r</sup> by Qr* = (*I* + *r J*−<sup>1</sup> *<sup>p</sup> A*)−<sup>1</sup> *(see, e.g., [33]).*

**Proposition 1.** *Ref. [30] Let A* : *<sup>E</sup>* <sup>→</sup> <sup>2</sup>*E*<sup>∗</sup> *be an operator satisfying the following range condition*

$$\mathbb{D}(A) \subset \mathbb{C} \subset J\_{\varphi}^{-1} \mathfrak{R}(J\_{\varphi} + \lambda A) \text{ for all } \lambda > 0.$$

*Define the ϕ-resolvent operator R<sup>ϕ</sup> <sup>λ</sup>* : *<sup>C</sup>* <sup>→</sup> <sup>2</sup>*<sup>E</sup> associated with operator A by*

$$R^{\mathfrak{g}}\_{\lambda}(\mathfrak{x}) = \{ z \in X : f\_{\mathfrak{g}}(\mathfrak{x}) \in (f\_{\mathfrak{g}} + \lambda A)z \}, \quad \mathfrak{x} \in \mathbb{C}.$$

*Then, for any u* ∈ *C and λ* > 0, *we see that*

$$\begin{aligned} 0 \in Au \quad \Leftrightarrow \quad & f\_{\varPhi}(u) \in (f\_{\varPhi} + \lambda A)u \\ \iff & u \in (f\_{\varPhi} + \lambda A)^{-1} f\_{\varPhi}(u) \\ \iff & u \in F(R\_{\varlambda}^{\operatorname{op}}). \end{aligned}$$

**Proposition 2.** *Ref. [30] Let C be a nonempty, closed, and convex subset of a reflexive, strictly convex Banach space <sup>E</sup> and let <sup>J</sup><sup>ϕ</sup>* : *<sup>E</sup>* <sup>→</sup> *<sup>E</sup>*<sup>∗</sup> *be the duality mapping with gauge <sup>ϕ</sup>*. *Let <sup>A</sup>* : *<sup>E</sup>* <sup>→</sup> <sup>2</sup>*E*<sup>∗</sup> *be a monotone operator satisfying the condition* <sup>D</sup> <sup>⊂</sup> *<sup>C</sup>* <sup>⊂</sup> *<sup>J</sup>*−<sup>1</sup> *<sup>ϕ</sup>* (*J<sup>ϕ</sup>* <sup>+</sup> *<sup>λ</sup>A*), *where <sup>λ</sup>* <sup>&</sup>gt; 0. *Let R<sup>ϕ</sup> <sup>λ</sup> be a resolvent operator of A*; *then,*


Let *E* be a uniformly convex and smooth Banach space. Let *A* be a monotone operator of *E* into 2*E*<sup>∗</sup> . From Browder [34], we know that *A* is maximal if and only if, for any *r* > 0,

$$
\Re(f\_\theta + rA) = E^\*.
$$

### **Remark 2.**

*(i) The smoothness and strict convexity of E ensures that Rϕ*,*<sup>A</sup> <sup>λ</sup> is single-valued. In addition, the range condition ensure that R<sup>ϕ</sup> <sup>λ</sup> single-valued operator from C into* <sup>D</sup>(*A*). *In other words,*

$$R\_{\lambda}(\mathfrak{x})^{\mathfrak{q}}(\mathfrak{x}) = (I\_{\mathfrak{q}} + \lambda A)^{-1} I\_{\mathfrak{q}}(\mathfrak{x}), \text{ for all } \mathfrak{x} \in \mathbb{C}.$$

*(ii) When A is maximal monotone, the range condition holds for C* = D(*A*).

In the sequel, we denote *R<sup>ϕ</sup> <sup>λ</sup>* by *<sup>R</sup><sup>λ</sup>* = (*Jp* + *<sup>λ</sup>A*)−<sup>1</sup> *Jp* for convenience.

Let *E* and *F* be real Banach spaces and let *T* : *E* → *F* be a bounded linear. The dual (adjoint) operator of *T*, denoted by *T*∗, is a bounded linear operator defined by *T*<sup>∗</sup> : *F*<sup>∗</sup> → *E*<sup>∗</sup>

$$\langle T^\*\vec{y}, x \rangle := \langle \vec{y}, Tx \rangle, \quad \text{for all } x \in E\_\prime \vec{y} \in F^\*$$

and the equalities *T*∗ = *T* and ℵ(*T*∗) = (*T*)<sup>⊥</sup> are valid, where (*T*)<sup>⊥</sup> := {*x*<sup>∗</sup> ∈ *F*<sup>∗</sup> : *x*∗, *u* = 0, for all *u* ∈ (*T*)} (see [35,36] for more details on bounded linear operators and their duals).

**Lemma 4.** *Ref. [9] Let E and F be uniformly convex and smooth Banach spaces, Let T* : *E* → *F be a bounded linear operator with the adjoint operator T*∗. *Let R<sup>λ</sup> be the resolvent operator associated with a maximal monotone operator A on E and let Qr be a metric resolvent associated with a maximal monotone operator B on F*. *Assume that A*−10 <sup>∩</sup> *<sup>T</sup>*−1(*B*−10) <sup>=</sup> <sup>∅</sup>. *Let <sup>λ</sup>*, *<sup>μ</sup>*,*<sup>r</sup>* <sup>&</sup>gt; <sup>0</sup> *and z* <sup>∈</sup> *<sup>E</sup>*. *Then, the following are equivalent:*

*(a) z* = *Rλ*(*J q <sup>E</sup>*<sup>∗</sup> (*J p <sup>E</sup>*(*z*) − *μT*<sup>∗</sup> *J p <sup>F</sup>* (*Tz* − *Q* − *rTz*))); *and (b) z* <sup>∈</sup> *<sup>A</sup>*−<sup>0</sup> <sup>∩</sup> *<sup>T</sup>*−1(*B*−10).

### **3. Main Results**

In this section, we present our algorithm and its convergence analysis. In the sequel, we assume that the following assumption hold.


In addition, we denote by *J p <sup>E</sup>*<sup>1</sup> and *J p <sup>E</sup>*<sup>2</sup> the duality mappings of *E*<sup>1</sup> and *E*2, respectively, while *J q E*∗ 1 is the duality mapping of *E*∗ <sup>1</sup> . It is worth mentioning that, when *E*<sup>∗</sup> <sup>1</sup> and *E*<sup>∗</sup> <sup>2</sup> are two *q*-uniformly smooth and uniformly convex Banach spaces, *J p <sup>E</sup>*<sup>1</sup> = (*<sup>J</sup> q E*∗ 1 )−<sup>1</sup> where 1 <sup>&</sup>lt; *<sup>q</sup>* <sup>≤</sup> <sup>2</sup> <sup>≤</sup> *<sup>p</sup>* <sup>&</sup>lt; <sup>+</sup><sup>∞</sup> with <sup>1</sup> *<sup>p</sup>* <sup>+</sup> <sup>1</sup> *<sup>q</sup>* = 1.

**Algorithm SASPM:** Given initial values *x*0, *x*<sup>1</sup> ∈ *C*<sup>1</sup> = *E*1, the sequence {*xn*} generated by the following iterative algorithm:

$$\begin{cases} w\_n = f\_{\mathcal{E}\_1^\*}^q \left[ f\_{\mathcal{E}\_1}^p x\_n + \theta\_n f\_{\mathcal{E}\_1}^p (x\_n - x\_{n-1}) \right], \\ z\_n = f\_{\mathcal{E}\_1^\*}^q \left[ f\_{\mathcal{E}\_1}^p (w\_n) - \rho\_n \frac{f^{p-1} (w\_n)}{\| \mathcal{T}^\* (f\_{\mathcal{E}\_2}^p (Tw\_n - Q\_{r\_n} Tw\_n)) \|^p} T^\* f\_{\mathcal{E}\_1}^p (Tw\_n - Q\_{r\_n} Tw\_n) \right], \\ y\_n = f\_{\mathcal{E}\_1}^q \left( a\_n f\_{\mathcal{E}\_1}^p z\_n + (1 - a\_n) f\_{\mathcal{E}\_1}^p R\_{\lambda\_n} z\_n \right), \\ \mathcal{C}\_{n+1} = \{ u \in \mathcal{C}\_n : \Delta\_p(y\_n, u) \le \Delta\_p(z\_n, u) \le \Delta\_p(w\_n, u) \}, \\ x\_{n+1} = \Pi\_{\mathcal{C}\_{n+1}} x\_0 \end{cases} \tag{20}$$

where {*rn*}, {*λn*} ⊂ (0, ∞), Π*Cn*+<sup>1</sup> is a Bregman projection of *E*<sup>1</sup> onto *Cn*+1, the sequence of real number {*αn*} ⊂ [*a*, *<sup>b</sup>*] <sup>⊂</sup> (0, 1) and {*θn*} ⊂ [*c*, *<sup>d</sup>*] <sup>⊂</sup> (−∞, <sup>+</sup>∞), *<sup>f</sup>*(*wn*) :<sup>=</sup> <sup>1</sup> *<sup>p</sup>* (*<sup>I</sup>* <sup>−</sup> *Qrn* )*Twnp*, and {*ρn*} ⊂ (0, +∞) satisfying

$$\liminf\_{n \to +\infty} \rho\_n \left( p - C\_q \frac{\rho\_n^{q-1}}{q} \right) > 0.$$

To prove the convergence analysis of Algorithm SASPM, we first prove some useful results.

**Lemma 5.** *Let E*<sup>1</sup> *be a p-uniformly convex and uniformly smooth real Banach space, and C*<sup>1</sup> = *E*1. *Then, for any sequence* {*yn*}, {*zn*} *and* {*wn*} *in E*1*, the set*

$$\mathcal{C}\_{n+1} = \{ \boldsymbol{\mu} \in \mathcal{C}\_n : \Delta\_p(\boldsymbol{y}\_{n\prime}\boldsymbol{\mu}) \le \Delta\_p(\boldsymbol{z}\_{n\prime}\boldsymbol{\mu}) \le \Delta\_p(\boldsymbol{w}\_{n\prime}\boldsymbol{\mu}) \}.$$

*is closed and convex for each n* ≥ 1.

**Proof.** First, since *C*<sup>1</sup> = *E*1, *C*<sup>1</sup> is closed and convex. Then, we assume that *Cn* is a closed and convex. For each *u* ∈ *Cn*, by the definition of the function Δ*p*, we have

$$\Delta\_p(y\_n, \mu) \le \Delta\_p(z\_n, \mu) \text{ if and only if } 2\langle J\_{E\_1}^p z\_n - J\_{E\_1}^p y\_n, \mu \rangle \le \frac{1}{q} (||z\_n||^p - ||y\_n||^p)\_+$$

and

$$\Delta\_{\mathcal{P}}(z\_{\boldsymbol{n}},\boldsymbol{u}) \le \Delta\_{\mathcal{P}}(w\_{\boldsymbol{n}},\boldsymbol{u}) \text{ if and only if } 2\langle \boldsymbol{J}\_{\mathcal{E}\_1}^{\boldsymbol{p}} w\_{\boldsymbol{n}} - \boldsymbol{J}\_{\mathcal{E}\_1}^{\boldsymbol{p}} z\_{\boldsymbol{n}}, \boldsymbol{u} \rangle \le \frac{1}{q} (||w\_{\boldsymbol{n}}||^{p} - ||z\_{\boldsymbol{n}}||^{p}).$$

Hence, we know that *Cn*+<sup>1</sup> is closed. In addition, we easily prove that *Cn*+<sup>1</sup> is convex. The proof is completed.

**Lemma 6.** *Let E*1, *E*2, *T T*∗ *A*, *B, and J p E*1 , *J p E*2 . *J p E*2 , *J q E*∗ 1 *be the same as above such that Conditions (1)–(4) are satisfied. If* <sup>Υ</sup> <sup>=</sup> {*<sup>z</sup>* : *<sup>z</sup>* <sup>∈</sup> *<sup>A</sup>*−10 <sup>∩</sup> *<sup>T</sup>*−1(*B*−10)}, *then* <sup>Υ</sup> <sup>⊆</sup> *Cn for any n* <sup>≥</sup> 1.

**Proof.** If Υ = ∅, it is obvious that Υ ⊆ *Cn*. Conversely, for any *z* ∈ Γ, according to Lemma 3 and using the fact that the resolvent *Rλ<sup>n</sup>* is non-expansive, we easily obtain

$$\begin{split} \Delta\_{\mathbb{P}}(y\_{\boldsymbol{n}}, \boldsymbol{z}) &= \quad \Delta\_{\mathbb{P}}(\boldsymbol{f}\_{\boldsymbol{E}\_{1}^{\boldsymbol{u}}}^{\boldsymbol{q}}(\boldsymbol{a}\_{\boldsymbol{n}} \boldsymbol{f}\_{\boldsymbol{E}\_{1}}^{\boldsymbol{p}} \boldsymbol{z}\_{\boldsymbol{n}} + (1 - \boldsymbol{a}\_{\boldsymbol{n}}) \boldsymbol{f}\_{\boldsymbol{E}\_{1}}^{\boldsymbol{p}} \boldsymbol{R}\_{\boldsymbol{\lambda}\_{\boldsymbol{n}}} \boldsymbol{z}\_{\boldsymbol{n}}), \boldsymbol{z}) \\ &\leq \quad \boldsymbol{a}\_{\boldsymbol{n}} \Delta\_{\mathbb{P}}(\boldsymbol{z}\_{\boldsymbol{n}}, \boldsymbol{z}) + (1 - \boldsymbol{a}\_{\boldsymbol{n}}) \Delta\_{\mathbb{P}}(\boldsymbol{R}\_{\boldsymbol{\lambda}\_{\boldsymbol{n}}} \boldsymbol{z}\_{\boldsymbol{n}}, \boldsymbol{z}) \\ &\leq \quad \Delta\_{\mathbb{P}}(\boldsymbol{z}\_{\boldsymbol{n}}, \boldsymbol{z}). \end{split} \tag{21}$$

From (20), let *un* = *J p E*1 (*wn*) <sup>−</sup> *<sup>ρ</sup><sup>n</sup> <sup>f</sup> <sup>p</sup>*−1(*wn*) *g*(*wn*)*<sup>p</sup> <sup>g</sup>*(*wn*) for all *<sup>n</sup>* <sup>≥</sup> 1, where *<sup>g</sup>*(*wn*) = *<sup>T</sup>*<sup>∗</sup> *<sup>J</sup> p E*1 (*Twn* − *QrnTwn*). We see from Lemma 1 that

$$\begin{split} \|u\_{n}\|\_{\dot{E}\_{1}^{\tau}}^{q} &= \quad \|\|l\_{\mathcal{E}\_{1}}^{p}(w\_{n}) - \rho\_{n}\frac{f^{p-1}(w\_{n})}{\|\mathcal{G}(w\_{n})\|^{p}}\mathcal{G}(w\_{n})\|\_{\dot{E}\_{1}^{\tau}}^{q} \\ &\leq \quad \|w\_{n}\|\|^{p} - q\rho\_{n}\frac{f^{p-1}(w\_{n})}{\|\mathcal{G}(w\_{n})\|\|^{p}}\langle w\_{n},\mathcal{G}(w\_{n})\rangle + c\_{q}\rho\_{n}^{q}\frac{f^{(p-1)q}(w\_{n})}{\|\mathcal{G}(w\_{n})\|\|^{pq}}\|\|\mathcal{G}(w\_{n})\|\|^{q} \\ &= \quad \|w\_{n}\|\|^{p} - q\rho\_{n}\frac{f^{p-1}(w\_{n})}{\|\mathcal{G}(w\_{n})\|\|^{p}}\langle w\_{n},\mathcal{G}(w\_{n})\rangle + c\_{q}\rho\_{n}^{q}\frac{f^{p}(w\_{n})}{\|\mathcal{G}(w\_{n})\|\|^{p}}. \end{split} \tag{22}$$

Then, by (16) and (22), we get

Δ*p*(*zn*, *z*) ≤ Δ*p*(*J p E*1 (*un*), *z*) <sup>=</sup> *z<sup>p</sup> p* + 1 *q J p E*1 (*un*)*<sup>p</sup>* − *z*, *<sup>u</sup>* <sup>=</sup> *z<sup>p</sup> p* + 1 *q un*(*q*−1)*<sup>p</sup>* − *z*, *un* <sup>=</sup> *z<sup>p</sup> p* + 1 *q un* (*q*−1) *<sup>q</sup>* (*q*−1) − *z*, *un* <sup>=</sup> *z<sup>p</sup> p* + 1 *q un<sup>q</sup>* − *z*, *un* <sup>=</sup> *z<sup>p</sup> p* + 1 *q un<sup>q</sup>* <sup>−</sup> *z*, *J p E*1 (*wn*) ! + *ρ<sup>n</sup> f <sup>p</sup>*−1(*wn*) *g*(*wn*)*<sup>p</sup> z*, *<sup>g</sup>*(*wn*) <sup>≤</sup> *z<sup>p</sup> p* + 1 *q wn<sup>p</sup>* <sup>−</sup> *<sup>q</sup>ρ<sup>n</sup> f <sup>p</sup>*−1(*wn*) *g*(*wn*)*<sup>p</sup> wn*, *<sup>g</sup>*(*wn*) <sup>+</sup> *cq<sup>ρ</sup> q n f <sup>p</sup>*(*wn*) *g*(*wn*)*<sup>p</sup>* − *z*, *J p E*1 (*wn*) ! + *ρ<sup>n</sup> f <sup>p</sup>*−1(*wn*) *g*(*wn*)*<sup>p</sup> z*, *<sup>g</sup>*(*wn*) <sup>=</sup> *z<sup>p</sup> p* <sup>+</sup> *wn<sup>p</sup> q* − *z*, *J p E*1 (*wn*) ! <sup>+</sup> *cq<sup>ρ</sup> q n q f <sup>p</sup>*(*wn*) *g*(*wn*)*<sup>p</sup>* <sup>+</sup> *<sup>ρ</sup><sup>n</sup> f <sup>p</sup>*−1(*wn*) *g*(*wn*)*<sup>p</sup> <sup>z</sup>* <sup>−</sup> *wn*, *<sup>g</sup>*(*wn*) <sup>=</sup> <sup>Δ</sup>*p*(*wn*, *<sup>z</sup>*) + *cq<sup>ρ</sup> q n q f <sup>p</sup>*(*wn*) *g*(*wn*)*<sup>p</sup>* <sup>+</sup> *<sup>ρ</sup><sup>n</sup> f <sup>p</sup>*−1(*wn*) *g*(*wn*)*<sup>p</sup> <sup>z</sup>* <sup>−</sup> *wn*, *<sup>g</sup>*(*wn*) (23)

On the other hand, observe that

$$\begin{aligned} \langle g(w\_n), z - w\_n \rangle &= \langle T^\* \rangle\_{\mathcal{E}\_2}^p \langle I - Q\_n T w\_n \rangle\_\ast z - w\_n \rangle \\ &= \langle \int\_{\mathcal{E}\_2}^p (I - Q\_{\mathcal{E}\_n} T w\_n), Tz - Tw\_n \rangle \\ &= \langle \int\_{\mathcal{E}\_2}^p (w\_n)(I - Q\_{\mathcal{E}\_n}) T w\_n, Q\_{\mathcal{E}\_n} T w\_n - T w\_n \rangle + \langle \int\_{\mathcal{E}\_2}^p (I - Q\_{\mathcal{E}\_n}) T w\_n, Tz - Q\_{\mathcal{E}\_n} T w\_n \rangle \\ &\le -\|(I - Q\_{\mathcal{E}\_n}) T w\_n\|^p = -p f(w\_n). \end{aligned} \tag{24}$$

By using (23) and (24), we get

$$
\Delta\_p(z\_n, z) \quad \le \quad \Delta\_p(w\_n, z) + \left(\frac{c\_q \rho\_n^q}{q} - \rho\_n p\right) \frac{f^p(w\_n)}{||\lg(w\_n)||^{p'}}\tag{25}
$$

which implies by our assumption that

$$
\Delta\_p(z\_n, z) \le \Delta\_p(w\_n, z). \tag{26}
$$

From (21) and (26), we have that *z* ∈ *Cn*+1, that is, Υ ⊆ *Cn*, for all *n* ≥ 1.

**Theorem 1.** *Let E*1, *E*<sup>2</sup> *T*, *T*∗, *A*, *B, and J p E*1 , *J p E*2 , *J q E*∗ 1 *be the same as above such that Conditions (1)–(4) are satisfied. If* <sup>Υ</sup> <sup>=</sup> {*<sup>z</sup>* : *<sup>z</sup>* <sup>∈</sup> *<sup>A</sup>*−10 <sup>∩</sup> *<sup>T</sup>*−1(*B*−10)} <sup>=</sup> <sup>∅</sup>, *then the sequence generated by Algorithm* (20) *converges strongly to a point z* = ΠΥ*x*<sup>0</sup> ∈ Υ.

**Proof.** By Lemmas 5 and 6, we know that Π*Cn*+<sup>1</sup> *x*<sup>0</sup> is well defined and Υ ⊂ *Cn*. According to Algorithm (20), we know that *xn* = Π*Cn x*<sup>0</sup> and *xn*+<sup>1</sup> = Π*Cn*+<sup>1</sup> *x*<sup>0</sup> for each *n* ≥ 1. Using Υ ⊂ *Cn* and (16), we have

$$
\Delta\_p(\mathbf{x}\_0, \mathbf{x}\_n) = \Delta\_p(\mathbf{x}\_0, \Pi\_{\mathbb{C}\_n} \mathbf{x}\_0) \le \Delta\_p(\mathbf{x}\_0, z) \; \; z \in \mathcal{Y}, \; \forall n \ge 1. \tag{27}
$$

It implies that {Δ*p*(*x*0, *xn*)} is bounded. Reusing (16), we also have

$$\begin{split} \Delta\_{\mathfrak{p}}(\mathbf{x}\_{n}, \mathbf{x}\_{n+1}) = \Delta\_{\mathfrak{p}}(\Pi\_{\mathbb{C}\_{n}} \mathbf{x}\_{0}, \mathbf{x}\_{n+1}) &\leq \quad \Delta\_{\mathfrak{p}}(\mathbf{x}\_{0}, \mathbf{x}\_{n+1}) - \Delta\_{\mathfrak{p}}(\mathbf{x}\_{0}, \Pi\_{\mathbb{C}\_{n}} \mathbf{x}\_{0}) \\ &= \quad \Delta\_{\mathfrak{p}}(\mathbf{x}\_{0}, \mathbf{x}\_{n+1}) - \Delta\_{\mathfrak{p}}(\mathbf{x}\_{0}, \mathbf{x}\_{n}). \end{split} \tag{28}$$

It follows that {Δ*p*(*x*0, *xn*+1)} is nondecreasing. Hence, the limit lim *<sup>n</sup>*→+∞Δ*p*(*x*0, *xn*) exists, and

$$\lim\_{n \to +\infty} \Delta\_p(\mathbf{x}\_{n\prime}, \mathbf{x}\_{n+1}) = 0 \tag{29}$$

It follows from (13) that

$$\lim\_{n \to +\infty} \|x\_{n+1} - x\_n\| = 0 \tag{30}$$

For some positive *m*, *n* with *m* ≥ *n*, we have *xm* = Π*Cm x*<sup>1</sup> ⊆ *Cn*. Using (16), we obtain

$$\begin{split} \Delta\_{\mathbb{P}}(\mathbf{x}\_{\mathsf{H}}, \mathbf{x}\_{\mathsf{H}}) &= \Delta\_{\mathbb{P}}(\Pi\_{\mathbb{C}\_{n}} \mathbf{x}\_{\mathsf{0}}, \mathbf{x}\_{\mathsf{m}}) \quad \leq \quad \Delta\_{\mathbb{P}}(\mathbf{x}\_{\mathsf{0}}, \mathbf{x}\_{\mathsf{m}}) - \Delta\_{\mathbb{P}}(\mathbf{x}\_{\mathsf{0}}, \Pi\_{\mathbb{C}\_{n}} \mathbf{x}\_{\mathsf{0}}) \\ &= \quad \Delta\_{\mathbb{P}}(\mathbf{x}\_{\mathsf{0}}, \mathbf{x}\_{\mathsf{m}}) - \Delta\_{\mathbb{P}}(\mathbf{x}\_{\mathsf{0}}, \mathbf{x}\_{\mathsf{m}}). \end{split} \tag{31}$$

Since the limit lim *<sup>n</sup>*→+∞Δ*p*(*x*0, *xn*) exists, it follows from (31) that lim *<sup>n</sup>*→+∞Δ*p*(*xn*, *xm*) = 0 and lim *<sup>n</sup>*→+∞*xn* <sup>−</sup> *xm* <sup>=</sup> 0. Therefore, {*xn*} is Cauchy sequence. Further, there exists a point *<sup>x</sup>*<sup>∗</sup> <sup>∈</sup> *<sup>C</sup>* such that *xn* → *x*∗.

From Algorithm (20), Definition 2, and Lemma 1, we have

<sup>Δ</sup>*p*(*wn*, *<sup>z</sup>*) = <sup>1</sup> *q J p E*∗ 1 (*J p E*1 *xn* + *θ<sup>n</sup> J p E*1 (*xn* <sup>−</sup> *xn*−1))*<sup>p</sup>* <sup>+</sup> 1 *p z<sup>p</sup>* −*J p E*1 *xn* + *θ<sup>n</sup> J p E*1 (*xn* − *xn*−1), *z* <sup>=</sup> <sup>1</sup> *q J p E*1 *xn* + *θ<sup>n</sup> J p E*1 (*xn* <sup>−</sup> *xn*−1)*<sup>q</sup>* <sup>+</sup> 1 *p z<sup>p</sup>* −*J p E*1 *xn* + *θ<sup>n</sup> J p E*1 (*xn* − *xn*−1), *z* ≤ 1 *q J p E*1 *xn<sup>q</sup>* <sup>+</sup> 1 *p z<sup>p</sup>* − *<sup>J</sup> p E*1 *xn*, *x*∗ − *θnJ p E*1 (*xn* − *xn*−1), *z* +*θnJ p E*1 (*xn* <sup>−</sup> *xn*−1), *xn* <sup>+</sup> *cq*(*θn*)*<sup>q</sup> <sup>q</sup> <sup>J</sup> p E*1 (*xn* <sup>−</sup> *xn*−1)*<sup>q</sup>* <sup>=</sup> <sup>1</sup> *q xn<sup>q</sup>* <sup>+</sup> 1 *p z<sup>p</sup>* − *<sup>J</sup> p E*1 *xn*, *x*∗ − *θnJ p E*1 (*xn* − *xn*−1), *z* +*θnJ p E*1 (*xn* <sup>−</sup> *xn*−1), *xn* <sup>+</sup> *cq*(*θn*)*<sup>q</sup> <sup>q</sup> <sup>J</sup> p E*1 (*xn* <sup>−</sup> *xn*−1)*<sup>q</sup>* = Δ*p*(*xn*, *z*) + *θnJ p E*1 (*xn* <sup>−</sup> *xn*−1), *xn* <sup>−</sup> *<sup>x</sup>*∗ <sup>+</sup> *cq*(*θn*)*<sup>q</sup> <sup>q</sup> xn* <sup>−</sup> *xn*−1*p*. (32)

By virtue of Remark 1 and the definition of *wn*, we know

$$\begin{split} \Delta\_{\mathcal{P}}(\boldsymbol{w}\_{\boldsymbol{n}}, \boldsymbol{z}) &= \quad \Delta\_{\mathcal{P}}(\boldsymbol{w}\_{\boldsymbol{n}}, \boldsymbol{x}\_{\boldsymbol{n}}) + \Delta\_{\mathcal{P}}(\boldsymbol{x}\_{\boldsymbol{n}}, \boldsymbol{z}) + \langle \mathbf{x}\_{\boldsymbol{n}} - \boldsymbol{z}, \boldsymbol{J}\_{\mathcal{E}\_{1}}^{\boldsymbol{p}} \boldsymbol{w}\_{\boldsymbol{n}} - \boldsymbol{J}\_{\mathcal{E}\_{1}}^{\boldsymbol{p}} \boldsymbol{x}\_{\boldsymbol{n}} \rangle \\ &= \quad \Delta\_{\mathcal{P}}(\boldsymbol{w}\_{\boldsymbol{n}}, \boldsymbol{x}\_{\boldsymbol{n}}) + \Delta\_{\mathcal{P}}(\boldsymbol{x}\_{\boldsymbol{n}}, \boldsymbol{z}) + \theta\_{\boldsymbol{n}} \langle \mathbf{x}\_{\boldsymbol{n}} - \boldsymbol{z}, \boldsymbol{J}\_{\mathcal{E}\_{1}}^{\boldsymbol{p}} \left(\mathbf{x}\_{\boldsymbol{n}} - \mathbf{x}\_{\boldsymbol{n}-1}\right) \rangle. \end{split} \tag{33}$$

By (32) and (33), we get <sup>Δ</sup>*p*(*wn*, *xn*) <sup>≤</sup> *cq*(*θn*)*<sup>q</sup> <sup>q</sup> xn* <sup>−</sup> *xn*−1*p*. Then, using (13) and (30) and the boundedness of the sequence {*θn*}, we can obtain

$$\lim\_{n \to +\infty} \|w\_n - \mathbf{x}\_n\| = 0.\tag{34}$$

Using a similar method, we can get

$$
\Delta\_p(w\_{n\prime}, \mathbf{x}\_{n+1}) = \Delta\_p(w\_{n\prime}, \mathbf{x}\_n) + \Delta\_p(\mathbf{x}\_n, \mathbf{x}\_{n+1}) + \langle \mathbf{x}\_n - \mathbf{x}\_{n+1}, f\_{\mathcal{E}\_1}^p w\_n - f\_{\mathcal{E}\_1}^p \mathbf{x}\_n \rangle.
$$

By setting *n* → +∞, we have

$$\lim\_{n \to +\infty} \|w\_n - x\_{n+1}\| = 0.\tag{35}$$

Since *xn*+<sup>1</sup> = Π*Cn*+<sup>1</sup> *x*<sup>0</sup> ∈ *Cn*+<sup>1</sup> ⊆ *Cn*, we have

$$
\Delta\_p(y\_{n\prime}, \mathbf{x}\_{n+1}) \le \Delta\_p(z\_{n\prime}, \mathbf{x}\_{n+1}) \le \Delta\_p(w\_{n\prime}, \mathbf{x}\_{n+1}) .
$$

According to (35), we obtain

$$\lim\_{n \to +\infty} \Delta\_{\mathcal{P}}(y\_n, \mathbf{x}\_{n+1}) = 0, \quad \lim\_{n \to +\infty} \Delta\_{\mathcal{P}}(z\_n, \mathbf{x}\_{n+1}) = 0,\tag{36}$$

which implies that lim *<sup>n</sup>*→+∞*yn* <sup>−</sup> *xn*+1 <sup>=</sup> 0, lim *<sup>n</sup>*→+∞*zn* <sup>−</sup> *xn*+1 <sup>=</sup> 0. Hence,

$$\|\|x\_n - z\_n\|\| \le \|\|x\_{n+1} - x\_n\|\| + \|\|x\_{n+1} - z\_n\|\| \to 0, \text{ as } n \to +\infty,\tag{37}$$

and

$$\|\|y\_n - z\_n\|\| \le \|\mathbf{x}\_{n+1} - y\_n\| + \|\mathbf{x}\_{n+1} - z\_n\| \to 0, \text{ as } n \to +\infty. \tag{38}$$

We also get from (34) and (37) that

$$\left|\left|w\_{n} - z\_{n}\right|\right| \le \left\|w\_{n} - \mathbf{x}\_{n}\right\| + \left\|\mathbf{x}\_{n} - z\_{n}\right\| \to 0, \text{ as } n \to \infty. \tag{39}$$

As *J p <sup>E</sup>*<sup>1</sup> is norm to norm uniformly continuous on a bounded subset of *E*1, we obtain

$$\|\|f\_{E\_1}^p(w\_n) - f\_{E\_1}^p(z\_n)\|\| \to 0, \text{ as } n \to +\infty. \tag{40}$$

Since *E*<sup>1</sup> is a *p*-uniformly convex and uniformly smooth real Banach space, then *J p <sup>E</sup>*<sup>1</sup> is uniformly norm-to-norm continuous. Thus, it follows from Algorithm (20) and real number sequence {*αn*} in [*a*, *b*] ⊂ (0, 1) that

$$\lim\_{n \to +\infty} \|f\_{E\_1}^p R\_{\lambda\_n} z\_n - f\_{E\_1}^p z\_n\| = 0 = \lim\_{n \to +\infty} \frac{1}{1 - \kappa\_n} \|f\_{E\_1}^p y\_n - f\_{E\_1}^p z\_n\| = 0\_\prime$$

which also implies that lim *<sup>n</sup>*→+∞*Rλ<sup>n</sup> zn* <sup>−</sup> *zn* <sup>=</sup> 0. From (25), and *<sup>z</sup>* being in <sup>Υ</sup>, we get

$$\begin{aligned} \Delta\_{\mathbb{P}}(z\_{n}, z) &\leq \quad \Delta\_{\mathbb{P}}(w\_{n}, z) + \rho\_{n} \left( \frac{c\_{q}\rho\_{n}^{q-1}}{q} - p \right) \frac{f^{p}(w\_{n})}{||g(w\_{n})||^{p}} \\ &= \quad \Delta\_{\mathbb{P}}(w\_{n}, z) - \rho\_{n} \left( p - \frac{c\_{q}\rho\_{n}^{q-1}}{q} \right) \frac{f^{p}(w\_{n})}{||g(w\_{n})||^{p}}. \end{aligned}$$

This implies that

$$\begin{split} \rho\_n \left( p - \frac{c\_q \rho\_n^{q-1}}{q} \right) \frac{f^p(w\_n)}{||g(w\_n)||^p} &\leq \quad \Delta\_p(w\_{n'}z) - \Delta\_p(z\_{n'}z) \\ &= \quad \frac{1}{q}||w\_n||^p - \frac{1}{q}||z\_n||^p - \langle f\_{E\_1}^p w\_n - f\_{E\_1}^p z\_{n'}z \rangle \\ &= \quad \Delta\_p(w\_{n'}z\_n) + \langle f\_{E\_1}^p w\_n - f\_{E\_1}^p z\_{n'}z\_n - z \rangle \\ &\leq \quad (||w\_n - z\_n|| + ||z\_n - z||) ||f\_{E\_1}^p w\_n - f\_{E\_1}^p z\_n||. \end{split}$$

By setting of *n* → +∞, the right-hand side of the last inequality tends to 0. This implies that

$$\rho\_n \left( p - \frac{c\_q \rho\_n^{q-1}}{q} \right) \frac{f^p(w\_n)}{\|g(w\_n)\|^p} \to 0, \ n \to +\infty. \tag{41}$$

Since lim inf *<sup>n</sup>*→+<sup>∞</sup> *<sup>ρ</sup><sup>n</sup> p* − *cq ρ q*−1 *n q* > 0, we get

$$\frac{f^p(w\_n)}{||g(w\_n)||^p} \to 0, \ n \to +\infty$$

and hence

$$\frac{f(w\_n)}{||g(w\_n)||^p} \to 0, \ n \to +\infty \tag{42}$$

Furthermore, since {*g*(*wn*)} is bounded, we obtain from (42) that

$$\begin{aligned} 0 \le \mathcal{g}(w\_n) &= \quad ||\mathcal{g}(w\_n)|| \frac{f(w\_n)}{||\mathcal{g}(w\_n)||} \\ &\le \quad M\_1 \frac{f(w\_n)}{||\mathcal{g}(w\_n)||} \to 0, \; n \to +\infty \end{aligned}$$

for some *M*<sup>1</sup> > 0. Therefore,

$$\lim\_{n \to +\infty} f(w\_n) = 0.$$

Hence,

$$\lim\_{n \to +\infty} \left\| \left( I - Q\_{r\_n} \right) T w\_n \right\| = 0.$$

In addition,

$$\|\|T^\*f\_{E\_2}^p(I-Q\_{r\_n})Tw\_n\|\| \le \|\|T\|\|\|(I-Q\_{r\_n})Tw\_n\|\| \to 0, \ n \to +\infty.$$

Since *xn* − *wn* → 0, as *n* → +∞, there exists a subsequence {*xnj* } of {*xn*} such that *xnj w* ∈ *E*1, as well as *xn* − *wn* → 0, as *n* → +∞ there exists a subsequence {*wnj* } of {*wn*} such that *wnj w* ∈ *E*1. From *Twn* − *QrnTwn* → 0 and by the boundedness and linearity of *T*, we have *Twnj Tw* and *Qrnj Twnj Tw*. Since *Qrn* is a metric resolvent on *B* for *rn* > 0, we have

$$\frac{J\_{E\_2}^p(Tw\_n - Q\_{r\_n}Tw\_n)}{r\_n} \in BQ\_{r\_n}Tw\_n$$

for all *<sup>n</sup>* <sup>∈</sup> <sup>N</sup>, thus we obtain

$$0 \le \left\langle v - Q\_{r\_{n\_j}} T w\_{n\_j} T w\_{n\_j}, v^\* - \frac{J\_{\mathbb{E}\_2}^p (T w\_{n\_j} - Q\_{r\_{n\_j}} T w\_{n\_j})}{r\_{n\_j}} \right\rangle\_0$$

for all (*v*, *v*∗) ∈ *B*. It follows that

$$0 \le \langle v - Tw, v^\* - 0 \rangle$$

for all (*v*, *<sup>v</sup>*∗) <sup>∈</sup> *<sup>B</sup>*. Since *<sup>B</sup>* is maximal monotone, *Tw* <sup>∈</sup> *<sup>B</sup>*−10 and hence *<sup>w</sup>* <sup>∈</sup> *<sup>T</sup>*−1(*B*−10). Let *bn* <sup>=</sup> *<sup>R</sup>λ<sup>n</sup> zn* and *kn* <sup>=</sup> *Twn* <sup>−</sup> *QrnTwn* <sup>∀</sup>*<sup>n</sup>* <sup>∈</sup> <sup>N</sup>

*bn* = *Jλ<sup>n</sup> J q E*∗ 1 *J p E*1 (*wn*) − *λnT*<sup>∗</sup> *J p E*2 (*kn*) ⇐⇒ *bn* = *J p <sup>E</sup>*<sup>1</sup> + *<sup>λ</sup>nA* −<sup>1</sup> *J p E*1 *J q E*∗ 1 (*J p E*1 (*wn*) − *λnT*<sup>∗</sup> *J p E*2 (*kn*)) ⇐⇒ *bn* = *J p <sup>E</sup>*<sup>1</sup> + *<sup>λ</sup>nA* −<sup>1</sup> *J p E*1 (*wn*) − *λnT*<sup>∗</sup> *J p E*2 (*kn*) ⇐⇒ *J p E*1 (*wn*) − *λnT*<sup>∗</sup> *J p E*2 (*kn*) ∈ *J p E*1 (*bn*) + *λnAbn* ⇐⇒ *J p E*1 (*wn*) − *J p E*1 (*bn*) *λn* − *T*<sup>∗</sup> *J p E*2 (*kn*) ∈ *Abn*.

Note that

$$\begin{aligned} \|\|f\_{\mathcal{E}\_1}^p(\boldsymbol{w}\_n) - f\_{\mathcal{E}\_1}^p(b\_n)\|\| &= \|\|f\_{\mathcal{E}\_1}^p(\boldsymbol{w}\_n) - f\_{\mathcal{E}\_1}^p(\mathcal{R}\_{\lambda\_n}\boldsymbol{z}\_n)\|\| \\ &\leq \|\|f\_{\mathcal{E}\_1}^p(\boldsymbol{w}\_n) - f\_{\mathcal{E}\_1}^p(\boldsymbol{z}\_n)\|\| + \|\|f\_{\mathcal{E}\_1}^p(\boldsymbol{z}\_n) - f\_{\mathcal{E}\_1}^p(\mathcal{R}\_{\lambda\_n}\boldsymbol{z}\_n)\|\| \to 0, \; n \to +\infty. \end{aligned} \tag{43}$$

By the monotonicity of *A*, it follows that

$$0 \le \left\langle v - b\_{\hbar \prime} v^\* - \frac{J\_{E\_1}^p(w\_n) - J\_{E\_1}^p(b\_{\hbar})}{\lambda\_n} + T^\* J\_{E\_2}^p(k\_{\hbar}) \right\rangle\_{\Omega}$$

for all (*v*, *v*∗) ∈ *A*. Then,

$$0 \le \left\langle v - b\_{n\_{\nu}} v^\* - \frac{J\_{E\_1}^p(w\_{n\_i}) - J\_{E\_1}^p(b\_{n\_i})}{\lambda\_{n\_i}} + T^\* J\_{E\_2}^p(k\_{n\_i}) \right\rangle.$$

Since *bni <sup>w</sup>*, (40) and (43), it follows that 0 ≤ *<sup>v</sup>* <sup>−</sup> *<sup>w</sup>*, *<sup>v</sup>*<sup>∗</sup> <sup>−</sup> <sup>0</sup> and hence *<sup>w</sup>* <sup>∈</sup> *<sup>A</sup>*−10. This concludes that *<sup>w</sup>* <sup>∈</sup> *<sup>A</sup>*−10 <sup>∩</sup> *<sup>T</sup>*−1(*B*−10). Then, from (28) and (20), we have

$$\langle f\_{\mathcal{E}\_1}^p \mathbf{x}\_0 - f\_{\mathcal{E}\_1}^p \mathbf{x}\_{\mathcal{U}}, p - \mathbf{x}\_{\mathcal{U}} \rangle, \text{ for all } p \in \mathcal{Y}. \tag{44}$$

By setting *n* → +∞ in (44), we obtain

$$\langle \langle f\_{E\_1}^p x\_0 - f\_{E\_1}^p x^\*, p - x^\* \rangle \le 0, \text{ for all } p \in \mathcal{Y}. \tag{45}$$

Again, from (15), we have *x*<sup>∗</sup> = ΠΥ*x*0. Definitely, we obtain that {*xn*} generated by Algorithm (20) strongly converges to *x*<sup>∗</sup> = ΠΥ*x*<sup>0</sup> ∈ Υ. The proof is completed.

As a corollary of Theorem 1, when *E*<sup>1</sup> and *E*<sup>2</sup> reduces to Hilbert spaces, the function Δ*<sup>p</sup>* is equal to <sup>1</sup> <sup>2</sup> *<sup>x</sup>* <sup>−</sup> *<sup>y</sup>*<sup>2</sup> and the Bregman projection <sup>Π</sup>*<sup>C</sup>* is equivalent to the metric projection *PC*. Then, we obtain the following result.

**Theorem 2.** *Let <sup>H</sup>*<sup>1</sup> *and <sup>H</sup>*<sup>2</sup> *be Hilbert spaces, <sup>A</sup>* : *<sup>H</sup>*<sup>1</sup> <sup>→</sup> <sup>2</sup>*H*<sup>1</sup> *and <sup>B</sup>* : *<sup>H</sup>*<sup>2</sup> <sup>→</sup> <sup>2</sup>*H*<sup>2</sup> *be maximal monotone operators, T* : *H*<sup>1</sup> → *H*<sup>2</sup> *be a bounded linear operator with T* = 0*, and T*<sup>∗</sup> : *H*<sup>2</sup> → *H*<sup>1</sup> *be the adjoint of T*. *Let R<sup>λ</sup> be the resolvent operator associated with a maximal monotone operator A on H*<sup>1</sup> *and Qr be metric resolvent associated with a maximal monotone operator <sup>B</sup> on <sup>H</sup>*2. *Suppose that* <sup>Υ</sup> <sup>=</sup> *<sup>A</sup>*−10 <sup>∩</sup> *<sup>T</sup>*−1(*B*−10) <sup>=</sup> <sup>∅</sup>. *For fixed <sup>x</sup>*<sup>0</sup> <sup>∈</sup> *<sup>H</sup>*1, *let* {*xn*}+<sup>∞</sup> *<sup>n</sup>*=<sup>0</sup> *be iteratively generated by x*<sup>1</sup> ∈ *H*<sup>1</sup> *and*

$$\begin{cases} w\_n = \mathbf{x}\_n + \theta\_n (\mathbf{x}\_n - \mathbf{x}\_{n-1}), \\ z\_n = w\_n - \rho\_n \frac{f(w\_n)}{\|T^\*(I - Q\_{r\_n})Tw\_n\|^2} [T^\*(I - Q\_{r\_n})Tw\_n] \\ y\_n = a\_n z\_n + (1 - a\_n) R\_{\lambda\_n} z\_n \\ \mathbb{C}\_{n+1} = \{\boldsymbol{\mu} \in \mathbb{C}\_n : ||y\_n - \boldsymbol{\mu}|| \le ||z\_n - \boldsymbol{\mu}|| \le ||w\_n - \boldsymbol{\mu}||\}, \\ x\_{n+1} = P\_{\mathbb{C}\_{n+1}} x\_{0\prime} \end{cases} \tag{46}$$

*where PCn*<sup>+</sup><sup>1</sup> *is the metric projection of H*<sup>1</sup> *onto Cn*+1, *the sequence of real numbers,* {*αn*} ⊂ [*a*, *b*] ⊂ (0, 1) *and* {*θn*} ⊂ [*c*, *<sup>d</sup>*] <sup>⊂</sup> (−∞, <sup>+</sup>∞). *<sup>f</sup>*(*wn*) :<sup>=</sup> <sup>1</sup> <sup>2</sup> (*<sup>I</sup>* <sup>−</sup> *Qrn* )*Twn*2*, and* {*ρn*} ∈ (0, 4). *Then, the sequence* {*xn*} *generated by* (46) *converges strongly to a point z*<sup>0</sup> = *P*Υ*x*<sup>0</sup> ∈ Υ.

### **4. Applications**

In this section, we provide some applications of our result to solving other nonlinear optimization problems.

### *4.1. Application to Minimization Problem*

First, we consider an application of our result to convex minimization problem in real Banach space *E*. Let *ϑ* : *E* → (−∞, +∞] be a proper, convex and lower semicontinuous function. The convex minimization problem is to find *x* ∈ *E* such that

$$
\vartheta(x) \le \vartheta(y), \quad \text{for all } y \in E.
$$

The set of minimizer of *ϑ* is denoted by *Argmin ϑ*. The subdifferential of *∂ϑ* of *ϑ* is defined as follows

$$\partial \theta(\mathfrak{u}) = \{ w \in E^\* : \theta(\mathfrak{u}) + \langle v - \mathfrak{u}, w \rangle \le \theta(\mathfrak{u}), \text{ for all } v \in E \},$$

for all *u* ∈ *E*. From Rockafellar [37], it is known that *∂ϑ* is a maximal monotone operator. Let *C* be a nonempty, closed, and convex subset of *E* and let *iC* be the indicator function of *C* i.e.,

$$i\_{\mathbb{C}}(u) = \begin{cases} 0, & u \in \mathbb{C} \\ \infty, & u \notin \mathbb{C}. \end{cases}$$

Then, *iC* is a proper, convex, and lower semicontinuous function on *E*. Thus, the subdifferential *∂iC* of *iC* is a maximal monotone operator. Then, we can define the resolvent *R<sup>λ</sup>* of *∂iC* for *λ* > 0 i.e.,

$$\mathcal{R}\_{\lambda}u = (I\_p + \lambda \partial\_{i\_C})^{-1} I\_{p}u$$

for all *u* ∈ *E* and *p* ∈ (1, +∞). We have that for any *x* ∈ *E* and *u* ∈ *C*

$$\begin{aligned} \boldsymbol{u} &= \mathcal{R}\_{\lambda}\mathbf{x} \quad \text{if and only if} \quad \boldsymbol{f}\_{\mathcal{P}}\mathbf{x} \in \boldsymbol{f}\_{\mathcal{P}}\boldsymbol{u} + \lambda \partial\_{\bar{\mathbf{c}}}\boldsymbol{u} \\ &\quad \text{if and only if} \quad \frac{1}{\lambda}(\boldsymbol{f}\_{\mathcal{P}}\mathbf{x} - \boldsymbol{f}\_{\mathcal{P}}\boldsymbol{u}) \in \partial\_{\bar{\mathbf{c}}}\boldsymbol{u} \\ &\quad \text{if and only if} \quad \boldsymbol{i}\_{\mathcal{C}}\boldsymbol{y} \ge \langle \boldsymbol{y} - \boldsymbol{u}, \frac{1}{\lambda}(\boldsymbol{f}\_{\mathcal{P}}\mathbf{x} - \boldsymbol{f}\_{\mathcal{P}}\boldsymbol{u}) \rangle + \boldsymbol{i}\_{\mathcal{C}}\boldsymbol{u} \text{ for all } \boldsymbol{y} \in \mathbb{C} \\ &\quad \text{if and only if} \quad \boldsymbol{0} \ge \langle \boldsymbol{y} - \boldsymbol{u}, \frac{1}{\lambda}(\boldsymbol{f}\_{\mathcal{P}}\mathbf{x} - \boldsymbol{f} - \boldsymbol{p}\boldsymbol{u}) \rangle, \text{ for all } \boldsymbol{y} \in \mathbb{C} \\ &\quad \text{if and only if} \quad \langle \boldsymbol{y} - \boldsymbol{u}, \boldsymbol{f}\_{\mathcal{P}}\mathbf{x} - \boldsymbol{f}\_{\mathcal{P}}\boldsymbol{u} \rangle \le 0, \text{ for all } \mathbf{x} \in \mathbb{C} \\ &\quad \text{if and only if} \quad \boldsymbol{u} = \boldsymbol{\Pi}\_{\mathcal{C}}\mathbf{x}. \end{aligned}$$

Let *E*<sup>1</sup> and *E*<sup>2</sup> be real Banach spaces and *ϑ* : *E*<sup>1</sup> → (−∞, +∞] and *ξ* : *E*<sup>2</sup> → (−∞, +∞] be proper, lower semicontinuous, and convex functions such that *Argminϑ* = ∅ and *Argminξ* = ∅. Consider the Split Proximal Feasibility Problem (SPFP) defined by: Find *x* ∈ *E*<sup>1</sup> such that

$$\mathbf{x} \in \operatorname{Argmin} \theta \quad \text{and} \quad \mathbf{A} \mathbf{x} \in \operatorname{Argmin} \mathfrak{J}\_{\mathbf{\uplor}} \tag{47}$$

where *Argmin ϑ* := {*x*¯ ∈ *E*<sup>1</sup> : *ϑ*(*x*¯) ≤ *ϑ*(*x*), for all *x* ∈ *E*1}, and *Argmin ξ* = {*y*¯ ∈ *E*<sup>2</sup> : *ξ*(*y*¯) ≤ *ξ*(*y*), for all *y* ∈ *E*2}. We denote the solution set of (47) by Ω. The PSFP is a generalization of the split feasibility problem and has been studied extensively by many authors in real Hilbert space (see, e.g., [38–42]).

By setting *A* = *∂ϑ* and *B* = *∂ξ*, we obtain a strong convergence result for solving (47) in real Banach spaces.

**Theorem 3.** *Let E*<sup>1</sup> *be a p-uniformly convex and uniformly smooth Banach space and E*<sup>2</sup> *be a uniformly convex smooth Banach space. Let ϑ and ξ be proper, lower semicontinuous, and convex functions of E*<sup>1</sup> *into* (−∞, +∞] *and <sup>E</sup>*<sup>2</sup> *into* (−∞, <sup>+</sup>∞] *such that* (*∂ϑ*)−10 <sup>=</sup> <sup>∅</sup> *and* (*∂ξ*)−10 <sup>=</sup> <sup>∅</sup>, *respectively. Let <sup>T</sup>* : *<sup>E</sup>*<sup>1</sup> <sup>→</sup> *<sup>E</sup>*<sup>2</sup> *be a bounded linear operator such that T* = 0 *and let T*<sup>∗</sup> *be the adjoint operator T*. *Suppose that* Ω = ∅. *For fixed x*<sup>0</sup> ∈ *E*1, *let* {*xn*}<sup>∞</sup> *<sup>n</sup>*=<sup>0</sup> *be iteratively generated by x*<sup>1</sup> ∈ *E*<sup>1</sup> *and*

$$\begin{cases} w\_{n} = \int\_{E\_{1}^{n}} \left[ \int\_{E\_{1}}^{p} x\_{n} + \theta\_{n} \right]\_{E\_{1}}^{p} \left( x\_{n} - x\_{n-1} \right) \right], \\ v\_{n} = \operatorname\*{arg\,min}\_{y \in E\_{2}} \left\{ \xi(y) + \frac{1}{\mu\_{n}} ||y||^{2} - \frac{1}{\mu\_{n}} \langle y, \int\_{E\_{2}}^{p} Tx\_{n} \rangle \right\} \\ z\_{n} = \int\_{E\_{1}^{n}}^{q} \left[ \int\_{E\_{1}}^{p} (w\_{n}) - \mu\_{n} \frac{f^{p-1}(w\_{n})}{\|\mathcal{T}^{\*} (\mu\_{n}^{p} - v\_{n})\|^{p}} T^{\*} f\_{E\_{1}}^{p} (Tw\_{n} - v\_{n}) \right], \\ u\_{n} = \operatorname\*{arg\,min}\_{y \in E\_{1}} \{ \xi(\mathbf{x}) + \frac{1}{\sigma\_{n}} ||\mathbf{x}||^{2} - \frac{1}{\sigma\_{n}} \langle \mathbf{x}, f\_{E\_{2}}^{p} z\_{n} \rangle \} \\ y\_{n} = f\_{E\_{1}}^{q} \left( \mathbf{a}\_{n} f\_{E\_{1}}^{p} z\_{n} + (1 - \mathbf{a}\_{n}) f\_{E\_{1}}^{p} u\_{n} \right), \\ \mathcal{C}\_{n+1} = \{ u \in \mathcal{C}\_{n} : \Delta\_{p}(y\_{n}, u) \le \Delta\_{p}(z\_{n}, u) \le \Delta\_{p}(w\_{n}, u) \}, \\ x\_{n+1} = \Pi\_{\mathcal{C}\_{n+1}} \mathbf{x}\_{0} \end{cases} \tag{48}$$

*where* {*σn*}, {*μn*} ⊂ (0, +∞), Π*Cn*+<sup>1</sup> *is a Bregman projection of E*<sup>1</sup> *onto Cn*+1, *the sequence of real number* {*αn*} ⊂ [*a*, *<sup>b</sup>*] <sup>⊂</sup> (0, 1) *and* {*θn*} ⊂ [*c*, *<sup>d</sup>*] <sup>⊂</sup> (−∞, <sup>+</sup>∞), *<sup>f</sup>*(*wn*) :<sup>=</sup> <sup>1</sup> *<sup>p</sup> Twn* <sup>−</sup> *vnp, and* {*ρn*} ⊂ (0, <sup>+</sup>∞) *satisfies*

$$\liminf\_{n \to +\infty} \rho\_n \left( p - \mathbb{C}\_q \frac{\rho\_n^{q-1}}{q} \right) > 0.$$

*where cq is the uniform smoothness coefficient of E*1*. Then, xn* <sup>→</sup> *<sup>z</sup>*<sup>0</sup> <sup>∈</sup> (*∂ϑ*)−10 <sup>∩</sup> *<sup>T</sup>*−1((*∂ξ*)−10), *where z*<sup>0</sup> :<sup>=</sup> <sup>Π</sup>(*∂ϑ*)−10∩*T*−1((*∂ξ*)−10)*x*<sup>0</sup>

*Axioms* **2020**, *9*, 140

**Proof.** We know from [43] that

$$\upsilon\_{\mathfrak{n}} = \arg\min\_{\mathcal{Y}\in E\_2} \{ \mathcal{\xi}(\mathcal{y}) + \frac{1}{2\mu\_n} \|\mathcal{y}\|^2 - \frac{1}{\mu\_n} \} \langle \mathcal{y}, J\_{E\_2}^p Tw\_{\mathfrak{n}} \rangle\_{\mathfrak{n}}$$

is equivalent to

$$0 \in (\partial \xi^\mathbf{r})\mathfrak{x}\_\mathbf{n} + \frac{1}{\mu\_\mathbf{n}} J\_{\to\_2}^p \mathfrak{x}\_\mathbf{n} - \frac{1}{\mu\_\mathbf{n}} J\_{\to\_2}^p Tw\_\mathbf{n}$$

From this, we have *J p E*2 *Twn* ∈ *J p E*2 *vn* + *μn*(*∂ξ*)*vn* i.e., *vn* = *QrnTwn*. Similarly, we have that

$$u\_{\mathfrak{n}} = \arg\min\_{\mathbf{x}\in E\_1} \{ \theta(\mathbf{x}) + \frac{1}{2\sigma\_{\mathfrak{n}}} \left\| \mathbf{x} \right\|^2 - \frac{1}{\sigma\_{\mathfrak{n}}} \langle \mathbf{x}, J\_{E\_1}^p z\_{\mathfrak{n}} \rangle \}.$$

is equivalent to *un* = *Rλ<sup>n</sup> zn*. Using Theorem 1, we get the conclusion.

### *4.2. Application to Equilibrium Problem*

Let *<sup>C</sup>* be a nonempty closed and convex subset of a Banach space *<sup>E</sup>* and let *<sup>G</sup>* : *<sup>C</sup>* <sup>×</sup> *<sup>C</sup>* <sup>→</sup> <sup>R</sup> be a bifunction. For solving the equilibrium problem, we assume that *G* satisfies the following conditions:


$$\limsup\_{t \to 0^+} G(tz + (1 - t)x, y) \le G(x, y).$$

(A4) *G*(*x*, 0) is convex and lower semicontinuous for each *x* ∈ *C*.

The equilibrium problem is to find *x*<sup>∗</sup> ∈ *C* such that

$$G(\mathfrak{x}^\*, y) \ge 0 \quad \text{for all } y \in \mathbb{C}.$$

The set of solution of this problem is denoted by *EP*(*G*).

**Lemma 7.** *[44] Let g* : *E* → (−∞, +∞] *be super coercive Legendre function, G be a bifunction of C* × *C into* <sup>R</sup> *satisfying Conditions (A1)–(A4), and x* <sup>∈</sup> *<sup>E</sup>*. *Define a mapping S<sup>g</sup> <sup>G</sup>* : *E* → *C as follows:*

$$\mathcal{S}\_{\mathbb{G}}^{\mathbb{g}}(\mathbf{x}) = \{ z \in \mathbb{C} : \mathbf{G}(z, y) + \langle y - z, \nabla \mathbf{g}(z) - \nabla \mathbf{g}(\mathbf{x}) \rangle \ge 0 \text{ for all } \mathbf{y} \in \mathbb{C} \}.$$

*Then,*


$$D\_{\mathcal{S}}(\boldsymbol{\mu}, \boldsymbol{S}\_{G}^{\mathcal{S}}(\boldsymbol{x})) + D\_{\mathcal{S}}(\boldsymbol{S}\_{G}^{\mathcal{S}}(\boldsymbol{x}), \boldsymbol{x}) \leq D\_{\mathcal{S}}(\boldsymbol{\mu}, \boldsymbol{x}).$$

**Proposition 3.** *[45] Let g* : *E* → (−∞, +∞] *be a super coercive Legendre Frécht differentiable and totally convex function. Let <sup>C</sup> be a closed and convex subset of <sup>E</sup> and assume that the bifunction <sup>G</sup>* : *<sup>C</sup>* <sup>×</sup> *<sup>C</sup>* <sup>→</sup> <sup>R</sup> *satisfies the Conditions (A1)–(A4). Let AG be a set-valued mapping of E into* 2*E*<sup>∗</sup> *defined by*

$$A\_G(\mathbf{x}) = \begin{cases} \{ \mathbf{z} \in E^\* : \mathbf{G}(\mathbf{x}, \mathbf{y}) \ge \langle \mathbf{y} - \mathbf{x}, \mathbf{z} \rangle \text{ for all } \mathbf{y} \in \mathbb{C} \}, & \mathbf{x} \in \mathbb{C} \\ \mathcal{O}, & \mathbf{x} \in E - \mathbb{C}. \end{cases}$$

*Then, AG is a maximal monotone operator, EP*(*G*) = *A*−<sup>1</sup> *<sup>G</sup>* (0) *and S<sup>g</sup> <sup>G</sup>* <sup>=</sup> *<sup>R</sup><sup>g</sup> AG* .

Let *E*<sup>1</sup> and *E*<sup>2</sup> real Banach spaces and *C* and *Q* be nonempty, closed, and convex subsets of *E*<sup>1</sup> and *<sup>E</sup>*2, respectively. Let *<sup>G</sup>*<sup>1</sup> : *<sup>C</sup>* <sup>×</sup> *<sup>C</sup>* <sup>→</sup> <sup>R</sup> and *<sup>G</sup>*<sup>2</sup> : *<sup>Q</sup>* <sup>×</sup> *<sup>Q</sup>* <sup>→</sup> <sup>R</sup> be bifunctions satisfying Conditions (A1)–(A4) and *T* : *E*<sup>1</sup> → *E*<sup>2</sup> be a bounded linear operator. We consider the Split Equilibrium Problem (SEP) defined by: Find *x* ∈ *C* such that

$$\mathbf{x} \in EP(\mathbf{G\_1}) \quad \text{and} \quad T\mathbf{x} \in EP(\mathbf{G\_2}). \tag{49}$$

The SEP was introduced by Moudafi [46] and has been studied by many authors for Hilbert and Banach spaces (see, e.g., [47–50]). We denote the set of solution of (49) by *SEP*(*G*1, *G*2).

Setting *A* = *AG*<sup>1</sup> and *B* = *AG*<sup>2</sup> in Algorithm (20), Lemma 7, and Proposition 3, we obtain a strong convergence result for solving SEP in real Banach spaces.

**Theorem 4.** *Let E*<sup>1</sup> *be a p-uniformly convex and uniformly smooth Banach space, E*<sup>2</sup> *be a uniformly smooth Banach space, and <sup>C</sup> and <sup>Q</sup> be nonempty closed subsets of <sup>E</sup>*<sup>1</sup> *and <sup>E</sup>*2*, respectively. Let <sup>G</sup>* : *<sup>C</sup>* <sup>×</sup> *<sup>C</sup>* <sup>→</sup> <sup>R</sup> *and <sup>H</sup>* : *<sup>Q</sup>* <sup>×</sup> *<sup>Q</sup>* <sup>→</sup> <sup>R</sup> *be bifunctions satisfying Conditions (A1)–(A4) and <sup>g</sup>* : *<sup>E</sup>*<sup>1</sup> <sup>→</sup> <sup>R</sup> *and <sup>h</sup>* : *<sup>E</sup>*<sup>2</sup> <sup>→</sup> <sup>R</sup> *be super coercive Legendre functions which are bounded, uniformly Frechet differentiable, and totally convex on bounded subset of E*2. *Let T* : *E*<sup>1</sup> → *E*<sup>2</sup> *be a bounded linear operator with T* = 0 *and T*<sup>∗</sup> : *E*<sup>∗</sup> <sup>2</sup> → *E*<sup>∗</sup> <sup>1</sup> *be the adjoint of T*. *Suppose that SEP*(*G*1, *<sup>G</sup>*2) <sup>=</sup> <sup>∅</sup> *for fixed x*<sup>0</sup> <sup>∈</sup> *<sup>E</sup>*1*, let* {*xn*}<sup>∞</sup> *<sup>n</sup>*=<sup>0</sup> *be iteratively generated by x*<sup>1</sup> ∈ *E*1*, and*

⎧ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨ ⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩ *wn* = *J q E*∗ 1 *J p E*1 *xn* + *θ<sup>n</sup> J p E*1 (*xn* − *xn*−1) , *zn* = *J q E*∗ 1 *J p E*1 (*wn*) <sup>−</sup> *<sup>ρ</sup><sup>n</sup> <sup>f</sup> <sup>p</sup>*−1(*wn*) *T*∗(*J p E*2 (*Twn*−*S<sup>h</sup> Hn Twn*)*p*) *T*∗ *J p E*1 (*Twn* <sup>−</sup> *<sup>S</sup><sup>h</sup> HnTwn*) , *yn* = *J q E*1 *α<sup>n</sup> J p E*1 *zn* + (1 − *αn*)*J p E*1 *Sg Gn zn* , *Cn*+<sup>1</sup> = {*u* ∈ *Cn* : Δ*p*(*yn*, *u*) ≤ Δ*p*(*zn*, *u*) ≤ Δ*p*(*wn*, *u*)}, *xn*+<sup>1</sup> = Π*Cn*+<sup>1</sup> *x*<sup>0</sup> (50)

where {*Hn*} and {*Gn*} ⊂ (0, <sup>+</sup>∞), *<sup>f</sup>*(*wn*) = <sup>1</sup> *<sup>p</sup>* (*<sup>I</sup>* <sup>−</sup> *<sup>S</sup><sup>h</sup> Hn* )*Tunp*, <sup>Π</sup>*Cn*+<sup>1</sup> is a Bregman projection of *<sup>E</sup>*<sup>1</sup> onto *Cn*+1, the sequence of real number {*αn*} ⊂ [*a*, *b*] ⊂ (0, 1) and {*θn*} ⊂ [*c*, *d*] ⊂ (−∞, +∞), and {*ρn*} ⊂ (0, +∞) satisfies

$$\lim\_{n \to +\infty} \inf\_{\theta} \rho\_n \left( p - \mathbb{C}\_q \frac{\rho\_n^{q-1}}{q} \right) > 0.$$

where *cq* is the uniform smoothness coefficient of *E*1. Then, *xn* → *z*<sup>0</sup> ∈ Π*SEP*(*G*1,*G*2)*x*0.

### **5. Conclusions**

In this paper, we introduce a new inertial shrinking projection method for solving the split common null point problem in uniformly convex and uniformly smooth real Banach spaces. The algorithm is designed such that its step size does not require prior knowledge of the norm of the bounded linear operator. A strong convergence result is also proved under some mild

conditions. We further provide some applications of our result to other nonlinear optimization problems. We highlight our contributions in this paper as follow:


**Author Contributions:** Conceptualization, C.C.O.; methodology, C.C.O. and L.O.J.; software, C.C.O. and L.O.J.; validation, C.C.O., L.O.J. and R.N.; formal analysis, C.C.O.and L.O.J.; writing–original draft preparation, C.C.O. and L.O.J.; writing–review and editing, C.C.O., L.O.J. and R.N.; supervision, L.O.J.; project administration, C.C.O.; funding acquisition, L.O.J. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research is funded by the Mathematical research fund at the Sefako Makgatho Health Sciences University.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Nonlinear Approximations to Critical and Relaxation Processes**

### **Simon Gluzman**

Materialica+ Research Group, Bathurst St. 3000, Apt. 606, Toronto, ON M6B 3B4, Canada; gluz@sympatico.ca

Received: 5 September 2020; Accepted: 22 October 2020; Published: 28 October 2020

**Abstract:** We develop nonlinear approximations to critical and relaxation phenomena, complemented by the optimization procedures. In the first part, we discuss general methods for calculation of critical indices and amplitudes from the perturbative expansions. Several important examples of the Stokes flow through 2D channels are brought up. Power series for the permeability derived for small values of amplitude are employed for calculation of various critical exponents in the regime of large amplitudes. Special nonlinear approximations valid for arbitrary values of the wave amplitude are derived from the expansions. In the second part, the technique developed for critical phenomena is applied to relaxation phenomena. The concept of time-translation invariance is discussed, and its spontaneous violation and restoration considered. Emerging probabilistic patterns correspond to a local breakdown of time-translation invariance. Their evolution leads to the time-translation invariance complete (or partial) restoration. We estimate the typical time extent, amplitude and direction for such a restorative process. The new technique is based on explicit introduction of origin in time as an optimization parameter. After some transformations, we arrive at the exponential and generalized exponential-type solutions (Gompertz approximants), with explicit finite time scale, which is only implicit in the initial parameterization with polynomial approximation. The concept of crash as a fast relaxation phenomenon, consisting of time-translation invariance breaking and restoration, is advanced. Several COVID-related crashes in the time series for Shanghai Composite and Dow Jones Industrial are discussed as an illustration.

**Keywords:** critical index; relaxation time; time-translation invariance breaking and restoration; market crash; COVID-19; Gompertz approximants

### **1. Introduction**

Let the function Φ(*x*) of a real variable *x* ∈ [0, ∞) be defined by some rather complicated problem. The variable *x* > 0 can represent, e.g., a coupling constant or concentration of particles. Of course, one should strive to find an exact solution to the problem [1,2]. Among such exact solutions one can find the solution to the celebrated Kondo problem and its thermodynamics. In a number of cases important for optical applications, such as Bessel beams and its generalizations [3], one can find an intriguing physics already within the linear wave equation. In optics, there are a variety of exact solutions: spatial, temporal, dark optical solitons and breathers all follow from the celebrated nonlinear Schrödinger equation and its modifications [4]. The so-called spatiotemporal *X*-waves, another type of the closed-form solutions are being studied as well (see, e.g., [5]).

What if such a problem does not allow for an explicit solution for the sought function? Let us assume that some kind of perturbation theory is still possible to develop, so that it generates formal power series about the point *x* = *x*<sup>0</sup> = 0, Φ(*x*) = ∑<sup>∞</sup> *<sup>n</sup>*=<sup>0</sup> *cnxn*, for the function in question [6]. The perturbation methods can generate the series (often slowly) convergent for all *x* smaller than the radius of convergence, or the series divergent for all *x*, except *x* = 0.

That is, for smooth function Φ(*x*) [7], we have the asymptotic power series [7,8],

$$\Phi(\mathbf{x}) \sim \sum\_{n=0}^{\infty} c\_n \mathbf{x}^n. \tag{1}$$

Our task is to recast the series (2) into some convergent expressions by means of a nonlinear analytical constructs, the so-called approximants. When literally all of the terms in divergent series are known, one can invoke Euler of Borel summation [8]. Even for convergent series there is still a problem of how to continue the expansion outside of radius of convergence [9], where the approximants could be useful.

However, in realistic problems, only a few terms on the RHS of (1) can be calculated, and applying various approximants is the only available analytical option for the truncated series (5) and (A6). The approximants are conditioned to be asymptotically equivalent to the series (1), truncated at some finite number *k*. However, the approximants are able to generate an additional infinite number of coefficients, approximating unknown exact coefficients. Determination of the best approximant is grounded solely on the empirical, numerical convergence [9], of the sequences of approximants.

One can always attempt to extrapolate the perturbative results by means of the Padé approximants *PM*,*<sup>N</sup>* (*x*) [6,9]. The Padé approximants *PM*,*<sup>N</sup>* can be understood as the ratio of two polynomials *PM*(*x*) and *QN*(*x*) of the order *M* and *N*, respectively. The diagonal Padé approximant of order *N* corresponds to the case of *M* = *N*. Conventionally, *QN*(0) = 1. The coefficients of the polynomials are derived directly from the asymptotic equivalence with the given power series for the sought function Φ(*x*). Sometimes, when there is a need to stress the role of Φ(*x*), we write *PadeApproximant* [Φ[*x*], *n*, *m*].

The Padé approximant might possess a pole associated with a finite critical point, but can only produce an integer critical index. While usually critical indices are not integers. The same concerns the large-variable behavior where the power of *x* produced from extrapolation with some form of Padé approximants is always an integer. Unfortunately, solutions to many problems exhibit irrational functional behavior. Such a behavior cannot be properly described by the standard rational Padé approximants. However, it would be highly desirable to modify somehow the familiar technique of Padé approximants in order to take into account the irrational behavior. Such modification can be performed by separating the sought modification of the Padé approximants into two factors [10]. The first factor is to be expressed as an iterated root or factor approximant [11,12]. It is specifically designed to take care of the irrational part of the solution. The second factor is simply a diagonal Padé approximant, and it is supposed to take care of the rational part of the solution. We arrive thus to the corrected Padé approximants. They appear to be applicable to a larger class of problems, even when the standard Padé technique is not applicable [11].

Many examples of application of the Padé approximants as well as their theoretical modifications, can be found in [13], including some important applications to aerodynamics and boundary layer problems [14]. The so-called two-point Padé is applied for interpolation, when in addition to the expansion about *x*<sup>0</sup> = 0, given by (1), additional information is available and contained in the asymptotic power series expansion about *<sup>x</sup>* <sup>=</sup> <sup>∞</sup>, <sup>Φ</sup>(*x*) <sup>∼</sup> <sup>∑</sup><sup>∞</sup> *<sup>n</sup>*=<sup>0</sup> *bnx*−*<sup>n</sup>* [8]. The two-point Padé approximant has the same form as the standard Padé approximant, but with the coefficients expressed through *cn* and *bn*.

The idea of combining information coming from the different limits appear to be fruitful and can be exploited for different types of approximants and various forms of asymptotic expansions [11,12,15,16]. Various self-similar approximants also allow extrapolating and interpolating between the small-variable and large-variable asymptotic expansions, as discussed recently in [16]. The key to the success is to introduce the so-called control functions to allow "to sew" the two limit-cases together in the form most natural for each concrete problem [11,12,15–17]. The example of such an approach is brought up in Appendix B. Although the expansions for small and large couplings are very bad, the resulting approximants are in a good agreement with the numerical data.

There are four main technical approaches to the approximants constructions, all aimed to optimize their performance. The first approach is conventional, also called accuracy-through order. It is based on progressive improvement of quality of approximants with adding new information through the higher-order coefficients, with the approximants becoming more and more complex. It is exemplified in construction of Pad'e and Euler super-exponential approximants [8], factor, root and additive approximants [11,12,16]. The latter "cluster" of approximations was derived based on the ideas of self-similar approximation theory, a close relative of the field-theoretic renormalization group [17]. The property of self-similarity is discussed in Section 3.1.

The second approach leads to corrected approximants. The idea is to ensure the correct form of the solution already in the starting approximation with some initial parameters. The initial parameters should be corrected by asymptotically matching with the truncated series/polynomial regressions in increasing orders. Thus, instead of increasing the order of approximation, one can correct the parameters of the initial approximation [11,12]. The form of the solution is not getting more complex, but the parameters take more and more complex form with increasing order.

In the third approach, predominantly adopted in Section 3, we keep the form and order of approximants the same in all orders, but let the series/regressions evolve into higher orders. Independent on the order of regression, we construct the same approximant, based on the first-order terms solely, only with the parameters changing with increasing order of regression. In the framework of such effective first-order theories, we employ exponential approximants and their extensions.

In the fourth approach, the critical index is treated as a vital part of optimization procedure. The critical index plays the role of a control parameter, to be determined from the optimization procedure described in Section 2.2, following Gluzman and Yukalov [18]. Different optimization techniques based on introduction of control parameters were proposed in [19,20].

The problems arising in approximation theory can vary. Note that, for a recovery problem, when measurements of the sought function are given for some finite set of points, there is Prony's method available, with the sought function represented as sums of polynomial or exponential functions combined with periodic functions [21]. For approximation of a continuous function on the interval *x* ∈ [0, 1], one can use Bernstein polynomials [22]. However, the two methods do not allow for inclusion of the asymptotic information.

Prony's and Bernstein methods are numerical and work only for interpolation problems. The latter method was further adapted to the region *x* ∈ [0, ∞), and applied in [23]. The technique of Cioslowski [23] allows for incorporation of the asymptotic information. The technique of self-similar roots [24] allows us to solve the same problems as in [23], but without resorting to fitting [23,24].

Our methods are analytical, user-friendly and applicable to the most difficult extrapolation problem [11,12,16], involving explicit calculation of various critical indices and amplitudes, with novel applications to finding relaxation times. However, our methods remain applicable also for various interpolation problems [11,12,16,24] (see also Appendix B).

It is likely impossible to find the same approximant to be the best for each and every realistic problem. Based on the same asymptotic information, such as series coefficients, thresholds, critical indices, and correction to scaling indices, one can construct not only Padé but quite a few different approximants, such as corrected Padé, additive, *DLog*-additive, etc. [12,16]. It is feasible that for each problem one can find an optimal different approximant. We think that the idea behind the method of corrected approximants [11,12,16], is the most progressive, since it allows to combine the strength of a few methods together and proceed, in the space of approximations, with piece-wise construction of the approximation sequences, as pointed out recently by Gluzman [16].In the following sections, we present a more expended description of the concept of approximants, applied now both to critical and relaxation phenomena, extending the earlier work of Chapter 1 of the book [12].

### **2. Critical Index and Relaxation Time**

The function Φ(*x*) of a real variable *x* exhibits critical behavior, with a critical index *α*, at a finite critical point *xc*, when

$$\Phi(\mathbf{x}) \cong A(\mathbf{x}\_{\mathfrak{c}} - \mathbf{x})^{a}, \text{ as } \mathbf{x} \to \mathbf{x}\_{\mathfrak{c}} - \mathbf{0} \text{ .} \tag{2}$$

The definition covers the case of negative index when function can tend to infinity, or the sought function can tend to zero if the index is positive. Sometimes, the values of critical index and critical point are known from some sources, and the problem consists in finding the critical amplitude *A*, as extensively exemplified in [11].

The case when critical behavior occur at infinity,

$$\Phi(\mathbf{x}) \simeq A\mathbf{x}^{\mathbf{a}}, \text{ as } \mathbf{x} \to \infty \text{ }, \tag{3}$$

can be analyzed similarly. It can be understood as the particular case with the critical point positioned at infinity.

Critical phenomena are ubiquitous [18], ranging from the field theory to hydrodynamics. It is vital to explain related critical indices theoretically. Regrettably, for realistic physical systems, one can as a rule learn only its behavior at small variable,

$$
\Phi(\mathbf{x}) \simeq \Phi\_k(\mathbf{x}), \text{ as } \mathbf{x} \to \mathbf{0} \text{ ,} \tag{4}
$$

which follows form some perturbation theory. The function Φ*k*(*x*) is approximated by an expansion

$$\Phi\_k(\mathbf{x}) = 1 + \sum\_{n=1}^k c\_n \mathbf{x}^n. \tag{5}$$

Most often one finds that such expansions give numerically divergent results, valid only for very small or very large *x* (see Appendix B). Constructively, the expansion is treated as a polynomial of the order *k*. Sometimes, theoretically, one even has a convergent series, resulting in a rather good numerically convergent, truncated polynomial approximations (A6). However, there is still a problem of extrapolating outside of the region of numerical convergence, where the critical behavior sets in. Three examples of such type are given in Appendix A, based on the results of Chapter 7 of the book [12].

The discussion below traces the basic ideas from Chapter 1 of the book [12]. One can always express the critical index directly by using its definition, and find it as the limit of explicitly expressed approximants. For instance, critical index can be estimated from a standard representation as the following derivative

$$\mathcal{B}\_{\mathfrak{a}}(\mathbf{x}) = \partial\_{\mathbf{x}} \log \left( \Phi(\mathbf{x}) \right) \simeq \frac{-\mathfrak{a}}{\mathfrak{x}\_{\mathfrak{c}} - \mathfrak{x}'} \tag{6}$$

as *x* → *xc*, thus defining the critical index as the residue in the corresponding single pole. The pole corresponds to the critical point *xc*. The critical index corresponds to the residue

$$\alpha = \lim\_{\mathbf{x} \to \mathbf{x}\_{\mathbf{c}}} (\mathbf{x} - \mathbf{x}\_{\mathbf{c}}) \mathcal{B}\_{\mathbf{d}} \ (\mathbf{x}) .$$

To the *DLog*-transformed series B*<sup>a</sup>* (*x*) one is bound to apply the Padé approximants [6]. Moreover, the whole table of Padé approximants can be constructed [9], That is, the *DLog* Padé method does not lead to a unique algorithm for finding critical indices. procedure. Basically, different values are produced by different Padé approximants. Then, it is not clear which of these estimates to prefer. The standard approach consists in applying a diagonal Padé approximants [6].

When a function, at asymptotically large variable, behaves as in (3), then the critical exponent can be defined similarly, by means of the *DLog* transformation. It is represented by the limit

$$\mathfrak{a} = \lim\_{\mathbf{x} \to \infty} \mathfrak{x} \mathcal{B}\_{\mathfrak{a}}(\mathbf{x}) \; . \tag{7}$$

Assume that the small-variable expansion for the function *Ba* (*x*) is given. In order for the critical index to be finite, it is necessary to take only the approximants behaving as *<sup>x</sup>*−<sup>1</sup> as *<sup>x</sup>* <sup>→</sup> <sup>∞</sup>. It leaves us no choice but to select the non-diagonal *Pn*,*n*+1(*x*) approximants, so that the corresponding approximation *α<sup>n</sup>* is finite. One can also apply, in place of Padé, some different approximants [12,16]. The examples of application of the *DLog* Padeé methods are given in Appendix A, based on the results first obtained in Chapter 1 of the book [12].

To simplify and standardize calculations different, and more powerful, approximants, called self-similar factor approximants, are introduced in [25]. The singular solutions emerging from factor approximants correspond to critical points and phase transitions [25], including also the case of singularity located at ∞. When the series is long, one would expect that the accuracy is going to improve with increasing numbers of terms. Sometimes, an optimum is achieved for some finite number of terms, reflecting the asymptotic nature of the underlying series. It is very difficult to improve the quality of results produced by the factor approximants, when the series are short. Some suggestions on such improvement were advanced by Gluzman [12].

In some simple but rather important cases of ODEs, the factor approximants allow to restore exact solutions, such a bell soliton, kink soliton, logistic equation solution and instanton-type solution [26]. However, as pointed out in the Introduction, such cases are quite special, and only an approximate solution could be found in many important cases [26,27]. More information about various methods of calculating critical index, amplitude and critical point can be found in [11,12,16].

### *2.1. Relaxation Time*

Consider the case of relaxation behavior when a function at asymptotically large variable decays as

$$\Phi(t) \simeq A \exp(\frac{t}{\tau}) \qquad (t \to \infty) \,, \tag{8}$$

with negative *τ*. Formally, the relaxation time is −*τ*. It can be found as the limit

$$\frac{1}{\tau} = \lim\_{t \to \infty} \frac{d}{dt} \,\, \ln \Phi(t) \,\,. \tag{9}$$

As in the case of critical behavior considered above, the small-variable expansion for the function is given by the sum Φ*k*(*t*). The effective relaxation time can be expressed in terms of the small-variable expansion as follows,

$$\frac{1}{\tau\_k(t)} = \frac{d}{dt} \,\,\ln \Phi\_k(t) \,\,. \tag{10}$$

It can be expanded in powers of *t*, leading to

$$\pi\_k(t) = \sum\_{n=0}^k b\_n t^n. \tag{11}$$

The coefficients *bn* are easily expressed through *cn* of the original series (1). Let us apply to the obtained expansion the self-similar or Padé approximants, That is, we have to derive an approximant *τ*∗ *<sup>k</sup>* (*t*) whose limit

$$
\pi\_k^\*(t) \to \text{const} \qquad (t \to \infty) \text{ ,}
$$

gives the relaxation time

$$\pi\_k^\* = \lim\_{t \to \infty} \pi\_k^\*(t) \,. \tag{12}$$

In such approach, the amplitude *A* does not enter the consideration. In practice, one can indeed construct the approximants with such required behavior. The complete approximant for the sought function Φ(*t*) denoted below as *E*(*t*,*r*), can be constructed as well. Even some *ad hoc* forms satisfying some general symmetry requirements can be suggested, as in Section 3.

As an illustration, let us find *τ*∗ *<sup>k</sup>* (*t*) in explicit form under some simple assumptions concerning its asymptotic behaviors. Assume simply that there are two distinct exponential behaviors for short and long times with two different *τ*1, *τ*2, and the transition from short to long time behavior also occurs at the duration of some third characteristic time *<sup>τ</sup>*<sup>3</sup> <sup>=</sup> <sup>−</sup>*β*−<sup>1</sup> <sup>3</sup> . The characteristic times can be found from the short-time expansion. The simple approximation to the effective relaxation time, expressed in second order of (12), can be written down in the spirit of Yukalov and Gluzman [28] as follows:

$$
\pi\_2 \mathbf{v}^\*(t)^{-1} = \beta\_2 + (\beta\_1 - \beta\_2) \exp\left(\beta\_3 t\right),
\tag{13}
$$

so that for negative *β*<sup>3</sup> we have *τ*∗ <sup>2</sup> (0)−<sup>1</sup> = *<sup>β</sup>*1, *<sup>τ</sup>*<sup>∗</sup> <sup>2</sup> (∞)−<sup>1</sup> = *<sup>β</sup>*2.

In the theory of reliability, the failure (hazard) rate or mortality force [29] is analogous to the inverse effective relaxation time, and the model of the type of formula (13) is known as the Gompertz–Makeham law of mortality.

The complete approximant corresponding to (13) is reconstructed after elementary integration

$$F(t) = A \exp\left(\frac{(\beta\_1 - \beta\_2) \exp(\beta\_3 t)}{\beta\_3} + \beta\_2 t\right),\tag{14}$$

with all unknown constituents of (13) expressed explicitly, from the asymptotic equivalence with the power-series,

$$\begin{array}{ll} \mathcal{A} = \mathfrak{c}\_{0} \exp\left(\frac{\left(\mathfrak{c}\_{1}^{2} - 2\mathfrak{c}\_{0}\mathfrak{c}\_{2}\right)^{3}}{4\left(3\mathfrak{c}\_{0}^{2}\mathfrak{c}\_{3} - 3\mathfrak{c}\_{0}\mathfrak{c}\_{1}\mathfrak{c}\_{2} + \mathfrak{c}\_{1}^{3}\right)^{2}}\right), & \mathcal{J}\_{1} = \frac{\mathfrak{c}\_{1}}{\mathfrak{c}\_{0}}, \quad \mathcal{J}\_{2} = \frac{6\mathfrak{c}\_{0}^{2}\mathfrak{c}\_{1}\mathfrak{c}\_{3} - 4\mathfrak{c}\_{0}^{2}\mathfrak{c}\_{2}^{2} - 2\mathfrak{c}\_{0}\mathfrak{c}\_{1}^{2}\mathfrak{c}\_{2} + \mathfrak{c}\_{1}^{4}}{2\mathfrak{c}\_{0}\left(3\mathfrak{c}\_{0}^{2}\mathfrak{c}\_{3} - 3\mathfrak{c}\_{0}\mathfrak{c}\_{1}\mathfrak{c}\_{2} + \mathfrak{c}\_{1}^{3}\right)}, \\ \mathcal{B}\_{3} = \frac{2\left(3\mathfrak{c}\_{0}^{2}\mathfrak{c}\_{3} - 3\mathfrak{c}\_{0}\mathfrak{c}\_{1}\mathfrak{c}\_{2} + \mathfrak{c}\_{1}^{3}\right)}{\mathfrak{c}\_{0}\left(2\mathfrak{c}\_{0}\mathfrak{c}\_{2} - \mathfrak{c}\_{1}^{2}\right)}. \end{array} \tag{15}$$

Most interesting, as *β*<sup>2</sup> = 0 the linear decay (growth) term in the formula for *F*(*t*) disappears, we arrive in different notations to the Gompertz function (54),

$$G(t) = A \exp\left(\frac{\beta\_1 \exp(\beta\_3 t)}{\beta\_3} t\right),\tag{16}$$

employed in calculations of Gluzman [30]. In this case, we have the effective relaxation time decaying (growing) exponentially with time.In Section 3, we apply this method of finding the effective relaxation time for time series.

### *2.2. Critical Index as Control Parameter. Optimization Technique*

The function's critical behavior follows from extrapolating the asymptotic expansion (1) to finite or large values of the variable. Such an extrapolation can be accomplished by means of a direct technique just discussed above. However, its successful application requires knowledge of a large number of terms in the expansion. However, it is also possible to obtain rather good estimates for the critical indices from a small number of terms in the asymptotic expansion [12,18]. To this end, we can employ the self-similar root approximants given by (17). The external power *mk* is to be determined here from additional conditions. More detailed explanations and more examples can be found in the book [12].

The self-similar root approximant has the following general form [15],

$$\mathcal{R}\_k^\*(\mathbf{x}, m\_k) = \left( \left( (1 + \mathcal{P}\_1 \mathbf{x})^{m\_1} + \mathcal{P}\_2 \mathbf{x}^2 \right)^{m\_2} + \dots + \mathcal{P}\_k \mathbf{x}^k \right)^{m\_k}.\tag{17}$$

In principle, all the parameters may be found from asymptotic equivalence with a given power series.

The large-variable power *α* in Equation (3) could be compared with the large-variable behavior of the root approximant (17),

$$\mathcal{R}\_k^\*(\mathbf{x}, m\_k) \simeq A\_k \mathbf{x}^{km\_k} \,, \tag{18}$$

where

$$A\_k = \left( \left( \left( \mathcal{P}\_1^{m\_1} + \mathcal{P}\_2 \right)^{m\_2} + \mathcal{P}\_3 \right)^{m\_3} + \dots + \mathcal{P}\_k \right)^{m\_k}.\tag{19}$$

This comparison yields the relation *kmk* = *α*, defining the external power *mk* = *<sup>α</sup> <sup>k</sup>* , when *α* is known. This way of defining the external power is used when the root approximants are applied for interpolation. The root approximants (17) are applied in Appendix B, in the context of interpolation problem, for construction of accurate formulas valid for all values of *x*.

Consider an exceptionally difficult situation of an extrapolation problem: the large-variable behavior of the function is not known and *α* is not given. In addition, the critical behavior can happen at a finite value *xc* of the variable *x*. The method for calculating the critical index *α* by employing the self-similar root approximants was developed by Gluzman and Yukalov [18].

In such approach, we construct several root approximants *R*∗ *<sup>k</sup>* (*x*, *mk*), and the external power *mk* plays the role of a control function. The sequence of approximants is considered as a trajectory of a dynamical system. The approximation order *k* plays the role of discrete time. A discrete-time dynamical system or the approximation cascade consists of the sequence of approximants. The cascade velocity is defined by Euler discretization formula [31–33]

$$\mathcal{V}\_{\mathbf{k}}(\mathbf{x}, m\_{\mathbf{k}}) = \mathcal{R}\_{\mathbf{k}+1}^{\*}(\mathbf{x}, m\_{\mathbf{k}}) - \mathcal{R}\_{\mathbf{k}}^{\*}(\mathbf{x}, m\_{\mathbf{k}}) + (m\_{\mathbf{k}+1} - m\_{\mathbf{k}}) \, \frac{\partial}{\partial m\_{\mathbf{k}}} \, \mathcal{R}\_{\mathbf{k}}^{\*}(\mathbf{x}, m\_{\mathbf{k}}) \,. \tag{20}$$

The effective limit of the sequence of approximants corresponds to the fixed point of the cascade. Based on just a few approximants, the cascade velocity has to decrease. In such a sense, the sequence appears to be convergent. The control functions *mk* = *mk*(*x*), have to minimize the absolute value of the cascade velocity

$$\left|V\_{k}(\mathbf{x},m\_{k}(\mathbf{x}))\right| = \min\_{m\_{k}} \left|V\_{k}(\mathbf{x},m\_{k})\right|.\tag{21}$$

A finite critical point *x<sup>c</sup> <sup>k</sup>*, in the *k*th approximation, is to be obtained from Equation (17) by imposing the condition on the critical behavior expressed by (2),

$$[\mathcal{R}\_k^\*(\mathfrak{x}\_{k'}^\varepsilon m\_k)]^{1/m\_k} = 0 \qquad (0 < \mathfrak{x}\_k^\varepsilon < \infty). \tag{22}$$

Its finite solution is denoted as *x<sup>c</sup> <sup>k</sup>* = *<sup>x</sup><sup>c</sup> <sup>k</sup>*(*mk*).

The critical index in the *k*th approximation is given by the limit

$$\mu\_k = \lim\_{\mathbf{x} \to \mathbf{x}\_k^{\varepsilon}} m\_k(\mathbf{x}).$$

In the case of the critical behavior at infinity, when *xc* ∼ ∞, the critical index is

$$\alpha = k \lim\_{\mathbf{x} \to \infty} m\_k(\mathbf{x})\_\prime \text{ as } \mathbf{x}\_\mathbf{c} \sim \infty \text{ .} \tag{23}$$

*Axioms* **2020**, *9*, 126

Thus, to find the critical indices, the control functions *mk*(*x*) have to be found. The minimization of the cascade velocity (50) is complicated. Equation (21) contains two control functions, *mk*<sup>+</sup><sup>1</sup> and *mk*. Nevertheless, the problem can be resolved.

This can be done in two ways. The first constructive approach notices that *mk*<sup>+</sup><sup>1</sup> should be close to *mk*. Then, we arrive to to the minimal difference condition

$$\min\_{m\_k} \left| \mathcal{R}\_{k+1}^\*(\mathbf{x}, m\_k) - \mathcal{R}\_k^\*(\mathbf{x}, m\_k) \right| \qquad (k = 1, 2, \dots) \,. \tag{24}$$

One should typically find a solution *mk* = *mk*(*x*) of the simpler equation

$$
\mathcal{R}\_{k+1}^\*(\mathbf{x}, m\_k) - \mathcal{R}\_k^\*(\mathbf{x}, m\_k) = 0 \; . \tag{25}
$$

The control functions *mk*, characterizing the critical behavior of Φ(*x*), become the numbers *mk*(*xc*). We simply write *mk* = *mk*(*xc*).

In the vicinity of a finite critical point, the function R<sup>∗</sup> *<sup>k</sup>* behaves as

$$\mathcal{R}\_k^\*(\mathbf{x}, m\_k) \simeq \left(1 - \frac{\mathbf{x}}{\mathbf{x}\_k^c}\right)^{m\_k}, \text{ as } \mathbf{x} \to \mathbf{x}\_k^c - \mathbf{0} \text{ .} \tag{26}$$

The condition (25) is expressed as follows,

$$\left(\mathbf{x}\_{k+1}^{\varepsilon}(m\_k) - \mathbf{x}\_k^{\varepsilon}(m\_k)\right) = 0 \qquad \left(0 < \mathbf{x}\_k^{\varepsilon} < \infty\right). \tag{27}$$

For the critical behavior at infinity, it is expedient to introduce the control function

$$\mathbf{s}\_k = km\_k \, . \tag{28}$$

The large-variable behavior reads as

$$\mathcal{R}\_k^\*(\mathbf{x}, \mathfrak{s}\_k) \simeq A\_k(\mathfrak{s}\_k) \mathfrak{x}^{\mathfrak{s}\_k}, \text{ as } \mathbf{x} \to \infty \text{ }. \tag{29}$$

As a result, the minimal difference condition is reduced to the equation

$$A\_{k+1}(\mathfrak{s}\_k) - A\_k(\mathfrak{s}\_k) = 0, \text{ as } \mathfrak{x}\_k^\varepsilon \sim \infty \text{ .} \tag{30}$$

The alternative equation for the control functions also follows from the minimal velocity condition (21), and is called the minimal derivative condition

$$\min\_{k} \left| \frac{\partial}{\partial m\_{k}} \, \mathcal{R}\_{k}^{\*} (\mathbf{x}, m\_{k}) \right| \qquad (k = 1, 2, \dots) \; , \tag{31}$$

In practice, we have to solve the equation

$$\frac{\partial}{\partial m\_k} \, \mathcal{R}\_k^\*(\mathfrak{x}, m\_k) = 0 \, . \tag{32}$$

To apply this condition, we have first to extract from the function its non-divergent parts. If the critical point is finite, one can study the residue of the function *∂* log R<sup>∗</sup> *<sup>k</sup>*/*∂mk*, expressed as

$$\lim\_{\mathbf{x}\to\mathbf{x}\_{k}^{c}} (\mathbf{x}\_{k}^{c} - \mathbf{x}) \, \frac{\partial}{\partial m\_{k}} \, \log \mathcal{R}\_{k}^{\*} (\mathbf{x}, m\_{k}) = m\_{k} \, \frac{\partial \mathbf{x}\_{k}^{c}}{\partial m\_{k}} \, \dots$$

Thus, from Equation (32), we arrive to the condition

$$\frac{\partial \mathbf{x}\_k^{\varepsilon}}{\partial m\_k} = 0 \qquad \left(0 < \mathbf{x}\_k^{\varepsilon} < \infty\right). \tag{33}$$

When the critical behavior occurs at infinity, then we can consider the limiting form of the amplitude and reduce Equation (32) to the form

$$\frac{\partial A\_k(\mathbf{s}\_k)}{\partial \mathbf{s}\_k} = 0, \text{ as } \mathbf{x}\_k^c \sim \infty \text{ .} \tag{34}$$

The final estimate for the critical index is given by a simple average of the minimal difference and minimal derivative results.

The technique reviewed in Section 2.2, following Chapter 1 of the book [12], turned out to be useful in calculating the critical properties of the classical analog of the graphene-type composites with varying concentration of vacancies [34].

In the next subsection, we give some examples, first presented in Chapter 1 of the book [12]. More information and details can also be found in Chapter 7 of the book [12].

### *2.3. Examples: Permeability in the Two-Dimensional Channels*

In the cases considered below, we deal with a unique theoretical opportunity to attack the problem of critical exponent and criticality in general, directly from the solution of the hydrodynamic Stokes problem. Let us consider as example the case of the two-dimensional channel bounded by the surfaces *z* = ±*b* (1 +  cos *x*) , as explained in Appendix A. Here, is termed *waviness*.

The permeability behaves critically [12], That is, it tends to zero as

$$K(\boldsymbol{\epsilon}) \sim (\boldsymbol{\epsilon}\_{\boldsymbol{\epsilon}} - \boldsymbol{\epsilon})^{\varkappa}, \text{ as } \boldsymbol{\epsilon} \to \boldsymbol{\epsilon}\_{\boldsymbol{\epsilon}} - \boldsymbol{0} \text{ , } \tag{35}$$

with *<sup>c</sup>* = 1 ,κ = <sup>5</sup> <sup>2</sup> . The permeability as a function of the waviness can be derived in the form of an expansion in powers of *-* [35]. In the particular case of *b* = 0.5, the permeability can be found explicitly as

$$K(\epsilon) \simeq 1 - 3.14963 \; \epsilon^2 + 4.08109 \; \epsilon^4, \; \text{as } \epsilon \to 0 \; . \tag{36}$$

By setting *<sup>c</sup>* = 1, and changing the variable *y* = *-*2 1−*-*<sup>2</sup> , one can move the critical point to infinity.

The critical index is calculated as explained above and in [18]. From the minimal-difference condition we find κ<sup>1</sup> = 2.184, with an error 12.6%. From the minimal derivative condition, we obtain κ<sup>2</sup> = 2.559, with an error 2.37%. The final answer κ<sup>∗</sup> is given by the average of two solutions <sup>κ</sup><sup>∗</sup> <sup>=</sup> 2.372 <sup>±</sup> 0.19 .

In another particular case considered in Chapter 1 of the book [12], for *b* = 0.25, the permeability expands as follows,

$$K(\epsilon) \simeq 1 - 3.03748 \; \epsilon^2 + 3.54570 \; \epsilon^4, \; \text{as } \epsilon \to 0 \; . \tag{37}$$

Setting *-* = 1, and using the same technique as above the approximations for critical index are found, so that <sup>κ</sup><sup>1</sup> <sup>=</sup> 2.342, and <sup>κ</sup><sup>2</sup> <sup>=</sup> 2.743. Finally, <sup>κ</sup><sup>∗</sup> <sup>=</sup> 2.543 <sup>±</sup> 0.2.

Let us also consider some examples of the numerical convergence of root approximants in high-orders, first presented in Chapter 1 of the book [12]. The technique is applied for calculating critical index κ. It seems instructive to consider the same two cases of permeability *K*(*-*), but with higher-order terms, up to 16th order inclusively.

The numerical form of the corresponding expansions can be found in Appendix A (see expansions (A8) and (A14)). Concretely, we construct the iterated root approximants

$$\mathcal{R}\_k^\*(y) = \left( \left( \left( (1 + \mathcal{P}\_1 y)^2 + \mathcal{P}\_2 y^2 \right)^{3/2} + \mathcal{P}\_3 y^3 \right)^{4/3} + \dots + \mathcal{P}\_k y^k \right)^{a/k} \,. \tag{38}$$

The parameters P*<sup>j</sup>* have to be found from the asymptotic equivalence with the expansions. The permeability has the required critical asymptotic forms

$$\mathcal{R}\_k^\*(y) \simeq A\_k y^\mu, \text{ as } y \to \infty \text{ .} \tag{39}$$

The amplitudes *Ak* = *Ak*(*αk*) are found explicitly as

$$A\_k = \left( \left( (\mathcal{P}\_1^2 + \mathcal{P}\_2)^{3/2} + \mathcal{P}\_3 \right)^{4/3} + \dots + \mathcal{P}\_k \right)^{a/k} . \tag{40}$$

To define the critical index *αk*, we analyze the differences

$$
\Delta\_{kn}(\mathfrak{a}\_k) = A\_k(\mathfrak{a}\_k) - A\_{\mathfrak{n}}(\mathfrak{a}\_k) \,. \tag{41}
$$

From the sequences Δ*kn* = 0, we find the related sequences of approximate values *α<sup>k</sup>* for the critical indices.

Although it is possible to investigate different sequences of the conditions Δ*kn* = 0, the most natural from is presented by the sequences of Δ*k*,*k*+<sup>1</sup> = 0 and Δ*k*<sup>8</sup> = 0, with *k* = 1, 2, 3, 4, 5, 6, 7.

The results for *b* = <sup>1</sup> <sup>2</sup> are shown in Table 1. We observe good numerical convergence of the approximations *<sup>α</sup><sup>k</sup>* <sup>≡</sup> <sup>κ</sup>*k*, to the value <sup>κ</sup> <sup>=</sup> <sup>5</sup> 2 .

Similar results, presented in Table 2 (for *b* = <sup>1</sup> <sup>4</sup> ), again demonstrate rather good numerical convergence of the approximate critical indices to the value κ = <sup>5</sup> 2 .

Comparison of the results for different parameters *b* allows us to think that the critical index does not depend on parameter *b*. In both examples considered above, the convergence sets in rather quickly.

The *DLog* Padé method appears to bring convergent sequences and consistent expressions for permeability as well. Further details can be found in Appendix A. The results obtained from the two different methods well agree with each other. A similar comparison was made by Gluzman and co-authors [34] for the effective conductivity of graphene-type composites.

**Table 1.** Walls can touch (*b* = 1/2). The problems described in Appendices A and A.1. Critical indices for the permeability κ*<sup>k</sup>* obtained from the optimization conditions (41). There is rather good numerical convergence to the number 5/2.



**Table 2.** Walls can touch (*b* = 1/4). The problems described in the Appendices A and A.2. Critical indices κ*<sup>k</sup>* are found from the optimization conditions (41). There is a good numerical convergence of the sequences to the value 5/2.

Consider a different case of permeability *K*(*-*) (see Appendixes A and A.3). The results were first obtained in Chapter 1 of the book [12]. For the parallel sinusoidal two-dimensional channel when the walls would not touch, the permeability remains finite. It is expected to decay as a power-law as the waviness becomes large,

$$K(\epsilon) \sim \epsilon^{\upsilon}, \text{ as } \epsilon \to \infty,$$

with negative index *ν*.

In the expansion of *K*(*-*) in small parameter *-*2, we retain the same number of terms as in the previous two examples. The numerical values of the corresponding coefficients can be found in Appendix A ( see expression (A16)). The results of calculations are presented in Table 3 (for *b* = <sup>1</sup> 2 ). They show rather good numerical convergence, especially in the last column, to the value −4. The sequence, based on the *DLog* Padé method, is convergent as well (see Appendixes A and A.3).

**Table 3.** Walls can not touch. Case of *b* = 1/2. Critical indices for the permeability for the problems in Appendixes A and A.3, obtained from the optimization conditions Δ*kn*(*νk*) = 0. The sequences demonstrate reasonably good numerical convergence to the value *ν* = −4.


More information on the problems of critical permeability, can be found in Appendix A. The three problems considered above are studied by applying the *DLog* Padé method of Section 2 to calculate the critical index for permeability. The computations complement and confirm the results for critical index, obtained above from the optimization technique. The optimization technique works better for short truncated series, converging more quickly, while the *DLog* Padé method is easier to apply for very long series. In addition, the *DLog* Padé method, as well as the Padé method, when its application is appropriate, allows us to compute the critical amplitudes.

### **3. Relaxation Phenomena in Time Series**

For the phenomenon to occur, the basic underlying symmetry must be broken. While studying the phenomenon it is important to distinguish between an explicit symmetry breaking when governing equations are not invariant under the desired symmetry and spontaneous symmetry breaking, without presence of any asymmetric cause [36]. When successful, the approach based on broken global symmetries leads to understanding of the key phenomena of magnetism, superconductivity and superfluidity. On the other hand, when some global inherent symmetry can be recognized in physical quantities, we arrive to the gloriously successful theory of critical phenomena and vital extensions of perturbation results in quantum field theories, jointly called renormalization group (RG) [17,37]. In a nutshell, we suggest below how to apply symmetry considerations and RG-inspired methods to the sharp moves which occur in time series, with the most notable examples given by stock market crashes.

Assume that numerical data on the time series variable (e.g., price) *s* is given for some time *t* segment. Typically, one considers *N* + 1 values *s*(*t*0),*s*(*t*1)... ,*s*(*tN*), for *N* + 1 given at equidistant successive moments in time *t* = *tj*, with *j* = 0, 1, 2 . . . , *N* [38].

In the study of time series, one is interested in the extrapolated to future value of *s*. In financial mathematics, one is particularly interested in the predicted value of log return [38,39],

$$R(t\_N + \delta t) = \ln \left( \frac{s(t\_N + \delta t)}{s(t\_N)} \right). \tag{42}$$

One can see from the definition that we are really interested in the quantity **S** = ln(*s*), to be called return. Let us place the origin at the very beginning of the time interval, setting also *t*<sup>0</sup> = 0. Naturally, one is interested in the value of **S**(*tN* + *δt*), allowing to find *R*(*tN* + *δt*) at a later time. Since the approach developed in [30,38] is invariant with regard to the time unit choice, we consider temporal points of the dataset as integer, while considering the actual time variable as continuous.

Modern physics when applied to financial theory is concerned with ergodicity violations [40–43]. Ergodicity violations may be understood as a manifestation of a non-stationarity, or violation of time-invariance of random process. Metastable phases in condensed matter also defy ergodicity over long observation timescales. In special quantum systems of ultracold atoms, spontaneous breaking of time-translation symmetry causes the formation of temporal crystalline structures [44]. The concept of a spontaneously broken time-translation invariance can be useful for time series in application to market dynamics, as first suggested in [38]. According to Andersen, Gluzman and Sornette [38], the window of forecasting of time series describing market evolution emerges due to a spontaneous breaking/restoration of the continuous time-translation invariance, dictated by relative probabilities of the evolution patterns [45]. In turn, the probabilities are derived from the stability considerations.

The notion of probability introduced in [45] is not based on the same conventional statistical ensemble probability for a collection of people, but it is closer to the time probability, concerned with a single person living through time (see Gell-Mann and Peters [42] and Taleb [43]). Probabilistic trading patterns correspond to local breakdown of time-translation invariance. Their evolution leads to the time-translation symmetry complete (or partial) restoration. We need to estimate typical time, amplitude and direction for such a restorative process. Thus, we are not confined to a binary outcome as in [38] but attempt to estimate also the magnitude of the event.

According to Hayek [46], markets are mechanisms for collecting vast amounts of information held by individuals and synthesizing it into a useful data point [46,47], e.g., price of the stock market index dependent on time. Conversely, consolidation of knowledge is done via prices and arbitrageurs (Taleb on Hayek).

A catastrophic downward acceleration regime in the time series is known as crash [48]. Time series representing market price dynamics in the vicinity of crisis (crash, melt-up), could be treated as a self-similar evolution, because of the prevalence of the collective coherent behavior of many trading, interacting agents [45,49], including humans and machine algorithms. The dominant

collective slow mode corresponding to such behavior, develops according to some law, formalized as a time-invariant, self-similar evolution. Away from crisis, there is a superposition of collective coherent mode (generalized trend) and of a stochastic incoherent behavior of the agents [39,45]. We do not attempt here to write down a generic evolution equation of behind the time series pertaining to market dynamics. Instead. we consider, locally in time, some trial functions—approximants—in the form inspired by the solutions to some well-known evolution equations. The approximants are designed to respect or violate self-similarity. If in physics the relation of phenomenon and symmetry violation is understood, in econophysics such connection is far from being clear. However, to realize the promise of econophysics [50], on a consistent basis and at par with physics achievements, one has to identify and study the phenomenon from the relevant symmetry viewpoint. Our primary goal here is not forecasting/timing the crash, but studying the crash as a particular phenomenon created by spontaneous, time-translation symmetry breaking/restoration.

Since the market dynamics is believed to be formed by a crowd (herd) behavior of many interacting agents, there are ongoing attempts to create empirical, binary-type prediction markets functioning on such principle, or mini Wall Streets [47]. Prediction markets often work pretty well, however there are many cases when they give wrong prediction or do not make any predictions at all. Such special set-ups are already very useful in reaching understanding that market crowds are correct only if they express a sufficient diversity of opinion. Otherwise, the market crowd can have a collective breakdown, i.e., is fallible, as expected by Soros [48]. In our understanding, such breakdowns amount to breaking of time-translation invariance. Restoration of the time-translation invariance—in theory—may be attributed to a small proportion of the traders having either superior information or market intellect [47]. Data from a survey conducted with high income and institutional investors show that they "generally exaggerated assessments of the risk of a stock market crash, and that these assessments are influenced by the news stories, especially front page stories, that they read" [51]. The division into two (at least) groups can be seen in the very parallel existence of future and spot markets for the same asset, such as S&P 500 index, with the futures market working 24 h. It is believed that a lot of the daily crashes, or melt-up days, start overnight. It is not that arbitrage is not effective, the spot market is just closed overnight, while the futures market operates in a discovery mode.

### *3.1. Self-Similarity and Time Translation Invariance*

According to Isaac Newton and Murray Gell-Mann, the laws of nature are somehow self-similar. The laws of Newtonian mechanics are invariant with respect to the Galilean group, expressing Galileo's principle of relativity [52]. The group includes time-translation invariance, or else the laws of classical mechanics are self-similar.

What should be the underlying symmetry for price dynamics? Mind that in normal times the average price trajectory is exponential, because of the compounding interests, and we enjoy an almost constant return (or price growth rate) [53]. Indeed, let *st*<sup>0</sup> be an underlying security (index) price at *t* = *t*0. Let *F<sup>P</sup> <sup>t</sup>* be the fair value of the future requiring a risk associated expected return *β* [43]. Then (see, e.g., [43]), expected forward price *F<sup>P</sup> <sup>t</sup>* = *st*<sup>0</sup> exp(*β*(*t* − *t*0). For example, a share of a stock would be correctly priced with the expected return calculated as the return of a risk-free money market fund minus the payout of the asset, being a continuous dividend for a stock [43]. Thus, rather simple and natural exponential estimates are constantly made for stocks and the alike. The formula for the forward price is self-similar, or time-translation invariant, as explained below.

However, as noted in [48,53], prices often significantly deviate from such a simple description. Bubbles can be formed, as well as other presumed patterns of technical analysis. Asset prices strongly deviate from the fundamental value over significant intervals of time. The fundamental value is not truly observable, making definition of such intervals somewhat elusive. There are some very real mechanisms in work, acting to increase and even accelerate the deviation from fundamental value. The causes of deviation could be "option hedging, portfolio insurance strategies, leveraging and margin requirements, imitation and herding behavior", as is the authoritative opinion expressed in [48,53].

Recall also that meaningful technical analysis starts from recasting the time series data using some polynomial representation to serve as the expansion [38]. The regression is constructed in standard fashion by minimizing mean-square deviation, with the effective result that the high-frequency component of the price is getting average out. Then, one can consider self-similarity in averages [49]. Indeed, the standard polynomial regressions are invariant under time-translation, retaining their form after arbitrary selection of origin of time with simple redefinition of all parameters. The position of origin in time can be explicitly introduced into the regression formula and included into the coefficients, but actual results of calculations with any arbitrary chosen origin will remain the same. Such property can be expressed as some symmetry.

We put forward the idea that it is the onset of broken time-translation invariance that signifies the birth of a bubble, or of some other temporal pattern preceding a crash. End of pattern corresponds to the restoration of time-translation invariance, partially or fully. Our task is to express this idea in quantitative terms by making explicit transformation from the regression-based technical analysis to the valuation formula in the exponential form, taking into account strong deviations from the standard valuation formulae.

Assume that a time series dynamics is predominantly governed by its own internal laws. This is the same as to write down a self-similar evolution for the marker price *s* [54], meaning that, for arbitrary shift *τ* , one can see that

$$s(t+\tau, a) = s(t, s(\tau, a)),\tag{43}$$

with the initial condition *s*(0, *a*) = *a* [55,56]. The value of the self-similar function *s* in the moment *t* + *τ* with given initial condition, is the same as in the moment *t*, with the initial condition shifted to the value of *s* in the moment *τ*.

When *t* stands for true time, the property of self-similarity means the time-translation invariance. Formally understood, Equation (43) gives a background for the field-theoretical RG, with addition of some perturbation expansion for the sought quantity, which should be resummed in accordance with self-similarity expressed in the form of ODE [55–57]. The time-translation invariance expressed by (43) means that the law for price evolution exists and remains unchanged with time, with proper transformation of the initial conditions [52]. The role of perturbation expansion when price dynamics is concerned, is accomplished by meaningful technical analysis, by recasting data in the form of some polynomial representation [38]. There is no formal difference in treating polynomials and expansions, as already mentioned in Section 2.

Consider first the simplest case of technical analysis. The linear function can be formally considered as the function of time and initial condition *a*, namely *s*1(*t*, *a*) = *a* + *bt*, and *s*1(0, *a*) = *a*. The linear function (regression) is self-similar, or time-translation invariant, as can be checked directly, by substitution into (43).

Through some standard procedure, let us obtain the linear regression on the data around the origin *t*<sup>0</sup> = 0, so that

$$s\_{0,1}(t) = a\_1 + b\_1 t\_-$$

Note that the position of origin is arbitrary, and it can be moved to arbitrary position given by real number *r*, so that

$$s\_{r,1}(t) = A\_1(r) + B\_1(r)(t - r)\_{\prime\prime}$$

with new and different coefficients. It turns out that the coefficients are related as follows

$$A\_1(r) = a\_1 + b\_1 r, \; B\_1(r) = b\_{1'} $$

so that

$$s\_{r,1}(t) \equiv s\_{0,1}(t).$$

By shifting the origin, we create an *r*-dependent form of the linear regression *sr*,1, which can be used constructively. Thus, instead of a single regression we have its *r*-replicas, equivalent to the original form of regression, and all replicas respect time-translation symmetry. In such a sense, one can speak about replica symmetry. Of course, we would like to avoid such redundancy in data parameterization and to find the origin(s) by imposing some optimal conditions (see Section 3.2).

The position of origin in time can be explicitly introduced into the regression formula and included into the coefficients, but actual results of calculations with any arbitrary chosen origin will remain the same. Such property can be expressed as some symmetry. However, intuitively, one would expect that the result of extrapolation with chosen predictors should be dependent on the point of origin *r*. Indeed, various patterns such as "heads and shoulders", "cup-with-handle"," hockey stick", etc., considered by technical analysts do depend on where the point of origin is placed. In physics, the point of origin (Big Bang) plays a fundamental role. We should find a way to break the replica symmetry.

As discussed above, it is exponential shapes that are natural in pricing. Exponential function

$$E(t, a) = a \exp(\beta t),$$

with initial condition *a* and arbitrary *β* satisfy functional self-similarity as well as the linear functions. It can be replicated as

$$\begin{aligned} E\_r(t) &= a(r) \exp(\beta(t - r)), \\ a(r) &= a \exp(\beta r). \end{aligned} \tag{44}$$

Having *β* dependent on *r* is going to *violate* the time-translation and replica symmetry. Instead of a global time-translation invariance, we have a set of *r* local "laws" near each point of origin. However, having *r* in Formula (44) fixed, by imposing some additional condition, or just being integrated out, should *restore* the global time-translation invariance completely as long as the exponential function is considered. Moreover, stability of the exponential function is measured by the exponential function with the same symmetry (see Formula (46)). Not only is exponential function time-translation invariant, but the expected return *β* has the same property. For exponential functions, the expected (predicted) value of return per unit time exactly equals *β*.

Another simple rational function, known as hyperbolic discounting function [58], *H*(*t*, *a*) = <sup>1</sup> *<sup>a</sup>*−1+*bt* , where *a* is the initial condition and *b* is arbitrary, is time-translation invariant. Note that shifted exponential function *Es*(*t*, *a*) = *c* + (*a* − *c*) exp(*bt*), with initial condition *a* and arbitrary *b* and *c*, is invariant under time-translation as well.

Another interesting symmetry is shape invariance [59], meaning

$$F\_{t+\tau}^P = \mathfrak{m} F\_t^P{}\_{\prime}$$

and an exponential function is shape invariant with m = exp (*βτ*), leaving the expected return unchanged. Keep in mind that our task is to calculate *β* from the time series. In principle, one can think about breaking/restoration of shape invariance, as a guide for construction of the concrete scheme for calculations.

A critical phenomenon, an underlying symmetry of the formula for the observable, is scaling

$$
\phi\_{\lambda t} = \Lambda \phi\_{t\_{\nu}}
$$

where Λ = *φλ*. The class of power laws, *φ<sup>t</sup>* = *t <sup>α</sup>*, with critical index *α*, is scaling-invariant. The central task is to calculate *c*. The statistical renormalization group formulated by Wilson [37] explains well the critical index in equilibrium statistical systems. When information on the critical index is encoded in some perturbation expansion, one can use resummation ideas to extract the index, even for short expansions and for non-equilibrium systems [11,12,18]. Some of the methods are discussed in the preceding section (see also [12,16]).

Working with power-law functions will not leave the return unchanged. However, one can envisage the scheme with broken scaling invariance, as an alternative to the former schemes. The log-periodic solutions extend the simple scaling [60] and are extensively employed in the form of a sophisticated seven-parametric fit to long historical dataset [53], as well as of its extensions [61]. The fit is tuned for prediction of the crossover point to a crash, understood as catastrophic downward acceleration regime [48]. However, one cannot exclude the possibility of the solutions with different time symmetries (scaling and time-invariance, for instance) competing to win over, or to coexist, all measured in terms of their stability characteristics.

Our primary concern is the crash per se, not the regime preceding it. We start analyzing crashes with the polynomial approximation that respects time-translation symmetry, have the symmetry broken, and then restored (completely or partially), by means of some optimization. Such sequence ends with a non-trivial outcome: *β* becomes renormalized *β*(*r*), with *r* being found using the optimization procedure(s) defined below. We discuss in Section 2.1 a general technique for correcting *β* directly, which accounts for higher order terms in regression, making it time-dependent.

In [38], the framework for technical analysis of time series was developed, based on second-degree regression and asymptotically equivalent exponential approximants, with some rudimentary, implicit breaking of the symmetry. We intend to go to higher-degree regressions and develop a consistent technique for explicit symmetry breaking with its subsequent restoration. According to textbooks, the fourth order should be considered as "high". Taleb (see footnote on p. 53 in [43]) also considered models with five parameters as more than sufficient.

### *3.2. Optimization, Approximants, Multipliers*

Higher-order regressions allow for replica symmetry. For instance, the quadratic regression *s*0,2(*t*) = *a*<sup>2</sup> + *b*2*t* + *c*2*t* <sup>2</sup> can be replicated as follows:

$$s\_{r,2}(t) = A\_2(r) + B\_2(r)(t - r) + C\_2(r)(t - r)^2,$$

with

$$A\_2(r) = a\_2 + b\_2r + c\_2r^2,\ \ B\_2(r) = b\_2 + 2c\_2r,\ \mathbb{C}\_2(r) = c\_2.$$

With such transformed parameters, we find that *sr*,2(*t*) ≡ *s*0,2(*t*). In fact, one can still formulate self-similarity analogous to (43), but in vector form with increased number of parameters/initial conditions in place of *a* [57]. However, if only the linear part of quadratic regression, or trend, is taken into account, we return to the conventional functional self-similarity ≡ time-translation invariance, discussed above extensively.

Such effective linear/trend approach to higher-order regressions allows applying the same idea at all orders and observe how the exponential structures change with increasing regression order. Note that, in the course of trading, a common pattern is trend following, which appears to be a collective, self-reinforcing motion that, intuitively, lends itself to a self-similar description. Indeed, some participants are waiting for a market confirmation of the trend before acting on it, which in turn acts as a confirmation for others. Having a universal model explaining this dynamics (if not predicting it) would be quite useful.

To take into account the dependence on origin, the replica symmetry has to be broken. Breaking of the symmetry means the dependence on origin of actual extrapolations with non-polynomial predictors. As the primary predictors, we suggest the simplest exponential approximants considered as the function of origin *r* and time,

$$E\_1^\*(t, r) = A(r) \exp\left(\frac{B(r)}{A(r)} (t - r)\right),\tag{45}$$

independent on the order of polynomial regression. The approximants (45) are constructed by requiring an asymptotic equivalence with the linear part of chosen polynomial regression. If the extrapolations *E*∗ <sup>1</sup> (*tN* + *δt*,*r*) are made by each of the approximants, they appear to be different for various *r*, meaning breaking of the replica symmetry and of the time-translation symmetry. Passage from polynomials to exponential functions leads to emergence of the continuous spectrum of relaxation (growth) times.

To compare the approximants quality, one can look at their stability. Stability of the approximants is characterized by the so-called multipliers defined as the variation derivative of the function with respect to some initial approximation function [45]. Following Yukalov and Gluzman [62], one can take the linear regression as zero approximation and find the multiplier

$$M\_1^\*(t, r) = \exp\left(\frac{B(r)}{A(r)}(t - r)\right). \tag{46}$$

The simple structure of multipliers (46) allows avoiding appearance of spurious zeroes which often complicate analysis with more complex approximants/multipliers.

Because of the multiplicity of solutions, embodied in their dependence on origin, it is both natural and expedient to introduce probability for each solution. As explained in [45], one can introduce

$$Probability \approx \left| M\_1^\*(t, r) \right|^{-1} \prime$$

with proper normalization, as shown below in Formula (48). Probability appears to be of a pure dynamic origin and is expressed only from the time series itself. When the approximants and multipliers of the first order are applied to the starting terms of the quadratic, third- or fourth-order regression, we are confined to *effective first-order models*, with velocity parameter from [38] dependent also on higher-order coefficients and origin.

To make extrapolation with approximants (45), one has still to know the origin. In other words, the time-translation symmetry has to be restored completely or partially, so that a specific predictor with specifically selected origin, or as close as possible to a time-translation invariant form, is devised. Fixing unique origin also selects unique relaxation (growth) time, during which the price is supposed to find a time-translation invariant state.

Exponential functions are chosen above because they are invariant under time translation. Any shift in origins is absorbed by the pre-exponential amplitude and does not influence the return *R*. A similar in spirit view that broken symmetries have to be restored in a correct theory was expressed by Duguet and Sadoudi [63].

In the approach predominantly adopted in this section, we keep the form and order of approximants the same in all orders, but let the series/regressions evolve into higher orders. Independent of the order of regression, we construct the same approximant, based only on the first-order terms, only with parameters changing with increasing order of regression. In the framework of the effective first-order theories, we employ exponential approximants.

Consider the value of origin as an optimization parameter [30]. To find it and restore the time-translation symmetry, we have to impose an additional condition directly on the exponential predictors with known last closing price,

$$E\_1^\*(t\_{N'}r) = s\_{N'} \, . \tag{47}$$

One has to solve the latter equation to find the particular origin(s) *r* = *r*∗. In this case, we consider a discrete spectrum of origins, consisting of several isolated values. To avoid double-counting when the last closing price enters both regression and optimization, one can determine the regression parameters in the segment limited from above by *tN*−1,*sN*−1. Alternatively, one can consider the two ways to define regression parameters and choose the one which leads to more stable solutions. Unless otherwise stated, we consider that such comparison was performed and the most stable way was selected.

The extrapolation for the price is simply *s*(*tN* + *δt*) = *E*<sup>∗</sup> <sup>1</sup> (*tN* + *δt*,*r*∗). The condition imposed by Equation (47) is natural, because then a first-order approximation to Formula (42), *<sup>R</sup>* <sup>≈</sup> *<sup>s</sup>*(*tN*+*δt*)−*s*(*tN*) *<sup>s</sup>*(*tN*) , is recovered (see, e.g., [39]), as one would expect intuitively.

The procedure embodied in (47), leads to a radical reduction of the set of *r*-predictors to just a few. Set of predictors and corresponding to each multiplier, define the probabilistic, poor man's order book. Instead of an unknown to us true numbers of buy and sell orders, we calculate a priori probabilities for the price going up or down and corresponding levels. Target price is estimated through weighted averaging developed in [45,62], in its concrete form (48) given below.

For the sake of uniqueness, one can simply choose the most stable result among such conditioned predictors. One can also consider extrapolation with a weighted average of all such selected solutions. With 1 ≤ *M* ≤ 6 solutions, their weighted average *E*<sup>1</sup> for the time *tN* + *δt* is given as follows,

$$\overline{E}\_{1}(t\_{N} + \delta t) = \frac{\sum\_{k=1}^{M} E\_{1}^{\*}(t\_{N} + \delta t, r\_{k}^{\*}) \left| M\_{1}^{\*}(t\_{N} + \delta t, r\_{k}^{\*}) \right|^{-1}}{\sum\_{k=1}^{M} \left| M\_{1}^{\*}(t\_{N} + \delta t, r\_{k}^{\*}) \right|^{-1}}. \tag{48}$$

Within the discrete spectrum, we can find solutions with varying degrees of adherence to the original data. They can follow data rather closely or be loosely defined by the parameters of regression. The former could be called "normal" solutions, and tend to be less stable, with multipliers ∼1, but the latter are "anomalous" solutions, since they cut through the data and typically are the most stable with small multipliers. Anomalous solutions are crashes (meltdowns) and melt-ups. The typical situation with the solutions in the discrete spectrum is presented in Figure 1. The novel feature introduced through (48) is that averaging is performed over all approximants of the same order, compatible with constraints expressed by (47).

**Figure 1.** All Gompertz approximants corresponding to the discrete spectrum, i.e., solutions to (56) are shown. The most stable downward and less stable upward solutions are shown with solid lines. Three additional solutions are shown as well. The solution shown with the dashed line is closest to the data. The "no-change", practically flat solution, is shown with a dot-dashed line. Another solution, corresponding to moderate growth, is shown with a dotted line. The level *s*<sup>16</sup> = 2746.61 is shown with black line. Several historical data points are shown as well.

One can also integrate out the dependence on origin *r*, considered as a continuous variable, by applying an averaging technique of weighted fixed points suggested in [45]. The dependence on origin enters the integration limit through parameter *T*. Integration can be performed numerically for the simplest exponential predictors according to the formula

$$I(t,T) = \frac{\int\_{t\_0 - T}^{t\_N + T} E\_1^\*(X, t) \left| M\_1^\*(X, t) \right|^{-1} dX}{\int\_{t\_0 - T}^{t\_N + T} \left| M\_1^\*(X, t) \right|^{-1} dX}. \tag{49}$$

To optimize the integral, we have to impose an additional condition on the weighted average/integral. It is natural to force it to pass precisely through the last historical point.

$$I(t\_N, T) = s(t\_N),\tag{50}$$

and solve the latter equation to find the integration limit *T* = *T*∗. The sought extrapolation value for the price *s* is simply *I*(*tN* + *δt*, *T*∗). We prefer to take into account the broadest possible region of integration. Under such conditions, if and when the solution to (50) exists, it is unique. The value of *sN* may enter the consideration twice: in the regression parameters and in the optimization condition (50). To avoid counting twice the last known value *sN*, one can slightly different definition

$$I(t,T) = \frac{\int\_{t\_0 - T}^{t\_{N-1} + T} E\_1^\*(X, t) \left| M\_1^\*(X, t) \right|^{-1} dX}{\int\_{t\_0 - T}^{t\_{N-1} + T} \left| M\_1^\*(X, t) \right|^{-1} dX}. \tag{51}$$

As an additional condition to find origin, one can also consider the minimal difference requirement on the lowest order predictors, as first suggested in [49]. Such approach is analogous to the technique discussed in Section 2.2. However, instead of a critical index, we calculate relaxation time. To this end, one has to construct the second order super-exponential approximant

$$\begin{split} E\_2^\*(t, r) &= A(r) \exp\left(\frac{B(r)(t - r) \exp\left(\frac{C(r)(t - r)\pi(r)}{B(r)}\right)}{A(r)}\right), \\ \tau(r) &= 1 - \frac{B(r)^2}{2A(r)C(r)}. \end{split} \tag{52}$$

and minimize its difference with the simplest exponential approximant in the time of interest *tN* + *δt*. Namely, one has to find all roots of the equation

$$\exp\left(\frac{\mathbb{C}(r)\tau(r)(t\_N + \delta t - r)}{B(r)}\right) = 1,\tag{53}$$

with respect to real variable *r*. Corresponding multiplier

$$M\_2^\*(t, r) = \frac{1}{B(r)} \frac{\partial E\_2^\*(t, r)}{\partial t} \,'$$

can be found as well.

The discrete spectrum optimization seems to be the most natural and transparent. Our goal is to find the approximants and probabilistic distributions in the last available historical point of time series. Crashes are attributed to the stable solutions with large negative *r*, meaning that the origin of time has to be moved to the deep past to explain the crash in near future. Preliminary results of Gluzman [30] suggest that, in the overwhelming majority of cases, a crash is preceded by similar, asymmetric probability pattern(s), of the type shown in the figures below. As noted in [51], Kahneman and Tversky explained that people tend to judge current events by their similarity to memories of representative events.

There are also additional solutions with multipliers of the order of unity, coming from the region of moderate *r*, and it is often possible to find some rather stable upward solution for large positive *r*. One can think that, for such stable time series as describing population dynamics, only the region of moderate *r* gives relevant solutions, while for time series describing price dynamics all types of solutions may exist simultaneously.

Within our approach to constructing approximants, one can also try to exploit the second order terms in regression. Instead of exponential approximants, one should try some other, higher-order approximants, but with time-translation invariance property. Such approximants are presented below. They are considered *ad hoc*, because they can be written in closed form only in special, low-order situations. It is not feasible to extend them systematically into arbitrary high order. Hence, our interest in special forms with desired symmetry. Sometimes, it is even not possible to find stable solutions with a single approximant, but it is still possible with corrected approximants.

*Axioms* **2020**, *9*, 126

Recall that exponential function can be obtained as the solution to simple linear first-order ODE. In the search for second-order approximants with time-translation invariance, we turned to some explicit formulas, emerging in the course of solving some first-order ODE with added nonlinear term with arbitrary positive power, which generalizes ODE for simple exponential growth. It is known as Bertalanffy–Richards (BR) growth model [64,65]. Among its solutions in the case of second-order nonlinear term, there is a celebrated logistic function [64],

$$L(t) = \frac{1}{q\_2 + \frac{(1 - q\_1 q\_2) \exp(-q\_0 t)}{q\_1}} \text{ }.$$

where *q*<sup>1</sup> is the initial condition. The logistic function is widely used to describe population growth phenomena and is also known to be the solution to the logistic equation of growth. The logistic function written in the form *L*(*t*, *q*1), dependent on the initial condition *L*(0, *q*1) = *q*1, with arbitrary *q*0, *q*2, is time-translation invariant. One can also introduce the second-order logistic approximant which generalizes logistic function [30]. In addition to describing situations with saturation at infinity, the logistic approximant include also the case of so-called finite-time singularity, which makes it redundant, since such solutions were axiomatically excluded from the price dynamics [38].

Another solution to the BR model in the case when the nonlinear term has power only slightly differing from unity, is known as Gompertz function [64],

$$G(t) = g\_0 \exp(\mathcal{g}\_1 \exp(\mathcal{g}\_2 t)),\tag{54}$$

used to describe growth (relaxation, decay) phenomena. However, as we demonstrate in Section 2.1, it is possible to explain *G*(*t*) directly from the resummation technique leading to Formula (16), without resorting to BR. Relaxation (growth) time behaves exponentially with time. The Gompertz function is *log*-time-translation invariant.

One can consider the second order Gompertz approximant. It simply generalizes the Gompertz function. Namely, one can find Gompertz approximant in the following form

$$\begin{array}{ll} G(t,r) = \operatorname{g}\_{0}(r) \exp(\operatorname{g}\_{1}(r) \exp(\operatorname{g}\_{2}(r)(t-r))),\\ \operatorname{g}\_{0}(r) = A(r)e^{-\operatorname{g}\_{1}(r)}, \quad \operatorname{g}\_{1}(r) = \frac{\operatorname{B}(r)}{A(r)\chi\_{2}(r)}, \quad \operatorname{g}\_{2}(r) = \frac{2A(r)\mathbb{C}(r) - \operatorname{B}(r)^{2}}{A(r)\mathbb{B}(r)},\end{array} \tag{55}$$

with the multiplier

$$M\_G(t,r) = \frac{g\_0(r)g\_1(r)g\_2(r)e^{(g\_1(r)e^{g\_2(r)(t-r)} + g\_2(r)(t-r))}}{B(r)}.$$

The Gompertz approximant, of course, is not limited to the situations with saturation at infinity, as it can also describe very fast decay (growth) at infinity.

With *r* to be found from some optimization procedure, the return *R* generated by Gompertz approximant is time-translation invariant and has a compact form

$$R(\delta t) = \mathcal{g}\_1(r) \exp(\mathcal{g}\_2(r)(t\_N - r))(\exp(\mathcal{g}\_2(r)\delta t) - 1).$$

For small *δt*, it becomes particularly transparent:

$$R(\delta t) \approx g\_1(r) g\_2(r) \exp(g\_2(r)(t\_N - r)) \times \delta t \equiv \frac{\delta t}{\pi (T\_{N\_I} r)^{\gamma}}$$

with the pre-factor giving the return per unit time. The inverse return per unit time has the physical meaning of the effective time for growth (relaxation)

$$\beta(t,r)^{-1} \equiv \tau(t,r) = \left(\mathfrak{g}\_1(r)\mathfrak{g}\_2(r)\right)^{-1} \exp\left(\mathfrak{g}\_2(r)(r-t)\right).$$

considered at the moment *t* = *TN*. Here, we employ the the effective relaxation (growth) time (see Section 2.1),

$$
\pi(t) = \left(\frac{d}{dt}\,\,\ln G(t)\right)^{-1}\,\,\,\,\,
$$

and replicate it. We find that the return for Gompertz approximant is solely determined by relaxation time

$$\mathbf{S}(t,r) = \frac{1}{\tau(t,r)'} $$

allowing to express the log return in a compact form

$$\mathcal{R}(\delta t) = \mathbf{S}(t\_N + \delta t, r) - \mathbf{S}(t\_{N'}r).$$

Thus, the return for Gompertz approximant appears as purely dynamic quantity, not involving any consent about equilibrium, fundamental value, etc. If relaxation time is found from the data to be very large as it should be close to equilibrium conditions [66], we have no potential for returns, i.e., near-equilibrium yields dull, everyday mundane events that are repetitive and lend themselves to statistical generalizations [48]. If relaxation time is anticipated to be very short, we have potentially huge returns. The far-from-equilibrium conditions give rise to unique, historic events [48], or to some very fast relaxation events/crashes. The latter condition makes real markets fragile [67].

Gompertz approximant can go at infinity faster or slower than exponential, and in some important examples such differences amounting to a few percent, can be detected. The function *g*0(*r*), could be called a gauge function for the price, expressing arbitrariness of choice of the price unit, as it does not enter the return. The time-translation invariance of return and gauge invariance for the price are considered very desirable in price model formulation [38], both properties are pertinent to exponential and Gompertz approximations for the price temporal dynamics.

We are interested in market prices on a daily level, and consider only significant market price drops/crashes with magnitude more than 5.5%. Such magnitude is selected to be comparable to the typical yearly return of Dow Jones Industrial Average index. Typically, a 2% daily move is considered as big, but not at the times of various turmoils.

It is widely accepted in practical finance that asset price moves in response to unexpected fundamental information. The information can be identified as well as the tone, positive versus negative. It is found that news arrival is concentrated among days with large return movements, positive or negative [68]. Spontaneously emerging narratives, a simple story or easily expressed explanation of events, might be considered as largely exogenous shocks to the aggregate economy [51]. Simply put, one should analyze what people are talking about in the search for the source of economic fluctuations. Moreover, as in true epidemics governed by evolutionary biology, mutations in narratives spring up randomly, and if contagious generate unpredictable changes in the economy [51]. As noted by Harmon et al. [69], panic on the market can be due to external shocks or self-generated nervousness.

It is argued [70] that cause and effect can be cleanly disentangled only in the case of exogenous shocks, as it is only needed to select some interesting set of shocks to which price is likely to respond. Effects of positive and negative oil price shocks on the stock price need not be symmetric. In macroeconomics, it is even accepted that only positive changes in the price of oil have important effects. Periods dominated by oil price shocks are reasonably easy to identify, and they can indeed be considered as exogenous as well as, often, strong, although difficult to model. Oil price shocks are the leading alternative to monetary shocks and may very well have similar effects [70].

Our goal here is not to forecast/timing the crash, but to study the crash as a particular phenomenon created by spontaneous, time-translation symmetry breaking/restoration. In essence, we ask the following questions:

1. What probabilistic pattern would an observer see the day before crash,

2. What would be the market reaction (expressed through the index), if we are aware that a Swan of some color has already arrived?

In our opinion, in the presence of a Swan, understood as a shock of unspecified strength, the problem simplifies, because of a reduced set of outcomes, dominated by the most extreme, very stable downward solution. Consider that, in natural sciences, most efforts are dedicated to creating a correct experimental setup. Studying reaction to shock is the only current viable substitute for clean experimental conditions.

### *3.3. Examples*

Consider as example a 7.72% drop in the value of Shanghai Composite index related to the first COVID-19 crash, which occurred on 3 February 2020. With *N* = 15, as recommended in [38], the following data points are available,

$$s\_0 = 3085.2, s\_1 = 3083.79, s\_2 = 3083.41, s\_3 = 3104.8, s\_4 = 3066.89, s\_5 = 3094.88,$$

$$s\_6 = 3092.29, s\_7 = 3115.57, s\_8 = 3106.82, s\_9 = 3090.04, s\_{10} = 3074.08, s\_{11} = 3075.5,$$

$$s\_{12} = 3095.79, s\_{13} = 3052.14, s\_{14} = 3060.75, s\_{15} = 2976.53.$$

The value of *s*<sup>16</sup> = 2746.61 is to be "predicted". From the whole set of daily data, we employ only several values of the closing price. Such coarse-grained description of the time series may be justified if one is interested in the phenomenon not dependent on the fine details, such as crash. In the examples presented below, we keep the number of data points per quartic regression parameter in the range 3–4. Lower order calculations can be found in [30]. Here, we show only the quartic regression

$$s\_{0,4}(t) = a\_4 + b\_4t + c\_4t^2 + d\_4t^3 + f\_4t^4,$$

and based on it optimize approximants and multipliers. It can be replicated as follows:

$$s\_{r,4}(t) = A\_4(r) + B\_4(r)(t - r) + \mathbb{C}\_4(r)(t - r)^2 + D\_4(r)(t - r)^3 + F\_4(r)(t - r)^4,$$

with

$$A\_4(r) = a\_4 + b\_4r + c\_4r^2 + d\_4r^3 + f\_4r^4,\ \ B\_4(r) = b\_4 + 2c\_4r + 3d\_4r^2 + 4f\_4r^3,$$

$$\mathbb{C}\_4(r) = c\_4 + 3d\_4r + 6f\_4r^2,\ \quad D\_4(r) = d\_4 + 4f\_4r,\ \quad \mathbb{F}\_4(r) = f\_4.$$

With such transformed parameters, we have *sr*,4(*t*) ≡ *s*0,4(*t*).

Within the data shown in Figure 2, one can discern competing trends. First, let us show the data compared to the regression. There are two obvious trends, "up" and "down", as can be seen in Figure 2.

**Figure 2.** COVID-19, Shanghai Composite, 3 February 2020. Fourth-order regression is shown against data points.

Our analysis indeed finds highly probable solutions of both types, with the downward trend developing into fast exponential decay. Let us analyze the typical approximant and multiplier dependencies on origin, for fixed time *t* = *tN*. The inverse multiplier is shown as a function of the origin *r* in Figure 3 as well as the first-order approximant.

**Figure 3.** Shanghai Composite, 3 February 2020. Calculations with fourth-order regression. The inverse multiplier is shown as a function of the origin *r* at *t* = *TN*, *N* = 15. The first-order approximant is shown in a separate figure. Level *s*<sup>15</sup> is shown as well, with dot-dashed line.

There are two uneven humps in the probabilistic inverse multiplier, suggesting that large negative and large positive *r* dominate, with more weight put on the negative region. Such dependence on *r* manifests the time-translation invariance violation, which should be lifted by finding appropriate origin. More details on the example can be found in [30]. Below, we discuss only the fourth-order calculations.

The results of extrapolation by method expressed by Equation (47) is given as

$$E\_1^\*(16) = 2804.32, \quad M\_1^\*(16) = 0.0113494,$$

with relative percentage error of 2.1%. There is also a less stable "upward" solution

$$E\_1^\*(16) = 3211.95, \quad M\_1^\*(16) = 0.0363796,$$

in agreement with intuitive picture based on naive data analysis. There are also two additional solutions in between with multipliers close to 1. They do not affect averages much, but in real time the metastable solutions, similar to the metastable phases in condensed matter, may show up under special conditions. Metastable solutions when realized violate the principle of maximal stability over the observation timescale, complicating or even negating a unique forecast, based on weighted averages or the most stable solution.

Calculation of the discrete spectrum can be extended to different approximants. For instance, one can also construct the second-order Gompertz approximant introduced above, and solve the following equation on origins:

$$G(t\_N, r) = s(t\_N) \,. \tag{56}$$

The most stable Gompertz approximant gives the most accurate estimate

$$G(16) = 2746.05, \quad M\_G(16) = 0.001539,$$

with a very small error of 0.02%. There are altogether five solutions to (56), in the discrete spectrum, as shown in Figure 1.

Thus, the Gompertz approximant of second order with log-time-translation invariance gives better results than symmetric exponential approximant *E*∗ <sup>1</sup> . Although Taleb's Black Swan did seem to materialize, the short-time stock market response was not different than in somewhat comparable instances of crashes brought up in [30], making it look like a Grey Swan. Indeed, it is plausible that the holiday season in China played the role here. It also helped our cause, effectively pinpointing the day for crash. One can think that all solutions, except the most extreme downward solution, were simply not considered.

Consider several most spectacular examples of crashes from the tumultuous spring and summer of 2020, caused by combination of economic causes such as oil anti-shock and COVID-19 related, enormous disruptions—a rare constellation of Two Swans of Gray coming together! There was a month long delay until DJ crashed. All three conspicuous crashes from March 2020 can be considered as an exponentially accelerated decay.

**Black Monday I.** Drop in DJ Industrial of 7.79% was caused by the shock from coronavirus, to the value of *s*<sup>19</sup> = 23,851, on 9 March 2020 (Black Monday I), as demonstrated in Figure 4. The data and the components defining spectrum of scenarios are presented.

Again, there are two asymmetric humps in the probabilistic space, and the region of large negative *r* dominates. The extrapolation by the most stable solution results in the following result,

$$E\_1^\*(19) = 24257.9, \ M\_1^\*(19) = 0.00629791,$$

of 1.7%. There is also less stable by order of magnitude "upward" solution, as well as four additional solutions in between with multipliers of the order of unity. Using the same methodology, we obtain Gompertz approximant, and find that it gives rather good extrapolation

$$G(19) = 23669.1, \ M\_G(19) = 0.000805813.$$

with a very small multiplier, and shows accuracy of 0.76%. There is also an upward solution, by order of magnitude less stable. Averaging the two solutions improves the estimate to the error of only 0.52%.

**Black Thursday.** Drop of 9.99%, to the level of *s*<sup>16</sup> = 21,200.6 on 12 March 2020 (Black Thursday), is also believed to be caused by the coronavirus-shock. In this case, we use the standard dataset with *N* = 15 and the third-order regression to see the typical pattern shown in Figure 5.

There is again a marked asymmetry on the graphs for the components in the probabilistic space, as the region of large negative *r* prevails. The extrapolation by the most stable solution gives

$$E\_1^\*(16) = 22\,237.1, \; M\_1^\*(16) = 0.0371606,$$

bringing the numerical error 4.89%. There is also a much less stable "upward" solution. Using the same methodology for finding the discrete spectrum, we obtain Gompertz approximant, and find that it gives rather good result

$$G(16) = 21,800.2, \ M\_G(16) = 0.00997846,$$

with a very small multiplier and an accuracy of 2.83%. There is also an additional solution, even slightly more stable, leading to a super-fast decay almost to zero. Such scenario, obviously, is absent in calculations with pure exponential approximants.

**Figure 4.** Black Monday I. Pattern in DJ Industrial index preceding 9 March 2020 Non-monotonous decay pattern reminds of a hockey stick. Fourth-order regression is shown against data points. The inverse multiplier is shown as a function of the origin *r* at *t* = *TN*, *N* = 18. The first-order approximant is shown in separate figures. Level *s*<sup>18</sup> = 25,864.8 is shown with a dot-dashed line.

**Figure 5.** Black Thursday. Pattern in DJ Industrial index preceding 12 March 2020. Monotonous decay pattern. Third-order regression is shown against data points. The inverse multiplier is shown as a function of the origin *r* at *t* = *TN*, *N* = 15. The first-order approximant is shown in separate figures. Level *s*<sup>15</sup> = 23,553.2 is shown with a dot-dashed line.

**Black Monday II.** Consider also the massive crash of 12.93%, to the value of *s*<sup>16</sup> = 20,188.5 on 16 March 2020 (Black Monday II), caused also by oil anti-shock. Because the USA is the largest producer of oil, the big drop in oil prices (anti-shock) caused an effect typically attributed to oil shock. In this case, we again use the dataset of standard length with *N* = 15, to see the typical pattern shown in Figure 6. It demonstrates the data, approximant and multiplier.

**Figure 6.** Black Monday II. Pattern in DJ Industrial index preceding 16 March 2020. Non-monotonous decay pattern. Fourth-order regression is shown against data points. The inverse multiplier is shown as a function of the origin *r* at *t* = *TN*, *N* = 15. The first-order approximant is shown in separate figures. Level *s*<sup>15</sup> = 23,185.6 is shown with a dot-dashed line.

There are two typical asymmetric humps in the probabilistic space, and the region of large negative *r* dominates. The extrapolation by the most stable solution gives the following values,

$$E\_1^\*(16) = 20,810.7, \; M\_1^\*(16) = 0.00777882,$$

bringing the numerical error of 3.08%. There is also much less stable "upward" solution,

$$E\_1^\*(16) = 27,387, \quad M\_1^\*(16) = 0.058839,$$

as well two additional solutions in between, with multipliers of the order of unity. Using the same optimization methodology, we obtain Gompertz approximant, and find extrapolations

$$G(16) = 19,987.4, \ M\_G(16) = 0.00100679,$$

with accuracy of 0.996%.

**Fear of second wave of coronavirus**.Bubble configuration corresponds to the price (index) going up monotonously, with rapid change of direction at some point, during the time scale of order of the time-series resolution. The growth finally becomes unsustainable. The crash of 11 June 2020 had started overnight. The index dropped to *s*<sup>17</sup> = 25,128.2, corresponding to a mini-crash of 6.9%. For the the dataset of length *N* = 16, we observe almost a perfect bubble, as shown in Figure 7. It demonstrates the data, approximant and multiplier as functions of origin.

**Figure 7.** Temporal bubble in Dow Jones Industrial index, preceding mini-crash of 11 June 2020. Fourth-order regression is shown against data points. The first-order approximant and multiplier are shown in separate figures. Level *s*<sup>16</sup> = 26,990 is shown with a dot-dashed line.

There is also a marked asymmetry in the probabilistic space, and the region of large negative *r* dominates. In the current case, the pattern appeared before the very day of crash and evolved into the mini-crash due to the overnight shock.

Extrapolation by the most stable solution results in

$$E\_1^\*(17) = 25,641, \ M\_1^\*(17) = 0.0124981,$$

bringing the error of 2.04%. There is also less stable "upward" solution,

$$E\_1^\*(17) = 28\,\text{814.7}, \quad M\_1^\*(17) = 0.0435021,$$

as well two additional solutions in between with multipliers of the order of unity.

Similar calculations with Gompertz approximant, give better estimate for the crash,

$$G(17) = 25,189.9, \quad M\_G(17) = 0.00169455,$$

with error of just 0.25%. One can think that fear of a second coronavirus wave leads to self-generated nervousness, leading to panic [69], having the net result of a shock. Bubbles are quite rare patterns in DJ index and more typical to Shanghai Composite [30].

### *3.4. Comments*

Many more examples of various notable crashes can be found in [30]. They were selected to exemplify market reaction to various shocks, including 9/11, Fukushima disaster, US entrance to the Great War, death of Chinese leader Deng Xiaoping, Friday the 13th, flash crash, etc. and to demonstrate similarity of early panics with coronavirus recession. Despite their different "geometry", different temporal patterns preceding crashes exhibit probabilistic distributions analogous in their main features, with significant difference only in the region of moderate *r*, but with analogous structure for large negative and positive origins. Crashes are attributed to the stable solutions with large negative *r*, meaning that the origin of time has to be moved to the deep past to explain the crash in the near future. Preliminary results of Gluzman [30] suggest that, in the overwhelming majority of cases, a crash is preceded by similar, asymmetric probability pattern(s), of the type shown in figures of this section.

Exponential and Gompertz approximants are found to work rather well, despite (or possibly due to) their simplicity. Unlike all other approximants, they give very clear graphic snapshots of the probabilistic space. Besides, their application is grounded in the exponential form of any future contract, with a transparent interpretation to the renormalized trend parameter *β*(*t*,*r*), as expected return per unit time, equivalent to inverse relaxation (growth) time.

Our theory explains or at least gives a hint why making predictions about the future is so notoriously difficult. Instead of a unique, ironclad solution to the problem, we advocate finding all solutions and interpreting them as bounds, as plainly illustrated in Figure 1. Bounds are given different strengths, a priori determined by multipliers. Reality is not completely confined to reaching the most stable bound, but various metastable bounds can be realized as well, blurring the picture and complicating emergent time dynamics.

After applying some arguments concerned with broken/restored time-invariance, we come to the exponential solution with explicit finite time scale, which was only implicit in initial parameterization with polynomial regressions. In condensed matter physics and field theory, there is a key Meissner–Higgs mechanism for generating mass or, equivalently, for creating some typical space scale from original fields through broken symmetry technique (see, e.g., [71]). Relatively recently, the concept was confirmed, culminating in the discovery of the Higgs boson. Our approach to market price evolution is by all means inspired by the Meissner–Higgs effect. However, instead of a mass of mind-boggling elementary particles, we have a mundane, but highly sought after return per unit time.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **Appendix A. Critical Index Calculations with Padé and** *DLog* **Pad'e Techniques**

For low Reynolds numbers *R*, the flow of a viscous fluid through a channel is described by the well-known Darcy's law. The Darcy law describes a linear relation between the average pressure gradient ∇*p* and the average velocity *u* along the pressure gradient [72]. It is given as follows,

$$|\overline{\nabla p}| = \frac{\eta}{K}\overline{u}\_{\prime} \tag{A1}$$

where *K* stands for the permeability and *η* is the dynamic viscosity of the fluid. The definition of permeability simply characterizes the amount of viscous fluid flow through a porous medium per unit time and unit area when a unit macroscopic pressure gradient is applied to the system [12]. The classical Poiseuille flow is a classic example, which yields the Darcy's law. It unfolds in the channel bounded by two parallel planes separated by a distance 2*b*, generated by an average pressure gradient ∇*p*. The flow profile is known to be parabolic when the Reynolds number is small. When the channel is "wavy", i.e., not straight and when the Reynolds number is not negligible, additional terms appear in this relation.

Darcy law holds in the interesting cases of the Stokes flow through a channel with two-dimensional and three-dimensional wavy walls. The enclosing wavy walls are described by the analytical expressions, including the amplitude of waviness. The amplitude is proportional to the mean clearance of the channel and is multiplied by the small dimensionless parameter *-*.

We briefly discuss below the main steps of the derivation leading to the expansions for permeability, as obtained by Mityushev, Malevich and Adler. In Ref. [35], a general asymptotic analysis was applied to a Stokes flow in curvilinear three-dimensional channel. It is bounded by walls of rather general shape described as follows

$$z = \mathcal{S}^+(\mathbf{x}\_1, \mathbf{x}\_2) \equiv b(1 + \epsilon T(\mathbf{x}\_1, \mathbf{x}\_2)),\tag{A2}$$

$$z = S^{-}(\mathbf{x}\_1, \mathbf{x}\_2) \equiv -b\left(1 + \epsilon B(\mathbf{x}\_1, \mathbf{x}\_2)\right). \tag{A3}$$

The formally small dimensionless parameter *-* ≥ 0 is considered. It is introduced in such a way to allow the general shape to be recast as the geometric perturbation around the straight channel. The expansion then is accomplished around the straight channel considered as zero-approximation. Such approach builds on an original work by Pozrikidis [73].

In [12,35], arbitrary profiles *S*±(*x*1, *x*2) were explored. It was assumed only that they satisfy some natural conditions, such as

$$|T(\mathbf{x}\_1, \mathbf{x}\_2)| \lessapprox 1 \quad \text{and} \quad |B(\mathbf{x}\_1, \mathbf{x}\_2)| \lessapprox 1. \tag{A4}$$

The infinite differentiability is assumed for the functions *T*(*x*1, *x*2) and *B*(*x*1, *x*2). Such assumption was made in order to calculate velocities and permeability, and to solve an emerging cascade of boundary value problems for the Stokes equations in a straight channel [35]. Influence of the curvilinear edges on flow is of significant theoretical interest. It illustrates the mechanism of viscous flow under different geometrical conditions.

To make our paper self-consistent, we bring below some general information about the mathematical formulation of the problem and some permeability definitions. Let **u** = **u**(*x*1, *x*2, *x*3) be the velocity vector, and *p* = *p*(*x*1, *x*2, *x*3) the pressure. The flow of a viscous fluid through a channel is considered under condition that the Reynolds number is small and the Stokes flow approximation is valid. The fluid is governed by the Stokes equations. The solution **u** of the Stokes equations is sought within the class of functions periodic with period 2*L* both in variable *x*<sup>1</sup> and in variable *x*2.

Let also *u* be the *x*-component of **u**. Let also an overall external gradient pressure ∇*p* be applied along the *x*1-direction. It corresponds to a constant jump 2*L*∇*p* along the *x*1-axis of the periodic cell. Then, the permeability of the channel in the *x*1-direction *Kx*<sup>1</sup> (*-*) is defined as the result of integration,

$$K\_{\mathbf{x}\_1}(\epsilon) = -\frac{\mu}{\overline{\nabla p}\,|\mathbf{r}\rangle} \int\_{-L}^{L} \int\_{-L}^{L} dx\_1 \, dx\_2 \int\_{-\left(\mathbf{x}\_1, \mathbf{x}\_2\right)}^{\mathbf{s}^+(\mathbf{x}\_1, \mathbf{x}\_2)} u(\mathbf{x}\_1, \mathbf{x}\_2, \mathbf{x}\_3) \, d\mathbf{x}\_3 \,. \tag{A5}$$

Here, |*τ*| stands for the volume of the unit cell *Q* of the channel. The sought *Kx*<sup>1</sup> (*-*) in (A5) is expressed explicitly as a function in *-*. More precisely, we are interested in the ratio *K* = *K*(*-*) of the dimensional permeability for the curvilinear channel and permeability of the Poiseuille flow.

Most important for our methodology, the formulae of Mityushev, Malevich and Adler [35] determine the coefficients of a Taylor series expansion for the permeability

$$K(\epsilon) = \sum\_{n=0}^{\infty} c\_n \epsilon^n,$$

with the normalization with respect to the dimensional permeability for the of the Poiseuille flow. In practical computations, *K*(*-*) is approximated by means of the truncation, leading to the Taylor polynomial of the order *k*

$$K\_k(\epsilon) = \sum\_{n=0}^k c\_n \epsilon^n. \tag{A6}$$

The domain of application of this formula appears to be restricted. The corresponding Taylor series are divergent for larger *-*.

### *Appendix A.1. Symmetric Sinusoidal Two-Dimensional Channel: Walls Can Touch*

Mityushev, Malevich and Adler [35] considered the following bounded two-dimensional channel

$$z = b(1 + \epsilon \cos x) \quad , \quad z = -b(1 + \epsilon \cos x) \,. \tag{A7}$$

The expansion for permeability was found up to *O*(*-*<sup>32</sup>), and for *b* = 0.5. This example is popular among the researchers, as is documented in [35]. The following truncated polynomial for the permeability as the function of "waviness" parameter was presented,

$$K\_{30}(\varepsilon) = \\\varepsilon^{1} - 3.14963\varepsilon^{2} + 4.08109\varepsilon^{4} - 3.48479\varepsilon^{6} + 2.93797\varepsilon^{8} - 2.56771\varepsilon^{10} + 2.21983\varepsilon^{12} - 1.93018\varepsilon^{16} + 1.67294\varepsilon^{16} - 1.45302\varepsilon^{18} + 1.26017\varepsilon^{20} - \\\ 1.09411\varepsilon^{22} + 0.949113\varepsilon^{24} - 0.823912\varepsilon^{26} + 0.714804\varepsilon^{28} - 0.620463\varepsilon^{20} \\\ + O(\varepsilon^{32}).$$

On the other hand, for larger *-*, a lubrication approximation *Kl* was discussed by Adler [72]. It is motivated by the solution in the case of two cylinders of different radii that are almost in contact with one another along a line. As *-* → *<sup>c</sup>* = 1, we arrive to the following power-law

$$K\_l \simeq \frac{8\sqrt{2}\sqrt{b^4}(\varepsilon-1)^{5/2}}{9\pi}.\tag{A9}$$

It has the general critical form, with the critical index for permeability κ = 5/2. The critical amplitude can be extracted as well, so that *A* = <sup>8</sup> <sup>√</sup>2*b*<sup>2</sup> <sup>9</sup>*<sup>π</sup>* . In the case under consideration, we calculate *A* = 0.100035.

The reasons for failure of lubrication approximation are explained in [35,72], as well as in [12]. In a nutshell, the main assumption of the lubrication approximation is that the velocity has a parabolic profile. Even for the plane channels [35], the lubrication approximation gives correct results only for channels in which the mean surface is sufficiently close to a plane and for small value of *-*.

In what follows, we completely avoid the lubrication approximation by following the approach of Gluzman [12] (Chapter 7). The technique of approximants allows approaching the critical region, when the walls nearly touch, only based on the expansion (A8).

As an input, we have the polynomial approximation (A8) of the function *K*(*-*). We intend to to calculate the critical index and amplitude(s) of the asymptotically equivalent approximants in the vicinity of the threshold *-* = *<sup>c</sup>* = 1. When such extrapolation problem is solved, one can proceed with an interpolation problem. In the latter case, assuming that the critical behavior is known in advance, one can derive the compact formula for all *-*(see Chapter 7, [12]).

Let us calculate the index and amplitude for the critical behavior written in general form

$$\mathcal{K}(\boldsymbol{\epsilon}) \simeq A(\boldsymbol{\epsilon}\_{\boldsymbol{\varepsilon}} - \boldsymbol{\epsilon})^{\varkappa}, \text{ as } \boldsymbol{\epsilon} \to \boldsymbol{\epsilon}\_{\boldsymbol{\varepsilon}} - \boldsymbol{0} \text{ , \tag{A10}$$

with unknown index and amplitude.

Let us first apply the transformation,

$$z = \frac{\epsilon}{1 - \epsilon} \iff \epsilon = \frac{z}{z + 1}\prime$$

to the series (A8). The transformation makes technical application of the different approximants more convenient.

To the transformed series *M*1(*z*), let us apply the *DLog* transformation and obtain the transformed series *M*(*z*). In terms of *M*(*z*) one can readily obtain the sequence of Padé approximations κ*<sup>n</sup>* for the critical index κ. Namely, we obtain the sequence of values

$$\varkappa\_n = -\lim\_{z \to \infty} (zPadeApproxinant[M[z], n, n+1])\,. \tag{A11}$$

as described in Section 2. The approximations for the critical index generated by the sequence of Padé approximants, converge nicely to the value 5/2, as shown below,

> κ<sup>1</sup> = 2.57972, κ<sup>2</sup> = 2.30995, κ<sup>3</sup> = 2.47451, κ<sup>4</sup> = 2.49689, κ<sup>5</sup> = 2.4959, κ<sup>6</sup> = 2.49791, κ<sup>7</sup> = 2.49923, κ<sup>8</sup> = 2.50113, κ<sup>9</sup> = 2.50028,κ<sup>10</sup> = 2.49783,κ<sup>11</sup> = 2.49778,κ<sup>12</sup> = 2.49829,κ<sup>13</sup> = 2.49836.

This result well agrees with estimates by the optimization technique of Section 2.3.

If *Bn*(*z*) = *PadeApproximant*[*M*[*z*], *n*, *n* + 1], then one can also find the approximation for permeability

$$K\_n^\*(\varepsilon) = \exp\left(\int\_0^{\frac{\varepsilon}{\varepsilon\_{\varepsilon}-\varepsilon}} B\_n(z) \, dz\right),\tag{A12}$$

and compute the corresponding amplitude

$$A\_{\mathfrak{n}} = \lim\_{\varepsilon \to \varepsilon\_{\varepsilon}} (\varepsilon\_{\mathfrak{c}} - \varepsilon)^{-\varkappa\_{\mathfrak{n}}} K\_{\mathfrak{n}}^{\*} (\varepsilon). \tag{A13}$$

The typical value of amplitude could be found as *A*<sup>9</sup> = 3.7758. It appears to be by order of magnitude larger than the value deduced from the lubrication approximation. Now, let us fix the critical index to a value of 5/2, obtained from the extrapolation procedure. Now, one can calculate *A* using the standard Padé technique, finding the value of 3.77188. The latter result turns out to be very close to the value just found above from the extrapolation.

It was illustrated by Gluzman [12] (Chapter 7) how the lubrication approximation approximation breaks down even in a close vicinity of *c*. The truncated polynomial is applicable only for small and moderately large *-*, breaking down for larger  in the vicinity of the critical point. But the final formula derived by means of factor approximant is qualitatively correct for all *-*. Obviously, the standard Padé approximants are not able to capture the non-trivial power-law in the vicinity of critical point *c*.

### *Appendix A.2. Symmetric Sinusoidal Two-Dimensional Channel: Example 2*

Let us again consider the channel bounded by the surfaces (A7), but with different parameter, *b* = 0.25. The truncated polynomial *K*(*-*) was obtained by Mityushev, Malevich and Adler [35] as well,

> *K*(*-*) = 1 − 3.03748*-*<sup>2</sup> + 3.54570*-*<sup>4</sup> <sup>−</sup> 2.33505*-*<sup>6</sup> + 1.35447*-*<sup>8</sup> <sup>−</sup> 0.83303*-*10 +0.49762*-*<sup>12</sup> <sup>−</sup> 0.30350*-*<sup>14</sup> + 0.18185*-*<sup>16</sup> <sup>−</sup> 0.11083*-*<sup>18</sup> + 0.06636*-*20 −0.04051*-*<sup>22</sup> + 0.02419*-*260.00880*-*<sup>28</sup> <sup>−</sup> 0.00544*-*<sup>30</sup>+ *O*(*-*<sup>32</sup>). (A14)

Again, as in the previous example, we follow Chapter 7 from the book [12], where the case was researched in great detail. Using Formula (A11), we found an excellent convergence in the sequence of estimates for the index,

$$\begin{aligned} \varkappa\_1 &= 2.64456, \quad \varkappa\_2 = 2.41346, \quad \varkappa\_3 = 2.49488, \quad \varkappa\_4 = 2.49992, \\\\ \varkappa\_5 &= 2.49991, \quad \varkappa\_6 = 2.50026, \quad \varkappa\_7 = 2.50068, \quad \varkappa\_8 = 2.50087, \end{aligned}$$

*Axioms* **2020**, *9*, 126

$$\varkappa\_{\theta} = 2.50086, \quad \varkappa\_{10} = 2.50063, \quad \varkappa\_{11} = 2.50063, \quad \varkappa\_{12} = 2.50086,$$

$$\varkappa\_{13} = 2.50087, \quad \varkappa\_{14} = 2.50068, \quad \varkappa\_{15} = 2.50026,$$

leading to the same value for the index as above, κ = 5/2. This result agrees with estimates by the optimization technique of Section 2.3. Clearly, the standard Padé technique fails.

The value of amplitude is estimated as well, as *A*<sup>15</sup> = 3.77362. Both amplitude and index appear to be independent on the parameter *b*, suggesting a universal regime in the vicinity of *c*.

Interpolating with the known critical index, one can calculate the amplitude *A*, using standard Padé technique, finding again the very close value of *A* ≈ 3.77316. As in the previous example, the lubrication approximation approximation breaks down even in a close vicinity of *<sup>c</sup>*. The truncated polynomial is applicable only for small and moderately large *-*, breaking down for larger  in the vicinity of the critical point. However, the final formula derived by means of factor approximant is qualitatively correct for all *-*(for more details, see Chapter 7, [12]).

The critical index, amplitude and overall behavior of permeability in the vicinity of *c*, practically do not depend on the parameter *b* [12].

### *Appendix A.3. Parallel Sinusoidal Two-Dimensional Channel. Walls Can Not Touch*

Let us proceed with the case principally different from the two cases just studied. Consider the channel bounded by the surfaces

$$z = b(1 + \epsilon \cos x) \,, \quad z = -b(1 - \epsilon \cos x) \,, \tag{A15}$$

with *b* = 0.5 [35]. It is not possible for the walls to touch, and permeability remains finite but expected to decay as a power-law as  becomes large. Instead of a critical transition from permeable to non-permeable phase, we have a non-critical transition, or crossover, as defined in [15]. The crossover is from high to low permeability and unravels with increasing parameter *-*. The crossover can still be characterized by the power-law, as one can study corresponding critical index at large *-*. Eddies are not expected in such channels even for very large *-*[35]. However, for large *b*, eddies are not excluded [35].

The truncated series expansion for the permeability were calculated up to *O*(*-*<sup>32</sup>),

$$\begin{array}{l} \mathrm{K}\_{30}(\varepsilon) = \\ \mathrm{I} - 2.53686 \times 10^{-1} \varepsilon^{2} + 4.28907 \times 10^{-2} \varepsilon^{4} - 5.46188 \times 10^{-3} \varepsilon^{6} \\ + 4.54695 \times 10^{-4} \varepsilon^{8} + 9.0656 \times 10^{-6} \varepsilon^{10} - 1.41572 \times 10^{-5} \varepsilon^{12} + 3.76584 \times 10^{-6} \varepsilon^{14} \\ - 6.72021 \times 10^{-7} \varepsilon^{16} + 7.58331 \times 10^{-8} \varepsilon^{18} + 2.34495 \times 10^{-9} \varepsilon^{20} - 4.59993 \times 10^{-9} \varepsilon^{22} \\ + 1.88446 \times 10^{-9} \varepsilon^{24} - 8.6005 \times 10^{-11} \varepsilon^{26} + 3.34156 \times 10^{-9} \varepsilon^{28} + 1.63748 \times 10^{-9} \varepsilon^{30}. \end{array} \tag{A16}$$

In this case, it is well understood that the velocity is analytic in  in the disk |*-*| < *-*0. Therefore, one can deduce that (A16) is valid for *-* < *-*0, where *-*<sup>0</sup> is of order <sup>1</sup> *<sup>b</sup><sup>χ</sup>* , with *χ* being the maximal wave number of *T*(*x*1, *x*2) and *B*(*x*1, *x*2). However, to extend *K*(*-*) for *- -*0, it was suggested to apply the Padé approximation to the polynomial (A16), which agrees with it up to *O*(*-*<sup>32</sup>).

The Padé approximant of the order (10, 20), denoted here as *K*10,20(*-*), was first developed by Mityushev, Malevich and Adler [35]. Its explicit expression can also be found in Chapter 7 of the book [12]. This approximant gives *K*10,20(*-*) ∼ *-*<sup>−</sup>10, as *-* → ∞. One can think then that the permeability decays as

$$\mathbb{K}(\epsilon) \simeq B\epsilon^{\nu},$$

as *-* → ∞, with the critical index *ν* different from the estimate given by *K*10,20(*-*). Calculation of the critical index *ν* was accomplished in Chapter 7 of the book [12].

Assuming that the small-variable expansion for the function is given by the truncated sum *K*30(*-*) (A16), we can find the corresponding small-variable expression for the effective critical exponent which

equals  *d d* log *K*30(*-*). By applying to the obtained series, the method of Padé approximants, as in two previous examples, the sought approximate expression for the critical exponent

$$\nu\_k = \lim\_{\epsilon \to \infty} \epsilon P\_{k,k+1}(\epsilon) \, , \tag{A17}$$

can be computed dependent on the approximation order *k*. Application of the method to the truncated power series (A16), is straightforward and suggests strongly the value of *ν* = −4, as can be seen in Figure A1. This result agrees with estimates by the optimization technique of Section 2.3. Clearly, the Padé estimate mentioned above fails. The amplitude *B*, corresponding to *k* = 14, is equal to 44.5872.

**Figure A1.** The index *ν* at infinity, is shown dependent on the approximation number *k*. The values found by computing (A17) are shown with black circles. They are compared with the most plausible value of −4 (shown with gray circles).

Assume now that *ν* = −4 and construct the sequence of Padé approximants *Pn*,*n*+<sup>4</sup> for the original truncated polynomial (A16). There is a convergence in the approximation sequence for the amplitude *B*. One can safely assume that it converges to the value of 43.2. The sequence is shown in Figure A2.

**Figure A2.** The amplitude *B* dependence on approximation number *k* is shown with black circles. One can see the convergence to the value of 43.2, shown with squares.

### **Appendix B. Example of Interpolation with Root Approximants: One-Dimensional Bose Gas**

Lieb and Liniger [74] considered a one-dimensional Bose gas with contact interactions. The ground-state energy of the gas can be written as a weak-coupling expansion, with respect to the coupling parameter *g* [75,76], as

$$E(\mathbf{g}) \simeq \mathbf{g} - \frac{4}{3\pi} \mathbf{g}^{3/2} + \frac{1.29}{2\pi^2} \mathbf{g}^2 - 0.017201 \mathbf{g}^{5/2} \,, \tag{A18}$$

as *g* → 0, In the strong-coupling limit, as *g* → ∞, we have the following expression [75,76]

$$E(\mathcal{g}) \simeq \frac{\pi^2}{3} (1 - \frac{4}{\mathcal{g}} + \frac{12}{\mathcal{g}^2}) \,. \tag{A19}$$

In what follows, *E*∗ <sup>3</sup>+3(*g*) assimilates the three coefficients from weak and strong coupling expansions, while *E*∗ <sup>4</sup>+3(*g*) is based on all four terms from the weak-coupling side.

The accuracy of the root approximants (17)

$$E\_{3+3}^{\*}(\mathbf{g}) = \frac{\pi^{2}}{3\sqrt[5]{\frac{385.383}{x^{5}} + \left(\frac{388.171}{x^{4}} + \left(\frac{164.914}{x^{3}} + \left(\frac{37.3454}{x^{2}} + \left(\frac{8.13698}{x} + 1\right)^{3/2}\right)^{5/4}\right)^{7/6}\right)^{9/8}} \tag{A20}$$

and

$$E\_{4+3}^{\*}(g) = \frac{\pi^{2}}{3\sqrt[6]{\frac{12\delta T\_{2}\delta g}{x^{6}} + \left(\frac{1546.85}{x^{3}} + \left(\frac{811.495}{x^{4}} + \left(\frac{254.69}{x^{3}} + \left(\frac{85.654}{x^{2}} + 1\right)^{3/2}\right)^{5/4}\right)^{7/6}\right)^{9/8}} \tag{A21}$$

turns out to be good. The approximants are constructed from "right-to-left". i.e., we self-similarly connect a known asymptotic expansion at the right boundary of the interval with a known asymptotic form at the left boundary.

In Table A1, they are compared to the extensive numerical data obtained by Dunjko and Olshanii *EDO* [77]. The Padé-estimates, *EP*, are also presented. The Padé approximant *P*3,5( √*g*) reads as follows:

$$P\_{3,5}(\sqrt{\mathbf{g}}) = \frac{\mathbf{g}\left(0.285957\mathbf{g}^{3/2} - 0.177533\mathbf{g} + 0.355474\sqrt{\mathbf{g}} + 1\right)}{0.458734\mathbf{g}^{3/2} + 0.0869206\mathbf{g}^{5/2} - 0.0599636\mathbf{g}^2 + 0.0881093\mathbf{g} + 0.77988\mathbf{Z}\sqrt{\mathbf{g}} + 1} \tag{A.22}$$

**Table A1.** Ground-state energy of Lieb-Liniger model, for the varying dimensionless parameter *g*, in different approximations: Root approximants *E*∗ <sup>3</sup>+3(*g*), *E*<sup>∗</sup> <sup>4</sup>+3(*g*), numerical data *EDO*, and the Padé approximant *EP*.


It should become completely clear from observing Figure A3 that the problem of interpolation is neither simple nor superficial. The asymptotic expressions for small and large couplings have little in common with each other. Although the expansions (A18) and (A19) appear to work only for very small and very large coupling constants, the deduced approximants work rather well. More examples of interpolation with various self-similar approximants can be found in [16].

**Figure A3.** The interpolation with root approximant (A20) is shown with solid line, while the Padé approximant is shown with dotted line. The weak- (dashed) and strong-coupling (dot-dashed) expansions are shown as well.

### **References**


[CrossRef]


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
