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;

A better implementation can be download at https://www.mathworks.com/matlabcentral/fileexchange/25634-smoothn

%--- 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