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