Open Lab/Local Modeling · Method

Locally Weighted Partial Least Squares Discriminant Analysis

LW-PLS-DA

Query-specific local PLS-DA: select n_local neighbors, apply Bevilacqua and Marini Table 1 distance weights through X_w = W X_l, and assign the class of the largest predicted dummy component.

ChemometricsPLS-DALW-PLS-DALocal ClassificationPythonMATLAB

LW-PLS-DA

What is LW-PLS-DA?

Bevilacqua and Marini 2014 extend locally weighted PLS to discrimination. For each unknown sample the algorithm finds the most similar calibration objects, optionally weights those neighbors by distance, and builds a PLS-DA model on that local set only. The next query requires a new local model.

This card implements that 2014 method. It is not the later Lesnoff / rchemo KNN-LWPLSR regression engine, and it is not Jchemo lwplsrda.

Why local classification?

The paper starts from the LWR idea that a globally nonlinear response can be treated as locally linear on a restricted domain. Global linear PLS-DA can fail under curvature, clustering, or masking of a middle class. Local PLS-DA is a strategy for those geometries. It does not prove that every nonlinear boundary is locally linear, and it does not guarantee improved prediction.

Query-specific modeling

Neighbors, weights, pretreatment, latent directions, and predicted dummy values can all change with the query. Fitting one global PLS-DA in fit and reusing it for every future sample is not LW-PLS-DA. The calibration database remains part of prediction.

PLS-DA response coding

Section 2.1 codes class membership as a binary dummy matrix $Y$ with one row per sample and one column per category. The row contains 1 in the column of the sample's class and 0 elsewhere. PLS predicts a real-valued $Y$ . The sample is assigned to the class of the highest predicted component. Those predicted values are not probabilities. This card does not use the Open Lab binary PLS-DA midpoint rule.

Finding the local neighborhood

Only the $n_{local}$ training objects closest to the query enter the local model. That count is a model parameter, selected together with PLS complexity by leave-one-out cross-validation. There is no universal rule such as 20 neighbors or a fraction of N. Dataset-specific grids in the paper (for example 50–300 on the sphere simulation) are experimental ranges, not defaults.

Distance and similarity

Similarity may be Euclidean or Mahalanobis, in the original variables or in PCA / PLS-DA score spaces (Section 2.2). The paper does not define LW-PLS-DA as Euclidean-only. Its reported experiments use Euclidean distance. Public V1 therefore uses original-space Euclidean distance (Eq. 1). Other geometries are source-allowed options, not implemented here.

d_{ik} = j = 1 \sum NV (x_{ij} - x_{k j})^{2} k \neq = i

NV is the number of predictors. Euclidean distance depends on variable scaling. Local pretreatment in this algorithm is applied after neighbor selection, so the V1 neighbor geometry is the stored calibration X, not the locally autoscaled X.

Distance-based weighting

Neighbor selection decides which rows enter. Table 1 then decides how strongly those rows scale the local predictor matrix. Before Table 1, local distances are mapped to [0, 1] by dividing by the farthest selected neighbor (Eq. 5). The Table 1 formulas were transcribed from the rendered PDF, not from OCR.

d_{m a x} = k = 1, \dots, n_{l} max d_{ik}

Scheme	$W_{k k}$	Note
Uniform	$1$	Locality from neighbors only
Triangular	$1 - d_{ik}$	0 at the farthest neighbor
Quadratic	$(1 - d_{ik}^{2})^{2}$	Bisquare on normalized d
Cubic	$(1 - d_{ik}^{3})^{3}$	Tricube on normalized d
Gaussian	$e^{- d_{ik}^{2}}$	Not zero at d = 1
Exponential	$e^{- ∣ d_{ik} ∣}$	Not zero at d = 1
Inverse	$1/ d_{ik}$	Singular at d = 0
Inverse quadratic	$1/ d_{ik}^{2}$	Singular at d = 0
Inverse shifted	$1/ (d_{ik} + 1)$	Finite at d = 0

The diagonal entries of W are predictor-row weights, not class probabilities. Section 2.3 says W elements are scaled between 0 and 1; in context that describes the normalized distances used in Table 1 and Fig. 1. Inverse and inverse quadratic of Table 1 can exceed 1. This implementation does not rewrite them as scikit-learn kernels and does not max-normalize W after Table 1. Different datasets in the paper favored different schemes, including Uniform. No scheme is universally best.

Uniform versus weighted

Uniform sets every selected neighbor to weight 1. Locality then comes only from who was selected. That is still local PLS-DA, not k-NN. Nonuniform Table 1 schemes further scale closer rows more strongly than farther selected rows, after d_max normalization.

Local PLS-DA construction

Order from Section 2.2: select neighbors; form W from Table 1; build $X_{w} = W X_{l}$ ; preprocess that local training matrix; apply the same X pretreatment to the query; fit PLS-DA if two or more classes are present. The query is not one of the weighted training rows. Y_l is not multiplied by W. This is not generic weighted least squares with sqrt(W) on both blocks.

Mathematics

Paper Eq. (2) is printed as $X_{w} = W X_{i}$ . Step c.4.1 and Section 2.3 name the same local predictor matrix $X_{l}$ . The public implementation uses

X_{w} = W X_{l}

with W diagonal of size n_l. Classification error for a held-out object is

E^{i} (n_{l}, A) = {10 c_{i, pred} \neq = c_{i, true} c_{i, pred} = c_{i, true}

and the leave-one-out percent error is

E_{CV} (n_{l}, A) = \frac{100}{N} i = 1 \sum N E^{i} (n_{l}, A)

Choose n_local and A_opt at a minimum of that surface. If several pairs share the exact minimum, the paper does not specify a breaker. OPEN LAB IMPLEMENTATION POLICY: smallest A, then smallest n_l.

Algorithm

Algorithm 1

LW-PLS-DA

InputCalibration X, y; query x_q; n_local; A_opt; Table 1 scheme; local X pretreatment

OutputPredicted class of x_q

01Require calibration $X, y$ , query $x_{q}$ , $n_{local}$ , $A_{opt}$ , Table 1 scheme, local pretreatment
02Encode global dummy $Y$ (one 1 per row, zeros elsewhere).
03Compute Euclidean distances (Eq. 1) from $x_{q}$ to calibration rows. Exclude the query from its own neighborhood.
04Keep the $n_{local}$ nearest rows as $X_{l}, Y_{l}$ .
05If one class only: assign that class and stop (no PLS-DA).
06Normalize distances by $d_{m a x}$ (Eq. 5). Build diagonal $W$ from Table 1.
07Form $X_{w} = W X_{l}$ (Eq. 2). Do not multiply $Y_{l}$ by $W$ .
08Fit local X pretreatment on $X_{w}$ and apply it to the query.
09Fit PLS-DA with $A_{opt}$ components. Assign $argmax$ of $y_{q}$ .
10Repeat independently for every new query.

Source steps: Bevilacqua and Marini 2014, Sections 2.1-2.3. Local n_local and A_opt are selected by the paper's LOO procedure when tuning is requested.

Multiclass classification

The method is explicitly multiclass. Global class columns are kept even if a neighborhood omits a class. Assignment is argmax of the G predicted dummy components. No one-versus-rest, pairwise voting, or midpoint threshold is used.

One-class local neighborhoods

If every selected neighbor belongs to one class, Y_l would be singular. The paper then assigns the query (or left-out object) directly to that class and does not fit PLS-DA. In LOO that 0/1 error is reused for every candidate A at that n_l.

Choosing the number of neighbors

Small n_local is more local. Large n_local uses more of the database. The paper tunes n_local jointly with A. Experimental grids are dataset-specific. Weighting-scheme choice, if data-driven, is another training-only hyperparameter.

Choosing latent variables

A_max is the largest local complexity examined. A_opt is taken with n_local at minimum E_CV. A cannot exceed the feasible local PLS rank. Impossible A values are skipped rather than factorized as invalid matrices.

Validation

SOURCE METHOD: for each training object i, leave i out, compute distances to the remaining N-1 objects, and never return i to its own neighborhood. Sweep n_l and A, accumulate Eqs. 3–4, then refit prediction models with the selected pair.

MODERN VALIDATION SAFEGUARD (Open Lab, not a 2014 quotation): any learned representation used for distance must be fitted inside the training fold. External-test labels must not choose n_local, A, pretreatment, or the Table 1 scheme. Use classification metrics from the Classification Metrics card; this page does not redefine them. The paper's reported accuracies are dataset-specific.

Prediction of a new sample

Repeat the local procedure for each query using the stored calibration database and the selected (or user-fixed) n_local, A, scheme, and pretreatment. Do not retune on the external test set.

LW-PLS-DA versus global PLS-DA

Global PLS-DA fits one discriminant model on all calibration samples. LW-PLS-DA builds a query-specific local PLS-DA. On the paper's four examples the local models were more accurate, including cases where global PLS-DA was near chance. That is not a universal superiority claim.

LW-PLS-DA versus k-NN

k-NN assigns a class from neighbor labels, typically by a vote. LW-PLS-DA uses neighbors as a local calibration set for PLS-DA. Using nearest neighbors does not make the methods the same. Kernel PLS-DA is a different nonlinear construction (one kernel latent model) and is not derived here.

Python

Fixed mode supplies n_neighbors and n_components. Tuned mode supplies neighbor_grid and n_components_max and runs the paper LOO search inside fit on training data only. Local PLS-DA after Eq. 2 uses the Open Lab educational NIPALS PLS2 engine, not observation-weighted rchemo plskern and not a midpoint classifier.

lw-pls-da

import numpy as np  TABLE1_SCHEMES = (    "uniform",    "triangular",    "quadratic",    "cubic",    "gaussian",    "exponential",    "inverse",    "inverse_quadratic",    "inverse_shifted",)  def encode_dummy(y):    y = np.asarray(y)    if y.ndim != 1:        raise ValueError("y must be a 1D array of class labels.")    if y.size == 0:        raise ValueError("y must contain at least one observation.")    classes = np.unique(y)    Y = np.zeros((y.shape[0], classes.shape[0]), dtype=float)    for g, lab in enumerate(classes):        Y[y == lab, g] = 1.0    return Y, classes  def assign_argmax(Yhat, classes):    Yhat = np.asarray(Yhat, dtype=float)    if Yhat.ndim == 1:        Yhat = Yhat.reshape(1, -1)    if Yhat.shape[1] != classes.shape[0]:        raise ValueError("Yhat must have one column per class.")    if not np.all(np.isfinite(Yhat)):        raise ValueError("Yhat must contain only finite values.")    # OPEN LAB IMPLEMENTATION DETAIL: exact ties keep the first class    # in np.unique order (numpy argmax).    return classes[np.argmax(Yhat, axis=1)]  def euclidean_distances_to_query(X, z):    X = np.asarray(X, dtype=float)    z = np.asarray(z, dtype=float).reshape(-1)    if z.shape[0] != X.shape[1]:        raise ValueError("Query z must have one value per predictor.")    if not np.all(np.isfinite(X)) or not np.all(np.isfinite(z)):        raise ValueError("X and z must contain only finite values.")    return np.sqrt(np.sum((X - z) ** 2, axis=1))  def select_local_neighbors(d, n_local):    d = np.asarray(d, dtype=float).reshape(-1)    n_local = int(n_local)    n = d.shape[0]    if n_local < 1 or n_local > n:        raise ValueError("n_local must be an integer between 1 and n inclusive.")    # OPEN LAB IMPLEMENTATION DETAIL: equal distances keep the smaller    # original index.    order = np.lexsort((np.arange(n), d))    idx = order[:n_local]    return idx, d[idx]  def normalize_local_distances(d_local):    d_local = np.asarray(d_local, dtype=float).reshape(-1)    if d_local.size == 0:        raise ValueError("Local distance vector must be non-empty.")    d_max = float(np.max(d_local))    if d_max == 0.0:        # OPEN LAB IMPLEMENTATION DETAIL: all selected neighbors are        # predictor-identical to the query. Eq. 5 is undefined.        return np.zeros_like(d_local), 0.0    return d_local / d_max, d_max  def table1_weights(d_norm, scheme):    d = np.asarray(d_norm, dtype=float).reshape(-1)    if scheme not in TABLE1_SCHEMES:        raise ValueError("Unknown Table 1 weighting scheme.")    if scheme == "uniform":        return np.ones_like(d)    if scheme == "triangular":        return 1.0 - d    if scheme == "quadratic":        return (1.0 - d ** 2) ** 2    if scheme == "cubic":        return (1.0 - d ** 3) ** 3    if scheme == "gaussian":        return np.exp(-(d ** 2))    if scheme == "exponential":        return np.exp(-np.abs(d))    if scheme == "inverse_shifted":        return 1.0 / (d + 1.0)    if scheme in ("inverse", "inverse_quadratic"):        w = np.empty_like(d)        pos = d > 0.0        if scheme == "inverse":            w[pos] = 1.0 / d[pos]        else:            w[pos] = 1.0 / (d[pos] ** 2)        zero = ~pos        if np.any(zero):            # OPEN LAB IMPLEMENTATION DETAIL: Table 1 Inverse and            # Inverse quadratic are singular at d = 0. Identical            # neighbors receive the largest finite weight from the            # same scheme in this neighborhood, or 1 if d_max = 0.            if np.any(pos):                w[zero] = np.max(w[pos])            else:                w[zero] = 1.0        return w    raise ValueError("Unknown Table 1 weighting scheme.")  def build_weighted_local_X(X_l, w):    X_l = np.asarray(X_l, dtype=float)    w = np.asarray(w, dtype=float).reshape(-1)    if w.shape[0] != X_l.shape[0]:        raise ValueError("Weight vector length must match n_local.")    return w[:, None] * X_l  def fit_local_preprocessing(Xw, method):    Xw = np.asarray(Xw, dtype=float)    method = str(method)    if method == "none":        mu = np.zeros(Xw.shape[1])        sigma = np.ones(Xw.shape[1])        return Xw.copy(), mu, sigma    mu = Xw.mean(axis=0)    if method == "center":        return Xw - mu, mu, np.ones(Xw.shape[1])    if method == "autoscale":        if Xw.shape[0] < 2:            raise ValueError("Autoscaling requires at least two local samples.")        sigma = Xw.std(axis=0, ddof=1)        # OPEN LAB IMPLEMENTATION DETAIL: zero-variance columns keep scale 1.        sigma = np.where(sigma == 0.0, 1.0, sigma)        return (Xw - mu) / sigma, mu, sigma    raise ValueError("preprocessing must be none, center, or autoscale.")  def apply_local_preprocessing(X, mu, sigma):    return (np.asarray(X, dtype=float) - mu) / sigma  def nipals_pls(X, Y, n_components, max_iter=500, tol=1e-6):    X = np.asarray(X, dtype=float)    Y = np.asarray(Y, dtype=float)    if Y.ndim == 1:        Y = Y.reshape(-1, 1)        one_d = True    else:        one_d = False    n_components = int(n_components)    if n_components < 1:        raise ValueError("n_components must be a positive integer.")    if n_components > min(X.shape):        raise ValueError("n_components cannot exceed min(n_local, n_features).")    n_samples, n_features = X.shape    n_targets = Y.shape[1]    x_mean = X.mean(axis=0)    y_mean = Y.mean(axis=0)    Xk = X - x_mean    Yk = Y - y_mean    eps = np.finfo(float).eps    x_weights = np.zeros((n_features, n_components))    x_scores = np.zeros((n_samples, n_components))    x_loadings = np.zeros((n_features, n_components))    y_loadings = np.zeros((n_targets, n_components))    for a in range(n_components):        if not np.any(np.abs(Xk) > eps):            raise ValueError("Local predictor residual has no variation.")        try:            u = next(col.copy() for col in Yk.T if np.any(np.abs(col) > eps))        except StopIteration:            raise ValueError("Y residual is constant.")        w_prev = None        for it in range(max_iter):            w = (Xk.T @ u) / (u @ u)            w = w / (np.sqrt(w @ w) + eps)            t = Xk @ w            c = (Yk.T @ t) / (t @ t)            u = (Yk @ c) / ((c @ c) + eps)            if w_prev is not None and (w - w_prev) @ (w - w_prev) < tol:                break            if n_targets == 1:                break            w_prev = w.copy()        t = Xk @ w        p_vec = (Xk.T @ t) / (t @ t)        q_vec = (Yk.T @ t) / (t @ t)        Xk = Xk - np.outer(t, p_vec)        Yk = Yk - np.outer(t, q_vec)        x_weights[:, a] = w        x_scores[:, a] = t        x_loadings[:, a] = p_vec        y_loadings[:, a] = q_vec    try:        x_rotations = x_weights @ np.linalg.pinv(x_loadings.T @ x_weights)    except np.linalg.LinAlgError as exc:        raise ValueError("Local PLS-DA factorization is infeasible.") from exc    if not np.all(np.isfinite(x_rotations)):        raise ValueError("Local PLS-DA factorization is infeasible.")    coef = (x_rotations @ y_loadings.T).T    return {        "coef": coef,        "intercept": y_mean.copy(),        "x_mean": x_mean,        "one_d": one_d,    }  def predict_nipals_pls(model, X_new):    X_new = np.asarray(X_new, dtype=float)    y_hat = (X_new - model["x_mean"]) @ model["coef"].T + model["intercept"]    if model["one_d"]:        return y_hat.ravel()    return y_hat  def classification_error(pred, true):    return 0 if pred == true else 1  def ecv_percent(errors):    errors = np.asarray(errors, dtype=float)    n = errors.size    if n == 0:        raise ValueError("Error vector must be non-empty.")    return 100.0 * np.sum(errors) / n  def loo_candidate_mask(n, i):    mask = np.ones(n, dtype=bool)    mask[i] = False    return mask  class LWPLSDA:    def __init__(        self,        n_neighbors=None,        n_components=None,        weighting="uniform",        preprocessing="center",        neighbor_grid=None,        n_components_max=None,    ):        self.n_neighbors = None if n_neighbors is None else int(n_neighbors)        self.n_components = None if n_components is None else int(n_components)        self.weighting = weighting        self.preprocessing = preprocessing        self.neighbor_grid = None if neighbor_grid is None else np.asarray(neighbor_grid, dtype=int)        self.n_components_max = None if n_components_max is None else int(n_components_max)        if self.weighting not in TABLE1_SCHEMES:            raise ValueError("weighting must be a Table 1 scheme name.")        if self.preprocessing not in ("none", "center", "autoscale"):            raise ValueError("preprocessing must be none, center, or autoscale.")     def fit(self, X, y):        X = np.asarray(X, dtype=float)        y = np.asarray(y)        if X.ndim != 2:            raise ValueError("X must be a 2D array of samples by variables.")        if X.shape[0] != y.shape[0]:            raise ValueError("X and y must have the same number of observations.")        if not np.all(np.isfinite(X)):            raise ValueError("X must contain only finite values.")        self.X_ = X        self.y_ = y        self.Y_, self.classes_ = encode_dummy(y)        tuned = self.neighbor_grid is not None and self.n_components_max is not None        fixed = self.n_neighbors is not None and self.n_components is not None        if tuned == fixed:            raise ValueError("Supply either fixed n_neighbors and n_components, or neighbor_grid and n_components_max.")        if tuned:            self.ecv_, self.n_neighbors, self.n_components = cross_validate_local_complexity(                X,                y,                self.neighbor_grid,                self.n_components_max,                self.weighting,                self.preprocessing,            )        return self     def query_local(self, z):        z = np.asarray(z, dtype=float).reshape(-1)        d = euclidean_distances_to_query(self.X_, z)        idx, dsel = select_local_neighbors(d, self.n_neighbors)        y_loc = self.y_[idx]        local_classes = np.unique(y_loc)        if local_classes.size == 1:            return {                "indices": idx,                "distances": dsel,                "d_max": float(np.max(dsel)),                "weights": None,                "one_class": True,                "yhat": None,                "pred": local_classes[0],            }        d_norm, d_max = normalize_local_distances(dsel)        if d_max == 0.0:            w = np.ones(idx.shape[0])        else:            w = table1_weights(d_norm, self.weighting)        X_w = build_weighted_local_X(self.X_[idx], w)        X_proc, mu, sigma = fit_local_preprocessing(X_w, self.preprocessing)        z_proc = apply_local_preprocessing(z.reshape(1, -1), mu, sigma)        Y_l = self.Y_[idx]        local = nipals_pls(X_proc, Y_l, self.n_components)        yhat = predict_nipals_pls(local, z_proc)        if yhat.ndim == 1:            yhat = yhat.reshape(1, -1)        pred = assign_argmax(yhat, self.classes_)[0]        return {            "indices": idx,            "distances": dsel,            "d_max": d_max,            "weights": w,            "one_class": False,            "yhat": yhat[0],            "pred": pred,        }     def predict(self, X_new):        X_new = np.asarray(X_new)        if X_new.ndim != 2:            raise ValueError("X_new must be a 2D array of samples by variables.")        return np.array([self.query_local(row)["pred"] for row in X_new])  def predict_one_held_out(X_train, y_train, Y_train, classes, z, y_true, n_local, A, weighting, preprocessing):    d = euclidean_distances_to_query(X_train, z)    idx, dsel = select_local_neighbors(d, n_local)    y_loc = y_train[idx]    local_classes = np.unique(y_loc)    if local_classes.size == 1:        pred = local_classes[0]        return classification_error(pred, y_true), pred    if A > min(idx.shape[0], X_train.shape[1]):        return 1, None    d_norm, d_max = normalize_local_distances(dsel)    w = np.ones(idx.shape[0]) if d_max == 0.0 else table1_weights(d_norm, weighting)    X_w = build_weighted_local_X(X_train[idx], w)    try:        X_proc, mu, sigma = fit_local_preprocessing(X_w, preprocessing)        z_proc = apply_local_preprocessing(np.asarray(z, dtype=float).reshape(1, -1), mu, sigma)        local = nipals_pls(X_proc, Y_train[idx], A)        yhat = predict_nipals_pls(local, z_proc)        if yhat.ndim == 1:            yhat = yhat.reshape(1, -1)        pred = assign_argmax(yhat, classes)[0]    except ValueError:        return 1, None    return classification_error(pred, y_true), pred  def cross_validate_local_complexity(X, y, neighbor_grid, n_components_max, weighting="uniform", preprocessing="center"):    X = np.asarray(X, dtype=float)    y = np.asarray(y)    Y, classes = encode_dummy(y)    n = X.shape[0]    neighbor_grid = np.asarray(neighbor_grid, dtype=int)    n_components_max = int(n_components_max)    A_grid = np.arange(1, n_components_max + 1)    ecv = np.full((neighbor_grid.size, A_grid.size), np.nan)    for a_idx, A in enumerate(A_grid):        for n_idx, n_l in enumerate(neighbor_grid):            errors = []            for i in range(n):                mask = loo_candidate_mask(n, i)                err, _ = predict_one_held_out(                    X[mask],                    y[mask],                    Y[mask],                    classes,                    X[i],                    y[i],                    int(n_l),                    int(A),                    weighting,                    preprocessing,                )                errors.append(err)            ecv[n_idx, a_idx] = ecv_percent(errors)    # OPEN LAB IMPLEMENTATION POLICY: if several (n_l, A) share the    # minimum ECV, keep the smallest A, then the smallest n_l.    finite = np.isfinite(ecv)    if not np.any(finite):        raise ValueError("No feasible (n_local, A) combination was found.")    best = np.min(ecv[finite])    candidates = np.argwhere(np.isclose(ecv, best) & finite)    A_vals = A_grid[candidates[:, 1]]    n_vals = neighbor_grid[candidates[:, 0]]    order = np.lexsort((n_vals, A_vals))    n_star = int(n_vals[order[0]])    A_star = int(A_vals[order[0]])    return ecv, n_star, A_star

MATLAB

MATLAB has no built-in LW-PLS-DA. The listing is the same 2014 algorithm. Do not pass row weights into plsregress and call the result Bevilacqua LW-PLS-DA. ORIGINAL MATLAB IMPLEMENTATION REFERENCED BY PAPER NOT RETRIEVED (Rome Chemometrics download last accessed May 2014 in the paper).

lw-pls-da

function model = fit_lwplsda(X, y, nNeighbors, nComponents, weighting, preprocessing)    if nargin < 5, weighting = "uniform"; end    if nargin < 6, preprocessing = "center"; end    y = y(:);    [Y, classes] = encodeDummy_lwplsda(y);    model.X = X;    model.y = y;    model.Y = Y;    model.classes = classes;    model.nNeighbors = double(nNeighbors);    model.nComponents = double(nComponents);    model.weighting = char(weighting);    model.preprocessing = char(preprocessing);end function pred = predict_lwplsda(model, Xnew)    m = size(Xnew, 1);    pred = repmat(model.classes(1), m, 1);    for i = 1:m        q = query_lwplsda(model, Xnew(i, :));        pred(i) = q.pred;    endend function q = query_lwplsda(model, z)    z = z(:)';    d = euclideanDistances_lwplsda(model.X, z);    [idx, dsel] = selectLocalNeighbors_lwplsda(d, model.nNeighbors);    yLoc = model.y(idx);    localClasses = unique(yLoc);    if numel(localClasses) == 1        q.indices = idx;        q.distances = dsel;        q.d_max = max(dsel);        q.weights = [];        q.one_class = true;        q.yhat = [];        q.pred = localClasses(1);        return    end    [dNorm, dMax] = normalizeLocalDistances_lwplsda(dsel);    if dMax == 0        w = ones(numel(idx), 1);    else        w = table1Weights_lwplsda(dNorm, model.weighting);    end    Xw = buildWeightedLocalX_lwplsda(model.X(idx, :), w);    [Xproc, mu, sigma] = fitLocalPreprocessing_lwplsda(Xw, model.preprocessing);    zproc = applyLocalPreprocessing_lwplsda(z, mu, sigma);    local = nipalsPls_lwplsda(Xproc, model.Y(idx, :), model.nComponents);    yhat = predictNipals_lwplsda(local, zproc);    q.indices = idx;    q.distances = dsel;    q.d_max = dMax;    q.weights = w;    q.one_class = false;    q.yhat = yhat;    q.pred = assignArgmax_lwplsda(yhat, model.classes);end function [Y, classes] = encodeDummy_lwplsda(y)    y = y(:);    classes = unique(y);    Y = zeros(numel(y), numel(classes));    for g = 1:numel(classes)        Y(:, g) = double(y == classes(g));    endend function pred = assignArgmax_lwplsda(Yhat, classes)    Yhat = Yhat(:)';    [~, idx] = max(Yhat);    pred = classes(idx);end function d = euclideanDistances_lwplsda(X, z)    d = sqrt(sum((X - z).^2, 2));end function [idx, dsel] = selectLocalNeighbors_lwplsda(d, nLocal)    n = numel(d);    [~, order] = sortrows([d(:), (1:n)']);    idx = order(1:nLocal);    dsel = d(idx);end function [dNorm, dMax] = normalizeLocalDistances_lwplsda(dLocal)    dLocal = dLocal(:);    dMax = max(dLocal);    if dMax == 0        dNorm = zeros(size(dLocal));    else        dNorm = dLocal / dMax;    endend function w = table1Weights_lwplsda(dNorm, scheme)    d = dNorm(:);    switch char(scheme)        case "uniform"            w = ones(size(d));        case "triangular"            w = 1 - d;        case "quadratic"            w = (1 - d.^2).^2;        case "cubic"            w = (1 - d.^3).^3;        case "gaussian"            w = exp(-(d.^2));        case "exponential"            w = exp(-abs(d));        case "inverse_shifted"            w = 1 ./ (d + 1);        case {"inverse", "inverse_quadratic"}            w = zeros(size(d));            pos = d > 0;            if strcmp(char(scheme), "inverse")                w(pos) = 1 ./ d(pos);            else                w(pos) = 1 ./ (d(pos).^2);            end            zero = ~pos;            if any(zero)                if any(pos)                    w(zero) = max(w(pos));                else                    w(zero) = 1;                end            end        otherwise            error("Unknown Table 1 weighting scheme.");    endend function Xw = buildWeightedLocalX_lwplsda(Xl, w)    Xw = w(:) .* Xl;end function [Xproc, mu, sigma] = fitLocalPreprocessing_lwplsda(Xw, method)    method = char(method);    p = size(Xw, 2);    if strcmp(method, "none")        mu = zeros(1, p);        sigma = ones(1, p);        Xproc = Xw;        return    end    mu = mean(Xw, 1);    if strcmp(method, "center")        sigma = ones(1, p);        Xproc = Xw - mu;        return    end    if strcmp(method, "autoscale")        sigma = std(Xw, 0, 1);        sigma(sigma == 0) = 1;        Xproc = (Xw - mu) ./ sigma;        return    end    error("preprocessing must be none, center, or autoscale.");end function Xout = applyLocalPreprocessing_lwplsda(X, mu, sigma)    Xout = (X - mu) ./ sigma;end function model = nipalsPls_lwplsda(X, Y, nComponents)    if isvector(Y)        Y = Y(:);    end    nComponents = double(nComponents);    if nComponents > min(size(X))        error("nComponents cannot exceed min(n_local, n_features).");    end    [nSamples, nFeatures] = size(X);    nTargets = size(Y, 2);    xMean = mean(X, 1);    yMean = mean(Y, 1);    Xk = X - xMean;    Yk = Y - yMean;    epsX = eps;    W = zeros(nFeatures, nComponents);    T = zeros(nSamples, nComponents);    P = zeros(nFeatures, nComponents);    Q = zeros(nTargets, nComponents);    maxIter = 500;    tol = 1e-6;    for a = 1:nComponents        if ~any(abs(Xk(:)) > epsX)            error("Local predictor residual has no variation.");        end        u = [];        for j = 1:nTargets            col = Yk(:, j);            if any(abs(col) > epsX)                u = col;                break            end        end        if isempty(u)            error("Y residual is constant.");        end        wPrev = [];        for it = 1:maxIter            w = (Xk' * u) / (u' * u);            w = w / (norm(w) + epsX);            t = Xk * w;            c = (Yk' * t) / (t' * t);            u = (Yk * c) / ((c' * c) + epsX);            if ~isempty(wPrev)                if (w - wPrev)' * (w - wPrev) < tol                    break                end            end            if nTargets == 1                break            end            wPrev = w;        end        t = Xk * w;        pVec = (Xk' * t) / (t' * t);        qVec = (Yk' * t) / (t' * t);        Xk = Xk - t * pVec';        Yk = Yk - t * qVec';        W(:, a) = w;        T(:, a) = t;        P(:, a) = pVec;        Q(:, a) = qVec;    end    R = W * pinv(P' * W);    model.coef = (R * Q')';    model.intercept = yMean;    model.x_mean = xMean;end function yhat = predictNipals_lwplsda(model, Xnew)    yhat = (Xnew - model.x_mean) * model.coef' + model.intercept;end

Practical notes

This is Bevilacqua and Marini 2014, not Lesnoff/rchemo KNN-LWPLSR.
Neighbor selection and Table 1 row weighting are different operations.
V1 distance is original-space Euclidean. Other spaces are named in the paper and not implemented here.
Weight X_l first, then fit local X pretreatment on X_w.
Predicted dummy components are not probabilities. Assign argmax.
n_local and A have no universal defaults. Tune them without external-test labels.
A one-class neighborhood is classified by that class, without PLS-DA.
Inverse and inverse quadratic are singular at zero distance; the safeguard is labelled as an Open Lab detail.
Query-specific models are heavier than one global PLS-DA.
Use Classification Metrics with the evaluation context stated.

References

1.
Bevilacqua, M., & Marini, F. (2014). Local classification: Locally weighted-partial least squares-discriminant analysis (LW-PLS-DA). Analytica Chimica Acta, 838, 20-30.
doi:10.1016/j.aca.2014.05.057
2.
Centner, V., & Massart, D. L. (1998). Optimization in locally weighted regression. Analytical Chemistry, 70(19), 4206-4211.
doi:10.1021/ac980208r
3.
Atkeson, C. G., Moore, A. W., & Schaal, S. (1997). Locally weighted learning. Artificial Intelligence Review, 11, 11-73.
doi:10.1023/A:1006559212014
4.
Barker, M., & Rayens, W. (2003). Partial least squares for discrimination. Journal of Chemometrics, 17(3), 166-173.
doi:10.1002/cem.785
5.
Stahle, L., & Wold, S. (1987). Partial least squares analysis with cross-validation for the two-class problem: a Monte Carlo study. Journal of Chemometrics, 1(3), 185-196.
doi:10.1002/cem.1180010306
6.
Cover, T. M., & Hart, P. E. (1967). Nearest neighbor pattern classification. IEEE Transactions on Information Theory, 13(1), 21-27.
doi:10.1109/TIT.1967.1053964

← Back to Open Lab

What is LW-PLS-DA?

Why local classification?

Query-specific modeling

PLS-DA response coding

Finding the local neighborhood

Distance and similarity

Distance-based weighting

Uniform versus weighted

Local PLS-DA construction

Mathematics

Algorithm

Multiclass classification

One-class local neighborhoods

Choosing the number of neighbors

Choosing latent variables

Validation

Prediction of a new sample

LW-PLS-DA versus global PLS-DA

LW-PLS-DA versus k-NN

Python

MATLAB

Practical notes

References

Related methods