基于Tikhonov正则化的图像去噪Python代码维度不匹配问题排查

问题描述 #

采用梯度下降法最小化Tikhonov正则化实现图像去噪的Python代码运行时持续出现维度不匹配错误，但相同逻辑的Matlab代码可正常执行，代码及报错信息如下。

报错信息 #

文件路径：C:\Users\angel\anaconda3\lib\site-packages\scipy\sparse_base.py，第559行，_mul_dispatch函数中
抛出错误：ValueError: 维度不匹配

错误原因分析 #

核心问题出在差分算子的维度设计与梯度计算逻辑：

1. 差分算子维度错误：原DiffOper(N)生成的是(N, N+1)的一维差分矩阵，但图像是N×N维度，展开为长度为N²的向量，二者维度完全不兼容，导致D @ x运算时出现维度不匹配。
2. 梯度公式错误：原代码中的梯度推导不符合Tikhonov正则化的损失函数定义，进一步加剧了运算维度冲突。

Matlab代码能正常运行，是因为其差分算子是针对二维图像结构构建的，与展开后的图像向量维度匹配，而Python代码仅实现了一维差分，未适配二维图像的向量展开形式。

修复方案 #

重构差分算子为适配二维图像的联合稀疏矩阵，并修正梯度计算逻辑，修改后的完整代码如下：

import numpy as np
import scipy.sparse as sp
import matplotlib.pyplot as plt

def psnr_fun(image1, image2):
    mse = np.mean((image1 - image2) ** 2)
    if mse == 0:
        return float('inf')
    max_pixel = np.max(image1)
    psnr = 20 * np.log10(max_pixel / np.sqrt(mse))
    return psnr

def DiffOper(N):
    # 构建水平差分算子（每行相邻像素差分）
    D_h = sp.diags([-np.ones(N), np.ones(N)], [0, 1], shape=(N, N+1), format='csc')
    D_h = D_h[:, :-1]  # 截断为(N, N)矩阵
    # 扩展为N×N图像的块对角矩阵，维度(N(N-1), N²)
    D_h_full = sp.block_diag([D_h]*N, format='csc')
    
    # 构建垂直差分算子（每列相邻像素差分）
    D_v = sp.diags([-np.ones(N*(N-1)), np.ones(N*(N-1))], [0, N], shape=(N*(N-1), N*N), format='csc')
    
    # 联合水平+垂直差分，总算子维度(2N(N-1), N²)
    D = sp.vstack([D_h_full, D_v], format='csc')
    return D

def tikhonov_gd(y, mu, lamda, nit, tol, n):
    x = y.copy()
    N = int(np.sqrt(n))
    D = DiffOper(N)
    RelErr = np.zeros(nit)

    for k in range(nit):
        x_old = x.copy()
        # 计算Tikhonov损失的梯度：∇L(x) = x - y + λD^T Dx
        Dx = D @ x
        grad = (x - y) + lamda * D.T @ Dx
        # 梯度下降更新：x = x - μ*grad
        x = x - mu * grad

        RelErr[k] = np.linalg.norm(x - x_old, 2) / np.linalg.norm(x, 2)

        if RelErr[k] < tol:
            RelErr = RelErr[:k+1]
            break

    return x, RelErr

# 主代码
N = 256
n = N * N
x = np.random.random(n)  # 示例原图
sigma = 0.09
y = x + sigma * np.max(x) * np.random.randn(n)  # 生成带噪图像

mu = 0.001  # 调整步长避免发散
nit = 1000
tol = 1e-10
lamda = 1
v, Err = tikhonov_gd(y, mu, lamda, nit, tol, n)

# 绘图展示
plt.figure(figsize=(10, 6))

plt.subplot(221)
plt.imshow(x.reshape(N, N), cmap='gray')
plt.axis('image')
plt.title('Original')

plt.subplot(222)
plt.imshow(y.reshape(N, N), cmap='gray')
plt.axis('image')
plt.title('Noisy')

plt.subplot(223)
plt.imshow(v.reshape(N, N), cmap='gray')
plt.axis('image')
psnr = psnr_fun(v, x)
plt.title(f'Tikhonov (PSNR = {psnr:.3f} dB)')

plt.subplot(224)
plt.semilogy(Err, linewidth=2.5, color='black')
plt.xlabel('Iterations (k)', fontsize=12)
plt.ylabel('Relative Error', fontsize=12)
plt.axis('tight')
plt.grid()
plt.legend(['Gradient Descent'], loc='best', fontsize=12)

plt.tight_layout()
plt.show()

关键修改说明 #

• 差分算子适配：重构后的DiffOper生成联合水平、垂直差分的稀疏矩阵，维度与展开后的图像向量完全匹配。
• 梯度公式修正：按照Tikhonov正则化损失函数的正确梯度推导更新迭代公式，保证逻辑正确性。
• 步长调整：将原步长0.05调整为0.001，避免迭代过程发散。

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