Improving Improved Kernel PLS Ole-Christian Galbo Engstrøma a FOSS Analytical A/S, Nils Foss Allé 1, Hillerød, 3400, Denmark
arXiv:2607.16138v1 [cs.LG] 17 Jul 2026
Abstract Improved Kernel Partial Least Squares (IKPLS) algorithms 1 and 2 are among the fastest PLS calibration algorithms. This article focuses on two shared steps, the computation of the X rotations, R, and the Y loadings, Q, and accelerates both. For R, term-by-term accumulation is replaced by a direct evaluation strategy that requires the same number of multiplications but parallelizes better on modern hardware. For Q, I identify — to the best of my knowledge, for the first time — equivalences showing that each Y loading is obtainable, up to explicitly derived constants, from quantities already computed earlier in the same iteration, and I exploit them in IKPLS to reduce the cost of each loading from Θ pKM q to Θ pM q operations whenever M “ 1 or 2 ď M ă K, with K predictor variables (number of columns in X) and M response variables (number of columns in Y). Both improvements provably yield exactly the same W, P, Q, R, and T as the original algorithms. Benchmarks with NumPy (CPU) and JAX (GPU) show speedups of up to two orders of magnitude for the isolated steps and of approximately 2ˆ (CPU) and 6ˆ (GPU) for entire fits. Both improvements are implemented in the free, open-source Python package ikpls. Keywords: partial least squares, computational complexity, algorithm design, chemometrics 1. Introduction Partial Least Squares (PLS, Wold 25) is a primary tool in chemometrics (Brereton et al. 7, Sørensen et al. 23). It is commonly used to analyze nearinfrared (NIR) spectra in both regression (Wold 25, Wold et al. 26) and classification (Sjöström et al. 22, Ståhle and Wold 24, Barker and Rayens 3, Brereton and Lloyd 8). Preprocessing is commonly applied to NIR spectra prior to PLS calibration. However, the best preprocessing is hard to determine ahead of modeling and thus typically requires experimentation to find (Rinnan et al. 20). Here, chemometricians should be concerned with the speed and stability of their chosen PLS algorithm. While the original NIPALS algorithm (Wold 25) is numerically stable (Andersson 2, Björck and Indahl 5), it is also slow (Alin 1). Improved Kernel PLS (IKPLS; Dayal and MacGregor
9) algorithms 1 and 2 are both fast PLS algorithms (Alin 1, Andersson 2). Algorithm 1 is also numerically stable (Andersson 2, Björck and Indahl 5), except on pathological data (Björck and Indahl 5), and algorithm 2 has been shown to produce equivalent results in practice (Engstrøm et al. 11). Faster computation was indeed a motivation for the development of IKPLS (Dayal and MacGregor 9). Thus, IKPLS is a strong choice of PLS algorithm for the chemometrician who requires both calibration precision and speed. In this article, I identify two steps in both IKPLS algorithms that can be further sped up. They concern the computations of R (X rotations) and Q (Y loadings), respectively. I prove the mathematical equivalence between the original IKPLS steps and the improvements introduced here, and I show that the computation of Q can be sped up by a factor of Θ pKq in most cases.
predictor variables matrix (N ˆ K). response variables matrix (N ˆ M ). PLS weights matrix for X (K ˆ A). PLS loadings matrix for X (K ˆ A). PLS loadings matrix for Y (M ˆ A). PLS weights matrix to compute scores T directly from original X (K ˆ A). T PLS scores matrix of X (N ˆ A). wa a column vector of W. pa a column vector of P. qa a column vector of Q. ra a column vector of R. ta a column vector of T. K number of X-variables. M number of Y-variables. N number of objects. A number of components in PLS model. a integer counter for latent variable dimension. This work will denote a matrix with subscript a as the state of the matrix at latent-variable dimension a. E.g., Pa is the PLS loadings matrix for X with the first a latent dimensions, and ` thus ˘ it has K rows and a columns. In contrast, XT Y a denotes the matrix product XT Y deflated with the first a ´ 1 components, ` and thus it ˘ has K `rows and M˘columns. Similarly, YT XXT Y a and XT YYT X a denote ` T ˘T ` T ˘ ` ˘ ` ˘T X Y a X Y a and XT Y a XT Y a , respectively. In addition to this notation, 0˚ denotes a length ˚ column vector of zeros, and δ is the Kronecker delta. Finally, PLS1 refers to the case when M “ 1 and PLS2 refers to the case when M ě 2, as is standard in the chemometric literature. X Y W P Q R
Although the improved computation of R has the same asymptotic runtime as the original computation, both improvements demonstrate improved runtime on datasets of varying sizes and shapes, illustrating their effect in the broader context of executing the full PLS algorithms. I have implemented the improved IKPLS algorithms with NumPy (Harris et al. 13) and JAX (Bradbury et al. 6) and released them in the free, open-source Python package ikpls (Engstrøm et al. 11), which also features sample-weighted PLS (Becker and Ismail 4, Engstrøm 10) and the combination of IKPLS with the fast cross-validation algorithms by Engstrøm and Jensen [12]. The rest of this article is organized as follows. Section 2 introduces the notation used and is consistent with the notation in the original IKPLS article (Dayal and MacGregor 9). Sections 3 and 4 introduce the two improvements to the IKPLS algorithms and prove both their correctness and, for the improvement in Section 4, superior asymptotic runtime. Section 5 briefly mentions additional efficient implementation details of IKPLS in ikpls (Engstrøm et al. 11). Section 6 exemplifies the faster practical runtime when applying the improvements introduced in Sections 3 and 4. Section 7 provides concluding remarks.
2. Nomenclature
3. Improved computation of R The original step 3 for computation of R in both IKPLS algorithms (Dayal and MacGregor 9) is due The notation in this work is identical to the one to Höskuldsson [14] and is given by the following reused by Dayal and MacGregor [9], so that readers currence relation. # familiar with that work may easily compare it to this wa a“1 work. For convenience, the nomenclature from Dayal (1) ra “ řa´1 ` T ˘ and MacGregor [9] is written below. wa ´ i“1 pi wa ri a ě 2 2
As written, Equation (1) invites evaluating the sum by term-by-term accumulation over the a ´ 1 preceding components, and this is how the implementation in the appendix of Dayal and MacGregor [9] proceeds. The accumulation is not necessary, however: de Jong [16, Equations 19–20], who attributes it to Höskuldsson [15], gives a direct computation of ra as a single matrix-vector product with an explicit K ˆ K matrix maintained across components by` a rank-one update. Their formulation requires ˘ Θ K 2 multiplications per component and an ad` ˘ ditional Θ K 2 memory, whereas Equation (1) requires just Θ pKpa ´ 1qq multiplications per component (Proposition 2). Since a ď K in PLS, and typically a ! K, this difference of a factor Θ pK{aq is substantial. Consider now an alternative, direct approach to computing ra that requires no such matrix but only the Ra´1 and Pa´1 that IKPLS accumulates in any case. It is given by
ra “
#
wa a“1 ˘ ` w a ě2 wa ´ Ra´1 PT a´1 a
fi pT 1 wa — pT ffi ` ˘ — 2 wa ffi w “ R Ra´1 PT ffi .. a´1 — a´1 a – fl . »
pT a´1 wa
“ “ r1 “ “
a´1 ÿ
i“1 a´1 ÿ i“1
r2
fi pT 1 wa T ffi ‰— — p2 w a ffi ¨ ¨ ¨ ra´1 — ffi .. – fl . »
pT a´1 wa
` ˘ ri pT i wa ` T ˘ pi wa ri
(3)
Thus, in the case of a ě 2, both Equation (1) and Equation (2) subtract the term in Equation (3) from wa to obtain ra , finalizing the proof that the two (2) equations are equivalent. 3.2. Runtime
The evaluation strategy of Equation (2) appears, without statement or proof, in the function kernelpls.fit of the R package pls [17, 19, version 2.9-0], an implementation of IKPLS algorithm 1, which uses it for a ą 5 and falls back to term-by-term accumulation for 2 ď a ď 5.
To analyze the computational complexity of computing r via Equation (1) and Equation (2), I ignore the summations and subtractions and focus solely on the number of multiplication operations required, as multiplication dominates both in number of operations and the practical constant of execution on modern hardware.
3.1. Correctness
Proposition 2. Equation (1) requires 2Kpa´1q multiplications to compute ra .
I will now prove the correctness of Equation (2).
Proof. When a “ 1, Equation (1) returns wa directly, requiring 2Kpa´1q “ 0 multiplications. When a ě 2, computing the inner vector product in the parentheses requires K multiplications, and so does the subsequent scalar vector multiplication, for a total of 2K multiplications. Summing over a ´ 1 products Proof. Recall from Dayal and ‰ MacGregor [9] means the total number of multiplications reaches “ and Pa “ 2Kpa ´ 1q. “ ‰ r1 r2 ¨ ¨ ¨ ra “that Ra p1 p2 ¨ ¨ ¨ pa . Then, in the case of a “ 1, the equations are identical. In the case of a ě 2, write Proposition 3. Equation (2) requires 2Kpa´1q mulout the matrix vector product to obtain tiplications to compute ra . Proposition 1 (Correctness, Equation (2)). Equation (1) and Equation (2) compute the same value of ra .
3
Proof. When a “ 1, Equation (2) returns wa directly, requiring 2Kpa´1q “ 0 multiplications. When a ě 2, the matrix-vector product in the parentheses must be evaluated first, as this yields the fewest total multiplications. PT a´1 wa is a multiplication between a matrix of size pa ´ 1q ˆ K and a vector of size K ˆ 1, Algorithm 1 Step 2 of IKPLS requiring K ˆ pa ´ 1q multiplications to yield a vec1: if M “ 1 then Ź PLS1 ` ˘ tor of size pa ´ 1q ˆ 1, which is then left-multiplied 2: w̃a Ð XT Y a by Ra´1 of size K ˆ pa ´ 1q requiring an additional 3: else Ź PLS2 K ˆ pa ´ 1q multiplications for a total of 2Kpa ´ 1q 4: if M ă K then Ź 2 ď M ă ˘K ` multiplications. 5: λa Ð largest eigenvalue of YT XXT Y a 6: q̃a Ð eigenvector corresponding to λa ` ˘ Proposition 4. Equation (1) and Equation (2) re7: w̃a Ð XT Y a q̃a quire the same number of multiplications to compute 8: else Ź 2 ď` M ^ K ď M ˘ ra . 9: λa Ð largest eigenvalue of XT YYT X a w̃a Ð eigenvector corresponding to λa Proof. This follows directly from Proposition 2 and 10: 11: end if Proposition 3. 12: end if Therefore, the gain from swapping the sequential 13: wa Ð w̃a kw̃a k2 Equation (1) for the direct Equation (2) is solely due to the matrix products being better suited for parallel execution on modern multi-processor hardware. 4. Improved computation of Q The improvement for computation of ra shown in Section 3 is always applicable when A ě 2. In contrast, the improvement in the computation of qa shown in this section applies only when M “ 1 (commonly referred to as PLS1) or 2 ď M ă K (some cases of PLS2). Thus, the improved computation of Algorithm 2 Step 4 of IKPLS qa does not apply when 2 ď M ^ K ď M , mean1: if algorithm “ 1 then Ź IKPLS algorithm 1 ing Y has at least two columns and at least as many 2: ta Ð Xra T columns as X (the remaining cases of PLS2). ta 3: pa Ð X kta k22 First, consider step 2 of the IKPLS algorithms as 4: else Ź IKPLS algorithm 2 given by Dayal and MacGregor [9], shown in Algo2 T T T T 5: kt k Ð r X Xr Ź rT a a a a X Xra “ ta ta 2 rithm 1. XT Xra 6: pa Ð kta k2 Then, step 3 consists of computing ra , which can 2 7: end if be done as shown in Section 3, and, in step 4, Dayal T T T pra pX Yqa q and MacGregor [9] compute pa and qa as shown in 8: qa Ð 2 kta k2 Algorithm 2. 2 Note that line 5 of Algorithm 2 produces kta k2 without computing ta . Substitute the definition of ta from line 2 into line 5 to see this. This allows reference to kta k22 regardless of whether IKPLS algorithm 1 or 2 is used. 4
Now, consider Algorithm 3, which is an improvement of Algorithm 2. The proportionality between qa and q̃a exploited by line 12 of Algorithm 3 is posed as a question— “is q proportional to q.a?”, with q denoting q̃a and q.a denoting qa — in a source comment in the function kernelpls.fit of the R package pls [17, 19, version 2.9-0]. However, it nevertheless computes qa as in line 8 of Algorithm 2. The comment has been present since version 2.0-0 (2006). Theorem 1 answers the question affirmatively and supplies the constant of proportionality. In the same function, the M “ 1 branch also computes qa as in line 8, although the divisor it uses to normalize w̃a into wa is exactly the numerator in line 9 of Algorithm 3.
` ˘ ` T ˘ 2 X Y a`1 “ XT Y a ´ pa qT a kta k2
Furthermore, it is assumed ` that ˘ the extraction of PLS components halts if XT Y a becomes a zero matrix or if ta “ 0N : in either case, all covariance between X and Y has been extracted, and extraction stops with a ´ 1 being the final component. The following lemmas can now be established. Lemma 1. PT a Ra “ I where I is the identity matrix with a rows and columns. Proof. Consider an index pi, jq into PT a Ra . It is given by ` T ˘ Pa Ra i,j “ pT i rj
Algorithm 3 Improved step 4 of IKPLS 1: if algorithm “ 1 then 2: ta Ð Xra
Ź IKPLS algorithm 1
T
ta pa Ð X kta k22 4: else Ź IKPLS algorithm 2 T T T 5: kta k22 Ð rT X Xr Ź rT a a a X Xra “ ta ta T a 6: pa Ð XktaXr k22 7: end if 8: if M “ 1 then Ź PLS1 w̃a k2 Ź k w̃ k from Algorithm 1. 9: qa Ð kkt 2 a 2 a k2 10: else Ź PLS2 11: if M ă K then Ź2ďM ăK ? 12: qa Ð kq̃a kλ2aktq̃aa k2 Ź λa and q̃a from 2 Algorithm 1. 13: else Ź 2ďM ^K ďM T T prT a pX Y qa q 14: qa Ð kta k22 15: end if 16: end if
3:
(4)
“
tT i Xrj kti k22
“
tT i tj kti k22
“
tT i tj tT i ti
“ δi,j . The last equality follows from the column vectors of T being mutually orthogonal (Höskuldsson 14), and kti k22 ą 0 holds for any extracted non-zero score, ti . Lemma 2. Consider an execution of IKPLS in which the qi consumed by each performed deflation ` ` T ˘ ˘T . Then at index i satisfies qi kti k22 “ rT i X Y i ` T ˘T X Y a ri “ 0M for every component a of the execution and all 1 ď i ď a ´ 1. Proof. Equivalently, the conclusion ` T by ˘ transposing, T states that rT i X Y a “ 0M for all 1 ď i ď a ´ 1. The supposition is in force throughout; the proof is by induction on a over the conclusion alone. Base step (a “ 1). The conclusion quantifies over the empty range 1 ď i ď 0 and is therefore vacuously true. Inductive step (a ě 2). Assume, as the inductive hypothesis, the conclusion for a ´ 1:
4.1. Correctness In this subsection, I will prove the correctness of Algorithm 3 by showing that it always computes the same qa as Algorithm 2. Consider step 5 of the IKPLS algorithms, as it is necessary to establish equivalence between ` Algo˘ rithm 2 and Algorithm 3. Step 5 deflates XT Y a and is given by Equation (4). 5
` T ˘ T rT i X Y a´1 “ 0M for all 1 ď i ď a ´ 2.
` ˘T ` T ˘T X Y 1 r1 “ XT Y 1 w1 .
(5)
For a ě 2, assume, as the inductive hypothesis, that the qi computed by Algorithm 3 satisfies To show the conclusion for a, fix i with 1 ď i ď Equation (7) at every index i with 1 ď i ď a ´ 1. T a ´ 1. Left-multiplying Equation (4) by ri yields The deflations performed before component a, at indices 1, . . . , a ´ 1, consumed precisely these qi , and ` T ˘ ` T ˘ T T 2 T p q kt k ´ r X Y “ r X Y rT multiplying Equation (7) by kti k22 ą 0 shows that a´1 a´1 i a´1 2 i i a´1 a ` ˘T ` ˘ T 2 T each satisfies qi kti k22 “ XT Y i ri . The execution “ rT ` T ˘ i X Y a´1 ´ δa´1,i qa´1 kta´1 k2 , (6) up to the formation of X Y a therefore satisfies the supposition of Lemma 2, whose conclusion yields ` ˘T where the second equality follows from Lemma 1. XT Y a ri “ 0M for all 1 ď i ď a ´ 1. Hence, Two cases remain. ` T ˘T If i ď a´2, then δa´1,i “ 0, and the remaining term X Y a ra is 0`T by the inductive hypothesis, Equation (5), so M ˘ ˘ ` ` ˘T ` T ˘T T T T ri X Y a “ 0M . “ X Y a wa ´ XT Y a Ra´1 PT a´1 wa If i “ a ´ 1, then δa´1,i “ 1, and the supposi` ˘T “ XT Y a wa , tion of Lemma 2, at `i “ a ´ 1 and transposed, gives ˘ T 2 T qT a´1 kta´1 k2 “ ra´1 X Y a´1 . Hence where the second equality holds by Lemma 2. Thus, for all a ě 1, ` T ˘ ` T ˘ ` T ˘ T T X Y ´ r X Y “ r X Y rT ` ˘T ` T ˘T a´1 a´1 a´1 a´1 a´1 a X Y a ra “ XT Y a wa . T “ 0M . Therefore, it suffices to show that each branch of ` T ˘ T Algorithm 3 reduces to “ 0 , completing the In both cases, rT X Y i M a ` T ˘T ` T ˘T induction. X Y a wa X Y a w̃a qa “ “ , (8) Theorem 1 (Correctness, Algorithm 3). The value kta k22 kw̃a k2 kta k22 of qa computed by Algorithm 3 equals the value of qa where the second equality comes from line 13 of computed by Algorithm 2. Algorithm 1. Under the halting assumptions stated (4), every extracted component has Proof. Algorithm 2 and Algorithm 3 differ only in below ` T Equation ˘ the computation of qa . Line 8 of Algorithm 2 always X Y a different from a zero matrix and ta ‰ 0N ; the former ensures kw̃a k2 ą 0 and the latter ensures computes qa as kta k22 ą 0, making Equation (8) well-defined. We ` T ` T ˘ ˘T ` T ˘T now consider the three cases: M “ 1, 2 ď M ă K, ra X Y a X Y a ra qa “ “ . (7) and 2 ď M ^ K ď M . 2 2 kta k2 kta k2 a) First, consider the case where M “ 1 (PLS1). I prove, by strong induction on a, that the qa Since ` ˘ M “ 1, line 2 of Algorithm 1 computes w̃a “ computed by Algorithm 3 also satisfies Equation (7) XT Y , so a for every a ě 1; the induction enters only through ` T ˘T Lemma 2, in relating ra to wa . Consider Equation (2) X Y a w̃a “ w̃aT w̃a “ kw̃a k22 , ` T ˘T for computing ra , and left-multiply it by X Y a to and Equation (8) becomes obtain the numerator in Equation (7). For a “ 1, Equation (2) contains no subtraction kw̃a k22 kw̃a k2 q “ “ , a term, so, with nothing assumed, kw̃a k2 kta k22 kta k22 6
which is line 9 of Algorithm 3. b) Now, consider the case where 2 ď M ă K (PLS2 with X larger than Y). Lines 5–7 of Algorithm 1 set λa to the largest ` ˘T ` ` ˘ ˘ eigenvalue of YT XXT Y a “ XT Y a XT Y a , set ` Tq̃a ˘to a corresponding eigenvector, and set w̃a “ X Y a q̃a . Then
Proof. Algorithm 2 and Algorithm 3 compute all quantities other than qa — in particular w̃a , wa , ra , pa , and kta k22 — identically, and qa influences later components only through the deflation in Equation (4). The run of IKPLS using Algorithm 3 therefore coincides with the run using Algorithm 2 component by component: ` ˘ the runs coincide up to the formation of XT Y 1 “ XT Y, and whenever they ` T ˘ ˘ ` ˘T ` ` T ˘T X Y a w̃a “ XT Y a XT Y a q̃a “ λa q̃a . (9) coincide up to the formation of X Y a , they compute identical w̃a , wa , ra , pa , and kta k22 , obtain iden` T ˘T ` T ˘ As X Y a X Y a is positive semi-definite, λa ě tical qa by Theorem ` T ˘1, and hence coincide up to the ` ˘ 0, and, for non-zero XT Y a , λa ą 0. Furthermore, formation of X Y a`1 through Equation (4). By induction on the components, the two runs coincide b and produce identical W, P, Q, R, and T. T T T kw̃a k2 “ q̃T a pX Yqa pX Yqa q̃a b b (10) λa kq̃a k22 “ λa q̃T a q̃a “ 4.2. Runtime a “ λa kq̃a k2 . In this subsection, I analyze the runtime of computing qa in Algorithms 2 and 3. I prove that the Substituting Equation (9) and Equation (10) into latter requires Θ pKq times fewer operations than the Equation (8) yields former, except when 2 ď M ^ K ď M , in which case ? the two coincide. λa q̃a λa q̃a “ , qa “ ? 2 2 kq̃a k2 kta k2 λa kq̃a k2 kta k2 Proposition 5. The computation of qa in line 8 Algorithm 2 requires Θ pKM q operations when which is line 12 of Algorithm 3. `of T ˘ c) Now, consider the final case where 2 ď M ^K ď X Y a has been precomputed. M (PLS2 with X smaller than Y). Line 14 of Algorithm 3 computes Proof. the vector-matrix product ` TComputing ˘ ` T ` T ˘ ˘T and its subsequent transposition rT a X Y a ra X Y a requires Θ pKM q operations. The subsequent qa Ð kta k22 element-wise division by kta k22 requires Θ pM q operations. Thus, the total number of required and is identical to line 8 of Algorithm 2. The nested else clauses partition the possibilities operations is Θ pKM q ` Θ pM q “ Θ pKM q. into M “ 1, 2 ď M ă K, and 2 ď M ^ K ď M , which are exhaustive and non-overlapping, and qa agrees with Algorithm 2 in all cases due to a), b), and c). This establishes Equation (7) at index a and completes the induction. Since line 8 of Algorithm 2 computes qa precisely by Equation (7), the qa computed by Algorithm 3 equals that of Algorithm 2 at every component a.
Proposition 6. The computation of qa in line 9 of Algorithm 3 requires ΘpM q “ Θ p1q operations. Proof. Scalar division between kw̃a k2 and kta k22 requires Θp1q operations. Since M “ 1 due to the satisfied condition in line 8 of Algorithm 3, the cost becomes ΘpM q “ Θ p1q.
Corollary 1. Executing IKPLS with Algorithm 3 in place of Algorithm 2 produces identical W, P, Q, R, Proposition 7. The computation of qa in line 12 of and T. Algorithm 3 requires ΘpM q operations. 7
? Proof. Computing λa requires Θp1q operations, since λa is already available from line 5 of Algorithm? 1. The subsequent vector-scalar product q̃a λa requires an additional Θ pM q operations. In the denominator, computing kq̃a k2 requires Θ pM q operations, and the subsequent scalar product kq̃a k2 kta k22 requires an additional Θ p1q operations. Finally, element-wise division of the numerator by the denominator requires Θ pM q operations. Thus, the total cost is 2Θp1q ` 3ΘpM q “ ΘpM q operations.
11) also optimizes the computation of Algorithm 1. These optimizations are not discussed here, as they rely on reusing precomputed matrix products and optimizing the order of multiplication in matrix products involving more than two matrices. For example, XT YYT X can be obtained quickly by using the preT T T computed ` TX Y ˘ and the equivalence X YY X “ T T X Y pX Yq , together with the fact that transposition is fast. These optimizations are not algorithmic improvements but efficient implementations. Therefore, they will not be discussed further, and the interested reader is referred to the source code of ikpls (Engstrøm et al. 11) which contains similar optimizations for many other steps in the IKPLS algorithms.
Proposition 8. The computation of qa in line 14 of Algorithm 3 requires Θ pKM q operations. Proof. This is identical to the proof of Proposition 5.
6. Benchmarks
Proposition 9. Algorithm 3 requires Θ pKq times fewer operations than Algorithm 2 to compute qa when M “ 1 (PLS1) or 2 ď M ă K (PLS2 with Y having fewer columns than X), and Algorithm 3 requires the same number of operations as Algorithm 2 when 2 ď M ^ K ď M (PLS2 with Y having at least as many columns as X).
This section presents benchmarks demonstrating how the improvements established in Sections 3 and 4 reduce the practical runtime of both IKPLS algorithms relative to the original formulation by Dayal and MacGregor [9]. All benchmarks are made using Python1 version 3.14 with the NumPy (Harris et al. 13, version 2.5.1) and JAX (Bradbury et al. 6, version 0.10.2, CUDA version 12.9) implementations of IKPLS in ikpls (Engstrøm et al. 11, version 6.1.2) with float64 precision. The NumPy implementations were executed on an AMD Ryzen 9 5950X processor utilizing all 16 cores (32 threads). The JAX implementations were executed on an Nvidia GeForce RTX 3090 Ti. For the JAX implementation, JIT-compiled variants were benchmarked, and compile-time overhead was not included. For both NumPy and JAX implementations, an initial few runs were executed and their results discarded to initialize the cache and ensure that the benchmarks executed first are not at a disadvantage due to a cold cache.
Proof. By Proposition 5, Algorithm 2 always requires Θ pKM q operations to compute qa . If M “ 1, then, by Proposition 6, Algorithm 3 requires Θ pM q operations to compute qa , which is K times fewer than Algorithm 2. If 2 ď M ă K, then, by Proposition 7, Algorithm 3 requires Θ pM q operations to compute qa , which is K times fewer than Algorithm 2. Otherwise, 2 ď M ^ K ď M and, then, by Proposition 8, Algorithm 3 requires Θ pKM q operations to compute qa , which is the same number as Algorithm 2.
6.1. Improvements in isolation The runtime of R is dependent on K and A as proven in Propositions 2 and 3. Therefore, Figure 1
5. Additional improvements of IKPLS In addition to the other improvements discussed in this article, the ikpls package (Engstrøm et al.
1 https://www.python.org/
8
shows practical runtimes for the computation of R with varying K and A. Likewise, the runtime of qa for any a is dependent only on K and M as proven by Propositions 5, 6, 7, and 8. Therefore, Figure 2 shows practical runtimes for the computation of Q with varying K and M . Each benchmark in Figure 1 and Figure 2 uses the same random seed and is the median of 30 to 1000 identical runs for each cell, to suppress variation from other factors. Faster cells receive more runs. Figure 1 reveals that the relative speedup increases with A and decreases with K but that the new evaluation strategy of R is always faster than the old one. My hypothesis is that as A grows, the original evaluation strategy has to do an increasing amount of sequential work that the new evaluation strategy can parallelize. And as K grows, the sequential steps themselves also grow, allowing better utilization of the parallel hardware within each step, thereby reducing the relative speedup. The relative speedups are generally most dramatic for the GPU, which is the more parallel of the two hardware types, corroborating the hypothesis. Figure 2 corroborates Proposition 5, Proposition 6, Proposition 7, and Proposition 8 nicely. For every value of K, the time required for computation of Q with Algorithm 3 is approximately constant when M “ 1 (« 8µs for the CPU and « 100µs for the GPU), and when 2 ď M ă K, the time consumption generally increases with M but is otherwise independent of K. When 2 ď M ^ K ď M , the time consumption increases with both K and M just like Algorithm 2 does for all values of K and M . For the GPU, when K is small, even though 2 ď M ^K ď M , the runtime does not seem to increase with increasing values of M . I hypothesize that this is an artifact of the GPU cores not being saturated at such smallsized operations. Sometimes, when 2 ď M ^ K ď M , Algorithm 3 seems slightly slower than Algorithm 2 when executed on a CPU. I attribute this to noise in the timing runs and, perhaps, to the fact that Algorithm 3 has to evaluate the conditionals on lines 8 and 11 before computing qa in the same way as Algorithm 2, which does not have to evaluate the conditionals. In practice, this tiny overhead is completely negligible. In
contrast, under JAX’s just-in-time compilation, M and K are static properties of the input shapes, so the conditionals on lines 8 and 11 of Algorithm 3 are resolved when the function is traced, and only the selected branch is compiled: the hot runs execute no conditionals, making Algorithm 2 and Algorithm 3 identical compiled programs when 2 ď M ^ K ď M . The dispatch that selects the compiled program from the input shapes at call time incurs the same cost for both algorithms, since compiled programs are specialized to the input shapes regardless of branching. 6.2. Improvements relative to full fit Figure 3 shows a representative subset of the practical speedup achieved by the improved R and Q computations in the context of entire IKPLS fits. The speedups are estimated by the medians and the first and third quartiles across 10 runs using the original and improved implementations of both IKPLS algorithms. For NumPy CPU implementations, the time spent on R and Q for the original and improved implementations is measured within the full fits. For the GPU runs, JAX compiles the entire fit routine into a fused XLA program, to the best of my knowledge, with no simple way to measure the wall-clock time spent on R and Q. So, the runtime for R and Q is not measured directly during the fit but is instead estimated using compiled standalone programs that compute the same operations required by R and Q, respectively. These isolated times are also each the median of 10 runs and include a roughly fixed CUDA kernel launch and JAX synchronization overhead that is independent of problem size and is not the in-fit contribution of R or Q; only their reduction is meaningful. The benchmarks are partitioned across N “ 103 as a surrogate for a reasonably sized dataset for PLS and N “ 200, which is relevant for either reasonably small datasets or local PLS modeling (Shenk et al. 21) on a larger dataset. In all cases, the fits were done with A “ 30 components. The regime with N ă K was not explored, as practitioners concerned with speed will, in these cases, likely prefer the parsimonious kernel PLS (Liland et al. 18) that computes the N ˆ N XXT matrix and operates with that, not unlike how IKPLS algorithm 2 operates on XT X. 9
62ms / 4.3ms
217ms / 15ms
11×
8.22×
14×
34×
15×
100
12ms / 363µs
50
3.1ms / 160µs
3.9ms / 427µs
15ms / 3.2ms
12×
7.54×
3.46×
16ms / 1.4ms
19×
30
1.1ms / 91µs
20
488µs / 58µs
9.24×
1.4ms / 189µs
8.46×
4.59×
4.47×
10
120µs / 26µs
5
30µs / 12µs
2
4.7µs / 3.9µs
4.7µs / 3.9µs
101
102
2.43×
1.22×
5.87×
617µs / 105µs
3.64×
122µs / 27µs
157µs / 43µs
2.42×
2.17×
30µs / 13µs
39µs / 18µs
1.21×
1.21×
K
60ms / 7.3ms
4.69×
10×
200
257ms / 2.8ms
456ms / 11ms
20×
5×
100
64ms / 970µs
47×
28×
10×
50
16ms / 529µs
17ms / 712µs
29ms / 1.4ms
17×
14×
17×
1×
3.06×
2.35×
403µs / 171µs
1.69×
101µs / 60µs
1.15×
17µs / 14µs
103
104
100×
1.6s / 11ms
1.6ms / 531µs
6.0µs / 5.0µs
54×
500
2×
3.6ms / 1.1ms
147×
20×
91×
66×
31×
30
6.2ms / 371µs
20
2.9ms / 242µs
23×
6.2ms / 438µs
12×
5.85×
745µs / 127µs
5
274µs / 97µs
2
112µs / 74µs
121µs / 84µs
101
102
2.83×
1.53×
10×
2.9ms / 282µs
6.43×
10
(a) R: speedup for K ˆ A (NumPy, CPU).
65ms / 1.4ms
5.06×
845µs / 131µs
850µs / 168µs
2.96×
2.87×
299µs / 101µs
301µs / 105µs
1.44×
1.48×
K
2.8s / 53ms
40×
115ms / 4.1ms
21×
11ms / 637µs
50×
5× 2× 1×
12×
4.9ms / 416µs
6.15×
1.3ms / 211µs
speedup (original / improved
200
1.3s / 127ms
A
10×
429ms / 22ms
speedup (original / improved
A
20×
500
2.87×
404µs / 141µs
1.52×
117µs / 79µs
128µs / 84µs
103
104
(b) R: speedup for K ˆ A (JAX, GPU).
Figure 1: Heatmap showing the speedup of the computation of R of shape K ˆA by using the evaluation strategy of Equation (2) as opposed to the term-by-term accumulation strategy of Equation (1). Since A is upper-bounded by minpN, Kq, the cells with K ă A are not computed.
101
100
1.00×
1.00×
14ms / 14ms
250ms / 250ms
2≤M∧K≤M
2≤M∧K≤M
2≤M∧K≤M
2≤M∧K≤M
0.99×
1.00×
1.00×
346×
35µs / 35µs
114µs / 115µs
473µs / 473µs
18ms / 52µs
2≤M∧K≤M
2≤M∧K≤M
2≤M∧K≤M
2≤M<K
0.99×
0.99×
8.24×
57×
21µs / 21µs
28µs / 28µs
131µs / 16µs
1.6ms / 28µs
2≤M∧K≤M
2≤M∧K≤M
2≤M<K
2≤M<K
0.99×
1.47×
2.48×
13×
20µs / 20µs
22µs / 15µs
37µs / 15µs
202µs / 15µs
2≤M∧K≤M
2≤M<K
2≤M<K
2≤M<K
1.83×
1.92×
1.97×
3.34×
15µs / 8.2µs
15µs / 8.0µs
17µs / 8.6µs
28µs / 8.3µs
M=1
M=1
M=1
M=1
101
102
103
104
K
300×
20× 10× 5× 2× 1×
1.00×
1.00×
1.00×
1.00×
102µs / 101µs
319µs / 319µs
1.2ms / 1.2ms
11ms / 11ms
2≤M∧K≤M
2≤M∧K≤M
2≤M∧K≤M
2≤M∧K≤M
1.00×
1.00×
1.00×
11×
103
101µs / 101µs
154µs / 154µs
350µs / 350µs
1.2ms / 116µs
2≤M∧K≤M
2≤M∧K≤M
2≤M∧K≤M
2≤M<K
1.00×
1.00×
1.46×
3.12×
102
102µs / 102µs
152µs / 152µs
149µs / 102µs
348µs / 112µs
2≤M∧K≤M
2≤M∧K≤M
2≤M<K
2≤M<K
1.00×
1.47×
1.39×
1.94×
101
101µs / 101µs
148µs / 101µs
140µs / 101µs
204µs / 105µs
2≤M∧K≤M
2≤M<K
2≤M<K
2≤M<K
1.02×
1.21×
1.21×
1.68×
101µs / 99µs
120µs / 99µs
120µs / 100µs
178µs / 106µs
M=1
M=1
M=1
M=1
101
102
103
104
104
100× 50×
100
(a) Q: speedup for K ˆ M, A “ 10 (NumPy, CPU).
K
10× 5×
2× 1×
speedup (origina / improved)
102
1.00× 403µs / 403µs
M
M
103
1.00× 160µs / 160µs
speedup (origina / improved)
104
(b) Q: speedup for K ˆ M, A “ 10 (JAX, GPU).
Figure 2: Heatmap showing the speedup of the computation of Q of shape M ˆ A with A “ 10 by using Algorithm 3 as opposed to Algorithm 2. Each qa takes time independent of a, and so using A “ 10 simply takes 10 times longer than computing any single qa .
Overall, the results are quite encouraging as the the common PLS1 case (M “ 1). But there are also improved versions are faster than the originals in ev- significant speedups achieved in the PLS2 case, most ery case. The relative speedup is particularly large for notably when 2 ď M ă K. Generally, it seems that 10
shape, they speed up entire PLS fits by up to « 2ˆ on a CPU and « 6ˆ on a GPU; and that, for realistic dataset sizes and component counts, the R improvement contributes the larger share of the speedup, although only the Q improvement lowers the asymptotic operation count. For IKPLS users and implementers, the results are twofold: Section 3 proves an evaluation strategy that is already used, without statement or proof, in the source code of the R package pls (Liland et al. 17, Mevik and Wehrens 19), and Theorem 1 answers, affirmatively, a question that its source has posed about PLS2 since 2006 while also covering the PLS1 case. As the replacement is drop-in, the improved computation of Q is an immediate candidate for adoption there and in other existing IKPLS im7. Conclusion plementations. Both improvements are implemented in the free, This article identified improvements to the computations of R (X rotations) and Q (Y loadings) open-source Python package ikpls [11, version in both IKPLS algorithms. The R improvement re- 6.1.2]. places term-by-term accumulation with a parallelfriendly evaluation strategy requiring the same number of multiplications and is applicable whenever A (the number of components) is at least 2. Correctness is proved in Proposition 1, and the multiplication counts are established in Propositions 2–4. The Q improvement (Algorithm 3) derives Y loadings from quantities already computed in step 2 of IKPLS and is applicable when M “ 1 (PLS1) or 2 ď M ă K (PLS2 with Y having fewer columns than X). In both cases, the improved computation requires a factor of Θ pKq fewer operations than the original computation. Correctness is proved in Theorem 1, Corollary 1 shows that the improvement is a drop-in replacement producing identical W, P, Q, R, and T, and the runtime analysis is given in Propositions 5–9. Finally, benchmarks showed that both improvements are realized in practice on both CPU and GPU implementations; that, depending on dataset size and the improvement in R is most significant, although it was solely a practical consideration, whereas the improvement in Q decreased the asymptotic runtime. This is the case both for the actual in-fit time measurements on the CPU and the estimated values on the GPU. It must be stated, however, that the relative speedups would be smaller for large N as none of the optimizations depend on N but other parts of the PLS fits do. Even when M becomes large, especially when M ě K, the relative speedups become vanishingly small as most of the time is then spent outside the computation of R and, if M ě K, Q falls back to the original implementation.
2 The observant reader may notice discrepancies between the relative speedups of R and Q in Figures 1 and 2 and those in Figure 3. This is likely because, within a full fit, the surrounding computation perturbs the cache residency and branch-predictor state under which R and Q execute.
11
A=30, N=103, K=103, (=1 [(=1]
A=30 N=103 K=103 M=1 [M=1]
R: 2.64ms 298µs Q: 98µs 48.5µs
original
R (isolated): 8.66ms 445µs Q (isolated): 282µs 189µs
original
→
→
→
improved 2.5 ms
A=30 N=103 K=103 M=10 [2
≤
0.0 s
M ( K] R: 2.59ms 297µs Q: 198µs 109µs
original
→
→
speedup (original / improved) 1.257× [1.180 1.360]
4.0 ms
A=30 N=103 K=103 M=103 [2
8.0 ms ≤
improved
improved
improved 1.5 ms
improved 3.0 ms
6.0 ms
A=30, N=2 102, K=102, (=10 [2
improved 1.5 ms
3.0 ms
R (isolated): 6.15ms 342µs Q (isolated): 319µs 194µs →
→
speedup (original / improved) 1.222 [1.222, 1.227]
improved 0.0 s
15.0 ms
30.0 ms
⋅
A=30, N=2 102, K=102, (=103 [2
R: 2.29ms → 282µs Q: 493µs → 495µs
original
R (isolated): 6.31ms → 334µs Q (isolated): 335µs → 336µs
speedup (original / improved) 1.053 [1.052, 1.053]
improved
40.0 ms
0.0 s
R
50.0 ms
100.0 ms
fit wall time
Q
Original ,KP./ fit
(a) IKPLS Algorithm 1, NumPy, CPU.
,mproved ,KP./ fit
(b) IKPLS Algorithm 1, JAX, GPU.
A=30 N=2 102 K=102 M=1 [M=1]
A=30, N=2 102, K=102, (=1
⋅
⋅
R: 1.21ms 107µs Q: 44µs 24µs
original
≤ M ∧ K ≤ M]
original
speedup (original / improved) 1.030× [1.026 1.044]
improved
M ) K]
original
≤ M ∧ K ≤ M]
T,e res) of ),e .K/L1 fi)
≤
⋅
→
→
N = 2 102
⋅
speedup (original / improved) 1.616× [1.614 1.627]
fi) wall )ime
→
→
speedup (original / improved) 5.541 [5.476, 5.590]
0.0 s
R: 1.29ms 117µs Q: 68µs 61.5µs
20.0 ms
R (isolated): 6.11ms 351µs Q (isolated): 256µs 189µs
original
→
M ( K]
⋅
N = 2 102
800.0 ms
→
3.0 ms
original
0.0 s
400.0 ms ⋅
164µs 39.5µs speedup (original / improved) 2.287× [2.160 2.295]
⋅
→
speedup (original / improved) 1.005 [1.005, 1.006]
A=30, N=2 102, K=102, (=1 [(=1]
R: 2.12ms Q: 77.5µs
A=30 N=2 102 K=102 M=103 [2
≤
R (isolated): 6.46ms 417µs Q (isolated): 911µs 913µs
0.0 s
original
0.0 s
∧
→
⋅
≤
M K M]
→
A=30 N=2 102 K=102 M=1 [M=1]
⋅
30.0 ms ≤
original
→
2.0 s
A=30 N=2 102 K=102 M=10 [2
15.0 ms
A=30, N=103, K=103, (=103 [2 1.3ms 16ms speedup (original / improved) 1.004× [0.999 1.008]
0.0 s
→
→
speedup (original / improved) 1.239 [1.232, 1.282]
0.0 s
R: 4.07ms Q: 15.3ms
1.0 s
R (isolated): 6.39ms 427µs Q (isolated): 309µs 192µs
original
≤
original
0.0 s
8.0 ms
M ) K]
≤
improved
M K M] ∧
4.0 ms
A=30, N=103, K=103, (=10 [2
improved 0.0 s
speedup (original / improved) 4.673 [4.606, 4.749]
improved
5.0 ms
N = 103
0.0 s
N = 103
→
speedup (original / improved) 1.587× [1.405 1.701]
(=1]
R (isolated): 6.19ms 361µs Q (isolated): 251µs 192µs
original
→
→
→
improved 800.0 µs ⋅
1.6 ms ≤
0.0 s
M ( K]
⋅
R: 1.24ms 112µs Q: 66µs 56.5µs
original
⋅
→
→
speedup (original / improved) 1.654× [1.648 1.658]
improved 0.0 s
1.0 ms
2.0 ms
⋅
A=30 N=2 102 K=102 M=103 [2
T,e res) of ),e .K/L1 fi)
→
→
improved 15.0 ms
⋅
≤ M ∧ K ≤ M] R (isolated): 6.15ms → 333µs Q (isolated): 329µs → 329µs
speedup (original / improved) 1.052× 1.051, 1.052]
improved
40.0 ms
0.0 s
R
30.0 ms
2
original
speedup (original / improved) 1.034× [1.018 1.049]
fi) wall )ime
M ) K]
speedup (original / improved) 1.222× 1.210, 1.233]
A=30, N=2 102, K=102, (=103
R: 2.31ms → 287µs Q: 493µs → 495µs
20.0 ms
≤
R (isolated): 6.12ms 325µs Q (isolated): 319µs 191µs
0.0 s
improved
6.0 ms
2
original
≤ M ∧ K ≤ M]
original
0.0 s
3.0 ms
A=30, N=2 102, K=102, (=10
⋅
A=30 N=2 102 K=102 M=10 [2
speedup (original / improved) 5.849× 5.782, 6.012]
improved
N = 2 102
0.0 s
N = 2 102
→
speedup (original / improved) 2.208× [2.193 2.221]
Q
50.0 ms
fit wall time
Original ,KP./ fit
(c) IKPLS Algorithm 2, NumPy, CPU.
100.0 ms
,mproved ,KP./ fit
(d) IKPLS Algorithm 2, JAX, GPU.
Figure 3: Time spent on full fits for the original and improved IKPLS algorithms. IKPLS algorithm 2 is only executed when N ą K.2
12
References
[9] Dayal, B.S., MacGregor, J.F., 1997. Improved PLS algorithms. Journal of Chemometrics 11, 73–85. doi:10.1002/(SICI)1099-128X(199701)11:1<73::AID-CEM435>
[1] Alin, A., 2009. Comparison of PLS algorithms when number of objects is much larger than number of variables. Statistical Papers 50, 711– [10] Engstrøm, O.C.G., 2025. Near-infrared hy720. doi:10.1007/s00362-009-0251-7. perspectral imaging applications in food analysis – improving algorithms and methodolo[2] Andersson, M., 2009. A comparison of nine gies. Ph.D. thesis. University of Copenhagen. PLS1 algorithms. Journal of Chemometrics 23, doi:10.48550/arXiv.2510.13452. 518–529. doi:10.1002/cem.1248.
[11] Engstrøm, O.C.G., Dreier, E.S., Jespersen, B.M., Pedersen, K.S., 2024. IKPLS: improved kernel partial least squares and fast crossvalidation algorithms for Python with CPU and GPU implementations using NumPy and [4] Becker, J.M., Ismail, I.R., 2016. AccountJAX. Journal of Open Source Software 9, 6533. ing for sampling weights in pls path moddoi:10.21105/joss.06533. eling: Simulations and empirical examples. European Management Journal 34, 606–617. [12] Engstrøm, O.C.G., Jensen, M.H., 2025. doi:10.1016/j.emj.2016.06.009. Fast partition-based cross-validation with
[3] Barker, M., Rayens, W., 2003. Partial least squares for discrimination. Journal of Chemometrics 17, 166–173. doi:10.1002/cem.785.
centering and scaling for XT X and XT Y. [5] Björck, Å., Indahl, U.G., 2017. Fast Journal of Chemometrics 39, e70008. and stable partial least squares modelling: doi:10.1002/cem.70008. A benchmark study with theoretical comments. Journal of Chemometrics 31, e2898. [13] Harris, C.R., Millman, K.J., van der Walt, S.J., doi:10.1002/cem.2898. e2898 cem.2898. Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N.J., [6] Bradbury, J., Frostig, R., Hawkins, P., JohnKern, R., Picus, M., Hoyer, S., van Kerkwijk, son, M.J., Leary, C., Maclaurin, D., Necula, G., M.H., Brett, M., Haldane, A., del Rı́o, J.F., Paszke, A., VanderPlas, J., Wanderman-Milne, Wiebe, M., Peterson, P., Gérard-Marchant, P., S., Zhang, Q., 2018. JAX: composable transSheppard, K., Reddy, T., Weckesser, W., Abformations of Python+NumPy programs. URL: basi, H., Gohlke, C., Oliphant, T.E., 2020. Array http://github.com/google/jax. programming with NumPy. Nature 585, 357– 362. doi:10.1038/s41586-020-2649-2. [7] Brereton, R.G., Jansen, J., Lopes, J., Marini, F., Pomerantsev, A., Rodionova, O., Roger, [14] Höskuldsson, A., 1988. PLS regression methJ.M., Walczak, B., Tauler, R., 2018. Chemoods. Journal of Chemometrics 2, 211–228. metrics in analytical chemistry—part II: moddoi:10.1002/cem.1180020306. eling, validation, and applications. Analytical The h-principle and Bioanalytical Chemistry 410, 6691–6704. [15] Höskuldsson, A., 1992. in modelling with applications to chemodoi:10.1007/s00216-018-1283-4. metrics. Chemometrics and Intelli[8] Brereton, R.G., Lloyd, G.R., 2014. Partial least gent Laboratory Systems 14, 139–153. squares discriminant analysis: taking the magic doi:10.1016/0169-7439(92)80099-P. proaway. Journal of Chemometrics 28, 213–225. ceedings of the 2nd Scandinavian Symposium doi:10.1002/cem.2609. on Chemometrics. 13
[16] de Jong, S., 1993. SIMPLS: An alterstudy. Journal of Chemometrics 1, 185–196. native approach to partial least squares doi:10.1002/cem.1180010306. regression. Chemometrics and Intelligent Laboratory Systems 18, 251–263. [25] Wold, H., 1966. Estimation of principal components and related models by iterative least doi:10.1016/0169-7439(93)85002-X. squares. Multivariate analysis , 391–420. [17] Liland, K.H., Mevik, B.H., Wehrens, R., 2026. pls: Partial least squares and [26] Wold, S., Sjöström, M., Eriksson, L., 2001. PLS-regression: a basic tool of principal component regression. URL: chemometrics. Chemometrics and Intelhttps://CRAN.R-project.org/package=pls. ligent Laboratory Systems 58, 109–130. r package version 2.9-0. doi:10.1016/S0169-7439(01)00155-1. [18] Liland, K.H., Stefansson, P., Indahl, U.G., 2020. Much faster cross-validation in PLSRmodelling by avoiding redundant calculations. Journal of Chemometrics 34, e3201. doi:10.1002/cem.3201. [19] Mevik, B.H., Wehrens, R., 2007. The PLS package: Principal component and partial least squares regression in r. Journal of Statistical Software 18, 1–23. doi:10.18637/jss.v018.i02. [20] Rinnan, Å., van den Berg, F., Engelsen, S.B., 2009. Review of the most common preprocessing techniques for near-infrared spectra. TrAC Trends in Analytical Chemistry 28, 1201– 1222. doi:10.1016/j.trac.2009.07.007. [21] Shenk, J.S., Westerhaus, M.O., Berzaghi, P., 1997. Investigation of a local calibration procedure for near infrared instruments. Journal of Near Infrared Spectroscopy 5, 223–232. doi:10.1255/jnirs.115. [22] Sjöström, M., Wold, S., Söderström, B., 1986. PLS discriminant plots. Elsevier. [23] Sørensen, K.M., van den Berg, F., Engelsen, S.B., 2021. NIR data exploration and regression by chemometrics—a primer. NearInfrared Spectroscopy: Theory, Spectral Analysis, Instrumentation, and Applications , 127– 189doi:10.1007/978-981-15-8648-4_7. [24] Ståhle, L., Wold, S., 1987. Partial least squares analysis with cross-validation for the two-class problem: a monte carlo 14