1
$\begingroup$

I would like to do a deconvolution of a noisy process. $$y_i = \sum_j k_{j-i} x_{j} + \nu$$ where $k$ is some well-behaved localized kernel (e.g. gaussian), and $\nu$ is gaussian noise with zero mean and known variance. This equation can be rewritten in a vector form as $$\vec{y} = K\vec{x} + \vec{\nu}$$ where $K$ is a Toeplitz matrix with the kernel propagating along the columns. Now, as far as I can tell, the MLE solution to this equation is the Moore-Penrose pseudoinverse

$$\vec{x} = M\vec{y} = (K^TK)^{-1}K^T\vec{y}$$

Unfortunately, $K$ happens to be singular for many interesting cases, so naive inversion is not possible. I have learned by googling that, among others, the Levinson recursion is used to find $\vec{x}$. Except I don't want $\vec{x}$. I want to get $M$, and explore how it looks like for different kernels.

I want a practical suggestion on how to compute $M$. I am ok with having to sacrifice a few values of $\vec{x}$ at the boundaries of the domain in order to get a well behaved result (in that case $M$ would be rectangular).

$\endgroup$
1
$\begingroup$

A practical (but not optimal) solution would be to exploit Tikhonov Regularization, that is (in the simplest case) accomplished with the help of diagonal loading in the following form:

$$\vec{x} = (K^TK + \lambda I)^{-1}K^T\vec{y} \tag{1}$$

Other high pass operators maybe used as a generalization of the above case, using a Tikhonov matrix $ \Gamma$ $$\vec{x} = (K^TK + \Gamma^T\Gamma )^{-1}K^T\vec{y} \tag{2}$$ Note that for $\Gamma = \sqrt{\lambda}I$, then $(2)$ becomes $(1)$. The above happens to be the solution of the problem $$\min_x \Vert y - Kx \Vert_2^2 + \Vert \Gamma x \Vert_2^2$$ where the second term is the regularization term (often referred to as penalty, or even a constraint). If you want an even more general solution, you could look at the Generalized Tikhonov regularization, i.e. the solution to the problem $$\hat{x} = \min_x \Vert y - Kx \Vert_P^2 + \Vert x - x_0 \Vert_Q^2 = x_0 + (K^T PK + Q)^{-1} (K^T P(y - Kx_0))\tag{3}$$ where $x_0$ could be seen as some prior info on $x$ and $P,Q$ are some co variances reflecting possible interactions between variables. Note that when $x_0 = 0$ and $P=I$ and $Q = \Gamma^T \Gamma$, then $(3)$ becomes $(2)$.

Now, another question is "how to find $\lambda$ in $(1)$". Some use heuristics (adhocs). Others resort to bayesian methods (as $\lambda = \frac{1}{\sigma^2}$ could be interpreted as a precision factor).

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks for your reply. This is kind of cool, for large $\lambda$ the matrix $M$ is diagonally-dominant. As I decrease $\lambda$, the diagonal spreads to be some kind of $sinc$ function that starts oscillating more and more wildly. Now I need to think if it makes sense to me :D $\endgroup$ – Aleksejs Fomins Sep 18 '19 at 12:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.