Skip to contents

Matrix Exponential Solving in rxode2

Source: vignettes/articles/rxode2-ind-lin.Rmd
rxode2-ind-lin.Rmd

Background

A linear system of ODEs

dy/dt = A . y

has a closed-form solution given by the matrix exponential,

y(t) = exp(A * t) . y0

which rxode2 evaluates through an eigenvalue decomposition,

y(t) = V . diag(exp(lambda * t)) . V^-1 . y0

where V are the eigenvectors and lambda the eigenvalues of A. For a system that really is linear this is exact – there is no discretization error at any step size, because there is no stepping.

Every standard compartmental PK model is of this form, which makes the matrix exponential a natural fit: the rate matrix is built once from the model’s rate constants and the solution is propagated across each interval in one operation.

rxSolve(..., method = "indLin") selects this integrator.

Three things sharing the indLin name

Name What it is What it does
indLin(model) An R function Rewrites a d/dt() model into matrix-exponential (matExp()) form and returns the model code as text
indLin(state) <- expr A model statement Adds an inhomogeneous forcing term to one compartment of a matExp() model
method = "indLin" An rxSolve() option Selects the matrix-exponential integrator instead of a time-stepping solver

What indLin() does

indLin() is the symbolic half. Give it a model written with d/dt() and it returns the same model written as a rate matrix:

library(rxode2)

oralOde <- function() {
  ini({
    ka <- 1
    cl <- 0.1
    v <- 1
  })
  model({
    d/dt(depot) <- -ka * depot
    d/dt(central) <- ka * depot - (cl / v) * central
    cp <- central / v
  })
}

cat(indLin(oralOde))

which prints

matExp()
cmt(depot)
cmt(central)
k_depot_central = ka
k_central_output = cl/v
cp=central/v

indLin() differentiated each right-hand side with respect to each state to get

A = [ -ka       0    ]
    [  ka   -cl/v    ]

and then reported A through its micro-constants rather than as a matrix: the off-diagonal A[central, depot] = ka becomes k_depot_central, and the column sum lost to the outside world, -(A[depot, depot] + A[central, depot]), becomes k_central_output. The diagonal is implied, so it never has to be written down. Nothing else changes – the output line cp = central/v is carried through untouched.

A nonlinear model converts too. The linear part goes into the rate matrix and the nonlinear remainder becomes an indLin() forcing:

mmOde <- function() {
  ini({
    ka <- 1
    km <- 0.5
    vmax <- 0.2
    v <- 1
  })
  model({
    d/dt(depot) <- -ka * depot
    cp <- central / v
    d/dt(central) <- ka * depot - vmax * cp / (km + cp)
  })
}

cat(indLin(mmOde))
matExp()
cmt(depot)
cmt(central)
k_depot_central = ka
indLin(central) <- -vmax*central/(v*(km + central/v))
cp=central/v

Absorption stays in the rate matrix, because ka is a genuine constant. Michaelis-Menten elimination does not: its rate depends on how much drug is in central. A rate constant that reads a compartment would make the rate matrix change with the solution, which is exactly what the matrix exponential assumes does not happen – so rxode2 rejects one:

rxode2("
matExp()
cmt(central)
vmax <- 10
km <- 5
k_central_output = vmax/(km + central)
")
#> Error: matrix exponential rate constant 'k_central_output' depends on the
#> compartment 'central'; rate constants must be constant in the states --
#> put the state-dependent part in 'indLin(central) <- ...' instead

The forcing is where a state-dependent term belongs, and it is what the solver iterates on.

You rarely have to call indLin() yourself: rxSolve(..., method = "indLin") on an ODE model runs this conversion internally. Calling it directly is how you inspect the factorization, or obtain a matExp() model you want to edit and keep.

How rxode2 solves it

Over each substep the rate matrix is evaluated once at the current state, the matrix exponential advances the states exactly for that frozen matrix, and the next substep starts from the updated state.

For a linear model this is exact and the substep length is irrelevant:

oralMe <- function() {
  ini({
    ka <- 1
    cl <- 0.1
    v <- 1
  })
  model({
    matExp()
    cmt(depot)
    cmt(central)
    k_depot_central <- ka
    k_central_output <- cl / v
    cp <- central / v
  })
}

ev <- et(amt = 100, ii = 12, until = 48) |> et(seq(0, 72, by = 0.5))

meSol <- rxSolve(oralMe, ev, method = "indLin")
odeSol <- rxSolve(oralOde, ev, method = "liblsoda")

## differences are at machine precision (~3e-13)
max(abs(meSol$cp - odeSol$cp))

For a model with a state-dependent indLin() forcing, the rate matrix is still constant but the forcing is not, so the substep is only an approximation. The solver therefore iterates the forcing to a fixed point within each substep and chooses the substep itself from a local error estimate – so atol and rtol control the accuracy, the same way they do for every other adaptive method:

solveMm <- function(...) {
  rxSolve(mmOde, et(amt = 3) |> et(c(0.1, 0.25, 0.5, 0.75, 1, 2, 4, 6, 8, 12, 16, 24, 30)), ...)
}

refMm <- solveMm(method = "liblsoda", atol = 1e-12, rtol = 1e-12)

maxRelErr <- function(sol) max(abs(sol$cp - refMm$cp) / refMm$cp)

maxRelErr(solveMm(method = "indLin"))
maxRelErr(solveMm(method = "indLin", atol = 1e-10, rtol = 1e-10))
atol/rtol Maximum relative error vs LSODA Steps ($counts$slvr)
1e-4 4.0e-03 148
1e-6 7.3e-05 1307
default (1e-8/1e-6) 7.6e-06 3212
1e-8 6.8e-07 12865
1e-10 6.4e-09 128442

The step is sized from a first-order error estimate, but each step advances on the average of two answers whose leading errors are equal and opposite – the one built at the step’s starting state and the converged one – so what is propagated is second order. The practical consequence is the scaling: the error falls roughly in proportion to the tolerance, not to its square root, while the step count still grows only as its square root. So the work needed for a given accuracy grows like 1/sqrt(error) rather than 1/error, which is worth a factor of tens to hundreds of steps at the accuracies people usually want.

For a tight tolerance, indLinRichardson = TRUE is worth trying. It runs each step once whole and twice at half length and extrapolates, which raises the step from second to third order at three times the cost per step. That is a losing trade at loose tolerances and a winning one at tight ones – on this model it is a wash at a relative error of 1e-4, about three times faster at 1e-6, and about six times faster at 1e-8:

maxRelErr(solveMm(method = "indLin", atol = 1e-8, rtol = 1e-8))
maxRelErr(solveMm(method = "indLin", atol = 1e-8, rtol = 1e-8,
                  indLinRichardson = TRUE))

hmax is now only an upper bound on the substep. Lowering it does not meaningfully improve the answer, because the solver is already choosing a shorter step than hmax allows:

maxRelErr(solveMm(method = "indLin", hmax = 1))    # 8.4e-06
maxRelErr(solveMm(method = "indLin", hmax = 0.01)) # 1.9e-06

This changed in rxode2 5.1.7. Previously the nonlinearity was folded into a rate constant, there was no iteration and no error control, and accuracy was governed entirely by hmax – the default gave a 71% error on the sampling grid above. If you set an explicit hmax to work around that, you can drop it and set atol/rtol instead.

Matrix exponential algorithms

Three algorithms are available via indLinMatExpType:

Value Algorithm Notes
"expokit" Sidje (1998) Krylov/Pade Default. Good for sparse or large matrices; controlled by indLinPhiTol and indLinPhiM
"Al-Mohy" Al-Mohy and Higham (2009) scaling and squaring Accurate; order controlled by indLinMatExpOrder
"arma" RcppArmadillo expmat Fastest for small dense matrices 
almohySol <- rxSolve(mmOde, ev, method = "indLin", hmax = 0.01,
                     indLinMatExpType = "Al-Mohy", indLinMatExpOrder = 8L)

armaSol <- rxSolve(mmOde, ev, method = "indLin", hmax = 0.01,
                   indLinMatExpType = "arma")

For the expokit method, two additional controls are available:

  • indLinPhiTol (default 1e-7): tolerance for the exponential-matrix action computation.
  • indLinPhiM (default 0L, meaning automatic): maximum Krylov basis size.

Multi-state terms

When an ODE term is the product of two or more states (y * z), no way of factoring it produces a legal rate constant: divide out y and the coefficient still contains z, and vice versa. Such a term therefore goes to the forcing whole. The van der Pol oscillator, whose mu * (1 - y^2) * dy term mixes both states, shows what that looks like.

vanPol <- function() {
  ini({
    mu <- 1
  })
  model({
    y(0) <- 2
    dy(0) <- 0
    d/dt(y) <- dy
    d/dt(dy) <- mu * (1 - y^2) * dy - y
  })
}

cat(indLin(vanPol))
matExp()
cmt(y)
cmt(dy)
k_y_dy = -1
k_y_output = 1
k_dy_y = 1
k_dy_output = -(1 + mu)
indLin(dy) <- -mu*dy*Rx_pow_di(y, 2)
y(0)=2
dy(0)=0

Expanding mu * (1 - y^2) * dy gives mu*dy - mu*y^2*dy. The first part is linear in dy and lands in the rate matrix; only the second, which mixes y and dy, becomes a forcing. The rate matrix stays constant in the states, which is what makes the matrix exponential valid, and the solver iterates the forcing.

rxIndLinStrategy() and rxIndLinState() used to choose which state a multi-state product was factored onto. They no longer affect the conversion – no choice of state yields a state-free coefficient – and are kept only so existing code keeps running.

Writing the rate matrix directly with matExp()

Everything above went through symengine to factor the ODEs symbolically. For complex expressions that step can fail or produce an unexpected factorization. matExp() skips it entirely: you write the rate constants, and rxode2 builds the rate matrix from them at compile time. This is:

  • The most stable path – no CAS factorization, no symengine dependency, no expression simplification.
  • Familiar to NONMEM users – it mirrors ADVAN5/7 K(i,j) notation with the same source-first index convention.
  • Self-documenting – every rate constant name says what it connects.

matExp() model syntax

Statement Meaning
matExp() Declare this model uses matrix exponential integration. Required. Cannot be combined with d/dt().
k_from_to <- expr First-order rate from compartment from to compartment to.
k.from.to <- expr Dot-separator alternative (identical meaning).
k_cmt_output <- expr Elimination from cmt to outside the system (like NONMEM’s K(i,0)).
k_from_to_nd <- expr Non-depleting transfer: only to gains; from is unchanged.
indLin(state) <- expr Add an inhomogeneous forcing term to state.
cmt(name) Optionally declare a compartment explicitly.
other <- expr Standard derived-variable assignment (output calculation).

Compartments are auto-detected from the k_from_to names, and the diagonal entries of the rate matrix are computed automatically from the outflow terms.

Declaring compartment order with cmt()

Without cmt() declarations, compartment order in the rate matrix – and therefore the CMT numbers used in the event table – follows the order in which names first appear in the k_from_to assignments. Use cmt() before the rate constants to pin the numbering, for example to match a NONMEM control stream where the depot is CMT=1:

mmComp <- function() {
  ini({
    ka <- 1
    km <- 0.5
    vmax <- 0.2
    v <- 1
  })
  model({
    matExp()
    cmt(depot)   ## CMT = 1
    cmt(central) ## CMT = 2
    cp <- central / v
    k_depot_central <- ka ## absorption: depot -> central
    ## MM elimination is state dependent, so it is a forcing, not a rate
    ## constant -- writing it as k_central_output would be an error.
    indLin(central) <- -vmax * cp / (km + cp)
  })
}

Compartments that appear in rate constants but have no cmt() declaration are appended after the declared ones in first-appearance order, so you can fix the positions that matter and let the rest be auto-ordered. cmt() declarations are purely for ordering; they do not change the model mathematics.

Rules and restrictions 

## ERROR: matExp() cannot be combined with d/dt()
rxode2({
  matExp()
  d/dt(depot) <- -ka * depot
})
#> :ERR: Matrix exponential models cannot be used with any ODEs

## ERROR: self-transfer is not allowed
rxode2({
  matExp()
  k_depot_depot <- 0.1
})
#> :ERR: transfer from a compartment to itself (e.g., k_cmt1_cmt1) is not allowed

## ERROR: indLin() without matExp()
rxode2({
  cmt(depot)
  indLin(depot) <- 0.5
})
#> :ERR: indLin() cannot be used without matExp() defined in the model

## ERROR: a rate constant that reads a compartment
rxode2({
  matExp()
  cmt(central)
  k_central_output <- vmax / (km + central)
})
#> :ERR: matrix exponential rate constant 'k_central_output' depends on the
#>       compartment 'central'; rate constants must be constant in the states --
#>       put the state-dependent part in 'indLin(central) <- ...' instead

State-dependent rate constants

Rate expressions can reference derived variables that depend on the current state, as in the Michaelis-Menten example above. The rate matrix is then no longer constant, and the substep argument from the accuracy section applies – set hmax:

library(ggplot2)

ev <- et(amt = 3) |> et(seq(0, 30, by = 0.25))

mmSol <- rxSolve(mmComp, ev, method = "indLin", hmax = 0.01)

ggplot(as.data.frame(mmSol), aes(time, cp)) +
  geom_line() +
  labs(x = "Time (h)", y = "Concentration (mg/L)",
       title = "Michaelis-Menten PK solved with matExp()")

rxode2 infers depot and central as compartments from the rate constant names. The output label is reserved for the elimination sink, analogous to compartment 0 in NONMEM.

Two-compartment PK with a peripheral compartment 

twoCmt <- function() {
  ini({
    ka <- 1
    cl <- 2
    vc <- 10
    q <- 1
    vp <- 20
  })
  model({
    matExp()
    k_depot_central <- ka   ## absorption
    k_central_output <- cl / vc ## central elimination
    k_central_periph <- q / vc  ## distribution to peripheral
    k_periph_central <- q / vp  ## redistribution from peripheral
    cp <- central / vc
  })
}

twoCmtSol <- rxSolve(twoCmt, et(amt = 100) |> et(seq(0, 48, by = 0.25)),
                     method = "indLin")

ggplot(as.data.frame(twoCmtSol), aes(time, cp)) +
  geom_line() +
  labs(x = "Time (h)", y = "Concentration (mg/L)")

Adding a forcing term with indLin()

indLin(state) <- expr supplies an inhomogeneous term: one that enters a compartment without being proportional to any state. This covers endogenous production, constant inputs that are not in the dosing event table, and steady-state baselines. A turnover response model is the smallest useful example – production kin is a forcing term, loss kout * resp is a rate constant:

turnover <- function() {
  ini({
    kin <- 5
    kout <- 0.1
  })
  model({
    matExp()
    cmt(resp)
    k_resp_output <- kout   ## first-order loss
    indLin(resp) <- kin     ## zero-order production
  })
}

turnoverSol <- rxSolve(turnover, et(seq(0, 100, by = 1)), method = "indLin")

## baseline approaches kin/kout = 50
ggplot(as.data.frame(turnoverSol), aes(time, resp)) +
  geom_line() +
  geom_hline(yintercept = 5 / 0.1, linetype = 2) +
  labs(x = "Time (h)", y = "Response")

The forcing is treated as constant across a substep, which makes it exact for a genuinely constant term like kin above.

Comparison with NONMEM ADVAN5/7

NONMEM’s general linear solver (ADVAN5 for non-steady-state, ADVAN7 for steady-state) specifies rate constants as K(i,j) scalars in $PK, where K(i,j) is the first-order transfer rate from compartment i to compartment j (source first, destination second). rxode2’s matExp() uses the same source-first convention: k_from_to.

NONMEM rxode2 matExp() Meaning
K12 = ka k_cmt1_cmt2 <- ka depot (1) to central (2)
K20 = ke k_cmt2_output <- ke central (2) elimination
K23 = Q/V2 k_cmt2_cmt3 <- q / v2 central to peripheral
K32 = Q/V3 k_cmt3_cmt2 <- q / v3 peripheral to central

Differences from NONMEM:

  • Compartment naming: NONMEM numbers compartments (1, 2, 3, …); rxode2 uses descriptive names (depot, central, periph).
  • Compartment declarations: NONMEM requires explicit COMP declarations; rxode2 auto-detects from rate constant names.
  • Elimination sink: NONMEM uses K(i,0) (compartment 0 = outside); rxode2 uses k_cmt_output (output is the reserved sink label).
  • Diagonal: both compute the diagonal (self-decay) entries automatically from the listed outflow rates.
  • Non-depleting inputs: rxode2 supports k_from_to_nd for zero-depletion transfers; NONMEM has no direct equivalent.
  • Forcing terms: rxode2’s indLin(state) <- expr adds an inhomogeneous term; in NONMEM this requires explicit ODEs (ADVAN6+).

When to use matExp() vs d/dt()

Situation Recommendation
Linear compartmental PK/PD with known rate constants matExp() – most stable, no symengine
Model migrated from NONMEM ADVAN5/7 matExp() – direct K(i,j) translation
Exploratory model with changing structure d/dt() – faster iteration
Nonlinear model d/dt() with a time-stepping solver, unless you have a specific reason not to

Use in nlmixr2

matExp() and indLin() models estimate in nlmixr2 across the focei, nlm and SAEM families and match the equivalent d/dt() model. Be aware of how that currently works: nlmixr2est materializes the rate constants back into ordinary ODEs before building each estimation model, so a fit solves the ODE form rather than propagating a matrix exponential. The syntax is therefore a modeling convenience there, not a change of integrator.

Full matrix exponential support during estimation is upcoming.

When to use method = "indLin"

Scenario Recommendation
Purely linear system method = "indLin" gives an exact, closed-form solution
Multi-subject population simulation, linear PK method = "liblsoda" (default) is typically fastest overall
State-dependent rate matrix method = "indLin" only with an explicit, verified hmax
Heavily nonlinear or chaotic ODE Time-stepping solvers (liblsoda, dop853)
Stiff system liblsoda or lsoda; the eigenvalue decomposition can struggle when eigenvalues span many orders of magnitude

Advantages

  1. Exact for linear systems: no discretization error at all, at any step size.
  2. Interval-based: the system is propagated across an interval rather than stepped through it, so no local error control is needed inside the interval.

Limitations

  1. The eigenvalue decomposition does not exist for all square matrices (defective matrices).
  2. For large systems the matrix operations become expensive; the decomposition scales as O(n^3).
  3. Accuracy on a state-dependent model depends entirely on hmax, which has no automatic error control behind it.
  4. The implementation is under active development; results may change between versions.

References

  • Moler C, Van Loan C (2003). Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later. SIAM Rev 45(1):3–49. https://doi.org/10.1137/S00361445024180

  • Al-Mohy AH, Higham NJ (2009). A new scaling and squaring algorithm for the matrix exponential. SIAM J Matrix Anal Appl 31(3):970–989.

  • Sidje RB (1998). Expokit: a software package for computing matrix exponentials. ACM Trans Math Softw 24(1):130–156.

  • Beal S, Sheiner LB, Boeckmann A, Bauer RJ (2009). NONMEM Users Guides (1989–2009). Icon Development Solutions, Ellicott City, MD. (ADVAN5/7 documentation for K(i,j) rate-constant specification.)