Suppose a (noisy) series y, of length m is given. The observations have been sampled at equal distances.
We want to fit a smooth series z to y. We then have to balance two conflicting goals: (1) fidelity to the data and (2) roughness of z. The smoother z is, the more it will deviate from y. Roughness of z can be expressed with differences.

To avoid a lot of tedious algebra, it is advantageous to introduce matrices and vectors
$$Q=|\mathbf{y}-\mathbf{z}|^{2}+\lambda|\mathbf{D z}|^{2}$$ (1)
where the first term $$|\mathbf{y}-\mathbf{z}|^{2}$$ indicates the quadratic norm of error, and $\mathbf{D}$ is a matrix such that $\mathbf{D} \mathbf{z}=\Delta \mathbf{z}$.

To solve the z, such that the $$Q=|\mathbf{y}-\mathbf{z}|^{2}+\lambda|\mathbf{D z}|^{2}$$ can be minimized, we should introduce the partial derivatives( the minimal or maximal is always at the point when the partial derivatives be zero).
$$\partial Q / \partial \mathbf{z}^{\prime}=-2(\mathbf{y}-\mathbf{z})+2 \lambda \mathbf{D}^{\prime} \mathbf{D} \mathbf{z}$$
and equating this to 0 we get the linear system of equations
$$\left(\mathbf{I}+\lambda \mathbf{D}^{\prime} \mathbf{D}\right) \mathbf{z}=\mathbf{y}$$ where I is the identity matrix.

The matrix $\mathbf{D}$ is actually an differential matrix. To make more sense, In the avbove, For example,when m=5, $\mathbf{D}$ would be
$$\mathbf{D}=\left[ \begin{array}{rrrrr}{-1} & {1} & {0} & {0} & {0} ;\ {0} & {-1} & {1} & {0} & {0} ;\ {0} & {0} & {-1} & {1} & {0}; \ {0} & {0} & {0} & {-1} & {1}\end{array}\right]$$

Implementation in Matlab is very straightforward, because we can use the built-in function diff() to compute $D$:

m = length(y);
E = eye(m);
D = diff(E);
z =  (E + lambda * D' * D)\y;


 %--- Example #1: smooth a curve ---
x = linspace(0,100,2^8);
y = cos(x/10)+(x/50).^2 + randn(size(x))/10;
y([70 75 80]) = [5.5 5 6];
z = smoothn(y); % Regular smoothing
zr = smoothn(y,'robust'); % Robust smoothing
subplot(121), plot(x,y,'r.',x,z,'k','LineWidth',2)
axis square, title('Regular smoothing')
subplot(122), plot(x,y,'r.',x,zr,'k','LineWidth',2)
axis square, title('Robust smoothing') Reach matcodelab@gmail.com