Computing $(A\otimes I + I \otimes A)^{-1} \text{vec}B$ Suppose I have two positive semi-definite $n$-by-$n$ matrices $A$, $B$ and an $n$-by-$n$ identity matrix $I$, and I'm looking for a way to compute, approximate or bound the following quantity:
$$(A\otimes I + I \otimes A)^{-1} \text{vec}B$$
Concretely I'm dealing with matrices with $n$ ranging from $100$ to $4000$, so $A$ is easy to invert, while $A \otimes I$ is too large, so need a way to compute this using operations on $n$-by-$n$ matrices
Additionally, I found the following to give a decent approximation when $A$, $B$ can be Kronecker-factored, wondering if there's a reason for that.
$$0.5 A^{-0.5} B  A^{-0.5}$$
Any tips or literature pointers are appreciated!
 A: Your equation is equivalent to solving the following equation for $P$
$$
A\,P + P\,A^\top = B. \tag{1}
$$
Such equation also solves for the infinite time (controllability) Gramian, which can also be calculated with
$$
P = -\int_0^\infty e^{A\,t}\,B\,e^{A^\top t}\,dt. \tag{2}
$$
This integral should also be of size $n \times n$. As time goes to infinity $e^{A\,t}$ is not well defined when $A$ is positive semi-definite. It can also be noted that in order to be able to solve for $P$ the spectra of $A$ and $-A$ need to be disjoint. In order words $A$ can't have an eigenvalue of zero, which would imply that $A$ should be positive definite instead of positive semi-definite. When multiplying both sides of $(1)$ by minus ones keeps it equivalent to the original problem. This is also equivalent to multiplying both $A$ and $B$ by minus one, which transforms $(2)$ into
$$
P = \int_0^\infty e^{-A\,t}\,B\,e^{-A^\top t}\,dt. \tag{3}
$$
Now the term $e^{-A\,t}$ does stay bounded and converges to zero as time goes to infinity. The rate of how fast this integral converges should be roughly inversely proportional to the smallest eigenvalue of $A$, however any finite time numerical approximation of this integral can be used as a lower bound for $P$.
A: Given symmetric matrices $\mathrm A \succ \mathrm O_n$ and $\mathrm B \succeq \mathrm O_n$, let
$$\mathrm X (t) := \int_0^t e^{-\tau \mathrm A} \mathrm B e^{-\tau \mathrm A} \,\mathrm d \tau$$
Using integration by parts,
$$\mathrm A \mathrm X (t) = \cdots = -\dot{\mathrm X} (t) + \mathrm B - \mathrm X (t) \mathrm A$$
and, thus, we obtain the following matrix differential equation
$$\boxed{\dot{\mathrm X} = - \mathrm A \mathrm X - \mathrm X \mathrm A + \mathrm B}$$
with initial condition $\mathrm X (0) = \mathrm O_n$. Assuming $\rm X$ converges to a steady-state, $\dot{\mathrm X} = \mathrm O_n$, we obtain the following linear matrix equation
$$\mathrm A \mathrm X + \mathrm X \mathrm A = \mathrm B$$
whose solution we denote by $\bar{\rm X}$. This is the desired solution.
Using a numerical ODE solver, we then integrate the matrix differential equation above and find an approximate value of matrix $\bar{\rm X}$. Since matrix $\rm X (t)$ is symmetric for all $t \geq 0$, we can half-vectorize both sides of the matrix differential equation in order to obtain a $\binom{n+1}{2}$-dimensional state vector (rather than a $n^2$-dimensional state vector, which would be wasteful).

Appendix
Note that the solution of the matrix differential equation is given by
$$\mathrm X (t) = \bar{\mathrm X} - e^{-t \mathrm A} \bar{\mathrm X} e^{-t \mathrm A}$$
Since $\mathrm A \succ \mathrm O_n$,
$$\lim_{t \to \infty} \mathrm X (t) = \bar{\mathrm X}$$
If matrix $\rm A$ were allowed to have an eigenvalue at zero, it would be unpleasant.
