Defines a new KeOps operator to compute a specific symbolic formula

This function is the core of the KeOps library, it allows you to create new operators based on kernel operation and matrix reduction described as a mathematic formula.

Usage

keops_kernel(formula, args, sum_scheme = "auto")

Arguments

formula

text string, an operator formula (see Details).

args

vector of text string, formula arguments (see Details).

sum_scheme

character string, method used to sum up results for reductions. This option may be changed only with reductions Sum, MaxSumShiftExp, LogSumExp, Max_SumShiftExpWeight, LogSumExpWeight, SumSoftMaxWeight. Default value "auto" will set this option to "block_sum" for these reductions. Possible values are:

"direct_sum": direct summation.
"block_sum": use an intermediate accumulator in each block before accumulating in the output. This improves accuracy for large sized data.
"kahan_scheme": use Kahan summation algorithm to compensate for round-off errors. This improves accuracy for large sized data.

Value

a function that can be used to compute the value of the symbolic formula on actual data. This function takes as input a list of data corresponding to the formula arguments and return the computed values (generally a vector or a matrix depending on the reduction). It has an additional character input parameter inner_dim indicating if the inner dimension (c.f. browseVignettes("rkeops")) corresponds to columns, i.e. inner_dim="col" (default), or rows, i.e. inner_dim="row", in the data.

Details

The use of the function keops_kernel is detailled in the vignettes, especially how to write formulae, specified input arguments, how to format data to apply the created operators, etc. Run browseVignettes("rkeops") to access the vignettes.

KeOps operators are defined thanks to formula, i.e. a text string describing the mathematical operations that you want to apply to your data, and a list defining the input arguments of your formula.

The function keops_kernel() compiles and imports a new operator that implements the formula given in input, it returns a function that can be used to compute the result of the formula on actual data.

The returned function expects a list of arguments, as data matrices, whose order corresponds to the order given in args to keops_kernel(). We use a list to avoid useless copies of data.

Important: Formula are assumed to include a reduction (e.g. Sum_Reduction(,<axis>)), with an <axis> parameter indicating on which dimension is done the reduction:

0 means a reduction over i (meaning that the result is a Vj variable).
1 means a reduction over j (meaning that the result is a Vi variable).

Note: Data are input as a list, because list are references and since argument passing is done by copy in R, it is better to copy a list of reference than the actual input data, especially for big matrices.

You should be careful with the input dimension of your data, so that it correspond to the input dimension specified in args (see inner or outer dimension in browseVignettes("rkeops").

It is possible to compute partial derivatives of user defined operators with the function keops_grad().

Author

Ghislain Durif

Examples

if (FALSE) {
set_rkeops_options()

## Example 1
# Defining a function that computes for each j the sum over i
# of the scalar products between `x_i` and `y_j` (both 3d vectors), 
# i.e. the sum over the rows of the result of the matrix product `X * t(Y)`
# where `x_i` and `y_j` are the respective rows of the matrices `X` and `Y`.
op <- keops_kernel(
    formula = "Sum_Reduction((x|y), 1)", args = c("x=Vi(3)", "y=Vj(3)"))
# data
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)
# compute the result (here, by default `inner_dim="col"` and 
# columns corresponds to the inner dimension)
res <- op(list(X,Y))

## Example 1 bis
# In example 1, the inner dimension (i.e. the common dimension of vectors 
# `x_i` and `y_j` corresponds to columns of the matrices `X` and `Y` resp.).
# We know consider the inner dimension to be the rows of the matrices `X` 
# and `Y`.

# data
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)
# compute the result (we specify `inner_dim="row"` to indicate that the rows 
# corresponds to the inner dimension)
res <- op(list(X,Y), inner_dim="row")

## Example 2
# Defining a function that computes the convolution with a Gaussian kernel 
# i.e. the sum over i of `e^(lambda * ||x_i - y_j||^2) * beta_j` where `x_i`, 
# `y_j` and `beta_j` are 3d vectors, and `lambda` is a scalar parameter.
op = keops_kernel(
    formula = "Sum_Reduction(Exp(lambda*SqNorm2(x-y))*beta, 1)",
    args = c("x=Vi(3)", "y=Vj(3)", "beta=Vj(3)", "lambda=Pm(1)"))

# data
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)
# beta_j = rows of the matrix beta
beta <- matrix(runif(ny*3), nrow=ny, ncol=3)
# !! important !! y and beta should have the same dimension

# parameter
lambda <- 0.25

# compute the result
res <- op(list(X, Y, beta, lambda))
}