In learning from demonstration, we assume access to an expert policy (e.g., a human teleoperator) that can provide actions for robot states. Behavior cloning treats imitation as a supervised learning problem: learn a mapping \( \pi_\theta: \mathcal{S} \rightarrow \mathcal{A} \) from states to actions that mimics the expert policy \( \pi_E \) on a dataset of expert trajectories.

Let the robot's state space be \( \mathcal{S} \subseteq \mathbb{R}^n \) and the action space \( \mathcal{A} \subseteq \mathbb{R}^m \). A trajectory of length \( T \) consists of states and actions \( \tau = (s_1, a_1, \dots, s_T, a_T) \), where \( a_t = \pi_E(s_t) \) during demonstration.

Behavior cloning typically solves

\[ \min_{\theta} \; \mathbb{E}_{(s,a) \sim \mathcal{D}_E} \left[ \ell\big(\pi_\theta(s), a\big) \right], \]

where \( \mathcal{D}_E \) is the empirical distribution of expert state-action pairs and \( \ell \) is a supervised loss (e.g., squared error for continuous torques or cross-entropy for discrete actions).

However, once the learned policy is deployed on the robot, its small prediction errors perturb the state distribution, causing the robot to visit states that are rare or unseen in the demonstration dataset. This covariate shift leads to compounding errors over time.

flowchart TD
  E["Expert policy pi_E"] --> D["Collect demos: (s_t, a_t)"]
  D --> SUP["Train pi_theta by supervised learning"]
  SUP --> DEP["Deploy pi_theta on robot"]
  DEP --> SD["New state distribution under pi_theta"]
  SD --> GAP["Distribution mismatch vs demos"]
  GAP --> ERR["Compounding control errors over horizon"]
        

We formalize the dynamics with a (possibly stochastic) system model \( P(s_{t+1} \mid s_t, a_t) \), e.g. the discrete-time dynamics of a manipulator under a low-level controller. Let \( d_t^{\pi}(s) \) denote the distribution of state \( s_t \) at time \( t \) when executing policy \( \pi \) from some initial distribution \( d_1(s) \).

Assume that behavior cloning achieves low classification/regression error on the expert distribution:

\[ \epsilon = \mathbb{E}_{s \sim d_E(s)} \Big[ \Pr\big( \pi_\theta(s) \neq \pi_E(s) \big) \Big], \]

where \( d_E \) is the marginal distribution of \( s \) under expert rollouts.

Let the per-step loss under policy \( \pi \) be \( c(s,\pi(s)) \in [0,1] \), measuring, for example, task failure probability or a normalized tracking error. Define the total expected loss over horizon \( T \) as

\[ J(\pi) = \sum_{t=1}^{T} \mathbb{E}_{s_t \sim d_t^{\pi}} \big[ c(s_t, \pi(s_t)) \big]. \]

Suppose that at each visited state, if the policy takes the expert action, the instantaneous cost is at most that of the expert plus a constant \( \Delta \), while taking a non-expert action can incur worst-case cost 1 and move to arbitrary states. Under suitable assumptions, one can show that behavior cloning may suffer \( \mathcal{O}(T^2 \epsilon) \) degradation compared to the expert.

Sketch of the compounding error bound.

1. At each time step, the probability that behavior cloning deviates from the expert is at most \( \epsilon \) under the expert state distribution. After one deviation, the state distribution can drift arbitrarily far from the expert distribution.
2. Let \( p_t \) be the probability that at least one mistake occurs up to time \( t \). A union bound yields \( p_t \leq t\epsilon \).
3. Conditioned on at least one past mistake, the cost may be as high as 1 per time step. Thus

   \[ J(\pi_\theta) \leq J(\pi_E) + \sum_{t=1}^{T} p_t \leq J(\pi_E) + \sum_{t=1}^{T} t\epsilon = J(\pi_E) + \mathcal{O}(T^2 \epsilon). \]

The key pathology is that behavior cloning minimizes one-step prediction error under the expert distribution, but is evaluated under the induced long-horizon distribution of its own actions, which may differ dramatically from the expert's.

Dataset Aggregation (DAgger) is an online imitation learning algorithm designed to mitigate covariate shift by actively querying the expert on states visited by the current policy. Instead of training once on a fixed dataset of expert rollouts, DAgger iteratively:

• rolls out the current policy,
• queries the expert for each visited state,
• aggregates the new state-action pairs into a dataset, and
• re-trains (or updates) the policy on the aggregated dataset.

The basic DAgger algorithm can be written as:

\[ \begin{aligned} &\text{Initialize dataset } \mathcal{D}_1 \text{ with expert rollouts}. \\ &\text{For } i = 1, \dots, N: \\ &\quad \pi_i = \text{TrainPolicy}(\mathcal{D}_i). \\ &\quad \text{Execute } \pi_i \text{ to obtain states } s_t. \\ &\quad \text{Query expert: } a_t^E = \pi_E(s_t). \\ &\quad \mathcal{D}_{i+1} = \mathcal{D}_i \cup \{(s_t, a_t^E)\}_t. \\ &\text{Return a suitable policy from } \{\pi_i\}_{i=1}^{N}. \end{aligned} \]

The crucial difference from behavior cloning is that the training distribution \( \mathcal{D}_i \) gradually becomes closer to the distribution of states encountered by the learned policies themselves.

flowchart TD
  D0["Initial expert dataset D_1"] --> T1["Train policy pi_1 on D_1"]
  T1 --> R1["Rollout pi_1 on robot/simulator"]
  R1 --> Q1["Query expert actions on visited states"]
  Q1 --> A1["Aggregate: D_2 = D_1 U new pairs"]
  A1 --> T2["Train pi_2 on D_2"]
  T2 --> LOOP["Repeat rollout, query, aggregate ..."]
  LOOP --> OUT["Final policy pi_N with improved robustness"]
        

DAgger can be analyzed via online learning and no-regret algorithms. At each iteration \( i \), the learner picks a policy \( \pi_i \) and observes a loss function \( \ell_i(\pi) \) defined by the expected surrogate loss over states visited by \( \pi_i \) and expert actions at those states:

\[ \ell_i(\pi) = \mathbb{E}_{s \sim d^{\pi_i}(s)} \Big[ \ell\big(\pi(s), \pi_E(s)\big) \Big]. \]

If the learner uses a no-regret online algorithm over the policy class \( \Pi \) (e.g., online gradient descent in a convex surrogate), then the average policy over iterations performs nearly as well as the best policy in hindsight on the induced losses.

Guarantee (simplified form). Suppose:

• We run DAgger for \( N \) iterations with horizon \( T \).
• The online learning algorithm has regret \( R_N = o(N) \) with respect to losses \( \ell_i \).

Then there exists some iteration \( i \in \{1, \dots, N\} \) such that the expected cost of \( \pi_i \) satisfies

\[ J(\pi_i) \leq J(\pi_E) + \mathcal{O}(T \epsilon_N), \]

where \( \epsilon_N \) is the average surrogate loss of the best policy in \( \Pi \) on the aggregated data. Notice the \( \mathcal{O}(T\epsilon_N) \) scaling, which is asymptotically better than the \( \mathcal{O}(T^2\epsilon) \) scaling of vanilla behavior cloning.

Proof sketch idea.

1. Interpret each DAgger iteration as an online learning round with loss \( \ell_i(\pi) \) defined on the states visited by the current policy.
2. No-regret implies

   \[ \frac{1}{N} \sum_{i=1}^{N} \ell_i(\pi_i) \leq \inf_{\pi \in \Pi} \frac{1}{N} \sum_{i=1}^{N} \ell_i(\pi) + o(1). \]

3. By relating the surrogate loss \( \ell_i \) to task cost under \( \pi_i \), one derives a bound of the form \( J(\pi_i) - J(\pi_E) \leq C T \ell_i(\pi_i) \) for some constant \( C \).
4. Combining these inequalities over \( i \) and using an averaging argument shows that at least one iterate satisfies the desired bound.

The key insight is that by training on the states induced by its own policies, the learner removes the mismatch between training and test distributions, so error does not compound as severely.

We demonstrate behavior cloning and DAgger in Python for a planar two-link manipulator tracking a desired end-effector position. We assume access to an expert controller (e.g., an inverse kinematics plus PD tracking law) implemented in a function expert_policy(state). We then train a neural network policy using PyTorch on aggregated data.

import numpy as np
import torch
import torch.nn as nn
import torch.optim as optim

# Simple planar arm state: [q1, q2, dq1, dq2]
STATE_DIM = 4
ACTION_DIM = 2   # joint torques tau1, tau2

class PolicyNet(nn.Module):
    def __init__(self):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(STATE_DIM, 64),
            nn.ReLU(),
            nn.Linear(64, 64),
            nn.ReLU(),
            nn.Linear(64, ACTION_DIM)
        )

    def forward(self, x):
        return self.net(x)

def expert_policy(state):
    """
    Placeholder for an expert controller, e.g.,
    a PD controller around a reference trajectory.
    state: np.array of shape (4,)
    returns: np.array of shape (2,)
    """
    # Here we just implement a dummy proportional controller
    q_des = np.array([0.5, -0.5])
    q = state[:2]
    dq = state[2:]
    Kp = 5.0
    Kd = 1.0
    tau = Kp * (q_des - q) - Kd * dq
    return tau

def simulate_dynamics(state, action, dt=0.02):
    """
    Very rough placeholder dynamics integrator.
    In practice, use a robotics simulator (e.g. PyBullet, Mujoco).
    """
    q = state[:2]
    dq = state[2:]
    # Simple double integrator approximation: ddq = action
    ddq = action
    dq_next = dq + ddq * dt
    q_next = q + dq_next * dt
    return np.concatenate([q_next, dq_next], axis=0)

def collect_expert_trajectory(T=200):
    states = []
    actions = []
    state = np.zeros(STATE_DIM, dtype=np.float32)
    for t in range(T):
        a = expert_policy(state)
        states.append(state.copy())
        actions.append(a.astype(np.float32))
        state = simulate_dynamics(state, a)
    return np.array(states), np.array(actions)

def train_bc(policy, dataset, epochs=20, batch_size=64, lr=1e-3):
    opt = optim.Adam(policy.parameters(), lr=lr)
    criterion = nn.MSELoss()
    states, actions = dataset
    states_t = torch.from_numpy(states)
    actions_t = torch.from_numpy(actions)
    N = states.shape[0]
    for ep in range(epochs):
        perm = torch.randperm(N)
        for i in range(0, N, batch_size):
            idx = perm[i:i+batch_size]
            s_batch = states_t[idx]
            a_batch = actions_t[idx]
            opt.zero_grad()
            pred = policy(s_batch)
            loss = criterion(pred, a_batch)
            loss.backward()
            opt.step()

def rollout(policy, T=200):
    states = []
    state = np.zeros(STATE_DIM, dtype=np.float32)
    for t in range(T):
        with torch.no_grad():
            s_t = torch.from_numpy(state).float().unsqueeze(0)
            a_t = policy(s_t).squeeze(0).numpy()
        states.append(state.copy())
        state = simulate_dynamics(state, a_t)
    return np.array(states)

# --- Behavior Cloning (single-shot) ---
expert_states, expert_actions = collect_expert_trajectory(T=200)
policy = PolicyNet()
train_bc(policy, (expert_states, expert_actions), epochs=50)

# --- DAgger loop ---
def dagger_training(num_iters=5, T=200):
    policy = PolicyNet()
    # Initial dataset from expert-only rollouts
    states, actions = collect_expert_trajectory(T)
    for i in range(num_iters):
        # Train policy on current dataset
        train_bc(policy, (states, actions), epochs=20)
        # Rollout current policy to explore its own state distribution
        roll_states = rollout(policy, T)
        # Query expert at visited states
        new_actions = np.array([expert_policy(s) for s in roll_states], dtype=np.float32)
        # Aggregate dataset
        states = np.concatenate([states, roll_states], axis=0)
        actions = np.concatenate([actions, new_actions], axis=0)
        print(f"DAgger iter {i+1}, dataset size = {states.shape[0]}")
    return policy, (states, actions)

dagger_policy, agg_dataset = dagger_training()
      

In a realistic robotics setting, simulate_dynamics would interface with a physics engine (e.g., PyBullet, MuJoCo, Gazebo) and expert_policy could be a teleoperation input or a pre-existing controller. The DAgger loop remains structurally identical.

In C++, behavior cloning and DAgger can be implemented in conjunction with ROS 2 and a simulation backend (e.g., Gazebo, Ignition) using Eigen for linear algebra. Below is a simplified sketch showing the basic structure of behavior cloning and data aggregation; actual code would integrate with ROS topics and services.

#include <iostream>
#include <vector>
#include <Eigen/Dense>

using State = Eigen::Vector4d;   // [q1, q2, dq1, dq2]
using Action = Eigen::Vector2d;  // [tau1, tau2]

struct Sample {
    State s;
    Action a;
};

class LinearPolicy {
public:
    Eigen::Matrix<double, 2, 4> W;
    Action b;

    LinearPolicy() {
        W.setZero();
        b.setZero();
    }

    Action operator()(const State &s) const {
        return W * s + b;
    }
};

// Dummy expert PD controller
Action expertPolicy(const State &s) {
    Eigen::Vector2d q_des(0.5, -0.5);
    Eigen::Vector2d q = s.segment<2>(0);
    Eigen::Vector2d dq = s.segment<2>(2);
    double Kp = 5.0, Kd = 1.0;
    return Kp * (q_des - q) - Kd * dq;
}

// Simple discrete-time dynamics for illustration
State simulate(const State &s, const Action &a, double dt = 0.02) {
    Eigen::Vector2d q = s.segment<2>(0);
    Eigen::Vector2d dq = s.segment<2>(2);
    Eigen::Vector2d ddq = a; // placeholder
    Eigen::Vector2d dq_next = dq + ddq * dt;
    Eigen::Vector2d q_next = q + dq_next * dt;
    State s_next;
    s_next << q_next, dq_next;
    return s_next;
}

// Least-squares fit of linear policy: minimize sum ||W s + b - a||^2
void trainBC(LinearPolicy &pi, const std::vector<Sample> &data, int iters = 100, double lr = 1e-3) {
    for (int k = 0; k < iters; ++k) {
        Eigen::Matrix<double, 2, 4> gradW = Eigen::Matrix<double, 2, 4>::Zero();
        Action gradb = Action::Zero();
        for (const auto &sample : data) {
            Action pred = pi(sample.s);
            Action err = pred - sample.a;
            gradW += err * sample.s.transpose();
            gradb += err;
        }
        gradW /= static_cast<double>(data.size());
        gradb /= static_cast<double>(data.size());
        pi.W -= lr * gradW;
        pi.b -= lr * gradb;
    }
}

std::vector<Sample> collectExpertRollout(int T) {
    std::vector<Sample> traj;
    State s = State::Zero();
    for (int t = 0; t < T; ++t) {
        Action a = expertPolicy(s);
        traj.push_back({s, a});
        s = simulate(s, a);
    }
    return traj;
}

// DAgger training loop
void daggerTrain(int numIters = 5, int T = 200) {
    LinearPolicy pi;
    std::vector<Sample> dataset = collectExpertRollout(T);

    for (int i = 0; i < numIters; ++i) {
        trainBC(pi, dataset, 200, 1e-3);

        // Rollout current policy
        State s = State::Zero();
        std::vector<Sample> newSamples;
        for (int t = 0; t < T; ++t) {
            Action a_pi = pi(s);
            // Query expert at visited state s
            Action a_exp = expertPolicy(s);
            newSamples.push_back({s, a_exp});
            s = simulate(s, a_pi);
        }

        // Aggregate
        dataset.insert(dataset.end(), newSamples.begin(), newSamples.end());
        std::cout << "DAgger iter " << (i+1)
                  << ", dataset size = " << dataset.size() << std::endl;
    }
}

int main() {
    daggerTrain();
    return 0;
}
      

In ROS-based systems, simulate and expertPolicy would interact with a robot simulation (or physical hardware), while states and actions would be passed through ROS message types. The learning logic remains identical.

Java is often used for robotics in contexts such as FIRST Robotics (WPILib), mobile robotics frameworks, or simulator front-ends. Below is a minimal Java sketch using simple arrays and gradient descent for behavior cloning and DAgger.

import java.util.ArrayList;
import java.util.List;

class Sample {
    double[] s; // length 4
    double[] a; // length 2
    Sample(double[] s, double[] a) {
        this.s = s;
        this.a = a;
    }
}

class LinearPolicy {
    // a = W s + b, W: 2x4, b: 2
    double[][] W = new double[2][4];
    double[] b = new double[2];

    double[] act(double[] s) {
        double[] a = new double[2];
        for (int i = 0; i < 2; i++) {
            a[i] = b[i];
            for (int j = 0; j < 4; j++) {
                a[i] += W[i][j] * s[j];
            }
        }
        return a;
    }
}

public class DAggerDemo {

    static double[] expertPolicy(double[] s) {
        double[] q_des = {0.5, -0.5};
        double[] q = {s[0], s[1]};
        double[] dq = {s[2], s[3]};
        double Kp = 5.0, Kd = 1.0;
        double[] tau = new double[2];
        for (int i = 0; i < 2; i++) {
            tau[i] = Kp * (q_des[i] - q[i]) - Kd * dq[i];
        }
        return tau;
    }

    static double[] simulate(double[] s, double[] a, double dt) {
        double q1 = s[0], q2 = s[1];
        double dq1 = s[2], dq2 = s[3];
        double ddq1 = a[0], ddq2 = a[1];
        double dq1n = dq1 + ddq1 * dt;
        double dq2n = dq2 + ddq2 * dt;
        double q1n = q1 + dq1n * dt;
        double q2n = q2 + dq2n * dt;
        return new double[]{q1n, q2n, dq1n, dq2n};
    }

    static List<Sample> collectExpertRollout(int T) {
        List<Sample> traj = new ArrayList<>();
        double[] s = {0, 0, 0, 0};
        for (int t = 0; t < T; t++) {
            double[] a = expertPolicy(s);
            traj.add(new Sample(s.clone(), a.clone()));
            s = simulate(s, a, 0.02);
        }
        return traj;
    }

    static void trainBC(LinearPolicy pi, List<Sample> data, int iters, double lr) {
        for (int it = 0; it < iters; it++) {
            double[][] gradW = new double[2][4];
            double[] gradb = new double[2];
            for (Sample sample : data) {
                double[] pred = pi.act(sample.s);
                double[] err = {pred[0] - sample.a[0], pred[1] - sample.a[1]};
                for (int i = 0; i < 2; i++) {
                    gradb[i] += err[i];
                    for (int j = 0; j < 4; j++) {
                        gradW[i][j] += err[i] * sample.s[j];
                    }
                }
            }
            int N = data.size();
            for (int i = 0; i < 2; i++) {
                gradb[i] /= N;
                for (int j = 0; j < 4; j++) gradW[i][j] /= N;
            }
            for (int i = 0; i < 2; i++) {
                pi.b[i] -= lr * gradb[i];
                for (int j = 0; j < 4; j++) {
                    pi.W[i][j] -= lr * gradW[i][j];
                }
            }
        }
    }

    static void daggerTrain() {
        int numIters = 5;
        int T = 200;
        LinearPolicy pi = new LinearPolicy();
        List<Sample> dataset = collectExpertRollout(T);

        for (int i = 0; i < numIters; i++) {
            trainBC(pi, dataset, 200, 1e-3);
            List<Sample> newSamples = new ArrayList<>();
            double[] s = {0, 0, 0, 0};
            for (int t = 0; t < T; t++) {
                double[] a_pi = pi.act(s);
                double[] a_exp = expertPolicy(s);
                newSamples.add(new Sample(s.clone(), a_exp.clone()));
                s = simulate(s, a_pi, 0.02);
            }
            dataset.addAll(newSamples);
            System.out.println("DAgger iter " + (i + 1) + ", dataset size = " + dataset.size());
        }
    }

    public static void main(String[] args) {
        daggerTrain();
    }
}
      

In a robotics codebase, the state vector s would typically include joint positions and velocities read from sensors, while actions correspond to torque or position commands sent to the robot controller.

MATLAB provides convenient tools for matrix-based policy representation and integration with Simulink models of robot dynamics. We sketch code that collects expert trajectories from a Simulink model, performs behavior cloning via linear regression, and describes how to embed the learned policy back into Simulink.

% State: [q1; q2; dq1; dq2], Action: [tau1; tau2]
STATE_DIM = 4;
ACTION_DIM = 2;

function tau = expert_policy(s)
    q_des = [0.5; -0.5];
    q = s(1:2);
    dq = s(3:4);
    Kp = 5.0; Kd = 1.0;
    tau = Kp * (q_des - q) - Kd * dq;
end

function s_next = simulate_step(s, tau, dt)
    q = s(1:2);
    dq = s(3:4);
    ddq = tau; % placeholder, use rigidBodyTree dynamics in practice
    dq_next = dq + ddq * dt;
    q_next = q + dq_next * dt;
    s_next = [q_next; dq_next];
end

function [S, A] = collect_expert_rollout(T)
    s = zeros(STATE_DIM, 1);
    S = zeros(STATE_DIM, T);
    A = zeros(ACTION_DIM, T);
    dt = 0.02;
    for t = 1:T
        tau = expert_policy(s);
        S(:, t) = s;
        A(:, t) = tau;
        s = simulate_step(s, tau, dt);
    end
end

% Behavior cloning via linear least squares: A ≈ W*S + b
T = 200;
[S, A] = collect_expert_rollout(T);
X = [S' ones(T,1)];          % [s', 1]
Y = A';                      % actions
theta = (X' * X) \ (X' * Y); % (STATE_DIM+1) x ACTION_DIM
W = theta(1:STATE_DIM, :)';  % 2 x 4
b = theta(STATE_DIM+1, :)';  % 2 x 1

% DAgger iterations
numIters = 5;
for k = 1:numIters
    % Rollout current policy
    dt = 0.02;
    s = zeros(STATE_DIM, 1);
    S_new = zeros(STATE_DIM, T);
    A_new = zeros(ACTION_DIM, T);
    for t = 1:T
        tau_pi = W * s + b;       % current policy
        tau_exp = expert_policy(s); % query expert at visited state
        S_new(:, t) = s;
        A_new(:, t) = tau_exp;
        s = simulate_step(s, tau_pi, dt);
    end
    % Aggregate and retrain
    S = [S, S_new];
    A = [A, A_new];
    X = [S' ones(size(S,2),1)];
    Y = A';
    theta = (X' * X) \ (X' * Y);
    W = theta(1:STATE_DIM, :)';
    b = theta(STATE_DIM+1, :)';
    fprintf('DAgger iter %d, dataset size = %d\n', k, size(S,2));
end

% To deploy in Simulink:
% 1. Create a MATLAB Function block implementing tau = W*s + b.
% 2. Feed joint states into this block and connect outputs to your torque inputs.
      

In a more detailed setup, the manipulator dynamics would be modeled as a rigidBodyTree in MATLAB with a Simulink block diagram, and the learned policy would be implemented as a controller subsystem driven by measured state signals.

Wolfram Mathematica offers symbolic and numeric tools useful for analyzing error propagation and fitting policies. Below is a sketch of behavior cloning and one DAgger iteration using linear regression constructs.

(* State and expert policy *)
stateDim = 4;
actionDim = 2;

expertPolicy[s_] := Module[
  {qDes, q, dq, Kp = 5.0, Kd = 1.0},
  qDes = {0.5, -0.5};
  q = s[[1 ;; 2]];
  dq = s[[3 ;; 4]];
  Kp (qDes - q) - Kd dq
]

simulateStep[s_, tau_, dt_:0.02] := Module[
  {q, dq, ddq, dqNext, qNext},
  q = s[[1 ;; 2]];
  dq = s[[3 ;; 4]];
  ddq = tau;
  dqNext = dq + ddq dt;
  qNext = q + dqNext dt;
  Join[qNext, dqNext]
]

collectExpertRollout[T_] := Module[
  {s, states, actions, tau},
  s = ConstantArray[0.0, stateDim];
  states = {}; actions = {};
  Do[
    tau = expertPolicy[s];
    AppendTo[states, s];
    AppendTo[actions, tau];
    s = simulateStep[s, tau],
    {t, 1, T}
  ];
  {states, actions}
]

{states, actions} = collectExpertRollout[200];

(* Behavior cloning via linear regression: actions ≈ W.states + b *)
X = Map[Append[#, 1.0] &, states];
Y = actions;
theta = LinearModelFit[
   Thread[X -> Y],
   Table[x[i], {i, 1, stateDim + 1}],
   Table[x[i], {i, 1, stateDim + 1}]
]["BestFitParameters"];

W = ArrayReshape[theta[[1 ;; stateDim*actionDim]], {actionDim, stateDim}];
b = theta[[stateDim*actionDim + 1 ;; stateDim*actionDim + actionDim]];

policy[s_] := W.s + b;

rollout[pi_, T_] := Module[
  {s, traj = {}, a},
  s = ConstantArray[0.0, stateDim];
  Do[
    a = pi[s];
    AppendTo[traj, s];
    s = simulateStep[s, a],
    {t, 1, T}
  ];
  traj
];

(* One DAgger iteration *)
T = 200;
trajStates = rollout[policy, T];
newActions = expertPolicy /@ trajStates;
statesAgg = Join[states, trajStates];
actionsAgg = Join[actions, newActions];
(* Retrain policy on aggregated data as above *)
      

Mathematica also allows symbolic derivation of error bounds and visualization of state distributions under different policies, which can be useful to illustrate covariate shift in a lecture setting.

Problem 1 (Behavior Cloning Objective): Consider a finite-horizon system with horizon \( T \) and expert policy \( \pi_E \). Show that if we assume the state distribution under the learned policy \( \pi_\theta \) equals the expert state distribution at each time step, i.e., \( d_t^{\pi_\theta} = d_t^{\pi_E} \) for all \( t \), then minimizing the behavior cloning loss \( \mathbb{E}_{(s,a)\sim\mathcal{D}_E}[\ell(\pi_\theta(s),a)] \) is equivalent to minimizing the task cost \( J(\pi_\theta) \) under a surrogate per-step loss.

Solution:

Under the assumption \( d_t^{\pi_\theta} = d_t^{\pi_E} \), we have

\[ J(\pi_\theta) = \sum_{t=1}^{T} \mathbb{E}_{s_t \sim d_t^{\pi_E}} \big[ c(s_t, \pi_\theta(s_t)) \big]. \]

If we choose the supervised loss \( \ell \) such that \( \ell(\pi_\theta(s), \pi_E(s)) \approx c(s, \pi_\theta(s)) - c(s, \pi_E(s)) \), then

\[ \mathbb{E}_{(s,a)\sim\mathcal{D}_E}[\ell(\pi_\theta(s),a)] \approx \frac{1}{T} \left( J(\pi_\theta) - J(\pi_E) \right), \]

since \( \mathcal{D}_E \) is formed from the union of expert states across time and expert actions \( a = \pi_E(s) \). Thus minimizing the behavior cloning objective approximately minimizes the difference in task cost between \( \pi_\theta \) and the expert. The approximation becomes exact if the surrogate loss is chosen to match cost increments exactly.

Problem 2 (Compounding Error Bound): Assume that behavior cloning yields a policy \( \pi_\theta \) with per-state error at most \( \epsilon \) under the expert distribution, i.e., \( \Pr(\pi_\theta(s) \neq \pi_E(s)) \leq \epsilon \) for \( s \sim d_t^{\pi_E} \) at each time step. Assume that a deviation at any time can cause worst-case loss 1 at all future time steps. Derive an upper bound on \( J(\pi_\theta) - J(\pi_E) \).

Solution:

Let \( E_t \) be the event that a deviation has occurred on or before time \( t \). Using a union bound over time steps,

\[ \Pr(E_t) \leq \sum_{k=1}^{t} \Pr(\text{deviation at time } k) \leq t \epsilon. \]

Denote the instantaneous cost under policy \( \pi \) as \( c_t(\pi) \). For each time step, decompose

\[ \mathbb{E}[c_t(\pi_\theta)] = \mathbb{E}[c_t(\pi_\theta) \mid E_t] \Pr(E_t) + \mathbb{E}[c_t(\pi_\theta) \mid \neg E_t] \Pr(\neg E_t). \]

If no deviation has occurred, \( \pi_\theta \) matches the expert and \( c_t(\pi_\theta) = c_t(\pi_E) \). If at least one deviation has occurred, the worst-case loss is bounded by 1. Therefore

\[ \mathbb{E}[c_t(\pi_\theta)] \leq \mathbb{E}[c_t(\pi_E)] + \Pr(E_t). \]

Summing over \( t = 1, \dots, T \) yields

\[ J(\pi_\theta) - J(\pi_E) \leq \sum_{t=1}^{T} \Pr(E_t) \leq \sum_{t=1}^{T} t \epsilon = \epsilon \frac{T(T+1)}{2} = \mathcal{O}(T^2 \epsilon). \]

This shows the quadratic dependence on horizon inherent in naïve behavior cloning.

Problem 3 (DAgger Loss Definition): Let \( \pi_1, \dots, \pi_N \) be the sequence of policies produced by DAgger. For each iteration, define the loss \( \ell_i(\pi) \) as in Section 4. Show that the empirical average of these losses approximates the average supervised error of \( \pi \) on the aggregated dataset.

Solution:

At iteration \( i \), DAgger produces a batch of samples \( \{(s_t^{(i)}, a_t^{E,(i)})\}_t \) where \( s_t^{(i)} \) is drawn from \( d_t^{\pi_i} \) and \( a_t^{E,(i)} = \pi_E(s_t^{(i)}) \). Using an empirical average,

\[ \widehat{\ell}_i(\pi) = \frac{1}{T} \sum_{t=1}^{T} \ell\big(\pi(s_t^{(i)}), a_t^{E,(i)}\big) \approx \mathbb{E}_{s \sim d^{\pi_i}} \big[ \ell(\pi(s), \pi_E(s)) \big] = \ell_i(\pi). \]

The aggregated dataset after \( N \) iterations is the union \( \bigcup_{i=1}^{N} \{(s_t^{(i)}, a_t^{E,(i)})\}_t \), so the empirical supervised error of \( \pi \) on this dataset is precisely the average of \( \widehat{\ell}_i(\pi) \) over \( i \), up to sampling noise. Thus optimizing a policy on the aggregated dataset approximates minimizing the average DAgger losses.

Problem 4 (No-Regret Condition): Suppose an online algorithm over the policy class \( \Pi \) satisfies for all \( \pi \in \Pi \)

\[ \sum_{i=1}^{N} \ell_i(\pi_i) - \sum_{i=1}^{N} \ell_i(\pi) \leq R_N, \]

with \( R_N / N \rightarrow 0 \). Show that there exists some iteration \( i \) such that

\[ \ell_i(\pi_i) \leq \inf_{\pi \in \Pi} \frac{1}{N} \sum_{j=1}^{N} \ell_j(\pi) + \frac{R_N}{N}. \]

Solution:

Divide the regret inequality by \( N \):

\[ \frac{1}{N} \sum_{i=1}^{N} \ell_i(\pi_i) \leq \inf_{\pi \in \Pi} \frac{1}{N} \sum_{i=1}^{N} \ell_i(\pi) + \frac{R_N}{N}. \]

If for all \( i \) we had \( \ell_i(\pi_i) > \inf_{\pi \in \Pi} \frac{1}{N} \sum_{j=1}^{N} \ell_j(\pi) + \frac{R_N}{N} \), then their average would also exceed the right-hand side, contradicting the inequality. Thus at least one \( i \) must satisfy the desired bound. This is the classic averaging argument used in online learning analysis.

In this lesson we formalized behavior cloning as supervised learning of a policy from expert demonstrations and showed how covariate shift leads to \( \mathcal{O}(T^2 \epsilon) \) compounding error. We then introduced Dataset Aggregation (DAgger), which iteratively collects expert labels on the states induced by the learner's own policy. Under a no-regret online learning framework, DAgger enjoys an improved \( \mathcal{O}(T \epsilon_N) \) performance guarantee.

We also illustrated concrete implementations in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica for simple manipulator-like systems, providing a bridge from theoretical guarantees to practical robotic control. These methods form the foundation for more advanced imitation and reinforcement learning techniques studied in subsequent lessons.