Inverting a finite difference operator plus a diagonal matrix with Fourier transform

125 Views Asked by At

I've been playing around with a 1-D smoothing problem using Tikhonov regularization in the objective function $J(x)$:

$$ J(x) = ||x - b||^2 + \lambda ||{\bf D}x||^2 $$

Where $b$ is data and ${\bf D}$ is a finite difference linear operator. This problem can be solved by:

$$ x = (I + \lambda {\bf D}^T {\bf D})^{-1}b $$

Where $I$ is the identity matrix.

When $b$ gets large, this can be solved using a sparse version of the linear algebra library in scipy. Alternatively, I could also solve it using FFT by doing the inversion as element-wise division in Fourier-space.

$$ {\bf H} = 1 + \lambda \cdot \mathcal{F}\{{\bf D}\}^* \mathcal{F}\{{\bf D}\} $$

$$ x = \mathcal{F}^{-1}\{ \mathcal{F}\{ b\} / {\bf H} \} $$

Where $\mathcal{F}$ and $\mathcal{F}^{-1}$ are the Fourier transform and its inverse, respectively. $\mathcal{F}\{x\}^*$ denotes complex conjugation of the Fourier-transformed $x$.

I don't have a strong math background and picked up this Fourier implementation from various image processing scripts. I'd like to learn more about why is it possible to replace the identity matrix with ones in the Fourier domain, or put in another way: is it possible to implement the same solution with Fourier transform when $I = W$, where $W$ is a diagonal matrix with non-zero, positive elements?