kappa(<singular>, exact=TRUE) |--> Inf "always"

git-svn-id: https://svn.r-project.org/R/trunk@87289 00db46b3-68df-0310-9c12-caf00c1e9a41
r-devel · Nov 6, 2024 · c66ee94 · c66ee94
1 parent cbb481a
commit c66ee94
Show file tree

Hide file tree

Showing 4 changed files with 38 additions and 5 deletions.
diff --git a/doc/NEWS.Rd b/doc/NEWS.Rd
@@ -352,6 +352,15 @@
   }
 }
 
+\section{\Rlogo CHANGES IN R 4.4.2 patched}{
+  \subsection{BUG FIXES}{
+    \itemize{
+      \item \code{kappa(A, exact=TRUE)} returns \code{Inf} more generally,
+      fixing  \PR{18817} reported by \I{Mikael Jagan}.
+    }
+  }
+}
+
 \section{\Rlogo CHANGES IN R 4.4.2}{
   \subsection{C-LEVEL FACILITIES}{
     \itemize{

diff --git a/src/library/base/R/kappa.R b/src/library/base/R/kappa.R
@@ -1,7 +1,7 @@
 #  File src/library/base/R/kappa.R
 #  Part of the R package, https://www.R-project.org
 #
-#  Copyright (C) 1998-2023 The R Core Team
+#  Copyright (C) 1998-2024 The R Core Team
 #  Copyright (C) 1998 B. D. Ripley
 #
 #  This program is free software; you can redistribute it and/or modify
@@ -73,8 +73,8 @@ kappa.default <- function(z, exact = FALSE,
     nNorm <- is.null(norm)
     if(exact) {
 	if(nNorm || norm == "2") {
-	    s <- svd(z, nu = 0L, nv = 0L)$d
-	    max(s)/min(s[s > 0])
+	    s <- svd(z, nu = 0L, nv = 0L)$d # decreasing, non-negative
+	    if(s[1]) s[1]/s[length(s)] else Inf # when s is all zero
 	}
 	else {
 	    if(nNorm) norm <- "1"
@@ -115,8 +115,9 @@ kappa.qr <- function(z, ...)
         return(0)
     if(exact) {
         if(is.null(norm) || identical("2", norm)) { # 2-norm : *not* assuming 'triangular'
+            ## identically to kappa.default(z, exact=TRUE) :
             s <- svd(z, nu = 0L, nv = 0L)$d
-            max(s)/min(s[s > 0]) ## <==> kappa.default(z, exact=TRUE)
+            if(s[1]) s[1]/s[length(s)] else Inf
         }
         else norm(z, type=norm) * norm(solve(z), type=norm) # == kappa.default(z, exact=TRUE, norm=norm,..)
     }

diff --git a/src/library/base/man/kappa.Rd b/src/library/base/man/kappa.Rd
@@ -1,6 +1,6 @@
 % File src/library/base/man/kappa.Rd
 % Part of the R package, https://www.R-project.org
-% Copyright 1995-2023 R Core Team
+% Copyright 1995-2024 R Core Team
 % Copyright 2008-2010 The R Foundation
 % Distributed under GPL 2 or later
 
@@ -203,5 +203,15 @@ iX <- pinv(m79)
 kappa(m79, exact=TRUE, norm="1", inv_z = iX)
 kappa(m79, exact=TRUE, norm="M", inv_z = iX)
 kappa(m79, exact=TRUE, norm="I", inv_z = iX)
+
+## Using a more "accurate" than default inv_z [example by Cleve Moler]:
+A <- rbind(c(4.1,   2.8),
+           c(9.676, 6.608))
+kappa(A) # -> Inf
+kappa(A, exact=TRUE) # 8.675057e+15 ( 2-norm )
+
+## now for the 1-norm :
+try(kappa(A, exact=TRUE, norm = "1")) #-> Error: computationally singular
+    kappa(A, exact=TRUE, norm = "1", inv_z = solve(A, tol = 1e-19)) ## 5.22057e16
 }
 \keyword{math}
diff --git a/tests/reg-tests-1e.R b/tests/reg-tests-1e.R
@@ -1540,6 +1540,19 @@ stopifnot(identical(ch1, ch2),
 ## error msg was  "'length = 2' in coercion to 'logical(1)'"
 
 
+## kappa(*, exact=TRUE)  for exactly singular cases - PR#18817
+for(x in list(x3 = {n <- 3L; x <- diag(n); x[n,n] <- 0; x},
+              z2 = rbind(1:2, 0),
+              D0 = diag(0, nrow = 3))) { print(x)
+  stopifnot(exprs = {
+    identical(Inf,      kappa(x, exact = TRUE))
+    identical(Inf,      kappa(x, exact = TRUE, norm = "2"))
+    identical(Inf, .kappa_tri(x, exact = TRUE, norm = "2"))
+  })
+}
+## kappa(..)  returned 1 or {0 with a warning} in R <= 4.4.2
+
+
 
 ## keep at end
 rbind(last =  proc.time() - .pt,