# SCM Repository

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

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

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