## 2 System Model and Problem Formulation  
### System Model  
As depicted in Figure 1, we consider a MIMO system consisting of a transmitter equipped with \(M_{t}\) antennas, a receiver equipped with \(M_{r}\) antennas, and a total of \(K\) interferers, where the \(k\)th interferer has \(M_{k}\) antennas, \(k=1,\cdots,K\). It is assumed that the background additive noise is much weaker than the interference, and hence we ignore the additive noise term for convenience (this is valid in practice as the noise power ranges from \(-192.5\) dBm\(/\)Hz to \(-174\) dBm\(/\)Hz, whereas the interference signal power ranges from \(-100\) dBm to \(-10\) dBm). In the channel training procedure, multiple training symbols are sent from the transmitter during \(L\) symbol times. The received signal matrix is then given by  
\[\mathbf{Y}=\mathbf{H}\mathbf{P}+\underbrace{\sum_{k=1}^{K}\mathbf{H}_{k} \mathbf{S}_{k}}_{=\mathbf{N}}=\mathbf{H}\mathbf{P}+\mathbf{N}, \tag{1}\]  
where \(\mathbf{P}\in\mathbb{C}^{M_{t}\times L}\) and \(\mathbf{S}_{k}\in\mathbb{C}^{M_{k}\times L}\) represent the training signal matrix sent from the transmitter and the interfering signal matrix sent from the \(k\)th interferer, respectively. We simply denote the total interference as \(\mathbf{N}\). The channel matrix between the transmitter and receiver is denoted by \(\mathbf{H}\in\mathbb{C}^{M_{t}\times M_{t}}\) and that between the \(k\)th interferer and receiver by \(\mathbf{H}_{k}\in\mathbb{C}^{M_{t}\times M_{t}}\). Without loss of generality, it is assumed that there exists spatial correlation at all nodes. According to the well-known Kronecker model [25], we represent the channels and interfering signals as follows (the analysis in our work is not confined to a specific distribution of the channel matrices, but valid for any distribution, as the proposed method requires only the knowledge of the covariance information of the channel matrices, not their probability distribution):  
\[\mathbf{H}=\mathbf{R}^{1/2}\mathbf{H}_{w}\mathbf{T}^{1/2}, \tag{2}\]  
\[\mathbf{H}_{k}=\mathbf{R}^{1/2}\mathbf{H}_{w\lambda}\mathbf{T}^{1/2}_{1,k}, \quad k=1,2,\cdots,K, \tag{3}\]  
\[\mathbf{S}_{k}=\mathbf{\Psi}^{1/2}_{s,k}\mathbf{S}_{w\lambda}\mathbf{\Psi}^{1/ 2}_{\tau,k},\quad k=1,2,\cdots,K, \tag{4}\]  
where the elements of the matrices \(\mathbf{H}_{w}\), \(\{\mathbf{H}_{w\lambda}\}_{k=1}^{K}\), and \(\{\mathbf{S}_{w\lambda}\}_{k=1}^{K}\) are independent and identically distributed (i.i.d.) as \(\mathcal{CN}(0,1)\). \(\mathbf{T}\succeq\mathbf{0}\), \(\mathbf{R}\succeq\mathbf{0}\), and \(\mathbf{T}_{1,k}\succeq\mathbf{0}\) are the spatial correlation matrices at the transmitter, receiver, and \(k\)th interferer, respectively. \(\mathbf{\Psi}_{\tau,k}\succeq\mathbf{0}\) and \(\mathbf{\Psi}_{s,k}\succeq\mathbf{0}\), respectively, represent the temporal and spatial correlations of the interfering signal. Taking the vectorizing operation on both sides of (1), the received signal can be rewritten in vector form as  
\[\text{vec}(\mathbf{Y})=\mathbf{y}=\left(\mathbf{P}^{T}\otimes\mathbf{I}_{M_{ r}}\right)\mathbf{h}+\mathbf{n}, \tag{5}\]  
where  
\[\mathbf{h}=\text{vec}(\mathbf{H})\]  
and  
\[\mathbf{n}=\text{vec}(\mathbf{N})=\sum_{k=1}^{K}(\mathbf{I}_{L}\otimes \mathbf{H}_{k})\text{vec}(\mathbf{S}_{k}).\]  
For the design purpose, we assume that the channel and interference covariance matrices are perfectly known at the receiver side (The issue with imperfect covariance information at the receiver side was studied in [20]. In [20] (Remark 2), it was shown that when the covariance information is imperfect at the receiver side, the resulting channel estimation MSE has a similar mathematical expression to (12), with an additional (negligible) loss term. In [20], it was also empirically demonstrated that the imperfect covariance information at the receiver side has a negligible impact on the performance, whereas that at the transmitter side has a much more significant impact).