[Question] matrix multiplication in rkeops

[yuhanH]

Hi, Thanks for your great package! I am trying to do a matrix multiplication in rkeops. But I cannot find an example... Could you show me an example?

[jeanfeydy]

Hi @yuhanH ,

Thanks for your interest in the library - and apologies for the delayed answer. As of today, RKeOps only provides support for the "low-level" Genred interface, without explicit support for LazyTensors. Implementing a matrix multiplication a = K * b in RKeOps is still do-able, but as a Sum_Reduction of a formula k(x_i, y_j) * b_j, where k(x_i, y_j) denotes the coefficients of your matrix and b_j refers to the "right-hand side" of your matrix product. An example is provided here. Does it fit your needs?

Best regards, Jean

[AmelieVernay]

Hello @yuhanH ,

I am currently co-working with @chloesrcb on LazyTensor support in RKeOps. Changes will be available in the coming weeks. Documentation and examples are coming as well, hopefully it will answer your question.

Best regards,

Amélie Vernay

[yuhanH]

HI @jeanfeydy and @adam-coogan Thanks for your reply! Actually, I have two following questions here.

1. matrix multiplication. Put a linear kernel between X and Y, and set X as an identity matrix to do the Y*B. ​Is this what you mean?

   
   formula = "Sum_Reduction( Sum(x* y) * b, 0)"
   args = c("x = Vi(2000)",     
        "y = Vj(2000)", 
        "b = Vj(3000)")      
   op <- keops_kernel(formula, args)

dim <- 2000 dim.b <- 3000 ny <- 5000

X <- diag(x = 1, nrow = dim) # Identity matrix Y <- matrix(runif(nydim), nrow=ny) # matrix ny x dim B <- matrix(runif(nydim.b), nrow=ny) # matrix ny x dim.b res <- op(list(X, Y, B))

If it is correct, I notice that when `dim.b` is big (such as 2000), this method works much faster than the normal matrix multiplication `t(Y)%*% B`. However, when `dim.b` is small (such as 20), this method works slower. Is it expected or anything I implemented incorrectly?

2. Sparse matrix. I think current rkeops doesn't support sparse matrix. I am wondering if you will consider a sparse matrix operation or your LazyTensor will handle sparse matrix directly. 
I am very looking forward to trying your LazyTensor work!

[jeanfeydy]

Hi @yuhanH ,

You’re very welcome! To answer your questions:

1. What is the exact computation that you are trying to perform? KeOps is meant to accelerate matrix operations with matrices that have a simple “symbolic” structure, i.e. whose coefficients are given by a mathematical formula that is evaluated on small vectors. (Think: a distance or kernel matrix.) If you want to accelerate a full, dense matrix-matrix product, the simplest option may be to directly link R to NVBLAS or rely on reticulate to make a call to a GPU-enabled Python framework such as cuPy or PyTorch for the compute-intensive parts of your program.
2. KeOps supports “block-sparse” matrices (as detailed here, these stream well on GPUs) but will probably never support “generic” sparse matrices: existing frameworks already support these objects very well.

What do you think? Best regards, Jean

[yuhanH]

Hi, @jeanfeydy Thanks for your reply! Both of them make sense. My exact need here is to use KNN graph/kernel to smooth a matrix B. I know ArgKMin_Reduction can calculate KNN neighbors index. Is it possible to implement KNN neighbors index as a kernel of X and Y in RKeOps? Thanks for your help!

[gdurif]

@yuhanH check this https://github.com/gdurif/optimizeR/blob/master/keops/tutorial/benchmarks/kNN.R for kNN implementation with RKeOps.

[yuhanH]

@gdurif Get it. Thanks for your example!