Nothing
context("Cholesky solves")
set.seed(1, kind="Mersenne-Twister", normal.kind="Inversion")
test_that("Cholesky for ddiMatrix works", {
n <- 100L
M <- Diagonal(x=runif(n, 0.1, 1))
y0 <- rnorm(n)
y1 <- matrix(rnorm(3*n), n, 3)
y2 <- rsparsematrix(n, 10, 0.1)
ch <- build_chol(M)
M <- Diagonal(x = 0.1 + 10*rnorm(n)^2)
ch$update(M)
expect_equal(ch$solve(y0), y0/M@x)
expect_equal(ch$solve(y1), y1/M@x)
expect_equal(ch$solve(y2), Diagonal(x=1/M@x) %*% y2)
expect_equal(ch$solve(y0, system="Lt"), y0/sqrt(M@x))
expect_equal(ch$solve(y0, system="Lt", systemP=TRUE), y0/sqrt(M@x))
expect_equal(ch$solve(y0, system="L"), y0/sqrt(M@x))
expect_equal(ch$solve(y1, system="Lt"), y1/sqrt(M@x))
expect_equal(ch$solve(y2, system="Lt"), Diagonal(x=1/sqrt(M@x)) %*% y2)
})
test_that("Cholesky for matrix works", {
n <- 12
M <- crossprod(matrix(rnorm(n^2), n))
diag(M) <- diag(M) + runif(n, 0.1, 1)
ch <- build_chol(M)
M <- crossprod(matrix(rnorm(n^2), n))
diag(M) <- diag(M) + runif(n, 0.1, 1)
ch$update(M)
y0 <- rnorm(n)
y1 <- matrix(rnorm(3*n), n, 3)
expect_equal(ch$solve(y0), solve(M, y0))
expect_equal(ch$solve(y1), solve(M, y1))
expect_equal(ch$solve(y0, system="Lt"), solve(chol(M), y0))
expect_equal(ch$solve(y0, system="L"), solve(t(chol(M)), y0))
expect_equal(ch$solve(y1, system="Lt"), solve(chol(M), y1))
x <- rnorm(n)
expect_equal(ch$Ltimes(x), ch$cholM %m*v% x)
})
test_that("Cholesky for dsCMatrix works", {
n <- 100
M <- crossprod(rsparsematrix(n, n, density=0.01)) + Diagonal(n)
ch <- build_chol(M)
# in-place update only works for same sparsity pattern!!
M <- M + Diagonal(n)
ch$update(M)
cholM <- Cholesky_dsC(M, perm=FALSE, super=NA)
y0 <- rnorm(n)
y1 <- matrix(rnorm(3*n), n, 3)
y2 <- rsparsematrix(n, 10, 0.1)
expect_equal(ch$solve(y0), as.vector(solve(cholM, y0)))
expect_equal(ch$solve(y1), array(solve(cholM, y1)@x, dim(y1)))
expect_equal(ch$solve(y2), solve(cholM, y2))
expect_equal(ch$solve(y0, system="Lt"), as.vector(solve(cholM, y0, system="Lt")))
expect_equal(ch$solve(y1, system="Lt"), array(solve(cholM, y1, system="Lt")@x, dim(y1)))
expect_equal(ch$solve(y2, system="Lt"), solve(cholM, y2, system="Lt"))
expect_equal(ch$solve(y0, system="L"), as.vector(solve(cholM, y0, system="L")))
# NB disable following test; currently cannot control cholmod ordering
#expect_equal(ch$solve(y0, system="L", systemP=TRUE), as.vector(solve(cholM, y0, system="L")))
expect_equal(ch$solve(y1, system="L"), array(solve(cholM, y1, system="L")@x, dim(y1)))
expect_equal(ch$solve(y2, system="L"), solve(cholM, y2, system="L"))
# with permutation
ch <- build_chol(M, control=chol_control(perm=TRUE))
M <- M + Diagonal(x=runif(n))
ch$update(M)
cholM <- Cholesky_dsC(M, perm=TRUE, super=NA)
expect_equal(ch$solve(y0), as.vector(solve(cholM, y0)))
expect_equal(ch$solve(y1), array(solve(cholM, y1)@x, dim(y1)))
expect_equal(ch$solve(y2), solve(cholM, y2))
expect_equal(ch$solve(y0, system="Lt"), as.vector(solve(cholM, y0, system="Lt")))
expect_equal(ch$solve(y1, system="Lt"), array(solve(cholM, y1, system="Lt")@x, dim(y1)))
expect_equal(ch$solve(y2, system="Lt"), solve(cholM, y2, system="Lt"))
expect_equal(ch$solve(y0, system="L"), as.vector(solve(cholM, y0, system="L")))
expect_equal(ch$solve(y1, system="L"), array(solve(cholM, y1, system="L")@x, dim(y1)))
expect_equal(ch$solve(y2, system="L"), solve(cholM, y2, system="L"))
})
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.