We consider a single-arm manipulation task (e.g., reaching to a pose, pushing an object). Let \( \mathbf{q}_t \in \mathbb{R}^n \) denote joint angles and \( \dot{\mathbf{q}}_t \) joint velocities at discrete time \( t \), and let the commanded action be joint velocities or torques \( \mathbf{a}_t \in \mathbb{R}^m \) (often \( m = n \)). A typical state representation for imitation is

\[ \mathbf{s}_t = \begin{bmatrix} \mathbf{q}_t \\ \dot{\mathbf{q}}_t \\ \mathbf{x}_t^{\text{goal}} - \mathbf{x}_t \end{bmatrix} \in \mathbb{R}^d, \quad \mathbf{a}_t = \boldsymbol{\tau}_t \ \text{or}\ \dot{\mathbf{q}}_t^{\text{cmd}}. \]

The lab pipeline is:

flowchart TD
  T["Choose manipulation task"] --> D["Collect expert demos (teleop / kinesthetic)"]
  D --> DS["Build dataset {(s_t, a_t)}"]
  DS --> PC["Select policy class pi_theta"]
  PC --> TR["Train policy on dataset"]
  TR --> EV["Evaluate in simulator or real robot"]
  EV --> IT["Iterate: tune features, add demos, refine policy"]
        

Mathematically, we will define a parameterized policy \( \pi_{\boldsymbol{\theta}} : \mathbb{R}^d \rightarrow \mathbb{R}^m \), a loss function \( \ell(\cdot,\cdot) \), and minimize empirical imitation loss over the demonstration dataset.

Suppose an expert provides \( K \) demonstrations, each a trajectory of length \( T_k \): \( \tau^{(k)} = (\mathbf{s}_1^{(k)}, \mathbf{a}_1^{(k)}, \dots, \mathbf{s}_{T_k}^{(k)}, \mathbf{a}_{T_k}^{(k)}) \). The full dataset is the union \( \mathcal{D} = \{(\mathbf{s}_n, \mathbf{a}_n)\}_{n=1}^N \), where \( N = \sum_{k=1}^K T_k \). For continuous actions (joint velocities or torques), a natural choice is squared error loss:

\[ \ell(\hat{\mathbf{a}}, \mathbf{a}) = \lVert \hat{\mathbf{a}} - \mathbf{a} \rVert_2^2 = (\hat{\mathbf{a}} - \mathbf{a})^{\top}(\hat{\mathbf{a}} - \mathbf{a}). \]

The empirical behavior cloning objective is then

\[ \mathcal{L}(\boldsymbol{\theta}) = \frac{1}{N} \sum_{n=1}^{N} \ell\bigl(\pi_{\boldsymbol{\theta}}(\mathbf{s}_n), \mathbf{a}_n\bigr) + \lambda \lVert \boldsymbol{\theta} \rVert_2^2, \]

where \( \lambda \ge 0 \) is a regularization parameter that penalizes overly large parameters and helps control overfitting.

Conceptually, behavior cloning assumes a supervised learning viewpoint: there exists an (unknown) expert mapping \( \pi^{\star}(\mathbf{s}) \) such that \( \mathbf{a}_n \approx \pi^{\star}(\mathbf{s}_n) \). We approximate \( \pi^{\star} \) inside some policy class (e.g., linear, neural network) using empirical risk minimization.

A useful baseline is a linear policy in features \( \boldsymbol{\phi}(\mathbf{s}) \in \mathbb{R}^d \):

\[ \pi_{\boldsymbol{\theta}}(\mathbf{s}) = \mathbf{W}\,\boldsymbol{\phi}(\mathbf{s}) + \mathbf{b}, \quad \mathbf{W} \in \mathbb{R}^{m \times d},\ \mathbf{b} \in \mathbb{R}^m, \]

with parameters \( \boldsymbol{\theta} = (\mathbf{W}, \mathbf{b}) \). Define augmented features

\[ \tilde{\boldsymbol{\phi}}(\mathbf{s}) = \begin{bmatrix}\boldsymbol{\phi}(\mathbf{s}) \\ 1\end{bmatrix} \in \mathbb{R}^{d+1}, \quad \tilde{\mathbf{W}} = \begin{bmatrix} \mathbf{W} & \mathbf{b} \end{bmatrix} \in \mathbb{R}^{m \times (d+1)}. \]

Stack the data into matrices \( \tilde{\mathbf{\Phi}} \in \mathbb{R}^{N \times (d+1)} \) and \( \mathbf{A} \in \mathbb{R}^{N \times m} \):

\[ \tilde{\mathbf{\Phi}} = \begin{bmatrix} \tilde{\boldsymbol{\phi}}(\mathbf{s}_1)^{\top} \\ \vdots \\ \tilde{\boldsymbol{\phi}}(\mathbf{s}_N)^{\top} \end{bmatrix}, \quad \mathbf{A} = \begin{bmatrix} \mathbf{a}_1^{\top} \\ \vdots \\ \mathbf{a}_N^{\top} \end{bmatrix}. \]

The regularized empirical loss can be written as a Frobenius norm:

\[ \mathcal{L}(\tilde{\mathbf{W}}) = \frac{1}{N}\bigl\lVert \tilde{\mathbf{\Phi}}\tilde{\mathbf{W}} - \mathbf{A} \bigr\rVert_F^2 + \lambda \lVert \tilde{\mathbf{W}} \rVert_F^2. \]

Taking the derivative with respect to \( \tilde{\mathbf{W}} \) and setting to zero:

\[ \frac{\partial \mathcal{L}}{\partial \tilde{\mathbf{W}}} = \frac{2}{N}\tilde{\mathbf{\Phi}}^{\top} (\tilde{\mathbf{\Phi}}\tilde{\mathbf{W}} - \mathbf{A}) + 2\lambda \tilde{\mathbf{W}} = \mathbf{0}. \]

Rearranging,

\[ \bigl(\tilde{\mathbf{\Phi}}^{\top}\tilde{\mathbf{\Phi}} + N\lambda \mathbf{I}\bigr)\tilde{\mathbf{W}} = \tilde{\mathbf{\Phi}}^{\top}\mathbf{A}, \]

and assuming \( \tilde{\mathbf{\Phi}}^{\top}\tilde{\mathbf{\Phi}} + N\lambda \mathbf{I} \) is invertible, the unique minimizer is

\[ \tilde{\mathbf{W}}^{\star} = \bigl(\tilde{\mathbf{\Phi}}^{\top}\tilde{\mathbf{\Phi}} + N\lambda \mathbf{I}\bigr)^{-1} \tilde{\mathbf{\Phi}}^{\top}\mathbf{A}. \]

This is just multivariate ridge regression, applied per joint dimension. It provides a non-iterative baseline for the lab, and is particularly convenient for C++/Java/MATLAB implementations using linear algebra libraries.

Assume a conditional Gaussian model for expert actions:

\[ p(\mathbf{a}\mid\mathbf{s};\boldsymbol{\theta}) = \mathcal{N}\bigl(\mathbf{a} \,\big|\, \pi_{\boldsymbol{\theta}}(\mathbf{s}),\ \sigma^2\mathbf{I}\bigr). \]

Then the log-likelihood of the dataset is

\[ \log p(\mathcal{D}\mid\boldsymbol{\theta}) = -\frac{1}{2\sigma^2}\sum_{n=1}^{N} \bigl\lVert \mathbf{a}_n - \pi_{\boldsymbol{\theta}}(\mathbf{s}_n) \bigr\rVert_2^2 + C, \]

where \( C \) does not depend on \( \boldsymbol{\theta} \). Maximizing this likelihood is equivalent to minimizing the squared error loss, so behavior cloning with MSE can be interpreted as maximum-likelihood estimation under Gaussian action noise.

A key subtlety is the state distribution. Let \( d_E(\mathbf{s}) \) denote the state distribution induced by the expert, and \( d_{\pi}(\mathbf{s}) \) the distribution induced when we execute the learned policy on the system. Behavior cloning minimizes

\[ \mathbb{E}_{\mathbf{s}\sim d_E} \bigl[\ell(\pi_{\boldsymbol{\theta}}(\mathbf{s}), \pi^{\star}(\mathbf{s}))\bigr], \]

but performance during execution depends on \( \mathbf{s} \sim d_{\pi} \), which can drift away from the expert distribution if small errors accumulate. Later in the chapter, dataset aggregation methods such as DAgger address this covariate shift; here we simply note the limitation for the lab experiments.

We assume a simulator or real robot interface providing:

• reset(): places the system in an initial state.
• get_state(): returns the current state \( \mathbf{s}_t \).
• apply_action(a): sends control command \( \mathbf{a}_t \).

For each demonstration episode, an expert (teleoperation, kinesthetic teaching, or scripted controller) provides an action a_demo for each observed state. The training loop alternates between sampling mini-batches of \( (\mathbf{s}_n,\mathbf{a}_n) \) and updating \( \boldsymbol{\theta} \) via gradient descent or the closed-form solution.

flowchart TD
  S["Start episode"] --> OBS["Observe state s_t"]
  OBS --> EXP["Expert chooses action a_t"]
  EXP --> REC["Record (s_t, a_t) in dataset"]
  REC --> NEXT["Step system, t := t+1"]
  NEXT -->|more steps| OBS
  NEXT -->|end| TRN["Train pi_theta on dataset"]
  TRN --> EVAL["Roll out pi_theta without expert"]
  EVAL --> UPD["Analyze errors, collect more demos if needed"]
        

In the concrete code below, we focus on the supervised training part (from dataset to policy) while assuming the demonstration collection is already done or provided as a separate tool.

We implement a small fully connected network that maps the state vector \( \mathbf{s}_t \) to joint velocity commands \( \dot{\mathbf{q}}_t^{\text{cmd}} \). The code uses numpy for data handling and torch for optimization. A typical robotics stack would embed this into a pybullet, MuJoCo, or ROS-based simulation node.

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

# -------------------------------------------------------
# 1. Load demonstrations: states S (N x d), actions A (N x m)
# -------------------------------------------------------
data = np.load("manip_demos.npz")
S = data["states"].astype(np.float32)   # shape (N, d)
A = data["actions"].astype(np.float32)  # shape (N, m)

N, d = S.shape
m = A.shape[1]

device = torch.device("cuda" if torch.cuda.is_available() else "cpu")

states = torch.from_numpy(S).to(device)
actions = torch.from_numpy(A).to(device)

# -------------------------------------------------------
# 2. Define neural policy pi_theta(s)
# -------------------------------------------------------
class MLPPolicy(nn.Module):
    def __init__(self, state_dim, action_dim, hidden_sizes=(128, 128)):
        super().__init__()
        layers = []
        last_dim = state_dim
        for h in hidden_sizes:
            layers.append(nn.Linear(last_dim, h))
            layers.append(nn.ReLU())
            last_dim = h
        layers.append(nn.Linear(last_dim, action_dim))
        self.net = nn.Sequential(*layers)

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

policy = MLPPolicy(d, m).to(device)

# -------------------------------------------------------
# 3. Optimization setup
# -------------------------------------------------------
optimizer = optim.Adam(policy.parameters(), lr=1e-3, weight_decay=1e-4)
loss_fn = nn.MSELoss()

batch_size = 256
num_epochs = 50

# -------------------------------------------------------
# 4. Training loop
# -------------------------------------------------------
perm = torch.randperm(N)
states = states[perm]
actions = actions[perm]

for epoch in range(num_epochs):
    epoch_loss = 0.0
    for start in range(0, N, batch_size):
        end = min(start + batch_size, N)
        batch_s = states[start:end]
        batch_a = actions[start:end]

        pred_a = policy(batch_s)
        loss = loss_fn(pred_a, batch_a)

        optimizer.zero_grad()
        loss.backward()
        optimizer.step()

        epoch_loss += loss.item() * (end - start)

    epoch_loss /= N
    print(f"Epoch {epoch+1:03d} | loss = {epoch_loss:.6f}")

# -------------------------------------------------------
# 5. Example rollout in a robotics environment (pseudo-code)
# -------------------------------------------------------
# import pybullet as p
# env = MyManipEnv()
# state = env.reset()
# for t in range(max_steps):
#     s_vec = torch.from_numpy(state.astype(np.float32)).unsqueeze(0).to(device)
#     with torch.no_grad():
#         a_vec = policy(s_vec).cpu().numpy()[0]
#     state, obs = env.apply_velocity_command(a_vec)
#     env.render()
      

In a ROS-based setup, the rollout loop would run inside a ROS node, subscribing to joint state topics and publishing velocity commands. The learned policy is purely kinematic/dynamic interface agnostic: only the state vector definition and low-level controller wrapper are platform-specific.

For real-time control on embedded controllers or within ROS, it is often convenient to use a linear policy with pre-computed \( \tilde{\mathbf{W}}^{\star} \). We demonstrate how to compute the ridge-regression solution using the Eigen library, assuming that feature and action matrices are loaded from disk.

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

// Load Phi (N x (d+1)) and A (N x m) from plain-text files
bool loadMatrix(const std::string& path, Eigen::MatrixXf& M, int rows, int cols) {
    std::ifstream in(path);
    if (!in) return false;
    M.resize(rows, cols);
    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < cols; ++j) {
            in >> M(i, j);
        }
    }
    return true;
}

int main() {
    int N = 10000;     // number of samples
    int d1 = 33;       // feature dimension + 1 for bias
    int m  = 7;        // number of joints

    Eigen::MatrixXf Phi_tilde;
    Eigen::MatrixXf A;
    if (!loadMatrix("Phi_tilde.txt", Phi_tilde, N, d1)) {
        std::cerr << "Could not load Phi_tilde\n";
        return 1;
    }
    if (!loadMatrix("A.txt", A, N, m)) {
        std::cerr << "Could not load A\n";
        return 1;
    }

    float lambda = 1e-3f;

    // Compute W_tilde* = (Phi^T Phi + N lambda I)^(-1) Phi^T A
    Eigen::MatrixXf I = Eigen::MatrixXf::Identity(d1, d1);
    Eigen::MatrixXf G = Phi_tilde.transpose() * Phi_tilde + N * lambda * I;
    Eigen::MatrixXf RHS = Phi_tilde.transpose() * A;

    Eigen::MatrixXf W_tilde = G.ldlt().solve(RHS);

    // Split into W (m x d) and b (m)
    int d = d1 - 1;
    Eigen::MatrixXf W = W_tilde.leftCols(d);
    Eigen::VectorXf b = W_tilde.col(d);

    // Example: evaluate policy for one state
    Eigen::VectorXf s(d);
    // ... fill s with features ...
    Eigen::VectorXf a_cmd = W * s + b;

    std::cout << "Commanded action:\n" << a_cmd.transpose() << std::endl;
    return 0;
}
      

In a ROS node, one would wrap the call W * s + b inside a subscriber callback that receives joint states and publishes either velocity or torque commands. Libraries like MoveIt can be used upstream to produce reference end-effector poses that are encoded in \( \boldsymbol{\phi}(\mathbf{s}) \).

Java is less common in manipulation, but with ROSJava and numerical libraries like EJML, it can be used for high-level control. The following code computes the same ridge-regression solution using EJML's SimpleMatrix API.

import org.ejml.simple.SimpleMatrix;

public class ImitationLinearPolicy {
    public static void main(String[] args) {
        int N = 10000;
        int d1 = 33;  // feature dim + 1
        int m  = 7;   // joints
        double lambda = 1e-3;

        // Assume these are loaded from files or constructed from data
        SimpleMatrix PhiTilde = loadMatrix("Phi_tilde.txt", N, d1);
        SimpleMatrix A = loadMatrix("A.txt", N, m);

        SimpleMatrix I = SimpleMatrix.identity(d1);
        SimpleMatrix G = PhiTilde.transpose().mult(PhiTilde)
                            .plus(I.scale(N * lambda));
        SimpleMatrix RHS = PhiTilde.transpose().mult(A);

        SimpleMatrix Wtilde = G.solve(RHS);

        int d = d1 - 1;
        SimpleMatrix W = Wtilde.rows(0, m).cols(0, d);
        SimpleMatrix b = Wtilde.cols(d, d1); // (m x 1)

        // Example evaluation for one state
        SimpleMatrix s = new SimpleMatrix(d, 1);
        // ... fill s with features ...
        SimpleMatrix aCmd = W.mult(s).plus(b);

        System.out.println("Action: " + aCmd.transpose());
    }

    private static SimpleMatrix loadMatrix(String path, int rows, int cols) {
        // Implement file parsing as needed
        SimpleMatrix M = new SimpleMatrix(rows, cols);
        // ...
        return M;
    }
}
      

In a Java-based robotics platform, aCmd would be translated into low-level controller messages. The same mathematical formulation as in Section 3 is reused across languages.

MATLAB is widely used in robotics for prototyping and is tightly integrated with Simulink for closed-loop simulation. We first compute the linear policy weights, then show how to integrate them into a Simulink model via a MATLAB Function block.

% Load features Phi_tilde (N x (d+1)) and actions A (N x m)
load("manip_demos.mat", "Phi_tilde", "A");  % variables: Phi_tilde, A

[N, d1] = size(Phi_tilde);
[~, m]  = size(A);
lambda  = 1e-3;

G   = Phi_tilde' * Phi_tilde + N * lambda * eye(d1);
RHS = Phi_tilde' * A;

W_tilde = G \ RHS;          % (d1 x m)
d = d1 - 1;
W = W_tilde(1:d, :);        % (d x m)
b = W_tilde(d1, :).';       % (m x 1)

% Save parameters for Simulink block
save("policy_params.mat", "W", "b");

% --- MATLAB Function block (inside Simulink) ---
% function a_cmd = imitation_policy(s)
% %#codegen
% load('policy_params.mat', 'W', 'b');
% a_cmd = W' * s + b;
      

In Simulink, the block imitation_policy takes a state vector \( \mathbf{s}_t \) (constructed from joint states and goal pose), and outputs joint velocity or torque commands, which are then passed to a joint-space controller block and ultimately to a simulated manipulator in the Robotics System Toolbox environment.

Mathematica is useful for symbolic inspection of the optimization problem and for rapid numerical experimentation. The snippet below sets up a linear policy and minimizes the empirical squared loss.

(* Assume data is imported: list of pairs {s_n, a_n} *)
data = Import["manip_demos.mx"];

states  = data[[All, 1]];  (* list of state vectors *)
actions = data[[All, 2]];  (* list of action vectors *)

d = Length[states[[1]]];
m = Length[actions[[1]]];
N = Length[states];
lambda = 10.^(-3);

(* Linear policy parameters: W (m x d), b (m) *)
W = Array[w, {m, d}];
b = Array[bv, m];

phi[s_] := s;  (* identity features; replace if needed *)

predict[s_] := W . phi[s] + b;

lossValue :=
  (1./N) *
    Sum[Norm[predict[states[[n]]] - actions[[n]]]^2, {n, 1, N}] +
  lambda * Total[Flatten[W]^2];

opt = NMinimize[lossValue, Flatten[{W, b}]];
Wopt = ArrayReshape[Flatten[{W}] /. opt[[2]], {m, d}];
bopt = Table[b[[i]] /. opt[[2]], {i, 1, m}];

(* Policy evaluation *)
imitationPolicy[s_] := Wopt . phi[s] + bopt;
      

Symbolic differentiation of lossValue with respect to W and b reproduces the matrix normal equation derived in Section 3, which is a useful consistency check between algebraic derivation and numerical implementation.

Problem 1 (Closed-form Linear Policy): Starting from the linear policy formulation in Section 3, derive the closed-form minimizer \( \tilde{\mathbf{W}}^{\star} \) of the regularized loss

\[ \mathcal{L}(\tilde{\mathbf{W}}) = \frac{1}{N}\bigl\lVert \tilde{\mathbf{\Phi}}\tilde{\mathbf{W}} - \mathbf{A} \bigr\rVert_F^2 + \lambda \lVert \tilde{\mathbf{W}} \rVert_F^2. \]

Solution: Expand the Frobenius norm as a trace:

\[ \bigl\lVert \tilde{\mathbf{\Phi}}\tilde{\mathbf{W}} - \mathbf{A} \bigr\rVert_F^2 = \operatorname{tr}\bigl( (\tilde{\mathbf{\Phi}}\tilde{\mathbf{W}} - \mathbf{A})^{\top} (\tilde{\mathbf{\Phi}}\tilde{\mathbf{W}} - \mathbf{A}) \bigr). \]

Use matrix calculus: \( \partial\lVert \mathbf{X}\mathbf{W} - \mathbf{Y} \rVert_F^2 / \partial \mathbf{W} = 2\mathbf{X}^{\top}(\mathbf{X}\mathbf{W} - \mathbf{Y}) \), and \( \partial \lVert \mathbf{W} \rVert_F^2 / \partial \mathbf{W} = 2\mathbf{W} \). Setting the gradient to zero yields

\[ \frac{2}{N}\tilde{\mathbf{\Phi}}^{\top} (\tilde{\mathbf{\Phi}}\tilde{\mathbf{W}} - \mathbf{A}) + 2\lambda \tilde{\mathbf{W}} = \mathbf{0}, \]

which simplifies to \( (\tilde{\mathbf{\Phi}}^{\top}\tilde{\mathbf{\Phi}} + N\lambda\mathbf{I})\tilde{\mathbf{W}} = \tilde{\mathbf{\Phi}}^{\top}\mathbf{A} \). Under invertibility, the unique solution is \( \tilde{\mathbf{W}}^{\star} = (\tilde{\mathbf{\Phi}}^{\top}\tilde{\mathbf{\Phi}} + N\lambda\mathbf{I})^{-1} \tilde{\mathbf{\Phi}}^{\top}\mathbf{A} \).

Problem 2 (Equivalence to Gaussian Maximum Likelihood): Assume the conditional Gaussian model from Section 4 with variance \( \sigma^2 \mathbf{I} \) and policy \( \pi_{\boldsymbol{\theta}}(\mathbf{s}) = \mathbf{W}\boldsymbol{\phi}(\mathbf{s}) + \mathbf{b} \). Show that maximizing the (unregularized) likelihood is equivalent to minimizing the empirical MSE loss.

Solution: The negative log-likelihood (up to additive constant) is

\[ -\log p(\mathcal{D}\mid\boldsymbol{\theta}) = \frac{1}{2\sigma^2}\sum_{n=1}^{N} \bigl\lVert \mathbf{a}_n - \pi_{\boldsymbol{\theta}}(\mathbf{s}_n) \bigr\rVert_2^2 + C'. \]

Since \( \sigma^2 \) and \( C' \) do not depend on \( \boldsymbol{\theta} \), minimizing the negative log-likelihood is equivalent to minimizing \( \sum_{n=1}^{N} \lVert \mathbf{a}_n - \pi_{\boldsymbol{\theta}}(\mathbf{s}_n) \rVert_2^2 \), i.e., the empirical MSE. Thus maximum-likelihood and least-squares behavior cloning coincide under this noise model.

Problem 3 (Feature Design for a 2-DOF Planar Arm): Consider a planar 2-DOF arm with joint angles \( \mathbf{q}_t = (q_{1,t}, q_{2,t})^{\top} \), joint velocities \( \dot{\mathbf{q}}_t \), and a goal end-effector position \( \mathbf{x}^{\text{goal}} \). Propose a feature vector \( \boldsymbol{\phi}(\mathbf{s}_t) \) suitable for behavior cloning of joint velocity commands, and justify your choice.

Solution: A simple and effective choice is

\[ \boldsymbol{\phi}(\mathbf{s}_t) = \begin{bmatrix} q_{1,t} \\ q_{2,t} \\ \dot{q}_{1,t} \\ \dot{q}_{2,t} \\ x^{\text{goal}}_x - x_t^{} \\ x^{\text{goal}}_y - y_t^{} \end{bmatrix}, \]

where \( (x_t^{}, y_t^{}) \) is the current end-effector position. Joint angles and velocities encode the configuration and motion of the arm; the Cartesian error \( \mathbf{x}^{\text{goal}} - \mathbf{x}_t \) encodes the task objective. These features are linear in the state variables and allow the linear policy to represent proportional-derivative-like behavior in joint or task space, consistent with classical control intuition.

Problem 4 (Effect of Regularization on Imitation Error): For the linear policy with loss \( \mathcal{L}(\tilde{\mathbf{W}}) \), discuss qualitatively how increasing \( \lambda \) affects: (i) the variance of the learned parameters, (ii) the training error, and (iii) the expected imitation error on new states drawn from \( d_E \).

Solution: (i) Larger \( \lambda \) shrinks \( \tilde{\mathbf{W}}^{\star} \) toward zero, reducing parameter variance across different training sets (ridge regression shrinks coefficients). (ii) Training error (empirical MSE) typically increases with \( \lambda \), since the solution is more constrained. (iii) Expected imitation error on unseen states may first decrease (due to reduced overfitting) and then increase (due to underfitting) as \( \lambda \) grows, producing the usual U-shaped bias-variance trade-off curve.

Problem 5 (Covariate Shift in a Reaching Task): In a planar reaching task, the expert never allows the arm to enter a near-singular configuration, but the learned policy makes slightly larger errors, eventually driving the system into that region. Explain in terms of \( d_E \) versus \( d_{\pi} \) why this can cause catastrophic failure, even if the supervised imitation loss on \( \mathcal{D} \) is small.

Solution: The supervised loss is computed on samples from \( d_E \), which by construction assigns negligible probability to near-singular states. The policy may thus behave poorly on such states without being penalized during training. During execution, small discrepancies cause the closed-loop trajectory distribution \( d_{\pi} \) to drift; once it reaches the near-singular region, the policy is effectively extrapolating and may issue large or unstable commands. This mismatch between training and execution distributions is covariate shift and is particularly severe in systems with unstable or highly nonlinear dynamics, such as manipulators near kinematic singularities.

In this lab you have:

• Formulated behavior cloning for manipulation as a supervised learning problem over state-action pairs derived from demonstrations.
• Derived the closed-form ridge-regression solution for linear imitation policies and connected it to maximum-likelihood estimation under Gaussian action noise.
• Implemented imitation policies in Python (neural network), C++ and Java (linear policies with Eigen/EJML), MATLAB/Simulink (for closed-loop simulation), and Wolfram Mathematica (for symbolic and numerical optimization).
• Analyzed regularization, feature design, and covariate shift effects that influence the success of imitation in manipulation tasks.

These tools form the foundation for the more advanced imitation and inverse reinforcement learning methods developed later in the course, and for integrating learned policies into full robotics stacks that include motion planning and perception.