Input matrix

You can use two type of input matrices with RKeOps:

• ones whose rows (or columns) are indexed by \(i=1,...,M\) such as \(\mathbf X = [x_{ik}]_{M \times D}\)
  

• others whose rows (or columns) are indexed by \(j=1,...,N\) such as \(\mathbf Y = [y_{ik'}]_{N \times D'}\)
  
  

The dimensions over indexes \(i\) or \(j\) are called the outer dimensions (i.e. \(M\) or \(N\)). The other dimensions (i.e. \(D\) or \(D'\)) are called the inner dimensions. These terms refer to the contiguity of the data in memory:

• Outer dimensions \(M\) and \(N\) (over indexes \(i\) and \(j\) respectively) can be very large, even to large for GPU memory.
  

• Inner dimensions \(D\) and \(D'\) should be small enough to fit in GPU memory, in particular to ensure data colocality and avoid useless memory transfers. Corresponding columns (or rows) should be contiguous in memory (this point is handled for you in RKeOps, see this section).
  

Note 1: The outer dimension can correspond to the rows or the columns of the input matrices (and vice-versa for the inner dimension). The optimal orientation of input matrices is discussed in this section .

Note 2: All matrices indexed by \(i\) should have the same outer dimension \(M\) over \(i\), same for all matrices indexed by \(j\) (outer dimension \(N\)). Only the inner dimensions \(D\) and \(D'\) should be known for the compilation of your operators. The respective outer dimensions \(M\) and \(N\) are set at runtime (and can change from one run to another).

Notations

Input arguments of the formula are defined by using keywords, they can be of different types:

You should provide a vector of text string specifying the name and the type of all arguments in your formula.

Each keyword takes as parameter the inner dimension of the corresponding object. For instance, to define an input variable indexed by \(i\) corresponding to a \(D\)-dimensional vector, you can use "Vi(D)", same for a \(D\)-dimensional variable indexed by \(j\) being "Vj(D)" or a \(D\)-dimensional parameter "Pm(D)".

The vector of arguments should be

args = c("<name1>=<type1>(dim1)", "<name2>=<type2>(dim2)", "<nameX>=<typeX>(dimX)")

where

• <nameX> is the name
• <type1> is the type (among Vi, Vj or Pm)
• <dimX> is the inner dimension

of the X\(^\text{th}\) variable in the formula.

Important: The names should correspond to the ones used in the formula. The input parameter order will be the one used when calling the compiled operator.

Example: We define the corresponding arguments of the previous formula, i.e. parameters or variables indexed by \(i\) or \(j\) with their corresponding inner dimensions:

args = c("x = Vi(3)",      # vector indexed by i (of dim 3)
         "y = Vj(3)",      # vector indexed by j (of dim 3)
         "b = Vj(6)",      # vector indexed by j (of dim 6)
         "s = Pm(1)")      # parameter (scalar) 

Compile options

• use_cuda_if_possible: by default, user-defined operators are compiled for GPU if CUDA is available (and compiled for CPU otherwise).

# enable compiling for GPU if available (not necessary if using default options)
compile4gpu()
# or equivalently
set_rkeops_option("use_cuda_if_possible", 1)
# disable compiling for GPU
set_rkeops_option("use_cuda_if_possible", 0)

• precision: by default, user-defined operators are compiled to use float 32bits for computations (faster than float 64bits or double, compensated sum is available to reduce errors inherent to float 32bits operations)

set_rkeops_option("precision", "float")    # float 32bits (default)
set_rkeops_option("precision", "double")   # float 64bits

You can directly change the precision used in compiled operators with the functions compile4float32 and compile4float64 which respectively enable float 32bits precision (default) and float 64bits (or double) precision.

• other compile options (including boolean value to enable verbosity or to add debugging flag), see ?compile_options

Choosing CPU or GPU computing at runtime

By default, RKeOps runs computations on CPU (even for GPU-compiled operators). To enable GPU computing, you can run (before calling your operator):

use_gpu()
# see `?runtime_options` for a more advanced use of GPU inside RKeOps

You can also specify the GPU id that you want to use, e.g. use_gpu(device=0) to use GPU 0 (default) for instance.

To deactivate GPU computations, you can run use_cpu().

In CPU mode, you can control the number of CPU cores used by RKeOps for computations, e.g. with use_cpu(ncore = 2) to run on 2 cores.

Other runtime options

• device_id: choose on which GPU the computations will be done, default is 0.

set_rkeops_option("device_id", 0)

Note: We recommend to handle GPU assignation outside RKeOps, for instance by setting the environment variable CUDA_VISIBLE_DEVICES. Thus, you can keep the default GPU device id = 0 in RKeOps.

• Other runtime options, see ?runtime_options

Precision

By default, RKeOps uses float 32bits precision for computations. Since R only considers 64bits floating point numbers, if you want to use float 32bits, input data and output results will be casted befors and after computations respectively in your RKeOps operator. If your application requires to use float 64bits (double) precision, keep in mind that you will suffer a performance loss (potentially not an issue on high-end GPUs). In any case, compensated summation reduction is available in KeOps to correct for 32bits floating point arithmetic errors.

Data storage orientation

In R, matrices are stored using a column-major order, meaning that a \(M \times D\) matrix is stored in memory as a succession of \(D\) vectors of length \(M\) representing each of its columns. A consequence is that two successive entries of a column are contiguous in memory, but two successive entries of a row are separated by \(M\) elements. See this page for more details.

For RKeOps to be computationnally efficient, it is important that elements of the input matrices are contiguous along the inner dimensions \(D\) (or \(D'\)). Thus, it is recommended to use input matrices where the outer dimension (i.e. indexes \(i\) or \(j\)) are the columns, and inner dimensions the rows, e.g. transpose matrices \(\mathbf X^{t} = [x_{ki}]_{D \times M}\) or \(\mathbf Y^{t} = [y_{k'i}]_{D' \times N}\).

Important: In machine learning and statistics, we generally use data matrices where each sample/observation/individual is a row, i.e. matrices where the outer dimensions correspond to rows, e.g. \(\mathbf X = [x_{ik}]_{M \times D}\), \(\mathbf Y = [y_{ik'}]_{N \times D'}\).

This is the default using case of RKeOps. RKeOps will then automatically convert your matrices to their transpose, where the outer dimensions correspond to columns.

If you want to use data where the inner dimension directly corresponds to rows of your matrices, i.e. \(\mathbf X^{t} = [x_{ki}]_{D \times M}\) or \(\mathbf Y^{t} = [y_{k'i}]_{D' \times N}\), you just need to specify the input parameter inner_dim=0 when calling your operator.

Example:

# standard column reduction of a matrix product
op <- keops_kernel(formula = "Sum_Reduction((x|y), 1)",
                   args = c("x=Vi(3)", "y=Vj(3)"))

# data (inner dimension = columns)
nx <- 10
ny <- 15
# x_i = rows of the matrix X
X <- matrix(runif(nx*3), nrow=nx, ncol=3)
# y_j = rows of the matrix Y
Y <- matrix(runif(ny*3), nrow=ny, ncol=3)
# computing the result (here, by default `inner_dim=1` and columns corresponds
# to the inner dimension)
res <- op(list(X,Y))

# data (inner dimension = rows)
nx <- 10
ny <- 15
# x_i = columns of the matrix X
X <- matrix(runif(nx*3), nrow=3, ncol=nx)
# y_j = columns of the matrix Y
Y <- matrix(runif(ny*3), nrow=3, ncol=ny)
# computing the result (we specify `inner_dim=0` to indicate that rows
# corresponds to the inner dimension)
res <- op(list(X,Y), inner_dim=0)

Compilation files and cleaning

The compilation of new operators produces shared library (or share object .so) files stored in a build sub-directory of the package installation directory, to be reused and avoid recompilation of already defined operators.

You can check where your compiled operators are stored by running get_build_dir(). To clean RKeOps install and remove all shared library files, you can run clean_rkeops().