R Under development (unstable) (2024-10-17 r87242 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 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 an HTML browser interface to help. Type 'q()' to quit R. > ## for R_DEFAULT_PACKAGES=NULL : > library(stats) > library(utils) > > library(Matrix) > > ### Matrix Products including cross products > > source(system.file("test-tools.R", package = "Matrix")) # is.EQ.mat(), dnIdentical() ..etc Loading required package: tools > doExtras [1] FALSE > options(warn=1, # show as they happen + Matrix.verbose = doExtras) > > ##' Check matrix multiplications with (unit) Diagonal matrices > chkDiagProd <- function(M) { + if(is.matrix(M)) { + noRN <- function(x) { + if(!is.null(dn <- dimnames(x))) + dimnames(x) <- c(list(NULL), dn[2L]) + x + } + noCN <- function(x) { + if(!is.null(dn <- dimnames(x))) + dimnames(x) <- c(dn[1L], list(NULL)) + x + } + } else if(is(M, "Matrix")) { + noRN <- function(x) { + x@Dimnames <- c(list(NULL), x@Dimnames[2L]) + x + } + noCN <- function(x) { + x@Dimnames <- c(x@Dimnames[1L], list(NULL)) + x + } + } else stop("'M' must be a [mM]atrix") + I.l <- Diagonal(nrow(M)) # I_n -- "unit" Diagonal + I.r <- Diagonal(ncol(M)) # I_d + D2.l <- Diagonal(nrow(M), x = 2) # D_n + D2.r <- Diagonal(ncol(M), x = 2) # I_d + stopifnot(is.EQ.mat(noCN(M), M %*% I.r), + is.EQ.mat(noRN(M), I.l %*% M), + is.EQ.mat(noCN(2*M), M %*% D2.r), + is.EQ.mat(noRN(M*2), D2.l %*% M), + ## crossprod + is.EQ.mat(noCN(t(M)), crossprod(M, I.l)), + is.EQ.mat(noRN(M), crossprod(I.l, M)), + is.EQ.mat(noCN(t(2*M)), crossprod(M, D2.l)), + is.EQ.mat(noRN(M*2) , crossprod(D2.l, M)), + ## tcrossprod + is.EQ.mat(noCN(M), tcrossprod(M, I.r)), + is.EQ.mat(noRN(t(M)), tcrossprod(I.r, M)), + is.EQ.mat(noCN(2*M), tcrossprod(M, D2.r)), + is.EQ.mat(noRN(t(M*2)), tcrossprod(D2.r, M))) + } > > ### dimnames -- notably for matrix products ---------------- > # > ##' Checks that matrix products are the same, including dimnames > ##' > ##' @param m matrix = traditional-R-matrix version of M > ##' @param M optional Matrix = "Matrix class version of m" > ##' @param browse > chkDnProd <- function(m = as(M, "matrix"), M = Matrix(m), browse=FALSE, warn.ok=FALSE) { + ## TODO: + ## if(browse) stopifnot <- f.unction(...) such that it enters browser() when it is not fulfilled + if(!warn.ok) { # NO warnings allowd + op <- options(warn = 2) + on.exit(options(op)) + } + stopifnot(is.matrix(m), is(M, "Matrix"), identical(dim(m), dim(M)), dnIdentical(m,M)) + ## m is n x d (say) + is.square <- nrow(m) == ncol(m) + + p1 <- (tm <- t(m)) %*% m ## d x d + p1. <- crossprod(m) + stopifnot(is.EQ.mat3(p1, p1., crossprod(m,m))) + + t1 <- m %*% tm ## n x n + t1. <- tcrossprod(m) + stopifnot(is.EQ.mat3(t1, t1., tcrossprod(m,m))) + + if(is.square) { + mm <- m %*% m + stopifnot(is.EQ.mat3(mm, crossprod(tm,m), tcrossprod(m,tm))) + } + + chkDiagProd(m)## was not ok in Matrix 1.2.0 + + ## Now the "Matrix" ones -- should match the "matrix" above + M0 <- M + cat("sparse: ") + for(sparse in c(TRUE, FALSE)) { + cat(sparse, "; ") + M <- as(M0, if(sparse) "sparseMatrix" else "denseMatrix") + P1 <- (tM <- t(M)) %*% M + P1. <- crossprod(M) + stopifnot(is.EQ.mat3(P1, P1., p1), + is.EQ.mat3(P1., crossprod(M,M), crossprod(M,m)), + is.EQ.mat (P1., crossprod(m,M))) + + ## P1. is "symmetricMatrix" -- semantically "must have" symm.dimnames + PP1 <- P1. %*% P1. ## still d x d + R <- triu(PP1); r <- as(R, "matrix") # upper - triangular + L <- tril(PP1); l <- as(L, "matrix") # lower - triangular + stopifnot(isSymmetric(P1.), isSymmetric(PP1), + isDiagonal(L) || is(L,"triangularMatrix"), + isDiagonal(R) || is(R,"triangularMatrix"), + isTriangular(L, upper=FALSE), isTriangular(R, upper=TRUE), + is.EQ.mat(PP1, (pp1 <- p1 %*% p1)), + dnIdentical(PP1, R), + dnIdentical(L, R)) + + T1 <- M %*% tM + T1. <- tcrossprod(M) + stopifnot(is.EQ.mat3(T1, T1., t1), + is.EQ.mat3(T1., tcrossprod(M,M), tcrossprod(M,m)), + is.EQ.mat (T1., tcrossprod(m,M)), + is.EQ.mat(tcrossprod(T1., tM), + tcrossprod(t1., tm)), + is.EQ.mat(crossprod(T1., M), + crossprod(t1., m))) + + ## Now, *mixing* Matrix x matrix: + stopifnot(is.EQ.mat3(tM %*% m, tm %*% M, tm %*% m)) + + if(is.square) + stopifnot(is.EQ.mat (mm, M %*% M), + is.EQ.mat3(mm, crossprod(tM,M), tcrossprod(M,tM)), + ## "mixing": + is.EQ.mat3(mm, crossprod(tm,M), tcrossprod(m,tM)), + is.EQ.mat3(mm, crossprod(tM,m), tcrossprod(M,tm))) + + ## Symmetric and Triangular + stopifnot(is.EQ.mat(PP1 %*% tM, pp1 %*% tm), + is.EQ.mat(R %*% tM, r %*% tm), + is.EQ.mat(L %*% tM, L %*% tm)) + + ## Diagonal : + chkDiagProd(M) + } + cat("\n") + invisible(TRUE) + } # chkDnProd() > > ## All these are ok {now, (2012-06-11) also for dense > (m <- matrix(c(0, 0, 2:0), 3, 5)) [,1] [,2] [,3] [,4] [,5] [1,] 0 1 0 0 2 [2,] 0 0 2 0 1 [3,] 2 0 1 0 0 > m00 <- m # *no* dimnames > dimnames(m) <- list(LETTERS[1:3], letters[1:5]) > (m.. <- m) # has *both* dimnames a b c d e A 0 1 0 0 2 B 0 0 2 0 1 C 2 0 1 0 0 > m0. <- m.0 <- m.. > dimnames(m0.)[1] <- list(NULL); m0. a b c d e [1,] 0 1 0 0 2 [2,] 0 0 2 0 1 [3,] 2 0 1 0 0 > dimnames(m.0)[2] <- list(NULL); m.0 [,1] [,2] [,3] [,4] [,5] A 0 1 0 0 2 B 0 0 2 0 1 C 2 0 1 0 0 > d <- diag(3); dimnames(d) <- list(c("u","v","w"), c("X","Y","Z")); d X Y Z u 1 0 0 v 0 1 0 w 0 0 1 > dU <- diagN2U(Matrix(d, doDiag = FALSE)) # unitriangular sparse > tU <- dU; tU[1,2:3] <- 3:4; tU[2,3] <- 7; tU # ditto "unitri" sparse 3 x 3 sparse Matrix of class "dtCMatrix" (unitriangular) X Y Z u I 3 4 v . I 7 w . . I > (T <- new("dtrMatrix", diag = "U", x= c(0,0,5,0), Dim= c(2L,2L), + Dimnames= list(paste0("r",1:2),paste0("C",1:2)))) # unitriangular dense 2 x 2 Matrix of class "dtrMatrix" (unitriangular) C1 C2 r1 1 5 r2 . 1 > pT <- pack(T)# ^^^^^^^^^^^^ > mt <- m[,2:3] %*% pT # deprecation warning in pre-1.5-4 > stopifnot(is(pT, "dtpMatrix"), validObject(pT), + validObject(mt), is(mt, "dgeMatrix"), + identical(as.matrix(mt), + array(c(1,0,0, 5,2,1), dim = 3:2, dimnames = list(c("A","B","C"), c("C1","C2")))) + ) > > A <- matrix(c(0.4, 0.1, 0, 0), 2) > B <- matrix(c(1.1, 0, 0, 0), 2); ABt <- tcrossprod(A, B) > ## Matrix version # both Matrix version: > class(AM <- as(A, "triangularMatrix")) # dtr [1] "dtrMatrix" attr(,"package") [1] "Matrix" > class(BM <- as(B, "Matrix")) # ddi [1] "ddiMatrix" attr(,"package") [1] "Matrix" > class(Bs <- as(as(B, "CsparseMatrix"), "symmetricMatrix")) # "dsC" [1] "dsCMatrix" attr(,"package") [1] "Matrix" > (ABst <- tcrossprod(A, Bs)) # was wrong since at least Matrix 1.3-3 (R-4.0.5) 2 x 2 Matrix of class "dgeMatrix" [,1] [,2] [1,] 0.44 0 [2,] 0.11 0 > assert.EQ.mat(ABst, ABt) > assert.EQ.mat(tcrossprod(AM, BM), ABt) > assert.EQ.mat(tcrossprod(AM, Bs), ABt) > assert.EQ.mat(tcrossprod(A , Bs), ABt) > > ## currently many warnings about sub-optimal matrix products : > chkDnProd(m..) sparse: TRUE ; FALSE ; > chkDnProd(m0.) sparse: TRUE ; FALSE ; > chkDnProd(m.0) sparse: TRUE ; FALSE ; > chkDnProd(m00) sparse: TRUE ; FALSE ; > chkDnProd(M = T) sparse: TRUE ; FALSE ; > chkDnProd(M = t(T)) sparse: TRUE ; FALSE ; > if(FALSE) { + ## FIXME: + ## Matrix() bug fix has revealed that diagonalMatrix product methods are not + ## handling (have never handled?) 'Dimnames' correctly, causing these to fail + chkDnProd(M = dU) + chkDnProd(M = t(dU)) + } > ## all the above failed in 1.2-0 and 1.1-5, 1.1-4 some even earlier > chkDnProd(M = tU) sparse: TRUE ; FALSE ; > chkDnProd(M = t(tU)) sparse: TRUE ; FALSE ; > ## the two above failed in Matrix <= 1.4-1 > chkDnProd(M = Diagonal(4)) sparse: TRUE ; FALSE ; > chkDnProd(diag(x=3:1)) sparse: TRUE ; FALSE ; > if(FALSE) { + ## FIXME: as for dU, t(dU) above + chkDnProd(d) + chkDnProd(M = as(d, "denseMatrix"))# M: dtrMatrix (diag = "N") + } > > > m5 <- 1 + as(diag(-1:4)[-5,], "generalMatrix") > ## named dimnames: > dimnames(m5) <- list(Rows= LETTERS[1:5], paste("C", 1:6, sep="")) > tr5 <- tril(m5[,-6]) > m. <- as(m5, "matrix") > m5.2 <- local({t5 <- as.matrix(tr5); t5 %*% t5}) > stopifnotValid(tr5, "dtrMatrix") > stopifnot(dim(m5) == 5:6, + class(cm5 <- crossprod(m5)) == "dpoMatrix") > assert.EQ.mat(t(m5) %*% m5, as(cm5, "matrix")) > assert.EQ.mat(tr5.2 <- tr5 %*% tr5, m5.2) > stopifnotValid(tr5.2, "dtrMatrix") > stopifnot(as.vector(rep(1,6) %*% cm5) == colSums(cm5), + as.vector(cm5 %*% rep(1,6)) == rowSums(cm5)) > ## uni-diagonal dtrMatrix with "wrong" entries in diagonal > ## {the diagonal can be anything: because of diag = "U" it should never be used}: > tru <- Diagonal(3, x=3); tru[i.lt <- lower.tri(tru, diag=FALSE)] <- c(2,-3,4) > tru@diag <- "U" ; stopifnot(diag(trm <- as.matrix(tru)) == 1) > ## TODO: Also add *upper-triangular* *packed* case !! > stopifnot((tru %*% tru)[i.lt] == + (trm %*% trm)[i.lt]) > > ## crossprod() with numeric vector RHS and LHS > i5 <- rep.int(1, 5) > isValid(S5 <- tcrossprod(tr5), "dpoMatrix")# and inherits from "dsy" [1] TRUE > G5 <- as(S5, "generalMatrix")# "dge" > assert.EQ.mat( crossprod(i5, m5), rbind( colSums(m5))) > assert.EQ.mat( crossprod(i5, m.), rbind( colSums(m5))) > assert.EQ.mat( crossprod(m5, i5), cbind( colSums(m5))) > assert.EQ.mat( crossprod(m., i5), cbind( colSums(m5))) > assert.EQ.mat( crossprod(i5, S5), rbind( colSums(S5))) # failed in Matrix 1.1.4 > ## tcrossprod() with numeric vector RHS and LHS : > stopifnot(identical(tcrossprod(i5, S5), # <- lost dimnames + tcrossprod(i5, G5) -> m51), + identical(dimnames(m51), list(NULL, Rows = LETTERS[1:5])) + ) > m51 <- m5[, 1, drop=FALSE] # [6 x 1] > m.1 <- m.[, 1, drop=FALSE] ; assert.EQ.mat(m51, m.1) > ## The only (M . v) case > assert.EQ.mat(tcrossprod(m51, 5:1), tcrossprod(m.1, 5:1)) > ## The two (v . M) cases: > assert.EQ.mat(tcrossprod(rep(1,6), m.), rbind( rowSums(m5)))# |v| = ncol(m) > assert.EQ.mat(tcrossprod(rep(1,3), m51), + tcrossprod(rep(1,3), m.1))# ncol(m) = 1 > > ## classes differ > tc.m5 <- m5 %*% t(m5) # "dge*", no dimnames (FIXME) > (tcm5 <- tcrossprod(m5)) # "dpo*" w/ dimnames 5 x 5 Matrix of class "dpoMatrix" Rows Rows A B C D E A 5 5 6 7 9 B 5 6 7 8 10 C 6 7 9 9 11 D 7 8 9 14 12 E 9 10 11 12 30 > assert.EQ.mat(tc.m5, mm5 <- as(tcm5, "matrix")) > ## tcrossprod(x,y) : > assert.EQ.mat(tcrossprod(m5, m5), mm5) > assert.EQ.mat(tcrossprod(m5, m.), mm5) > assert.EQ.mat(tcrossprod(m., m5), mm5) > > M50 <- m5[,FALSE, drop=FALSE] > M05 <- t(M50) > s05 <- as(M05, "sparseMatrix") > s50 <- t(s05) > assert.EQ.mat(M05, matrix(1, 0,5)) > assert.EQ.mat(M50, matrix(1, 5,0)) > assert.EQ.mat(tcrossprod(M50), tcrossprod(as(M50, "matrix"))) > assert.EQ.mat(tcrossprod(s50), tcrossprod(as(s50, "matrix"))) > assert.EQ.mat( crossprod(s50), crossprod(as(s50, "matrix"))) > stopifnot(identical( crossprod(s50), tcrossprod(s05)), + identical( crossprod(s05), tcrossprod(s50))) > (M00 <- crossprod(M50))## used to fail -> .Call(dgeMatrix_crossprod, x, FALSE) 0 x 0 Matrix of class "dpoMatrix" <0 x 0 matrix> > stopifnot(identical(M00, tcrossprod(M05)), + all(M00 == t(M50) %*% M50), dim(M00) == 0) > > ## simple cases with 'scalars' treated as 1x1 matrices: > d <- Matrix(1:5) > d %*% 2 5 x 1 Matrix of class "dgeMatrix" [,1] [1,] 2 [2,] 4 [3,] 6 [4,] 8 [5,] 10 > 10 %*% t(d) 1 x 5 Matrix of class "dgeMatrix" [,1] [,2] [,3] [,4] [,5] [1,] 10 20 30 40 50 > assertError(3 %*% d) # must give an error , similar to > assertError(5 %*% as.matrix(d)) # -> error > > ## right and left "numeric" and "matrix" multiplication: > (p1 <- m5 %*% c(10, 2:6)) 5 x 1 Matrix of class "dgeMatrix" Rows [,1] A 20 B 30 C 33 D 38 E 54 > (p2 <- c(10, 2:5) %*% m5) 1 x 6 Matrix of class "dgeMatrix" C1 C2 C3 C4 C5 C6 [1,] 14 24 27 32 24 44 > (pd1 <- m5 %*% diag(1:6)) 5 x 6 Matrix of class "dgeMatrix" Rows [,1] [,2] [,3] [,4] [,5] [,6] A 0 2 3 4 5 6 B 1 2 3 4 5 6 C 1 2 6 4 5 6 D 1 2 3 12 5 6 E 1 2 3 4 5 30 > (pd. <- m5 %*% Diagonal(x = 1:6)) 5 x 6 Matrix of class "dgeMatrix" Rows [,1] [,2] [,3] [,4] [,5] [,6] A 0 2 3 4 5 6 B 1 2 3 4 5 6 C 1 2 6 4 5 6 D 1 2 3 12 5 6 E 1 2 3 4 5 30 > (pd2 <- diag (10:6) %*% m5) 5 x 6 Matrix of class "dgeMatrix" C1 C2 C3 C4 C5 C6 [1,] 0 10 10 10 10 10 [2,] 9 9 9 9 9 9 [3,] 8 8 16 8 8 8 [4,] 7 7 7 21 7 7 [5,] 6 6 6 6 6 30 > (pd..<- Diagonal(x = 10:6) %*% m5) 5 x 6 Matrix of class "dgeMatrix" C1 C2 C3 C4 C5 C6 [1,] 0 10 10 10 10 10 [2,] 9 9 9 9 9 9 [3,] 8 8 16 8 8 8 [4,] 7 7 7 21 7 7 [5,] 6 6 6 6 6 30 > stopifnot(exprs = { + dim(crossprod(t(m5))) == c(5,5) + c(class(p1),class(p2),class(pd1),class(pd2), class(pd.),class(pd..)) == "dgeMatrix" + identical(dimnames(pd.), c(dimnames(m5)[1L], list(NULL))) + identical(dimnames(pd..), c(list(NULL), dimnames(m5)[2L])) + }) > > assert.EQ.mat(p1, cbind(c(20,30,33,38,54))) > assert.EQ.mat(pd1, m. %*% diag(1:6)) > assert.EQ.mat(pd2, diag(10:6) %*% m.) > assert.EQ.mat(pd., as(pd1,"matrix")) > assert.EQ.mat(pd..,as(pd2,"matrix")) > > ## check that 'solve' and '%*%' are inverses > suppressWarnings(RNGversion("3.5.0")); set.seed(1) > A <- Matrix(rnorm(25), ncol = 5) > y <- rnorm(5) > all.equal((A %*% solve(A, y))@x, y) [1] TRUE > Atr <- new("dtrMatrix", Dim = A@Dim, x = A@x, uplo = "U") > all.equal((Atr %*% solve(Atr, y))@x, y) [1] TRUE > > ## R-forge bug 5933 by Sebastian Herbert, > ## https://r-forge.r-project.org/tracker/index.php?func=detail&aid=5933&group_id=61&atid=294 > mLeft <- matrix(data = double(0), nrow = 3, ncol = 0) > mRight <- matrix(data = double(0), nrow = 0, ncol = 4) > MLeft <- Matrix(data = double(0), nrow = 3, ncol = 0) > MRight <- Matrix(data = double(0), nrow = 0, ncol = 4) > stopifnot(exprs = { + class(mLeft) == class(mRight) + class(MLeft) == class(MRight) + class(MLeft) == "dgeMatrix" + }) > > Qidentical3 <- function(a,b,c) Q.eq(a,b) && Q.eq(b,c) > Qidentical4 <- function(a,b,c,d) Q.eq(a,b) && Q.eq(b,c) && Q.eq(c,d) > chkP <- function(mLeft, mRight, MLeft, MRight, cl = class(MLeft)) { + ident4 <- + if(extends(cl, "generalMatrix")) { + function(a,b,c,d) { + r <- eval(Matrix:::.as.via.virtual( + class(a)[1L], cl[1], quote(a))) + identical4(r,b,c,d) + } + } else { + function(a,b,c,d) { + assert.EQ.mat(M=b, m=a, tol=0) + Qidentical3(b,c,d) + } + } + mm <- mLeft %*% mRight # ok + m.m <- crossprod(mRight) + mm. <- tcrossprod(mLeft, mLeft) + stopifnot(mm == 0, + ident4(mm, + mLeft %*% MRight, + MLeft %*% mRight, + MLeft %*% MRight),# now ok + m.m == 0, identical(m.m, crossprod(mRight, mRight)), + mm. == 0, identical(mm., tcrossprod(mLeft, mLeft)), allow.logical0 = TRUE) + stopifnot(ident4(m.m, + crossprod(MRight, MRight), + crossprod(MRight, mRight), + crossprod(mRight, MRight))) + stopifnot(ident4(mm., + tcrossprod(mLeft, MLeft), + tcrossprod(MLeft, MLeft), + tcrossprod(MLeft, mLeft))) + } > > chkP(mLeft, mRight, MLeft, MRight, "dgeMatrix") > m0 <- mLeft[FALSE,] # 0 x 0 > for(cls in c("triangularMatrix", "symmetricMatrix")) { + cat(cls,": "); stopifnotValid(ML0 <- as(MLeft[FALSE,], cls), cls) + chkP(m0, mRight, ML0, MRight, class(ML0)) + chkP(mLeft, m0 , MLeft, ML0 , class(ML0)) + chkP(m0, m0 , ML0, ML0 , class(ML0)); cat("\n") + } triangularMatrix : symmetricMatrix : > > > ## New in R 3.2.0 -- for traditional matrix m and vector v > for(spV in c(FALSE,TRUE)) { + cat("sparseV:", spV, "\n") + v <- if(spV) as(1:3, "sparseVector") else 1:3 + stopifnot(identical(class(v2 <- v[1:2]), class(v))) + assertError(crossprod(v, v2)) ; assertError(v %*% v2) + assertError(crossprod(v, 1:2)); assertError(v %*% 1:2) + assertError(crossprod(v, 2)) ; assertError(v %*% 2) + assertError(crossprod(1:2, v)); assertError(1:2 %*% v) + cat("doing vec x vec ..\n") + stopifnot(identical(crossprod(2, v), t(2) %*% v), + identical(5 %*% v, 5 %*% t(v))) + for(sp in c(FALSE, TRUE)) { + m <- Matrix(1:2, 1,2, sparse=sp) + cat(sprintf("class(m): '%s'\n", class(m))) + stopifnot(identical( crossprod(m, v), t(m) %*% v), # m'v gave *outer* prod wrongly! + identical(tcrossprod(m, v2), m %*% v2)) + assert.EQ.Mat(m %*% v2, m %*% 1:2, tol=0) + } + ## gave error "non-conformable arguments" + } sparseV: FALSE doing vec x vec .. class(m): 'dgeMatrix' class(m): 'dgCMatrix' sparseV: TRUE doing vec x vec .. class(m): 'dgeMatrix' class(m): 'dgCMatrix' > ## crossprod(m, v) t(1 x 2) * 3 ==> (2 x 1) * (1 x 3) ==> 2 x 3 > ## tcrossprod(m,v2) 1 x 2 * 2 ==> (1 x 2) * t(1 x 2) ==> 1 x 1 > > ### ------ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - > ### Sparse Matrix products > ### ------ > ## solve() for dtC* > mc <- round(chol(crossprod(A)), 2) > B <- A[1:3,] # non-square on purpose > stopifnot(all.equal(sum(rowSums(B %*% mc)), 5.82424475145)) > assert.EQ.mat(tcrossprod(B, mc), as.matrix(t(tcrossprod(mc, B)))) > > m <- kronecker(Diagonal(2), mc) > stopifnot(is(mc, "dtrMatrix"), + is(m, "sparseMatrix")) > im <- solve(m) > round(im, 3) 10 x 10 sparse Matrix of class "dtCMatrix" [1,] 0.515 -0.161 0.412 -0.095 2.126 . . . . . [2,] . 0.763 0.896 -0.473 -9.344 . . . . . [3,] . . 0.645 -0.028 -10.444 . . . . . [4,] . . . 0.862 1.533 . . . . . [5,] . . . . 11.111 . . . . . [6,] . . . . . 0.515 -0.161 0.412 -0.095 2.126 [7,] . . . . . . 0.763 0.896 -0.473 -9.344 [8,] . . . . . . . 0.645 -0.028 -10.444 [9,] . . . . . . . . 0.862 1.533 [10,] . . . . . . . . . 11.111 > itm <- solve(t(m)) > iim <- solve(im) # should be ~= 'm' of course > iitm <- solve(itm) > I <- Diagonal(nrow(m)) > (del <- c(mean(abs(as.numeric(im %*% m - I))), + mean(abs(as.numeric(m %*% im - I))), + mean(abs(as.numeric(im - t(itm)))), + mean(abs(as.numeric( m - iim))), + mean(abs(as.numeric(t(m)- iitm))))) [1] 2.081668e-17 2.220446e-18 1.168510e-16 2.942091e-17 2.442491e-17 > stopifnot(is(m, "triangularMatrix"), is(m, "sparseMatrix"), + is(im, "dtCMatrix"), is(itm, "dtCMatrix"), is(iitm, "dtCMatrix"), + del < 1e-15) > > ## crossprod(.,.) & tcrossprod(), mixing dense & sparse > v <- c(0,0,2:0) > mv <- as.matrix(v) ## == cbind(v) > (V <- Matrix(v, 5,1, sparse=TRUE)) 5 x 1 sparse Matrix of class "dgCMatrix" [1,] . [2,] . [3,] 2 [4,] 1 [5,] . > sv <- as(v, "sparseVector") > a <- as.matrix(A) > cav <- crossprod(a,v) > tva <- tcrossprod(v,a) > assert.EQ.mat(crossprod(A, V), cav) # gave infinite recursion > assert.EQ.mat(crossprod(A,sv), cav) > assert.EQ.mat(tcrossprod( sv, A), tva) > assert.EQ.mat(tcrossprod(t(V),A), tva) > > ## [t]crossprod() for . incl. one arg.: > stopifnotValid(s.s <- crossprod(sv,sv), "Matrix") > stopifnotValid(ss. <- tcrossprod(sv,sv), "sparseMatrix") > stopifnot(identical(as(s.s, "symmetricMatrix"), crossprod(sv)), + identical(as(ss., "symmetricMatrix"), tcrossprod(sv))) > assert.EQ.mat(s.s, crossprod(v,v)) > assert.EQ.mat(ss., tcrossprod(v,v)) > > dm <- Matrix(v, sparse=FALSE) > sm <- Matrix(v, sparse=TRUE) > > vvt <- tcrossprod(v, v) > validObject(d.vvt <- as(as(vvt, "unpackedMatrix"), "generalMatrix")) [1] TRUE > validObject(s.vvt <- as(as(vvt, "CsparseMatrix"), "generalMatrix")) [1] TRUE > > stopifnot( + identical4(tcrossprod(v, v), tcrossprod(mv, v), + tcrossprod(v,mv), tcrossprod(mv,mv))## (base R) + , + identical4(d.vvt, tcrossprod(dm, v), + tcrossprod(v,dm), tcrossprod(dm,dm)) ## (v, dm) failed + , + identical(s.vvt, tcrossprod(sm,sm)) + , + identical3(d.vvt, tcrossprod(sm, v), + tcrossprod(v,sm)) ## both (sm,v) and (v,sm) failed + ) > assert.EQ.mat(d.vvt, vvt) > assert.EQ.mat(s.vvt, vvt) > > M <- Matrix(0:5, 2,3) ; sM <- as(M, "sparseMatrix"); m <- as(M, "matrix") > v <- 1:3; v2 <- 2:1 > sv <- as( v, "sparseVector") > sv2 <- as(v2, "sparseVector") > tvm <- tcrossprod(v, m) > assert.EQ.mat(tcrossprod( v, M), tvm) > assert.EQ.mat(tcrossprod( v,sM), tvm) > assert.EQ.mat(tcrossprod(sv,sM), tvm) > assert.EQ.mat(tcrossprod(sv, M), tvm) > assert.EQ.mat(crossprod(M, sv2), crossprod(m, v2)) > stopifnot(identical(tcrossprod(v, M), v %*% t(M)), + identical(tcrossprod(v,sM), v %*% t(sM)), + identical(tcrossprod(v, M), crossprod(v, t(M))), + identical(tcrossprod(sv,sM), sv %*% t(sM)), + identical(crossprod(sM, sv2), t(sM) %*% sv2), + identical(crossprod(M, v2), t(M) %*% v2)) > > ## *unit* triangular : > t1 <- new("dtTMatrix", x= c(3,7), i= 0:1, j=3:2, Dim= as.integer(c(4,4))) > ## from 0-diagonal to unit-diagonal {low-level step}: > tu <- t1 ; tu@diag <- "U" > cu <- as(tu, "CsparseMatrix") > cl <- t(cu) # unit lower-triangular > cl10 <- cl %*% Diagonal(4, x=10) > assert.EQ.mat(cl10, as(cl, "matrix") %*% diag(4, x=10)) > stopifnot(is(cl,"dtCMatrix"), cl@diag == "U") > (cu2 <- cu %*% cu) 4 x 4 sparse Matrix of class "dtCMatrix" (unitriangular) [1,] I . . 6 [2,] . I 14 . [3,] . . I . [4,] . . . I > cl2 <- cl %*% cl > validObject(cl2) [1] TRUE > > cu3 <- tu[-1,-1] > assert.EQ.mat(crossprod(tru, cu3),## <- "FIXME" should return triangular ... + crossprod(trm, as.matrix(cu3))) > > cl2 4 x 4 sparse Matrix of class "dtCMatrix" (unitriangular) [1,] I . . . [2,] . I . . [3,] . 14 I . [4,] 6 . . I > mcu <- as.matrix(cu) > cu2. <- Diagonal(4) + Matrix(c(rep(0,9),14,0,0,6,0,0,0), 4,4) > D4 <- Diagonal(4, x=10:7); d4 <- as(D4, "matrix") > D.D4 <- crossprod(D4); assert.EQ.mat(D.D4, crossprod(d4)) > stopifnotValid(D.D4, "ddiMatrix") > stopifnotValid(su <- crossprod(cu), "dsCMatrix") > asGe <- function(x) as(as(x, "unpackedMatrix"), "generalMatrix") > stopifnot(exprs = { + all(cu2 == cu2.) + ## was wrong for ver. <= 0.999375-4 + identical(D.D4, tcrossprod(D4)) + identical4(crossprod(d4, D4), crossprod(D4, d4), tcrossprod(d4, D4), asGe(D.D4)) + is(cu2, "dtCMatrix") # triangularity preserved + is(cl2, "dtCMatrix") + cu2@diag == "U" # unit triangularity preserved + cl2@diag == "U" + all.equal(D4 %*% mcu, asGe(D4 %*% cu)) + all.equal(mcu %*% D4, asGe(cu %*% D4)) + all(D4 %*% su == D4 %*% as.mat(su)) + all(su %*% D4 == as.mat(su) %*% D4) + identical(t(cl2), cu2) + ## !! + identical(crossprod(mcu, D4), asGe(crossprod(cu, D4))) + identical4(asGe(tcrossprod(cu, D4)), tcrossprod(mcu, D4), + asGe(cu %*% D4), mcu %*% D4) + identical4(asGe(tcrossprod(D4, cu)), tcrossprod(D4,mcu), + asGe(D4 %*% t(cu)), D4 %*% t(mcu)) + identical( crossprod(cu), Matrix( crossprod(mcu),sparse=TRUE)) + identical(tcrossprod(cu), Matrix(tcrossprod(mcu),sparse=TRUE)) + }) > assert.EQ.mat( crossprod(cu, D4), crossprod(mcu, d4)) > assert.EQ.mat(tcrossprod(cu, D4), tcrossprod(mcu, d4)) > tr8 <- kronecker(Matrix(rbind(c(2,0),c(1,4))), cl2) > T8 <- tr8 %*% (tr8/2) # triangularity preserved? > T8.2 <- (T8 %*% T8) / 4 > stopifnot(is(T8, "triangularMatrix"), T8@uplo == "L", is(T8.2, "dtCMatrix")) > mr8 <- as(tr8,"matrix") > m8. <- (mr8 %*% mr8 %*% mr8 %*% mr8)/16 > assert.EQ.mat(T8.2, m8.) > > data(KNex, package = "Matrix") > mm <- KNex$mm > M <- mm[1:500, 1:200] > MT <- as(M, "TsparseMatrix") > cpr <- t(mm) %*% mm > cpr. <- crossprod(mm) > cpr.. <- crossprod(mm, mm) > stopifnot(is(cpr., "symmetricMatrix"), + identical3(cpr, as(cpr., "generalMatrix"), cpr..)) > ## with dimnames: > m <- Matrix(c(0, 0, 2:0), 3, 5) > dimnames(m) <- list(LETTERS[1:3], letters[1:5]) > m 3 x 5 sparse Matrix of class "dgCMatrix" a b c d e A . 1 . . 2 B . . 2 . 1 C 2 . 1 . . > p1 <- t(m) %*% m > (p1. <- crossprod(m)) 5 x 5 sparse Matrix of class "dsCMatrix" a b c d e a 4 . 2 . . b . 1 . . 2 c 2 . 5 . 2 d . . . . . e . 2 2 . 5 > t1 <- m %*% t(m) > (t1. <- tcrossprod(m)) 3 x 3 sparse Matrix of class "dsCMatrix" A B C A 5 2 . B 2 5 2 C . 2 5 > stopifnot(isSymmetric(p1.), + isSymmetric(t1.), + identical(p1, as(p1., "generalMatrix")), + identical(t1, as(t1., "generalMatrix")), + identical(dimnames(p1), dimnames(p1.)), + identical(dimnames(p1), list(colnames(m), colnames(m))), + identical(dimnames(t1), dimnames(t1.)) + ) > > showMethods("%*%", classes=class(M)) Function: %*% (package base) x="ddiMatrix", y="dgCMatrix" (inherited from: x="diagonalMatrix", y="CsparseMatrix") x="dgCMatrix", y="ddiMatrix" (inherited from: x="CsparseMatrix", y="diagonalMatrix") x="dgCMatrix", y="dgCMatrix" (inherited from: x="CsparseMatrix", y="CsparseMatrix") x="dgCMatrix", y="dgRMatrix" (inherited from: x="CsparseMatrix", y="RsparseMatrix") x="dgCMatrix", y="dtCMatrix" (inherited from: x="CsparseMatrix", y="CsparseMatrix") x="dgCMatrix", y="integer" (inherited from: x="CsparseMatrix", y="vector") x="dgCMatrix", y="isparseVector" (inherited from: x="Matrix", y="sparseVector") x="dgCMatrix", y="matrix" (inherited from: x="CsparseMatrix", y="matrix") x="dgRMatrix", y="dgCMatrix" (inherited from: x="RsparseMatrix", y="CsparseMatrix") x="dgeMatrix", y="dgCMatrix" (inherited from: x="denseMatrix", y="CsparseMatrix") x="dtCMatrix", y="dgCMatrix" (inherited from: x="CsparseMatrix", y="CsparseMatrix") x="integer", y="dgCMatrix" (inherited from: x="vector", y="CsparseMatrix") x="isparseVector", y="dgCMatrix" (inherited from: x="sparseVector", y="Matrix") x="matrix", y="dgCMatrix" (inherited from: x="matrix", y="CsparseMatrix") > > v1 <- rep(1, ncol(M)) > str(r <- M %*% Matrix(v1)) Formal class 'dgeMatrix' [package "Matrix"] with 4 slots ..@ Dim : int [1:2] 500 1 ..@ Dimnames:List of 2 .. ..$ : NULL .. ..$ : NULL ..@ x : num [1:500] 0.277 0.5 0.277 0.5 0.5 ... ..@ factors : list() > str(rT <- MT %*% Matrix(v1)) Formal class 'dgeMatrix' [package "Matrix"] with 4 slots ..@ Dim : int [1:2] 500 1 ..@ Dimnames:List of 2 .. ..$ : NULL .. ..$ : NULL ..@ x : num [1:500] 0.277 0.5 0.277 0.5 0.5 ... ..@ factors : list() > stopifnot(identical(r, rT)) > str(r. <- M %*% as.matrix(v1)) Formal class 'dgeMatrix' [package "Matrix"] with 4 slots ..@ Dim : int [1:2] 500 1 ..@ Dimnames:List of 2 .. ..$ : NULL .. ..$ : NULL ..@ x : num [1:500] 0.277 0.5 0.277 0.5 0.5 ... ..@ factors : list() > stopifnot(identical4(r, r., rT, M %*% as(v1, "matrix"))) > > v2 <- rep(1,nrow(M)) > r2 <- t(Matrix(v2)) %*% M > r2T <- v2 %*% MT > str(r2. <- v2 %*% M) Formal class 'dgeMatrix' [package "Matrix"] with 4 slots ..@ Dim : int [1:2] 1 200 ..@ Dimnames:List of 2 .. ..$ : NULL .. ..$ : NULL ..@ x : num [1:200] 1.66 2 1.33 1.15 1.79 ... ..@ factors : list() > stopifnot(identical3(r2, r2., t(as(v2, "matrix")) %*% M)) > > > ###------------------------------------------------------------------ > ### Close to singular matrix W > ### (from micEconAids/tests/aids.R ... priceIndex = "S" ) > (load(system.file("external", "symW.rda", package="Matrix"))) # "symW" [1] "symW" > stopifnot(is(symW, "symmetricMatrix")) > n <- nrow(symW) > I <- .sparseDiagonal(n, shape="g") > S <- as(symW, "matrix") > sis <- solve(S, S) > ## solve(, ) when Cholmod fails was buggy for *long*: > o. <- options(Matrix.verbose = 2) # <-- showing Cholmod error & warning now > SIS <- solve(symW, symW) Warning in .local(A, ...) : CHOLMOD warning 'not positive definite' at file 'Cholesky/t_cholmod_rowfac_worker.c', line 449 > iw <- solve(symW) ## << TODO: LU *not* saved in @factors Warning in .local(A, ...) : CHOLMOD warning 'not positive definite' at file 'Cholesky/t_cholmod_rowfac_worker.c', line 449 > iss <- iw %*% symW > ## nb-mm3 openBLAS (Avi A.) > assert.EQ.(I, drop0(sis), tol = 1e-8)# 2.6e-10; 7.96e-9 Mean relative difference: 2.602095e-10 > assert.EQ.(I, SIS, tol = 1e-7)# 8.2e-9 Mean relative difference: 3.044365e-09 > assert.EQ.(I, iss, tol = 4e-4)# 3.3e-5 Mean relative difference: 3.255303e-05 > ## solve(, ) : > I <- diag(nrow=n) > SIS <- solve(symW, as(symW,"denseMatrix")) Warning in .local(A, ...) : CHOLMOD warning 'not positive definite' at file 'Cholesky/t_cholmod_rowfac_worker.c', line 449 > iw <- solve(symW, I) Warning in .local(A, ...) : CHOLMOD warning 'not positive definite' at file 'Cholesky/t_cholmod_rowfac_worker.c', line 449 > iss <- iw %*% symW > assert.EQ.mat(SIS, I, tol = 1e-7, giveRE=TRUE) Mean relative difference: 4.079679e-09 > assert.EQ.mat(iss, I, tol = 4e-4, giveRE=TRUE) Mean relative difference: 3.398617e-05 > rm(SIS,iss) > > WW <- as(symW, "generalMatrix") # the one that gave problems > IW <- solve(WW) > class(I1 <- IW %*% WW)# "dge" or "dgC" (!) [1] "dgCMatrix" attr(,"package") [1] "Matrix" > class(I2 <- WW %*% IW) [1] "dgCMatrix" attr(,"package") [1] "Matrix" > ## these two were wrong for for M.._1.0-13: > assert.EQ.(as(I1,"matrix"), I, tol = 1e-4) Mean relative difference: 3.255197e-05 > assert.EQ.(as(I2,"matrix"), I, tol = 7e-7) Mean relative difference: 5.434198e-09 > > ## now slightly perturb WW (and hence break exact symmetry > set.seed(131); ii <- sample(length(WW), size= 100) > WW[ii] <- WW[ii] * (1 + 1e-7*runif(100)) replCmat[x,i,j,..,val] : nargs()=3; missing (i,j) = (0,1) diagnosing replTmat(x,i,j,v): nargs()= 3; missing (i,j) = (0,1) > SW. <- symmpart(WW) > SW2 <- forceSymmetric(WW) > stopifnot(all.equal(as(SW.,"matrix"), + as(SW2,"matrix"), tolerance = 1e-7)) > (ch <- all.equal(WW, as(SW., "generalMatrix"), tolerance = 0)) [1] "Mean relative difference: 1.746244e-08" > stopifnot(is.character(ch), length(ch) == 1)## had length(.) 2 previously > IW <- solve(WW) # ( => stores in WW@factors !) > class(I1 <- IW %*% WW)# "dge" or "dgC" (!) [1] "dgCMatrix" attr(,"package") [1] "Matrix" > class(I2 <- WW %*% IW) [1] "dgCMatrix" attr(,"package") [1] "Matrix" > I <- diag(nrow=nrow(WW)) > stopifnot(all.equal(as(I1,"matrix"), I, check.attributes=FALSE, tolerance = 1e-4), + ## "Mean relative difference: 3.296549e-05" (or "1.999949" for Matrix_1.0-13 !!!) + all.equal(as(I2,"matrix"), I, check.attributes=FALSE)) #default tol gives "1" for M.._1.0-13 > options(o.) # revert to less Matrix.verbose > > if(doExtras) { + print(kappa(WW)) ## [1] 5.129463e+12 + print(rcond(WW)) ## [1] 6.216103e-14 + ## Warning message: rcond(.) via sparse -> dense coercion + } > > class(Iw. <- solve(SW.))# FIXME? should be "symmetric" but is not Warning in .local(A, ...) : CHOLMOD warning 'not positive definite' at file 'Cholesky/t_cholmod_rowfac_worker.c', line 449 [1] "dgCMatrix" attr(,"package") [1] "Matrix" > class(Iw2 <- solve(SW2))# FIXME? should be "symmetric" but is not Warning in .local(A, ...) : CHOLMOD warning 'not positive definite' at file 'Cholesky/t_cholmod_rowfac_worker.c', line 449 [1] "dgCMatrix" attr(,"package") [1] "Matrix" > class(IW. <- as(Iw., "denseMatrix")) [1] "dgeMatrix" attr(,"package") [1] "Matrix" > class(IW2 <- as(Iw2, "denseMatrix")) [1] "dgeMatrix" attr(,"package") [1] "Matrix" > > ### The next two were wrong for very long, too > assert.EQ.(I, as.matrix(IW. %*% SW.), tol= 4e-4) Mean relative difference: 8.17382e-05 > assert.EQ.(I, as.matrix(IW2 %*% SW2), tol= 4e-4) Mean relative difference: 8.776952e-05 > dIW <- as(IW, "denseMatrix") > assert.EQ.(dIW, IW., tol= 4e-4) Mean relative difference: 5.254132e-05 > assert.EQ.(dIW, IW2, tol= 8e-4) Mean relative difference: 0.0002166145 > > ##------------------------------------------------------------------ > > > > ## Sparse Cov.matrices from Harri Kiiveri @ CSIRO > a <- matrix(0,5,5) > a[1,2] <- a[2,3] <- a[3,4] <- a[4,5] <- 1 > a <- a + t(a) + 2*diag(5) > b <- as(a, "CsparseMatrix") ## ok, but we recommend to use Matrix() ``almost always'' : > (b. <- Matrix(a, sparse = TRUE)) 5 x 5 sparse Matrix of class "dsCMatrix" [1,] 2 1 . . . [2,] 1 2 1 . . [3,] . 1 2 1 . [4,] . . 1 2 1 [5,] . . . 1 2 > stopifnot(identical(b, b.)) > > ## calculate conditional variance matrix ( vars 3 4 5 given 1 2 ) > (B2 <- b[1:2, 1:2]) 2 x 2 sparse Matrix of class "dsCMatrix" [1,] 2 1 [2,] 1 2 > bb <- b[1:2, 3:5] > stopifnot(is(B2, "dsCMatrix"), # symmetric indexing keeps symmetry + identical(as.mat(bb), rbind(0, c(1,0,0))), + ## TODO: use fully-sparse cholmod_spsolve() based solution : + is(z.s <- solve(B2, bb), "sparseMatrix")) > assert.EQ.mat(B2 %*% z.s, as(bb, "matrix")) > ## -> dense RHS and dense result > z. <- solve(as(B2, "generalMatrix"), bb)# now *sparse* > z <- solve( B2, as(bb, "denseMatrix")) > stopifnot(is(z., "sparseMatrix"), + all.equal(z, as(z.,"denseMatrix"))) > ## finish calculating conditional variance matrix > v <- b[3:5,3:5] - crossprod(bb,z) > stopifnot(all.equal(as.mat(v), + matrix(c(4/3, 1:0, 1,2,1, 0:2), 3), tolerance = 1e-14)) > > > ###--- "logical" Matrices : --------------------- > > ##__ FIXME __ now works for lsparse* and nsparse* but not yet for lge* and nge* ! > > ## Robert's Example, a bit more readable > fromTo <- rbind(c(2,10), + c(3, 9)) > N <- 10 > nrFT <- nrow(fromTo) > rowi <- rep.int(1:nrFT, fromTo[,2]-fromTo[,1] + 1) - 1:1 > coli <- unlist(lapply(1:nrFT, function(x) fromTo[x,1]:fromTo[x,2])) - 1:1 > > # ## "n" --- nonzero pattern Matrices > > chk.ngMatrix <- function(M, verbose = TRUE) { + if(!(is(M, "nsparseMatrix") && length(d <- dim(M)) == 2 && d[1] == d[2])) + stop("'M' must be a square sparse [patter]n Matrix") + if(verbose) + show(M) + m <- as(M, "matrix") + + ## Part I : matrix products of pattern Matrices + ## ------ For now [by default]: *pattern* <==> boolean arithmetic + ## ==> FIXME ??: warning that this will change? + MM <- M %*% M # numeric (dgC) + if(verbose) { cat("M %*% M:\n"); show(MM) } + assert.EQ.mat(MM, m %*% m) + assert.EQ.mat(t(M) %&% M, + (t(m) %*% m) > 0, tol=0) + cM <- crossprod(M) # pattern {FIXME ?? warning ...} + tM <- tcrossprod(M) # pattern {FIXME ?? warning ...} + if(verbose) {cat( "crossprod(M):\n"); show(cM) } + if(verbose) {cat("tcrossprod(M):\n"); show(tM) } + stopifnot(is(cM,"symmetricMatrix"), is(tM,"symmetricMatrix"), + identical(as(as(cM, "CsparseMatrix"), "generalMatrix"), + t(M) %&% M), + identical(as(as(tM, "CsparseMatrix"), "generalMatrix"), + M %&% t(M))) + assert.EQ.mat(cM, crossprod(m) > 0) + assert.EQ.mat(tM, as(tcrossprod(m), "matrix") > 0) + + ## Part II : matrix products pattern Matrices with numeric: + ## + ## "n" x "d" (and "d" x "n") --> "d", i.e. numeric in any case + dM <- as(M, "dMatrix") + stopifnot(exprs = { + ## dense ones: + identical(M %*% m, m %*% M -> Mm) + ## sparse ones : + identical3(M %*% dM, dM %*% M -> sMM, + eval(Matrix:::.as.via.virtual( + "matrix", class(sMM), quote(m %*% m)))) + }) + if(verbose) {cat( "M %*% m:\n"); show(Mm) } + stopifnotValid(Mm, "dMatrix") # not "n.." + stopifnotValid(sMM, "dMatrix") # not "n.." + stopifnotValid(cdM <- crossprod(dM, M), "CsparseMatrix") + stopifnotValid(tdM <- tcrossprod(dM, M), "CsparseMatrix") + assert.EQ.mat (cdM, crossprod(m)) + assert.EQ.mat (tdM, tcrossprod(m)) + stopifnot(identical( crossprod(dM), as(cdM, "symmetricMatrix"))) + stopifnot(identical(tcrossprod(dM), as(tdM, "symmetricMatrix"))) + invisible(TRUE) + } > > sM <- new("ngTMatrix", i = rowi, j=coli, Dim=as.integer(c(N,N))) > chk.ngMatrix(sM) # "ngTMatrix" 10 x 10 sparse Matrix of class "ngTMatrix" [1,] . | | | | | | | | | [2,] . . | | | | | | | . [3,] . . . . . . . . . . [4,] . . . . . . . . . . [5,] . . . . . . . . . . [6,] . . . . . . . . . . [7,] . . . . . . . . . . [8,] . . . . . . . . . . [9,] . . . . . . . . . . [10,] . . . . . . . . . . M %*% M: 10 x 10 sparse Matrix of class "dgCMatrix" [1,] . . 1 1 1 1 1 1 1 . [2,] . . . . . . . . . . [3,] . . . . . . . . . . [4,] . . . . . . . . . . [5,] . . . . . . . . . . [6,] . . . . . . . . . . [7,] . . . . . . . . . . [8,] . . . . . . . . . . [9,] . . . . . . . . . . [10,] . . . . . . . . . . crossprod(M): 10 x 10 sparse Matrix of class "nsCMatrix" [1,] . . . . . . . . . . [2,] . | | | | | | | | | [3,] . | | | | | | | | | [4,] . | | | | | | | | | [5,] . | | | | | | | | | [6,] . | | | | | | | | | [7,] . | | | | | | | | | [8,] . | | | | | | | | | [9,] . | | | | | | | | | [10,] . | | | | | | | | | tcrossprod(M): 10 x 10 sparse Matrix of class "nsCMatrix" [1,] | | . . . . . . . . [2,] | | . . . . . . . . [3,] . . . . . . . . . . [4,] . . . . . . . . . . [5,] . . . . . . . . . . [6,] . . . . . . . . . . [7,] . . . . . . . . . . [8,] . . . . . . . . . . [9,] . . . . . . . . . . [10,] . . . . . . . . . . M %*% m: 10 x 10 Matrix of class "dgeMatrix" [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 0 0 1 1 1 1 1 1 1 0 [2,] 0 0 0 0 0 0 0 0 0 0 [3,] 0 0 0 0 0 0 0 0 0 0 [4,] 0 0 0 0 0 0 0 0 0 0 [5,] 0 0 0 0 0 0 0 0 0 0 [6,] 0 0 0 0 0 0 0 0 0 0 [7,] 0 0 0 0 0 0 0 0 0 0 [8,] 0 0 0 0 0 0 0 0 0 0 [9,] 0 0 0 0 0 0 0 0 0 0 [10,] 0 0 0 0 0 0 0 0 0 0 > chk.ngMatrix(tsM <- as(sM, "triangularMatrix")) # ntT 10 x 10 sparse Matrix of class "ntTMatrix" [1,] . | | | | | | | | | [2,] . . | | | | | | | . [3,] . . . . . . . . . . [4,] . . . . . . . . . . [5,] . . . . . . . . . . [6,] . . . . . . . . . . [7,] . . . . . . . . . . [8,] . . . . . . . . . . [9,] . . . . . . . . . . [10,] . . . . . . . . . . M %*% M: 10 x 10 sparse Matrix of class "dtCMatrix" [1,] . . 1 1 1 1 1 1 1 . [2,] . . . . . . . . . . [3,] . . . . . . . . . . [4,] . . . . . . . . . . [5,] . . . . . . . . . . [6,] . . . . . . . . . . [7,] . . . . . . . . . . [8,] . . . . . . . . . . [9,] . . . . . . . . . . [10,] . . . . . . . . . . crossprod(M): 10 x 10 sparse Matrix of class "nsCMatrix" [1,] . . . . . . . . . . [2,] . | | | | | | | | | [3,] . | | | | | | | | | [4,] . | | | | | | | | | [5,] . | | | | | | | | | [6,] . | | | | | | | | | [7,] . | | | | | | | | | [8,] . | | | | | | | | | [9,] . | | | | | | | | | [10,] . | | | | | | | | | tcrossprod(M): 10 x 10 sparse Matrix of class "nsCMatrix" [1,] | | . . . . . . . . [2,] | | . . . . . . . . [3,] . . . . . . . . . . [4,] . . . . . . . . . . [5,] . . . . . . . . . . [6,] . . . . . . . . . . [7,] . . . . . . . . . . [8,] . . . . . . . . . . [9,] . . . . . . . . . . [10,] . . . . . . . . . . M %*% m: 10 x 10 Matrix of class "dgeMatrix" [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 0 0 1 1 1 1 1 1 1 0 [2,] 0 0 0 0 0 0 0 0 0 0 [3,] 0 0 0 0 0 0 0 0 0 0 [4,] 0 0 0 0 0 0 0 0 0 0 [5,] 0 0 0 0 0 0 0 0 0 0 [6,] 0 0 0 0 0 0 0 0 0 0 [7,] 0 0 0 0 0 0 0 0 0 0 [8,] 0 0 0 0 0 0 0 0 0 0 [9,] 0 0 0 0 0 0 0 0 0 0 [10,] 0 0 0 0 0 0 0 0 0 0 > chk.ngMatrix(as( sM, "CsparseMatrix")) # ngC 10 x 10 sparse Matrix of class "ngCMatrix" [1,] . | | | | | | | | | [2,] . . | | | | | | | . [3,] . . . . . . . . . . [4,] . . . . . . . . . . [5,] . . . . . . . . . . [6,] . . . . . . . . . . [7,] . . . . . . . . . . [8,] . . . . . . . . . . [9,] . . . . . . . . . . [10,] . . . . . . . . . . M %*% M: 10 x 10 sparse Matrix of class "dgCMatrix" [1,] . . 1 1 1 1 1 1 1 . [2,] . . . . . . . . . . [3,] . . . . . . . . . . [4,] . . . . . . . . . . [5,] . . . . . . . . . . [6,] . . . . . . . . . . [7,] . . . . . . . . . . [8,] . . . . . . . . . . [9,] . . . . . . . . . . [10,] . . . . . . . . . . crossprod(M): 10 x 10 sparse Matrix of class "nsCMatrix" [1,] . . . . . . . . . . [2,] . | | | | | | | | | [3,] . | | | | | | | | | [4,] . | | | | | | | | | [5,] . | | | | | | | | | [6,] . | | | | | | | | | [7,] . | | | | | | | | | [8,] . | | | | | | | | | [9,] . | | | | | | | | | [10,] . | | | | | | | | | tcrossprod(M): 10 x 10 sparse Matrix of class "nsCMatrix" [1,] | | . . . . . . . . [2,] | | . . . . . . . . [3,] . . . . . . . . . . [4,] . . . . . . . . . . [5,] . . . . . . . . . . [6,] . . . . . . . . . . [7,] . . . . . . . . . . [8,] . . . . . . . . . . [9,] . . . . . . . . . . [10,] . . . . . . . . . . M %*% m: 10 x 10 Matrix of class "dgeMatrix" [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 0 0 1 1 1 1 1 1 1 0 [2,] 0 0 0 0 0 0 0 0 0 0 [3,] 0 0 0 0 0 0 0 0 0 0 [4,] 0 0 0 0 0 0 0 0 0 0 [5,] 0 0 0 0 0 0 0 0 0 0 [6,] 0 0 0 0 0 0 0 0 0 0 [7,] 0 0 0 0 0 0 0 0 0 0 [8,] 0 0 0 0 0 0 0 0 0 0 [9,] 0 0 0 0 0 0 0 0 0 0 [10,] 0 0 0 0 0 0 0 0 0 0 > chk.ngMatrix(as(tsM, "CsparseMatrix")) # ntC 10 x 10 sparse Matrix of class "ntCMatrix" [1,] . | | | | | | | | | [2,] . . | | | | | | | . [3,] . . . . . . . . . . [4,] . . . . . . . . . . [5,] . . . . . . . . . . [6,] . . . . . . . . . . [7,] . . . . . . . . . . [8,] . . . . . . . . . . [9,] . . . . . . . . . . [10,] . . . . . . . . . . M %*% M: 10 x 10 sparse Matrix of class "dtCMatrix" [1,] . . 1 1 1 1 1 1 1 . [2,] . . . . . . . . . . [3,] . . . . . . . . . . [4,] . . . . . . . . . . [5,] . . . . . . . . . . [6,] . . . . . . . . . . [7,] . . . . . . . . . . [8,] . . . . . . . . . . [9,] . . . . . . . . . . [10,] . . . . . . . . . . crossprod(M): 10 x 10 sparse Matrix of class "nsCMatrix" [1,] . . . . . . . . . . [2,] . | | | | | | | | | [3,] . | | | | | | | | | [4,] . | | | | | | | | | [5,] . | | | | | | | | | [6,] . | | | | | | | | | [7,] . | | | | | | | | | [8,] . | | | | | | | | | [9,] . | | | | | | | | | [10,] . | | | | | | | | | tcrossprod(M): 10 x 10 sparse Matrix of class "nsCMatrix" [1,] | | . . . . . . . . [2,] | | . . . . . . . . [3,] . . . . . . . . . . [4,] . . . . . . . . . . [5,] . . . . . . . . . . [6,] . . . . . . . . . . [7,] . . . . . . . . . . [8,] . . . . . . . . . . [9,] . . . . . . . . . . [10,] . . . . . . . . . . M %*% m: 10 x 10 Matrix of class "dgeMatrix" [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 0 0 1 1 1 1 1 1 1 0 [2,] 0 0 0 0 0 0 0 0 0 0 [3,] 0 0 0 0 0 0 0 0 0 0 [4,] 0 0 0 0 0 0 0 0 0 0 [5,] 0 0 0 0 0 0 0 0 0 0 [6,] 0 0 0 0 0 0 0 0 0 0 [7,] 0 0 0 0 0 0 0 0 0 0 [8,] 0 0 0 0 0 0 0 0 0 0 [9,] 0 0 0 0 0 0 0 0 0 0 [10,] 0 0 0 0 0 0 0 0 0 0 > > ## "l" --- logical Matrices -- use usual 0/1 arithmetic > nsM <- sM > sM <- as(sM, "lMatrix") > sm <- as(sM, "matrix") > stopifnot(identical(sm, as.matrix(nsM))) > stopifnotValid(sMM <- sM %*% sM, "dsparseMatrix") > assert.EQ.mat (sMM, sm %*% sm) > assert.EQ.mat(t(sM) %*% sM, + t(sm) %*% sm, tol=0) > stopifnotValid(cM <- crossprod(sM), "dsCMatrix") > stopifnotValid(tM <- tcrossprod(sM), "dsCMatrix") > stopifnot(identical(cM, as(t(sM) %*% sM, "symmetricMatrix")), + identical(tM, forceSymmetric(sM %*% t(sM)))) > assert.EQ.mat( cM, crossprod(sm)) > assert.EQ.mat( tM, as(tcrossprod(sm),"matrix")) > dm <- as(sM, "denseMatrix") > ## the following 6 products (dm o sM) all failed up to 2013-09-03 > stopifnotValid(dm %*% sM, "denseMatrix")## failed {missing coercion} > stopifnotValid(crossprod (dm , sM),"denseMatrix") > stopifnotValid(tcrossprod(dm , sM),"denseMatrix") > dm[2,1] <- TRUE # no longer triangular > stopifnotValid( dm %*% sM, "denseMatrix") > stopifnotValid(crossprod (dm , sM),"denseMatrix") > stopifnotValid(tcrossprod(dm , sM),"denseMatrix") > > ## A sparse example - with *integer* matrix: > M <- Matrix(cbind(c(1,0,-2,0,0,0,0,0,2.2,0), + c(2,0,0,1,0), 0, 0, c(0,0,8,0,0),0)) > t(M) 6 x 10 sparse Matrix of class "dgCMatrix" [1,] 1 . -2 . . . . . 2.2 . [2,] 2 . . 1 . 2 . . 1.0 . [3,] . . . . . . . . . . [4,] . . . . . . . . . . [5,] . . 8 . . . . 8 . . [6,] . . . . . . . . . . > (-4:5) %*% M 1 x 6 Matrix of class "dgeMatrix" [,1] [,2] [,3] [,4] [,5] [,6] [1,] 8.8 -3 0 0 8 0 > stopifnot(as.vector(print(t(M %*% 1:6))) == + c(as(M,"matrix") %*% 1:6)) 1 x 10 Matrix of class "dgeMatrix" [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [1,] 5 0 38 2 0 4 0 40 4.2 0 > (M.M <- crossprod(M)) 6 x 6 sparse Matrix of class "dsCMatrix" [1,] 9.84 4.2 . . -16 . [2,] 4.20 10.0 . . . . [3,] . . . . . . [4,] . . . . . . [5,] -16.00 . . . 128 . [6,] . . . . . . > MM. <- tcrossprod(M) > stopifnot(class(MM.) == "dsCMatrix", + class(M.M) == "dsCMatrix") > > M3 <- Matrix(c(rep(c(2,0),4),3), 3,3, sparse=TRUE) > I3 <- as(Diagonal(3), "CsparseMatrix") > m3 <- as.matrix(M3) > iM3 <- solve(m3) > stopifnot(all.equal(unname(iM3), matrix(c(3/2,0,-1,0,1/2,0,-1,0,1), 3))) > assert.EQ.mat(solve(as(M3, "sparseMatrix")), iM3) > assert.EQ.mat(solve(I3,I3), diag(3)) > assert.EQ.mat(solve(M3, I3), iM3)# was wrong because I3 is unit-diagonal > assert.EQ.mat(solve(m3, I3), iM3)# gave infinite recursion in (<=) 0.999375-10 > > stopifnot(identical(ttI3 <- crossprod(tru, I3), t(tru) %*% I3), + identical(tI3t <- crossprod(I3, tru), t(I3) %*% tru), + identical(I3tt <- tcrossprod(I3, tru), I3 %*% t(tru))) > I3@uplo # U pper triangular [1] "U" > tru@uplo# L ower triangular [1] "L" > ## "FIXME": These are all FALSE now; the first one *is* ok (L o U); the others *not* > isValid(tru %*% I3, "triangularMatrix") [1] FALSE > isValid(ttI3, "triangularMatrix") [1] TRUE > isValid(tI3t, "triangularMatrix") [1] TRUE > isValid(I3tt, "triangularMatrix") [1] TRUE > > ## even simpler > m <- matrix(0, 4,7); m[c(1, 3, 6, 9, 11, 22, 27)] <- 1 > (mm <- Matrix(m)) 4 x 7 sparse Matrix of class "dgCMatrix" [1,] 1 . 1 . . . . [2,] . 1 . . . 1 . [3,] 1 . 1 . . . 1 [4,] . . . . . . . > (cm <- Matrix(crossprod(m))) 7 x 7 sparse Matrix of class "dsCMatrix" [1,] 2 . 2 . . . 1 [2,] . 1 . . . 1 . [3,] 2 . 2 . . . 1 [4,] . . . . . . . [5,] . . . . . . . [6,] . 1 . . . 1 . [7,] 1 . 1 . . . 1 > stopifnot(identical(crossprod(mm), cm)) > (tm1 <- Matrix(tcrossprod(m))) #-> had bug in 'Matrix()' ! 4 x 4 sparse Matrix of class "dsCMatrix" [1,] 2 . 2 . [2,] . 2 . . [3,] 2 . 3 . [4,] . . . . > (tm2 <- tcrossprod(mm)) 4 x 4 sparse Matrix of class "dsCMatrix" [1,] 2 . 2 . [2,] . 2 . . [3,] 2 . 3 . [4,] . . . . > Im2 <- solve(tm2[-4,-4]) > P <- as(as.integer(c(4,1,3,2)),"pMatrix") > p <- as(P, "matrix") > P %*% mm 4 x 7 sparse Matrix of class "dgCMatrix" [1,] . . . . . . . [2,] 1 . 1 . . . . [3,] 1 . 1 . . . 1 [4,] . 1 . . . 1 . > assertError(mm %*% P) # dimension mismatch > assertError(m %*% P) # ditto > assertError(crossprod(t(mm), P)) # ditto > stopifnotValid(tm1, "dsCMatrix") > stopifnot(exprs = { + all.equal(tm1, tm2, tolerance = 1e-15) + all.equal(Im2 %*% tm2[1:3,], Matrix(cbind(diag(3), 0))) + identical(p, as.matrix(P)) + all(P %*% m == as.matrix(P) %*% m) + all(P %*% mm == P %*% m) + all(P %*% mm - P %*% m == 0) + all(t(mm) %*% P == t(m) %*% P) + all(crossprod(m, P) == crossprod(mm, P)) + }) > > d <- function(m) as(m,"dsparseMatrix") > IM1 <- as(c(3,1,2), "indMatrix") > IM2 <- as(c(1,2,1), "indMatrix") > assert.EQ.Mat(crossprod( IM1, IM2), + crossprod(d(IM1),d(IM2)), tol=0)# failed at first > iM <- as(cbind2(IM2, 0), "indMatrix") > assert.EQ.Mat(crossprod(iM), Diagonal(x = 2:0)) > assert.EQ.Mat(crossprod(iM, iM), Diagonal(x = 2:0)) > > N3 <- Diagonal(x=1:3) > U3 <- Diagonal(3) # unit diagonal (@diag = "U") > C3 <- as(N3, "CsparseMatrix") > lM <- as(IM2, "lMatrix") > nM <- as(IM2, "nMatrix") > nCM <- as(nM, "CsparseMatrix") > NM <- N3 %*% IM2 > NM. <- C3 %*% IM2 > stopifnot(Q.C.identical(NM, ## <- failed + d(N3) %*% d(IM2), checkClass=FALSE), + identical(NM, N3 %*% lM), + identical(NM, N3 %*% nM) + , ## all these "work" (but partly wrongly gave non-numeric Matrix: + Q.C.identical(NM, NM., checkClass=FALSE) + , + mQidentical(as.matrix(NM.), array(c(1, 0, 3, 0, 2, 0), dim=3:2)) + , + identical(NM., C3 %*% lM) + , + identical(NM., C3 %*% nM) # wrongly gave n*Matrix + , + isValid(U3 %*% IM2, "dsparseMatrix")# was l* + , + isValid(U3 %*% lM, "dsparseMatrix")# was l* + , + isValid(U3 %*% nM, "dsparseMatrix")# was n* + , + identical(C3 %*% nM -> C3n, # wrongly gave ngCMatrix + C3 %*% nCM) + , + isValid(C3n, "dgCMatrix") + , + identical3(U3 %*% IM2, # wrongly gave lgTMatrix + U3 %*% lM -> U3l, # ditto + U3 %*% nM) # wrongly gave ngTMatrix + , + isValid(U3l, "dgRMatrix") + ) > > selectMethod("%*%", c("dtCMatrix", "ngTMatrix")) # x %*% .T.2.C(y) --> Method Definition: function (x, y) .Call(R_sparse_matmult, x, y, FALSE, FALSE, FALSE, FALSE) Signatures: x y target "dtCMatrix" "ngTMatrix" defined "CsparseMatrix" "TsparseMatrix" > selectMethod("%*%", c("dtCMatrix", "ngCMatrix")) # .Call(Csparse_Csparse_prod, x, y) Method Definition: function (x, y) .Call(R_sparse_matmult, x, y, FALSE, FALSE, FALSE, FALSE) Signatures: x y target "dtCMatrix" "ngCMatrix" defined "CsparseMatrix" "CsparseMatrix" > selectMethod("%*%", c("ddiMatrix", "indMatrix")) # x %*% as(y, "lMatrix") -> Method Definition: function (x, y) { if (y@margin == 1L) return(x %*% .M2kind(y, "d")) matmultDim(x@Dim, y@Dim, type = 1L) r <- .M2kind(x[, y@perm, drop = FALSE], ",") r@Dimnames <- matmultDN(dimnames(x), y@Dimnames, type = 1L) r } Signatures: x y target "ddiMatrix" "indMatrix" defined "Matrix" "indMatrix" > selectMethod("%*%", c("ddiMatrix", "lgTMatrix")) # diagCspprod(as(x, "Csp.."), y) Method Definition: function (x, y) .Call(R_diagonal_matmult, x, y, FALSE, FALSE, FALSE) Signatures: x y target "ddiMatrix" "lgTMatrix" defined "diagonalMatrix" "TsparseMatrix" > selectMethod("%*%", c("ddiMatrix", "ngTMatrix")) # (ditto) Method Definition: function (x, y) .Call(R_diagonal_matmult, x, y, FALSE, FALSE, FALSE) Signatures: x y target "ddiMatrix" "ngTMatrix" defined "diagonalMatrix" "TsparseMatrix" > > stopifnot( + isValid(show(crossprod(C3, nM)), "dgCMatrix"), # wrongly gave ngCMatrix + identical3(## the next 4 should give the same (since C3 and U3 are symmetric): + show(crossprod(U3, IM2)),# wrongly gave ngCMatrix + crossprod(U3, nM), # ditto + crossprod(U3, lM))) # wrongly gave lgCMatrix 3 x 2 sparse Matrix of class "dgCMatrix" [1,] 1 . [2,] . 2 [3,] 3 . 3 x 2 sparse Matrix of class "dgRMatrix" [1,] 1 . [2,] . 1 [3,] 1 . > > > set.seed(123) > for(n in 1:250) { + n1 <- 2 + rpois(1, 10) + n2 <- 2 + rpois(1, 10) + N <- rpois(1, 25) + ii <- seq_len(N + min(n1,n2)) + IM1 <- as(c(sample(n1), sample(n1, N, replace=TRUE))[ii], "indMatrix") + IM2 <- as(c(sample(n2), sample(n2, N, replace=TRUE))[ii], "indMatrix") + ## stopifnot(identical(crossprod( IM1, IM2), + ## crossprod(d(IM1), d(IM2)))) + if(!identical(C1 <- crossprod( IM1, IM2 ), + CC <- crossprod(d(IM1), d(IM2))) && + !all(C1 == CC)) { + cat("The two crossprod()s differ: C1 - CC =\n") + print(C1 - CC) + stop("The two crossprod()s differ!") + } else if(n %% 25 == 0) cat(n, " ") + }; cat("\n") 25 50 75 100 125 150 175 200 225 250 > > ## two with an empty column --- these failed till 2014-06-14 > X <- as(c(1,3,4,5,3), "indMatrix") > Y <- as(c(2,3,4,2,2), "indMatrix") > > ## kronecker: > stopifnot(identical(X %x% Y, + as(as.matrix(X) %x% as.matrix(Y), "indMatrix"))) > ## crossprod: > (XtY <- crossprod(X, Y))# gave warning in Matrix 1.1-3 5 x 4 sparse Matrix of class "dgTMatrix" [1,] . 1 . . [2,] . . . . [3,] . 1 1 . [4,] . . . 1 [5,] . 1 . . > XtY_ok <- as(crossprod(as.matrix(X), as.matrix(Y)), "TsparseMatrix") > assert.EQ.Mat(XtY, XtY_ok) # not true, previously > > ###------- %&% -------- Boolean Arithmetic Matrix products > > x5 <- c(2,0,0,1,4) > D5 <- Diagonal(x=x5) > N5 <- as(D5 != 0, "nMatrix") ## an "ndiMatrix" > D. <- Diagonal(x=c(TRUE,FALSE,TRUE,TRUE,TRUE)) > stopifnot(identical(D5 %&% D., N5)) > stopifnot(identical(D5 %&% as(D.,"CsparseMatrix"), + as(N5,"CsparseMatrix"))) > > set.seed(7) > L <- Matrix(rnorm(20) > 1, 4,5) > (N <- as(L, "nMatrix")) 4 x 5 sparse Matrix of class "ngCMatrix" [1,] | . . | . [2,] . . | . . [3,] . . . | . [4,] . . | . . > D <- Matrix(round(rnorm(30)), 5,6) # "dge", values in -1:1 (for this seed) > L %&% D 4 x 6 sparse Matrix of class "ngCMatrix" [1,] | | | | | . [2,] | | . | . | [3,] | | . . | . [4,] | | . | . | > stopifnot(identical(L %&% D, N %&% D), + all(L %&% D == as((L %*% abs(D)) > 0, "sparseMatrix"))) > stopifnotValid(show(crossprod(N )) , "nsCMatrix") # (TRUE/FALSE : boolean arithmetic) 5 x 5 sparse Matrix of class "nsCMatrix" [1,] | . . | . [2,] . . . . . [3,] . . | . . [4,] | . . | . [5,] . . . . . > stopifnotValid(show(crossprod(N +0)) -> cN0, "dsCMatrix") # -> numeric Matrix (with same "pattern") 5 x 5 sparse Matrix of class "dsCMatrix" [1,] 1 . . 1 . [2,] . . . . . [3,] . . 2 . . [4,] 1 . . 2 . [5,] . . . . . > stopifnot(all(crossprod(N) == t(N) %&% N), + identical(crossprod(N, boolArith=TRUE) -> cN., + as(cN0 != 0, "nMatrix")), + identical (cN., crossprod(L, boolArith=TRUE)), + identical3(cN0, crossprod(L), crossprod(L, boolArith=FALSE)) + ) > stopifnotValid(cD <- crossprod(D, boolArith = TRUE), "nsCMatrix") # sparse: "for now" > ## another slightly differing test "series" > L.L <- crossprod(L) > (NN <- as(L.L > 0,"nMatrix")) 5 x 5 sparse Matrix of class "nsCMatrix" [1,] | . . | . [2,] . . . . . [3,] . . | . . [4,] | . . | . [5,] . . . . . > nsy <- as(NN,"denseMatrix") > stopifnot(identical(NN, crossprod(NN)))# here > stopifnotValid(csy <- crossprod(nsy), "nsCMatrix") > stopifnotValid(csy. <- crossprod(nsy, boolArith=TRUE),"nsCMatrix") > stopifnot(all((csy > 0) == csy.), + all(csy. == (nsy %&% nsy))) > > ## for "many" more seeds: > set.seed(7); for(nn in 1:256) { + L <- Matrix(rnorm(20) > 1, 4,5) + D <- Matrix(round(rnorm(30)), 5,6) + stopifnot(all(L %&% D == as((L %*% abs(D)) > 0, "sparseMatrix"))) + } > > ## [Diagonal] o [0-rows/colums] : > m20 <- matrix(nrow = 2, ncol = 0); m02 <- t(m20) > M20 <- Matrix(nrow = 2, ncol = 0); M02 <- t(M20) > stopifnot(identical(dim(Diagonal(x=c(1,2)) %*% m20), c(2L, 0L)), + identical(dim(Diagonal(2) %*% M20), c(2L, 0L)), + identical(dim(Diagonal(x=2:1) %*% M20), c(2L, 0L))) > stopifnot(identical(dim(m02 %*% Diagonal(x=c(1,2))), c(0L, 2L)), + identical(dim(M02 %*% Diagonal(2) ), c(0L, 2L)), + identical(dim(M02 %*% Diagonal(x=2:1) ), c(0L, 2L))) > > ## RsparseMatrix --- Arko Bose (Jan.2022): "Method for %*% " > (m <- Matrix(c(0,0,2:0), 3,5)) 3 x 5 sparse Matrix of class "dgCMatrix" [1,] . 1 . . 2 [2,] . . 2 . 1 [3,] 2 . 1 . . > stopifnotValid(R <- as(m, "RsparseMatrix"), "RsparseMatrix") > stopifnotValid(T <- as(m, "TsparseMatrix"), "TsparseMatrix") > stopifnot(exprs = { + all.equal(as(t(R) %*% R, "symmetricMatrix"), crossprod(R) -> cR) + all.equal(as(R %*% t(R), "symmetricMatrix"), tcrossprod(R) -> tR) + all.equal(as(R %*% t(m), "symmetricMatrix"), as(tcrossprod(m), "RsparseMatrix")) + all.equal(as(m %*% t(R), "symmetricMatrix"), as(tcrossprod(m), "CsparseMatrix")) + ## failed in Matrix <= 1.4.1 (since 1.2.0, when 'boolArith' was introduced): + all.equal(as(cR, "RsparseMatrix"), as( crossprod(R, T), "symmetricMatrix")) + all.equal(as(cR, "CsparseMatrix"), as( crossprod(T, R), "symmetricMatrix")) + all.equal(as(tR, "RsparseMatrix"), as(tcrossprod(R, T), "symmetricMatrix")) + all.equal(as(tR, "CsparseMatrix"), as(tcrossprod(T, R), "symmetricMatrix")) + }) > > ## More for kronecker() ------------------------------------------------ > > checkKronecker <- function(X, Y, ...) { + k1 <- as(k0 <- kronecker(X, Y, ...), "matrix") + k2 <- kronecker(as(X, "matrix"), as(Y, "matrix"), ...) + cldX <- getClassDef(class(X)) + cldY <- getClassDef(class(Y)) + if(extends(cldX, "indMatrix") && + extends(cldY, "indMatrix")) { + stopifnot(is(k0, "indMatrix")) + storage.mode(k2) <- "logical" + } + if(extends(cldX, "triangularMatrix") && + extends(cldY, "triangularMatrix") && X@uplo == Y@uplo) + stopifnot(is(k0, "triangularMatrix"), + k0@uplo == X@uplo, + k0@diag == (if(X@diag == "N" || Y@diag == "N") "N" else "U")) + else if(extends(cldX, "symmetricMatrix") && + extends(cldY, "symmetricMatrix")) + stopifnot(is(k0, "symmetricMatrix"), k0@uplo == X@uplo) + ## could test for more special cases + if(!isTRUE(ae <- all.equal(k1, k2))) + stop(sprintf("checkKronecker(<%s>, <%s>):\n %s", + class(X), class(Y), paste0(ae, collapse = "\n "))) + } > > set.seed(145133) > lXY <- lapply(rpois(2L, 10), + function(n) { + x <- rsparsematrix(n, n, 0.6) + dn <- replicate(2L, if(sample(c(FALSE, TRUE), 1L)) + sample(letters, n, TRUE), + simplify = FALSE) + x@Dimnames <- dn + x4 <- list(as(as(x, "dMatrix"), "denseMatrix"), + as(as(x, "dMatrix"), "CsparseMatrix"), + as(as(x, "lMatrix"), "RsparseMatrix"), + as(as(x, "nMatrix"), "TsparseMatrix")) + + mkList <- function(y) + list(x, + triu(x), + tril(x), + { z <- triu(x, 1L); z@diag <- "U"; z }, + { z <- tril(x, -1L); z@diag <- "U"; z }, + forceSymmetric(x, "U"), + forceSymmetric(x, "L"), + { z <- Diagonal(x = diag(x)); z@Dimnames <- dn; z }, + as(sample.int(n, replace = TRUE), "indMatrix"), + as(x, "matrix")) + unlist(lapply(x4, mkList), FALSE, FALSE) + }) > lX <- lXY[[1L]] > lY <- lXY[[2L]] > for(i in seq_along(lX)) + for(j in seq_along(lY)) + checkKronecker(lX[[i]], lY[[j]], make.dimnames = TRUE) > > cat('Time elapsed: ', proc.time(),'\n') # for ``statistical reasons'' Time elapsed: 7.14 0.71 8.34 NA NA > > proc.time() user system elapsed 7.14 0.71 8.34