Baseline Correction with AsLS and LAsLS

Baseline drift has been a persistent headache in spectroscopy for decades. In Raman spectroscopy, fluorescence interference produces broad, curved backgrounds that dwarf the actual Raman peaks. In infrared and NIR spectroscopy, scattering and instrumental effects introduce sloping baselines that shift from sample to sample. Before the mid-2000s, correcting these baselines often meant either fitting polynomials by hand (selecting baseline points manually, then hoping the polynomial behaved between them) or using rigid automated methods that struggled with the diversity of real-world spectra. Neither approach scaled well to large datasets or complex baselines.

In 2005, Paul H. C. Eilers and Hans F. M. Boelens at the Leiden University Medical Centre proposed a remarkably elegant solution in a technical report titled “Baseline Estimation by Quantile Regression.” Their method, now widely known as Asymmetric Least Squares (AsLS), built directly on Eilers’ earlier work on the Whittaker smoother (2003), which itself revived and modernized a penalized smoothing idea from 1923. The key insight was deceptively simple: take the Whittaker smoother’s objective function and make the weights asymmetric. Points lying above the current fitted curve (likely peaks) receive very small weights, while points below it (likely baseline) receive large weights. Through a few iterations of this reweighting scheme, the fitted curve is pulled down beneath the peaks and settles onto the baseline.

What made AsLS particularly appealing was its minimal parameter set: a smoothness parameter $\lambda$ (inherited from the Whittaker smoother) and an asymmetry parameter $p$ controlling how aggressively peaks are ignored. No manual point selection, no polynomial degree to guess, and no assumptions about peak shapes. Despite being published only as an unpublished technical report rather than a peer-reviewed journal article, AsLS quickly became one of the most widely cited and implemented baseline correction methods in analytical chemistry, spawning a family of variants including airPLS (Zhang et al., 2010) and drPLS (Xu et al., 2019).

The baseline problem in spectroscopy

In spectroscopy, baseline drift is one of the most common—and frustrating—preprocessing challenges. Your spectrum might show:

• A sloping baseline from scattering
• A curved baseline from fluorescence interference
• Instrumental drift over time
• Sample-specific background contributions

The fundamental challenge: you want to fit the baseline (the slowly varying background) while ignoring the peaks (your actual signal of interest). Traditional smoothing methods like Whittaker can’t distinguish between baseline and peaks—they’ll smooth everything equally, giving you neither good baseline estimation nor preserved peak shapes.

From smoothing to baseline correction

Whittaker smoothing provides an excellent foundation. It minimizes:

$$S=\sum _{i=1}^n(y_i-z_i{)}^2+\lambda \sum _{i=1}^{n-d}({\Delta }^d{z}_i{)}^2$$

This works great for smoothing, where we want z to be close to y. But for baseline correction, we need something clever: asymmetric weighting.

The key insight? Points below the fitted curve are likely baseline (high weight, fit closely), while points above the curve are likely peaks (low weight, ignore them).

This simple but brilliant idea leads to Asymmetric Least Squares (AsLS) and its more advanced cousin, Local Asymmetric Least Squares (LAsLS).

---

Asymmetric Least Squares (AsLS)

AsLS extends Whittaker smoothing with asymmetric weights that force the fitted curve to “hug” the baseline from below while ignoring peaks above.

The AsLS objective function

We modify the Whittaker objective by adding weights w_i:

$$S=\sum _{i=1}^n{w}_i(y_i-z_i{)}^2+\lambda \sum _{i=1}^{n-d}({\Delta }^d{z}_i{)}^2$$

The weights are updated iteratively based on the current fit:

$$w_i=\{\begin{array}{ll}p& \text{if }{y}_i>z_i\text{ (peak)}\\ 1-p& \text{if }{y}_i\le z_i\text{ (baseline)}\end{array}$$

Where:

• p is the asymmetry parameter (typically 0.001 to 0.1)
• Small p → more asymmetric → ignores peaks more strongly
• Large p → less asymmetric → closer to symmetric smoothing

The algorithm

AsLS is an iterative method:

1. Initialize: Start with all weights w_i = 1
2. Solve: Solve the weighted system (W + λD^TD)z = Wy for z
3. Update weights: Set w_i = p if y_i > z_i, else w_i = 1-p
4. Repeat steps 2-3 until convergence (typically 10 iterations)

The result is a smooth baseline that passes below the peaks.

Parameters

λ (lambda): Controls smoothness of the baseline

• Larger λ → smoother baseline (less flexible)
• Typical range: 10⁵ to 10⁹ (much larger than for smoothing!)
• Start with λ = 10⁶ and adjust
• If baseline is too rough → increase λ
• If baseline doesn’t follow curvature → decrease λ

p: Controls asymmetry (how much to ignore peaks)

• Smaller p → more asymmetric, ignores peaks more
• Typical range: 0.001 to 0.1
• Start with p = 0.01 (1%)
• If baseline cuts into peaks → decrease p (e.g., 0.001)
• If baseline is similar height to peaks → increase p (e.g., 0.05)

d: Difference order

• Usually d = 2 (same as Whittaker smoothing)
• Penalizes curvature changes

niter: Number of iterations

• Typically 10 is sufficient
• Algorithm usually converges in 5-10 iterations

Quick start parameters

Recommended starting point:

• λ = 10⁶
• p = 0.01
• d = 2
• niter = 10

Then adjust λ and p based on visual inspection:

1. Plot your spectrum with the fitted baseline
2. If baseline is too rough → increase λ
3. If baseline cuts into peaks → decrease p
4. If baseline doesn’t follow baseline curves → decrease λ

Interactive visualization

Let’s explore how λ and p affect baseline correction. The visualization below shows a Raman spectrum with fluorescence baseline. You can adjust the parameters and see the results in real-time:

Try this:

1. Start with the default parameters (λ = 10⁶, p = 0.01)
2. Increase λ and watch the baseline become smoother
3. Decrease λ and see it follow the baseline curvature more closely
4. Decrease p to make it more asymmetric (ignore peaks more)
5. Increase p and observe how the baseline starts cutting into peaks
6. Toggle between “None” and “AsLS” to compare original vs corrected

Notice how you can find parameter values that fit the baseline well while preserving the peak shapes!

Code examples: AsLS baseline correction

Python

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 correction.

    Parameters:
    -----------
    y : array-like
        Input spectrum
    lam : float
        Smoothness parameter (lambda). Larger = smoother baseline.
        Typical range: 1e5 to 1e9
    p : float
        Asymmetry parameter. Smaller = more asymmetric.
        Typical range: 0.001 to 0.1
    niter : int
        Number of iterations (typically 10)
    d : int
        Order of differences (typically 2)

    Returns:
    --------
    z : array
        Fitted baseline
    """
    n = len(y)

    # Difference matrix
    D = sparse.diags([1, -2, 1], [0, 1, 2], shape=(n-2, n), format='csc')
    for _ in range(d - 2):
        D = D[1:] - D[:-1]

    # Precompute D'D
    DTD = D.T @ D

    # Initialize weights
    w = np.ones(n)

    # Iterative reweighting
    for _ in range(niter):
        W = sparse.diags(w, 0, shape=(n, n), format='csc')
        Z = W + lam * DTD
        z = spsolve(Z, w * y)

        # Update weights asymmetrically
        w = p * (y > z) + (1 - p) * (y <= z)

    return z

# Example: Raman spectrum with fluorescence baseline
np.random.seed(42)
wavelength = np.linspace(400, 1800, 700)

# True Raman peaks
peaks = (0.8 * np.exp(-((wavelength - 600)**2) / 500) +
         0.6 * np.exp(-((wavelength - 900)**2) / 400) +
         0.5 * np.exp(-((wavelength - 1200)**2) / 600))

# Fluorescence baseline (curved)
baseline = 3 + 0.002 * wavelength - 1e-6 * wavelength**2

# Noisy spectrum
noise = np.random.normal(0, 0.05, len(wavelength))
spectrum = baseline + peaks + noise

# Apply AsLS baseline correction
fitted_baseline = baseline_asls(spectrum, lam=1e6, p=0.01, niter=10)
corrected = spectrum - fitted_baseline

# Plot
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8))

# Original spectrum with fitted baseline
ax1.plot(wavelength, spectrum, label='Spectrum', linewidth=1.5)
ax1.plot(wavelength, fitted_baseline, label='Fitted baseline',
         linewidth=2, color='red')
ax1.set_xlabel('Raman shift (cm⁻¹)')
ax1.set_ylabel('Intensity')
ax1.set_title('AsLS Baseline Correction')
ax1.legend()
ax1.grid(True, alpha=0.3)

# Baseline-corrected spectrum
ax2.plot(wavelength, corrected, linewidth=1.5, color='green')
ax2.axhline(0, color='black', linestyle='--', linewidth=1, alpha=0.5)
ax2.set_xlabel('Raman shift (cm⁻¹)')
ax2.set_ylabel('Intensity (corrected)')
ax2.set_title('After Baseline Correction')
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

MATLAB

function z = baseline_asls(y, lam, p, niter, d)
    % Asymmetric Least Squares baseline correction
    %
    % Parameters:
    %   y - Input spectrum (vector)
    %   lam - Smoothness (default: 1e6)
    %   p - Asymmetry (default: 0.01)
    %   niter - Iterations (default: 10)
    %   d - Difference order (default: 2)

    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

    n = length(y);

    % Difference matrix
    D = spdiags(ones(n-2, 1) * [1 -2 1], 0:2, n-2, n);
    for i = 1:(d-2)
        D = diff(D, 1, 1);
    end

    % Precompute D'D
    DTD = D' * D;

    % Initialize weights
    w = ones(n, 1);

    % Iterative reweighting
    for iter = 1:niter
        W = spdiags(w, 0, n, n);
        Z = W + lam * DTD;
        z = Z \ (w .* y(:));

        % Update weights asymmetrically
        w = p * (y(:) > z) + (1 - p) * (y(:) <= z);
    end
end

% Example
rng(42);
wavelength = linspace(400, 1800, 700);

% Raman peaks + fluorescence baseline
peaks = 0.8 * exp(-((wavelength - 600).^2) / 500) + ...
        0.6 * exp(-((wavelength - 900).^2) / 400) + ...
        0.5 * exp(-((wavelength - 1200).^2) / 600);
baseline = 3 + 0.002 * wavelength - 1e-6 * wavelength.^2;
spectrum = baseline + peaks + 0.05 * randn(size(wavelength));

% Apply AsLS
fitted_baseline = baseline_asls(spectrum, 1e6, 0.01, 10, 2);
corrected = spectrum - fitted_baseline;

% Plot
figure;
subplot(2,1,1);
plot(wavelength, spectrum, 'LineWidth', 1.5);
hold on;
plot(wavelength, fitted_baseline, 'r', 'LineWidth', 2);
xlabel('Raman shift (cm^{-1})');
ylabel('Intensity');
title('AsLS Baseline Correction');
legend('Spectrum', 'Fitted baseline');
grid on;

subplot(2,1,2);
plot(wavelength, corrected, 'g', 'LineWidth', 1.5);
hold on;
plot(wavelength, zeros(size(wavelength)), 'k--');
xlabel('Raman shift (cm^{-1})');
ylabel('Intensity (corrected)');
title('After Baseline Correction');
grid on;

R

baseline_asls <- function(y, lambda = 1e6, p = 0.01, niter = 10, d = 2) {
  # Asymmetric Least Squares baseline correction
  #
  # Parameters:
  #   y - Input spectrum
  #   lambda - Smoothness (default: 1e6)
  #   p - Asymmetry (default: 0.01)
  #   niter - Iterations (default: 10)
  #   d - Difference order (default: 2)

  library(Matrix)

  n <- length(y)

  # Difference matrix
  D <- diff(Diagonal(n), differences = d)
  DTD <- t(D) %*% D

  # Initialize weights
  w <- rep(1, n)

  # Iterative reweighting
  for (i in 1:niter) {
    W <- Diagonal(x = w)
    Z <- W + lambda * DTD
    z <- as.vector(solve(Z, w * y))

    # Update weights asymmetrically
    w <- p * (y > z) + (1 - p) * (y <= z)
  }

  return(z)
}

# Example
set.seed(42)
wavelength <- seq(400, 1800, length.out = 700)

# Raman peaks + fluorescence baseline
peaks <- 0.8 * exp(-((wavelength - 600)^2) / 500) +
         0.6 * exp(-((wavelength - 900)^2) / 400) +
         0.5 * exp(-((wavelength - 1200)^2) / 600)
baseline <- 3 + 0.002 * wavelength - 1e-6 * wavelength^2
spectrum <- baseline + peaks + rnorm(length(wavelength), 0, 0.05)

# Apply AsLS
fitted_baseline <- baseline_asls(spectrum, lambda = 1e6, p = 0.01,
                                  niter = 10, d = 2)
corrected <- spectrum - fitted_baseline

# Plot
par(mfrow = c(2, 1))

plot(wavelength, spectrum, type = 'l', lwd = 1.5,
     xlab = 'Raman shift (cm⁻¹)', ylab = 'Intensity',
     main = 'AsLS Baseline Correction')
lines(wavelength, fitted_baseline, col = 'red', lwd = 2)
legend('topright', legend = c('Spectrum', 'Fitted baseline'),
       col = c('black', 'red'), lwd = c(1.5, 2))
grid()

plot(wavelength, corrected, type = 'l', col = 'green', lwd = 1.5,
     xlab = 'Raman shift (cm⁻¹)', ylab = 'Intensity (corrected)',
     main = 'After Baseline Correction')
abline(h = 0, lty = 2)
grid()

---

Local Asymmetric Least Squares (LAsLS): Adaptive baseline correction

AsLS works remarkably well for many spectroscopic baselines, but it has a limitation: it uses global parameters (λ and p) for the entire spectrum. This means:

• The same smoothness is applied everywhere
• The same asymmetry is applied everywhere

In real spectra, baseline complexity often varies locally:

• Some regions might need more smoothness (gentle baseline drift)
• Other regions might need more flexibility (rapid baseline changes)
• Peak heights may vary dramatically (requiring different asymmetry)

LAsLS (Local Asymmetric Least Squares) solves this by making the parameters adaptive—they adjust automatically to local signal characteristics.

The core innovation: local parameter adaptation

Instead of using global λ and p for all points, LAsLS computes local versions for each spectral region:

1. Local smoothness λ_i adapts to the local baseline curvature
2. Local asymmetry p_i adapts to the local peak height

This allows the algorithm to be:

• Very smooth in flat regions
• More flexible in regions with rapid baseline changes
• More asymmetric where peaks are tall
• Less asymmetric where peaks are small

The LAsLS objective function

The mathematical formulation extends AsLS:

$$S=\sum _{i=1}^n{w}_i(y_i-z_i{)}^2+\sum _{i=1}^{n-d}{\lambda }_i({\Delta }^d{z}_i{)}^2$$

Notice that λ is now λ_i—a different smoothness parameter for each difference term.

Local parameter computation

At each iteration, LAsLS computes local parameters in a window around each point:

Local smoothness λ_i:

$${\lambda }_i=\frac{{\lambda }_{\text{global}}}{1+\alpha \cdot \frac{|{\Delta }^d{z}_i|}{s_d}}$$

Where:

• s_d = robust estimate of local curvature (median absolute deviation)
• α = adaptation strength (typically 0.5 to 2)

This means: in regions with high curvature (rapid baseline changes), λ_i decreases → more flexibility.

Local asymmetry p_i:

$$p_i=\frac{p_{\text{global}}}{1+\beta \cdot \frac{\max (0,y_i-z_i)}{s_y}}$$

Where:

• s_y = robust estimate of signal scale (MAD of positive residuals)
• β = adaptation strength (typically 0.5 to 2)

This means: where peaks are tall (large positive residuals), p_i decreases → more asymmetric (ignore peaks more).

The LAsLS algorithm

1. Initialize: Compute initial baseline (use AsLS or Whittaker)
2. Compute local parameters:
   • Calculate local curvature and update λ_i for all i
   • Calculate local peak heights and update p_i for all i
3. Update weights: w_i = p_i if y_i > z_i, else w_i = 1-p_i
4. Solve: (W + Λ_localD^TD)z = Wy where Λ_local = diag(λ_i)
5. Repeat steps 2-4 until convergence (typically 10-15 iterations)

Computational cost

LAsLS is more computationally expensive than AsLS because:

1. It requires computing local statistics at each iteration
2. The system matrix changes each iteration (can’t precompute)
3. More iterations may be needed for convergence

However, for typical spectroscopic data (<10,000 points), the difference is negligible on modern computers. The improved baseline fit is usually worth it.

Parameters

LAsLS has more parameters than AsLS, but most can use default values:

Global parameters (same as AsLS):

• λ_global: Initial smoothness (start with 10⁵ to 10⁷)
• p_global: Initial asymmetry (start with 0.01)
• d: Difference order (typically 2)
• niter: Iterations (15-20 for LAsLS, more than AsLS)

Adaptation parameters (new):

• α: Smoothness adaptation strength (default: 1.0)
  • Larger α → more local smoothness variation
  • Typical range: 0.5 to 2.0
• β: Asymmetry adaptation strength (default: 1.0)
  • Larger β → more local asymmetry variation
  • Typical range: 0.5 to 2.0
• window_size: Window for computing local statistics (default: 10-50 points)
  • Depends on your spectral resolution
  • Should be larger than typical peak width

Start simple, adapt when needed

Recommended workflow:

1. Try AsLS first with λ = 10⁶, p = 0.01
2. If AsLS works well, you’re done (don’t overcomplicate)
3. If baseline fit is poor in some regions but good in others → try LAsLS
4. Start LAsLS with same λ and p as AsLS, use α = β = 1.0
5. Adjust α if smoothness needs to vary more
6. Adjust β if peak heights vary dramatically

LAsLS shines when you have heterogeneous baselines—different behavior in different spectral regions.

When to use LAsLS vs AsLS

Use AsLS when:

• Baseline character is similar across the spectrum
• You want a simple, fast method
• Parameters are easy to tune visually
• Your baseline is smoothly varying

Use LAsLS when:

• Baseline complexity varies significantly across the spectrum
• You have regions with very different curvature
• Peak heights vary dramatically (some tall, some small)
• AsLS either over-corrects some regions or under-corrects others
• You’re willing to invest time in parameter tuning

Real-world examples favoring LAsLS:

• Raman spectra with both sharp peaks and broad bands
• NIR spectra with variable scattering + absorption
• Chromatograms with baseline drift + very tall peaks
• Fluorescence with multiple emission bands of different intensities

Code examples: LAsLS baseline correction

Python

import numpy as np
from scipy import sparse
from scipy.sparse.linalg import spsolve
from scipy.ndimage import uniform_filter1d

def baseline_lasls(y, lam=1e6, p=0.01, alpha=1.0, beta=1.0,
                   niter=15, d=2, window=20):
    """
    Local Asymmetric Least Squares baseline correction.

    Parameters:
    -----------
    y : array-like
        Input spectrum
    lam : float
        Global smoothness parameter
    p : float
        Global asymmetry parameter
    alpha : float
        Smoothness adaptation strength (default: 1.0)
    beta : float
        Asymmetry adaptation strength (default: 1.0)
    niter : int
        Number of iterations (default: 15)
    d : int
        Difference order (default: 2)
    window : int
        Window size for local statistics (default: 20)

    Returns:
    --------
    z : array
        Fitted baseline
    """
    n = len(y)

    # Difference matrix
    D = sparse.diags([1, -2, 1], [0, 1, 2], shape=(n-2, n), format='csc')
    for _ in range(d - 2):
        D = D[1:] - D[:-1]

    DTD = D.T @ D

    # Initialize with AsLS
    w = np.ones(n)
    for _ in range(3):
        W = sparse.diags(w, 0, shape=(n, n), format='csc')
        Z = W + lam * DTD
        z = spsolve(Z, w * y)
        w = p * (y > z) + (1 - p) * (y <= z)

    # LAsLS iterations
    for iteration in range(niter):
        # Compute residuals
        residuals = y - z

        # Local smoothness adaptation
        # Compute local curvature
        if d == 2:
            curvature = np.abs(np.diff(z, n=2))
        else:
            curvature = np.abs(np.diff(z, n=d))

        # Pad to match length
        curvature = np.pad(curvature, (d//2, n - len(curvature) - d//2),
                          mode='edge')

        # Smooth curvature estimate
        local_curve = uniform_filter1d(curvature, size=window, mode='nearest')

        # Robust scale estimate
        mad_curve = np.median(np.abs(curvature - np.median(curvature)))
        if mad_curve == 0:
            mad_curve = 1.0

        # Adaptive lambda (smaller where curvature is high)
        lambda_local = lam / (1 + alpha * local_curve / mad_curve)

        # Local asymmetry adaptation
        # Compute local peak heights
        peak_heights = np.maximum(0, residuals)
        local_peaks = uniform_filter1d(peak_heights, size=window, mode='nearest')

        # Robust scale estimate
        mad_peaks = np.median(np.abs(peak_heights[peak_heights > 0]))
        if mad_peaks == 0:
            mad_peaks = 1.0

        # Adaptive p (smaller where peaks are tall)
        p_local = p / (1 + beta * local_peaks / mad_peaks)

        # Update weights with local asymmetry
        w = p_local * (y > z) + (1 - p_local) * (y <= z)

        # Build local smoothness matrix
        lambda_diff = lambda_local[:-d] if d > 0 else lambda_local
        Lambda = sparse.diags(lambda_diff, 0, format='csc')

        # Solve with local parameters
        W = sparse.diags(w, 0, shape=(n, n), format='csc')
        Z = W + D.T @ Lambda @ D
        z_new = spsolve(Z, w * y)

        # Check convergence
        if np.max(np.abs(z_new - z)) < 1e-6:
            break

        z = z_new

    return z

# Example: Complex baseline with variable features
np.random.seed(42)
wavelength = np.linspace(400, 2000, 800)

# Multiple peaks of different heights
peaks = (2.0 * np.exp(-((wavelength - 600)**2) / 400) +    # Tall peak
         0.3 * np.exp(-((wavelength - 900)**2) / 300) +    # Small peak
         1.5 * np.exp(-((wavelength - 1200)**2) / 500) +   # Medium peak
         0.4 * np.exp(-((wavelength - 1500)**2) / 350))    # Small peak

# Complex baseline: smooth in some regions, curved in others
baseline = (2 + 0.003 * wavelength - 2e-6 * wavelength**2 +
            0.5 * np.exp(-((wavelength - 1000)**2) / 50000))  # Local bump

# Noisy spectrum
noise = np.random.normal(0, 0.03, len(wavelength))
spectrum = baseline + peaks + noise

# Compare AsLS vs LAsLS
baseline_asls_fit = baseline_asls(spectrum, lam=1e6, p=0.01, niter=10)
baseline_lasls_fit = baseline_lasls(spectrum, lam=1e6, p=0.01,
                                     alpha=1.0, beta=1.0, niter=15)

corrected_asls = spectrum - baseline_asls_fit
corrected_lasls = spectrum - baseline_lasls_fit

# Plot comparison
fig, axes = plt.subplots(3, 1, figsize=(12, 10))

# Original with fitted baselines
axes[0].plot(wavelength, spectrum, label='Spectrum', linewidth=1.5, alpha=0.7)
axes[0].plot(wavelength, baseline_asls_fit, label='AsLS baseline',
             linewidth=2, color='red')
axes[0].plot(wavelength, baseline_lasls_fit, label='LAsLS baseline',
             linewidth=2, color='purple', linestyle='--')
axes[0].set_xlabel('Wavenumber (cm⁻¹)')
axes[0].set_ylabel('Intensity')
axes[0].set_title('Baseline Correction Comparison')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# AsLS corrected
axes[1].plot(wavelength, corrected_asls, linewidth=1.5, color='red')
axes[1].axhline(0, color='black', linestyle='--', linewidth=1, alpha=0.5)
axes[1].set_xlabel('Wavenumber (cm⁻¹)')
axes[1].set_ylabel('Intensity')
axes[1].set_title('AsLS Corrected')
axes[1].grid(True, alpha=0.3)

# LAsLS corrected
axes[2].plot(wavelength, corrected_lasls, linewidth=1.5, color='purple')
axes[2].axhline(0, color='black', linestyle='--', linewidth=1, alpha=0.5)
axes[2].set_xlabel('Wavenumber (cm⁻¹)')
axes[2].set_ylabel('Intensity')
axes[2].set_title('LAsLS Corrected')
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Print statistics
print(f"AsLS - Baseline residuals std: {np.std(corrected_asls):.4f}")
print(f"LAsLS - Baseline residuals std: {np.std(corrected_lasls):.4f}")

MATLAB

function z = baseline_lasls(y, lam, p, alpha, beta, niter, d, window)
    % Local Asymmetric Least Squares baseline correction
    %
    % Parameters:
    %   y - Input spectrum
    %   lam - Global smoothness (default: 1e6)
    %   p - Global asymmetry (default: 0.01)
    %   alpha - Smoothness adaptation (default: 1.0)
    %   beta - Asymmetry adaptation (default: 1.0)
    %   niter - Iterations (default: 15)
    %   d - Difference order (default: 2)
    %   window - Window size (default: 20)

    if nargin < 8, window = 20; end
    if nargin < 7, d = 2; end
    if nargin < 6, niter = 15; end
    if nargin < 5, beta = 1.0; end
    if nargin < 4, alpha = 1.0; end
    if nargin < 3, p = 0.01; end
    if nargin < 2, lam = 1e6; end

    n = length(y);
    y = y(:);

    % Difference matrix
    D = spdiags(ones(n-2, 1) * [1 -2 1], 0:2, n-2, n);
    for i = 1:(d-2)
        D = diff(D, 1, 1);
    end
    DTD = D' * D;

    % Initialize with AsLS
    w = ones(n, 1);
    for iter = 1:3
        W = spdiags(w, 0, n, n);
        Z = W + lam * DTD;
        z = Z \ (w .* y);
        w = p * (y > z) + (1 - p) * (y <= z);
    end

    % LAsLS iterations
    for iter = 1:niter
        residuals = y - z;

        % Local smoothness adaptation
        if d == 2
            curvature = abs(diff(z, 2));
        else
            curvature = abs(diff(z, d));
        end
        curvature = [curvature(1)*ones(floor(d/2),1); curvature; ...
                     curvature(end)*ones(n-length(curvature)-floor(d/2),1)];

        % Smooth curvature
        local_curve = movmean(curvature, window);
        mad_curve = median(abs(curvature - median(curvature)));
        if mad_curve == 0, mad_curve = 1.0; end

        % Adaptive lambda
        lambda_local = lam ./ (1 + alpha * local_curve / mad_curve);

        % Local asymmetry adaptation
        peak_heights = max(0, residuals);
        local_peaks = movmean(peak_heights, window);
        mad_peaks = median(abs(peak_heights(peak_heights > 0)));
        if isempty(mad_peaks) || mad_peaks == 0, mad_peaks = 1.0; end

        % Adaptive p
        p_local = p ./ (1 + beta * local_peaks / mad_peaks);

        % Update weights
        w = p_local .* (y > z) + (1 - p_local) .* (y <= z);

        % Local lambda matrix
        lambda_diff = lambda_local(1:end-d);
        Lambda = spdiags(lambda_diff, 0, length(lambda_diff), length(lambda_diff));

        % Solve with local parameters
        W = spdiags(w, 0, n, n);
        Z = W + D' * Lambda * D;
        z_new = Z \ (w .* y);

        % Check convergence
        if max(abs(z_new - z)) < 1e-6
            break;
        end
        z = z_new;
    end
end

% Example usage - see Python example for detailed context

R

baseline_lasls <- function(y, lambda = 1e6, p = 0.01, alpha = 1.0,
                           beta = 1.0, niter = 15, d = 2, window = 20) {
  # Local Asymmetric Least Squares baseline correction

  library(Matrix)
  library(zoo)  # For rollmean

  n <- length(y)

  # Difference matrix
  D <- diff(Diagonal(n), differences = d)
  DTD <- t(D) %*% D

  # Initialize with AsLS
  w <- rep(1, n)
  for (i in 1:3) {
    W <- Diagonal(x = w)
    Z <- W + lambda * DTD
    z <- as.vector(solve(Z, w * y))
    w <- p * (y > z) + (1 - p) * (y <= z)
  }

  # LAsLS iterations
  for (iter in 1:niter) {
    residuals <- y - z

    # Local smoothness adaptation
    curvature <- abs(diff(z, differences = d))
    curvature <- c(rep(curvature[1], floor(d/2)), curvature,
                   rep(curvature[length(curvature)],
                       n - length(curvature) - floor(d/2)))

    # Smooth curvature
    local_curve <- rollmean(curvature, k = window, fill = "extend")
    mad_curve <- median(abs(curvature - median(curvature)))
    if (mad_curve == 0) mad_curve <- 1.0

    # Adaptive lambda
    lambda_local <- lambda / (1 + alpha * local_curve / mad_curve)

    # Local asymmetry adaptation
    peak_heights <- pmax(0, residuals)
    local_peaks <- rollmean(peak_heights, k = window, fill = "extend")
    mad_peaks <- median(abs(peak_heights[peak_heights > 0]))
    if (length(mad_peaks) == 0 || mad_peaks == 0) mad_peaks <- 1.0

    # Adaptive p
    p_local <- p / (1 + beta * local_peaks / mad_peaks)

    # Update weights
    w <- p_local * (y > z) + (1 - p_local) * (y <= z)

    # Local lambda matrix
    lambda_diff <- lambda_local[1:(n-d)]
    Lambda <- Diagonal(x = lambda_diff)

    # Solve with local parameters
    W <- Diagonal(x = w)
    Z <- W + t(D) %*% Lambda %*% D
    z_new <- as.vector(solve(Z, w * y))

    # Check convergence
    if (max(abs(z_new - z)) < 1e-6) break
    z <- z_new
  }

  return(z)
}

# Example usage - see Python example for detailed context

---

Practical guidelines

Choosing between AsLS and LAsLS

Start with AsLS for most applications:

• Simpler with fewer parameters
• Faster computation
• Works well for uniform baselines
• Easier to understand and tune

Upgrade to LAsLS when AsLS fails:

• Baseline varies significantly across spectrum
• Some regions need different smoothness
• Peak heights vary dramatically
• AsLS over/under-corrects in different regions

Common applications

AsLS is ideal for:

• Standard Raman spectra with fluorescence
• NIR spectra with multiplicative scatter
• FTIR spectra with sloping baselines
• Chromatography with smooth baseline drift

LAsLS excels at:

• Complex Raman with varying peak heights
• Mixed spectral regions (some flat, some curved)
• Spectra with localized baseline artifacts
• Data with heterogeneous noise levels

Parameter tuning workflow

1. Start with defaults: λ = 10⁶, p = 0.01
2. Visual inspection: Plot spectrum + fitted baseline
3. Adjust smoothness (λ):
   • Too wiggly? Increase λ
   • Missing curves? Decrease λ
4. Adjust asymmetry (p):
   • Cuts into peaks? Decrease p
   • Doesn’t reach baseline? Increase p
5. Check residuals: Corrected spectrum should have flat baseline around zero

Consistency for quantitative analysis

Critical: If building calibration models, use the same parameters for all samples (training and test). Inconsistent baseline correction will hurt your predictions more than imperfect baseline removal.

Summary: AsLS vs LAsLS

Feature	AsLS	LAsLS
Complexity	Low (2-3 parameters)	High (4-6 parameters)
Speed	Fast (10 iterations)	Slower (15-20 iterations)
Best for	Uniform baselines	Variable baselines
Tuning	Easy (visual)	Moderate (requires testing)
Robustness	High	Moderate (can overfit)
Typical λ	105 – 109	105 – 107
Typical p	0.001 – 0.1	0.001 – 0.1

Decision rule:

1. Try AsLS first
2. If it works → stop (don’t overcomplicate)
3. If baseline fit is inconsistent across spectrum → try LAsLS
4. If LAsLS doesn’t help → your baseline may need manual correction or different preprocessing

Further reading

AsLS baseline correction:

• Eilers, P. H. C., & Boelens, H. F. M. (2005). “Baseline correction with asymmetric least squares smoothing”. Unpublished manuscript.
  • The original AsLS paper

LAsLS and adaptive variants:

• Zhang, Z.-M., Chen, S., & Liang, Y.-Z. (2010). “Baseline correction using adaptive iteratively reweighted penalized least squares”. Analyst, 135(5), 1138-1146.
  • airPLS method with adaptive parameters
• Xu, D., et al. (2019). “Baseline correction method based on doubly reweighted penalized least squares”. Applied Optics, 58(14), 3913-3920.
  • drPLS extension with robust weighting

Related methods:

• Whittaker smoothing (foundation for these methods)
• Polynomial baseline fitting
• Rubberband baseline correction
• Morphological baseline estimation

These methods are among the most powerful and widely-used baseline correction techniques in modern spectroscopy. Master them, and you’ll be equipped to handle the vast majority of baseline challenges in chemometric applications.