dlode45

Deep learning solution of nonstiff ordinary differential equation (ODE)

Since R2021b

Syntax

Y = dlode45(net,tspan,Y0)

Y = dlode45(odefun,tspan,Y0,theta)

Y = dlode45(___,DataFormat=FMT)

Y = dlode45(___,Name=Value)

Description

The neural ordinary differential equation (ODE) operation returns the solution of a specified ODE.

The dlode45 function applies the neural ODE operation to dlarray data. Using dlarray objects makes working with high dimensional data easier by allowing you to label the dimensions. For example, you can label which dimensions correspond to spatial, time, channel, and batch dimensions using the "S", "T", "C", and "B" labels, respectively. For unspecified and other dimensions, use the "U" label. For dlarray object functions that operate over particular dimensions, you can specify the dimension labels by formatting the dlarray object directly, or by using the DataFormat option.

Note

This function applies the neural ODE operation to dlarray data in deep learning models defined as functions or in custom layer functions. If you want to apply the neural ODE operation within a dlnetwork object or Layer array, use neuralODELayer. To solve ODEs for other workflows, use ode45.

Y = dlode45(net,tspan,Y0) integrates the system of ODEs characterized by the dlnetwork object net on the time interval defined by the first and last elements of tspan and with the initial conditions Y0. (since R2023a)

Y = dlode45(odefun,tspan,Y0,theta) integrates the system of ODEs given by the function handle odefun on the time interval defined by the first and last elements of tspan, with the initial conditions Y0 and parameters theta.

Y = dlode45(___,DataFormat=FMT) specifies the data format for the unformatted initial conditions Y0. The format must contain "S" (spatial), "C" (channel), and "B" (batch) dimension labels only.

Y = dlode45(___,Name=Value) specifies additional options using one or more name-value arguments. For example, Y = dlode45(odefun,tspan,Y0,theta,GradientMode="adjoint") integrates the system of ODEs given by odefun and computes gradients by solving the associated adjoint ODE system.

Examples

Apply Neural ODE Operation

For the initial conditions, create a formatted dlarray object containing 128 observations of numeric feature data with 4 channels. Specify the format "CB" (channel, batch).

miniBatchSize = 128;
numChannels = 4;
Y0 = rand(numChannels,miniBatchSize);
Y0 = dlarray(Y0,"CB");

View the size and format of the initial conditions.

size(Y0)

ans = 1×2

     4   128

dims(Y0)

ans = 
'CB'

Create a dlnetwork object that characterizes the ODE function given by two fully connect operations with a tanh activation between them.

layers = [
    featureInputLayer(4)
    fullyConnectedLayer(10)
    tanhLayer
    fullyConnectedLayer(4)];
net = dlnetwork(layers);

Specify an interval of integration of [0 0.1].

tspan = [0 0.1];

Apply the neural ODE operation.

Y = dlode45(net,tspan,Y0);

View the size and format of the output.

size(Y)

ans = 1×2

     4   128

dims(Y)

ans = 
'CB'

Specify Custom Neural ODE Function

For the initial conditions, create a formatted dlarray object containing a batch of 128 28-by-28 images with 64 channels. Specify the format "SSCB" (spatial, spatial, channel, batch).

miniBatchSize = 128;
inputSize = [28 28];
numChannels = 64;
Y0 = rand(inputSize(1),inputSize(2),numChannels,miniBatchSize);
Y0 = dlarray(Y0,"SSCB");

View the size and format of the initial conditions.

size(Y0)

ans = 1×4

    28    28    64   128

dims(Y0)

ans = 
'SSCB'

Specify the ODE function. Define the function odeModel, listed in the ODE Function section of the example, which applies a convolution operation followed by a hyperbolic tangent operation to the input data.

odefun = @odeModel;

Initialize the parameters for the convolution operation in the ODE function. The output size of the ODE function must match the size of the initial conditions, so specify the same number of filters as the number of input channels.

filterSize = [3 3];
numFilters = numChannels;

parameters = struct;
parameters.Weights = dlarray(rand(filterSize(1),filterSize(2),numChannels,numFilters));
parameters.Bias = dlarray(zeros(1,numFilters));

Specify an interval of integration of [0 0.1].

tspan = [0 0.1];

Apply the neural ODE operation.

Y = dlode45(odefun,tspan,Y0,parameters);

View the size and format of the output.

size(Y)

ans = 1×4

    28    28    64   128

dims(Y)

ans = 
'SSCB'

ODE Function

The ODE function odeModel takes as input the function inputs t (unused) and y, and the ODE function parameters p containing the convolution weights and biases, and returns the output of the convolution-tanh block operation. The convolution operation applies padding such that the output size matches the input size.

function z = odeModel(t,y,p)

weights = p.Weights;
bias = p.Bias;
z = dlconv(y,weights,bias,Padding="same");
z = tanh(z);

end

Input Arguments

net — Neural network characterizing neural ODE function

initialized dlnetwork object

Since R2023a

Neural network characterizing neural ODE function, specified as an initialized dlnetwork object.

If net has one input, then predict(net,Y) defines the ODE system. If net has two inputs, then predict(net,T,Y) defines the ODE system, where T is a time step repeated over the batch dimension.

When the first argument is a dlnetwork object and GradientMode is "adjoint", the network State property must be empty. To use a network with a nonempty State property, set GradientMode to "direct".

odefun — Function to solve

function handle

Function to solve, specified as a function handle that defines the function to integrate.

Specify odefun as a function handle with syntax z = fcn(t,y,p), where t is a scalar, y is a dlarray, and p is a set of parameters given by theta. The function returns a dlarray with the same size and format as y. The function must accept all three input arguments t, y, and p, even if not all the arguments are used in the function. The size of the ODE function output z must match the size of the initial conditions.

For example, specify the ODE function that applies a convolution operation followed by a tanh operation.

function z = dlconvtanh(t,y,p)

weights = p.Weights;
bias = p.Bias;
z = dlconv(y,weights,bias,Padding="same");
z = tanh(z);

end

Note here that the t argument is unused.

Data Types: function_handle

tspan — Interval of integration

numeric vector | unformatted dlarray vector

Interval of integration, specified as a numeric vector or an unformatted dlarray vector with two or more elements. The elements in tspan must be all increasing or all decreasing.

The solver imposes the initial conditions given by Y0 at the initial time tspan(1), then integrates the ODE function from tspan(1) to tspan(end).

• If tspan has two elements, [t0 tf], then the solver returns the solution evaluated at point tf.

• If tspan has more than two elements, [t0 t1 t2 ... tf], then the solver returns the solution evaluated at the given points [t1 t2 ... tf]. The solver does not step precisely to each point specified in tspan. Instead, the solver uses its own internal steps to compute the solution, then evaluates the solution at the points specified in tspan. The solutions produced at the specified points are of the same order of accuracy as the solutions computed at each internal step.

  Specifying several intermediate points has little effect on the efficiency of computation, but for large systems it can affect memory management.

Note

The behavior of the dlode45 function differs from the ode45 function.

If InitialStep or MaxStep is [], then the software uses the values of tspan to initialize the values.

• If InitialStep is [], then the software uses the elements of tspan as an indication of the scale of the task. When you specify tspan with different numbers of elements, the solution of the solver can change.

• If MaxStep is [], then the software calculates the maximum step size using the first and last elements of tspan. When you change the initial or final values of tspan, the solution of the solver can change because the solver uses a different step sequence.

Y0 — Initial conditions

dlarray object

Initial conditions, specified as a formatted or unformatted dlarray object.

If Y0 is an unformatted dlarray, then you must specify the format using the DataFormat option.

For neural ODE operations, the data format must contain "S", "C", and "B" dimension labels only. The initial conditions must not have a "T" or "U" dimension.

theta — Parameters of ODE function

dlnetwork object (since R2023a) | dlarray object | structure | cell array | table

Parameters of ODE function, specified as one of these values:

• dlnetwork object (since R2023a)

• dlarray object

• Structure or nested structures of one or more of these values:

  • dlnetwork object (since R2023a)

  • dlarray object

  • Any MATLAB^® data type (since R2023b)

• Cell array of one or more of these values:

  • dlnetwork object (since R2023a)

  • dlarray object

  • Any MATLAB data type (since R2023b)

• Table with these variables:

  • Layer — Layer name

  • Parameter — Parameter name

  • Value — Parameter value, specified as a one of these values:

    • dlarray object

    • Any MATLAB data type (since R2023b)

For automatic differentiation workflows (for example, when you define a custom layer or use the dlgradient and dlfeval functions), the software only computes gradients with respect to the dlarray objects contained in the object specified by theta.

The software passes the parameters to odefun as its third input argument.

Name-Value Arguments

Specify optional pairs of arguments as Name1=Value1,...,NameN=ValueN, where Name is the argument name and Value is the corresponding value. Name-value arguments must appear after other arguments, but the order of the pairs does not matter.

Example: Y = dlode45(odefun,tspan,Y0,theta,GradientMode="adjoint") integrates the system of ODEs given by odefun and computes gradients by solving the associated adjoint ODE system.

DataFormat — Description of data dimensions

character vector | string scalar

Description of the data dimensions, specified as a character vector or string scalar.

A data format is a string of characters, where each character describes the type of the corresponding data dimension.

The characters are:

• "S" — Spatial

• "C" — Channel

• "B" — Batch

• "T" — Time

• "U" — Unspecified

For example, consider an array containing a batch of sequences where the first, second, and third dimensions correspond to channels, observations, and time steps, respectively. You can specify that this array has the format "CBT" (channel, batch, time).

For neural ODE operations, the data format must contain "S", "C", and "B" dimension labels only. The initial conditions must not have a "T" or "U" dimension.

If Y0 is not a formatted dlarray object, then you must specify the DataFormat option.

For more information, see Deep Learning Data Formats.

Example: DataFormat="SSCB"

Data Types: char | string

GradientMode — Method to compute gradients

"direct" (default) | "adjoint"

Method to compute gradients with respect to the initial conditions and parameters when using the dlgradient function, specified as one of these values:

• "direct" — Compute gradients by backpropagating through the operations undertaken by the numerical solver. This option best suits large mini-batch sizes or when tspan contains many values.

• "adjoint" — Compute gradients by solving the associated adjoint ODE system. This option best suits small mini-batch sizes or when tspan contains a small number of values.

• "adjoint-seminorm" (since R2023b) — Compute gradients using adjoint seminorm computations [4]. This option typically uses fewer terms in the approximation error computations and can speed up some neural ODE models.

When the first argument is a dlnetwork object and GradientMode is "adjoint", the network State property must be empty. To use a network with a nonempty State property, set GradientMode to "direct".

The dlaccelerate function does not support accelerating the dlode45 function when the GradientMode option is "direct". To accelerate the code that calls the dlode45 function, set the GradientMode option to "adjoint" or accelerate parts of your code that do not call the dlode45 function with the GradientMode option set to "direct".

Warning

When GradientMode is "adjoint", odefun must support function acceleration. Otherwise, the function can return unexpected results.

When GradientMode is "adjoint", the software traces the ODE function input to determine the computation graph used for automatic differentiation. This tracing process can take some time and can end up recomputing the same trace. By optimizing, caching, and reusing the traces, the software can speed up the gradient computation.

Because of the nature of caching traces, not all functions support acceleration.

The caching process can cache values that you might expect to change or that depend on external factors. You must take care when you accelerate functions that:

• have inputs with random or frequently changing values

• have outputs with frequently changing values

• generate random numbers

• use if statements and while loops with conditions that depend on the values of dlarray objects

• have inputs that are handles or that depend on handles

• Read data from external sources (for example, by using a datastore or a minibatchqueue object)

For more information on deep learning function acceleration, see Deep Learning Function Acceleration for Custom Training Loops.

InitialStepSize — Initial step size

[] (default) | positive scalar

Initial step size, specified as a positive scalar or [].

If InitialStepSize is [], then the function automatically determines the initial step size based on the interval of integration and the output of the ODE function corresponding to the initial conditions.

Data Types: single | double | int8 | int16 | int32 | int64 | uint8 | uint16 | uint32 | uint64

MaxStepSize — Maximum step size

[] (default) | positive scalar

Maximum step size, specified as a positive scalar or [].

If MaxStepSize is [], then the function uses a tenth of the interval of integration size.

Data Types: single | double | int8 | int16 | int32 | int64 | uint8 | uint16 | uint32 | uint64

RelativeTolerance — Relative error tolerance

1e-3 (default) | positive scalar

Relative error tolerance, specified as a positive scalar. The relative tolerance applies to all components of the solution.

Data Types: single | double | int8 | int16 | int32 | int64 | uint8 | uint16 | uint32 | uint64

AbsoluteTolerance — Absolute error tolerance

1e-6 (default) | positive scalar

Absolute error tolerance, specified as a positive scalar. The absolute tolerance applies to all components of the solution.

Data Types: single | double | int8 | int16 | int32 | int64 | uint8 | uint16 | uint32 | uint64

Output Arguments

Y — Solution of neural ODE

dlarray object

Solution of the neural ODE at the times given by tspan(2:end), returned as a dlarray object with the same underlying data type as Y0.

If Y0 is a formatted dlarray and tspan contains exactly two elements, then Y has the same format as Y0. If Y0 is not a formatted dlarray and tspan contains exactly two elements, then Y is an unformatted dlarray with the same dimension order as Y0.

If Y0 is a formatted dlarray and tspan contains more than two elements, then Y has the same format as Y0 with an additional appended "T" (time) dimension. If Y0 is not a formatted dlarray and tspan contains more than two elements, then Y is an unformatted dlarray with the same dimension order as Y0 with an additional appended dimension corresponding to time.

Algorithms

Neural Ordinary Differential Equation

The neural ordinary differential equation (ODE) operation returns the solution of a specified ODE. In particular, given an input, a neural ODE operation outputs the numerical solution of the ODE $$y\prime =f(t,y,\theta )$$ for the time horizon (t₀,t₁) and with the initial condition y(t₀) = y₀, where t and y denote the ODE function inputs and θ is a set of learnable parameters. Typically, the initial condition y₀ is either the network input or the output of another deep learning operation.

To apply the operation, dlode45 uses the ode45 function, which is based on an explicit Runge-Kutta (4,5) formula, the Dormand-Prince pair. It is a single-step solver—in computing y(t_n), it needs only the solution at the immediately preceding time point, y(t_n-1) [2] [3].

Deep Learning Array Formats

Most deep learning networks and functions operate on different dimensions of the input data in different ways.

For example, an LSTM operation iterates over the time dimension of the input data, and a batch normalization operation normalizes over the batch dimension of the input data.

To provide input data with labeled dimensions or input data with additional layout information, you can use data formats.

A data format is a string of characters, where each character describes the type of the corresponding data dimension.

The characters are:

• "S" — Spatial

• "C" — Channel

• "B" — Batch

• "T" — Time

• "U" — Unspecified

For example, consider an array containing a batch of sequences where the first, second, and third dimensions correspond to channels, observations, and time steps, respectively. You can specify that this array has the format "CBT" (channel, batch, time).

To create formatted input data, create a dlarray object and specify the format using the second argument.

To provide additional layout information with unformatted data, specify the format using the DataFormat argument.

For more information, see Deep Learning Data Formats.

References

[1] Chen, Ricky T. Q., Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. “Neural Ordinary Differential Equations.” Preprint, submitted June 19, 2018. https://arxiv.org/abs/1806.07366.

[2] Dormand, J. R., and P. J. Prince. “A Family of Embedded Runge-Kutta Formulae.” Journal of Computational and Applied Mathematics 6, no. 1 (March 1980): 19–26. https://doi.org/10.1016/0771-050X(80)90013-3.

[3] Shampine, Lawrence F., and Mark W. Reichelt. “The MATLAB ODE Suite.” SIAM Journal on Scientific Computing 18, no. 1 (January 1997): 1–22. https://doi.org/10.1137/S1064827594276424.

[4] Kidger, Patrick, Ricky T. Q. Chen, and Terry Lyons. “‘Hey, That’s Not an ODE’: Faster ODE Adjoints via Seminorms.” arXiv, May 10, 2021. https://doi.org/10.48550/arXiv.2009.09457.

Extended Capabilities

GPU Arrays

Accelerate code by running on a graphics processing unit (GPU) using Parallel Computing Toolbox™.

The dlode45 function supports GPU array input with these usage notes and limitations:

• This function runs on the GPU when Y0 has underlying data of type gpuArray.

For more information, see Run MATLAB Functions on a GPU (Parallel Computing Toolbox).

Version History

Introduced in R2021b

R2023b: Apply neural ODE operations using adjoint seminorm computations

To apply neural ODE operations using adjoint seminorm computations, set the GradientMode option of the dlode45 function to "adjoint-seminorm".

R2023b: dlode45 supports more data types as parameters

The dlode45 function now supports more data types as parameters for the ODE function, such as structures, cell arrays, and tables that contain numeric arrays and gridded interpolant objects. Use these parameters to define neural ODE functions with greater flexibility.

R2023a: Specify system of ODEs or learnable parameters as a neural network

Specify a system of ODEs as a neural network by specifying a dlnetwork object as the first input argument net.

Specify learnable parameters theta as one of these values:

• dlnetwork object

• Cell array of dlnetwork and dlarray objects

• Structure of dlnetwork and dlarray objects or nested structures of dlnetwork and dlarray objects