The R Project SVN R

Rev

Rev 89125 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 89125 Rev 89181
Line 52... Line 52...
52 
  is true that \code{t(Q[, oo]) \%*\% Q[, oo]} equals \code{x},
 52 
  is true that \code{t(Q[, oo]) \%*\% Q[, oo]} equals \code{x},
53 
  or, alternatively, \code{t(Q) \%*\% Q} equals \code{x[pivot,
 53 
  or, alternatively, \code{t(Q) \%*\% Q} equals \code{x[pivot,
54 
    pivot]}. Notice also, that \code{Q[,oo]} is typically not triangular.
 54 
    pivot]}. Notice also, that \code{Q[,oo]} is typically not triangular.
55 
   See the examples.
 55 
   See the examples.
56 
 
 56 
 
57 
  NOTE: For versions of R up to 4.5.2, the above wasn't quite true.
 57 
  NOTE: For versions of R up to 4.5.x, the above wasn't quite true,
58 
  You also needed \code{rank <- attr(Q, "rank")} and
58 
  since LAPACK does not zero out the rows of \code{Q} beyond the rank.
59 
  \code{t(Q[1:rank, oo, drop=FALSE]) \%*\% Q[1:rank, oo, drop=FALSE]}
59 
  This is now done in R, but if you want code to be compatible with
60 
  and similar for the alternative version.
60 
  earlier versions of R, you may want to ensure that you only use
- 
 
 61 
  the first \code{rank} rows.
61
62
62 
  The value of \code{tol} is passed to LAPACK, with negative values
 63 
  The value of \code{tol} is passed to LAPACK, with negative values
63 
  selecting the default tolerance of (usually) \code{nrow(x) *
 64 
  selecting the default tolerance of (usually) \code{nrow(x) *
64 
  .Machine$double.neg.eps * max(diag(x))}.  The algorithm terminates once
 65 
  .Machine$double.neg.eps * max(diag(x))}.  The algorithm terminates once
65 
  the pivot is less than \code{tol}.
 66 
  the pivot is less than \code{tol}.
Line 125... Line 126...
125 
try(chol(m))  # fails
 126 
try(chol(m))  # fails
126 
(Q <- chol(m, pivot = TRUE)) # warning
 127 
(Q <- chol(m, pivot = TRUE)) # warning
127 
crossprod(Q)  # not equal to m
 128 
crossprod(Q)  # not equal to m
128 
 
 129 
 
129 
## another example, this time with positive, negative, and zero eigenvalues
 130 
## another example, this time with positive, negative, and zero eigenvalues
130 
## notice that this (and the previous example) gets the rank wrong
 131 
## notice that this (and the previous example) gets the rank wrong - it
- 
 
 132 
## returns the rank of an approximating positive semidefinite matrix,
- 
 
 133 
## effectively discounting negative eigenvalues.
131 
(m <- matrix(c(1,1,2, 1,1,2, 2,2,1), 3, 3))
 134 
(m <- matrix(c(1,1,2, 1,1,2, 2,2,1), 3, 3))
132 
chol(m, pivot=TRUE)
 135 
(Q <- chol(m, pivot=TRUE))
133 
eigen(m)$values
 136 
eigen(m)$values
- 
 
 137 
 
- 
 
 138 
## conjecture: indefinite cases should be detectable by comparing the
- 
 
 139 
## diagonal of the reconstructed matrix to the original
- 
 
 140 
piv <- attr(Q, "pivot")
- 
 
 141 
(D <- diag(crossprod(Q))) ## faster: colSums(Q^2)
- 
 
 142 
diag(m)[piv] - D
- 
 
 143 
 
134 
}
 144 
}
135 
 
 145 
 
136 
\keyword{algebra}
 146 
\keyword{algebra}
137 
\keyword{array}
 147 
\keyword{array}