18

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.

Lagerbaer
  • 487
  • 4
  • 9
  • This is the question you asked in math.SE a few months before, no? Is $A$ square or rectangular? – J. M. Nov 30 '11 at 07:01
  • 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)

     – Lagerbaer Nov 30 '11 at 19:22
  • Just curious, what is the application you want to use this for? – Hjulle Dec 03 '11 at 20:34
  • 1
    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. – Lagerbaer Dec 03 '11 at 21:04
  • You might reconsider whether this question is suitable for cstheory.stackexchange.com – shuhalo Dec 08 '11 at 16:16
  • 1
    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]: http://mh2009.ulb.ac.be/ResActiv.pdf – user189035 Jul 04 '12 at 09:46

2 Answers2

9

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

David Ketcheson
  • 16,522
  • 4
  • 54
  • 105
6

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.

Jack Poulson
  • 7,599
  • 32
  • 40