Baseline Correction with AsLS and LAsLS

Baseline drift has been a persistent headache in spectroscopy for decades. In Raman spectroscopy, fluorescence from the sample (or from impurities) produces a broad, curved background that can be far larger than the Raman peaks themselves. In infrared and near-infrared spectroscopy, scattering and instrumental effects add sloping backgrounds that change from sample to sample. For a long time, the usual fix was to pick a few "baseline points" by hand and draw a polynomial through them, hoping it behaved in between. That works for one spectrum on a quiet afternoon. It does not work for a thousand spectra, or for backgrounds with an awkward shape.

In 2005, Paul H. C. Eilers (Leiden University Medical Centre) and Hans F. M. Boelens (University of Amsterdam) described a remarkably elegant alternative in a short report titled Baseline Correction with Asymmetric Least Squares Smoothing [1]. Their method, now known as Asymmetric Least Squares (AsLS), builds directly on the Whittaker smoother that Eilers had revived two years earlier [2],[3]. The idea fits in one sentence: take a smoother, and give the points above the curve almost no say in where the curve goes. A few rounds of this, and the smooth curve sinks below the peaks and settles onto the background.

The name comes from statistics, where "asymmetric least squares" had already been used to estimate expectiles, a close cousin of quantiles [4]. The report was never published in a journal, yet AsLS became one of the most widely used baseline methods in analytical chemistry, and the starting point of a whole family of variants such as airPLS and arPLS [5],[6]. Later in this article we also look at Local Asymmetric Least Squares (LAsLS), an extension that lets the AsLS settings change from one region of the signal to another [8].

A measured spectrum is usually the sum of two very different things:

The background can come from many places:

• a curved hump from fluorescence in Raman spectra,
• a slope from scattering in diffuse reflectance,
• slow instrumental drift between measurements,
• broad contributions from the sample holder, solvent or matrix.

The peaks are what we care about. The background is not, and if it changes from sample to sample it can easily dominate a PCA or a PLS model.

So we want a curve that follows the slow background but ignores the peaks. An ordinary smoother cannot do that: it treats every point the same way, so it is pulled up by the peaks and ends up running through their middle. We need a smoother that knows the difference between "a point on the background" and "a point on a peak". AsLS gets that knowledge from a simple observation: in the signals it is designed for, peaks stick up above the background, never below it. (For signals with negative peaks, such as transmittance spectra, flip the sign of the signal first.)

The Whittaker smoother looks for a smooth curve $\mathbf{z}$ that stays close to the data $\mathbf{y}$ [2],[3]. It does that by minimizing a sum of two terms:

The first term is a familiar sum of squared residuals. The second term measures roughness. The symbol ${\Delta }^2{z}_i$ is a second difference:

It is the discrete version of a second derivative, so it is large where the curve bends sharply and zero where the curve is a straight line. The parameter $\lambda$ sets the balance between the two wishes. A small $\lambda$ lets the curve follow the data closely; a large $\lambda$ forces it to be stiff and gently curved. (The Whittaker smoothing article explores this balance in detail.)

This is excellent for smoothing, but for baseline correction it has exactly the wrong symmetry: a point that is far above the curve pulls just as hard as a point that is far below it.

AsLS adds one ingredient: a weight $w_i$ for every point.

A point with a large weight pulls the curve strongly towards itself. A point with a tiny weight is almost ignored. The trick is in how the weights are chosen [1]:

with $p$ a small number, typically between 0.001 and 0.1. With $p=0.01$, a point below the curve counts 99 times more than a point above it.

There is a chicken-and-egg problem here: to set the weights we need the curve, and to find the curve we need the weights. AsLS solves it the simplest possible way, by alternating.

1. Start with all weights equal to 1. The first fit is therefore an ordinary Whittaker smooth, which runs through the middle of the peaks.
2. Fit the smooth curve with the current weights.
3. Reweight: points above the curve get weight $p$, points below get $1-p$.
4. Repeat steps 2 and 3. After each round the peaks lose more influence and the curve sinks lower. Usually 5 to 10 rounds are enough, and the weights stop changing. More iterations are not always better, though: with a very small $p$ and a stiff baseline, the curve can keep sinking slowly below the background.

Each fit in step 2 is a linear least squares problem, just like the ones in the least squares section. Write the weights on the diagonal of a matrix $\mathbf{W}$ and collect the second differences in a matrix $\mathbf{D}$, so that $\mathbf{Dz}$ is the vector of all ${\Delta }^2{z}_i$. Setting the gradient of $S$ to zero gives:

That is a set of normal equations. It looks large (one unknown per data point), but the matrix is banded: each row only has five non-zero entries, so sparse solvers handle thousands of points in milliseconds [2].

The spectrum below is simulated, so we know the true background (green, dashed) and the true peaks. Two numbers at the bottom tell you how well the correction worked: the error of the baseline, and how much of the true peak area survives in the corrected spectrum. Start with the defaults, then:

1. Move Iterations to 1. That is a plain Whittaker smooth: the curve runs straight through the peaks and the peak area recovered drops to about half.
2. Step the iterations up one by one and watch the curve sink below the peaks as the weights are updated.
3. Lower $\lambda$ to about $10^2$. The baseline becomes so flexible that it climbs into the peaks and eats a large part of them.
4. Raise $\lambda$ to $10^9$. The baseline is now so stiff that it cannot follow the hump, and a lot of background is left behind (the "recovered" area goes far above 100 %).
5. Push $p$ up to 0.5. The weights become symmetric and AsLS turns back into an ordinary smoother.
6. Try the other background shapes. The straight slope survives any stiffness, even $\lambda =10^9$. The wavy drift has two full waves, and at the default $10^5$ the baseline is already too stiff to follow them: lower $\lambda$ to about $10^4$.

Notice that the error is small over a fairly wide middle range of both parameters (here roughly $10^4$ to $10^6$ for $\lambda$). That is typical, and it is part of why AsLS is popular: it is forgiving.

AsLS has two parameters that matter and two that rarely need attention.

$\lambda$, the smoothness. Larger values give a stiffer baseline. Eilers and Boelens suggest searching between $10^2$ and $10^9$ [1], which is a much larger range than for ordinary smoothing because the baseline has to be much smoother than the peaks. A good starting point is $10^5$ to $10^6$.

• If the baseline climbs into the peaks or follows broad bands, increase $\lambda$.
• If it cannot follow the curvature of the background, decrease $\lambda$.

Keep in mind that the penalty is a sum over neighbouring points. If you measure the same spectrum with twice as many points, the same physical curvature needs a $\lambda$ about $2^{2d}=16$ times larger (for $d=2$). The reason: halving the point spacing $h$ doubles the number of terms in the fit, while each second difference shrinks by a factor of four, so $\lambda$ scales with $1/h^4$. Parameters tuned on one instrument do not automatically transfer to another with a different resolution.

$p$, the asymmetry. Smaller values ignore the peaks more strongly. The usual range is 0.001 to 0.1 [1].

• If the baseline is pulled up into the peaks, decrease $p$.
• If the baseline sits clearly below the noisy background, increase $p$ a little.

The difference order $d$ is almost always 2, as above. The number of iterations is usually 10; the weights typically stop changing well before that.

If the corrected spectra feed a calibration or classification model, you can also let the model choose: try a grid of $\lambda$ and $p$ values and keep the pair that gives the lowest cross-validated prediction error, as proposed by Liland, Almøy and Mevik [9].

The implementations below follow the original report [1] and use sparse matrices, so they are fast even for long spectra.

import numpy as np
from scipy import sparse
from scipy.sparse.linalg import spsolve
import matplotlib.pyplot as plt

def baseline_asls(y, lam=1e6, p=0.01, niter=10, d=2):
    """
    Asymmetric Least Squares baseline (Eilers and Boelens, 2005).

    y     : 1D spectrum
    lam   : smoothness (lambda). Typical range 1e2 to 1e9.
    p     : asymmetry. Typical range 0.001 to 0.1.
    niter : number of reweighting iterations.
    d     : order of the differences (usually 2).
    Returns the fitted baseline z.
    """
    y = np.asarray(y, dtype=float)
    n = len(y)

    # Difference matrix of order d, shape (n - d, n)
    D = sparse.eye(n, format='csr')
    for _ in range(d):
        D = D[1:] - D[:-1]
    DTD = lam * (D.T @ D)

    w = np.ones(n)
    for _ in range(niter):
        W = sparse.diags(w, 0, shape=(n, n), format='csc')
        z = spsolve(W + DTD, w * y)
        w = np.where(y > z, p, 1 - p)   # asymmetric reweighting
    return z

# Example: Raman-like spectrum on a fluorescence background
rng = np.random.default_rng(42)
shift = np.linspace(400, 1800, 700)
peaks = (0.8 * np.exp(-((shift - 600) ** 2) / 500) +
         0.6 * np.exp(-((shift - 900) ** 2) / 400) +
         0.5 * np.exp(-((shift - 1200) ** 2) / 600))
background = 3 + 0.002 * shift - 1e-6 * shift ** 2
spectrum = background + peaks + rng.normal(0, 0.02, len(shift))

baseline = baseline_asls(spectrum, lam=1e6, p=0.01)
corrected = spectrum - baseline

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(9, 7), sharex=True)
ax1.plot(shift, spectrum, label='Spectrum')
ax1.plot(shift, baseline, label='AsLS baseline', lw=2)
ax1.legend()
ax2.plot(shift, corrected, color='green')
ax2.axhline(0, color='k', ls='--', lw=1)
ax2.set_xlabel('Raman shift (cm$^{-1}$)')
plt.tight_layout()
plt.show()

% Example
rng(42);
shift = linspace(400, 1800, 700)';
peaks = 0.8 * exp(-((shift - 600).^2) / 500) + ...
        0.6 * exp(-((shift - 900).^2) / 400) + ...
        0.5 * exp(-((shift - 1200).^2) / 600);
background = 3 + 0.002 * shift - 1e-6 * shift.^2;
spectrum = background + peaks + 0.02 * randn(size(shift));

baseline = baseline_asls(spectrum, 1e6, 0.01, 10, 2);
corrected = spectrum - baseline;

subplot(2, 1, 1); plot(shift, spectrum, shift, baseline, 'LineWidth', 1.5);
legend('Spectrum', 'AsLS baseline');
subplot(2, 1, 2); plot(shift, corrected, 'LineWidth', 1.5); yline(0, '--');
xlabel('Raman shift (cm^{-1})');

function z = baseline_asls(y, lam, p, niter, d)
    % Asymmetric Least Squares baseline (Eilers and Boelens, 2005)
    if nargin < 5, d = 2; end
    if nargin < 4, niter = 10; end
    if nargin < 3, p = 0.01; end
    if nargin < 2, lam = 1e6; end

    y = y(:);
    n = length(y);
    D = diff(speye(n), d);          % (n - d) x n difference matrix
    DTD = lam * (D' * D);

    w = ones(n, 1);
    for it = 1:niter
        W = spdiags(w, 0, n, n);
        z = (W + DTD) \ (w .* y);
        w = p * (y > z) + (1 - p) * (y <= z);
    end
end

library(Matrix)

baseline_asls <- function(y, lambda = 1e6, p = 0.01, niter = 10, d = 2) {
  # Asymmetric Least Squares baseline (Eilers and Boelens, 2005)
  n <- length(y)
  D <- Matrix(diff(diag(n), differences = d), sparse = TRUE)
  DTD <- lambda * crossprod(D)

  w <- rep(1, n)
  for (it in 1:niter) {
    W <- Diagonal(x = w)
    z <- as.vector(solve(W + DTD, w * y))
    w <- ifelse(y > z, p, 1 - p)
  }
  z
}

# Example
set.seed(42)
shift <- seq(400, 1800, length.out = 700)
peaks <- 0.8 * exp(-((shift - 600)^2) / 500) +
         0.6 * exp(-((shift - 900)^2) / 400) +
         0.5 * exp(-((shift - 1200)^2) / 600)
background <- 3 + 0.002 * shift - 1e-6 * shift^2
spectrum <- background + peaks + rnorm(length(shift), 0, 0.02)

baseline <- baseline_asls(spectrum, lambda = 1e6, p = 0.01)
corrected <- spectrum - baseline

par(mfrow = c(2, 1))
plot(shift, spectrum, type = "l"); lines(shift, baseline, col = "red", lwd = 2)
plot(shift, corrected, type = "l", col = "darkgreen"); abline(h = 0, lty = 2)

AsLS is simple and robust, but its weights are crude: every point is either "peak" (weight $p$) or "background" (weight $1-p$), with nothing in between. Two consequences follow.

• Noise is treated as peak. On a noisy background, about half of the points sit above the true baseline just because of noise. AsLS down-weights them, so the fitted curve drifts towards the lower edge of the noise band and ends up too low. With $p=0.01$, the baseline typically sits one to two noise standard deviations $\sigma$ below the true background. This is where the expectiles mentioned at the start come in: for pure Gaussian noise and a very stiff baseline, the asymmetric weights estimate the 1 % expectile of the noise, about $1.7\sigma$ below its centre (about $0.9\sigma$ for $p=0.1$ and $2.4\sigma$ for $p=0.001$).
• You must tune $p$. The right value depends on how tall the peaks are and how noisy the spectrum is.

Several later methods keep the same penalized least squares engine and only change the way the weights are computed. They are often called penalized least squares (PLS) baseline methods, not to be confused with partial least squares regression.

airPLS (adaptive iteratively reweighted penalized least squares) gives zero weight to every point above the curve and lets the weight of points below it grow with their distance from the curve, more strongly at each iteration [5]. There is no $p$ to tune, only $\lambda$.

arPLS (asymmetrically reweighted penalized least squares) estimates the noise level from the points below the curve and uses a smooth, S-shaped (logistic) weight function [6]. Points that are above the curve by only about the noise level keep a reasonable weight, so the baseline runs through the middle of the noise instead of along its bottom edge.

Many more variants exist, for example drPLS [7] or improved AsLS (IAsLS) [10], but they all share the structure you have just learned: a Whittaker smoother inside a loop that keeps updating the weights.

You rarely need to write these yourself. In Python, the pybaselines package implements AsLS and most of its relatives (asls, airpls, arpls, drpls, iasls and others). In R, the baseline package offers AsLS as baseline.als [9]. The code in this article is still worth reading, because it shows exactly what those functions do.

Raise the noise and compare the error of each method. With clean spectra and a moderate $\lambda$ the three are almost indistinguishable. With noisy spectra AsLS settles too low (at the default noise level, by a bit more than one noise standard deviation), and airPLS, which gives zero weight to every point above the curve, sinks with it; arPLS stays close to the true background. Now lower the noise to the minimum and raise $\lambda$: airPLS and arPLS are more sensitive than AsLS to a $\lambda$ that is too stiff for the background, and above about $10^{6.5}$ they collapse on this example while AsLS degrades more gradually. Real data are rarely this tidy, so in practice it is worth trying two or three of these methods on your own spectra and looking at the result.

Everything so far uses one $\lambda$ and one $p$ for the whole signal. That is fine when the background looks roughly the same everywhere. But many real signals are not like that. Picture a chromatogram that starts with a nearly flat, noisy stretch with a few small peaks, and later climbs steeply while tall and small peaks sit on top of the rise.

Now try to choose a single setting:

• The steep rise needs a fairly flexible baseline (small $\lambda$), but the flat region would be better served by a stiff one.
• The tall peaks need a very small $p$ so they cannot lift the baseline. In the noisy region, the same small $p$ treats every noise spike above the curve as a peak and pushes the baseline down to the lower edge of the noise.

Whatever you pick is a compromise, and the compromise is wrong somewhere.

Local Asymmetric Least Squares (LAsLS), published in 2026 by Gómez-Sánchez, Torres-Cobos and Rocha de Oliveira, removes that compromise [8]. (Its official code is published by Lovelace's Square, the project behind this Library.) It keeps exactly the same penalized least squares engine as AsLS, but lets the asymmetric weights and the smoothing penalty change across user-defined regions of the signal. The authors present it for chromatography, Raman and mid-infrared signals.

With LAsLS you mark one or more intervals on the signal, typically the regions that contain peaks. Each interval $j$ gets its own pair of AsLS settings:

• an asymmetry $p_j$, used for the weights of the points inside that interval,
• a smoothness ${\lambda }_j$, used for the penalty inside that interval.

Outside every interval there are, by assumption, no peaks. There, LAsLS uses symmetric weights (every point has weight 1) and a smoothness ${\lambda }_{\text{Whit}}$. In other words, outside the intervals the baseline is an ordinary Whittaker smoother, which runs through the middle of the noise instead of along its bottom edge. Finally, one extra parameter $\mu$ adds a gentle "tension" to the whole baseline.

Two special cases make the relationship with what you already know clear. With no intervals at all (and $\mu =0$), LAsLS is a Whittaker smoother. With a single interval that covers the whole signal (and $\mu =0$), it is AsLS with $\lambda =\max ({\lambda }_{\text{Whit}},{\lambda }_1)$, apart from the stopping rule (LAsLS stops when the baseline stops changing, while the AsLS code above runs a fixed number of iterations).

LAsLS minimizes a sum of three terms:

Compare it with AsLS. The first term is the same weighted fit to the data. The second term is the same second-difference penalty (written here centred on point $i$), except that $\lambda$ now carries an index: each point has its own amount of smoothing. The third term is new. It uses first differences, $z_{i+1}-z_i$, so it penalizes slope rather than curvature and pulls the baseline towards flatter shapes.

The weights follow the AsLS rule inside the intervals and are fixed to 1 outside them:

The local smoothness is ${\lambda }_{\text{Whit}}$ everywhere, raised to ${\lambda }_j$ inside interval $j$:

Two details of the official implementation are worth knowing. Because of the maximum, a ${\lambda }_j$ smaller than ${\lambda }_{\text{Whit}}$ has no effect: inside an interval the baseline is never more flexible than outside. And when intervals overlap, the larger ${\lambda }_j$ wins, while the $p_j$ of the interval listed last is used.

Setting the gradient of $S$ to zero gives a linear system that looks just like the AsLS one, with two extra pieces:

Here $\mathbf{W}$ holds the weights on its diagonal, $\mathbf{D}$ is the second-difference matrix from before, $\mathbf{\Lambda }$ is a diagonal matrix with the local ${\lambda }_i$ values, and $\mathbf{L}$ is the first-difference matrix. The system matrix is still banded, so it is solved as fast as plain AsLS.

The iterations are the same iteratively reweighted least squares (IRLS) loop as before: solve, update the weights, solve again.

1. Start with all weights equal to 1.
2. Solve the system above for $\mathbf{z}$.
3. Check convergence: stop when the baseline barely changes, that is when $\Vert {\mathbf{z}}_{\text{new}}-{\mathbf{z}}_{\text{old}}\Vert /\Vert {\mathbf{z}}_{\text{old}}\Vert$ falls below a tolerance, or after a maximum number of iterations.
4. Reweight the points inside the intervals with their own $p_j$ (points outside keep weight 1) and go back to step 2.

The signal below is a simulated chromatogram with exactly the problem described above: a noisy, nearly flat start and a steep rise later on, with tall and small peaks. Plot 1a uses global AsLS, plot 1b uses LAsLS with two intervals (shaded). The bottom panel shows how far each baseline is from the true one.

1. Start with the defaults and compare the two errors at the bottom.
2. Try to make global AsLS as good as LAsLS. A small $p$ protects the tall peaks but sinks the baseline below the noise at the start; a $p$ large enough to keep it in the middle of the noise lets the tall peaks lift it. $\lambda$ has the same kind of conflict: flexible enough for the rise, or stiff enough to stay out of the peaks. In fact the AsLS defaults are already close to the best single setting we could find by searching the whole grid (error ≈ 0.031). LAsLS gets to ≈ 0.012, a ratio of about 2.7, inside the 2.6 to 5.7 range reported in the paper.
3. Lower $p_1$ to about 0.001. Inside interval 1 the LAsLS baseline now also sinks to the lower edge of the noise, just like AsLS.
4. Move the start of interval 1 to about 3 minutes, so that the first small peak falls outside. The symmetric smoother there treats the peak as background and bends up into it. The overall error hardly changes, so look at the error plot around 2.6 minutes.
5. Change ${\lambda }_2$ and watch the steep region respond: too large and the baseline can no longer follow the rise. Then lower it step by step. Once it is at or below ${\lambda }_{\text{Whit}}$, nothing changes any more, because of the maximum in the formula above.
6. Raise $\mu$. The tension pulls the whole baseline towards straight lines. Here the background really does rise steeply, so larger values soon stop the baseline from following it.

Use AsLS when the background behaves the same way across the whole signal: similar noise, similar curvature, peaks of comparable size. It has fewer parameters, needs no intervals, and is easy to apply to large datasets.

Consider LAsLS when the characteristics of the baseline are clearly non-uniform: noisy regions next to clean ones, flat stretches next to steep ones, tall peaks next to small ones. The original paper tested this situation on simulated chromatogram-like signals with heterogeneous baselines [8]. Across four cases, the relative RMSE of the AsLS baseline was 2.6 to 5.7 times higher than that of LAsLS, and the relative RMSE of peak quantification was 1.4 to 3.5 times higher. In a controlled PLSR example, LAsLS also reduced the external prediction error (relative RMSE of 11.3% against 17.9% with AsLS).

The price is that you must say where the regions are. LAsLS does not find the intervals for you, so it is most useful when you know your signal: where the peaks elute or appear, and where the background is free of them.

The advice below follows from how the method is built; it is not a set of values taken from the paper. Treat it as a starting point and, as with AsLS, judge the result by plotting the baselines.

Intervals. Put every region that contains peaks inside an interval. Any peak left outside is fitted by the symmetric smoother, which will bend the baseline up into it. Leave peak-free stretches outside: there the symmetric weights give an unbiased baseline through the noise. If peak positions shift a little between samples, make the intervals wide enough to contain the peaks in every sample, and use the same intervals for the whole dataset. In the official function the intervals are given as start and end indices of the signal.

$p_j$. It plays the same role as $p$ in AsLS, region by region. Where peaks are tall compared with the noise, a small value (for example 0.001 to 0.01) keeps them from lifting the baseline. Where the noise is large compared with the peaks, a larger value keeps the baseline from sinking below the noise.

${\lambda }_j$. The stiffness inside interval $j$. Use a larger value where the baseline under the peaks is flat or gently curved, and a smaller one where it bends quickly. Remember that only values above ${\lambda }_{\text{Whit}}$ make a difference.

${\lambda }_{\text{Whit}}$. The smoothness of the baseline outside the intervals, and the minimum smoothness everywhere. It has to be large enough to average out the noise in the peak-free stretches and small enough to follow the background there.

$\mu$. The global first-difference penalty. With $\mu =0$ the term disappears. Larger values pull the whole baseline towards flatter, straighter shapes, which can steady it but also stop it from following real slopes.

Iterations and tolerance. The official function stops when the relative change of the baseline drops below a tolerance $\text{tol}$ or after $\text{maxIter}$ iterations; its documentation example uses 50 iterations and $10^{-6}$.

The official implementation is the MATLAB function LASLS_CL, released under the MIT license in the lasls-cl repository on GitHub, and also listed on The Square. The Python and R versions below are line-by-line ports of the same algorithm (same operators, same interval rules, same stopping criterion). The example signal is a chromatogram-like trace like the one in the demo.

import warnings
import numpy as np
from scipy import sparse
from scipy.sparse.linalg import spsolve
import matplotlib.pyplot as plt

def lasls(y, intervals, p_vals, lambdas_asym, lambda_whit, mu,
          max_iter=50, tol=1e-6):
    """
    Local Asymmetric Least Squares (LAsLS) baseline.
    Python port of LASLS_CL.m (github.com/LovelaceSquare/lasls-cl).

    y            : 1D signal (n points)
    intervals    : list of (start, end) index pairs, 0-based and inclusive
    p_vals       : asymmetry p_j for each interval
    lambdas_asym : smoothness lambda_j for each interval
    lambda_whit  : smoothness outside the intervals (symmetric weights)
    mu           : global first-difference ("tension") penalty
    max_iter     : maximum number of IRLS iterations
    tol          : stop when the relative change of the baseline < tol
    Returns (baseline, weights).
    """
    y = np.asarray(y, dtype=float).ravel()
    n = y.size

    # Local parameters: p_vec = -1 marks points outside every interval
    p_vec = -np.ones(n)
    lam_vec = np.full(n - 2, float(lambda_whit))
    for (start, end), p, lam in zip(intervals, p_vals, lambdas_asym):
        if start < 0 or end > n - 1 or start > end:
            warnings.warn(f"Skipping invalid interval [{start}, {end}].")
            continue
        p_vec[start:end + 1] = p                    # later intervals overwrite p
        r0, r1 = max(0, start - 1), min(n - 3, end - 1)
        lam_vec[r0:r1 + 1] = np.maximum(lam_vec[r0:r1 + 1], lam)  # overlaps: max

    # Second-difference (n-2 x n) and first-difference (n-1 x n) operators
    D = sparse.diags([1.0, -2.0, 1.0], [0, 1, 2], shape=(n - 2, n))
    L = sparse.diags([-1.0, 1.0], [0, 1], shape=(n - 1, n))
    P = D.T @ sparse.diags(lam_vec) @ D + mu * (L.T @ L)

    inside = p_vec > 0
    w = np.ones(n)
    z = np.zeros(n)
    for _ in range(max_iter):
        A = (sparse.diags(w) + P).tocsc()
        z_new = spsolve(A, w * y)
        rel_change = np.linalg.norm(z_new - z) / (np.linalg.norm(z) + np.finfo(float).eps)
        z = z_new
        if rel_change < tol:
            break
        # Asymmetric weights inside the intervals, weight 1 outside
        w = np.where(inside, np.where(y > z, p_vec, 1 - p_vec), 1.0)
    return z, w

# Example: chromatogram-like signal with a heterogeneous baseline
rng = np.random.default_rng(1)
t = np.linspace(0, 25, 600)                       # time (min)
true_baseline = 0.08 + 0.012 * t + 0.9 / (1 + np.exp(-(t - 16) / 1.1))

def peak(height, centre, width):
    return height * np.exp(-0.5 * ((t - centre) / width) ** 2)

peaks = (peak(0.22, 2.6, 0.12) + peak(0.30, 6.1, 0.13) + peak(0.12, 8.2, 0.16) +
         peak(2.2, 14.6, 0.16) + peak(0.9, 17.6, 0.18) + peak(0.2, 20.4, 0.15))
noise = np.where(t < 11, 0.035, 0.008) * rng.standard_normal(t.size)
y = true_baseline + peaks + noise

def idx(minutes):
    """Index of the time point closest to a given time."""
    return int(np.argmin(np.abs(t - minutes)))

# Interval 1: noisy, almost flat region with small peaks
# Interval 2: steep region with tall and small peaks
intervals = [(idx(1.5), idx(9.5)), (idx(13.5), idx(23.5))]
baseline, weights = lasls(y, intervals,
                          p_vals=[0.1, 0.005],
                          lambdas_asym=[1e6, 1e4],
                          lambda_whit=1e3,
                          mu=0.0)
corrected = y - baseline

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(9, 7), sharex=True)
ax1.plot(t, y, color='0.5', label='Signal')
ax1.plot(t, baseline, lw=2, label='LAsLS baseline')
for start, end in intervals:
    ax1.axvspan(t[start], t[end], alpha=0.1)
ax1.legend()
ax2.plot(t, corrected, color='green')
ax2.axhline(0, color='k', ls='--', lw=1)
ax2.set_xlabel('Time (min)')
plt.tight_layout()
plt.show()

% Official implementation: LASLS_CL.m from https://github.com/LovelaceSquare/lasls-cl
% addpath('path/to/lasls-cl');   % folder that contains LASLS_CL.m

% Example: chromatogram-like signal with a heterogeneous baseline
rng(1);
t = linspace(0, 25, 600)';                     % time (min), column vector
trueBaseline = 0.08 + 0.012 * t + 0.9 ./ (1 + exp(-(t - 16) / 1.1));
peak = @(h, c, s) h * exp(-0.5 * ((t - c) / s).^2);
peaks = peak(0.22, 2.6, 0.12) + peak(0.30, 6.1, 0.13) + peak(0.12, 8.2, 0.16) + ...
        peak(2.2, 14.6, 0.16) + peak(0.9, 17.6, 0.18) + peak(0.2, 20.4, 0.15);
noise = (0.035 * (t < 11) + 0.008 * (t >= 11)) .* randn(size(t));
y = trueBaseline + peaks + noise;              % y must be a column vector

% Intervals as [startIdx, endIdx] (1-based indices, one row per interval)
idx = @(minutes) find(abs(t - minutes) == min(abs(t - minutes)), 1);
intervals   = [idx(1.5),  idx(9.5); ...      % noisy, almost flat region
               idx(13.5), idx(23.5)];        % steep region, tall and small peaks
pVals       = [0.1; 0.005];                  % p_j for each interval
lambdasAsym = [1e6; 1e4];                    % lambda_j for each interval
lambdaWhit  = 1e3;                           % smoothness outside the intervals
mu          = 0;                             % global first-difference penalty
maxIter     = 50;
tol         = 1e-6;

[baseline, weights] = LASLS_CL(y, intervals, pVals, lambdasAsym, ...
                               lambdaWhit, mu, maxIter, tol);
corrected = y - baseline;

subplot(2, 1, 1); plot(t, y, t, baseline, 'LineWidth', 1.5);
legend('Signal', 'LAsLS baseline');
subplot(2, 1, 2); plot(t, corrected, 'LineWidth', 1.5); yline(0, '--');
xlabel('Time (min)');

library(Matrix)

lasls <- function(y, intervals, p_vals, lambdas_asym, lambda_whit, mu,
                  max_iter = 50, tol = 1e-6) {
  # Local Asymmetric Least Squares (LAsLS) baseline.
  # R port of LASLS_CL.m (github.com/LovelaceSquare/lasls-cl).
  # intervals: matrix with one row per interval, c(start, end),
  #            1-based inclusive indices (as in MATLAB)
  y <- as.numeric(y)
  n <- length(y)

  # Local parameters: p_vec = -1 marks points outside every interval
  p_vec <- rep(-1, n)
  lam_vec <- rep(lambda_whit, n - 2)
  for (j in seq_along(p_vals)) {
    i_start <- intervals[j, 1]
    i_end <- intervals[j, 2]
    if (i_start < 1 || i_end > n || i_start > i_end) {
      warning(sprintf("Skipping invalid interval [%d, %d].", i_start, i_end))
      next
    }
    p_vec[i_start:i_end] <- p_vals[j]               # later intervals overwrite p
    rows <- max(1, i_start - 1):min(n - 2, i_end - 1)
    lam_vec[rows] <- pmax(lam_vec[rows], lambdas_asym[j])  # overlaps: max
  }

  # Second-difference (n-2 x n) and first-difference (n-1 x n) operators
  D <- Matrix(diff(diag(n), differences = 2), sparse = TRUE)
  L <- Matrix(diff(diag(n), differences = 1), sparse = TRUE)
  P <- t(D) %*% Diagonal(x = lam_vec) %*% D + mu * crossprod(L)

  w <- rep(1, n)
  z <- rep(0, n)
  for (iter in seq_len(max_iter)) {
    A <- Diagonal(x = w) + P
    z_new <- as.vector(solve(A, w * y))
    rel_change <- sqrt(sum((z_new - z)^2)) / (sqrt(sum(z^2)) + .Machine$double.eps)
    z <- z_new
    if (rel_change < tol) break
    # Asymmetric weights inside the intervals, weight 1 outside
    w <- ifelse(p_vec > 0, ifelse(y > z, p_vec, 1 - p_vec), 1)
  }
  list(baseline = z, weights = w)
}

# Example: chromatogram-like signal with a heterogeneous baseline
set.seed(1)
tt <- seq(0, 25, length.out = 600)              # time (min); 't' is R's transpose
true_baseline <- 0.08 + 0.012 * tt + 0.9 / (1 + exp(-(tt - 16) / 1.1))
peak <- function(h, c, s) h * exp(-0.5 * ((tt - c) / s)^2)
peaks <- peak(0.22, 2.6, 0.12) + peak(0.30, 6.1, 0.13) + peak(0.12, 8.2, 0.16) +
         peak(2.2, 14.6, 0.16) + peak(0.9, 17.6, 0.18) + peak(0.2, 20.4, 0.15)
noise <- ifelse(tt < 11, 0.035, 0.008) * rnorm(length(tt))
y <- true_baseline + peaks + noise

idx <- function(minutes) which.min(abs(tt - minutes))
intervals <- rbind(c(idx(1.5), idx(9.5)),       # noisy, almost flat region
                   c(idx(13.5), idx(23.5)))     # steep region, tall and small peaks

fit <- lasls(y, intervals, p_vals = c(0.1, 0.005),
             lambdas_asym = c(1e6, 1e4), lambda_whit = 1e3, mu = 0)
corrected <- y - fit$baseline

par(mfrow = c(2, 1))
plot(tt, y, type = "l", col = "grey40")
lines(tt, fit$baseline, col = "blue", lwd = 2)
plot(tt, corrected, type = "l", col = "darkgreen", xlab = "Time (min)")
abline(h = 0, lty = 2)

Use the same parameters for every spectrum. If you are building a calibration or classification model, the training, validation and test spectra must all be corrected with the same $\lambda$ and $p$. Tuning each spectrum separately creates differences that have nothing to do with chemistry.

Look at the baselines, not only the corrected spectra. Plot a handful of spectra with their baselines on top, including the most extreme ones (highest fluorescence, strongest peaks). A baseline that cuts into a broad band is easy to miss in the corrected spectrum and very visible in the overlay.

Broad bands are the hard case. AsLS separates "narrow" from "broad" through $\lambda$. If a real chemical band is nearly as broad as the background, no value of $\lambda$ will separate them perfectly, and the method will remove part of the band. Keep this in mind for water bands in NIR or for amorphous materials in Raman.

Corrected spectra can go slightly negative. Because of noise, a well-placed baseline leaves some points below zero. That is expected and not a sign of failure.

AsLS is a Whittaker smoother with lopsided weights. Points above the curve are assumed to be peaks and are almost ignored; points below it are assumed to be background and are followed closely. Repeating the fit a few times lets the curve sink onto the background.

Two parameters control it: $\lambda$ for the stiffness of the baseline and $p$ for how strongly the peaks are ignored. Each fit is a banded linear least squares problem, which is why it is fast. And when noise or tuning become a problem, the same engine with smarter weights (airPLS, arPLS and friends) is usually the next step.

When the background changes character along the signal, so that no single pair of settings fits everywhere, LAsLS keeps the same engine but lets you assign $p$ and $\lambda$ region by region, with a symmetric smoother in between.