For R Users¶

torch-openreml can be used in R through the reticulate package. We recommend installing torch-openreml inside a dedicated conda environment.

Installation¶

1. Install Conda¶

You can skip this step if conda is already installed on your system.

if (is.null(reticulate:::find_conda()[[1]])) {
    reticulate::install_miniconda()
}

2. Create a Conda Environment¶

if (reticulate::condaenv_exists("torch-openreml")) {
    reticulate::conda_remove("torch-openreml")
}

reticulate::conda_create("torch-openreml",
                         python_version = "3.12.12")

3. Install torch-openreml¶

reticulate::conda_install("torch-openreml",
                          pip = TRUE,
                          packages = c("torch-openreml")))

Usage¶

The following example mirrors the Getting Started example. In R, it is often more convenient to construct design matrices directly using model.matrix().

library(reticulate)
use_condaenv("torch-openreml")

openreml <- import("torch_openreml", convert = FALSE)
torch <- import("torch", convert = FALSE)

ScalarMatrix <- openreml$covariance$ScalarMatrix
IdentityMatrix <- openreml$covariance$IdentityMatrix
KroneckerProduct <- openreml$covariance$KroneckerProduct
CovariancePropagation <- openreml$covariance$CovariancePropagation
Sum <- openreml$covariance$Sum
REML <- openreml$REML

data <- agridat::john.alpha

y <- torch$tensor(data$yield, dtype = torch$float32)
X <- model.matrix(~ rep, data = data) |>
    torch$tensor(dtype = torch$float32)

Z_gen <- model.matrix(~ gen - 1, data = data) |>
    torch$tensor(dtype = torch$float32)
Z_rep_block <- model.matrix(~ interaction(rep, block, lex.order = TRUE) - 1, data = data) |>
    torch$tensor(dtype = torch$float32)

G_gen <- ScalarMatrix(length(unique(data$gen)))
G_rep_block <- KroneckerProduct(IdentityMatrix(length(unique(data$rep))),
                                ScalarMatrix(length(unique(data$block))))
V <- Sum(CovariancePropagation(Z_gen, G_gen),
         CovariancePropagation(Z_rep_block, G_rep_block),
         ScalarMatrix(nrow(data)))

fit_openreml <- REML(V)
result <- fit_openreml$optimize(y, X, torch$zeros(3L), verbose = 2L)

print(py_to_r(fit_openreml$get_theta()$numpy()))
print(py_to_r(V$build_params(fit_openreml$get_theta())$numpy()))
print(py_to_r(fit_openreml$get_beta()$numpy()))

∥∇∥:      41.8197, ∥Δ∥: 8.9438, η: 1.00, ∥Δᶜ∥: 8.9438, log 𝓛: -34.3129
∥∇∥:  315227.5000, ∥Δ∥: 0.8838, η: 1.00, ∥Δᶜ∥: 0.8838, log 𝓛: -201185.4531 (-201151.1402)
∥∇∥:  119231.8828, ∥Δ∥: 0.8520, η: 1.00, ∥Δᶜ∥: 0.8520, log 𝓛: -71845.0547 (+129340.3984)
∥∇∥:   44274.7461, ∥Δ∥: 0.8143, η: 1.00, ∥Δᶜ∥: 0.8143, log 𝓛: -26065.9766 (+45779.0781)
∥∇∥:   16330.8594, ∥Δ∥: 0.7587, η: 1.00, ∥Δᶜ∥: 0.7587, log 𝓛: -9457.3301 (+16608.6465)
∥∇∥:    6001.5029, ∥Δ∥: 0.7142, η: 1.00, ∥Δᶜ∥: 0.7142, log 𝓛: -3392.9321 (+6064.3979)
∥∇∥:    2194.9402, ∥Δ∥: 0.7031, η: 1.00, ∥Δᶜ∥: 0.7031, log 𝓛: -1181.3041 (+2211.6281)
∥∇∥:     793.7997, ∥Δ∥: 0.6952, η: 1.00, ∥Δᶜ∥: 0.6952, log 𝓛: -382.4992 (+798.8049)
∥∇∥:     278.6937, ∥Δ∥: 0.6715, η: 1.00, ∥Δᶜ∥: 0.6715, log 𝓛: -102.3017 (+280.1975)
∥∇∥:      90.4135, ∥Δ∥: 0.5957, η: 1.00, ∥Δᶜ∥: 0.5957, log 𝓛: -11.4286 (+90.8731)
∥∇∥:      23.9223, ∥Δ∥: 0.4006, η: 1.00, ∥Δᶜ∥: 0.4006, log 𝓛:  12.6798 (+24.1084)
∥∇∥:       3.8675, ∥Δ∥: 0.1392, η: 1.00, ∥Δᶜ∥: 0.1392, log 𝓛:  16.5910 (+3.9112)
∥∇∥:       0.2448, ∥Δ∥: 0.0157, η: 1.00, ∥Δᶜ∥: 0.0157, log 𝓛:  16.8083 (+0.2173)
∥∇∥:       0.0103, ∥Δ∥: 0.0007, η: 1.00, ∥Δᶜ∥: 0.0007, log 𝓛:  16.8099 (+0.0017)
∥∇∥:       0.0004, ∥Δ∥: 0.0000, η: 1.00, ∥Δᶜ∥: 0.0000, log 𝓛:  16.8096 (-0.0004)
∥∇∥:       0.0000, ∥Δ∥: 0.0000, η: 1.00, ∥Δᶜ∥: 0.0000, log 𝓛:  16.8103 (+0.0007)
∥∇∥:       0.0000, ∥Δ∥: 0.0000, η: 1.00, ∥Δᶜ∥: 0.0000, log 𝓛:  16.8101 (-0.0002)
∥∇∥:       0.0000, ∥Δ∥: 0.0000, η: 1.00, ∥Δᶜ∥: 0.0000, log 𝓛:  16.8100 (-0.0001)

[∇: score, Δ: 𝐉⁻¹∇, η: learning rate, Δᶜ: clip(𝛉 + ηΔ, lb, ub) - 𝛉, 𝓛: restricted likelihood]

✓ Converged at iteration 18
[1] -0.9727983 -1.3280727 -1.2528644
[1] 0.14290193 0.07021837 0.08161610
[1]  4.5182509  0.2978451 -0.4140463