Description of Changes

Motivating Summary

This PR dramatically improves the speed of PLS regression by implementing the Dayal-MacGregor modified kernel algorithm (Dayal-MacGregor 1997). This has significant performance improvements over the current default, NIPALS, without sacrificing accuracy or numerical stability.

The Dayal-MacGregor algorithm accurately computes the weight and rotation matrices for the maximal covariance projection(s) of X without explicitly calculating the X score matrix and without deflating the Y target(s) at all. Thus it is (much) faster than NIPALS in the base case of univariate Y, and becomes asymptotically faster as the number of Y targets increases. Testing shows speed improvements of over 40% in the base case, and rising to over 95% when Y is multivariate.

In other words, this PR makes PLS regression 2x - 20x faster for equivalent results. The paper contains detailed numerical stability characteristics and accuracy results which are confirmed by my empirical tests below.

Specific Changes

• the parameter constraints on the algorithm parameter for the _PLS base class were modified so that a user can select the new algorithm via algorithm="dayalmacgregor" or algorithm="kernel" - see comments below, I'm not attached to having both arguments refer to the new algorithm, but I've implemented it this way for now
• the PLSRegression child class was modified to allow selection of the algorithm parameter during instantiation (previously a user could not select the algorithm, and the regression child class implicitly took the default of the base class). this argument defaults to algorithm="nipals" for backwards compatibility
• the fit method was augmented to include the implementation of Dayal-MacGregor, within a conditional branch on the value of the algorithm attribute
• Added a new test to cover the new algorithm

Performance & Accuracy Results

Small X (20 x 10) and Univariate Y (20 x 1) with 5 Components

NIPALS: 344 µs ± 4.77 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

Dayal-MacGregor: 198 µs ± 3.87 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

The output predictions agree to 6 decimal places.

43% faster

Small X (20 x 10) and Multivariate Y (20 x 2) with 5 Components

NIPALS: 700 µs ± 13.7 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

Dayal-MacGregor: 268 µs ± 4.92 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

The output predictions are equal to 4 decimal places.

62% faster

Moderate X (200 x 100) and Univariate Y (200 x 1) with 50 Components

NIPALS: 32.6 ms ± 2.16 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

Dayal-MacGregor: 15 ms ± 470 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

The output predictions agree to 6 decimal places.

54% faster

Moderate X (200 x 100) and Multivariate Y (200 x 2) with 50 Components

NIPALS: 226 ms ± 21.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Dayal-MacGregor: 15.8 ms ± 280 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

The output predictions are equal to 6 decimal places.

94% faster

Moderate X (200 x 100) and Multivariate Y (200 x 3) with 50 Components

NIPALS: 346 ms ± 29.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Dayal-MacGregor: 17 ms ± 1.48 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)

The output predictions are equal to 5 decimal places.

95% faster

Large X (2000 x 1000) and Univariate Y (2000 x 1) with 100 Components

NIPALS: 859 ms ± 27.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Dayal-MacGregor: 271 ms ± 20.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

The output predictions are equal to 6 decimal places.

59% faster

Large X (2000 x 1000) and Multivariate Y (2000 x 2) with 100 Components

NIPALS: 5.71 s ± 390 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Dayal-MacGregor: 233 ms ± 10.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

The output predictions are equal to 6 decimal places.

96% faster

Large X (2000 x 1000) and Multivariate Y (2000 x 3) with 100 Components

NIPALS: 6.43 s ± 398 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Dayal-MacGregor: 236 ms ± 13.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

The output predictions agree to 6 decimal places.

97% faster

Testing

Dev Environment & Architecture

My testing environment is as follows:

• Hardware: MacBook Pro (16-inc, 2021) with Apple M1 Max Chip and 64GB RAM
• OS: macOS Monterey v12.4
• Python: 3.10.5
• Installer/Package Manager: Anaconda3 (for Apple Silicon/arm64)
• Building Compiler: llvm-openmp

I don't think this will materially impact the results, but wanted to note it.

Testing Script

I put together a performance and accuracy testing script.

After running conda create -n sklearn-dev -c conda-forge python numpy scipy cython joblib threadpoolctl pytest compilers llvm-openmp, conda activate sklearn-dev and pip install --verbose --no-build-isolation --editable . in the directory, the following testing script should validate the foregoing performance and accuracy assessment. It randomly generates an m x n matrix X and m x n_targets matrix Y, and then fits PLS regression models estimated with each of the NIPALS and Dayal-MacGregor algorithms, then assesses performance and accuracy.

import numpy as np
from numpy.testing import assert_array_almost_equal
from sklearn.cross_decomposition import PLSRegression

#%% Config
m = 1000
n = 200
n_targets = 2

X = np.random.rand(m, n)
y = np.random.rand(m, n_targets)
n_components = 50
scale = True

####
# Performance Testing
####

#%% NIPALS
%%timeit
nipals = PLSRegression(n_components=n_components, scale=scale, algorithm='nipals')
nipals.fit(X=X, Y=y)
nipals.predict(X)

#%% DayalMacGregor Kernel
%%timeit
dayal = PLSRegression(n_components=n_components, scale=scale, algorithm='dayalmacgregor')
dayal.fit(X=X, Y=y)
dayal.predict(X)

#%%
####
# Accuracy/Stability Testing
####

control = PLSRegression(n_components=n_components, scale=scale, algorithm='nipals')
test = PLSRegression(n_components=n_components, scale=scale, algorithm='dayalmacgregor')
control.fit(X, y)
test.fit(X, y)
assert_array_almost_equal(test.predict(X), control.predict(X), decimal=4)

Comments

• I kept NIPALS as the default algorithm. If this PR is accepted, we might want to discuss whether or not NIPALS should remain the default given that Dayal-MacGregor is strictly superior for PLSRegression. Maybe keep NIPALS as the default, mention the new algorithm prominently in documentation and issue a warning that the new algorithm will become the default in a release or two?

• Given that Dayal-MacGregor doesn't deflate the Y target(s), there are no y_scores_, y_weights_ and y_rotations_ class attributes. This means that Dayal-MacGregor is only implemented for PLS Regression (not PLS Canonical), and a ValueError check is implemented to that effect. Moreover, the PLSRegression class documentation was edited to reflect the fact that these three attributes are only given for NIPALS/SVD.

• Should this algorithm be available via both algorithm="dayalmacgregor" and algorithm="kernel", or is it better to choose one? I currently have both implemented. I can see arguments either way but I'm not strongly opinionated. The Dayal-MacGregor algorithm is a kernel algorithm, but it's not the only one that exists in the literature (e.g. SIMPLS is also a kernel algorithm, albeit numerically unstable).