xcp ref implementation #2895
Open
xcp ref implementation #2895
Conversation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This routine computes the matrix of cross product of data stored in column major format, in batches. For matrix X of dimensions p x n, the i,j th entry of the cross product matrix is C_ij = \sum_k (x_ik-\mu_i) (x_jk-\mu_k) where x_ij is the jth element of the ith row, of the matrix X. Implementation uses the BLAS routine GEMM. Signed-off-by: Dhanus M Lal <[email protected]>
DhanusML
requested review from
Alexsandruss,
samir-nasibli and
Alexandr-Solovev
as code owners
Signed-off-by: Dhanus M Lal <[email protected]>
| double alpha; | ||
| if (accWtOld != 0) | ||
| { | ||
| double * sumOld = daal::services::internal::service_malloc<double, cpu>(nFeatures, sizeof(double)); |
There was a problem hiding this comment.
Suggested change
| double * sumOld = daal::services::internal::service_malloc<double, cpu>(nFeatures, sizeof(double)); | |
| double* const sumOld = daal::services::internal::service_malloc<double, cpu>(nFeatures, sizeof(double)); |
| sumOld[i] = sum[i]; | ||
| } | ||
| // S_old S_old^t/accWtOld | ||
| alpha = 1.0 / accWtOld; |
There was a problem hiding this comment.
Any checks for overflow?
Comment on lines
+223
to
+225
| transb = 'T'; | ||
| alpha = 1.0; | ||
| beta = 1.0; |
There was a problem hiding this comment.
I would recommend not to use the same variables. It's confusing and also can be dangerous if someone will forgot to change them.
Comment on lines
+223
to
+227
| transb = 'T'; | ||
| alpha = 1.0; | ||
| beta = 1.0; | ||
| blasInst.xgemm(&transa, &transb, &nFeatures, &nFeatures, &nVectors, &alpha, data, &nFeatures, data, &nFeatures, &beta, crossProduct, | ||
| &nFeatures); |
There was a problem hiding this comment.
Suggested change
| transb = 'T'; | |
| alpha = 1.0; | |
| beta = 1.0; | |
| blasInst.xgemm(&transa, &transb, &nFeatures, &nFeatures, &nVectors, &alpha, data, &nFeatures, data, &nFeatures, &beta, crossProduct, | |
| &nFeatures); | |
| { | |
| constexpr char transb = 'T'; | |
| constexpr double alpha = 1.0; | |
| constexpr double beta = 1.0; | |
| blasInst.xgemm(&transa, &transb, &nFeatures, &nFeatures, &nVectors, &alpha, data, &nFeatures, data, &nFeatures, &beta, crossProduct, | |
| &nFeatures); | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
Notation
Data is a matrix$X\in\mathbb{R}^{p\times n}$ . Each column is a $p$ -dimensional vector sampled independently. The matrix $X$ is assumed to be stored in column-major fashion.
xcp
Matrix of cross products is a$p\times p$ matrix whose $(i,j)$th entry is given by
$$C_{ij} = \sum_{k=1}^n(x_{ik}-\mu_i)(x_{jk}-\mu_j).$$ $x_{ik}$ is the $i$th component of the $k$th random vector.
Here
The implementation requires this to be computed in batches of$X$ , where the raw sum is passed between the batch computations. This can be done using the following:
Let$\mu'$ be the mean of the data until the previous batch and let $\mu$ be the mean including the data from current batch. Let $n'$ be the total number of data points until the previous batch and $S'$ be their sum (i.e., $\mu' = S'/n'$ and $\mu = S/n$ ).
Then the$(i,j)$ entry of $C$ can be computed using
$$C_{ij}\leftarrow C_{ij}+(\mu_j'-\mu_j)S_i' + (\mu_i'-\mu_i)S_j' + (\mu_i\mu_j-\mu_i'\mu_j')n' + \sum_{k=n'+1}^n(x_{ik}-\mu_i)(x_{jk}-\mu_j).$$
This can be simplified to
$$C_{ij} \leftarrow C_{ij} + \frac{S_i' S_j'}{n'} - \frac{S_i S_j}{n} + \sum_{k=n'+1}^n x_{ik}x_{jk}$$
$$C \leftarrow C + \frac{S'S'^t}{n'} - \frac{SS^t}{n} + XX^t$$
or
The level3 BLAS routine GEMM is used to compute the matrix products.