NIPALS

NIPALS (Nonlinear Iterative Partial Least Squares) is an iterative computational procedure for extracting latent components in two-block PLS models. Geladi and Kowalski present an algorithm for a predictive PLS in that sense. The name refers to the algorithm, not to the PLS modeling framework itself. The inner- loop algebra on this page follows scikit-learn 1.6.1 _pls.py, not a transcription of the Geladi and Kowalski listing.

Partial Least Squares Regression is the two-block calibration method. NIPALS is one way to compute a PLS solution. SIMPLS is another. The Open Lab PLSR page treats the method. This page treats one NIPALS formulation in detail.

Components are extracted sequentially. Each component is obtained by iterating weight and score updates on the current residual matrices. After the inner iteration converges, those residuals are deflated and the next component is computed from what remains.

At component $a$, the inner loop estimates the first left and right singular vectors of the current cross-covariance $X_a^{\mathsf{T}}{Y}_a$. Those vectors are the X-weight $w$ and the Y-weight $c$. They are chosen so that the covariance between the projected blocks, $\mathrm{Cov}(X_{a}w,Y_{a}c)$, is maximized subject to the normalization of this formulation: $\Vert w\Vert =1$, while $c$ is not unit-normalized.

A full SVD of $X_a^{\mathsf{T}}{Y}_a$ would return those same leading singular vectors. The NIPALS option replaces that SVD by a power method: alternate between the two blocks until the X-weight stabilizes. Iteration is required because each block update uses the current score from the other block. The loop is the numerical search for those singular vectors, not a search for an undefined "best component".

Column-center $X$ and $Y$. The educational implementation does not scale columns. scikit-learn PLSRegression scales both blocks by default.

For each component, initialize the Y-score $u$ as the first non-constant column of the current $Y$ residual. Then iterate: regress the current $X$ residual on $u$ to obtain $w$, normalize $w$ to unit Euclidean length, form the X-score $t=X_{a}w$, regress the current $Y$ residual on $t$ to obtain $c$, and update $u=Y_{a}c/(c^{\mathsf{T}}c)$.

Stop when the squared Euclidean change in successive X-weights is below a tolerance, or after one pass if $Y$ has a single column. Then compute loadings by regressing the current $X$ and $Y$ residuals on $t$, and subtract those rank-one fits. The next component uses the deflated matrices.

This is PLS regression NIPALS: Mode A weights and regression-mode Y deflation on X scores. It covers PLS2 in general. PLS1 is the single-response case of the same procedure. It is not PLS canonical deflation, which residualizes $Y$ on Y scores, and it is not SIMPLS.

Samples are rows. Predictors occupy the columns of $X\in {\mathbb{R}}^{n\times p}$. Responses occupy the columns of $Y\in {\mathbb{R}}^{n\times m}$. Notation is held fixed below. scikit-learn documents the X-weight as the left singular vector $u_k$. That library symbol is this page's $w$, not the Y-score $u$.

Algorithm 1 is the two-block PLS regression NIPALS used by the educational code: Mode A, Y-weights not unit-normalized, regression deflation, column centering, no scaling.

PLS1 is this algorithm with a single response column. PLS2 is the same algorithm with two or more response columns. A scalar-y shortcut is not a different model class here. It only removes the need to iterate the inner loop.

NIPALS and SIMPLS both compute PLS regression models. They are not the same algorithm. NIPALS builds successive residual matrices. SIMPLS forms factors from the original variables and avoids that NIPALS-style deflation. MATLAB plsregress uses SIMPLS. Current scikit-learn PLSRegression uses this NIPALS power-method route internally.

The inner power method estimates the leading singular vectors of $X_a^{\mathsf{T}}{Y}_a$. That is related to SVD, and scikit-learn also offers an SVD solver for the same inner vectors. NIPALS is not a synonym for SVD. A single SVD of the original $X^{\mathsf{T}}Y$ without sequential deflation is a different estimator (PLSSVD).

The Python listing implements Algorithm 1 with NumPy. It does not call scikit-learn inside the educational fit. The sklearn tab is the current library estimator, which uses the same Mode A regression-NIPALS route and, by default, also scales both blocks. Compare numerical intermediates only with scale=False. Latent vectors are defined up to sign.

MATLAB plsregress is SIMPLS. The MATLAB tab is an explicit NIPALS implementation of the same formulation. Do not label plsregress output as NIPALS.

import numpy as np  def _validate_xy(X, Y):    X = np.asarray(X, dtype=float)    Y = np.asarray(Y, dtype=float)    if X.ndim != 2:        raise ValueError("X must be a 2D array of samples by variables.")    if Y.ndim == 1:        Y = Y.reshape(-1, 1)        one_d = True    elif Y.ndim == 2:        one_d = False    else:        raise ValueError("Y must be a 1D vector or a 2D array of samples by responses.")    if X.shape[0] != Y.shape[0]:        raise ValueError("X and Y must have the same number of observations.")    if min(X.shape) == 0 or Y.shape[1] < 1:        raise ValueError("X and Y must have at least one row and one column.")    if not np.all(np.isfinite(X)) or not np.all(np.isfinite(Y)):        raise ValueError("X and Y must contain only finite values.")    return X, Y, one_d  def fit_nipals_pls(X, Y, n_components, tol=1e-6, max_iter=500):    X, Y, one_d = _validate_xy(X, Y)    n_components = int(n_components)    max_iter = int(max_iter)    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_samples, n_features).")    if max_iter < 1:        raise ValueError("max_iter must be a positive integer.")    if tol < 0:        raise ValueError("tol must be nonnegative.")     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))    y_weights = np.zeros((n_targets, n_components))    x_scores = np.zeros((n_samples, n_components))    y_scores = np.zeros((n_samples, n_components))    x_loadings = np.zeros((n_features, n_components))    y_loadings = np.zeros((n_targets, n_components))    n_iter = []     for a in range(n_components):        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()        n_iter.append(it + 1)         t = Xk @ w        u_score = (Yk @ c) / (c @ c)        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        y_weights[:, a] = c        x_scores[:, a] = t        y_scores[:, a] = u_score        x_loadings[:, a] = p_vec        y_loadings[:, a] = q_vec     x_rotations = x_weights @ np.linalg.pinv(x_loadings.T @ x_weights)    coef = (x_rotations @ y_loadings.T).T    return {        "x_weights": x_weights,        "y_weights": y_weights,        "x_scores": x_scores,        "y_scores": y_scores,        "x_loadings": x_loadings,        "y_loadings": y_loadings,        "x_rotations": x_rotations,        "coef": coef,        "intercept": y_mean.copy(),        "x_mean": x_mean,        "y_mean": y_mean,        "n_iter": n_iter,        "one_d": one_d,    }  def predict_nipals_pls(model, X_new):    X_new = np.asarray(X_new, dtype=float)    if X_new.ndim != 2:        raise ValueError("X_new must be a 2D array of samples by variables.")    if not np.all(np.isfinite(X_new)):        raise ValueError("X_new must contain only finite values.")    y_hat = (X_new - model["x_mean"]) @ model["coef"].T + model["intercept"]    if model["one_d"]:        return y_hat.ravel()    return y_hat 

1. 1.

   Geladi, P., & Kowalski, B. R. (1986). Partial least-squares regression: a tutorial. Analytica Chimica Acta, 185, 1-17.

   doi:10.1016/0003-2670(86)80028-9
2. 2.

   Wegelin, J. A. (2000). A Survey of Partial Least Squares (PLS) Methods, with Emphasis on the Two-Block Case. University of Washington, Department of Statistics, Technical Report 371.

3. 3.

   de Jong, S. (1993). SIMPLS: an alternative approach to partial least squares regression. Chemometrics and Intelligent Laboratory Systems, 18(3), 251-263.

   doi:10.1016/0169-7439(93)85002-X
4. 4.

   scikit-learn Developers (1.6.1 / 1.9 user guide). PLSRegression, Cross decomposition, and sklearn.cross_decomposition._pls. scikit-learn documentation and source.

5. 5.

   The MathWorks, Inc. (n.d.). plsregress. MATLAB documentation. States that plsregress uses SIMPLS.