Solving powered matrix using LU decomposition

Question

My task is to solve using $LU$ decomposition $A^k x = b$ where $k$ is a positive integer.

I can raise the matrix $A$ to the power of $k$ and use $LU$ decomposition to solve the linear system $A^k x = b$. While this solution is correct, I was told that I should try to use only one $LU$ decomposition along with the identity $A^k = LU \cdots LU$. I understand this identity but do not know how to proceed.

Answer 1

Note that the system $A^k x = b$ can be written $LUA^{k-1}x=b$ by factoring the first matrix.

1. Let $y_1 \equiv UA^{k-1} x$ so that $Ly_1 = b$. Solve for $y_1$ using a forward solve.
2. Now you are left with the system $UA^{k-1} x = y_1$. Let $z_1 \equiv A^{k-1} x$ so that $U z_1 = y_1$. Solve for $z_1$ using a backward solve.
3. Now you are left with the system $A^{k-1} x = z_1$. If $k = 1$, you are done and your answer is $x = z_1$. Otherwise, rewrite the system as $LUA^{k-2} x = z_1$ and repeat the above procedure.

Let $n$ be the order of the matrix $A$.

Runtime. Since a forward/backward solve costs $O(n^2)$, this algorithm requires $O(n^3)$ FLOPs to perform the $LU$ factorization and $O(n^2 k)$ FLOPs to perform the $k$ forward/backward solve pairs.

Naive algorithm. Alternatively, we could have computed $A^k$ and solved the system $A^k x = b$ directly (i.e., Gaussian elimination). Taking the power naively by the divide and conquer approach requires $O(n^3 \lg k)$ FLOPs. $O(n^3)$ FLOPs are needed for Gaussian elimination.

• The method of repeated forward/backward solves is desirable over the naive algorithm when $k$ is strictly sublinear in $n$. Even when $k$ is linear in $n$, the algorithm still has a bit of an edge over the naive algorithm (due to the $\lg k$ term).
• If $k$ is strictly superlinear in $n$, the naive algorithm is the more desirable of the two.