Asymmetric least squares smoothing

Average rating
(0 votes)

// by <a href="mailto:tony.withers@uwo.ca" rel="nofollow">tony.withers@uwo.ca</a>, 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

Back to top