# Asymmetric least squares smoothing

// by tony.withers@uwo.ca, using method of Eilers, PHC and Boelens, HFM
// (2005) Baseline correction with asymmetric least squares smoothing.

// Creates (and overwrites) w_base, a baseline estimate for w_data. The
// asymmetry parameter (Eilers and Boelens' p) generally takes values
// between 0.001 and 0.1. Try varying lambda in orders of magnitude
// between 10^2 and 10^9. Not efficient for large N, try it for w_data
// with fewer than 1000 points.
function ALS(w_data, lambda, asymmetry)
wave w_data
variable lambda, asymmetry

variable i, N=numpnts(w_data), rms=inf
variable maxIts=20

matrixOp /free  D = identity(N)
differentiate /EP=1/METH=2/DIM=0 D
differentiate /EP=1/METH=2/DIM=0 D

// this step (specifically the matrix multiplication) is slow:
matrixOp /free H = lambda * (D^t x D)

duplicate /o/free w_data w, w_new
w=1

for (i=0;i<maxIts;i+=1)
matrixOp /o/free  C = chol(diagRC(w, N, N)+H)
matrixOp /o w_base = backwardSub(C,(forwardSub(C^t, w * w_data)))
w_new = asymmetry * (w_data>w_base) + (1-asymmetry) * (w_data<w_base)

// convergence test
w-=w_new
wavestats /Q w
if (v_rms>=rms)
return i+1
else
rms=v_rms
w=w_new
endif
endfor
return 0
end

I changed  the original code from tony a little bit to be able to use it on larger datasets:

// Originally developed by tony withers:
//
// by tony.withers@uwo.ca, using method of Eilers, PHC and Boelens, HFM
// (2005) Baseline correction with asymmetric least squares smoothing.

// Creates (and overwrites) w_base, a baseline estimate for w_data. The
// asymmetry parameter (Eilers and Boelens' p) generally takes values
// between 0.001 and 0.1. Try varying lambda in orders of magnitude
// between 10^2 and 10^9. Not efficient for large N, try it for w_data
// with fewer than 1000 points.
//
//
// I just changed the code to avoid the slow matrix multiplication.
// The H-matrix is now constructed "manually". This saves time and memory
// allows larger datasets.
// (kmichel@wzw.tum.de)

function ALS(w_data, lambda, asymmetry)
wave w_data
variable lambda, asymmetry
variable i, N=numpnts(w_data), rms=inf
variable maxIts=20

//    matrixOp /free  D = identity(N)
//    differentiate /EP=1/METH=2/DIM=0 D
//    differentiate /EP=1/METH=2/DIM=0 D

//      this step (specifically the matrix multiplication) is slow:
//      matrixOp /o H = lambda * (D^t x D)

Make /O/N=(N) diag0 = 6
wave diag0
diag0=1
diag0=5
diag0[N-2]=5
diag0[N-1]=1

Make /O/N=(N-1) diag1 = -4
wave diag1
diag1=-2
diag1[N-2]=-2

Make /O/N=(N-2) diag2 = 1
Make /O/N=(N,N) H
wave  H
matrixoP/o H = setoffdiag(H,0,diag0)
matrixop/o H = setoffdiag(H,-1,diag1)
matrixop/o H = setoffdiag(H,1,diag1)
matrixop/o H = setoffdiag(H,-2,diag2)
matrixop/o H = setoffdiag(H,2,diag2)
matrixop/o H = lambda * H
killwaves diag0, diag1, diag2
duplicate /o/free w_data w, w_new
w=1

for (i=0;i<maxIts;i+=1)
matrixOp /o /free  C = chol(diagRC(w, N, N)+H)
matrixOp /o w_base = backwardSub(C,(forwardSub(C^t, w * w_data)))
w_new = asymmetry * (w_data>w_base) + (1-asymmetry) * (w_data<w_base)

// convergence test
w-=w_new
wavestats /Q w
if (v_rms>=rms)
killwaves H
return i+1
else
rms=v_rms
w=w_new
endif

endfor

return 0
end Forum Support Gallery