Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

What is the most efficient way to make this linear algebra computation? I am interested in computing a vector $y^{(k)}$ that updates as shown below. $$y^{(k)} = A^k B A^k x$$ where the matrices $A,B \in \mathbb{R}^{n \times n}$ are tridiagonal matrices and the vectors $x,y^{(k)} \in \mathbb{R}^{n \times 1}$. The matrices $A$, $B$ and $x$ are known.

For instance, $y^{(0)} = Bx$, $y^{(1)} = ABAx$, $y^{(2)} = A^2 B A^2x$ and so on.

Recall that a tri-diagonal matrix vector product costs $\mathcal{O}(n)$.

A naive operation would cost $\mathcal{O}((2k+1)n)$ at step $k$. It is possible to cut down the cost by $1/2$, if we define $z^{(k)} = A^{(k)}x$ and at each step compute $z^{(k)}$ by the recurrence $z^{(k+1)} = A z^{(k)}$, which costs $\mathcal{O}(n)$ and then compute $y^{(k)}$ as $y^{(k)} = A^k B z^{k}$, which costs $\mathcal{O}((k+1)n)$ and hence the total cost would be $\mathcal{O}((k+2)n)$ as opposed to $\mathcal{O}((2k+1)n)$.

But I want to do better than the above, by cutting down the dependence of $k$ at each time step.

Obviously I do not want to compute and store $A^k B A^k$ at each time step since the matrix will keep increasing its bandwidth making the computation and storage only more and more expensive.

Ideally, what I am looking for is to make use of some recurrence for $y^{(k)}$'s so that the update at each step $k$ can be done in $\mathcal{O}(n)$ instead of being a function of $k$ as well. But constructing such a nice recurrence, which can be computationally advantageous has eluded me. (The main reason why I asking is that $k$ also gets large as I proceed.)

share|cite|improve this question

Let $b_i$ be the $i$th column vector of $B$. Then, write $b_i = x + \hat{b}_i$. This is very simple to compute, particularly when $B$ is tri-diagonal.

Then, $$A^kB = A^k \begin{pmatrix} x+\hat{b}_1 & x+\hat{b}_2 & \cdots & x+\hat{b}_n \end{pmatrix} = A^k X+A^k\hat{B},$$ where $X$ is a matrix whose column vectors are all $x$, and $\hat{B}$ is a matrix whose column vectors are $\hat{b}_i$.

Then, you compute $A^kx$ once, and you can obtain $A^kX$ for "free." Finally, you can store each element-wise product of each row of $A^kx$, and compute only the three relevant products in $A^k\hat{B}$.

For $A^k$ full, $A^kx$ has $n^2$ multiplications. Store these. Then, for $B$ tri-diagonal, $\hat{b}_i$ has at most three elements that are distinct from elements of $-x$. Therefore, it should require no more than $3n$ new multiplications to compute $A^k\hat{B}$. Finally, the matrix vector product $(A^kX+A^kB)(A^kx)$ has $n^2$ multiplications, making the total op-count less than $2n^2+3n$. This should be better than 2 matrix-matrix multiplications --assuming you can compute $A^k$ efficiently at each step.

share|cite|improve this answer

You might try,

$y^{(k)} = M^{(k)} x$

where $M^{(k)} = A M^{(k-1)} A$ and $M^{(0)} = B$. This way, each successive step only costs one matrix-vector multiplication and 2 matrix-matrix multiplications.

share|cite|improve this answer
Thanks for the reply. Ideally, I do not want to perform these matrix matrix operations. I just want to use these matrices as a black-box for matrix vector products. If nothing works out though, I will probably resort to something like this. The problem with this is that the bandwidth of $M^{(k)}$ keeps increasing with $k$. Hence, the matrix vector product again scales only as $\mathcal{O}(kn)$. – user17762 Nov 4 '12 at 3:13

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.