##-*- mode: org -*-
–> see SM (3e6 x 1023) ex. in tests/Simple.R
Matrix-Bugs [#6729] 2021-06 ./tests/AAA_latest.R: provide .KhatriRao() for “general” (notably dense,complex,..) matrices
see FIXME in symmetricMatrix_validate() in src/dsyMatrix.c —> see TODO below: `forceCspSymmetric()`
`unique()` and `duplicated()` methods for “sparseVector” & “(sparse)Matrix”; have “*.matrix” S3 methods
–> tests/Simple.R {‘max.print’}
Why change? Originally “ddiMatrix” etc extended denseMatrix but now sparseMatrix Currently ‘doDiag’s documentation starts as ‘doDiag: only when ‘sparse = FALSE’, ....’ This would change, and doDiag would be active also for ‘sparse = TRUE’
also mentioning forceSymmetric(); maybe that wd “inherit” arg. ‘symDimnames’ in {T,F,NA} from forceCspSymmetric().
https://r-forge.r-project.org/tracker/index.php?func=detail&aid=1330&group_id=61&atid=294
sparse.model.matrix(.) bug with NA’s & na.action = “na.pass” => ~/R/MM/Pkg-ex/Matrix/sparse.model-bugs_EB.R
sparse.model.matrix(~ poly(x,3)) gives wrong column names => ~/R/MM/Pkg-ex/Matrix/sparse.model-bugs_EB.R
R/
------- e.g., for src/Csparse.c, Csp_dense_products(…) around lines 600
Definitely should use it for solve(<dgCMatrix>) {it seems MATLAB does}; alternatively also is_sym() [in src/cs_utils.c], see below.
“ddiMatrix” does not inherit from “dsparseMatrix”, nor does “ldiMatrix” from “lparseMatrix”. Seems an undesirable inconsistency. Try changing setClass(“ddiMatrix”, contains = c(“diagonalMatrix”, “dMatrix”)) to setClass(“ddiMatrix”, contains = c(“diagonalMatrix”, “dsparseMatrix”))
–> ~/R/MM/Pkg-ex/Matrix/sparseOrderedLogit.R
so people like RP may stop whining about its non-availability in “base R” (2015-07-09)
level e.g. for solve(<dsyMatrix>). NB: is a generalized LDL’ [with pivoting!]. Should define expand() method or similar, see man/BunchKaufman-methods.Rd and R/dsyMatrix.R (at end).
leave the for loops, as soon as is_lower == is_upper == 0.
when one of the arguments is diagonal
it first copies the matrix to a dgeMatrix {re-allocating!}, then goes to sparse via cholmod_(l_)dense_to_sparse. ==>
matrix_to_Csparse() plus .m2dgC() etc R short cuts
extend C’s matrix_to_Csparse() to optionally check for diagonal, (upper or lower) triangular, and/or symmetric case
a pattern, i.e., nsparseMatrix as well because cholmod_ssmult() just does that even if only one of the two matrices is a pattern matrix. The latter case is really wrong. The above behavior seems many years old.. and sometimes is interesting and good, using Boolean arithmetic: T+T := T|T = T
For 1.2-0, changed the result to return numeric when one of the two matrices is not nsparse. ==> Provide the previous functionality via a Matrix package R function: ==> We’ve introduced ‘%&%’ for Matrix 1.2-0 and ‘boolArith = TRUE’ for crossprod/tcrossprod.
behave the same as for ‘ndenseMatrix’ or ‘lsparseMatrix’ or traditional logical matrices?
see misc/products-Mv.R and *.Rout – now tests/matprod.R (“3.2.0”)
qr.coef() has wrong (column)names, even in full-rank case: see man/qr-methods.Rd (“FIXME”); maybe related to
drop0(R. <- qr.R(qx), tol=1e-15) # columns are int b1 c1 c2 b2 c3 {in man/qr-methods.Rd}
Bigger drawbacks: Really I have to define Matrix methods for base functions that just worked perfectly via as.matrix
==> I’ve introduced “Matrix” S4 methods (and hence made eigen() S4 generic)
{as eigen()} would also need “Matrix” S4 methods
and now warns, also, e.g., from rcond(<non-square dgeMatrix>)
2) base :: lower.tri() and upper.tri() also use as.matrix() but are not generic => would need to make them S4 genric
for use in diverse ./tests/*.R.
---- give substitute .asmatrix(), but that’s not faster; simply calls as(*,”matrix”)
of some samples, i.e., columns of all 0 s. It’s mathematically more natural –> typically will be useful.
-> ~/R/MM/Pkg-ex/Matrix/cor_sparse-propos.R <- http://stackoverflow.com/questions/5888287/ -> ~/R/MM/Pkg-ex/Matrix/cor_cos.R and ~/R/MM/Pkg-ex/Matrix/cor_cos_testing Provide cor.sparse() and other association measures for sparse matrices.
typically sparse ones: Have rsparseMatrix() [exported] already; then rspMat(), rUnitTri(), mkLDL() [!] in inst/test-tools-Matrix.R ; then, e.g., rBlockTri() in man/bdiag.Rd. (man/* ?; tests/* )
from Jeremy D Silver, per posts to R-devel on Aug.26,2011 {MM: ~/R/MM/Pkg-ex/Matrix/bands-Jeremy_Silver-ex.R }
in the sense that if ‘dg.mat’ is a (“LINPACK/EISPACK”-format) dense (n x K) matrix containing K diagonals, and BS <- bandSparse(.., diagonals=dg.mat); dg.m <- getbands(BS,..) would exactly return the ‘dg.mat’ matrix.
happen in C, see R/bind2.R and src/Csparse.c (Csparse_horzcat etc).
third calls in the Synopsis should be to ldexpf and ldexpl.
2012-07: DONE with Ops, etc, also pack() / unpack(); not yet: “Math”
(well, this is “related to” the fact that we do not have too many packed matrix methods).
solution wrt special forms by using a ‘right’ argument and a ‘classed’ argument. [done with dgeMatrix_matrix_mm(); not yet for other classes; and for _crossprod()]
via .set.factors()
depends on (perm, LDL, super) arguments:
Can I find a uniform language ‘1-based indexing’ or ‘0-origin indexing’ ?
-> src/Mutils.c
when appropriate: particularly desirable for “Math2”: round(), signif()
(Golub+Van Loan 1996, 3.1.8, p.93)”
“d” -> “l” : drops the ‘x’ slot “l” -> “d” : construct an ‘x’ slot of all ‘1’ We currently have many of these conversions explicitly, e.g. setAs(“dsTMatrix”, “lsTMatrix”, function(from) new(“lsTMatrix”, i = from@i, j = from@j, uplo = from@uplo, Dim = from@Dim, Dimnames = from@Dimnames)) but I would rather want to automatically construct all these coercion methods at once by a “method constructor”, i.e., for all “dsparse*” -> “lsparse*” and vice versa. How can one do this {in a documented way} ?
basically enable all “sensible” as(fromMatrix, toMatrix) calls, possibly using canCoerce(.)
==> implement setAs(“dspMatrix”, “sparseMatrix”) setAs(“dppMatrix”, “sparseMatrix”) setAs(“dtpMatrix”, “sparseMatrix”) and the same for “lsp” , “ltp” and “nsp” , “ntp” !
calling Lapack’s DGEMM for “dense” [2005-12-xx: done for dgeMatrix at least]
both with ‘tol = 0.’ to drop “values that happen to be 0” and for zapsmall() methods for Csparse*
in src/CHOLMOD/Include/cholmod_matrixops.h : for another function specifically for multiplying a cholmod_sparse object by a diagonal matrix. Use it in %*% and [t]crossprod methods.
e.g. zapsmall(<pMatrix>) fails “badly”
returns <dense>, since –> Csparse_dense_prod –> cholmod_sdmult and that does only return dense. When the sparse matrix is very sparse, i.e. has many rows with only zero entries, it would make much sense to return sparse.
!M where M is “sparseMatrix”, currently always gives dense. This only makes sense when M is “really sparse”.
This (`[<-` method) is now “smart” for diagonalMatrix, but needs also to be for triangularMatrix, and probably also “dense*general*Matrix” since the above currently goes via “matrix” and back instead of using the ‘x’ slot directly; in particular, the triangular* “class property” is lost! [current ??]
when the replacement region is large, or also when do many millions of one-element assignments (say in a 100’000^2 Matrix).
(mostly done ?)
shows that chol() does not seem to make use of an already present factorization and rather uses one with more ‘0’ in x slot.
probably rather in man/CHMfactor-class.Rd than man/Cholesky.Rd
{diagonal entries of triangular Csparse} –> see comment in determinant(<dsC>) in R/dsCMatrix.R, will give faster determinant
we really mean to factorize AA’ ( + beta * I)
with non-symmetric A and then AA’ + mult * I is really meant. .updateCHMfactor() ## allows that already(?)
Cholesky(......, Imult), also tests for hidden {hence no examples} ldetL2up() { R/CHMfactor.R }; see ex in man/wrld_1deg.Rd MM: See e.g. ~/R/MM/Pkg-ex/Matrix/CholUpdate.R – for solve(<CHM>, <type>)
src/Csparse.c’s diag_tC_ptr() .
- diag_tC_ptr() functionality now exported via R/dsCMatrix.R .diag.dsC() : the name is silly, but functionality nice. See (hidden) example in man/Cholesky.Rd
but should give the symbolic factorization; similarly Cholesky(.) is not enabled
frequently in some computations {conditional normal distr. …}. Since this can be done faster than by sum(diag(A %*% B)) even for traditional matrices, e.g. sum(A * t(B)) or {sometimes even faster for “full” mat} crossprod(as.vector(A), as.vector(t(B))) and even more so for, e.g. <sparse> %*% <dense> {used in Soeren’s ‘gR’ computations}, we should also provide a generic and methods.
but as the above tricks show, is not really. Still, it’s well worth to provide diag.prod(A, B):
Well, if A %*% B is square, diag(A %*% B) === colSums(t(A) * B) and we should probably teach people about that !
but also for symmetric -> dsyMatrix [LAPACK dsyev() uses UPLO !], but also simply for dgeMatrix (without going via tradition matrices). What about Sparse? There’s fill-in, but it may still be sensible, e.g. mlist <- list(1, 2:3, diag(x=5:3), 27, cbind(1,3:6), 100:101) ee <- eigen(tcrossprod(bdiag(lapply(mlist, as.matrix)))) Matrix( signif(ee$vectors, 3) )
sparse or dense: “qr.R(<sparse>) may differ from qr.R(<dense>) because of permutations”
Users are still a bit frustrated and it currently influences rcond() as well.
can we go via qr() in any case? In some cases, e.g. lmer()’s “Lambda” (block triangular, small blocks) rcond(L) := 1 / (norm(L) * norm(solve(L))) is simple {and remains sparse, as solve(L) is still block triangular}
(instead of the infamous qr.qy, qr.qty,.... functions)
and could be made more efficient. Consider going via asTuniq() or something very close to .Arith.Csparse() in R/Ops.R For a traditional “matrix” object, we should speedup, using C code ..
replCmat() uses as_CspClass() and drop0(.) which itself call as_CspClass() quite a bit. –> try to replace these by as(*, “CsparseMatrix”); forceSymmetric, etc.
assumed to be a string, i.e. cannot be a connection. An R (instead of C) version should be pretty simple, and would work with connections automatically [“lsparse” become either “real” or “pattern”, “depending if they have NAs or not].
could return <diagMatrix> in the same cases where <diagMatrix> o <numeric> does.
our sparse matrices. Maybe this needs to be re-implemented using CHOLMOD routines.
as per Csparse_validate() and the tests in tests/validObj.R. This is as in CHOLMOD/CSparse, where nzmax (>= .@p[n]) corresponds to length(.@i), and makes sense e.g. for M[.,.] <- v assignments which could allocate in chunks and would not need to re-allocate anything in many cases. HOWEVER, replCmat() in R/Csparse.R is still far from making use of that.
the intermediate ‘M1 == M2’ typically is dense, hence potentially using humongous amount of memory. We should/could devise something like allCompare(M1, M2, `==`) which would remain sparse in all its computations.
packages. It may be better simply to add all the src/<nm>/Include directories to the include path for all compilations. I don’t think there is a big overhead. Right now we need to modify the include file src/SPQR/Include/SuiteSparseQR_C.h so that it does not expect to have src/UFsparse and src/CHOLMOD/Include on the include path. Maybe just those two should be added to the include path.
0-dimensional input themselves; LAPACK gives an integer error code)
[ in tests/indexing.R ]. Still would be nice to be able to use abIndex (see replTmat in R/Tsparse.R)
Sept. 2009: Subject: chol2inv() |-> solve(<CHMfactor>)
when testing and documenting chol2inv(), I found that it’s pretty simple to also define a method for “CHMfactor” objects, namely simply the solve(*, Diagonal(.) “A”) method. This is not particularly exciting, and also does not, I think help for defining a chol2inv() method for sparse (upper) triangular matrices.
Note that defining xtfrm() does not “help” (as sort() then goes via dense index). See “mean” in R/Matrix.R
correct and equivalent to src/CHOLMOD/Include/cholmod_core.h and siblings ??? {currently need to do this manually (Emacs M-x compare-windows) for the typedefs}
ensuring that the size allocated cannot grow with the vector/matrix/nnzero sizes of the input. [see the change needed in svn r2770 in src/dtCMatrix.c !]