Matrix exponentiation and more.

Description

Methods for computing the solution of the ODE

w'(t) = x w(t) + u

with initial condition w(0) = v at one or more time points.

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
  ## S4 method for signature 'matrix,vector'
expv(x, v, t = 1, u = NULL,
    Markov = FALSE, transpose = Markov, ...)

  ## S4 method for signature 'CsparseMatrix,vector'
expv(x, v, t = 1,
    u = NULL, Markov = FALSE, transpose = Markov, ...)

  ## S4 method for signature 'dgCMatrix,vector'
expv(x, v, t = 1, u = NULL,
    Markov = FALSE, transpose = Markov, ...)

  ## S4 method for signature 'TsparseMatrix,vector'
expv(x, v, t = 1,
    u = NULL, Markov = FALSE, transpose = Markov, ...)

  ## S4 method for signature 'matrix.csc,vector'
expv(x, v, t = 1, u = NULL,
    Markov = FALSE, transpose = Markov, ...)

  ## S4 method for signature 'matrix.csr,vector'
expv(x, v, t = 1, u = NULL,
    Markov = FALSE, transpose = Markov, ...)

  ## S4 method for signature 'matrix.coo,vector'
expv(x, v, t = 1, u = NULL,
    Markov = FALSE, transpose = Markov, ...)

Arguments

x

a matrix.

v

a numeric vector. The initial value.

t

a numeric vector of time points at which the solution is computed.

u

a numeric vector. Default NULL.

Markov

logical. If TRUE the matrix is taken to be an rate matrix and steps are taken to ensure that the computed result is a probability vector. Default FALSE.

transpose

logical. If TRUE transpose the matrix before the solution is computed. Default equals Markov.

...

other arguments passed to Rexpv.

Details

Analytically the solution is given as

w(t) = exp(tx)v + tphi(tx)u

with phi(z) = (exp(z)-1)/z. For large matrices x the computation of the full matrices exp(tx) and phi(tx) is slow. An alternative is to compute the solution w directly for a given initial value v using Krylov subspace methods. This is, in particular, efficient for large sparse matrices.

Note that if Q is a rate matrix for a homogeneous continuous time Markov process (non-negative off-diagonals and row-sums 0) and v is a probability vector (the initial distribution at time point 0), then the distribution at time point t solves

w'(t) = Q^T w(t).

In this case we want to take x to be Q^T to get the desired solution.

The solution is computed using the Fortran package Expokit. Methods are available for matrix classes implemented in the Matrix package as well as the SparseM package. The implementation avoids the computation of the full matrix exponential of tx and the approach is advantageous when we want to compute w(t) for one or a few initial values v. The full matrix exponential should not be computed this way by looping over n different initial values.

Though there is a method implemented for ordinary (dense) matrices, such a matrix is simply coerced into a CsparseMatrix before the solution is computed. It is recommended that large sparse matrices are stored and handled as such, e.g. using the classes and methods from the Matrix package. Dense intermediates should be avoided.

The x matrix is allowed to be a dense complex matrix in which case v and u are also allowed to be complex.

Value

An n by k matrix with k the length of t and n the dimension of the matrix x. The i'th column contains the solution of the ODE at time point i.

Author(s)

Niels Richard Hansen Niels.R.Hansen@math.ku.dk

See Also

Rexpv, expm, expm

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
### Small 4 by 4 example.
x <- matrix(c(-1, 0, 1, 0,
              0, -2, 0, -3,
              1, 0, 0, 0,
              0, 2, -1, 0),
            4, 4)
v <- c(1, 1, 1, 1)

require(Matrix)
require(SparseM)

w <- cbind(padm(x) %*% v,
           expv(x, v),
           expv(Matrix(x, sparse = TRUE), v),
           expv(as.matrix.coo(x), v),
           expv(as.matrix.csr(x), v),
           expv(as.matrix.csc(x), v)
           )

stopifnot(all.equal(w[, 1], w[, 2]),
          all.equal(w[, 1], w[, 3]),
          all.equal(w[, 1], w[, 4]),
          all.equal(w[, 1], w[, 5]),
          all.equal(w[, 1], w[, 6]))

u <- c(2, 0, 1, 1)
ex <- padm(x)
w <- cbind(ex %*% v + (ex - diag(1, 4)) %*% solve(x, u),
           expv(x, v, u = u),
           expv(Matrix(x, sparse = TRUE), v, u = u),
           expv(as.matrix.coo(x), v, u = u),
           expv(as.matrix.csr(x), v, u = u),
           expv(as.matrix.csc(x), v, u = u)
           )

stopifnot(all.equal(w[, 1], w[, 2]),
          all.equal(w[, 1], w[, 3]),
          all.equal(w[, 1], w[, 4]),
          all.equal(w[, 1], w[, 5]),
          all.equal(w[, 1], w[, 6]))

############################################################
### Linear birth-death Markov process with immigration
############################################################

alpha <- 2.1  ## Death rate per individual
beta <- 2     ## Birth rate per individual
delta <- 20   ## Immigration rate

n <- 500L     ## state space {0, ..., n-1}
i <- seq(1, n)
rates <- c(alpha * i[-1], ## subdiagonal
           -(c(0, alpha * i[-1]) +
             c(delta, beta * i[-c(1,n)] + delta, 0)), ## diagonal
           c(delta, beta * i[-c(1, n)] + delta)) ## superdiagonal
j <- c(i[-n], i, i[-1])
i <- c(i[-1], i, i[-n])

## Sparse rate matrix constructed without dense intermediate
require(Matrix)
Q <- sparseMatrix(i = i, j = j, x = rates, dims = c(n, n))

## Evolution of uniform initial distribution
p0 <- rep(1, n)/n
time <- seq(0, 10, 0.2)
Pt <- expv(Q, p0, t = time, Markov = TRUE)

## Not run: 
matplot(1:n, Pt, type = "l")
image(time, 0:(n-1), -t(Pt), col = terrain.colors(100))

## End(Not run)