SCM

SCM Repository

[matrix] Diff of /pkg/inst/doc/Comparisons.Rnw
ViewVC logotype

Diff of /pkg/inst/doc/Comparisons.Rnw

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 188, Mon May 31 18:50:01 2004 UTC revision 342, Mon Nov 15 13:55:52 2004 UTC
# Line 23  Line 23 
23  \end{abstract}  \end{abstract}
24  <<preliminaries, echo=FALSE>>=  <<preliminaries, echo=FALSE>>=
25  options(width=75)  options(width=75)
 sysgc.time <- function(expr, gcFirst = TRUE)  
 {  
     if (!exists("proc.time"))  
         return(rep(as.numeric(NA), 5))  
     loc.frame <- parent.frame()  
     if (gcFirst)  
         gc(FALSE)  
     on.exit(cat("Timing stopped at:", proc.time() - time, "\n"))  
     expr <- substitute(expr)  
     time <- proc.time()  
     eval(expr, envir = loc.frame)  
     new.time <- proc.time()  
     on.exit()  
     if (length(new.time) == 3)  
         new.time <- c(new.time, 0, 0)  
     if (length(time) == 3)  
         time <- c(time, 0, 0)  
     new.time - time  
 }  
26  @  @
27    
28  \section{Linear least squares calculations}  \section{Linear least squares calculations}
# Line 88  Line 69 
69  On modern computers this calculation is performed so quickly that it cannot  On modern computers this calculation is performed so quickly that it cannot
70  be timed accurately in \RR{}  be timed accurately in \RR{}
71  <<timedNaive>>=  <<timedNaive>>=
72  sysgc.time(solve(t(mm) %*% mm) %*% t(mm) %*% y)  system.time(solve(t(mm) %*% mm) %*% t(mm) %*% y, gc = TRUE)
73  @  @
74  and it provides essentially the same results as the standard  and it provides essentially the same results as the standard
75  \code{lm.fit} function that is called by \code{lm}.  \code{lm.fit} function that is called by \code{lm}.
# Line 96  Line 77 
77  dput(c(solve(t(mm) %*% mm) %*% t(mm) %*% y))  dput(c(solve(t(mm) %*% mm) %*% t(mm) %*% y))
78  dput(lm.fit(mm, y)$coefficients)  dput(lm.fit(mm, y)$coefficients)
79  @  @
80    %$
81    
82  \subsection{A large example}  \subsection{A large example}
83  \label{sec:largeLSQ}  \label{sec:largeLSQ}
# Line 108  Line 89 
89  library(Matrix)  library(Matrix)
90  data(mm, package = "Matrix")  data(mm, package = "Matrix")
91  data(y, package = "Matrix")  data(y, package = "Matrix")
92  mm = as(mm, "matrix")  mm <- as(mm, "matrix")
93  dim(mm)  dim(mm)
94  sysgc.time(naive.sol <- solve(t(mm) %*% mm) %*% t(mm) %*% y)  system.time(naive.sol <- solve(t(mm) %*% mm) %*% t(mm) %*% y, gc = TRUE)
95  @  @
96    
97  Because the calculation of a ``cross-product'' matrix, such as  Because the calculation of a ``cross-product'' matrix, such as
# Line 133  Line 114 
114    
115  Combining these optimizations we obtain  Combining these optimizations we obtain
116  <<crossKoenNg>>=  <<crossKoenNg>>=
117  sysgc.time(cpod.sol <- solve(crossprod(mm), crossprod(mm,y)))  system.time(cpod.sol <- solve(crossprod(mm), crossprod(mm,y)), gc = TRUE)
118  all.equal(naive.sol, cpod.sol)  all.equal(naive.sol, cpod.sol)
119  @  @
120    
# Line 142  Line 123 
123  naive calculation.  In fact, the entire crossprod solution is  naive calculation.  In fact, the entire crossprod solution is
124  faster than simply calculating $\bX\trans\bX$ the naive way.  faster than simply calculating $\bX\trans\bX$ the naive way.
125  <<xpxKoenNg>>=  <<xpxKoenNg>>=
126  sysgc.time(t(mm) %*% mm)  system.time(t(mm) %*% mm, gc = TRUE)
127  @  @
128    
129  \subsection{Least squares calculations with Matrix classes}  \subsection{Least squares calculations with Matrix classes}
# Line 157  Line 138 
138  use a Cholesky decomposition.  The Cholesky decomposition could be used  use a Cholesky decomposition.  The Cholesky decomposition could be used
139  but it is rather awkward  but it is rather awkward
140  <<naiveChol>>=  <<naiveChol>>=
141  sysgc.time(ch <- chol(crossprod(mm)))  system.time(ch <- chol(crossprod(mm)), gc = TRUE)
142  sysgc.time(chol.sol <- backsolve(ch,  system.time(chol.sol <- backsolve(ch,
143     forwardsolve(ch, crossprod(mm, y), upper = TRUE, trans = TRUE)))     forwardsolve(ch, crossprod(mm, y), upper = TRUE, trans = TRUE)),
144       gc = TRUE)
145  all.equal(chol.sol, naive.sol)  all.equal(chol.sol, naive.sol)
146  @  @
147    
# Line 172  Line 154 
154  Cholesky decomposition.  Cholesky decomposition.
155  <<MatrixKoenNg>>=  <<MatrixKoenNg>>=
156  data(mm, package = "Matrix")  data(mm, package = "Matrix")
157  mm = as(mm, "geMatrix")  mm <- as(mm, "geMatrix")
158  class(crossprod(mm))  class(crossprod(mm))
159  sysgc.time(Mat.sol <- solve(crossprod(mm), crossprod(mm, y)))  system.time(Mat.sol <- solve(crossprod(mm), crossprod(mm, y)), gc=TRUE)
160  all.equal(naive.sol, as(Mat.sol, "matrix"))  all.equal(naive.sol, as(Mat.sol, "matrix"))
161  @  @
162    
# Line 182  Line 164 
164  decomposition or factorization stores the resulting factorization with  decomposition or factorization stores the resulting factorization with
165  the original object so that it can be reused without recalculation.  the original object so that it can be reused without recalculation.
166  <<saveFactor>>=  <<saveFactor>>=
167  xpx = crossprod(mm)  xpx <- crossprod(mm)
168  xpy = crossprod(mm, y)  xpy <- crossprod(mm, y)
169  sysgc.time(solve(xpx, xpy))  system.time(solve(xpx, xpy), gc=TRUE)
170  sysgc.time(solve(xpx, xpy))  system.time(solve(xpx, xpy), gc=TRUE)
171  @  @
172    
173  The model matrix \code{mm} is sparse; that is, most of the elements of  The model matrix \code{mm} is sparse; that is, most of the elements of
# Line 195  Line 177 
177  <<SparseKoenNg>>=  <<SparseKoenNg>>=
178  data(mm, package = "Matrix")  data(mm, package = "Matrix")
179  class(mm)  class(mm)
180  sysgc.time(sparse.sol <- solve(crossprod(mm), crossprod(mm, y)))  system.time(sparse.sol <- solve(crossprod(mm), crossprod(mm, y)), gc=TRUE)
181  all.equal(naive.sol, as(sparse.sol, "matrix"))  all.equal(naive.sol, as(sparse.sol, "matrix"))
182  @  @
183    
# Line 205  Line 187 
187  difficult to determine the difference in the solution times.  difficult to determine the difference in the solution times.
188    
189  <<SparseSaveFactor>>=  <<SparseSaveFactor>>=
190  xpx = crossprod(mm)  xpx <- crossprod(mm)
191  xpy = crossprod(mm, y)  xpy <- crossprod(mm, y)
192  sysgc.time(solve(xpx, xpy))  system.time(solve(xpx, xpy), gc=TRUE)
193  sysgc.time(solve(xpx, xpy))  system.time(solve(xpx, xpy), gc=TRUE)
194  @  @
195    
196  \bibliography{Matrix}  \bibliography{Matrix}

Legend:
Removed from v.188  
changed lines
  Added in v.342

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge