Cross-Validation Algorithms IKPLS, Improved Kernel Partial Least Squares;
NIPALS, Nonlinear Iterative Partial Least Squares; PLS, Partial Least Squares;
PLS-DA, Partial Least Squares Discriminant Analysis; PLS-R, Partial Least
Squares Regression; SG, Savitzky-Golay; SNV, Standard Normal Variate. Ole-
Christian Galbo Engstrøm, Department of Computer Science, University of
Copenhagen Faculty of Science, Universitetsparken 1, 2100 Copenhagen,
Denmmark<EMAIL_ADDRESS>R&D, Foss Analytical A/S, Nils Foss Allé 1, 3400
Hillerød, Denmark. FOSS Analytical A/S; Innovation Fund Denmark, Grant Number:
1044-00108B. Research Article

# Shortcutting Cross-Validation: Efficiently Deriving Column-Wise Centered and
Scaled Training Set $\mathbf{X}^{\mathbf{T}}\mathbf{X}$ and
$\mathbf{X}^{\mathbf{T}}\mathbf{Y}$ Without Full Recomputation of Matrix
Products or Statistical Moments

Ole-Christian Galbo Engstrøm  FOSS Analytical A/S, Denmark.  Department of
Computer Science (DIKU), University of Copenhagen, Denmark.  Department of
Food Science (UCPH FOOD), University of Copenhagen, Denmark.

###### Abstract

Cross-validation is a widely used technique for assessing the performance of
predictive models on unseen data. Many predictive models, such as Kernel-Based
Partial Least-Squares (PLS) models, require the computation of
$\mathbf{X}^{\mathbf{T}}\mathbf{X}$ and $\mathbf{X}^{\mathbf{T}}\mathbf{Y}$
using only training set samples from the input and output matrices,
$\mathbf{X}$ and $\mathbf{Y}$, respectively. In this work, we present three
algorithms that efficiently compute these matrices. The first one allows no
column-wise preprocessing. The second one allows column-wise centering around
the training set means. The third one allows column-wise centering and column-
wise scaling around the training set means and standard deviations.
Demonstrating correctness and superior computational complexity, they offer
significant cross-validation speedup compared with straight-forward cross-
validation and previous work on fast cross-validation - all without data
leakage. Their suitability for parallelization is highlighted with an open-
source Python implementation combining our algorithms with Improved Kernel
PLS.

Keywords — fast, cross-validation, cross validation, partial least squares,
autoscaling.

##  1 Introduction

In predictive modeling, the task is often to find a relation from $\mathbf{X}$
to $\mathbf{Y}$, where both are matrices with $N$ rows, each representing a
sample. Each sample is associated with $K$ predictors, making up the columns
of $\mathbf{X}$, and $M$ observations making up the columns of $\mathbf{Y}$.
Instead of directly using $\mathbf{X}$ and $\mathbf{Y}$, some models use one
or both of the matrix products, $\mathbf{X}^{\mathbf{T}}\mathbf{X}$ and
$\mathbf{X}^{\mathbf{T}}\mathbf{Y}$.

We present three novel algorithms for substantially speeding up cross-
validation requiring the computation of $\mathbf{X}^{\mathbf{T}}\mathbf{X}$
and $\mathbf{X}^{\mathbf{T}}\mathbf{Y}$ per dataset partition. $\mathbf{X}$
and $\mathbf{Y}$ are matrices with $N$ rows, each representing a sample.
$\mathbf{X}$ has $K$ columns representing predictors of the samples, and
$\mathbf{Y}$ has $M$ columns representing observations of the samples. Our
algorithms include the possibility of using the variants of
$\mathbf{X}^{\mathbf{T}}\mathbf{X}$ and $\mathbf{X}^{\mathbf{T}}\mathbf{Y}$
where $\mathbf{X}$ and $\mathbf{Y}$ have been mean-centered and standard
deviation-scaled on a per training set basis, avoiding data leakage from the
corresponding validation sets. Our algorithms do not require recomputing the
full $\mathbf{X}^{\mathbf{T}}\mathbf{X}$ and
$\mathbf{X}^{\mathbf{T}}\mathbf{Y}$ for each training set in the cross-
validation scheme, nor do they require recomputation of the statistical
moments.

Our algorithms can rapidly compute partition-wise, potentially mean-centered,
and optionally sample standard deviation-scaled
$\mathbf{X}^{\mathbf{T}}\mathbf{X}$ and $\mathbf{X}^{\mathbf{T}}\mathbf{Y}$.
These algorithms find valuable applications across various PLS-R [18, 19] and
PLS-DA [3] algorithms. Particularly noteworthy is their seamless integration
with the IKPLS algorithms [6], known for their speed [1] and numerical
stability [2]. Leveraging our cross-validation algorithms with IKPLS
facilitates swift testing of diverse preprocessing methods within a condensed
timeframe. Selecting the optimal preprocessing technique requires model
validation [14, 16] and is imperative for achieving peak performance [7].

Our algorithms have the same asymptotic runtime of $\Theta(NK(K+M))$ as the
cross-validation algorithm proposed in [13] in the simple case with no column-
wise preprocessing and column-wise mean-centering case. In the case of column-
wise standard deviation-scaling, our algorithm retains this runtime while the
one proposed in [13] increases in runtime. Additionally, we use the training
sets’ statistical moments to avoid data leakage from the validation sets where
[13] uses statistical moments computed on the entire dataset. Additionally,
the space complexity of our algorithms is asymptotically lower than [13] by a
factor $P$ where $P$ is the number of cross-validation partitions. Our
algorithms’ runtime and space complexity are entirely independent of the
number of cross-validation partitions. Furthermore, we show how to derive the
sample mean and standard deviation of the training sets from the sample mean
and standard deviation of the entire dataset without requiring full
recomputation.

We begin with a comparison between our own and related work in Section 2
followed by incremental introductions and analyses of our algorithms and the
naive alternatives in Section 3, Section 4, and Section 5, a reference to an
open-source parallel implementation of all our algorithms in Section 6, and a
conclusion of our contributions in Section 7.

##  2 Related Work

Several related works in efficient cross-validation exist. The most relevant,
perhaps, is the open source implementation [8] of all combinations of all our
algorithms and both IKPLS Algorithm 1 and IKPLS Algorithm 2 [6]. The open-
source implementation also handles multiprocessing for parallelization of the
cross-validation. It is elaborated upon in Section 6.

The primary, recurring idea behind all our algorithms is that matrix products
$\mathbf{X}^{\mathbf{T}}\mathbf{X}$ and $\mathbf{X}^{\mathbf{T}}\mathbf{Y}$
can be computed once using all samples. Then, for each cross-validation
partition, we remove the contribution of the validation set. The same idea
goes behind the computation of statistical moments. The idea behind computing
$\mathbf{X}^{\mathbf{T}}\mathbf{X}$ and $\mathbf{X}^{\mathbf{T}}\mathbf{Y}$
for the validation part only can be traced back to [12] with improvements by
[13]. Here, however, [13] uses the dataset-wide standard deviation for
scaling. Thus, if one wishes to use the approach presented by [13] to scale by
only the training set standard deviation, then recomputing the standard
deviation for each training set is required. However, this has devastating
effects on the time complexity, as evidenced by Proposition 5.10. Furthermore,
although the approach in [13] has many similarities with ours, the scaling by
the standard deviation is done on a per cross-validation partition basis by
matrix multiplication requiring an additional $\Theta(PK^{3})$ operations to
scale $\mathbf{X}^{\mathbf{T}}\mathbf{X}$ and $\Theta(PKM^{2})$ operations to
scale $\mathbf{X}^{\mathbf{T}}\mathbf{Y}$. The interested reader is referred
to equation (33) in [13] to see why this is true.

Like our algorithms, the algorithm in [17] computes
$\mathbf{X}^{\mathbf{T}}\mathbf{X}$ and $\mathbf{X}^{\mathbf{T}}\mathbf{Y}$
for the entire dataset, and then, for each cross-validation partition, they
subtract the contribution of the validation set, leaving them with the
training set. The difference is that they remove features (columns) of
$\mathbf{X}$ while we remove samples (rows) of $\mathbf{X}$ and $\mathbf{Y}$.
Although we drew inspiration from [17], our algorithm is vastly different.

In the same spirit as ours, [11] shows how to compute training set-wise
$\mathbf{X}\mathbf{X}^{\mathbf{T}}$ from the dataset-wise
$\mathbf{X}\mathbf{X}^{\mathbf{T}}$. They also show how to derive it based on
the centered $\mathbf{X}$. Their cross-validation algorithm is helpful for PLS
algorithms that use $\mathbf{X}\mathbf{X}^{\mathbf{T}}$. However, as datasets
grow ever larger, $\mathbf{X}\mathbf{X}^{\mathbf{T}}$ will be much larger than
$\mathbf{X}^{\mathbf{T}}\mathbf{X}$, the size of the first being quadratic in
the dataset size and the size of the latter being completely independent of
the dataset size. Additionally, [11] states that if scaling by the standard
deviation is required, then the dataset-wise standard deviation must be used,
or full recomputation must be performed. As previously stated, our algorithms
handle scaling by the training set-wise standard deviation without fully
recomputing matrix products or statistical moments.

In Table 1, we show a comparison of time and space complexities for algorithms
computing cross-validation partitions of $\mathbf{X}^{\mathbf{T}}\mathbf{X}$
and $\mathbf{X}^{\mathbf{T}}\mathbf{Y}$ where the partitioning is based on
sample selection. The comparison includes our algorithms, the algorithms
proposed in [13], and straight-forward, naive approaches for cross-validation
as described in [5, 9].

Work |  No column-wise preprocessing |  Column-wise centering |  Column-wise scaling
---|---|---|---
[5, 9] (naive cross-validation) |  Time: $\Theta(PNK(K+M))$ Space: $\Theta\mathbf{(}(K+N)(K+M))$ |  Time: $\Theta\mathbf{(}PNK(K+M))$ Space: $\Theta((K+N)(K+M))$ |  Time: $\Theta(PNK(K+M))$ Space: $\Theta((K+N)(K+M))$
[13]* |  Time: $\Theta(NK(K+M))$ Space: $\Theta\mathbf{(}(PK+N)(K+M))$ |  Time: $\Theta\mathbf{(}NK(K+M))$ Space: $\Theta((PK+N)(K+M))$ |  Time: $\Theta\mathbf{(}NK(K+M)+PK^{3}+PKM^{2})$ Space: $\Theta((PK+N)(K+M))$
Ours |  Time: $\bm{\Theta}\mathbf{(NK(K+M))}$ Space: $\bm{\Theta}\mathbf{((K+N)(K+M))}$ |  Time: $\bm{\Theta}\mathbf{(NK(K+M))}$ Space: $\bm{\Theta}\mathbf{((K+N)(K+M))}$ |  Time: $\bm{\Theta}\mathbf{(NK(K+M))}$ Space: $\bm{\Theta}\mathbf{((K+N)(K+M))}$
Table 1: Comparison of time and space complexities as functions of the input
sizes for algorithms computing $\mathbf{X}^{\mathbf{T}}\mathbf{X}$ and
$\mathbf{X}^{\mathbf{T}}\mathbf{Y}$ for each training set across all cross-
validation partitions. [5] is the straight-forward cross-validation algorithm
where the explicit computation of all training set matrices and statistical
moments must occur. *[13] use the dataset statistics when centering around the
mean and scaling around the standard deviation. Using dataset statistics
introduces data leakage between training and validation sets. Our algorithms
use the statistical moments from the training sets only, preventing data
leakage.

As a final note in this section, we wish to state that other preprocessing
methods operating independently on each sample can trivially be applied
directly to any of the algorithms computing
$\mathbf{X}^{\mathbf{T}}\mathbf{X}$ and $\mathbf{X}^{\mathbf{T}}\mathbf{Y}$
for each cross-validation training set, including our algorithms, without
changing their correctness. Such preprocessing, however, would have to be
applied prior to the potential column-wise centering and scaling. Common
examples of preprocessing methods that operate on each sample independently
include SNV [4], Detrend [4], and SG filters [15].

##  3 Cross-Validation Algorithm Without Column-Wise Preprocessing

In all our analyses, we assume for simplicity that the matrix multiplication
operation is the naive one that computes the dot product of each row vector of
the left matrix with each column vector of the right matrix - i.e., the
iterative algorithm that requires $ABC$ multiplications to multiply two
matrices with dimensions $A\times B$ and $B\times C$. We know faster
algorithms exist [10] \- especially for symmetric matrices such as
$\mathbf{X}^{\mathbf{T}}\mathbf{X}$, where only the upper or lower triangular
matrix needs to be computed. While the matrix product involves both
multiplication and addition, the former dominates the latter in terms of
amount and the practical constant required for multiplication. Moreover,
hardware supporting efficient fused multiply–add (such as GPUs) collapse
multiplication and addition into a single operation. Therefore, we distinguish
between the types of operations in our analyses.

An entry $(i,j)$ in a matrix $\mathbf{M}$ is denoted $\mathbf{M}_{ij}$. Let
$R=\\{1,\ldots,N\\}$ be the set of indices (rows) and $R^{\prime}\subseteq R$
then we denote by $\mathbf{X}_{R^{\prime}}$ ($\mathbf{Y}_{R^{\prime}}$) the
submatrix of $\mathbf{X}$ ($\mathbf{Y}$) containing each index (row) in
$R^{\prime}$. We write $\mathbf{X}^{\mathbf{T}}_{R^{\prime}}$ to mean
$(\mathbf{X}_{R^{\prime}})^{\mathbf{T}}$ and omit parentheses for brevity.

######  Definition 3.1 (partition, training set, validation set, cross-
validation).

Let $2\leq P\leq N$ be the number of cross-validation partitions. Each index
$n\in R$ is associated with a (cross-validation) _partition_
$p\in\\{1,\ldots,P\\}$, indicating which partition $n$ belongs to. The
_validation set_ $V_{p}\subseteq R$ denotes the set of indices of partition
$p$, and $T_{p}=R\setminus V_{p}$ the _training set_ of partition $p$. Observe
that $T_{p}\cap V_{p}=\emptyset$ and $T_{p}\cup V_{p}=R$. A _cross-validation_
is the fitting of models using $\mathbf{X}_{T_{p}}$ and $\mathbf{Y}_{T_{p}}$
and evaluating their performance using $\mathbf{X}_{V_{p}}$ and
$\mathbf{Y}_{V_{p}}$, for each partition $p\in\\{1,\ldots,P\\}$.

Observe that this definition gives us the property that
$\sum_{p=1}^{P}|V_{p}|=N$ since all $V_{p}$ are pairwise disjoint. Commonly,
cross-validation is identified by the number of partitions $P$ and named
$P$-fold cross-validation, with $P=N$ being leave-one-out cross-validation.
Note that we do not assume the partitions are balanced (containing roughly the
same amount of samples). We consider Algorithm 1 the immediate method for
computing cross-validation as given in Definition 3.1, and present Algorithm 2
as the more efficient alternative. Efficiency is achieved by computing matrix
products only for the validation set rather than the larger training set for
each partition.

Algorithm 1 Naive Cross-Validation Algorithm

1:Determine $T_{p}$ and $V_{p}$ for each partition $p\in\\{1,\ldots,P\\}$ as
per Definition 3.1.

2:for each partition $p\in\\{1,\ldots,P\\}$ do

3: Compute $\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{X}_{T_{p}}$ and
$\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{Y}_{T_{p}}$.

4:end for

Algorithm 2 Cross-Validation Algorithm

1:Determine $T_{p}$ and $V_{p}$ for each partition $p\in\\{1,\ldots,P\\}$ as
per Definition 3.1.

2:Compute $\mathbf{X}^{\mathbf{T}}\mathbf{X}$ and
$\mathbf{X}^{\mathbf{T}}\mathbf{Y}$.

3:for each partition $p\in\\{1,\ldots,P\\}$ do

4: Compute $\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{X}_{V_{p}}$ and
$\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{Y}_{V_{p}}$.

5: Compute
$\mathbf{X}^{\mathbf{T}}\mathbf{X}-\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{X}_{V_{p}}$
to obtain $\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{X}_{T_{p}}$, and
$\mathbf{X}^{\mathbf{T}}\mathbf{Y}-\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{Y}_{V_{p}}$
to obtain $\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{Y}_{T_{p}}$.

6:end for

###  3.1 Correctness

A useful relationship between dataset-wise $\mathbf{X}^{\mathbf{T}}\mathbf{X}$
($\mathbf{X}^{\mathbf{T}}\mathbf{Y}$) and submatrices formed by training and
validation sets $\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{X}_{T_{p}}$ and
$\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{X}_{V_{p}}$
($\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{Y}_{T_{p}}$ and
$\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{Y}_{V_{p}}$) is given by Lemma 3.1
(Lemma 3.2).

######  Lemma 3.1.

Let $p$ be a partition, $T_{p}$ a training set and $V_{p}$ a validation set,
then
$\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{X}_{T_{p}}=\mathbf{X}^{\mathbf{T}}\mathbf{X}-\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{X}_{V_{p}}$.

###### Proof.

We show the equivalent statement
$\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{X}_{T_{p}}+\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{X}_{V_{p}}=\mathbf{X}^{\mathbf{T}}\mathbf{X}$.
Recall that each of $\mathbf{X}^{\mathbf{T}}\mathbf{X}$,
$\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{X}_{T_{p}}$ and
$\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{X}_{V_{p}}$ are $K\times K$ matrices,
and consider any entry $(i,j)$ in $\mathbf{X}^{\mathbf{T}}\mathbf{X}$ for
$i\in\\{1,\ldots,K\\}$ and $j\in\\{1,\ldots,K\\}$ given by,

$\left(\mathbf{X}^{\mathbf{T}}\mathbf{X}\right)_{ij}=\sum_{n\in
R}\mathbf{X}^{\mathbf{T}}_{in}\mathbf{X}_{nj}$ (1)

Similarly, we have

$\left(\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{X}_{T_{p}}\right)_{ij}=\sum_{n\in
T_{p}}\mathbf{X}^{\mathbf{T}}_{in}\mathbf{X}_{nj}\ \ \text{ and }\ \
\left(\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{X}_{V_{p}}\right)_{ij}=\sum_{n\in
V_{p}}\mathbf{X}^{\mathbf{T}}_{in}\mathbf{X}_{nj}$ (2)

for $\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{X}_{T_{p}}$ respectively
$\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{X}_{V_{p}}$. By assumption $T_{p}\cap
V_{p}=\emptyset$ and $T_{p}\cup V_{p}=R$ so with commutativity of addition we
have

$\left(\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{X}_{T_{p}}\right)_{ij}+\left(\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{X}_{V_{p}}\right)_{ij}=\sum_{n\in
T_{p}}\mathbf{X}^{\mathbf{T}}_{in}\mathbf{X}_{nj}+\sum_{n\in
V_{p}}\mathbf{X}^{\mathbf{T}}_{in}\mathbf{X}_{nj}=\sum_{n\in T_{p}\cup
V_{p}}\mathbf{X}^{\mathbf{T}}_{in}\mathbf{X}_{nj}=\sum_{n\in
R}\mathbf{X}^{\mathbf{T}}_{in}\mathbf{X}_{nj}=\left(\mathbf{X}^{\mathbf{T}}\mathbf{X}\right)_{ij}$
(3)

∎

######  Lemma 3.2.

Let $p$ be a partition, $T_{p}$ a training set and $V_{p}$ a validation set,
then
$\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{Y}_{T_{p}}=\mathbf{X}^{\mathbf{T}}\mathbf{Y}-\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{Y}_{V_{p}}$.

###### Proof.

As proof of Lemma 3.1 with $j\in\\{1,\ldots,M\\}$. ∎

######  Proposition 3.3 (Correctness of Algorithm 2).

$\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{X}_{V_{p}}$ and
$\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{Y}_{V_{p}}$ computed in step 5 of
Algorithm 2 are identical to those computed in step 3 of Algorithm 1.

###### Proof.

Due to the selection in step 1, it follows from Lemma 3.1 and Lemma 3.2 that
$\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{X}_{T_{p}}$ and
$\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{Y}_{T_{p}}$ obtained in step 5 of
Algorithm 2 is identical to those in step 3 of Algorithm 1. ∎

###  3.2 Computational Complexity

We analyze the computational complexity of the algorithms using $O$-notation
and $\Theta$-notation.

######  Proposition 3.4.

Algorithm 1 requires a total of $\Theta(PNK(K+M))$ multiplication operations.

###### Proof.

For some $p$, step 3 requires $\Theta(|T_{p}|K(K+M))$ operations to compute
$\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{X}_{V_{p}}$ and
$\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{Y}_{V_{p}}$. Iterating over partitions
in step 2 implies all matrix products require
$\Theta\left(\sum_{p=1}^{P}|T_{p}|K(K+M)\right)=\Theta\left(K(K+M)\cdot\sum_{p=1}^{P}|T_{p}|\right)$
operations. By definition $|T_{p}|=N-|V_{p}|$ and

$\sum_{p=1}^{P}|V_{p}|=N$ (4)

so we have,

$\sum_{p=1}^{P}|T_{p}|=\sum_{p=1}^{P}N-|V_{p}|=\sum_{p=1}^{P}N-\sum_{p=1}^{P}|V_{p}|=PN-N=N(P-1)$
(5)

showing a total of $\Theta(PNK(K+M))$ operations are required. ∎

######  Proposition 3.5.

Algorithm 2 requires a total of $\Theta(NK(K+M))$ multiplication operations.

###### Proof.

Step 2 computes $\mathbf{X}^{\mathbf{T}}\mathbf{X}$ with iterative matrix
multiplication requiring $\Theta(KNK)=\Theta(NK^{2})$ operations and
$\mathbf{X}^{\mathbf{T}}\mathbf{Y}$ requiring $\Theta(KNM)$ operations, for a
total of $\Theta(NK(K+M))$ operations. For some $p$, step 4 requires
$\Theta\left(|V_{p}|K(K+M)\right)$ operations to compute
$\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{X}_{V_{p}}$ and
$\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{Y}_{V_{p}}$. Iterating over partitions
in step 3 implies all matrix products require
$\Theta\left(\sum_{p=1}^{P}|V_{p}|K(K+M)\right)=\Theta\left(K(K+M)\cdot\sum_{p=1}^{P}|V_{p}|\right)$
operations. Since all $V_{p}$ are pairwise disjoint we have
$\sum_{p=1}^{P}|V_{p}|=N$ and so steps 2 and 4 require a total of
$\Theta(NK(K+M))$ multiplication operations. ∎

######  Proposition 3.6.

Algorithm 2 requires $\Theta(P)$ fewer multiplication operations than
Algorithm 1.

###### Proof.

By Proposition 3.4 and Proposition 3.5 the ratio of multiplication operations
is $\frac{PNK(K+M)}{NK(K+M)}=\Theta(P)$ ∎

However, we must consider that there is an additional cost to Algorithm 2 from
the required matrix subtractions.

######  Proposition 3.7.

Algorithm 2 requires a total of $\Theta(PK(K+M))=O(NK(K+M))$ subtraction
operations.

###### Proof.

Matrix subtraction operations are determined by the dimensions of the matrices
involved, which are independent of the validation set sizes of a particular
partition. Due to step 3 we perform $P$ such subtractions in step 5, which
requires a total of $\Theta(P(K^{2}+KM))=\Theta(PK(K+M))$ operations which,
due to $P\leq N$, is $O(NK(K+M))$. ∎

Thus, Algorithm 1 requires $\Theta(PNK(K+M))$ total operations, and Algorithm
2 requires $\Theta(NK(K+M))+O(NK(K+M))=\Theta(NK(K+M))$ total operations.
Therefore, Algorithm 2 is faster than Algorithm 1, even when it comes to total
operations.

The space complexities of Algorithm 1 and Algorithm 2 are analyzed in
Proposition 3.9 and Proposition 3.8, respectively. A comparison between the
two follows.

######  Proposition 3.8.

Algorithm 1 requires storing $NK+NM+N+K^{2}+KM=\Theta((K+N)(K+M))$ matrix
entries.

###### Proof.

Storing $\mathbf{X}$ and $\mathbf{Y}$ requires $NK$ and $NM$ entries,
respectively. Both must be stored to compute $\mathbf{X}_{T_{p}}$ and
$\mathbf{Y}_{T_{p}}$ in step 3. $T_{p}$ and $V_{p}$ can be stored in a single
vector of length $N$, identifying the partition for which any given sample is
part of the validation set. Additionally, storing
$\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{X}_{T_{p}}$ and
$\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{Y}_{T_{p}}$ in step 3 require $K^{2}$
and $KM$ entries each. Summing all these requirements yields a total space
complexity of

$NK+NM+N+K^{2}+KM=\Theta((K+N)(K+M))$ (6)

∎

######  Proposition 3.9.

Algorithm 2 requires storage of $\Theta((K+N)(K+M))$ matrix entries.

###### Proof.

In addition to storing the same matrices as Algorithm 1, Algorithm 2 stores
$\mathbf{X}^{\mathbf{T}}\mathbf{X}$, $\mathbf{X}^{\mathbf{T}}\mathbf{Y}$,
$\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{X}_{T_{p}}$, and
$\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{Y}_{V_{p}}$. The latter two can share
memory with $\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{X}_{T_{p}}$, and
$\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{Y}_{T_{p}}$. Thus, the additional
storage requirements are $K^{2}+KM$, yielding a total space complexity of

$\Theta((K+N)(K+M))+K^{2}+KM=\Theta((K+N)(K+M))$ (7)

∎

Thus, Algorithm 1 and Algorithm 2 have identical space complexities.

As a final note in this section, we wish to address the algorithm proposed in
[12] and its successor [13]. They propose to, instead of precomputing
$\mathbf{X}^{\mathbf{T}}\mathbf{X}$ and $\mathbf{X}^{\mathbf{T}}\mathbf{Y}$,
independently compute $\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{X}_{V_{p}}$ and
$\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{Y}_{V_{p}}$ for all
$p\in\\{1,\ldots,P\\}$ requiring a total of $NK(K+M)$ multiplications. Then we
can add these matrix products together for each partition,
$p\in\\{1,\ldots,P\\}$, to form
$\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{X}_{T_{p}}$ and
$\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{Y}_{T_{p}}$. This modification would
speed up Algorithm 2 by a constant factor of $2$ \- meaning asymptotically no
change. However, it would require us to explicitly store each
$\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{X}_{V_{p}}$ and
$\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{Y}_{V_{p}}$. This additional storage
requirement increases the storage complexity to
$NK+NM+PK^{2}+PKM=(PK+N)(K+M)=\Theta((PK+N)(K+M))$. Therefore, we opted not to
pursue this modification of Algorithm 2.

##  4 Cross-Validation Algorithm With Column-Wise Centering

Algorithm 2 can be extended to account for training set-wise, column-wise
mean-centering of $\mathbf{X}_{T_{p}}$ and $\mathbf{Y}_{T_{p}}$ directly in
$\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{X}_{T_{p}}$ and
$\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{Y}_{T_{p}}$. We show how to compute
the mean-centering only once for the entire dataset and then, for each
partition, $p$, remove the contribution of any sample $n\in V_{p}$.

Let $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{R^{\prime}}}$ be the mean row of $\mathbf{X}_{R^{\prime}}$,
defined by

$\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{R^{\prime}}}=\frac{1}{|R^{\prime}|}\sum_{n\in
R^{\prime}}\mathbf{X}_{n}.$ (8)

Equivalently, let $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}$ be the mean row of $\mathbf{Y}_{R^{\prime}}$,
defined by

$\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}=\frac{1}{|R^{\prime}|}\sum_{n\in
R^{\prime}}\mathbf{Y}_{n}.$ (9)

Stacking these row vectors to obtain matrices with similar dimensions as
$\mathbf{X}_{R^{\prime}}$ and $\mathbf{Y}_{R^{\prime}}$, we get

$\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}_{R^{\prime}}}=\begin{bmatrix}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{R^{\prime}}}\\\ \vdots\\\ \macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{R^{\prime}}}\end{bmatrix}$ (10)

and

$\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{R^{\prime}}}=\begin{bmatrix}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}\\\ \vdots\\\ \macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}\end{bmatrix},$ (11)

where the dimensions of $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}_{R^{\prime}}}$ and $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{R^{\prime}}}$ are $|R^{\prime}|\times K$ and
$|R^{\prime}|\times M$, respectively. We take $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}^{\mathbf{T}}_{R^{\prime}}}$ and $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}$ to mean $\left(\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}_{R^{\prime}}}\right)^{\mathbf{T}}$ and $\left(\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{R^{\prime}}}\right)^{\mathbf{T}}$, respectively, and omit
parentheses for brevity.

Then, the mean centering of $\mathbf{X}_{R^{\prime}}$ and
$\mathbf{Y}_{R^{\prime}}$ can be defined by
$\mathbf{X}_{R^{\prime}\mathbf{c}}$ and $\mathbf{Y}_{R^{\prime}\mathbf{c}}$ as

$\mathbf{X}_{R^{\prime}\mathbf{c}}=\mathbf{X}_{R^{\prime}}-\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}_{R^{\prime}}}$ (12)

and

$\mathbf{Y}_{R^{\prime}\mathbf{c}}=\mathbf{Y}_{R^{\prime}}-\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{R^{\prime}}},$ (13)

respectively.

Now, we present Algorithm 3, which we consider the immediate extension of the
naive Algorithm 1 to account for column-wise mean-centering. Additionally, we
present Algorithm 4, which extends Algorithm 2 to include mean-centering of
$\mathbf{X}_{T_{p}}$ and $\mathbf{Y}_{T_{p}}$ without requiring full
recomputation of training set matrix products or mean vectors.

Algorithm 3 Naive Cross-Validation Algorithm with Mean-Centering

1:Determine $T_{p}$ and $V_{p}$ for each partition $p\in\\{1,\ldots,P\\}$ as
per Definition 3.1.

2:for each partition $p\in\\{1,\ldots,P\\}$ do

3: Compute $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}_{T_{p}}}$ and $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{T_{p}}}$.

4: Compute $\mathbf{X}_{T_{p}}-\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}_{T_{p}}}$ to obtain $\mathbf{X}_{T_{p}\mathbf{c}}$ and
$\mathbf{Y}_{T_{p}}-\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{T_{p}}}$ to obtain $\mathbf{Y}_{T_{p}\mathbf{c}}$.

5: Compute
$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{c}}\mathbf{X}_{T_{p}\mathbf{c}}$ and
$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{c}}\mathbf{Y}_{T_{p}\mathbf{c}}$.

6:end for

Algorithm 4 Cross-Validation Algorithm with Mean-Centering

1:Determine $T_{p}$ and $V_{p}$ for each partition $p\in\\{1,\ldots,P\\}$ as
per Definition 3.1.

2:Compute $\mathbf{X}^{\mathbf{T}}\mathbf{X}$,
$\mathbf{X}^{\mathbf{T}}\mathbf{Y}$, $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}}$, and $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}}$.

3:for each partition $p\in\\{1,\ldots,P\\}$ do

4: Compute $\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{X}_{V_{p}}$ and
$\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{Y}_{V_{p}}$.

5: Compute
$\mathbf{X}^{\mathbf{T}}\mathbf{X}-\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{X}_{V_{p}}$
to obtain $\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{X}_{T_{p}}$, and
$\mathbf{X}^{\mathbf{T}}\mathbf{Y}-\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{Y}_{V_{p}}$
to obtain $\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{Y}_{T_{p}}$.

6: Compute $\frac{|V_{p}|}{|T_{p}|}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{V_{p}}}$, $\frac{|V_{p}|}{|T_{p}|}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{V_{p}}}$, $\frac{N}{|T_{p}|}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}}$, and $\frac{N}{|T_{p}|}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}}$.

7: Compute $\frac{N}{|T_{p}|}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}}-\frac{|V_{p}|}{|T_{p}|}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{V_{p}}}$ to obtain $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{T_{p}}}$ and $\frac{N}{|T_{p}|}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}}-\frac{|V_{p}|}{|T_{p}|}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{V_{p}}}$ to obtain $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{T_{p}}}$.

8: Compute $|T_{p}|\left(\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{T_{p}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{T_{p}}}\right)$ and $|T_{p}|\left(\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{T_{p}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{T_{p}}}\right)$.

9: Compute
$\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{X}_{T_{p}}-|T_{p}|\left(\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{T_{p}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{T_{p}}}\right)$ to obtain
$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{c}}\mathbf{X}_{T_{p}\mathbf{c}}$ and
$\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{Y}_{T_{p}}-|T_{p}|\left(\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{T_{p}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{T_{p}}}\right)$ to obtain
$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{c}}\mathbf{Y}_{T_{p}\mathbf{c}}$.

10:end for

###  4.1 Correctness

Proposition 4.1, Proposition 4.2, Proposition 4.3, and Proposition 4.4 yield
equalities that show us how to compute the mean-centered
$\mathbf{X}^{\mathbf{T}}_{R^{\prime}\mathbf{c}}\mathbf{X}_{R^{\prime}\mathbf{c}}$
and
$\mathbf{X}^{\mathbf{T}}_{R^{\prime}\mathbf{c}}\mathbf{Y}_{R^{\prime}\mathbf{c}}$
given $\mathbf{X}^{\mathbf{T}}_{R^{\prime}}\mathbf{X}_{R^{\prime}}$ and
$\mathbf{X}^{\mathbf{T}}_{R^{\prime}}\mathbf{Y}_{R^{\prime}}$ with the only
multiplication requirement being an outer vector product. Proposition 4.5 and
Proposition 4.6 show how to efficiently compute the training set means
$\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{T_{p}}}$ and $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{T_{p}}}$ from the dataset means $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}}$ and $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}}$. We use these propositions to prove the correctness of
Algorithm 4 in Proposition 4.7.

######  Proposition 4.1.

The matrix products $\mathbf{X}^{\mathbf{T}}_{R^{\prime}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{R^{\prime}}}$, $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}^{\mathbf{T}}_{R^{\prime}}}\mathbf{Y}_{R^{\prime}}$ and
$\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}^{\mathbf{T}}_{R^{\prime}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{R^{\prime}}}$ are all equal to $|R^{\prime}|\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}$. That is,
$\mathbf{X}^{\mathbf{T}}_{R^{\prime}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{R^{\prime}}}=\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}^{\mathbf{T}}_{R^{\prime}}}\mathbf{Y}_{R^{\prime}}=\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}^{\mathbf{T}}_{R^{\prime}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{R^{\prime}}}=|R^{\prime}|\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}$

###### Proof.

We start by deriving $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}^{\mathbf{T}}_{R^{\prime}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{R^{\prime}}}=|R^{\prime}|\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}$ for this proof. Then, we build upon this
equality to derive the other two equalities.

Consider any entry $(i,j)$ into $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}^{\mathbf{T}}_{R^{\prime}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{R^{\prime}}}$. It is defined by

$\left(\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}^{\mathbf{T}}_{R^{\prime}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{R^{\prime}}}\right)_{ij}=\sum_{n\in
R^{\prime}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}^{\mathbf{T}}_{R^{\prime}}}_{{}_{{}_{{}_{in}}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{R^{\prime}}}_{{}_{nj}}.$ (14)

By definitions in Equation (10) and Equation (11), any column, $n$, in
$\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}^{\mathbf{T}}_{R^{\prime}}}$ is $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}$ and any row, $n$, in
$\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{R^{\prime}}}$ is $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}$. With these definitions and using the
distributivity of multiplication over addition, the right-hand side simplifies
to

$\sum_{n\in R^{\prime}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}_{{}_{{}_{{}_{i}}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}_{{}_{j}}=\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}_{{}_{{}_{{}_{i}}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}_{{}_{j}}\sum_{n\in
R^{\prime}}1=|R^{\prime}|\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}_{{}_{{}_{{}_{i}}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}_{{}_{j}}.$ (15)

This is the definition of $|R^{\prime}|\left(\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}\right)_{ij}$, which concludes the proof of
$\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}^{\mathbf{T}}_{R^{\prime}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{R^{\prime}}}=|R^{\prime}|\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}$.

Next, consider expanding the right-hand side of Equation (15) by using the
definition of $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}$ in Equation (8) and commutativity
of scalar multiplication to get

$|R^{\prime}|\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}_{{}_{{}_{{}_{i}}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}_{{}_{j}}=|R^{\prime}|\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}_{{}_{j}}\frac{1}{|R^{\prime}|}\sum_{n\in
R^{\prime}}\mathbf{X}^{\mathbf{T}}_{R^{\prime}_{in}}=\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}_{{}_{j}}\sum_{n\in
R^{\prime}}\mathbf{X}^{\mathbf{T}}_{R^{\prime}_{in}}$ (16)

We can put $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}_{{}_{j}}$ inside the sum using the distributivity
of multiplication over addition and then, using Equation (11) to see that
$\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}=\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{R^{\prime}}}_{{}_{n}}$ for any row, $n$, in $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{R^{\prime}}}$, we get

$\sum_{n\in
R^{\prime}}\mathbf{X}^{\mathbf{T}}_{R^{\prime}_{in}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}_{{}_{j}}=\sum_{n\in
R^{\prime}}\mathbf{X}^{\mathbf{T}}_{R^{\prime}_{in}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{R^{\prime}}}_{{}_{nj}}.$ (17)

This is the definition of
$\left(\mathbf{X}^{\mathbf{T}}_{R^{\prime}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{R^{\prime}}}\right)_{ij}$, proving that
$\mathbf{X}^{\mathbf{T}}_{R^{\prime}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{R^{\prime}}}=\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}$.

Now, we only need to prove equality to $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}^{\mathbf{T}}_{R^{\prime}}}\mathbf{Y}_{R^{\prime}}$. To this
end, go back to Equation (15) and expand the right-hand side by using the
definition of $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}_{{}_{j}}$ in Equation (9) to get

$|R^{\prime}|\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}_{{}_{{}_{{}_{i}}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}_{{}_{j}}=|R^{\prime}|\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}_{{}_{{}_{{}_{i}}}}\frac{1}{|R^{\prime}|}\sum_{n\in
R^{\prime}}\mathbf{Y}_{R^{\prime}_{nj}}=\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}_{{}_{{}_{{}_{i}}}}\sum_{n\in
R^{\prime}}\mathbf{Y}_{R^{\prime}_{nj}}.$ (18)

From the definition in Equation (10) that $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{R^{\prime}}}=\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}_{R^{\prime}}}_{n}$ for any row, $n$, in $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}_{R^{\prime}}}$, we have have the equivalent statement that any
column, $n$, in $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}^{\mathbf{T}}_{R^{\prime}}}$ is equal to $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}$. Using this and the distributivity
of multiplication over addition, we can get

$\sum_{n\in R^{\prime}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}_{{}_{{}_{{}_{i}}}}\mathbf{Y}_{R^{\prime}_{nj}}=\sum_{n\in
R^{\prime}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}^{\mathbf{T}}_{R^{\prime}}}_{{}_{{}_{{}_{in}}}}\mathbf{Y}_{R^{\prime}_{nj}},$
(19)

which is the definition of $\left(\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}^{\mathbf{T}}_{R^{\prime}}}\mathbf{Y}_{R^{\prime}}\right)_{ij}$,
proving that $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}^{\mathbf{T}}_{R^{\prime}}}\mathbf{Y}_{R^{\prime}}=|R^{\prime}|\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}$ which was the only remaining equality to prove.
∎

Note that $|R^{\prime}|\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}=|R^{\prime}|\left(\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}\right)$. When executed using floating point
precision on computers, the latter is more numerically stable. The more stable
approach comes at the price of scaling each entry in the $K\times M$ matrix
$\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}$ instead of just the vector with the fewest
entries of $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}$ or $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}$. We opt for the numerically stable approach,
and, as we will see later, this does not affect the asymptotic runtime of
Equation (4) in terms of number of multiplications required.

######  Proposition 4.2.

The matrix products $\mathbf{X}^{\mathbf{T}}_{R^{\prime}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}_{R^{\prime}}}$, $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}^{\mathbf{T}}_{R^{\prime}}}\mathbf{X}_{R^{\prime}}$ and
$\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}^{\mathbf{T}}_{R^{\prime}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}_{R^{\prime}}}$ are all equal to $|R^{\prime}|\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{R^{\prime}}}$. That is,
$\mathbf{X}^{\mathbf{T}}_{R^{\prime}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}_{R^{\prime}}}=\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}^{\mathbf{T}}_{R^{\prime}}}\mathbf{X}_{R^{\prime}}=\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}^{\mathbf{T}}_{R^{\prime}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}_{R^{\prime}}}=|R^{\prime}|\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{R^{\prime}}}$

###### Proof.

As of proof of Proposition 4.1 with
$\mathbf{Y}_{R^{\prime}}=\mathbf{X}_{R^{\prime}}$. ∎

######  Proposition 4.3.

Let $\mathbf{X}_{R^{\prime}\mathbf{c}}$ and
$\mathbf{Y}_{R^{\prime}\mathbf{c}}$ be defined by Equation (12) and Equation
(13), respectively. Then
$\mathbf{X}^{\mathbf{T}}_{R^{\prime}\mathbf{c}}\mathbf{Y}_{R^{\prime}\mathbf{c}}=\mathbf{X}^{\mathbf{T}}_{R^{\prime}}\mathbf{Y}_{R^{\prime}}-|R^{\prime}|\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}$

###### Proof.

Using the definitions in Equation (12) and Equation (13), we have

$\mathbf{X}^{\mathbf{T}}_{R^{\prime}\mathbf{c}}\mathbf{Y}_{R^{\prime}\mathbf{c}}=\left(\mathbf{X}_{R^{\prime}}-\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}_{R^{\prime}}}\right)^{\mathbf{T}}\left(\mathbf{Y}_{R^{\prime}}-\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{R^{\prime}}}\right).$ (20)

Expanding the right-hand side and using Proposition 4.1, this simplifies to

$\mathbf{X}^{\mathbf{T}}_{R^{\prime}}\mathbf{Y}_{R^{\prime}}-\mathbf{X}^{\mathbf{T}}_{R^{\prime}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{R^{\prime}}}-\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}^{\mathbf{T}}_{R^{\prime}}}\mathbf{Y}_{R^{\prime}}+\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}^{\mathbf{T}}_{R^{\prime}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{R^{\prime}}}=\mathbf{X}^{\mathbf{T}}_{R^{\prime}}\mathbf{Y}_{R^{\prime}}-|R^{\prime}|\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}.$ (21)

∎

######  Proposition 4.4.

Let $\mathbf{X}_{R^{\prime}\mathbf{c}}$ be defined by Equation (12). Then
$\mathbf{X}^{\mathbf{T}}_{R^{\prime}\mathbf{c}}\mathbf{X}_{R^{\prime}\mathbf{c}}=\mathbf{X}^{\mathbf{T}}_{R^{\prime}}\mathbf{X}_{R^{\prime}}-|R^{\prime}|\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{R^{\prime}}}$

###### Proof.

As of proof of Proposition 4.3 with
$\mathbf{Y}_{R^{\prime}}=\mathbf{X}_{R^{\prime}}$ and using proof of
Proposition 4.2 instead of proof of Proposition 4.1. ∎

######  Proposition 4.5.

Given $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}}$, we can compute the mean row of $\mathbf{X}_{T_{p}}$ as
$\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{T_{p}}}=\frac{N}{|T_{p}|}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}}-\frac{|V_{p}|}{|T_{p}|}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{V_{p}}}$ where $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{V_{p}}}$ is the mean row of $\mathbf{X}_{V_{p}}$.

###### Proof.

Using the definitions of $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}}$ and $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{V_{p}}}$ from Equation (8), the right-hand side can be written
as

$\frac{N}{|T_{p}|}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}}-\frac{|V_{p}|}{|T_{p}|}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{V_{p}}}=\frac{N}{|T_{p}|}\frac{1}{N}\sum_{n\in
R}\mathbf{X}_{n}-\frac{|V_{p}|}{|T_{p}|}\frac{1}{|V_{p}|}\sum_{n\in
V_{p}}\mathbf{X}_{n}=\frac{1}{|T_{p}|}\sum_{n\in
R}\mathbf{X}_{n}-\frac{1}{|T_{p}|}\sum_{n\in V_{p}}\mathbf{X}_{n}$ (22)

By Definition 3.1 and distribution of multiplication over addition, we can
expand the right-hand side to get

$\frac{1}{|T_{p}|}\sum_{n\in T_{p}}\mathbf{X}_{n}+\frac{1}{|T_{p}|}\sum_{n\in
V_{p}}\mathbf{X}_{n}-\frac{1}{|T_{p}|}\sum_{n\in
V_{p}}\mathbf{X}_{n}=\frac{1}{|T_{p}|}\sum_{n\in T_{p}}\mathbf{X}_{n},$ (23)

which is the definition of $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{T_{p}}}$

∎

######  Proposition 4.6.

Given $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}}$, we can compute the mean row of $\mathbf{Y}_{T_{p}}$ as
$\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{T_{p}}}=\frac{N}{|T_{p}|}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}}-\frac{|V_{p}|}{|T_{p}|}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{V_{p}}}$.

###### Proof.

As of proof of Proposition 4.5 with $\mathbf{X}_{T_{p}}=\mathbf{Y}_{T_{p}}$. ∎

######  Proposition 4.7 (Correctness of Algorithm 4).

$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{c}}\mathbf{X}_{T_{p}\mathbf{c}}$ and
$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{c}}\mathbf{Y}_{T_{p}\mathbf{c}}$
computed in step 9 of Algorithm 4 are identical to those computed in step 5 of
Algorithm 3.

###### Proof.

Due to the selection in step 1, correctness of
$\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{X}_{T_{p}}$ and
$\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{Y}_{T_{p}}$ obtained from step 5 of
Algorithm 4 has already been proven in Proposition 3.3. Correctness of
$\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{T_{p}}}$ and $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{T_{p}}}$ computed in step 7 of Algorithm 4 follows from
Proposition 4.5 and Proposition 4.6, respectively. Finally, Proposition 4.4
and Proposition 4.3 ensure that
$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{c}}\mathbf{X}_{T_{p}\mathbf{c}}$ and
$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{c}}\mathbf{Y}_{T_{p}\mathbf{c}}$
computed in step 9 of Algorithm 4 are identical to those obtained in step 5 of
Algorithm 3. ∎

###  4.2 Computational Complexity

Here, we analyze the runtime and space complexities of Algorithm 3 and
Algorithm 4. As noted previously, we ignore the additions required for matrix
multiplication as the multiplication operations dominate in terms of quantity
and practical constant required for computation. However, we account for
additions, which are not part of matrix multiplication operations.

######  Proposition 4.8.

Algorithm 3 requires a total of $\Theta(PN(K+M))$ additions.

###### Proof.

For some $p$, step 3 requires $\Theta\left(|T_{p}|(K+M)\right)$ additions to
compute $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{T_{p}}}$ and $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{T_{p}}}$ and thus construct $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}_{T_{p}}}$ and $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{T_{p}}}$. Iterating over partitions in step 2, we use Equation
(5) to derive a total of $\Theta(PN(K+M))$ additions required. ∎

######  Proposition 4.9.

Algorithm 3 requires a total of $\Theta(PN(K+M))$ subtractions.

###### Proof.

For some $p$, step 4 requires $\Theta\left(|T_{p}|(K+M)\right)$ subtractions
to compute $\mathbf{X}_{T_{p}\mathbf{c}}$ and $\mathbf{Y}_{T_{p}\mathbf{c}}$.
Iterating over partitions in step 2, we use Equation (5) to derive a total of
$\Theta(PN(K+M))$ subtractions required. ∎

######  Proposition 4.10.

Algorithm 3 requires a total of $\Theta(PNK(K+M))$ multiplication operations.

###### Proof.

For any partition, $p$, the matrix multiplications in step 5 of Algorithm 3
require the same amount of multiplications as those in step 3 of Algorithm 1.
Thus, the derivation in Proposition 3.4 also applies here yielding
$\Theta(PNK(K+M))$. For any partition, $p$, step 3 of Algorithm 3 requires
$K+M$ multiplications. Thus, the total number of multiplications is

$\Theta(PNK(K+M))+\Theta(P(K+M))=\Theta(PNK(K+M)).$ (24)

∎

######  Proposition 4.11.

Algorithm 3 requires a total of $\Theta(P)$ division operations.

###### Proof.

For any partition, $p$, the computation of $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}_{T_{p}}}$ and $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{T_{p}}}$ require $2$ divisions. Iterating over all partitions,
we get a total of $2P=\Theta(P)$ division operations. ∎

######  Proposition 4.12.

Algorithm 4 requires a total of $\Theta(N(K+M))$ additions.

###### Proof.

Step 2 requires $2N(K+M)=\Theta(N(K+M))$ additions to compute
$\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}}$ and $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}}$.

Step 6 requires $\Theta\left(|V_{p}|(K+M)\right)$ additions for each
partition. So, with Equation (4), step 6 requires a total of $\Theta(N(K+M))$
additions over all partitions.

Thus, the total amount of addition operations are
$\Theta(N(K+M))+\Theta(N(K+M))=\Theta(N(K+M))$ ∎

######  Proposition 4.13.

Algorithm 4 requires a total of $\Theta(PK(K+M))=O(NK(K+M))$ subtractions.

###### Proof.

It follows from Proposition 3.7 that, iterating over all partitions, step 5
requires $\Theta(PK(K+M))=O(NK(K+M))$ subtractions.

Step 7 requires $K+M$ subtractions for each partition for a total of
$\Theta(P(K+M))=O(N(K+M))$ subtractions over all partitions.

Step 9 requires $K^{2}+KM$ subtractions for each partition for a total of
$\Theta(PK(K+M))=O(NK(K+M))$ subtractions over all partitions.

Thus, the total amount of subtraction operations are

$\Theta(PK(K+M))+\Theta(P(K+M))+\Theta(PK(K+M))=\Theta(PK(K+M))=O(NK(K+M))$
(25)

∎

######  Proposition 4.14.

Algorithm 4 requires a total of $\Theta(NK(K+M))$ multiplication operations.

###### Proof.

Step 2 requires $\Theta(NK(K+M))$ multiplications to compute
$\mathbf{X}^{\mathbf{T}}\mathbf{X}$ and $\mathbf{X}^{\mathbf{T}}\mathbf{Y}$.
It requires an additional $K+M$ multiplications to compute $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}}$ and $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}}$ for a total of $\Theta(NK(K+M))$ multiplications in step 2.

It follows from Proposition 3.5 that, iterating over all partitions, step 4
requires $\Theta(NK(K+M))$ multiplication operations.

Step 6 requires $2K+2M$ multiplications for each partition for
$\Theta(PKM)=O(NKM)$ multiplications over all partitions.

Step 8 requires $2K^{2}+2KM$ multiplications for each partition for a total of
$\Theta(PK(K+M))=O(NK(K+M))$ multiplications over all partitions.

The total number of multiplications is then

$\Theta(NK(K+M))+\Theta(NK(K+M))+O(NKM)+O(NK(K+M))=\Theta(NK(K+M)),$ (26)

showing a total of $\Theta(NK(K+M))$ multiplications over all partitions. ∎

######  Proposition 4.15.

Algorithm 4 requires a total of $\Theta(P)=O(N)$ division operations.

###### Proof.

Step 2 requires $2$ division operations to compute $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}}$ and $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}}$.

For any partition, $p$, step 6 requires $6$ divisions. Iterating over all
partitions, we get a total of $6P$ division operations. So we have
$2+6P=\Theta(P)=O(N)$ division operations. ∎

It follows that Algorithm 3 requires

$\Theta(PN(K+M))+\Theta(PK(K+M))+\Theta(PNK(K+M))+\Theta(P)=\Theta(PNK(K+M))$
(27)

total operations. It also follows that Algorithm 4 requires

$\Theta(N(K+M))+O(NK(K+M))+\Theta(NK(K+M))+O(N)=\Theta(NK(K+M))$ (28)

total operations. Therefore, Algorithm 4 is faster than Algorithm 3.

######  Proposition 4.16.

Algorithm 3 requires storing $\Theta((K+N)(K+M))$ matrix entries.

###### Proof.

Storing $\mathbf{X}$ and $\mathbf{Y}$ requires $N(K+M)$ entries. Storing all
$V_{p}$ and $T_{p}$ requires $N$ entries. In step 3, storing $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}_{T_{p}}}$ and $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{T_{p}}}$ can be accomplished by only storing $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{T_{p}}}$ and $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{T_{p}}}$, requiring $K+M$ entries.. In step 4, storing
$\mathbf{X}_{T_{p}\mathbf{c}}$ and $\mathbf{Y}_{T_{p}\mathbf{c}}$ requires
$|T_{p}|(K+M)=O(N(K+M))$ entries. In step 5, storing
$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{c}}\mathbf{X}_{T_{p}\mathbf{c}}$
requires $K(K+M)$ entries. Thus, the total storage requirement is

$N(K+M)+N+K+M+O(N(K+M))+K(K+M)=\Theta((K+N)(K+M)).$ (29)

∎

######  Proposition 4.17.

Algorithm 4 requires storing $\Theta((K+N)(K+M))$ matrix entries.

###### Proof.

Storing $\mathbf{X}$ and $\mathbf{Y}$ requires $N(K+M)$ entries. Storing all
$V_{p}$ and $T_{p}$ requires $N$ entries. In step 2, storing of
$\mathbf{X}^{\mathbf{T}}\mathbf{X}$, $\mathbf{X}^{\mathbf{T}}\mathbf{Y}$,
$\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{x}$,
and $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a 111{y}$
requires $K(K+M)+K+M$ entries. In step 4, storing
$\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{X}_{V_{p}}$ and
$\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{Y}_{V_{p}}$ requires $K(K+M)$ entries.
In step 5, $\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{X}_{T_{p}}$ and
$\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{Y}_{V_{p}}$ can overwrite the location
of $\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{X}_{V_{p}}$ and
$\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{Y}_{V_{p}}$ for no additional storage
cost. In step 6, storing $\frac{|V_{p}|}{|T_{p}|}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{V_{p}}}$, $\frac{|V_{p}|}{|T_{p}|}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{V_{p}}}$, $\frac{N}{|T_{p}|}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}}$, and $\frac{N}{|T_{p}|}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}}$ require $2K+2M$ entries. In step 7, $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{T_{p}}}$ and $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{T_{p}}}$ can overwrite the locations of
$\frac{|V_{p}|}{|T_{p}|}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{V_{p}}}$ and $\frac{|V_{p}|}{|T_{p}|}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{V_{p}}}$ for no additional storage requirements. In step 8
storing $|T_{p}|\left(\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{T_{p}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{T_{p}}}\right)$ and $|T_{p}|\left(\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{T_{p}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{T_{p}}}\right)$ require $K(K+M)$ entries. In step 9,
$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{c}}\mathbf{X}_{T_{p}\mathbf{c}}$ and
$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{c}}\mathbf{Y}_{T_{p}\mathbf{c}}$ can
overwrite the locations of $\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{X}_{T_{p}}$
and $\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{Y}_{T_{p}}$ for no additional
storage requirements. Thus, the total storage requirement is

$N(K+M)+N+3K(K+M)+3K+3M=\Theta((K+N)(K+M))$ (30)

∎

Thus, Algorithm 3 and Algorithm 4 have identical asymptotic space
complexities.

##  5 Cross-Validation Algorithm With Column-Wise Centering and Scaling

Algorithm 4 can be extended to account for a subsequent training set-wise,
column-wise sample standard deviation-scaling of
$\mathbf{X}_{T_{p}\mathbf{c}}$ and $\mathbf{Y}_{T_{p}\mathbf{c}}$ directly in
$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{c}}\mathbf{X}_{T_{p}\mathbf{c}}$ and
$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{c}}\mathbf{Y}_{T_{p}}\mathbf{c}$. We
show how to compute the standard deviation-scaling only once for the entire
dataset and then for each partition, $p$, remove the contribution of any
sample $n\in V_{p}$.

Let $\widehat{\mathbf{x}_{R^{\prime}}}$ be the row vector of sample standard
deviations of $\mathbf{X}_{R^{\prime}}$, defined using Hadamard exponentiation
(∘) by

$\widehat{\mathbf{x}_{R^{\prime}}}=\left(\frac{1}{|R^{\prime}|-1}\sum_{n\in
R^{\prime}}\left(\left(\mathbf{X}_{n}-\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{R^{\prime}}}\right)^{\circ
2}\right)\right)^{\circ\frac{1}{2}}.$ (31)

Likewise, let Let $\widehat{\mathbf{y}}$ be the row vector of sample standard
deviations of $\mathbf{Y}_{R^{\prime}}$, defined by

$\widehat{\mathbf{y}_{R^{\prime}}}=\left(\frac{1}{|R^{\prime}|-1}\sum_{n\in
R^{\prime}}\left(\left(\mathbf{Y}_{n}-\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}\right)^{\circ
2}\right)\right)^{\circ\frac{1}{2}}.$ (32)

Note that when using $\widehat{\mathbf{x}_{R^{\prime}}}$ and
$\widehat{\mathbf{y}_{R^{\prime}}}$ for scaling, replacing any $0$-entries
with $1$-entries is standard practice to avoid division by zero. Also, note
that the sample standard deviation is only defined when $|R^{\prime}|>1$.When
scaling, the sample standard deviation must be computed on the training set to
avoid data leakage. So, the training set sample standard deviation is only
undefined when $|T_{p}|=1$, which is rarely the case.

Stacking these row vectors to obtain matrices with similar dimensions as
$\mathbf{X}_{R^{\prime}}$ and $\mathbf{Y}_{R^{\prime}}$, we get

$\widehat{\mathbf{X}_{R^{\prime}}}=\begin{bmatrix}\widehat{\mathbf{x}_{R^{\prime}}}\\\
\vdots\\\ \widehat{\mathbf{x}_{R^{\prime}}}\end{bmatrix}$ (33)

and

$\widehat{\mathbf{Y}_{R^{\prime}}}=\begin{bmatrix}\widehat{\mathbf{y}_{R^{\prime}}}\\\
\vdots\\\ \widehat{\mathbf{y}_{R^{\prime}}}\end{bmatrix}$ (34)

,

where the dimensions of $\widehat{\mathbf{X}_{R^{\prime}}}$ and
$\widehat{\mathbf{Y}_{R^{\prime}}}$ are $|R^{\prime}|\times K$ and
$|R^{\prime}|\times M$, respectively. We take
$\widehat{\mathbf{X}^{\mathbf{T}}_{R^{\prime}}}$ and
$\widehat{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}$ to mean
$\left(\widehat{\mathbf{X}_{R^{\prime}}}\right)^{\mathbf{T}}$ and
$\left(\widehat{\mathbf{x}_{R^{\prime}}}\right)^{\mathbf{T}}$, respectively,
and omit parentheses for brevity.

Then, the mean-centering and standard deviation-scaling of
$\mathbf{X}_{R^{\prime}}$ and $\mathbf{Y}_{R^{\prime}}$ can be defined using
Hadamard division ($\oslash$) by $\mathbf{X}_{R^{\prime}\mathbf{cs}}$ and
$\mathbf{Y}_{R^{\prime}\mathbf{cs}}$

as

$\mathbf{X}_{R^{\prime}\mathbf{cs}}=\mathbf{X}_{R^{\prime}\mathbf{c}}\oslash\widehat{\mathbf{X}_{R^{\prime}}}$
(35)

and

$\mathbf{Y}_{R^{\prime}\mathbf{cs}}=\mathbf{Y}_{R^{\prime}\mathbf{c}}\oslash\widehat{\mathbf{Y}_{R^{\prime}}},$
(36)

respectively.

Now, we present Algorithm 5, which we consider the immediate extension of the
naive Algorithm 3 to account for column-wise standard deviation-scaling.
Additionally, we present Algorithm 6, which extends Algorithm 4 to include
standard deviation scaling of $\mathbf{X}_{T_{p}\mathbf{c}}$ and
$\mathbf{Y}_{T_{p}\mathbf{c}}$ without requiring full recomputation of
training set matrix products or standard deviation vectors. Note that in some
steps, Algorithm 6 uses Hadamard multiplication ($\odot$).

Algorithm 5 Naive Cross-Validation Algorithm with Mean-Centering and Sample
Standard Deviation-Scaling

1:Determine $T_{p}$ and $V_{p}$ for each partition $p\in\\{1,\ldots,P\\}$ as
per Definition 3.1.

2:for each partition $p\in\\{1,\ldots,P\\}$ do

3: Compute $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}_{T_{p}}}$ and $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{T_{p}}}$.

4: Compute $\mathbf{X}_{T_{p}}-\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}_{T_{p}}}$ to obtain $\mathbf{X}_{T_{p}\mathbf{c}}$ and
$\mathbf{Y}_{T_{p}}-\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{T_{p}}}$ to obtain $\mathbf{Y}_{T_{p}\mathbf{c}}$.

5: Compute $\widehat{\mathbf{X}_{T_{p}}}$ and $\widehat{\mathbf{Y}_{T_{p}}}$.

6: Replace with $1$ any entry in $\widehat{\mathbf{X}_{T_{p}}}$ and
$\widehat{\mathbf{Y}_{T_{p}}}$ that is $0$.

7: Compute $\mathbf{X}_{T_{p}\mathbf{c}}\oslash\widehat{\mathbf{X}_{T_{p}}}$
to obtain $\mathbf{X}_{T_{p}\mathbf{cs}}$ and
$\mathbf{Y}_{T_{p}\mathbf{c}}\oslash\widehat{\mathbf{Y}_{T_{p}}}$ to obtain
$\mathbf{Y}_{T_{p}\mathbf{cs}}$.

8: Compute
$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{cs}}\mathbf{X}_{T_{p}\mathbf{cs}}$ and
$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{cs}}\mathbf{Y}_{T_{p}\mathbf{cs}}$.

9:end for

Algorithm 6 Cross-Validation Algorithm with Mean-Centering and Sample Standard
Deviation-Scaling

1:determine training indices $T_{p}$ and validation indices $V_{p}$ as per
Definition 3.1.

2:Compute $\mathbf{X}^{\mathbf{T}}\mathbf{X}$,
$\mathbf{X}^{\mathbf{T}}\mathbf{Y}$, $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}}$, $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}}$, $\sum_{n\in R}\mathbf{X}_{n}$, $\sum_{n\in
R}\mathbf{Y}_{n}$, $\sum_{n\in R}\left(\mathbf{X}_{n}^{\circ 2}\right)$, and
$\sum_{n\in R}\left(\mathbf{Y}_{n}^{\circ 2}\right)$.

3:for each partition $p\in\\{1,\ldots,P\\}$ do

4: Compute $\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{X}_{V_{p}}$ and
$\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{Y}_{V_{p}}$.

5: Compute
$\mathbf{X}^{\mathbf{T}}\mathbf{X}-\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{X}_{V_{p}}$
to obtain $\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{X}_{T_{p}}$, and
$\mathbf{X}^{\mathbf{T}}\mathbf{Y}-\mathbf{X}^{\mathbf{T}}_{V_{p}}\mathbf{Y}_{V_{p}}$
to obtain $\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{Y}_{T_{p}}$.

6: Compute $\frac{|V_{p}|}{|T_{p}|}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{V_{p}}}$, $\frac{|V_{p}|}{|T_{p}|}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{V_{p}}}$, $\frac{N}{|T_{p}|}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}}$, and $\frac{N}{|T_{p}|}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}}$.

7: Compute $\frac{N}{|T_{p}|}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}}-\frac{|V_{p}|}{|T_{p}|}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{V_{p}}}$ to obtain $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{T_{p}}}$ and $\frac{N}{|T_{p}|}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}}-\frac{|V_{p}|}{|T_{p}|}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{V_{p}}}$ to obtain $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{T_{p}}}$.

8: Compute $|T_{p}|\left(\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{T_{p}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{T_{p}}}\right)$ and $|T_{p}|\left(\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{T_{p}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{T_{p}}}\right)$.

9: Compute
$\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{X}_{T_{p}}-|T_{p}|\left(\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{T_{p}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{T_{p}}}\right)$ to obtain
$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{c}}\mathbf{X}_{T_{p}\mathbf{c}}$ and
$\mathbf{X}^{\mathbf{T}}_{T_{p}}\mathbf{Y}_{T_{p}}-|T_{p}|\left(\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}^{\mathbf{T}}_{T_{p}}}\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{T_{p}}}\right)$ to obtain
$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{c}}\mathbf{Y}_{T_{p}\mathbf{c}}$.

10: Compute $\sum_{n\in V_{p}}\mathbf{X}_{n}$ and $\sum_{n\in
V_{p}}\mathbf{Y}_{n}$.

11: Compute $\sum_{n\in V_{p}}\left(\mathbf{X}_{n}^{\circ 2}\right)$ and
$\sum_{n\in V_{p}}\left(\mathbf{Y}_{n}^{\circ 2}\right)$.

12: Compute $\sum_{n\in R}\mathbf{X}_{n}-\sum_{n\in V_{p}}\mathbf{X}_{n}$ to
obtain $\sum_{n\in T_{p}}\mathbf{X}_{n}$ and $\sum_{n\in
R}\mathbf{Y}_{n}-\sum_{n\in V_{p}}\mathbf{Y}_{n}$ to obtain $\sum_{n\in
T_{p}}\mathbf{Y}_{n}$.

13: Compute $\sum_{n\in R}\left(\mathbf{X}_{n}^{\circ 2}\right)-\sum_{n\in
V_{p}}\left(\mathbf{X}_{n}^{\circ 2}\right)$ to obtain $\sum_{n\in
T_{p}}\left(\mathbf{X}_{n}^{\circ 2}\right)$ and $\sum_{n\in
R}\left(\mathbf{Y}_{n}^{\circ 2}\right)-\sum_{n\in
V_{p}}\left(\mathbf{Y}_{n}^{\circ 2}\right)$ to obtain $\sum_{n\in
T_{p}}\left(\mathbf{Y}_{n}^{\circ 2}\right)$

14: Compute
$\widehat{\mathbf{x}_{T_{p}}}=\left(\frac{1}{|T_{p}|-1}\left(-2\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{T_{p}}}\odot\left(\sum_{n\in
T_{p}}\mathbf{X}_{n}\right)+|T_{p}|\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{T_{p}}}^{\circ 2}+\sum_{n\in T_{p}}\left(\mathbf{X}_{n}^{\circ
2}\right)\right)\right)^{\circ\frac{1}{2}}$.

15: Compute
$\widehat{\mathbf{y}_{T_{p}}}=\left(\frac{1}{|T_{p}|-1}\left(-2\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{T_{p}}}\odot\left(\sum_{n\in
T_{p}}\mathbf{Y}_{n}\right)+|T_{p}|\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{T_{p}}}^{\circ 2}+\sum_{n\in T_{p}}\left(\mathbf{Y}_{n}^{\circ
2}\right)\right)\right)^{\circ\frac{1}{2}}$.

16: Replace with $1$ any entry in $\widehat{\mathbf{x}_{T_{p}}}$ that is $0$.

17: Replace with $1$ any entry in $\widehat{\mathbf{y}_{T_{p}}}$ that is $0$.

18: Compute
$\widehat{\mathbf{x}^{\mathbf{T}}_{T_{p}}}\widehat{\mathbf{x}_{T_{p}}}$ and
$\widehat{\mathbf{x}^{\mathbf{T}}_{T_{p}}}\widehat{\mathbf{y}_{T_{p}}}$.

19: Compute
$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{c}}\mathbf{X}_{T_{p}\mathbf{c}}\oslash\left(\widehat{\mathbf{x}^{\mathbf{T}}_{T_{p}}}\widehat{\mathbf{x}_{T_{p}}}\right)$
to obtain
$\mathbf{X}^{\mathbf{T}}_{\mathbf{cs}T_{p}}\mathbf{X}_{\mathbf{cs}T_{p}}$ and
$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{c}}\mathbf{Y}_{T_{p}\mathbf{c}}\oslash\left(\widehat{\mathbf{x}^{\mathbf{T}}_{T_{p}}}\widehat{\mathbf{y}_{T_{p}}}\right)$
to obtain
$\mathbf{X}^{\mathbf{T}}_{\mathbf{cs}T_{p}}\mathbf{Y}_{\mathbf{cs}T_{p}}$.

20:end for

###  5.1 Correctness

Proposition 5.1 and Proposition 5.2 prove equalities that show us how to
efficiently compute mean-centered and standard deviation-scaled
$\mathbf{X}^{\mathbf{T}}_{R^{\prime}\mathbf{cs}}\mathbf{X}_{R^{\prime}\mathbf{cs}}$
and
$\mathbf{X}^{\mathbf{T}}_{R^{\prime}\mathbf{cs}}\mathbf{Y}_{R^{\prime}\mathbf{cs}}$
given
$\mathbf{X}^{\mathbf{T}}_{R^{\prime}\mathbf{c}}\mathbf{Y}_{R^{\prime}\mathbf{c}}$
and
$\mathbf{X}^{\mathbf{T}}_{R^{\prime}\mathbf{c}}\mathbf{Y}_{R^{\prime}\mathbf{c}}$.
Proposition 5.3 and Proposition 5.4 show how to efficiently derive the
training set standard deviations $\widehat{\mathbf{x}^{\mathbf{T}}_{T_{p}}}$
and $\widehat{\mathbf{y}_{T_{p}}}$. Proposition 5.5 and Proposition 5.6 show
additional equalities that are relevant to the derivations of
$\widehat{\mathbf{x}^{\mathbf{T}}_{T_{p}}}$ and
$\widehat{\mathbf{y}_{T_{p}}}$.

######  Proposition 5.1.

Let $\mathbf{X}_{R^{\prime}\mathbf{cs}}$ and
$\mathbf{Y}_{R^{\prime}\mathbf{cs}}$ be defined by Equation (35) and Equation
(36), respectively. Then
$\mathbf{X}^{\mathbf{T}}_{R^{\prime}\mathbf{cs}}\mathbf{Y}_{R^{\prime}\mathbf{cs}}=\left(\mathbf{X}^{\mathbf{T}}_{R^{\prime}\mathbf{c}}\mathbf{Y}_{R^{\prime}\mathbf{c}}\right)\oslash\left(\widehat{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}\widehat{\mathbf{y}_{R^{\prime}}}\right)$

###### Proof.

Consider an entry $(i,j)$ into
$\mathbf{X}^{\mathbf{T}}_{R^{\prime}\mathbf{cs}}\mathbf{Y}_{R^{\prime}\mathbf{cs}}$,
it is given by

$\left(\mathbf{X}^{\mathbf{T}}_{R^{\prime}\mathbf{cs}}\mathbf{Y}_{R^{\prime}\mathbf{cs}}\right)_{ij}=\left(\left(\mathbf{X}_{R^{\prime}\mathbf{c}}\oslash\widehat{\mathbf{X}_{R^{\prime}}}\right)^{\mathbf{T}}\left(\mathbf{Y}_{R^{\prime}\mathbf{c}}\oslash\widehat{\mathbf{Y}_{R^{\prime}}}\right)\right)_{ij}.$
(37)

Writing out the sums yields

$\sum_{n\in
R^{\prime}}\frac{\mathbf{X}^{\mathbf{T}}_{R^{\prime}\mathbf{c}_{{}_{in}}}}{\widehat{\mathbf{X}^{\mathbf{T}}_{R^{\prime}}}_{{}_{in}}}\frac{\mathbf{Y}_{R^{\prime}\mathbf{c}_{{}_{nj}}}}{\widehat{\mathbf{Y}_{R^{\prime}}}_{{}_{nj}}}=\sum_{n\in
R^{\prime}}\frac{\mathbf{X}^{\mathbf{T}}_{R^{\prime}\mathbf{c}_{{}_{in}}}\mathbf{Y}_{R^{\prime}\mathbf{c}_{{}_{nj}}}}{\widehat{\mathbf{X}^{\mathbf{T}}_{R^{\prime}}}_{{}_{in}}\widehat{\mathbf{Y}_{R^{\prime}}}_{{}_{nj}}}$
(38)

By Equation (10) and Equation (11), we have that every column, $n$, in
$\widehat{\mathbf{X}^{\mathbf{T}}_{R^{\prime}}}$ is
$\widehat{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}$ and every row, $n$, in
$\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{R^{\prime}}}$ is $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}$. Using these definitions and the distributivity
of division over summation, we get

$\sum_{n\in
R^{\prime}}\frac{\mathbf{X}^{\mathbf{T}}_{R^{\prime}\mathbf{c}_{{}_{in}}}\mathbf{Y}_{R^{\prime}\mathbf{c}_{{}_{nj}}}}{\widehat{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}_{{}_{i}}\widehat{\mathbf{y}_{R^{\prime}}}_{{}_{j}}}=\frac{\sum_{n\in
R^{\prime}}\mathbf{X}^{\mathbf{T}}_{R^{\prime}\mathbf{c}_{{}_{in}}}\mathbf{Y}_{R^{\prime}\mathbf{c}_{{}_{nj}}}}{\widehat{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}_{{}_{i}}\widehat{\mathbf{y}_{R^{\prime}}}_{{}_{j}}}.$
(39)

By the definitions of matrix products, outer vector product, and Hadamard
division, we get

$\frac{\left(\mathbf{X}^{\mathbf{T}}_{R^{\prime}\mathbf{c}}\mathbf{Y}_{R^{\prime}\mathbf{c}}\right)_{ij}}{\left(\widehat{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}\widehat{\mathbf{y}_{R^{\prime}}}\right)_{ij}}=\left(\left(\mathbf{X}^{\mathbf{T}}_{R^{\prime}\mathbf{c}}\mathbf{Y}_{R^{\prime}\mathbf{c}}\right)\oslash\left(\widehat{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}\widehat{\mathbf{y}_{R^{\prime}}}\right)\right)_{ij}.$
(40)

∎

######  Proposition 5.2.

Let $\mathbf{X}_{R^{\prime}\mathbf{cs}}$ be defined by Equation (35). Then
$\mathbf{X}^{\mathbf{T}}_{R^{\prime}\mathbf{cs}}\mathbf{X}_{R^{\prime}\mathbf{cs}}=\left(\mathbf{X}^{\mathbf{T}}_{R^{\prime}\mathbf{c}}\mathbf{X}_{R^{\prime}\mathbf{c}}\right)\oslash\left(\widehat{\mathbf{x}^{\mathbf{T}}_{R^{\prime}}}\widehat{\mathbf{x}_{R^{\prime}}}\right)$

###### Proof.

As of proof of Proposition 5.2 with
$\mathbf{Y}_{R^{\prime}}=\mathbf{X}_{R^{\prime}}$. ∎

######  Proposition 5.3.

Given the definition of $\widehat{\mathbf{x}_{R^{\prime}}}$ in Equation (31),
an equivalent definition is

$\widehat{\mathbf{x}_{R^{\prime}}}=\left(\frac{1}{|R^{\prime}|-1}\left(-2\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{R^{\prime}}}\odot\left(\sum_{n\in
R^{\prime}}\mathbf{X}_{n}\right)+|R^{\prime}|\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{R^{\prime}}}^{\circ 2}+\sum_{n\in
R^{\prime}}\left(\mathbf{X}_{n}^{\circ
2}\right)\right)\right)^{\circ\frac{1}{2}}$ (41)

###### Proof.

Consider expanding the Hadamard square in the definition of
$\widehat{\mathbf{x}_{R^{\prime}}}$ from Equation (31) to get

$\left(\frac{1}{|R^{\prime}|-1}\sum_{n\in
R^{\prime}}\left(\left(\mathbf{X}_{n}-\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{R^{\prime}}}\right)^{\circ
2}\right)\right)^{\circ\frac{1}{2}}=\left(\frac{1}{|R^{\prime}|-1}\sum_{n\in
R^{\prime}}\left(\mathbf{X}_{n}^{\circ 2}+\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{R^{\prime}}}^{\circ 2}-2\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{R^{\prime}}}\odot\mathbf{X}_{n}\right)\right)^{\circ\frac{1}{2}}.$
(42)

By commutativity of addition, the distributivity of Hadamard multiplication
over addition, computing the middle sum, and rearranging terms, we get

$\left(\frac{1}{|R^{\prime}|-1}\left(\sum_{n\in
R^{\prime}}\left(\mathbf{X}_{n}^{\circ 2}\right)+\sum_{n\in
R^{\prime}}\left(\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{R^{\prime}}}^{\circ 2}\right)-\sum_{n\in
R^{\prime}}2\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{R^{\prime}}}\odot\mathbf{X}_{n}\right)\right)^{\circ\frac{1}{2}}=\left(\frac{1}{|R^{\prime}|-1}\left(-2\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{R^{\prime}}}\odot\left(\sum_{n\in
R^{\prime}}\mathbf{X}_{n}\right)+|R^{\prime}|\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{R^{\prime}}}^{\circ 2}+\sum_{n\in
R^{\prime}}\left(\mathbf{X}_{n}^{\circ
2}\right)\right)\right)^{\circ\frac{1}{2}}$ (43)

∎

######  Proposition 5.4.

Given the definition of $\widehat{\mathbf{y}_{R^{\prime}}}$ in Equation (32),
an equivalent definition is

$\widehat{\mathbf{y}_{R^{\prime}}}=\left(\frac{1}{|R^{\prime}|-1}\left(-2\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}\odot\left(\sum_{n\in
R^{\prime}}\mathbf{Y}_{n}\right)+|R^{\prime}|\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{R^{\prime}}}^{\circ 2}+\sum_{n\in
R^{\prime}}\left(\mathbf{Y}_{n}^{\circ
2}\right)\right)\right)^{\circ\frac{1}{2}}$ (44)

###### Proof.

As proof of Proposition 5.3 with
$\widehat{\mathbf{x}_{R^{\prime}}}=\widehat{\mathbf{y}_{R^{\prime}}}$. ∎

######  Proposition 5.5.

Given Definition 3.1, we have

$\sum_{n\in T_{p}}\mathbf{X}_{n}=\sum_{n\in R}\mathbf{X}_{n}-\sum_{n\in
V_{p}}\mathbf{X}_{n}.$ (45)

###### Proof.

We show the equivalent statement $\sum_{n\in T_{p}}\mathbf{X}_{n}+\sum_{n\in
V_{p}}=\sum_{n\in R}\mathbf{X}_{n}$. By Definition 3.1, we have that
$T_{p}\cap V_{p}=\emptyset$ and $T_{p}\cup V_{p}=R$. Therefore, we can write

$\sum_{n\in T_{p}}\mathbf{X}_{n}+\sum_{n\in V_{p}}\mathbf{X}_{n}=\sum_{n\in
T_{p}\cup V_{p}}\mathbf{X}_{n}=\sum_{n\in R}\mathbf{X}_{n}.$ (46)

∎

######  Proposition 5.6.

Given Definition 3.1, we have

$\sum_{n\in T_{p}}\mathbf{Y}_{n}=\sum_{n\in R}\mathbf{Y}_{n}-\sum_{n\in
V_{p}}\mathbf{Y}_{n},$ (47)

$\sum_{n\in T_{p}}\left(\mathbf{X}_{n}^{\circ 2}\right)=\sum_{n\in
R}\left(\mathbf{X}_{n}^{\circ 2}\right)-\sum_{n\in
V_{p}}\left(\mathbf{X}_{n}^{\circ 2}\right),$ (48)

and

$\sum_{n\in T_{p}}\left(\mathbf{Y}_{n}^{\circ 2}\right)=\sum_{n\in
R}\left(\mathbf{Y}_{n}^{\circ 2}\right)-\sum_{n\in
V_{p}}\left(\mathbf{Y}_{n}^{\circ 2}\right),$ (49)

.

###### Proof.

As of proof of Proposition 5.5 but with $\mathbf{X}=\mathbf{Y}$,
$\mathbf{X}=\mathbf{X}^{\circ 2}$, and $\mathbf{X}=\mathbf{Y}^{\circ 2}$,
respectively.

∎

######  Proposition 5.7 (Correctness of Algorithm 6).

$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{cs}}\mathbf{X}_{T_{p}\mathbf{cs}}$ and
$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{cs}}\mathbf{Y}_{T_{p}\mathbf{cs}}$
computed in step 19 of Algorithm 6 are identical to those computed in step 8
of Algorithm 5.

###### Proof.

Except for some additional computations in step 1, steps 1-9 of Algorithm 6
are identical to steps 1-9 of Algorithm 4 and so correctness of these steps
follows from Proposition 4.7. Correctness of $\sum_{n\in
T_{p}}\mathbf{X}_{n}$, $\sum_{n\in T_{p}}\mathbf{Y}_{n}$, $\sum_{n\in
T_{p}}\left(\mathbf{X}_{n}^{\circ 2}\right)$, and $\sum_{n\in
T_{p}}\left(\mathbf{Y}_{n}^{\circ 2}\right)$ computed in steps 12-13 follows
from Proposition 5.5 and Proposition 5.6. Correctness of
$\widehat{\mathbf{x}_{T_{p}}}$ and $\widehat{\mathbf{y}_{T_{p}}}$ computed in
steps 14-17 follows from Proposition 5.3, Proposition 5.4 and the standard
practice of replacing with $1$ any $0$-entries in these vectors, which is also
done in step 6 of Algorithm 5. Finally, correctness of
$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{cs}}\mathbf{X}_{T_{p}\mathbf{cs}}$ and
$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{cs}}\mathbf{Y}_{T_{p}\mathbf{cs}}$
computed in step 19 follows from Proposition 5.2 and Proposition 5.1,
respectively. ∎

###  5.2 Computational Complexity

Here, we analyze the runtime and space complexities of Algorithm 5 and
Algorithm 6. Again, we disregard any additions part of a matrix multiplication
operation.

######  Proposition 5.8.

Algorithm 5 requires a total of $\Theta(PN(K+M))$ additions.

###### Proof.

It follows from Proposition 4.8 that steps 2-4 require $\Theta(PN(K+M))$
additions.

For some $p$, step 5 requires $|T_{p}|(K+M)$ additions to compute
$\widehat{\mathbf{X}_{T_{p}}}$ and $\widehat{\mathbf{Y}_{T_{p}}}$ as
$\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{T_{p}}}$ and $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{T_{p}}}$ have already been computed in step 3. Iterating over
all partitions in step 2, we use Equation (5) to derive a total of
$\Theta(PN(K+M))$ additions in step 5.

Thus, the total number of additions is

$2\Theta(PN(K+M))=\Theta(PN(K+M)).$ (50)

∎

######  Proposition 5.9.

Algorithm 5 requires a total of $\Theta(PN(K+M))$ subtractions.

###### Proof.

It follows from Proposition 4.9 that steps 2-4 require $\Theta(PN(K+M))$
subtractions.

For some $p$, step 5 requires $|T_{p}|(K+M)+2$ subtractions to compute
$\widehat{\mathbf{X}_{T_{p}}}$ and $\widehat{\mathbf{Y}_{T_{p}}}$. Iterating
over all partitions in step 2, we use Equation (5) to derive a total of
$\Theta(PN(K+M))+2P$ subtractions in step 5.

Thus, the total number of subtractions is

$2\Theta(PN(K+M))+2P=\Theta(PN(K+M)).$ (51)

∎

######  Proposition 5.10.

Algorithm 5 requires a total of $\Theta(PNK(K+M))$ multiplications.

###### Proof.

Algorithm 5 computes the $\Theta(PNK(K+M))$ also required by Algorithm 3 given
by Proposition 4.10. The multiplications covered are those for all iterations
of steps 3 and 8 in Algorithm 5.

Additionally, for some $p$, step 5 of Algorithm 5 requires $|T_{p}|(K+M)+K+M$
multiplications to compute $\widehat{\mathbf{X}_{T_{p}}}$ and
$\widehat{\mathbf{Y}_{T_{p}}}$ as $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{T_{p}}}$ and $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{T_{p}}}$ have already been computed in step 3. These
multiplications arise from the need to compute the Hadamard squares in
Equation (31) and Equation (32). Iterating over all partitions in step 2, we
use Equation (5) to derive a total of $\Theta(PN(K+M))$ multiplications in
step 5.

Thus, the total number of multiplications is

$\Theta(PNK(K+M))+\Theta(PN(K+M))=\Theta(PNK(K+M)).$ (52)

∎

######  Proposition 5.11.

Algorithm 5 requires a total of $\Theta(PK(K+M))$ divisions.

###### Proof.

It follows from Proposition 4.11 that iterating over all partitions, step 3
requires $\Theta(P)$ divisions.

Iterating over all partitions, step 5 requires $2P=\Theta(P)$ divisions.

Step 7 requires $|T_{p}|(K+M)$ divisions. Iterating over all partitions in
step 2, we use Equation (5) to derive a total of $\Theta(PN(K+M))$ divisions
in step 7.

Thus, the total number of divisions is

$2\Theta(P)+PN(K+M)=\Theta(PN(K+M)).$ (53)

∎

######  Proposition 5.12.

Algorithm 5 requires a total of $\Theta(P(K+M))$ square roots.

###### Proof.

The square roots arise from the need to compute the Hadamard exponents of
$\frac{1}{2}$ in Equation (31) and Equation (32). These are required for step
5. For any partition, $p$, there are a total of $K+M$ square roots to compute.
Iterating over all partitions, we get $\Theta(P(K+M))$ square roots. ∎

######  Proposition 5.13.

Algorithm 5 requires a total of $O(P(K+M))$ replacements.

###### Proof.

Replacements arise from the need to replace zero-entries in
$\widehat{\mathbf{X}_{T_{p}}}$ and $\widehat{\mathbf{Y}_{T_{p}}}$ with ones.
$\widehat{\mathbf{X}_{T_{p}}}$ and $\widehat{\mathbf{Y}_{T_{p}}}$ are uniquely
defined by $\widehat{\mathbf{x}_{T_{p}}}$ and $\widehat{\mathbf{y}_{T_{p}}}$,
respectively. Thus, any efficient implementation need only replacements in
$\widehat{\mathbf{x}_{T_{p}}}$ and $\widehat{\mathbf{y}_{T_{p}}}$ to have them
take effect in $\widehat{\mathbf{X}_{T_{p}}}$ and
$\widehat{\mathbf{Y}_{T_{p}}}$. At most, $K+M$ entries are to be replaced in
each iteration. Iterating over all partitions gives an upper bound of
$O(P(K+M))$ replacements. ∎

######  Proposition 5.14.

Algorithm 6 requires a total of $\Theta(N(K+M))$ additions.

###### Proof.

Except for some additional computations in step 2, steps 1-9 of Algorithm 6
are identical to those of Algorithm 4 where a total of $\Theta(N(K+M))$
additions follows from Proposition 4.12. The additional addition operations in
step 2 of Algorithm 6 are due to the computations of $\sum_{n\in
R}\mathbf{X}_{n}$, $\sum_{n\in R}\mathbf{Y}_{n}$, $\sum_{n\in
R}\left(\mathbf{X}_{n}^{\circ 2}\right)$, and $\sum_{n\in
R}\left(\mathbf{Y}_{n}^{\circ 2}\right)$. These require $2N(K+M)$ additions.

For any partition, $p$, steps 10-11 require $2|V_{p}|(K+M)$ additions.
Iterating over all partitions in step 3, we use Equation (4) to derive a total
of $\Theta(N(K+M))$ additions in steps 10-11.

The additions in steps 12-13 have already been computed in previous steps.

The additions required to compute the terms in steps 14-15 have already been
computed in previous steps. Thus, we only need to sum the terms, requiring
$2K+2M$ additions per iteration. Iterating over all partitions, this gives a
total of $\Theta(P(K+M))=O(N(K+M))$ additions for steps 14-15.

Thus, the total number of additions is

$2\Theta(N(K+M))+2N(K+M)+O(N(K+M))=\Theta(N(K+M)).$ (54)

∎

######  Proposition 5.15.

Algorithm 6 requires a total of $\Theta(N(K+M))$ subtractions.

###### Proof.

Except for some additional computations in step 2, steps 1-9 of Algorithm 6
are identical to those of Algorithm 4 where a total of $\Theta(N(K+M))$
subtractions follows from Proposition 4.13. The additional operations in step
2 of Algorithm 6 require no subtraction operations.

For any partition, $p$, steps 12-13 require $K+M$ subtractions. Iterating over
all partitions in step 3 gives us a total of $P(K+M)=O(N(K+M))$ subtractions
in steps 12-13.

Steps 14-15 each have only one subtraction operation. Iterating over all
partitions in step 3 gives us a total of $2P=O(N)$ subtractions in steps
14-15.

Thus, the total number of subtractions is

$\Theta(N(K+M))+2O(N(K+M))+O(N)=\Theta(N(K+M))$ (55)

∎

######  Proposition 5.16.

Algorithm 6 requires a total of $\Theta(NK(K+M))$ multiplications.

###### Proof.

Except for some additional computations in step 2, steps 1-9 of Algorithm 6
are identical to those of Algorithm 4 where a total of $\Theta(NK(K+M))$
multiplications follows from Proposition 4.14. The additional operations that
require multiplications in step 2 of Algorithm 6 are $\sum_{n\in
R}\left(\mathbf{X}_{n}^{\circ 2}\right)$, and $\sum_{n\in
R}\left(\mathbf{Y}_{n}^{\circ 2}\right)$. These require $N(K+M)$
multiplications.

For any partition, $p$, step 11 requires $|V_{p}|(K+M)$ multiplications.
Iterating over all partitions in step 3, we use Equation (4) to derive a total
of $N(K+M)$ multiplications for step 11.

All the multiplications in step 13 have been computed in previous steps.

Most expressions in steps 14-15 have been computed in previous steps. The only
remaining multiplication operations are the scaling with
$\frac{1}{|T_{p}|-1}$, the scaling with $-2$ in the first summand, the
Hadamard product in the first summand, the scaling by $|T_{p}|$ in the second
summand, and the Hadamard square in the second summand. Each of these requires
$K$ multiplications in step 14 and $M$ multiplications in step 15 for a total
of $5(K+M)$ multiplications. Iterating over all partitions, we get a total of
$5P(K+M)=O(N(K+M))$ multiplications in steps 14-15.

For any partition, $p$, step 18 requires $K(K+M)$ multiplications. Iterating
over all partitions in step 3, we get a total of $PK(K+M)=O(NK(K+M))$
multiplications in step 18.

$\Theta(NK(K+M))+2N(K+M)+O(N(K+M))+O(NK(K+M))=\Theta(NK(K+M))$ (56)

∎

######  Proposition 5.17.

Algorithm 6 requires a total of $O(NK(K+M))$ divisions.

###### Proof.

Except for some additional computations in step 2, steps 1-9 of Algorithm 6
are identical to those of Algorithm 4 where a total of $O(N)$ multiplications
follows from Proposition 4.15. The additional operations in step 2 of
Algorithm 6 require no division operations.

For any partition, $p$, steps 14-15 require $2$ division operations. Iterating
over all partitions, we get a total of $2P=O(N)$ divisions in steps 14-15.

For any partition, $p$, step 19 requires $K(K+M)$ division operations.
Iterating over all partitions, we get a total of $PK(K+M)=O(NK(K+M))$
divisions in step 19.

$2O(N)+O(NK(K+M))=O(NK(K+M))$ (57)

∎

######  Proposition 5.18.

Algorithm 6 requires a total of $\Theta(P(K+M)=O(N(K+M))$ square roots.

###### Proof.

The square roots arise from the need to compute the Hadamard exponents of
$\frac{1}{2}$ in Equation (31) and Equation (32). These are required for steps
14-15. For any partition, $p$, there are a total of $K+M$ square roots to
compute. Iterating over all partitions, we get a total of
$\Theta(P(K+M))=O(N(K+M))$ square roots. ∎

######  Proposition 5.19.

Algorithm 6 requires a total of $O(P(K+M))=O(N(K+M))$ replacements.

###### Proof.

Replacements arise from the need to replace zero-entries in
$\widehat{\mathbf{x}_{T_{p}}}$ and $\widehat{\mathbf{y}_{T_{p}}}$ with ones.
At most, $K+M$ entries are to be replaced in each iteration. Iterating over
all partitions, this gives an upper bound of $O(P(K+M))=O(N(K+M))$
replacements. ∎

It follows that Algorithm 5 requires

$3\Theta(PN(K+M))+\Theta(PNK(K+M))+\Theta(P(K+M))+O(P(K+M))=\Theta(PNK(K+M))$
(58)

total operations. It also follows that Algorithm 6 requires

$2\Theta(N(K+M))+\Theta(NK(K+M))+O(NK(K+M))+2O(N(K+M))=\Theta(NK(K+M))$ (59)

total operations. Therefore Algorithm 6 is faster than Algorithm 5.

######  Proposition 5.20.

Algorithm 5 requires storing $\Theta((K+N)(K+M))$ matrix entries.

###### Proof.

Storing $\mathbf{X}$ and $\mathbf{Y}$ requires $N(K+M)$ entries. Storing all
$V_{p}$ and $T_{p}$ requires $N$ entries. In step 3, storing $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{X}_{T_{p}}}$ and $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{Y}_{T_{p}}}$ can be accomplished by only storing $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{x}_{T_{p}}}$ and $\macc@depth\char
1\relax\frozen@everymath{\macc@group}\macc@set@skewchar\macc@nested@a
111{\mathbf{y}_{T_{p}}}$, requiring $K+M$ entries. In step 4,
$\mathbf{X}_{T_{p}}$ and $\mathbf{X}_{T_{p}\mathbf{c}}$ can share the same
memory, as can $\mathbf{Y}_{T_{p}}$ and $\mathbf{Y}_{T_{p}\mathbf{c}}$. Thus,
the storage requirement for step 4 is $N(K+M)$. In step 5, storing
$\widehat{\mathbf{X}_{T_{p}}}$ and $\widehat{\mathbf{Y}_{T_{p}}}$ can be
accomplished by only storing $\widehat{\mathbf{x}_{T_{p}}}$ and
$\widehat{\mathbf{y}_{T_{p}}}$, requiring $K+M$ entries. In step 7,
$\mathbf{X}_{T_{p}\mathbf{cs}}$ and $\mathbf{Y}_{T_{p}\mathbf{cs}}$ can
overwrite the memory locations of $\mathbf{X}_{T_{p}\mathbf{c}}$ and
$\mathbf{Y}_{T_{p}\mathbf{c}}$, respectively, for no additional storage
requirement. In step 8, storing
$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{cs}}\mathbf{X}_{T_{p}\mathbf{cs}}$ and
$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{cs}}\mathbf{Y}_{T_{p}\mathbf{cs}}$
requires $K(K+M)$ entries. Thus, the total storage requirement is

$2N(K+M)+N+2K+2M+K(K+M)=\Theta((K+N)(K+M)).$ (60)

∎

######  Proposition 5.21.

Algorithm 6 requires storing $\Theta((K+N)(K+M))$ matrix entries.

###### Proof.

Except for some computation of some additional terms in step 2, steps 1-9 of
Algorithm 6 are identical to those of Algorithm 4 where a total storage
requirement of $\Theta((K+N)(K+M))$ matrix entries follows from Proposition
4.17. The additional terms that must be stored in step 2 are $\sum_{n\in
R}\mathbf{X}_{n}$, $\sum_{n\in R}\mathbf{Y}_{n}$, $\sum_{n\in
R}\left(\mathbf{X}_{n}^{\circ 2}\right)$, and $\sum_{n\in
R}\left(\mathbf{Y}_{n}^{\circ 2}\right)$. They require storage of $2(K+M)$
matrix entries. Steps 10-11 require storage of $2(K+M)$ matrix entries for
$\sum_{n\in V_{p}}\mathbf{X}_{n}$, $\sum_{n\in V_{p}}\mathbf{Y}_{n}$,
$\sum_{n\in V_{p}}\left(\mathbf{X}_{n}^{\circ 2}\right)$, and $\sum_{n\in
V_{p}}\left(\mathbf{Y}_{n}^{\circ 2}\right)$. In steps 12-13, $\sum_{n\in
T_{p}}\mathbf{X}_{n}$, $\sum_{n\in T_{p}}\mathbf{Y}_{n}$, $\sum_{n\in
T_{p}}\left(\mathbf{X}_{n}^{\circ 2}\right)$, and $\sum_{n\in
T_{p}}\left(\mathbf{Y}_{n}^{\circ 2}\right)$ can overwrite the memory
locations of $\sum_{n\in V_{p}}\mathbf{X}_{n}$, $\sum_{n\in
V_{p}}\mathbf{Y}_{n}$, $\sum_{n\in V_{p}}\left(\mathbf{X}_{n}^{\circ
2}\right)$, and $\sum_{n\in V_{p}}\left(\mathbf{Y}_{n}^{\circ 2}\right)$,
respectively, for no additional storage requirement. Steps 14-17 require
storing $K+M$ matrix entries for $\widehat{\mathbf{x}_{T_{p}}}$ and
$\widehat{\mathbf{y}_{T_{p}}}$. Step 18 requires storing $K(K+M)$ matrix
entries for
$\widehat{\mathbf{x}^{\mathbf{T}}_{T_{p}}}\widehat{\mathbf{x}_{T_{p}}}$ and
$\widehat{\mathbf{x}^{\mathbf{T}}_{T_{p}}}\widehat{\mathbf{y}_{T_{p}}}$. Step
19 requires no additional storage as
$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{cs}}\mathbf{X}_{T_{p}\mathbf{cs}}$ and
$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{cs}}\mathbf{Y}_{T_{p}\mathbf{cs}}$ can
overwrite the memory locations of either
$\widehat{\mathbf{x}^{\mathbf{T}}_{T_{p}}}\widehat{\mathbf{x}_{T_{p}}}$ and
$\widehat{\mathbf{x}^{\mathbf{T}}_{T_{p}}}\widehat{\mathbf{y}_{T_{p}}}$,
respectively, or
$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{c}}\mathbf{X}_{T_{p}\mathbf{c}}$ and
$\mathbf{X}^{\mathbf{T}}_{T_{p}\mathbf{c}}\mathbf{Y}_{T_{p}\mathbf{c}}$,
respectively. Thus, the total storage requirement is

$\Theta((K+N)(K+M))+5(K+M)+K(K+M)=\Theta((K+N)(K+M)).$ (61)

∎

Thus, Algorithm 5 and Algorithm 6 have identical asymptotic space
complexities.

##  6 Parallel algorithms

Algorithm 2, Algorithm 4, and Algorithm 6 can all have their loops
parallelized for speeding up computation in a multi-processor setting. Let $C$
be the number of parallel processes used; then the storage requirement
increases from $\Theta((K+N)(K+M))$ to $\Theta((CK+N)(K+M))$ while reducing
the runtime of the loop from $\Theta(NK(K+M))$ to $\frac{\Theta(NK(K+M))}{C}$.
Regarding the number of operations required, the total runtime, including step
$1$, will still be $\Theta(NK(K+M))$. In practice, multiprocessing can greatly
speed up cross-validation, although the asymptotic runtime is unchanged from
the sequential alternative. With the possibility of parallelizing the loops,
all three algorithms are implemented in [8]. The implementations are available
for installation on the Python Package Index (PyPI) under the name ikpls.

##  7 Conclusion

In this work, we derive three algorithms that can efficiently, in both space
and time, compute training set-wise $\mathbf{X}^{\mathbf{T}}\mathbf{X}$ and
$\mathbf{X}^{\mathbf{T}}\mathbf{Y}$ as required in a cross-validation setting.
Applying any preprocessing that only depends on the sample on which it is
applied can be performed before initializing our algorithms and does not
change the correctness of our algorithms. However, when centering and scaling
the columns of $\mathbf{X}$ and $\mathbf{Y}$, the centering and scaling depend
on dataset statistics. These dataset statistics must only be computed on the
training set during cross-validation. Our algorithms provide elegant solutions
for computing training set-wise $\mathbf{X}^{\mathbf{T}}\mathbf{X}$ and
$\mathbf{X}^{\mathbf{T}}\mathbf{Y}$ where $\mathbf{X}$ and $\mathbf{Y}$ have
possibly been centered around the training set-wise means and scaled around
the training set-wise sample standard deviations - without fully recomputing
any training set matrices or statistical moments. The algorithms work on the
principle of computing $\mathbf{X}^{\mathbf{T}}\mathbf{X}$,
$\mathbf{X}^{\mathbf{T}}\mathbf{Y}$, means, and standard deviations only once
for the entire dataset in the beginning, and then, during each cross-
validation partition, remove the contribution of the validation set. We prove
this is equivalent to but much faster than directly computing each partition’s
training set-wise matrices and statistical moments. Additionally, we prove
that our algorithms require the same amount of storage as the standard cross-
validation algorithms. Finally, we note that our algorithms are well-suited
for parallelization and provide a reference to a working Python implementation
for all algorithms using the extremely fast Improved Kernel PLS Algorithms.

While we have shown how to efficiently compute possibly mean-centered and
possibly standard deviation-scaled partition-wise
$\mathbf{X}^{\mathbf{T}}\mathbf{X}$ and $\mathbf{X}^{\mathbf{T}}\mathbf{Y}$,
subsequent training of the PLS model to derive the PLS regression matrix is
still required for each cross-validation partition. We wonder if our approach
can be extended to derive the partition-wise PLS regression coefficient
matrices directly to provide an algorithm for cross-validation PLS
calibrations completely independent of the number of cross-validation
partitions.

## References

  * [1] A. Alin.  Comparison of pls algorithms when number of objects is much larger than number of variables.  Statistical papers, 50:711–720, 2009.
  * [2] M. Andersson.  A comparison of nine pls1 algorithms.  Journal of Chemometrics: A Journal of the Chemometrics Society, 23(10):518–529, 2009.
  * [3] M. Barker and W. Rayens.  Partial least squares for discrimination.  Journal of Chemometrics: A Journal of the Chemometrics Society, 17(3):166–173, 2003.
  * [4] R. J. Barnes, M. S. Dhanoa, and S. J. Lister.  Standard normal variate transformation and de-trending of near-infrared diffuse reflectance spectra.  Applied spectroscopy, 43(5):772–777, 1989.
  * [5] Daniel Berrar et al.  Cross-validation., 2019.
  * [6] B. S. Dayal and J. F. MacGregor.  Improved pls algorithms.  Journal of Chemometrics: A Journal of the Chemometrics Society, 11(1):73–85, 1997.
  * [7] Z. Du, W. Tian, M. Tilley, D. Wang, G. Zhang, and Y. Li.  Quantitative assessment of wheat quality using near-infrared spectroscopy: A comprehensive review.  Comprehensive Reviews in Food Science and Food Safety, 21(3):2956–3009, 2022.
  * [8] O.-C. G. Engstrøm, E.S. Dreier, B.M. Jespersen, and K.S. Pedersen.  Fast cross-validation algorithm with improved kernel partial least squares for python: Fast cpu and gpu implementations with numpy and jax (under review).  The Journal of Open Source Software, 2024.
  * [9] T. Hastie, R. Tibshirani, and J. H. Friedman.  The elements of statistical learning: data mining, inference, and prediction, volume 2.  Springer, 2009.
  * [10] F. Le Gall.  Powers of tensors and fast matrix multiplication.  In Proceedings of the 39th international symposium on symbolic and algebraic computation, pages 296–303, 2014.
  * [11] K. H. Liland, P. Stefansson, and U. G. Indahl.  Much faster cross-validation in plsr-modelling by avoiding redundant calculations.  Journal of Chemometrics, 34(3):e3201, 2020.
  * [12] F. Lindgren, P. Geladi, and S. Wold.  The kernel algorithm for pls.  Journal of Chemometrics, 7(1):45–59, 1993.
  * [13] F. Lindgren, P. Geladi, and S. Wold.  Kernel-based pls regression; cross-validation and applications to spectral data.  Journal of chemometrics, 8(6):377–389, 1994.
  * [14] Å. Rinnan, F. van den Berg, and S. B. Engelsen.  Review of the most common pre-processing techniques for near-infrared spectra.  TrAC Trends in Analytical Chemistry, 28(10):1201–1222, 2009.
  * [15] A. Savitzky and M. J. E. Golay.  Smoothing and differentiation of data by simplified least squares procedures.  Analytical chemistry, 36(8):1627–1639, 1964.
  * [16] K. M. Sørensen, F. van den Berg, and S. B. Engelsen.  Nir data exploration and regression by chemometrics—a primer.  Near-Infrared Spectroscopy: Theory, Spectral Analysis, Instrumentation, and Applications, pages 127–189, 2021.
  * [17] P. Stefansson, U. G. Indahl, K. H. Liland, and I. Burud.  Orders of magnitude speed increase in partial least squares feature selection with new simple indexing technique for very tall data sets.  Journal of Chemometrics, 33(11):e3141, 2019.
  * [18] H. Wold.  Estimation of principal components and related models by iterative least squares.  Multivariate analysis, pages 391–420, 1966.
  * [19] S. Wold, M. Sjöström, and L. Eriksson.  Pls-regression: a basic tool of chemometrics.  Chemometrics and intelligent laboratory systems, 58(2):109–130, 2001.