The R Project SVN R-packages

R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.0 alpha (2005-03-31), ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for a HTML browser interface to help.
Type 'q()' to quit R.

> ### Testing positive definite matrices
> 
> library(Matrix)
> 
> 
> h9 <- Hilbert(9)
> stopifnot(c(0,0) == dim(Hilbert(0)),
+           c(9,9) == h9)
> str(h9)
Formal class 'dpoMatrix' [package "Matrix"] with 6 slots
  ..@ uplo    : chr "U"
  ..@ rcond   : num(0) 
  ..@ factors : list()
  ..@ x       : num [1:81] 1.000 0.500 0.333 0.250 0.200 ...
  ..@ Dim     : int [1:2] 9 9
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL
> all.equal(determinant(h9)$modulus, -96.7369450737858, tol= 1e-15)
[1] TRUE
> ##-> TRUE only on some platforms; seen relative difference of 10^-8
> stopifnot(0 == length(h9@factors))# nothing yet
> str(f9 <- as(chol(h9), "dtrMatrix"))
Formal class 'dtrMatrix' [package "Matrix"] with 7 slots
  ..@ uplo    : chr "U"
  ..@ diag    : chr "N"
  ..@ rcond   : num(0) 
  ..@ factors : list()
  ..@ x       : num [1:81] 1 0 0 0 0 0 0 0 0 0.5 ...
  ..@ Dim     : int [1:2] 9 9
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL
> ## h9 now has factorization
> stopifnot(names(h9@factors) == "Cholesky",
+           all.equal(rcond(h9), 9.0938e-13),
+           all.equal(rcond(f9), 9.1272e-7, tol = 1e-6))# more precision fails
> str(h9)# has 'rcond' and 'factors'
Formal class 'dpoMatrix' [package "Matrix"] with 6 slots
  ..@ uplo    : chr "U"
  ..@ rcond   : Named num 9.1e-13
  .. ..- attr(*, "names")= chr "O"
  ..@ factors :List of 1
  .. ..$ Cholesky:Formal class 'Cholesky' [package "Matrix"] with 7 slots
  .. .. .. ..@ uplo    : chr "U"
  .. .. .. ..@ diag    : chr "N"
  .. .. .. ..@ rcond   : num(0) 
  .. .. .. ..@ factors : list()
  .. .. .. ..@ x       : num [1:81] 1 0 0 0 0 0 0 0 0 0.5 ...
  .. .. .. ..@ Dim     : int [1:2] 9 9
  .. .. .. ..@ Dimnames:List of 2
  .. .. .. .. ..$ : NULL
  .. .. .. .. ..$ : NULL
  ..@ x       : num [1:81] 1.000 0.500 0.333 0.250 0.200 ...
  ..@ Dim     : int [1:2] 9 9
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL
> options(digits=4)
> (cf9 <- crossprod(f9))# looks the same as  h9 :
9 x 9 Matrix of class "dpoMatrix"
        [,1]   [,2]    [,3]    [,4]    [,5]    [,6]    [,7]    [,8]    [,9]
 [1,] 1.0000 0.5000 0.33333 0.25000 0.20000 0.16667 0.14286 0.12500 0.11111
 [2,] 0.5000 0.3333 0.25000 0.20000 0.16667 0.14286 0.12500 0.11111 0.10000
 [3,] 0.3333 0.2500 0.20000 0.16667 0.14286 0.12500 0.11111 0.10000 0.09091
 [4,] 0.2500 0.2000 0.16667 0.14286 0.12500 0.11111 0.10000 0.09091 0.08333
 [5,] 0.2000 0.1667 0.14286 0.12500 0.11111 0.10000 0.09091 0.08333 0.07692
 [6,] 0.1667 0.1429 0.12500 0.11111 0.10000 0.09091 0.08333 0.07692 0.07143
 [7,] 0.1429 0.1250 0.11111 0.10000 0.09091 0.08333 0.07692 0.07143 0.06667
 [8,] 0.1250 0.1111 0.10000 0.09091 0.08333 0.07692 0.07143 0.06667 0.06250
 [9,] 0.1111 0.1000 0.09091 0.08333 0.07692 0.07143 0.06667 0.06250 0.05882
> stopifnot(all.equal(as.matrix(h9),
+                     as.matrix(cf9), tol= 1e-15))
> 
> str(hp9 <- as(h9, "dppMatrix"))
Formal class 'dppMatrix' [package "Matrix"] with 6 slots
  ..@ uplo    : chr "U"
  ..@ rcond   : Named num 9.1e-13
  .. ..- attr(*, "names")= chr "O"
  ..@ factors :List of 1
  .. ..$ pCholesky:Formal class 'pCholesky' [package "Matrix"] with 7 slots
  .. .. .. ..@ uplo    : chr "U"
  .. .. .. ..@ diag    : chr "N"
  .. .. .. ..@ rcond   : num(0) 
  .. .. .. ..@ factors : list()
  .. .. .. ..@ x       : num [1:45] 1.000 0.500 0.289 0.333 0.289 ...
  .. .. .. ..@ Dim     : int [1:2] 9 9
  .. .. .. ..@ Dimnames:List of 2
  .. .. .. .. ..$ : NULL
  .. .. .. .. ..$ : NULL
  ..@ x       : num [1:45] 1.000 0.500 0.333 0.333 0.250 ...
  ..@ Dim     : int [1:2] 9 9
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL
> 
> s9 <- solve(hp9, seq(nrow(hp9)))
> signif(t(s9)/10000, 4)# only rounded numbers are platform-independent
1 x 9 Matrix of class "dgeMatrix"
       [,1]  [,2]  [,3]   [,4] [,5]  [,6] [,7]  [,8] [,9]
[1,] 0.0729 -5.76 109.5 -864.9 3468 -7668 9459 -6095 1597
> hp9 %*% s9
9 x 1 Matrix of class "dgeMatrix"
      [,1]
 [1,]    1
 [2,]    2
 [3,]    3
 [4,]    4
 [5,]    5
 [6,]    6
 [7,]    7
 [8,]    8
 [9,]    9
> ## Works in 2.1.0, but not earlier:
> ## stopifnot(all.equal(cbind(1:9), as.matrix(hp9 %*% s9)))
> stopifnot(all.equal(1:9, (hp9 %*% s9)@x))
>