If you are still interested in smoothing with L0-penalties I would give a look to the following reference: "Visualization of Genomic Changes by Segmented Smoothing Using an L0 Penalty" - DOI: 10.1371/journal.pone.0038230 (a nice intro to the Whittaker smoother can be found in P. Eilers paper "A perfect smoother" - DOI: 10.1021/ac034173t). Of course, in order to achieve your objective you have to work a bit around the method.
In principle, you need 3 ingredients:
- The smoother - I would use the Whittaker smoother. Also, I will use matrix augmentation (see Eilers and Marx, 1996 - "Flexible Smoothing with B-splines
and Penalties", p.101).
- Quantile regression - I will use the R package quantreg (rho = 0.5) for laziness :-)
- L0-penalty - I will follow the mentioned "Visualization of Genomic Changes by Segmented Smoothing Using an L0 Penalty" - DOI: 10.1371/journal.pone.0038230
Of course, you would need also a way to select the optimal amount of smoothing. This is done by my carpenter eyes for this example. You could use the criteria in DOI: 10.1371/journal.pone.0038230 (pg. 5, but I did not try it on your example).
You will find a small code below. I left some comments as guide through it.
# Cross Validated example
rm(list = ls()); graphics.off(); cat("\014")
library(splines)
library(Matrix)
library(quantreg)
# The data
set.seed(20181118)
n = 400
x = 1:n
true_fct = stepfun(c(100, 200, 250), c(200, 250, 300, 250))
y = true_fct(x) + rt(length(x), df = 1)
# Prepare bases - Identity matrix (Whittaker)
# Can be changed for B-splines
B = diag(1, n, n)
# Prepare penalty - lambda parameter fix
nb = ncol(B)
D = diff(diag(1, nb, nb), diff = 1)
lambda = 1e2
# Solve standard Whittaker - for initial values
a = solve(t(B) %*% B + crossprod(D), t(B) %*% y, tol = 1e-50)
# est. loop with L0-Diff penalty as in DOI: 10.1371/journal.pone.0038230
p = 1e-6
nit = 100
beta = 1e-5
for (it in 1:nit) {
ao = a
# Penalty weights
w = (c(D %*% a) ^ 2 + beta ^ 2) ^ ((p - 2)/2)
W = diag(c(w))
# Matrix augmentation
cD = lambda * sqrt(W) %*% D
Bp = rbind(B, cD)
yp = c(y, 1:nrow(cD)*0)
# Update coefficients - rq.fit from quantreg
a = rq.fit(Bp, yp, tau = 0.5)$coef
# Check convergence and update
da = max(abs((a - ao)/ao))
cat(it, da, '\n')
if (da < 1e-6) break
}
# Fit
v = B %*% a
# Show results
plot(x, y, pch = 16, cex = 0.5)
lines(x, y, col = 8, lwd = 0.5)
lines(x, v, col = 'blue', lwd = 2)
lines(x, true_fct(x), col = 'red', lty = 2, lwd = 2)
legend("topright", legend = c("True Signal", "Smoothed signal"),
col = c("red", "blue"), lty = c(2, 1))
PS. This is my first answer on Cross Validated. I hope it is useful and clear enough :-)