Golub and Van Loan 3rd Ed Basic Algos Getting Started

Question

I am trying to get my brain back into matrix algebra, and I am starting right from the beginning of the standard Matrix Computations, by Golub and Van Loan ( 3rd Edition ). A couple of the practice problems in chapter 1.1 have to do with writing algorithms for matrices raised to a power. I don't see these cases explicitly described in the text, and I'm having trouble connecting the dots here. Specifically, problem P1.1.3:

P1.1.3 Give an algorithm for computing $C = (xy^T)^k$ where $x$ and $y$ are $n$-vectors.

So I understand the $xy^T$ is an outer product, so $C$ is going to be an $n$-by-$n$ matrix. I can easily use one of the matrix multiplication update algorithms for the case where $k=2$. For example, for the a algorithm

$$C=AB + C$$

where $A,B = xy^T$, then $A(i,j),B(i,j)=x^{}_{i}y^{}_{j}$. So the matrix update algorithm would be:

for i = 1:n:
    for j = 1:n:
        for p = 1:n:
            C(i,j) = x[i]*y[p]*x[p]*y[j]

I used the variable $p$ here since $k$ has a different meaning in the question. There are a few other variations of the matrix multiplication algorithm that could be substituted. But the big thing I'm not getting is how to expand this algorithm to account for $k>2$ in $(xy^T)^k$? It seems like there needs to be another loop in the algorithm to account for the $k$. But also I have this gut feeling that I'm supposed to be using outer product updated algorithm, which takes the form $A=A+xy^T$.

Big thing I'm missing is seeing how these update algorithms correspond to raising matrices to powers, or multiplying many times. Any help is greatly appreciated, thanks.

UPDATE

After considering the suggestion in the comments, I think I have it figured out. So as in the comment, $$C=(xy^T)^k$$ can be re-written as $$C=(y^Tx)^{k-1}(xy^T)$$ where $y^Tx$ is a dot product, so this represents a scalar. $y^Tx$ is an outer product matrix, so this operation simplifies to a scalar-matrix multiplication. From the text we could use a simple matrix-scalar algorithm to accomplish this.

Let $${\alpha} = (y^Tx)^{k-1}$$ $$A = xy^T$$ Then $$C={\alpha}A {\implies} c^{}_{ij}={\alpha}a^{}_{ij}$$

Answer 1

I think the purpose of the question is to let us take advantage of associativity of matrix product.

A naive way of compute this would be to compute $xy^T$ explicitly. This involves complexity $O(n^2)$ and then multiply $xy^T$ to the power $k$ directly, $O(n^3)$ each time, we are doing this $k$ times, hence the overall complexity would be $O(kn^3)$.

A smarter way has been pointed out in the comment, that is first we compute $y^Tx$, first, complexity $O(n)$, then raise to the power $k-1$, and then multiply it to the product of $xy^T$. Overall complexity would be $O(\max(n^2,k))$.

Remark about your attempt:

It is not clear to me what do you mean by $A,B = xy^T$ and in your for loop, only the last $p$ matters because you are overwriting it. I can't tell what is going on though.