The Tensor Pade´-Type Approximant with Application in Computing Tensor Exponential Function

Abstract

Tensor exponential function is an important function that is widely used. In this paper, tensor Pad$\acute{e}$-type approximant (TPTA) is defined by introducing a generalized linear functional for the first time. The expression of TPTA is provided with the generating function form. Moreover, by means of formal orthogonal polynomials, we propose an efficient algorithm for computing TPTA. As an application, the TPTA for computing the tensor exponential function is presented. Numerical examples are given to demonstrate the efficiency of the proposed algorithm.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Tensor exponential function is an important function that is widely used, owing to its key role in the solution of tensor differential equations [1–4]. For instance, Markovian master equation can be written as tensor differential equations $\left(\partial /\partial t\right)\mathbf{P}(t)=\mathbf{A}·\mathbf{P}(t)$, where the probabilities tensor $\mathbf{P}(t)\in {\mathbb{R}}^{n_{\mathrm{1}}\times n_{\mathrm{2}}\times \cdots \times n_d}$ [5]. Consider the initial value problem defined by the tensor ordinary differential equation [6, 7]$$\begin{array}{c}\text{(1)}& \dot{\mathcal{Y}}\left(t\right)=\mathcal{A}\mathcal{Y}\left(t\right),\mathcal{Y}\left(t_{\mathrm{0}}\right)={\mathcal{Y}}_{\mathrm{0}},\end{array}$$where the superimposed dot denotes differentiation with respect to $t$ and $\mathcal{A}$ and ${\mathcal{Y}}_{\mathrm{0}}$ are given constant tensors. The solution to system (1) is $\mathcal{Y}(t)=\mathrm{e}\mathrm{x}\mathrm{p}[(t-t_{\mathrm{0}})\mathcal{A}]{\mathcal{Y}}_{\mathrm{0}}$, and $\mathrm{e}\mathrm{x}\mathrm{p}(·)$ is the tensor exponential function.

To solve (1), we need to calculate an exponential function about tensor $\mathcal{A}$. Recently, tensor computation, especially eigenvalues of tensors, has attracted attention of many scholars; some important results have been found in the current literature [8–15]. For instance, computation of $\mathrm{e}\mathrm{x}\mathrm{p}(\mathcal{A}t)$ is required in applications such as finite strain hyperelastic-based multiplicative plasticity models [7, 16–19]. Explicitly, for a generic tensor $\mathcal{A}$, the tensor exponential can be expressed by means of its series representation [7] (see p.749).$$\begin{array}{}\text{(2)}& \mathrm{e}\mathrm{x}\mathrm{p}\left(\mathcal{A}\right)={\displaystyle \sum }_{n=\mathrm{0}}^{\infty }\frac{\mathrm{1}}{n!}{\mathcal{A}}^n,\end{array}$$and the preceding series is absolutely convergent for any argument $\mathcal{A}$ and, as its scalar counterpart, can be used to evaluate the tensor exponential function to any prescribed degree of accuracy [16]. The computation of (2) is carried out by simply truncating the infinite series with $n_{max}$ terms: $$\begin{array}{}\text{(3)}& \mathrm{e}\mathrm{x}\mathrm{p}\left(\mathcal{A}\right)={\displaystyle \sum }_{n=\mathrm{0}}^{n_{max}}\frac{\mathrm{1}}{n!}{\mathcal{A}}^n,\end{array}$$ with $n_{max}$ being such that $\left(\mathrm{1}/n_{max}\right)‖{\mathcal{A}}^{n_{max}}‖\le {ϵ}_{tol}$.

However, the accuracy and effectiveness of the preceding algorithm is limited by round-off and choice of termination criterion [16]. Pad$\acute{e}$ approximant has become by far the most widely used one in calculation of exponential function or formal power series due to the following reasons: first, the series may converge too slowly to be of any use and the approximation can accelerate its convergence; second, only few coefficients of the series may be known and a good approximation to the series is needed to obtain properties of the function that it represents [20]. For instance, matrix Pad$\acute{e}$-type approximant (MPTA) [21] can be used to simplify the high degree multivariable system by approximating transfer function matrix $G(s)$ that can be expanded into a power series with matrix coefficients, i.e., $G(s)={\sum }_{i=\mathrm{0}}^{\infty }{c}_i{s}^i$, where $c_i\in {\mathbb{C}}^{s\times t}$. The key to construct TPTA is to maintain the same order of tensor $\mathcal{A}$ with different powers. For this issue, we introduce t-product [22, 23] of two tensors to solve it. In addition, in order to give the definition of TPTA, we introduce a generalized linear functional in the tensor space for the first time.

This paper is organized as follows. In Section 2, we provide some preliminaries. First, we introduce the t-product of two tensors; then, we show the definitions of tensor exponential function and the Frobenius norm of a tensor. In Section 3.1, we define the tensor Pad$\acute{e}$-type approximant by using generalized linear functional; the expression of TPTA is of the form of tensor numerator and scalar denominator; and then we introduce the definition of orthogonal polynomial with respect to generalized linear functional and sketch an algorithm to compute the TPTA. Numerical examples are given and analyzed in Section 4. Finally, we finish the paper with concluding remarks in Section 5.

2. Preliminaries

There arise mainly a problem for approximating tensor exponential function. That is how to expand $e^{\mathcal{A}t}$ into the power series for order-$p$ ($p\ge \mathrm{3}$) tensors. For a Symmetric and second-order tensor $A$, higher powers of $A$ can be computed by the Cayley-Hamilton theorem [24], but it fails for the order-$p$ ($p\ge \mathrm{3}$) tensors. Therefore, we shall utilize the t-product to obtain higher powers of order-$p$ ($p\ge \mathrm{3}$) tensors in this section. Firstly, we introduce some notations and basic definitions which will be used in the sequel. Throughout this paper tensors are denoted by calligraphic letters (e.g., $\mathcal{A}$, $\mathcal{B}$), while capital letters represent matrices, and lowercase letters refer to scalars.

An order-$p$ tensor, $\mathcal{A}$, can be written as $$\begin{array}{}\text{(4)}& \mathcal{A}=\left(a_{i_{\mathrm{1}}{i}_{\mathrm{2}}\cdots i_p}\right)\in {\mathbb{R}}^{n_{\mathrm{1}}\times n_{\mathrm{2}}\times \cdots \times n_p}.\end{array}$$

Thus, a matrix is considered a second-order tensor, and a vector is a first-order tensor [22], for $i=\mathrm{1},\dots ,n_p$, denoted by ${\mathcal{A}}_i\in {\mathbb{R}}^{n_{\mathrm{1}}\times n_{\mathrm{2}}\times \cdots \times n_{p-\mathrm{1}}}$, the tensor whose order is $p-\mathrm{1}$ and is created by holding the $p$th index of $\mathcal{A}$ fixed at $i$. For example, consider a third-order tensor, $\mathcal{A}=(a_{i_{\mathrm{1}}{i}_{\mathrm{2}}{i}_{\mathrm{3}}})\in {\mathbb{R}}^{\mathrm{3}\times \mathrm{3}\times \mathrm{3}}.$ Fixing the $\mathrm{3}$rd index of $\mathcal{A}$, we can get three $\mathrm{3}\times \mathrm{3}$ matrices, namely, 2-order tensor, which are $A_{\mathrm{1}}$, $A_{\mathrm{2}}$, and $A_{\mathrm{3}}$ and with elements $$\begin{array}{c}\text{(5)}& A_{\mathrm{1}}:\begin{array}{ccccccccc}{a}_{\mathrm{111}}& a_{\mathrm{121}}& a_{\mathrm{131}}& a_{\mathrm{211}}& a_{\mathrm{221}}& a_{\mathrm{231}}& a_{\mathrm{311}}& a_{\mathrm{321}}& a_{\mathrm{331}}\end{array},A_{\mathrm{2}}:\begin{array}{ccccccccc}{a}_{\mathrm{112}}& a_{\mathrm{122}}& a_{\mathrm{132}}& a_{\mathrm{212}}& a_{\mathrm{222}}& a_{\mathrm{232}}& a_{\mathrm{312}}& a_{\mathrm{322}}& a_{\mathrm{332}}\end{array},\\ A_{\mathrm{3}}:\begin{array}{ccccccccc}{a}_{\mathrm{113}}& a_{\mathrm{123}}& a_{\mathrm{133}}& a_{\mathrm{213}}& a_{\mathrm{223}}& a_{\mathrm{233}}& a_{\mathrm{313}}& a_{\mathrm{323}}& a_{\mathrm{333}}\end{array},\end{array}$$ respectively.

Now, we will define the t-product of two tensors.

Definition 1 (see [22, 23]).

Let $\mathcal{A}\in {\mathbb{R}}^{n_{\mathrm{1}}\times n_{\mathrm{2}}\times \cdots \times n_p}$. Then the block circulant pattern tensor of $\mathcal{A}$ is denoted by $$\begin{array}{}\text{(6)}& circ\left(\mathcal{A}\right)=\left(\begin{array}{ccccc}{\mathcal{A}}_{\mathrm{1}}& {\mathcal{A}}_{n_p}& {\mathcal{A}}_{n_p-\mathrm{1}}& \cdots & {\mathcal{A}}_{\mathrm{2}}\\ {\mathcal{A}}_{\mathrm{2}}& {\mathcal{A}}_{\mathrm{1}}& {\mathcal{A}}_{n_p}& \cdots & {\mathcal{A}}_{\mathrm{3}}\\ ⋮& \ddots & \ddots & \ddots & ⋮\\ {\mathcal{A}}_{n_p}& {\mathcal{A}}_{n_p-\mathrm{1}}& \cdots & {\mathcal{A}}_{\mathrm{2}}& {\mathcal{A}}_{\mathrm{1}}\end{array}\right),\end{array}$$ where ${\mathcal{A}}_i$$\in {\mathbb{R}}^{n_{\mathrm{1}}\times n_{\mathrm{2}}\times \cdots \times n_{p-\mathrm{1}}}$, $i=\mathrm{1,2},\dots ,n_p.$

Define unfold$(·)$ to an $n_{\mathrm{1}}\times n_{\mathrm{2}}\times \cdots \times n_p$ tensor by an $n_{\mathrm{1}}{n}_p\times n_{\mathrm{2}}\times \cdots \times n_{p-\mathrm{1}}$ block tensor in the following way: $$\begin{array}{}\text{(7)}& unfold\left(\mathcal{A}\right)=\left(\begin{array}{c}{\mathcal{A}}_{\mathrm{1}}\\ {\mathcal{A}}_{\mathrm{2}}\\ ⋮\\ {\mathcal{A}}_{n_p}\end{array}\right).\end{array}$$ If $\mathcal{A}$ is order-3 tensor, then $unfold(\mathcal{A})$ is a block vector. Similarly, define fold$(·)$ as the inverse operation, which takes an $n_{\mathrm{1}}{n}_p\times n_{\mathrm{2}}\times \cdots \times n_{p-\mathrm{1}}$ block tensor and returns an $n_{\mathrm{1}}\times n_{\mathrm{2}}\times \cdots \times n_p$ tensor; then $$\begin{array}{}\text{(8)}& fold\left(unfold\left(\mathcal{A}\right)\right)=\mathcal{A}.\end{array}$$

Definition 2 (see [23]).

Let $\mathcal{A}$ be ${\mathbb{R}}^{n_{\mathrm{1}}\times n_{\mathrm{2}}\times \cdots \times n_p}$ and $\mathcal{B}$ be ${\mathbb{R}}^{n_{\mathrm{2}}\times l\times n_{\mathrm{3}}\cdots \times n_p}$. Then the t-product $\mathcal{A}*\mathcal{B}$ is the $n_{\mathrm{1}}\times l\times n_{\mathrm{3}}\times \cdots \times n_p$ tensor defined recursively as $$\begin{array}{}\text{(9)}& \mathcal{A}*\mathcal{B}=fold\left(circ\left(\mathcal{A}\right)*unfold\left(\mathcal{B}\right)\right).\end{array}$$

Remark 3.

If $\mathcal{A}$ and $\mathcal{B}$ are order-2 tensors, then the product “$\ast$” can be replaced by standard matrix multiplication.

Remark 4.

The $k$ times power of $\mathcal{A}$ is defined as ${\mathcal{A}}^k=\mathcal{A}*\mathcal{A}*\cdots *\mathcal{A}$ ($k$ times); “$\ast$” denotes the t-product.

Example 5.

Letting$$\begin{array}{}\text{(10)}& \mathcal{A}=\left(\begin{array}{cccccc}\mathrm{1}& \mathrm{0}& -\mathrm{1}& \mathrm{3}& \mathrm{1}& \mathrm{0}\\ \mathrm{2}& \mathrm{1}& \mathrm{0}& \mathrm{2}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -\mathrm{1}& \mathrm{1}& \mathrm{1}& -\mathrm{1}& \mathrm{2}\end{array}\right)\in {\mathbb{R}}^{\mathrm{3}\times \mathrm{3}\times \mathrm{2}},\end{array}$$then, from Definition 2, we have $$\begin{array}{}\text{(11)}& {\mathcal{A}}^{\mathrm{4}}=\mathcal{A}*\mathcal{A}*\mathcal{A}*\mathcal{A}=\left(\begin{array}{cccccc}\mathrm{271}& \mathrm{113}& -\mathrm{98}& \mathrm{237}& \mathrm{104}& -\mathrm{119}\\ \mathrm{242}& \mathrm{97}& -\mathrm{96}& \mathrm{242}& \mathrm{96}& -\mathrm{96}\\ -\mathrm{73}& -\mathrm{46}& \mathrm{56}& -\mathrm{94}& -\mathrm{52}& \mathrm{43}\end{array}\right).\end{array}$$

Remark 6.

One of the characteristic features of t-product is that it ensures that the order of multiplication result of two tensors does not change, whereas other tensor multiplications do not have the feature; that is why we chose the t-product as the multiplication of tensors.

The tensor exponential function is a tensor function on tensors analogous to the ordinary exponential function, which can be defined as follows.

Definition 7.

Let $\mathcal{A}$ be an $n_{\mathrm{1}}\times n_{\mathrm{2}}\times \cdots \times n_p$ real or complex tensor. The tensor exponential function of $t$, denoted by $e^{\mathcal{A}t}$ or exp$(\mathcal{A}t)$, is the $n_{\mathrm{1}}\times n_{\mathrm{2}}\times \cdots \times n_p$ tensor given by the power series $$\begin{array}{}\text{(12)}& e^{\mathcal{A}t}={\displaystyle \sum }_{k=\mathrm{0}}^{\infty }\frac{\mathrm{1}}{k!}{\left(\mathcal{A}t\right)}^k,\end{array}$$ where ${\mathcal{A}}^{\mathrm{0}}$ is defined to be the identity tensor $\mathcal{I}$ (see Definition 8) with the same orders as $\mathcal{A}$.

Definition 8 (see [23]).

The $n\times n\times l_{\mathrm{1}}\times \cdots \times l_{p-\mathrm{2}}$ order-$p$ identity tensor ($p>\mathrm{3}$) $\mathcal{I}$ is the tensor such that ${\mathcal{I}}_{\mathrm{1}}$ is the $n\times n\times l_{\mathrm{1}}\times \cdots \times l_{p-\mathrm{3}}$ order-$(p-\mathrm{1})$ identity tensor, and ${\mathcal{I}}_j$ is the order-$(p-\mathrm{1})$ zero tensor, for $j=\mathrm{2,3},\dots ,l_{p-\mathrm{2}}$.

By Definition 8, we can define tensor inverse, transpose, and orthogonality. However, we do not discuss these works here, as it is beyond the scope of the present work. For the details of these definitions of tensor, we refer to reader to [22, 23, 25] and the references therein.

Let $\mathcal{A}$$\in {\mathbb{R}}^{n_{\mathrm{1}}\times n_{\mathrm{2}}\times \cdots \times n_p}$; then the norm of a tensor is the square root of the sum of the squares of all its entries [25]; i.e.,$$\begin{array}{}\text{(13)}& ‖\mathcal{A}‖=\sqrt{{\displaystyle \sum }_{i_{\mathrm{1}}=\mathrm{1}}^{n_{\mathrm{1}}}{\displaystyle \sum }_{i_{\mathrm{2}}=\mathrm{1}}^{n_{\mathrm{2}}}\hspace{0.17em}\cdots {\displaystyle \sum }_{i_p=\mathrm{1}}^{n_p}{a}_{i_{\mathrm{1}}{i}_{\mathrm{2}}\cdots i_p}^{\mathrm{2}}}.\end{array}$$

This is analogous to the matrix Frobenius norm. The inner product of two same-sized tensors $\mathcal{A}$, $\mathcal{B}\in {\mathbb{R}}^{n_{\mathrm{1}}\times n_{\mathrm{2}}\times \cdots \times n_p}$ is the sum of the products of their elements [25]; i.e.,$$\begin{array}{}\text{(14)}& \left(\mathcal{A},\mathcal{B}\right)={\displaystyle \sum }_{i_{\mathrm{1}}=\mathrm{1}}^{n_{\mathrm{1}}}{\displaystyle \sum }_{i_{\mathrm{2}}=\mathrm{1}}^{n_{\mathrm{2}}}\hspace{0.17em}\cdots {\displaystyle \sum }_{i_p=\mathrm{1}}^{n_p}{a}_{i_{\mathrm{1}}{i}_{\mathrm{2}}\cdots i_p}{b}_{i_{\mathrm{1}}{i}_{\mathrm{2}}\cdots i_p}.\end{array}$$It follows immediately that $(\mathcal{A},\mathcal{A})={‖\mathcal{A}‖}^2.$

3. Tensor Pad$\acute{e}$-Type Approximant

Let $f(x)$ be a given power series with tensor coefficients; i.e.,$$\begin{array}{}\text{(15)}& f\left(x\right)={\mathcal{A}}_{\mathrm{0}}+{\mathcal{A}}_{\mathrm{1}}x+{\mathcal{A}}_{\mathrm{2}}{x}^{\mathrm{2}}+\cdots +{\mathcal{A}}_n{x}^n+\cdots ,\hspace{1em}{\mathcal{A}}_i\in {\mathbb{R}}^{n_{\mathrm{1}}\times n_{\mathrm{2}}\times \cdots \times n_p},\hspace{0.17em}\hspace{0.17em}x\in \mathbb{R}.\end{array}$$

Let $\mathbf{P}$ denote the set of scalar polynomials in one real variable whose coefficients belong to the real field $\mathbb{R}$ and ${\mathbf{P}}_k$ denote the set of elements of $\mathbf{P}$ of degree less than or equal to $k$.

Let $\varphi :$$\mathbf{P}$$\to$${\mathbb{R}}^{n_{\mathrm{1}}\times n_{\mathrm{2}}\times \cdots \times n_p}$ be a linear functional on $\mathbf{P}$. Let it act on $t$ by$$\begin{array}{}\text{(16)}& \varphi \left(t^i\right)={\mathcal{A}}_i,\hspace{1em}i=\mathrm{0,1},\mathrm{2},\dots .\end{array}$$Then, by the linearity of $\varphi$, we have $$\begin{array}{}\text{(17)}& \varphi \left({\left(\mathrm{1}-tx\right)}^{-\mathrm{1}}\right)=\varphi \left(\mathrm{1}+tx+{\left(tx\right)}^{\mathrm{2}}+\cdots +{\left(tx\right)}^n+\cdots \right)={\mathcal{A}}_{\mathrm{0}}+{\mathcal{A}}_{\mathrm{1}}x+{\mathcal{A}}_{\mathrm{2}}{x}^{\mathrm{2}}+\cdots +{\mathcal{A}}_n{x}^n+\cdots =f\left(x\right).\end{array}$$

3.1. Definition of Tensor Pad$\acute{e}$-Type Approximant

Let $v_n\left(x\right)=b_{\mathrm{0}}+b_{\mathrm{1}}x+\cdots +b_n{x}^n\hspace{0.17em}\hspace{0.17em}(b_n\ne \mathrm{0})$ be a scalar polynomial of ${\mathbf{P}}_n$ of exact degree $n$. In this case, $v_n(x)$ is said to be quasi-monic. Define the tensor polynomial ${\mathcal{W}}_{n-\mathrm{1}}(x)$ associated with $v_n(x)$ with tensor-valued coefficients, by$$\begin{array}{}\text{(18)}& {\mathcal{W}}_{n-\mathrm{1}}\left(x\right)=\varphi \left(\frac{v_n\left(x\right)-v_n\left(t\right)}{x-t}\right).\end{array}$$It is easily seen that ${\mathcal{W}}_{n-\mathrm{1}}(x)$ is a tensor polynomial of exact degree $n-\mathrm{1}$ in $x$. Set$$\begin{array}{}\text{(19)}& {\tilde{v}}_n\left(x\right)=x^n{v}_n\left(x^{-\mathrm{1}}\right),\text{(20)}& {\tilde{\mathcal{W}}}_{n-\mathrm{1}}\left(x\right)=x^{n-\mathrm{1}}{\mathcal{W}}_{n-\mathrm{1}}\left(x^{-\mathrm{1}}\right).\end{array}$$Then, the polynomials ${\tilde{v}}_n(x)$ and ${\tilde{\mathcal{W}}}_{n-\mathrm{1}}(x)$ are obtained from $v_n(x)$ and ${\mathcal{W}}_{n-\mathrm{1}}(x)$, respectively, by reversing the numbering of the coefficients. By the procedure given above, the following conclusion is obtained.

Theorem 9.

Let ${\tilde{v}}_n(\mathrm{0})\ne \mathrm{0}$; then ${\tilde{\mathcal{W}}}_{n-\mathrm{1}}(x)/{\tilde{v}}_n(x)-f(x)=O(x^n).$

Proof.

Expanding $\left(v_n(x)-v_n(t)\right)/\left(x-t\right)$ in (18) and applying $\varphi$ yields that $$\begin{array}{}\text{(21)}& {\tilde{\mathcal{W}}}_{n-\mathrm{1}}\left(x\right)={\displaystyle \sum }_{l=\mathrm{0}}^{n-\mathrm{1}}\left({\displaystyle \sum }_{i=\mathrm{0}}^{n-l-\mathrm{1}}{b}_{l+i+\mathrm{1}}{\mathcal{A}}_i\right)x^l={\displaystyle \sum }_{l=\mathrm{0}}^{n-\mathrm{1}}\left({\displaystyle \sum }_{i=\mathrm{0}}^l{b}_{n-l+i}{\mathcal{A}}_i\right)x^l.\end{array}$$ Computing ${\tilde{v}}_n(x)f(x)$, we get $$\begin{array}{}\text{(22)}& {\tilde{v}}_n\left(x\right)f\left(x\right)=\left({\displaystyle \sum }_{j=\mathrm{0}}^n{b}_{n-j}{x}^j\right)\left({\displaystyle \sum }_{i=\mathrm{0}}^{\infty }{\mathcal{A}}_i{x}^i\right)={\displaystyle \sum }_{l=\mathrm{0}}^{\infty }\left({\displaystyle \sum }_{i=\mathrm{0}}^l{b}_{n-l+i}{\mathcal{A}}_i\right)x^l.\end{array}$$ Thus, $$\begin{array}{}\text{(23)}& {\tilde{v}}_n\left(x\right)f\left(x\right)-{\tilde{\mathcal{W}}}_{n-\mathrm{1}}\left(x\right)={\displaystyle \sum }_{l=n}^{\infty }\left({\displaystyle \sum }_{i=\mathrm{0}}^l{b}_{n-l+i}{\mathcal{A}}_i\right)x^l=O\left(x^n\right).\end{array}$$

Definition 10.

$R_{n-\mathrm{1},n}(x)={\tilde{\mathcal{W}}}_{n-\mathrm{1}}(x)/{\tilde{v}}_n(x)$ is called a tensor Pad$\acute{e}$-type approximant with order $n$ for the given power series (15) and is denoted by $(n-\mathrm{1}/n{)}_f(x).$

Remark 11.

The polynomial $v_n(x)$, called the generating polynomial of $(n-\mathrm{1}/n{)}_f$ with respect to power series $f(x)$, can be arbitrarily chosen.

Remark 12.

The tensor Pad$\acute{e}$-type approximant ${(n-\mathrm{1}/n)}_f$ possesses the degree constraint, which is caused by its construction process. The constraint implies that the method does not construct tensor Pad$\acute{e}$-type approximant of type $(m/n)$ when $m$ is different from $n-\mathrm{1}$.

To fill this gap, we define a new tensor Pad$\acute{e}$-type approximant by introducing a generalized linear functional.

Let ${\varphi }^{(q)}:$$\mathbf{P}$$\to$${\mathbb{R}}^{n_{\mathrm{1}}\times n_{\mathrm{2}}\times \cdots \times n_p}$ be a generalized linear functional on $\mathbf{P}$. Let it act on $t$ by$$\begin{array}{}\text{(24)}& {\varphi }^{\left(q\right)}\left(t^i\right)={\mathcal{A}}_{q+i},\hspace{1em}i=\mathrm{0,1},\mathrm{2},\dots .\end{array}$$Similarly to what was done for ${\mathcal{W}}_{n-\mathrm{1}}(x)$, we consider the polynomial ${\mathcal{W}}_q(x)$ associated with $v_n(x)$, and defined by$$\begin{array}{}\text{(25)}& {\mathcal{W}}_q\left(x\right)={\varphi }^{\left(q\right)}\left(\frac{v_n\left(x\right)-v_n\left(t\right)}{x-t}\right),\hspace{1em}q=m-n+\mathrm{1}.\end{array}$$Set$$\begin{array}{}\text{(26)}& {\tilde{\mathcal{W}}}_q\left(x\right)=x^{n-\mathrm{1}}{\mathcal{W}}_q\left(x^{-\mathrm{1}}\right),\end{array}$$and define$$\begin{array}{}\text{(27)}& {\mathcal{P}}_{mn}\left(x\right)={\tilde{v}}_n\left(x\right){\displaystyle \sum }_{i=\mathrm{0}}^{m-n}{\mathcal{A}}_i{x}^i+x^{m-n+\mathrm{1}}{\tilde{\mathcal{W}}}_q\left(x\right),\hspace{1em}m\ge n.\end{array}$$Then we have the following conclusion.

Theorem 13.

Let ${\tilde{v}}_n(\mathrm{0})\ne \mathrm{0}$; then ${\mathcal{P}}_{mn}(x)/{\tilde{v}}_n(x)-f(x)=O(x^{m+\mathrm{1}}).$

Proof.

Let $f_{m-n+\mathrm{1}}$ be the formal power series $$\begin{array}{}\text{(28)}& f_{m-n+\mathrm{1}}\left(x\right)={\displaystyle \sum }_{j=\mathrm{0}}^{\infty }{\mathcal{A}}_{m-n+\mathrm{1}+j}{x}^j.\end{array}$$ Then$$\begin{array}{}\text{(29)}& x^{m-n+\mathrm{1}}{f}_{m-n+\mathrm{1}}\left(x\right)=f\left(x\right)-{\displaystyle \sum }_{i=\mathrm{0}}^{m-n}{\mathcal{A}}_i{x}^i.\end{array}$$Expanding (25) and using (26) we obtain $$\begin{array}{}\text{(30)}& {\tilde{\mathcal{W}}}_q\left(x\right)={\displaystyle \sum }_{q=\mathrm{0}}^{n-\mathrm{1}}\left({\displaystyle \sum }_{i=\mathrm{0}}^q{b}_{n-q+i}{\mathcal{A}}_{m-n+\mathrm{1}+i}\right)x^q.\end{array}$$ Computing the product ${\tilde{v}}_n(x)f_{m-n+\mathrm{1}}(x)$, we find that $$\begin{array}{}\text{(31)}& {\tilde{v}}_n\left(x\right)f_{m-n+\mathrm{1}}\left(x\right)={\displaystyle \sum }_{q=\mathrm{0}}^{\infty }\left({\displaystyle \sum }_{i=\mathrm{0}}^q{b}_{n-q+i}{\mathcal{A}}_{m-n+\mathrm{1}+i}\right)x^q.\end{array}$$ By Theorem 9, one has $$\begin{array}{}\text{(32)}& \frac{{\tilde{\mathcal{W}}}_q\left(x\right)}{{\tilde{v}}_n\left(x\right)}={\left(n-\frac{\mathrm{1}}{n}\right)}_{f_{m-n+\mathrm{1}}}\left(x\right).\end{array}$$ Then, for $m\ge n$ we deduce from (27) and (29) that $$\begin{array}{}\text{(33)}& {\tilde{v}}_n\left(x\right)f\left(x\right)-{\mathcal{P}}_{mn}\left(x\right)={\tilde{v}}_n\left(x\right)\left({\displaystyle \sum }_{i=\mathrm{0}}^{m-n}{\mathcal{A}}_i{x}^i+x^{m-n+\mathrm{1}}{f}_{m-n+\mathrm{1}}\left(x\right)\right)-\tilde{v}\left(x\right){\displaystyle \sum }_{i=\mathrm{0}}^{m-n}{\mathcal{A}}_i{x}_i+x^{m-n+\mathrm{1}}{\tilde{\mathcal{W}}}_q\left(x\right)=x^{m-n+\mathrm{1}}\left({\tilde{v}}_n\left(z\right)f_{m-n+\mathrm{1}}-{\tilde{\mathcal{W}}}_q\left(x\right)\right)=x^{m-n+\mathrm{1}}\left\{{\displaystyle \sum }_{q=n}^{\infty }\left({\displaystyle \sum }_{i=\mathrm{0}}^q{b}_{n-q+i}{\mathcal{A}}_{m-n+\mathrm{1}+i}\right)x^q\right\}=O\left(x^{m+\mathrm{1}}\right).\end{array}$$

Now, we can achieve $(m/n{)}_f(x)$ by the above procedure, and it will be denoted by ${\mathcal{P}}_{mn}(x)/{\tilde{v}}_n(x)$.

Definition 14.

$R_{m,n}(x)={\mathcal{P}}_{mn}(x)/{\tilde{v}}_n(x)$ is called TPTA with order $m+\mathrm{1}$ and is denoted by $(m/n{)}_f(x).$

Algorithm 15 (compute $R_{m,n}(x)={\mathcal{P}}_{mn}(x)/{\tilde{v}}_n(x)$ with $v_n(x)$ being arbitrarily chosen).

(1)

Set $q=m-n+\mathrm{1}$ and chose a quasi-monic polynomial $v_n(x)$.

(2)

Use (19) to compute ${\tilde{v}}_n(x)$.

(3)

Compute ${\mathcal{W}}_q(x)$ and ${\tilde{\mathcal{W}}}_q(x)$ by (25) and (26), respectively.

(4)

Substitute ${\tilde{v}}_n(x)$ and ${\tilde{\mathcal{W}}}_q(x)$ into (27) to compute ${\mathcal{P}}_{mn}(x)$.

(5)

Set $R_{m,n}(x)={\mathcal{P}}_{mn}(x)/{\tilde{v}}_n(x)$.

Example 16.

Let $$\begin{array}{}\text{(34)}& f\left(x\right)=\left(\begin{array}{cccccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{1}& -\mathrm{1}\\ \mathrm{0}& \mathrm{2}& -\mathrm{1}& \mathrm{2}& \mathrm{2}& \mathrm{1}\end{array}\right)+\left(\begin{array}{cccccc}\mathrm{1}& \mathrm{0}& \mathrm{2}& \mathrm{1}& \mathrm{1}& -\mathrm{2}\\ \mathrm{1}& \mathrm{2}& -\mathrm{1}& \mathrm{2}& \mathrm{3}& \mathrm{1}\end{array}\right)x+\left(\begin{array}{cccccc}\mathrm{1}& \mathrm{2}& \mathrm{1}& \mathrm{1}& \mathrm{0}& -\mathrm{1}\\ \mathrm{0}& \mathrm{2}& -\mathrm{1}& \mathrm{3}& \mathrm{2}& \mathrm{0}\end{array}\right)x^{\mathrm{2}}+\cdots ={\mathcal{A}}_{\mathrm{0}}+{\mathcal{A}}_{\mathrm{1}}x+{\mathcal{A}}_{\mathrm{2}}{x}^{\mathrm{2}}+\cdots .\end{array}$$

Now we apply Algorithm 15 to compute TPTA of type $(\mathrm{2}/\mathrm{2})$ for this example.

$(1)$

Chose $v_{\mathrm{2}}(x)=x^{\mathrm{2}}-\mathrm{2}x+\mathrm{4}$, $q=m-n+\mathrm{1}=\mathrm{1}.$

$(2)$

Use (19) to compute ${\tilde{v}}_{\mathrm{2}}(x)$: $$\begin{array}{}\text{(35)}& {\tilde{v}}_{\mathrm{2}}\left(x\right)=\mathrm{1}-\mathrm{2}x+\mathrm{4}{x}^{\mathrm{2}}.\end{array}$$

$(3)$

By using (25) and (26) we get$$\begin{array}{}\text{(36)}& {\mathcal{W}}_{\mathrm{1}}\left(x\right)=\left(\begin{array}{cccccc}x-\mathrm{1}& \mathrm{2}& \mathrm{2}x+\mathrm{1}& x-\mathrm{1}& x-\mathrm{2}& \mathrm{1}-\mathrm{2}x\\ x& \mathrm{2}x-\mathrm{2}& -x+\mathrm{1}& \mathrm{2}x-\mathrm{1}& \mathrm{3}x-\mathrm{2}& x-\mathrm{2}\end{array}\right),\end{array}$$

 

and $$\begin{array}{}\text{(37)}& \widetilde{{\mathcal{W}}_{\mathrm{1}}}\left(x\right)=\left(\begin{array}{cccccc}\mathrm{1}-x& \mathrm{2}x& x& \mathrm{1}-x& \mathrm{1}-\mathrm{2}x& x-\mathrm{1}\\ \mathrm{0}& \mathrm{2}-\mathrm{2}x& x-\mathrm{1}& \mathrm{2}-x& \mathrm{2}-\mathrm{2}x& \mathrm{1}-\mathrm{2}x\end{array}\right).\end{array}$$

$(4)$

Substituting ${\tilde{v}}_{\mathrm{2}}(x)$ and ${\tilde{\mathcal{W}}}_{\mathrm{1}}(x)$ into (27), we obtain$$\begin{array}{}\text{(38)}& {\mathcal{P}}_{\mathrm{22}}\left(x\right)={\tilde{v}}_{\mathrm{2}}\left(x\right){\mathcal{A}}_{\mathrm{0}}+x\widetilde{{\mathcal{W}}_{\mathrm{1}}}\left(x\right)=\left(\begin{array}{cccccc}\mathrm{3}{x}^{\mathrm{2}}-x+\mathrm{1}& \mathrm{2}{x}^{\mathrm{2}}& -\mathrm{3}{x}^{\mathrm{2}}+\mathrm{2}x& \mathrm{3}{x}^{\mathrm{2}}-x+\mathrm{1}& \mathrm{2}{x}^{\mathrm{2}}-x+\mathrm{1}& -x^{\mathrm{2}}-\mathrm{1}\\ -\mathrm{2}{x}^{\mathrm{2}}+x& \mathrm{6}{x}^{\mathrm{2}}-\mathrm{2}x+\mathrm{2}& -\mathrm{3}{x}^{\mathrm{2}}+x-\mathrm{1}& \mathrm{7}{x}^{\mathrm{2}}-\mathrm{2}x+\mathrm{2}& \mathrm{4}{x}^{\mathrm{2}}-x+\mathrm{2}& \mathrm{2}{x}^{\mathrm{2}}-x+\mathrm{1}\end{array}\right).\end{array}$$

$(5)$

Set $R_{\mathrm{22}}(x)={\mathcal{P}}_{\mathrm{22}}(x)/{\tilde{v}}_{\mathrm{2}}(x)$. It is easy to verify that $$\begin{array}{}\text{(39)}& \frac{{\mathcal{P}}_{\mathrm{22}}\left(x\right)}{{\tilde{v}}_{\mathrm{2}}\left(x\right)}-f\left(x\right)=O\left(x^{\mathrm{3}}\right).\end{array}$$

3.2. Algorithm for Computing TPTA

Generally, the precision of TPTA is limited, since the denominator polynomials of TPTA are arbitrarily prescribed. In this subsection, in order to improve the precision of approximation, we propose an algorithm for computing the denominator polynomials and illustrate the efficiency of this algorithm in next section.

First, we give the following conclusion.

Theorem 17 (error formula).

$$\begin{array}{}\text{(40)}& f\left(x\right)-{\left(n-\frac{\mathrm{1}}{n}\right)}_f\left(x\right)=\frac{x^n}{{\tilde{v}}_n\left(x\right)}\varphi \left(\frac{v_n\left(t\right)}{\mathrm{1}-tx}\right).\end{array}$$

Proof.

Note that $\varphi$ is a linear functional on $\mathbf{P}$, only acting on $t$. From (18) and (20) we deduced that $$\begin{array}{}\text{(41)}& {\tilde{\mathcal{W}}}_{n-\mathrm{1}}\left(x\right)=x^{n-\mathrm{1}}{\mathcal{W}}_{n-\mathrm{1}}\left(x^{-\mathrm{1}}\right)=x^{n-\mathrm{1}}\varphi \left(\frac{x{v}_n\left(x^{-\mathrm{1}}\right)-x{v}_n\left(t\right)}{\mathrm{1}-tx}\right)=\varphi \left(\frac{x^n{v}_n\left(x^{-\mathrm{1}}\right)-x^n{v}_n\left(t\right)}{\mathrm{1}-tx}\right)=\varphi \left(\frac{x^n{v}_n\left(x^{-\mathrm{1}}\right)}{\mathrm{1}-tx}\right)-x^n\varphi \left(\frac{v_n\left(t\right)}{\mathrm{1}-tx}\right)={\tilde{v}}_n\left(x\right)f\left(x\right)-x^n\varphi \left(\frac{v_n\left(t\right)\text{)}}{\mathrm{1}-tx}\right),\end{array}$$ and then this error formula holds.

In terms of the error formula, it holds that$$\begin{array}{}\text{(42)}& f\left(x\right)-{\left(n-\frac{\mathrm{1}}{n}\right)}_f\left(x\right)=\frac{x^n}{{\tilde{v}}_n\left(x\right)}\varphi \left(\frac{v_n\left(t\right)}{\mathrm{1}-tx}\right)=\frac{x^n}{{\tilde{v}}_n\left(x\right)}\varphi \left(v_n\left(t\right)+v_n\left(t\right)tx+v_n\left(t\right)t^{\mathrm{2}}{x}^{\mathrm{2}}+\cdots \right)=\frac{x^n}{{\tilde{v}}_n\left(x\right)}\left(\varphi \left(v_n\left(t\right)\right)+\varphi \left(v_n\left(t\right)t\right)x+\varphi \left(v_n\left(t\right)t^{\mathrm{2}}\right)x^{\mathrm{2}}+\cdots \right).\end{array}$$If we impose that $v_n(t)$ satisfies the condition $\varphi (v_n(t))=\mathrm{0}$, then the first term of (42) disappears, and the order of approximation becomes $n+\mathrm{1}$. If, in addition, we also impose the condition $\varphi (t{v}_n(t))=\mathrm{0}$, the second term in the expansion of the error also disappears, and the order of approximation becomes $n+\mathrm{2}$, and so on. We indicate that $v_n(x)$ depends on $n+\mathrm{1}$ arbitrary constants; however, on the other side, a rational function is defined apart from a multiplying factor in its numerator and its denominator. It implies that $(n-\mathrm{1}/n{)}_f(x)$ depends on $n$ arbitrary constants. So let us take $v_n(t)$ such that$$\begin{array}{}\text{(43)}& \varphi \left(v_n\left(t\right)t^k\right)=\mathrm{0},\hspace{1em}k=\mathrm{0,1},\mathrm{2},\dots ,n-\mathrm{1}.\end{array}$$

Definition 18.

$v_n(t)$ in (43) is called an orthogonal polynomial with respect to the linear functional $\varphi$ and $(n-\mathrm{1}/n{)}_f(x)$ in (42) is also called a TPTA for the given power series (15) when (43) is satisfied.

From (43) we obtain$$\begin{array}{}\text{(44)}& \varphi \left(v_n\left(t\right)t^k\right)={\displaystyle \sum }_{i=\mathrm{0}}^n{b}_i{\mathcal{A}}_{i+k}=\mathrm{0},\hspace{1em}k=\mathrm{0,1},\mathrm{2},\dots ,n-\mathrm{1}.\end{array}$$Let $b_n=\mathrm{1}$ in (44); then it follows that$$\begin{array}{}\text{(45)}& {\displaystyle \sum }_{i=\mathrm{0}}^{n-\mathrm{1}}{b}_i{\mathcal{A}}_{i+k}=-{\mathcal{A}}_{k+n},\hspace{1em}k=\mathrm{0,1},\mathrm{2},\dots ,n-\mathrm{1}.\end{array}$$Forming the scalar product of both sides of (45) with ${\mathcal{A}}_{\mathrm{0}},{\mathcal{A}}_{\mathrm{1}},\dots ,{\mathcal{A}}_{n-\mathrm{1}}$, respectively, we get$$\begin{array}{}\text{(46)}& {\displaystyle \sum }_{i=\mathrm{0}}^{n-\mathrm{1}}{b}_i\left({\mathcal{A}}_{i+k},{\mathcal{A}}_k\right)=-\left({\mathcal{A}}_{k+n},{\mathcal{A}}_k\right),\hspace{1em}k=\mathrm{0,1},\mathrm{2},\dots ,n-\mathrm{1}.\end{array}$$Denote$$\begin{array}{c}\text{(47)}& H_n\left({\mathcal{A}}_{\mathrm{0}}\right)=\left(\begin{array}{cccc}\left({\mathcal{A}}_{\mathrm{0}},{\mathcal{A}}_{\mathrm{0}}\right)& \left({\mathcal{A}}_{\mathrm{1}},{\mathcal{A}}_{\mathrm{0}}\right)& \cdots & \left({\mathcal{A}}_{n-\mathrm{1}},{\mathcal{A}}_{\mathrm{0}}\right)\\ \left({\mathcal{A}}_{\mathrm{1}},{\mathcal{A}}_{\mathrm{1}}\right)& \left({\mathcal{A}}_{\mathrm{2}},{\mathcal{A}}_{\mathrm{1}}\right)& \cdots & \left({\mathcal{A}}_n,{\mathcal{A}}_{\mathrm{1}}\right)\\ ⋮& ⋮& \ddots & ⋮\\ \left({\mathcal{A}}_{n-\mathrm{1}},{\mathcal{A}}_{n-\mathrm{1}}\right)& \left({\mathcal{A}}_n,{\mathcal{A}}_{n-\mathrm{1}}\right)& \cdots & \left({\mathcal{A}}_{\mathrm{2}n-\mathrm{2}},{\mathcal{A}}_{n-\mathrm{1}}\right)\end{array}\right),{\overrightarrow{h}}_n=\left(\begin{array}{c}-\left({\mathcal{A}}_n,{\mathcal{A}}_{\mathrm{0}}\right)\\ -\left({\mathcal{A}}_{n+\mathrm{1}},{\mathcal{A}}_{\mathrm{1}}\right)\\ ⋮\\ -\left({\mathcal{A}}_{\mathrm{2}n-\mathrm{1}},{\mathcal{A}}_{n-\mathrm{1}}\right)\end{array}\right),\end{array}$$and call det$(H_n({\mathcal{A}}_{\mathrm{0}}))$ the Hankel determinant of $f(x)$ with respect to the coefficients ${\mathcal{A}}_{\mathrm{0}},{\mathcal{A}}_{\mathrm{1}},\dots ,{\mathcal{A}}_{n-\mathrm{1}}$.

Then (46) is converted into$$\begin{array}{}\text{(48)}& H_n\left({\mathcal{A}}_{\mathrm{0}}\right)x={\overrightarrow{h}}_n,\hspace{1em}x={\left(b_{\mathrm{0}},b_{\mathrm{1}},\dots ,b_{n-\mathrm{1}}\right)}^T.\end{array}$$

In the case of TPTA, $v_n(x)$ is not arbitrarily chosen any more but is determined by the preceding system. The choice of $v_n(x)$ can help to improve the accuracy of approximation, but unfortunately we have not been able to guarantee that the solution of system (48) comes into existence, so far. We only give the following basic theorem about system (48) on the basis of linear algebra.

Theorem 19.

The solution of (48) exists if and only if $rank(H_n({\mathcal{A}}_{\mathrm{0}}))=rank(H_n({\mathcal{A}}_{\mathrm{0}}):{\overrightarrow{h}}_n)$. Moreover, the solution is unique if det$(H_n({\mathcal{A}}_{\mathrm{0}}))\ne \mathrm{0}$.

Proof.

The proof of the assertion follows from the simple fact that, for a system of linear equations, described by $Ax=b$, where $A$ is matrix and $x$, $b$ are vectors, the solution of system (48) comes into existence for $x$ if and only if $rank(A)=rank(A:b)$; i.e., the right-hand vector must be in the vector space spanned by the columns of the coefficient matrix $A$. Moreover, if det$(H_n({\mathcal{A}}_{\mathrm{0}}))\ne \mathrm{0}$, according to Cramer’s rule, the solution is unique.

Theorem 20 (existence).

Let $f(x)$ be given power series (15); then ${(n-\mathrm{1}/n)}_f$ and ${(m/n)}_f$ exist and are unique if and only if $det(H_n({\mathcal{A}}_{\mathrm{0}}))\ne \mathrm{0}$.

Proof.

“$\Leftarrow$” By Theorem 19, if $det(H_n({\mathcal{A}}_{\mathrm{0}}))\ne \mathrm{0}$, it means that nonhomogeneous equation (48) exists as a unique solution $b_{\mathrm{0}},b_{\mathrm{1}},\dots ,b_{n-\mathrm{1}}$ for $v_n(x)$. From (19), it also means that ${\tilde{v}}_n(x)$ and ${\tilde{\mathcal{W}}}_{n-\mathrm{1}}(x)$ exist. Hence, by the construction of $(n-\mathrm{1}/n{)}_f$, existence holds.

“$\Rightarrow$” Let $(n-\mathrm{1}/n{)}_f={\tilde{\mathcal{W}}}_{n-\mathrm{1}}(x)/{\tilde{v}}_{n-\mathrm{1}}(x)$ exist and be unique; then it implies that $$\begin{array}{}\text{(49)}& rank\left(H_n\left({\mathcal{A}}_{\mathrm{0}}\right)\right)=rank\left(H_n\left({\mathcal{A}}_{\mathrm{0}}\right):{\overrightarrow{h}}_n\right),\end{array}$$ and if $det(H_n({\mathcal{A}}_{\mathrm{0}}))=\mathrm{0}$, then $rank(H_n({\mathcal{A}}_{\mathrm{0}}):{\overrightarrow{h}}_n)<n$; the fact that equation (46) has $n-rank(H_n({\mathcal{A}}_{\mathrm{0}}):{\overrightarrow{h}}_n)$ solutions, namely, that we can construct $n-rank(H_n({\mathcal{A}}_{\mathrm{0}}):{\overrightarrow{h}}_n)$ generating polynomials, which is contradictory to the uniqueness of $(n-\mathrm{1}/n{)}_f$, holds.

The proof of existence and uniqueness of $(m/n{)}_f$ is similar to the preceding process.

Theorem 21.

Let det$(H_n({\mathcal{A}}_{\mathrm{0}}))\ne \mathrm{0}$; then $(n-\mathrm{1}/n{)}_f(x)={\tilde{\mathcal{W}}}_{n-\mathrm{1}}(x)/{\tilde{v}}_n(x)$, where the generating polynomial $v_n(x)$ is given by$$\begin{array}{}\text{(50)}& v_n\left(x\right)=det\left(\begin{array}{cc}{H}_n\left({\mathcal{A}}_{\mathrm{0}}\right)& \overrightarrow{\gamma }\\ {\overrightarrow{\delta }}^T& x^n\end{array}\right),\end{array}$$where $\overrightarrow{\gamma }={\left(\left({\mathcal{A}}_n,{\mathcal{A}}_{\mathrm{0}}\right),\left({\mathcal{A}}_{n+\mathrm{1}},{\mathcal{A}}_{\mathrm{1}}\right),\dots ,\left({\mathcal{A}}_{\mathrm{2}n-\mathrm{1}},{\mathcal{A}}_{n-\mathrm{1}}\right)\right)}^T$, ${\overrightarrow{\delta }}^T=(\mathrm{1},x,\dots ,x^{n-\mathrm{1}})$, and ${\tilde{v}}_n(x)$ and ${\tilde{\mathcal{W}}}_{n-\mathrm{1}}(x)$ are given by (19) and (20), respectively.

Now, we can derive an algorithm to calculate $(m/n{)}_f(x)$ using (26), (27), and (50).

Algorithm 22 (compute $(m/n{)}_f(x)={\mathcal{P}}_{mn}(x)/{\tilde{v}}_n(x)$).

(1)

Use (14) to calculate $H_n({\mathcal{A}}_{\mathrm{0}})$ and $\overrightarrow{\gamma }$.

(2)

Use (50) and (19) to compute $v_n(x)$ and ${\tilde{v}}_n(x)$, respectively.

(3)

Set $q=m-n+\mathrm{1}$ and compute ${\mathcal{W}}_q(x)$ and ${\tilde{\mathcal{W}}}_q(x)$ by $$\begin{array}{c}\text{(51)}& {\mathcal{W}}_q\left(x\right)={\varphi }^{\left(m-n+\mathrm{1}\right)}\left(\frac{v_n\left(x\right)-v_n\left(t\right)}{x-t}\right),{\tilde{\mathcal{W}}}_q\left(x\right)=x^{n-\mathrm{1}}{\mathcal{W}}_q\left(x^{-\mathrm{1}}\right).\end{array}$$

(4)

Compute the numerator of TPTA by $$\begin{array}{}\text{(52)}& {\mathcal{P}}_{mn}\left(x\right)={\tilde{v}}_n\left(x\right){\displaystyle \sum }_{i=\mathrm{0}}^{m-n}{\mathcal{A}}_i{x}^i+x^{m-n+\mathrm{1}}{\tilde{\mathcal{W}}}_q\left(x\right).\end{array}$$

(5)

Obtain $(m/n{)}_f(x)={\mathcal{P}}_{mn}(x)/{\tilde{v}}_n(x)$.

4. Application for Computing the Tensor Exponential Function

The method of truncated infinite series has abroad applications in finite single crystal plasticity for computing tensor exponential function [16]. However, the accuracy and effectiveness of such algorithm are limited by round-off and choice of termination criterion. In this section, we will utilize the method of TPTA to compute tensor exponential function. We start by briefly reviewing some basic equations that model the behaviour of single crystals in the finite strain range [16].

Consider a single crystal model$$\begin{array}{}\text{(53)}& F=F^e{F}^p,\end{array}$$where $F^e$ and $F^p$ denote elastic part and plastic part, respectively.

For a single crystal with a total number $n_{syst}$ of slip systems, the evolution of the inelastic deformation gradient, $F^p$, is defined by means of the following rate form:$$\begin{array}{}\text{(54)}& {\dot{F}}^p{F}^{p-\mathrm{1}}={\displaystyle \sum }_{\alpha =\mathrm{1}}^{n_{syst}}{\dot{\gamma }}^{\alpha }{s}_{\mathrm{0}}^{\alpha }\otimes m_{\mathrm{0}}^{\alpha },\end{array}$$where ${\dot{\gamma }}^{\alpha }$ denotes the contribution of slip system $\alpha$ to the total inelastic rate of deformation. The vectors $s_{\mathrm{0}}^{\alpha }$ and $m_{\mathrm{0}}^{\alpha }$ denote, respectively, the slip direction and normal direction of slip system $\alpha$.

The above tensor differential equation can be discretized in an implicit fashion with use of the tensor exponential function. The implicit exponential approximation to the inelastic flow equation results in the following discrete form:$$\begin{array}{}\text{(55)}& F_{n+\mathrm{1}}^p=\mathrm{e}\mathrm{x}\mathrm{p}\left({\displaystyle \sum }_{\alpha =\mathrm{1}}^{n_{syst}}\Delta {\gamma }^{\alpha }{s}_{\mathrm{0}}^{\alpha }\otimes m_{\mathrm{0}}^{\alpha }\right)F_n^p.\end{array}$$

The above formula is analogous to the exact solution of initial value problem (1) and it is necessary to calculate$$\begin{array}{}\text{(56)}& \mathrm{e}\mathrm{x}\mathrm{p}\left({\displaystyle \sum }_{\alpha =\mathrm{1}}^{n_{syst}}\Delta {\gamma }^{\alpha }{s}_{\mathrm{0}}^{\alpha }\otimes m_{\mathrm{0}}^{\alpha }\right).\end{array}$$

In [7], the author used Algorithm 23 to calculate (56).

Algorithm 23 (truncated infinite series method [7] (p.749)).

(1)

Given tensor $\mathcal{X}$, initialise $n=\mathrm{0}$ and $\mathrm{e}\mathrm{x}\mathrm{p}\left(\mathcal{X}\right)≔\mathcal{I}$.

(2)

Increment counter $n≔n+\mathrm{1}$.

(3)

Compute $n!$ and ${\mathcal{X}}^n$.

(4)

Add new term to the series $\mathrm{e}\mathrm{x}\mathrm{p}\left(\mathcal{X}\right)≔\mathrm{e}\mathrm{x}\mathrm{p}(\mathcal{X})+\left(\mathrm{1}/n!\right){\mathcal{X}}^n$.

(5)

Check convergence, if $$\begin{array}{}\text{(57)}& \frac{‖{\mathcal{X}}^n‖}{n!}<{ϵ}_{tol},\hspace{1em}then\hspace{0.17em}\hspace{0.17em}exit,\hspace{0.17em}\hspace{0.17em}else\hspace{0.17em}\hspace{0.17em}goto\hspace{0.17em}\hspace{0.17em}\left(\mathrm{1}\right).\end{array}$$

Example 24.

Consider a tensor exponential function $\mathrm{e}\mathrm{x}\mathrm{p}(\mathcal{A}x)$; the entries of $\mathcal{A}$ are $a_{\mathrm{121}}=\mathrm{1}$, $a_{\mathrm{221}}=-\mathrm{2}$, $a_{\mathrm{122}}=\mathrm{2}$, $a_{\mathrm{222}}=-\mathrm{1}$, and zero elsewhere.

To find a tensor Pad$\acute{e}$-type approximation of type $(\mathrm{3}/\mathrm{3})$ for the tensor exponential function, first we should expand $\mathrm{e}\mathrm{x}\mathrm{p}(\mathcal{A}x)$ into power series by means of Definition 7. We can obtain$$\begin{array}{}\text{(58)}& \mathrm{e}\mathrm{x}\mathrm{p}\left(\mathcal{A}x\right)=\left(\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\end{array}\right)+\left(\begin{array}{cccc}\mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{2}\\ \mathrm{0}& -\mathrm{2}& \mathrm{0}& -\mathrm{1}\end{array}\right)x+\left(\begin{array}{cccc}\mathrm{0}& -\mathrm{2}& \mathrm{0}& -\frac{\mathrm{5}}{\mathrm{2}}\\ \mathrm{0}& \frac{\mathrm{5}}{\mathrm{2}}& \mathrm{0}& \mathrm{2}\end{array}\right)x^{\mathrm{2}}+\left(\begin{array}{cccc}\mathrm{0}& \frac{\mathrm{13}}{\mathrm{6}}& \mathrm{0}& \frac{\mathrm{14}}{\mathrm{6}}\\ \mathrm{0}& -\frac{\mathrm{14}}{\mathrm{6}}& \mathrm{0}& -\frac{\mathrm{13}}{\mathrm{6}}\end{array}\right)x^{\mathrm{3}}+\left(\begin{array}{cccc}\mathrm{0}& -\frac{\mathrm{40}}{\mathrm{24}}& \mathrm{0}& -\frac{\mathrm{41}}{\mathrm{24}}\\ \mathrm{0}& \frac{\mathrm{41}}{\mathrm{24}}& \mathrm{0}& \frac{\mathrm{40}}{\mathrm{24}}\end{array}\right)x^{\mathrm{4}}+\cdots ={\mathcal{A}}_{\mathrm{0}}+{\mathcal{A}}_{\mathrm{1}}x+{\mathcal{A}}_{\mathrm{2}}{x}^{\mathrm{2}}+{\mathcal{A}}_{\mathrm{3}}{x}^{\mathrm{3}}+{\mathcal{A}}_{\mathrm{4}}{x}^{\mathrm{4}}+\cdots .\end{array}$$ By Algorithm 22, the following can be done.

$(1)$

Use (14) to compute $H_{\mathrm{3}}({\mathcal{A}}_{\mathrm{0}})$ and $\overrightarrow{\gamma }$: $$\begin{array}{c}\text{(59)}& H_{\mathrm{3}}\left({\mathcal{A}}_{\mathrm{0}}\right)=\left(\begin{array}{ccc}\mathrm{2}& -\mathrm{2}& \frac{\mathrm{5}}{\mathrm{2}}\\ \mathrm{10}& -\mathrm{14}& \frac{\mathrm{41}}{\mathrm{3}}\\ \frac{\mathrm{41}}{\mathrm{2}}& -\frac{\mathrm{61}}{\mathrm{3}}& \frac{\mathrm{365}}{\mathrm{24}}\end{array}\right),\overrightarrow{\gamma }=\left(\begin{array}{c}-\frac{\mathrm{7}}{\mathrm{3}}\\ -\frac{\mathrm{61}}{\mathrm{6}}\\ -\frac{\mathrm{547}}{\mathrm{60}}\end{array}\right).\end{array}$$

$(2)$

Use (50) to calculate $v_{\mathrm{3}}(x)$: $$\begin{array}{}\text{(60)}& v_{\mathrm{3}}\left(x\right)=\frac{\mathrm{15041}}{\mathrm{1080}}+\frac{\mathrm{3947}}{\mathrm{60}}x+\frac{\mathrm{1189}}{\mathrm{10}}{x}^{\mathrm{2}}+\frac{\mathrm{1493}}{\mathrm{18}}{x}^{\mathrm{3}},\end{array}$$

and  

compute $\widetilde{v_{\mathrm{3}}}(x)$ by (19), so we get $$\begin{array}{}\text{(61)}& {\tilde{v}}_{\mathrm{3}}\left(x\right)=\frac{\mathrm{15041}}{\mathrm{1080}}{x}^{\mathrm{3}}+\frac{\mathrm{3947}}{\mathrm{60}}{x}^{\mathrm{2}}+\frac{\mathrm{1189}}{\mathrm{10}}x+\frac{\mathrm{1493}}{\mathrm{18}}.\end{array}$$

$(3)$

Set $q=m-n+\mathrm{1}=\mathrm{1}$ and compute ${\mathcal{W}}_{\mathrm{1}}(x)$, ${\tilde{\mathcal{W}}}_{\mathrm{1}}(x)$: $$\begin{array}{}\text{(62)}& {\mathcal{W}}_{\mathrm{1}}\left(x\right)=\left(\begin{array}{cccc}\mathrm{0}& \frac{\mathrm{1493}}{\mathrm{18}}{x}^{\mathrm{2}}-\frac{\mathrm{4229}}{\mathrm{90}}x+\frac{\mathrm{1039}}{\mathrm{135}}& \mathrm{0}& \frac{\mathrm{1493}}{\mathrm{9}}{x}^{\mathrm{2}}+\frac{\mathrm{5479}}{\mathrm{180}}x+\frac{\mathrm{15041}}{\mathrm{540}}\\ \mathrm{0}& -\frac{\mathrm{1493}}{\mathrm{9}}{x}^{\mathrm{2}}-\frac{\mathrm{5479}}{\mathrm{180}}x-\frac{\mathrm{15041}}{\mathrm{540}}& \mathrm{0}& -\frac{\mathrm{1493}}{\mathrm{18}}{x}^{\mathrm{2}}+\frac{\mathrm{4229}}{\mathrm{90}}x-\frac{\mathrm{1039}}{\mathrm{135}}\end{array}\right),\end{array}$$

 

and$$\begin{array}{}\text{(63)}& {\tilde{\mathcal{W}}}_{\mathrm{1}}\left(x\right)=\left(\begin{array}{cccc}\mathrm{0}& \frac{\mathrm{1493}}{\mathrm{18}}-\frac{\mathrm{4229}}{\mathrm{90}}x+\frac{\mathrm{1039}}{\mathrm{135}}{x}^{\mathrm{2}}& \mathrm{0}& \frac{\mathrm{1493}}{\mathrm{9}}+\frac{\mathrm{5479}}{\mathrm{180}}x+\frac{\mathrm{15041}}{\mathrm{540}}{x}^{\mathrm{2}}\\ \mathrm{0}& -\frac{\mathrm{1493}}{\mathrm{9}}-\frac{\mathrm{5479}}{\mathrm{180}}x-\frac{\mathrm{15041}}{\mathrm{540}}{x}^{\mathrm{2}}& \mathrm{0}& -\frac{\mathrm{1493}}{\mathrm{18}}+\frac{\mathrm{4229}}{\mathrm{90}}x-\frac{\mathrm{1039}}{\mathrm{135}}{x}^{\mathrm{2}}\end{array}\right).\end{array}$$

$(4)$

Substitute ${\tilde{v}}_{\mathrm{3}}(x)$ and ${\tilde{\mathcal{W}}}_{\mathrm{1}}(x)$ into (27) to compute ${\mathcal{P}}_{\mathrm{33}}(x)$: $$\begin{array}{}\text{(64)}& {\mathcal{P}}_{\mathrm{33}}\left(x\right)=\left(\begin{array}{cccc}\frac{\mathrm{15041}}{\mathrm{1080}}{x}^{\mathrm{3}}+\frac{\mathrm{3947}}{\mathrm{60}}{x}^{\mathrm{2}}+\frac{\mathrm{1189}}{\mathrm{10}}x+\frac{\mathrm{1493}}{\mathrm{18}}& \frac{\mathrm{1039}}{\mathrm{135}}{x}^{\mathrm{3}}-\frac{\mathrm{4229}}{\mathrm{90}}{x}^{\mathrm{2}}+\frac{\mathrm{1493}}{\mathrm{18}}x& \mathrm{0}& \frac{\mathrm{15041}}{\mathrm{540}}{x}^{\mathrm{3}}+\frac{\mathrm{5479}}{\mathrm{180}}{x}^{\mathrm{2}}+\frac{\mathrm{1493}}{\mathrm{9}}x\\ \mathrm{0}& -\frac{\mathrm{15041}}{\mathrm{1080}}{x}^{\mathrm{3}}+\frac{\mathrm{3181}}{\mathrm{90}}{x}^{\mathrm{2}}-\frac{\mathrm{4229}}{\mathrm{90}}x+\frac{\mathrm{1493}}{\mathrm{18}}& \mathrm{0}& -\frac{\mathrm{1039}}{\mathrm{135}}{x}^{\mathrm{3}}+\frac{\mathrm{4229}}{\mathrm{90}}{x}^{\mathrm{2}}-\frac{\mathrm{1493}}{\mathrm{18}}x\end{array}\right).\end{array}$$

$(5)$

Obtain $(m/n{)}_f(x)={\mathcal{P}}_{\mathrm{33}}(x)/{\tilde{v}}_{\mathrm{3}}(x)$.

In Table 1 we compare the number of exact figures given by the method of TPTA of type $(\mathrm{3}/\mathrm{3})$ with corresponding exact value of $\mathrm{e}\mathrm{x}\mathrm{p}(\mathcal{A}x)$ referring to the entries of $(\mathrm{1,2},\mathrm{1})$, $(\mathrm{2,2},\mathrm{1})$, $(\mathrm{1,2},\mathrm{2})$, and $(\mathrm{2,2},\mathrm{2})$. We also compute the norm of absolute residual tensor (denoted by $Res$). Here, $$\begin{array}{}\text{(65)}& Res\left(x_i\right)={‖\exp \left(\mathcal{A}{x}_i\right)-{\left(\frac{\mathrm{3}}{\mathrm{3}}\right)}_{e^{\mathcal{A}x}}\left(x_i\right)‖}^2,\end{array}$$ where the operation $‖·‖$ is defined by (13).

Table 1

Numerical results of Example 24 at different points by using Algorithm 22.

$$x$$		$(\mathrm{1,2},\mathrm{1})$	$(\mathrm{2,2},\mathrm{1})$	$$(\mathrm{1,2},\mathrm{2})$$	$(\mathrm{2,2},\mathrm{2})$	$$Res$$
$$\mathrm{0.1}$$	$$e^{\mathcal{A}x}$$	0.08200959	0.82282781	0.17717281	-0.08200959	8.34e-12
	$(\mathrm{3}/\mathrm{3}{)}_{e^{\mathcal{A}x}}(x)$	0.08200778	0.82282688	0.17717311	-0.08200778
0.2	$$e^{\mathcal{A}x}$$	0.13495955	0.68377119	0.31622880	-0.13495955	1.62e-9
	$$(\mathrm{3}/\mathrm{3}{)}_{e^{\mathcal{A}x}}(x)$$	0.13493452	0.68375764	0.31624235	-0.13493452
$\mathrm{0.3}$	$$e^{\mathcal{A}x}$$	0.16712428	0.57369394	0.42630605	-0.16712428	3.14e-8
	$$(\mathrm{3}/\mathrm{3}{)}_{e^{\mathcal{A}x}}(x)$$	0.16701602	0.57363058	0.42636941	-0.16701602
$\mathrm{0.4}$	$$e^{\mathcal{A}x}$$	0.18456291	0.48575712	0.51424287	-0.18456291	2.38e-7
	$$(\mathrm{3}/\mathrm{3}{)}_{e^{\mathcal{A}x}}(x)$$	0.18427224	0.48557038	0.51442961	-0.18427224

From Table 1, it is observed that the estimates from TPTA can reach the desired accuracy.

Example 25.

Let $f(x)$ be given by Example 24.

By Algorithm 22 for preceding example again, we calculate $(m/m{)}_f(\mathrm{1}),m=\mathrm{1,2},\mathrm{3,4},\mathrm{5}$. The exact value and approximant value associated with the entries of $(\mathrm{1,2},\mathrm{1})$, $(\mathrm{2,2},\mathrm{1})$, $(\mathrm{1,2},\mathrm{2})$, and $(\mathrm{2,2},\mathrm{2})$ are listed in Tables 2 and 3, respectively, where $x=\mathrm{1}$.

Table 2

The exact value of $\mathrm{e}\mathrm{x}\mathrm{p}(\mathcal{A})$.

	$$(\mathrm{1,2},\mathrm{1})$$	$$(\mathrm{2,2},\mathrm{1})$$	$$(\mathrm{1,2},\mathrm{2})$$	$$(\mathrm{2,2},\mathrm{2})$$
$$\mathrm{e}\mathrm{x}\mathrm{p}(\mathcal{A})$$	$$\mathrm{0.15904705}$$	$$\mathrm{0.20883238}$$	$$\mathrm{0.79116761}$$	$-\mathrm{0.15904705}$

Table 3

Numerical approximations of $\exp \left(\mathcal{A}\right)$ using Algorithm 22 for Example 25.

$$m$$	$$(\mathrm{1,2},\mathrm{1})$$	$$(\mathrm{2,2},\mathrm{1})$$	$$(\mathrm{1,2},\mathrm{2})$$	$$(\mathrm{2,2},\mathrm{2})$$	$$Res$$
$$\mathrm{1}$$	$$\mathrm{0.50000000}$$	$$\mathrm{0.00000000}$$	$$\mathrm{1.00000000}$$	$$-\mathrm{0.50000000}$$	3.19e-1
$$\mathrm{2}$$	$$-\mathrm{2.12500000}$$	$$\mathrm{0.75000000}$$	$$\mathrm{0.25000000}$$	$$\mathrm{2.12500000}$$	1.10e+1
$$\mathrm{3}$$	$$\mathrm{0.15503865}$$	$$\mathrm{0.20377270}$$	$$\mathrm{0.79622729}$$	$$-\mathrm{0.15503865}$$	8.33e-5
$$\mathrm{4}$$	$$\mathrm{0.17454584}$$	$$\mathrm{0.22682585}$$	$$\mathrm{0.77317481}$$	$$-\mathrm{0.17454584}$$	1.12e-3
$$\mathrm{5}$$	$$\mathrm{0.17625313}$$	$$\mathrm{0.19112365}$$	$$\mathrm{0.80887636}$$	$$-\mathrm{0.17625313}$$	1.21e-3

From Table 3, we can see that ${(\mathrm{3}/\mathrm{3})}_{e^{\mathcal{A}x}}$ has the best approximation for this example. We also compute $\mathrm{e}\mathrm{x}\mathrm{p}(\mathcal{A})$ by using Algorithm 23, and the corresponding numerical results are listed in Table 4. By comparison of Table 3 with Table 4, we find that it requires at most 6 coefficients (since $m=\mathrm{3}$) of power series expansion of $\mathrm{e}\mathrm{x}\mathrm{p}(\mathcal{A}x)$ to achieve an error of $\mathrm{1}{\mathrm{0}}^{-\mathrm{5}}$ in Algorithm 22, while requiring 11 coefficients in Algorithm 23. It is straightforward to understand that Algorithm 23 is more expensive than Algorithm 22 especially for higher order tensor exponential function. In practical applications, only few coefficients of the series may be known, so, we may get the desired results by means of TPTA. Thus the effectiveness of the proposed Algorithm 22 is verified.

Table 4

Numerical results of Example 25 by using Algorithm 23.

$$n_{max}$$	$$(\mathrm{1,2},\mathrm{1})$$	$$(\mathrm{2,2},\mathrm{1})$$	$$(\mathrm{1,2},\mathrm{2})$$	$$(\mathrm{2,2},\mathrm{2})$$	$$Res$$
$$\mathrm{1}$$	$$\mathrm{1}$$	$$-\mathrm{1}$$	$$\mathrm{2}$$	$$\mathrm{2}$$	4.33
$$\mathrm{2}$$	$$-\mathrm{1}$$	$$\mathrm{1.5}$$	$$-\mathrm{0.5}$$	$$\mathrm{1}$$	6.02
$$\mathrm{3}$$	$$\mathrm{1.16666666}$$	$$-\mathrm{0.83333333}$$	$$\mathrm{1.83333333}$$	$$-\mathrm{1.16666666}$$	4.20
$$\mathrm{4}$$	$$-\mathrm{0.50000000}$$	$$\mathrm{0.87500000}$$	$$\mathrm{0.12500000}$$	$$\mathrm{0.50000000}$$	1.75
$$\mathrm{5}$$	$$\mathrm{0.50833333}$$	$$-\mathrm{0.14166666}$$	$$\mathrm{1.14166666}$$	$$-\mathrm{0.50833333}$$	4.89e-1
$$\mathrm{6}$$	$$\mathrm{0.00277777}$$	$$\mathrm{0.36527777}$$	$$\mathrm{0.63472222}$$	$$-\mathrm{0.00277777}$$	9.77e-2
$$\mathrm{7}$$	$$\mathrm{0.21964285}$$	$$\mathrm{0.14821428}$$	$$\mathrm{0.85178571}$$	$$-\mathrm{0.21964285}$$	1.47e-2
$$\mathrm{8}$$	$$\mathrm{0.13829365}$$	$$\mathrm{0.22958829}$$	$$\mathrm{0.77041170}$$	$$-\mathrm{0.13829365}$$	1.72e-3
$$\mathrm{9}$$	$$\mathrm{0.16541280}$$	$$\mathrm{0.20246638}$$	$$\mathrm{0.79753361}$$	$$-\mathrm{0.16541280}$$	1.62e-4
$$\mathrm{10}$$	$$\mathrm{0.15735119}$$	$$\mathrm{0.21060267}$$	$$\mathrm{0.78939732}$$	$$-\mathrm{0.15735119}$$	1.22e-5
$$\mathrm{11}$$	$$\mathrm{0.15957013}$$	$$\mathrm{0.20838371}$$	$$\mathrm{0.79161628}$$	$$-\mathrm{0.15957013}$$	8.77e-7

5. Conclusion

In this paper, we presented tensor Pad$\acute{e}$-type approximant method for computing tensor exponential function; the expression of TPTA is of the form of tensor numerator and scalar denominator. In order to have a tensor Pad$\acute{e}$-type approximant with the higher possible precision of approximation, we proposed an algorithm for computing denominator polynomials of TPTA, and its effectiveness has been investigated in one example of tensor exponential function. The key to the TPTA to be applied to the tensor exponential function is that it can be expanded into power series with the same order tensors coefficients by means of t-product. Of course, there are several ways to multiply tensors [26–30], but the order of the resulting tensor may be changed. For example, if $\mathcal{A}$ is $n_{\mathrm{1}}\times n_{\mathrm{2}}\times n_{\mathrm{3}}$ and $\mathcal{B}$ is $n_{\mathrm{1}}\times m_{\mathrm{2}}\times m_{\mathrm{3}}$, then the contracted product [26] of $\mathcal{A}$ and $\mathcal{B}$ is $n_{\mathrm{2}}\times n_{\mathrm{3}}\times m_{\mathrm{2}}\times m_{\mathrm{3}}$. So, the choice of the multiplication of two tensors is an open question for expanding tensor exponential function, and the corresponding tensor Pad$\acute{e}$ approximant theoretic is a subject of further research.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The work is supported by National Natural Science Foundation (11371243) and Key Disciplines of Shanghai Municipality (S30104).

References

[1] S. Dolgov, B. Khoromskij, "Simultaneous state-time approximation of the chemical master equation using tensor product formats," Numerical Linear Algebra with Applications, vol. 22 no. 2, pp. 197-219, DOI: 10.1002/nla.1942, 2015.

[2] N. G. van Kampen, Stochastic Process in Physics and Chemistry, 2007.

[3] P. GelB, S. Matera, C. Schütte, "Solving the master equation without kinetic Monte Carlo: tensor train approximations for a CO oxidation model," Journal of Computational Physics, vol. 314, pp. 489-502, DOI: 10.1016/j.jcp.2016.03.025, 2016.

[4] W. Ding, K. Liu, E. Belyaev, F. Cheng, "Tensor-based linear dynamical systems for action recognition from 3D skeletons," Pattern Recognition, pp. 75-86, DOI: 10.1016/j.patcog.2017.12.004, 2017.

[5] P. Gelß, S. Klus, S. Matera, C. Schütte, "Nearest-neighbor interaction systems in the tensor-train format," Journal of Computational Physics, vol. 341, pp. 140-162, DOI: 10.1016/j.jcp.2017.04.007, 2017.

[6] C. Gu, Y. Z. Huang, Z. B. Chen, Continued Fractional Recurrence Algorithm for Generalized Inverse Tensor Padé Approximation, 2018.

[7] E. A. de Souza, D. Peric', D. R. J. Owen, Computational methods for plasticity: Theory and Applications, 2008.

[8] H. Chen, Y. Chen, G. Li, L. Qi, "A semidefinite program approach for computing the maximum eigenvalue of a class of structured tensors and its applications in hypergraphs and copositivity test," Numerical Linear Algebra with Applications, vol. 25 no. 1,DOI: 10.1002/nla.2125, 2018.

[9] G. Zhou, G. Wang, L. Qi, M. Alqahtani, "A fast algorithm for the spectral radii of weakly reducible nonnegative tensors," Numerical Linear Algebra with Applications, vol. 25 no. 2,DOI: 10.1002/nla.2134, 2018.

[10] H. Chen, Y. Wang, "On computing minimal H -eigenvalue of sign-structured tensors," Frontiers of Mathematics in China, vol. 12 no. 6, pp. 1289-1302, DOI: 10.1007/s11464-017-0645-0, 2017.

[11] G. Wang, G. Zhou, L. Caccetta, "Z -eigencvalue inclusion theorems for tensors," Discrete and Continuous Dynamical Systems - Series B, vol. 22 no. 1, pp. 187-198, DOI: 10.3934/dcdsb.2017009, 2017.

[12] K. Zhang, Y. Wang, "An H -tensor based iterative scheme for identifying the positive definiteness of multivariate homogeneous forms," Journal of Computational and Applied Mathematics, vol. 305,DOI: 10.1016/j.cam.2016.03.025, 2016.

[13] Y. Wang, K. Zhang, H. Sun, "Criteria for strong H-tensors," Frontiers of Mathematics in China, vol. 11 no. 3, pp. 577-592, DOI: 10.1007/s11464-016-0525-z, 2016.

[14] H. Chen, L. Qi, Y. Song, "Column sufficient tensors and tensor complementarity problems," Frontiers of Mathematics in China, vol. 13 no. 2, pp. 255-276, DOI: 10.1007/s11464-018-0681-4, 2018.

[15] Y. Wang, L. Caccetta, G. Zhou, "Convergence analysis of a block improvement method for polynomial optimization over unit spheres," Numerical Linear Algebra with Applications, vol. 22 no. 6, pp. 1059-1076, DOI: 10.1002/nla.1996, 2015.

[16] E. A. de Souza Neto, "The exact derivative of the exponential of an unsymmetric tensor," Computer Methods Applied Mechanics and Engineering, vol. 190 no. 18-19, pp. 2377-2383, DOI: 10.1016/S0045-7825(00)00241-3, 2001.

[17] A. Cuitino, M. Ortiz, "A material-independent method for extending stress update algorithms from small-strain plasticity to finite plasticity with multiplicative kinematics," Engineering Computations, vol. 9 no. 4, pp. 437-451, DOI: 10.1108/eb023876, 1992.

[18] A. L. Eterovic, K. Bathe, "A hyperelastic‐based large strain elasto‐plastic constitutive formulation with combined isotropic‐kinematic hardening using the logarithmic stress and strain measures," International Journal for Numerical Methods in Engineering, vol. 30 no. 6, pp. 1099-1114, DOI: 10.1002/nme.1620300602, 1990.

[19] J. C. Simo, "Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping schemes of the infinitesimal theory," Computer Methods Applied Mechanics and Engineering, vol. 99 no. 1, pp. 61-112, DOI: 10.1016/0045-7825(92)90123-2, 1992.

[20] C. Brezinski, Padé-type approximation and general orthogonal polynomials, vol. 50, 1980.

[21] C. Gu, "Matrix Padé-type approximant and directional matrix Padé approximant in the inner product space," Journal of Computational and Applied Mathematics, vol. 164, pp. 365-385, DOI: 10.1016/S0377-0427(03)00487-4, 2004.

[22] M. E. Kilmer, K. Braman, N. Hao, R. C. Hoover, "Third-order tensors as operators on matrices: a theoretical and computational framework with applications in imaging," SIAM Journal on Matrix Analysis and Applications, vol. 34 no. 1, pp. 148-172, DOI: 10.1137/110837711, 2013.

[23] C. D. Martin, R. Shafer, B. Larue, "An order- p tensor factorization with applications in imaging," SIAM Journal on Scientific Computing, vol. 35 no. 1, pp. A474-A490, DOI: 10.1137/110841229, 2013.

[24] M. Itskov, "Computation of the exponential and other isotropic tensor functions and their derivatives," Computer Methods Applied Mechanics and Engineering, vol. 192 no. 35-36, pp. 3985-3999, DOI: 10.1016/S0045-7825(03)00397-9, 2003.

[25] T. G. Kolda, B. W. Bader, "Tensor decompositions and applications," SIAM Review, vol. 51 no. 3, pp. 455-500, DOI: 10.1137/07070111X, 2009.

[26] M. E. Kilmer, C. D. Martin, "Factorization strategies for third-order tensors," Linear Algebra and its Applications, vol. 435 no. 3, pp. 641-658, DOI: 10.1016/j.laa.2010.09.020, 2011.

[27] B. W. Bader, T. G. Kolda, "Algorithm 862: {MATLAB} tensor classes for fast algorithm prototyping," ACM Transactions on Mathematical Software, vol. 32 no. 4, pp. 635-653, DOI: 10.1145/1186785.1186794, 2006.

[28] H. A. L. Kiers, "Towards a standardized notation and terminology in multiway analysis," Journal of Chemometrics, vol. 14 no. 3, pp. 105-122, DOI: 10.1002/1099-128x(200005/06)14:3<105::aid-cem582>3.0.co;2-i, 2000.

[29] L. De Lathauwer, B. De Moor, J. Vandewalle, "A multilinear singular value decomposition," SIAM Journal on Matrix Analysis and Applications, vol. 21 no. 4, pp. 1253-1278, DOI: 10.1137/S0895479896305696, 2000.

[30] D. Liu, W. Li, S.-W. Vong, "The tensor splitting with application to solve multi-linear systems," Journal of Computational and Applied Mathematics, vol. 330, pp. 75-94, DOI: 10.1016/j.cam.2017.08.009, 2018.