[Math] Connection between SVD and Discrete Fourier Transform for Denoising

fourier analysisfourier transformmatricesmatrix decompositionsignal processing

Denoising signals (in particular, 2D arrays, such as images) can be done by removing the high frequency components of the discrete Fourier transform (which is related to convolution with a Gaussian kernel) or by removing the smallest singular values.
I was wondering if there is a known, specific mathematical connection between these two approaches.

I've seen a little discussion on the topic here and here, but I didn't really glean much specifically except the mention of circulant matrices.

Best Answer

There is an entire textbook written about this "Introduction to Inverse Problems and Imaging" by Bertaro and Boccacci. Chapter 4 sets up the formalism which intertwines inverse problems and fourier transforms. Chapter 5 shows how tikhonov regularization is related to windowing functions in fourier transforms. Here, I'll flush out one connection, that may appease your curiosity.

------------- THE CONNECTION -----------------------------------

Specifically, both the SVD and the fourier transform are characterized by some unitary transform (i.e. a change of basis) (1). Once in some new basis, you can make a choice to remove parts of the basis which are sensitive to small changes in the data (2). For both SVD and DFT these components are removed by some "windowing function" which select a viable window of either singular values or frequencies.

Assertion (1): For the SVD a matrix $A=U \Sigma V^T$ where $U$ and $V$ are unitary transforms into either the row or column space of $A$. The Fourier transform is also a unitary transform into a new basis "from the point of view of the linear canonical transformation, the Fourier transform is rotation by 90° in the time–frequency domain, and preserves the symplectic form." (https://en.wikipedia.org/wiki/Fourier_transform#Uncertainty_principle)

Assertion (2): For the SVD "parts of the basis sensitive to small changes" means the singular vectors (basis vectors) which correspond to small singular values. For Fourier Transform, I mean the Hilbert space basis vectors $\exp[i ( t \cdot f)]$ with large frequency $f$. The small singular values raise the condition number of the matrix and increase the upper bound on the signal to noise ratio and basis vectors with large frequency do the same.

------------- THE MATH -----------------------------------

Consider the a integral transform with kernel $K$ $$g(t) = \int dt' K(t,t') h(t') $$

%%%%% SINGULAR VALUE DECOMPOSITION %%%%%

Discretize the integral domain into $N \times N$ intervals. By the SVD, $K = U \Sigma V^T$ we can write the ordinary least squares (OLS) solution, $$\vec{h} = V \Sigma^{-1} U^T \vec{g} $$ or equivalently, $$\vec{h} = \sum_{i=1}^N \frac{1}{\sigma_i} \vec{v}_i \vec{u}_i^T \vec{g} $$ If we cutout $\sigma_i$'s based on some tolerance, $\epsilon$ then we have a windowing function, $$ W^{cut}_i(\sigma) = \frac{1}{\sigma_i} \Theta(\sigma_i - \epsilon)$$ Where $\Theta$ is the step function so that if $\sigma_i > \epsilon$ then $W_i=\sigma_i^{-1}$ and 0 otherwise. We can write the cutoff solution as $$\vec{h}^{cut} = \sum_{i=1}^N W^{cut}_i(\sigma_i) \vec{v}_i \vec{u}_i^T \vec{g} $$ Many different windowing functions can be concocted, for example: $$ W^{Tik}_i(\sigma_i) = \frac{\sigma_i}{\sigma_i^2+\mu}$$ This is the Tikhonov regularization solution (take the lime $\mu \rightarrow 0$, you'll recover the original $\frac{1}{\sigma_i}$ from the OLS solution ).

%%%%%NOTE ON LINEAR EXPANSIONS%%%%%

I have assumed the definition of a linear expansion and it's important to understanding the answer so I'll flush it out.

Consider, $\vec{h} = \sum_{i=1}^N \frac{1}{\sigma_i} \vec{v}_i \vec{u}_i^T \vec{g} $. In accordance with the singular value decomposition (wiki it), the $\vec{u}$ $\vec{v}$ vectors serve as basis vectors of the four fundamental sub-spaces of linear algebra.

Note $\vec{u}_i^T \vec{g} $ and $\sigma_i$ are both scalars. Thus these terms only serve as coefficients of a linear expansion of the solution in the column space of K. This looks like $$\vec{h} = \sum_{i=1}^N a_i \vec{v}_i$$ here I have grouped $a_i = \frac{\vec{u}_i^T \vec{g}}{\sigma_i}$ so the windowing function is reweighting the basis components, $\vec{v}_i$ in the expansion.


%%%% FOURIER TRANSFORM %%%%

If we fourier transform this problem into frequency space the convolution turns into a multiplication (See E.O. Brigham's ``Introduction to FFT'').

Thus, $g(f) = K(f,f') h(f') $ and consequently $$h(f') = g(f)/K(f,f')$$ We can fourier transform over frequency space $\mathcal{F}$ to arrive at a solution. $$h(t') = \int_{\mathcal{F}} df \frac{g(f')}{K(f,f')}\exp[i ( t \cdot f)]$$ or equivalently, $$h(t') = \int_{\mathcal{F}} df \frac{\exp[i ( t \cdot f)]}{K(f,f')} g(f')$$ If you want to introduce some cutoff on high frequency data then you can introduce a windowing function as before: $$ W^{cut}(f) = \frac{\exp[i ( t \cdot f)]}{K(f,f')} \Theta(f - \epsilon)$$ So that $$h(t') = \int_{\mathcal{F}} df W^{cut}(f) g(f')$$

You could introduce the same Tikhonov regularization to the frequency basis vectors by means of $$ W^{Tik}(f) = \exp[i ( t \cdot f)] \frac{K^*(f,f')}{\vert K(f,f') \vert^2 + \mu}$$ Or you could select some other windowing function more broadly used by the DFT community (e.g. Hamming function)


Related Question