cdeterman / gpuR

R interface to use GPU's
241 stars 26 forks source link

cholesky decomposition and chol2inv #143

Open vmoprojs opened 5 years ago

vmoprojs commented 5 years ago

Hi,

I would like to use some gpuR features in likelihood functions optimizations. To this end, I have two particular questions:

rm(list=ls())
library(gpuR)
library(GeoModels)
################################################################
### Example 1. Spatial covariance matrix associated to
### a Matern correlation model
###############################################################
# "GeoModels" package may be installed from github as:
# install_github("vmoprojs/GeoModels")
# library(GeoModels)

# Define the spatial-coordinates of the points:

NN=5000
x <- runif(NN, 0, 1)
y <- runif(NN, 0, 1)
coords <- cbind(x,y)

# Correlation Parameters for Matern model 
CorrParam("Matern")

# Matern Parameters 
param=list(smooth=0.5,sill=1,scale=0.2,nugget=0)

cc <- GeoCovmatrix(coordx=coords, corrmodel="Matern", param=param)$covmatrix

gpd1=vclMatrix(cc,type = "double")

set.seed(5)
vv=rnorm(NN)

vv1=vclVector(vv,type="double")

######## Chol Decom %*% Vector
system.time(chol(cc)%*%vv)
# user  system elapsed 
# 17.968   0.174  18.270
system.time((chol(gpd1)%*%vv1)[])
# user   system  elapsed 
# 1200.162    9.631  205.755 

######## Inverse
system.time(solve(cc))
# user   system  elapsed 
# 1147.829    7.406  241.336 
system.time(solve(gpd1)[])
# user   system  elapsed 
# 1514.267   29.650  244.032

### Quadratic form:
system.time(vv%*%solve(cc)%*%vv)
# user  system elapsed 
# 135.174   0.848 136.465
system.time((vv1%*%solve(gpd1)%*%vv1)[])
# user   system  elapsed 
# 1553.757   31.814  251.499 

Best,

VMO

cdeterman commented 5 years ago

@vmoprojs thanks for the report. There currently is no chol2inv function. I will need to add that to the list of functions to implement.

Regarding performance, cholesky decomposition is experimental. As such, the performance is not were I would like it to be. I will try to work on improving it but I can't promise a timeline. Regarding solve, that is being addressed in #127.

A couple of quick checks would be to see what GPU hardware you have and to make sure you are using the GPU. Can you provide the output of listContexts()? By default you will be using the first device in the returned list.