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.

  • $\begingroup$ This is the question you asked in math.SE a few months before, no? Is $A$ square or rectangular? $\endgroup$
    – J. M.
    Nov 30, 2011 at 7:01
  • $\begingroup$ I recall that someone in the comments at math.SE suggested I ask this here once the beta is online :) In my special case, A is rectangular. The recursive equations correspond to a hierarchical set of equations, and the number of quantities grows with $n$. In my case, the dimension of A[n] is n x (n-1) $\endgroup$
    – Lagerbaer
    Nov 30, 2011 at 19:22
  • $\begingroup$ Just curious, what is the application you want to use this for? $\endgroup$
    – Hjulle
    Dec 3, 2011 at 20:34
  • 1
    $\begingroup$ Very briefly, using Dyson's identity for a particular Hamiltonian generates Green's functions that I can label with a certain index $N$. Collecting all functions with the same index into a vector $V_N$ allows me to write $V_N = \alpha_N V_{N-1} + \beta_N V_{N+1}$ using Dyson's identity and a suitable approximation. Using a cut-off so that $V_N = 0$ for all $n \ge N$ allows me to find matrices $A_n$ so that $V_n = A_n V_{n-1}$ and these matrices are given by my continued-fraction style equation. This technique can, for example, compute lattice Green's functions for tight-binding models. $\endgroup$
    – Lagerbaer
    Dec 3, 2011 at 21:04
  • 1
    $\begingroup$ It's not my field, but i was a while back at a seminar where something pertaining to this problem was presented. [Here][1] is the only trace i could find of it online. I really don't know whether it helps. [1]: mh2009.ulb.ac.be/ResActiv.pdf $\endgroup$
    – user189035
    Jul 4, 2012 at 9:46

2 Answers 2


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}$


$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$.



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$.



  • $\begingroup$ This looks very interesting. I will see if I can apply it to my specific problem but it answers the question since my b[n] * A[n+1] is a square matrix $\endgroup$
    – Lagerbaer
    Feb 29, 2012 at 3:53
  • $\begingroup$ Ah, but I just noticed that the matrix $X$ is the same everywhere in your solution, but mine isn't necessarily. $\endgroup$
    – Lagerbaer
    Feb 29, 2012 at 3:54
  • $\begingroup$ Okay, I have generalized it. $\endgroup$ Feb 29, 2012 at 14:07

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.