# Find effective inverse of Toeplitz matrix

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).

$$\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).
• 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 – Aleksejs Fomins Sep 18 '19 at 12:15