Covariance estimation for matrix-valued data has been a hot topic in applications. A new method was developed that doesn't rely on specific assumptions about the data's distribution or size. By breaking down the original covariance matrix into smaller ones, the method can estimate the variability in different directions. The method was shown to be optimal in certain situations and can handle heavy-tailed data. It was tested on temperature and stock data, showing better performance than existing methods.