2. Minimal PH representation problem: continuous-time case

For the continuous-time case, in view of the work by Asmussen and Bladt [Reference Asmussen, Bladt, Alfa and Chakravarthy29] and O’Cinneide [Reference O’Cinneide30], for any probability distribution $F(\!\cdot\!)$ on $[0, \infty)$ with continuous positive density on $(0, \infty)$ , its LST is given by the rational function

(1) \begin{equation} {\mathcal{L}}(s)=\int^\infty_{0-}\mathrm{e}^{-st}{\sf d} F(t)=\frac{p(s)}{q(s)}=\frac{p_{n-1}s^{n-1}+p_{n-2}s^{n-2}+\cdots+p_0}{s^n+q_{n-1}s^{n-1}+\cdots+q_0}, \quad {\mathrm{Re}}(s)\geq 0.\end{equation}

$F(\!\cdot\!)$ is a PH distribution characterized by the absorbing time of the continuous-time Markov chain on the states ${\mathcal{S}}=\{1,\ldots, m,m+1\}$ with $m\geq n$ , if ${\mathcal{L}}(s)$ satisfies the following:

• (A1) its coefficients are real numbers;

• (A2) p(s) and q(s) are coprime;

• (A3) ${\mathcal{L}}(0)=1$ and ${\mathcal{L}}(s)$ has a unique pole of maximal real part.

Specifically, the Markov chain’s first m states are all transient while the last state is absorbing, the initial state distribution is $\boldsymbol{\alpha}=(\alpha_1,\ldots,\alpha_m)$ with $\alpha{\mathbf{1}}^\top=1$ , and the infinitesimal generator is

\[\left[\begin{array}{c@{\quad}c}A&-A {\mathbf{1}}^\top\\{\mathbf{0}}& 0\end{array}\right],\]

where the $m\times m$ matrix $A=(a_{ij})_{m\times m}$ satisfies $a_{ii}<0$ for $1\leq i\leq m$ , and $a_{ij}\geq 0$ for $i\neq j$ . In particular, $F(t)=1-\alpha \mathrm{e}^{At} {\mathbf{1}}^\top$ for $t\geq 0$ . Then ${\mathcal{L}}(\!\cdot\!)$ can be written as

(2) \begin{equation}{\mathcal{L}}(s)=\int_{0}^{\infty}-\mathrm{e}^{-st}\alpha A\mathrm{e}^{At}{\mathbf{1}}^\top \ {\sf d}t=-\alpha A(sI-A)^{-1}{\mathbf{1}}^\top=\frac{-\alpha A\cdot {\mathrm{adj}}(sI-A){\mathbf{1}}^\top}{{\mathrm{det}}(sI-A)}\end{equation}

for $\mathrm{Re}(s)\geq 0$ . We call such $(\alpha, A)$ an m-order continuous-time PH representation of $F(\!\cdot\!)$ with LST ${\mathcal{L}}(s)$ , and n the algebraic degree of the PH distribution $F(\!\cdot\!)$ . In the following, for the continuous-time case, we start with LST ${\mathcal{L}}(s)$ given in (1) satisfying (A1)–(A3) to identify when there exists such a Markov chain with $(n+1)$ states, that is, when there exists a PH representation with the same order as its algebraic degree. Furthermore, if it exists, we develop an algorithm to determine the initial distribution $\alpha=(\alpha_1,\ldots,\alpha_n)$ and the infinitesimal subgenerator $A=(a_{ij})_{n\times n}$ .

Before developing our analysis on the minimal representation, we comment on ${\mathcal{L}}(s)$ given in (1). In (1), ${\mathcal{L}}(s)$ assumes that there is no probability mass at zero. That is, the Markov chain we identify based on ${\mathcal{L}}(s)$ initially starts at the transient states with probability 1. This assumption can be dropped by directly adding a constant term (the probability of initially starting with the absorbing state) to ${\mathcal{L}}(s)$ . Our following analysis and results on ${\mathcal{L}}(s)$ given in (1) can be directly adopted to the case with positive probability mass at zero by just setting the initial probability at the absorbing state to be this positive probability mass at zero for the Markov chain we have identified. Without loss of generality, we only focus on ${\mathcal{L}}(s)$ given in (1) for the rest of the paper.

3. AM algorithm

Consider the rational function ${\mathcal{L}}(s)$ as given in (1) with (A1)–(A3) holding. Suppose that all distinct roots of q(s) in ${\mathcal{L}}(s)$ are

(3) \begin{equation}\{\lambda_1, \lambda_2,\ldots,\lambda_{\ell}, \mu_{\ell+1}\pm\mathrm{i}\omega_{\ell+1}, \mu_{\ell+2}\pm \mathrm{i}\omega_{\ell+2},\ldots,\mu_{k}\pm \mathrm{i}\omega_{k} \},\end{equation}

with $\lambda_i,\mu_j,\omega_j\in \mathbb{R}$ , $\omega_j\neq 0$ for $1\leq i\leq \ell, \ell+1\leq j \leq k$ , and their multiplicities are $n_1,n_2,\ldots,n_k$ , respectively. That is,

(4) \begin{align}q(s)&=\prod_{i=1}^{\ell}(s-\lambda_{i})^{n_{i}}\cdot\prod_{j=\ell+1}^{k}((s-\mu_{j}-\mathrm{i} \omega_{j})(s-\mu_{j}+\mathrm{i} \omega_{j}))^{n_{j}}\nonumber\\&=\prod_{i=1}^{\ell}(s-\lambda_{i})^{n_{i}}\cdot\prod_{j=\ell+1}^{k}((s-\mu_{j})^2+\omega_{j}^{2})^{n_{j}}.\end{align}

Let $\lambda_1=\arg\max_{1\leq i\leq \ell}\{\lambda_i\}$ . From (A3), we have

\begin{align*}\lambda_{1} > \max_{\substack{2 \leq i \leq \ell \\ \ell + 1 \leq j \leq k}} \{\lambda_{i},\mu_{j}\}.\end{align*}

In view of (2), our problem can be equivalently rephrased as follows. Finding $\alpha \in \mathbb{R}_+^n$ with $\alpha{\mathbf{1}}^\top=1$ , and the non-singular matrix $A=(a_{ij})_{n\times n}\in \mathbb{R}^{n\times n}$ with $a_{ii}<0$ for $1\leq i\leq n$ , $a_{ij}\geq 0$ for $i\neq j$ , and $\smash{\sum_{j=1}^na_{ij}\leq 0}$ for $1\leq i\leq n$ such that

(5) \begin{equation}{\mathcal{L}}(s)=\frac{p(s)}{q(s)}=-\alpha A(sI-A)^{-1}{\mathbf{1}}^\top=\frac{-\alpha A\cdot \mathrm{adj}(sI-A){\mathbf{1}}^\top}{\mathrm{det}(sI-A)}.\end{equation}

Recall (4) and note that $\mathrm{det}(sI-A)$ is a polynomial of nth degree. To ensure (5) holds, a necessary condition is that $q(s)=\mathrm{det}(sI-A)$ . That is, the roots (including multiplicity) of q(s) must be the spectrum of A. Let J be the Jordan form of A. Then there exists a non-singular matrix P such that

(6) \begin{equation}A=P^{-1}JP.\end{equation}

Thus,

(7) \begin{equation}-\alpha A(sI-A)^{-1}{\mathbf{1}}^\top=-\alpha P^{-1} J(sI-J)^{-1}P{\mathbf{1}}^\top.\end{equation}

In the following, based on (6) and (7), we first explore some properties on the Jordan form J and matrix P. Using these properties, we transform our problem into the problem on how to identify whether there exists a feasible solution to the set of linear and quadratic equations. Furthermore, with the help of the non-convex optimization theory, the problem to identify a feasible solution is formulated into a non-convex optimization problem.

3.1. Some properties and an alternative expression

Following the previous discussion, for each real root $\lambda_i$ of q(s), let $J^{(i)}(\lambda_i)$ be the corresponding Jordan block for $1\leq i\leq\ell$ , and for each pair of conjugate complex roots $\mu_{j}\pm\mathrm{i}\omega_{j}$ , let $J^{(j)}(\mu_{j},\omega_{j})$ be the corresponding Jordan block for $\ell+1\leq j\leq k$ . Then

(8) \begin{equation}J={\mathrm{diag}} (J^{(1)}(\lambda_1),\ldots,J^{(\ell)}(\lambda_{\ell}),J^{(\ell+1)}(\mu_{\ell+1},\omega_{\ell+1}),\ldots,J^{(k)}(\mu_{k},\omega_{k})).\end{equation}

For any positive integer m, let

\begin{align*}{\mathcal{J}}_m(\lambda)=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}\lambda& & & \\1 &\ddots & & \\&\ddots &\ddots & \\& &1 &\lambda\end{array}\right]_{m\times m};\end{align*}

\begin{align*}{\mathcal{J}}_m(\mu,\omega)=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}\Theta(\mu,\omega)& & & \\I &\ddots & & \\&\ddots &\ddots & \\& &I &\Theta(\mu,\omega)\end{array}\right]_{(2m)\times (2m)} \quad \mbox{with} \ \Theta(\mu,\omega)=\left[\begin{array}{c@{\quad}c}\mu&-\omega\\\omega& \mu\end{array}\right].\end{align*}

Due to the multiplicity, the form of J is generally not unique. If we consider J with (5) together, however, we are able to prove the uniqueness of the Jordan form J in the following proposition, whose proof is postponed to Appendix A.

Proposition 1. For the rational function $ {\mathcal{L}}(s)$ given in (1), suppose that there exists an nth-order PH representation $(\alpha, A)$ such that (5) holds. Then each distinct eigenvalue of A has exactly one Jordan block, that is, for $1\leq i\leq \ell$ , $J^{(i)}(\lambda_i)={\mathcal{J}}_{n_i}(\lambda_i)$ and, for $\ell+1\leq j\leq k$ , $J^{(j)}(\mu_j,\omega_{j})={\mathcal{J}}_{n_j}(\mu_j,\omega_{j})$ . In other words, the Jordan block has a unique form, i.e.

(9) \begin{equation}{\mathrm{diag}} ( {\mathcal{J}}_{n_1}(\lambda_1),\ldots,{\mathcal{J}}_{n_\ell}(\lambda_{\ell}),{\mathcal{J}}_{n_{\ell+1}}(\mu_{\ell+1},\omega_{\ell+1}),\ldots,{\mathcal{J}}_{n_k}(\mu_{k},\omega_{k}) ),\end{equation}

which is denoted by ${\mathcal{J}}$ in the following analysis.

Next we provide a property of P that helps us to further simplify (7), and whose proof is given in Appendix B.

Proposition 2. For the rational function ${\mathcal{L}}(s)$ given in (1), suppose that there exists an nth-order PH representation $(\alpha,A)$ such that (5) holds. Then there exists a non-singular matrix P such that $A=P^{-1}{\mathcal{J}}P$ and $P{\mathbf{1}}^\top={\mathbf{1}}^\top$ , where ${\mathcal{J}}$ is given in (9) of Proposition 1.

Now we perform a partial-fraction expansion for ${\mathcal{L}}(s)$ according to the factorization of q(s) given in (4). That is,

\begin{equation*}{\mathcal{L}}(s)=\sum_{i=1}^\ell \sum_{r=1}^{n_i}\frac{ c_{i,r}}{(s-\lambda_i)^r}+\sum_{j=\ell+1}^k \sum_{r=1}^{n_j}\biggr(\frac{c_{j,r}}{(s-\mu_j-\mathrm{i}\omega_j)^r}+\frac{\overline{c}_{j,r}}{(s-\mu_j+\mathrm{i}\omega_j)^r}\biggr),\end{equation*}

where $ c_{i,r}\in \mathbb{R}$ and $ c_{i,n_i}\neq 0$ for $1\leq r\leq n_i$ and $1\leq i\leq\ell$ , $ c_{j,r}\in \mathbb{C}$ , $ \overline{c}_{j,r}$ is the conjugate of $ c_{j,r}$ , and $ c_{j,n_j}\neq 0$ for $1\leq r\leq n_j$ and $\ell+1\leq j\leq k$ . Let $\beta$ be the solution to

(10) \begin{align}&\sum_{i=1}^\ell \sum_{r=1}^{n_i}\frac{c_{i,r}}{(s-\lambda_i)^r}+\sum_{j=\ell+1}^k \sum_{r=1}^{n_j} \biggl(\frac{c_{j,r}}{(s-\mu_j-\mathrm{i}\omega_j)^r}+\frac{\overline{c}_{j,r}}{(s-\mu_j+\mathrm{i}\omega_j)^r}\biggr) \nonumber\\& \qquad =-\beta {\mathcal{J}}(sI-{\mathcal{J}})^{-1}{\mathbf{1}}^\top,\end{align}

where ${\mathcal{J}}$ is given in (9) of Proposition 1. To be consistent with the block structure of ${\mathcal{J}}$ , we take the blocks of $\beta$ in the following way:

\begin{gather*}\beta =(\beta(\lambda_1),\ldots,\beta(\lambda_\ell),\beta(\mu_{\ell+1},\omega_{\ell+1}),\ldots,\beta(\mu_k,\omega_k));\\\beta(\lambda_i) =(\beta_1(\lambda_i),\ldots,\beta_{n_i}(\lambda_i)) \quad \mathrm{for } 1\leq i\leq \ell;\\\beta(\mu_{j},\omega_j)=(\beta^{(1)}_1(\mu_{j},\omega_j),\beta^{(1)}_2(\mu_{j},\omega_j),\ldots,\beta^{(n_j)}_1(\mu_{j},\omega_j),\beta^{(n_j)}_2(\mu_{j},\omega_j)), \quad \ell+1\leq j\leq k.\end{gather*}

The following proposition shows that $\beta$ is well defined.

Proposition 3. For the rational function ${\mathcal{L}}(s)$ given in (1), $\beta$ is well defined by (10). Furthermore, $\beta$ is an n-dimensional real vector, $\beta_{n_i}(\lambda_i)\neq 0$ for $1\leq i\leq \ell$ , and $\smash{(\beta^{(n_j)}_{1}(\mu_j,\omega_j))^2}$ $+(\beta^{(n_j)}_2(\mu_j,\omega_j))^2\neq 0$ for $\ell+1 \leq j\leq k$ .

Proof. Note that

(11) \begin{equation}(sI-{\mathcal{J}}_{m}(\lambda))^{-1}=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}(s-\lambda)^{-1}& & & \\(s-\lambda)^{-2}&(s-\lambda)^{-1}& &\\\vdots&\ddots &\ddots &\\(s-\lambda)^{-m}&\cdots&(s-\lambda)^{-2}&(s-\lambda)^{-1}\end{array}\right],\end{equation}

(12) \begin{align}&(sI-{\mathcal{J}}_{m}(\mu,\omega))^{-1}\nonumber\\&\quad=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}(sI-\Theta(\mu,\omega))^{-1}& & &\\(sI-\Theta(\mu,\omega))^{-2}&(sI-\Theta(\mu,\omega))^{-1}& &\\\vdots&\ddots &\ddots &\\(sI-\Theta(\mu,\omega))^{-m}&\cdots&(sI-\Theta(\mu,\omega))^{-2}&(sI-\Theta(\mu,\omega))^{-1}\end{array}\right],\\[6pt]\nonumber\end{align}

where, for any positive integer m, $(sI-\Theta(\mu,\omega))^{-m}=((sI-\Theta(\mu,\omega))^{-1})^m$ and

\begin{equation*}(sI-\Theta(\mu,\omega))^{-1}=\left[\begin{array}{c@{\quad}c}s-\mu &\omega\\-\omega& s-\mu\end{array}\right]^{-1}=\frac{1}{(s-\mu)^2+\omega^2}\left[\begin{array}{c@{\quad}c}s-\mu &-\omega\\\omega& s-\mu\end{array}\right].\end{equation*}

We first prove the part corresponding to the real pole $\lambda_i$ . By (9) and (11),

\[{\mathcal{J}}_{n_i}(\lambda_i)(sI- {\mathcal{J}}_{n_i}(\lambda_i))^{-1}{\mathbf{1}}^\top=\left(\begin{array}{c@{\quad}cc}\frac{\lambda_i}{s-\lambda_i}\\\frac{\lambda_i}{(s-\lambda_i)^2}+\frac{\lambda_i+1}{s-\lambda_i}\\\vdots\\\frac{\lambda_i}{(s-\lambda_i)^{n_i}}+\sum_{r=1}^{n_i-1}\frac{\lambda_i+1}{(s-\lambda_i)^r}\end{array}\right).\]

Then $\beta(\lambda_i)=\big(\beta_{1}(\lambda_i),\ldots,\beta_{n_i}(\lambda_i)\big)$ can be uniquely determined by (10) as follows:

(13) \begin{align}\beta_{n_i}(\lambda_i) = -\frac{ c_{i,n_i}}{\lambda_i}; \end{align}

(14) \begin{align}\beta_{r}(\lambda_i) = -\frac{1}{\lambda_i}\biggl( c_{i,r}+(\lambda_i+1)\sum_{l=1}^{n_i-r}\biggl(-\frac{1}{\lambda_i}\biggr)^l c_{i,r+l}\biggr) \quad\mathrm{for } 1\leq r\leq n_i-1. \end{align}

The proposition for $\beta(\lambda_i)$ with $1\leq i\leq \ell$ follows from (13) and (14).

Now consider the proposition corresponding to the pair of conjugate complex poles $\mu_{j}\pm \mathrm{i}\omega_{j}$ with $\ell+1\leq j\leq k$ . Let

\begin{equation*}U_{n_j}(\text{i},-\text{i})=\frac{1}{\sqrt{2}}\begin{bmatrix} 1&\quad\text{i}& & & & \\ & & \quad\cdots&\quad\cdots& & \\ & & & & \quad1&\quad\text{i}\\ 1&\quad-\text{i}& & & & \\ & &\quad\cdots&\quad\cdots& & \\ & & & & \quad1&\quad-\text{i} \end{bmatrix}_{2n_j\times2n_j}, \end{equation*}

and $U_{n_j}^ H(\mathrm{i},-\mathrm{i})$ denote its conjugate transpose. Note that $U_{n_j}(\mathrm{i},-\mathrm{i})U_{n_j}^ H(\mathrm{i},-\mathrm{i})=I$ and

(15) \begin{align}{\mathcal{J}}_{n_j}(\mu_j,\omega_j)=U_{n_j}^H(\mathrm{i},-\mathrm{i}) \widetilde{\mathcal{J}}_{n_j}(\mu_j,\omega_j)U_{n_j}(\mathrm{i},-\mathrm{i}),\end{align}

where

$$\widetilde{\mathcal{J}}_{n_j}(\mu_j,\omega_j)=\begin{bmatrix}\widetilde{\mathcal{J}}_{n_j}(\mu_j+\mathrm{i}\omega_j)&\\& \widetilde{\mathcal{J}}_{n_j}(\mu_j-\mathrm{i}\omega_j)\end{bmatrix}$$

with $ \widetilde{\mathcal{J}}_{n_j}(\mu_j+\mathrm{i}\omega_j)$ and $\widetilde{\mathcal{J}}_{n_j}(\mu_j-\mathrm{i}\omega_j)$ being given by

\begin{align*}\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}\mu_j+\mathrm{i}\omega_j& & &\\1 &\ddots & &\\&\ddots &\ddots &\\& &1 &\mu_j+\mathrm{i}\omega_j\end{array}\right]_{n_j\times n_j} \quad \mbox{and} \quad\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}\mu_j-\mathrm{i}\omega_j& & &\\1 &\ddots & &\\&\ddots &\ddots &\\& &1 &\mu_j-\mathrm{i}\omega_j\end{array}\right]_{n_j\times n_j},\end{align*}

respectively. From (15),

\begin{align*}&\beta(\mu_j,\omega_j) {\mathcal{J}}_{n_j}(\mu_j,\omega_j) \Big(sI- {\mathcal{J}}_{n_j}(\mu_j,\omega_j)\Big)^{-1}{\mathbf{1}}^\top\\& \quad =\beta(\mu_j,\omega_j)U_{n_j}^ H(\mathrm{i},-\mathrm{i}) \widetilde{\mathcal{J}}_{n_j}(\mu_j,\omega_j)\Big(sI- \widetilde{\mathcal{J}}_{n_j}(\mu_j,\omega_j)\Big)^{-1}U_{n_j}(\mathrm{i},-\mathrm{i}){\mathbf{1}}^\top.\end{align*}

Observe that

\begin{align*}U_{n_j}(\mathrm{i},-\mathrm{i}){\mathbf{1}}^\top=\frac{1}{\sqrt{2}}(\underbrace{1+\mathrm{i},1+\mathrm{i},\ldots,1+\mathrm{i}}_{n_j},\underbrace{1-\mathrm{i},1-\mathrm{i},\ldots,1-\mathrm{i}}_{n_j})^\top.\end{align*}

Then

\begin{align*}& \widetilde{\mathcal{J}}_{n_j}(\mu_j,\omega_j)(sI- \widetilde{\mathcal{J}}_{n_j}(\mu_j,\omega_j))^{-1}U_{n_j}(\mathrm{i},-\mathrm{i}){\mathbf{1}}^\top\\ & \ \quad=\left(\begin{array}{cc}\frac{1+\mathrm{i}}{\sqrt{2}}\widetilde{\mathcal{J}}_{n_j}(\mu_j+\mathrm{i}\omega_j)\Big(sI- \widetilde{\mathcal{J}}_{n_j}(\mu_j+\mathrm{i}\omega_j)\Big)^{-1}{\mathbf{1}}^\top\\\frac{1-\mathrm{i}}{\sqrt{2}}\widetilde{\mathcal{J}}_{n_j}(\mu_j-\mathrm{i}\omega_j)\Big(sI- \widetilde{\mathcal{J}}_{n_j}(\mu_j-\mathrm{i}\omega_j)\Big)^{-1}{\mathbf{1}}^\top\end{array}\right).\end{align*}

Similar to the case of real poles, by (12), we have

\begin{align*}\widetilde{\mathcal{J}}_{n_j}(\mu_j+\mathrm{i}\omega_j)(sI- \widetilde{\mathcal{J}}_{n_j}(\mu_j+\mathrm{i}\omega_j))^{-1}{\mathbf{1}}^\top=\left(\begin{array}{c@{\quad}c@{\quad}c}\frac{\mu_j+ \mathrm{i}\omega_j}{s-(\mu_j+\mathrm{i}\omega_j)}\\\frac{\mu_j+\mathrm{i}\omega_j}{(s-(\mu_j+\mathrm{i}\omega_j))^2}+\frac{1+\mu_j+\mathrm{i}\omega_j}{s-(\mu_j+\mathrm{i}\omega_j)}\\ \vdots\\\frac{\mu_j+\mathrm{i}\omega_j}{(s-(\mu_j+\mathrm{i}\omega_j))^{n_j}}+\sum_{r=1}^{n_j-1}\frac{1+\mu_j+\mathrm{i}\omega_j}{(s-(\mu_j+\mathrm{i}\omega_j))^r}\end{array}\right)\end{align*}

and

\begin{align*}\widetilde{\mathcal{J}}_{n_j}(\mu_j-\mathrm{i}\omega_j)(sI- \widetilde{\mathcal{J}}_{n_j}(\mu_j-\mathrm{i}\omega_j))^{-1}{\mathbf{1}}^\top=\left(\begin{array}{ccc}\frac{\mu_j- \mathrm{i}\omega_j}{s-(\mu_j-\mathrm{i}\omega_j)}\\\frac{\mu_j-\mathrm{i}\omega_j}{(s-(\mu_j-\mathrm{i}\omega_j))^2}+\frac{1+\mu_j-\mathrm{i}\omega_j}{s-(\mu_j-\mathrm{i}\omega_j)}\\ \vdots\\\frac{\mu_j-\mathrm{i}\omega_j}{(s-(\mu_j-\mathrm{i}\omega_j))^{n_j}}+\sum_{r=1}^{n_j-1}\frac{1+\mu_j-\mathrm{i}\omega_j}{(s-(\mu_j-\mathrm{i}\omega_j))^r}\end{array}\right).\end{align*}

Note that

\begin{align*}\beta(\mu_j,\omega_j)U_{n_j}^H(\mathrm{i},-\mathrm{i})=(&\beta^{(1)}_1(\mu_j,\omega_j)-\mathrm{i}\beta^{(1)}_2(\mu_j,\omega_j),\beta^{(2)}_1(\mu_j,\omega_j)-\mathrm{i}\beta^{(2)}_2(\mu_j,\omega_j),\ldots,\\&\beta^{(n_j)}_1(\mu_j,\omega_j)-\mathrm{i} \beta^{(n_j)}_2(\mu_j,\omega_j),\beta^{(1)}_1(\mu_j,\omega_j)+\mathrm{i} \beta^{(1)}_2(\mu_j,\omega_j),\\&\beta^{(2)}_1(\mu_j,\omega_j)+\mathrm{i} \beta^{(2)}_2(\mu_j,\omega_j),\ldots,\beta^{(n_j)}_1(\mu_j,\omega_j)+\mathrm{i}\beta^{(n_j)}_2(\mu_j,\omega_j)).\end{align*}

To prove the proposition, by (9) and with consideration of (13) and (14) for $1\leq i\leq \ell$ , it suffices to prove that, for $j=\ell+1,\ldots,k$ ,

\[\beta(\mu_j,\omega_j)=(\beta^{(1)}_1(\mu_j,\omega_j),\beta^{(1)}_2(\mu_j,\omega_j),\ldots,\beta^{(n_j)}_1(\mu_j,\omega_j),\beta^{(n_j)}_2(\mu_j,\omega_j))\]

is uniquely determined by the following equation:

\begin{align*}&\sum_{r=1}^{n_j} \biggl(\frac{ c_{j,r}}{(s-\mu_j-\mathrm{i}\omega_j)^r} +\frac{\overline{c}_{j,r}}{(s-\mu_j+\mathrm{i}\omega_j)^r}\biggr)\nonumber\\ & \quad = -\beta(\mu_j,\omega_j)U_{n_j}^ H(\mathrm{i},-\mathrm{i})\left(\begin{array}{cc} \frac{1+\mathrm{i}}{\sqrt{2}} \widetilde{\mathcal{J}}_{n_j}(\mu_j+\mathrm{i}\omega_j)(sI- \widetilde{\mathcal{J}}_{n_j}(\mu_j+\mathrm{i}\omega_j))^{-1}{\mathbf{1}}^\top\\[7pt] \frac{1-\mathrm{i}}{\sqrt{2}} \widetilde{\mathcal{J}}_{n_j}(\mu_j-\mathrm{i}\omega_j)(sI- \widetilde{\mathcal{J}}_{n_j}(\mu_j-\mathrm{i}\omega_j))^{-1}{\mathbf{1}}^\top\end{array}\right).\end{align*}

The two terms on the left-hand side of the previous equation can be expressed as

\begin{align*}\sum_{r=1}^{n_j} \frac{ c_{j,r}}{(s-\mu_j-\mathrm{i}\omega_j)^r}\nonumber= &-(\beta^{(1)}_1(\mu_j,\omega_j)-\mathrm{i}\beta^{(1)}_2(\mu_j,\omega_j),\beta^{(2)}_1(\mu_j,\omega_j)-\mathrm{i}\beta^{(2)}_2(\mu_j,\omega_j),\ldots,\\&\quad\beta^{(n_j)}_1(\mu_j,\omega_j)-\mathrm{i} \beta^{(n_j)}_2(\mu_j,\omega_j))\\&\quad \quad \times \frac{1+\mathrm{i}}{\sqrt{2}}\widetilde{\mathcal{J}}_{n_j}(\mu_j+\mathrm{i}\omega_j)(sI- \widetilde{\mathcal{J}}_{n_j}(\mu_j+\mathrm{i}\omega_j))^{-1}{\mathbf{1}}^\top;\\\sum_{r=1}^{n_j} \frac{ \overline{c}_{j,r}}{(s-\mu_j+\mathrm{i}\omega_j)^r}=& -(\beta^{(1)}_1(\mu_j,\omega_j)+\mathrm{i}\beta^{(1)}_2(\mu_j,\omega_j),\beta^{(2)}_1(\mu_j,\omega_j)+\mathrm{i}\beta^{(2)}_2(\mu_j,\omega_j),\ldots,\\&\quad \quad \beta^{(n_j)}_1(\mu_j,\omega_j)+\mathrm{i}\beta^{(n_j)}_2(\mu_j,\omega_j) )\\&\quad \quad \times \frac{1-\mathrm{i}}{\sqrt{2}} \widetilde{\mathcal{J}}_{n_j}(\mu_j-\mathrm{i}\omega_j)(sI- \widetilde{\mathcal{J}}_{n_j}(\mu_j-\mathrm{i}\omega_j))^{-1}{\mathbf{1}}^\top.\end{align*}

Thus, we have

\begin{align*}&\beta^{(n_j)}_1(\mu_j,\omega_j)-\mathrm{i} \beta^{(n_j)}_2(\mu_j,\omega_j)=-\frac{(1-\mathrm{i}) c_{j,n_j}}{\sqrt{2}(\mu_j+\mathrm{i}\omega_j)};\\& \beta^{(r)}_1(\mu_j,\omega_j)-\mathrm{i}\beta^{(r)}_2(\mu_j,\omega_j)\\&\quad =-\frac{1-\mathrm{i}}{\sqrt{2}(\mu_j+\mathrm{i}\omega_j)}\biggl(c_{j,r}+(1+\mu_j+\mathrm{i}\omega_j)\sum_{l=1}^{n_j-r}\biggl(-\frac{1}{\mu_j+\mathrm{i}\omega_j}\biggr)^l c_{j,r+l}\biggr), \quad 1\leq r\leq n_j-1; \end{align*}

\begin{align*}&\beta^{(n_j)}_1(\mu_j,\omega_j)+\mathrm{i}\beta^{(n_j)}_2(\mu_j,\omega_j)=-\frac{(1+\mathrm{i})\overline{c}_{j,n_j}}{\sqrt{2}(\mu_j-\mathrm{i}\omega_j)};\\& \beta^{(r)}_1(\mu_j,\omega_j)+\mathrm{i} \beta^{(r)}_2(\mu_j,\omega_j)\\&\quad =-\frac{1+\mathrm{i}}{\sqrt{2}(\mu_j-\mathrm{i}\omega_j)}\biggl(\overline{c}_{j,r}+(1+\mu_j-\mathrm{i}\omega_j)\sum_{l=1}^{n_j-r}\biggl(-\frac{1}{\mu_j-\mathrm{i}\omega_j} \biggr)^l \overline{c}_{j,r+l}\biggr),\quad 1\leq r\leq n_j-1. \end{align*}

Therefore, $\beta(\mu_j,\omega_j)$ can be uniquely determined as follows: for $ 1\leq r\leq n_j-1$ ,

(16) \begin{align}\beta^{(r)}_1(\mu_j,\omega_j)&=-{\mathrm{Re}}\biggl(\frac{1-\mathrm{i}}{\sqrt{2}(\mu_j+\mathrm{i}\omega_j)}(c_{j,r}+(1+\mu_j+\mathrm{i}\omega_j)\sum_{l=1}^{n_j-r}\biggl(-\frac{1}{\mu_j+\mathrm{i}\omega_j}\biggr)^l c_{j,r+l}\biggr)\biggr), \end{align}

(17) \begin{align}\beta^{(r)}_2(\mu_j,\omega_j)&=\mathrm{Im}\biggl(\frac{1-\mathrm{i}}{\sqrt{2}(\mu_j+\mathrm{i}\omega_j)}(c_{j,r}+(1+\mu_j+\mathrm{i}\omega_j)\sum_{l=1}^{n_j-r}\biggl(-\frac{1}{\mu_j+\mathrm{i}\omega_j}\biggr)^l c_{j,r+l}\biggr)\biggr), \end{align}

(18) \begin{align}\beta^{(n_j)}_1(\mu_j,\omega_j)&=-{\mathrm{Re}}\biggl(\frac{(1-\mathrm{i}) c_{j,n_j}}{\sqrt{2}(\mu_j+\mathrm{i}\omega_j)}\biggr),\end{align}

(19) \begin{align}\beta^{(n_j)}_2(\mu_j,\omega_j)&=\mathrm{Im}\biggl(\frac{(1-\mathrm{i})c_{j,n_j}}{\sqrt{2}(\mu_j+\mathrm{i}\omega_j)}\biggr).\\[6pt]\nonumber\end{align}

The proposition corresponding to the pair of conjugate complex poles $\mu_{j}\pm \mathrm{i}\omega_{j}$ with $\ell+1\leq j\leq k$ follows from (16)–(19). Therefore, $\beta$ is well defined by (10) and is an n-dimensional real vector.

Now we show that $\beta_{n_i}(\lambda_i)\neq 0$ for $1\leq i\leq \ell$ , and $\big(\beta^{(n_j)}_{1}(\mu_j,\omega_j)\big)^2$ $+\big(\beta^{(n_j)}_2(\mu_j,\omega_j)\big)^2\neq 0$ for $\ell+1 \leq j\leq k$ . Note that

\begin{align*}&\beta {\mathcal{J}}(sI-{\mathcal{J}})^{-1}{\mathbf{1}}^\top= \beta (sI-{\mathcal{J}})^{-1}{\mathcal{J}}{\mathbf{1}}^\top\\&\qquad= (\beta(\lambda_1)(sI-{\mathcal{J}}_{n_1}(\lambda_1)\big)^{-1},...,\beta(\lambda_{\ell})(sI- {\mathcal{J}}_{n_{\ell}}(\lambda_{\ell}))^{-1},\\&\hskip11mm\beta(\mu_{\ell+1},\omega_{\ell+1})(sI- {\mathcal{J}}_{n_{\ell+1}}(\mu_{\ell+1},\omega_{\ell+1}))^{-1},\ldots,\beta(\mu_k,\omega_k)(sI- {\mathcal{J}}_{n_k}(\mu_k,\omega_k))^{-1}){\mathcal{J}}{\mathbf{1}}^\top.\end{align*}

Recall (11) and (12), $(s-\lambda_i)^{-n_i}$ only appears in the $n_i$ th row of $(sI- {\mathcal{J}}_{n_i}(\lambda_i))^{-1}$ and $((s-\mu_j)^2+\omega_j^2)^{-n_j}$ only appears in the $(2n_j-1)$ th and $2n_j$ th rows of $(sI-{\mathcal{J}}_{n_j}(\mu_j\pm \mathrm{i} \omega_j))^{-1}$ . By (10), since $ c_{i,n_i}\neq 0$ and $ c_{j,n_j}\neq 0$ , we have $\beta_{n_i}(\lambda_i)\neq 0$ for $1\leq i\leq \ell$ , and $(\beta^{(n_j)}_{1}(\mu_j,\omega_j))^2$ $+(\beta^{(n_j)}_2(\mu_j,\omega_j))^2$ $\neq 0$ for $\ell+1 \leq j\leq k$ . Hence, the proof is complete.

For any rational function $ {\mathcal{L}}(s)$ given in (1), from its poles given in (3), we can determine $ {\mathcal{J}}_{n_i}(\lambda_i)$ for $1\leq i\leq \ell$ , and $ {\mathcal{J}}_{n_j}(\mu_j,\omega_j)$ for $\ell+1\leq j\leq k$ , then J can be constructed as (9). Furthermore, the n-dimensional vector $\beta$ can be determined in (10) from Proposition 3. With the help of Propositions 1 and 2, we can transform the question of finding an nth-order PH representation $(\alpha,A)$ such that (5) holds to the problem of identifying a feasible solution of the set of linear and quadratic equations. To that end, let

(20) \begin{align}\xi=-\biggl(\sum_{i=1}^\ell n_i\lambda_i+\sum_{j=\ell+1}^k2n_j\mu_j\biggr).\end{align}

For the sake of notation simplicity, for any matrix $H=(h_{ij})_{n\times n}\in \mathbb{R}^{n\times n}$ , its following row-blocked form is called a ${\mathcal{J}}$ -blocked partition with respect to the block structure of ${\mathcal{J}}$ given in (9):

\begin{equation*}H=\left[\begin{array}{cccc}H(\lambda_1)\\\vdots\\H(\lambda_\ell)\\H(\mu_{\ell+1},\omega_{\ell+1})\\\vdots\\H(\mu_k,\omega_k)\end{array}\right] \quad \mbox{with} \ H(\lambda_i)=\left[\begin{array}{cccc}H_1(\lambda_i)\\\vdots\\H_{n_i}(\lambda_i)\end{array}\right]\end{equation*}

and

\begin{equation*}H(\mu_j,\omega_j)=\left[\begin{array}{cccc}H^{(1)}(\mu_j,\omega_j)\\\vdots\\H^{(n_j)}(\mu_j,\omega_j)\end{array}\right], \quad H^{(r)}(\mu_j,\omega_j)=\left[\begin{array}{cc}H^{(r)}_1(\mu_j,\omega_j)\\H^{(r)}_2(\mu_j,\omega_j)\end{array}\right], \quad 1\leq r\leq n_j.\end{equation*}

Here $H(\lambda_1)$ is the first $n_i$ rows of H and, in general, $H(\lambda_i)$ is the rows of H through $(n_{[0,i-1]}+1)$ to $n_{[0,i]}$ ; and $H(\mu_j,\omega_j)$ is the rows of H through $(n_{[0,j-1]}+1)$ to $n_{[0, j]}$ with the convention

\[n_{[0,r]}=\left\{\begin{array}{lll}0 &\quad\mbox{if $r=0$},\\n_1+\cdots+n_r &\quad\mbox{if $1\leq r\leq \ell$},\\n_1+\cdots+n_\ell+2(n_{\ell+1}+\cdots+n_{r}) &\quad\mbox{if $\ell+1\leq r\leq k$}.\end{array}\right.\]

Clearly, $H(\lambda_i)\in \mathbb{R}^{n_i\times n}$ for $1\leq i\leq \ell$ , and $H(\mu_j,\omega_j)\in \mathbb{R}^{2n_j\times n}$ for $\ell+1\leq j\leq k$ .

Theorem 1. For the rational function ${\mathcal{L}}(s)$ given in (1), there exists an nth-order continuous-time PH representation $(\alpha,A)$ such that (5) holds if and only if there exist $P=(p_{ij})_{n\times n}$ and $A=(a_{ij})_{n\times n}$ such that the following (in)equalities hold:

(21) \begin{gather}\begin{gathered}PA={\mathcal{J}}P; \quad P {\mathbf{1}}^\top={\mathbf{1}}^\top; \quad \beta P\geq {\mathbf{0}};\quad A {\mathbf{1}}^\top \leq {\mathbf{0}};\\-\xi\leq a_{ii}\leq 0 \quad \mathrm{for } 1\leq i\leq n; \quad 0\leq a_{ij}\leq \xi \quad \mathrm{for } i\neq j, 1\leq i,j\leq n,\end{gathered}\end{gather}

where ${\mathcal{J}}$ is given in (9) of Proposition 1, and $\beta$ is given by Proposition 3. Moreover, the feasible solution (A,P) to (21) directly provides an nth-order continuous-time PH representation $(\alpha,A)$ with $\alpha=\beta P$ .

Proof. First we prove the necessity. From Proposition 2, there exists non-singular P such that the first two equations hold in (21). That is, $A=P^{-1}{\mathcal{J}}P$ and $P {\mathbf{1}}^\top={\mathbf{1}}^\top$ . Thus, from (5),

\begin{align*}{\mathcal{L}}(s)=\frac{p(s)}{q(s)}=-\alpha A(sI-A)^{-1}{\mathbf{1}}^\top= -\alpha P^{-1} {\mathcal{J}}(sI-{\mathcal{J}})^{-1}{\mathbf{1}}^\top.\end{align*}

Then the uniqueness of $\beta$ in Proposition 3 implies that $\alpha P^{-1}=\beta$ . This gives $\beta P=\alpha\geq {\mathbf{0}}$ . Again from $PA={\mathcal{J}}P$ and the non-singularity of P, the spectrum of A is the same as that of ${\mathcal{J}}$ . Then

\begin{align*}\sum_{i=1}^na_{ii}={\sf tr}(A)=\sum_{i=1}^\ell n_i\lambda_i+\sum_{j=\ell+1}^k2n_j\mu_j= -\xi.\end{align*}

As $a_{ii}\leq0$ for each i, we have $-\xi\leq a_{ii}\leq 0$ for $1\leq i\leq n$ . By $A {\mathbf{1}}^\top \leq {\mathbf{0}}$ and $ a_{ij}\geq 0$ for $i\neq j$ , we obtain $a_{ij}\leq-a_{ii}\leq \xi$ for $i\neq j, 1\leq i,j\leq n$ .

Now consider the sufficiency. Make a ${\mathcal{J}}$ -blocked partition for matrix P,

(22) \begin{align}P=\left[\begin{array}{ccc}P(\lambda_1)\\\vdots\\P(\lambda_\ell)\\P(\mu_{\ell+1},\omega_{\ell+1})\\\vdots\\P(\mu_{k},\omega_k)\\\end{array}\right].\end{align}

From the first relation in (21), we have $P(\lambda_i)A= {\mathcal{J}}_{n_i}(\lambda_i)P(\lambda_i)$ and $P(\mu_{j}, \omega_j) A= {\mathcal{J}}_{n_j}(\mu_{j},\omega_j)P(\mu_{j},\omega_j)$ . Recall the form of $ {\mathcal{J}}_{n_i}(\lambda_i)$ . We have

(23) \begin{equation}\begin{gathered}P_{1}(\lambda_i)(A-\lambda_iI)=0,\\P_{r}(\lambda_i)(A-\lambda_iI)=P_{r-1}(\lambda_i), \quad 2\leq r\leq n_i.\end{gathered}\end{equation}

Recall the structure of ${\mathcal{J}}_{n_j}(\mu_{j}, \omega_j)$ . We have

\begin{gather*}P^{(1)}_1(\mu_{j},\omega_j)A=\mu_jP^{(1)}_1(\mu_{j},\omega_j)-\omega_jP^{(1)}_2(\mu_{j},\omega_j),\\P^{(1)}_2(\mu_{j},\omega_j)A=\mu_jP^{(1)}_2(\mu_{j},\omega_j)+\omega_jP^{(1)}_1(\mu_{j},\omega_j),\\P^{(r)}_1(\mu_{j},\omega_j)A=P^{(r-1)}_1(\mu_{j},\omega_j)+\mu_jP^{(r)}_1(\mu_{j},\omega_j)-\omega_jP^{(r)}_2(\mu_{j},\omega_j),\quad2\leq r\leq n_j,\\P^{(r)}_2(\mu_{j},\omega_j)A=P^{(r-1)}_2(\mu_{j},\omega_j)+\mu_jP^{(r)}_2(\mu_{j},\omega_j)+\omega_jP^{(r)}_1(\mu_{j},\omega_j),\quad2\leq r\leq n_j.\end{gather*}

Equivalently,

(24) \begin{align}\begin{split}&(P^{(1)}_1(\mu_{j},\omega_j)+\mathrm{i} P^{(1)}_2(\mu_{j},\omega_j))(A-(\mu_j+\mathrm{i} \omega_j)I)=0,\\&(P^{(1)}_1(\mu_{j},\omega_j)-\mathrm{i}P^{(1)}_2(\mu_{j},\omega_j))(A-(\mu_j-\mathrm{i} \omega_j)I)=0,\\&(P^{(r)}_1(\mu_{j},\omega_j)+\mathrm{i}P^{(r)}_2(\mu_{j},\omega_j))(A-(\mu_j+\mathrm{i} \omega_j)I)\\&\ =P^{(r-1)}_1(\mu_{j},\omega_j)+\mathrm{i}P^{(r-1)}_2(\mu_{j},\omega_j),\quad 2\leq r\leq n_j,\\&(P^{(r)}_1(\mu_{j},\omega_j)-\mathrm{i}P^{(r)}_2(\mu_{j},\omega_j))(A-(\mu_j-\mathrm{i} \omega_j)I)\\&\ =P^{(r-1)}_1(\mu_{j},\omega_j) -\mathrm{i}P^{(r-1)}_2(\mu_{j},\omega_j),\quad 2\leq r\leq n_j.\end{split}\end{align}

It follows from the second equation in (21) that each row vector of P is not zero. By the Exercise in Section 6.3 of [Reference Lancaster and Tismenetsky31], (23) and (24) imply that P is the (generalized) eigenvector matrix of A so it is non-singular. Then the non-singularity of ${\mathcal{J}}$ means that A is non-singular. Let $\alpha=\beta P$ . By the definition of $\beta$ , we have

\begin{align*}{\mathcal{L}}(s)=-\beta {\mathcal{J}}(sI-{\mathcal{J}})^{-1}{\mathbf{1}}^\top=-\beta PP^{-1}{\mathcal{J}}(sI-{\mathcal{J}})^{-1}P{\mathbf{1}}^\top=-\alpha A(sI-A)^{-1}{\mathbf{1}}^\top.\end{align*}

It remains to show that $\alpha{\mathbf{1}}^\top=1$ . This is straightforward by ${\mathcal{L}}(0)=1$ . Therefore, the sufficiency is proved.

With the help of Propositions 1–3, Theorem 1 equivalently transforms the problem of finding the minimal representation $(\alpha, A)$ with the corresponding LST ( ${\mathcal{L}}(s)$ ) given in (1) into the solution existence problem for a set of linear and quadratic equations given in (21). The transformation is based on vector $\beta$ and the Jordan form ${\mathcal{J}}$ . Both of them are uniquely determined from ${\mathcal{L}}(s)$ although ${\mathcal{J}}$ is the Jordan form of A. Telek and Horváth [Reference Telek and Horváth24] also provided a representation transformation but by considering the Jordan form of $(-A)^{-1}$ rather than A itself. We know that the algorithm involving the computation of the matrix inverse usually leads to numerical instability and high computational complexity. Compared with Telek and Horváth [Reference Telek and Horváth24], moreover, the properties of $\beta$ obtained in Proposition 3 can help us to establish the convergence of the AM algorithm proposed in the next subsection for identifying the feasible solution to (21).

Next, we equivalently transform the problem of identifying the existence of a feasible solution to (21) into an optimization problem.

3.2. Non-convex optimization problem

We first formulate condition (21) as a non-convex optimization problem, then use the AM approach developed in non-convex optimization theory to solve it. For the AM approach, see [Reference Bertsekas and Tsitsiklis32].

To solve (21), we reformulate it as an optimization problem:

(25) \begin{eqnarray}{\mathrm{(OP):}}\left\{\begin{array}{ll}\min_{P,A} \|PA-{\mathcal{J}}P \|^2_F \\ \mathrm{such that } \left\{\begin{array}{llll} P {\mathbf{1}}^\top={\mathbf{1}}^\top;\\\beta P\geq {\mathbf{0}};\\A {\mathbf{1}}^\top \leq {\mathbf{0}};\\ -\xi\leq a_{ii}\leq 0 \quad \mbox{for} \ 1\leq i\leq n;\\0\leq a_{ij}\leq \xi \quad\ \mbox{for} \ i\neq j, 1\leq i,j\leq n.\end{array}\right.\end{array}\right. \end{eqnarray}

Note that (21) is solvable if and only if the above optimization problem (OP) has a solution such that the optimal objective value is zero. We use the AM approach to explore the optimal solution of the problem (OP).

Step 1. At the beginning, we pick up $A_0$ and solve the following optimization problem with $A=A_0$ :

(26) \begin{eqnarray}({\mathrm{OP}}[A]): \left \{\begin{array}{ll}\min_{P} \|P A-{\mathcal{J}}P \|^2_F\\ \mathrm{such that } \left\{\begin{array}{ll}P {\mathbf{1}}^\top ={\mathbf{1}}^\top;\\\beta P\geq {\mathbf{0}}.\end{array}\right.\end{array}\right. \end{eqnarray}

Step 2. Let $P_0$ denote the solution of the optimization problem $({\mathrm{OP}}[A])$ given in (26) with $A=A_0$ . Then we solve another optimization problem with $P=P_0$ :

(27) \begin{eqnarray}{\mathrm{(OP[}P\mathrm{])}}: \left \{\begin{array}{ll}\min_{A} \|P A-{\mathcal{J}} P \|^2_F \\\mathrm{such that } \left\{\begin{array}{ll}A {\mathbf{1}}^\top \leq {\mathbf{0}};\\-\xi\leq a_{ii}\leq 0 \quad \mbox{for} \ 1\leq i\leq n;\\0\leq a_{ij}\leq \xi \quad \mbox{for} \ i\neq j, 1\leq i,j\leq n.\end{array}\right.\end{array}\right. \end{eqnarray}

Step 3. Let $A_1$ be the solution to (27). Then solve the optimization problem $({\mathrm{OP}}[A])$ given in (26) by replacing A with $A_1$ , and its solution is denoted by $P_1$ . Again, we solve the optimization problem (OP[P]) given in (27) by replacing P with $P_1$ , and denote its solution by $A_2$ . Repeating step 1 and step 2, we generate a sequence $\{(A_l,P_l): l\geq 0\}$ .

Take P into the ${\mathcal{J}}$ -blocked partition, and rewrite $\|PA-{\mathcal{J}}P\|^2_F$ as

(28) \begin{align}\begin{aligned}\|PA-{\mathcal{J}}P \|^2_F=&\sum_{i=1}^\ell\|P(\lambda_i)A- {\mathcal{J}}_{n_i}(\lambda_i)P(\lambda_i) \|^2_F\\&+\sum_{j=\ell+1}^k\|P(\mu_{j},\omega_j)A- {\mathcal{J}}_{n_j}(\mu_{j},\omega_j)P(\mu_{j},\omega_j) \|^2_F.\end{aligned}\end{align}

To make the above steps workable, we establish the following two lemmas whose proofs are relegated to Appendices C and D. To get an expression for $\|PA-{\mathcal{J}}P\|^2_F$ , let

(29) \begin{gather}\overline{P}(\lambda_i)=(P_1(\lambda_i),P_2(\lambda_i),\ldots,P_{n_i}(\lambda_i));\\[-30pt]\nonumber\end{gather}

(30) \begin{gather}\overline{P}(\mu_{j},\omega_j)=(P^{(1)}_1(\mu_{j},\omega_j),P^{(1)}_2(\mu_{j},\omega_j),\ldots,P^{(n_j)}_1(\mu_{j},\omega_j),P^{(n_j)}_2 (\mu_{j},\omega_j)).\end{gather}

That is, $\overline{P}(\lambda_i)$ is an $(n\times n_i)$ -dimensional row vector containing all rows of $P(\lambda_i)$ . Similarly, $\overline{P}(\mu_{j},\omega_j)$ is an $(n\times 2n_j)$ -dimensional row vector containing all row vectors of $P(\mu_{j},\omega_j)$ . The first lemma is about the convexity on the objective functions in step 1 and step 2.

Lemma 1. For any given $A=(a_{ij})_{n\times n}\in \mathbb{R}^{n\times n}$ , $\|PA-{\mathcal{J}}P \|^2_F$ is convex with respect to $n^2$ variables $\{p_{ij}\colon 1\leq i,j\leq n\}$ given by the entries of matrix $P=(p_{ij})_{n\times n}$ . Similarly, for any given $P=(p_{ij})_{n\times n}\in \mathbb{R}^{n\times n}$ , $\|PA-{\mathcal{J}}P \|^2_F$ is convex with respect to $n^2$ variables $\{a_{ij}\colon 1\leq i,j\leq n\}$ given by the entries of matrix $A=(a_{ij})_{n\times n}$ . Moreover,

(31) \begin{gather}\Big\|P(\lambda_i)A- {\mathcal{J}}_{n_i}(\lambda_i)P(\lambda_i)\Big\|^2_F=\overline{P}(\lambda_i)B_{\lambda_i}B_{\lambda_i}^\top\overline{P}(\lambda_i)^\top,\\[-30pt]\nonumber\end{gather}

(32) \begin{gather}\Big\|P(\mu_{j},\omega_j)A- {\mathcal{J}}_{n_j}(\mu_{j},\omega_j)P(\mu_{j},\omega_j)\Big\|^2_F=\overline{P}(\mu_{j},\omega_j)B_{(\mu_{j},\omega_j)}B_{(\mu_{j},\omega_j)}^\top\overline{P}(\mu_{j},\omega_j)^\top,\end{gather}

where $B_{\lambda_i}$ is an $((nn_i)\times (nn_i))$ matrix, i.e.

(33) \begin{align}B_{\lambda_i}=\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}A-\lambda_iI&-I & &\\& A-\lambda_iI &\ddots &\\& & \ddots & -I\\& & &A-\lambda_iI\end{array}\right]\end{align}

and $B_{(\mu_{j},\omega_j)}$ is a $((2nn_j)\times (2nn_j))$ matrix, i.e.

(34) \begin{align}\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}A_{(\mu_j,\omega_j)}&-I & &\\& A_{(\mu_j,\omega_j)} &\ddots &\\& & \ddots & -I\\& & &A_{(\mu_j,\omega_j)}\end{array}\right]\quad \mathrm{with } A_{(\mu_j,\omega_j)}=\left[\begin{array}{c@{\quad}c}A-\mu_jI&-\omega_jI\\\omega_j I& A-\mu_jI\end{array}\right].\end{align}

The next lemma shows us that both the convex optimization problem $({\mathrm{OP}}[A])$ given in (26) and the convex optimization problem (OP[P]) given in (27) are well posed.

Lemma 2. For any given $A\in \mathbb{R}^{n\times n}$ , the convex optimization problem $({\mathrm{OP}}[A])$ given in (26) has an optimal solution if the set of the feasible solutions to its constraints is nonempty $;$ For any given $P\in \mathbb{R}^{n\times n}$ , the convex optimization problem ( $\mathrm{OP}[P]$ ) given in (27) has an optimal solution.

Note that if $\beta$ , given by Proposition 3, leads to no feasible solution to the convex optimization problem $({\mathrm{OP}}[A])$ given in (26), then the (in)equalities in (21) for such $\beta$ in Theorem 1 have no feasible solution. This means that the associated rational function $ {\mathcal{L}}(s)$ given in (1) has no nth-order continuous-time PH representation.

To verify the AM approach workable to the non-convex optimization problem (OP), following the procedure of the AM, we need to establish the convergence of $\{(A_l,P_l)\colon l\geq 0\}$ and that the corresponding limit is a critical point. Here a critical point for problem $\min\{g(x),x\in X\}$ is a point $x\in X$ such that $\bigtriangledown g(x)^\top(y-x)\geq 0$ for each $y\in X$ , where $\bigtriangledown g(x)$ denotes the gradient of g at x.

To this end, we first establish the boundedness of the optimal solution to $({\mathrm{OP}}[A])$ given in (26) for any given $A \in\mathbb{R}^{n\times n}$ satisfying the constraints in (27).

Proposition 4 For any given $A \in \mathbb{R}^{n\times n}$ satisfying the constraints in (27), the optimal solution to $({\mathrm{OP}}[A])$ given in (26) is uniformly bounded.

Proof. Denote by Sol $({\mathrm{OP}}[A])$ the solution set of $({\mathrm{OP}}[A])$ . We first show that fixed A, Sol $({\mathrm{OP}}[A])$ is bounded. From Lemma 1, $({\mathrm{OP}}[A])$ is a quadratic convex optimization. Thus, its optimal solution set is the same as its local optimal solution set. If Sol $({\mathrm{OP}}[A])$ is unbounded, by Theorem 4.2 in [Reference Lee, Tam and Yen33], there exist $P\in \mathbb{R}^{n\times n}$ and $Q\in \mathbb{R}^{n\times n}$ with $Q\neq {\mathbf{0}}$ such that $\{P+t Q, t\geq0\} \in$ Sol $({\mathrm{OP}}[A])$ . Then P and Q must satisfy

(35) \begin{equation}P {\mathbf{1}}^\top={\mathbf{1}}^\top, \qquad (P+t Q){\mathbf{1}}^\top={\mathbf{1}}^\top,\quad\beta P\geq {\mathbf{0}}, \qquad \beta(P+t Q)\geq {\mathbf{0}},\end{equation}

and

(36) \begin{align}\|PA-{\mathcal{J}}P\|^2_F = \|(P+t Q)A-{\mathcal{J}}(P+t Q)\|^2_F\end{align}

for $t\geq 0$ . Since t is an arbitrary non-negative real number, (35) implies that $Q{\mathbf{1}}^\top={\mathbf{0}}$ and $\beta Q\geq {\mathbf{0}}$ . Thus, we have $\beta Q{\mathbf{1}}^\top=0$ . This together with $\beta Q\geq{\mathbf{0}}$ gives us $\beta Q={\mathbf{0}}$ . We make a ${\mathcal{J}}$ -blocked partition for Q, and similar to (29) and (30), let

\begin{gather*}\overline{Q}(\lambda_i)=(Q_1(\lambda_i),Q_2(\lambda_i),\ldots,Q_{n_i}(\lambda_i));\\\overline{Q}(\mu_{j},\omega_j)=(Q^{(1)}_1(\mu_{j},\omega_j),Q^{(1)}_2(\mu_{j},\omega_j),\ldots,Q^{(n_j)}_1(\mu_{j},\omega_j),Q^{(n_j)}_2(\mu_{j},\omega_j)).\end{gather*}

By (28)–(32), (36) gives, for $t\geq 0$ ,

\begin{align*}&\sum_{i=1}^\ell\overline{P}(\lambda_i)B_{\lambda_i}B_{\lambda_i}^\top\overline{P}(\lambda_i)^\top+\sum_{j=\ell+1}^k\overline{P}(\mu_{j},\omega_j)B_{(\mu_{j},\omega_j)}B_{(\mu_{j},\omega_j)}^\top\overline{P}(\mu_{j},\omega_j)^\top\\&\quad=\sum_{i=1}^\ell(\overline{P}(\lambda_i)+t\overline{Q}(\lambda_i))B_{\lambda_i}B_{\lambda_i}^\top(\overline{P}(\lambda_i)+t\overline{Q}(\lambda_i))^\top\\&\hskip15mm +\sum_{j=\ell+1}^k(\overline{P}(\mu_{j},\omega_j)+t\overline{Q}(\mu_{j},\omega_j))B_{(\mu_{j},\omega_j)}B_{(\mu_{j},\omega_j)}^\top(\overline{P}(\mu_{j},\omega_j)+t\overline{Q}(\mu_{j},\omega_j))^\top.\end{align*}

It can be directly calculated that

\begin{align*}&t^2\biggl(\sum_{i=1}^\ell\overline{Q}(\lambda_i)B_{\lambda_i}B_{\lambda_i}^\top\overline{Q}(\lambda_i)^\top+\sum_{j=\ell+1}^k\overline{Q}(\mu_{j},\omega_j)B_{(\mu_{j},\omega_j)}B_{(\mu_{j},\omega_j)}^\top\overline{Q}(\mu_{j},\omega_j)^\top\biggr)\\&+2t\biggl(\sum_{i=1}^\ell\overline{P}(\lambda_i)B_{\lambda_i}B_{\lambda_i}^\top\overline{Q}(\lambda_i)^\top+\sum_{j=\ell+1}^k\overline{P}(\mu_{j},\omega_j)B_{(\mu_{j},\omega_j)}B_{(\mu_{j},\omega_j)}^\top\overline{Q}(\mu_{j},\omega_j)^\top\biggr)\\&\quad=0\end{align*}

for $t\geq 0$ . Since t is an arbitrary non-negative real number, we have

\begin{align*}\sum_{i=1}^\ell\overline{Q}(\lambda_i)B_{\lambda_i}B_{\lambda_i}^\top\overline{Q}(\lambda_i)^\top+\sum_{j=\ell+1}^k\overline{Q}(\mu_{j},\omega_j)B_{(\mu_{j},\omega_j)}B_{(\mu_{j},\omega_j)}^\top\overline{Q}(\mu_{j},\omega_j)^\top=0.\end{align*}

Note that each term is non-negative. Thus, for $1\leq i\leq \ell$ and $\ell+1\leq j\leq k$ ,

\begin{gather*}\overline{Q}(\lambda_i)B_{\lambda_i}B_{\lambda_i}^\top\overline{Q}(\lambda_i)^\top=0,\\\overline{Q}(\mu_{j},\omega_j)B_{(\mu_{j},\omega_j)}B_{(\mu_{j},\omega_j)}^\top\overline{Q}(\mu_{j},\omega_j)^\top=0.\end{gather*}

Equivalently,

(37) \begin{gather}\overline{Q}(\lambda_i)B_{\lambda_i}={\mathbf{0}},\\[-30pt]\nonumber\end{gather}

(38) \begin{gather}\overline{Q}(\mu_{j},\omega_j)B_{(\mu_{j},\omega_j)}={\mathbf{0}}.\end{gather}

We first focus on (37), by the definition of $B_{\lambda_i}$ in (33),

\begin{align*}\overline{Q}(\lambda_i)B_{\lambda_i}=(Q_1(\lambda_i),Q_2(\lambda_i),\ldots,Q_{n_i}(\lambda_i))\cdot\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}A-\lambda_iI&-I & &\\& A-\lambda_iI &\ddots &\\& & \ddots & -I\\& & &A-\lambda_iI\end{array}\right]= {\mathbf{0}}.\end{align*}

This gives

\begin{gather*}Q_1(\lambda_i)\cdot(A-\lambda_iI)={\mathbf{0}},\\Q_r(\lambda_i)\cdot(A-\lambda_iI)=Q_{r-1}(\lambda_i), \quad 2\leq r\leq n_i.\end{gather*}

Note that if there exists r such that $Q_r(\lambda_i)={\mathbf{0}}$ , then $Q_l(\lambda_i)={\mathbf{0}}$ for all $l\leq r$ . If $Q_{n_i}(\lambda_i)\neq {\mathbf{0}}$ , let $ \sigma_i=\min\{r:Q_r(\lambda_i)\neq {\mathbf{0}}\}$ . Otherwise, let $\sigma_i=n_i+1$ . Hence, when $\sigma_i\leq n_i$ , we have $Q_l(\lambda_i)={\mathbf{0}}$ for all $l< \sigma_i$ and $Q_l(\lambda_i)\neq{\mathbf{0}}$ for all $\sigma_i\leq l\leq n_i$ .

For (38), by the definition of $B_{(\mu_j,\omega_j)}$ in (34),

\begin{align*}\overline{Q}(\mu_{j},\omega_j)B_{(\mu_{j},\omega_j)}&=(Q^{(1)}_1(\mu_{j},\omega_j),Q^{(1)}_2(\mu_{j},\omega_j), \ldots,Q^{(n_j)}_1(\mu_{j},\omega_j)_{(1)},Q^{(n_j)}_2(\mu_{j},\omega_j)_{(2)})\\&\hskip4.5mm\times\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c}A_{(\mu_j,\omega_j)}&-I & &\\& A_{(\mu_j,\omega_j)} &\ddots &\\& & \ddots & -I\\& & &A_{(\mu_j,\omega_j)}\end{array}\right]\\&= {\mathbf{0}}.\end{align*}

Equivalently,

(39a) \begin{gather}Q^{(1)}_1(\mu_{j},\omega_j)\cdot(A-\mu_jI)+\omega_jQ^{(1)}_2(\mu_{j},\omega_j)={\mathbf{0}},\\[-30pt]\nonumber\end{gather}

(39b) \begin{gather}Q^{(1)}_2(\mu_{j},\omega_j)\cdot(A-\mu_jI)-\omega_jQ^{(1)}_1(\mu_{j},\omega_j)={\mathbf{0}},\\[-30pt]\nonumber\end{gather}

(39c) \begin{gather}Q^{(r)}_1(\mu_{j},\omega_j)\cdot(A-\mu_jI)+\omega_jQ^{(r)}_2(\mu_{j},\omega_j) =Q^{(r-1)}_1(\mu_{j},\omega_j),\quad 2\leq r\leq n_j,\\[-30pt]\nonumber\end{gather}

(39d) \begin{gather}Q^{(r)}_2(\mu_{j},\omega_j)\cdot(A-\mu_jI)-\omega_jQ^{(r)}_1(\mu_{j},\omega_j)=Q^{(r-1)}_2(\mu_{j},\omega_j),\quad 2\leq r\leq n_j.\end{gather}

By multiplying (39b) by i, then adding (39a), we obtain the following first equation. By multiplying (39b) by i, then taking it away from (39a), we obtain the following second equation. Similarly, for $2\leq r\leq n_j$ , by multiplying (39d) by i, then adding (39c), we have the following third equation. Multiplying (39d) by i, then taking it away from (39c), we have the following last equation:

(40a) \begin{align}(Q^{(1)}_1(\mu_{j},\omega_j)+\mathrm{i} Q^{(1)}_2(\mu_{j},\omega_j))\cdot(A-\mu_jI-\mathrm{i} \omega_jI)={\mathbf{0}},\\[-30pt]\nonumber\end{align}

(40b) \begin{align}(Q^{(1)}_1(\mu_{j},\omega_j)-\mathrm{i}Q^{(1)}_2(\mu_{j},\omega_j))\cdot(A-\mu_jI+\mathrm{i} \omega_jI)={\mathbf{0}},\\[-30pt]\nonumber\end{align}

(40c) \begin{align}&(Q^{(r)}_1(\mu_{j},\omega_j)+\mathrm{i}Q^{(r)}_2(\mu_{j},\omega_j))\cdot(A-\mu_jI-\mathrm{i} \omega_jI)\\&\qquad=Q^{(r-1)}_1(\mu_{j},\omega_j)+\mathrm{i}Q^{(r-1)}_2(\mu_{j},\omega_j), \quad 2\leq r\leq n_j\nonumber,\\[-30pt]\nonumber\end{align}

(40d) \begin{align}&(Q^{(r)}_1(\mu_{j},\omega_j)-\mathrm{i}Q^{(r)}_2(\mu_{j},\omega_j))\cdot(A-\mu_jI+\mathrm{i} \omega_jI)\\&\qquad=Q^{(r-1)}_1(\mu_{j},\omega_j)-\mathrm{i}Q^{(r-1)}_2(\mu_{j},\omega_j),\quad 2\leq r\leq n_j\nonumber.\end{align}

If $Q^{(n_j)}_l(\mu_{j},\omega_j)\neq {\mathbf{0}}$ for $l=1,2$ , let

$$\zeta_{j}=\min\{r\colon Q^{(r)}_1(\mu_{j},\omega_j)\neq {\mathbf{0}}\},\quad\xi_j=\min\{r\colon Q^{(r)}_2(\mu_{j},\omega_j)\neq {\mathbf{0}}\}.$$

Otherwise, if $\smash{Q^{(n_j)}_1(\mu_{j},\omega_j)= {\mathbf{0}}}$ , let $\zeta_{j}=n_j+1$ , and if $\smash{Q^{(n_j)}_2(\mu_{j},\omega_j)= {\mathbf{0}}}$ , let $\xi_{j}=n_j+1$ . We first claim that $\zeta_{j}=\xi_j$ . If it is not true, assume that $\zeta_{j}>\xi_j$ . Using (40a)–(40d) with $r=\xi_j$ , and by $\smash{Q^{(\xi_j)}_1(\mu_{j},\omega_j)=0}$ as $\zeta_{j}>\xi_j$ ,

\begin{equation*}\mathrm{i} Q^{(\xi_j)}_2(\mu_{j},\omega_j)\cdot(A-\mu_jI-\mathrm{i} \omega_jI)={\mathbf{0}},\qquad -\mathrm{i} Q^{(\xi_j)}_2(\mu_{j},\omega_j)\cdot(A-\mu_jI+\mathrm{i}\omega_jI)={\mathbf{0}}.\end{equation*}

Adding these two equations together, $\omega_j\neq 0$ implies that $Q^{(\xi_j)}_2(\mu_{j},\omega_j)={\mathbf{0}}$ . This contradicts the definition of $\xi_j$ . Conversely, if $\zeta_{j}<\xi_j$ , the proof is similar.

Now we show that $Q^{(r)}_1(\mu_{j},\omega_j)=Q^{(r)}_2(\mu_{j},\omega_j)={\mathbf{0}}$ for $r<\zeta_j$ , and $Q^{(r)}_1(\mu_{j},\omega_j)\neq{\mathbf{0}}$ , $Q^{(r)}_2(\mu_{j},\omega_j)\neq{\mathbf{0}}$ for $\zeta_j\leq r\leq n_j$ . The former is obvious by the definition of $\xi_j$ . We only show the latter. If there exists $\zeta_j< r\leq n_j$ such that $Q^{(r)}_1(\mu_{j},\omega_j)={\mathbf{0}}$ or $Q^{(r)}_2(\mu_{j},\omega_j)={\mathbf{0}}$ , let

\begin{equation*}\zeta_{j}^0=\min\{l>\zeta_j\colon Q^{(l)}_1(\mu_{j},\omega_j)= {\mathbf{0}}\},\quad\xi_j^0=\min\{l>\zeta_j\colon Q^{(l)}_2(\mu_{j},\omega_j)= {\mathbf{0}}\}.\end{equation*}

We first show that $\zeta_{j}^0=\xi_j^0$ . Otherwise, if $\zeta_{j}^0>\xi_j^0$ , from (40a)–(40d),

\begin{align*}&Q^{(\xi_j^0)}_1(\mu_{j},\omega_j)\cdot(A-\mu_jI-\mathrm{i} \omega_jI)=Q^{(\xi_j^0-1)}_1(\mu_{j},\omega_j)+\mathrm{i}Q^{(\xi_j^0-1)}_2(\mu_{j},\omega_j),\\&Q^{(\xi_j^0)}_1(\mu_{j},\omega_j)\cdot(A-\mu_jI+\mathrm{i}\omega_jI)=Q^{(\xi_j^0-1)}_1(\mu_{j},\omega_j)-\mathrm{i}Q^{(\xi_j^0-1)}_2(\mu_{j},\omega_j),\\&(Q^{(\xi_j^0-1)}_1(\mu_{j},\omega_j)+\mathrm{i}Q^{(\xi_j^0-1)}_2(\mu_{j},\omega_j))\cdot(A-\mu_jI-\mathrm{i}\omega_jI)^{\xi_j^0-1-\zeta_j}\\&\qquad=Q^{(\zeta_j)}_1(\mu_{j},\omega_j)+\mathrm{i}Q^{(\zeta_j)}_2(\mu_{j},\omega_j),\\&(Q^{(\xi_j^0-1)}_1(\mu_{j},\omega_j)-\mathrm{i}Q^{(\xi_j^0-1)}_2(\mu_{j},\omega_j))\cdot(A-\mu_jI+\mathrm{i}\omega_jI)^{\xi_j^0-1-\zeta_j}\\&\qquad=Q^{(\zeta_j)}_1(\mu_{j},\omega_j)-\mathrm{i}Q^{(\zeta_j)}_2(\mu_{j},\omega_j),\\&(Q^{(\zeta_j)}_1(\mu_{j},\omega_j)+\mathrm{i}Q^{(\zeta_j)}_2(\mu_{j},\omega_j))\cdot(A-\mu_jI-\mathrm{i} \omega_jI)={\mathbf{0}},\\&(Q^{(\zeta_j)}_1(\mu_{j},\omega_j)-\mathrm{i}Q^{(\zeta_j)}_2(\mu_{j},\omega_j))\cdot(A-\mu_jI+\mathrm{i} \omega_jI)={\mathbf{0}}.\end{align*}

Note that $Q^{(l)}_1(\mu_{j},\omega_j)\neq {\mathbf{0}}$ and $Q^{(l)}_2(\mu_{j},\omega_j)\neq {\mathbf{0}}$ for $\zeta_j\leq l\leq\xi_j^0-1$ . This implies that $Q^{(\xi_j^0)}_1(\mu_{j},\omega_j)$ is the generalized eigenvector corresponding to both $\mu_j+\mathrm{i} \omega_j$ and $\mu_j-\mathrm{i} \omega_j$ . This leads to a contradiction. Conversely, if $\zeta_{j}^0<\xi_j^0$ , the proof is similar.

Now we have $\zeta_{j}^0=\xi_j^0$ . This together with (40a)–(40d) gives

$$Q^{(r)}_1(\mu_{j},\omega_j)=Q^{(r)}_2(\mu_{j},\omega_j)={\mathbf{0}},\quad r<\zeta_j^0.$$

This leads to a contradiction. Then we obtain $Q^{(r)}_1(\mu_{j},\omega_j)\neq{\mathbf{0}}$ , $Q^{(r)}_2(\mu_{j},\omega_j)\neq{\mathbf{0}}$ for $\zeta_j\leq r\leq n_j$ .

Let ${\mathcal{N}}=\{1\leq i \leq \ell\colon \sigma_i\leq n_i\}\cup\{\ell+ 1\leq j\leq k\colon \zeta_j\leq n_j\}.$ It follows from $Q\neq {\mathbf{0}}$ that ${\mathcal{N}}$ is not empty. In summary, for $i,j\in {\mathcal{N}}$ ,

\begin{gather*}Q_{\sigma_i}(\lambda_i)\cdot(A-\lambda_iI)={\mathbf{0}},\\Q_r(\lambda_i)\cdot(A-\lambda_iI)=Q_{r-1}(\lambda_i), \quad \sigma_i+1\leq r\leq n_i\end{gather*}

and

\begin{align*}&(Q^{(\zeta_j)}_1(\mu_{j},\omega_j)+\mathrm{i} Q^{(\zeta_j)}_2(\mu_{j},\omega_j))\cdot(A-\mu_jI-\mathrm{i} \omega_jI)={\mathbf{0}},\\&(Q^{(\zeta_j)}_1(\mu_{j},\omega_j)-\mathrm{i}Q^{(\zeta_j)}_2(\mu_{j},\omega_j))\cdot(A-\mu_jI+\mathrm{i} \omega_jI)={\mathbf{0}},\\&(Q^{(r)}_1(\mu_{j},\omega_j)+\mathrm{i}Q^{(r)}_2(\mu_{j},\omega_j))\cdot(A-\mu_jI-\mathrm{i} \omega_jI)\\&\qquad=Q^{(r-1)}_1(\mu_{j},\omega_j)+\mathrm{i} Q^{(r-1)}_2(\mu_{j},\omega_j),\quad\zeta_j+1\leq r\leq n_j,\\&(Q^{(r)}_1(\mu_{j},\omega_j)-\mathrm{i}Q^{(r)}_2(\mu_{j},\omega_j))\cdot(A-\mu_jI+\mathrm{i} \omega_jI)\\&\qquad=Q^{(r-1)}_1(\mu_{j},\omega_j) -\mathrm{i} Q^{(r-1)}_2(\mu_{j},\omega_j) ,\quad \zeta_j+1\leq r\leq n_j.\end{align*}

It follows that, for matrix A, $Q_{\sigma_i}(\lambda_i)$ is the eigenvector related to $\lambda_i$ and $Q_r(\lambda_i)$ is the generalized eigenvector related to $\lambda_i$ for $\sigma_i+1\leq r\leq n_i$ ; $Q^{(\zeta_j)}_1(\mu_{j},\omega_j)\pm \mathrm{i}Q^{(\zeta_j)}_2(\mu_{j},\omega_j)$ are the eigenvectors related to $\mu_j\pm \mathrm{i}\omega_j$ and $Q^{(r)}_1(\mu_{j},\omega_j)\pm \mathrm{i}Q^{(r)}_2(\mu_{j},\omega_j)$ are the generalized eigenvectors related to $\mu_j\pm \mathrm{i} \omega_j$ for $\zeta_j+1\leq r\leq n_j$ . Thus, they are linearly independent. This contradicts $\beta Q={\mathbf{0}}$ from Proposition 3. Therefore, such Q does not exist. This implies that the optimal solution set to $({\mathrm{OP}}[A])$ is bounded.

Now we show that Sol $({\mathrm{OP}}[A])$ is uniformly bounded for A satisfying constraints in (27). If it is not true, there must exist $\{A_l\colon l\geq 1\}$ such that

(41) \begin{align}{\sup}_{P_l\in {\mathrm{Sol}} \mbox{{(OP[$A_l$])} }} \|P_l \|_F^2\rightarrow \infty.\end{align}

Note that $\{A_l\colon l\geq 1\}$ belongs to a compact set. Without loss of generality, we assume that $\{A_l\colon l\geq 1\}$ converges to $A_\infty$ . By (28)–(32), for each (OP[ $A_l$ ]) with $1\leq l\leq \infty$ , its objective function is quadratic and the corresponding Hessian matrix is positive semi-definite. Consider the following affine variational inequality problem denoted by (AVI[ $A_l$ ]):

\begin{eqnarray*}\left\{\begin{array}{lll}&\mbox{Find $P\in \mathbb{R}^{n\times n}$ such that}\\&\mbox{for all $U\in \mathbb{R}^{n\times n}$ with $U{\mathbf{1}}^{\top}={\mathbf{1}}^{\top}$ and $\beta U\geq {\mathbf{0}}$, we have}\\&\sum_{i=1}^\ell\langle B_{\lambda_i}B^\top_{\lambda_i}\overline{P}^\top(\lambda_i),\overline{U}^\top(\lambda_i)-\overline{P}^\top(\lambda_i)\rangle\\&+\sum_{j=\ell+1}^k\langle B_{(\mu_{j},\omega_j)}B^\top_{(\mu_{j},\omega_j)}\overline{P}^\top(\mu_{j},\omega_j),\overline{U}^\top(\mu_{j},\omega_j)-\overline{P}^\top(\mu_{j},\omega_j)\rangle\geq 0,\\&\mathrm{where}\ \overline{P}(\lambda_i), \overline{U}(\lambda_i), \overline{P}(\mu_{j},\omega_j)\ \mathrm{and}\ \overline{U}(\mu_{j},\omega_j)\ \mathrm{are \, given \, by \, the}\ {\mathcal{J}}\mathrm{-blocked}\\&\mathrm{ partitions \, of}\ P\ \mathrm{and}\ U\ \mathrm{as \, ({29}) \, and \, ({30}), and}\ B_{\lambda_i}\ \mathrm{and}\ B_{(\mu_j,\omega_j)}\ \mathrm{are \, given \, by}\\ &\mathrm{({33}) \, and \, ({34}) by \, replacing}\ A\ \mathrm{with}\ A_l.\end{array}\right.\end{eqnarray*}

Using Corollary 5.2 in [Reference Lee, Tam and Yen33], the solution set of (AVI[ $A_l$ ]) is the same as the Karush–Kuhn–Tucker (KKT) solution set of problem (OP[ $A_l$ ]), and so it is the same as the optimal solution set of problem (OP[ $A_l$ ]) as it is convex. However, by Theorem 7.6 in [Reference Lee, Tam and Yen33], we know the solution sets of (AVI[ $A_l$ ]) for $l\geq 1$ are uniformly bounded. Therefore, the optimal solution sets for Sol (OP[ $A_l$ ]) for $l\geq 1$ are uniformly bounded. This implies that (41) is not true. Thus, we complete the proof of this proposition.

Now we can obtain the effectiveness of the AM for the non-convex optimization problem (OP) given in (25). For simplicity, let ${\mathcal{F}} (A,P)=\|PA-{\mathcal{J}}P \|^2_F$ .

Theorem 2. The algorithm (given in step 1 and step 2) can guarantee the following.

1. (i) The sequence $\{{\mathcal{F}}(A_l,P_l)\colon l\geq 0\}$ has a finite limit denoted by ${\mathcal{F}}_\infty$ .

2. (ii) The sequence $\{(A_l, P_l)\colon l\geq 0\}$ admits at least one limit point, and its any limit point, say $( {A}_\infty,{P}_\infty)$ , is a critical point of the non-convex optimization problem $(\mathrm{OP})$ given in (25). Furthermore, ${\mathcal{F}}( {A}_\infty,{P}_\infty)= {\mathcal{F}}_\infty$ .

Proof. We first prove that $\{{\mathcal{F}}(A_l, P_l)\colon l\geq 0\}$ is monotonically decreasing. Recalling the procedure to generate $\{(A_l, P_l)\colon l\geq 0\}$ in step 1 and step 2, we first arbitrarily pick up $A_0$ that satisfies the constraints in (27), and that:

• • $P_l$ is a solution for the convex optimization problem (OP[A]) given in (26) with $A=A_{l}$ for $l\geq0$ ;

• • $A_l$ is a solution for the convex optimization problem (OP[P]) given in (27) with $P=P_{l-1}$ for $l\geq 1.$

Hence, we have, for $l\geq 0$ ,

\begin{equation*}{\mathcal{F}}(A_l, P_l)\geq {\mathcal{F}}(A_{l+1}, P_l) \geq {\mathcal{F}}(A_{l+1},P_{l+1}).\end{equation*}

Since ${\mathcal{F}}(A_l, P_l) \geq 0$ for each l, there exists a non-negative real number $ {\mathcal{F}}_\infty$ such that

\begin{equation*}\lim_{l \rightarrow \infty} {\mathcal{F}}(A_l, P_l) ={\mathcal{F}}_\infty.\end{equation*}

This gives the first part of the theorem.

Now we prove the second part of the theorem. The existence of the limit point admitted by the sequence $\{(A_l, P_l)\colon l\geq 0\}$ follows directly from Proposition 4. Furthermore, from the first part of the theorem, we know that the objective value given by any limit point of $\{(A_l, P_l)\colon l\geq 0\}$ must be ${\mathcal{F}}_\infty$ . Thus, the remainder of the proof of the second part is to show that any limit point of $\{(A_l, P_l)\colon l\geq 0\}$ must be a critical point of the non-convex optimization problem $(\mathrm{OP})$ given in (25). From Proposition 4, this follows from Corollary 2 in [Reference Grippo and Sciandrone34]. Then we complete this theorem.

Telek and Horváth [Reference Telek and Horváth24] developed the MM method to determine an nth-order PH representation of the rational function ${\mathcal{L}}(s)$ given in (1). To make their algorithm possibly converge to a PH distribution, continuously differentiable functions are used to penalize the negative off-diagonal terms that must be properly constructed and implemented in each step of the algorithm. When their algorithm converges by keeping properly tuning the penalization functions, if it terminates with a PH representation, then it successfully finds an nth-order PH representation of ${\mathcal{L}}(s)$ . Otherwise, if it terminates without a PH representation, it cannot tell us whether ${\mathcal{L}}(s)$ has an nth-order PH representation (a more detailed discussion is given in Section 5). Here, Theorem 2, however, shows that the AM algorithm converges, and the convergence limit point is a critical point without keeping properly constructing the penalization functions. Similar to the MM method of [Reference Telek and Horváth24], on the other hand, as the outcome objective value obtained from the AM algorithm may not be the global minimum, when the outcome objective value computed by the AM algorithm turns out to be non-zero, it does not imply that ${\mathcal{L}}(s)$ has no nth-order PH representation; see Example 3 in Section 5.

Next we show that the previous analysis and results for the continuous-time minimal PH representation problem can be used to analyse the discrete-time case.

4. Minimal discrete-time PH representation

Parallel to the continuous-time case, consider a discrete non-negative-integer-value random variable with the generating function

(42) \begin{align}{\mathcal{G}}(z)=\frac{\widetilde{p}(z)}{\widetilde{q}(z)}=\frac{\widetilde{p}_nz^n+\widetilde{p}_{n-1}z^{n-1}+\cdots+\widetilde{p}_1z}{\widetilde{q}_nz^n+\widetilde{q}_{n-1}z^{n-1}+\cdots+1}\quad \mbox{for $|z|\leq 1$}.\end{align}

From Theorem 1.2 of [Reference O’Cinneide30], we know that this discrete non-negative-integer-value random variable is given by an absorbing time of a discrete-time Markov chain if ${\mathcal{G}}(z)$ satisfies the following:

• ( $\widetilde{A1}$ ) its coefficients are real numbers;

• ( $\widetilde{A2}$ ) $\widetilde{p}_n^2+\widetilde{q}_n^2\neq 0$ , $\widetilde{p}(s)$ and $\widetilde{q}(s)$ are coprime;

• ( $\widetilde{A3}$ ) ${\mathcal{G}}(1)=1$ and the root with minimal modulus of $\widetilde{q}(z)$ is real and strictly greater than 1.

Similar to the continuous-time case, starting with the moment generating function ${\mathcal {G}}(z)$ given in (42) with ( $\widetilde{A1}$ )–( $\widetilde{A3}$ ) holding, our problem can be phrased as finding $\widetilde{\alpha} \in \mathbb{R}_+^n$ with $\widetilde\alpha{\mathbf{1}}^\top=1$ , and a substochastic matrix $\widetilde{A}=(\widetilde{a}_{ij})_{n\times n}\in\mathbb{R}_+^{n\times n}$ with $(I-\widetilde{A})$ being non-singular such that

(43) \begin{align}{\mathcal{G}}(z)&=\sum_{l=1}^{\infty}\widetilde\alpha\widetilde{A}^{l-1}(I-\widetilde{A}){\mathbf{1}}^\top z^l\nonumber\\&=z\widetilde\alpha (I-\widetilde{A})(I-z\widetilde{A})^{-1}{\mathbf{1}}^\top\nonumber\\&=\frac{z\widetilde\alpha (I-\widetilde{A})\cdot \mathrm{adj}(I-z\widetilde{A}){\mathbf{1}}^\top}{\mathrm{det}(I-z\widetilde{A})} \quad\mbox{for $|z|\leq 1$}. \end{align}

Note that the form of (42) with ( $\widetilde{A1}$ )–( $\widetilde{A3}$ ) holding actually includes the generating functions of all discrete-time PH distributions that assign no mass at zero and have no periodicity. Here if a discrete-time PH distribution has no periodicity it means that it has a representation whose corresponding Markov chain has no periodic states. To see this, suppose that $(\widetilde\alpha,\widetilde{A})$ is such a representation from an aperiodic and $(m+1)$ -state Markov chain in which (i) all the states are transient except one state is absorbing, (ii) the initial state distribution given by $\widetilde\alpha$ has zero mass at the absorbing state, and (iii) the transition probabilities among the transient states is $\widetilde{A}$ . We claim that $z\widetilde\alpha (I-\widetilde{A})(I-z\widetilde{A})^{-1}{\mathbf{1}}^\top$ can be written in the form of (42) and satisfies ( $\widetilde{A1}$ )–( $\widetilde{A3}$ ). Note first that the numerator polynomial in $z\widetilde\alpha(I-\widetilde{A})(I-z\widetilde{A})^{-1}{\mathbf{1}}^\top$ has a factor z but the denominator polynomial does not have it. Therefore, after some simplification, when the numerator polynomial and denominator polynomial become coprime, the numerator polynomial has no constant term and so $z\widetilde\alpha (I-\widetilde{A})(I-z\widetilde{A})^{-1}{\mathbf{1}}^\top$ has the form of (42). Furthermore, note that

\[(z\widetilde\alpha (I-\widetilde{A})(I-z\widetilde{A})^{-1} {\mathbf{1}}^\top)_{z=1}=1.\]

Hence, if $z\widetilde\alpha (I-\widetilde{A})(I-z\widetilde{A})^{-1}{\mathbf{1}}^\top$ is a polynomial, then it satisfies ( $\widetilde{A1}$ )–( $\widetilde{A3}$ ). Otherwise, it is straightforward to see that $z\widetilde\alpha (I-\widetilde{A})(I-z\widetilde{A})^{-1}{\mathbf{1}}^\top$ satisfies ( $\widetilde{A1}$ ) and ( $\widetilde{A2}$ ). To see ( $\widetilde{A3}$ ), it suffices to show that its pole with minimal modulus is real and strictly larger than 1. From the proof of the necessity in Theorem 1.2 of [Reference O’Cinneide30, p.50], its pole with minimal modulus is real and positive. To see that it is strictly larger than 1, note that its poles belong to the set of reciprocals of the non-zero eigenvalues of $\widetilde{A}$ . As $\widetilde{A}$ is substochastic and $(I-\widetilde{A})$ is invertible, the spectral radius of $\widetilde{A}$ is strictly less than 1. Hence, the pole of $z\widetilde\alpha (I-\widetilde{A})(I-z\widetilde{A})^{-1}{\mathbf{1}}^\top$ with minimal modulus is strictly larger than 1, and ( $\widetilde{A3}$ ) is satisfied.

Although ${\mathcal{G}}(z)$ in (42) assumes that there is no probability mass at zero, similar to the continuous-time case ${\mathcal{L}}(s)$ given in (1), our following analysis and results on ${\mathcal{G}}(z)$ in (42) can be directly adopted to the case with positive probability mass at zero just by setting the initial probability at the absorbing state to be this positive probability mass at zero for the Markov chain we identified. Without loss of generality, here we only focus on ${\mathcal{G}}(z)$ in (42).

First, we establish the following theorem to transform this discrete-time PH representation problem into a continuous-time PH representation problem.

Proof. Suppose that ${\mathcal{G}}(z)$ has an nth-order discrete-time PH representation $(\widetilde\alpha,\widetilde{A})$ , by (43),

\begin{align*}{\mathcal{G}}(z)=z\widetilde\alpha(I-\widetilde{A}) (I-z\widetilde{A})^{-1}{\mathbf{1}}^\top=\widetilde\alpha(I-\widetilde{A}) (z^{-1}I-\widetilde{A})^{-1}{\mathbf{1}}^\top.\end{align*}

Note that both sides of this equation are analytical in the region $\{z\in\mathbb{C}\colon|z|\leq 1\}$ and when $s\in \mathbb{C}$ and ${\mathrm{Re}}(s)\geq 0$ , $|1/(s+1)|\leq 1$ . Then, for $s\in \mathbb{C}$ and ${\mathrm{Re}}(s)\geq 0$ ,

\begin{align*}{\mathcal{G}}\biggl(\frac{1}{s+1}\biggr)=\widetilde\alpha(I-\widetilde{A})((s+1)I-\widetilde{A})^{-1}{\mathbf{1}}^\top=-\widetilde \alpha(\widetilde{A}-I)(sI-(\widetilde{A}-I))^{-1}{\mathbf{1}}^\top.\end{align*}

Let $A=\widetilde{A}-I$ and write $A=(a_{ij})_{n\times n}$ . In view of the property of $\widetilde{A}$ , we know that A is non-singular, $A{\mathbf{1}}^\top\leq {\mathbf{0}}$ , $a_{ii}\geq -1$ , and $a_{ij}\geq 0$ for $i\neq j,1\leq i,j\leq n$ . Therefore, $(\widetilde \alpha, A)$ is an nth-order continuous-time PH representation.

Conversely, suppose that ${\mathcal{G}}(1/(s+1))$ has an nth-order continuous-time representation $(\alpha,A)$ with the diagonal entries of A larger than or equal to $-1$ . That is,

\begin{align*}{\mathcal{G}}\biggl(\frac{1}{s+1}\biggr)=-\alpha A (sI-A)^{-1}{\mathbf{1}}^\top=-\alpha A((s+1)I-(A+I))^{-1}{\mathbf{1}}^\top\end{align*}

for $s\in \mathbb{C}$ and ${\mathrm{Re}}(s)\geq 0$ . Then, for $z\in \{z\in \mathbb{C}\colon{\mathrm{Re}}(z)>0,|z|\leq 1\}$ ,

\begin{align*}{\mathcal{G}} (z)=-\alpha A (z^{-1}I-(A+I))^{-1}{\mathbf{1}}^\top.\end{align*}

Observe that both sides of the above equation are analytical in the region $\{z\in \mathbb{C}\colon |z|\leq 1\}$ . Using the uniqueness of the analytical function, the equation also holds in the region $\{z\in\mathbb{C}\colon |z|\leq 1\}$ . Let $\widetilde{A}=A+I$ . We have

\begin{align*}{\mathcal{G}}(z)=\alpha(I-\widetilde{A}) (z^{-1}I-\widetilde{A})^{-1} {\mathbf{1}}^\top=z\alpha(I-\widetilde{A}) (I-z\widetilde{A})^{-1}{\mathbf{1}}^\top.\end{align*}

In view of the property of A and $a_{ii}\geq -1$ , we know $(\alpha,\widetilde{A})$ is an nth-order discrete-time PH representation.

Maier [Reference Maier35] showed that ${\mathcal{G}}(c/(s+c))$ with any $c>0$ is the LST of a continuous-time PH distribution if ${\mathcal{G}}(z)$ is the generating function of a discrete-time PH distribution. This result does not show (i) the minimal representation order relation between these two continuous-time and discrete-time PH distributions and (ii) whether ${\mathcal{G}}(z)$ is the generating function of a discrete-time PH distribution given that ${\mathcal {G}}(c/(s+c))$ with any $c>0$ is the LST of a continuous-time PH distribution. To completely adopt the results developed in Section 3, we need the equivalent relation given by Theorem 3 between the minimal representations for the continuous-time and discrete-time PH distributions. To answer the minimal representation problem for the discrete-time PH distribution specified by (43), the equivalent relation obtained by Theorem 3 implies that it is sufficient to consider the continuous-time case with the LST given by ${\mathcal{G}}(1/(s+1))$ . To apply what we have developed for the continuous-time case in Section 3, we have to show that ${\mathcal{G}}(1/(s+1))$ has the canonical form of (1).

Proof. Items (i) and (ii) are straightforward by replacing z with $1/(s+1)$ in (42). Note that ${\mathcal{L}}_0(0)={\mathcal{G}}(1)=1$ . Furthermore, we can factor $\widetilde{q}(z)$ as the following form:

\begin{align*}\widetilde{q}(z)=(1-\gamma_1z)^{n_1}(1-\gamma_2z)^{n_2}\cdots(1-\gamma_bz)^{n_b}.\end{align*}

Here $\gamma_i\neq 0$ . Let $n_{[1,b]}=\sum_{l=1}^bn_l$ . If $\widetilde{q}_n\neq 0$ , then $n_{[1,b]}=n$ ; otherwise, $n_{[1,b]}<n$ . That is, the set of all distinct roots of $\widetilde{q}(z)$ is $\{1/\gamma_1,\ldots,1/\gamma_b\}$ with $1>\gamma_1>|\gamma_i|$ for $2\leq i\leq b$ from ( $\widetilde{A3}$ ). Then

\begin{align*}q_0(s)=(s+1)^{n- n_{[1,b]}}(s-(\gamma_1-1))^{n_1}(s-(\gamma_2-1))^{n_2}\cdots(s-(\gamma_b-1))^{n_b}.\end{align*}

It follows that $\{\gamma_1-1,\gamma_2-1,\ldots,\gamma_b-1\}$ are all distinct roots of q(s), and $-1$ is also a root when $n_{[1,b]}<n$ . Note that ${\mathrm{Re}}(\gamma_i-1)<\gamma_1-1<0$ for $2\leq i\leq b$ and $\gamma_1-1>-1$ . We find that the root with maximal real part of $q_0(s)$ is real and negative. Hence, we have part (iv).

Now it remains to show part (iii). If it is not true, suppose that $\gamma_i-1$ is also a root of $p_0(s)$ . That is,

\[0=p_0(\gamma_i-1)=\gamma_i^n\widetilde{p}\biggl(\frac{1}{\gamma_i}\biggr).\]

Here $\gamma_i\neq 0$ implies that $\widetilde{p}(1/\gamma_i)=0$ . This contradicts the fact that $\widetilde{p}(z)$ and $\widetilde{q}(z)$ are coprime from ( $\widetilde{A2} $ ). Then $\gamma_i-1$ cannot be the root of $p_0(s)$ . Thus, if $n_{[1,b]}=n$ , we complete the proof of part (iii).

If $n_{[1,b]}<n$ , we need to show that $-1$ is not the root of $p_0(s)$ . Note that $n_{[1,b]}<n$ means that $\widetilde{q}_n= 0$ . This tells us that $\widetilde{p}_n\neq 0$ . Then $-1$ is not the root of $p_0(s)$ and part (iii) still holds. Hence, part (iii) is proved.

With the help of Theorems 3 and 4, the discrete-time case with ${\mathcal{G}}(z)$ given in (42) is equivalently transformed into the continuous-time case with ${\mathcal{L}}_0(s)$ , which has the standard form defined by (1) and satisfies (A1)–(A3). Thus, the discrete-time case can be solved by formulating the optimization problem (OP) given in (25) by modifying $\xi$ given in (20) to be 1.

5. Numerical experiments

In this section, we use the AM algorithm to conduct the numerical experiments and report on the corresponding numerical performance of the minimal representations for both continuous-time and discrete-time cases. The numerical experiments lead to the following conclusions: (i) the AM algorithm produces A and P such that $ {\mathcal{F}} (A,P)=0$ and they satisfy the constraints in the non-convex optimization problem (OP) given in (25), then we find a minimal PH representation; (ii) when the AM algorithm produces a non-zero value, it often does not have a minimal PH representation with the same order as its algebraic degree. Thus, the AM algorithm may provide an efficient approach to simultaneously solve the problem of whether there exists an algebraic degree minimal PH representation, and the problem of how to find the corresponding representation for both continuous-time and discrete-time cases; (iii) choosing different initial $A_0$ , we may obtain different minimal PH representations; and (iv) the AM algorithm becomes more efficient and effective when ${\mathcal{L}}(s)$ only has real poles in terms of the computation time and whether a minimal PH representation can be successfully found. At the same time, we make a comparison between the AM algorithm and the MM method in [Reference Mészáros, Horváth and Telek23, Reference Telek and Horváth24] by extensive numerical experiments.

First, we provide an example to show how to adjust the initial $A_0$ to generate different minimal PH representations.

Example 1. Consider

\begin{align*}{\mathcal{L}}(s)=\frac{1.161}{s+1}-\frac{0.23}{s+2.8+0.4\mathrm{i}}-\frac{0.23}{s+2.8-0.4\mathrm{i}}.\end{align*}

If it has a third-order PH representation $(\alpha, A)$ , the spectrum of A must be $\{-1,-2.8\pm 0.4\mathrm{i}\}$ . Then, by Proposition 1, the Jordan form ${\mathcal{J}}$ of A is

\begin{align*}\begin{bmatrix}-1&0 & 0\\0&-2.8&0.4\\0&-0.4&-2.8\\\end{bmatrix}.\end{align*}

From (10), (20) and Theorem 1, we can calculate that

\begin{align*}\beta=(1.161,-0.092,-0.069),\end{align*}

and $\xi=6.6$ . Now we use the AM algorithm to solve the optimization problem ( $\mathrm{OP}$ ) given in (25). Let $A_0={\mathcal{J}}$ . The outcome of the AM algorithm is

\begin{align*}P_{\infty} \approx\left[\begin{array}{c@{\quad}c@{\quad}c}0.5347 & 0.1988 & 0.2665\\-0.1520 & 0.2325 & 0.9195\\-0.1251 & -0.0858 & 1.2109\end{array}\right], \quad A_{\infty} \approx\left[\begin{array}{c@{\quad}c@{\quad}c}-1.1120 & 0.5827 & 0.0033\\0.0000 & -2.5666 & 2.5666\\0.2246 & 0.0000 & -2.9216\end{array}\right],\end{align*}

and ${\mathcal{F}}(A_{\infty},P_{\infty}) \approx 3.6262\times 10^{-9}$ . By Theorem 1, let $\alpha=\beta P_{\infty} \approx(0.6434,0.2154,0.1412)$ . We obtain a minimal PH representation $(\alpha,A_{\infty})$ of ${\mathcal{L}}(s)$ .

If we choose

\begin{align*}\widehat{A}_0=\left[\begin{array}{c@{\quad}c@{\quad}c}-1&1&0\\0&0&0\\0&0&0\end{array}\right],\end{align*}

the outcome of the AM algorithm is

\begin{align*}\widehat{P}_{\infty} \approx\left[\begin{array}{c@{\quad}c@{\quad}c}0.3438 & 0.3351 & 0.3211\\0.3568 & 0.8141 &-0.1709\\0.0427 & 1.0944 &-0.1371\end{array}\right], \quad\widehat{A}_{\infty} \approx\left[\begin{array}{c@{\quad}c@{\quad}c}-2.1582 & 2.0974 & 0.0000\\0.0000 &-3.1762 & 0.2547\\1.2403 & 0.0253 &-1.2658\end{array}\right],\end{align*}

and ${\mathcal{F}}(\widehat{A}_{\infty},\widehat{P}_{\infty}) \approx 4.0136\times 10^{-9}$ . Let $\widehat{\alpha}=\beta \widehat{P}_{\infty} \approx(0.3634,0.2386,0.3980)$ . We obtain another minimal PH representation $(\widehat{\alpha},\widehat{A}_{\infty})$ of ${\mathcal{L}}(s)$ .

The following example shows that when the AM algorithm produces a non-zero value, quite often the minimal algebraic degree PH representation does not exist.

Example 2. Let

\begin{align*}{\mathcal{G}} (z)=\frac{0.8854z^3-0.5948z^2+0.0294z}{(1-0.5z)(1-0.2z)^2}.\end{align*}

We want to identify whether it has a third-order discrete-time PH representation. We can convert it into a continuous-time case. From Theorem 3, let

\begin{align*}{\mathcal{L}}_0(s)={\mathcal{G}}\biggl(\frac{1}{s+1}\biggr)=\frac{0.0294s^2-0.536s+0.32}{(s+0.5)(s+0.8)^2}.\end{align*}

Then we only need to verify whether ${\mathcal{L}}_0(s)$ has a third-order PH representation. Here ${\mathcal{J}}$ is

\begin{align*}\left[\begin{array}{c@{\quad}c@{\quad}c}-0.5&0 & 0\\0&-0.8&0\\0&1&-0.8\\\end{array}\right].\end{align*}

From (10), we can calculate that

\begin{align*}\beta=(13.23,-9.0316,-3.1984).\end{align*}

Now we use the AM algorithm to solve the optimization problem ( $\mathrm{OP}$ ) given in (25). Regardless of what we choose as the initial $A_0$ , the AM algorithm never produces an optimal zero-value result.

Actually, for this specific numerical example, we can also theoretically prove that ${\mathcal{G}}(z)$ has no third-order discrete-time PH representation.

The next example shows that sometimes the initial $A_0$ is critical to the efficiency of the AM algorithm (steps 1–3).

Example 3. Consider

\begin{align*}{\mathcal{L}}(s)=\frac{0.9421}{s+1}-\frac{0.0942}{(s+1.2)^2}+\frac{0.1884}{s+1.3}-\frac{0.0471}{(s+1.3)^3}.\end{align*}

If it has a sixth-order continuous-time PH representation $(\alpha, A)$ , the spectrum of A must be $\{-1,-1.2,-1.3\}$ , and their multiplicities are 1,2,3, respectively. Then the Jordan form ${\mathcal{J}}$ of A is

\begin{align*}\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}-1& & & & &\\0& -1.2& & & &\\& 1&-1.2& & &\\& &0 & -1.3& &\\& & & 1& -1.3 &\\& & & &1&-1.3\end{array}\right].\end{align*}

From (10) and (20) and Theorem 1, we can calculate that

\begin{align*}\beta=(0.9421,0.1798,-0.0904,-0.0391,0.0080,-0.0004),\end{align*}

and $\xi=1+2\times 1.2+3\times 1.3=7.3$ . Now we use the AM algorithm to solve the optimization problem ( $\mathrm{OP}$ ) given in (25). Let $A_0={\mathcal{J}}$ . The outcome of the AM algorithm is

\begin{gather*}P_{\infty} \approx \left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}0.5918 & 0.2821 & 0.0559 & 0.0516 & 0.0148 & 0.0038\\-0.0018 & 0.5788 & 0.0001 & 0.3819 & 0.0376 & 0.0035\\-0.0097 & 4.0031 & 0.5785 & -3.6285 & 0.0188 & 0.0378\\-0.0020 & -0.0037 & 0.0000 & 1.0066 & -0.0008 & 0.0000\\-0.0079 & -0.0199 & -0.0021 & 0.0248 & 1.0065 & -0.0013\\-0.0279 & 0.0157& -0.0221 & 0.0036 & 0.0239 & 1.0069\end{array}\right],\\[5pt]A_{\infty} \approx\left[\begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}-1.0004 & 0.0005 & 0.0188 & 0.0006 & 0.0008 & 0.0008\\0.0002 & -1.2005 & 0.0000 & 0.0004 & 0.0004 & 0.0006\\0.0024 & 1.0059 & -1.1995 & 0.0026 & 0.0004 & 0.0085\\0.0006 & 0.0004 & 0.0000 &-1.2992 & 0.0000 & 0.0000\\0.0003 & 0.0003 & 0.0003 & 1.0001 &-1.2995 & 0.0000\\0.0005 & 0.0008 & 0.0006 & 0.0009 & 0.9997 &-1.3011\\\end{array}\right],\end{gather*}

and ${\mathcal{F}}(A_{\infty},P_{\infty}) \approx 6.2781 \times 10^{-10}$ . By Theorem 1, let

$$\alpha=\beta P_{\infty} \approx(0.5581,0.0079,0.0004,0.4061,0.0271,0.0004).$$

We obtain a minimal PH representation $(\alpha,A_{\infty})$ of ${\mathcal{L}}(s)$ .

However, if we choose $A_0=-\xi I$ , the outcome is that ${\mathcal{F}}_{\infty} \approx0.2114$ . In this case, we cannot find a sixth-order PH representation.

Now we demonstrate the effectiveness and efficiency of the AM algorithm through a series of numerical experiments using the following procedure.

1. (P1) Generate an nth-order continuous-time PH distribution $(\alpha, A)$ randomly and compute its LST ${\mathcal{L}}(s)$ .

2. (P2) Select $(\alpha, A)$ from (P1) such that its algebraic degree is also n.

3. (P3) Calculate ${\mathcal{J}}$ from (9) in Proposition 1, $\beta$ from (10) in Proposition 3, and $\xi$ by $\xi=-\sum_{i=1}^n{\mathcal{J}}_{ii}$ .

4. (P4) Implement the AM algorithm to solve the optimization problem (OP) in (25) for ${\mathcal{J}}$ , $\beta$ , and $\xi$ .

In (P1) the choice of $\alpha=(\alpha_1,\alpha_2,\ldots,\alpha_n)$ follows the rule: $\alpha_1$ is sampled from the uniform distribution ${\mathcal{U}}(0,1)$ ; $\alpha_i$ is sampled from ${\mathcal{U}}(0,1-\sum_{j=1}^{i-1}\alpha_j)$ for $2\leq i\leq n-1$ ; and $\alpha_n=1-\sum_{j=1}^{n-1}\alpha_j$ . For the choice of $A=(a_{ij})_{n\times n}$ , the off-diagonal elements of A are sampled from independent uniform distributions ${\mathcal{U}}(0,c)$ for some positive number c; the diagonal element $a_{ii}=-\sum_{j=1,j\neq i}^n a_{ij}-\theta$ , and $\theta$ is also sampled from ${\mathcal{U}}(0,c)$ . Hence, we randomly generate a continuous-time PH distribution $(\alpha, A)$ with the algebraic degree n. For each n with $3\leq n\leq 7$ , we first repeat the procedure to randomly generate 1000 continuous-time PH distributions $(\alpha, A)$ with the algebraic degree n.

We then employ the AM algorithm by setting $A_0={\mathcal{J}} + {\mathbf{1}}^\top {\mathbf{1}}- I$ . The iterations will be terminated at step $k$ if $|{\mathcal{F}}(A_{k+1},P_{k})-{\mathcal{F}}(A_{k},P_{k})|\leq 10^{-13}$ . If the outcome objective value ${\mathcal{F}}(A_{k},P_{k})$ is less than $n^2\times10^{-10}$ , then we consider (21) to be solvable and the AM algorithm provides an nth-order PH representation. Here we solve the convex optimization problems (OP[P]) and $({\mathrm{OP}}[A])$ in each step using the ‘quadprog’ function in MATLAB. Table 1 summarizes the performance over randomly generated 1000 PH distributions for each n. These distributions are divided into two groups based on whether $\mathcal{L}(s)$ has complex poles: the ‘real case’ for distributions without complex poles and the ‘complex case’ for those with complex poles. The table shows the number of successful experiments and the average number of iterations for both cases, denoted by ‘Real Succ.’, ‘Real Ave. Iter.’, ‘Complex Succ.’, and ‘Complex Ave. Iter.’, respectively, for brevity. In addition, the table provides the total number of successful experiments, the average number of iterations, and the average runtime, denoted by ‘Total Succ.’, ‘Ave. Iter.’, and ‘Ave. Runtime’.

Table 1.

Outcome of the AM algorithm

A few of the following observations are made from Table 1. (a) In each experiment, simply keeping $A_0$ fixed at $A_0={\mathcal{J}} + {\mathbf{1}}^\top {\mathbf{1}}- I$ , the AM algorithm can effectively obtain a minimal PH representation in most cases. (b) It can be seen that all failures occur in the complex case, and the number of iterations in the complex case is significantly higher than in the real case. Therefore, the AM algorithm is more efficient for the real case compared with the complex case. (c) As n increases, the average number of iterations and the average runtime tend to increase.

Finally, we compare the AM algorithm with the MM method proposed in [Reference Mészáros, Horváth and Telek23, Reference Telek and Horváth24]. For a given ${\mathcal{L}}(s)$ defined in (1), their MM method starts with a matrix representation (v, H) of ${\mathcal{L}}(s)$ , where $v\in \mathbb{R}^n$ and $H\in \mathbb{R}^{n\times n}$ . We say (v, H) is a matrix representation of ${\mathcal{L}}(s)$ if $-vH(sI-H)^{-1}{\mathbf{1}}^\top={\mathcal{L}}(s)$ . Note that the components of v and the off-diagonal elements of H may not be non-negative, and the diagonal elements of H may not be negative. To address this incorrectness issue for H to be a generator, they introduce different parameterized error functions to penalize the terms with incorrect positive or negative signs. In each step, they choose an elementary transformation that minimizes the error function. By applying such elementary transformations to the initial matrix representation in each step, their MM method may or may not generate a PH representation. To make these two comparable, the initial matrix representation for the MM method is also set to $(\beta, {\mathcal{J}})$ defined by (10). From (25), this initial setting makes both algorithms utilize the same information of ${\mathcal{L}}(s)$ , namely $\beta$ and ${\mathcal{J}}$ .

To evaluate their performance, we follow the procedure (P1)–(P4). In particular, to generate the following three cases studied in [Reference Mészáros, Horváth and Telek23, Reference Telek and Horváth24], the choices for the off-diagonal elements of $A=(a_{ij})_{n\times n}$ in (P1) are modified by the following.

• • Balanced case (BC): its off-diagonal elements are sampled from independent uniform distributions ${\mathcal{U}}(0,c)$ .

• • Sparse case (SP): with probability p, its off-diagonal elements are set to zero, and with probability $1-p$ , they are sampled from independent uniform distributions ${\mathcal{U}}(0,c)$ .

• • Stiff case (ST): with probability p, its off-diagonal elements are sampled from independent uniform distributions ${\mathcal{U}}(0,1000c)$ , and with probability $1-p$ , they are sampled from independent uniform distributions ${\mathcal{U}}(0,c)$ .

For the sparse case, we consider $p=0.5, 0.7$ , and $0.9$ and, for the stiff case, we consider $p=0.3, 0.5$ , and $0.7$ . In summary, three cases with seven scenarios are considered. For each scenario, 100 continuous-time PH distributions $(\alpha, A)$ with the algebraic degree n are randomly generated for $3\leq n\leq 6$ . For the AM algorithm starting with $A_0={\mathcal{J}} +{\mathbf{1}}^\top {\mathbf{1}}- I$ , the iterations are terminated at step $k$ if $|{\mathcal{F}}(A_{k+1},P_{k})-{\mathcal{F}}(A_{k},P_{k})|\leq 10^{-13}$ . If the outcome objective value ${\mathcal{F}}(A_{k},P_{k})$ is less than $n^2\times10^{-10}$ , then we consider (21) to be solvable and the AM algorithm provides an nth-order PH representation. For the MM method starting with $A_0={\mathcal{J}} + {\mathbf{1}}^\top {\mathbf{1}}- I$ , its several modified versions are developed in [Reference Mészáros, Horváth and Telek23, Reference Telek and Horváth24]. Here we report the best results over all its modified versions. When executing the MM method, the number of iterations should be determined to terminate the computation, and the error functions may be adjusted to improve the performance. To ensure that the two algorithms are compared within the same error range, the tolerance is set to $10^{-5}$ . That is, if the error in the incorrect positive or negative signs of the outcome is within the tolerance, we consider that the MM method has successfully provided a PH representation.

Table 2 summarizes the number of successful experiments separately obtained by the AM algorithm and MM method for the randomly generated 100 PH distributions for each n and each scenario. ‘Total’ represents the number of successes, where either the AM algorithm or the MM method can successfully provide an nth-order PH representation. The table indicates that the two algorithms exhibit almost similar performance. For the BC, both algorithms perform well. When the generator becomes sparser or contains elements with a big difference, however, neither of them performs well.

Table 2.

Outcome of the AM algorithm and MM method

Actually, the two algorithms operate within two fundamentally different frameworks. In the AM algorithm, at each step $k$ , the constraints ensure that $\beta P_k$ remains a distribution and $A_k$ is a generator while minimizing the objective function aims to enhance the similarity between $A_k$ and $\mathcal{J}$ , thereby aligning the Laplace transform of $(\beta P_k, A_k)$ as closely as possible with $\mathcal{L}(s)$ . In contrast, for the MM method, the elementary transformation at step $k$ with output $(\alpha_k, A_k)$ keeps the similarity between $A_k$ and $\mathcal{J}$ , ensuring that the Laplace transform of $(\alpha_k, A_k)$ is $\mathcal{L}(s)$ , while the error functions attempt to ensure that $\alpha_k$ is an initial distribution and $A_k$ is a generator. When both of them yield a solution, the AM algorithm guarantees that the outcome $(\beta P_{\infty}, A_{\infty})$ is a PH representation with any error reflected in the difference between the Laplace transform of $(\beta P_{\infty}, A_{\infty})$ and $\mathcal{L}(s)$ . In contrast, the outcome $(\alpha_{\infty}, A_{\infty})$ produced by the MM method implies that its Laplace transform is $\mathcal{L}(s)$ ; however, due to the presence of a tolerance, $(\alpha_{\infty}, A_{\infty})$ may not be a PH representation. Moreover, from the robustness and simplicity perspective, the AM algorithm proves to be more robust and simpler across different problems, as there is no need to properly and dynamically tune the error functions and set the step length to determine the similarity transformation matrix, which are critical to the MM method. The following example studied in [Reference Telek and Horváth24] illustrates this in a more concrete way.

Example 4. Consider

\begin{align*}\beta=(0.2,0.3,0.5) \quad \mbox{and} \quad {\mathcal{J}}=\begin{bmatrix}-1&0&0\\0&-3&h\\0&-h&-3\end{bmatrix}.\end{align*}

Let ${\mathcal{L}}(s)=-\beta {\mathcal{J}}(sI-{\mathcal{J}})^{-1}{\mathbf{1}}^\top$ . In [Reference Telek and Horváth24] it was pointed out that if $h\leq 0.552\,748\,375$ , ${\mathcal{L}}(s)$ has a third-order continuous-time PH representation. As h tends to $0.552\,748\,375$ , when using the MM algorithm, modifications to the parameters of the error function are necessary to obtain a third-order PH representation. However, with the AM algorithm, we can always choose $A_0={\mathcal{J}} + {\mathbf{1}}^\top{\mathbf{1}}- I$ to obtain a third-order PH representation. For instance, if $h=0.545$ , the outcome of the AM algorithm is

\begin{gather*}P_{\infty}\approx\left[\begin{array}{c@{\quad}c@{\quad}c}0.395\, 172\,218\,1 & 0.353\,666\,305\,1 & 0.251\,161\,476\,7\\-0.112\,701\,865\,8 & 0.868\,412\,056\,8 & 0.244\,289\,808\,9\\-0.090\,447\,764\,6 & 1.254\,052\,302\,1 & -0.163\,604\,537\,4\end{array}\right],\\[3pt]A_{\infty}\approx\left[\begin{array}{c@{\quad}c@{\quad}c}-1.162\,767\,578\,5 & 0.000\,174\,459\,2 & 1.162\,324\,760\,7\\0.181\,546\,631\,1 & -3.008\,462\,097\,6 & 0.000\,011\,359\,2\\0.000\,448\,511\,8 & 2.827\,890\,616\,1 & -2.828\,786\,558\,6\end{array}\right],\end{gather*}

and ${\mathcal{F}}(A_{\infty},P_{\infty})\approx1.5411\times 10^{-10}$ . By Theorem 1, $(\alpha,A_\infty)$ is a third-order PH distribution with $\alpha=\beta P_{\infty}\approx(0.000\,000\,001\,6,0.958\,283\,029\,1,0.041\,716\,969\,3)$ . If $h=0.552$ , the outcome of the AM algorithm is

\begin{gather*}P_{\infty}\approx\left[\begin{array}{c@{\quad}c@{\quad}c}0.395\,535\,909\,5 & 0.243\,160\,386\,5 & 0.361\,303\,703\,9\\-0.112\,989\,907\,0 & 0.228\,432\,582\,3 & 0.884\,557\,324\,7\\-0.090\,420\,403\,0 & -0.193\,095\,478\,9 & 1.283\,515\,881\,9\end{array}\right],\\[3pt]A_{\infty}\approx\left[\begin{array}{c@{\quad}c@{\quad}c}-1.163\,003\,988\,2& 1.162\,797\,680\,2 & 0.000\,066\,524\,2\\0.000\,537\,465\,5 & -2.891\,464\,676\,7 & 2.890\,692\,306\,2\\0.178\,085\,932\,4 & 0.000\,003\,703\,6 & -2.945\,532\,324\,8\end{array}\right],\end{gather*}

and ${\mathcal{F}}(A_{\infty},P_{\infty})\approx2.0198\times 10^{-12}$ . By Theorem 1, $(\alpha,A_\infty)$ is a third-order PH distribution with $\alpha=\beta P_{\infty}\approx(0.000\,000\,008\,3,0.020\,614\,112\,6,0.979\,385\,879\,1)$ .

6. Concluding remarks

In this paper, we have investigated the minimal representation problem of PH distributions starting with their LST. The if-and-only-if condition has been established to verify whether a PH distribution has the same order representation as its algebraic degree. To this end, we have developed a unique expression for the LST in terms of the Jordan form. The unique expression helps us equivalently transform the minimal representation problem into the solution existence problem for the set of linear and quadratic equations. The latter problem is further equivalently transformed into the non-convex optimization problem, and the AM algorithm is then developed to solve this optimization problem, the proposed algorithm is proved to be convergent for any initial value, and numerically demonstrated to be quite effective. Moreover, the approach we have developed for the continuous-time PH distributions can be applied to the discrete-time PH distributions by establishing the equivalence between the LST of continuous-time PH distributions and the generating function of discrete-time distributions.

The next research topic will be to consider the case in which the order of its PH representation is not the same as its algebraic degree. The challenge for this case is that from the LST’s poles, how can one identify the other poles. This is very subtle from the discussions conducted in [Reference O’Cinneide30, Reference Dehon and Latouche36].