3.1 Matrix quadratics, pencils, and doubling

We begin by expressing the UQME in (5) as a subspace problem via the first companion linearization of the matrix quadratic problem (Hammarling et al. Reference Hammarling, Munro and Tisseur2013)

(6) \begin{align} \mathscr{A} \mathscr{X} = \mathscr{B} \mathscr{X} \mathscr{M} \end{align}

with

\begin{align*} \mathscr{X}= \begin{pmatrix} I \\ P \end{pmatrix} && \mathscr{A}= \begin{pmatrix} 0 &\quad I \\ C &\quad B \end{pmatrix}, && \mathscr{B}= \begin{pmatrix} I &\quad 0\\ 0 &\quad -A \end{pmatrix}, && \mathscr{M}=P \end{align*}

Clearly, any $P$ satisfying (5) is a solution of (6). Further, the eigenvalues of $\mathscr{M}$ are a subset of the generalized eigenvalues of the matrix pencil $\mathscr{A}-\lambda \mathscr{B}$ , i.e., $\operatorname {eig}(\mathscr{M})\subset \operatorname {eig}(\mathscr{A},\mathscr{B})$ .

We make the standard assumption in the literature for the existence of a unique stable solvent $P$ , i.e., Blanchard and Kahn (Reference Blanchard and Kahn1980) celebrated rank and order conditions.

The quadratic convergence of our SDAs requires the existence of a unique solvent to the reverse problem, which follows directly from Assumption 1.Footnote ⁹

Corollary 1 (Solvent of the Dual Matrix Equation). Let Assumption 1 hold. There exists a $P_d\in \mathbb{R}^{n_y \times n_y}$ with $\rho (P_d)\lt 1$ that satisfies the dual matrix quadratic

\begin{equation*}0=C P_d^2+B P_d+A.\end{equation*}

The problem in (6) is an eigenvalue problem and can thus be solved numerically using the QZ or generalized Schur decomposition of Moler and Stewart (Reference Moler and Stewart1973). In the online appendix, we derive the solution with QZ by working directly with the linear algebraic problem. This should link the more familiar QZ with the doubling algorithms we will subsequently present more clearly to readers looking for additional intuition.

We will approach the matrix pencil problem here via a doubling approach. In detail we will try to generate sequences $\{\mathscr{A}_k\}_{k=1}^N$ and $\{\mathscr{B}_k\}_{k=1}^N$ such that

(7) \begin{align} \mathscr{A}_k \mathscr{X} = \mathscr{B}_k \mathscr{X} \mathscr{M\,}^{2^k} \end{align}

squaring the matrix $\mathscr{M}$ in each iteration. At the same time, we will seek to preserve the structure of the pencil, e.g., a special block structure of $\mathscr{A}_k$ and $\mathscr{B}_k$ . This will enable us to express the algorithm in terms of submatrices of $\mathscr{A}_k$ and $\mathscr{B}_k$ instead of the matrices in their entirety and makes sure that we can repeat the doubling step in the next iteration.

Following Guo et al. (Reference Guo, Lin and Xu2006), a transformation $\,\,\,\,\widehat{\!\!\!\!\mathscr{A}}-\lambda \,\,\widehat {\!\!\mathscr{B}}$ of a pencil $\mathscr{A}-\lambda \mathscr{B}$ is called a doubling transform if

(8) \begin{align} \,\,\,\,\widehat{\!\!\!\!\mathscr{A}}=\overline {\mathscr{A}}\mathscr{A},\;\;\;\;\,\,\widehat {\!\!\mathscr{B}}=\overline {\mathscr{B}}\mathscr{B} \end{align}

for $\overline {\mathscr{A}}$ and $\overline {\mathscr{B}}$ that satisfy

(9) \begin{align} \text{rank}\left (\begin{bmatrix} \overline {\mathscr{A}} &\quad \overline {\mathscr{B}}\end{bmatrix}\right )=2n_y,\;\;\;\;\begin{bmatrix} \overline {\mathscr{A}}&\quad \overline {\mathscr{B}}\end{bmatrix}\begin{bmatrix} \mathscr{B}\\-\mathscr{A}\end{bmatrix}=0 \end{align}

That is, we will be transforming the problem here as in QZ, seeking a structure that is amenable to doubling instead of the upper triangularity sought there.

Such a doubling transformation is eigenspace preserving and eigenvalue squaring (“doubling”) following Guo et al. (Reference Guo, Lin and Xu2006, Theorem 2.1) or Huang et al. (Reference Huang, Li and Lin2018, Theorem 3.1) that we repeat here adapted to our problem

Theorem 1 (Doubling Pencil). Suppose $\,\,\,\,\widehat{\!\!\!\!\mathscr{A}}-\lambda \,\,\widehat {\!\!\mathscr{B}}$ is a doubling transformation of the pencil $\mathscr{A}-\lambda \mathscr{B}$ . Then, as $\mathscr{A} \mathscr{X} = \mathscr{B} \mathscr{X} \mathscr{M}$ it follows from ( 6 ) that

(10) \begin{align} \widehat{\!\!\!\!\mathscr{A}} \mathscr{X} = \,\,\widehat {\!\!\mathscr{B}} \mathscr{X} \mathscr{M\,}^2 \end{align}

Following this theorem, we obtain a doubling algorithm for (6) by iterating on

(11) \begin{align} \underbrace {\,\,\,\,\widehat{\!\!\!\!\mathscr{A}}_{k}}_{\mathscr{A}_{k+1}}=\overline {\mathscr{A}}_{k}\mathscr{A}_{k},\;\;\;\;\underbrace {\,\,\widehat {\!\!\mathscr{B}}_{k}}_{\mathscr{B}_{k+1}}=\overline {\mathscr{B}}_{k}\mathscr{B}_{k} \end{align}

initializing with $\mathscr{A}$ and $\mathscr{B}$ , as $\mathscr{A}\mathscr{X}=\mathscr{B}\mathscr{X}\mathscr{M}$ then gives

(12) \begin{align} \mathscr{A}_{1}\mathscr{X}&=\mathscr{B}_{1}\mathscr{X}\mathscr{M\,}^2&&\rightarrow & \mathscr{A}_{2}\mathscr{X}&=\mathscr{B}_{2}\mathscr{X}\mathscr{M\,}^4&&\rightarrow & \mathscr{A}_{k}\mathscr{X}&=\mathscr{B}_{k}\mathscr{X}\mathscr{M\,}^{2^k} \end{align}

We seek a structure in $\mathscr{A}$ and $\mathscr{B}$ such that we calculate $\mathscr{A}_{k}\rightarrow \mathscr{A}_{k+1}$ and $\mathscr{B}_{k}\rightarrow \mathscr{B}_{k+1}$ by recursions in elements (here we will settle for submatrices). In the following subsections, we provide two recursions of exactly this structure. Both require that we rearrange our pencil $\mathscr{A}-\lambda \mathscr{B}$ to conform to a specific structure of each recursion, as we now show.

3.2 First Standard Form

To obtain the first SDA, we will transform the pencil $\mathscr{A}-\lambda \mathscr{B}$ . Assuming that $B$ is non-singular, we receive the primal problem in SF1

(13) \begin{align} \mathscr{A}_0 \mathscr{X} = \mathscr{B}_0 \mathscr{X} \mathscr{M} \end{align}

by multiplying (6) from the left by $S$ , where

\begin{align*} \mathscr{A}_0 &= \begin{pmatrix} E_0 &\quad 0 \\ -X_0 &\quad I \end{pmatrix}\,:\!=\,S\mathscr{A}, &\quad \mathscr{B}_0 &\quad = \begin{pmatrix} I &\quad -Y_0 \\ 0 &\quad F_0 \end{pmatrix}\,:\!=\,S\mathscr{B}, &\quad S &\quad = \begin{pmatrix} I &\quad -B^{-1}\\ 0 &\quad B^{-1} \end{pmatrix} \end{align*}

Note that (6) and (13) are equivalent in the sense that the pencils $\mathscr{A}-\lambda \mathscr{B}$ and $\mathscr{A}_0-\lambda \mathscr{B}_0$ share the same set of generalized eigenvalues, i.e., $\operatorname {eig}(\mathscr{A},\mathscr{B})=\operatorname {eig}(\mathscr{A}_0,\mathscr{B}_0)$ . We now use doubling transformations to arrive at a SDA for SF1, computing $\{{\mathscr{A}}_{k}\}_{k=0}^{\infty }$ , $\{{\mathscr{B}}_{k}\}_{k=0}^{\infty }$ with

(14) \begin{align} \mathscr{A}_{k+1}&= \begin{pmatrix} E_{k+1} &\quad 0 \\ -X_{k+1} &\quad I \end{pmatrix} \,:\!=\, \begin{pmatrix} E_k\left (I-Y_k X_k \right )^{-1} &\quad 0\\ -F_k\left (I-X_k Y_k\right )^{-1} X_k &\quad I \end{pmatrix} \mathscr{A}_{k} \\[-12pt] \nonumber \end{align}

(15) \begin{align} \mathscr{B}_{k+1}&= \begin{pmatrix} I &\quad -Y_{k+1} \\ 0 &\quad F_{k+1} \end{pmatrix} \,:\!=\, \begin{pmatrix} I &\quad -E_k\left (I-Y_k X_k \right )^{-1} Y_k \\ 0 &\quad F_k\left (I-X_k Y_k\right )^{-1} \end{pmatrix} \mathscr{B}_{k} \\[6pt] \nonumber \end{align}

such that

(16) \begin{align} \mathscr{A}_k \mathscr{X} = \mathscr{B}_k \mathscr{X} \mathscr{M\,}^{2^k} \end{align}

A key feature here is that (16) retains SF1 for all $k\in \mathbb{N}$ . Algorithm 1 with an initial guess $P_0=0$ summarizes this SDA for the SF1. The intuition here is that under some fairly general preconditions, e.g., Assumption 1, the term $\mathscr{B}_k \mathscr{X} \mathscr{M\,}^{2^k}$ on the right-hand side of (16) converges to zero for $k\rightarrow \infty$ , so that consequently $X_k$ converges to $P$ .

A disadvantage of SDAs is that numerical inaccuracies can propagate through iteration if, say, the matrices to be inverted are not well-conditioned. In the context of discrete algebraic Riccati equations (DAREs), Mehrmann and Tan (Reference Mehrmann and Tan1988) show that such a defection of the approximate solution again satisfies a DARE. As a result, after solving a DARE, one can solve the associated DAREs of the approximation error to increase the overall accuracy. Following this, Bini and Meini (Reference Bini and Meini2023) show how to incorporate an initial guess to SDA for the SF1 by means of an equivalence transformation of the pencil $\mathscr{A}-\lambda \mathscr{B}$ . Huang et al. (Reference Huang, Li and Lin2018) discuss similar transformations for algebraic Riccati equations (AREs) in general. We introduce such an initial guess by transforming the eigenvalue problem (6) to

(17) \begin{align} \,\,\,\,\widehat{\!\!\!\!\mathscr{A}}\, \,\,\,\widehat {\!\!\!\mathscr{X}}\, = \,\,\widehat {\!\!\mathscr{B}}\, \,\,\,\widehat {\!\!\!\mathscr{X}}\, \mathscr{M} \end{align}

with $P_0$ as the initial guess such that $P=\widehat {P}+P_0$ and

\begin{align*} \,\,\,\widehat {\!\!\!\mathscr{X}}\, \,:\!=\, \begin{pmatrix} I \\ \hat {P} \end{pmatrix} && \,\,\,\,\widehat{\!\!\!\!\mathscr{A}}\, \,:\!=\, \mathscr{A} \begin{pmatrix} I &\quad 0 \\ P_0 &\quad I \end{pmatrix}, && \,\,\widehat {\!\!\mathscr{B}}\, \,:\!=\, \mathscr{B} \begin{pmatrix} I &\quad 0 \\ P_0 &\quad I \end{pmatrix} \end{align*}

Assuming that $B+AP_0$ has full rank, we multiply $\,\,\,\,\widehat{\!\!\!\!\mathscr{A}}$ and $\,\,\widehat {\!\!\mathscr{B}}$ by

\begin{align*} \widehat {S}= \begin{pmatrix} \left (B+AP_0\right )^{-1}B &\quad -\left (B+AP_0\right )^{-1}\\ I-\left (B+AP_0\right )^{-1}B &\quad\left (B+AP_0\right )^{-1} \end{pmatrix} \end{align*}

to receive the corresponding problem in SF1, i.e.,

(18) \begin{align} \,\,\,\,\widehat{\!\!\!\!\mathscr{A}}_0\, \,\,\,\widehat {\!\!\!\mathscr{X}}\, = \,\,\widehat {\!\!\mathscr{B}}_0\, \,\,\,\widehat {\!\!\!\mathscr{X}}\, \mathscr{M} \end{align}

with

\begin{align*} \,\,\,\,\widehat{\!\!\!\!\mathscr{A}}_0\, &= \begin{pmatrix} \widehat {E}_0 &\quad 0 \\ -\widehat {X}_0 &\quad I \end{pmatrix} \,:\!=\, \widehat {S}\,\,\,\,\widehat{\!\!\!\!\mathscr{A}}, &\quad\hat {\mathscr{B}}_0\, &= \begin{pmatrix} I &\quad -\widehat {Y}_0 \\ 0 &\quad \widehat {F}_0 \end{pmatrix} \,:\!=\, \widehat {S}\,\,\widehat {\!\!\mathscr{B}} \end{align*}

Algorithm 1 summarizes the SDA for SF1 based on an initial $P_0$ . Note that a good $P_0$ should increase the probability that $B+AP_0$ is well-conditioned, since under Assumption 1 we know that $B+AP$ has full rank (see Lan and Meyer-Gohde, Reference Lan and Meyer-Gohde2014). In comparison to $P_0=0$ above which is only applicable for non-singular $B$ , Algorithm 1 can handle a singular $B$ , provided we have a $P_0$ such that $B+AP_0$ is non-singular.Footnote ¹⁰

Algorithm 1:

Structure-Preserving Doubling Algorithm (SF1)

3.3 Second Standard Form

Alternatively, Chiang et al. (Reference Chiang, Chu, Guo, Huang, Lin and Xu2009) use a doubling algorithm on the eigenvalue problem

(19) \begin{align} {\mathscr{A\,\,}}^{\dagger }_0 {\mathscr{X\,\,}}^{\dagger } = {\mathscr{B\,}}^{\dagger }_0 {\mathscr{X\,\,}}^{\dagger } \mathscr{M} \end{align}

with

\begin{align*} {\mathscr{X\,\,}}^{\dagger }&= \begin{pmatrix} I \\ AP \end{pmatrix}, & {\mathscr{A\,\,}}^{\dagger }_0 &= \begin{pmatrix} {E}^{\dagger }_0 &\quad 0 \\ -{X}^{\dagger }_0 &\quad I \end{pmatrix}\,:\!=\, \begin{pmatrix} -C &\quad 0 \\ 0 &\quad I \end{pmatrix}, &{\mathscr{B\,}}^{\dagger }_0&= \begin{pmatrix} -{Y}^{\dagger }_0 &\quad I \\ -{F}^{\dagger }_0 &\quad 0 \end{pmatrix}\,:\!=\, \begin{pmatrix} B &\quad I \\ A &\quad 0 \end{pmatrix}, \mathscr{M}=P \end{align*}

which satisfies the so-called Second Standard Form (SF2). Again any $P$ satisfying (5) is also a solution of (19). Note that $\mathscr{A}-\lambda \mathscr{B}$ and ${\mathscr{A\,\,}}^{\dagger }_0-\lambda {\mathscr{B\,}}^{\dagger }_0$ again share the same set of generalized eigenvalues, i.e., $\operatorname {eig}(\mathscr{A},\mathscr{B})=\operatorname {eig}({\mathscr{A\,\,}}^{\dagger }_0,{\mathscr{B\,}}^{\dagger }_0)$ .

Similarly to SF1 above, the SDA for SF2 computes sequences $\{{\mathscr{A\,\,}}^{\dagger }_{k}\}_{k=0}^{\infty }$ , $\{{\mathscr{B\,}}^{\dagger }_{k}\}_{k=0}^{\infty }$ with

(20) \begin{align} {\mathscr{A\,\,}}^{\dagger }_{k+1}&= \begin{pmatrix} {E}^{\dagger }_{k+1} &\quad 0 \\ -{X}^{\dagger }_{k+1} &\quad I \end{pmatrix} \,:\!=\, \begin{pmatrix} {E}^{\dagger }_k\left ( {X}^{\dagger }_k - {Y}^{\dagger }_k \right )^{-1} &\quad 0 \\ {F}^{\dagger }_k\left ( {X}^{\dagger }_k - {Y}^{\dagger }_k \right )^{-1} &\quad I \end{pmatrix} \mathscr{A}_{k} \end{align}

(21) \begin{align} {\mathscr{B\,}}^{\dagger }_{k+1}&= \begin{pmatrix} -{Y}^{\dagger }_{k+1} &I \\ -{F}^{\dagger }_{k+1} &0 \end{pmatrix} \,:\!=\, \begin{pmatrix} I & {E}^{\dagger }_k\left ( {X}^{\dagger }_k - {Y}^{\dagger }_k \right )^{-1} \\ 0 & {F}^{\dagger }_k\left ( {X}^{\dagger }_k - {Y}^{\dagger }_k \right )^{-1} \end{pmatrix} \mathscr{B}_{k} \\[6pt] \nonumber \end{align}

such that

(22) \begin{align} {\mathscr{A\,\,}}^{\dagger }_k {\mathscr{X\,\,}}^{\dagger } = {\mathscr{B\,}}^{\dagger }_k {\mathscr{X\,\,}}^{\dagger } \mathscr{M\,}^{2^k} \end{align}

Note that compared to the SDA for SF1, the sequence $\{{X}^{\dagger }_{k}\}_{k=0}^{\infty }$ now converges to $AP$ instead of $P$ . However, in case of convergence we may obtain a $P$ as $-({X}^{\dagger }_k+B)^{-1}C$ .

Following the idea of Bini and Meini (Reference Bini and Meini2023) and applying it to the SDA for SF2, we take again $P_0$ as the initial guess such that $P=\widehat {P}+P_0$ and insert this into the UQME (5)

(23) \begin{align} 0&=A \left (\widehat {P}+P_0\right )^2+B \left (\widehat {P}+P_0\right )+C=A\widehat {P}^2+A \widehat {P}P_0 +\left (AP_0+B\right )\widehat {P}+AP_0^2+B P_0+C \end{align}

This can be rewritten as (19) with

\begin{align*} \,\,\,\widehat {\!\!\!\mathscr{X\,\,}}^{\dagger }&= \begin{pmatrix} I \\ A\widehat {P} \end{pmatrix}, & \,\,\,\,\widehat{\!\!\!\!\mathscr{A\,\,\,}}^{\dagger }_0&= \begin{pmatrix} \widehat {E\,}^{\dagger }_0 &\ 0 \\ -\widehat {X}^{\dagger }_0 &\ I \end{pmatrix}\,:\!=\, \begin{pmatrix} -C &\ 0 \\ A P_0 &\ I \end{pmatrix}, &\,\,\widehat {\!\!\mathscr{B\,}}^{\dagger }_0 &= \begin{pmatrix} -\widehat {Y}^{\dagger }_0 &\ I \\ -\widehat {F}^{\dagger }_0 &\ 0 \end{pmatrix}\,:\!=\, \begin{pmatrix} AP_0+B &\ I \\ A &\ 0 \end{pmatrix}, \mathscr{M}=P \end{align*}

Algorithm 2 summarizes the SDA for SF2 based on an initial guess $P_0$ .

While this is a straightforward application of the initial guess approach of Algorithm 1 to the SDA for SF2, the algorithm will deliver the same approximation to $P$ regardless of $P_0$ . We summarize this in the following (Huang et al. Reference Huang, Li and Lin2018, Theorem 3.32).

Theorem 2. Let ${X}^{\dagger }_{t}$ , ${Y}^{\dagger }_{t}$ , ${E}^{\dagger }_{t}$ , and ${F}^{\dagger }_{t}$ denote the quantities of the Structure-Preserving Doubling Algorithm (SF2) and let $\widehat {X}^{\dagger }_{t}$ , $\widehat {Y}^{\dagger }_{t}$ , $\widehat {E\,}^{\dagger }_{t}$ , and $\widehat {F}^{\dagger }_{t}$ be the quantities of the Structure-Preserving Doubling Algorithm (SF2) with an initial guess $P_0$ . Then we have that

\begin{align*} \widehat {X}^{\dagger }_{t}&={X}^{\dagger }_{t}-AP_0,& \widehat {Y}^{\dagger }_{t}&={Y}^{\dagger }_{t}-AP_0\\ \widehat {E\,}^{\dagger }_{t}&={E}^{\dagger }_{t},& \widehat {F}^{\dagger }_{t}&={F}^{\dagger }_{t} \end{align*}

Hence, both algorithms will eventually return the same approximation for $P$ .

We see that it is not trivial to generate versions of these algorithms that enable refinement of arbitrary initializations of the solution. The results above are known for Riccati equations, we have extended them here to the specific case of our UQME. That we can provide a version of SF1 in Algorithm 1 which operates on arbitrary initializations is all the more interesting as it can potentially profit in terms of increased accuracy and reduced computation time from initializations that are close to the stable solvent $P$ being sought. We will confirm this and, following Theorem 2, that the Algorithm 2 does not.

Algorithm 2:

Structure-Preserving Doubling Algorithm (SF2)

4.1 Cyclic and Logarithmic Reduction

SDAs generate a sequence of matrix pencils, in each step squaring the corresponding eigenvalues. Similarly, Bini and Meini’s (Reference Bini and Meini1996) cyclic reduction and Latouche and Ramaswami’s (Reference Latouche and Ramaswami1993) logarithmic reduction, generate sequences of matrix polynomials whose eigenvalues are squared in each step. In the following we outline the idea of the cyclic SDA of SF2.Footnote ¹¹

Cyclic reduction computes the sequences $\{A_k\}_{k=0}^{\infty }$ , $\{B_k\}_{k=0}^{\infty }$ , $\{C_k\}_{k=0}^{\infty }$ , and $\{\widehat {B}_k\}_{k=0}^{\infty }$ with

(24) \begin{align} A_{k+1} &= - A_{k}\,B_{k}^{-1}\,A_{k}, &A_{0}&=A \\[-10pt] \nonumber\end{align}

(25) \begin{align} B_{k+1} &= B_{k} - A_{k}\,B_{k}^{-1}\,C_{k} - C_{k}\,B_{k}^{-1}\,A_{k}, &B_{0}&=B \\[-10pt] \nonumber \end{align}

(26) \begin{align} C_{k+1} &= - C_{k}\,B_{k}^{-1}\,C_{k}, &C_{0}&=C \\[-10pt] \nonumber \end{align}

(27) \begin{align} \widehat {B}_{k+1}&= \widehat {B}_{k} - A_{k}\,B_{k}^{-1}\,C_{k} &\widehat {B}_{0}&=B \\[9pt] \nonumber \end{align}

Using divide-and-conquer, cyclic reduction defines a sequence of UQMEs with

(28) \begin{align} 0 &= A_{k} \mathscr{M}_{k}^{\,2} + B_{k} \mathscr{M}_{k} + C_{k}, & 0 &= A_{k} \mathscr{M}_{k}\mathscr{M} + \widehat {B}_{k} \mathscr{M} + C, \end{align}

where $\mathscr{M}_{k}=\mathscr{M\,}^{2^{k}}$ . Furthermore, assuming that $\lim \limits _{k \rightarrow \infty }\widehat {B}_{k}^{-1}$ exists and that $\lim \limits _{k \rightarrow \infty }\widehat {B}_{k}^{-1}A_{k}\mathscr{M}_{k}\mathscr{M}=0$ we can express the solvent $P$ as $P =- \lim \limits _{k \rightarrow \infty } \widehat {B}_{k}^{-1}C$ .

As pointed out by Bini et al. (Reference Bini, Iannazzo and Meini2011, pp. 167) and Chiang et al. (Reference Chiang, Chu, Guo, Huang, Lin and Xu2009, pp. 236), the SDA for SF2 is connected to cyclic reduction. In detail, it follows via induction that

\begin{align*} \widehat {B}_{k} &= B+{X}^{\dagger }_k, &B_{k} &={X}^{\dagger }_k-{Y}^{\dagger }_k, &C_{k} &=-{E}^{\dagger }_k, &A_{k} &= -{F}^{\dagger }_k, & \forall k \geq 0 \end{align*}

Hence, the SDA for SF2 with $P_0=0$ and cyclic reduction are theoretically equivalent.

This equivalence with a special case ( $P_0=0$ ) of an SDA is significant as this cyclic reduction is included yet only mildly documented in Dynare. The same is true for the logarithmic reduction of Latouche and Ramaswami (Reference Latouche and Ramaswami1993), which uses the same divide-and-conquer strategy to obtain $\{H_k\}_{k=0}^{\infty }$ , $\{L_k\}_{k=0}^{\infty }$ , $\{\widehat {H}_k\}_{k=0}^{\infty }$ , and $\{\widehat {L}_k\}_{k=0}^{\infty }$ with

(29) \begin{align} H_{k+1} &= \left (I-H_{k}L_{k}-L_{k}H_{k}\right )^{-1}H_{k}^2, &H_{0}= -B^{-1}A \\[-8pt] \nonumber\end{align}

(30) \begin{align} L_{k+1} &= \left (I-H_{k}L_{k}-L_{k}H_{k}\right )^{-1}L_{k}^2, &L_{0}= -B^{-1}C, \\[-8pt] \nonumber\end{align}

(31) \begin{align} \widehat {H}_{k+1} &= \widehat {H}_{k}H_{k+1}, &\widehat {H}_{0} = -B^{-1}A, \\[-8pt] \nonumber\end{align}

(32) \begin{align} \widehat {L}_{k+1} &= \widehat {L}_{k} +\widehat {H}_{k}L_{k+1}, &\widehat {L}_{0} = -B^{-1}C, \\[8pt] \nonumber \end{align}

that define a sequence of UQMEs similarly to the cyclic reduction above with

(33) \begin{align} 0&= H_{k} \mathscr{M}_{k}^{\,2} - \mathscr{M}_{k} + L_{k}. & 0&= \widehat {H}_{k} \mathscr{M}_{k+1} - \mathscr{M} + \widehat {L}_{k} \end{align}

As pointed out by Bini et al. (Reference Bini, Latouche and Meini2005, Theorem 7.5), it follows via induction that for all $k\geq 0$

\begin{align*} H_{k} &= -B_{k}^{-1}A_{k}, &L_{k} &= -B_{k}^{-1}C_{k}. \end{align*}

Consequently, there is also a link between the SDA for SF2 with $P_0=0$ and the logarithmic reduction. In detail, we get for all $k \geq 0$ that

\begin{align*} H_{k}&=\left ({X}^{\dagger }_k-{Y}^{\dagger }_k\right )^{-1}{F}^{\dagger }_k &L_{k}&=\left ({X}^{\dagger }_k-{Y}^{\dagger }_k\right )^{-1}{E}^{\dagger }_k \end{align*}

Similar to the cyclic reduction assuming here $\lim \limits _{k \rightarrow \infty } \widehat {H}_{k} \mathscr{M}_{k+1}=0$ the stable solvent is

(34) \begin{align} P &= \lim \limits _{k \rightarrow \infty } \widehat {L}_{k}, \end{align}

Pseudocode of both reductions in our notation are in the online appendix. Compared to the doubling algorithm, the reductions follow a similar idea by squaring the eigenvalues of matrix polynomials. Beyond that, abstracting from numerical inaccuracies, we find that the SDA for SF2 with $P_0=0$ and the cyclic reduction will deliver identical approximation to $P$ for a given number of iterations. While the two reduction generate interchangeable sequences of UQMEs (the left equations in (28) and (33)), the two algorithms differ in the way in which they recover the approximation to $P$ (the right equations in (28) and (33)). Hence, although we can link some quantities computed by the logarithmic reduction to the quantities of the cyclic reduction / SF2, we will in general receive distinct solutions $P$ .

4.2 Convergence

A major advantage of SDAs—regardless of a particular starting value—is that they provide quadratic convergence at relatively low computational cost per iteration. We establish sufficient conditions for quadratic convergence in the following theorem.

Theorem 3 also applies to the cyclic reduction due to its equivalence to the SDA for SF2. Consequently, under Assumption 1, all algorithms presented above, from SDA for SF1 to logarithmic reduction, will converge quadratically to the unique (semi) stable solution $P$ .

4.3 Accuracy

To compare the different algorithms we need a measure of accuracy, that is, how far the solution produced by a particular method, $\hat {P}$ , is from the true solution, $P$ . As the latter is not known we will appeal to two measures, one from numerical stability analysis and another by comparing to arbitrary precision. For the numerical measure, we use the tighter forward error bound of Meyer-Gohde (Reference Meyer-Gohde2023)

(35) \begin{align} \underbrace {\frac {\left \lVert P-\hat {P}\right \rVert _F}{\left \lVert P\right \rVert _F}}_{\text{Forward Error}}\le \underbrace {\frac {\left \lVert H_{\hat {P}}^{-1}\text{vec}\left (R_{\hat {P}}\right )\right \rVert _2}{\left \lVert \hat {P}\right \rVert _F}}_{\text{Forward Error Bound}} \end{align}

where the residual of the UQME is $R_{\hat {P}}=A\hat {P}^2+B\hat {P}+C$ and $H_{\hat {P}}=I_{n_y}\otimes (A\hat {P}+B)+\hat {P}'\otimes A$ , and the forward error relative to a high precision solution of the problem, $\hat {P}_{HP}$ ,

(36) \begin{align} \underbrace {\frac {\left \lVert \hat {P}_{HP}-\hat {P}\right \rVert _F}{\left \lVert \hat {P}_{HP}\right \rVert _F}}_{\text{High Precision Forward Error}} \end{align}

The right-hand side of (35) can be rearranged so the numerator solves a Sylvester equation and we use the Advanpix Multiprecision Toolbox (Advanpix LLC., 2025) for $\hat {P}_{HP}$ .Footnote ¹² We will see in the applications that the two measures provide essentially the same answer, which favors (35) that operates directly with the numerical solution, $\hat {P}$ . Nonetheless, having both measures should support belay skepticism in either particular measure.

5.1 Smets and Wouters (2007) model

We start with Smets and Wouters’ (Reference Smets and Wouters2007) medium scale, estimated DSGE model that is a benchmark for policy analysis. They estimate and analyze a New Keynesian model with US data featuring the usual frictions, sticky prices and wages, inflation indexation, consumption habit formation as well as production frictions concerning investment, capital and fixed costs. For our purposes, the monetary policy rule is particularly important, as we will compare the accuracy of different solution methods when solving under alternate, but nearby parameterizations of the following Taylor rule

(37) \begin{gather} r_t=\rho r_{t-1} + (1-\rho ) \big(r_{\pi } \pi _t + r_{Y} (y_t - y_t^p) \big) + r_{\Delta y} \big(\big(y_t-y_t^p\big)-\big(y_{t-1}-y_{t-1}^p\big) \big) + \varepsilon ^r_t, \end{gather}

where $r_t$ is the interest rate, $\pi _t$ inflation, and $(y_t - y_t\,^p)$ the current output gap. The parameters $r_{\pi }$ , $r_{Y}$ , and $r_{\Delta y}$ describing the sensitivity of the interest rate to each of these variables, and also the change in the output gap, and $\rho$ measures interest rate smoothing.

The results for the posterior mode of Smets and Wouters (Reference Smets and Wouters2007) can be found in Table 1. Compared with the QZ algorithm, both the doubling algorithms, SF1 and SF2, along with the reduction algorithms already implemented in Dynare, perform very comparably. This is in contrast to the baseline Newton algorithm that fails to converge to the stable solvent (there is no guarantee that this algorithm will converge to a particular solvent though extensions of this baseline algorithm perform more favorably). The baseline Bernoulli algorithm gives a trade-off repeated throughout that analysis: generally more accurate, but often about an order of magnitude slower, than QZ. This is not the case with our doubling algorithms, they are as fast or faster than QZ and provide an order of magnitude more of accuracy. All four of the algorithms explored here, doubling and reduction, require about 10 iterations to converge and provide solutions that do not differ in an economic sense from the solution provided by QZ. Note that both of the measures of accuracy, the forward error bound in (35) and the error relative to a high precision solution (36), give similar results with the exception being the baseline Newton algorithm that solves the matrix quadratic, the forward error bound of (35) is small, but with a different solution than the rest of the algorithms—i.e., it converged to a solution containing unstable eigenvalues.

Table 1.

Smets and Wouters (Reference Smets and Wouters2007), posterior mode

Run Time for Dynare (QZ) in seconds, for all others, run time relative to Dynare. Max Abs. Diff. measures the largest abs. difference from the $P$ produced by Dynare.

In Table 2 we examine the different methods as solution refinement techniques, by parameterizing the model of Smets and Wouters (Reference Smets and Wouters2007) to a numerical instability following Meyer-Gohde (Reference Meyer-Gohde2023), see the online appendix for this parameterization, and initializing the different methods at the Dynare QZ solution. Note that the reduction methods cannot work with an initial value for the solution, so we compare the solution provided by Dynare’s QZ with the doubling refinement versions of SF1 and SF2 as well as the Newton and the Bernoulli method. First, consistent with theorem 2 and, hence, the irrelevance of initial conditions for that algorithm, the SF2 algorithm does not appear to be comparable with the Newton or SF1 algorithm, converging to a solution that has a much less substantial improvement in accuracy. As above, the two different measures of accuracy we examine agree on the relative ordering of the different algorithms. The Newton algorithm provides substantially more accuracy than the Bernoulli algorithm but at a higher additional computational cost, echoing a trade-off mentioned above. The SF1 algorithm breaks through this barrier, providing roughly the same level of accuracy as the Newton algorithm at less computational costs. That the Newton method is so time consuming despite the relatively few number of iterations performed emphasizes the computational intensity of the Newton step that is not shared by the doubling algorithms—the former involves solving structured linear (Sylvester) equations whereas the latter solve standard systems of linear equations. The resulting implied variances of inflation, $\pi _t$ , differ on an economically relevant scale. The more accurate methods all agree that the QZ solution understates the variance of inflation.Footnote ¹⁶

Table 2.

Smets and Wouters (Reference Smets and Wouters2007), numerically problematic parameterization

Run Time for Dynare in seconds, for all others, run time relative to Dynare. Variance $\pi _t$ gives the associated value for the population or theoretical variance of inflation.

We now take this refinement perspective and apply the different algorithms to solve iteratively for different parameterizations of the Taylor rule in the model of Smets and Wouters (Reference Smets and Wouters2007). We would like to establish whether solutions from previous, nearby parameterizations can be used to efficiently initialize the doubling methods similarly to the experiment above with the QZ solution as the initial guess. To that end, the experiment iterates through an evenly spaced $10\times 10$ grid for $r_\pi$ and $r_y$ , the responses to inflation and the output gap in the Taylor rule (37), holding all other parameters at their posterior mode values. We iterate through the grid taking the solution from the previous parameterization as the initialization for the next iteration, that is $P_{0}|_{r_\pi =r_{\pi ,j},r_y=r_{y,1}}=\hat {P}|_{r_\pi =r_{\pi ,j-1},r_y=r_{y,1}}$ initializing with $P_{0}|_{r_\pi =r_{\pi ,1},r_y=r_{y,1}}=\hat {P}^{QZ}|_{r_\pi =r_{\pi ,1},r_y=r_{y,1}}$ for the first row of the grid in the different values of $r_\pi$ and then $P_{0}|_{r_\pi =r_{\pi ,j},r_y=r_{y,j}}=\hat {P}|_{r_\pi =r_{\pi ,j},r_y=r_{y,j-1}}$ to fill out the grid in the values of $r_y$ . Thus we initialize all algorithms at the QZ solution and, as we move to nearby parameterizations, we measure the time and accuracy of the different methods. We calculate the median over the whole grid and then repeat the experiment with the grid points closer together by varying the size of the interval in the grid—setting $r_\pi \in [1.5,1.5 \times (1+10^{-x})]$ and $r_Y \in [0.125,0.125 \times (1+10^{-x})]$ , where $x \in [-1,8]$ . An increase in $x$ is a decrease in the spacing between the 100 grid points that thus increases the precision of the initialization.

Figure 1.

Forward errors and computation time per grid point for different parameterizations of the model by Smets and Wouters (Reference Smets and Wouters2007), log10 scale all axes.

Figure 1 summarizes the experiment graphically. Firstly, the left panel shows that the accuracy of the algorithms is independent of the grid spacing (and, hence, how close the parameter steps are from each other), this is expected and confirms the consistency of the algorithms and our measures of their accuracy.Footnote ¹⁷ The Bernoulli method is the exception and displays a significant drop in forward errors that coincides with the drop in computation time in the right figure—at a close enough parameterization, the Bernoulli algorithm starts with a guess from the previous parameterization that is accurate beyond the convergence threshold and the single iteration that is performed provides substantial (relative to more widely spaced grids) accuracy gains. In terms of their relative accuracies, the Newton is the most and the Bernoulli is generally (that is, except at the very closely spaced grids) the least accurate. All of the doubling and reduction algorithms are more accurate than QZ, with the cyclic reduction and SF2 algorithms roughly the same (as we would expect from Section 4.1) and close to Newton, followed by the log reduction, and finally SF1 doubling. As to be expected from Theorem 2, the SF2 doubling algorithm does not systematically profit in terms of reduced computation time as the parameter initializations improve (due to closer grid spacing) as the Newton and Bernoulli algorithms do. The SF1 doubling algorithm, in contrast, does profit with a downward trend in the computation time, although the Newton and Bernoulli algorithms clearly profit much more substantially with more significant downward trends.

5.2 MMB Suite Comparison

To guard against our comparison of methods from being specific to a particular model, we also adopt the model-robust approach inherit in the MMB. By abstracting from individual model characteristics we can gain insights into the algorithms that are robust to the particular features of a specific model. Besides drawing then general conclusions, we show an advantage for the SDAs with increasing model size. We implement version 3.3 of the MMB with the comparison repository Rep-MMB,Footnote ¹⁸ which contains over 160 different models, ranging from small scale to large-scale models and including estimated models of the US, EU, and multi-country economies. We test the algorithm on a subset of models appropriate for linear, rational expectations DSGE soution Footnote ¹⁹ —a histogram of the differing sizes can be found in Figure 2. With this large set of models, we feel confident in asserting that our assessment of the algorithms is robust to a representative sample of models relevant to policy settings.

Figure 2.

Histogram over the number of state variables for the 142 MMB models.

We solve each model in the MMB first initializing with the zero matrix and then with the QZ solution provided by Dynare. To get a reliable measure of the computation times, we solve each model 100 times and take the average computation time of the middle three quintiles to minimize the effects of outliers.Footnote ²⁰

Table 3.

MMB results

Run time and forward errors relative to Dynare, number of models that converge to the stable solution and number of iterations in absolute terms.

We begin with the initialization at the zero matrix and the top half of Table 3 summarizes the results. As noted in our previous studies, the Newton method is somewhat slower but more accurate than QZ when it converges to the stable solution—yet it does this in its baseline form only slightly more than half the time—and the baseline Bernoulli algorithm (almost) always converges to the unique stable solution—but does so more slowly than QZ at about the same level of accuracy. The cyclic reduction method is an intermediate case, converging to the stable solution for about 3/4 of the models, at generally comparable computation time, and somewhat higher accuracy than QZ. The logarithmic reduction is an improvement, converging for more models, and doing so at about the same speed and somewhat more accurately than the cyclic reduction. In terms of our doubling algorithms, they converge more frequently than even the logarithmic reduction (but still not for all and when, this failure is due to a breakdown of nonsingularity of the coefficient matrices in their recursions). As we will see in the next experiment, this can be overcome at least for the SF1 doubling algorithm by an appropriate (re)initialization of the algorithm.Footnote ²¹ The SF2 doubling algorithm on average outperforms QZ both in terms of speed and accuracy, although not for every model as the max or worst case shows (and again with the caveat that it does not successfully converge for 10 of the 142 models). The SF1 algorithm is not quite as fast at the median but has a far lower worst case computation time. This average performance of SF2 being faster than SF1 is consistent with the former inverting only one matrix, ${X}^{\dagger }_k - {Y}^{\dagger }_k$ —see Algorithm 2, whereas the latter needs to invert two, $I-{Y}_k {X}_k$ and $I-{X}_k{Y}_k$ —see Algorithm 1. Comparing the errors of the two doubling algorithms, SF2 has the lower worst case and mean forward error of the two and SF1 has the most accurate best case model and second best worst case model. Again, the information provided by the two error measures roughly coincides. In sum, the doubling algorithms (less so the reduction algorithms) perform favorably relative to QZ.

Figure 3.

Distribution of forward error bound for the MMB.

A graphical overview of the entire distribution of forward errors is plotted in Figure 3, the left relative to QZ and the right in absolute terms. The Newton algorithm is the most accurate with a clearly left shifted distribution relative to Dynare—with a mode improvement of about one order of magnitude—and the Bernoulli method the least with a clearly right shifted distribution relative to QZ. All of the methods here, reduction and doubling, lie in between but with modes and medians all to the left of QZ. While they are all very comparable, the SF2 algorithm has almost its entire mass to the left of Dynare and the cyclic reduction algorithm is right skewed with a substantial mass to the right of Dynare. The absolute values show that almost all of the models have errors less than 1e-10 and the accuracy of the Newton method is also clear with a substantial mass of forward error bounds inside machine precision and all mass essentially below e-12. The distributions of forward errors relative to the high precision solution (36) give virtually the same picture, as we show in Figure 4c. In any case, it appears that the models in the MMB were solved with acceptable accuracy and the SF2 doubling or logarithmic reduction is to be preferred among the algorithms here.

Figure 4.

Forward errors and computation time, log10 scales, for the MMB.

Figure 4 provides a model-by-model comparison of the different algorithms’ performance relative to QZ with each dot being the result of a model. All four panels express computation times and forward errors relative to Dynare on a log scale—hence a negative value in a dimension means the algorithm is more accurate or computationally efficient than QZ. In the top left panel, all of the methods are plotted using the forward error bound. First, recognize that the Bernoulli method is sometimes more and sometimes less accurate than QZ but usually slower (higher computation time), whereas the baseline Newton method is generally more accurate but usually slightly slower than QZ. The doubling and reduction algorithms require similar computation times as QZ (almost always within one half an order of magnitude slower/faster) but are generally more accurate. The doubling algorithms are more accurate than the reduction algorithms (notice the outlier along the x axis of the reduction algorithms with more than three orders of magnitude higher forward errors than QZ) and the SF2 doubling algorithm is the most accurate of the doubling algorithms as can be seen by comparing the lower two panels—visually, this algorithm is on average slightly faster than QZ and provides an order of magnitude more accuracy. The upper right panel plots the forward error relative to a high precision solution to the numerical forward error bound—as all points, hence the solutions from all algorithms for all models, lie along the 45 degree line, we conclude these two measures provide the same information. This provides strong evidence in favor of Meyer-Gohde (Reference Meyer-Gohde2023) practical forward error bound (35) that calculates this bound without having to calculate an additional, computationally expensive higher precision solution, but rather operates through the numerical stability analysis of the solution provided by the method in question only.

Figure 5.

Forward errors, computation time, and state variables for MMB.

Figure 5 continues the model-by-model comparison of the different algorithms’ performance relative to QZ, but now with a focus on the effect of model size, measured by the number of state or backward-looking variables. The top left panel gives the accuracy of the different methods for the different models plotted against the dimension of the endogenous state. There is not a visually compelling correlation—although for very large models, the reduction and doubling algorithms like the Newton algorithm appear to perform better than QZ. For computational times the story is different: the right panel plots the computational time against the number of state variables and a clear downward trend or negative correlation in particular for the doubling algorithms is obvious. In sum, for the largest models in the MMB, our doubling algorithms provide about two orders of magnitude more accuracy at about one tenth the computation time.

To assess the potential for improving on solutions, we repeat the exercise, but now initialize with the solution provided by QZ, see the bottom half of Table 3 for the results. The cyclic and logarithmic reductions are no longer included as they cannot operate with an arbitrary initialization. The baseline Newton method significantly improves the number of models it converges for. Note especially that, in contrast to the top half of the table when initialized at the zero matrix, now the SF1 doubling algorithm alongside the baseline Bernoulli method converges successfully to the unique stable solution in all models. Consistent with theorem 2, the SF2 doubling algorithm, however, does not and continues to converge or not for all but one of the same models as under the initialization with the zero matrix.

Both the Bernoulli and Newton methods run just a few iterations, and the later generally achieves a greater increase in accuracy albeit at a higher additional computation cost due to the computational intensity of the Newton step. The SF1 and SF2 doubling algorithms run through multiple iterations, ending up being faster and slightly more accurate than the Newton method at the median. Far more impressive here are the max or worst-case outcomes, with SF1 doubling at worst adding on an additional 7.5 times QZ computation cost and 4.9 times the forward error. The worst case Newton computation costs are an additional 840 times and Bernoulli 50+ thousand times the QZ initial time and the SF2 doubling algorithm has at worst a forward error bound about 300 times and Newton 1e9 times that of QZ. The results also highlight the effects of finite precision: numerical calculations generally introduce nonzero numerical errors, e.g. rounding errors, and if QZ’s solution is already accurate up to machine precision, additional calculations with any algorithm are unlikely to increase and more likely to decrease accuracy. We conclude that the SF1 doubling algorithm ought to be preferred as a solution refinement method, usually providing significant accuracy gains at low additional computation costs that is robust even in the worst case relative to alternatives.

Figure 6, like Figure 3 but now initialized at the QZ solution, provides an overview of the entire distribution of forward errors, the left relative to those from Dynare’s QZ method and the middle in absolute terms. We display the results from the forward error bound and, again, the forward error relative to a high precision solution gives the same picture, as confirmed by the figure on the right. It is apparent that the baseline Bernoulli algorithm provides only a marginal improvement on the QZ solution. While this should be tempered with the observation that only two Bernoulli iterations were performed at the median, the Newton algorithm also only runs a few iterations. The doubling algorithms perform favorably, with both displaying left shifted distributions of error bounds relative to QZ. Taken together with the results from Table 3, the SF1 doubling algorithm can be considered an ideal solution refinement algorithm across a wide variety of models.

Figure 6.

Distribution of FE bound for the MMB, initialized at dynare (QZ).