如何复原被两个循环矩阵模糊的图像矩阵？

Alright, let's break down how to recover your original image $X$ from the blurred observation $Y = \Phi_1 X \Phi_2$. The key here is leaning into the diagonalization property of cyclic matrices in the Fourier domain—this turns a messy matrix multiplication problem into simple element-wise operations, which is a game-changer.

1. First, Map the Blur to the Fourier Domain #

Cyclic matrices commute with the Discrete Fourier Transform (DFT) matrix. That means we can rewrite $\Phi_1$ and $\Phi_2$ as:
$$\Phi_1 = F_1^\ast \Lambda_1 F_1, \quad \Phi_2 = F_2^\ast \Lambda_2 F_2$$
where $F_1$, $F_2$ are the 1D DFT matrices for sizes $n_1$ and $n_2$, and $\Lambda_1$, $\Lambda_2$ are diagonal matrices. The diagonal entries of $\Lambda_1$ are just the DFT of the first row of $\Phi_1$, and same for $\Lambda_2$ and $\Phi_2$'s first row.

Applying the 2D DFT to both sides of your blur equation gives us a much simpler relationship:
$$\hat{Y} = \Lambda_1 \hat{X} \Lambda_2$$
Here, $\hat{Y}$ is the 2D FFT of $Y$, $\hat{X}$ is the 2D FFT of your original $X$, and the right-hand side is just element-wise multiplication (since diagonal matrices multiply element-wise when acting on a matrix).

2. Basic Restoration: Inverse Filtering (Noiseless Case) #

If your observation $Y$ is completely noise-free (a rare scenario, but let's cover it first), you can directly solve for $\hat{X}$:
$$\hat{X}[u,v] = \frac{\hat{Y}[u,v]}{\Lambda_1[u,u] \cdot \Lambda_2[v,v]}$$
Where $\hat{X}[u,v]$ is the $(u,v)$-th element of the 2D FFT result, and $\Lambda_1[u,u]$, $\Lambda_2[v,v]$ are the $u$-th and $v$-th diagonal entries of $\Lambda_1$ and $\Lambda_2$.

Then just take the inverse 2D FFT of $\hat{X}$ to get your restored $X$. Warning: If $\Lambda_1[u,u] \cdot \Lambda_2[v,v]$ is close to zero, this will blow up numerically—so you need regularization for practical cases.

3. Practical Restoration: Regularized Methods (For Noise/Ill-Conditioned Cases) #

Even if you don't explicitly have noise, the blur itself can make the inverse problem ill-posed. Given that your cyclic matrices have Gaussian-distributed non-zero elements, statistical regularization makes sense here. Two go-to methods:

3.1 Wiener Filtering #

Wiener filtering minimizes the mean squared error between your restored image and the true $X$. In the Fourier domain, the solution is:
$$\hat{X}[u,v] = \frac{\overline{\Lambda_1[u,u]} \cdot \overline{\Lambda_2[v,v]}}{|\Lambda_1[u,u]|^2 |\Lambda_2[v,v]|^2 + \frac{\sigma_n^2}{\sigma_x2}} \hat{Y}[u,v]$$

• $\sigma_n^2$: Variance of any observation noise (if you don't have explicit noise, you can estimate this using the given $\sigma$ of your cyclic matrix elements, since the blur acts like a form of structured noise)
• $\sigma_x^2$: Variance of your original image's pixels—you can estimate this from local patches of $Y$ (assuming your image is stationary) or use a prior if you know something about the image.

3.2 Tikhonov Regularization (Regularized Least Squares) #

If you prefer a deterministic approach, set up a least squares problem with a regularization term to penalize overly "rough" restored images:
$$\min_X | Y - \Phi_1 X \Phi_2 |_F^2 + \lambda | X |_F^2$$
Where $\lambda$ is a regularization parameter, and $|\cdot|_F$ is the Frobenius norm. Again, using Fourier domain simplification, the solution is:
$$\hat{X}[u,v] = \frac{\overline{\Lambda_1[u,u]} \cdot \overline{\Lambda_2[v,v]}}{|\Lambda_1[u,u]|^2 |\Lambda_2[v,v]|^2 + \lambda} \hat{Y}[u,v]$$

To pick $\lambda$, you can use cross-validation (like leave-one-out) or derive it using the known $\sigma$: since each cyclic matrix has $2K+1$ non-zero elements per row, each distributed as $N(0,\sigma^2)$, you can compute the Frobenius norm of $\Phi_1$ and $\Phi_2$ to set a reasonable $\lambda$ that balances fitting the observation and keeping the solution stable.

4. Optimize Using the Bandwidth Constraint #

Since $\Phi_1$ and $\Phi_2$ have bandwidth $2K$, their first rows only have $2K+1$ non-zero elements. You don't need to compute the full DFT of the entire row—just compute the DFT of those non-zero positions, which saves computation time, especially for large images.

For implementation, use fast FFT libraries (like numpy.fft in Python or fft2 in MATLAB) to handle all the Fourier transforms efficiently. After computing $\hat{X}$, take the inverse FFT and discard any tiny imaginary parts (they're just numerical artifacts).

Quick Step-by-Step Recap #

1. Compute the 1D FFT of the first row of $\Phi_1$ and $\Phi_2$ to get the diagonal entries of $\Lambda_1$ and $\Lambda_2$.
2. Compute the 2D FFT of your blurred image $Y$ to get $\hat{Y}$.
3. Choose a restoration method (inverse filtering for noiseless, Wiener/Tikhonov for practical cases) to compute $\hat{X}$ element-wise.
4. Compute the inverse 2D FFT of $\hat{X}$, take the real part, and that's your restored image $X$.

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