<p><div><div><div><div data-width="100" data-show-faces="false" data-send="true" data-layout="button_count" data-height="21" data-href="https://www.r-bloggers.com/linear-algebra-do-not-scale-your-matrix/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://privefl.github.io/blog/(Linear-Algebra)-Do-not-scale-your-matrix/"> Florian Privé</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<section>
<p>In this post, I will show you that you generally don’t need to explicitly scale a matrix. Maybe you wanted to know more about WHY matrices should be scaled when doing linear algebra. I will remind about that in the beginning but the rest will focus on HOW to not explicitly scale matrices. We will apply our findings to the computation of Principal Component Analysis (PCA) and then Pearson correlation at the end.</p>
<div>
<h2>WHY scaling matrices?</h2>
<p>Generally, if you don’t center columns of a matrix before PCA, you end up with the loadings of PC1 being the column means, which is not of must interest.</p>
<div>
<pre><code>n &lt;-<span> </span><span>100</span>; m &lt;-<span> </span><span>10</span>
a &lt;-<span> </span><span>matrix</span>(<span>0</span>, n, m); a[] &lt;-<span> </span><span>rnorm</span>(<span>length</span>(a))
a &lt;-<span> </span><span>sweep</span>(a, <span>2</span>, <span>1</span>:m, <span>'+'</span>)
<span>colMeans</span>(a)</code></pre>
</div>
<pre><code>##  [1]  0.9969826  1.9922219  3.0894905  3.9680835  4.9433477  6.0007352
##  [7]  6.9524794  8.0898546  8.9506657 10.0354630</code></pre>
<div>
<pre><code>pca &lt;-<span> </span><span>prcomp</span>(a, <span>center =</span> <span>FALSE</span>)
<span>cor</span>(pca$rotation[, <span>1</span>], <span>colMeans</span>(a))</code></pre>
</div>
<pre><code>## [1] 0.9999991</code></pre>
<p>Now, say you have centered column or your matrix, do you also need to scale them? That is to say, do they need to have the same norm or standard deviation?</p>
<p>PCA consists in finding an orthogonal basis that maximizes the variation in the data you analyze. So, if there is a column with much more variation than the others, it will probably end up being PC1, which is not of must interest.</p>
<div>
<pre><code>n &lt;-<span> </span><span>100</span>; m &lt;-<span> </span><span>10</span>
a &lt;-<span> </span><span>matrix</span>(<span>0</span>, n, m); a[] &lt;-<span> </span><span>rnorm</span>(<span>length</span>(a))
a[, <span>1</span>] &lt;-<span> </span><span>100</span> *<span> </span>a[, <span>1</span>]
<span>apply</span>(a, <span>2</span>, sd)</code></pre>
</div>
<pre><code>##  [1] 90.8562836  0.8949268  0.9514208  1.0437561  0.9995103  1.0475651
##  [7]  0.9531466  0.9651383  0.9804789  0.9605121</code></pre>
<div>
<pre><code>pca &lt;-<span> </span><span>prcomp</span>(a, <span>center =</span> <span>TRUE</span>)
pca$rotation[, <span>1</span>]</code></pre>
</div>
<pre><code>##  [1]  9.999939e-01  6.982439e-04 -5.037773e-04 -9.553513e-05  8.937869e-04
##  [6] -1.893732e-03 -3.053232e-04  1.811950e-03 -1.658014e-03  9.937442e-04</code></pre>
<p>Hope I convinced you on WHY it is important to scale matrix columns before doing PCA. I will now show you <strong>HOW not to do it explictly</strong>.</p>
</div>
<div>
<h2>Reimplementation of PCA</h2>
<p>In this part, I will show you the basic code if you want to reimplement PCA yourself. It will serve as a basis to show you how to do linear algebra on a scaled matrix, without explicitly scaling the matrix.</p>
<div>
<pre><code><span># True one</span>
n &lt;-<span> </span><span>100</span>; m &lt;-<span> </span><span>10</span>
a &lt;-<span> </span><span>matrix</span>(<span>0</span>, n, m); a[] &lt;-<span> </span><span>rnorm</span>(<span>length</span>(a))
pca &lt;-<span> </span><span>prcomp</span>(a, <span>center =</span> <span>TRUE</span>, <span>scale. =</span> <span>TRUE</span>)
<span># DIY</span>
a.scaled &lt;-<span> </span><span>scale</span>(a, <span>center =</span> <span>TRUE</span>, <span>scale =</span> <span>TRUE</span>)
K &lt;-<span> </span><span>crossprod</span>(a.scaled)
K.eigs &lt;-<span> </span><span>eigen</span>(K, <span>symmetric =</span> <span>TRUE</span>)
v &lt;-<span> </span>K.eigs$vectors
PCs &lt;-<span> </span>a.scaled %*%<span> </span>v
<span># Verif, recall that PCs can be opposites between runs</span>
<span>plot</span>(v, pca$rotation)</code></pre>
</div>
<p><img src="https://i2.wp.com/privefl.github.io/blog/knitr_files/post-scale-matrix_files/figure-html/unnamed-chunk-3-1.png?w=456&ssl=1" data-recalc-dims="1"></p>
<div>
<pre><code><span>all.equal</span>(<span>sqrt</span>(K.eigs$values), <span>sqrt</span>(n -<span> </span><span>1</span>) *<span> </span>pca$sdev)</code></pre>
</div>
<pre><code>## [1] TRUE</code></pre>
<div>
<pre><code><span>plot</span>(PCs, pca$x)</code></pre>
</div>
<p><img src="https://i0.wp.com/privefl.github.io/blog/knitr_files/post-scale-matrix_files/figure-html/unnamed-chunk-3-2.png?w=456&ssl=1" data-recalc-dims="1"></p>
</div>
<div>
<h2>Linear algebra behind the previous implementation</h2>
<p>Suppose <span>\(m &lt; n\)</span> (<span>\(m\)</span> is the number of columns and <span>\(n\)</span> is the number of rows). Let us denote <span>\(\tilde{X}\)</span> the scaled matrix. A partial singular value decomposition of <span>\(\tilde{X}\)</span> is <span>\(\tilde{X} \approx U \Delta V^T\)</span> where <span>\(U\)</span> is an <span>\(n \times K\)</span> matrix such that <span>\(U^T U = I_K\)</span>, <span>\(\Delta\)</span> is a <span>\(K \times K\)</span> diagonal matrix and <span>\(V\)</span> is an <span>\(m \times K\)</span> matrix such that <span>\(V^T V = I_K\)</span>. Taking <span>\(K = m\)</span>, you end up with <span>\(\tilde{X} = U \Delta V^T\)</span>.</p>
<p><span>\(U \Delta\)</span> are the scores (PCs) of the PCA and <span>\(V\)</span> are the loadings (rotation coefficients). <span>\(K = \tilde{X}^T \tilde{X} = (U \Delta V^T)^T \cdot U \Delta V^T = V \Delta U^T U \Delta V^T = V \Delta^2 V^T\)</span>. So, when doing the eigen decomposition of K, you get <span>\(V\)</span> and <span>\(\Delta^2\)</span> because <span>\(K V = V \Delta^2\)</span>. For getting the scores, you then compute <span>\(\tilde{X} V = U \Delta\)</span>.</p>
<p>These are exactly the steps implemented above.</p>
</div>
<div>
<h2>Implicit scaling of the matrix</h2>
<p>Do you know the matrix formulation of column scaling? <span>\(\tilde{X} = C_n X S\)</span> where <span>\(C_n = I_n – \frac{1}{n} 1_n 1_n^T\)</span> is the <a target="_blank" rel="nofollow" href="https://en.wikipedia.org/wiki/Centering_matrix">centering matrix</a> and <span>\(S\)</span> is an <span>\(m \times m\)</span> diagonal matrix with the scaling coefficients (typically, <span>\(S_{j,j} = 1 / \text{sd}_j\)</span>).</p>
<div>
<pre><code><span># Let's verify</span>
sds &lt;-<span> </span><span>apply</span>(a, <span>2</span>, sd)
a.scaled2 &lt;-<span> </span>(<span>diag</span>(n) -<span> </span><span>tcrossprod</span>(<span>rep</span>(<span>1</span>, n)) /<span> </span>n) %*%<span> </span>a %*%<span> </span><span>diag</span>(<span>1</span> /<span> </span>sds)
<span>all.equal</span>(a.scaled2, a.scaled, <span>check.attributes =</span> <span>FALSE</span>)</code></pre>
</div>
<pre><code>## [1] TRUE</code></pre>
<p>In our previous implementation, we computed <span>\(\tilde{X}^T \tilde{X}\)</span> and <span>\(\tilde{X} V\)</span>. We are going to compute them again, without explicitly scaling the matrix.</p>
<div>
<h3>Product</h3>
<p>Let us begin by something easy: <span>\(\tilde{X} V = C_n X S V = C_n (X (S V))\)</span>. So, you can compute <span>\(\tilde{X} V\)</span> without explicitly scaling <span>\(X\)</span>. Let us verify:</p>
<div>
<pre><code>SV &lt;-<span> </span>v /<span> </span>sds
XSV &lt;-<span> </span>a %*%<span> </span>SV
CXSV &lt;-<span> </span><span>sweep</span>(XSV, <span>2</span>, <span>colMeans</span>(XSV), <span>'-'</span>)
<span>all.equal</span>(CXSV, PCs)</code></pre>
</div>
<pre><code>## [1] TRUE</code></pre>
</div>
<div>
<h3>Self cross-product</h3>
<p>A little more tricky: <span>\(\tilde{X}^T \tilde{X} = (C_n X S)^T \cdot C_n X S = S^T X^T C_n X S\)</span> (<span>\(C_n^2 = C_n\)</span> is intuitive because centering an already centered matrix doesn’t change it).</p>
<p><span>\(\tilde{X}^T \tilde{X} = S^T X^T (I_n – \frac{1}{n} 1_n 1_n^T) X S = S^T (X^T X – X^T (\frac{1}{n} 1_n 1_n^T) X) S = S^T (X^T X – \frac{1}{n} s_X * s_X^T) S\)</span> where <span>\(s_X\)</span> is the vector of column sums of X.</p>
<p>Let us verify with a rough implementation:</p>
<div>
<pre><code>sx &lt;-<span> </span><span>colSums</span>(a)
K2 &lt;-<span> </span>(<span>crossprod</span>(a) -<span> </span><span>tcrossprod</span>(sx) /<span> </span>n) /<span> </span><span>tcrossprod</span>(sds)
<span>all.equal</span>(K, K2)</code></pre>
</div>
<pre><code>## [1] TRUE</code></pre>
</div>
</div>
<div>
<h2>Conclusion</h2>
<p>We have recalled some steps about the computation of Principal Component Analysis on a scaled matrix. We have seen how to compute the different steps of the implementation without having to explicitly scale the matrix. This “implicit” scaling can be quite useful if you manipulate very large matrices because you are not copying the matrix nor making useless computation. In my next post, I will present to you my new package that uses this trick to make a lightning partial SVD on very large matrices.</p>
</div>
<div>
<h2>Appendix: application to Pearson correlation</h2>
<p>Pearson correlation is merely a self cross-product on a centered and normalized (columns with unit norm) matrix. Let us just implement that with our new trick.</p>
<div>
<pre><code><span>#include &lt;Rcpp.h></span>
<span>using</span> <span>namespace</span> Rcpp;

<span>// [[Rcpp::export]]</span>
NumericMatrix&amp; correlize(NumericMatrix&amp; mat,
                         <span>const</span> NumericVector&amp; shift,
                         <span>const</span> NumericVector&amp; scale) {
  
  <span>int</span> n = mat.nrow();
  <span>int</span> i, j;
  
  <span>for</span> (j = <span>0</span>; j &lt; n; j++) {
    <span>for</span> (i = <span>0</span>; i &lt; n; i++) {
      <span>// corresponds to &quot;- \frac{1}{n} s_X * s_X^T&quot;</span>
      mat(i, j) -= shift(i) * shift(j);
      <span>// corresponds to &quot;S^T (...) S&quot;</span>
      mat(i, j) /= scale(i) * scale(j);
    }
  }
  
  <span>return</span> mat;
}</code></pre>
</div>
<div>
<pre><code>cor3 &lt;-<span> </span>function(mat) {
  sums &lt;-<span> </span><span>colSums</span>(mat) /<span> </span><span>sqrt</span>(<span>nrow</span>(mat))
  corr &lt;-<span> </span><span>crossprod</span>(mat)
  diags &lt;-<span> </span><span>sqrt</span>(<span>diag</span>(corr) -<span> </span>sums^<span>2</span>)
  <span>correlize</span>(corr, <span>shift =</span> sums, <span>scale =</span> diags)
}

a &lt;-<span> </span><span>matrix</span>(<span>0</span>, <span>1000</span>, <span>1000</span>); a[] &lt;-<span> </span><span>rnorm</span>(<span>length</span>(a))
<span>all.equal</span>(<span>cor3</span>(a), <span>cor</span>(a))</code></pre>
</div>
<pre><code>## [1] TRUE</code></pre>
<div>
<pre><code><span>library</span>(microbenchmark)
<span>microbenchmark</span>(
  <span>cor3</span>(a),
  <span>cor</span>(a),
  <span>times =</span> <span>20</span>
)</code></pre>
</div>
<pre><code>## Unit: milliseconds
##     expr       min        lq      mean    median        uq       max neval
##  cor3(a)  39.38138  39.67276  40.68596  40.09893  40.65623  46.46785    20
##   cor(a) 635.74350 637.33605 639.34810 638.09980 639.36110 651.61876    20</code></pre>
</div>
</section>
	<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="https://privefl.github.io/blog/(Linear-Algebra)-Do-not-scale-your-matrix/"> Florian Privé</a></strong>.</div>
<hr>
<a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers.com</a> offers <strong><a rel="nofollow" href="https://feedburner.google.com/fb/a/mailverify?uri=RBloggers">daily e-mail updates</a></strong> about <a title="The R Project for Statistical Computing" rel="nofollow" href="https://www.r-project.org/">R</a> news and <a title="R tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/tutorial">tutorials</a> on topics such as: <a title="Data science" rel="nofollow" href="https://www.r-bloggers.com/search/data%20science">Data science</a>, <a title="Big Data" rel="nofollow" href="https://www.r-bloggers.com/search/Big%20Data">Big Data, <a title="R jobs" rel="nofollow" href="https://www.r-users.com/">R jobs</a>, visualization (<a title="ggplot and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/ggplot2">ggplot2</a>, <a title="Boxplots using lattice and ggplot2 tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/boxplot">Boxplots</a>, <a title="Maps and gis" rel="nofollow" href="https://www.r-bloggers.com/search/map">maps</a>, <a title="Animation in R" rel="nofollow" href="https://www.r-bloggers.com/search/animation">animation</a>), programming (<a title="RStudio IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/RStudio">RStudio</a>, <a title="Sweave and literate programming" rel="nofollow" href="https://www.r-bloggers.com/search/sweave">Sweave</a>, <a title="LaTeX in R" rel="nofollow" href="https://www.r-bloggers.com/search/LaTeX">LaTeX</a>, <a title="SQL and databases" rel="nofollow" href="https://www.r-bloggers.com/search/SQL">SQL</a>, <a title="Eclipse IDE for R" rel="nofollow" href="https://www.r-bloggers.com/search/eclipse">Eclipse</a>, <a title="git and github, Version Control System" rel="nofollow" href="https://www.r-bloggers.com/search/git">git</a>, <a title="Large data in R using Hadoop" rel="nofollow" href="https://www.r-bloggers.com/search/hadoop">hadoop</a>, <a title="Web Scraping of google, facebook, yahoo, twitter and more using R" rel="nofollow" href="https://www.r-bloggers.com/search/Web+Scraping">Web Scraping</a>) statistics (<a title="Regressions and ANOVA analysis tutorials" rel="nofollow" href="https://www.r-bloggers.com/search/regression">regression</a>, <a title="principal component analysis tutorial" rel="nofollow" href="https://www.r-bloggers.com/search/PCA">PCA</a>, <a title="Time series" rel="nofollow" href="https://www.r-bloggers.com/search/time+series">time series</a>, <a title="finance trading" rel="nofollow" href="https://www.r-bloggers.com/search/trading">trading</a>) and more...

	<img height="1" alt="" width="1" src="http://feeds.feedburner.com/~r/RBloggers/~4/5-f6GoZv23w"></a></div>