Open Lab/Multiblock Data Analysis · Method

Multiblock Partial Least Squares

MB-PLS

Supervised multiblock PLS with multiple predictor blocks and one response block, in the Westerhuis 1998 super-score-deflation formulation equivalent to standard PLS on a sqrt(J_b)-weighted superblock.

ChemometricsMultiblockPLSData FusionRegressionPythonMATLAB

MB-PLS

What is MB-PLS?

This card teaches one supervised multiblock method: Multiblock Partial Least Squares in the formulation analysed by Westerhuis, Kourti and MacGregor, with super-score deflation of the predictor blocks and of the response. Several predictor blocks share the same observations. One response block $Y$ guides the latent directions. The names MB-PLS, hierarchical PLS, SO-PLS, PO-PLS, OnPLS, and N-PLS are not synonyms.

The PRIMARY IMPLEMENTED FORMULATION is that Westerhuis 1998 super-score-deflation MB-PLS. After the selected block scaling, its global predictive solution is the Open Lab NIPALS PLS model of the weighted superblock. Block weights, block scores, super weights, and block loadings are then recovered from the original block structure. That is not the claim that every historical MB-PLS algorithm equals ordinary PLS.

Wangen and Kowalski 1989 is the historical origin of a general multiblock PLS construction for complex chemical systems. Their pathway algorithm, and later block-score deflation, are HISTORICAL FORMULATIONS. They are described, not implemented. Westerhuis and Smilde 2001 discuss a Y-only deflation alternative. That DEFLATION ALTERNATIVE is described from the 2001 abstract. Its regression-vector recursion is not implemented, because the 2001 body was not opened for this card.

The supervised multiblock problem

MB-PCA finds variation in predictor blocks without using a response. MB-PLS is supervised: latent directions are determined so that combined predictor information is related to $Y$ . It is not MB-PCA followed by regression, and it is not a set of independent PLS models, one per block.

Keeping blocks separate is a scientific choice. Different instruments, process units, or variable groups remain labelled after the global model is fit. The global prediction can still be a PLS model of a concatenated weighted matrix. The block labels are what make block-level weights, scores, and loadings interpretable.

SO-PLS is a different supervised multiblock method. In the selected MB-PLS formulation the predictor blocks enter jointly through the weighted superblock. SO-PLS treats predictor blocks sequentially with orthogonalization. They do not solve the same algebraic problem. SO-PLS equations are not derived here. OnPLS is a symmetric multi-block extension of O2PLS. Its use of the word predictive does not make it an ordinary external-Y replacement for this MB-PLS. N-PLS operates on a genuine multiway array. MB-PLS uses multiple linked matrices. Multiblock and multiway are not the same structure.

Predictor blocks and response block

For $B$ predictor blocks measured on the same $I$ observations in matching row order,

X_{b} \in R^{I \times J_{b}}, b = 1, \dots, B

The response is one block $Y \in R^{I \times M}$ . If $M = 1$ , the global model is PLS1. If $M > 1$ , it is PLS2 under the selected NIPALS convention. This V1 card does not develop multiple Y-block path models, L-shaped structures, or the general Wangen network of predictor and predicted blocks. Those are advanced extensions of the 1989 algorithm, not this implementation.

Code refuses silently mismatched rows. Blocks are not reordered or truncated to force alignment.

Preprocessing and block scaling

Within-block preprocessing and between-block scaling are different operations. Training column means of each $X_{b}$ are stored and reused. Optional autoscaling uses the training sample standard deviation with denominator $I - 1$ , matching the Normalization & Scaling card. Autoscaling is not prescribed for every application. A zero-variance column makes it undefined.

The NIPALS core column-centers $Y$ and stores the training response mean as the intercept. Educational Y is not autoscaled. Do not inherit a software package default that scales $Y$ unless that choice is stated.

Whole-block scaling changes the model. MB-PLS does not automatically balance all blocks. In the Westerhuis 1998 formulation analysed here, block scores divide by $m_{X_{b}}$ , where $m_{X_{b}}$ is their number of variables in block $b$ . This card writes that count as $J_{b}$ . The matching standard-PLS superblock therefore uses

α_{b} = 1/ J_{b}, Z_{b} = α_{b} X_{b}

That $1/ J_{b}$ rule belongs to the selected historical/equivalent formulation. It is not a universal optimum. The option none sets $α_{b} = 1$ and matches their discussion of MB-PLS without that block-score scaling. Other block weights exist. Changing them changes the effective objective.

From blocks to the predictor superblock

After within-block preprocessing and the chosen $α_{b}$ ,

Z = [Z_{1} ∣ Z_{2} ∣ \dots ∣ Z_{B}] \in R^{I \times J}, J = b = 1 \sum B J_{b}

Concatenation is the algebraic device for the global PLS model. Column slices that belong to each block are stored. The concatenated array is not a PARAFAC tensor, and it is not an excuse to forget which variables came from which block.

MB-PLS latent structure

Each component has a global X super score, a Y score, PLS weights and loadings on the superblock, and block-level weights, scores, and loadings. Y influences the latent directions through the PLS covariance criterion of the NIPALS card. The super score is a linear combination of block scores. It contains information from all predictor blocks under the selected scaling.

Weights and loadings are not interchangeable. A PLS weight is the direction used to form a score from the current residual matrix. A PLS loading is the regression of that residual matrix on the score used for deflation. The NIPALS card uses $w$ and $c$ for inner-loop weights and $p$ and $q$ for loadings. This card keeps that split.

Block weights and block scores

In the HISTORICAL Westerhuis 1998 iteration, the current residual block (after the same within-block preprocessing) is used, not a silently concatenated residual. The unnormalized block variable weight is

w_{b} = X_{b}^{⊤} u / (u^{⊤} u)

Each $w_{b}$ is then normalized to unit Euclidean length. The block score under the default Westerhuis block scaling is

t_{b} = X_{b} w_{b} / J_{b} = α_{b} X_{b} w_{b}

Order matters: normalize the weight first, then form the score, then apply the block factor already present in $α_{b}$ . The block-score matrix for one component is

T = [t_{1} ∣ \dots ∣ t_{B}] \in R^{I \times B}

A block score summarizes sample variation from one predictor block along that block direction. Its later components depend on the deflation convention. The educational recovery uses residual processed $X_{b}$ at the current component, which is the super-score-deflation residual, not the original raw block after the first component.

Super weights and super scores

The super weight combines block scores into the global X super score. Under the same historical identities,

w_{sup} = T^{⊤} u / (u^{⊤} u)

after which $w_{sup}$ is normalized to unit length and

t_{sup} = T w_{sup}

The super weight reflects how block scores enter that component under the selected scaling. It is not causal importance, not predictive necessity, and not proof of unique information. Other blocks may carry correlated predictive variation. Do not call a partitioned regression coefficient a block-importance statistic unless a separate, sourced definition is given. This card does not introduce a generic BIP formula.

Relationship with standard PLS

Westerhuis, Kourti and MacGregor show that MB-PLS using super-score deflation can be calculated from ordinary PLS when the predictor variables are given the corresponding scaling. With the $J_{b}$ convention, that ordinary PLS is fit to $Z$ . The MB-PLS super scores then match the standard PLS X scores. Predictions match. The equivalence is formulation-specific, scaling-dependent, and tied to super-score deflation. It is not the sentence "all MB-PLS is ordinary PLS."

MB-PLS still adds scientifically meaningful block organization, explicit block scaling, and block-level quantities. Hierarchical PLS is a different method in the same 1998 paper: different normalization, different block-score construction, and a different role for Y. Westerhuis et al. do not give a standard-PLS equivalence for HPLS. HPLS is not taught here.

The 1998 comparison assumes matching variable scaling and no missing values. Missing data can break the simple concatenated-PLS calculation they describe. Missing blocks and missing entries are outside this educational implementation.

Mathematics

The global model is the NIPALS PLS1/PLS2 formulation already validated on the Open Lab NIPALS card: Mode A weights, regression deflation of $Y$ on the X score, inner-loop Y-weight $c$ not forced to unit length. For $M = 1$ the inner loop is a single pass. For $M > 1$ it iterates until the squared change in $w$ is below the documented tolerance. NIPALS PLS2 and SIMPLS are not assumed to return identical raw latent vectors when $M > 1$ .

Westerhuis et al. state that this MB-PLS has the same objective as the corresponding standard PLS: covariance between the X super score and the Y score is maximized under that PLS convention. This card does not display a separate sample-space eigenproblem. The operational criterion is the NIPALS inner loop of the weighted superblock.

Selected Westerhuis 1998 MB-PLS identities Z = [\frac{X _{1}}{J _{1}} ∣ \dots ∣ \frac{X _{B}}{J _{B}}] (1) w_{b} = \frac{X _{b}^{⊤} u}{u ^{⊤} u}, w_{b} \leftarrow w_{b} / ∥ w_{b} ∥_{2} (2) t_{b} = \frac{X _{b} w _{b}}{J _{b}} (3) T = [t_{1} ∣ \dots ∣ t_{B}] (4) w_{sup} = \frac{T ^{⊤} u}{u ^{⊤} u}, w_{sup} \leftarrow w_{sup} / ∥ w_{sup} ∥_{2} (5) t_{sup} = T w_{sup} (6) c = \frac{Y ^{⊤} t _{sup}}{t _{sup}^{⊤} t _{sup}}, u = \frac{Yc}{c ^{⊤} c} (7) p_{b} = \frac{X _{b}^{⊤} t _{sup}}{t _{sup}^{⊤} t _{sup}}, X_{b} \leftarrow X_{b} - t_{sup} p_{b}^{⊤} (8) Y = Z B^{⊤} + 1 μ_{Y}^{⊤} (9) X_{b} Processed residual of predictor block b at the current component, I by J_b. Z Weighted predictor superblock used by the global NIPALS PLS model. w_{b} Unit-length block variable weight. Not a loading and not a regression vector. t_{sup} X super score. Equal to the NIPALS X-score of Z under this formulation. B NIPALS regression coefficient matrix of shape M by J, applying to the column-centered superblock Z. The intercept is the training mean of Y. Interpretation Equations (1) to (8) are the selected Westerhuis 1998 super-score-deflation identities, with J_{b} in place of their m_{X_{b}} . Equation (9) is the NIPALS prediction convention of the Open Lab PLSR/NIPALS cards. Inner-loop c is a Y-weight. The Y-loading used for deflation is q = Y^{⊤} t_{sup} / (t_{sup}^{⊤} t_{sup}), matching the NIPALS card.

Algorithm

Algorithm 1

Multiblock PLS

InputPredictor blocks X_1,...,X_B, response block Y, component count A, within-block preprocessing, block scaling.

OutputGlobal NIPALS PLS of the weighted superblock, recovered block-level weights/scores/loadings, stored training parameters, predicted Y.

01Verify that every predictor block and Y have I rows in the same order.
02Fit training column means, and optional training column scales, on each X_b and on Y.
03 $X_{b, c} \leftarrow X_{b} - 1 μ_{b}^{⊤}$
04Optionally autoscale predictor columns with the training sample standard deviation.
05 $α_{b} \leftarrow 1/ J_{b} (or α_{b} = 1 if unweighted)$
06 $Z_{b} \leftarrow α_{b} X_{b, p}, Z \leftarrow [Z_{1} ∣ \dots ∣ Z_{B}]$
07Fit Open Lab NIPALS PLS (Mode A, regression deflation) of Z to Y.
08Read global super scores t and Y scores u from that PLS solution.
09Recover block weights, block scores, super weights, and block loadings from residual processed X_b using the Westerhuis 1998 identities.
10Keep block column slices for interpretation. Predict Y from the global PLS model of Z.
11Select A by training-only validation. Do not inspect the external test set.

PRIMARY IMPLEMENTED FORMULATION: Westerhuis 1998 MB-PLS with super-score deflation, computed as Open Lab NIPALS PLS on the weighted superblock plus block recovery. Historical iterative MB-PLS and its deflation variants are discussed in Why deflation matters.

Algorithm 2

Historical MB-PLS iteration from Westerhuis et al. 1998

InputCurrent residual blocks X_b and Y, block factors alpha_b.

OutputOne component: block weights, block scores, super weight, super score, Y-weight, updated u, super-score deflation.

01Initialize u as a non-constant column of the current Y residual.
02 $w_{b} \leftarrow X_{b}^{⊤} u / (u^{⊤} u)$
03Normalize each w_b to unit Euclidean length.
04 $t_{b} \leftarrow α_{b} X_{b} w_{b}$
05 $T \leftarrow [t_{1} ∣ \dots ∣ t_{B}]$
06 $w_{sup} \leftarrow T^{⊤} u / (u^{⊤} u)$
07Normalize w_sup to unit Euclidean length.
08 $t_{sup} \leftarrow T w_{sup}$
09 $c \leftarrow Y^{⊤} t_{sup} / (t_{sup}^{⊤} t_{sup})$
10 $u \leftarrow Yc / (c^{⊤} c)$
11Repeat until the selected score/weight criterion converges. The educational NIPALS core uses squared X-weight change below 1e-6, or one pass for a single response. That tolerance is an implementation parameter.
12 $p_{b} \leftarrow X_{b}^{⊤} t_{sup} / (t_{sup}^{⊤} t_{sup})$
13 $X_{b} \leftarrow X_{b} - t_{sup} p_{b}^{⊤}$
14Deflate Y on the super score using the NIPALS Y-loading, then extract the next component from the residuals.

HISTORICAL FORMULATION shown for explanation. The public code does not run this nested loop as the computational engine. It recovers the same super-score-deflation quantities from NIPALS of Z.

Why deflation matters

Deflation is not an implementation detail. Westerhuis and Smilde 2001 treat three strategies for multiblock PLS. The 2001 body was not opened. The comparison below uses their abstract together with the 1998 super-score versus block-score discussion. No regression-vector recursion $r^{(a)}$ is displayed or coded.

Strategy	What is deflated	Predictive behavior	Block interpretation
Block-score X deflation	Each X_b using its own block score t_b.	Inferior prediction of Y relative to the other two strategies studied in 2001. Variation can be removed from a block even though that entire block-score direction was not used to predict Y.	Associated with the original Wangen-type / block-score pathway. Not recommended as the primary V1 strategy. Block scores can be orthogonal while super scores need not be, under the 1998 discussion of this option.
Super-score X deflation	Each X_b, and Y, using the super score t_sup.	Same predictions as standard PLS with all variables in one large X-block, under matching scaling. This is the PRIMARY IMPLEMENTED FORMULATION.	Later block scores can contain information influenced by other blocks, because X residualization uses the global super score. Super scores are orthogonal; block scores need not be.
Y-only deflation	Y using the super score. X is not residualized.	2001 abstract: the prediction problems of block-score X deflation, and the mixing problem of super-score X deflation, disappear under the study they report.	X remains the original processed blocks, so later block scores can stay functions of their own X_b. DEFLATION ALTERNATIVE only. Not implemented. No r, w, p recursion is copied from memory.

Do not copy PCA orthogonality language onto MB-PLS. Different deflation strategies produce different score properties. This card never claims that all block scores are orthogonal.

Block interpretation

Interpret the global latent variables, the block-level quantities, the relation to $Y$ , and the validated prediction together. Weights without the preprocessing and block scaling context are not interpretable.

Block loadings $p_{b}$ describe how the residual processed block relates to the super score used for deflation. They are not the unit-length block weights $w_{b}$ . Partitioned columns of $B$ are regression coefficients on the block-scaled superblock. They are not a block-importance index.

A large super weight means that, on the processed and block-scaled scale, that block's score contributes strongly to the super score of that component. It does not prove unique information. Adding predictor blocks can add relevant signal, redundant signal, irrelevant variation, or noise. MB-PLS does not guarantee improvement over the best individual block. MB-PLS does not establish causality.

Choosing the number of components

The number of latent variables $A$ is a model-complexity choice. Do not use a universal 95% variance rule, a fixed $A = 2$ , or inspection of the external test RMSEP to select $A$ . For predictive MB-PLS, select $A$ with the Cross-Validation card. Every learned operation, including column means, optional autoscaling, block factors that depend on training statistics, and the PLS fit, must be reconstructed inside each training fold. Metrics follow the Regression Metrics card. This page does not redefine RMSECV, RMSEP, MAE, R², or bias.

Prediction

For new observations, apply stored training means, stored column scales, and stored $α_{b}$ to every block, in the training block order, then concatenate and apply the trained NIPALS coefficients. Do not refit latent variables on the new samples. Equation (9) applies to the processed superblock $Z$ . Coefficients are not automatically in raw original units if autoscaling or block scaling was used.

Python

The listing reuses the exact educational fit_nipals_pls from the NIPALS card (tolerance 1e-6, at most 500 inner iterations). It does not invent a second PLS algorithm. scikit-learn PLSRegression is an independent global cross-check on already processed $Z$ with scale=False. sklearn documentation for 1.9 still defaults to scale=True. The numerical cross-check in this repository used the installed scikit-learn 1.6.1 with the same scale=False setting. . sklearn does not know predictor block boundaries.

The Python package mbpls (Baum and Vermue 2019) is older scientific software with methods NIPALS, UNIPALS, SIMPLS, and KERNEL and default standardize=True. It is not this card's definition. If it is installed, it may be used only as a labelled cross-check under matched options. If it is missing or incompatible, the transparent implementation remains primary. No runtime is invented.

multiblock-pls

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  def _validate_blocks(blocks, Y, label="blocks"):    if not isinstance(blocks, (list, tuple)) or len(blocks) < 1:        raise ValueError("%s must be a non-empty list of 2D arrays." % label)    out = []    n_samples = None    for index, block in enumerate(blocks):        block = np.asarray(block, dtype=float)        if block.ndim != 2:            raise ValueError("Block %d must be a 2D array." % (index + 1))        if min(block.shape) < 1:            raise ValueError("Block %d must have positive dimensions." % (index + 1))        if not np.all(np.isfinite(block)):            raise ValueError("Block %d must contain only finite values." % (index + 1))        if n_samples is None:            n_samples = block.shape[0]        elif block.shape[0] != n_samples:            raise ValueError("All predictor blocks must share the same number of samples.")        out.append(block)    Y = np.asarray(Y, dtype=float)    if Y.ndim == 1:        Y_work = Y.reshape(-1, 1)        one_d = True    elif Y.ndim == 2:        Y_work = Y        one_d = False    else:        raise ValueError("Y must be a 1D vector or a 2D array of samples by responses.")    if Y_work.shape[0] != n_samples:        raise ValueError("Every predictor block and Y must have the same number of observations.")    if Y_work.shape[1] < 1 or not np.all(np.isfinite(Y_work)):        raise ValueError("Y must contain only finite values and at least one response.")    return out, Y_work, one_d  def _block_scale_factors(n_variables, block_scale, processed, tol):    n_blocks = len(n_variables)    factors = np.ones(n_blocks, dtype=float)    if block_scale == "sqrt_nvars":        for index, width in enumerate(n_variables):            factors[index] = 1.0 / np.sqrt(float(width))    elif block_scale == "none":        pass    else:        raise ValueError("block_scale must be 'sqrt_nvars' or 'none'.")    for index, block in enumerate(processed):        if not np.isfinite(factors[index]) or (            block_scale == "sqrt_nvars" and np.linalg.norm(block) <= tol and factors[index] == 0        ):            raise ValueError("Block scaling is undefined for block %d." % (index + 1))    return factors  def _preprocess_blocks(blocks, means, variable_scales, block_scale_factors):    processed = []    scaled = []    for block, mean, scale, factor in zip(blocks, means, variable_scales, block_scale_factors):        centered = (block - mean) / scale        processed.append(centered)        scaled.append(centered * factor)    return processed, scaled  def fit_multiblock_pls(    blocks,    Y,    n_components,    autoscale=False,    block_scale="sqrt_nvars",    tol=1e-6,    max_iter=500,    recovery_tol=1e-12,):    blocks, Y_work, one_d = _validate_blocks(blocks, Y)    n_samples = blocks[0].shape[0]    n_blocks = len(blocks)    n_variables = [block.shape[1] for block in blocks]    n_components = int(n_components)    n_total = sum(n_variables)    if n_components < 1 or n_components > min(n_samples, n_total):        raise ValueError(            "n_components must be between 1 and min(n_samples, n_variables)."        )    if autoscale and n_samples < 2:        raise ValueError("Autoscaling requires at least two observations.")     means = [block.mean(axis=0) for block in blocks]    centered = [block - mean for block, mean in zip(blocks, means)]    variable_scales = []    processed = []    for block in centered:        if autoscale:            scale = block.std(axis=0, ddof=1)            if np.any(~np.isfinite(scale) | (scale == 0)):                raise ValueError(                    "Scaling cannot be applied to a zero variance variable."                )        else:            scale = np.ones(block.shape[1], dtype=float)        variable_scales.append(scale)        processed.append(block / scale)     block_scale_factors = _block_scale_factors(        n_variables, block_scale, processed, recovery_tol    )    scaled = [block * factor for block, factor in zip(processed, block_scale_factors)]    superblock = np.concatenate(scaled, axis=1)    pls = fit_nipals_pls(superblock, Y_work, n_components, tol=tol, max_iter=max_iter)     residuals = [block.copy() for block in processed]    block_weights = [np.zeros((width, n_components)) for width in n_variables]    block_scores = [np.zeros((n_samples, n_components)) for _ in range(n_blocks)]    block_loadings = [np.zeros((width, n_components)) for width in n_variables]    super_weights = np.zeros((n_blocks, n_components))    reconstructed_super = np.zeros((n_samples, n_components))    eps = np.finfo(float).eps     for a in range(n_components):        t = pls["x_scores"][:, a]        u = pls["y_scores"][:, a]        scores = []        for b, block in enumerate(residuals):            raw = (block.T @ u) / (u @ u)            length = np.sqrt(raw @ raw)            if length <= recovery_tol:                weight = np.zeros(block.shape[1], dtype=float)                score = np.zeros(n_samples, dtype=float)            else:                weight = raw / (length + eps)                score = (block @ weight) * block_scale_factors[b]            block_weights[b][:, a] = weight            block_scores[b][:, a] = score            scores.append(score)        T = np.column_stack(scores)        raw_super = (T.T @ u) / (u @ u)        super_length = np.sqrt(raw_super @ raw_super)        if super_length <= recovery_tol:            w_super = np.zeros(n_blocks, dtype=float)            t_super = np.zeros(n_samples, dtype=float)        else:            w_super = raw_super / (super_length + eps)            t_super = T @ w_super        super_weights[:, a] = w_super        reconstructed_super[:, a] = t_super        denom = float(t @ t)        for b, block in enumerate(residuals):            loading = (block.T @ t) / denom            block_loadings[b][:, a] = loading            residuals[b] = block - np.outer(t, loading)     starts = np.cumsum([0] + n_variables)    coef_by_block = [        pls["coef"][:, starts[b] : starts[b + 1]] for b in range(n_blocks)    ]    return {        "n_samples": n_samples,        "n_blocks": n_blocks,        "n_variables": n_variables,        "n_components": n_components,        "n_responses": Y_work.shape[1],        "one_d": one_d,        "pls1": Y_work.shape[1] == 1,        "autoscale": autoscale,        "block_scale": block_scale,        "deflation": "super-score X and Y",        "means": means,        "variable_scales": variable_scales,        "block_scale_factors": block_scale_factors,        "processed_blocks": processed,        "superblock": superblock,        "pls": pls,        "super_scores": pls["x_scores"],        "y_scores": pls["y_scores"],        "y_weights": pls["y_weights"],        "y_loadings": pls["y_loadings"],        "superblock_weights": pls["x_weights"],        "superblock_loadings": pls["x_loadings"],        "block_weights": block_weights,        "block_scores": block_scores,        "block_loadings": block_loadings,        "super_weights": super_weights,        "reconstructed_super_scores": reconstructed_super,        "coef": pls["coef"],        "coef_by_block": coef_by_block,        "intercept": pls["intercept"],        "tol": tol,        "max_iter": max_iter,    }  def predict_multiblock_pls(model, blocks_new):    blocks_new, _, _ = _validate_blocks(blocks_new, np.zeros((blocks_new[0].shape[0], 1)))    if len(blocks_new) != model["n_blocks"]:        raise ValueError("The number of new blocks must match the trained model.")    for index, block in enumerate(blocks_new):        if block.shape[1] != model["n_variables"][index]:            raise ValueError(                "New block %d must have the same number of variables as training."                % (index + 1)            )    _, scaled = _preprocess_blocks(        blocks_new,        model["means"],        model["variable_scales"],        model["block_scale_factors"],    )    superblock = np.concatenate(scaled, axis=1)    y_hat = predict_nipals_pls(model["pls"], superblock)    return {        "superblock": superblock,        "y_hat": y_hat,    }  def mbpls_cv_rmse(    blocks,    Y,    n_components,    n_folds=5,    autoscale=False,    block_scale="sqrt_nvars",):    blocks, Y_work, one_d = _validate_blocks(blocks, Y)    n_samples = Y_work.shape[0]    n_folds = int(n_folds)    if n_folds < 2 or n_folds > n_samples:        raise ValueError("n_folds must be between 2 and n_samples.")    fold = np.array_split(np.arange(n_samples), n_folds)    errors = []    for hold in fold:        train = np.setdiff1d(np.arange(n_samples), hold, assume_unique=False)        model = fit_multiblock_pls(            [block[train] for block in blocks],            Y_work[train],            n_components,            autoscale=autoscale,            block_scale=block_scale,        )        y_hat = predict_multiblock_pls(model, [block[hold] for block in blocks])["y_hat"]        if one_d:            y_hat = np.asarray(y_hat).reshape(-1, 1)        errors.append(np.mean((Y_work[hold] - y_hat) ** 2))    return float(np.sqrt(np.mean(errors)))

MATLAB

The MATLAB listing follows the same primary formulation: preprocess, apply $α_{b}$ , concatenate, run the educational NIPALS core, recover block quantities. It does not depend on plsregress. MathWorks documents plsregress as SIMPLS: it centers and does not rescale columns. de Jong 1993 shows SIMPLS equivalent to PLS1 for univariate y, and not strictly identical to NIPALS-PLS2 for multivariate Y. For PLS1, predictions may be compared after matched preprocessing and sign alignment. For PLS2, compare predictions or the predictive subspace, not raw latent vectors.

multiblock-pls

function model = fit_multiblock_pls(blocks, Y, nComponents, autoscale, blockScale)    if nargin < 4, autoscale = false; end    if nargin < 5, blockScale = "sqrt_nvars"; end    [blocks, Y, oneD] = validateBlocks_mbpls(blocks, Y);    nSamples = size(blocks{1}, 1);    nBlocks = numel(blocks);    nVariables = zeros(1, nBlocks);    for b = 1:nBlocks        nVariables(b) = size(blocks{b}, 2);    end    nTotal = sum(nVariables);    if nComponents < 1 || nComponents > min(nSamples, nTotal)        error('nComponents must be between 1 and min(nSamples, nVariables).');    end    if autoscale && nSamples < 2        error('Autoscaling requires at least two observations.');    end    means = cell(1, nBlocks);    variableScales = cell(1, nBlocks);    processed = cell(1, nBlocks);    blockScaleFactors = ones(1, nBlocks);    for b = 1:nBlocks        means{b} = mean(blocks{b}, 1);        centered = blocks{b} - means{b};        if autoscale            scale = std(centered, 0, 1);            if any(~isfinite(scale) | scale == 0)                error('Scaling cannot be applied to a zero variance variable.');            end        else            scale = ones(1, nVariables(b));        end        variableScales{b} = scale;        processed{b} = centered ./ scale;        if strcmp(string(blockScale), "sqrt_nvars")            blockScaleFactors(b) = 1 / sqrt(nVariables(b));        elseif ~strcmp(string(blockScale), "none")            error('blockScale must be "sqrt_nvars" or "none".');        end    end    scaled = cell(1, nBlocks);    for b = 1:nBlocks        scaled{b} = processed{b} * blockScaleFactors(b);    end    superblock = [scaled{:}];    pls = fitNipalsPls_mbpls(superblock, Y, nComponents);    residuals = processed;    blockWeights = cell(1, nBlocks);    blockScores = cell(1, nBlocks);    blockLoadings = cell(1, nBlocks);    for b = 1:nBlocks        blockWeights{b} = zeros(nVariables(b), nComponents);        blockScores{b} = zeros(nSamples, nComponents);        blockLoadings{b} = zeros(nVariables(b), nComponents);    end    superWeights = zeros(nBlocks, nComponents);    reconstructedSuper = zeros(nSamples, nComponents);    recoveryTol = 1e-12;    for a = 1:nComponents        t = pls.x_scores(:, a);        u = pls.y_scores(:, a);        T = zeros(nSamples, nBlocks);        for b = 1:nBlocks            raw = (residuals{b}' * u) / (u' * u);            lengthW = sqrt(raw' * raw);            if lengthW <= recoveryTol                weight = zeros(nVariables(b), 1);                score = zeros(nSamples, 1);            else                weight = raw / lengthW;                score = (residuals{b} * weight) * blockScaleFactors(b);            end            blockWeights{b}(:, a) = weight;            blockScores{b}(:, a) = score;            T(:, b) = score;        end        rawSuper = (T' * u) / (u' * u);        lengthS = sqrt(rawSuper' * rawSuper);        if lengthS <= recoveryTol            wSuper = zeros(nBlocks, 1);            tSuper = zeros(nSamples, 1);        else            wSuper = rawSuper / lengthS;            tSuper = T * wSuper;        end        superWeights(:, a) = wSuper;        reconstructedSuper(:, a) = tSuper;        denom = t' * t;        for b = 1:nBlocks            loading = (residuals{b}' * t) / denom;            blockLoadings{b}(:, a) = loading;            residuals{b} = residuals{b} - t * loading';        end    end    coefByBlock = cell(1, nBlocks);    start = 1;    for b = 1:nBlocks        stop = start + nVariables(b) - 1;        coefByBlock{b} = pls.coef(:, start:stop);        start = stop + 1;    end    model.nSamples = nSamples;    model.nBlocks = nBlocks;    model.nVariables = nVariables;    model.nComponents = nComponents;    model.oneD = oneD;    model.pls1 = size(Y, 2) == 1;    model.autoscale = autoscale;    model.blockScale = blockScale;    model.deflation = "super-score X and Y";    model.means = means;    model.variableScales = variableScales;    model.blockScaleFactors = blockScaleFactors;    model.processedBlocks = processed;    model.superblock = superblock;    model.pls = pls;    model.superScores = pls.x_scores;    model.yScores = pls.y_scores;    model.blockWeights = blockWeights;    model.blockScores = blockScores;    model.blockLoadings = blockLoadings;    model.superWeights = superWeights;    model.reconstructedSuperScores = reconstructedSuper;    model.coef = pls.coef;    model.coefByBlock = coefByBlock;    model.intercept = pls.intercept;end function out = predict_multiblock_pls(model, blocksNew)    dummyY = zeros(size(blocksNew{1}, 1), 1);    [blocksNew, ~, ~] = validateBlocks_mbpls(blocksNew, dummyY);    if numel(blocksNew) ~= model.nBlocks        error('The number of new blocks must match the trained model.');    end    scaled = cell(1, model.nBlocks);    for b = 1:model.nBlocks        if size(blocksNew{b}, 2) ~= model.nVariables(b)            error('New block %d must have the same number of variables as training.', b);        end        centered = (blocksNew{b} - model.means{b}) ./ model.variableScales{b};        scaled{b} = centered * model.blockScaleFactors(b);    end    superblock = [scaled{:}];    out.superblock = superblock;    out.yHat = (superblock - model.pls.x_mean) * model.pls.coef' + model.pls.intercept;end function [blocks, Y, oneD] = validateBlocks_mbpls(blocks, Y)    if ~iscell(blocks) || isempty(blocks)        error('blocks must be a non-empty cell array of 2D arrays.');    end    nSamples = [];    for b = 1:numel(blocks)        block = double(blocks{b});        if ~ismatrix(block) || any(size(block) < 1)            error('Block %d must be a 2D array with positive dimensions.', b);        end        if any(~isfinite(block(:)))            error('Block %d must contain only finite values.', b);        end        if isempty(nSamples)            nSamples = size(block, 1);        elseif size(block, 1) ~= nSamples            error('All predictor blocks must share the same number of samples.');        end        blocks{b} = block;    end    Y = double(Y);    if isvector(Y)        Y = Y(:);        oneD = true;    else        oneD = false;    end    if size(Y, 1) ~= nSamples        error('Every predictor block and Y must have the same number of observations.');    end    if any(~isfinite(Y(:)))        error('Y must contain only finite values.');    endend function model = fitNipalsPls_mbpls(X, Y, nComponents)    if isvector(Y)        Y = Y(:);    end    tol = 1e-6;    maxIter = 500;    [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);    C = zeros(nTargets, nComponents);    T = zeros(nSamples, nComponents);    U = zeros(nSamples, nComponents);    P = zeros(nFeatures, nComponents);    Q = zeros(nTargets, nComponents);    nIter = zeros(1, nComponents);    for a = 1:nComponents        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)                dw2 = (w - wPrev)' * (w - wPrev);                if dw2 < tol                    break                end            end            if nTargets == 1                break            end            wPrev = w;        end        nIter(a) = it;        t = Xk * w;        uScore = (Yk * c) / (c' * c);        pVec = (Xk' * t) / (t' * t);        qVec = (Yk' * t) / (t' * t);        Xk = Xk - t * pVec';        Yk = Yk - t * qVec';        W(:, a) = w;        C(:, a) = c;        T(:, a) = t;        U(:, a) = uScore;        P(:, a) = pVec;        Q(:, a) = qVec;    end    R = W * pinv(P' * W);    coef = (R * Q')';    model.x_weights = W;    model.y_weights = C;    model.x_scores = T;    model.y_scores = U;    model.x_loadings = P;    model.y_loadings = Q;    model.x_rotations = R;    model.coef = coef;    model.intercept = yMean;    model.x_mean = xMean;    model.y_mean = yMean;    model.n_iter = nIter;end

Practical notes

Report the exact MB-PLS formulation. Super-score deflation, block-score deflation, Y-only deflation, and HPLS are not interchangeable.
Blocks must be row-aligned with Y. Shared sample mode is required for this card.
Within-block preprocessing and whole-block scaling are different operations. 1/sqrt(J_b) is the selected Westerhuis convention, not a universal optimum.
The global model of this card is NIPALS PLS of the weighted superblock. MB-PLS still keeps block structure for interpretation.
Weights and loadings are different quantities. Super weights are not causal importance.
Do not recommend block-score X deflation as the V1 default. It can remove unused block variation and harm prediction.
Super-score X deflation matches standard PLS prediction under the studied scaling, and it can mix block information in later X residuals.
Y-only deflation is described from Westerhuis and Smilde 2001, not implemented.
Block scores are not automatically orthogonal.
PLS1 and PLS2 must be distinguished. NIPALS and SIMPLS must be distinguished when M > 1.
Select A inside training folds. Do not scale on the full data and do not tune on the test set.
More blocks do not guarantee better prediction. A large block contribution does not prove unique information.
MB-PLS is not MB-PCA plus regression, not SO-PLS, not OnPLS, and not N-PLS.
scikit-learn PLSRegression is not MB-PLS. The Python mbpls package is an independent library with its own defaults.

References

1.
Wangen, L. E., & Kowalski, B. R. (1989). A multiblock partial least squares algorithm for investigating complex chemical systems. Journal of Chemometrics, 3(1), 3-20.
doi:10.1002/cem.1180030104
2.
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
3.
Westerhuis, J. A., Kourti, T., & MacGregor, J. F. (1998). Analysis of multiblock and hierarchical PCA and PLS models. Journal of Chemometrics, 12(5), 301-321.
doi:10.1002/(SICI)1099-128X(199809/10)12:5<301::AID-CEM515>3.0.CO;2-S
4.
Westerhuis, J. A., & Smilde, A. K. (2001). Deflation in multiblock PLS. Journal of Chemometrics, 15(5), 485-493.
doi:10.1002/cem.652
5.
Wold, S., Sjostrom, M., & Eriksson, L. (2001). PLS-regression: a basic tool of chemometrics. Chemometrics and Intelligent Laboratory Systems, 58(2), 109-130.
doi:10.1016/S0169-7439(01)00155-1
6.
Baum, A., & Vermue, L. (2019). Multiblock PLS: Block dependent prediction modeling for Python. Journal of Open Source Software, 4(34), 1190.
doi:10.21105/joss.01190
7.
Mishra, P., Roger, J. M., Jouan-Rimbaud-Bouveresse, D., Biancolillo, A., Marini, F., Nordon, A., & Rutledge, D. N. (2021). Recent trends in multi-block data analysis in chemometrics for multi-source data integration. TrAC Trends in Analytical Chemistry, 137, 116206.
doi:10.1016/j.trac.2021.116206
8.
Smilde, A. K., Næs, T., & Liland, K. H. (2022). Multiblock Data Fusion in Statistics and Machine Learning: Applications in the Natural and Life Sciences. Wiley.
doi:10.1002/9781119600978

PLSR

Partial Least Squares Regression

Open

Nonlinear Iterative Partial Least Squares

NIPALS

Open

Multiblock Data Analysis & Data Fusion

Open

MB-PCA

Multiblock PCA / Consensus PCA

Open

Normalization & Scaling

Open

Cross-Validation

Open

Metrics

Regression Metrics

Open

SO-PLS

Sequential and Orthogonalized Partial Least Squares

Open

ComDim

Common Components and Specific Weights Analysis

Open

HPLS

Hierarchical PLS

Coming soon

OnPLS

Coming soon

N-PLS

Coming soon

PO-PLS

Coming soon