Is there an efficient algorithm for matrix-valued continued fractions?

Question

Suppose I have a matrix equation recursively defined as

A[n] = inverse([1 - b[n]A[n+1]]) * a[n]

Then the equation for A[1] looks similar to a continued fraction, for which there are some highly efficient methods that avoid tedious recalculation (See "Numerical Recipes" for some examples).

However, I wonder if there are analogous methods that allow the coefficients b[n] and a[n] to be matrices, with the only constraint that b[n] A[n+1] be a square matrix so that the matrix

1 - b[n]A[n+1]

is actually invertible.

Answer 1

The following two methods are given in Functions of Matrices: Theory and Computation by Nicholas Higham, on page 81. These formulas evaluate

$$r(X) = b_0 + \frac{a_1 X}{b_1+\frac{a_2 X}{b_2+\cdots + \frac{a_{2m-1} X}{b_{2m-1} + \frac{a_{2m}X}{b_{2m}}}}} $$ where $X$ is a square matrix.

---

Top-down method:

$P_{−1}=I, Q_{−1}=0,P_0 =b_0 I,Q_0 =I$

for j=1:2m

$P_j =b_j P_{j−1}+a_jXP_{j−2}$

$Q_j =b_j Q_{j−1}+a_j X Q_{j−2}$

end

$r_m =P_{2m} Q_{2m}^{-1}$

---

Bottom-up method:

$Y_{2m} = (a_{2m}/b_{2m})X$

for j=2m−1:−1:1

Solve $(b_jI+Y_{j+1})Y_j =a_jX$ for $Y_j$.

end

$r_m=b_0I+Y_1$

---

The question asks for evaluation of the more general form

$$b_0 + \frac{a_1 X_1}{b_1+\frac{a_2 X_2}{b_2+\cdots + \frac{a_{2m-1} X_{2m-1}}{b_{2m-1} + \frac{a_{2m}X_{2m}}{b_{2m}}}}}$$

This can be evaluated by a simple generalization of the formulas above; for instance the bottom-up method becomes

$Y_{2m} = (a_{2m}/b_{2m})X_{2m}$

for j=2m−1:−1:1

Solve $(b_jI+Y_{j+1})Y_j =a_j X_j$ for $Y_j$.

end

$r_m=b_0I+Y_1$.

Answer 2

I know that this answer makes a lot of assumptions, but it at least generalizes your algorithm:

Suppose that $\{A_n\}$, $\{B_n\}$, and the seeding matrix, $V_N$, all form a commuting family of normal matrices, where the eigenvalue decompositions of $\{A_n\}$ and $\{B_n\}$ are known a priori, say $U' V_N U = \Lambda_N$, $U' A_n U = \Omega_n$, and $U' B_n U = \Delta_n$, where $U$ is unitary and $\Lambda_N$, $\{\Omega_n\}$, and $\{\Delta_n\}$ are complex-valued diagonal matrices.

Once we have said decomposition, by induction,

$$V_n = (I - B_n V_{n+1})^{-1} A_n = (I - U \Delta_n U' U\Lambda_{n+1} U')^{-1} U \Omega_n U',$$

which can be rearranged into the form

$$V_n = U (I - \Delta_n \Lambda_{n+1})^{-1} \Omega_n U' \equiv U \Lambda_n U',$$

where $\Lambda_n$ is of course still diagonal, so the entire family $\{V_n\}$ will necessarily commute with the other operators, and we have shown that the diagonal values of each $\Lambda_n$ are decoupled, so your fast scalar recursion formula can be applied independently on the eigenvalues of $V_N$ and the coefficient matrices.

Note that a special case is when $A_n \equiv \alpha_n I$ and $B_n \equiv \beta_n I$, so that the only requirement is that $V_N$ be a normal matrix.