# Advantage of diagonal “jitter” for numerical stability?

In a machine learning code, that computes optimum parameters $$\theta _{MLE}$$ of a linear regression model, by maximum likelihood estimation:

$$\boldsymbol \theta^\text{ML} = (\boldsymbol\Phi^T\boldsymbol\Phi )^{-1}\boldsymbol\Phi^T\boldsymbol y$$

Where $$y$$ is the target vector and $$\Phi$$ is the polynomial feature matrix. In the linked notebook we can find:

For reasons of numerical stability, we often add a small diagonal "jitter" $$\kappa$$ to $$\boldsymbol\Phi^T\boldsymbol\Phi$$ so that we can invert the matrix without significant problems so that the maximum likelihood estimate becomes $$\boldsymbol \theta^\text{ML} = (\boldsymbol\Phi^T\boldsymbol\Phi + \kappa\boldsymbol I)^{-1}\boldsymbol\Phi^T\boldsymbol y$$

In the code, $$\kappa$$ is very small value of 1e-08.

So, how does the diagonal "jitter" $$\kappa$$ affects stability?

So, you want to invert your matrix $$A=\Phi^T\Phi$$. For $$A$$ to be invertible it must not have zero eigenvalues. We can show that $$A$$ is positive semi-definite as follows. Positive semi-definite means that the eigenvalues of $$A$$ are $$\geq 0$$. This is equivalent to showing $$y^TAy \geq 0, \forall y \neq 0$$.

$$y^TAy = y^T\Phi^T\Phi{y}=(\Phi{y})^T(\Phi{y}) \geq 0$$

So we have proved that $$A$$ is positive semi-definite. Therefore, it could have eigenvalues which are zero and will therefore render it non-invertible. So, we replace $$A$$ with $$A+\kappa{I}$$, where $$\kappa > 0$$ can be chosen to render $$A+\kappa{I}$$ positive definite and therefore invertible.

$$y^T(A+\kappa{I})y = y^T{A}y + \kappa{y}^Ty$$

Since $$y^T{A}y \geq 0$$ and $$y^Ty>0$$, choosing a small positive $$\kappa$$ renders $$y^T(A+\kappa{I})y > 0 \,\,\, \forall y\neq0$$ and therefore invertible.

Look up something on Tikhonov regularization, also known as ridge regression in machine learning. This is a standard technique (but I agree that the explanation in that notebook is somewhat poor).

Technically speaking, it does not affect the numerical stability of that algorithm, but it modifies the problem to a more well-conditioned one, from $$\min \|\Phi \theta - y\|^2$$ to $$\min \|\Phi \theta - y\|^2 + \kappa \|\theta\|^2.$$

• oops thanks for the fix @cfdlab. I am not used to this notation, where I work $y$ is a perfectly good name for an unknown and $\theta$ for a known parameter. :p – Federico Poloni Nov 22 '20 at 15:25
• Hello Federico!, I don't know if Tikhonov regularization is the same as $l_2$ regularization or not, but if it is, there is already one deliverable in the notebook to fit a regularized linear regression model (through maximizing a posteriori method), and as you can see $\kappa$ is set to be very much close to zero: 1e-08, so it's useless as a regularization coefficient. I believe it's intended for some inverse calculation stability for numpy, or may be something close to ConvexHull's answer. – Algo Nov 22 '20 at 16:31
• @Algo I think there are several ways how to intepret such modifications. Federico added a reference on a specific method which is frequently used in machine learning. Here you can justify it from a Bayesian viewpoint. – ConvexHull Nov 22 '20 at 16:56
• @Algo yes, it is the same as l_2 regularization. I don't know what it is "intended" for, but in practice that is Tikhonov regularization. It modifies the problem: what you return is not a solution to the original problem, but to the regularized one. You may think it is for stability, but that does not change what it does. (Since we are speaking about stability, by the way, that method to solve least-squares problem is not backward stable, and you should probably think about switching to SVD.) – Federico Poloni Nov 22 '20 at 18:32
• Also, do not think that $\kappa$ is useless as a regularization coefficient: the cases in which it matters are precisely those in which there are singular values of the order of $10^{-8}$, and hence that little modification changes the exact solution significantly. (I assume that $\|\Psi\|$ is significantly larger than $10^{-8}$, here; proper practice would suggest using a relative threshold like $\kappa \|\Psi\|$.) – Federico Poloni Nov 22 '20 at 18:39

Think of the simplest case when $$\Phi$$ is a scalar value.

Not well defined:

$$\boldsymbol \theta^\text{ML} = (0^T 0)^{-1}0^T ~ y = \frac{1}{0} 0~y= \frac{0}{0}$$

Well defined:

$$\boldsymbol \theta^\text{ML} = (0^T 0 + \kappa)^{-1}0^T~y =\frac{1}{\kappa} 0 ~y= 0$$