Open Lab/Multiblock Data Analysis · Method

Multiblock PCA / Consensus PCA

MB-PCA

Unsupervised multiblock PCA in the CPCA-W / SUM-PCA equivalent formulation: processed superblock PCA with block-level contributions, loadings, and scores.

ChemometricsMultiblockPCAData FusionSVDPythonMATLAB

MB-PCA

01

What is Multiblock PCA?

This card teaches unsupervised multiblock principal component analysis in one exact formulation: Consensus PCA with unit-length block loadings, written CPCA-W by Smilde, Westerhuis and de Jong, and shown by them to be identical to SUM-PCA. SUM-PCA is PCA of the concatenated processed superblock. The names Multiblock PCA, Consensus PCA, CPCA, CPCA-W, Hierarchical PCA, and SUM-PCA are not universal synonyms.

The scientific purpose is not a different global latent space from PCA of that processed superblock. Westerhuis, Kourti and MacGregor showed that results of the Consensus PCA formulation they studied can be calculated from standard PCA when the same variable scaling is used. Qin, Valle and Piovoso state that both scores and loadings of the Consensus PCA they analysed can be obtained from regular PCA. Those equivalences are formulation-specific. The multiblock representation keeps block identity so that contributions, block loadings, and block scores can be recovered after the global PCA.

Original Consensus PCA of Wold and co-workers is a different iterative algorithm. Hierarchical PCA is a different algorithm. Neither is implemented here. No external response enters the objective.

02

Data structure

The formulation uses shared sample mode. For blocks measured on the same observations in the same row order,

This is the shared-object layout of Smilde, Westerhuis and de Jong. It is not every multiblock topology. Independent PCA of each block separately is not this model: those analyses do not share one global component space obtained from the joint superblock.

03

From blocks to a superblock

After within-block preprocessing and the chosen whole-block scaling, concatenate along the variable mode

Concatenation is an algebraic device. The scientific blocking remains: which columns belong to which block. The concatenated array is not a PARAFAC tensor.

04

Why block scaling matters

Within-block preprocessing and between-block scaling are different operations. Centering, and optional autoscaling to unit sample standard deviation, act on variables inside a block. Whole-block scaling multiplies an entire processed block by a factor and thereby changes its contribution to the SUM-PCA criterion.

If and block is multiplied by , then becomes . Its weight in therefore changes by . CPCA-W / SUM-PCA does not automatically make every block equally important.

Smilde, Westerhuis and de Jong discuss prescaling blocks to equal total variance as one way to introduce fairness before SUM-PCA. The educational option frobenius divides a nonzero processed block by its Frobenius norm so that . That is one source-supported choice. It is not a universal MB-PCA rule. Unscaled blocks are not automatically wrong. A zero block makes unit-sum-of-squares scaling undefined.

05

The SUM-PCA / CPCA-W model

Smilde, Westerhuis and de Jong write SUM-PCA as PCA of the superblock. Each processed block is modelled with the same orthonormal super-score matrix and its own block loadings

subject to . This is their equation (1), with their written here as so that it is not confused with conventional PCA scores. Equivalently, . CPCA-W solves the same problem. Original CPCA does not.

06

Global scores and loadings

Let the processed superblock have the thin SVD . For retained components this page uses two sample-space matrices that must not be given the same name.

contains the orthonormal SUM-PCA super-score directions: . contains ordinary chemometric PCA scores. contains ordinary PCA loadings. They satisfy and . The SUM-PCA block loading in the least-squares model is , which is the corresponding block rows of , not the unit-norm CPCA-W loading defined below.

07

Block loadings and block contributions

Partition the conventional loading matrix by variable boundaries

Those rows are partitioned PCA loadings. They are not the CPCA-W unit-norm block loadings.

For an orthonormal direction Smilde, Westerhuis and de Jong define the explained sum of squares of block as with . Equivalently

Here is the SUM-PCA sample-space eigenvalue, equal to the squared singular value of the processed superblock. It is not the sample-covariance eigenvalue used on the PCA card unless that denominator is introduced separately. The derived share is a diagnostic for that component. It is not the historical CPCA weight vector, which Smilde et al. show is different for original CPCA and for CPCA-W.

08

Block scores

For CPCA-W, Smilde, Westerhuis and de Jong normalize the unnormalized block loading and then form the block score. For component and a nonzero map ,

These are CPCA-W block loadings and block scores. They are not the global scores or . If , the normalized loading is undefined. The code then returns a zero vector and marks the quantity as undefined rather than dividing by zero.

09

Mathematics

SUM-PCA / CPCA-W on the processed superblock

(1)
(2)
(3)
(4)
(5)
(6)
processed block b used in the superblock
orthonormal SUM-PCA super-score direction (Smilde T_sup column)
conventional PCA scores Q Sigma
explained sum of squares of block b along q_a

Interpretation

Equations (1) to (5) are SUM-PCA / CPCA-W. Equation (6) is the CPCA-W block loading and block score construction. Do not insert conventional PCA scores T into (4) without removing the singular-value scale.

For retained orthonormal directions, the global subspace reconstruction of a block is , and . The ratio of that quantity to is a block explained sum-of-squares fraction on the processed scale. It is not an ordinary coefficient of determination unless those conditions are stated.

Super-score directions are orthogonal. Smilde, Westerhuis and de Jong note that block-level quantities do not inherit all of those orthogonality properties. Block scores from different global components need not be orthogonal.

10

Algorithm

Algorithm 1

Multiblock PCA via CPCA-W / SUM-PCA

InputBlocks sharing observations, component count A, within-block preprocessing, whole-block scaling.

OutputGlobal PCA of the processed superblock, partitioned loadings, block contributions, optional CPCA-W block scores, stored training parameters.

  1. 01Validate that all blocks share the sample count I.
  2. 02
  3. 03
  4. 04Optionally autoscale columns with the training sample standard deviation.
  5. 05Optionally scale each processed block to unit Frobenius norm.
  6. 06
  7. 07
  8. 08
  9. 09Partition P by block variable boundaries.
  10. 10
  11. 11If CPCA-W block quantities are requested, normalize X_b^T q_a and form t_ba = X_b p_ba.

This is PCA of the processed concatenated superblock plus block-level recovery. It is not the historical iterative CPCA algorithm of Wold et al.

11

CPCA, CPCA-W and HPCA

Wold, Hellberg, Lundstedt, Sjöström and Wold 1987 is the historical Consensus PCA source cited by Smilde et al. The proceedings text was not opened for this card. Mathematical statements about original CPCA come from Smilde, Westerhuis and de Jong 2003. That CPCA uses an iterative consensus loop. After convergence the super-score is proportional to the mean of the block scores. The implied eigenproblem weights each by , so blocks that contribute strongly are downweighted. That is their built-in fairness tendency. The weights depend on the unknown solution. There is no closed SUM-PCA-style eigenproblem. Historical implementations had convergence problems that motivated normalization changes. Original CPCA is described here only. It is not the educational code.

CPCA-W normalizes block loadings to length one. Smilde et al. prove that CPCA-W is exactly SUM-PCA and therefore can be computed by PCA of the processed superblock. Westerhuis et al. 1998 and Qin et al. 2001 established, for the Consensus PCA formulations they analysed, that scores and, in Qin et al., loadings can be obtained from regular PCA under matching scaling. Those statements are not a licence to call every historical CPCA identical to PCA.

Hierarchical PCA of Wold, Kettaneh-Wold and Tjessem 1996 is a different hierarchical construction. Smilde et al. show that HPCA weights blocks that already contribute strongly even more, acting as a block selector, and that it lacks a closed solution of the SUM-PCA type. HPCA is not taught here.

Tchandao Mangamana, Cariou, Vigneau, Glèlè Kakaï and Qannari 2019 call MB-PCA Consensus PCA in a unified unsupervised framework of global and block components. Their eigenanalysis recovers the first standardized principal component of the concatenated matrix. That modern naming is not used here to erase the CPCA versus CPCA-W distinction of Smilde 2003.

12

Choosing the number of components

SUM-PCA uses one common for all blocks. Blocks do not choose independent ranks in this formulation. There is no universal rule such as 95% explained variation, eigenvalue greater than 1, or a fixed component count. Useful evidence may include the global sum of squares , block-wise fractions on the processed scale, stability, and interpretability. If predictive claims are made about later supervised use of the scores, rank decisions belong inside the training structure of the Cross-Validation card.

13

Interpretation

Global scores describe observations in the latent space of the processed combined blocks. They depend on within-block preprocessing, whole-block scaling, and the variance structure of all blocks. They should not be interpreted without reporting those choices. Component signs are arbitrary, as on the PCA card.

A large means that, on the processed scale, block contributes substantial explained sum of squares along . That is not automatic scientific importance. A global component can dominate the superblock while representing blocks unequally. That is why block-wise quantities exist.

Partitioned PCA loadings show how variables of a block enter the global loading vector. CPCA-W block loadings are unit vectors used to construct block scores. Do not mix those two definitions in one sentence of interpretation. Block scores show how observations project in that block-specific direction. Similar block-score patterns do not prove a causal relationship between blocks.

Mishra et al. discuss multi-source spectroscopic tables on the same samples as a chemometric motive for multiblock methods. This page does not attach claims to a particular unpublished sensor pair. New observations are projected only when every trained block is present, in the training variable order, using stored means, stored column scales, and stored block factors. Missing values and missing blocks are outside this educational SVD.

14

Python

The listing uses numpy.linalg.svd with full_matrices=False. Conventional scores are . Contributions use . Autoscaling uses the sample standard deviation with denominator , matching the Normalization & Scaling card. scikit-learn PCA may be used as an independent check on the processed superblock. It does not know block structure and it centers columns. Do not double-scale. Block logic is calculated here, not by sklearn.

import numpy as np  def _validate_blocks(blocks, 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 blocks must share the same number of samples.")        out.append(block)    return out  def _preprocess_block(block, mean, variable_scale, block_scale_factor):    processed = block - mean    processed = processed / variable_scale    return processed * block_scale_factor  def fit_multiblock_pca(    blocks,    n_components,    autoscale=False,    block_scale="none",    tol=1e-12,):    blocks = _validate_blocks(blocks)    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 block_scale not in ("none", "frobenius"):        raise ValueError("block_scale must be 'none' or 'frobenius'.")    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 = np.ones(n_blocks, dtype=float)    scaled = []    for index, block in enumerate(processed):        if block_scale == "frobenius":            norm = np.linalg.norm(block, ord="fro")            if not np.isfinite(norm) or norm <= tol:                raise ValueError(                    "Unit-sum-of-squares block scaling is undefined for a zero block."                )            block_scale_factors[index] = 1.0 / norm        scaled.append(block * block_scale_factors[index])     superblock = np.concatenate(scaled, axis=1)    U, singular_values, Vt = np.linalg.svd(superblock, full_matrices=False)    Q = U[:, :n_components]    sigma = singular_values[:n_components]    loadings = Vt[:n_components, :].T    scores = Q * sigma     pca_loadings_by_block = []    start = 0    for width in n_variables:        pca_loadings_by_block.append(loadings[start : start + width, :])        start += width     contributions = np.zeros((n_blocks, n_components))    cpca_block_loadings = []    cpca_block_scores = []    defined = np.ones((n_blocks, n_components), dtype=bool)    for b, block in enumerate(scaled):        block_loadings = np.zeros((block.shape[1], n_components))        block_scores = np.zeros((n_samples, n_components))        for a in range(n_components):            direction = Q[:, a]            residual_map = block.T @ direction            value = float(np.dot(residual_map, residual_map))            contributions[b, a] = value            length = np.sqrt(value)            if length <= tol:                defined[b, a] = False            else:                loading = residual_map / length                block_loadings[:, a] = loading                block_scores[:, a] = block @ loading        cpca_block_loadings.append(block_loadings)        cpca_block_scores.append(block_scores)     totals = contributions.sum(axis=0)    share = np.divide(        contributions,        totals,        out=np.zeros_like(contributions),        where=totals > tol,    )    block_ss = np.array(        [float(np.linalg.norm(block, ord="fro") ** 2) for block in scaled]    )    block_explained_ss = contributions.sum(axis=1)    block_explained_fraction = np.divide(        block_explained_ss,        block_ss,        out=np.zeros(n_blocks),        where=block_ss > tol,    )     return {        "n_samples": n_samples,        "n_blocks": n_blocks,        "n_variables": n_variables,        "n_components": n_components,        "autoscale": autoscale,        "block_scale": block_scale,        "tol": tol,        "means": means,        "variable_scales": variable_scales,        "block_scale_factors": block_scale_factors,        "processed_blocks": scaled,        "superblock": superblock,        "Q": Q,        "scores": scores,        "loadings": loadings,        "singular_values": sigma,        "pca_loadings_by_block": pca_loadings_by_block,        "block_contributions": contributions,        "block_contribution_share": share,        "block_ss": block_ss,        "block_explained_ss": block_explained_ss,        "block_explained_fraction": block_explained_fraction,        "cpca_block_loadings": cpca_block_loadings,        "cpca_block_scores": cpca_block_scores,        "cpca_block_loading_defined": defined,    }  def transform_multiblock_pca(model, blocks_new):    blocks_new = _validate_blocks(blocks_new, label="blocks_new")    if len(blocks_new) != model["n_blocks"]:        raise ValueError("The number of new blocks must match the trained model.")    processed = []    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)            )        processed.append(            _preprocess_block(                block,                model["means"][index],                model["variable_scales"][index],                model["block_scale_factors"][index],            )        )    superblock = np.concatenate(processed, axis=1)    return {        "processed_blocks": processed,        "superblock": superblock,        "scores": superblock @ model["loadings"],    } 
15

MATLAB

The MATLAB listing uses base svd(..., "econ") so that the Statistics and Machine Learning Toolbox is not required. Blocks are a cell array. The algebra matches the Python listing, including CPCA-W block quantities and stored training parameters for projection. Signs of singular vectors may differ from NumPy. Align signs before comparing implementations.

function model = fit_multiblock_pca(blocks, nComponents, autoscale, blockScale, tol)    if nargin < 3, autoscale = false; end    if nargin < 4, blockScale = "none"; end    if nargin < 5, tol = 1e-12; end    blocks = validateBlocks_mbpca(blocks, "blocks");    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);    scaled = 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 = centered ./ scale;        if strcmp(string(blockScale), "frobenius")            nrm = norm(processed, "fro");            if ~isfinite(nrm) || nrm <= tol                error('Unit-sum-of-squares block scaling is undefined for a zero block.');            end            blockScaleFactors(b) = 1 / nrm;        elseif ~strcmp(string(blockScale), "none")            error('blockScale must be "none" or "frobenius".');        end        scaled{b} = processed * blockScaleFactors(b);    end     superblock = [scaled{:}];    [U, S, V] = svd(superblock, "econ");    Q = U(:, 1:nComponents);    sigma = diag(S);    sigma = sigma(1:nComponents);    loadings = V(:, 1:nComponents);    scores = Q * diag(sigma);     pcaLoadingsByBlock = cell(1, nBlocks);    start = 1;    for b = 1:nBlocks        stop = start + nVariables(b) - 1;        pcaLoadingsByBlock{b} = loadings(start:stop, :);        start = stop + 1;    end     contributions = zeros(nBlocks, nComponents);    cpcaBlockLoadings = cell(1, nBlocks);    cpcaBlockScores = cell(1, nBlocks);    defined = true(nBlocks, nComponents);    for b = 1:nBlocks        block = scaled{b};        blockLoadings = zeros(nVariables(b), nComponents);        blockScores = zeros(nSamples, nComponents);        for a = 1:nComponents            direction = Q(:, a);            residualMap = block' * direction;            value = residualMap' * residualMap;            contributions(b, a) = value;              mapNorm = sqrt(value);              if mapNorm <= tol                defined(b, a) = false;              else                loading = residualMap / mapNorm;                blockLoadings(:, a) = loading;                blockScores(:, a) = block * loading;            end        end        cpcaBlockLoadings{b} = blockLoadings;        cpcaBlockScores{b} = blockScores;    end     totals = sum(contributions, 1);    share = zeros(size(contributions));    mask = totals > tol;    share(:, mask) = contributions(:, mask) ./ totals(mask);    blockSS = zeros(1, nBlocks);    for b = 1:nBlocks        blockSS(b) = norm(scaled{b}, "fro")^2;    end    blockExplainedSS = sum(contributions, 2)';    blockExplainedFraction = zeros(1, nBlocks);    positive = blockSS > tol;    blockExplainedFraction(positive) = blockExplainedSS(positive) ./ blockSS(positive);     model.nSamples = nSamples;    model.nBlocks = nBlocks;    model.nVariables = nVariables;    model.nComponents = nComponents;    model.autoscale = autoscale;    model.blockScale = blockScale;    model.tol = tol;    model.means = means;    model.variableScales = variableScales;    model.blockScaleFactors = blockScaleFactors;    model.processedBlocks = scaled;    model.superblock = superblock;    model.Q = Q;    model.scores = scores;    model.loadings = loadings;    model.singularValues = sigma;    model.pcaLoadingsByBlock = pcaLoadingsByBlock;    model.blockContributions = contributions;    model.blockContributionShare = share;    model.blockSS = blockSS;    model.blockExplainedSS = blockExplainedSS;    model.blockExplainedFraction = blockExplainedFraction;    model.cpcaBlockLoadings = cpcaBlockLoadings;    model.cpcaBlockScores = cpcaBlockScores;    model.cpcaBlockLoadingDefined = defined;end function out = transform_multiblock_pca(model, blocksNew)    blocksNew = validateBlocks_mbpca(blocksNew, "blocksNew");    if numel(blocksNew) ~= model.nBlocks        error('The number of new blocks must match the trained model.');    end    processed = 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};        processed{b} = (centered ./ model.variableScales{b}) * model.blockScaleFactors(b);    end    superblock = [processed{:}];    out.processedBlocks = processed;    out.superblock = superblock;    out.scores = superblock * model.loadings;end function blocks = validateBlocks_mbpca(blocks, label)    if ~iscell(blocks) || isempty(blocks)        error('%s must be a non-empty cell array of 2D arrays.', label);    end    nSamples = [];    for b = 1:numel(blocks)        block = 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 blocks must share the same number of samples.');        end        blocks{b} = double(block);    endend 
16

Practical notes

  • Blocks must be scientifically meaningful and row-aligned. Shared sample mode is required for this formulation.
  • Report within-block preprocessing and whole-block scaling. They are different operations.
  • CPCA-W / SUM-PCA is PCA of the processed superblock. Original CPCA is not that algorithm.
  • Unit Frobenius block scaling is optional. Equal block sum of squares does not imply equal contribution to every component.
  • Inspect block contributions and block explained fractions. Global variance can hide unequal block representation.
  • Do not confuse partitioned PCA loadings with CPCA-W unit-norm block loadings, or global scores with block scores.
  • Global super-score directions are orthogonal. Block scores need not be.
  • A large statistical contribution is not automatic scientific importance.
  • Do not choose A by a universal percentage or eigenvalue rule.
  • Project new samples with training means, column scales, and block factors. Do not refit on the new data.
  • This method is unsupervised. MB-PLS uses a response and belongs on a separate card.
  • PARAFAC models a tensor. Concatenated matrices are not a PARAFAC decomposition.
  • HPCA is a different hierarchical method. It is not implemented here.
17

References

  1. 1.

    Wold, S., Hellberg, S., Lundstedt, T., Sjöström, M., & Wold, H. (1987). PLS modeling with latent variables in two or more dimensions. Proceedings of the PLS Meeting, Frankfurt, Germany.

  2. 2.

    Wold, S., Kettaneh-Wold, N., & Tjessem, K. (1996). Hierarchical multiblock PLS and PC models for easier model interpretation and as an alternative to variable selection. Journal of Chemometrics, 10(5-6), 463-482.

    doi:10.1002/(SICI)1099-128X(199609)10:5/6<463::AID-CEM445>3.0.CO;2-L
  3. 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. 4.

    Qin, S. J., Valle, S., & Piovoso, M. J. (2001). On unifying multiblock analysis with application to decentralized process monitoring. Journal of Chemometrics, 15(9), 715-742.

    doi:10.1002/cem.667
  5. 5.

    Smilde, A. K., Westerhuis, J. A., & de Jong, S. (2003). A framework for sequential multiblock component methods. Journal of Chemometrics, 17(6), 323-337.

    doi:10.1002/cem.811
  6. 6.

    Tchandao Mangamana, E., Cariou, V., Vigneau, E., Glèlè Kakaï, R. L., & Qannari, E. M. (2019). Unsupervised multiblock data analysis: A unified approach and extensions. Chemometrics and Intelligent Laboratory Systems, 194, 103856.

    doi:10.1016/j.chemolab.2019.103856
  7. 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. 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