I am trying to solve a system of equations which can be represented in the well-known form of $\bf Ax = b$, where $\bf A$ is a large sparse matrix consisting of known values, $\bf x$ is a vector of unknowns and $\bf b$ is a vector of known values. The matrix $\bf A$ is obtained as follows
$$ {\bf A} := \left( {\bf I}^\top \otimes {\bf A}_1 \right) + \left( {\bf I}^\top \otimes {\bf A}_2 \right) + \operatorname{diag} ( \operatorname{vec} \left({\bf A}_3\right)) \left( {\bf I}^\top \otimes {\bf A}_4 \right)$$
where ${\bf A}_1$, ${\bf A}_2$, ${\bf A}_3$ and ${\bf A}_4$ are all known matrices, and $\bf I$ is the identity matrix of conforming dimensions. Now, computing the Kronecker product while viable for small matrices(size in the order of a few hundreds), is obviously not feasible for very large systems.
The only approach I could think of is inverting this equation to have ${\bf x} = {\bf A}^{-1} {\bf b}$ and then applying the inverse individually to each term in the aforementioned equation and then basically computing the Kronecker-vector product without explicitly computing the Kronecker product itself. However, for the matrices I have, they are singular or nearly singular in nature.
Is there a known set of methods or literature that can be used to solve these types of problems without computing the Kronecker product explicitly?
Best Answer
Many iterative solvers (splitting-based methods such as power method, Jacobi method, and Gauss-Seidel method, Krylov subspace methods, multi-level methods etc.) can be modified so that they do not require explicit generation of the Kronecker products. Application of these solution techniques to Markov chains discussed in various articles and books. Analyzing Markov chains using Kronecker products: theory and applications (Dayar, 2012) is a good introduction book for the topic of solving Markov chains using Kronecker products. Therein, Kronecker representation and details of the solvers are discussed, and many references to related research are given.
In most of these solvers, vector-matrix multiplication is replaced by vector-Kronecker product multiplication so that elements of the large matrix are not stored. In vector-Kronecker product multiplication algorithms,the multiplication of matrix $A = \otimes_{j = 1}^{J} A_j$ and vector $n$, is computed without explicitly storing $A$. A detailed discussion of these algorithms can be found in Vector-Kronecker product multiplication As a result, vector-matrix multiplication in a solver can be replaced with vector-Kronecker product multiplication where the Kronkecker product is not explicitly generated.
Jacobi method Let $A = D + R$ where $D$ and $R$ respectively include diagonal and off-diagonal elements of $A$. Note that, the inverse of $D$, that is $D^{-1}$, is given by replacing the main diagonal elements of the matrix with their reciprocals. Since $D^{-1}$ is used in the iterations, it is more efficient to store elements of $D^{-1}$ instead of elements of $D$.
The solution for $Ax = b$ can be obtained iteratively using Jacobi iterations $$ x^{(k)} = D^{-1} \left( b - R x^{(k-1)} \right) ~~~~~~ \text{ for } k = 1,2,\ldots $$ where $x^{(0)}$ is the initial guess and $x^{(k)}$ is the approximation at $k$th iteration for $k = 1,2,\ldots$. The vector $r^{(k-1)} = R x^{(k-1)}$ is computed using a vector-matrix multiplication algorithm. Then, elementwise multiplication of vectors including $D^{-1}$ and $(b - r^{(k-1)})$ gives $x^{(k)}$.
Now, assume that $R$ is represented as sums of Kronecker products. Then, in a Kronecker-based Jacobi solver, the vector $r^{(k-1)} = R x^{(k-1)}$ is computed using a vector-Kronecker multiplication algorithm instead of a vector-matrix multiplication algorithm.
The approach is similar in Krylov subspace method, that is the vector-matrix multiplication is computed without explicitly storing the large matrix by employing a Kronecker-based vector-matrix multiplication algorithm.
Note In the matrix you have asked, there is a diagonal matrix in the expression. The solution methods needs to be modified so that that diagonal matrix is taken into account. You may simply store the elements of the elements of that matrix in a vector and element-wise multiply with the intermediate vector when necessary.