如何在MIMO检测中用共轭梯度法替换LMMSE均衡器？

核心原理回顾 #

LMMSE均衡器的本质是求解线性方程组：
$$(H^H H + \sigma^2 I) \mathbf{x} = H^H \mathbf{y}$$
其中：

• $H \in \mathbb{C}^{N_r \times N_t}$：MIMO信道矩阵（$N_r$接收天线数，$N_t$发射天线数）
• $\mathbf{y} \in \mathbb{C}^{N_r \times 1}$：接收信号向量
• $\sigma^2$：接收端噪声方差
• $\mathbf{x} \in \mathbb{C}^{N_t \times 1}$：待检测的发射符号向量

共轭梯度(CG)法专门用于求解Hermitian正定矩阵的线性方程组，而$A = H^H H + \sigma^2 I$恰好满足正定条件，因此完全适配该场景。

具体实现步骤 #

1. 初始化系统与变量 #

• 确定MIMO系统参数：$N_t$（发射天线数）、$N_r$（接收天线数），生成信道矩阵$H$，接收信号$\mathbf{y}$，噪声方差$\sigma^2$。
• 构造系数矩阵$A$和右端向量$\mathbf{b}$：

  A = H.conj().T @ H + sigma2 * np.eye(Nt)
  b = H.conj().T @ y
  

• 初始化CG迭代变量：
  • 初始解$\mathbf{x}_0$：通常取全零向量（$\mathbf{x}_0 = \mathbf{0}$）
  • 初始残差$\mathbf{r}_0 = \mathbf{b} - A\mathbf{x}_0$
  • 初始搜索方向$\mathbf{p}_0 = \mathbf{r}_0$
  • 迭代计数器$k=0$，设置收敛阈值（如$10^{-6}$）和最大迭代次数（最多$N_t$次，CG理论上$N_t$次收敛）

2. 迭代求解循环 #

重复以下步骤直到残差小于阈值或达到最大迭代次数：

1. 计算步长$\alpha_k$：
   $$\alpha_k = \frac{\mathbf{r}_k^H \mathbf{r}_k}{\mathbf{p}_k^H A \mathbf{p}_k}$$
2. 更新解向量：
   $$\mathbf{x}_{k+1} = \mathbf{x}_k + \alpha_k \mathbf{p}_k$$
3. 更新残差：
   $$\mathbf{r}_{k+1} = \mathbf{r}_k - \alpha_k A \mathbf{p}_k$$
4. 计算方向更新系数$\beta_k$：
   $$\beta_k = \frac{\mathbf{r}{k+1}^H \mathbf{r}{k+1}}{\mathbf{r}_k^H \mathbf{r}_k}$$
5. 更新搜索方向：
   $$\mathbf{p}{k+1} = \mathbf{r}{k+1} + \beta_k \mathbf{p}_k$$
6. 迭代计数器加1，检查终止条件。

Python代码示例 #

import numpy as np

def cg_mimo_detector(H, y, sigma2, tol=1e-6, max_iter=None):
    Nr, Nt = H.shape
    if max_iter is None:
        max_iter = Nt  # CG最多Nt次收敛
    
    # 构造A和b
    A = H.conj().T @ H + sigma2 * np.eye(Nt)
    b = H.conj().T @ y
    
    # 初始化变量
    x = np.zeros((Nt, 1), dtype=np.complex128)
    r = b - A @ x
    p = r.copy()
    r_norm_sq = np.abs(r.conj().T @ r)[0, 0]
    
    for k in range(max_iter):
        if r_norm_sq < tol:
            break
        
        # 计算Ap
        Ap = A @ p
        # 步长alpha
        alpha = r_norm_sq / np.abs(p.conj().T @ Ap)[0, 0]
        # 更新解
        x += alpha * p
        # 更新残差
        r_new = r - alpha * Ap
        # 计算beta
        r_new_norm_sq = np.abs(r_new.conj().T @ r_new)[0, 0]
        beta = r_new_norm_sq / r_norm_sq
        # 更新搜索方向
        p = r_new + beta * p
        
        # 更新残差和其模平方
        r = r_new
        r_norm_sq = r_new_norm_sq
    
    return x

# 测试示例
if __name__ == "__main__":
    # 设置MIMO参数
    Nt = 4  # 发射天线数
    Nr = 6  # 接收天线数
    sigma2 = 0.1  # 噪声方差
    
    # 生成信道矩阵和发送符号
    H = np.random.randn(Nr, Nt) + 1j * np.random.randn(Nr, Nt)
    x_true = np.random.randint(0, 2, (Nt, 1)) + 1j * np.random.randint(0, 2, (Nt, 1))  # QPSK符号
    # 生成接收信号
    noise = np.sqrt(sigma2/2) * (np.random.randn(Nr, 1) + 1j * np.random.randn(Nr, 1))
    y = H @ x_true + noise
    
    # 用CG求解
    x_cg = cg_mimo_detector(H, y, sigma2)
    # 对比LMMSE解
    x_lmmse = np.linalg.inv(H.conj().T @ H + sigma2 * np.eye(Nt)) @ (H.conj().T @ y)
    
    print("CG解与LMMSE解的均方误差：", np.mean(np.abs(x_cg - x_lmmse)**2))

关键注意事项 #

• 收敛特性：在无数值误差的情况下，CG最多$N_t$次迭代就能得到精确解，实际中可根据系统要求设置提前终止阈值，减少计算量。
• 复杂度优势：直接求LMMSE的复杂度为$O(N_t^3)$，而CG每次迭代仅需$O(N_r N_t)$的运算量，当$N_t$较大时，CG的效率远高于直接求逆。
• 数值稳定性：由于$A$是正定矩阵，CG迭代过程中不会出现发散问题，无需额外的正则化操作（除非噪声方差极小，但此时$\sigma^2 I$已经提供了正则化）。

内容的提问来源于stack exchange，提问作者Chaima Hmidi