如何在Python中由模糊函数生成矩阵以构建Ax=y逆问题？

你需要构建形如 Ax=y 的逆问题，其中A是高斯模糊算子的矩阵形式，作用于扁平化的图像数组x，y为输出数据。当前代码中A = gaussian_filter(x, sigma=1)得到的是模糊后的图像，而非算子矩阵，因此无法用A(x)调用。以下是解决方法：

核心思路 #

高斯模糊本质是局部加权求和，将该操作转化为矩阵形式时：

• 若图像尺寸为M×N，扁平化后是M*N×1的向量，A则是M*N×M*N的矩阵
• 矩阵A的每一行对应输出图像中一个像素的加权系数，仅与高斯核覆盖的输入像素有非零值，因此用稀疏矩阵存储可避免内存爆炸

完整修改代码 #

# import library
import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import gaussian_filter
from PIL import Image
from scipy.sparse import lil_matrix

# 生成高斯核
def create_gaussian_kernel(sigma, kernel_size=None):
    if kernel_size is None:
        kernel_size = int(6 * sigma + 1)
        if kernel_size % 2 == 0:
            kernel_size += 1
    ax = np.linspace(-(kernel_size//2), kernel_size//2, kernel_size)
    xx, yy = np.meshgrid(ax, ax)
    kernel = np.exp(-(xx**2 + yy**2)/(2*sigma**2))
    kernel /= np.sum(kernel)
    return kernel

# 构建模糊算子的稀疏矩阵形式
def build_blur_matrix(image_shape, kernel):
    rows, cols = image_shape
    total_pixels = rows * cols
    kernel_rows, kernel_cols = kernel.shape
    pad_rows = kernel_rows // 2
    pad_cols = kernel_cols // 2

    # 使用LIL格式构建稀疏矩阵，便于逐行赋值
    A = lil_matrix((total_pixels, total_pixels))
    
    for i in range(rows):
        for j in range(cols):
            idx = i * cols + j  # 当前输出像素的扁平化索引
            # 遍历高斯核的每个元素
            for ki in range(kernel_rows):
                for kj in range(kernel_cols):
                    ni = i + ki - pad_rows
                    nj = j + kj - pad_cols
                    # 仅处理图像范围内的输入像素
                    if 0 <= ni < rows and 0 <= nj < cols:
                        nidx = ni * cols + nj
                        A[idx, nidx] = kernel[ki, kj]
    return A.tocsr()  # 转换为CSR格式，提升矩阵乘法效率

# import image
img = Image.open("cameraman.jpg")
# img.show()

# 扁平化图像
x = np.squeeze(np.asarray(img))
image_shape = x.shape

# 构建矩阵形式的A
sigma = 1
kernel = create_gaussian_kernel(sigma)
A = build_blur_matrix(image_shape, kernel)

# 计算y = Ax：先扁平化x，矩阵乘法后恢复形状
x_flat = x.flatten()
y_flat = A @ x_flat
y = y_flat.reshape(image_shape)

# 验证：与直接调用gaussian_filter的结果对比
y_direct = gaussian_filter(x, sigma=sigma)
print("矩阵乘法与直接滤波的最大误差：", np.max(np.abs(y - y_direct)))

关键说明 #

1. 稀疏矩阵优势：对于512×512的图像，稠密矩阵需占用约524GB内存，而稀疏矩阵仅存储非零元素，内存占用大幅降低
2. 边界处理：代码中默认忽略核超出图像边界的部分，若需与gaussian_filter默认的0填充边界完全匹配，可先对图像做边界填充再构建矩阵
3. 正则化操作：稀疏矩阵的转置可直接通过A.T获取，便于后续正则化计算（如Tikhonov正则化中的A.T @ A）

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