chol: zero trailing with pivot=TRUE + update help page

git-svn-id: https://svn.r-project.org/R/trunk@89125 00db46b3-68df-0310-9c12-caf00c1e9a41

Author: pd

doc/NEWS.Rd

@@ -306,6 +306,11 @@
 
   \subsection{BUG FIXES}{
     \itemize{
+      \item \code{chol(x, pivot=TRUE)} now zeroes the trailing submatrix when
+      the algorithm stops early. Following the help page instructions on
+      how to reconstruct the matrix in postive semidefinite cases gave
+      incorrect results if the rank was 2 or more less than the dimension.
+
       \item Profiling \code{"mle"} objects (\pkg{stats4} package) could fail
       because profile MLE used starting value from full MLE. Now starts from profile MLE
       from previous step.

src/library/base/man/chol.Rd

@@ -42,15 +42,22 @@ chol(x, \dots)
   error occurs.  If \code{x} is positive semi-definite (i.e., some zero
   eigenvalues) an error will also occur as a numerical tolerance is used.
 
-  If \code{pivot = TRUE}, then the Cholesky decomposition of a positive
+  If \code{pivot = TRUE}, then the Cholesky decomposition of a (rearranged)
+  positive
   semi-definite \code{x} can be computed.  The rank of \code{x} is
   returned as \code{attr(Q, "rank")}, subject to numerical errors.
   The pivot is returned as \code{attr(Q, "pivot")}.  It is no longer
   the case that \code{t(Q) \%*\% Q} equals \code{x}.  However, setting
   \code{pivot <- attr(Q, "pivot")} and \code{oo <- order(pivot)}, it
   is true that \code{t(Q[, oo]) \%*\% Q[, oo]} equals \code{x},
   or, alternatively, \code{t(Q) \%*\% Q} equals \code{x[pivot,
-    pivot]}.  See the examples.
+    pivot]}. Notice also, that \code{Q[,oo]} is typically not triangular.
+   See the examples.
+
+  NOTE: For versions of R up to 4.5.2, the above wasn't quite true.
+  You also needed \code{rank <- attr(Q, "rank")} and 
+  \code{t(Q[1:rank, oo, drop=FALSE]) \%*\% Q[1:rank, oo, drop=FALSE]} 
+  and similar for the alternative version.    
   
   The value of \code{tol} is passed to LAPACK, with negative values
   selecting the default tolerance of (usually) \code{nrow(x) *
@@ -99,8 +106,8 @@ crossprod(cm)  #-- = 'm'
 
 # now for something positive semi-definite
 x <- matrix(c(1:5, (1:5)^2), 5, 2)
-x <- cbind(x, x[, 1] + 3*x[, 2])
-colnames(x) <- letters[20:22]
+x <- cbind(x, x[, 1] + 3*x[, 2], 2*x[, 1] + x[, 2])
+colnames(x) <- letters[20:23]
 m <- crossprod(x)
 qr(m)$rank # is 2, as it should be
 
@@ -118,6 +125,13 @@ crossprod(Q[, order(pivot)]) # recover m
 try(chol(m))  # fails
 (Q <- chol(m, pivot = TRUE)) # warning
 crossprod(Q)  # not equal to m
+
+## another example, this time with positive, negative, and zero eigenvalues
+## notice that this (and the previous example) gets the rank wrong
+(m <- matrix(c(1,1,2, 1,1,2, 2,2,1), 3, 3))
+chol(m, pivot=TRUE)
+eigen(m)$values
 }
+
 \keyword{algebra}
 \keyword{array}

src/modules/lapack/Lapack.c

@@ -1118,6 +1118,10 @@ static SEXP La_chol(SEXP A, SEXP pivot, SEXP stol)
 		error(_("argument %d of Lapack routine %s had invalid value"),
 		      -info, "dpstrf");
 	}
+ 	/* zero remaining upper triangle if rank < m */
+	for (int j = rank ; j < m ; j++)
+            for (int i = rank ; i <= j ; i++)
+                REAL(ans) [i + N * j] = 0. ;
 	setAttrib(ans, install("pivot"), piv);
 	SEXP s_rank = install("rank");
 	setAttrib(ans, s_rank, ScalarInteger(rank));