含特殊稀疏矩阵P的矩阵运算计算复杂度求解

先明确矩阵P的特殊结构 #

P是M×N稀疏矩阵，由N阶单位矩阵的M个不同行采样得到，因此：

• P的每行恰好有一个1，其余元素为0，1的位置对应采样的行索引picks；
• 转置矩阵P^T是N×M矩阵，每列恰好有一个1，对应picks中的索引位置。

a) 计算 $H = \beta I + P^{T}P$ 的复杂度 #

常规N×M与M×N矩阵乘法的复杂度是O(N²M)，但利用P的结构可以大幅简化：

1. P^T P的本质：P^T的第k列是单位向量$e_{picks[k]}$，因此 $P^T P = \sum_{k=1}^M e_{picks[k]} e_{picks[k]}^T$，结果是一个对角矩阵——对角线上，索引在picks中的位置值为1，其余位置为0。
2. 构造H的步骤：
   • $\beta I$是对角矩阵，对角线元素全为$\beta$，构造其对角线的复杂度为O(N)；
   • 给picks对应的对角线位置加1，复杂度为O(M)；
3. 整体复杂度：O(N + M)，远低于常规矩阵乘法的O(N²M)。

对应代码中的get_H_hat_diag就是利用这个逻辑：通过P_diag标记picks位置，再用克罗内克积np.kron构造二维形式的对角元素数组，避免了矩阵乘法。

b) 计算 $a = F^{-1}Z + P^{T}B$ 的复杂度 #

第一项：$F^{-1}Z$（逆快速傅里叶变换）

Z是N×N矩阵，二维逆FFT的复杂度为O(N² log N)（对每行做O(N log N)的FFT，共N行；再对每列做同样操作，总复杂度为O(N² log N)）。

第二项：$P^{T}B$

常规矩阵乘法复杂度是O(N²M)，但利用P的结构简化：

• P^T乘以B的结果是N×N矩阵，其中：
  • 若行索引i在picks中，对应B中第k行（k为i在picks中的位置）直接复制到结果的第i行；
  • 若行索引i不在picks中，结果的第i行为全0。
• 构造这个矩阵的复杂度：初始化全0矩阵为O(N²)，复制B的M行到对应位置为O(M×N)，总复杂度为O(N² + M×N)。

整体复杂度

由于M≤N（P从N阶单位矩阵采样M行），M×N≤N²，而N² log N远大于N²，因此整体复杂度由逆FFT主导，为O(N² log N)，而非你之前认为的O(N²M)。

对应代码中a_hat[np.ix_(picks, picks)] = B就是直接将B放入对应位置，实现P^T B的二维形式（注：代码中的$\hat{a}$对应问题中的a，是二维矩阵的形式）。

附带代码（原代码） #

def get_I_matrix(N):
    return np.diag(np.ones(N))
  
def get_P_matrix(picks, N):
    return get_I_matrix(N)[np.ix_(picks),:]
  
def get_H_hat_diag(N, beta, picks):
    H_hat_diag = beta*np.ones(N*N)
    P_diag = np.zeros(N)
    P_diag[picks] += 1
    H_hat_diag = H_hat_diag + np.kron(P_diag, P_diag)
    return H_hat_diag
  
def get_a_hat_vec(a_hat, z, beta):
    return (a_hat + ifftw2d(beta*z)).flatten('F') # column first

a_hat = np.zeros([N, N])
a_hat = a_hat.astype(complex)
a_hat[np.ix_(picks, picks)] = B
    
H_hat_diag = get_H_hat_diag(N, beta, picks)
a_hat_vec = get_a_hat_vec(y_hat, z, beta)

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