Solving powered matrix using LU decomposition 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.
 A: Note that the system $A^k x = b$ can be written $LUA^{k-1}x=b$ by factoring the first matrix.


*

*Let $y_1 \equiv UA^{k-1} x$ so that $Ly_1 = b$.
Solve for $y_1$ using a forward solve.

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

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

