<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/the-statisticians-apprentice-an-introduction-to-the-swp-operator/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="http://localhost:4000/2016/12/13/swp.html"> MeanMean</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>

<p></p>
<p>
The sweep operator as defined in (Dempster, 1969), commonly referred to as the SWP operator, is a useful tool for a computational statistician working with covariance matrices. In particular, the SWP operator allows a statistician to quickly regress all variables against one specified variable, obtaining OLS estimates for regression coefficients and variances in a single application. Subsequent applications of the SWP operator allows for regressing against more variables.
</p>
<p>
In this blog post, I will define the sweep operator, provide an application of the sweep data to simulated data, and provide some references for further study. Examples will be provided in R through the <code>ISR3</code> package (Lisic, 2016).
</p>
<p>
The basic SWP operator is parameterized by a matrix and a set of indices associated with rows and columns of the matrix being swept. E.g. to sweep the matrix <span></span> by the <span></span> index, such that <span></span>, the operator has the form <span></span>. The operator may also be written more descriptively by listing the elements of <span></span> as <span></span>. The application of the SWP operator to the <span></span> row and column of <span></span> produces a new matrix with elements
</p>
<div><span></span></div>

<p>
As an example, consider the application of the SWP operator to the first row and column of the <span></span> covariance matrix for the random variables <span></span>,
</p>
<div><span></span></div>
<p>resulting in </p>
<div><span></span></div>

<p>
How is this useful? Well, if we substitute the <span></span> covariance matrix <span></span> with the sample covariance matrix, where <span></span> is an estimator of <span></span> and
</p>
<div><span></span></div>
<p> is an estimator of <span></span>, then by sweeping by the <span></span> index, <span></span>, we can get the simple linear regression (SLR) estimates for <span></span> for all <span></span>:</p>
<div><span></span></div>
<p>where <span></span> is the simple linear model estimator of the regression coefficient for the linear model <span></span> for <span></span>, <span></span> is the inverse of the covariance matrix of <span></span> unscaled by <span></span>, and <span></span> is the sample covariance of the residuals from the regression models on <span></span> unscaled by <span></span>.</p>
<h2>Example</h2>
<p>
Lets try this with some real data. I have used the SWP function included with the <code>ISR3</code> package.
</p>
<div>
<div>
<div>
<pre><span>library</span><span>(ISR3)</span>

<span>set.seed</span><span>(</span><span>12</span><span>)</span>
<span>n</span> <span>&lt;-</span> <span>100</span>
<span>p</span> <span>&lt;-</span> <span>3</span>

<span># generate a positive definite matrix for covariance</span>
<span>Sigma</span> <span>&lt;-</span> <span>rWishart</span><span>(</span><span>1</span><span>,p</span><span>+</span><span>1</span><span>,</span><span>diag</span><span>(p))[,,</span><span>1</span><span>]</span>
<span>Sigma_inv</span> <span>&lt;-</span> <span>chol2inv</span><span>(</span><span>chol</span><span>(Sigma))</span>
<span>Sigma_inv_chol</span> <span>&lt;-</span> <span>chol</span><span>(Sigma_inv)</span>

<span># generate 'n' multivate normal deviates</span>
<span>X</span> <span>&lt;-</span> <span>Sigma_inv_chol</span> <span>%*%</span> <span>matrix</span><span>(</span><span>rnorm</span><span>(n</span><span>*</span><span>p),</span><span>nrow</span><span>=p)</span>
<span>X</span> <span>&lt;-</span> <span>t</span><span>(X)</span>
<span>colnames</span><span>(X)</span> <span>&lt;-</span> <span>sprintf</span><span>(</span><span>&quot;X_%d&quot;</span><span>,</span><span>1</span><span>:</span><span>p)</span>

<span>XX</span> <span>&lt;-</span> <span>t</span><span>(X)</span> <span>%*%</span> <span>X</span>

<span>#Sweep by the first row/column of XX</span>
<span>SWP_1</span> <span>&lt;-</span> <span>SWP</span><span>(XX,</span><span>1</span><span>)</span>

<span>SWP_1</span>
</pre>
</div>
<div>
<pre>##              X_1        X_2         X_3
## X_1 -0.002381092  0.1375626 -0.01674788
## X_2  0.137562633 50.0987932 -4.95181624
## X_3 -0.016747875 -4.9518162 41.73979419
</pre>
</div>
</div>
</div>
<p>We can compare this output with the results from the simple linear models <span></span> and <span></span>: </p>
<div>
<div>
<div>
<pre><span># linear models</span>
<span>fit_21</span> <span>&lt;-</span> <span>lm</span><span>(X_2</span> <span>~ -</span><span>1</span> <span>+</span> <span>X_1,</span> <span>data</span><span>=</span><span>as.data.frame</span><span>(X))</span>
<span>fit_31</span> <span>&lt;-</span> <span>lm</span><span>(X_3</span> <span>~ -</span><span>1</span> <span>+</span> <span>X_1,</span> <span>data</span><span>=</span><span>as.data.frame</span><span>(X))</span>

<span># SLR coefficents for X_1</span>
<span>print</span><span>(fit_21</span><span>$</span><span>coefficients)</span>
</pre>
</div>
<div>
<pre>##       X_1 
## 0.1375626
</pre>
</div>
<div>
<pre><span>print</span><span>(fit_31</span><span>$</span><span>coefficients)</span>
</pre>
</div>
<div>
<pre>##         X_1 
## -0.01674788
</pre>
</div>
<div>
<pre><span># sum of squares</span>
<span>print</span><span>(</span><span>sum</span><span>(fit_21</span><span>$</span><span>residuals</span><span>^</span><span>2</span><span>))</span>
</pre>
</div>
<div>
<pre>## [1] 50.09879
</pre>
</div>
<div>
<pre><span>print</span><span>(</span><span>sum</span><span>(fit_31</span><span>$</span><span>residuals</span><span>^</span><span>2</span><span>))</span>
</pre>
</div>
<div>
<pre>## [1] 41.73979
</pre>
</div>
<div>
<pre><span># covariance between X_2 and X_3 given X_1</span>
<span>print</span><span>(</span><span>sum</span><span>(fit_21</span><span>$</span><span>residuals</span> <span>*</span> <span>fit_31</span><span>$</span><span>residuals))</span>
</pre>
</div>
<div>
<pre>## [1] -4.951816
</pre>
</div>
<div>
<pre><span># negative inverse of the covariance matrix of X_1</span>
<span>print</span><span>(</span> <span>-</span><span>1</span><span>/</span><span>((n</span><span>-</span><span>1</span><span>)</span><span>*</span><span>var</span><span>(X[,</span><span>1</span><span>])) )</span>
</pre>
</div>
<div>
<pre>## [1] -0.002381293
</pre>
</div>
</div>
</div>
<p>Here we see that the coefficients from <code>lm</code> match up to our swept results. Likewise, the diagonal terms for the second and third column are equivalent to the sum of squares or the residuals. The remaining unmatched term is just the negative inverse of the covariance matrix of <span></span> unscaled by <span></span>. </p>
<p>
A subsequent sweep on the second rows and columns produces the matrix,
</p>
<div><span></span></div>
<p>where <span></span>, <span></span> is now the inverse covariance matrix of <span></span> and <span></span> and <span></span> is the matrix </p>
<div><span></span></div>
<p>Again, we can quickly check the results by comparing the output from <code>lm</code> to our SWP operator. </p>
<p>
A final sweep by the third row and column produces <span></span>, the negative inverse of the unscaled covariance matrix of <span></span>. This can be quickly verified through <code>R</code>:
</p>
<div>
<div>
<div>
<pre><span>SWP</span><span>(XX,</span><span>1</span><span>:</span><span>3</span><span>)</span> <span>+</span> <span>chol2inv</span><span>(</span><span>chol</span><span>(XX))</span>
</pre>
</div>
<div>
<pre>##              X_1          X_2          X_3
## X_1 0.000000e+00 0.000000e+00 1.355253e-20
## X_2 0.000000e+00 3.469447e-18 8.673617e-19
## X_3 1.355253e-20 8.673617e-19 0.000000e+00
</pre>
</div>
</div>
</div>
<h2>Undo!</h2>
<p>
We can also undo sweeps with the reverse sweep (RSWP) operator through the following element-wise operations
</p>
<div><span></span></div>
<p>Caution should be used when reverse sweeping, as in application <span></span> is not necessarily equal to <span></span> due to floating point error. In particular when the matrix is ill-conditioned there may be considerable departure between the calculated identity and the original matrix. </p>
<h2>What Next?</h2>
<p>
A few good references on the subject beyond (Dempster, 1969) can be found in (Goodknight, 1979) and (Little and Rubin, 2014). There are a few slight differences between these implementations, but all generally do the same thing.
</p>
<p>One function that doesn’t do the same thing is the <code>sweep</code> function in <code>R</code>. This function can be used to write a SWP implementation, but doesn’t sweep by itself.
</p>
<p>
Code examples can be found in the <code>ISR3</code> package on <code>CRAN</code> and in lots of other <a target="_blank" rel="nofollow" href="https://www.rdocumentation.org/search?q=sweep">packages</a>.<br>
Examples in <code>C</code> can be found in <code>ISR3</code> and as an internal function in the <code>MNP</code> package (Imai and Van Dyk, 2013).
</p>
<p>
Enjoy your sweeping, next time I’ll go over sweeping with the generalized inverse!
</p>
<p><img title="broom by  ericlemerdy ( https://openclipart.org/user-detail/ericlemerdy )" alt="broom" src="https://openclipart.org/image/300px/svg_to_png/223800/1438128529.png&amp;disposition=attachment"></p>
<h2>References</h2>
<p>
Dempster, A.P. (1969). <i>Elements of continuous multivariate analysis</i>. Reading, MA: Addison-Wesley.
</p>
<p>
Goodnight, J. H. (1979). A tutorial on the SWEEP operator. The American Statistician, 33(3), 149-158.
</p>
<p>
Imai, K., and Van Dyk, D. A. (2013). MNP: R Package for fitting the Multinomial Probit Model, R package version 2.6-4, <a target="_blank" rel="nofollow" href="https://cran.r-project.org/package=MNP">https://CRAN.R-project.org/package=MNP</a>
</p>
<p>
Lisic, J. J. (2016). ISR3: Iterative Sequential Regression, R package version 0.98, <a target="_blank" rel="nofollow" href="https://cran.r-project.org/package=ISR3">https://CRAN.R-project.org/package=ISR3</a>
</p>
<p>
Little, R. J., and Rubin, D. B. (2014). <i>Statistical analysis with missing data</i>. John Wiley &amp; Sons.</p>
<p>
<div>
<div>To <strong>leave a comment</strong> for the author, please follow the link and comment on their blog: <strong><a href="http://localhost:4000/2016/12/13/swp.html"> MeanMean</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/sjafjbGOw7I"></a></div>