Skip to content

Instantly share code, notes, and snippets.

@perimosocordiae
Created September 25, 2015 18:03
Show Gist options
  • Save perimosocordiae/efabc30c4b2c9afd8a83 to your computer and use it in GitHub Desktop.
Save perimosocordiae/efabc30c4b2c9afd8a83 to your computer and use it in GitHub Desktop.
Asymmetric Least Squares
import numpy as np
from scipy.linalg import solveh_banded
def als_baseline(intensities, asymmetry_param=0.05, smoothness_param=1e6,
max_iters=10, conv_thresh=1e-5, verbose=False):
'''Computes the asymmetric least squares baseline.
* http://www.science.uva.nl/~hboelens/publications/draftpub/Eilers_2005.pdf
smoothness_param: Relative importance of smoothness of the predicted response.
asymmetry_param (p): if y > z, w = p, otherwise w = 1-p.
Setting p=1 is effectively a hinge loss.
'''
smoother = WhittakerSmoother(intensities, smoothness_param, deriv_order=2)
# Rename p for concision.
p = asymmetry_param
# Initialize weights.
w = np.ones(intensities.shape[0])
for i in xrange(max_iters):
z = smoother.smooth(w)
mask = intensities > z
new_w = p*mask + (1-p)*(~mask)
conv = np.linalg.norm(new_w - w)
if verbose:
print i+1, conv
if conv < conv_thresh:
break
w = new_w
else:
print 'ALS did not converge in %d iterations' % max_iters
return z
class WhittakerSmoother(object):
def __init__(self, signal, smoothness_param, deriv_order=1):
self.y = signal
assert deriv_order > 0, 'deriv_order must be an int > 0'
# Compute the fixed derivative of identity (D).
d = np.zeros(deriv_order*2 + 1, dtype=int)
d[deriv_order] = 1
d = np.diff(d, n=deriv_order)
n = self.y.shape[0]
k = len(d)
s = float(smoothness_param)
# Here be dragons: essentially we're faking a big banded matrix D,
# doing s * D.T.dot(D) with it, then taking the upper triangular bands.
diag_sums = np.vstack([
np.pad(s*np.cumsum(d[-i:]*d[:i]), ((k-i,0),), 'constant')
for i in xrange(1, k+1)])
upper_bands = np.tile(diag_sums[:,-1:], n)
upper_bands[:,:k] = diag_sums
for i,ds in enumerate(diag_sums):
upper_bands[i,-i-1:] = ds[::-1][:i+1]
self.upper_bands = upper_bands
def smooth(self, w):
foo = self.upper_bands.copy()
foo[-1] += w # last row is the diagonal
return solveh_banded(foo, w * self.y, overwrite_ab=True, overwrite_b=True)
@zhenyou1
Copy link

zhenyou1 commented Apr 4, 2019

It would be great to give a data of intensities. People then can run it and gain a first experience.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment