Understanding an algorithm for computing a matrix polynomial

I'm trying to understand this algorithm by Charles Van Loan for evaluating a matrix polynomial $p(\mathbf A)=\sum\limits_{k=0}^q b_k \mathbf A^k=b_0\mathbf I+b_1\mathbf A+\cdots$ (where $\mathbf A$ is $n\times n$):

$\displaystyle s=\lfloor q\rfloor, \quad r=\lfloor q/s\rfloor$
$\displaystyle \mathbf Y=\mathbf A^s$
$\displaystyle \text{for}\quad j=1,\dots,n$
$\displaystyle\quad\quad \mathbf y_0^{(j)}=\mathbf e_j$ ($\mathbf e_j$ is the $j$th column of the identity matrix)
$\displaystyle\quad\quad \text{for}\quad k=1,\dots,s-1$
$\displaystyle\quad\quad\quad \mathbf y_k^{(j)}=\mathbf A\mathbf y_{k-1}^{(j)}$
$\displaystyle\quad\quad \mathbf f_0^{(j)}=\sum_{k=q}^{rs}b_k\mathbf y_{k-rs}^{(j)}$
$\displaystyle\quad\quad \text{for}\quad k=1,\dots,r$
$\displaystyle\quad\quad\quad \mathbf f_k^{(j)}=\mathbf Y\mathbf p_{k-1}^{(j)}+\sum_{i=0}^{s-1}b_{s(r-k)+i}\mathbf y_{i}^{(j)}$

I've tried going through it a number of times, but it seems that the quantity (vector?) $\mathbf p_{k-1}^{(j)}$ was never defined anywhere in the paper. What might it be? Also, the paper claims (or more probably it's how I (mis)understood the paper) that only three matrices need to be stored for the evaluation: $\mathbf A$, $\mathbf Y$, and the final $p(\mathbf A)$. However, I can't see how to properly organize things so that an extra array for storing the $\mathbf y_k^{(j)}$ is not needed (otherwise, four matrices are now needed instead of the three claimed in the paper).

Could anyone enlighten me on how to implement this algorithm properly? Better yet, is there any code available anywhere that implements this algorithm (I'm language-agnostic, so whatever code you post, I can likely translate into the language I'm currently using).

Thanks for any help!

-
For those who can't access IEEE, here is a mirror on Van Loan's home page. – Jerry Jun 30 '11 at 8:22
Since in the paper is said $p(A) = [f_r^{(1)}|\cdots|f_r^{(n)}]$, I would guess $p_{k-1}^{(j)}$ should be $f_{k-1}^{(j)}$. Otherwise why bother to compute $f_k^{(j)}$ for $k \neq r$? – Giacomo d'Antonio Jun 30 '11 at 9:57
And you're probably right on the fourth matrix (but I also didn't read the paper carefully). I guess the author is very good in "selling smoke". – Giacomo d'Antonio Jun 30 '11 at 10:04
@Giacomo: I don't want to judge Van Loan too harshly; after all he is the co-author of a certain useful reference for matrix computations... – Jerry Jul 1 '11 at 2:16

The storage required for the $y^{(j)}_k$ is less than a "full" matrix in proportion as $s \ll n$. That is, $s$ such vectors are required for each step $j$, and the storage for one step may be reused in the next step.

The Question states $s$ to be the floor of $q$, the degree of the polynomial, but this is mistaken. Per the discussion of the paper preceding Algorithm I, $s$ can be any integer strictly between 1 and $q$. While this does not guarantee $s \ll n$, it allows for it. See the discussion of applying Algorithm II to the Taylor polynomial for $e^A$, where $s=4$ and $q=16$ were chosen. [Note also the discussion of minimizing "work count" in Algorithm II by choosing $s$ around $\sqrt{q/2}$.]

That leaves still the mysterious $p^{(j)}_k$. My guess is it is a typo, that $f^{(j)}_k$ is being computed by adjusting $Yf^{(j)}_{k-1}$, after initializing $f^{(j)}_0$, but I need to do a bit of scribbling to confirm that works out.

Added: My guess above is correct, that we should read $f^{(j)}_{k-1}$ instead of $p^{(j)}_{k-1}$ at one point in Algorithm II where the latter appears.

The algorithm may be understood as combining two ideas. One is that the matrix polynomial $p(A)$ will be evaluated column by column, i.e. by applying it successively to unit vectors $e_j$, $j=1,..,n$.

The other idea is simple though a bit cumbersome to spell out. Recall that $1 \lt s \lt q$ where $q$ is the degree of $p$. Let the polynomial:

$$p(x) = \sum^q_{i=0} b_i x^i$$

be grouped according to the largest power of $x^s$ that divides each term and rewrite $p(x)$ in an expansion involving powers of $x^s$ and "coefficients" that are polynomials of degree at most $s-1$. Taking $r = \lfloor q/s \rfloor$, and assuming any indicated coefficients $b_m$ where $m \gt q$ are zero:

$$p(x) = \sum^r_{k=0} x^{sk} \beta_k(x)$$

where $\beta_k(x) = \sum^{s-1}_{i=0} b_{sk+i} x^i$ are the polynomial "coefficients" as above.

Assume that $Y = A^s$ has already been computed, e.g. by binary exponentiation or other scheme.

The outermost loop of Algorithm II is then iteratively applying $p(A)$ to $n$ unit vectors $v = e_j$ as explained above. For simplicity we just describe computing $p(A)v$ for an arbitrary vector $v$ of compatible dimension, a column of length $n$ where $A$ is $n \times n$. In effect we thus omit superscript $(j)$ from the remainder of the discussion.

First is a loop to initialize $y_0 = v$ and $y_i = A y_{i-1}$ for $i = 1,..,s-1$.

Second is a loop that essentially applies Horner's method to get $p(A)v$ using the expansion in powers of $A^s$ and "coefficients" $\beta_k(A)v$. Note that the latter are column vectors, and for $k = 0,..,r$ we evaluate $\beta_k(A)v$ by taking the appropriate linear combination of previously initialized $y_i$. That is:

$$f_0 = \beta_r(A) v$$ $$f_k = Y f_{k-1} + \beta_{r-k}(A) v$$

for $k = 1,..,r$. This shows us that the term $p^{(j)}_{k-1}$ is a misprint for $f^{(j)}_{k-1}$ in Algorithm II.

-
This is interesting... if you tried it out in some computing environment, would you mind sharing the code? – J. M. Aug 11 '11 at 21:30
@J. M.: I'd be glad to share. I suspect for many readers code would make Van Loan's idea appear in a more natural light. At first I was thinking just a vanilla C implementation (for the sake of being self-contained), but it occurs to me an Octave implementation would be more concise, if perhaps not ideal for checking running times. – hardmath Aug 12 '11 at 20:52
Yes, Octave code (which is essentially MATLAB ready) would be lovely. No rush, though. :) Thanks in advance! – J. M. Aug 13 '11 at 6:51
Hey, just wondering: did you ever manage to come up with neat code for Van Loan's method? – J. M. Oct 30 '11 at 16:10