Fitting a piecewise constant function to noisy data

f[x_]:=Piecewise[{{3/4,0≤x≤1.5},{2,1.5<x≤3},{4/3,3<x≤5}}]
original=Plot[f[x],{x,0,5},PlotRange->{{0,5.1},{0,2.5}}]

dx=1/10;
Npoints=5/dx+1;

y=Table[f[k dx],{k,0,Npoints-1}]+RandomVariate[NormalDistribution[0,0.1],Npoints]
pts=Table[{k dx,y[[k+1]]},{k,0,Npoints-1}]
Show[original,ListPlot[pts]]

Original function and noisy data

DM=SparseArray[{Band[{1,1}]->-1,Band[{1,2}]->1},{Npoints-1,Npoints}];DM//MatrixForm

xhat=Array[x,Npoints]

DM.xhat

Fit using 2-norm penalty term (weak and strong strengths)

gamma=1;
soln=Minimize[Norm[y-xhat,2] +gamma Norm[DM.xhat,2],xhat];
Show[original,ListPlot[Table[{k dx,x[k+1]}/.soln[[2]],{k,0,Npoints-1}]]]

gamma=2;
soln=Minimize[Norm[y-xhat,2] +gamma Norm[DM.xhat,2],xhat];
Show[original,ListPlot[Table[{k dx,x[k+1]}/.soln[[2]],{k,0,Npoints-1}]]]

Fit using 1-norm penalty term of different strengths

gamma=.1;
soln=Minimize[Norm[y-xhat,2] +gamma Norm[DM.xhat,1],xhat];
Show[original,ListPlot[Table[{k dx,x[k+1]}/.soln[[2]],{k,0,Npoints-1}]]]

gamma=.4;
soln=Minimize[Norm[y-xhat,2] +gamma Norm[DM.xhat,1],xhat];
Show[original,ListPlot[Table[{k dx,x[k+1]}/.soln[[2]],{k,0,Npoints-1}]]]

gamma=2;
soln=Minimize[Norm[y-xhat,2] +gamma Norm[DM.xhat,1],xhat];
Show[original,ListPlot[Table[{k dx,x[k+1]}/.soln[[2]],{k,0,Npoints-1}]]]

Fit using 1-norm inequality constraint (very slow)

gamma=.1;
soln=Minimize[{Norm[y-xhat,2] , Norm[DM.xhat,1]≤2},xhat];
Show[original,ListPlot[Table[{k dx,x[k+1]}/.soln[[2]],{k,0,Npoints-1}]]]