yixuan / RSpectra

R Interface to the Spectra Library for Large Scale Eigenvalue and SVD Problems
http://cran.r-project.org/package=RSpectra
80 stars 12 forks source link

Support dsCMatrix/dsRMatrix #15

Closed flying-sheep closed 5 years ago

flying-sheep commented 5 years ago

The dsCMatrix and dsRMatrix classes are the only useful ones you don’t support yet. They’re like dgCMatrix/dgRMatrix except that only one triangle is stored and a uplo ∈ {'U', 'L'} slot exists to signify which.

The dgRMatrix docs say “These "..RMatrix" classes are currently still mostly unimplemented!”, but it still can’t hurt to implement a 2-line method for them.

PS: Thank you for the amazing package, I wish it had been there from the beginning!

yixuan commented 5 years ago

Yeah, adding new matrix classes should be a relatively easy task. However, the downside for these two is that they do not have a counterpart in the Eigen library. It may take some effort to implement all the matrix operations, but I'll gradually add them if I get some time. Also you are more than welcome to contribute. 😃

flying-sheep commented 5 years ago

Well, eigs_sym only uses a triangle, right? so you could just do this, right?

##' @rdname eigs
##' @export
eigs.dsCMatrix <- function(A, k, which = "LM", opts = list(), ...)
{
    if (! which %in% c("LM", "SM", "LR", "SR"))
        stop('Symmetric matrices can only support `which %in% c("LM", "SM", "LR", "SR")`')
    if (which == "LR") which <- "LA"
    if (which == "SR") which <- "SA"
    # I think you can just use it as a dgCMatrix without adding a mattype here
    eigs_real_sym(A, nrow(A), k, which, sigma, opts, mattype = "sym_dgCMatrix",
                  extra_args = list(use_lower = A@uplo == "L"))
}

Also you are more than welcome to contribute. :smiley:

With pleasure! If above code is about correct, I’ll happily send a PR.

yixuan commented 5 years ago

The R code is good. What I mean is that the C++ part needs to correctly handle the dsCMatrix object, whose internal structure may be different from dgCMatrix.

privefl commented 5 years ago

You could always use the function interface, right?

library(Matrix)
mat <- sparseMatrix(i = 1:2, j = c(2, 2), x = 1:2, dims = c(3, 3))
(mat2 <- forceSymmetric(mat))
RSpectra::eigs_sym(mat2, k = 1)                                # NOK
RSpectra::eigs_sym(function(x, args) as.vector(mat2 %*% x), 
                   k = 1, n = nrow(mat2))                      # OK

Would this be slower to do it like this?

flying-sheep commented 5 years ago

@yixuan look at the comment in the R code: I think it’s the exact same thing except for the uplo property and the fact that only a triangle is specified. Which matches exactly your usage of sym_dgCMatrix. (I think you could just rename that to dsCMatrix)

@privefl as said: since I think that dsCMatrix is exactly what @yixuan calls sym_dgCMatrix (with an added uplo attribute), it’s just much more elegant to directly use that.

yixuan commented 5 years ago

Great. I'll take a look later.

yixuan commented 5 years ago

I've checked the structure of dsCMatrix and can confirm that it is compatible with dgCMatrix. So, yes, go ahead and submit a PR. :smiley:

flying-sheep commented 5 years ago

I would do the following:

  1. Check if you can somehow coerce a dgCMatrix to a dsCMatrix while setting the uplo slot and not copying the data
  2. Check if that works if more than that triangle is filled in the dgCmatrix
  3. Change the internal names sym_dgCMatrixdsCMatrix and the lower parameter to use the uplo slot instead
  4. change the R methods so that eigs.dgCMatrix delegates to eigs_sym.dsCMatrix if isSymmetric(A) and so on
alexpghayes commented 5 years ago

Thanks for starting the PR, this will be very useful for a project I'm currently working on!