This study explores the estimation of parameters in a matrix-valued linear regression model, where the $T$ responses $(Y_t)_{t=1}^T \in \mathbb{R}^{n \times p}$ and predictors $(X_t)_{t=1}^T \in \mathbb{R}^{m \times q}$ satisfy the relationship $Y_t = A^* X_t B^* + E_t$ for all $t = 1, \ldots, T$. In this model, $A^* \in \mathbb{R}_+^{n \times m}$ has $L_1$-normalized rows, $B^* \in \mathbb{R}^{q \times p}$, and $(E_t)_{t=1}^T$ are independent noise matrices following a matrix Gaussian distribution. The primary objective is to estimate the unknown parameters $A^*$ and $B^*$ efficiently. We propose explicit optimization-free estimators and establish non-asymptotic convergence rates to quantify their performance. Additionally, we extend our analysis to scenarios where $A^*$ and $B^*$ exhibit sparse structures. To support our theoretical findings, we conduct numerical simulations that confirm the behavior of the estimators, particularly with respect to the impact of the dimensions $n, m, p, q$, and the sample size $T$ on finite-sample performances. We complete the simulations by investigating the denoising performances of our estimators on noisy real-world images.