Open Lab/Multiblock Data Analysis · Method

Sequential and Orthogonalized Partial Least Squares

SO-PLS

Supervised sequential multiblock PLS: each later predictor block is orthogonalized to previously extracted PLS scores, then used to model the current response residual.

ChemometricsMultiblockPLSData FusionRegressionPythonMATLAB

SO-PLS

What is SO-PLS?

This card teaches one supervised sequential multiblock method: Sequential and Orthogonalized Partial Least Squares in the standard regression formulation of Biancolillo (2016) section 2.3.3, the same two-block algorithm published as Biancolillo and Næs (2019) DHST chapter 6, generalized to $B$ predictor blocks as in Næs et al. (2021). Predictor blocks are entered in a chosen order. Each later block is orthogonalized with respect to the PLS scores already extracted from earlier blocks, then used to model the current response residual.

Næs et al. (2013) is the primary methodological SO-PLS paper cited here as the definitional method reference. Its body was not opened for this card. The equations displayed follow Biancolillo (2016), Næs et al. (2021), and the inspected educational code. SO-PLS is not MB-PLS, not PO-PLS, not OnPLS, and not N-PLS. The names are not synonyms.

Jørgensen, Mevik and Næs (2007) describe LS-PLS. That paper does not name SO-PLS. Later SO-PLS literature identifies the sequential orthogonalization tradition as a precursor and foundation. It is not called here the original SO-PLS paper.

The sequential multiblock regression problem

Several predictor blocks share the same $N$ observations in matching row order. One response block $Y \in R^{N \times M}$ is to be predicted. A compact two-block notation is

Y = XB + ZC + E

with $X \in R^{N \times J}$ , $Z \in R^{N \times L}$ , and residual $E$ . SO-PLS does not estimate ordinary least squares $B$ and $C$ in general. It builds a sequential PLS model with orthogonalization between stages.

MB-PCA finds variation without using $Y$ . SO-PLS is supervised. MB-PLS enters blocks jointly through a weighted superblock. SO-PLS enters blocks sequentially. They do not solve the same algebraic problem. OnPLS belongs to a different family. N-PLS operates on multiway arrays. SO-N-PLS is an extension, not this card.

Why block order matters

SO-PLS is order dependent. The first block models $Y$ directly. Each later block models what remains after earlier blocks, conditional on the selected component counts and preprocessing. Reordering blocks changes the partition of fitted variation. Neither order is automatically correct without domain context.

Næs et al. (2021) note that order is sometimes obvious from the study design, and otherwise must be chosen. Experience in that source suggests order can matter more for interpretation than for prediction. Campos et al. (2018) discuss optimal block order, natural order versus no obvious order, and the combinatorial cost of exhaustive permutation when $B > 3$ . Stepwise SO-PLS is an extension not implemented in V1.

A natural order can arise when blocks follow industrial process stages. That is a literature-supported example, not proof of causality from order alone. Process monitoring applications in Næs et al. (2021) illustrate ordered blocks. No causal claim follows from sequential entry.

First block: PLS regression

Stage 1 is standard PLS regression of the first processed block to $Y$ . Use the Open Lab PLSR and NIPALS cards for PLS theory. This card does not duplicate that material. With $A_{X}$ retained components,

Y_{X} = T_{X} Q_{X}^{⊤}

$T_{X}$ and $Q_{X}$ come from NIPALS PLS on the first block. If $A_{X} = 0$ , the first block contributes nothing at that stage. Zero components are allowed in the educational code, consistent with the R multiblock::sopls documentation.

Response residuals

After the first block,

E_{Y} = Y - T_{X} Q_{X}^{⊤}

$E_{Y}$ is not pure noise. It may contain systematic information not captured by the selected first-block model, information associated with later blocks, and residual error. The next stage uses $E_{Y}$ as the current response to be modeled.

Orthogonalizing the next block

Before the second block enters, its processed matrix is orthogonalized with respect to the selected PLS score space of the first block, not automatically the full column space of raw $X$ :

Z_{orth} = Z - T_{X} (T_{X}^{⊤} T_{X})^{- 1} T_{X}^{⊤} Z

The educational code uses least squares, not an explicit matrix inverse. Numerically, $T_{X}^{⊤} Z_{orth} \approx 0$ . Orthogonal here means zero inner product with the selected score columns. It does not mean statistically independent variables or blocks.

If $A_{X}$ components are used, orthogonalization is with respect to the span of those scores. The component count of earlier blocks affects later blocks.

Modeling additional information

Stage 2 fits PLS of $Z_{orth}$ to $E_{Y}$ :

Y_{Z ∣ X} = T_{Z_{orth}} Q_{Z_{orth}}^{⊤}

The full two-block prediction is

Y = T_{X} Q_{X}^{⊤} + T_{Z_{orth}} Q_{Z_{orth}}^{⊤}

The additional information contributed by $Z$ is additional given $X$ , the chosen order, the selected $A_{X}$ , and the preprocessing. It is not an intrinsic property of $Z$ alone.

Generalization to multiple blocks

For $B$ blocks, any further block is orthogonalized with respect to the scores of all preceding blocks. With $T_{prev} = [T_{1} ∣ \dots ∣ T_{b - 1}]$ ,

X_{b, orth} = [I - T_{prev} (T_{prev}^{⊤} T_{prev})^{- 1} T_{prev}^{⊤}] X_{b}

when $T_{prev}$ has full column rank. Næs et al. (2021) describe the same sequential workflow: PLS on the current block, orthogonalize remaining blocks with respect to previous PLS scores, deflate $Y$ and fit the orthogonalized next block to the deflated response. Each block may retain a different number of components $A_{1}, \dots, A_{B}$ .

Mathematics

All blocks and $Y$ are centered throughout, following Næs et al. (2021). V1 stores training column means for each $X_{b}$ and for $Y$ . Optional within-block autoscaling uses the training sample standard deviation with denominator $N - 1$ in Python and std(...,0,1) in MATLAB. It is not mandatory.

Between-block relative scaling is a separate question. Næs et al. (2021) report invariance, due to the orthogonalization, to between-block scaling, but not to within-block scaling because PLS is used. Biancolillo (2016) gives a scale-invariant sequential formulation on blocks orthogonalized with respect to previous scores. That scope is whole-block relative scaling, not the claim that scaling never matters.

Standard SO-PLS regression (Biancolillo 2016; Næs et al. 2021) Y = XB + ZC + E (1) Y_{X} = T_{X} Q_{X}^{⊤} (2) E_{Y} = Y - T_{X} Q_{X}^{⊤} (3) Z_{orth} = Z - T_{X} (T_{X}^{⊤} T_{X})^{- 1} T_{X}^{⊤} Z (4) Y = T_{X} Q_{X}^{⊤} + T_{Z_{orth}} Q_{Z_{orth}}^{⊤} (5) X_{b, orth} = [I - T_{prev} (T_{prev}^{⊤} T_{prev})^{- 1} T_{prev}^{⊤}] X_{b} (6) X_{b, orth, new} = X_{b, new} - T_{prev, new} C_{b} (7) Y_{new} = 1 μ_{Y}^{⊤} + b \sum T_{b, new} Q_{b}^{⊤} (8) X, Z Two predictor blocks in the canonical two-block notation. Generalized to X_1,...,X_B. T_{X}, Q_{X} PLS scores and Y-loadings from the first block stage. Span defines the orthogonalization projector. Z_{orth} Second block after removal of its least-squares fit to the selected earlier score space. C_{b} Training-only projection coefficients from lstsq(T_prev, X_b). Shape A_prev by J_b. μ_{Y} Training response mean restored at prediction. Interpretation Equations (1) to (5) are the inspected two-block SO-PLS regression workflow. Equation (6) generalizes orthogonalization to B blocks. Equations (7) and (8) are the new-sample prediction convention of Biancolillo (2016), using training-only C_{b} and NIPALS rotations. SO-PLS is not equivalent to estimating B and C by OLS on raw blocks.

Algorithm

Algorithm 1

Sequential orthogonalized PLS

InputPredictor blocks X_1,...,X_B, response Y, component counts A_1,...,A_B, block order, within-block preprocessing.

OutputSequential PLS stages, orthogonalized blocks, stored projection coefficients, block contributions, predicted Y.

01Verify that every predictor block and Y have N rows in the same order.
02Fit training column means of each X_b and of Y. Optionally autoscale predictor columns with the training sample standard deviation.
03 $Y_{res} \leftarrow Y - 1 μ_{Y}^{⊤}$
04Initialize T_prev as an empty score matrix.
05For block b = 1,...,B in the chosen order:
06If T_prev is empty, set X_use to the processed block. Otherwise orthogonalize the processed block against T_prev by least squares.
07 $X_{b, orth} = X_{b, p} - T_{prev} C_{b}, C_{b} = lstsq (T_{prev}, X_{b, p})$
08If A_b = 0, skip PLS for this block and leave its contribution at zero.
09Otherwise fit Open Lab NIPALS PLS of X_use to the current Y_res with A_b components.
10 $Y_{b} = T_{b} Q_{b}^{⊤}, Y_{res} \leftarrow Y_{res} - Y_{b}$
11Append T_b to T_prev. Store projection coefficients, PLS parameters, and the block contribution.
12 $Y = 1 μ_{Y}^{⊤} + \sum_{b = 1}^{B} Y_{b}$
13Select A_1,...,A_B and block order inside training-only validation. Do not inspect the external test set.

PRIMARY IMPLEMENTED FORMULATION: standard SO-PLS regression matching the educational Python and MATLAB listings. PLS substeps follow Open Lab NIPALS.

Interpreting block contributions

Each block contribution $Y_{b}$ is the part of the fitted response allocated at that sequential stage, given earlier blocks, order, preprocessing, and selected $A_{b}$ . It is additional information conditional on blocks already entered. It is not proof of causal importance and not an intrinsic property of the block alone.

A smaller contribution after adding a block does not by itself prove the block is unimportant. A larger contribution does not prove the block carries information unavailable elsewhere. Correlated blocks can exchange roles under reordering. Report order, preprocessing, and component counts together with any block-level summary.

Component selection and validation

Each block may retain a different $A_{b}$ . Biancolillo (2016) and Næs et al. (2021) discuss global versus sequential component-search strategies. The R multiblock::sopls package defaults to global search with sequential=FALSE. Neither strategy is universally superior. Exhaustive order search is costly when $B > 3$ (Campos et al. 2018).

Do not use a universal 95% variance rule or inspection of the external test set to choose complexity. Select $A_{b}$ and, if data-driven, block order with training cross-validation. Every learned operation, including column means, optional autoscaling, component counts, and order when searched, must be reconstructed inside each training fold. Metrics follow the Regression Metrics card.

SO-PLS and Type I ANOVA logic

SO-PLS allocates variation sequentially, analogous to Type I sums of squares in ANOVA (Næs 2010 presentation; Næs et al. 2011 interaction paper abstract; Biancolillo incremental formulation). SO-PLS is not an ANOVA algorithm. It does not produce ANOVA tables or F-tests on training residuals in this V1 card.

Næs et al. (2021) mention cross-validated ANOVA ideas, including paired t-tests on cross-validated residuals for assessing significance of new blocks. That concept is noted here only. No p-values are implemented. Do not treat sequential fit reduction as formal hypothesis testing.

SO-PLS versus MB-PLS

SO-PLS and MB-PLS address the same broad supervised multiblock setting with different algebra. MB-PLS, in the Westerhuis 1998 formulation taught on the MB-PLS card, jointly uses blocks through a weighted superblock and super-score deflation. SO-PLS uses sequential PLS stages with orthogonalization of later blocks against earlier scores. Neither method is universally better.

Topic	SO-PLS	MB-PLS
Block entry	Sequential, fixed or chosen order.	Joint through weighted superblock.
Later-block role	Models response residual after orthogonalization to earlier scores.	Enters global latent structure with all blocks each component.
Order dependence	Yes. Partition of fitted variation depends on order.	Block order in concatenation matters for interpretation of column slices, not sequential allocation.
Relative block scaling	SO-PLS literature reports between-block scale invariance under the sequential orthogonalized formulation. Within-block scaling still matters because PLS is used.	Sensitive to relative block weighting through the superblock construction.
Versus PLSR	PLSR uses one $X$ block. Neither SO-PLS nor MB-PLS equals PLS on a silently concatenated superblock without the corresponding preprocessing and algorithm.
Versus PO-PLS	Næs (2013) and Mishra et al. (2021) distinguish SO-PLS incremental sequential modeling from PO-PLS common and distinct decomposition via a different orthogonalization. PO-PLS is not derived here.

MB-PCA remains unsupervised. OnPLS is mentioned only as a different family. N-PLS and SO-N-PLS belong to multiway extensions outside this matrix-based V1 card.

Prediction

For new observations, never apply a training $N_{train} \times N_{train}$ projector built from training scores to $X_{new}$ . Store training projection coefficients $C_{b} = lstsq (T_{prev, train}, X_{b, train})$ . Process each new block with stored means and scales. Form $T_{prev, new}$ sequentially from new scores. Orthogonalize with equation (7). Compute block scores from stored NIPALS rotations. Sum block contributions and add the training response mean as in equation (8).

Coefficients and scores are not refit on new samples. Order, preprocessing, and retained components must match the trained model.

Python

The listing reuses the exact educational fit_nipals_pls from the NIPALS card. scikit-learn PLSRegression is a PLS building block only, not SO-PLS. If used as a PLS check, set scale=False.

The R package multiblock provides sopls (Kristian Hovde Liland). Its documentation describes sequential PLS, orthogonalization of remaining blocks on extracted components, fit method PKPLS, default scale=FALSE, global versus sequential component search, and allowance of zero components. The public Open Lab code is Python and MATLAB. R is an independent cross-check. R was not runtime tested in this repository because Rscript is not installed.

sequential-orthogonalized-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 orthogonalize_against_scores(X, T):    X = np.asarray(X, dtype=float)    T = np.asarray(T, dtype=float)    if X.ndim != 2:        raise ValueError("X must be a 2D array.")    if T.ndim != 2:        raise ValueError("T must be a 2D score matrix.")    if T.shape[0] != X.shape[0]:        raise ValueError("T and X must have the same number of observations.")    if not np.all(np.isfinite(X)) or not np.all(np.isfinite(T)):        raise ValueError("T and X must contain only finite values.")    if T.shape[1] == 0:        coef = np.zeros((0, X.shape[1]))        return X.copy(), coef    coef, *_ = np.linalg.lstsq(T, X, rcond=None)    return X - T @ coef, coef  def fit_so_pls(blocks, Y, n_components, autoscale=False, tol=1e-6, max_iter=500):    blocks, Y_work, one_d = _validate_blocks(blocks, Y)    n_blocks = len(blocks)    n_components = np.asarray(n_components, dtype=int).reshape(-1)    if n_components.size != n_blocks:        raise ValueError("n_components must contain one nonnegative integer per block.")    if np.any(n_components < 0):        raise ValueError("n_components must be nonnegative integers.")    n_samples = blocks[0].shape[0]    n_variables = [block.shape[1] for block in blocks]    if autoscale and n_samples < 2:        raise ValueError("Autoscaling requires at least two observations.")     y_mean = Y_work.mean(axis=0)    Y_c = Y_work - y_mean    means = []    variable_scales = []    processed = []    for block in blocks:        mean = block.mean(axis=0)        centered = block - mean        if autoscale:            scale = centered.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)        means.append(mean)        variable_scales.append(scale)        processed.append(centered / scale)     Y_res = Y_c.copy()    T_prev = np.zeros((n_samples, 0))    stages = []    contributions = []    for index in range(n_blocks):        Xb = processed[index]        if T_prev.shape[1] == 0:            X_use = Xb.copy()            projection_coef = np.zeros((0, n_variables[index]))        else:            X_use, projection_coef = orthogonalize_against_scores(Xb, T_prev)        A = int(n_components[index])        if A == 0:            pls = None            T_b = np.zeros((n_samples, 0))            Q_b = np.zeros((Y_work.shape[1], 0))            y_hat_b = np.zeros_like(Y_c)        else:            if A > min(X_use.shape):                raise ValueError(                    "Block %d n_components cannot exceed min(n_samples, n_features)."                    % (index + 1)                )            pls = fit_nipals_pls(X_use, Y_res, A, tol=tol, max_iter=max_iter)            T_b = pls["x_scores"]            Q_b = pls["y_loadings"]            y_hat_b = T_b @ Q_b.T            Y_res = Y_res - y_hat_b        contributions.append(y_hat_b)        stages.append(            {                "n_components": A,                "projection_coef": projection_coef,                "X_orth": X_use,                "pls": pls,                "x_scores": T_b,                "y_loadings": Q_b,                "y_hat": y_hat_b,            }        )        if T_b.shape[1] > 0:            T_prev = np.hstack([T_prev, T_b])     y_hat_c = np.sum(contributions, axis=0)    y_hat = y_hat_c + y_mean    residual = Y_work - y_hat    return {        "n_samples": n_samples,        "n_blocks": n_blocks,        "n_variables": n_variables,        "n_components": n_components.copy(),        "one_d": one_d,        "pls1": Y_work.shape[1] == 1,        "autoscale": autoscale,        "means": means,        "variable_scales": variable_scales,        "y_mean": y_mean,        "processed_blocks": processed,        "stages": stages,        "T_prev": T_prev,        "contributions": contributions,        "y_hat_centered": y_hat_c,        "y_hat": y_hat.ravel() if one_d else y_hat,        "residual": residual.ravel() if one_d else residual,        "Y_centered": Y_c,    }  def predict_so_pls(model, blocks_new):    dummy = np.zeros((np.asarray(blocks_new[0]).shape[0], 1))    blocks_new, _, _ = _validate_blocks(blocks_new, dummy, label="blocks_new")    if len(blocks_new) != model["n_blocks"]:        raise ValueError("The number of new blocks must match the trained model.")    n_new = blocks_new[0].shape[0]    T_prev = np.zeros((n_new, 0))    y_hat_c = np.zeros((n_new, model["y_mean"].shape[0]))    stage_scores = []    stage_orth = []    contributions = []    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)            )        Xb = (block - model["means"][index]) / model["variable_scales"][index]        stage = model["stages"][index]        coef = stage["projection_coef"]        if T_prev.shape[1] == 0:            X_use = Xb        else:            X_use = Xb - T_prev @ coef        A = int(stage["n_components"])        if A == 0:            T_b = np.zeros((n_new, 0))            y_hat_b = np.zeros_like(y_hat_c)        else:            pls = stage["pls"]            T_b = (X_use - pls["x_mean"]) @ pls["x_rotations"]            y_hat_b = T_b @ pls["y_loadings"].T        contributions.append(y_hat_b)        y_hat_c = y_hat_c + y_hat_b        stage_scores.append(T_b)        stage_orth.append(X_use)        if T_b.shape[1] > 0:            T_prev = np.hstack([T_prev, T_b])    y_hat = y_hat_c + model["y_mean"]    if model["one_d"]:        y_hat = y_hat.ravel()        y_hat_c = y_hat_c.ravel()    return {        "y_hat": y_hat,        "y_hat_centered": y_hat_c,        "contributions": contributions,        "stage_scores": stage_scores,        "X_orth": stage_orth,    }  def sopls_cv_rmse(blocks, Y, n_components, n_folds=5, autoscale=False):    blocks, Y_work, one_d = _validate_blocks(blocks, Y)    n_samples = blocks[0].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_so_pls(            [block[train] for block in blocks],            Y_work[train],            n_components,            autoscale=autoscale,        )        y_hat = predict_so_pls(model, [block[hold] for block in blocks])["y_hat"]        y_hat = np.asarray(y_hat).reshape(-1, Y_work.shape[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 sequential workflow with the educational NIPALS core. It does not use plsregress. MathWorks documents plsregress as SIMPLS. Orthogonalization uses mldivide(T, X), the documented MathWorks least-squares solver equivalent to T \\ X. Sample standard deviation uses std(X,0,1).

sequential-orthogonalized-pls

function model = fit_so_pls(blocks, Y, nComponents, autoscale)    if nargin < 4, autoscale = false; end    [blocks, Y, oneD] = validateBlocks_sopls(blocks, Y);    nBlocks = numel(blocks);    nComponents = nComponents(:)';    if numel(nComponents) ~= nBlocks        error('nComponents must contain one nonnegative integer per block.');    end    if any(nComponents < 0 | nComponents ~= floor(nComponents))        error('nComponents must be nonnegative integers.');    end    nSamples = size(blocks{1}, 1);    nVariables = zeros(1, nBlocks);    for b = 1:nBlocks        nVariables(b) = size(blocks{b}, 2);    end    if autoscale && nSamples < 2        error('Autoscaling requires at least two observations.');    end    yMean = mean(Y, 1);    Yc = Y - yMean;    means = cell(1, nBlocks);    variableScales = cell(1, nBlocks);    processed = cell(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;    end    Yres = Yc;    Tprev = zeros(nSamples, 0);    stages = cell(1, nBlocks);    contributions = cell(1, nBlocks);    for b = 1:nBlocks        Xb = processed{b};        if size(Tprev, 2) == 0            Xuse = Xb;            projectionCoef = zeros(0, nVariables(b));        else            [Xuse, projectionCoef] = orthogonalize_against_scores(Xb, Tprev);        end        A = nComponents(b);        if A == 0            Tb = zeros(nSamples, 0);            Qb = zeros(size(Y, 2), 0);            yHatB = zeros(size(Yc));            pls = [];        else            if A > min(size(Xuse))                error('Block %d nComponents cannot exceed min(nSamples, nFeatures).', b);            end            pls = fitNipalsPls_sopls(Xuse, Yres, A);            Tb = pls.x_scores;            Qb = pls.y_loadings;            yHatB = Tb * Qb';            Yres = Yres - yHatB;        end        contributions{b} = yHatB;        stage.n_components = A;        stage.projection_coef = projectionCoef;        stage.X_orth = Xuse;        stage.pls = pls;        stage.x_scores = Tb;        stage.y_loadings = Qb;        stage.y_hat = yHatB;        stages{b} = stage;        if size(Tb, 2) > 0            Tprev = [Tprev, Tb];        end    end    yHatC = zeros(size(Yc));    for b = 1:nBlocks        yHatC = yHatC + contributions{b};    end    yHat = yHatC + yMean;    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.means = means;    model.variableScales = variableScales;    model.yMean = yMean;    model.processedBlocks = processed;    model.stages = stages;    model.Tprev = Tprev;    model.contributions = contributions;    model.yHatCentered = yHatC;    model.yHat = yHat;    model.residual = Y - yHat;end function out = predict_so_pls(model, blocksNew)    dummyY = zeros(size(blocksNew{1}, 1), 1);    [blocksNew, ~, ~] = validateBlocks_sopls(blocksNew, dummyY);    if numel(blocksNew) ~= model.nBlocks        error('The number of new blocks must match the trained model.');    end    nNew = size(blocksNew{1}, 1);    Tprev = zeros(nNew, 0);    yHatC = zeros(nNew, numel(model.yMean));    stageScores = cell(1, model.nBlocks);    stageOrth = cell(1, model.nBlocks);    contributions = 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        Xb = (blocksNew{b} - model.means{b}) ./ model.variableScales{b};        stage = model.stages{b};        if size(Tprev, 2) == 0            Xuse = Xb;        else            Xuse = Xb - Tprev * stage.projection_coef;        end        A = stage.n_components;        if A == 0            Tb = zeros(nNew, 0);            yHatB = zeros(size(yHatC));        else            pls = stage.pls;            Tb = (Xuse - pls.x_mean) * pls.x_rotations;            yHatB = Tb * pls.y_loadings';        end        contributions{b} = yHatB;        yHatC = yHatC + yHatB;        stageScores{b} = Tb;        stageOrth{b} = Xuse;        if size(Tb, 2) > 0            Tprev = [Tprev, Tb];        end    end    out.yHat = yHatC + model.yMean;    out.yHatCentered = yHatC;    out.contributions = contributions;    out.stageScores = stageScores;    out.Xorth = stageOrth;end function [Xorth, coef] = orthogonalize_against_scores(X, T)    if size(T, 2) == 0        coef = zeros(0, size(X, 2));        Xorth = X;        return    end    coef = mldivide(T, X);    Xorth = X - T * coef;end function [blocks, Y, oneD] = validateBlocks_sopls(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_sopls(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 block order, preprocessing, and A_1,...,A_B. SO-PLS results are conditional on all three.
Blocks must be row-aligned with Y. Shared sample mode is required.
Within-block preprocessing and between-block relative scaling are different operations. Scaling still matters inside blocks because PLS is used.
Orthogonalized blocks are orthogonal to selected earlier score spaces, not automatically statistically independent.
Later-block contribution is additional information conditional on blocks already entered. Do not call it unique.
SO-PLS is not OLS on raw blocks and not PLS on a silently concatenated superblock.
More blocks do not guarantee better prediction. Fit reduction is not proof of block importance.
Order can matter more for interpretation than prediction. Stepwise SO-PLS is not implemented in V1.
Zero components per block are allowed and yield zero contribution at that stage.
Select complexity and, if searched, order inside training folds only.
Variable selection, VIP, SR, stepwise extensions, interactions, SO-N-PLS, and SO-PLS-LDA are extensions not in V1.
SO-PLS is not MB-PLS, PO-PLS, OnPLS, MB-PCA, or N-PLS.
sklearn PLSRegression is not SO-PLS. R multiblock::sopls is an independent cross-check, not runtime tested here.

References

1.
Jørgensen, K., Mevik, B.-H., & Næs, T. (2007). Combining designed experiments with several blocks of spectroscopic data. Chemometrics and Intelligent Laboratory Systems, 88(2), 154-166.
doi:10.1016/j.chemolab.2007.04.002
2.
Næs, T., Tomic, O., Afseth, N. K., Segtnan, V., & Måge, I. (2013). Multi-block regression based on combinations of orthogonalisation, PLS-regression and canonical correlation analysis. Chemometrics and Intelligent Laboratory Systems, 124, 32-42.
doi:10.1016/j.chemolab.2013.03.006
3.
Biancolillo, A., & Næs, T. (2019). The Sequential and Orthogonalized PLS Regression for Multiblock Regression: Theory, Examples, and Extensions. In Cocchi, M. (Ed.), Data Fusion Methodology and Applications. Data Handling in Science and Technology, 31, 157-177.
doi:10.1016/B978-0-444-63984-4.00006-5
4.
Biancolillo, A. (2016). Method development in the area of multi-block analysis focused on food analysis (PhD thesis). University of Copenhagen.
5.
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
6.
Campos, M. P., Sousa, R., & Reis, M. S. (2018). Establishing the optimal blocks' order in SO-PLS: Stepwise SO-PLS and alternative formulations. Journal of Chemometrics, 32(8), e3032.
doi:10.1002/cem.3032
7.
Næs, T., Romano, R., Tomic, O., Måge, I., Smilde, A., & Liland, K. H. (2021). Sequential and orthogonalized PLS (SO-PLS) regression for path analysis: Order of blocks and relations between effects. Journal of Chemometrics, 35(10), e3243.
doi:10.1002/cem.3243
8.
Wold, S., Sjöström, 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
9.
Næs, T., Måge, I., & Segtnan, V. H. (2011). Incorporating interactions in multi-block sequential and orthogonalised partial least squares regression. Journal of Chemometrics, 25(11), 601-609.
doi:10.1002/cem.1406
10.
Biancolillo, A., Liland, K. H., Måge, I., Næs, T., & Bro, R. (2016). Variable selection in multi-block regression. Chemometrics and Intelligent Laboratory Systems, 156, 89-101.
doi:10.1016/j.chemolab.2016.05.016
11.
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
12.
Liland, K. H. (n.d.). multiblock: Multivariate block methods. R package documentation, function sopls. CRAN / khliland.github.io.

PLSR

Partial Least Squares Regression

Open

Nonlinear Iterative Partial Least Squares

NIPALS

Open

MB-PLS

Multiblock Partial Least Squares

Open

Multiblock Data Analysis & Data Fusion

Open

MB-PCA

Multiblock PCA / Consensus PCA

Open

Normalization & Scaling

Open

Cross-Validation

Open

Metrics

Regression Metrics

Open

ComDim

Common Components and Specific Weights Analysis

Open

PO-PLS

Coming soon

OnPLS

Coming soon

N-PLS

N-PLS / SO-N-PLS

Coming soon

Stepwise SO-PLS

Coming soon

ROSA

Coming soon